Certified Matrix Decompositions

BallArithmetic.jl provides rigorous certification for several matrix decompositions. Each method computes both a numerical approximation and a rigorous error bound.

Overview

Rigorous vs Non-Rigorous Results

Each decomposition method returns:

• Non-rigorous residual: The computed ‖A - UΣVᴴ‖_F using standard floating-point
• Rigorous bound: A mathematically guaranteed upper bound on the true error

The rigorous bound accounts for:

• Rounding errors in all computations
• Potential instabilities in the algorithm
• Accumulated error from multiple operations

For certification purposes, always use the rigorous bound.

Other Matrix Decompositions (100×100 Matrix)

Beyond SVD, BallArithmetic provides verified versions of standard matrix decompositions.

Standard (Float64 → BigFloat)

These use Float64 computation with BigFloat certification:

GLA-Based (Native BigFloat)

With GenericLinearAlgebra, we achieve 60+ orders of magnitude better residuals:

GLA functions:

using BallArithmetic, GenericLinearAlgebra

result_lu = verified_lu(A)              # in BigFloat, the default
result_qr = verified_qr(A)
result_chol = verified_cholesky(A_spd)
result_svd = verified_svd(A)            # a VerifySvdAllResult
result_fast = verified_lu(Float64, A)
result_d64 = verified_lu(A; oracle = Double64)   # starting factorisation in Double64

Usage examples:

using BallArithmetic
using LinearAlgebra

A = randn(100, 100) + 5I

# LU decomposition
result_lu = verified_lu(A)
# result_lu.L, result_lu.U are BallMatrix enclosures
# result_lu.residual_norm is the rigorous error bound

# QR decomposition
result_qr = verified_qr(A)
# result_qr.Q, result_qr.R are BallMatrix enclosures
# result_qr.orthogonality_defect bounds ‖QᴴQ - I‖

# Cholesky (for symmetric positive definite)
A_spd = A'A + 100I
result_chol = verified_cholesky(A_spd)
# result_chol.G encloses the upper triangular factor, A_spd = G'G

# Polar decomposition A = QP
result_polar = verified_polar(A)
# result_polar.Q (unitary), result_polar.P (Hermitian positive definite)

# Takagi factorization (for complex symmetric)
A_symm = A + transpose(A)
result_takagi = verified_takagi(A_symm)
# A = U Σ Uᵀ (not Uᴴ!)

Schur Complement Bounds (Oishi 2023 / Rump-Oishi 2024)

For block-structured matrices, the Schur complement method provides efficient σ_min bounds:

result = rump_oishi_2024_sigma_min_bound(G; block_size=m)

This implements:

• Oishi 2023: Base Schur complement algorithm
• Rump-Oishi 2024: Improved ψ(N) formula that works for ‖N‖ ≥ 1

Key advantages:

• Works on matrices too large for full SVD
• Exploits block diagonal dominance structure
• Automatic optimal block size selection

References

• Ogita, T. & Aishima, K. (2020), "Iterative refinement for singular value decomposition based on matrix multiplication"
• Rump, S.M. & Ogita, S. (2024), "A Note on Oishi's Lower Bound for the Smallest Singular Value of Linearized Galerkin Equations"
• Oishi, S. (2023), "Lower bounds for the smallest singular values of generalized asymptotic diagonal dominant matrices"