API - Matrix Decompositions

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

Inclusions of the factors of the LU decomposition of the m × n matrix A, by Section 3 of Rump and Ogita (2024); the method is described at _rumpogita2024_lu. Returns a VerifiedLUResult: when certified is true, A[p, q] = LU with L of size m × min(m, n) unit lower triangular and U of size min(m, n) × n upper triangular, both in the ball matrices returned. p comes from partial pivoting; q is the identity for m ≥ n, and for m < n the column permutation of Section 3.4 of the paper.

Arguments

• T: the arithmetic of the enclosure, BigFloat (the default) or Float64, as a leading argument. In BigFloat the computation runs at precision_bits and the factors are enclosed to about that precision; in Float64 the radii are of the order of the rounding unit times the size of the factors. Float64 with a BigFloat input raises.
• precision_bits: the BigFloat precision. A BigFloat input with more digits is rounded to it and the rounding is enclosed.
• oracle: a floating-point type (Double64, a MultiFloat) in which the approximate factorisation the method starts from is computed, in place of the working arithmetic; nothing, the default, means the working arithmetic. It has no part in the verification: the factors are enclosed from residuals in the working arithmetic whatever their accuracy, which only decides how small the radii come out. The same keyword is taken by verified_cholesky, verified_qr, verified_polar and verified_takagi.

Example

A = randn(100, 100)
r = verified_lu(Float64, A)
r.certified      # A[r.p, r.q] = L U with L in r.L and U in r.U

Reference

Section 3 of [15].

VerifiedLUResult{LM, UM, RT}

Result of verified_lu.

Fields

• L::LM: ball matrix containing the unit lower triangular (trapezoidal) factor
• U::UM: ball matrix containing the upper triangular (trapezoidal) factor
• p::Vector{Int}: row permutation
• certified::Bool: whether the verification succeeded; when false the radii are infinite and nothing is proved
• residual_norm::RT: an upper bound of ‖L̃Ũ − A[p, q]‖_∞ / ‖A‖_∞ over all L̃ ∈ L, Ũ ∈ U
• q::Vector{Int}: column permutation, the identity unless the matrix has more columns than rows

Statement

When certified is true, A[p, q] has a unique decomposition LU with L unit lower and U upper triangular, and these factors lie in the ball matrices L and U.

Reference

Section 3 of [15].

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.

_rumpogita2024_alg3_2(E::BallMatrix) -> (; LE, UE, LinvE, UinvE) or nothing

Algorithm 3.2 (LU_E) of Rump and Ogita (2024): inclusions of the factors of the LU decomposition, without pivoting, of a perturbed identity matrix, and of their inverses.

For E of size m × n, p = min(m, n) and every Ẽ of the ball E, the matrix I_{m,n} + Ẽ has a unique decomposition LU with L of size m × p unit lower triangular (trapezoidal) and U of size p × n upper triangular (trapezoidal), and

L ∈ I_{m,p} + LE,      U ∈ I_{p,n} + UE,

LE strictly lower and UE upper triangular. LinvE is returned when L is square (m ≤ n), with L⁻¹ ∈ I + LinvE, and UinvE when U is square (m ≥ n), with U⁻¹ ∈ I + UinvE; otherwise the field is nothing. The cost is O(mn) operations.

The bounds are (3.1)-(3.7) of the paper. With Δ_L := sum(|E^{[ℓ]}|, 2) max(|E|) / (1 − ‖E‖_∞), strictly lower part, |L − I − E^{[ℓ]}| ≤ Δ_L (3.1); with G_L = |E^{[ℓ]}| + Δ_L, |L⁻¹ − I + E^{[ℓ]}| ≤ Δ_L + sum(G_L, 2) max(G_L)/(1 − ‖G_L‖_∞) (3.4); with B equal to |E| with its strictly lower part replaced by Δ_L, |U − I − E^{[u]}| ≤ Δ_U := sum(G_L, 2) max(B) / (1 − ‖G_L‖_∞), upper part (3.5); and with G_U = |E^{[u]}| + Δ_U, |U⁻¹ − I + E^{[u]}| ≤ Δ_U + sum(G_U, 2) max(G_U)/(1 − ‖G_U‖_∞) (3.7). Here sum(·, 2) is the column of row sums, max(·) the row of column maxima, and the products are outer products.

Returns nothing when one of the three norms ‖E‖_∞, ‖G_L‖_∞, ‖G_U‖_∞ is not below one, which the paper's listing assumes. Every bound is rounded up.

Reference

Section 3.1 and Algorithm 3.2 of [15].

verified_qr([T = BigFloat,] A::AbstractMatrix; compute_full_Q = false, precision_bits = 256)

Inclusions of the factors of the QR decomposition of the m × n matrix A, real or complex, by Section 5 of Rump and Ogita (2024); the method is described at _rumpogita2024_qr. Returns a VerifiedQRResult: when certified is true, A has full rank and A = QR with R upper triangular with positive diagonal entries.

For m ≥ n the decomposition is the economy-size one, Q of size m × n and R of size n × n, which is unique. With compute_full_Q = true, Q is m × m unitary and R is m × n; the last m − n columns of Q are an orthonormal basis of the orthogonal complement of the range, which is not unique, and the statement is that some such basis lies in the ball. For m < n, Q is m × m, from the leading square block, and R = Q*A is m × n.

Arguments

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

Example

A = randn(100, 50)
r = verified_qr(Float64, A)
r.certified      # A = Q R with Q (100 × 50) in r.Q and R (50 × 50) in r.R

Reference

Section 5 and Lemma 5.1 of [15].

VerifiedQRResult{QM, RM, RT}

Result of verified_qr.

Fields

• Q::QM: ball matrix containing the factor with orthonormal columns
• R::RM: ball matrix containing the upper triangular factor
• certified::Bool: whether the verification succeeded; when false the radii are infinite and nothing is proved
• residual_norm::RT: an upper bound of ‖Q̃R̃ − A‖_∞ / ‖A‖_∞ over all Q̃ ∈ Q, R̃ ∈ R
• orthogonality_defect::RT: an upper bound of ‖Q̃*Q̃ − I‖₂ over all Q̃ ∈ Q; the factor itself has orthonormal columns, and this measures the width of the enclosure

Statement

When certified is true, A has full rank and A = QR with Q having orthonormal columns and R upper triangular with positive diagonal entries, Q and R in the ball matrices returned.

Reference

Section 5 of [15].

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.

_rumpogita2024_qr(A::BallMatrix; full = false) -> (; Q, R) or nothing

Section 5 of Rump and Ogita (2024): inclusions of the factors of the QR decomposition. For every matrix Ã of the m × n ball A, Ã is proved to have full rank and à = QR with R upper triangular with positive diagonal entries, Q and R in the ball matrices returned. Returns nothing when the verification fails.

For m ≥ n and full = false the decomposition is the economy-size one, Q of size m × n with orthonormal columns and R of size n × n, which is unique. With full = true, Q is m × m unitary and R is m × n; its last m − n columns are an orthonormal basis of the orthogonal complement of the range, which is not unique, and the statement is that some such basis lies in the ball. For m < n the decomposition is the full one of the leading square block, Q of size m × m, and R = Q*Ã of size m × n.

The method

Q̃R̃ is a floating-point decomposition of the midpoint, R̃ made to have a real positive diagonal, and X_R ≈ R̃⁻¹. C := Ã X_R, kept as an unevaluated sum of two terms, is close to a matrix with orthonormal columns, and C*C = I + E is enclosed in two-fold precision without forming Ã*Ã. The Cholesky factor of I + E is G_E = D^{1/2} L_E* by Algorithm 3.2 as in Section 4, and

R = G_E X_R⁻¹,        Q₁ = Ã X_R G_E⁻¹,

the second being the possibility the paper finds better for Q₁. R is the Cholesky factor of Ã*Ã, upper triangular with positive diagonal since those of G_E and of X_R are.

For the orthogonal complement, with Q̃₂ its floating-point approximation, the square system (5.1), [Ã*; N] Y = [0; αI] with N = fl(α Q̃₂*), has a solution whose columns span the orthogonal complement of the range of Ã. An inclusion Q̃₂ + Δ of it comes from a verified solve for the correction, with the residual in two-fold precision, and Lemma 5.1 gives an orthonormal basis Q₂ of that complement with

‖Q₂ − Q̃₂‖₂ ≤ ‖I − Q̃₂*Q̃₂‖₂ + √2 ‖Δ‖₂,

which bounds every entry of Q₂ − Q̃₂.

Reference

Section 5 and Lemma 5.1 of [15].

cholesky_qr2(X) -> (Z, U)

The CholeskyQR2 algorithm of Yamamoto, Nakatsukasa, Yanagisawa and Fukaya (2015) for the thin QR decomposition X ≈ Z U of a matrix with at least as many rows as columns (Section 2.2):

A = X*X,   R = chol(A),   Y = X R⁻¹,
B = Y*Y,   S = chol(B),   Z = Y S⁻¹,   U = S R.

Under their assumption (2.1), 8 κ₂(X) √(mnu + n(n+1)u) ≤ 1, so κ₂(X) at most of order u^(-1/2), the paper proves that the loss of orthogonality of Z and the relative residual are of order u. For a larger condition number the Cholesky factorisation can break down, and the exception of cholesky is then thrown.

The result is an approximation with no error bound; see verified_qr.

Reference

The roundoff analysis is [18], whose standing assumption (2.1) is the inequality above.

verified_polar([T = BigFloat,] A::AbstractMatrix; right = true, kappa = 0, precision_bits = 256)

Inclusions of the factors of the polar decomposition of the m × n matrix A, by Section 7 of Rump and Ogita (2024), from the inclusions of the singular value decomposition of Rump and Lange (2023); the method is described at _rumpogita2024_polar. Returns a VerifiedPolarResult.

With right = true, which needs m ≥ n, the decomposition is A = QP with Q of size m × n with orthonormal columns and P of size n × n Hermitian positive definite. With right = false, which needs m ≤ n, it is A = PQ with P of size m × m and Q of size m × n with orthonormal rows, obtained from the right decomposition of A*. certified is false when a singular value is not separated from zero.

Arguments

• kappa: the clustering threshold of verifysvdall.
• T: the arithmetic of the enclosure, BigFloat (the default), at precision_bits, or Float64, as a leading argument. Float64 with a BigFloat input raises.
• precision_bits: the BigFloat precision. A BigFloat input with more digits is rounded to it and the rounding is enclosed.

Example

A = randn(ComplexF64, 50, 50)
r = verified_polar(Float64, A)
r.certified      # A = Q P with Q in r.Q and P in r.P

Reference

Section 7 of [15], over the singular value enclosures of [19].

VerifiedPolarResult{QM, PM, RT}

Result of verified_polar.

Fields

• Q::QM: ball matrix containing the factor with orthonormal columns (rows, for the left decomposition)
• P::PM: ball matrix containing the Hermitian positive definite factor
• is_right::Bool: true for A = QP, false for A = PQ
• certified::Bool: whether the verification succeeded; when false the radii are infinite and nothing is proved
• residual_norm::RT: an upper bound of ‖Q̃P̃ − A‖_∞ / ‖A‖_∞ (of ‖P̃Q̃ − A‖_∞ / ‖A‖_∞ for the left decomposition) over all Q̃ ∈ Q, P̃ ∈ P

Statement

When certified is true, A has full rank and A = QP (or A = PQ) with P Hermitian positive definite and Q with orthonormal columns (rows), Q and P in the ball matrices returned.

Reference

Section 7 of [15].

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.

_rumpogita2024_polar(A::BallMatrix; kappa = 0) -> (; Q, P, svd) or nothing

Section 7 of Rump and Ogita (2024): inclusions of the factors of the polar decomposition. For every matrix Ã of the m × n ball A, m ≥ n, Ã is proved to have full rank and à = QP with Q of size m × n with orthonormal columns and P Hermitian positive definite, Q and P in the ball matrices returned. Returns nothing when a singular value is not separated from zero. svd is the result of verifysvdall the inclusions come from, and kappa its threshold.

The method

The paper obtains the factors from the singular value decomposition à = UΣV* as Q = UV* and P = VΣV*, with the inclusions of Rump and Lange (2023). Those are, for each cluster μ of singular values, inclusions of orthonormal bases of the left and of the right singular subspace, and the two bases of a cluster are not matched to each other (for one singular value, the two vectors are determined up to a factor of modulus one each). So the factors are written with the right bases alone: with V_μ any orthonormal basis of the right singular subspace of the cluster and H_μ := V_μ*Ã*ÃV_μ, whose eigenvalues are the squares of the singular values of the cluster,

P = Σ_μ V_μ H_μ^{1/2} V_μ*,        Q = Ã Σ_μ V_μ H_μ^{-1/2} V_μ*.

The singular values of the cluster lie in an interval [a, b] with a > 0, so H_μ^{1/2} is within (b − a)/2 of ((a + b)/2) I in the spectral norm, hence entrywise, and H_μ^{-1/2} within half the width of [1/b, 1/a] of its midpoint times I. For a cluster of one singular value these are P = Σ σ_j v_j v_j* and Q = Σ Ãv_j v_j*/σ_j.

Reference

Section 7 of [15], over the singular value enclosures of [19]; see verifysvdall.

PolarIterationResult

Result of polar_newton, polar_newton_schulz and polar_qdwh: the computed factors of A ≈ U H, U with orthonormal columns and H Hermitian, the number of iterations, and whether the stopping test of the algorithm was met. The factors are floating-point approximations.

polar_newton(A; maxiter = 100, tol = nothing) -> PolarIterationResult

Algorithm Polar of Higham (1986) for a nonsingular square matrix A, giving A ≈ U H:

1. X₀ = A;
2. repeat Y_k = X_k⁻¹; γ_k = 1 when close to convergence, otherwise γ_k = √(β_k/α_k) with α_k = √(‖X_k‖₁‖X_k‖_∞) and β_k = √(‖Y_k‖₁‖Y_k‖_∞); X_{k+1} = (γ_k X_k + Y_k*/γ_k)/2; until converged;
3. U = X_{k+1}, H = (U*A + (U*A)*)/2.

The convergence test is (3.17), ‖X_{k+1} − X_k‖₁ ≤ δ_n ‖X_k‖₁, with δ_n = tol; the paper asks for a small multiple of the unit roundoff depending on n, and the default is 8n times it. "Close to convergence" is taken as a relative change below 10⁻², the value at which the paper reverts to the unaccelerated iteration. A matrix with more rows than columns is reduced to its triangular factor first (Section 3.1 of the paper).

The result is an approximation with no error bound; see verified_polar.

Reference

The iteration, its acceleration parameters and a way of estimating them are [20], which also proves the convergence quadratic. Theorem 1.1 there is the existence statement: for an m × n matrix with m ≥ n there are a matrix with orthonormal columns and a unique Hermitian positive semidefinite factor whose product is it, and the factor is positive definite and the first matrix unique when the rank is full.

polar_newton_schulz(A; maxiter = 100, tol = nothing) -> PolarIterationResult

The Newton-Schulz iteration (8.20) of Higham (2008),

X_{k+1} = X_k (3I − X_k* X_k)/2,      X₀ = A,

which converges quadratically to the unitary polar factor when every singular value of A lies in (0, √3), and applies to rectangular A of full column rank. Throws when the computed ‖A‖₂ is not below √3. The iteration stops when the relative change of the iterate in the 1-norm is at most tol (default 8n times the unit roundoff), the test of polar_newton; the book states no test for this iteration.

The result is an approximation with no error bound; see verified_polar.

Reference

Chapter 8 of [21], on the polar decomposition.

polar_qdwh(A; maxiter = 20) -> PolarIterationResult

The QDWH algorithm of Nakatsukasa, Bai and Gygi (2010), the QR-based dynamically weighted Halley iteration, for a nonsingular square A:

1. X₀ = A/α, ℓ₀ = β/α;
2. repeat: a_k = h(ℓ_k), b_k = (a_k − 1)²/4, c_k = a_k + b_k − 1, with h of (3.15); the QR decomposition [√c_k X_k; I] = [Q₁; Q₂] R; X_{k+1} = (b_k/c_k) X_k + (1/√c_k)(a_k − b_k/c_k) Q₁ Q₂*; ℓ_{k+1} = ℓ_k (a_k + b_k ℓ_k²)/(1 + c_k ℓ_k²); until convergence;
3. U = X_k.

The estimates are those of Section 4 of the paper: α = ‖A‖_F, and β = 1/(γ√n) (4.7) with γ the LAPACK estimate of ‖A⁻¹‖₁ (for element types outside the BLAS, the 1-norm of the computed inverse). The stopping test is (3.21), ‖X_k − X_{k−1}‖_F ≤ (4ε)^{1/3}, ε the machine epsilon. The paper bounds the number of iterations by six for condition numbers up to 10¹⁶ in double precision. A matrix with more rows than columns is reduced to its triangular factor first, as in polar_newton.

The result is an approximation with no error bound; see verified_polar.

Reference

The dynamically weighted Halley iteration and its QR-based variant are [22], which chooses the weights by a max-min problem from α = ‖A‖₂ and a lower bound of β = σ_min(A), and proves the convergence global and asymptotically cubic, with at most six iterations reaching 10⁻¹⁶ in exact arithmetic when κ₂(A) ≤ 10¹⁶.

verified_takagi([T = BigFloat,] A::AbstractMatrix{<:Complex}; precision_bits = 256)

Inclusions of the factors of the Takagi decomposition A = UΣUᵀ of a complex symmetric matrix, by the third method of Section 9 of Rump and Ogita (2024), the one through the real symmetric matrix [E F; F −E] for A = E + iF; the method is described at _rumpogita2024_takagi. Returns a VerifiedTakagiResult.

certified is false when the singular values are not proved positive and simple, which includes every matrix with a multiple singular value, the identity for one: the decomposition is not unique there and the method does not apply.

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

Arguments

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

Example

B = randn(ComplexF64, 20, 20); A = B + transpose(B)
r = verified_takagi(Float64, A)
r.certified      # A = U Σ Uᵀ with U in r.U and the diagonal of Σ in r.Σ

Reference

Section 9 of [15].

VerifiedTakagiResult{UM, SV, RT}

Result of verified_takagi.

Fields

• U::UM: ball matrix containing the unitary factor
• Σ::SV: balls containing the diagonal of Σ, the singular values of A, in decreasing order
• certified::Bool: whether the verification succeeded; when false the radii are infinite and nothing is proved
• residual_norm::RT: an upper bound of ‖ŨΣ̃Ũᵀ − A‖_∞ / ‖A‖_∞ over all Ũ ∈ U, Σ̃ ∈ Σ

Statement

When certified is true, the symmetric part (A + Aᵀ)/2 of the input, which is the input when it is symmetric (Aᵀ = A, not A* = A), is nonsingular with simple singular values and equals UΣUᵀ with U unitary and Σ diagonal with positive entries, U and the diagonal of Σ in the balls returned. U is determined up to the sign of each column.

Reference

Section 9 of [15].

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.

_rumpogita2024_takagi(A::BallMatrix) -> (; U, sigma) or nothing

Section 9 of Rump and Ogita (2024): inclusions of the factors of the Takagi decomposition of a complex symmetric matrix. For every symmetric matrix à = Ãᵀ of the ball A, Ã is proved nonsingular with simple singular values and à = UΣUᵀ with U unitary and Σ diagonal with positive entries in decreasing order, U in the ball matrix U and the diagonal of Σ in the balls sigma. U is determined up to the sign of each column. Returns nothing when the singular values are not proved positive and simple; the decomposition is discontinuous at a singular matrix.

The method

The third method of the paper, which it finds best. With à = E + iF, E and F real symmetric, the real symmetric matrix M = [E F; F −E] has the eigenvalues ±σ_j; if [x; y] is a unit eigenvector for σ_j > 0, then u = x + iy is a unit vector with à ū = σ_j u, and these u are the columns of U. The eigenvalues and eigenvectors of M are enclosed by _rumplange2023_eig; each positive eigenvalue is required to be a cluster of its own, so that its unit eigenvector is determined up to sign. The bound on the eigenvector is a bound of the norm of the error of [x; y], hence of the modulus of each entry of x + iy.

Reference

Section 9 of [15].

_rumpogita2024_schur(A::BallMatrix; kwargs...) -> (; Q, T, X, D) or nothing

Section 8 of Rump and Ogita (2024): inclusions of the factors of a complex Schur decomposition of a matrix with simple eigenvalues. For every matrix Ã of the ball A there are a unitary Q and an upper triangular T with à = QTQ*, Q and T in the ball matrices returned; the diagonal of T holds the eigenvalues in the order of D. Returns nothing when the verification fails, which includes every matrix for which the eigenvalues are not proved simple: the paper restricts the method to that case, the decomposition being discontinuous at a multiple eigenvalue.

The method

(8.1) of the paper: with à = XDX⁻¹ an eigendecomposition and X = QR the QR decomposition of X, à = QTQ* with T = RDR⁻¹. Inclusions of X and of the diagonal D come from verifyeigall (Section 6 refers to Rump (2022) for a general matrix), each cluster being required to be a single certified eigenvalue; kwargs are passed to it. The inclusions of Q and R are those of _rumpogita2024_qr for the ball matrix X, valid for the true eigenvector matrix in it; and T is the solution of the linear system TR = RD with R and D replaced by their inclusions, enclosed by a verified solve. Its diagonal is D and its entries below the diagonal are exact zeros. The paper notes that the replacement of X, R and D by inclusions is a source of overestimation that the other decompositions do not have.

Reference

Section 8 of [15].

verified_svd([T = BigFloat,] A::AbstractMatrix; precision_bits = 256, kappa = 0)
    -> VerifySvdAllResult

verifysvdall for a matrix with exact entries, in the arithmetic T: enclosures of all the singular values and of bases of the singular subspaces, by the method of Rump and Lange (2023). T is BigFloat, at precision_bits, or Float64; Float64 with a BigFloat input raises. kappa is the clustering threshold of verifysvdall. The arguments are those of verified_lu and the other verified_* decompositions.

Reference

The enclosures are those of [19]; see verifysvdall.

SVDMethod

Which enclosure of the singular values a certification uses, as a type, so that the choice is made by dispatch and the routine that implements a theorem is named after it: MiyajimaM1 is Theorem 7 of [23], MiyajimaM3 Theorem 10, MiyajimaM4 Theorem 11, and MiyajimaAuto is not a theorem but the rule that picks between the first two.

rigorous_svd, certify_svd, rigorous_svd_gpu and adaptive_ogita_svd take one of these as their method keyword. svdbox, which returns the values and nothing else, selects by a Symbol instead and reaches three further enclosures that these types do not name.

MiyajimaAuto <: SVDMethod

The default of rigorous_svd: MiyajimaM3, Theorem 10, for an exact input, and MiyajimaM1, Theorem 7, for one carrying a radius.

This is a selection rule of this package and not a result of [23]; the measurements behind it are in _svd_auto_theorem. Pass method explicitly to override it.

On an input carrying a radius this is not what svdbox does. svdbox returns values only, so it can certify mid(A) by Theorem 10 and widen each interval by ‖rad(A)‖₂ (_miyajima2014_thm10_weyl), which is tighter; rigorous_svd also returns singular-vector bounds and the residuals E, F, G measured against the whole ball, which a midpoint-only certification does not produce, so it stays on Theorem 7. For the values alone, prefer svdbox.

MiyajimaM1 <: SVDMethod

Theorem 7 of [23], the economy enclosure the paper itself calls svdbox and its numerical section labels M1.

With Û, V̂ the first q = min(m, n) columns of the approximate frames, Σ̂ the approximate singular values, and the residuals Ê = ÛΣ̂V̂ᵀ − A, F̂ = V̂ᵀV̂ − I, Ĝ = ÛᵀÛ − I measured against A itself, the theorem gives, when ‖F̂‖₂ < 1 and ‖Ĝ‖₂ < 1,

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

the singular values in decreasing order on both sides, at a cost of 24Qq² + 12q³ operations with Q = max(m, n). Theorem 8 of the paper proves these bounds never looser than those of its Theorem 4, the square-frame enclosure of Oishi (_miyajima2014_thm4).

MiyajimaM3, Theorem 10, uses the same economy frames and is cheaper and at least as tight, so rigorous_svd prefers it on an exact input and keeps this one on a ball, where Theorem 10's residual would propagate the input radius through a product.

MiyajimaM3 <: SVDMethod

Theorem 10 of [23], the algorithm his numerical section labels M3, and what rigorous_svd uses on an exact input.

The economy frames of Theorem 7 with a one-sided residual: with Λᵢᵢ and Λ̄ᵢᵢ the scalings of Σ̂ᵢᵢ by √((1∓‖Ĝ‖)/(1±‖F̂‖)) and ρ = ‖AV̂ − ÛΣ̂‖₂/√(1−‖F̂‖₂) for m ≥ n, with F̂ and Ĝ exchanged otherwise,

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

Cheaper than Theorem 7, 22Qq² + 12q³ against 24Qq² + 12q³, since the residual is m × q rather than m × n, and at least as tight in the paper's Tables 2 and 3.

MiyajimaM4 <: SVDMethod

Theorem 11 of [23], the algorithm his numerical section labels M4.

One frame instead of two: V̂ approximates the eigenvectors of AᵀA (of AAᵀ when m < n), whose eigenvalues are the squares of the singular values, and the theorem isolates those eigenvalues in the frame (AV̂)ᵀ(AV̂) = D̂ + Ê with D̂ diagonal. The statement, its hypotheses and the branch it takes on a cluster of indices are at _miyajima2014_thm11. It says nothing about the left singular vectors and nothing about a residual ÛΣ̂V̂ᵀ − A.

This method does not go through certify_svd, whose two frames and residual Ê the theorem does not use: rigorous_svd(A; method = MiyajimaM4()) is routed to rigorous_svd_m4, certify_svd with this method raises an ArgumentError, and svdbox(A; method = :miyajima2014_thm11) gives the values alone.

certify_svd(A::BallMatrix, svdA, method::SVDMethod; apply_vbd = true) -> RigorousSVDResult

Enclosures of the singular values of every matrix in A from an approximate singular value decomposition svdA that the caller supplies: an SVD object or a named tuple (U, S, V, Vt), in any precision. This is the certification step of rigorous_svd, which computes the approximate decomposition itself; it is called directly where the approximation comes from elsewhere, a refinement in higher precision or the decomposition of a nearby matrix.

Returns a RigorousSVDResult. The residual U Σ Vᵀ − A and the defects VᵀV − I, UᵀU − I are enclosed in ball arithmetic and the bounds of method, one of the types listed at SVDMethod, are applied to upper bounds of their spectral norms; the enclosures returned hold for every matrix of A. A method of MiyajimaAuto is resolved first, to Theorem 10 on an exact input and Theorem 7 on one carrying a radius, and MiyajimaM4 raises an ArgumentError, Theorem 11 not being stated for these frames.

Where one of the two defects is not below 1, which every one of these theorems assumes, a warning is issued and the result proves nothing: the singular values carry infinite radius, residual_norm is Inf, and block_diagonalisation is nothing. With apply_vbd = true, the default, that last field holds the refinement schur_gershgorin_enclosure produces from Σ'Σ, the diagonal matrix of the squared enclosures; it is nothing when the flag is false and nothing in BigFloat whatever the flag.

svd_bounds(A::BallMatrix; method = :auto) -> (lo, hi)

Certified bounds lo[i] ≤ σᵢ(X) ≤ hi[i] for every X in the ball A and every i = 1:min(size(A)...), singular values in decreasing order, read off svdbox with the given method and rounded outward. Where the enclosure fails, or the input carries a non-finite entry, the bounds are the trivial lo = 0, hi = Inf, so the result is always a valid bound.

A caller whose numbers must not change with a change of default names the theorem, for instance method = :miyajima2014_thm10_weyl, rather than relying on :auto.

svd_lower_bound_sigma_min(A::BallMatrix; method = :auto) -> T

Certified lower bound on the smallest singular value σ_q, q = min(size(A)...), of every matrix in the ball, from svd_bounds; zero when the enclosure fails.

rigorous_svd_m4(A::BallMatrix; apply_vbd = true) -> RigorousSVDResult

The singular values of A enclosed by Theorem 11 of [23], the algorithm his numerical section labels M4, as a RigorousSVDResult; the theorem and its hypotheses are stated at _miyajima2014_thm11, which this routine calls.

The theorem uses one frame, the eigenvectors of A'A, and says nothing about left singular vectors or about a residual UΣV' − A. The fields U and residual of the result are therefore nothing, and residual_norm and left_orthogonality_defect are Inf. V is the approximate frame as a point matrix, with right_orthogonality_defect an upper bound of ‖V'V − I‖₂.

Reference

Theorem 11 of [23].

rigorous_svd_gpu(A::BallMatrix{Float64,Float64}; method = MiyajimaM1(),
                 apply_vbd = true, seed_on = :cpu)

GPU-accelerated rigorous SVD of a (CPU-resident) Float64 ball matrix, returning the same RigorousSVDResult as rigorous_svd.

The certification math is identical to rigorous_svd — only the O(n³) certification products (UΣVᵀ, VᵀV, UᵀU, and the residual widening term) are evaluated on the GPU through the rigorous INT8-Ozaki MMul4 dispatch, while the O(n²) rigorous norm bounds and singular-value enclosures are finished on the CPU with directed rounding. The seed factorisation defaults to CPU LAPACK (seed_on = :cpu, the hybrid configuration, which benchmarks fastest on consumer/A40 cards because their Float64 throughput is throttled); pass seed_on = :gpu to seed with cuSOLVER instead.

This method is only available when the CUDA extension is loaded (using CUDA) and a functional GPU is present. The enclosures it returns overlap those of rigorous_svd; the GPU midpoint truncation makes the singular-value radii modestly looser (still rigorous, typically ~10⁻¹⁰).

Reference

The enclosure is the one method names, by default Theorem 7 of [23]; see SVDMethod for the others.

svd_bigfloat(A::AbstractMatrix{T}) where T<:Union{BigFloat, Complex{BigFloat}}

A floating-point singular value decomposition of a BigFloat matrix, by the svd of GenericLinearAlgebra.jl, at the current precision. An approximation with no error bound: it is the starting point that rigorous_svd certifies for a BigFloat ball matrix. The method is defined in the GenericLinearAlgebra extension.

_rigorous_svd_bigfloat(A, method; apply_vbd = true, use_cache = true)
    -> RigorousSVDResult

The BigFloat path of rigorous_svd: an approximate decomposition of mid(A) is computed at the current BigFloat precision by svd_bigfloat, and certify_svd encloses the singular values of every matrix of A from it with the bounds of method. apply_vbd is passed on, and has no effect here, the block refinement being skipped in BigFloat.

With use_cache = true, the default, the approximate factors are kept in the module-level cache (set_svd_cache!, svd_cache_stats, clear_svd_cache!) and reused when the next call has a midpoint of the same hash and factors of matching size. The cache holds an approximation and takes no part in what is proved: a stale or wrong entry can only widen the enclosure, since the residuals are measured against A.

The refinement path of ogita_svd_refine and adaptive_ogita_svd is not used here and remains available on its own.

_miyajima2014_thm10_bounds(S, normR, normF, normG, tall, ::Type{T}) -> (lo, hi)

The bounds of Theorem 10 from the quantities it needs: the economy singular values S, the one-sided residual norm normR (‖AV̂ − ÛΣ̂‖₂ when m ≥ n, else ‖ÛᵀA − Σ̂V̂ᵀ‖₂), the two orthogonality defects, and tall = m ≥ n. Held separately so that _miyajima2014_thm10, which builds its own frames, and rigorous_svd, which already has them, share one copy of the formula.

ogita_iterations_for_precision(target_bits::Int; start_bits::Int = 53) -> Int

How many steps of Algorithm 1 (RefSVD) of [6] ogita_svd_refine runs when it is asked to carry an approximate decomposition from start_bits correct bits, the 53 of Float64 by default, to target_bits: ceil(log2(target_bits / start_bits)), and zero when target_bits ≤ start_bits.

The count comes from the quadratic convergence of that algorithm, each step roughly doubling the number of correct bits; the paper's Theorem 1 gives the quadratic rate under a closeness condition on the initial guess relative to the gaps between the singular values. It is a count of steps to run and not a bound on anything: this function checks no hypothesis and the accuracy actually reached is not certified by it. What is certified is the enclosure certify_svd returns from the refined factors, whatever their accuracy.

From 53 bits: 3 steps for 256, 4 for 512, 5 for 1024.

_ogita_svd_refine_impl(A, U, Σ, V, max_iterations, check_convergence)
    -> OgitaSVDRefinementResult

Algorithm 1 (RefSVD) of [6] applied to the approximate decomposition A ≈ UΣVᵀ. The exact factors are written Û(I + F) and V̂(I + G); the orthogonality defects I − ÛᵀÛ and I − V̂ᵀV̂ together with the projected matrix ÛᵀAV̂ give linear equations for the corrections, and the factors are updated with their solution. The step is a Newton step, so it converges quadratically from a modestly accurate start, and all of its cost sits in matrix products, which is what lets the refinement be run at a higher precision than the certification.

The paper states the algorithm for a real m × n matrix with m ≥ n, simple singular values and the full m × m left factor. A matrix with fewer rows than columns is sent to the same routine as A* with the two frames exchanged, and a thin left factor is completed by an approximately orthonormal basis of the orthogonal complement before the first step and cut back afterwards. max_iterations steps are run, or fewer when check_convergence is set and the iterate stops moving.

Returns floating-point factors with no error bound, and a residual_norm which is the estimate of _spectral_norm_bound and not a bound either. An enclosure is obtained by handing the refined factors to certify_svd, which is what adaptive_ogita_svd does at each precision level.

_spectral_norm_bound(A::AbstractMatrix)

The floating-point value of √(‖A‖∞ ‖A‖₁), with ‖A‖∞ the largest row sum of |A| and ‖A‖₁ the largest column sum. In exact arithmetic that quantity is an upper bound of ‖A‖₂; here the sums, the product and the root are rounded to nearest and no rounding is accounted for, so the value returned is an estimate and not a bound. It is used for the stopping test of the refinement and for the diagnostic residual_norm of its result, and for nothing that is stated. An upper bound of the spectral norm is upper_bound_L2_opnorm on a ball matrix.

block_enclosure(vbd) -> Vector{NamedTuple}

Block-disc enclosure of σ(A) from a VBD result (SchurGershgorinResult or Miyajima2014aSchurNewtonResult): one disc per cluster, (center = cᵢ, radius = nᵢ + rᵢ, mult = |Cᵢ|) with the barycentre cᵢ (block_centers), within-block non-normality nᵢ (block_nonnormality) and localized coupling rᵢ (block_coupling).

From the barycentre floor σ_min(Pᵢ − zI) ≥ |z − cᵢ| − nᵢ and the block Gershgorin region {z : σ_min(Pᵢ − zI) ≤ rᵢ}, every eigenvalue of A lies in one of the discs. mult is the count that Theorem 4 of [24] gives for a union of discs disjoint from the others; see schur_gershgorin_enclosure for the reference and for what of that count is not written out for these radii.

overlap_components(balls) -> Vector{Vector{Int}}

The connected components of the graph on balls (real intervals or complex discs) in which two balls are joined when they may meet, as sorted vectors of indices. Two balls are joined unless they are shown disjoint: a lower bound of the distance of the centres (dist_down) exceeding an upper bound of the sum of the radii, or intersect_ball returning nothing for real intervals. So balls in different components are disjoint, and the union of the balls of one component is disjoint from the others.

sep_clusters(ints, compA, compB)

A lower bound of the distance between the union of the balls ints[compA] and the union of the balls ints[compB], zero where it cannot be shown positive. Each union is enclosed in one ball by ball_hull, and for two balls with centres c_A, c_B and radii r_A, r_B the distance is |c_A − c_B| − r_A − r_B, evaluated here from below with dist_down and downward subtractions.

Discs of the complex plane: what _spectrum_kind asks for when nothing in the type of the argument makes the spectrum real.

Intervals of the real line: what _spectrum_kind asks for when the type of the argument guarantees a real spectrum, a HermitianBallMatrix or a real ball matrix whose midpoint and radius are both Diagonal.

_vbd_solve(Bc::Matrix{CT}; sep = -1, maxsteps = 6, mode = :entrywise) -> (W, clusters)

The candidate frame of miyajima2014a_schurnewton, in floating point, on the complexified point matrix Bc: a Schur form; the eigenvalues clustered by their mutual distances at the level sep, √eps · ‖Bc‖_∞ when it is not given; then, with mode = :entrywise, up to maxsteps sweeps of the Newton transform Xᵢⱼ = −Aᵢⱼ/(dᵢ − dⱼ) on the cross-cluster pairs whose gap dominates their off-diagonal entry, or, with mode = :block, one exact block-Sylvester elimination per block boundary, merging two blocks where the transform would be too large to survive the certification; and a QR of each block, which leaves W orthonormal within each block but not unitary overall, so consumers invert it rather than taking its adjoint (see _vbd_unitary_basis). clusters is a Vector{Vector{Int}} of index sets, in general not contiguous.

Nothing here is verified and nothing here needs to be. W is a candidate; every claim rests on the residuals _certify_ball forms against the input ball, so the clustering, the merge rule and the Newton steps can be changed without touching what is proved.

_certify_ball(A::BallMatrix, W; B = nothing, kappa_mode = :cheap)
    -> (transformed, discs, nrmR2, beta, kappa, FB)

The certification step of miyajima2014a_schurnewton, over the candidate frame W, which is not assumed nonsingular. FB is the frame the approximate inverse Y = inv(mid(FB)) is formed of: W for Ax = λx, and B·W for the pencil Ax = λBx, which is the whole of the difference between the two problems. The routine forms the residual R₂ = Y·FB − I and the collapsed enclosure transformed = Y·A·W, a ball product against the input ball, so the radii of A propagate into it; ‖R₂‖_∞ < 1 is required and proves FB and Y nonsingular rather than assuming it.

discs are the Gershgorin discs of transformed, each widened by its own row of the slack, βᵢ = ‖rowᵢ(R₂)‖₁‖transformed‖_∞/(1 − ‖R₂‖_∞), which charges the error Y ≠ FB⁻¹ to disc i alone instead of to every disc at the worst row's rate. nrmR2 is ‖R₂‖_∞, beta the largest βᵢ, and kappa an estimate of the condition number of FB, reported and taking no part in the enclosure.

Throws when ‖R₂‖_∞ ≥ 1: the frame is then not proved nonsingular at this precision and there is nothing to return.

Reference

The two residuals and the Neumann slack are Corollary 3.2 of [25], which states the radii as u + ⟨u⟩_t t with u = |R₁|1, t = |R₂|1 and ⟨v⟩_w = maxᵢ|vᵢ|/(1 − |wᵢ|); the slack is charged per row here rather than once for the whole matrix.

triangular_eigenvectors(L::AbstractMatrix{T}; tol=nothing) where T

Compute eigenvectors and their inverse for a triangular matrix in O(n²).

For a lower triangular matrix L, the eigenvector for eigenvalue λ_i = L[i,i] has v[j] = 0 for j < i, v[i] = 1, and entries for j > i by forward substitution. The resulting V is unit lower triangular, so W = V⁻¹ is also unit lower triangular.

For an upper triangular matrix U, the eigenvector for eigenvalue λ_i = U[i,i] has v[j] = 0 for j > i, v[i] = 1, and entries for j < i by back-substitution. The resulting V is unit upper triangular, so W = V⁻¹ is also unit upper triangular.

Fallback

If min_{i≠j} |L[i,i] - L[j,j]| ≤ tol, falls back to eigen() + inv().

Arguments

• L::AbstractMatrix: A triangular matrix (upper or lower).
• tol: Tolerance for detecting repeated eigenvalues. Defaults to n * eps(real(eltype(L))).

Returns

(V, W, eigenvalues) where V is the eigenvector matrix, W = V⁻¹, and eigenvalues = diag(L), all complex and in floating point: they are approximations, with no error bound. A matrix that is not triangular is an ArgumentError.