API - Linear Systems

backward_substitution(A::BallMatrix, v::BallVector)
backward_substitution(A::BallMatrix, B::BallMatrix)

Enclose the solution of A x = v, or of A X = B column by column, for an upper triangular A: the rows are eliminated from the last to the first, each step a ball subtraction and one ball division by the diagonal entry, so the result encloses the solution of every real or complex system whose matrix and right hand side lie in A and in v.

Both the midpoints and the radii of A are required to be upper triangular, since a nonzero radius below the diagonal would describe a matrix that is not; an ArgumentError is raised otherwise, and a DimensionMismatch when A is not square. A diagonal ball containing zero raises an ArgumentError from the division rather than returning an unbounded enclosure, so a matrix whose enclosure admits a singular member is rejected and not reported as solved.

forward_substitution(A::BallMatrix, v::BallVector)
forward_substitution(A::BallMatrix, B::BallMatrix)

The same as backward_substitution for a lower triangular A, eliminating the rows from the first to the last; midpoints and radii are both required to be lower triangular.

verified_cholesky([T = BigFloat,] A::AbstractMatrix; precision_bits = 256)

An inclusion of the Cholesky factor of A, by Section 4 of Rump and Ogita (2024); the method is described at _rumpogita2024_cholesky. Returns a VerifiedCholeskyResult: when certified is true, A is proved positive definite and A = G*G with G upper triangular with positive diagonal, in the ball matrix returned. Positive definiteness is a conclusion and not an assumption.

When A is not Hermitian a warning is issued (above a rounding-level threshold) and the statement is for its Hermitian part (A + A*)/2, enclosed in ball arithmetic.

Arguments

• T: the arithmetic of the enclosure, BigFloat (the default) or Float64, as a leading argument. In BigFloat the computation runs at precision_bits and the factor is enclosed to about that precision. Float64 with a BigFloat input raises.
• precision_bits: the BigFloat precision. A BigFloat input with more digits is rounded to it and the rounding is enclosed.

Example

B = randn(100, 100); A = B' * B + I
r = verified_cholesky(Float64, A)
r.certified      # A is positive definite and its Cholesky factor lies in r.G

Reference

S. M. Rump and T. Ogita, Verified Error Bounds for Matrix Decompositions, SIAM J. Matrix Anal. Appl. 45(4) (2024), 2155-2183, doi 10.1137/24M165096X, Section 4.

VerifiedCholeskyResult{GM, RT}

Result of verified_cholesky.

Fields

• G::GM: ball matrix containing the upper triangular Cholesky factor, A = G*G
• certified::Bool: whether the verification succeeded; when false the radii are infinite and nothing is proved, in particular not that A fails to be positive definite
• residual_norm::RT: an upper bound of ‖G̃*G̃ − A‖_∞ / ‖A‖_∞ over all G̃ ∈ G

Statement

When certified is true, the Hermitian part (A + A*)/2 of the input, which is the input when it is Hermitian, is proved positive definite, and its Cholesky factor, upper triangular with positive diagonal, lies in the ball matrix G.

Reference

S. M. Rump and T. Ogita, Verified Error Bounds for Matrix Decompositions, SIAM J. Matrix Anal. Appl. 45(4) (2024), 2155-2183, doi 10.1137/24M165096X, Section 4.

The flag is called certified since version 0.4 and was called success in 0.3.0; the old name still resolves, with a deprecation warning, and iscertified reads it in either version.

GaussianEliminationResult{VT, MT}

The result of interval_gaussian_elimination.

Fields

• solution::VT: a ball vector containing the solution of  x = b̂ for every Â in A and b̂ in b; its radii are infinite when not certified
• certified::Bool: whether the elimination ran to the end with no pivot containing zero. Then every matrix of A is nonsingular. false says only that the elimination stopped: it is not a proof that A contains a singular matrix
• L::MT, U::MT: ball matrices, unit lower triangular and upper triangular. When certified, for every Â in A the rows of Â taken in the order p factor as Â[p, :] = L̂ Û with L̂ in L and Û in U
• p::Vector{Int}: the order of the rows

The flag is called certified since version 0.4 and was called success in 0.3.0; the old name still resolves, with a deprecation warning, and iscertified reads it in either version.

interval_gaussian_elimination(A::BallMatrix, b::BallVector; pivoting = :partial)
    -> GaussianEliminationResult

Gaussian elimination and backward substitution in ball arithmetic for the square real interval system A x = b: Algorithms 5.9 and 5.10 of Horáček (2019).

Elimination, for k = 1, …, n − 1: a row j ≥ k whose entry in column k does not contain zero is brought to position k; for every row j > k the entry a_jk is set to zero and

a_j,k+1:n := a_j,k+1:n − (a_jk / a_kk) a_k,k+1:n,        b_j := b_j − (a_jk / a_kk) b_k.

If no such row exists the elimination stops, and the result is not certified. Backward substitution: x_i := (b_i − Σ_{j>i} a_ij x_j) / a_ii for i = n, …, 1.

Why the result is an enclosure: for a fixed Â in A and b̂ in b, Gaussian elimination on them with the same order of the rows evaluates the same formulae on numbers contained in the balls, and every operation on balls contains the result of the operation on their points; a pivot that does not contain zero is not zero for Â. So Â is nonsingular, its solution is in solution, and its factors are in L and U.

pivoting = :partial chooses, among the rows whose entry does not contain zero, the one whose entry is farthest from zero (the largest lower bound of its absolute value). Step 2 of Algorithm 5.9 asks only for a row whose entry does not contain zero and leaves the choice free; the thesis then says that it selects "a pivot with the smallest mignitude". pivoting = :none keeps the rows in their order and stops at the first pivot that contains zero.

The enclosure may be much wider than the hull of the solution set: the radii grow with the number of operations, as Example 5.11 of the thesis shows. The thesis recommends preconditioning by the inverse of the midpoint (compute_preconditioner).

Reference

J. Horáček, "Interval linear and nonlinear systems", PhD thesis, Charles University, Prague, 2019, Section 5.6.1.

interval_gaussian_elimination_det(A::BallMatrix) -> Ball

A ball containing the determinant of every matrix of the real interval matrix A, by Section 8.4.1 of Horáček (2019): the product of the diagonal of the row echelon form computed by interval_gaussian_elimination, with the sign of the permutation of the rows. The ball has infinite radius when the elimination stops at a pivot that contains zero; the elimination is then inconclusive, and the determinant may or may not vanish on A.

is_regular_gaussian_elimination(A::BallMatrix) -> Bool

true when interval_gaussian_elimination runs to the end on the real interval matrix A, which proves that every matrix of A is nonsingular. false is inconclusive.

IterativeResult{T, VT}

Result of interval_jacobi and interval_gauss_seidel.

Fields

• solution: a box containing the part of the solution set that lies in the initial box
• converged: whether the stopping criterion of Section 5.7.2 was met within max_iterations
• iterations: number of sweeps
• final_width: the largest radius of solution
• convergence_rate: the ratio of the largest radii of the last two iterates
• empty: true when an intersection was empty, which proves that the solution set does not meet the initial box; solution is then the last nonempty iterate and is not an enclosure

interval_initial_enclosure(A::BallMatrix, b::BallVector; C = inv(mid(A)), method = :auto)
    -> BallVector

A box that contains the solution set of the square interval system A x = b, by Section 5.7.1 of Horáček (2019), with C a real matrix (any; the default is the floating-point inverse of the midpoint):

• method = :proposition, Proposition 5.14: if ⟨CA⟩u ≥ v > 0 for some u ≥ 0, then Σ ⊆ ‖Cb‖_v [−u, u], with ‖x‖_v = max_i |x_i|/v_i. u is an approximate solution of ⟨CA⟩u = e and v a lower bound of ⟨CA⟩u, as the thesis suggests.
• method = :norm, formula (5.4): if ‖I − CA‖_∞ < 1, then Σ ⊆ ‖Cb‖_∞/(1 − ‖I − CA‖_∞) [−e, e].
• method = :auto: the first, and the second where the first does not apply.

CA and Cb are formed in ball arithmetic and every bound is rounded outwards. Where the condition of the method is not established the box returned has infinite radii.

interval_gauss_seidel(A::BallMatrix, b::BallVector; x0 = nothing, preconditioner = nothing,
    max_iterations = 20, tol = nothing) -> IterativeResult

Algorithm 5.17 of Horáček (2019), the Gauss-Seidel method: the Jacobi method (interval_jacobi, with the same arguments and the same statement) in which each component uses the components already updated in the sweep,

y_i⁽ᵏ⁺¹⁾ = (b_i − Σ_{j<i} a_ij x_j⁽ᵏ⁺¹⁾ − Σ_{j>i} a_ij x_j⁽ᵏ⁾)/a_ii.

interval_jacobi(A::BallMatrix, b::BallVector; x0 = nothing, preconditioner = nothing,
    max_iterations = 20, tol = nothing) -> IterativeResult

Algorithm 5.15 of Horáček (2019), the Jacobi method for a real square interval system A x = b: from a box x⁽⁰⁾ containing the solution set,

y_i⁽ᵏ⁺¹⁾ = (b_i − Σ_{j≠i} a_ij x_j⁽ᵏ⁾)/a_ii,        x⁽ᵏ⁺¹⁾ = x⁽ᵏ⁾ ∩ y⁽ᵏ⁺¹⁾,

until the stopping criterion of Section 5.7.2: both ends of every component move by less than tol, or max_iterations sweeps (the thesis uses 20). tol = nothing takes (5.5), min_ij wid(A_ij) × 10⁻⁵.

• x0: the initial box. When it contains the solution set, so does the result; when it contains only a part of it, the result contains that part. nothing takes interval_initial_enclosure, and where that gives none the result has infinite radii and converged = false.
• preconditioner: a real matrix C; the iteration is then applied to (CA) x = Cb, formed in ball arithmetic, whose solution set contains that of the original system (Section 5.4).

A diagonal entry that contains zero leaves its component unchanged. An empty intersection is reported in the field empty.

HBRResult{VT}

Result of hbr_method.

Fields

• solution: a box containing the solution set when certified is true; infinite radii otherwise
• certified: whether the hypotheses of Theorem 5.12 were verified and the enclosure computed

The flag is called certified since version 0.4 and was called success in 0.3.0; the old name still resolves, with a deprecation warning, and iscertified reads it in either version.

hbr_method(A::BallMatrix, b::BallVector; preconditioner = nothing) -> HBRResult

Theorem 5.12 of Horáček (2019). Let A x = b be a real square interval system with A an H-matrix, ⟨A⟩ its comparison matrix,

u = ⟨A⟩⁻¹ mag(b),    d_i = (⟨A⟩⁻¹)_ii,    α_i = ⟨A⟩_ii − 1/d_i,    β_i = u_i/d_i − mag(b_i).

Then the solution set is contained in the box with components

x_i = (b_i + [−β_i, β_i]) / (A_ii + [−α_i, α_i]).

When the midpoint of A is diagonal the box is the hull (the thesis, after the theorem).

The theorem is applied to the preconditioned system (CA) x = Cb, formed in ball arithmetic, whose solution set contains that of A x = b (Section 5.4). preconditioner = nothing takes C the floating-point inverse of the midpoint; a real matrix is used as given; hbr_method_simple applies the theorem to A x = b itself.

How it is verified: that the matrix is an H-matrix by is_H_matrix; ⟨A⟩⁻¹ is enclosed by verifylss, as the thesis suggests with the ε-inflation method; α_i and β_i are replaced by upper bounds, which enlarges the box, since the quotient grows with both. certified is false where one of these steps fails.

hbr_method_simple(A::BallMatrix, b::BallVector) -> HBRResult

Theorem 5.12 of Horáček (2019) applied to A x = b without preconditioning; A must then be an H-matrix itself. See hbr_method.

ShavingResult{T, VT}

Result of interval_shaving.

Fields

• solution: the shaved box; it contains every point of the solution set that the initial box contains
• shaved_amount: the sum of the widths removed
• iterations: number of passes over the variables
• components_shaved: number of slices removed
• empty: true when a slice as wide as the box was removed, which proves that the initial box contains no point of the solution set; solution is then meaningless

interval_shaving(A::BallMatrix, b::BallVector, x0::BallVector; passes = 1,
    improvements = 3, min_improvement = 1e-6, lp_solver = nothing) -> ShavingResult

The shaving method of Section 5.9 of Horáček (2019) for a real square interval system A x = b: slices of the box x0 that contain no point of the solution set are cut off. With x0 an enclosure of the solution set, so is the result.

For the variable i and the upper end, the slice of width α is the box with x_i restricted to [x̄_i − α, x̄_i] (5.8), and similarly at the lower end. For each of the 2n slices:

1. A certificate for the midpoint system (Section 5.9.1). By Lemma 5.21, A x = b, x ∈ x has no solution exactly when

   Aᵀw + y − z = 0,    bᵀw + x̄ᵀy − x̲ᵀz = −1,    y, z ≥ 0                        (5.9)

   is solvable. The linear programme min b_cᵀw + x̄ᵀy − x̲ᵀz subject to A_cᵀw + y − z = 0, −e ≤ w ≤ e, y, z ≥ 0 is solved for the slice of width zero. If its value is not negative the slice is kept.

2. A square interval system. The entries of y and z that vanish at the optimum are fixed to zero; the remaining system (5.12) is completed to a square one, C v = d, C ∈ C, with an orthonormal basis of the null space of its midpoint, the right-hand sides chosen so that the normalised optimum solves it (5.13). Its solution set is enclosed by verifylss. If the components u of the enclosure that come from y and z are nonnegative, (5.9) is solvable for every A ∈ A, b ∈ b, and the slice contains no point of the solution set.

3. The width (Section 5.9.2). The width α enters one coefficient, C + α e_p e_qᵀ, and by Lemma 5.22 the solution becomes v − α v_q/(1 + α C⁻¹_qp) · C⁻¹_*p. With c an enclosure of C⁻¹_*p, the solution set of C u = e_p, the largest α allowed is given by α < −1/c̲_q when c̲_q < 0 (5.15) and, for every component k of u other than q, α ≤ u̲_k / (sup(u_q c_k) − u̲_k c̲_q) when the denominator is positive (5.16), all rounded downwards.

4. Iterative improvement (Section 5.9.3), improvements times: the coefficient is moved by the width found, the system enclosed again, c updated by c/(1 + α c_q), with its q-th component computed as 1/(α + 1/c_q), and a further width added.

passes is the number of sweeps over the variables; a slice narrower than min_improvement times the width of its variable is not cut.

lp_solver(c, Aeq, beq, lb, ub) must return a minimiser of c'v subject to Aeq v = beq, lb ≤ v ≤ ub, or nothing. It is a heuristic: whatever it returns, a slice is cut only after step 2 has verified the certificate. nothing uses JuMP with HiGHS, which must then be loaded (using JuMP, HiGHS); see _shaving_lp.

Reference

J. Horáček, Interval linear and nonlinear systems, doctoral thesis, Charles University, Department of Applied Mathematics, Prague, 2019 ([26]), Section 5.9, Lemma 5.21 and Lemma 5.22.

sherman_morrison_inverse_update(A_inv::Matrix{T}, u::Vector{T}, v::Vector{T}) where {T}

Compute (A + uv^T)^(-1) using Sherman-Morrison formula.

Sherman-Morrison Formula

(A + uv^T)^(-1) = A^(-1) - (A^(-1)uv^T A^(-1)) / (1 + v^T A^(-1)u)

This allows O(n²) update of inverse for rank-1 perturbation instead of O(n³) recomputation.

Arguments

• A_inv: Inverse of base matrix A (n×n)
• u: First vector of rank-1 update (n)
• v: Second vector of rank-1 update (n)

Returns

(A + uv^T)^(-1)

Example

A = [3.0 1.0; 1.0 2.0]
A_inv = inv(A)
u = [1.0, 0.0]
v = [0.0, 1.0]

# Efficiently compute inv(A + u*v')
A_updated_inv = sherman_morrison_inverse_update(A_inv, u, v)

Notes

• O(n²) complexity vs O(n³) for full inverse
• Numerically stable if |1 + v^T A^(-1)u| is not too small
• A floating-point formula: the result is an approximation, with no error bound

_shaving_lp(c, Aeq, beq, lb, ub) -> Vector or nothing

A minimiser of c'v subject to Aeq v = beq, lb ≤ v ≤ ub, in floating point, or nothing.

The function has no method in the package itself: one is added by the extension loaded with JuMP and HiGHS, and interval_shaving raises an ArgumentError naming them when neither that extension nor an lp_solver of its own is available. Whatever the solver returns is a guess and is not verified: a slice is cut only after the square interval system of step 2 of interval_shaving has established the certificate for the whole ball.

PreconditionerType

Enumeration of available preconditioning strategies.

• :midpoint: C = A_c^(-1) (inverse of midpoint)
• :lu: C from LU factorization of midpoint
• :ldlt: C from LDLT factorization of midpoint (symmetric matrices)
• :identity: C = I (no preconditioning)

PreconditionerType: C = inv(mid(A)) in floating point.

PreconditionerType: C the inverse of the LU factorisation of mid(A), in floating point.

PreconditionerType: C the inverse of the LDLᵀ factorisation of the symmetric mid(A) with the pivoting of Bunch and Kaufman (bunchkaufman), in floating point.

PreconditionerType: C = I.

PreconditionerResult{CT, T, F}

Result from preconditioner computation. Nothing in it is verified: the preconditioner is an approximate inverse and the condition number an estimate.

Fields

• preconditioner::Matrix{CT}: Approximate inverse C, real or complex as the midpoint is
• method::PreconditionerType: Method used
• condition_number::T: Estimated condition number, a real number
• success::Bool: Whether computation succeeded
• factorization: the stored factorization (for LU/LDLT), nothing otherwise

compute_preconditioner(A::BallMatrix{T};
                      method::Symbol=:midpoint,
                      check_conditioning::Bool=true) where {T}

Compute preconditioner C ≈ A^(-1) for interval matrix A.

Methods

• :midpoint: Direct inverse of A_c (default, simplest)
• :lu: LU factorization of A_c (more stable)
• :ldlt: LDLT factorization of A_c (for symmetric A)
• :identity: Identity matrix (no preconditioning)

Arguments

• A: Interval matrix to precondition
• method: Preconditioning strategy
• check_conditioning: Compute and warn about condition number

Returns

PreconditionerResult with approximate inverse and diagnostic info.

Example

A = BallMatrix([4.0 1.0; 1.0 3.0], fill(0.1, 2, 2))

# Midpoint inverse (default)
prec = compute_preconditioner(A)
C = prec.preconditioner

# LU factorization (more stable)
prec_lu = compute_preconditioner(A, method=:lu)

Notes

• Midpoint inverse: O(n³), simplest but may be numerically unstable
• LU factorization: O(n³) once, O(n²) per solve, more stable
• LDLT factorization: O(n³/6) once, O(n²) per solve, only for symmetric
• Identity: O(1), use for well-conditioned or diagonally dominant matrices

apply_preconditioner(prec::PreconditionerResult{T}, v::Vector{T}) where {T}

Apply preconditioner to vector: compute C * v.

If factorization is available, uses it for efficiency.

Arguments

• prec: Preconditioner result from compute_preconditioner
• v: Vector to multiply

Returns

C * v

Example

A = BallMatrix([3.0 1.0; 1.0 2.0], fill(0.1, 2, 2))
prec = compute_preconditioner(A, method=:lu)

b = [5.0, 4.0]
Cb = apply_preconditioner(prec, b)

apply_preconditioner(prec::PreconditionerResult{T}, M::Matrix{T}) where {T}

Apply preconditioner to matrix: compute C * M.

Arguments

• prec: Preconditioner result
• M: Matrix to multiply

Returns

C * M

apply_preconditioner(prec::PreconditionerResult, A::BallMatrix) -> BallMatrix
apply_preconditioner(prec::PreconditionerResult, b::BallVector) -> BallVector

The preconditioned data C A and C b of Section 5.4 of Horáček (2019), formed in ball arithmetic with C = prec.preconditioner, so that the solution set of (CA) x = Cb contains that of A x = b whatever the accuracy of C.

is_well_preconditioned(A::BallMatrix, prec::PreconditionerResult; threshold = 0.5) -> Bool

true when ‖I − C A‖_∞ < threshold is proved for every matrix A of the ball, with C = prec.preconditioner: the product is formed in ball arithmetic and its norm bounded from above. With threshold ≤ 1 this makes CA strictly diagonally dominant (Example 4.25 of Horáček (2019)) and gives the initial enclosure (5.4). false means that the bound was not established.

OverdeterminedResult{VT}

Result of the routines for overdetermined interval systems.

Fields

• solution: a box containing the solution set (for interval_least_squares, the set of least squares solutions) when certified is true; infinite radii otherwise
• certified: whether a finite enclosure was obtained
• unsolvable: true when the solution set is proved empty, by an empty intersection of enclosures; solution is then meaningless
• method: the method used
• subsystems_checked: number of square subsystems for which an enclosure was obtained

select_subsquares(m, n; overlap = round(Int, n/3), rng = nothing) -> Vector{Vector{Int}}

Algorithm 6.14 of Horáček (2019): row index sets of square subsystems covering the m equations of an m × n system. The first is chosen at random; each following one takes overlap equations already covered and n − overlap not yet covered; the last takes the remaining ones and completes with covered ones. Their number is (6.5), 1 + ⌈(m − n)/(n − overlap)⌉. The thesis reports overlap ≈ n/3 as a reasonable choice.

subsquares_method(A::BallMatrix, b::BallVector; selection = :covering, k = 10,
    overlap = round(Int, n/3), solver = :hbr, rng = nothing) -> OverdeterminedResult

Algorithm 6.12 of Horáček (2019), the simple subsquares algorithm, for a real overdetermined interval system. The solution set lies in the solution set of every square subsystem (Proposition 6.5), so:

1. select square subsystems: selection = :random, k random ones, as in the algorithm; :covering, those of select_subsquares (Algorithm 6.14); :all, all binomial(m, n);
2. enclose the solution set of each by a method for square systems, solver = :hbr (hbr_method), :verifylss or :gaussian_elimination;
3. return the intersection of the enclosures;
4. if the intersection is empty, the system has no solution: unsolvable = true.

A subsquare for which the solver gives no enclosure is skipped. certified is false when none gives one.

multi_jacobi_method(A::BallMatrix, b::BallVector; x0 = nothing, overlap = round(Int, n/3),
    max_iterations = 20, tol = nothing, rng = nothing) -> OverdeterminedResult

Section 6.8.3 of Horáček (2019), the multi-Jacobi method: subsquares are selected by Algorithm 6.14 (select_subsquares); then, instead of solving each, one iteration of the Jacobi method (Algorithm 5.15) is applied to every subsquare from the current global box, the global box is intersected with all the results, and this is repeated until the stopping criterion of Section 5.7.2 or max_iterations.

• x0: a box containing the solution set. nothing takes the HBR enclosure of a subsquare, as in Example 6.15 of the thesis (the first selected subsquare for which it exists).
• Each subsquare is preconditioned by the floating-point inverse of its midpoint, the products formed in ball arithmetic. The thesis does not say how the diagonal of a subsquare is kept free of zero; this is the preconditioning of Section 5.4 applied to each subsquare.

An empty intersection proves that the system has no solution (unsolvable = true). certified is false when no initial box is found.

interval_least_squares(A::BallMatrix, b::BallVector; method = :supersquare,
    solver = :verifylss) -> OverdeterminedResult

Section 6.2 of Horáček (2019): an enclosure of the set of least squares solutions (Definition 6.3)

Σ_lsq = {x : AᵀA x = Aᵀb for some A ∈ A, b ∈ b},

whose hull contains that of the solution set of A x = b; the result is nonempty also when the system has no solution.

• method = :supersquare: the square system of order m + n

  [I A; Aᵀ 0] [y; x] = [b; 0],

  enclosed by solver, a method for square systems (:verifylss, :hbr or :gaussian_elimination); the part x is returned. The two occurrences of each entry of A are treated as independent, which the thesis notes as a source of overestimation.

• method = :normal_equations: the interval normal equation AᵀA x = Aᵀb, with the products formed in ball arithmetic. The thesis warns that the overestimation of the product may be huge.

certified is false when the solver gives no enclosure.

VerifiedLinearSystemResult

The result of verified_linear_solve_hmatrix.

• solution: the ball vector with midpoint x_approx and radius error_bound; when certified it contains the solution of every system with data in the given balls;
• x_approx: the approximate solution x̃;
• error_bound: ε with |Â⁻¹ b̂ − x̃| ≤ ε entrywise, infinite when not certified;
• certified: whether the hypothesis of the theorems was proved;
• method: the bound that was computed;
• R: the approximate inverse used;
• v: the positive vector with ⟨C⟩ v > 0, empty when none was found.

The flag is called certified since version 0.4 and was called verified in 0.3.0; the old name still resolves, with a deprecation warning, and iscertified reads it in either version.

verified_linear_solve_hmatrix(A::BallMatrix, b::BallVector;
    method = :theorem_3_3, m = 1, R = nothing, x_approx = nothing)

An entrywise bound ε of the error of an approximate solution x̃ of a real linear system: |Â⁻¹ b̂ − x̃| ≤ ε for every matrix Â in A and every vector b̂ in b, together with the proof that every such Â is nonsingular. Returns a VerifiedLinearSystemResult; certified is false, and the bound infinite, when the hypothesis could not be proved.

The hypothesis, common to all four bounds: with R an approximate inverse and C a ball matrix containing R Â for every Â in A, there is a vector v > 0 with u := ⟨C⟩ v > 0, where ⟨C⟩ is the comparison matrix (see comparison_matrix). Then C is an H-matrix, Â and R are nonsingular, and |(R Â)⁻¹| ≤ ⟨C⟩⁻¹ (Theorems 2.1 and 2.3, Remark 1). Write ⟨C⟩ = D − E with D the diagonal part, G := E D⁻¹, and c for an upper bound of |R (b̂ − Â x̃)|.

method selects the bound:

• :theorem_3_1 (Rump): ε = (D⁻¹ + v wᵀ) c with w_k := max_i G_ik / u_i, equation (3);
• :theorem_3_2 (Minamihata, Sekine, Ogita, Rump and Oishi 2015): ε = (D⁻¹ + v wᵀ) (I − Σ_{k=1}^m (Δ^{2k−1} − Δ^{2k})) c with Δ := u wᵀ − G, equation (4);
• :theorem_3_3 (the default): ε = min_{1 ≤ k ≤ m} ε^(k), ε^(k) := D⁻¹ (c^(0) + … + c^(k−1)) + max_i (c^(k)_i / u_i) v, with c^(0) := c and c^(k) := G c^(k−1), equations (5) and (6);
• :corollary_3_2: ε = D⁻¹ r + α v with r := ⟨mid C⟩ |ỹ| + |mid(C) ỹ − R (b̂ − Â x̃)|, α := max_i (G r)_i / u_i and ỹ an approximate solution of mid(C) y = R (b − A x̃), equation (17).

m is the number of terms in Theorems 3.2 and 3.3. R defaults to the floating-point inverse of mid(A) and x_approx to R mid(b).

How it is computed. C and the residual R (b − A x̃) are products and sums of balls, so they contain what they should for every point of A and b; the residual is formed in working precision. v is the floating-point solution of ⟨C⟩ v = (1, …, 1)ᵀ, the first of the two choices named in Section 3 of the paper, and u is then enclosed, its lower end being the one that has to be positive. The bounds of Theorems 3.1, 3.3 and of Corollary 3.2 are sums and products of nonnegative numbers, rounded upwards, with u replaced by its lower end where it divides. The bound of Theorem 3.2 has terms of both signs and is evaluated in balls, w being any vector with G ≤ u wᵀ, which is all its proof uses; the upper end is returned.

What is not reproduced. The paper distinguishes two enclosures of R A, Method (a) from two products with directed rounding and Method (b) from one product and the a priori bound 2 n u |R| |A|, and proposes Corollary 3.2 for the second. Here there is one enclosure, the ball product of the package, and each of the four bounds is applied to it (Remark 5 of the paper allows Corollary 3.2 on any enclosure); the comparison of costs between the two methods does not carry over.

Reference

A. Minamihata, T. Ogita, S. M. Rump and S. Oishi, "Modified error bounds for approximate solutions of dense linear systems", Journal of Computational and Applied Mathematics 369 (2020), 112546, doi 10.1016/j.cam.2019.112546.

comparison_matrix(C::BallMatrix)

The comparison matrix ⟨C⟩ of Definition 2.7 of Minamihata, Ogita, Rump and Oishi (2020), in floating point: mig of the entry on the diagonal, minus mag of the entry off it. Both are rounded so that the matrix returned is entrywise below the comparison matrix of C; a positive vector that it maps to a positive vector therefore proves that C is an H-matrix (Remark 1 there).

is_H_matrix(A::BallMatrix) -> Bool

true when every matrix of the ball A, real or complex, is proved to be a nonsingular H-matrix: its comparison matrix M(A), with |a_ii| on the diagonal and −|a_ij| off it, is a nonsingular M-matrix (Varga (2004), Definition C.6), which makes A nonsingular (Proposition C.7 there). The test is that of is_M_matrix on the comparison matrices of the ball; false means that this was not proved.

Reference

R. S. Varga, Geršgorin and His Circles, Springer Series in Computational Mathematics 36, Springer, Berlin, 2004, doi 10.1007/978-3-642-17798-9 ([27]), Appendix C, (C.4), Definition C.6 and Proposition C.7.

krawczyk_linear_system(A::BallMatrix, b::BallVector; max_iterations = 20) -> KrawczykResult

A verified enclosure of the solution of A x = b, for every matrix of the ball A and every vector of the ball b: verifylss, which is Algorithm 10.7 of Rump (2010), with its result in the fields of a KrawczykResult.

With R a floating-point inverse of the midpoint of A and x̃ = R mid(b), the algorithm iterates X ← Z + C·Y, Z ∋ R(b − Ax̃), C ∋ I − RA, Y an inflation of X, all in ball arithmetic, and stops when X is in the interior of Y; then A and R are nonsingular and the solution is in x̃ + X (Theorem 10.8 there). This is the iteration of Krawczyk (1969) for a linear system with the inflation of Rump.

certified is the success of that test. iterations is the number of iterations, contraction_factor an upper bound of ‖I − RA‖₂ (a diagnostic: the test may succeed when it is not below one), residual_norm an upper bound of ‖A x̂ − b‖₂ over the two balls for x̂ the midpoint of the returned enclosure. When certified is false the enclosure proves nothing.

The preconditioner and the approximate solution are those of the algorithm; the keywords R, x_approx and expansion_factor of earlier versions are gone with the code that used them, which formed I − RA in floating point and left the residual out of the inclusion test.

References

S. M. Rump, Verification methods: rigorous results using floating-point arithmetic, Acta Numerica 19 (2010), 287-449, doi 10.1017/S096249291000005X ([28]), Algorithm 10.7 and Theorem 10.8.

R. Krawczyk, Newton-Algorithmen zur Bestimmung von Nullstellen mit Fehlerschranken, Computing 4(3) (1969), 187-201, doi 10.1007/BF02234767 ([29]).

KrawczykResult{T, VT}

Result from Krawczyk verification method.

The flag is called certified since version 0.4 and was called verified in 0.3.0; the old name still resolves, with a deprecation warning, and iscertified reads it in either version.

verifylss(A::BallMatrix, b; method = :rump2010, r = 0.1, ϵ = 1e-20, iter_max = 20)
    -> VerifyLssResult

Verified enclosure of the solution of A x = b, with b a BallVector or a BallMatrix. When the returned certified is true, every matrix in the ball A is nonsingular and the enclosure contains the solution for every A ∈ A and b ∈ b; nonsingularity is a conclusion of the method, not a hypothesis.

This function selects an algorithm and does nothing else. Each algorithm is a separate unexported function named for the paper it implements:

Reference

S. M. Rump, Verification methods: rigorous results using floating-point arithmetic, Acta Numerica 19 (2010) 287-449, doi 10.1017/S096249291000005X. Algorithm 10.7 and Theorem 10.8; Rump attributes the algorithm to Rump (1980).

VerifyLssResult{T, VT}

Result from epsilon-inflation linear system verification.

Fields

• solution::VT: Enclosure of the solution (BallVector or BallMatrix)
• certified::Bool: Whether the solution is mathematically certified
• iterations::Int: Number of iterations performed
• spectral_radius_bound::T: Upper bound on ‖I - RA‖₂ (convergence requires < 1)
• condition_number::T: Approximate condition number of A

Convergence Diagnostics

If certified == false, check:

• spectral_radius_bound ≥ 1: Iteration matrix doesn't contract, method cannot converge
• condition_number > 1/eps(T): Matrix is numerically singular
• iterations == iter_max: May need more iterations or smaller inflation parameters

sylvester_miyajima_enclosure(A, B, C, X̃)

An entrywise enclosure of the solution X* of A X + X B = C, as a BallMatrix with midpoint X̃ (an approximate solution, in floating point) and radius X^ε with |X̃ − X*| ≤ X^ε. A (m×m), B (n×n) and C may be matrices or BallMatrix; with balls, the enclosure holds for the solution of every equation with data in them.

Theorems 1 and 2 of

S. Miyajima, Fast enclosure for solutions of Sylvester equations, Linear Algebra Appl. 439 (2013), 856–878, doi:10.1016/j.laa.2012.07.001 ([3]),

Section 2.1 (spectral decomposition). The block-diagonalisation variant of Section 2.2 (Theorems 3 and 4), the accelerations of Section 3 and the refinement of Section 4 are not implemented.

The statement

Let Ṽ_A, W_A (m×m), Ṽ_B, W_B (n×n) be any matrices and D̃_A, D̃_B diagonal; in practice A ≈ Ṽ_A D̃_A Ṽ_A⁻¹, Bᵀ ≈ Ṽ_B D̃_B Ṽ_B⁻¹, W_A ≈ Ṽ_A⁻¹, W_B ≈ Ṽ_B⁻¹, all computed in floating point. E is the m×n matrix of ones, ./ and |·| are entrywise, ‖·‖_∞ is the maximum row sum and ‖·‖_M the largest entry in modulus.

R_A = W_A(Ṽ_A D̃_A − A Ṽ_A),     R_B = W_B(Ṽ_B D̃_B − Bᵀ Ṽ_B),
S_A = I − W_A Ṽ_A,              S_B = I − W_B Ṽ_B,
T_A = |R_A| + ‖R_A‖_∞/(1 − ‖S_A‖_∞) |S_A|,     T_B likewise,
T = T_A E + E T_Bᵀ,     D̃ = D̃_A E + E D̃_B,     T_D = T ./ |D̃|.

Theorem 1. If ‖S_A‖_∞ < 1, ‖S_B‖_∞ < 1, D̃ has no zero entry and ‖T_D‖_M < 1, the equation has a unique solution X*.

With R = A X̃ + X̃ B − C, R_W = W_A R W_Bᵀ, and ρ_i(M), κ_j(M) the largest entry in modulus of row i and of column j,

R_W⁽¹⁾ = |R_W| + 1/(1 − ‖S_B‖_∞) · diag(ρ(R_W)) E |S_B|ᵀ,
R_V⁽¹⁾ = R_W⁽¹⁾ + 1/(1 − ‖S_A‖_∞) · |S_A| E diag(κ(R_W⁽¹⁾)),
R_W⁽²⁾ = |R_W| + 1/(1 − ‖S_A‖_∞) · |S_A| E diag(κ(R_W)),
R_V⁽²⁾ = R_W⁽²⁾ + 1/(1 − ‖S_B‖_∞) · diag(ρ(R_W⁽²⁾)) E |S_B|ᵀ,
R_V = min(R_V⁽¹⁾, R_V⁽²⁾),     R_D = R_V ./ |D̃|,
U = R_D + ‖R_D‖_M/(1 − ‖T_D‖_M) · T_D.

Theorem 2. Under the hypotheses of Theorem 1, |X̃ − X*| ≤ X^ε := |Ṽ_A| U |Ṽ_B|ᵀ.

R_V bounds |Ṽ_A⁻¹ R Ṽ_B⁻ᵀ|, the residual in the exact inverse frames, from R_W, the one in the computed W_A, W_B.

How it is evaluated

The products of matrices are ball products, of which the entrywise upper bound is taken. The products with E are row sums ((M E)_ij is the sum of row i of M), accumulated rounded up; |D̃_ij| is bounded from below; every other operation is rounded up. Throws ArgumentError when a hypothesis of Theorem 1 is not proved.

verified_sylvester_enclosure(A, B, C; X̃ = nothing, prefer_complex_schur = true)

sylvester_miyajima_enclosure with the approximate solution computed by schur_sylvester_midpoint when X̃ is not given: an entrywise enclosure of the solution of A X + X B = C for every datum of the balls, by Theorems 1 and 2 of [3].

It throws ArgumentError when the hypotheses of that enclosure are not proved, as happens for a defective A or B; the message names the condition that failed, one of ‖S_A‖_∞ < 1, ‖S_B‖_∞ < 1, a zero spectral gap or ‖T_D‖_M < 1. The candidate X̃ enters the enclosure only through the residual, so its accuracy is not a hypothesis of the theorems and a X̃ of one's own may be passed.

triangular_sylvester_miyajima_enclosure(T, k; sylvester_fallback = :direct)

For the upper triangular T = [T₁₁ T₁₂; 0 T₂₂] with T₁₁ of size k×k, an enclosure of the solution Y of

T₂₂* Y − Y T₁₁* = T₁₂*,

that is of A Y + Y B = C with A = T₂₂*, B = −T₁₁*, C = T₁₂*. T is a matrix or a BallMatrix; with a ball the enclosure holds for every matrix in it, the radii entering the residuals (no perturbation expansion is used).

An approximate solution is computed column by column, and enclosed by sylvester_miyajima_enclosure. When that does not apply (its hypotheses are not proved, typically for large triangular matrices whose eigenvector matrix is ill conditioned), the enclosure comes from the triangular structure, which requires the radii of T below the diagonal to be zero:

• sylvester_fallback = :direct: the column recurrence solved by forward substitution in ball arithmetic;
• sylvester_fallback = :residual: the residual of the approximate solution in ball arithmetic and, column by column from the last, (A + b_jj I) e_j = R_j − Σ_{l>j} e_l b_lj for the error e_j of column j, so ‖e_j‖₂ ≤ ‖(A + b_jj I)⁻¹‖₂ (‖R_j‖₂ + Σ_{l>j} ‖e_l‖₂ |b_lj|), the norm of the inverse of the triangular matrix bounded as in triangular_inverse_two_norm_bound. The radius is uniform in each column.

SylvesterResolventResult{T}

The quantities of the bound of parametric_resolvent_bound that do not depend on z, for the split T = [T₁₁ T₁₂; 0 T₂₂] at index k and an approximate solution X of T₁₁X − XT₂₂ = −T₁₂.

Fields

• residual_norm: upper bound of ‖R‖₂, R = T₁₂ + T₁₁X − XT₂₂
• coupling_norm: upper bound of ‖T₁₂‖₂
• similarity_cond: upper bound of κ₂(S(X))
• reduction_factor, net_penalty: the diagnostics residual_norm/coupling_norm and its product with similarity_cond, in plain floating point; they enter no bound
• k, n: the split index and the size of T
• certified: true when the bounds were obtained; false, with failure_reason, when the structure is not the required one or X is not finite
• X: the matrix X used
• R: a BallMatrix containing R

The flag is called certified since version 0.4 and was called precomputation_success in 0.3.0; the old name still resolves, with a deprecation warning, and iscertified reads it in either version.

parametric_resolvent_bound(precomp, T, z, config = config_v1(); svd_warm_start = nothing)

An upper bound of ‖(zI − T)⁻¹‖₂ for the block upper triangular matrix T with exact entries, through the similarity S(X) = [I X; 0 I] of precomp = sylvester_resolvent_precompute(T, k).

The bound

With S(X)⁻¹TS(X) = [T₁₁ R; 0 T₂₂], R = T₁₂ + T₁₁X − XT₂₂ (an identity for every X), A_z = zI − T₁₁ and D_z = zI − T₂₂,

(zI − T)⁻¹ = S(X) [A_z⁻¹  A_z⁻¹RD_z⁻¹; 0  D_z⁻¹] S(X)⁻¹,

‖(zI − T)⁻¹‖₂ ≤ κ₂(S(X)) ( ‖A_z⁻¹‖₂ + ‖D_z⁻¹‖₂ + ‖A_z⁻¹RD_z⁻¹‖₂ ).

It is the bound for a block triangular matrix after a similarity, so it carries the condition number κ₂(S(X)) = ψ(‖X‖₂)² (psi_squared) where the bound computed directly on T carries none; in exchange the certified singular value computation is on the k×k block only.

The terms

1. ‖A_z⁻¹‖₂ ≤ 1/σ, σ a certified lower bound of the smallest singular value of the ball matrix zI − T₁₁ (rigorous_svd with the method of config).
2. ‖D_z⁻¹‖₂ by the DInverseEstimator of config.
3. ‖A_z⁻¹RD_z⁻¹‖₂ by the CouplingEstimator of config.

The moduli |z − t_ii| are rounded down, R and the residuals of the approximate solves are ball matrices, and the scalar operations are rounded up.

Returns a ParametricResolventResult. T must be the matrix precomp was computed from.

parametric_resolvent_bound(precomp, T, z_list::AbstractVector, config = config_v1())

The bound at each point of z_list.

parametric_resolvent_bound(T, k::Int, z, config = config_v1(); X_oracle = nothing)

Precompute for the split k and evaluate at z (a point or a vector of points). Returns (precomp, result).

ParametricResolventResult{T}

The bound of parametric_resolvent_bound at a point z and its terms: resolvent_bound (upper bound of ‖(zI − T)⁻¹‖₂), K_S (of κ₂(S(X))), M_A (of ‖(zI − T₁₁)⁻¹‖₂), M_D (of ‖(zI − T₂₂)⁻¹‖₂), r (of ‖R‖₂) and coupling_term (of ‖(zI − T₁₁)⁻¹R(zI − T₂₂)⁻¹‖₂), with the config used. When certified is false the numbers are Inf and failure_reason says which step did not apply.

The flag is called certified since version 0.4 and was called success in 0.3.0; the old name still resolves, with a deprecation warning, and iscertified reads it in either version.

schur_sylvester_midpoint(A, B, C; prefer_complex_schur = true)

An approximate solution of A X + X B = C in floating point: Schur forms A = Q_A T_A Q_A* and B = Q_B T_B Q_B* reduce the equation to T_A Y + Y T_B = Q_A* C Q_B, which is solved by LinearAlgebra.sylvester, and the solution is transformed back as Q_A Y Q_B*. This is the construction of Bartels and Stewart.

prefer_complex_schur = true takes the Schur forms of the complexified data, which makes T_A and T_B triangular rather than quasi-triangular.

No bound comes with the result, and nothing here is verified: the reduction, the triangular solve and the two products are all in floating point. It is the default candidate of verified_sylvester_enclosure, which is where a bound is attached to it.

Reference

R. H. Bartels and G. W. Stewart, Algorithm 432: Solution of the matrix equation AX + XB = C, Comm. ACM 15(9) (1972), 820-826, doi 10.1145/361573.361582 ([30]).

_sylvester_triangular_columns(A, B, C)

An approximate solution of A X + X B = C in floating point, for A and B lower triangular, computed column by column in O(m²k) operations with m = size(A, 1) and k = size(B, 1).

Sweeping j = k, k−1, …, 1, the j-th column solves

(A + B[j,j] I) x_j = c_j − Σ_{l>j} x_l B[l,j],

whose coefficient matrix is lower triangular, by _forward_substitution. No bound comes with the result: it is the midpoint that the enclosures of this file take as their candidate, and _sylvester_triangular_direct_ball is the same recurrence carried in ball arithmetic.

_forward_substitution(L, b)

The solution of L x = b for a lower triangular L, by forward substitution in floating point.

Nothing is bounded here: the rounding of the substitution is not accounted for and L[i, i] is not tested, so a zero on the diagonal gives a non-finite result rather than an error. The verified counterpart, which returns an enclosure for a lower triangular BallMatrix, is forward_substitution.

_sylvester_triangular_direct_ball(A, B, C)

An entrywise enclosure of the solution X of A X + X B = C, for A and B lower triangular, holding for every datum of the balls A, B and C.

It is the column recurrence of _sylvester_triangular_columns carried in ball arithmetic throughout: the right-hand side c_j − Σ_{l>j} x_l B[l,j] is formed from the enclosures of the columns already computed, the diagonal of A + B[j,j] I is a sum of balls, so its rounding goes into the radius, and the column is enclosed by forward_substitution on the resulting lower triangular BallMatrix.

Nothing is assumed about the eigenvectors of A or B, which is why this route remains available where the hypotheses of sylvester_miyajima_enclosure are not proved. Where a ball on the diagonal of A + B[j,j] I contains zero no enclosure exists, and forward_substitution raises an ArgumentError naming that ball.

_has_distinct_diagonal(A::AbstractMatrix, tol)

true when |A[i,i] − A[j,j]| > tol for every i ≠ j.

The test is run in floating point on the entries as they stand, with no rounding accounted for, and tol is a threshold rather than a bound: true is not a proof that the diagonal of any nearby matrix is distinct. It is used only to decide whether the O(n²) substitution route of triangular_eigenvectors is taken instead of a dense eigen, a choice of method that does not enter any bound.

NormEstimator

The upper bound of the spectral norm used for the matrices of the approximate solves.

• OneInfNorm: √(‖·‖₁‖·‖_∞)
• FrobeniusNorm: ‖·‖_F

NormEstimator: the spectral norm bounded by √(‖·‖₁‖·‖_∞).

NormEstimator: the spectral norm bounded by the Frobenius norm.

DInverseEstimator

How ‖D_z⁻¹‖₂ is bounded for the upper triangular D_z = zI − T₂₂.

• TriBacksub: the recursion of triangular_inverse_two_norm_bound
• NeumannOneInf: D_z = Δ(I + N) with Δ the diagonal of D_z and N = Δ⁻¹(D_z − Δ), so that ‖D_z⁻¹‖_p ≤ ‖Δ⁻¹‖_p/(1 − ‖N‖_p) when ‖N‖_p < 1, for p = 1 and p = ∞, combined by ‖·‖₂ ≤ √(‖·‖₁‖·‖_∞)
• NeumannCollatz2: the same factorisation in the spectral norm, ‖N‖₂ bounded by collatz_upper_bound_L2_opnorm applied to an entrywise upper bound of |N|

DInverseEstimator: ‖D_z⁻¹‖₂ by the recursion of triangular_inverse_two_norm_bound.

DInverseEstimator: ‖D_z⁻¹‖₂ by the Neumann bound for D_z = Δ(I + N) in the norms 1 and ∞, combined by ‖·‖₂ ≤ √(‖·‖₁‖·‖_∞).

DInverseEstimator: ‖D_z⁻¹‖₂ by the Neumann bound for D_z = Δ(I + N) in the spectral norm, ‖N‖₂ bounded by collatz_upper_bound_L2_opnorm.

CouplingEstimator

How the off-diagonal term ‖A_z⁻¹ R D_z⁻¹‖₂ is bounded.

• CouplingNone: ‖A_z⁻¹‖₂‖R‖₂‖D_z⁻¹‖₂
• CouplingARSolve: ‖A_z⁻¹R‖₂‖D_z⁻¹‖₂, with A_z⁻¹R = Y + A_z⁻¹(R − A_zY) for an approximate solution Y of A_zY = R, so ‖A_z⁻¹R‖₂ ≤ ‖Y‖₂ + ‖A_z⁻¹‖₂‖R − A_zY‖₂
• CouplingOffDirect: A_z⁻¹RD_z⁻¹ = W + A_z⁻¹(R − A_zWD_z)D_z⁻¹ for an approximate W, so the term is at most ‖W‖₂ + ‖A_z⁻¹‖₂‖R − A_zWD_z‖₂‖D_z⁻¹‖₂

In the last two the residual is computed in ball arithmetic; when the approximate solve does not return finite numbers the first bound is used.

CouplingEstimator: ‖A_z⁻¹RD_z⁻¹‖₂ ≤ ‖A_z⁻¹‖₂‖R‖₂‖D_z⁻¹‖₂.

CouplingEstimator: ‖A_z⁻¹RD_z⁻¹‖₂ ≤ ‖A_z⁻¹R‖₂‖D_z⁻¹‖₂, with ‖A_z⁻¹R‖₂ bounded through an approximate solution of A_zY = R and its residual in ball arithmetic.

CouplingEstimator: ‖A_z⁻¹RD_z⁻¹‖₂ bounded through an approximate W of A_zWD_z = R and the residual R − A_zWD_z in ball arithmetic.

RegularityResult{T}

Result of a regularity test.

Fields

• certified::Bool: true when regularity is proved; false when the test is inconclusive
• method::Symbol: the condition that was tried last
• certificate::T: the quantity compared in that condition, as described by each routine

The conditions are sufficient: certified = false does not say that the matrix is singular.

The flag is called certified since version 0.4 and was called is_regular in 0.3.0; the old name still resolves, with a deprecation warning, and iscertified reads it in either version.

is_regular(A::BallMatrix; methods = [:diagonal_dominance, :gershgorin, :h_matrix,
    :sufficient_condition, :spectral_radius], verbose = false) -> RegularityResult

Try the sufficient conditions for regularity in the given order and return the first that proves it, or the last result when none does. :h_matrix is is_H_matrix (an H-matrix is regular, Section 4.4 of Horáček (2019)); the others are the routines of this file. Every one of them is sufficient only, so a certified = false result from all of them leaves the question open.

Reference

J. Horáček, Interval linear and nonlinear systems, doctoral thesis, Charles University, Department of Applied Mathematics, Prague, 2019 ([26]), Section 4.4.

is_regular_sufficient_condition(A::BallMatrix; condition = 1) -> RegularityResult

Theorem 11.12 of Horáček (2019): the interval matrix A is regular if

1. λ_max(A_Δᵀ A_Δ) < λ_min(A_cᵀ A_c) (condition = 1), or
2. A_cᵀ A_c − ‖A_Δᵀ A_Δ‖ I is positive definite for some consistent matrix norm (condition = 2, evaluated with the ∞-norm).

Both are evaluated through singular values: λ_min(A_cᵀ A_c) = σ_min(A_c)² ≥ 1/b² with b ≥ ‖A_c⁻¹‖₂ from svd_bound_L2_opnorm_inverse, and λ_max(A_Δᵀ A_Δ) ≤ ν² with ν ≥ ‖A_Δ‖₂ from upper_bound_L2_opnorm, or, for the second condition, an upper bound of ‖A_Δᵀ A_Δ‖_∞ computed in ball arithmetic. The certificate is a lower bound of the difference of the two sides, and regularity is proved when it is positive.

The thesis states the theorem for real interval matrices. For a complex ball matrix the first condition still implies regularity, since σ_min(A_c + E) ≥ σ_min(A_c) − ‖E‖₂ and ‖E‖₂ ≤ ‖A_Δ‖₂ for |E| ≤ A_Δ; the second is weaker than the first.

Reference

J. Horáček, Interval linear and nonlinear systems, doctoral thesis, Charles University, Department of Applied Mathematics, Prague, 2019 ([26]), Theorem 11.12.

is_regular_gershgorin(A::BallMatrix) -> RegularityResult

Regularity from the Geršgorin discs of the rows or of the columns: if zero lies outside every row disc of every matrix of the ball, or outside every column disc, the matrix is regular. For the rows this is strict diagonal dominance (is_regular_diagonal_dominance); for the columns it is the same condition on the transpose. The certificate is the larger of the two margins, and the condition is sufficient only.

Reference

J. Horáček, Interval linear and nonlinear systems, doctoral thesis, Charles University, Department of Applied Mathematics, Prague, 2019 ([26]), Definition 4.24 and Section 4.5.

is_regular_diagonal_dominance(A::BallMatrix) -> RegularityResult

Definition 4.24 of Horáček (2019): A is strictly diagonally dominant when

mig(a_ii) > Σ_{k ≠ i} mag(a_ik)        for every i,

with mig and mag the smallest and the largest modulus over an entry. Such a matrix is an H-matrix and therefore regular (Section 4.5 of the thesis). The certificate is a lower bound of the smallest difference of the two sides over the rows, and the condition is sufficient only.

Reference

J. Horáček, Interval linear and nonlinear systems, doctoral thesis, Charles University, Department of Applied Mathematics, Prague, 2019 ([26]), Definition 4.24 and Section 4.5.

is_regular_spectral_radius(A::BallMatrix; R = nothing) -> RegularityResult

Theorem 4.2 of Horáček (2019): a square interval matrix A is regular if for some real matrix R

ϱ(|I − R A_c| + |R| A_Δ) < 1.

R = nothing takes the floating-point inverse of the midpoint, which need not be exact: the theorem holds for any R. The matrix |I − R A_c| + |R| A_Δ is bounded above entrywise in ball arithmetic and its spectral radius by collatz_upper_bound; the certificate is that bound, and regularity is proved when it is below one. With R = A_c⁻¹ the condition is ϱ(|A_c⁻¹| A_Δ) < 1, strong regularity (Theorem 4.33 of the thesis). The condition is sufficient only: certified = false does not say that the ball contains a singular matrix.

Reference

J. Horáček, Interval linear and nonlinear systems, doctoral thesis, Charles University, Department of Applied Mathematics, Prague, 2019 ([26]), Theorem 4.2 and Theorem 4.33.

is_singular_sufficient_condition(A::BallMatrix; conditions = 1:4) -> Bool

Theorem 11.13 of Horáček (2019): a real interval matrix A is singular, that is, it contains a singular matrix, if at least one of the following holds:

1. λ_max(A_cᵀ A_c) ≤ λ_min(A_Δᵀ A_Δ);
2. max_j (|A_c⁻¹| A_Δ)_jj ≥ 1;
3. (A_Δ − |A_c|)⁻¹ ≥ 0;
4. A_Δᵀ A_Δ − A_cᵀ A_c is positive semidefinite.

Returns true when one of the conditions listed in conditions is proved, false otherwise, which does not mean that the matrix is regular. How each is evaluated:

1. as ‖A_c‖₂ ‖A_Δ⁻¹‖₂ ≤ 1 with upper bounds of the two norms;
2. with an enclosure of A_c⁻¹ by verifylss, the moduli of its entries bounded from below and the products with A_Δ rounded down;
3. with an enclosure of the inverse of the ball matrix A_Δ − |A_c| (the subtraction is rounded) by verifylss, all of whose lower ends must be nonnegative; an inverse with an entry equal to zero is therefore not recognised;
4. with enclosures of the eigenvalues of the symmetric ball matrix A_Δᵀ A_Δ − A_cᵀ A_c by verifyeigall on that HermitianBallMatrix (Rump and Lange 2023), all of whose lower ends must be nonnegative.

Reference

J. Horáček, Interval linear and nonlinear systems, doctoral thesis, Charles University, Department of Applied Mathematics, Prague, 2019 ([26]), Theorem 11.13.

DeterminantResult

Result of a determinant enclosure.

Fields

• determinant::Ball: a ball containing det(A) for every matrix A of the input; radius Inf where the method gives no bound
• method::Symbol: the method used

interval_det(A::BallMatrix; method = :auto) -> DeterminantResult

An enclosure of the determinant of the interval matrix A by

• :gaussian_elimination: Section 8.4.1 of Horáček (2019), the product of the diagonal of the row echelon form (interval_gaussian_elimination_det);
• :gershgorin: det_gershgorin;
• :hadamard: det_hadamard;
• :cramer: det_cramer;
• :cofactor: det_cofactor;
• :auto: the cofactor expansion for orders up to three; above, Gaussian elimination for a real matrix and :cramer for a complex one, the elimination being for real matrices.

det_hadamard(A::BallMatrix) -> DeterminantResult

Section 8.4.3 of Horáček (2019): by Hadamard's inequality, for every matrix of the interval matrix A,

|det(A)| ≤ d = ∏_i ( Σ_j mag(A_ji)² )^(1/2),

the product of the largest Euclidean norms of the columns. The same holds with the rows, and the smaller of the two bounds is returned, as the ball of centre zero and radius d. The sign of the determinant cannot be detected in this way. Real or complex input.

det_gershgorin(A::BallMatrix) -> DeterminantResult

Section 8.4.2 of Horáček (2019), for a real interval matrix. Every eigenvalue of every matrix of A lies in an interval Gerschgorin disc, of centre mid(A_ii) and radius rad(A_ii) + Σ_{j≠i} mag(A_ij), and a group of k discs whose union is connected (a "bunch") contains k eigenvalues. With [a, b] the projection of a bunch on the real axis, the product of its k eigenvalues lies in the interval of Figures 8.3 and 8.4 of the thesis:

• [a, b] for one disc;
• [a^k, b^k] when a > 0; when b < 0, [a^k, b^k] for odd k and [b^k, a^k] for even k;
• when a ≤ 0 ≤ b and |a| ≤ |b|, [a b^(k−1), b^k];
• when a ≤ 0 ≤ b and |a| > |b|, [a^(k−1) b, a^k] for even k and [a^k, a^(k−1) b] for odd k.

Complex eigenvalues come in conjugate pairs and the product of a pair, the square of its modulus, lies between the squares of the smallest and the largest modulus of [a, b] (Figure 8.5), so the same intervals hold. The determinant is enclosed by the product of the intervals of the bunches. The ends of the discs and the powers are rounded outwards.

The thesis remarks that the overestimation may be large, and that the method suits tight intervals or a matrix close to a diagonal one. Throws on complex input, where the discs are not centred on the real axis.

det_cramer(A::BallMatrix) -> DeterminantResult

Section 8.4.4 of Horáček (2019). By Cramer's rule the first coefficient of the solution of A x = e₁ is x₁ = det(A_{2:n})/det(A), with A_{2:n} the matrix without its first row and column, so that (8.3)

det(A) ⊆ det(A_{2:n}) / x₁

for an enclosure x₁ of that coefficient over the interval matrix. The rule is applied recursively down to order one. The enclosure of x₁ is that of verifylss. The method works when no enclosure of x₁ in the recursion contains zero and every solve is certified; otherwise the returned ball has radius Inf.

det_cofactor(A::BallMatrix) -> DeterminantResult

The Laplace expansion of the determinant along the first row, evaluated in ball arithmetic. It is not a method of the thesis; its cost is n!, and it is used by interval_det for orders up to three. For order two each entry occurs once in the formula, so the only overestimation is that of the ball product.

contains_zero(det_result::DeterminantResult) -> Bool

Whether zero lies in the enclosure of the determinant. false proves that the interval matrix is regular; true proves nothing.