API - Pseudospectra

pseudospectrum_cells(lower, upper, lo, hi; ε, min_halfdiag, max_evaluations = 100_000)
    -> PseudospectrumCells

Divide the rectangle with opposite corners lo and hi into cells that lie outside the set {z : σ_min(A − zI) ≤ ε}, cells that lie inside {z : σ_min(A − zI) < ε}, and cells left undecided.

lower(z) and upper(z) are functions returning certified bounds lower(z) ≤ σ_min(A − zI) ≤ upper(z); upper may return Inf and lower zero. For instance z -> sigma_min_floor(f, z) with f from block_resolvent_floor or svd_frame_floor, or the larger of several such, and z -> sigma_min_upper(A, z).

What is proved

z ↦ σ_min(A − zI) is 1-Lipschitz, since |σ_min(X + E) − σ_min(X)| ≤ ‖E‖₂ and ‖(z − w)I‖₂ = |z − w|. For a cell with centre c whose points are within r of c:

• if lower(c) − r > ε then σ_min(A − zI) > ε on the cell, which is excluded;
• if upper(c) + r < ε then σ_min(A − zI) < ε on the cell, which is included.

A cell that is neither is cut into four, down to cells whose r is below min_halfdiag, which are returned as undecided; so are the cells not yet examined when max_evaluations centres have been used, the cells being examined from the largest to the smallest. r is rounded up and carries a margin for the rounding of the centres of the parts, so that the discs of radius r about the centres of the returned cells cover the rectangle.

The cells that are cut are those the bounds do not decide. With bounds within a factor of σ_min these lie along the level curve σ_min = ε; where lower returns zero and σ_min > ε, they fill that region down to min_halfdiag.

No use is made of the maximum principle for the resolvent norm on cells free of eigenvalues, which would let the bound on the boundary of such a cell serve for its interior.

PseudospectrumCells{T}

Returned by pseudospectrum_cells: excluded, cells on which σ_min(A − zI) > ε at every point; included, cells on which σ_min(A − zI) < ε at every point; undecided, the others; epsilon, the ε the three lists are for, converted to T; evaluations, the number of cell centres at which the bounds were computed.

The three lists together cover the rectangle that was given, so the ε-pseudospectrum within it meets no cell of excluded and contains every cell of included, and what is left undetermined is contained in the union of undecided.

PseudospectrumCell{T}

The closed rectangle {z : |Re(z − centre)| ≤ hx, |Im(z − centre)| ≤ hy}, with fields centre, of type Complex{T}, and the two half sizes hx and hy. A cell carries no bound of its own: which of the three lists of a PseudospectrumCells it is in is what says what was proved on it.

sigma_min_upper(A::BallMatrix, z)

An upper bound of σ_min(M − zI) for every matrix M of the ball A: with v a floating-point right singular vector of the midpoint of A − zI for its smallest singular value, σ_min(M − zI) ≤ ‖(M − zI)v‖₂/‖v‖₂, the numerator bounded above and the denominator below in ball arithmetic. Any nonzero v gives a valid bound; the singular vector is the one that makes it close. One singular value decomposition in floating point for each z.

svd_frame_floor(A::BallMatrix) -> SVDFrameFloor

Prepare, from one singular value decomposition of the midpoint of A, a lower bound of σ_min(A − zI) that costs O(n²) operations at each z. For a ball matrix A the bound holds for every matrix of the ball.

The bound

Let V be any invertible matrix and N(z) := V*(A − zI)*(A − zI)V, Hermitian. Expanding,

N(z) = H − z̄ K − z K* + |z|² G,      H = V*A*AV,   K = V*AV,   G = V*V,

so N(z) is known entrywise from three matrices that do not depend on z. For y = Vx, ‖(A − zI)y‖² = x*N(z)x ≥ λ_min(N(z))‖x‖² and ‖y‖² ≤ ‖V‖₂²‖x‖², hence, when λ_min(N(z)) ≥ 0,

σ_min(A − zI)² ≥ λ_min(N(z)) / ‖V‖₂².

λ_min(N(z)) is bounded below by Gershgorin's theorem for a Hermitian matrix,

λ_min(N(z)) ≥ γ(z) := min_i ( N_ii(z) − Σ_{j≠i} |N_ij(z)| ).

V is the matrix of right singular vectors of the midpoint of A, computed in floating point: then H is close to the diagonal matrix of the squared singular values, G to the identity, and the off-diagonal part of N(z) is −z̄K − zK* up to rounding, which vanishes for a normal matrix with distinct singular values and grows with the departure from normality. Nothing is assumed of V beyond what is checked: H, K, G are computed in ball arithmetic from A and V, ‖V‖₂² = ‖G‖₂ ≤ 1 + ‖G − I‖₂, and V is invertible when ‖G − I‖₂ < 1. The left singular vectors are not used.

The statement is the theorem "Certified σ_min surface from one verified SVD" of I. Nisoli, Enclosure of eigenvalues via approximate diagonalization (unpublished note, 2026), with the diagonal part of the term linear in z kept in the Gershgorin centre and with the deviation of the computed decomposition from an exact one carried by the radii of H, K, G in place of a separate term; the argument is the one above.

SVDFrameFloor{T}

The quantities of svd_frame_floor that do not depend on z: ball matrices H ∋ V*A*AV, K ∋ V*AV and G ∋ V*V for a floating-point matrix V, and gram_norm ≥ ‖V‖₂² (Inf when V is not proved invertible). Evaluate with sigma_min_floor or resolvent_bound.

BlockResolventFloor{T}

The quantities of block_resolvent_floor that do not depend on z: per block the barycentre centres[j], nonnormality[j] ≥ ‖P_j − c_j I‖₂, the coupling radius coupling[j] ≥ ‖(M − Λ)[C_j, :]‖₂ and the ball matrix blocks[j] containing P_j; the slack slack ≥ ‖M − Λ‖₂; and kappa ≥ ‖W‖₂‖W⁻¹‖₂. Evaluate with sigma_min_floor or resolvent_bound.

sigma_min_floor(f::BlockResolventFloor, z; near = false)

A lower bound of σ_min(A − zI): the larger of the two bounds of block_resolvent_floor, or zero when neither is positive. With near = false the blocks enter through |z − c_j| − n_j only, at O(number of blocks) operations. With near = true a certified smallest singular value of P_j − zI is computed for each block of size larger than one, which is what gives a positive bound inside the disc |z − c_j| ≤ n_j of a non-normal block. Every operation is rounded toward a smaller result.

sigma_min_floor(f::SVDFrameFloor, z)

A lower bound of σ_min(A − zI): √(γ(z)/‖V‖₂²) with the quantities of svd_frame_floor, every operation rounded toward a smaller result; zero when γ(z) is not positive or V was not proved invertible.

resolvent_bound(f::BlockResolventFloor, z; near = false)

An upper bound of ‖(A − zI)⁻¹‖₂: the reciprocal, rounded up, of sigma_min_floor; Inf where that is zero.

resolvent_bound(f::SVDFrameFloor, z)

An upper bound of ‖(A − zI)⁻¹‖₂: the reciprocal, rounded up, of sigma_min_floor; Inf where that is zero.

Enclosure

A chain of discs with a certified bound on each.

• λ: the point the chain was built around, as a ComplexF64
• points[i], radii[i]: centre and radius of disc i
• bounds[i]: a ball containing σ_min(T − zI) for every z of disc i
• loop_closure: true when every disc is proved to overlap the next and the last the first. The closed polygon through the centres is then contained in the union of the discs, since a segment between two centres at distance less than the sum of the radii is covered by the two discs. That the polygon winds around λ is not checked.

The bounds are for the matrix T the chain was built from.

compute_enclosure(A::BallMatrix, r1, r2, ϵ; max_initial_newton = 30,
    max_steps = Int64(ceil(256 * π)), rel_steps = 16,
    svd_method = MiyajimaM1(), apply_vbd = false)

Chains of discs around the spectrum of the Schur factor T of the midpoint of A (computed in Float64), each returned as an Enclosure:

• one around each diagonal entry λ of T with r1 < |λ| < r2, following the level curve σ_min(T − zI) = ϵ;
• one along the circle |z| = r1 when T has diagonal entries of modulus below r1;
• one along the circle |z| = r2 when it has diagonal entries of modulus above r2.

max_initial_newton bounds the Newton steps used to reach a level curve, max_steps the length of a chain, and rel_steps sets the step relative to the distance from the nearest diagonal entry.

A level curve is followed by a predictor-corrector continuation, an Euler step along the curve and then one Newton step back to it, both directions read from the singular vectors belonging to the smallest singular value of T − zI; this is the continuation of [35], carried out here in floating point. Nothing about the curve is certified, and the bounds the chain carries do not depend on how well it was followed: the disc of a step is enclosed for the T it was built from, so a step that leaves the level curve widens the bound on that disc and proves nothing false. What the step sizes do govern is whether consecutive discs overlap, which is what loop_closure records.

The bounds are for T, not for A: the radius of A and the defects of the Schur factorisation are not used. CertifScripts.run_certification certifies a circle for A itself. bound_resolvent gives the bound on the union of the discs of a chain, and its loop_closure says whether the chain is proved closed.

check_enclosure(E::Enclosure) -> Bool

Whether every disc of E overlaps the next and the last overlaps the first, recomputed from the centres and radii with directed rounding: the distances between centres rounded up and the sums of radii rounded down, so that a true is proved and a false leaves the question open. This is the test whose outcome E.loop_closure holds; it is recomputed here rather than read back, and the bounds E.bounds are not used.

bound_resolvent(E::Enclosure)

An upper bound of ‖(zI − T)⁻¹‖₂ for every z in the union of the discs of E: the reciprocal of the smallest lower end of the bounds, rounded up; Inf when one of them is not positive.

sylvester_resolvent_precompute(T, k; X_oracle = nothing)

The part of the bound of parametric_resolvent_bound that does not depend on z, for the split of T at index k: T₁₁ = T[1:k, 1:k], T₂₂ = T[k+1:end, k+1:end].

T must have exact entries, be zero in the block [k+1:end, 1:k], and have T₂₂ upper triangular. X_oracle is the matrix X of the similarity; when it is not given, solve_sylvester_oracle computes one. Any X gives a valid bound: the residual R = T₁₂ + T₁₁X − XT₂₂ is enclosed in ball arithmetic for the X used, and its norm and the norm of X are bounded from above (upper_bound_L2_opnorm).

Returns a SylvesterResolventResult; certified is false, with the reason, when the structure is not the required one or X is not finite.

solve_sylvester_oracle(T11::AbstractMatrix, T12::AbstractMatrix, T22::AbstractMatrix)

An approximate solution X of the Sylvester equation T₁₁X − XT₂₂ = −T₁₂, in floating point: no bound comes with it, and none is needed. It is computed by LAPACK, through LinearAlgebra.sylvester, in Float64; for other element types the data are rounded to Float64 and the solution converted back, so for a BigFloat matrix the solution is no more accurate than a Float64 one.

What X is used for is the similarity S(X) = [I X; 0 I] of sylvester_resolvent_precompute, and the identity S(X)⁻¹TS(X) = [T₁₁ R; 0 T₂₂] with R = T₁₂ + T₁₁X − XT₂₂ holds for every X; the residual R is then enclosed in ball arithmetic for the X actually used. An inaccurate X leaves a larger ‖R‖₂ and a larger κ₂(S(X)), and so a weaker bound, and nothing else. A caller with a better X passes it as the X_oracle of the precomputation instead.

SVDWarmStart(U, Σ, V)

An approximate singular value decomposition of zI − T₁₁ at a nearby point, given to parametric_resolvent_bound in place of a fresh decomposition at the current point: U, Σ and V are passed to ogita_svd_refine, the refinement of [6], for two iterations, and the refined triple is then certified by certify_svd.

It changes the approximation the certificate is computed from and nothing else. The certified enclosure of the singular values is of the ball matrix zI − T₁₁, whatever triple is supplied, so a warm start that has drifted from the current point gives a wider enclosure of the same singular values, and parametric_resolvent_bound then declines when the resulting σ_min is not proved positive.

find_optimal_split(T, z; k_range = 2:min(size(T, 1) - 2, 50), config = config_v1())

The split index in k_range with the smallest bound at z, as (k, precomp, result); nothing when the bound applies for no index of the range.

Every index gives a valid bound on ‖(zI − T)⁻¹‖₂, since the identity behind parametric_resolvent_bound holds at each split, so the smallest of them is valid too and the search costs nothing in rigour. It costs one sylvester_resolvent_precompute, one Sylvester solve among it, and one certified singular value computation of size k for each index of the range; the index is chosen for one z, and a precomp is for the split and not for the point, so the returned precomp can be reused at other points.

compare_all_configs(T, k, z; miyajima_method = :M1, power_iterations = 3)

The bound at z for config_v1, config_v2, config_v2p5 and config_v3 on one precomputation, with miyajima_method and power_iterations substituted into each. Returns (; precomp, results, bounds, best, z): results holds the ParametricResolventResult of each configuration and bounds the bound it reached, both indexed by "V1", "V2", "V2.5" and "V3", with Inf where a configuration declined, and best is the key of the smallest entry of bounds.

All four are bounds on the same ‖(zI − T)⁻¹‖₂, so bounds[best] is a bound whenever it is finite; it is Inf when none of the four applied, and the failure_reason of each entry of results says at which step it stopped.

ResolventBoundConfig(norm_estimator, d_inv_estimator, coupling_estimator, miyajima_method,
                     power_iterations, fallback_to_tri)

The choices of parametric_resolvent_bound: a NormEstimator, a DInverseEstimator, a CouplingEstimator, the method for the certified singular values of the small block (:M1 or :M4), the number of iterations of the Collatz bound, and whether TriBacksub replaces a Neumann bound that does not apply (‖N‖ ≥ 1).

config_v1, config_v2, config_v2p5 and config_v3 are the combinations in use; every one of them is a valid bound, so the smallest of several may be taken, which is what compare_all_configs does.

config_v1() -> ResolventBoundConfig

TriBacksub for ‖D_z⁻¹‖₂ and CouplingNone, the product ‖A_z⁻¹‖₂‖R‖₂‖D_z⁻¹‖₂, for the coupling term; OneInfNorm, :M1, three Collatz iterations, and the fallback to TriBacksub enabled. This is the default of parametric_resolvent_bound and the one that needs no approximate solve at each z.

config_v2() -> ResolventBoundConfig

config_v1 with CouplingARSolve in place of CouplingNone: the coupling term is bounded through an approximate solution of A_zY = R and its residual, which costs one solve of size k at each z.

config_v2p5() -> ResolventBoundConfig

config_v1 with CouplingOffDirect: the whole off-diagonal block is approximated by a W with A_zWD_z ≈ R and the residual R − A_zWD_z is enclosed, which costs a solve on each side at each z.

config_v3() -> ResolventBoundConfig

config_v2 with NeumannCollatz2 in place of TriBacksub: ‖D_z⁻¹‖₂ comes from the factorisation D_z = Δ(I + N) with ‖N‖₂ bounded by collatz_upper_bound_L2_opnorm, which applies only where that bound is below one; since fallback_to_tri is set, TriBacksub is used where it is not.

estimate_2norm(M, method::NormEstimator)

An upper bound of ‖M‖₂, holding for every matrix of the ball when M is a BallMatrix and for M itself when it is a matrix with exact entries, by the chosen NormEstimator: √(‖M‖₁‖M‖_∞) for OneInfNorm, which bounds ‖M‖₂ since ‖M‖₂² = ρ(M*M) ≤ ‖M*M‖_∞ ≤ ‖M‖₁‖M‖_∞, and ‖M‖_F for FrobeniusNorm. Neither is always the smaller: ‖M‖_F is the norm itself for a matrix of rank one, while for the identity of size n it is √n and the first bound is one. The operand norms are computed in ball arithmetic and the square root is rounded up.

The bounds of parametric_resolvent_bound use this for the norms of the approximate solutions and of their residuals only, never for a quantity a certified singular value computation already covers.