API - Eigenvalues & SVD

RigorousEigenvaluesResult

Container returned by rigorous_eigenvalues bundling the midpoint eigenvector factorisation, the certified eigenvalue enclosures, 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)

Compute rigorous enclosures for the eigenvalues of A, following Ref. [7]. The returned RigorousEigenvaluesResult exposes both the interval enclosures and the norm bounds used during certification.

TODO: Using Miyajima's algorithm is overkill, may be worth using

References

• [7] Miyajima, JCAM 246, 9 (2012)

evbox(A::BallMatrix{T})

Backward-compatible wrapper returning only the vector of eigenvalue enclosures produced by rigorous_eigenvalues.

Rump2022aResult

Container for Rump2022a eigenvalue and eigenvector error bounds.

Extends the standard eigenvalue result with:

• Individual eigenvector error bounds
• Condition number estimates
• Residual-based refinements

rump_2022a_eigenvalue_bounds(A::BallMatrix; method=:standard, hermitian=false)

Compute verified error bounds for all eigenvalues and eigenvectors following Rump (2022a).

This method provides:

1. Individual eigenvalue enclosures with guaranteed containment
2. Eigenvector error bounds for each eigenpair
3. Condition number estimates for stability assessment
4. Residual-based verification

Arguments

• A::BallMatrix: Square matrix for eigenvalue problem
• method::Symbol: Verification method
  • :standard - Standard residual-based bounds (default)
  • :refined - Refined bounds using Gershgorin + residuals
  • :krawczyk - Krawczyk operator for sharper enclosures
• hermitian::Bool: Whether A is Hermitian (enables tighter bounds)

Method Description

Standard method:

For each eigenpair (λᵢ, vᵢ), computes:

1. Residual: rᵢ = Avᵢ - λᵢvᵢ
2. Eigenvalue bound: |λ̃ᵢ - λᵢ| ≤ ‖rᵢ‖/(1 - κᵢ*‖rᵢ‖)
3. Eigenvector bound: ‖ṽᵢ - vᵢ‖ ≤ κᵢ‖rᵢ‖/(1 - κᵢ‖rᵢ‖) where κᵢ is the condition number

Refined method:

Combines Gershgorin discs with residual bounds for tighter enclosures, especially effective when eigenvalues are clustered.

Krawczyk method:

Uses interval Newton-Krawczyk operator for quadratic convergence in eigenvector refinement.

Returns

Rump2022aResult containing verified bounds.

Example

A = BallMatrix([2.0 1.0; 1.0 2.0])
result = rump_2022a_eigenvalue_bounds(A; hermitian=true)

# Access results
λ = result.eigenvalues  # Eigenvalue balls
κ = result.condition_numbers  # Condition numbers
err = result.eigenvector_errors  # Eigenvector error bounds

Reference

• Rump, S.M. (2022), "Verified Error Bounds for All Eigenvalues and Eigenvectors of a Matrix", SIAM J. Matrix Anal. Appl. 43(4):1736–1754. DOI: 10.1137/21M1451440

RumpLange2023Result

Container for RumpLange2023 eigenvalue cluster bounds.

Emphasizes fast computation with clustering information:

• Cluster structure identification
• Per-cluster error bounds
• Fast bounds optimized for clustered spectra

Fields

• eigenvectors::VT: Approximate eigenvectors (as ball matrix)
• eigenvalues::ΛT: Certified eigenvalue enclosures
• cluster_assignments::Vector{Int}: Cluster assignments (cluster index for each eigenvalue)
• cluster_bounds::Vector{Ball{T, T}}: Cluster bounds (interval enclosure for each cluster)
• num_clusters::Int: Number of clusters identified
• cluster_residuals::Vector{T}: Per-cluster residual norms
• cluster_separations::Vector{T}: Per-cluster separation gaps
• cluster_sizes::Vector{Int}: Cluster sizes
• verified::Bool: Overall verification status

rump_lange_2023_cluster_bounds(A::BallMatrix; hermitian=false, cluster_tol=1e-6, fast=true)

Compute fast error bounds for eigenvalues with emphasis on cluster structure, following Rump & Lange (2023).

This method excels when eigenvalues form clusters, providing:

1. Fast identification of eigenvalue clusters
2. Tight bounds within each cluster
3. Optimized computation exploiting cluster structure
4. Scaling to large matrices via cluster-wise processing

Arguments

• A::BallMatrix: Square matrix for eigenvalue problem
• hermitian::Bool: Whether A is Hermitian (enables faster algorithms)
• cluster_tol::Real: Tolerance for cluster identification (default: 1e-6)
• fast::Bool: Use fast approximations vs. rigorous bounds (default: true)

Method Description

Cluster identification:

Uses Gershgorin discs with connectivity analysis to identify clusters of eigenvalues that are close together. Two eigenvalues belong to the same cluster if their Gershgorin discs overlap.

Per-cluster bounds:

For each cluster Ck with eigenvalues {λᵢ}ᵢ∈Ck:

1. Compute cluster interval: [min λᵢ - δᵢ, max λᵢ + δᵢ]
2. Refine using projected residuals within cluster
3. Apply separation bounds between clusters

Fast mode:

When fast=true, uses:

• Single power iteration for norms (vs. convergence)
• Simplified residual bounds
• Cluster-level (vs. individual) error propagation

Results in ~10x speedup with typically <2x looser bounds.

Returns

RumpLange2023Result containing cluster structure and bounds.

Example

Example usage for a matrix with two eigenvalue clusters:

• Create interval matrix with clustered spectrum
• Call rump_lange_2023_cluster_bounds(A; hermitian=true)
• Result contains cluster assignments and per-cluster bounds

Performance Notes

• For n×n matrix: O(n²) flops in fast mode, O(n³) in rigorous mode
• Cluster count << n gives significant speedup
• Most effective when cluster separation >> cluster width

Reference

• Rump, S.M. & Lange, M. (2023), "Fast Computation of Error Bounds...", SIAM J. Matrix Anal. Appl., to appear

refine_cluster_bounds(result::RumpLange2023Result, A::BallMatrix; iterations=1)

Refine cluster bounds using iterative residual computation.

Takes an existing RumpLange2023Result and performs additional refinement iterations to tighten the bounds, particularly for well-separated clusters.

Arguments

• result: Initial cluster bound result
• A: Original ball matrix
• iterations: Number of refinement iterations (default: 1)

Returns

New RumpLange2023Result with refined bounds.

RigorousGeneralizedEigenvaluesResult

Container returned by rigorous_generalized_eigenvalues bundling the midpoint eigenvector factorisation, the certified eigenvalue enclosures, 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 rigorous enclosures for the eigenvalues of the generalised problem A * x = λ * B * x, following Ref. [7]. The returned RigorousGeneralizedEigenvaluesResult exposes both the interval enclosures and the norm bounds used during certification.

References

• [7] Miyajima, JCAM 246, 9 (2012)

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

Backward-compatible wrapper returning only the vector of eigenvalue enclosures produced by rigorous_generalized_eigenvalues.

GEVResult

Result structure for verified generalized eigenvalue problem.

All numeric fields use Float64 precision. This struct is not currently parametric; extension to other numeric types would require making it GEVResult{T}.

Fields

• success::Bool: Whether verification succeeded
• eigenvalue_intervals::Vector{Tuple{Float64, Float64}}: Verified intervals [λ̃ᵢ - ηᵢ, λ̃ᵢ + ηᵢ] for each eigenvalue
• eigenvector_centers::Matrix{Float64}: Approximate eigenvectors (centers)
• eigenvector_radii::Vector{Float64}: Verified radii ξᵢ for eigenvector balls
• beta::Float64: Preconditioning factor β ≥ √‖B⁻¹‖₂
• global_bound::Float64: Global eigenvalue bound δ̂ (Theorem 4)
• individual_bounds::Vector{Float64}: Individual eigenvalue bounds ε (Theorem 5)
• separation_bounds::Vector{Float64}: Separation bounds η (Lemma 2)
• residual_norm::Float64: Norm of residual matrix ‖Rg‖₂
• message::String: Diagnostic message (especially if success = false)

Interpretation

• If success = true: All eigenvalue intervals are guaranteed to contain exactly one true eigenvalue, and all eigenvector balls contain the corresponding normalized eigenvector, rigorously accounting for all matrices in the input intervals [A] and [B].
• If success = false: Check message for diagnostic information. Common failures:
  • Approximate eigenvectors not sufficiently orthogonal (‖I - Gg‖₂ >= 1)
  • Eigenvalues too clustered to separate
  • B not positive definite

verify_generalized_eigenpairs(A::BallMatrix, B::BallMatrix, X̃::Matrix, λ̃::Vector) -> GEVResult

Verify all eigenpairs of the generalized eigenvalue problem Ax = λBx.

Implements Algorithm 1 from Miyajima et al. (2010).

Numeric Type Support

Currently supports Float64 only. All computations and error bounds are performed using Float64 arithmetic with IEEE 754 double precision. The implementation uses Float64-specific rounding error constants for rigorous verification.

Extension to BigFloat would require:

• Parametric GEVResult{T} struct
• Type-dependent unit roundoff (eps(T))
• Modified error analysis for arbitrary precision

The mathematical algorithms (Theorems 4, 5, 7, 10) are precision-independent, but this implementation is optimized for Float64 hardware arithmetic.

Arguments

• A::BallMatrix: Symmetric interval matrix (n×n, Float64 elements)
• B::BallMatrix: Symmetric positive definite interval matrix (n×n, Float64 elements)
• X̃::Matrix: Approximate eigenvectors (n×n), typically from eigen(A.c, B.c)
• λ̃::Vector: Approximate eigenvalues (n), assumed sorted

Returns

• GEVResult with verified eigenvalue intervals and eigenvector balls

Algorithm (4 steps)

1. Compute β ≥ √‖B⁻¹‖₂ using Theorem 10
2. Compute global bound δ̂ and individual bounds ε using Theorems 4, 5
3. Determine separation bounds η using Lemma 2
4. Compute eigenvector bounds ξ using Theorem 7

Complexity

O(12n³) dominated by matrix multiplications with interval arithmetic

Verification Guarantees

When success = true:

• Each interval [λ̃ᵢ - ηᵢ, λ̃ᵢ + ηᵢ] contains exactly one true eigenvalue
• Each ball B(x̃⁽ⁱ⁾, ξᵢ) contains the normalized true eigenvector
• Results are rigorous for ALL matrices in the intervals [A] and [B]

Example

using LinearAlgebra

A = BallMatrix([4.0 1.0; 1.0 3.0], fill(1e-10, 2, 2))
B = BallMatrix([2.0 0.5; 0.5 2.0], fill(1e-10, 2, 2))

F = eigen(Symmetric(A.c), Symmetric(B.c))
result = verify_generalized_eigenpairs(A, B, F.vectors, F.values)

if result.success
    println("Eigenvalue 1 ∈ ", result.eigenvalue_intervals[1])
    println("Eigenvalue 2 ∈ ", result.eigenvalue_intervals[2])
    println("Eigenvector radii: ", result.eigenvector_radii)
end

References

Miyajima, S., Ogita, T., Rump, S. M., Oishi, S. (2010). "Fast Verification for All Eigenpairs in Symmetric Positive Definite Generalized Eigenvalue Problems". Reliable Computing 14, pp. 24-45.

compute_beta_bound(B::BallMatrix) -> Float64

Compute verified upper bound β ≥ √‖B⁻¹‖₂ using Theorem 10.

This function efficiently computes a bound on the square root of the 2-norm of B⁻¹ using Cholesky factorization and an approximate inverse.

Numeric Type Support

Currently supports Float64 only. The error analysis uses Float64-specific rounding error constants (eps(Float64)). Extension to BigFloat would require:

• Type-parametric unit roundoff: eps(T) instead of eps(Float64)
• Appropriate error analysis for higher precision
• Parametric GEVResult struct

The mathematical algorithm itself is not Float64-specific, but the rigorous error bounds in the implementation assume IEEE 754 double precision arithmetic.

Arguments

• B::BallMatrix: Symmetric positive definite interval matrix (Float64 elements)

Returns

• β::Float64: Upper bound on √‖B⁻¹‖₂

Algorithm (Theorem 10)

1. Compute Cholesky factorization B ≈ LL^T
2. Compute approximate inverse X_L ≈ L⁻¹
3. Use interval arithmetic to bound error
4. Return β = √((α₁α∞)/(1 - α₁α∞αC))

Complexity

O(n³) for Cholesky and inverse

Example

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

RigorousSVDResult

Container returned by rigorous_svd bundling the midpoint factorisation, the certified singular-value enclosures, and the block-diagonal refinement obtained from miyajima_vbd. 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 miyajima_vbd 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 = MiyajimaM1(), apply_vbd = true)

Backward-compatible wrapper returning only the vector of singular-value enclosures produced by rigorous_svd. New code should prefer rigorous_svd directly to access the additional certification data. The optional method and apply_vbd flags mirror those in rigorous_svd.

MiyajimaVBDResult

Container returned by miyajima_vbd encapsulating the data produced by the verified block diagonalisation (VBD) step. The fields contain the basis that block diagonalises the midpoint matrix, the transformed enclosure, its block-diagonal truncation, the rigorous remainder, and the Gershgorin clusters that certify how the spectrum groups together.

miyajima_vbd(A::BallMatrix; hermitian = false)

Perform Miyajima's verified block diagonalisation (VBD) on the square ball matrix A. The midpoint matrix is reduced either by an eigenvalue decomposition (when hermitian = true) or by a unitary Schur form (for the general case). The enclosure is transported to that basis, the Gershgorin discs are clustered, and a block-diagonal truncation together with a rigorous remainder is produced.

Overlapping discs are grouped via their connectivity graph so that each cluster becomes contiguous after a basis permutation. The remainder bound combines the classical Collatz estimate with a block-separation bound that exploits the verified gaps between clusters.

When hermitian = true the routine expects A to be Hermitian and the resulting eigenvalues and intervals are real. Otherwise the Schur form is used and the clusters are discs in the complex plane.

refine_svd_bounds_with_vbd(result::RigorousSVDResult)

Attempt to refine singular value bounds using VBD isolation information.

For singular values whose squared values fall in isolated Gershgorin clusters, we can potentially tighten the bounds using Miyajima's Theorem 11.

Returns a new RigorousSVDResult with potentially tighter bounds, or the original result if no refinement is possible.

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 with Float64 precision
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

Arguments

• A: Input ball matrix
• tolerance: Target tolerance for max(rad(σᵢ))
• method: SVD certification method (MiyajimaM1, MiyajimaM4, RumpOriginal)
• 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.

RigorousSpectralProjectorsResult

Container returned by miyajima_spectral_projectors encapsulating the rigorously computed spectral projectors for each eigenvalue cluster identified by verified block diagonalization (VBD).

Each projector P_k is a ball matrix satisfying:

• P_k^2 ≈ P_k (idempotency)
• ∑_k P_k ≈ I (resolution of identity)
• P_i * P_j ≈ 0 for i ≠ j (orthogonality)
• A * P_k ≈ P_k * A * P_k (invariance)

The projectors are constructed from the basis that block-diagonalizes the matrix, restricted to each spectral cluster.

miyajima_spectral_projectors(A::BallMatrix; hermitian=false, verify_invariance=true)

Compute rigorous enclosures for spectral projectors corresponding to each eigenvalue cluster identified by Miyajima's verified block diagonalization (VBD).

The method follows the approach from Ref. [8]:

1. Apply VBD to obtain basis V that block-diagonalizes A
2. For each cluster k, extract columns V[:, cluster_k]
3. Construct projector P_k = V[:, cluster_k] * V[:, cluster_k]' as ball matrix
4. Verify idempotency, orthogonality, and resolution of identity

When hermitian = true, the basis is computed via eigendecomposition and projectors are Hermitian. Otherwise, the Schur basis is used.

When verify_invariance = true, additionally verifies that A * P_k ≈ P_k * A * P_k for each projector, confirming that the columns of P_k span an invariant subspace.

Arguments

• A::BallMatrix: Square ball matrix whose spectral projectors to compute
• hermitian::Bool = false: Whether to assume A is Hermitian
• verify_invariance::Bool = true: Whether to verify invariant subspace property

Returns

RigorousSpectralProjectorsResult containing the projectors and verification data.

Example

using BallArithmetic, LinearAlgebra

# Create a matrix with clustered eigenvalues
A = BallMatrix(Diagonal([1.0, 1.1, 5.0, 5.1]))

# Compute projectors
result = miyajima_spectral_projectors(A; hermitian=true)

# Access projectors
P1 = result[1]  # Projector for first cluster (eigenvalues ≈ 1.0, 1.1)
P2 = result[2]  # Projector for second cluster (eigenvalues ≈ 5.0, 5.1)

# Verify properties
@assert result.idempotency_defect < 1e-10
@assert result.orthogonality_defect < 1e-10

References

• [8] Miyajima, SIAM J. Matrix Anal. Appl. 35, 1205–1225 (2014)

compute_invariant_subspace_basis(proj_result::RigorousSpectralProjectorsResult, k::Int)

Extract an orthonormal basis for the invariant subspace corresponding to cluster k from the spectral projector result.

Returns a BallMatrix whose columns span the invariant subspace associated with the k-th eigenvalue cluster.

verify_projector_properties(proj_result::RigorousSpectralProjectorsResult; tol=1e-10)

Verify that all projector properties hold within the specified tolerance. Returns true if all properties are satisfied, false otherwise.

Checks:

1. Idempotency: ‖Pk^2 - Pk‖₂ < tol for all k
2. Orthogonality: ‖Pi * Pj‖₂ < tol for all i ≠ j
3. Resolution: ‖∑k Pk - I‖₂ < tol
4. Invariance: ‖APk - PkA*P_k‖₂ < tol for all k (if computed)

projector_condition_number(proj_result::RigorousSpectralProjectorsResult, k::Int)

Estimate the condition number of the k-th spectral projector based on the gap between eigenvalue clusters.

A small gap indicates potential ill-conditioning of the projector.

RigorousBlockSchurResult

Container returned by rigorous_block_schur encapsulating a verified block Schur decomposition A ≈ Q * T * Q' where:

• Q is an approximately orthogonal/unitary transformation (as ball matrix)
• T is in block upper quasi-triangular form
• Diagonal blocks correspond to eigenvalue clusters
• Off-diagonal blocks are rigorously bounded

The decomposition satisfies:

• A ≈ Q * T * Q' with rigorous residual bounds
• Q' * Q ≈ I with rigorous orthogonality defect bounds
• Each diagonal block T[cluster_k, cluster_k] contains eigenvalues from cluster k

rigorous_block_schur(A::BallMatrix; hermitian=false, block_structure=:quasi_triangular)

Compute a rigorous block Schur decomposition A ≈ Q * T * Q' where Q is orthogonal/unitary and T is in block form determined by eigenvalue clustering.

The method follows Miyajima's VBD framework:

1. Apply verified block diagonalization to identify eigenvalue clusters
2. Construct orthogonal basis Q from the diagonalizing transformation
3. Transform matrix to block form T = Q' * A * Q
4. Verify orthogonality of Q and residual bounds

Arguments

• A::BallMatrix: Square ball matrix to decompose
• hermitian::Bool = false: Whether to assume A is Hermitian
• block_structure::Symbol = :quasi_triangular: Block structure to compute
  • :diagonal: Keep only diagonal blocks (same as VBD)
  • :quasi_triangular: Keep upper triangular block structure
  • :full: Keep all blocks (no truncation)

Returns

RigorousBlockSchurResult containing the decomposition and verification data.

Example

using BallArithmetic, LinearAlgebra

# Create a matrix with clustered eigenvalues
A = BallMatrix([2.0 0.1 0.05 0.02;
                0.1 2.1 0.03 0.01;
                0.05 0.03 5.0 0.15;
                0.02 0.01 0.15 5.1])

# Compute block Schur form
result = rigorous_block_schur(A; hermitian=true)

# Access components
Q = result.Q
T = result.T

# Verify decomposition
@assert result.residual_norm < 1e-10
@assert result.orthogonality_defect < 1e-10

References

• [8] Miyajima, SIAM J. Matrix Anal. Appl. 35, 1205–1225 (2014)
• [11] Miyajima, Japan J. Indust. Appl. Math. 31, 513–539 (2014)

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)

Verify that the block Schur decomposition satisfies all required properties within the specified tolerance.

Checks:

1. Residual: ‖A - QTQ'‖₂ < tol
2. Orthogonality: ‖Q'*Q - I‖₂ < tol
3. Block structure preserved (if applicable)

Returns true if all properties are satisfied, false otherwise.

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.

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

Refine the (i,j) off-diagonal block by solving the block Sylvester equation with Miyajima's verified solver.

For i < j, solves T_ii' * Y - Y * T_jj' = T_ij' to obtain a refined enclosure for the (i,j) off-diagonal block.

Returns a BallMatrix with rigorous enclosure for the refined block.

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

For clusters i < j, compute the right-hand side for the block Sylvester equation that would refine the (i,j) off-diagonal block.

Given T_ii * X + X * T_jj = C, returns the matrix C that should equal T[cluster_i, cluster_j] if the block Schur form were exact.

SchurSpectralProjectorResult{T}

Result of verified spectral projector computation from Schur decomposition.

Fields

• projector::BallMatrix{T}: Rigorous enclosure of the spectral projector P
• schur_projector::BallMatrix{T}: Projector in Schur coordinates
• coupling_matrix::BallMatrix{T}: Solution Y to the Sylvester equation T₁₁Y - YT₂₂ = T₁₂
• eigenvalue_separation::T: Lower bound on min|λᵢ - λⱼ| for i ∈ S, j ∉ S
• projector_norm::T: Upper bound on ‖P‖₂
• idempotency_defect::T: Upper bound on ‖P² - P‖₂
• schur_basis::Matrix{T}: Unitary matrix Q from A = QTQ^†
• schur_form::Matrix{T}: Upper triangular T from A = QTQ^†
• cluster_indices::UnitRange{Int}: Indices in Schur form corresponding to projected eigenvalues

compute_spectral_projector_schur(A::BallMatrix, cluster_indices::UnitRange{Int};
                                 verify_idempotency::Bool=true)

Compute rigorous enclosure of the spectral projector for eigenvalues at specified Schur form indices using verified Sylvester equation solver.

Mathematical Background

Given a matrix A with Schur decomposition A = QTQ^†, where T is upper triangular, partition T into blocks:

T = [T₁₁  T₁₂]
    [0    T₂₂]

where T₁₁ contains the eigenvalues we want to project onto (diagonal entries at positions given by cluster_indices).

The spectral projector (Riesz projector) onto the eigenspace corresponding to eigenvalues of T₁₁ is:

P = Q * P_Schur * Q^†

where P_Schur is the projector in Schur coordinates:

P_Schur = [I   Y]
          [0   0]

and Y solves the Sylvester equation:

T₁₁*Y - Y*T₂₂ = T₁₂

This is solved rigorously using triangular_sylvester_miyajima_enclosure.

Algorithm

1. Compute Schur decomposition: A = QTQ^† (using approximate Schur for A.c)
2. Extract blocks T₁₁, T₁₂, T₂₂ based on cluster_indices
3. Solve Sylvester equation rigorously: T₁₁Y - YT₂₂ = T₁₂
4. Construct P_Schur = [I Y; 0 0] with interval arithmetic
5. Transform back: P = Q * P_Schur * Q^† with interval arithmetic
6. Verify idempotency: ‖P² - P‖₂ < tol

Arguments

• A::BallMatrix: Square interval matrix (n × n)
• cluster_indices::UnitRange{Int}: Indices of eigenvalues to project onto (corresponds to positions in Schur form, typically 1:k for first k eigenvalues)
• verify_idempotency::Bool=true: Check that ‖P² - P‖₂ is small

Returns

SchurSpectralProjectorResult containing the verified projector and diagnostics.

Examples

Project onto first two eigenvalues

using BallArithmetic, LinearAlgebra

# Matrix with small uncertainties
A = BallMatrix([1.0 2.0 0.0; 0.0 3.0 1.0; 0.0 0.0 5.0], fill(1e-10, 3, 3))

# Compute projector onto eigenvalues 1.0 and 3.0 (first two in Schur form)
result = compute_spectral_projector_schur(A, 1:2)

P = result.projector
@show result.eigenvalue_separation  # Gap to third eigenvalue
@show result.idempotency_defect     # ‖P² - P‖₂

# Project a vector
v = BallVector([1.0, 2.0, 3.0])
v_projected = P * v

For a clustered matrix

# Upper triangular matrix with eigenvalues [1.0, 1.1, 5.0, 5.1]
A = BallMatrix(triu(randn(4,4)) .+ Diagonal([1.0, 1.1, 5.0, 5.1]))

# After sorting eigenvalues by Schur decomposition, project onto first cluster
result = compute_spectral_projector_schur(A, 1:2)

Notes

• Requires non-zero eigenvalue separation: min|λᵢ - λⱼ| for i ∈ cluster, j ∉ cluster
• The Sylvester equation solver may fail if separation is too small
• Complexity: O(n³) for Schur decomposition + O(k²(n-k)²) for Sylvester solver
• For hermitian matrices, use compute_spectral_projector_hermitian instead (more efficient)
• The projector P satisfies: P² ≈ P and PA ≈ AP (modulo idempotency defect)
• Caveat: The Schur decomposition of the midpoint matrix is not itself rigorously enclosed (Q and T carry zero radii). The Sylvester coupling Y is rigorously bounded, but the final projector P = Q·PSchur·Q† does not account for the Schur approximation error. The idempotency and invariance checks verify self-consistency of the result, but do not guarantee that the true spectral projector lies within the computed ball. For a fully rigorous enclosure, use [`miyajimaspectral_projectors`](@ref) instead.

References

• Kato, T. "Perturbation Theory for Linear Operators" (1995), Chapter II.4
• Stewart, G. W., Sun, J. "Matrix Perturbation Theory" (1990), Chapter V
• Miyajima, S. "Fast enclosure for all eigenvalues in generalized eigenvalue problems" SIAM J. Matrix Anal. Appl. 35, 1205–1225 (2014)

See Also

• triangular_sylvester_miyajima_enclosure: Verified Sylvester solver used internally
• project_vector_spectral: Convenient interface to project vectors using this result

compute_spectral_projector_schur(A::BallMatrix, target_indices::AbstractVector{Int};
                                 verify_idempotency::Bool=true,
                                 schur_data=nothing)

Compute rigorous spectral projector for eigenvalues at arbitrary Schur form positions.

Unlike the UnitRange method, this accepts any subset of eigenvalue indices. Internally, the Schur form is reordered (via ordschur) so the target eigenvalues occupy positions 1:k, then the standard Sylvester-based projector pipeline is applied.

Arguments

• A::BallMatrix: Square interval matrix (n × n)
• target_indices::AbstractVector{Int}: Positions of eigenvalues in Schur form to project onto
• verify_idempotency::Bool=true: Check ‖P² - P‖₂
• schur_data=nothing: Optional pre-computed (Q, T) Schur factors to skip internal Schur

Returns

SchurSpectralProjectorResult with the reordered Schur factors and verified projector.

Example

A = BallMatrix(randn(5, 5))
# Project onto 2nd and 4th eigenvalues (arbitrary positions)
result = compute_spectral_projector_schur(A, [2, 4])

compute_spectral_projector_hermitian(A::BallMatrix, cluster_indices::UnitRange{Int})

Compute spectral projector for Hermitian matrix (simplified, no Sylvester equation needed).

For Hermitian matrices, eigenvectors are orthogonal, so the spectral projector is simply:

P = V_S * V_S^†

where V_S are the eigenvectors corresponding to the selected eigenvalues.

This is more efficient than the general Schur-based method since no Sylvester equation needs to be solved.

Arguments

• A::BallMatrix: Hermitian interval matrix
• cluster_indices::UnitRange{Int}: Indices of eigenvalues to project onto (after sorting)

Returns

SchurSpectralProjectorResult with the verified projector.

Example

# Symmetric matrix
A = BallMatrix([4.0 1.0; 1.0 3.0], fill(1e-10, 2, 2))

# Project onto first eigenvalue
result = compute_spectral_projector_hermitian(A, 1:1)
P = result.projector

project_vector_spectral(v::BallVector, result::SchurSpectralProjectorResult)

Project interval vector onto eigenspace using precomputed spectral projector.

Arguments

• v::BallVector: Interval vector to project
• result::SchurSpectralProjectorResult: Precomputed projector from compute_spectral_projector_schur or compute_spectral_projector_hermitian

Returns

BallVector containing P*v with rigorous error bounds.

Example

A = BallMatrix([1.0 2.0; 0.0 3.0], fill(1e-10, 2, 2))
result = compute_spectral_projector_schur(A, 1:1)

v = BallVector([1.0, 2.0], [1e-10, 1e-10])
v_projected = project_vector_spectral(v, result)

project_vector_spectral(v::AbstractVector, result::SchurSpectralProjectorResult)

Project standard vector onto eigenspace (non-verified version).

Arguments

• v::AbstractVector: Vector to project
• result::SchurSpectralProjectorResult: Precomputed projector

Returns

Projected vector P*v (center value only, no error bounds).

verify_spectral_projector_properties(result::SchurSpectralProjectorResult, A::BallMatrix;
                                     tol::Real=1e-10)

Verify that the computed spectral projector satisfies expected mathematical properties.

Checks

1. Idempotency: ‖P² - P‖₂ < tol
2. Bounded norm: ‖P‖₂ < ∞ (typically ‖P‖₂ ≤ κ(V) where κ is condition number of eigenvectors)
3. Commutation (if A normal): ‖AP - PA‖₂ < tol
4. Eigenvalue separation: min|λᵢ - λⱼ| > 0 for i ∈ S, j ∉ S

Arguments

• result::SchurSpectralProjectorResult: Computed projector
• A::BallMatrix: Original matrix
• tol::Real=1e-10: Tolerance for property verification

Returns

• true if all properties satisfied, false otherwise

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

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

refine_schur_decomposition(A::AbstractMatrix{T}, Q0::Matrix, T0::Matrix;
                            target_precision::Int=256,
                            max_iterations::Int=20,
                            tol::Real=0.0) where T

Refine an approximate Schur decomposition A ≈ Q₀ T₀ Q₀^H to higher precision.

Starting from an approximate Schur decomposition computed in Float64, this function iteratively refines Q and T to achieve the accuracy of target_precision bits. The algorithm uses Newton-like iterations that converge quadratically.

Algorithm (from Bujanović et al. 2022, Algorithm 4)

Given A = Q₀ T₀ Q₀^H + E where E is the initial error:

1. Pre-computation (2 high-precision matrix multiplications):

   • Compute AQ = A * Q and Q^HA = Q^H * A in target precision

2. Iteration (4 high-precision matrix multiplications per step):

   • Compute residual R = Q^H * A_Q - T
   • Extract diagonal correction: T_new = T + diag(R)
   • Solve Sylvester equations for off-diagonal corrections
   • Apply Newton-Schulz to re-orthogonalize Q
   • Repeat until convergence

The method achieves 10-20× speedup over computing Schur directly in high precision.

Arguments

• A: Original matrix (will be converted to target precision)
• Q0: Approximate orthogonal/unitary factor from Schur decomposition
• T0: Approximate upper triangular factor from Schur decomposition
• target_precision: Target precision in bits (default: 256 for BigFloat)
• max_iterations: Maximum refinement iterations (default: 20)
• tol: Convergence tolerance (default: machine epsilon of target precision)

Returns

SchurRefinementResult containing refined Q, T and convergence information.

Example

using BallArithmetic, LinearAlgebra

# Compute Schur in Float64
A = randn(100, 100)
F = schur(A)
Q0, T0 = F.Z, F.T

# Refine to BigFloat precision
result = refine_schur_decomposition(A, Q0, T0; target_precision=256)

# Check residual
@show result.residual_norm      # Should be ≈ 10^-77 for 256-bit precision
@show result.orthogonality_defect

References

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

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

Compute rigorous Schur decomposition with BigFloat precision using iterative refinement.

This combines:

1. Fast approximate Schur decomposition in Float64 (or use schur_seed)
2. Iterative refinement to BigFloat precision
3. Rigorous error certification via 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}: Rigorous enclosure of Schur basis
• T_ball::BallMatrix{BigFloat}: Rigorous enclosure of Schur form
• result::SchurRefinementResult: Refinement diagnostics

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

• [9] for the iterative refinement algorithm
• [12] for rigorous error certification

certify_schur_decomposition(A::BallMatrix, result::SchurRefinementResult)

Wrap a SchurRefinementResult in BallMatrix form with rigorous error bounds.

Given an input ball matrix A and a SchurRefinementResult from refine_schur_decomposition, construct BallMatrix enclosures of the Schur factors Q and T. The radii are derived from:

• The orthogonality defect ‖Q^H Q - I‖₂ (contributes to Q radius)
• The residual norm ‖stril(Q^H A Q)‖₂ / ‖A‖_F (contributes to T radius)
• The input uncertainty rad(A) (propagated to both)

The BigFloat precision is automatically matched to the precision of result.Q.

Arguments

• A::BallMatrix: The original input matrix (provides rad(A) for uncertainty propagation)
• result::SchurRefinementResult: Refinement output with Q, T, and diagnostics

Returns

(Q_ball, T_ball) — BallMatrix enclosures of the Schur factors.

Example

A = BallMatrix(randn(ComplexF64, 5, 5), fill(1e-10, 5, 5))
F = schur(mid(A))
result = refine_schur_decomposition(mid(A), F.Z, F.T; target_precision=256)
Q_ball, T_ball = certify_schur_decomposition(A, result)

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)

This converges quadratically if ‖Q^H Q - I‖₂ < 1.

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

References

• [9] Section 2.2.2 (Newton-Schulz iteration)

NOTE: This is an ORACLE computation - convergence is checked via the defect

(‖Q^H Q - I‖_F < tol), which serves as an A POSTERIORI verification. The

intermediate computations don't need directed rounding because the final

defect measurement validates the orthogonalization quality.

project_onto_eigenspace(v::AbstractVector, V::AbstractMatrix, indices::UnitRange{Int};
                        hermitian::Bool=false, left_eigenvectors::Union{Nothing, AbstractMatrix}=nothing)

Project vector v onto the eigenspace spanned by eigenvectors at the specified indices.

Mathematical Background

For a diagonalizable matrix A = VΛV⁻¹, the Riesz projection onto the eigenspace corresponding to eigenvalues λ_k (k ∈ indices) is:

Hermitian/Normal case (hermitian=true or V orthogonal):

P_S v = V_S (V_S^† V_S)⁻¹ V_S^† v = V_S V_S^† v

where V_S = V[:, indices] and V^† is the conjugate transpose.

Non-normal case (requires left eigenvectors):

P_S v = V_S W_S^† v

where W_S contains the corresponding left eigenvectors (rows of V⁻¹).

Arguments

• v::AbstractVector: Vector to project (length n)
• V::AbstractMatrix: Matrix of right eigenvectors (n × n)
• indices::UnitRange{Int}: Column indices of eigenvectors to project onto
• hermitian::Bool=false: If true, assumes V is orthogonal (V^† V = I)
• left_eigenvectors::Union{Nothing, AbstractMatrix}=nothing: Left eigenvectors W = V⁻¹ (required for non-normal matrices when hermitian=false)

Returns

Projected vector P_S v of length n.

Examples

Symmetric/Hermitian case

using LinearAlgebra, BallArithmetic

# Symmetric matrix with clustered eigenvalues
A = [4.0 1.0 0.0; 1.0 4.0 0.0; 0.0 0.0 10.0]
F = eigen(Symmetric(A))

# Project onto first two eigenvectors (eigenvalues ≈ 3, 5)
v = [1.0, 2.0, 3.0]
v_proj = project_onto_eigenspace(v, F.vectors, 1:2; hermitian=true)

Non-normal case (requires left eigenvectors)

# Non-normal matrix
A = [1.0 1.0; 0.0 2.0]  # Upper triangular, non-symmetric
F = eigen(A)

# Compute left eigenvectors
V_inv = inv(F.vectors)  # Rows are left eigenvectors

# Project onto first eigenspace
v = [1.0, 1.0]
v_proj = project_onto_eigenspace(v, F.vectors, 1:1;
                                  hermitian=false, left_eigenvectors=V_inv)

Notes

• For hermitian=true, uses the simple formula P = VS VS^†
• For hermitian=false without left eigenvectors, attempts V⁻¹ computation (may fail)
• For non-normal matrices, providing left_eigenvectors avoids computing V⁻¹
• The projection is exact in exact arithmetic but subject to rounding errors
• For verified computation, use verified_project_onto_eigenspace with BallMatrix

References

• Kato, T. "Perturbation theory for linear operators" (1995), Chapter II
• Stewart, G. W., Sun, J. "Matrix Perturbation Theory" (1990), Section V.3

project_onto_schur_subspace(v::AbstractVector, Q::AbstractMatrix, indices::UnitRange{Int})

Project vector v onto the invariant subspace spanned by Schur vectors at specified indices.

Mathematical Background

For a Schur decomposition A = QTQ^†, the Schur vectors (columns of Q) span nested invariant subspaces. The projection onto the subspace spanned by Q[:, indices] is:

P_S v = Q_S Q_S^† v

where Q_S = Q[:, indices]. Since Q is unitary (Q^† Q = I), this is always well-defined.

Arguments

• v::AbstractVector: Vector to project (length n)
• Q::AbstractMatrix: Unitary Schur matrix (n × n, with Q^† Q = I)
• indices::UnitRange{Int}: Column indices of Schur vectors to project onto

Returns

Projected vector P_S v of length n.

Example

using LinearAlgebra, BallArithmetic

# Matrix with complex eigenvalues
A = [0.0 1.0; -1.0 0.0]  # Rotation matrix
F = schur(A)

# Project onto first Schur vector
v = [1.0, 1.0]
v_proj = project_onto_schur_subspace(v, F.Z, 1:1)

Notes

• Schur vectors are always orthonormal, so projection is straightforward
• The subspace spanned by Q[:, 1:k] is an invariant subspace of A
• For real Schur form, complex eigenvalue pairs share a 2D invariant subspace
• This is simpler than eigenspace projection for non-normal matrices

References

• Golub, G. H., Van Loan, C. F. "Matrix Computations" (2013), Section 7.6

verified_project_onto_eigenspace(v::BallVector, V::BallMatrix, indices::UnitRange{Int};
                                 hermitian::Bool=false)

Rigorously verified projection of interval vector onto eigenspace with error bounds.

Similar to project_onto_eigenspace but uses interval arithmetic to provide rigorous bounds on the projection result accounting for uncertainties in both the vector and eigenvector matrix.

Arguments

• v::BallVector: Interval vector to project
• V::BallMatrix: Interval matrix of eigenvectors
• indices::UnitRange{Int}: Column indices of eigenvectors
• hermitian::Bool=false: If true, assumes V is orthogonal

Returns

BallVector containing the projected vector with rigorous error bounds.

Example

using BallArithmetic, LinearAlgebra

# Matrix with small uncertainties
A = BallMatrix([4.0 1.0; 1.0 3.0], fill(1e-10, 2, 2))
F = eigen(Symmetric(A.c))

# Vector with uncertainties
v = BallVector([1.0, 2.0], [1e-10, 1e-10])

# Verified projection
V_ball = BallMatrix(F.vectors)
v_proj = verified_project_onto_eigenspace(v, V_ball, 1:1; hermitian=true)

Notes

• Uses interval arithmetic throughout for rigorous error bounds
• The Gram system G·c = V_S†·v is solved rigorously via krawczyk_linear_system; if Krawczyk does not verify, a Neumann-series residual bound is used as fallback
• Result contains all possible projections for vectors/eigenvectors in input intervals
• More expensive than non-verified version due to interval arithmetic overhead
• Currently only supports hermitian case (non-normal requires Sylvester solver)

compute_eigenspace_projector(V::AbstractMatrix, indices::UnitRange{Int};
                            hermitian::Bool=false,
                            left_eigenvectors::Union{Nothing, AbstractMatrix}=nothing)

Compute the projection matrix P_S onto the eigenspace spanned by specified eigenvectors.

Returns the n × n projection matrix PS such that PS v projects any vector v onto the specified eigenspace.

Arguments

• V::AbstractMatrix: Matrix of right eigenvectors (n × n)
• indices::UnitRange{Int}: Column indices of eigenvectors defining the subspace
• hermitian::Bool=false: If true, assumes V is orthogonal
• left_eigenvectors::Union{Nothing, AbstractMatrix}=nothing: Left eigenvectors W = V⁻¹

Returns

n × n projection matrix PS with PS² = P_S.

Example

F = eigen(A)
P = compute_eigenspace_projector(F.vectors, 1:2; hermitian=true)

# Now can project multiple vectors
v1_proj = P * v1
v2_proj = P * v2

See Also

• project_onto_eigenspace: Project single vector without forming matrix

compute_schur_projector(Q::AbstractMatrix, indices::UnitRange{Int})

Compute the projection matrix onto the invariant subspace spanned by Schur vectors.

Returns the n × n projection matrix PS = QS QS^† where QS = Q[:, indices].

Arguments

• Q::AbstractMatrix: Unitary Schur matrix (n × n)
• indices::UnitRange{Int}: Column indices of Schur vectors

Returns

n × n projection matrix PS with PS² = P_S.

Example

F = schur(A)
P = compute_schur_projector(F.Z, 1:3)  # Project onto first 3-dimensional invariant subspace

qi_intervals(A::BallMatrix)

Qi (1984, Theorem 2) intervals for singular values. Returns the intervals Bᵢ as a Vector{Ball}. See Ref. [13].

qi_sqrt_intervals(A::BallMatrix)

Sharper square-root intervals for the singular values (Qi 1984, Theorem 3). See Ref. [13].