API - Eigenvalues & SVD

RigorousEigenvaluesResult

Container returned by rigorous_eigenvalues bundling the midpoint eigenvector factorisation, the balls whose union contains the eigenvalues, and the norm bounds underpinning their verification. Besides behaving like the underlying vector of eigenvalue balls, the struct exposes the interval residual, the projected residual used in the Miyajima certification, and the inverse defect (inverse * vectors - I).

rigorous_eigenvalues(A::BallMatrix)
rigorous_eigenvalues(H::Hermitian)

Compute a set of balls whose union contains the eigenvalues of A, following Theorem 1 of Miyajima (2010), below, with B = I. The returned RigorousEigenvaluesResult exposes both the balls and the norm bounds used during certification.

What is certified. The returned balls have the computed eigenvalues as centres and one common radius, and every eigenvalue of every matrix in the input ball lies in their union. A single ball need not contain an eigenvalue, and the number of eigenvalues in a ball or in a group of overlapping balls is not certified. For sets with a certified number of eigenvalues use verifyeigall.

The type of the argument decides the kind of ball:

• for a BallMatrix the balls are discs of the complex plane, with complex centres, also when the matrix is real and the computed eigenvalues are real: matrices in the input ball may have complex eigenvalues, which lie in the discs;
• for a HermitianBallMatrix, the Hermitian matrices within the tolerances, and for a point matrix that is Hermitian by its type (Hermitian(M), or Symmetric(M) with M real), the eigenvalues are real and the balls are real intervals: every eigenvalue of every Hermitian matrix within the tolerances lies in their union.

References

S. Miyajima, Fast enclosure for all eigenvalues in generalized eigenvalue problems, J. Comput. Appl. Math. 233(11) (2010) 2994-3004, doi 10.1016/j.cam.2009.11.048, Theorem 1 ([5]): for an arbitrary Y, with R₁ = Y(AX̃ − BX̃D̃) and R₂ = YBX̃ − I, if ‖R₂‖_∞ < 1 then every eigenvalue is within ε = ‖R₁‖_∞/(1 − ‖R₂‖_∞) of one of the diagonal entries of D̃.

evbox(A::BallMatrix)
evbox(H::Hermitian)

Backward-compatible wrapper returning only the vector of balls produced by rigorous_eigenvalues: their union contains the eigenvalues, a single ball need not contain one.

VerifyEigAllResult{T, CT}

Outcome of verifyeigall.

Fields

• clusters::Vector{Vector{Int}}: the partition mu of 1:n the algorithm worked with, one entry per cluster of approximate eigenvalues; a cluster of length > 1 is a multiple or nearly multiple eigenvalue whose individual eigenvectors are ill-posed.
• certified::Vector{Bool}: which clusters satisfied the self-mapping test (2.10).
• centers::Vector{CT}: the approximate eigenvalue lambda_i of each cluster.
• radii::Vector{T}: for a certified cluster rho(|V_i' Z V_i|) of Theorem 2.2; for one where the self-mapping test (2.10) declined, the smallest disc about centers[i] containing the Gershgorin connected component of the transformed matrix that the cluster's rows fall in, which is a valid enclosure but proves no Jordan structure. Either way every eigenvalue of cluster i lies in centers[i] ± radii[i], and certified[i] says which of the two claims is being made.
• subspaces::Vector{BallMatrix}: Inf unless the cluster is certified, since no subspace follows from the Gershgorin fallback. For each certified cluster, an enclosure of a basis of the corresponding invariant subspace of the INPUT matrix, n by |mu_i|, satisfying B * subspaces[i] = subspaces[i] * blocks[i]; radius Inf where not certified.
• blocks::Vector{BallMatrix}: the enclosed Jordan block Mhat_i of each certified cluster.
• basis::Matrix{CT}: the approximate eigenvector matrix W the transformation used.
• spectrum_covered::Bool: true when every cluster is certified, in which case the union of the discs contains the whole spectrum of the input.
• iterations::Int, transform_defect::T: the number of interval iterations, and the certified ‖I − RW‖ of the transformation, which must be below one.
• gershgorin_centers::Vector{CT}, gershgorin_radii::Vector{T}: the Gershgorin row discs of the transformed ball matrix, which encloses W⁻¹BW and so carries the eigenvalues of B; those of B itself when the transformation failed; empty when the result comes from another method. Every eigenvalue lies in their union, and a union of discs disjoint from the others holds as many eigenvalues as it has discs.
• similarity::BallMatrix, transformed::BallMatrix: a ball matrix containing the similarity S the algorithm ended with (W, multiplied by the transformations of the recursion of step 6 when it ran) and the ball matrix containing S⁻¹BS to which Theorem 2.2 was applied, for every matrix B of the input. When the transformation failed they are the identity and B. For :miyajima2014a they are X̃ and a ball matrix containing X̃⁻¹B⁻¹AX̃ (X̃⁻¹AX̃ for one matrix), empty when the method did not prove its basis nonsingular. block_resolvent_floor reads them.

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

Verified inclusions of all eigenvalues and invariant subspaces of B. Returns a VerifyEigAllResult, whose spectrum_covered says whether the union of the returned discs is proved to contain the whole spectrum. Where a cluster is not certified, its radius is a Gershgorin bound and its subspace and block are Inf; certified says which is which. When the transformation itself cannot be certified, every eigenvalue index gets the Gershgorin disc of B about its diagonal entry and nothing is certified.

maxlevels bounds the recursion of step 6: where a pass leaves clusters uncertified, those columns are transformed again by the eigenvectors of their diagonal block and Theorem 2.2 is applied to the result. The paper describes this step in one sentence and gives no code; the implementation is the algorithm's own first step applied to those columns. maxlevels = 0 switches it off.

This function selects an algorithm and, where it leaves the spectrum uncovered, runs the one named by fallback: the result of method is returned when its spectrum_covered is true; otherwise fallback is run with the same keywords and its result is returned when it covers the spectrum or certifies more columns, the result of method when it does neither. Which of the two was returned can be read from similarity. Called without method, the function runs :rump2022a with the fallback :rump2022aschurstep6; called with a method, it runs that algorithm alone unless a fallback is given, so that a result asked for by name is that algorithm's. fallback = nothing runs :rump2022a alone, and no fallback is run after :miyajima2014a, whose spectrum_covered is another condition (below). Each algorithm is a separate unexported function named for the paper it implements, so that a deviation from a paper is visible in the name rather than buried in a docstring:

For a Hermitian matrix the algorithm of Rump and Lange (2023) is reached by the type of the argument, verifyeigall(A) with A a HermitianBallMatrix or Hermitian(M); it is not a method of a plain ball matrix.

The Neumann variant is kept because it does not need the solve to succeed, so it still returns a result where the solve declines; where both succeed, the faithful one is at least as tight. :miyajima2014a called with one matrix solves A x = λ x; for a pencil see verifyeigall(A, B).

Note that spectrum_covered is not the same condition in the two families. Rump's needs the self-mapping test on every cluster together with ρ(Z) < 1; Miyajima's needs only ‖t‖_∞ < 1, after which the discs cover the spectrum whether or not the per-cluster tests succeed. In both, certified[i] is the per-cluster claim about the invariant subspace.

Example

r = verifyeigall(BallMatrix(randn(50, 50)))
if r.spectrum_covered
    for i in eachindex(r.clusters)
        println(r.centers[i], " ± ", r.radii[i])     # eigenvalue disc
        Y = r.subspaces[i]                            # basis of the invariant subspace of B
    end
end

Reference

S. M. Rump, Verified error bounds for all eigenvalues and eigenvectors of a matrix, SIAM J. Matrix Anal. Appl. 43(4):1736-1754, 2022, doi 10.1137/21M1451440.

verifyeigall(A::HermitianBallMatrix; kappa = 0) -> VerifyEigAllResult
verifyeigall(H::Hermitian; kappa = 0)

Verified inclusions of all eigenvalues and of orthonormal bases of the invariant subspaces of a Hermitian matrix, by the algorithm verifyeigall of Rump and Lange (2023); see _rumplange2023 for the fields of the result. For a HermitianBallMatrix the statements hold for every Hermitian matrix within the tolerances, as in the paper ("the results are true for every symmetric or Hermitian A ∈ A"), and for no other matrix of the entrywise ball. kappa is the clustering threshold.

Reference

S. M. Rump and M. Lange, Fast computation of error bounds for all eigenpairs of a Hermitian and all singular pairs of a rectangular matrix with emphasis on eigen- and singular value clusters, J. Comput. Appl. Math. 434 (2023) 115332, doi 10.1016/j.cam.2023.115332.

verifyeigall(A::BallMatrix, B::BallMatrix; method = :miyajima2014a, kwargs...)
    -> VerifyEigAllResult

Verified inclusions of all eigenvalues and invariant subspaces of the pencil A x = λ B x.

B is not assumed nonsingular: every method proves it, and returns an uncertified result when it cannot. B⁻¹A is formed by none of them. :rump2022aneumann is not implemented for a pencil. No fallback is run. The keywords maxiter, inflate and maxlevels are those of the Rump methods for one matrix; :miyajima2014a takes none.

inclusion_discs(r::VerifyEigAllResult) -> (centers, radii, counts, ok)

Discs holding all the eigenvalues of every matrix of the ball passed to verifyeigall, each with the number of eigenvalues it holds when it is disjoint from the others.

When every cluster is certified (r.spectrum_covered) the discs are those of Theorem 2.2 of [31], centers[i] ± radii[i], and the i-th holds as many eigenvalues as the i-th cluster has indices: that theorem gives one Jordan block of that size per cluster, and together the blocks carry the whole spectrum with algebraic multiplicity.

Otherwise the discs are the Gershgorin row discs of the transformed matrix W⁻¹BW, each holding one eigenvalue. Every eigenvalue is then in their union, and a union of discs disjoint from the remaining ones holds exactly as many eigenvalues as it has discs, by Theorem 1.6 of [27], applied to each matrix of the ball W⁻¹BW separately.

ok is false, with centers and radii empty, when the clusters are not all certified and no Gershgorin discs were recorded either.

eigencount_outside(r::VerifyEigAllResult, R) -> (count, ok)

Number of eigenvalues of modulus greater than R, with algebraic multiplicity, from the inclusion_discs of verifyeigall: the sum of the counts of the discs outside the circle, when no disc meets it. Each disc is then inside or outside, so the discs outside are a union of components disjoint from the rest and hold exactly that many eigenvalues. ok is false, and count zero, when no discs are available or some disc meets the circle |z| = R.

eigencount_in_disc(r::VerifyEigAllResult, c, ρ) -> (count, ok)

Number of eigenvalues in the open disc |z − c| < ρ, with algebraic multiplicity, from the inclusion_discs of verifyeigall, as in eigencount_outside; ok is false when no discs are available or some disc meets the circle |z − c| = ρ.

AlmostInvariantBasis{T, CT}

An orthonormal basis for one certified cluster of a VerifyEigAllResult, with certified defects. Produced by orthonormal_invariant_basis.

Fields

• cluster::Vector{Int}: the indices of the cluster this basis belongs to.
• basis::BallMatrix: Q, n × k, orthonormal to within orthogonality_defect.
• block::BallMatrix: the Rayleigh block H, k × k, as an exact floating-point matrix with zero radius. It is a candidate, not an enclosure of anything: the certified statement is invariance_defect, a bound on ‖BQ − QH‖₂ for this H.
• orthogonality_defect::T: a rigorous bound on ‖Q*Q − I‖₂.
• invariance_defect::T: a rigorous bound on ‖BQ − QH‖₂. This is the measure of how far Q is from spanning an invariant subspace of B; it is zero exactly when Q spans one.

orthonormal_invariant_basis(B::BallMatrix, r::VerifyEigAllResult)
    -> Vector{AlmostInvariantBasis}

For every certified cluster of r, an orthonormal basis of its subspace together with certified bounds on how far that basis is from orthonormal and from invariant. One entry per certified cluster, in cluster order; clusters where the self-mapping test declined are skipped, since there is no subspace to orthonormalise.

Theorem 2.2 of [31] returns the basis in the frozen-rows normalisation Vᵢᵀ Ŷᵢ = I_k, whose columns are not orthonormal; measured on three semisimple doubles, cond(mid Y) was 2.6e3, 2.6e3 and 1.2e4. Consumers of the basis — a spectral projector, the smallest singular value of a block, a resolvent bound through the block — degrade with that conditioning, while the eigenvalue radius does not depend on it at all.

The orthonormal Q is obtained by a floating-point QR of mid(Ŷᵢ) and is a candidate: nothing here claims it spans the invariant subspace. What is rigorous is the triple

‖Q*Q − I‖₂,        H = fl(Q*BQ),        ‖BQ − QH‖₂,

the last computed by prodK, the error-free-transformation residual, against the input B, and it is zero exactly when Q spans an invariant subspace of B. Taking H as a floating-point candidate rather than a ball matters: a ball H carries a radius of order u‖B‖, which enters the defect multiplied by ‖Q‖ = 1 and dominated it, overstating the defect by 22 to 78 times for a simple eigenvalue. With prodK and a float H the bound agrees with a 256-bit reference to three digits. Turning invariance_defect into a distance to a true invariant subspace needs a separation between the cluster and the rest of the spectrum as well, and is not done here.

Example

r = verifyeigall(B)
for b in orthonormal_invariant_basis(B, r)
    println(length(b.cluster), "  ", b.invariance_defect)
end

NKEigenpairResult{T, CT}

Result of Newton-Kantorovich certification for a single eigenpair.

Fields

• certified::Bool: whether the eigenpair was rigorously certified
• eigenvalue::Ball{T,CT}: certified eigenvalue enclosure
• eigenvector::Vector{Ball{T,CT}}: certified eigenvector enclosure
• enclosure_radius::T: NK ball radius in the product space
• defect_q::T: q₀ = ‖I - C·J‖₂
• C_norm::T: ‖C‖₂
• residual_y::T: y = ‖C·F(x₀)‖₂
• discriminant::T: (1-q₀)² - 4·‖C‖₂·y

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.

NKEigenpairsResult{T, CT}

Returned by certify_eigenpairs: one NKEigenpairResult per eigenpair attempted, in the order of the indices passed, with the number of them that were certified.

Fields

• results::Vector{NKEigenpairResult{T,CT}}: the individual results, certified or not
• n_certified::Int: how many of them have certified == true
• n_total::Int: how many were attempted, that is length(results)

length, getindex, iterate, firstindex and lastindex are forwarded to results, so the object can be indexed and iterated in place of the vector.

The count is called n_certified since version 0.4 and was called n_verified in 0.3.0; the old name still resolves, with a deprecation warning.

certify_eigenpair(A, λ_approx, v_approx; u_approx=nothing, norm_method=:fast)

An enclosure of one eigenpair of A near the approximate pair (λ̃, ṽ), valid for every matrix of A when it is a ball matrix.

Only a simple eigenvalue can be certified here, and that is a property of the statement rather than of the implementation: simplicity is part of the conclusion below, obtained from the invertibility of F′ at the zero, so at a multiple eigenvalue q₀ < 1 cannot hold and the result comes back uncertified however accurate the pair handed in. For a spectrum with clusters use verifyeigall, which encloses each cluster with its multiplicity rather than each eigenpair, and reports spectrum_covered when the union of its discs holds the whole spectrum.

What is proved

Let u be a vector with uᴴṽ = 1 up to rounding (an approximate left eigenvector), and F(z, w) := (A w − z w, uᴴ w − 1) on ℂ × ℂᴺ with the Euclidean norm. With x₀ := (λ̃, ṽ), J := F′(x₀), a matrix C (the floating-point inverse of the midpoint of J) and

y ≥ ‖C F(x₀)‖₂,        q₀ ≥ ‖I − C J‖₂,        c ≥ ‖C‖₂,

write d := (1 − q₀)² − 4 c y. If q₀ < 1 and d ≥ 0, then with

r := 2y / ((1 − q₀) + √d)

there is exactly one (z, w) with ‖(z − λ̃, w − ṽ)‖₂ ≤ r, A w = z w and uᴴ w = 1; and z is a simple eigenvalue of A. When certified, the result has eigenvalue = ⟨λ̃, r⟩ and each component of eigenvector of radius r about that of ṽ.

That radius is Yamamoto's: [32], Theorem 1, writes it 2‖ε‖/(1 − ‖K‖ + √t) with the same discriminant. Until 2026-11-07 this returned 2y/(1 − q₀), the same numerator over a smaller denominator and hence up to twice as large; both are valid, and the proof below covers either.

Proof. F is quadratic: F(x₀ + h) = F(x₀) + J h + Q(h) with Q(h) = (−h_z h_w, 0), so ‖Q(h)‖₂ = |h_z| ‖h_w‖₂ ≤ ‖h‖₂²/2 and ‖Q′(h)‖₂ ≤ ‖h‖₂. For T(x) := x − C F(x),

‖T(x₀ + h) − x₀‖₂ ≤ y + q₀ ‖h‖₂ + c ‖h‖₂²/2,        ‖T′(x₀ + h)‖₂ ≤ q₀ + c ‖h‖₂.

The first bound is at most ρ at ‖h‖₂ ≤ ρ exactly when (c/2)ρ² − (1 − q₀)ρ + y ≤ 0, so every ρ between the two roots of that quadratic serves, and both r above and the larger 2y/(1 − q₀) lie between them: the larger because 4cy ≤ (1 − q₀)² gives (c/2)(2y/(1−q₀))² − y = y(2cy/(1−q₀)² − 1) ≤ −y/2 ≤ 0, and r ≤ 2y/(1 − q₀) because √d ≥ 0, while r is at least the smaller root because d understates the quadratic's own discriminant (1 − q₀)² − 2cy. For either, c ρ ≤ 2cy/(1 − q₀) ≤ (1 − q₀)/2, so the second bound is at most (1 + q₀)/2 < 1. So T maps the ball of radius r into itself and is a contraction on it; its fixed point is the only zero of C F there, and C is invertible because ‖I − CJ‖₂ < 1. At the zero the derivative of F is invertible, by the second bound, which excludes a second independent eigenvector w′ (then (0, w′ − (uᴴw′) w) would be in its kernel) and a generalised eigenvector w′ with (A − z) w′ = w (then (1, w′ − (uᴴw′) w) would be).

y, q₀ and c are computed in ball arithmetic over the ball A, so the conclusion holds for each of its matrices, with an eigenpair that depends on the matrix.

References

The construction is Yamamoto's, and the statement above is self-contained, so neither paper is needed to check it.

• [33] writes the eigenpair problem as a nonlinear system and applies an existence theorem for its zeros, for a real matrix and a real pair, normalising by (1 − ‖x‖₂²)/2. Its Lemma 2, credited to Wilkinson, is that an eigenvalue is simple exactly when the bordered Jacobian is nonsingular, which is where the restriction above comes from. Its Remark 1 sets aside the rounding committed in forming I − CJ and C F(x₀), which is what ball arithmetic supplies here.
• [32] drops the realness, normalises by fixing the largest coordinate of the computed eigenvector rather than by the Euclidean length, and carries the discriminant and the enclosure radius used here. Its bounds are componentwise in the maximum norm, where this returns one Euclidean radius for the whole pair.

Arguments

• A: square matrix (plain or BallMatrix)
• λ_approx::Number: approximate eigenvalue
• v_approx::AbstractVector: approximate eigenvector
• u_approx: optional left eigenvector (computed automatically if nothing)
• norm_method::Symbol: :fast uses upper_bound_L2_opnorm, :svd uses svd_bound_L2_opnorm

The enclosed eigenvector is the one normalised by uᴴ w = 1, with u the vector given, or computed, divided by the conjugate of uᴴ ṽ.

Returns

NKEigenpairResult with certification status and enclosure.

certify_eigenpairs(A; indices=nothing, norm_method=:fast)

Certify all (or selected) eigenpairs of matrix A.

Arguments

• A: square matrix (plain or BallMatrix)
• indices: optional vector of eigenvalue indices to certify
• norm_method::Symbol: :fast or :svd

The approximate eigenpairs come from the eigendecomposition of the midpoint, the Hermitian one when the midpoint is Hermitian; each is then certified by certify_eigenpair, for every matrix of the ball.

Returns

NKEigenpairsResult containing individual results plus summary counts.

VerifySvdAllResult{T, NT}

Outcome of verifysvdall for an m × n matrix, p = min(m, n).

• values::Vector{Ball{T, T}}: p intervals; there is a numbering of the singular values with σ_j ∈ values[j].
• clusters::Vector{Vector{Int}}: a partition of 1:p; the unions of the intervals of the clusters are mutually disjoint and the union for a cluster holds exactly as many singular values as the cluster has indices.
• U::BallMatrix (m × p), V::BallMatrix (n × p): for each cluster μ, U[:, μ] and V[:, μ] contain orthonormal bases of the left and of the right singular subspace of the singular values of μ. A radius is Inf where this was not proved.

verifysvdall(A::BallMatrix; kappa = 0) -> VerifySvdAllResult

Verified inclusions of all singular values and of the left and right singular subspaces of an m × n matrix: Algorithm verifysvdall of Rump and Lange (2023), with its refinement refinesvd and the singular vector bounds of its Section 7. The statements, Theorem 7.3 there, are those of VerifySvdAllResult and hold for every matrix of the ball A. kappa ≥ 0 is the paper's threshold: singular value inclusions closer than kappa relative to their size are put in one cluster, which widens their intervals and tightens their subspaces.

The algorithm

For m ≥ n, with X̃ Σ̃ Ỹ* a floating-point economy-size decomposition of the midpoint and B = [0 A*; A 0], whose eigenvalues are ±σ_j and m − n zeros:

1. σ̃_j is the midpoint of an enclosure of the Rayleigh quotient of B at [ỹ_j; x̃_j], 2 Re(x̃_j*Aỹ_j)/(‖x̃_j‖² + ‖ỹ_j‖²), with the products in twice the working precision.
2. E ∋ AỸ − X̃Σ̃, F ∋ A*X̃ − ỸΣ̃, and δ_j = √(‖E e_j‖² + ‖F e_j‖²)/‖ỹ_j‖: the interval max(0, σ̃_j ± δ_j) holds a singular value (Theorem 7.1).
3. The intervals are grouped, with kappa, into connected components, and for a component μ of more than one index δ_j becomes √(‖E(:, μ)‖² + ‖F(:, μ)‖²)/σ_min(Ỹ(:, μ)), until the partition no longer changes.
4. For a cluster of one index the interval is intersected with the one of Theorem 5.1 applied to B, the gap including the distance to the eigenvalues −σ_i and, for m > n, to zero.
5. For a cluster μ with gap ε to the other intervals, the left and the right singular subspaces are within τ = √(‖E(:, μ)‖² + ‖F(:, μ)‖²)/ε of X̃(:, μ) and Ỹ(:, μ) (Lemma 6.1 applied to B; Lemma 7.2 for the cluster of the smallest singular values), and Lemma 3.1 gives orthonormal bases. For m > n the gap for the left subspace of the cluster of the smallest singular values also includes its distance to zero.

For m < n the algorithm is applied to A* and the two sides are exchanged.

The departures from the paper's listing are those of _rumplange2023_eig: the stopping test of step 3, the centres and the single ε of a cluster in step 5, and the cluster of all singular values, for which the listing returns radius zero. The treatment of the left null space for m > n, which the paper describes and leaves out of its listing, is not implemented: U has min(m, n) columns.

Reference

S. M. Rump and M. Lange, Fast computation of error bounds for all eigenpairs of a Hermitian and all singular pairs of a rectangular matrix with emphasis on eigen- and singular value clusters, J. Comput. Appl. Math. 434 (2023) 115332, doi 10.1016/j.cam.2023.115332, Section 7.

_rumplange2023_eig(A::BallMatrix; kappa = 0)
    -> (; lo, hi, clusters, vectors, basis, centers, radii)

Algorithm verifyeigall of Rump and Lange (2023) for a Hermitian matrix, with the refinement refineeig of its Section 5 and the eigenvector bounds of its Section 6. For every Hermitian matrix Ã of the ball A:

• Theorem 4.2. There is a numbering λ_1, …, λ_n of the eigenvalues of Ã with λ_j ∈ [lo[j], hi[j]]. clusters is a partition of 1:n; the unions of the intervals of the clusters are mutually disjoint, and the union for a cluster holds exactly as many eigenvalues as the cluster has indices.
• Theorem 6.2. For each cluster μ, the ball matrix vectors[:, μ] contains an orthonormal basis of the invariant subspace of the eigenvalues in the intervals of μ.

basis is the floating-point eigenvector matrix X̃, centers the approximate eigenvalues λ̃ (midpoints of the Rayleigh quotients) and radii the δ of the algorithm before refinement.

The algorithm

1. X̃ are the eigenvectors of (mid(A)* + mid(A))/2; the Rayleigh quotients ϱ_j = x̃_j*Ãx̃_j / x̃_j*x̃_j are enclosed with the products accumulated to twice the working precision, and λ̃_j is the midpoint of the enclosure (Section 5).
2. E ∋ ÃX̃ − X̃Λ̃, and δ_j = ‖E e_j‖/‖x̃_j‖, so that λ̃_j ± δ_j holds an eigenvalue (Theorem 4.1 with k = 1).
3. The intervals are grouped into the connected components of "these two intersect", or, with the threshold kappa > 0 of Section 6, "these two are closer than kappa relative to their ends": eigenvalues that the algorithm could separate are then kept as one cluster, whose invariant subspace is enclosed much more tightly than the single eigenvectors would be. For a component μ with more than one index, δ_j is replaced for j ∈ μ by ‖E(:, μ)‖/σ_min(X̃(:, μ)) (Theorem 4.1 with k = |μ|), the intervals keeping their own midpoints, and the grouping is repeated until the partition no longer changes.
4. For a cluster of one index, with e_j the gap of its interval to the others, |λ_j − ϱ_j| ≤ ‖Ãx̃_j − ϱ_j x̃_j‖²/(e_j‖x̃_j‖²) (Theorem 5.1), and the interval is intersected with that one.
5. For a cluster μ, with ε a lower bound of the distance from the λ̃_j, j ∈ μ, to the intervals outside μ, there is Y with columns in the invariant subspace and ‖X̃(:, μ) − Y‖ ≤ τ := ‖E(:, μ)‖/ε (Lemma 6.1), and then an orthonormal basis Q of that subspace with ‖Q − X̃(:, μ)‖ ≤ ‖I − X̃(:, μ)*X̃(:, μ)‖ + √2 τ (Lemma 3.1), which bounds every entry.

Where the code departs from the paper's listing

Three places, each toward the statement of the lemma used:

• the loop of step 3 stops when the partition is the one for which the δ were computed, where the listing compares the number of components;
• in step 5 the distance is measured from λ̃_j, the value in the residual E, where the listing takes the midpoint of the refined interval; and one ε is used for the whole cluster, as (6.3) and the sentence after (6.4) say, where the listing divides column by column;
• when all eigenvalues form one cluster the listing returns radius zero for X̃; here the radius is ‖I − X̃*X̃‖, Lemma 3.1 with Y = X̃, since X̃ is orthonormal only up to rounding.

Reference

S. M. Rump and M. Lange, Fast computation of error bounds for all eigenpairs of a Hermitian and all singular pairs of a rectangular matrix with emphasis on eigen- and singular value clusters, J. Comput. Appl. Math. 434 (2023) 115332, doi 10.1016/j.cam.2023.115332, Sections 3 to 6.

RigorousGeneralizedEigenvaluesResult

Container returned by rigorous_generalized_eigenvalues bundling the midpoint eigenvector factorisation, the balls whose union contains the eigenvalues, and the norm bounds underpinning their verification. Besides behaving like the underlying vector of eigenvalue balls, the struct exposes the interval residual, the projected residual used in the Miyajima certification, and the coupling defect of the left action (left_action * B * right_vectors - I).

rigorous_generalized_eigenvalues(A::BallMatrix, B::BallMatrix)

Compute a set of balls whose union contains the eigenvalues of the generalised problem A * x = λ * B * x, by Theorem 1 of Miyajima (2010), below. The returned RigorousGeneralizedEigenvaluesResult exposes both the balls and the norm bounds used during certification.

What is certified. The returned balls have the computed eigenvalues as centres and one common radius, and every eigenvalue of every pair of matrices in the input balls lies in their union. A single ball need not contain an eigenvalue, and the number of eigenvalues in a ball or in a group of overlapping balls is not certified. For sets with a certified number of eigenvalues use verifyeigall.

The balls are discs of the complex plane, with complex centres, also when A and B are real and the computed eigenvalues are real: pairs of matrices in the input balls may have complex eigenvalues, which lie in the discs.

References

S. Miyajima, Fast enclosure for all eigenvalues in generalized eigenvalue problems, J. Comput. Appl. Math. 233(11) (2010) 2994-3004, doi 10.1016/j.cam.2009.11.048, Theorem 1 ([5]): for an arbitrary Y, with R₁ = Y(AX̃ − BX̃D̃) and R₂ = YBX̃ − I, if ‖R₂‖_∞ < 1 then every eigenvalue is within ε = ‖R₁‖_∞/(1 − ‖R₂‖_∞) of one of the diagonal entries of D̃.

gevbox(A::BallMatrix{T}, B::BallMatrix{T})

Backward-compatible wrapper returning only the vector of balls produced by rigorous_generalized_eigenvalues: their union contains the eigenvalues, a single ball need not contain one.

GEVResult{T}

The result of verify_generalized_eigenpairs for the problem A x = λ B x with A symmetric and B symmetric positive definite, of order n. The eigenvalues λ₁ ≤ … ≤ λₙ and their approximations λ̃₁ ≤ … ≤ λ̃ₙ are both in increasing order.

Fields

• certified::Bool: whether the bounds on the eigenvalues were proved. When it is false every bound below is infinite and message says which condition failed.
• eigenvalue_intervals::Vector{Tuple{T, T}}: (λ̃ᵢ − ηᵢ, λ̃ᵢ + ηᵢ), the ends rounded outwards; the i-th interval contains λᵢ. The intervals of a cluster may overlap.
• eigenvector_centers::Matrix{T}: the approximate eigenvectors x̃⁽ⁱ⁾, as columns.
• eigenvector_radii::Vector{T}: ξᵢ with ‖x̂⁽ⁱ⁾ − x̃⁽ⁱ⁾‖₂ ≤ ξᵢ for some eigenvector x̂⁽ⁱ⁾ of λᵢ (Theorem 7); Inf for the indices where the hypotheses of the theorem were not proved, which includes every eigenvalue that is not separated from its neighbours.
• beta::T: β ≥ √‖B⁻¹‖₂.
• global_bound::T: δ with |λᵢ − λ̃ᵢ| ≤ δ for all i (Theorem 4).
• individual_bounds::Vector{T}: εᵢ with min_j |λⱼ − λ̃ᵢ| ≤ εᵢ (Theorem 5).
• separation_bounds::Vector{T}: the radii ηᵢ of the intervals, by (39): min(εᵢ, δ) where the eigenvalue nearest to λ̃ᵢ was proved to be λᵢ, and δ elsewhere.
• residual_norm::T: an upper bound of ‖A X̃ − B X̃ D̃‖₂.
• message::String: empty on success.

All the statements hold for every symmetric A and every symmetric B within the tolerances of the arguments.

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.

verify_generalized_eigenpairs(A, B, X̃, λ̃; method) -> GEVResult

Error bounds for all the approximate eigenpairs (λ̃ᵢ, x̃⁽ⁱ⁾) of A x = λ B x, with A real symmetric and B real symmetric positive definite: Algorithm 1 of Miyajima, Ogita, Rump and Oishi (2010). X̃ is the matrix of the approximate eigenvectors and λ̃ the vector of the approximate eigenvalues, in increasing order, as returned by eigen(Symmetric(mid(A)), Symmetric(mid(B))). Positive definiteness of B is proved, not assumed.

That A and B are symmetric is the hypothesis of the theorem, and it is stated by the types of the arguments: each is a point matrix that is Symmetric (or real Hermitian) by its type, or a HermitianBallMatrix, and the bounds then hold for every symmetric A and every symmetric B within the tolerances. A plain matrix or a plain ball matrix is refused: a ball matrix with positive radii contains matrices that are not symmetric. For a pencil without symmetry see miyajima_gev_enclosure and verifyeigall.

With R_g := A X̃ − B X̃ D̃, D̃ := diag(λ̃), and G_g := X̃ᵀ B X̃:

1. β ≥ √‖B⁻¹‖₂, by compute_beta_bound.
2. Theorem 4: if ‖I − G_g‖₂ < 1 then B is positive definite and |λᵢ − λ̃ᵢ| ≤ δ := β ‖R_g‖₂ / (1 − ‖I − G_g‖₂) for all i. Theorem 5: min_j |λⱼ − λ̃ᵢ| ≤ εᵢ := β ‖r⁽ⁱ⁾‖₂ / √gᵢ, with r⁽ⁱ⁾ the i-th column of R_g and gᵢ the i-th diagonal entry of G_g.
3. The radii ηᵢ of (39): min(εᵢ, δ) where Lemma 2 or Lemma 3 proves that the eigenvalue nearest to λ̃ᵢ is λᵢ, and δ elsewhere.
4. Theorem 7: with ρᵢ the smaller of λ̃ᵢ − λ̃ᵢ₋₁ − ηᵢ₋₁ and λ̃ᵢ₊₁ − λ̃ᵢ − ηᵢ₊₁ (the one that exists for i = 1 and i = n) and ξ̂ᵢ := β ‖r⁽ⁱ⁾‖₂ / ρᵢ: if ρᵢ > 0 and ξ̂ᵢ < √gᵢ, there is an eigenvector x̂⁽ⁱ⁾ of λᵢ with ‖x̂⁽ⁱ⁾ − x̃⁽ⁱ⁾‖₂ ≤ ξᵢ := β ξ̂ᵢ.

The result is certified when steps 1 and 2 succeed; the eigenvalue bounds then hold for all i, clusters included. The eigenvector bound is given index by index, and is Inf where the hypotheses of Theorem 7 are not proved.

Every scalar operation is rounded in the direction that keeps the bound.

The two versions

method = :paper is Algorithm 1 as the paper computes it, for A and B in floating point (Float32 or Float64, radius zero); it is the default when both have radius zero in one of those types:

• β by Theorem 10, from a floating-point Cholesky factor L̃ of B, an approximate inverse X_L of it obtained by substitution, and the a priori bounds of Theorems 8 and 9;
• Z̲ = fl▽(B X̃), Z̄ = fl△(B X̃) and the same two products for A X̃, with the rounding mode downward and upward (Technique 1);
• ‖I − G_g‖∞ by the a priori estimate (48), with Z_c and Z_r of (49) and (50), from one product in rounding to nearest (Technique 2);
• δ := β √(‖R_g‖₁ ‖R_g‖∞) / (1 − ‖I − G_g‖∞), formula (26);
• r⁽ⁱ⁾, gᵢ from the enclosures of R_g and of B X̃ already computed (Techniques 3, 4).

One constant differs from the paper's text: its (51) states the backward error of the Cholesky factorisation with γₙ, citing Higham, whose Theorem 10.3 has γₙ₊₁; α_C is computed with γₙ₊₁. The paper does not write the formula of the enclosure of R_g; here column j of B X̃ D̃ is enclosed by the products of λ̃ⱼ with Z̲ and Z̄, rounded outwards.

method = :ball, the default otherwise, gives the same four steps for matrices with tolerances and for any floating-point type: R_g and G_g are enclosed by ball products, the two norms of Theorem 4 are bounded by upper_bound_L2_opnorm, and β comes from compute_beta_bound (Remark 2 of the paper leaves the method for β open). It does not use Technique 2 or Theorem 10.

The input is checked: a DimensionMismatch or an ArgumentError for wrong sizes, complex data or arguments that are not of the symmetric kinds; a result with certified = false when λ̃ is not sorted.

Reference

S. Miyajima, T. Ogita, S. M. Rump and S. Oishi, "Fast verification for all eigenpairs in symmetric positive definite generalized eigenvalue problems", Reliable Computing 14 (2010), 24-45.

compute_beta_bound(B::BallMatrix{T}) -> T

Compute a verified upper bound β ≥ √‖B⁻¹‖₂ via the Miyajima–Rump inversion bound, valid for every matrix in the ball B.

Arguments

• B::BallMatrix: Symmetric positive definite interval matrix

Returns

• β: Upper bound on √‖B⁻¹‖₂ valid for every matrix in the ball B, or Inf when the residual condition fails (B not certifiably nonsingular).

Algorithm (Miyajima–Rump inversion bound)

Take any floating-point Y ≈ B⁻¹ and form the residual R = I - Y·B in ball arithmetic. If ‖R‖₂ < 1 then every matrix in the ball is nonsingular and

‖B⁻¹‖₂ ≤ ‖Y‖₂ / (1 - ‖R‖₂),

so β = √(‖Y‖₂/(1-‖R‖₂)), rounded outward. This is the Neumann/approximate- inverse lemma used throughout the package (see verified_cholesky) and is the same two-residual device as Theorem 1 of Miyajima (2010): Y is free, its quality entering only through ‖R‖, and no inverse is ever formed.

Because R is a ball product it absorbs the radius of B automatically, so the bound covers the whole input interval. The previous implementation followed Theorem 10 with hand-rolled γₙ constants applied to B.c alone: it ignored B.r entirely (returning the same β for B ± 0 and for a ball containing singular matrices) and, when its conditions failed, fell back to sqrt(cond(B.c)) — a non-rigorous estimate of the wrong quantity, since cond = ‖B‖·‖B⁻¹‖ undershoots √‖B⁻¹‖ by √‖B‖ whenever ‖B‖ < 1.

Carrying no unit-roundoff constants, this version is also type-generic rather than Float64-only.

Complexity

O(n³) for the approximate inverse and the ball product.

Example

B = BallMatrix([2.0 0.5; 0.5 2.0], fill(1e-10, 2, 2))
β = compute_beta_bound(B)

miyajima_gev_enclosure(A::BallMatrix, B::BallMatrix, X̃, λ̃; Y = nothing)

Enclose all eigenvalues of the generalized problem Ax = λBx by Theorem 1 of [5].

Given approximate eigenpairs X̃, λ̃ (so that AX̃ ≈ BX̃D̃, D̃ = diag(λ̃)) and an arbitrary matrix Y, set

R₁ := Y(AX̃ - BX̃D̃),    R₂ := YBX̃ - I .

If ‖R₂‖_∞ < 1 then B, X̃ and Y are nonsingular and every eigenvalue λ satisfies minᵢ |λ - λ̃ᵢ| ≤ ε with

ε = ‖R₁‖_∞ / (1 - ‖R₂‖_∞) .

The returned radii are sharper: Corollary 3.2 of [25] replaces the two global norms by row sums u = |R₁|𝟙, t = |R₂|𝟙, giving

radii = u + ‖u‖_t · t,    ‖v‖_w := maxᵢ vᵢ/(1 - wᵢ),

which satisfies radii[i] ≤ ε for every i — usually with room to spare, since a row with a small residual is no longer charged the worst row's.

Y defaults to a floating-point (B X̃)⁻¹ (Remark 2 of the paper). Its accuracy affects only how small ε comes out, never correctness — the same "Y is free" device as the Miyajima–Rump eigenvalue certification: no inverse is ever formed, all conditioning is folded into the single residual ‖R₂‖.

Both residuals are evaluated in ball arithmetic, so the enclosure covers every matrix in the input balls A and B, including all rounding errors.

Unlike verify_generalized_eigenpairs — which implements the symmetric-definite algorithm and returns per-eigenvalue intervals plus eigenvector bounds — this needs no symmetry, no positive definiteness and no √‖B⁻¹‖₂ prefactor; it returns one common radius for all eigenvalues.

Example

A = BallMatrix([4.0 1.0; 1.0 3.0])
B = BallMatrix([2.0 0.5; 0.5 2.0])
F = eigen(Symmetric(mid(A)), Symmetric(mid(B)))
res = miyajima_gev_enclosure(A, B, Matrix(F.vectors), collect(F.values))
res.certified && (res.radius)

MiyajimaGEVResult

Result of miyajima_gev_enclosure.

Fields

• certified::Bool: whether ‖R₂‖_∞ < 1 held, which is the hypothesis of both statements below and without which nothing is proved
• centers::Vector: the approximate eigenvalues λ̃
• radius: one radius good for every eigenvalue, the smaller of ‖R₁‖_∞/(1-‖R₂‖_∞) from Theorem 1 of [5] and maximum(radii)
• radii: the radii of Corollary 3.2 of [25], one for each λ̃ᵢ; every eigenvalue of the pencil lies in ⋃ᵢ {z : |z - λ̃ᵢ| ≤ radii[i]}, for every pair of matrices of the balls A and B, and radii[i] ≤ radius. Infinite when not certified
• nrmR1, nrmR2: upper bounds of the two residual norms ‖R₁‖_∞, ‖R₂‖_∞, computed in ball arithmetic
• message::String: what happened, for a run that was not certified

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.

RigorousSVDResult

Container returned by rigorous_svd bundling the midpoint factorisation, the certified singular-value enclosures, and the block-diagonal refinement obtained from schur_gershgorin_enclosure. Besides the singular values themselves the struct exposes the residual and orthogonality defect bounds that underpin the certification.

rigorous_svd(A::BallMatrix; apply_vbd = true)

Compute a rigorous singular value decomposition of the ball matrix A. The midpoint SVD is certified following Theorem 3.1 of Ref. [10]; optionally, the resulting singular-value enclosure can be refined by applying schur_gershgorin_enclosure to Σ'Σ, yielding a block-diagonal structure with a rigorously bounded remainder.

The returned RigorousSVDResult exposes both the enclosures and the intermediate norm bounds that justify them. When apply_vbd is set to false, the block_diagonalisation field is nothing.

References

• Miyajima S. (2014), "Verified bounds for all the singular values of matrix", Japan J. Indust. Appl. Math. 31, 513–539.
• Rump S.M. (2011), "Verified bounds for singular values", BIT 51, 367–384.

svdbox(A::BallMatrix; method = :auto) -> Vector{Ball}

Verified enclosures of all the singular values of A, one ball each, in the order the chosen theorem produces them.

svdbox is the name Miyajima (2014) gives his Theorem 7. This function selects an enclosure and does nothing else; each enclosure is a separate unexported function named for the theorem it implements, so a deviation from a paper is visible in the name:

These are the paper's own algorithms M3, M1, O and R1 in its numerical section, which labels M1 = Theorem 7, M2 = Theorem 9, M3 = Theorem 10, M4 = Theorem 11, O = Theorem 4, R1 = Theorem 5 and R2 = Theorem 6. Theorem 9, the sharpening of Rump's Theorem 5, has no economy form: Remark 5 shows the upper half of it fails on an economy frame and leaves the improvement as an open challenge.

Prefer the default. Both enclosures it chooses between are anchored to Σ̂ᵢᵢ of the SVD, so the i-th interval encloses σᵢ(A) with no permutation caveat, unlike :rump2011_thm3_1, and neither degrades with conditioning, unlike :miyajima2014_thm11.

For the singular vectors, the residual and the block refinement as well as the values, use rigorous_svd.

Reference

S. Miyajima, Verified bounds for all the singular values of matrix, Japan J. Indust. Appl. Math. 31 (2014) 513-539, doi 10.1007/s13160-014-0145-5.

_rump2011_thm3_1(A::BallMatrix) -> Vector{Ball}

Theorem 3.1 of Rump (2011), cited as Theorem 5 by Miyajima (2014). With U, V square frames from a numerical SVD of A, Σ = UᴴAV = D + E with D diagonal, α = ‖I − UᴴU‖₂ and β = ‖I − VᴴV‖₂ both below one,

(|Dᵢᵢ| − ‖E‖₂)/√((1+α)(1+β))  ≤  σ_ν(i)(A)  ≤  (|Dᵢᵢ| + ‖E‖₂)/√((1−α)(1−β)).

The numbering matters. The theorem asserts the inequalities for some permutation ν which it does not determine, so what is proved relates the two multisets, not the i-th interval to σᵢ. The intervals are returned in decreasing order of |Dᵢᵢ|; reading the i-th of them as an enclosure of σᵢ(A) is justified only when they are pairwise disjoint, which forces the numbering. Where they overlap, what holds is that a union of k overlapping intervals contains k singular values counted with multiplicity. _miyajima2014_thm7 has no such caveat, since its bounds are anchored to Σ̂ᵢᵢ of the SVD itself.

A must be square: Miyajima's Remark 3 shows this split loses the upper bound on an economy frame. Returns Inf radii where the orthogonality condition fails.

_miyajima2014_thm4(A::BallMatrix) -> Vector{Ball}

Theorem 4 of Miyajima (2014), which that paper attributes to Oishi: with U, Σ, V a numerical SVD of A, E = UΣVᵀ − A, F = VᵀV − Iₙ and G = UᵀU − Iₘ, and provided ‖F‖₂ < 1 and ‖G‖₂ < 1,

Σᵢᵢ − δᵢ ≤ σᵢ(A) ≤ Σᵢᵢ + δᵢ,     δᵢ = Σᵢᵢ max(‖F‖₂, ‖G‖₂) + ‖E‖₂.

Both frames are square, so this is the expensive enclosure; Theorem 8 of the same paper proves _miyajima2014_thm7 is never worse. Kept because it is the baseline the paper compares against. Returns Inf radii where the orthogonality condition fails.

_miyajima2014_thm7(A::BallMatrix) -> Vector{Ball}

Theorem 7 of Miyajima (2014), the economy enclosure that the paper itself calls svdbox. With Û, V̂ the first q = min(m,n) columns of a numerical SVD, Σ̂ the q × q diagonal, Ê = Û Σ̂ V̂ᵀ − A, F̂ = V̂ᵀV̂ − I_q, Ĝ = ÛᵀÛ − I_q, and provided ‖F̂‖₂ < 1 and ‖Ĝ‖₂ < 1,

Σ̂ᵢᵢ √((1−‖F̂‖₂)(1−‖Ĝ‖₂)) − ‖Ê‖₂  ≤  σᵢ(A)  ≤  Σ̂ᵢᵢ √((1+‖F̂‖₂)(1+‖Ĝ‖₂)) + ‖Ê‖₂.

Theorem 8 proves these bounds are at least as tight as Theorem 4's, and the residual is measured against A rather than through Rump's split UᵀAV = D + E, which Remark 3 shows can lose the upper bound on an economy frame: for m = 2n, A = [Iₙ; 0], U the block swap and V = Iₙ every truncated residual vanishes while every singular value is 1.

No longer reached by svdbox's default, which routes an exact input to _miyajima2014_thm10 and a ball to _miyajima2014_thm10_weyl; ask for it by name. It remains the default of rigorous_svd on a ball, which also returns vector bounds measured against that ball. Returns Inf radii where the orthogonality condition fails.

_miyajima2014_thm10(A::BallMatrix) -> Vector{Ball}

Theorem 10 of Miyajima (2014), the algorithm his numerical section labels M3, and what svdbox uses on an exact input. It uses the same economy frames as _miyajima2014_thm7 but a one-sided residual. With Û, Σ̂, V̂ the economy SVD, F̂ = V̂ᵀV̂ − I_q and Ĝ = ÛᵀÛ − I_q, and ‖F̂‖₂ < 1, ‖Ĝ‖₂ < 1, define for m ≥ n

Λᵢᵢ = √((1 − ‖Ĝ‖₂)/(1 + ‖F̂‖₂)) Σ̂ᵢᵢ,   Λ̄ᵢᵢ = √((1 + ‖Ĝ‖₂)/(1 − ‖F̂‖₂)) Σ̂ᵢᵢ,
ρ   = ‖A V̂ − Û Σ̂‖₂ / √(1 − ‖F̂‖₂),

and for m < n the same with the roles of F̂ and Ĝ exchanged and ρ = ‖Ûᵀ A − Σ̂ V̂ᵀ‖₂ / √(1 − ‖Ĝ‖₂). Then

Λᵢᵢ − ρ ≤ σᵢ(A) ≤ Λ̄ᵢᵢ + ρ.

The orientation follows the proof: σᵢ(ÛΣ̂V̂⁻¹) ≤ ‖Û‖₂σᵢ(Σ̂)‖V̂⁻¹‖₂, with ‖Û‖₂ ≤ √(1+‖Ĝ‖₂) from ÛᵀÛ = I + Ĝ and ‖V̂⁻¹‖₂ ≤ 1/√(1−‖F̂‖₂); V̂ is square exactly when m ≥ n, which is why the two cases exchange the frames.

Why this is preferred to Theorem 7. The residual is A V̂ − Û Σ̂, which is m × q, where Theorem 7's Ê = Û Σ̂ V̂ᵀ − A is m × n: cheaper, 22Qq² + 12q³ against 24Qq² + 12q³, and in the paper's own measurements at least as tight. Table 2 gives M3 against M1 as 3.0e-14 against 3.1e-14 at cnd = 1, 3.7e-14 against 3.6e-14 at 1e8 and 5.3e-14 against 5.7e-14 at 1e16, with neither degrading as the conditioning worsens, unlike M4 (Theorem 11), which goes to 9.4e-9; Table 3 gives M3 9.5e-15, 4.5e-15, 4.5e-15 against M1's 1.1e-14, 6.0e-15, 6.0e-15.

Returns Inf radii where the orthogonality condition fails.

_miyajima2014_thm10_weyl(A::BallMatrix) -> Vector{Ball}

Theorem 10 of Miyajima (2014) applied to the midpoint of A alone, each interval then widened by one perturbation term that covers the whole ball. A variation on the paper, hence the name: Theorem 10 is used exactly as stated, on an exact matrix, and the radius is handled outside it.

Write A = Ã ± R with à = mid(A) and R = rad(A) ≥ 0 entrywise. For every X ∈ A and every i,

|σᵢ(X) − σᵢ(Ã)| ≤ ‖X − Ã‖₂ ≤ ‖R‖₂.

The first inequality is Weyl's monotonicity theorem applied to the Hermitian dilation H(M) = [0 M; Mᴴ 0], whose eigenvalues are ±σᵢ(M): H is linear, so H(X) − H(Ã) = H(X−Ã), and ‖H(E)‖₂ = ‖E‖₂, so each eigenvalue of H(Ã) moves by at most ‖X−Ã‖₂ and the singular values with them. The second holds because |X − Ã| ≤ R entrywise gives |X−Ã|ᵀ|X−Ã| ≤ RᵀR entrywise, whose spectral radii are then ordered by Perron-Frobenius monotonicity on non-negative matrices. Both inequalities index the singular values in decreasing order, which is the order Theorem 10 anchors its intervals to, so the two compose index by index. ‖R‖₂ is bounded by upper_bound_L2_opnorm.

Mirsky (1960), doi 10.1093/qmath/11.1.50, proves the same inequality for every unitarily invariant norm. Only the spectral-norm case is used here, and that case needs nothing beyond Weyl, so the name says Weyl.

Why this is not the same as handing A to _miyajima2014_thm10. There the residual is AV̂ − ÛΣ̂, so R reaches the bound as the ball product R·V̂, whose radius is bounded entrywise by R|V̂|; taking ‖R|V̂|‖₂ discards the cancellation that makes ‖RV̂‖₂ ≤ ‖R‖₂, and costs a factor growing like √n. Here R never multiplies anything.

The lower bounds are not clamped at zero, so a rank-deficient A returns intervals whose lower end is negative, as _miyajima2014_thm7 and _miyajima2014_thm10 do. Clamping would be sound, since singular values are non-negative, but it would put σ = 0 on the boundary of its interval instead of inside it, which is the signal that the matrix is rank deficient. Returns Inf radii wherever Theorem 10 does on the midpoint.

_miyajima2014_thm11(A::BallMatrix) -> Vector{Ball}

Theorem 11 of Miyajima (2014): the enclosure from one frame only, the eigenvectors V of AᵀA when m ≥ n and of AAᵀ otherwise. With F = VᵀV − I, D̂ + Ê = (AV)ᵀ(AV) and D̂ diagonal, f = |Ê|𝟙 the Gershgorin row sums of the off-diagonal part, and ‖F‖₂ < 1:

for an i whose interval ⟨D̂ᵢᵢ, fᵢ⟩ is isolated from the others,

ρᵢ = min_{j≠i}(|D̂ᵢᵢ − D̂ⱼⱼ| − fⱼ),   gᵢ = ‖Ê eᵢ‖₂²/ρᵢ,   hᵢ = min(fᵢ, gᵢ),
√((D̂ᵢᵢ − hᵢ)/(1+‖F‖₂)) ≤ σᵢ(A) ≤ √((D̂ᵢᵢ + hᵢ)/(1−‖F‖₂)),

the lower bound being 0 when D̂ᵢᵢ < hᵢ; and for an i inside a connected component {i₁,…,i_k} of overlapping intervals,

ξᵢ = max(min_{j∈cluster}(D̂ⱼⱼ − fⱼ), D̂ᵢᵢ − ‖Ê‖_∞),
ξ̄ᵢ = min(max_{j∈cluster}(D̂ⱼⱼ + fⱼ), D̂ᵢᵢ + ‖Ê‖_∞),

with the same division by 1 ± ‖F‖₂ under the square root.

gᵢ is Parlett's bound, Theorem 3 of the paper, |θ − α| ≤ ‖Ay − θy‖₂²/gap(A,θ), so its denominator is ρᵢ and not 2ρᵢ. The cluster branch is what makes the clustered case sound: an individual eigenvalue of D̂ + Ê need not lie in its own Gershgorin interval, only the cluster's eigenvalues in the union of the cluster's intervals, so hᵢ = fᵢ is not available there.

Only V is needed, so this costs less than the two-frame enclosures and says nothing about the singular vectors. Returns Inf radii where ‖F‖₂ < 1 fails.

_svd_auto_theorem(A::BallMatrix) -> Symbol

Theorem 10 either way: on A itself when the input is exact, and through _miyajima2014_thm10_weyl when it carries a radius.

This is not a result of the paper; it is a selection rule. Theorem 10 is the tighter of the two economy enclosures on an exact matrix and the cheaper of the two always, but its residual AV̂ − ÛΣ̂ multiplies A by V̂, so an input radius R reaches the bound as a ball product, bounded entrywise by R|V̂|; ‖R|V̂|‖₂ discards the cancellation that makes ‖RV̂‖₂ ≤ ‖R‖₂ and costs a factor growing like √n. That is a reason to keep R out of the theorem, not a reason to use Theorem 7: the Weyl route runs Theorem 10 on mid(A) and widens by ‖R‖₂ once. Maximum radii, randn matrices at MersenneTwister(2026+n):

A graded matrix with σ₁/σ_q = 10⁸ and radius 1e-14 is the one case measured where the Weyl route is the wider of the two, 2.076e-13 against Theorem 7's 2.068e-13, four parts in a thousand; the absolute widening ‖R‖₂ is felt at the small singular values, where Theorem 7's ‖Ê‖ is felt equally. No threshold in rad(A) is fitted, since there is no longer a factor to avoid on either side of one.

The test is iszero(rad(A)), which is maximum(rad(A)) == 0 since radii are non-negative; it is spelled with iszero only because that stops at the first nonzero instead of scanning the whole matrix. A single radius of 1e-30 therefore sends the matrix through the Weyl route, which adds ‖R‖₂ = 1e-30 and changes nothing. Pass method explicitly to override the rule.

Miyajima2014aSchurNewtonResult

Container returned by miyajima2014a_schurnewton. Field names are duck-type compatible with SchurGershgorinResult so the same downstream consumers (block Schur, spectral projectors) work unchanged, with three extra certification scalars (nrmR2, beta, kappa).

The basis W is block-orthonormal, not globally unitary, so projector / similarity consumers must use inv(W) (not adjoint(W)); see _vbd_unitary_basis.

miyajima2014a_schurnewton(A::BallMatrix; sep = -1, maxsteps = 6, kappa_mode = :cheap)

Verified block diagonalisation of the square ball matrix A via the O(n³) Schur + Newton route, with a rigorous Miyajima two-residual certification.

The midpoint is reduced to a block-orthonormal basis W (Schur frame, distance clustering, between-cluster Newton decoupling, per-block QR). The enclosure is transported to that basis and the β-inflated Gershgorin discs of à = inv(W)·A·W rigorously enclose σ(M) for every M ∈ A. The discs are re-clustered by overlap (the rigorous arbiter — an over-optimistic Newton split whose discs still overlap is re-merged) so that each cluster is a contiguous range holding exactly its eigenvalue count.

Unlike the removed NJD route this is rigorous (it charges the Y ≠ W⁻¹ error through β) and refuses (throws) when the basis fails the Neumann condition ‖R₂‖_∞ < 1, i.e. when mid(A) is numerically defective at the working precision. kappa_mode selects the κ₂ diagnostic: :cheap (Collatz, keeps O(n³)) or :svdbox (one verified SVD of W).

miyajima2014a_schurnewton(A::BallMatrix, B::BallMatrix; sep = -1, maxsteps = 6,
                 kappa_mode = :cheap)

Verified block diagonalisation of the pencil Ax = λBx, by the same Schur + Newton route as the one-argument method.

The change of variable x = Wy turns the pencil into (YAW)y = λ(YBW)y, so with an approximate inverse Y ≈ (BW)⁻¹ and the residual R₂ = YBW − I the collapsed matrix is again M = (I+R₂)⁻¹Ã with à = YAW. That is structurally identical to the standard problem, so the per-row certification slack, the β-inflated Gershgorin discs, the overlap re-clustering and the block enclosure carry over unchanged — the only substitution is W ⟶ B·W wherever the frame is paired with Y or with the candidate Λ (R₁ = Y(AW − BWΛ)).

B is not assumed Hermitian or positive definite; its nonsingularity is not assumed either but proved, by ‖R₂‖_∞ < 1 (which certifies B, W and Y nonsingular at once). This follows Miyajima (2014), Fast enclosure for all eigenvalues and invariant subspaces in generalized eigenvalue problems, SIAM J. Matrix Anal. Appl. 35(3), 1205–1225, whose block-diagonalisation algorithm (§4) is obtained here through the Schur–Newton frame rather than the Kronecker fixed point of that paper.

The candidate basis is produced in floating point from mid(B) \ mid(A); being a candidate it needs no verification of its own, all correctness resting on the final certification.

ordschur_bigfloat(T::AbstractMatrix, Q::AbstractMatrix,
                  select::AbstractVector{Bool})

Reorder a Schur decomposition A = Q T Q^H so that the eigenvalues corresponding to select[i] == true appear in the top-left block.

The reordering is ordschur of GenericSchur, applied after the arguments are wrapped in a Schur object, and it runs in the arithmetic of T and Q, which may be BigFloat or Complex{BigFloat} where LAPACK offers only Float64. The eigenvalues are taken as diag(T), so T is expected triangular.

It is floating point throughout and carries no bound: the rotations are not accounted anywhere, so the returned pair is an approximate Schur pair of whatever Q T Qᴴ the arguments came from, and nothing is claimed about any matrix. ordschur_ball is the routine that reorders and then measures what the reordered pair satisfies.

Arguments

• T: Upper triangular Schur form (n × n)
• Q: Unitary Schur basis (n × n)
• select: Boolean mask of length n — true for eigenvalues to move to top-left

Returns

(Q_ord, T_ord, values) — reordered Schur basis, Schur form, and eigenvalues.

ordschur_ball(Q_ball::BallMatrix, T_ball::BallMatrix, select; A = nothing)

Reorder an approximate Schur pair so that the diagonal entries with select[i] == true come first, and bound the defects of the reordered pair.

The reordering itself is the floating-point one: ordschur on the midpoint of T_ball returns the reordered matrix and the accumulated rotation G. No error is propagated through the rotations. What the reordered pair satisfies is measured afterwards, in ball arithmetic.

What is returned

A NamedTuple with

• G: the accumulated rotation, a matrix with exact entries (it is only approximately unitary);
• T: T̃, the upper triangle of the reordered midpoint, as a BallMatrix of radius zero. It is exactly upper triangular, and whatever the floating-point reordering left below the diagonal is in the residuals below;
• Q: the ball matrix Q_ball·G;
• values: the diagonal of T̃;
• rotation_orth_defect: upper bound of ‖G*G − I‖₂;
• rotation_residual: upper bound of ‖TG − GT̃‖₂ for every T in T_ball;
• orth_defect, fact_defect: when the matrix A is given (a BallMatrix), upper bounds of ‖Q*Q − I‖₂ and of ‖AQ − QT̃‖₂ for every Q in the returned ball and every A in the ball A; nothing otherwise.

These are the two numbers spectral_projector_error_bound takes. Without A, they follow from the defects δ₀ ≥ ‖Q₀*Q₀ − I‖₂ and ε₀ ≥ ‖AQ₀ − Q₀T‖₂ of the pair that was given, since AQ₀G − Q₀GT̃ = (AQ₀ − Q₀T)G + Q₀(TG − GT̃):

‖AQ − QT̃‖₂ ≤ ε₀ √(1 + δ_G) + √(1 + δ₀) ρ,

with δ_G and ρ the two rotation bounds.

T_ball is expected to be upper triangular (a complex Schur form). The cost is that of ordschur plus three products of n×n ball matrices, and three more when A is given.

spectral_projector_error_bound(; resolvent_bound_A, contour_radius, orth_defect, fact_defect)

An upper bound of ‖P_A − P_c‖₂, where P_A is the spectral projector of A for the part of its spectrum inside a contour Γ, and P_c = Q P_T̃ Q* with P_T̃ the spectral projector of a matrix T̃ for the part of its spectrum inside the same contour. Q and T̃ are any pair with δ = orth_defect ≥ ‖I − Q*Q‖₂ < 1 and ε = fact_defect ≥ ‖AQ − QT̃‖₂, for instance the one returned by ordschur_ball with A given.

The bound

Both projectors are contour integrals, P_A = (1/2πi)∮_Γ (zI − A)⁻¹ dz and P_c = (1/2πi)∮_Γ Q(zI − T̃)⁻¹Q* dz. With R = AQ − QT̃, so that (zI − A)Q = Q(zI − T̃) − R,

(zI − A)⁻¹ − Q(zI − T̃)⁻¹Q* = (zI − A)⁻¹ [ (I − QQ*) + R (zI − T̃)⁻¹ Q* ].

With M_A ≥ ‖(zI − A)⁻¹‖₂ on Γ (resolvent_bound_A), M_T̃ ≥ ‖(zI − T̃)⁻¹‖₂ on Γ, ‖I − QQ*‖₂ = ‖I − Q*Q‖₂ ≤ δ for a square Q, ‖Q‖₂ ≤ √(1 + δ), and a contour of length at most 2πr (contour_radius; a circle of radius r or a polygon inscribed in it),

‖P_A − P_c‖₂ ≤ r · M_A · ( δ + ε · M_T̃ · √(1 + δ) ).

M_T̃ is obtained from M_A: with E = I − Q*Q and F = Q*(zI − A)Q, (I − E)(zI − T̃) = F + Q*R, ‖F⁻¹‖₂ ≤ M_A/(1 − δ), so with γ = M_A √(1 + δ) ε/(1 − δ) < 1,

M_T̃ ≤ M_A (1 + δ) / ((1 − δ)(1 − γ)).

The hypothesis that T̃ has no spectrum on Γ follows from γ < 1. Every operation is rounded up. Returns Inf when δ ≥ 1 or γ ≥ 1.

Reference

The comparison of two Riesz projections through the resolvent on a contour is the lemma "Resolvent bridge" of I. Nisoli, Certified spectral approximation of transfer operators and the Gauss map, arXiv:2602.19435 (2026), in its section on spectral stability and projections; the form above, with a pair Q, T̃ that is neither exactly unitary nor an exact factorisation, is derived in this docstring.

SchurGershgorinResult

Container returned by schur_gershgorin_enclosure. The fields hold the unitary basis (Schur, or eigenvectors when the input is Hermitian), the enclosure of A transported to it, the block-diagonal part of that enclosure, the rigorous remainder, and the clusters read off the inflated Gershgorin discs.

The basis is not a block-diagonalising one: it is the unitary basis of the midpoint matrix, permuted so the clusters are contiguous, so remainder is the off-block part of Z*AZ and is not small. Field names are duck-type compatible with Miyajima2014aSchurNewtonResult, whose basis does decouple the blocks, so the same downstream consumers accept either.

schur_gershgorin_enclosure(A::BallMatrix)
schur_gershgorin_enclosure(H::Hermitian)

Enclose σ(M) for every M in the square ball matrix A by inflated Gershgorin discs in a unitary basis of mid(A), grouped into blocks.

The type of the argument decides the kind of enclosure. For a BallMatrix the clusters are discs of the complex plane: with positive radii it is an open set of matrices, which contains matrices that are not Hermitian, whatever its midpoint. For a HermitianBallMatrix, the Hermitian matrices within the tolerances, and for a point matrix that is Hermitian by its type (Hermitian(M), or Symmetric(M) with M real), the eigenvalues are real and the clusters are real intervals; the statement is then for those Hermitian matrices only. The same holds for a real ball matrix stored as Diagonal in midpoint and radius, the type of the matrix of singular value enclosures of rigorous_svd: every matrix in it is a real diagonal matrix, whose eigenvalues are its diagonal entries.

The midpoint matrix is reduced by a unitary Schur form, or by an eigenvalue decomposition when the midpoint is Hermitian, which only chooses the basis and does not change the kind of enclosure; the enclosure is transported to that basis as à = Z*AZ, a ball product against the input ball so input radii propagate. Z is unitary only up to rounding, so the residual R₂ = Z*Z − I is formed and ‖R₂‖_∞ < 1 required, and the per-row slack βᵢ = ‖row_i(R₂)‖₁‖Ã‖_∞/(1−‖R₂‖_∞) is added to disc i before anything is read off it. The inflated discs are then clustered by their connectivity graph and the basis permuted so each cluster is contiguous; the block-diagonal part of Ã and its rigorous remainder follow, with remainder_norm a rigorous bound on ‖Ã − block‖₂ (the best of the Collatz and ‖·‖₁‖·‖∞ estimates).

Each cluster Cᵢ also carries a disc D(cᵢ, nᵢ + rᵢ), with barycentre cᵢ = mean(diag Pᵢ), within-block non-normality nᵢ = ‖Pᵢ − cᵢI‖₂, and coupling rᵢ = min(Σ_{j≠i}‖Ñ[Cᵢ,Cⱼ]‖₂, √(q−1)‖Ñ[Cᵢ,:]‖₂) + β₂. Every eigenvalue lies in one of these discs: by Theorem 2 of the reference every eigenvalue z has σ_min(Pᵢ − zI) ≤ rᵢ for some i, its proof needing only a bound of ‖Σ_{j≠i} Ñ[Cᵢ,Cⱼ] w_j‖₂ for ‖w_j‖₂ ≤ 1, which both terms of the minimum are (the second is not in the reference: the q − 1 blocks w_j have ‖(w_j)_j‖₂ ≤ √(q−1)); and σ_min(Pᵢ − zI) ≥ |z − cᵢ| − nᵢ. Theorem 4 of the reference counts the eigenvalues in a union of its sets disjoint from the others, by a continuity argument written for the radii Σ_{j≠i}‖A_ij‖; for the radii rᵢ above, with the strip form and β₂, that argument has not been written out here.

This is not either algorithm of Miyajima (2014), whose name the routine formerly carried. That paper certifies clusters by a Brouwer fixed-point test on a Newton operator, over an approximate eigendecomposition (§3, radii u + ⟨u⟩_t t from infinity-norm row sums) or a genuine block diagonalisation of B⁻¹A (§4); neither the fixed-point test nor a diagonalising transformation appears here, and the basis is only permuted. For the block-diagonalising route see miyajima2014a_schurnewton.

Reference

D. G. Feingold and R. S. Varga, Block diagonally dominant matrices and generalizations of the Gerschgorin circle theorem, Pacific J. Math. 12(4) (1962), 1241-1250, doi 10.2140/pjm.1962.12.1241, Theorems 2 and 4.

OgitaSVDRefinementResult

Result from Ogita's iterative SVD refinement.

AdaptiveSVDResult

Result from adaptive precision SVD computation.

ogita_svd_refine(A::AbstractMatrix{T}, U, Σ, V;
                 max_iterations=10, precision_bits=256,
                 check_convergence=false,
                 use_optimal_iterations=false) where {T<:AbstractFloat}

Refine an approximate SVD using Ogita's iterative method (RefSVD algorithm).

Arguments

• A: Original matrix (in higher precision if needed)
• U: Initial left singular vectors
• Σ: Initial singular values
• V: Initial right singular vectors
• max_iterations: Maximum number of iterations to run (default: 10)
• precision_bits: Working precision in bits (default: 256)
• check_convergence: If true, check convergence via spectral norm and stop early. If false (default), run fixed number of iterations.
• use_optimal_iterations: If true, compute optimal iteration count from quadratic convergence theory and ignore max_iterations. This is the fastest option for known target precision.

Iteration count (quadratic convergence from Float64)

Based on theory, starting from Float64 (~53 bits), each iteration doubles precision:

• 3 iterations: sufficient for 256-bit (~77 decimal digits)
• 4 iterations: sufficient for 512-bit (~154 decimal digits)
• 5 iterations: sufficient for 1024-bit (~308 decimal digits)

Returns

• OgitaSVDRefinementResult containing refined SVD

References

• [6] Ogita, T. & Aishima, K. (2020), "Iterative refinement for singular value decomposition based on matrix multiplication", J. Comput. Appl. Math. 369, 112512.

adaptive_ogita_svd(A::BallMatrix{T};
                   tolerance=1e-10,
                   method::SVDMethod=MiyajimaM1(),
                   apply_vbd=true,
                   max_precision_bits=1024,
                   max_refinement_iterations=5) where {T}

Compute rigorous SVD bounds with adaptive precision using Ogita's refinement.

Algorithm

1. Start in the arithmetic of A: Float64, or BigFloat at the precision of the midpoints of A. The result for a BigFloat matrix is in BigFloat; for a Float64 matrix it is in Float64 when the tolerance is met at the first level and in BigFloat otherwise. Where a level works at a lower precision than an entry of A has, the entry is rounded and the rounding error goes into the radii, so that what is certified holds for A.
2. Compute rigorous SVD bounds using rigorous_svd
3. Check if max(rad(σᵢ)) < tolerance
4. If not satisfied, refine using Ogita's method with doubled precision
5. Repeat until tolerance met or max precision reached

A matrix with non-zero radii

The enclosures must contain the singular values of every matrix of the ball, so their radii cannot be brought below what the radii of A impose, at any precision. The singular values are Lipschitz with constant 1 in the 2-norm, |σᵢ(X) − σᵢ(M)| ≤ ‖X − M‖₂, and ‖X − M‖₂ ≤ ‖R‖₂ when |X − M| ≤ R entrywise; both are proved in the docstring of _miyajima2014_thm10_weyl, which uses the same step. The precision is therefore adapted on the midpoint of A, and tolerance, radii_history and tolerance_achieved are about the enclosures of the singular values of the midpoint. In rigorous_result each of those enclosures is widened by an upper bound of ‖rad(A)‖₂, so that it holds for every matrix of the ball; the residual there is U Σ Vᴴ − A with that Σ.

Arguments

• A: Input ball matrix
• tolerance: Target tolerance for max(rad(σᵢ)), of the midpoint when A has radii
• method: SVD certification method (MiyajimaM1 for Theorem 7, MiyajimaM4 for Theorem 11)
• apply_vbd: Whether to apply verified block diagonalization
• max_precision_bits: Maximum precision to use (bits)
• max_refinement_iterations: Maximum number of refinement steps

Returns

• AdaptiveSVDResult containing final rigorous result and adaptation history

Example

A = BallMatrix([3.0 1.0; 1.0 2.0], fill(1e-8, 2, 2))
result = adaptive_ogita_svd(A; tolerance=1e-12)

# Access final result
σ = result.rigorous_result.singular_values
println("Final precision: ", result.final_precision, " bits")
println("Max radius: ", maximum(rad.(σ)))

clear_svd_cache!()

Clear the BigFloat SVD cache used for warm-starting Ogita refinement.

svd_cache_stats()

Return statistics about the BigFloat SVD cache usage.

set_svd_cache!(U, S, V, A_hash)

Set the SVD cache with the given factors and matrix hash. Used for warm-starting Ogita refinement on nearby matrices.

RigorousBlockSchurResult

Returned by rigorous_block_schur. For the matrix A, a basis Q, a matrix Q_inv ≈ Q⁻¹ computed in floating point, and a block matrix T, it holds the two residuals

R₁ = Q_inv (A Q − Q T),        R₂ = Q_inv Q − I,

through upper bounds of their spectral norms, projected_residual_norm and inverse_defect_norm. When inverse_defect_norm < 1, Q is invertible, Q⁻¹ = (I + R₂)⁻¹ Q_inv, and

Q⁻¹ A Q = T + S,        S = (I + R₂)⁻¹ R₁,        ‖S‖₂ ≤ ‖R₁‖₂ / (1 − ‖R₂‖₂),

an identity; the last bound is perturbation_norm (Inf when inverse_defect_norm ≥ 1). So A is similar to a matrix within perturbation_norm of T, and no inverse is enclosed. The statements hold for every A, Q and T in the respective balls.

S itself is enclosed through a floating-point solution S̃ of (I + R₂)S̃ = R₁ (perturbation_approx) and the residual V = (I + R₂)S̃ − R₁, computed in ball arithmetic:

S − S̃ = −(I + R₂)⁻¹V,        ‖S − S̃‖_p ≤ ‖V‖_p / (1 − ‖R₂‖_p),        p ∈ {2, ∞},

and an entry of a matrix is at most its norm in either, so every entry of S − S̃ is at most perturbation_error, the smaller of the two bounds (Inf when neither ‖R₂‖_p is below 1). similar = T + S̃ ± perturbation_error is then a ball matrix containing Q⁻¹AQ.

Q_inv is the adjoint of Q when the basis is unitary and inv of its midpoint otherwise; it is NOT Q' in general.

rigorous_block_schur(A::BallMatrix; block_structure = :quasi_triangular, vbd_method = :auto)
rigorous_block_schur(H::Hermitian; block_structure = :quasi_triangular, vbd_method = :auto)

A block form T of A in the basis Q of a verified block diagonalisation, with the residuals that bound its distance from a similarity. See RigorousBlockSchurResult for the statement that the result carries.

Steps

1. The verified block diagonalisation gives the clusters, their intervals and the basis: miyajima2014a_schurnewton (Schur form followed by a Newton step for the decoupling) for a ball matrix, schur_gershgorin_enclosure for a HermitianBallMatrix or a point matrix that is Hermitian by its type (vbd_method = :schur_newton or :nsd; :auto takes :nsd when the midpoint is Hermitian, which is a choice of basis and not a statement about the ball). The clusters and intervals of the result are the certified ones of that routine: discs of the complex plane for a ball matrix, which with positive radii contains matrices that are not Hermitian, and real intervals for a Hermitian ball matrix, where the statement is for the Hermitian matrices within the tolerances.
2. Q is that basis and Q_inv an approximate inverse in floating point (the adjoint for a unitary basis).
3. T = Q_inv A Q in ball arithmetic, with the blocks that block_structure drops set to zero: :diagonal keeps the diagonal blocks, :quasi_triangular the blocks on and above the diagonal, :full everything.
4. R₁ = Q_inv (A Q − Q T) and R₂ = Q_inv Q − I are computed in ball arithmetic and their spectral norms bounded from above; the dropped blocks are therefore in R₁.
5. (I + R₂)S̃ = R₁ is solved in floating point and the residual of that solve bounded in ball arithmetic, which gives the ball matrix similar containing Q⁻¹AQ. The dropped blocks are in S̃, so similar does not depend on block_structure beyond rounding.

The device of an arbitrary Y ≈ X⁻¹ with the two residuals Y(AX − XD) and YX − I in place of an enclosed inverse is that of Theorem 2 of

S. Miyajima, Fast enclosure for all eigenvalues and invariant subspaces in generalized eigenvalue problems, SIAM J. Matrix Anal. Appl. 35 (2014), 1205–1225, doi:10.1137/140953150 ([25]),

used here with a block matrix T in place of the diagonal one. The enclosure of S by a computed S̃ and the residual of its solve is the lemma "Residual control of S" of I. Nisoli, Enclosure of eigenvalues via approximate diagonalization (unpublished note, 2026); the argument is the identity S − S̃ = −(I + R₂)⁻¹((I + R₂)S̃ − R₁) stated in RigorousBlockSchurResult.

extract_cluster_block(result::RigorousBlockSchurResult, i::Int, j::Int)

Extract the (i,j)-th block from the block Schur form T. Returns a BallMatrix corresponding to T[cluster_i, cluster_j].

verify_block_schur_properties(result::RigorousBlockSchurResult; tol = 1e-10)

Whether projected_residual_norm and inverse_defect_norm are both below tol.

estimate_block_separation(result::RigorousBlockSchurResult, i::Int, j::Int)

Estimate the spectral separation between clusters i and j. A small separation indicates potential numerical difficulties in separating the corresponding invariant subspaces.

SchurSpectralProjectorResult{T, NT}

Returned by compute_spectral_projector_schur and compute_spectral_projector_hermitian.

Fields

• projector: when projector_error_bound is finite, an enclosure of the spectral projector P_A of every matrix of the ball A; otherwise an enclosure of the computed P_c = Q P_T̃ Q* only, which is not P_A
• schur_projector: enclosure of P_T̃ = [I Y; 0 0], the spectral projector of T̃ for its first k diagonal entries
• coupling_matrix: enclosure of Y, the solution of T₁₁Y − YT₂₂ = T₁₂
• eigenvalue_separation: min |t_ii − t_jj| over i ≤ k < j, in floating point (a diagnostic)
• projector_norm: upper bound of ‖P‖₂ over the returned ball
• idempotency_defect: upper bound of ‖P² − P‖₂ over the returned ball (zero when not asked)
• schur_basis, schur_form: the pair Q, T̃ used, the latter exactly upper triangular
• cluster_indices: 1:k, the positions of the selected eigenvalues in T̃
• orth_defect: upper bound of ‖Q*Q − I‖₂
• fact_defect: upper bound of ‖AQ − QT̃‖₂ over the ball A
• projector_error_bound: upper bound of ‖P_A − P_c‖₂ from spectral_projector_error_bound; Inf when no resolvent bound was given

compute_spectral_projector_schur(A::BallMatrix, cluster_indices;
    resolvent_bound = nothing, contour_radius = nothing, verify_idempotency = true,
    schur_data = nothing, ordered_schur_data = nothing, sylvester_fallback = :direct)

The spectral projector of A for the eigenvalues at the positions cluster_indices (a range or a vector of distinct indices) of the diagonal of its Schur form.

What is computed, and what is proved

1. In floating point: a Schur pair of the midpoint of A (or schur_data = (Q, T)), reordered so that the selected eigenvalues come first (or ordered_schur_data = (Q, T, k), already in that order). T̃ is its upper triangle.
2. In ball arithmetic: the defects δ ≥ ‖Q*Q − I‖₂ and ε ≥ ‖AQ − QT̃‖₂, the second over the whole ball A.
3. P_T̃ = [I Y; 0 0], the spectral projector of the triangular T̃ for its first k diagonal entries, with Y the solution of T₁₁Y − YT₂₂ = T₁₂, enclosed by triangular_sylvester_miyajima_enclosure; and P_c = Q P_T̃ Q* in ball arithmetic.
4. With resolvent_bound ≥ ‖(zI − A)⁻¹‖₂ on a contour of length at most 2π·contour_radius that has the selected eigenvalues of A inside and the others outside (for instance from CertifScripts.run_certification), spectral_projector_error_bound gives b ≥ ‖P_A − P_c‖₂, and the radius of every entry of the returned projector is increased by b: the result is then an enclosure of P_A for every matrix of the ball A.

Without the two keywords step 4 is not done, projector_error_bound is Inf, and the returned ball encloses P_c only. P_c is the projector of a nearby matrix; its distance from P_A is exactly what step 4 bounds, and ball arithmetic on Q and T̃ alone says nothing about it.

That the contour separates the spectrum as required is the caller's hypothesis; the bound γ < 1 inside spectral_projector_error_bound shows that T̃ has no eigenvalue on it.

verify_idempotency adds the bound of ‖P² − P‖₂ over the returned ball.

compute_spectral_projector_hermitian(A::BallMatrix, cluster_indices::UnitRange{Int};
    resolvent_bound = nothing, contour_radius = nothing)

compute_spectral_projector_schur with the eigenvector matrix of the Hermitian midpoint of A as Q and the diagonal matrix of its eigenvalues (in increasing order) as T̃; cluster_indices are positions in that order. The coupling Y is zero, and the statement and the role of the two keywords are those of the general routine.

project_vector_spectral(v, result::SchurSpectralProjectorResult)

The product of the returned projector ball with v (a BallVector, or a vector with exact entries), in ball arithmetic. It encloses P_A v when result.projector_error_bound is finite, and P_c v otherwise.

verify_spectral_projector_properties(result, A::BallMatrix; tol = 1e-10,
                                     check_commutation = false)

Whether result.idempotency_defect < tol, result.projector_norm is finite, result.eigenvalue_separation is positive and, when check_commutation is set, the upper bound of ‖AP − PA‖₂ over the two balls is below tol. These are consistency checks of the returned ball; they do not replace projector_error_bound.

SpectralCoefficientResult{T, NT}

Returned by compute_spectral_coefficient.

Fields

• coefficients: when coefficient_error_bound is finite, an enclosure of Q₁* P_A v, with Q₁ the first k columns of schur_basis, for every matrix of the ball A; otherwise an enclosure of the computed [I Y] Q* v only
• left_eigenvector_schur: enclosure of the k×n block [I Y]
• coupling_matrix: enclosure of Y
• eigenvalue_separation: min |t_ii − t_jj| over i ≤ k < j, in floating point
• schur_basis, schur_form: the pair Q, T̃ used
• orth_defect, fact_defect: upper bounds of ‖Q*Q − I‖₂ and ‖AQ − QT̃‖₂
• coefficient_error_bound: what was added to the radius of every coefficient; Inf when no resolvent bound was given

compute_spectral_coefficient(A::BallMatrix, v, cluster_indices;
    resolvent_bound = nothing, contour_radius = nothing,
    schur_data = nothing, ordered_schur_data = nothing, sylvester_fallback = :direct)

The coordinates of the spectral projection of v in the first k Schur vectors, without forming the n×n projector: with the ordered pair (Q, T̃) and the coupling Y of compute_spectral_projector_schur,

c = [I Y] Q* v,        so that        P_c v = Q₁ c,

Q₁ being the first k columns of Q. v is a BallVector or a vector with exact entries.

With resolvent_bound and contour_radius (as in compute_spectral_projector_schur), b ≥ ‖P_A − P_c‖₂ and

‖Q₁* P_A v − c‖₂ ≤ ‖Q₁‖₂ ‖(P_A − P_c) v‖₂ + ‖(Q₁*Q₁ − I) c‖₂ ≤ √(1 + δ) · b · ‖v‖₂ + δ ‖c‖₂,

which is added to the radius of every coefficient: the result then encloses Q₁* P_A v. Without them coefficient_error_bound is Inf and the ball encloses the computed c only.

SchurRefinementResult{T, RT}

Result from iterative Schur refinement to higher precision.

Fields

• Q::Matrix{T}: Refined orthogonal/unitary Schur basis
• T::Matrix{T}: Refined upper triangular Schur form
• iterations::Int: Number of refinement iterations performed
• residual_norm::RT: Final residual ‖stril(Q^H A Q)‖₂ / ‖A‖_F — rigorous upper bound (real type)
• orthogonality_defect::RT: Final ‖Q^H Q - I‖₂ — rigorous upper bound (real type)
• converged::Bool: Whether refinement converged to desired tolerance

References

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

refine_schur_decomposition(A, Q0, T0; method = :bujanovicetal2023,
    target_precision = 256, max_iterations = 20, tol = 0, threshold = 1e-5, rng = nothing)

Refines an approximate Schur decomposition A ≈ Q0 T0 Q0ᴴ computed in Float64, T0 upper triangular, to BigFloat at target_precision bits. Returns a SchurRefinementResult with the factors in BigFloat; its residual_norm and orthogonality_defect are floating-point diagnostics, not bounds.

Each algorithm is an unexported function named for the paper it implements, a deviation being visible in the name:

The default is the algorithm as published. It and the truncated one reorder the eigenvalues along a random line (rng gives the generator), so the diagonal of the result is in another order than that of T0, and they need Q0, T0 in a LAPACK type. The callers in this package that order the Schur form before refining it and rely on that order pass method = :bujanovicetal2023unordered.

tol is the tolerance of the stopping test, by default 100 rounding units of the target precision. A T0 that is not upper triangular is an ArgumentError: for a real matrix with complex eigenvalues pass the complex Schur form, schur(complex(A)); the paper mentions an extension to the real Schur form and does not give it.

Example

A = randn(100, 100)
F = schur(complex(A))
r = refine_schur_decomposition(A, F.Z, F.T; method = :bujanovicetal2023)
r.residual_norm          # near 1e-77 at 256 bits

Reference

Z. Bujanović, D. Kressner and C. Schröder, "Iterative refinement of Schur decompositions", Numerical Algorithms 92 (2023), 247-267, doi 10.1007/s11075-022-01327-6.

rigorous_schur_bigfloat(A::BallMatrix{T}; target_precision::Int=256,
                        max_iterations::Int=20,
                        schur_seed=nothing) where T

An approximate Schur pair of the midpoint of A at BigFloat precision, as ball matrices of radius zero:

1. approximate Schur decomposition in Float64 (or schur_seed);
2. iterative refinement to BigFloat precision (refine_schur_decomposition);
3. the pair wrapped by certify_schur_decomposition.

The name is historical: the pair is not an enclosure of a Schur decomposition of A, and the rigorous statements about it are the defects measured by ordschur_ball and compute_spectral_projector_schur. See certify_schur_decomposition.

Arguments

• A::BallMatrix: Input ball matrix
• target_precision: Target BigFloat precision in bits (default: 256)
• max_iterations: Maximum refinement iterations (default: 20)
• schur_seed: Optional (Q0, T0) tuple — external Schur seed to use instead of computing a Float64 Schur internally. Useful when an accurate decomposition is already available (e.g. from GenericSchur.schur in BigFloat).

Returns

A tuple (Q_ball, T_ball, result) where:

• Q_ball::BallMatrix{BigFloat}: the refined, approximately unitary factor, radius zero
• T_ball::BallMatrix{BigFloat}: the refined upper triangular factor, radius zero
• result::SchurRefinementResult: floating-point diagnostics of the refinement

Example

A = BallMatrix(randn(10, 10), fill(1e-10, 10, 10))
Q, T, result = rigorous_schur_bigfloat(A; target_precision=512)
@show result.converged

# With external seed:
F = schur(Complex{BigFloat}.(mid(A)))
Q2, T2, r2 = rigorous_schur_bigfloat(A; schur_seed=(F.Z, F.T))

References

• [7] for the iterative refinement algorithm

certify_schur_decomposition(A::BallMatrix, result::SchurRefinementResult) -> (Q_ball, T_ball)

The refined pair Q, T of a SchurRefinementResult as ball matrices of radius zero, T with its part below the diagonal set to zero.

They are not enclosures of a Schur decomposition of A. Q is an approximately unitary matrix and T an upper triangular one with A Q ≈ Q T; what can be proved about them is the size of the two defects ‖Q*Q − I‖₂ and ‖AQ − QT‖₂ over the ball A, and those are measured in ball arithmetic by ordschur_ball (with A given) and by compute_spectral_projector_schur, which take this pair. The fields residual_norm and orthogonality_defect of result are floating-point diagnostics of the refinement.

Until 2026-10-11 the radii were the two diagnostics plus the norm of rad(A), added uniformly to every entry, and the result was described as an enclosure of the Schur factors; the audit found an eigenvalue at distance 1.65e-71 from a diagonal entry given the radius 6.4e-77. For inclusions of the factors of a Schur decomposition, see _rumpogita2024_schur (Rump and Ogita (2024), Section 8, for simple eigenvalues).

newton_schulz_orthogonalize!(Q::Matrix{T}; max_iter=10, tol=nothing) where T

Apply Newton-Schulz iteration to orthogonalize Q in-place.

The iteration is: Q{k+1} = (1/2) Qk (3I - Qk^H Qk)

It converges quadratically to the unitary factor of the polar decomposition of Q when Q is invertible and ‖Q‖₂ < √3 (Section 2.2.2 of the reference); ‖Q^H Q - I‖₂ < 1 implies both.

Arguments

• Q: Matrix to orthogonalize (modified in place)
• max_iter: Maximum iterations (default: 10)
• tol: Convergence tolerance (default: machine epsilon of T)

Returns

• Number of iterations performed
• Final orthogonality defect ‖Q^H Q - I‖_F

The iteration and the defect are computed in floating point; the result is an approximation and nothing about it is certified.

Reference

Z. Bujanović, D. Kressner and C. Schröder, Iterative refinement of Schur decompositions, Numer. Algorithms 92(1) (2023), 247-267, doi 10.1007/s11075-022-01327-6, Section 2.2.2.

solve_triangular_matrix_equation(T_mat::Matrix{T}, E::Matrix{T};
                                 use_recursive::Bool=true, nmin::Int=32, max_entry::Real=1e5) where T

An approximate solution, in floating point, of the triangular matrix equation stril(TL - LT) = -E for a strictly lower triangular L, by Algorithm 2 (successive substitution in bottom-to-top columnwise order) or Algorithm 3 (recursive, by 2 × 2 blocks) of [7]. For an upper triangular T with pairwise distinct diagonal entries and a strictly lower triangular E the equation has a unique solution; every entry of the computed L whose modulus exceeds max_entry is set to zero when it is formed, so the result need not satisfy the equation. No bound comes with it.

Arguments

• T_mat: Upper triangular matrix (n × n) with pairwise distinct diagonal entries
• E: Strictly lower triangular matrix (n × n)
• use_recursive: If true, use recursive block algorithm (default: true)
• nmin: Minimum block size for recursive algorithm (default: 32)
• max_entry: entries of L of larger modulus are set to zero (default: 1e5). The paper proposes this in the discussion of its Example 4, with the threshold written there as 1e-5, to keep large entries from propagating when eigenvalues are clustered; the default here is ten orders looser.

Returns

The strictly lower triangular matrix L.

Reference

[7], Section 2.1, equation (6) for the equation and Algorithms 2 and 3 for the two solvers.

_bujanovicetal2023_alg4(A, Q0, T0, max_iterations, tol; hp = BigFloat, rng = nothing)

Algorithm 4 of Bujanović, Kressner and Schröder (2023), the mixed-precision refinement of a Schur decomposition, from line 1's ordering on: Q0, T0 are the Schur factors of A in low precision (Float64), T0 upper triangular.

1. The Schur form is reordered so that eigenvalues close to each other are neighbours on the diagonal: by the order of their projections on a random line through the origin.
2. Q̂ is converted to high precision and Q ← ½ Q̂ (3I − Q̂ᴴQ̂), one Newton-Schulz step.
3. Until E is small: T̂ ← QᴴAQ in high precision; E ← stril(T̂), T ← T̂ − E; the equation stril(T L − L T) = −E is solved for the strictly lower triangular L in low precision, by the recursive Algorithm 3; W = L − Lᴴ, Y = QᴴQ − I; the products YW, W², W³ in low precision; Σ = 2I + 2W − Y − YW + W² + W³ and Q ← ½ Q Σ in high precision.

hp is the high precision: BigFloat at the current precision, or a type such as Double64. The iteration stops when ‖E‖/‖A‖ ≤ tol. The result is in floating point: its two norms are diagnostics, not bounds. The paper reports that the iteration can fail when eigenvalues are poorly separated (its Example 3), the entries of L growing until they are not numbers.

Not the paper's: nmin, the size at which Algorithm 3 stops recursing, is 32 here; the paper's experiments use 4 for the complex Schur form.

_bujanovicetal2023_alg4_truncated(A, Q0, T0, max_iterations, tol;
    hp = BigFloat, rng = nothing, threshold = 1e-5)

_bujanovicetal2023_alg4 with the modification of Example 4 of the paper: "once an element of L is computed (Line 4 of Algorithm 2, Line 6 of Algorithm 3), if its absolute value turns out to be larger than, e.g., 10⁻⁵, such element is immediately set to zero. The idea is to prevent propagation of large elements in L, at the cost of obtaining a less accurate solution of the triangular equation." The paper adds that "there is little theoretical support that such a modification will lead to convergence"; it resolved the failures of its Table 3 and not the one of its Example 3.

_bujanovicetal2023_alg4_unordered(A, Q0, T0, max_iterations, tol; hp = BigFloat)

The iteration of _bujanovicetal2023_alg4 without the reordering of its line 1, so that the eigenvalues stay on the diagonal in the order of T0: a deviation, named for it. It is for the callers that have ordered the Schur form themselves and use that order afterwards. Three further differences from the paper's listing: Q0, T0 are returned as they are when both defects are already below the tolerance; the initial orthogonalisation takes up to three Newton-Schulz steps; and the iteration stops when both ‖E‖/‖A‖ and ‖QᴴQ − I‖ are below the tolerance. Q0 and T0 may be in any precision.

SymmetricEigenRefinementResult{T, RT}

Result from RefSyEv iterative refinement of symmetric eigenvalue decomposition.

Fields

• Q::Matrix{T}: Refined orthogonal eigenvector matrix
• λ::Vector{RT}: Refined eigenvalues (real)
• iterations::Int: Number of refinement iterations performed
• residual_norm::RT: Final ‖A - QΛQ^T‖F / ‖A‖F
• orthogonality_defect::RT: Final ‖Q^T Q - I‖_F
• converged::Bool: Whether refinement converged to desired tolerance

References

• [34] Ogita, T. & Aishima, K., "Iterative refinement for symmetric eigenvalue decomposition", Japan J. Indust. Appl. Math. 35, 1007–1035 (2018).

refine_symmetric_eigen(A::AbstractMatrix{T}, Q0::Matrix, λ0::AbstractVector;
                       target_precision::Int=256,
                       max_iterations::Int=20,
                       tol::Real=0.0) where T

Refine an approximate symmetric eigenvalue decomposition A ≈ Q₀ Λ₀ Q₀^T to higher precision using the RefSyEv algorithm from Ogita & Aishima (2018).

This algorithm is specifically designed for symmetric (or Hermitian) matrices and handles multiple eigenvalues gracefully by using a special correction formula that avoids division by nearly-zero eigenvalue differences.

Algorithm (RefSyEv from Ogita & Aishima 2018)

For symmetric A with approximate decomposition A ≈ Q̂ D̂ Q̂^T:

1. Compute residuals: R = I - Q̂^T Q̂, S = Q̂^T A Q̂
2. Update eigenvalues: λ̃i = sii / (1 - r_ii) (Rayleigh quotient)
3. Compute threshold: δ = 2(‖S - D̃‖ + ‖A‖ · ‖R‖)
4. Compute correction matrix Ẽ:
   • For distinct eigenvalues (|λ̃i - λ̃j| > δ): ẽij = (sij + λ̃j · rij) / (λ̃j - λ̃i)
   • For multiple eigenvalues (|λ̃i - λ̃j| ≤ δ): ẽij = rij / 2
5. Update: Q' = Q̂(I + Ẽ)

The method converges quadratically, even for matrices with multiple eigenvalues.

Arguments

• A: Symmetric matrix (will be converted to target precision)
• Q0: Approximate orthogonal eigenvector matrix
• λ0: Approximate eigenvalues
• target_precision: Target precision in bits (default: 256)
• max_iterations: Maximum refinement iterations (default: 20)
• tol: Convergence tolerance (default: machine epsilon of target precision)

Returns

SymmetricEigenRefinementResult containing refined Q, λ and convergence info.

Example

using BallArithmetic, LinearAlgebra

# Create symmetric matrix
A = randn(10, 10); A = A + A'

# Compute eigen decomposition in Float64
F = eigen(Symmetric(A))
Q0, λ0 = F.vectors, F.values

# Refine to BigFloat precision
result = refine_symmetric_eigen(A, Q0, λ0; target_precision=256)

@show result.residual_norm       # Should be ≈ 10^-77
@show result.orthogonality_defect
@show result.converged

Notes

• For general (non-symmetric) matrices, use refine_schur_decomposition
• This algorithm is more efficient than general Schur refinement for symmetric matrices
• Handles multiple eigenvalues without instability

References

• [34] Ogita, T. & Aishima, K., "Iterative refinement for symmetric eigenvalue decomposition", Japan J. Indust. Appl. Math. 35, 1007–1035 (2018). Algorithm 1: RefSyEv.