API - Algorithm Steps

_rump2022a(B::BallMatrix; maxiter = 20, inflate = 0.1, maxlevels = 3)
    -> VerifyEigAllResult

Theorem 2.2 of Rump (2022) with the paper's own transformation, a verified linear solve. Not exported; reached through verifyeigall with method = :rump2022a.

Identical to _rump2022aneumann in every step of the theorem, and differing only in _rump2022a_transform: the correction Δ to the approximate eigenvalue matrix is enclosed by solving WΔ = BW − WX with _rump2010_alg10_7, which is the Krawczyk iteration Rump's verifylss performs, instead of being bounded through an explicit inverse. The solve proves W nonsingular as a by-product, so no separate nonsingularity argument is needed.

_rump2022aneumann(B::BallMatrix; maxiter = 20, inflate = 0.1, maxlevels = 3)
    -> VerifyEigAllResult

Theorem 2.2 of Rump (2022), with the transformation of §2.3 replaced by an explicit inverse and a uniform Neumann bound. Not exported; reached through verifyeigall with method = :rump2022aneumann.

Everything the theorem requires is here: the clustering of the diagonal by connected components, Rtilde of (2.8)-(2.9) as an enclosure, the iteration Y = X_O X_D - E - E X_O with Z = Rtilde .* Y, the self-mapping test (2.10) per cluster with its left side rounded up, and the condition rho(Z) < 1 that Theorem 2.2 needs in addition to (2.10) before the union of the discs may be called the spectrum (Remark 2.4). spectrum_covered reports the conjunction of the two.

Where it departs from the paper. Rump transforms by a verified linear solve, verifylss(W, B*W), which encloses W^{-1}BW and proves W nonsingular at the same time; _rump2022aneumann_transform instead forms R = inv(W) in floating point, bounds R*W - I in the spectral norm, and charges that defect to every entry of the residual enclosure uniformly. The name records the substitution. Its cost is quantitative: tight clusters of size three and above are declined where Rump's Table 4 reports no failures, at 27 of 30 clusters certified for k = 3, because for three eigenvalues separated by u^(1/3) the entries of Rtilde reach 1.6e5 against transformed off-diagonals of order u*cond(W), and the quadratic term of the iteration then puts Z outside X. That figure was measured before the step 6 recursion transformed the uncertified columns again.

No numerical Jordan decomposition is formed. The accuracy attainable is governed by the Jordan structure: an eigenvalue whose largest Jordan block has size k cannot be enclosed more tightly than about u^(1/k) in floating-point arithmetic, so a triple eigenvalue is limited to about 1e-5, and the method declines rather than return a bound it cannot justify. Uncertified clusters carry a Gershgorin radius and Inf in subspace and block.

When a pass leaves clusters uncertified, those columns are transformed again and the pass is repeated, up to maxlevels times, as for _rump2022a. Soundness does not rest on that refinement: Theorem 2.1's assertions hold "for any quality" of the approximation, and the certification is (2.10) alone.

_rump2022a_conditioned(B::BallMatrix; kappa_max = 1/√eps, maxiter = 20, inflate = 0.1,
                       maxlevels = 3) -> VerifyEigAllResult

Theorem 2.2 of Rump (2022) with the frame chosen by the condition number of the eigenvector matrix W of the midpoint: W when cond(W) ≤ kappa_max, the unitary factor of a Schur decomposition otherwise. The same rule is applied in the recursion of step 6: the columns left uncertified are transformed by the eigenvectors of their diagonal block when the condition number of that eigenvector matrix is at most kappa_max, and are left as they are otherwise. A deviation from Rump's algorithm, named for it; reached through verifyeigall with method = :rump2022aconditioned.

cond is computed in floating point and only selects the frame: every claim of the result is that of Theorem 2.2 on the enclosure of the transformed matrix, whatever frame was selected. The effect of the rule is that the condition number of the similarity returned, which block_resolvent_floor divides by, is a product of factors at most kappa_max; the price is that clusters which only an ill-conditioned frame separates stay uncertified.

_rump2022a_discclusters(B::BallMatrix; maxiter = 20, inflate = 0.1, maxlevels = 3)
    -> VerifyEigAllResult

Theorem 2.2 of Rump (2022) with the paper's transformation, but clustering the diagonal at √eps·‖A‖_∞ instead of the 1e-14·‖A‖_∞ of its step 2. A deviation from Rump's algorithm, named for it; reached through verifyeigall with method = :rump2022adiscclusters.

The reason is that the fixed threshold cannot group a multiple eigenvalue. A k-fold eigenvalue computed in floating point spreads by about u^(1/k), and the paper states the case: "this attempt to construct a matrix with 3-fold eigenvalue generates a matrix with a cluster of radius 1e-5". Against a threshold of 1e-14‖A‖ those three eigenvalues stay three singletons, and a singleton cannot satisfy (2.10) for a defective eigenvalue, so the algorithm declines on exactly the case Theorem 2.2 was written to handle. Miyajima (2014a) groups on overlap of the Gershgorin discs of the rows, whose radii follow the residual; this variant does not have the row radii and widens the diagonal discs by the k = 2 sensitivity √eps·‖A‖ instead.

Soundness does not depend on the choice. The partition in Theorem 2.2 is arbitrary and (2.10) is what proves the inclusion, so a different clustering changes how often the test succeeds and never whether a success is valid. The cost is that well-separated eigenvalues closer than √eps‖A‖ are merged into one block, which widens their discs to the block's.

_rump2022a_pencil(A::BallMatrix, B::BallMatrix; frame = :eigen,
                  cluster_rule = _rump2022a_clusters, maxiter = 20, inflate = 0.1,
                  maxlevels = 3) -> VerifyEigAllResult

Theorem 2.2 of Rump (2022) for the pencil A x = λ B x: the theorem is applied to a ball matrix containing W⁻¹B⁻¹AW, which has the eigenvalues of the pencil. An extension of Rump's algorithm, which is stated for one matrix; reached through verifyeigall(A, B; method) with method one of :rump2022a, :rump2022adiscclusters, :rump2022aschur, :rump2022aschurstep6.

The steps after the transformation are those of _rump2022a, the recursion of step 6 included, since they act on the transformed matrix alone. The transformation is _rump2022a_pencil_transform: the residual AW − BWD is enclosed with error-free transformations and (BW)Δ = AW − BWD by a verified solve with the ball matrix BW, which proves B nonsingular. B⁻¹A is not formed. W is the matrix of generalised eigenvectors (frame = :eigen) or the factor Z of a generalised Schur decomposition Q*AZ = S, Q*BZ = T (frame = :schur), for which Z⁻¹B⁻¹AZ = T⁻¹S is triangular and Z unitary.

In the result, subspaces[i] and blocks[i] of a certified cluster satisfy A * subspaces[i] = B * subspaces[i] * blocks[i]; similarity and transformed are W and the ball matrix containing W⁻¹B⁻¹AW; the Gershgorin fields are those of that matrix. When the transformation is not certified (a singular or nearly singular B among the causes) nothing is certified and every radius is Inf.

_rump2022a_schur(B::BallMatrix; maxiter = 20, inflate = 0.1) -> VerifyEigAllResult

Theorem 2.2 of Rump (2022) with the paper's transformation and clustering, applied in the frame of a Schur decomposition mid(B) ≈ QTQ* in place of the eigenvector matrix of the paper's step 1. A deviation from Rump's algorithm, named for it; reached through verifyeigall with method = :rump2022aschur.

The transformation encloses Q⁻¹BQ by the same verified solve, which proves Q nonsingular and does not use that it is unitary. The enclosed matrix is then close to the triangular T and not to a diagonal matrix: its part above the diagonal enters E of Theorem 2.2 at full size. The theorem does not ask E to be small, its conclusions hold for the clusters on which the self-mapping test (2.10) succeeds, so every claim of the result is the theorem's; what changes is how often (2.10) succeeds. The recursion of step 6 is not run: it would transform the uncertified columns by the eigenvectors of their block, which is the frame this variant avoids.

The reason for the variant is the similarity it returns. κ₂(Q) is one up to rounding where the eigenvector matrix of a non-normal matrix is ill conditioned, and block_resolvent_floor divides by that condition number.

_rump2022a_schur_step6(B::BallMatrix; maxiter = 20, inflate = 0.1, maxlevels = 3)
    -> VerifyEigAllResult

_rump2022a_schur followed by the recursion of step 6 as _rump2022a runs it: where a pass leaves clusters uncertified, their columns are transformed by the eigenvectors of their diagonal block and Theorem 2.2 is applied again, up to maxlevels times and for as long as the number of certified columns grows. A deviation from Rump's algorithm, named for it; reached through verifyeigall with method = :rump2022aschurstep6.

The similarity returned is then Q multiplied by those transformations, the identity outside the columns that were transformed again; it is unitary up to rounding only when the recursion did not run, and its condition number is what block_resolvent_floor bounds from similarity.

Enclosure of B*W - W*X with B a ball matrix and W, X floating point, the midpoints in doubled precision by _accurate_product_sum.

It stands where Rump's Algorithm transform calls prodK(B, W, -W, X), described there as computing "an accurate approximation of the residual BW - WX using error-free transformations", with a second output giving the radius. The paper attributes prodK to no numbered reference; it is a routine of INTLAB, whose algorithm is the scheme of [16], a manuscript. For Float64 data that scheme is used here, _lange_rump_two_terms in src/types/MMul/lange_rump.jl; for other types, and where it declines, the accumulation is Dot2 of [11].

_rump2022a_prodK_pencil(A::BallMatrix, B::BallMatrix, W::Matrix, D::Vector)

Enclose the residual A W − B W D of a pencil, D the diagonal of the approximate eigenvalues, as _rump2022a_prodK encloses B W − W X.

The product of three factors is brought to products of two by error-free transformations of W D: _lr_diagonal_pieces! returns four matrices whose sum is W D, exact unless a product underflows, that being two_product of the real and imaginary parts. The residual is then a sum of five products, enclosed by _accurate_product_sum; the interval parts rad(A)|W| and rad(B)|W||D| are added, together with |B| times eight times the smallest subnormal, for the underflow of the four two_product of an entry of W D.

_rump2022a_transform(B::BallMatrix, W::Matrix, X0::Matrix)

Enclose W⁻¹ B W as [31] writes the transformation, returning the enclosure and the bound on the spectral radius that the inner solve reports, or nothing and that bound when the solve does not certify.

The correction is enclosed by a verified linear solve and not bounded through an explicit inverse: with X the approximate eigenvalue matrix,

W⁻¹ B W  =  X + W⁻¹(B W − W X),      so      W Δ = B W − W X,

and verifylss(BallMatrix(W), B W − W X) encloses Δ while proving W nonsingular, which is Theorem 10.8 of [28] applied to the point matrix W. Without that proof there is no enclosure of W⁻¹ B W, so the transformation declines rather than returning a bound. The right hand side is the _rump2022a_prodK enclosure, so the two terms that nearly cancel are handled by the error-free transformations and not by ball arithmetic, and the Newton step on X is taken in floating point, since it has only to be a better approximation.

_rump2022aneumann_transform(B::BallMatrix, W::Matrix, X0::Matrix)

Enclose W⁻¹ B W through an explicit approximate inverse of W, returning the enclosure and the defect ‖R W − I‖₂, or nothing and that defect when it is not below one.

What is enclosed is the correction and not the product: with X the approximate eigenvalue matrix,

W⁻¹ B W  =  X + W⁻¹(B W − W X),

and the residual B W − W X is small, so an enclosure of W⁻¹ times it is small as well, while an enclosure of W⁻¹ B W formed directly would carry the radius of the product. One Newton step on X, taken in floating point since nothing in it is certified, makes the residual smaller still; the residual itself is enclosed by _rump2022a_prodK, through the error-free transformations rather than in ball arithmetic, since its two terms nearly cancel.

With R an approximate inverse of W one has W⁻¹ B W = (R W)⁻¹(R B W), so an enclosure of R(B W − W X) together with ‖R W − I‖ < 1 gives one of the correction by a Neumann series, the extra radius being ‖R W − I‖/(1 − ‖R W − I‖) times the norm of that enclosure. The factor multiplies the size of the residual and not the size of B, which is what makes the radius small. _rump2022a_transform replaces the inverse and this series by a verified linear solve, as [31] writes the transformation.

_rump2022a_pencil_transform(A::BallMatrix, B::BallMatrix, W::Matrix, X0::Matrix)

The transformation of _rump2022a_transform for a pencil A x = λ B x, returning the enclosure of W⁻¹ B⁻¹ A W and the bound the inner solve reports, or nothing and that bound when the solve does not certify.

With D the diagonal matrix of the approximate eigenvalues, read off X0,

W⁻¹ B⁻¹ A W  =  D + (B W)⁻¹(A W − B W D),      so      (B W) Δ = A W − B W D,

and verifylss applied to the ball matrix B W encloses Δ while proving B W, hence both B and W, nonsingular, which is Theorem 10.8 of [28] for an interval matrix. The residual is enclosed by _rump2022a_prodK_pencil.

The Newton step that _rump2022a_transform takes on the eigenvalue matrix is not taken here: it would make the matrix in the third factor full, and the residual is brought to products of two factors only for a diagonal one.

_vbd_unitary_basis(result) -> Bool

Trait distinguishing VBD results whose basis is genuinely unitary (so that adjoint(basis) is its inverse) from those whose basis is only block-orthonormal (so that consumers must use inv(basis)). SchurGershgorinResult (NSD / Schur / Hermitian eigenvectors) is unitary; Miyajima2014aSchurNewtonResult is not.

_rumplange2023(A::BallMatrix; kappa = 0) -> VerifyEigAllResult

_rumplange2023_eig with its result in the fields of a VerifyEigAllResult; reached through verifyeigall for a HermitianBallMatrix, of which A is the dense ball matrix that contains it. The statements are for the Hermitian matrices of A.

• clusters is the partition of Theorem 4.2; centers[i] and radii[i] describe the smallest interval containing the intervals of cluster i, which holds exactly length(clusters[i]) eigenvalues; certified is true where that interval is bounded, and spectrum_covered when all are.
• gershgorin_centers and gershgorin_radii hold the n intervals λ_j of the theorem, one for each eigenvalue: every eigenvalue lies in their union and a union of k of them that is disjoint from the others holds k eigenvalues, which is the contract of those fields.
• subspaces[i] contains an orthonormal basis Q_i of the invariant subspace of cluster i (Theorem 6.2), and blocks[i] the Hermitian matrix Q_i*AQ_i, with A Q_i = Q_i blocks[i].
• similarity contains the unitary matrix Q made of those bases, and transformed the matrix Q*AQ: block diagonal on the clusters, exact zeros elsewhere.

_miyajima2014a_alg1(A::BallMatrix, B::BallMatrix) -> VerifyEigAllResult

Algorithms 1 and 2 of Miyajima (2014): verified inclusions of all eigenvalues and invariant subspaces of the pencil A x = λ B x. Not exported; reached through verifyeigall with method = :miyajima2014a.

B is not assumed nonsingular; its nonsingularity, and that of X̃ and Y, is proved by ‖t‖_∞ < 1 with t = |I − YBX̃|𝟙 (Theorem 3.1, Corollary 3.2). Once that holds, every eigenvalue of the pencil lies in one of the n discs ⟨λ̃ᵢ, rᵢ⟩ with r = u + ⟨u⟩_t t, so spectrum_covered reports that single condition and not the per-cluster tests: a cluster of k connected discs contains precisely k eigenvalues by Remark 3.3 whether or not Algorithm 2 then succeeds on it.

certified[i] is the stronger, per-cluster claim. For a singleton it is Theorem 3.5: the disc holds exactly one eigenvalue, of geometric multiplicity one, and subspaces[i] encloses an eigenvector for it. For a cluster of size k it is Theorem 3.11 via Algorithm 2: radii[i] is ρ(V ᵀP̄̄_ε), subspaces[i] encloses a basis of the invariant subspace, and blocks[i] encloses the block Λ ∈ ⟨λ̃I_k, V ᵀP̄̄_ε⟩.

Where Algorithm 2 declines, radii[i] falls back to the smallest disc about λ̃ containing the union of the cluster's Gershgorin discs, which is still a rigorous enclosure of those k eigenvalues, while subspaces[i] and blocks[i] carry Inf, since no subspace was proved.

The result also carries similarity, the matrix X̃, and transformed, a ball matrix containing X̃⁻¹B⁻¹AX̃ = D̃ + S. Here S = (I + R₂)⁻¹R₁ with the paper's two residuals R₁ = Y(AX̃ − BX̃D̃) and R₂ = YBX̃ − I, and it is enclosed through a floating-point solution of (I + R₂)S̃ = R₁ and the residual of that solve. This is an addition to the paper's algorithm, which works with the norms of the rows of R₁; the discs and the clusters above do not use it.

Reference

S. Miyajima, Fast Enclosure for All Eigenvalues and Invariant Subspaces in Generalized Eigenvalue Problems, SIAM J. Matrix Anal. Appl. 35(3):1205-1225, 2014, doi 10.1137/140953150.

miyajima2010_algorithm1_steps1to3(A::BallMatrix, B::BallMatrix)

Steps 1 to 3 of Algorithm 1 of [5] for the pencil A x = λ B x, with the approximate eigenvectors X̃ and eigenvalues taken from eigen(mid(A), mid(B)):

• Step 1: Z = ⟨Z_c, Z_r⟩ with B X̃ ∈ Z, by Algorithm 2 (cprod, _cprod in src/types/MMul/oishi_mmul.jl) for the bounds of the real and of the imaginary part and Algorithm 3 (ccr, _ccr in the same file) for the ball;
• Step 2: Y ≈ Z_c⁻¹, in floating point;
• Step 3: ε₂ ≥ ‖R₂‖_∞, R₂ = Y B X̃ − I, by inequality (8) of Corollary 1 (bdd_R2).

Returns (X, Y, Z, eigenvalues, coupling_bound) with coupling_bound = ε₂. Step 4 of the algorithm is the test ε₂ < 1, which proves that B, X̃ and Y are nonsingular; Steps 5 and 6, the bound ε₁ ≥ ‖R₁‖_∞ by Procedure 1 or 2 of the paper and ε = ε₁/(1 − ε₂), are not done here. (This function was named miyajima_algorithm_1_procedure_1; Procedure 1 of the paper is the enclosure of R₁, which it does not compute.)

The paper states the algorithm for point matrices. Here B may have a positive radius: the bounds of Step 1 are then those of ⟨mid(B), rad(B)⟩ X̃, by the kernel ccrprod′ of Procedure 2 (_ccrprod_prime, same file), so that B̃ X̃ ∈ Z and ‖Y B̃ X̃ − I‖_∞ ≤ ε₂ for every B̃ of the ball, which is our use of the paper's kernel and not a statement of the paper. A enters through the approximate eigenvectors only.

bdd_R2(Y::AbstractMatrix, Z::BallMatrix)

An upper bound of ‖Y Z̃ − I‖_∞ over the matrices Z̃ of the ball Z = ⟨Z_c, Z_r⟩, for a floating-point matrix Y: inequality (8) of Corollary 1 of [5],

‖R₂‖_∞ ≤ ‖fl(Y Z_c − I)‖_∞ + ‖(|Y_r| + |Y_i|) Z_r e‖_∞ + γ′ₙ₋₁ ‖|Y| |Z_c| e‖_∞
          + u (‖fl(Y Z_c)‖_∞ + 1) + 4 n² (1 + γₙ₋₁) u̲,

with Y_r, Y_i the real and imaginary parts of Y, e the vector of ones, u the unit roundoff, u̲ the underflow unit, γₙ = n u/(1 − n u) and γ′ₙ₋₁ = √5 u + γₙ₋₁ (1 + √5 u). The one product Y Z_c is formed in round-to-nearest, as the corollary asks, and its rounding is what the last three terms pay for; it holds in the presence of underflow.

Differences from the paper, all on the side of a larger bound: eps is used for u, which is twice the unit roundoff, and γₙ for γₙ₋₁. The two terms that the corollary states as exact quantities, ‖(|Y_r| + |Y_i|) Z_r e‖_∞ and ‖|Y| |Z_c| e‖_∞, are evaluated with every operation rounded up.

_up_bound_Linf_opnorm(A::AbstractMatrix{T}) where T

An upper bound of the ∞-norm of the matrix A, of type float(real(T)): the largest row sum of the moduli, each modulus and each sum rounded up.

_up_bound_Linf_norm(v::AbstractVector{T}) where T

An upper bound of the ∞-norm of the vector v, of type float(real(T)): the largest modulus, rounded up.

The default u of certify_eigenpair when the caller passes none: the eigenvector of A_mid' whose eigenvalue is closest to conj(λ), from eigen of the adjoint in floating point.

Nothing is verified here and nothing need be. u enters the certificate only through the normalisation uᴴw = 1 and the bordered Jacobian, and the statement holds for any u with uᴴṽ ≠ 0; a poor one makes q₀ and y larger, and a u whose |uᴴṽ| falls below 100 eps makes certify_eigenpair decline. The approximate eigenvector argument is not read.

The (N+1) × (N+1) ball matrix standing for J = F′(x₀) of certify_eigenpair, from which that routine takes q₀ ≥ ‖I − CJ‖₂:

J = [ −ṽ   A − λ̃ I ]
    [  0     uᴴ     ]

the bordered Jacobian of the eigenpair system, whose nonsingularity at a zero is the simplicity of the eigenvalue there (Lemma 2 of [33], credited there to Wilkinson). The first column and the last row are exact, so their radii are zero, and the A − λ̃ I block carries the radii of A_ball and nothing else: its midpoints are mid(A_ball) − λ̃ I in round-to-nearest, so the rounding of the N diagonal subtractions, at most eps(|aᵢᵢ − λ̃|)/2 each, is not in the radii.

An upper bound of ‖M̃‖₂ over the matrices M̃ of the ball M, by the method certify_eigenpair passes in its norm_method: :svd goes through the O(n³) verified factorisation, anything else through the cheaper bound built from the row and column sums. Both q₀ and ‖C‖₂ of the certificate are taken here, so a looser method costs a larger radius and never correctness.

The NKEigenpairResult returned when certify_eigenpair declines before any bound has been computed: certified = false, the eigenvalue and every eigenvector component a ball of infinite radius about the pair handed in, each of enclosure_radius, defect_q, C_norm and residual_y infinite and discriminant = -1.

The two cases that reach it are |uᴴṽ| < 100 eps, where the normalisation cannot be imposed, and a midpoint Jacobian that inv could not factor. Nothing is claimed about the spectrum of A.

The NKEigenpairResult returned when certify_eigenpair has computed q₀, ‖C‖₂ and y and one of its two tests then failed: q₀ ≥ 1, or q₀ < 1 with a negative discriminant. The three numbers are kept, so that a caller can see which test failed and by how much; certified is false, enclosure_radius infinite and the eigenvalue and eigenvector balls have infinite radius about the pair handed in.

discriminant is recomputed here by _nk_discriminant, a lower bound of (1 − q₀)² − 4‖C‖₂y; it is a diagnostic only, and for q₀ ≥ 1 it bounds nothing.

_refine_schur_impl(A, Q0, T0, max_iterations, tol)

The iteration behind refine_schur_decomposition: Algorithm 4 of [7], the mixed-precision refinement, returning a SchurRefinementResult. It throws an ArgumentError on a T0 that is not upper triangular, the real Schur form of a real matrix with a complex pair among them, for which the iteration diverges.

The arithmetic is floating point in the high precision hp, so the two numbers reported are measured on the computed defects and bound nothing about A itself: residual_norm and orthogonality_defect are the sizes of the computed stril(QᴴAQ) and QᴴQ − I, divided in the first case by a Frobenius norm of A, each size being an upper bound of the spectral norm for BigFloat and a Frobenius norm from _frobenius_norm otherwise. A rigorous statement about the refined pair comes from certify_schur_decomposition.

The steps:

1. Initial orthogonalization: Apply Newton-Schulz to Q₀ converted to high precision
2. Main iteration loop: a. Compute T̂ = Q^H A Q in high precision b. Extract E = stril(T̂), T = T̂ - E (upper triangular part) c. Solve triangular matrix equation stril(TL - LT) = -E for L d. Construct skew-Hermitian correction W = L - L^H e. Compute orthogonality defect Y = Q^H Q - I f. Update Q with combined Newton-Schulz step: Q_new = ½ Q (2I + 2W - Y - YW + W² + W³)
3. Convergence: Stop when ‖E‖F and ‖Y‖F are below tolerance

The algorithm achieves quadratic convergence and typically needs only 3-4 iterations to refine from double to quadruple precision.

References

• [7] Bujanović, Kressner & Schröder, "Iterative refinement of Schur decompositions", Numer. Algorithms 92, 247–267 (2023). Algorithm 4: Mixed-precision Schur refinement.

_solve_triangular_equation_direct!(L::Matrix{T}, T_mat::Matrix{T}, E::Matrix{T};
                                   max_entry::Real=1e5) where T

Solve the triangular matrix equation stril(TL - LT) = -E for strictly lower triangular L.

This implements Algorithm 2 from Bujanović et al. (2022): successive substitution with bottom-to-top columnwise order.

Arguments

• L: Output matrix (will contain strictly lower triangular solution)
• T_mat: Upper triangular matrix with pairwise distinct diagonal entries
• E: Strictly lower triangular right-hand side matrix
• max_entry: Maximum allowed absolute value for entries of L (default: 1e5). Larger values are truncated to zero to prevent instability from clustered eigenvalues.

Details

For 1 ≤ j < i ≤ n, the entry L[i,j] is computed as:

L[i,j] = -1/(T[i,i] - T[j,j]) * (E[i,j] + Σ_{k=i+1}^n T[i,k]*L[k,j] - Σ_{k=1}^{j-1} L[i,k]*T[k,j])

where entries to the left or below (i,j) must be computed first.

Notes

When eigenvalues are clustered (nearly equal), the denominator T[i,i] - T[j,j] becomes very small, causing large entries in L that can destabilize the algorithm. Following the suggestion in Bujanović et al. (Example 4), entries larger than max_entry are set to zero.

References

• [7] Algorithm 2: Successive substitution

_solve_triangular_equation_recursive!(L::Matrix{T}, T_mat::Matrix{T}, E::Matrix{T};
                                      nmin::Int=32, max_entry::Real=1e5) where T

Solve stril(TL - LT) = -E using the recursive block algorithm (Algorithm 3).

This achieves better cache efficiency for large matrices by recursively partitioning into 2×2 blocks and using Sylvester equation solvers for off-diagonal blocks.

Arguments

• L: Output matrix (will contain strictly lower triangular solution)
• T_mat: Upper triangular matrix with pairwise distinct diagonal entries
• E: Strictly lower triangular right-hand side matrix
• nmin: Minimum block size before switching to direct algorithm (default: 32)
• max_entry: Maximum allowed absolute value for entries of L (default: 1e5)

Algorithm

Partition T, E, L into 2×2 blocks: T = [T11 T12] E = [E11 0 ] L = [L11 0 ] [ 0 T22] [E21 E22] [L21 L22]

1. Solve Sylvester equation: T22 * L21 - L21 * T11 = -E21
2. Update: E11 ← E11 + stril(T12 * L21)
3. Update: E22 ← E22 - stril(L21 * T12)
4. Recursively solve for L11 from stril(T11L11 - L11T11) = -E11
5. Recursively solve for L22 from stril(T22L22 - L22T22) = -E22

References

• [7] Algorithm 3: Recursive block algorithm

_solve_sylvester_triangular!(X::Matrix{T}, T22::Matrix{T}, T11::Matrix{T}, C::Matrix{T}) where T

Solve the triangular Sylvester equation T22 * X - X * T11 = C where T11 and T22 are upper triangular (from Schur form).

This corresponds to equation (15) in Bujanović et al. (2022): T22 * L21 - L21 * T11 = -E21

Arguments

• X: Output matrix (will contain solution), size (m × k)
• T22: Upper triangular matrix (m × m)
• T11: Upper triangular matrix (k × k)
• C: Right-hand side matrix (m × k)

References

• [7] Section 2.1, equation (15)

_stril(A::AbstractMatrix)

A new matrix of the size of A whose entries below the diagonal are those of A and whose other entries are zero: the operator written stril throughout [7], in which the refinement step is the equation stril(TL − LT) = −E. Entries are copied, not rounded.

_frobenius_norm(A::AbstractMatrix)

The Frobenius norm of A, accumulated in real(eltype(A)) with round-to-nearest arithmetic, so the value is an approximation and not a bound in either direction. It serves the stopping tests of the refinement iterations in this file, where the arithmetic may be BigFloat or a MultiFloat and the rigorous norms of norm_bounds/ are unavailable.

_refine_symmetric_eigen_impl(A, Q0, λ0, max_iterations, tol)

The iteration behind refine_symmetric_eigen: Algorithm 1 of [34], named RefSyEv there, run at the precision in force, which the caller sets.

A, Q0 and λ0 are converted to hp (BigFloat at the current precision, or the type given), and the steps R = I − QᵀQ, S = QᵀAQ, the Rayleigh quotients λᵢ = sᵢᵢ/(1 − rᵢᵢ), the threshold δ, the correction Ẽ and the update Q(I + Ẽ) are taken in that arithmetic. It returns a SymmetricEigenRefinementResult.

Everything here is floating point in the high precision, with no directed rounding: the residual_norm and orthogonality_defect it reports are Frobenius norms computed by _frobenius_norm and are approximations, not bounds. A refined pair that needs a bound is handed to certify_eigenpair, which takes an approximate pair and returns an enclosure.

_apply_block_structure(T::BallMatrix, clusters, structure::Symbol)

The block matrix that rigorous_block_schur reports as its T, obtained from the transformed matrix by keeping the blocks the requested structure asks for, where clusters are the index ranges of the blocks:

• :diagonal keeps the blocks (i, i) and sets every other entry, midpoint and radius, to zero;
• :quasi_triangular keeps the blocks (i, j) with i ≤ j and zeroes the rest;
• :full returns the argument unchanged.

Discarding a block is not an approximation that is left unaccounted: what A is similar to is measured against the matrix returned here, through the residual R₁ = Q_inv(AQ − QT) of RigorousBlockSchurResult, so a discarded block reappears in projected_residual_norm and in perturbation_norm.

_compute_off_diagonal_norm(T::BallMatrix, clusters)

The largest of the rigorous upper bounds of ‖T̃[μᵢ, μⱼ]‖₂, i ≠ j, over the matrices T̃ of the ball T, with μᵢ the index ranges in clusters; zero when there is one block only. Each block bound is collatz_upper_bound_L2_opnorm, and it is the maximum of the bounds rather than a bound on the maximum of the norms, which for a single matrix are the same thing.

This is off_diagonal_norm of RigorousBlockSchurResult, a measure of how far the block form is from block diagonal, and it is read off the structured T of _apply_block_structure: under :diagonal those blocks have been set to zero, so it comes back zero and says nothing about the matrix before the truncation.

_rump2010_alg10_7(A::BallMatrix{T}, b::BallVector{T};
                  r=0.1, ϵ=1e-20, iter_max=20) where {T<:AbstractFloat}

Algorithm 10.7 of Rump (2010), the algorithm behind INTLAB's verifylss: an enclosure of the solution of the square linear system $Ax = b$ by epsilon-inflation. Not exported; reached through verifylss with method = :rump2010.

Input

• A – square matrix of size n × n
• b – vector of length n or matrix of size n × m
• r – relative inflation, default 10%
• ϵ – absolute inflation, default 1e-20
• iter_max – maximum number of iterations

Output

Returns an VerifyLssResult containing:

• solution – enclosure of the solution of the linear system
• certified – Boolean flag, if true, then solution is certified to contain the true solution; if false, certification failed (check diagnostics)
• iterations – number of iterations performed
• spectral_radius_bound – upper bound on ‖I - RA‖₂ (must be < 1 for convergence)
• condition_number – approximate condition number of A

The result can be destructured as (x, cert) = verifylss(A, b).

Algorithm

This is Algorithm 10.7 of

S. M. Rump, Verification methods: rigorous results using floating-point arithmetic, Acta Numerica 19 (2010) 287-449, doi 10.1017/S096249291000005X,

which Rump attributes to Rump (1980) and which is the algorithm behind INTLAB's verifylss, transcribed step for step:

R  = inv(mid(A))                     # approximate inverse
xs = R * mid(b)                      # approximate solution
C  = I - R * A                       # iteration matrix, R against the BALL A
Z  = R * (b - A * xs)
X  = Z
for iter in 1:15
    Y = X * [0.9, 1.1] + [-1e-20, 1e-20]
    X = Z + C * Y
    in0(X, Y) && return xs + X
end

Two points decide whether it proves anything, and both are easy to get wrong.

C is I minus R times the ball A, not I - R*mid(A) in floating point: the product R*mid(A) is itself rounded, and an enclosure of I - RA has to carry that rounding. On a 50x50 randn matrix with an exact midpoint, the float form BallMatrix(I - R*mid(A), abs.(R)*rad(A)) claims radius 0.0 while the entries of I - RA differ from its midpoint by up to 4.9e-15, so it is not an enclosure and nothing downstream of it is proved. This routine previously did that.

The containment test is in0, strict interior containment, per Rump's footnote 19: in0(X, Y) checks X ⊂ int(Y) componentwise. in, which allows equality, does not give the contraction the argument rests on.

What success proves

Theorem 10.8: if the iteration ends successfully then every matrix in the ball A is nonsingular, and the returned enclosure contains $Σ(A, b) = {x : Ax = b for some A ∈ A, b ∈ b}$. Nonsingularity is a conclusion, not an assumption, which is why Rump uses this to enclose W⁻¹BW while proving W invertible.

The note following Theorem 10.8 records that the result holds for a matrix right-hand side $b ∈ IR^{n×k}$ as well, which is the method taken by the BallMatrix right-hand side below.

Convergence

An inclusion is computed essentially if and only if $ρ(|I − RA|) < 1$; spectral_radius_bound reports $‖I − RA‖₂$, an upper bound for it, and a value $≥ 1$ means the iteration cannot contract.

Example

A = BallMatrix(randn(5, 5))
b = BallVector(randn(5))
result = verifylss(A, b)
if result.certified
    println("Certified solution: ", result.solution)
else
    println("Verification failed:")
    println("  Spectral radius bound: ", result.spectral_radius_bound)
    println("  Condition number: ", result.condition_number)
end

_rumpogita2024_cholesky(A::BallMatrix) -> (; G) or nothing

Section 4 of Rump and Ogita (2024): an inclusion of the Cholesky factor. For every Hermitian matrix Ã of the ball A, Ã is proved positive definite and its Cholesky factor, the upper triangular G with positive diagonal and G*G = Ã, lies in the ball matrix G returned. Positive definiteness is a conclusion and not an assumption. Returns nothing when the verification fails. The midpoint of A must be Hermitian and its radii symmetric.

The method

G̃ is a floating-point Cholesky factor of the midpoint and X_G ≈ G̃⁻¹, upper triangular. I_E := X_G* Ã X_G is a perturbed identity, enclosed with à X_G kept as an unevaluated sum of two terms; Algorithm 3.2 encloses its factors, I_E = L_E U_E. I_E being Hermitian, U_E = D L_E* with D the diagonal of U_E, which is real; when D > 0 is proved, G_E = D^{1/2} L_E* is the Cholesky factor of I_E, and by uniqueness, (4.1) of the paper,

G = G_E X_G⁻¹ = D^{1/2} L_E* X_G⁻¹.

X_G⁻¹ is enclosed through G̃ X_G = I + F in two-fold precision and Algorithm 3.2 for (I + F)⁻¹, as in _rumpogita2024_lu.

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.

_rumpogita2024_lu(A::BallMatrix) -> (; L, U, p, q) or nothing

Section 3 of Rump and Ogita (2024): inclusions of the factors of the LU decomposition of an m × n matrix. For every matrix Ã of the ball A, the permuted matrix Ã[p, q] has a unique decomposition LU with L of size m × min(m, n) unit lower triangular (trapezoidal) and U of size min(m, n) × n upper triangular (trapezoidal), L and U in the ball matrices returned. Returns nothing when the verification fails; nothing is then proved.

p is the row permutation of a floating-point decomposition with partial pivoting. q is the identity for m ≥ n; for m < n it is the column permutation of a decomposition with partial pivoting of the transpose, as in Section 3.4, after which the left square block has a decomposition.

The method

With p = min(m, n), B the leading p × p block of Ã[p, q], L̃, Ũ its floating-point factors, X_L ≈ L̃⁻¹ unit lower and X_U ≈ Ũ⁻¹ upper triangular,

I_E := X_L B X_U = L_E U_E

is a perturbed identity, enclosed with the product B X_U kept as an unevaluated sum of two terms and the second product in two-fold precision; Algorithm 3.2 encloses L_E, U_E and their inverses. Then, by uniqueness of the decomposition,

• for m ≥ n, (3.8) and (3.9): U = U_E X_U⁻¹ and L = Ã[p, :] X_U U_E⁻¹, the latter from the two-term product of the whole matrix with X_U;
• for m < n, (3.10): L = X_L⁻¹ L_E and U = L_E⁻¹ X_L Ã[p, q], with X_L Ã[p, q] kept in two terms.

These are the methods the paper finds best in each case (the blue curves of its Figures 3, 5 and 7), at O(max(m, n) min(m, n)²) operations. The inverses X_U⁻¹ and X_L⁻¹, which the paper leaves to the implementation, are enclosed through X_U Ũ = I + F, respectively X_L L̃ = I + F, in two-fold precision and Algorithm 3.2 for (I + F)⁻¹, which also proves X_U and X_L nonsingular. The entries that are zero or one in the factors are exact.

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, Sections 3.2, 3.3 and 3.4.

block_resolvent_floor(vbd) -> BlockResolventFloor

Prepare, from a verified block diagonalisation vbd of A (the result of miyajima2014a_schurnewton or of schur_gershgorin_enclosure, for the standard eigenvalue problem; for a result of verifyeigall see the method below), a lower bound of σ_min(A − zI) that costs O(number of blocks) at each z. For a ball matrix A the bound holds for every matrix of the ball.

Setting

W is the basis of vbd, M := W⁻¹AW, and Λ = blockdiag(P_j) its block diagonal part, with blocks on the index sets C_j. The result carries β_Λ ≥ ‖M − Λ‖₂ (block_residual_norm, from M − Λ = (I + R₂)⁻¹R₁, R₁ = Y(AW − WΛ), R₂ = YW − I), the radii r_j ≥ ‖(M − Λ)[C_j, :]‖₂ (block_coupling), the barycentres c_j and n_j ≥ ‖P_j − c_j I‖₂.

The two bounds

With s_j(z) ≤ σ_min(P_j − zI) and s(z) = min_j s_j(z):

1. σ_min(M − zI) ≥ s(z) − β_Λ, since σ_min(Λ − zI) = min_j σ_min(P_j − zI) and σ_min is 1-Lipschitz in the spectral norm (|σ_min(X + E) − σ_min(X)| ≤ ‖E‖₂).
2. If every s_j(z) > 0 and τ(z)² := Σ_j (r_j / s_j(z))² < 1, then σ_min(M − zI) ≥ s(z)(1 − τ(z)): factor M − zI = (Λ − zI)(I + G), G = (Λ − zI)⁻¹(M − Λ); the block row j of G is (P_j − zI)⁻¹(M − Λ)[C_j, :], of norm at most r_j/s_j, and a matrix whose block rows have norms g_j has norm at most √(Σ g_j²). Each block is charged its own coupling, so this bound is the larger one where the coupling is uneven across the blocks.

Both are carried to A by A − zI = W(M − zI)W⁻¹:

σ_min(A − zI) ≥ σ_min(M − zI) / κ₂(W).

κ₂(W) is bounded here, not taken from vbd.kappa (a diagnostic there): with Y the floating-point inverse of W and ρ ≥ ‖YW − I‖₂ < 1, ‖W⁻¹‖₂ ≤ ‖Y‖₂/(1 − ρ).

s_j(z) is |z − c_j| − n_j (from P_j − zI = (c_j − z)I + (P_j − c_j I)), or, on request, the larger of that and a certified smallest singular value of the ball P_j − zI.

The statements are Lemma "Block resolvent floor" and Proposition "Localized block-dominance floor" of I. Nisoli, A certified resolvent bound via verified block diagonalization (unpublished draft, 2026); the arguments are the ones above.

block_resolvent_floor(r::VerifyEigAllResult) -> BlockResolventFloor

The same bound from a result of verifyeigall, with the frame W := S, the similarity the algorithm ended with, and the blocks on its clusters.

S is a product of floating-point matrices, contained in the ball matrix r.similarity, and r.transformed is a ball matrix containing M = S⁻¹AS, so the quantities of the bound are read from it and no residual is formed: P_j is its block on the cluster C_j, c_j the barycentre of the midpoints of the diagonal of that block, n_j ≥ ‖P_j − c_j I‖₂, r_j ≥ ‖(M − Λ)[C_j, :]‖₂ the norm of the block row with the block itself set to zero, and β_Λ ≥ ‖M − Λ‖₂ the norm of the ball matrix with every diagonal block set to zero. κ₂(S) is bounded over the ball r.similarity as above.

The bound uses the enclosure of S⁻¹AS and a partition, and neither the self-mapping test (R̃ ∘ Y)V_i ⊆ int(XV_i) of Theorem 2.2 of [31], which is what proves the Jordan blocks and the invariant subspaces, nor r.spectrum_covered: it holds whether or not the clusters are certified. When the transformation failed, r.transformed is A and the frame is the identity.

For method = :miyajima2014a the matrix r.transformed is the one of _miyajima2014a_alg1. For a result of the pencil A x = λ B x it contains S⁻¹B⁻¹AS, and the bound is then a lower bound of σ_min(B⁻¹A − zI).

The enclosure of S⁻¹AS and the clustering are those of [31], where S⁻¹AS is enclosed by a verified linear solve that also proves S nonsingular; the bound read off it is the one of the method above.

is_M_matrix(A::BallMatrix) -> Bool

true when every matrix of the ball A is proved to be a nonsingular M-matrix, by the definition of Varga (2004), Definition 5.4 (repeated as Definition C.3): a real matrix A with a_ij ≤ 0 for all i ≠ j, written A = μI − B with μ real and B ≥ 0 entrywise, is a nonsingular M-matrix if ρ(B) < μ. Then A⁻¹ ≥ 0 entrywise (Proposition C.4 there).

The two conditions are checked for the whole ball: the upper ends of the off-diagonal entries are not positive, and with μ the largest upper end of a diagonal entry, an upper bound of the spectral radius of a matrix that dominates every B entrywise is below μ (collatz_upper_bound). false means that this was not proved, and is returned for a complex matrix, an M-matrix being real.

For the property that needs no sign condition, that the comparison matrix is a nonsingular M-matrix, see is_H_matrix.

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]), Definition 5.4 and Appendix C.

collatz_upper_bound(A::BallMatrix; iterates = 10)

A rigorous upper bound on the spectral radius of |M| for every M in the ball A, hence on the spectral radius of every such M, by Collatz's inclusion: for a nonnegative matrix B,

ρ(B) ≤ max_i (Bx)_i / x_i      for every positive vector x,

applied to B =upper_abs(A) with x chosen by iterates power iterations from the vector of ones. Any number of iterations gives a valid bound; the result is Inf, not NaN, when the quotient is not finite (the vector is kept strictly positive, so a zero row or a nilpotent matrix is handled).

References

L. Collatz, Einschließungssatz für die charakteristischen Zahlen von Matrizen, Math. Z. 48 (1942) 221-226, doi 10.1007/BF01180013.

S. M. Rump, Verified bounds for singular values, in particular for the spectral norm of a matrix and its inverse, BIT Numer. Math. 51 (2011) 367-384, doi 10.1007/s10543-010-0294-0, where the inclusion is used in the form (3.3) for BᵀB; see collatz_upper_bound_L2_opnorm.