API - CertifScripts

CertificationCircle(center, radius; samples = 256)

Discretisation of a circle with centre center, radius radius, and samples equally spaced points used for certification runs.

GramTransform{MT, RT}

Verified data for a Gram-weighted change of variables, as produced by gram_transform.

Fields

• factor::MT: enclosure of a factor L with G = L*L, so ‖x‖_G = ‖Lx‖₂.
• factor_inv::MT: enclosure of L⁻¹.
• cond_factor::RT: rigorous upper bound on κ₂(L) = √κ₂(G).
• cond_gram::RT: rigorous upper bound on κ₂(G). This is the factor the cheap Cbound = √κ₂(G) route would have paid at every point of the contour; the change of variables avoids it entirely.
• gram_residual::Union{Nothing, RT}: rigorous upper bound on ‖L*L - G‖₂ / ‖G‖₂, or nothing when no gram was supplied (then L defines the norm and there is nothing to compare against).
• source::Symbol: :cholesky, :diagonal or :user_factor — how the factor, and hence the certified norm, was obtained.

_bf_ogita_cache_stats()

Return statistics about the BigFloat Ogita cache usage.

_clear_bf_ogita_cache!()

Clear the worker-local BigFloat SVD cache used for Ogita optimization.

_clear_ogita_cache!()

Clear the worker-local SVD cache used for Ogita optimization.

_clear_parametric_cache!()

Clear the worker-local parametric certifier cache.

_compute_schur_bigfloat_direct(A; polynomial = nothing)

Compute BigFloat Schur decomposition directly via GenericSchur.jl (no Float64 seed). Works for matrices whose eigenvalues span many orders of magnitude.

_diagonal_gram_factor(G)

Diagonal fast path: the exact Cholesky factor is L = diag(√dᵢ), enclosed entrywise by the rigorous ball sqrt, and L⁻¹ = diag(1/√dᵢ) by the ball reciprocal. Avoids both the Cholesky and the O(n³) inverse, which matters because a diagonal weight is the common case (orthogonal-polynomial and Fourier bases). Both are stored as Diagonal, so the products of the change of variables cost one directed product per entry of A.

_diagonal_inverse(L::BallMatrix)

Enclosure of L⁻¹ for a factor stored as Diagonal: the ball reciprocal of each diagonal entry, which raises when an entry is not proved to be nonzero. Diagonal storage, O(n).

_evaluate_sample_ogita_bigfloat(T_matrix::BallMatrix, z::Number, idx::Int;
                                 max_iterations::Int=4,
                                 target_precision::Int=256,
                                 distance_threshold::Real=1e-4,
                                 use_cache::Bool=true)

Evaluate the smallest singular value at z using Ogita refinement with caching.

This function uses a worker-local cache to store the last refined BigFloat SVD. For nearby points (distance < distance_threshold), it uses the cached SVD as a starting point for Ogita refinement, which requires fewer iterations than starting from Float64.

Cache Strategy

• If cache hit: use cached BigFloat SVD → 1-2 Ogita iterations
• If cache miss: use Float64 SVD → 3-4 Ogita iterations

Arguments

• T_matrix: The Schur matrix T (BallMatrix)
• z: Sample point on the circle
• idx: Sample index for logging
• max_iterations: Number of Ogita iterations for cache miss (3-4 for 256-bit)
• target_precision: Target precision in bits
• distance_threshold: Maximum distance for cache reuse
• use_cache: Whether to use caching (default: true)

_evaluate_sample_parametric(T, z, idx; use_warm_start=true, distance_threshold=1e-4)

Evaluate sample using the parametric Sylvester-based certifier.

Returns a result compatible with the standard certification format.

_evaluate_sample_with_ogita_cache(T, z, idx; ogita_distance_threshold, ogita_quality_threshold)

Evaluate sample with Ogita optimization: try to refine from cached SVD if available.

If the cached SVD is from a nearby point (distance < ogitadistancethreshold), attempt Ogita refinement. If the refined SVD has acceptable quality (residual < ogitaqualitythreshold), use it; otherwise fall back to full SVD.

_gram_convert(M::BallMatrix, T2)

Re-type the enclosure M so its radii have type T2, rounding outward so that the result still encloses everything M did. This is what lets the Cholesky run in BigFloat — where the factor enclosure is far tighter — while the certifier still receives a Float64 matrix.

_gram_promote(A::BallMatrix, T2)

Widen A to radius type T2. The conversion is exact only when T2 carries at least the precision of T (and, for BigFloat, the exponent range), so a narrowing raises: an enclosure is narrowed with the outward rounding of _gram_convert, never here.

_gram_residual_bound(G, L)

Rigorous upper bound on ‖L*L - G‖₂ / ‖G‖₂. The numerator is bounded with upper_bound_L2_opnorm over the ball product; the denominator needs a rigorous lower bound on ‖G‖₂, for which maxᵢ Gᵢᵢ = maxᵢ eᵢ*Geᵢ serves for Hermitian positive definite G.

The lower end of a real ball, rounded down; negative when the ball straddles zero.

_ogita_cache_stats()

Return statistics about the Ogita cache usage.

_parametric_cache_stats()

Return statistics about the parametric cache usage.

_prepare_gram(A, gram, factor, factor_inv, schur_data; kwargs...)

Shared front end for the certificators. Returns (A′, gram_info) where A′ is either A unchanged (no Gram matrix requested, gram_info === nothing) or the transformed matrix à = L A L⁻¹ together with its GramTransform.

_set_center_svd_cache!(U, S, V, z)

Set the center SVD cache. This should be called once at the start of certification with the SVD computed at the circle center.

_triangular_inverse(L::BallMatrix)

Enclosure of L⁻¹ for upper (or lower) triangular L, by rigorous back- (resp. forward-) substitution against the identity. Tighter than _verify_inverse, but O(n³) in scalar ball operations rather than BLAS-speed.

An upper bound of |x| over the ball: the modulus of the centre bounded above, plus the radius.

_verify_inverse(F::BallMatrix, X::BallMatrix)

Turn an approximate inverse X ≈ F⁻¹ into a rigorous enclosure of F⁻¹.

With R = I - X F we have F⁻¹ = (I - R)⁻¹ X, hence ‖F⁻¹ - X‖₂ ≤ ‖R‖₂‖X‖₂ / (1 - ‖R‖₂) whenever ‖R‖₂ < 1. Every entry of the deviation is bounded by that quantity, so inflating X uniformly by it yields a valid enclosure — and ‖R‖₂ < 1 also proves F invertible. Looser than a triangular back-substitution, but it runs at BLAS speed and applies to an arbitrary invertible factor.

adaptive_arcs!(arcs, cache, pending, η; kwargs...)

Drive the adaptive refinement routine. When job channels are provided the refinement uses asynchronous workers; otherwise the evaluation is carried out serially using the supplied evaluator.

Returns a named tuple (; ok, reason, z, processed, evaluations, remaining). ok is false when the refinement stopped without covering the contour, with reason either :nonpositive, the evaluator returned a bound that was zero, negative or not finite at z, which no amount of bisection repairs, or :cap, the optional maxarcs ceiling on the number of evaluations was reached. maxarcs = 0, the default, leaves the refinement uncapped.

Bisection takes the midpoint of the chord, so after the first refinement the samples lie on the polygon inscribed in the circle rather than on the circle itself. That is the contour the run certifies, which is a closed contour like any other; a caller checking a region condition should use the inradius r·cos(Δθ/2) and not r.

apply_gram_transform(A::BallMatrix, gt::GramTransform) -> BallMatrix

Return the enclosure of à = L A L⁻¹, whose 2-norm resolvent is the G-norm resolvent of A:

‖(A - z)⁻¹‖_G = ‖(Ã - z)⁻¹‖₂    for every z ∉ σ(A) = σ(Ã)

The product is formed in the wider of the factor's type and A's, to which both are widened exactly, and then rounded outward to A's radius type: a BigFloat factor still yields a Float64 enclosure, and a BigFloat A with a Float64 factor is not narrowed.

bound_res_original(l2pseudo, η, norm_Z, norm_Z_inv, errF, errT, N; zmax, Cbound = 1.0,
                   method = :best)

An upper bound of ‖(zI − A)⁻¹‖₂ on the certified contour for the original matrix A, from the samples of its Schur factor T. bound_resolvent_schur gives the bound for T on the polygon, and it is carried to A with errF ≥ ‖Z*Z − I‖ (which equals ‖I − Z Z*‖ for a square Z) and errT ≥ ‖Z T Z* − A‖, by

• method = :bntc: schur_to_original_resolvent, Lemma 5.1 of Blumenthal, Nisoli and Taylor-Crush, with ϵ = max(errF, errT, ‖Z‖ − 1, ‖Z⁻¹‖ − 1);
• method = :defects: schur_to_original_resolvent_defects, the two defects kept separate, sharper by about a factor 2;
• method = :best (default): the smaller of the two, both being valid bounds.

zmax must bound |z| on the contour; for a circle it is |centre| + radius. Inf when no method applies. N is not used.

Reference

A. Blumenthal, I. Nisoli and T. Taylor-Crush, A pseudospectral approach to rigorous numerical estimation of resonances of transfer operators, Found. Comput. Math. (accepted; arXiv:2507.09021), Lemma 5.1.

bound_resolvent_schur(l2pseudo, η)

An upper bound of ‖(zI − T)⁻¹‖₂ over the whole certified contour, from its samples.

The contour is the closed polygon through the samples: adaptive_arcs! starts from the vertices of the regular polygon inscribed in the circle and bisects SIDES, so every later sample is the midpoint of a side. A side [z_a, z_b] is accepted when |z_b − z_a| / σ_a ≤ η, with σ_a a certified lower bound of σ_min(T − z_a I). Since σ_min(T − zI) is 1-Lipschitz in z (Weyl's inequality for singular values: |σ_min(M + E) − σ_min(M)| ≤ ‖E‖₂, here with E = (z_a − z)I), every z on the side has

σ_min(T − zI) ≥ σ_a − |z − z_a| ≥ σ_a − |z_b − z_a| ≥ σ_a (1 − η),

so ‖(zI − T)⁻¹‖₂ ≤ (1/σ_a)/(1 − η). With l2pseudo the largest 1/σ_a over the samples this is l2pseudo/(1 − η), computed with the numerator rounded up and the denominator down. It is Inf for η ≥ 1.

The bound is for the polygon, not for the circle it is inscribed in: between two vertices the circle lies outside the polygon, by up to r(1 − cos(π/N)) for N initial samples, and a point of the spectrum in that sliver, or just outside the circle, can make the resolvent on the circle far larger than on the polygon. Any statement drawn from the bound (a Riesz projector, a count of eigenvalues, a region free of spectrum) is therefore one about the polygon and its interior.

choose_snapshot_to_load(basepath)

Return the most recent valid snapshot stored at basepath.

circle_resolvent_bound(M, circle::CertificationCircle)

Given M ≥ ‖(zI − T)⁻¹‖₂ for every z on the certified polygon (the output of bound_resolvent_schur), an upper bound of ‖(zI − T)⁻¹‖₂ for every z on the CIRCLE circle, or Inf when the argument below does not apply.

Every point z of the circle is within h =polygon_sagitta(circle) of a point z' of the polygon, and σ_min(T − zI) is 1-Lipschitz in z (Weyl's inequality for singular values, |σ_min(M₁ + E) − σ_min(M₁)| ≤ ‖E‖₂, with E = (z' − z)I). Hence

σ_min(T − zI) ≥ σ_min(T − z'I) − h ≥ 1/M − h,

and for M h < 1 the resolvent on the circle is at most M/(1 − M h). The same holds for any matrix in place of T, so it applies equally to a bound for the original matrix. Nothing is assumed about where the spectrum is: when M h ≥ 1 the polygon bound does not exclude a point of the spectrum between the polygon and the circle, and the result is Inf.

compute_schur_and_error(A; polynomial = nothing)

Compute a Schur decomposition of the midpoint of A and certified bounds for its defects: returns (S, errF, errT, norm_Z, norm_Z_inv) with errF ≥ ‖Z*Z − I‖₂, errT ≥ ‖Z T Z* − A‖₂ for every matrix of the ball A, norm_Z ≥ ‖Z‖₂ and norm_Z_inv ≥ ‖Z⁻¹‖₂. When polynomial is provided (coefficients in ascending order), errT bounds ‖Z p(T) Z* − p(A)‖₂ instead.

The precision follows the element type of A. For BigFloat the decomposition is computed directly in that precision by GenericSchur (_compute_schur_bigfloat_direct).

configure_certification!(; job_channel, result_channel, certification_log, snapshot, io)

Cache common resources used by the certification helpers. The stored values are used as defaults by adaptive_arcs! and save_snapshot!. Any keyword may be omitted to keep its previous value.

dowork(jobs, results)

Process tasks received on jobs, computing the SVD certification routine for T - zI. The Schur factor must have been registered in advance with set_schur_matrix!.

dowork_ogita(jobs, results; ogita_distance_threshold=1e-4, ogita_quality_threshold=1e-10)

Process tasks with Ogita optimization enabled. Similar to dowork but tries to use cached SVD from previous evaluations to speed up nearby points.

The worker maintains a cache of the last computed SVD. When a new point z arrives, if it's within ogita_distance_threshold of the cached point, Ogita refinement is attempted. If the refined SVD has acceptable quality (relative residual < ogita_quality_threshold), it's used; otherwise a full SVD is computed.

This is beneficial when the adaptive bisection sends consecutive jobs for nearby points, which happens naturally as arcs get smaller during refinement.

dowork_ogita_bigfloat(jobs, results; target_precision=256, max_ogita_iterations=3,
                      distance_threshold=1e-4)

Process tasks with BigFloat precision using Ogita SVD refinement with caching.

For each job (id, z), this worker:

1. Checks if a cached BigFloat SVD exists from a nearby point
2. If cache hit: refines from cached SVD (1-2 iterations)
3. If cache miss: computes Float64 SVD and refines to BigFloat (3-4 iterations)
4. Certifies with Miyajima bounds
5. Caches the result for future reuse

This is the distributed equivalent of run_certification_ogita for parallel execution. The Schur matrix must be a BigFloat BallMatrix registered with set_schur_matrix!.

Performance

• Cache miss: Ogita from Float64 requires ~4 iterations (10^-16 → 10^-64)
• Cache hit: Ogita from BigFloat requires ~1-2 iterations (already at high precision)
• For adaptive bisection with nearby points, expect 50-90% cache hit rate

dowork_parametric(jobs, results; distance_threshold=1e-4)

Process tasks using the parametric Sylvester-based certifier.

The Schur factor and parametric config must have been registered in advance with set_schur_matrix! and set_parametric_config!.

gram_transform(gram; kwargs...) -> GramTransform

Build the verified data for certifying resolvent bounds in the weighted norm ‖x‖_G = √(x*Gx) induced by the Hermitian positive definite Gram matrix gram, via the change of variables à = L A L⁻¹ with G = L*L.

Pass the result to apply_gram_transform, or — more conveniently — pass gram directly as a keyword to run_certification, run_certification_ogita or run_certification_parametric, which do both steps for you.

The factor comes from verified_cholesky, so positive definiteness of gram is proven, not assumed; the call raises if it cannot be established.

Keyword Arguments

• factor = nothing: use this factor instead of computing a Cholesky, given as a plain matrix or a BallMatrix enclosure. Worth supplying when it is analytically known (diagonal weights, banded factors). Note that with factor the supplied factor defines the certified norm ‖Lx‖₂; it is not checked against gram, and gram may be omitted entirely. A Diagonal, UpperTriangular, LowerTriangular or Symmetric factor keeps its storage, and so do the factor of a diagonal gram (Diagonal) and the verified Cholesky factor (triangular), so that the products of the change of variables go to the kernels of that type; the inverse is stored the same way. Any other plain matrix is converted to a dense one.
• factor_inv = nothing: an approximate inverse of the factor. It is verified before use — an unverified inverse would silently break rigor, so a bad one raises rather than being trusted. Supplying it replaces the O(n³) scalar back-substitution with a BLAS-speed inverse plus verification. It is not used when the factor is a Diagonal: then L⁻¹ is enclosed by the ball reciprocals of its entries. Only its midpoints are used, converted to the factor's type.
• radius_type = Float64: radius type of the enclosures handed to the certifier.
• arithmetic = BigFloat: the arithmetic of the verified Cholesky, BigFloat or Float64, forwarded to verified_cholesky. The BigFloat factor enclosure is many orders of magnitude tighter (radii ~1e-31 against ~1e-14 for a Float64 Gram matrix), and that tightness propagates straight into Ã, so it is worth the cost; the transformed matrix is rounded outward to radius_type afterwards. Float64 trades enclosure width for speed.
• precision_bits = 256: BigFloat precision when arithmetic = BigFloat.
• inverse_method = :auto: :backsub for the tight O(n³) scalar triangular substitution, :verify for the BLAS-speed approximate-inverse-plus- verification route. :auto picks :backsub for triangular factors with n ≤ 256 and :verify otherwise.
• hermitian_tol = 1e-10: relative tolerance of the Hermitian check on gram.

Example

G = collect(Diagonal([1.0, 4.0, 9.0]))     # weighted ℓ² norm
gt = gram_transform(G)
à = apply_gram_transform(BallMatrix(A), gt)

points_on(circle)

Return the discretisation of circle used for certification.

poly_from_roots(roots::AbstractVector)

Given a list of roots r₁, r₂, …, rₙ, returns the coefficients [a₀, a₁, …, aₙ] of the monic polynomial p(x) = (x - r₁)(x - r₂)…(x - rₙ) so that p(x) = a₀ + a₁x + a₂x² + … + aₙ*xⁿ.

polyconv(a, b)

The coefficients of the product of the polynomials with coefficients a and b, both in ascending order: the discrete convolution c[i + j − 1] = Σ a[i] b[j], in the arithmetic of the entries (rounded to nearest for floating-point entries; no bound comes with it).

polygon_sagitta(circle::CertificationCircle)

An upper bound of the largest distance from a point of the circle to the polygon the certification runs on, the regular polygon with circle.samples vertices inscribed in it: the sagitta r(1 − cos(π/N)), bounded by r(π/N)²/2 since 1 − cos x ≤ x²/2, evaluated with every operation rounded upward and π replaced by an upper bound. Bisection adds points on the sides only, so the refined polygon has the same distance to the circle.

run_certification(A, circle; schur_data = nothing, polynomial = nothing, kwargs...)

Run the adaptive certification routine on circle using a serial evaluator.

Arguments

• A: matrix to certify. Converted to BallMatrix when required.
• circle: CertificationCircle describing the contour used for the adaptive refinement.

Keyword Arguments

• schur_data = nothing: pre-computed Schur data as the 5-tuple (S, errF, errT, norm_Z, norm_Z_inv) returned by compute_schur_and_error. When provided, the expensive Schur computation is skipped entirely. This is useful for reusing the same Schur decomposition across multiple circles or when the default Float64-seeded refinement fails (e.g., matrices with eigenvalues below ~10⁻¹⁶).
• polynomial = nothing: optional coefficients (ascending order) describing a polynomial p. When provided the certification is carried out on p(T) and the returned error corresponds to the reconstruction error of p(A).
• η = 0.5: admissible threshold for the adaptive refinement. Must lie in the open unit interval.
• check_interval = 100: number of processed arcs between progress reports and consistency checks.
• log_io = stdout: destination IO for log messages.
• Cbound = 1.0: constant used by bound_res_original when lifting resolvent bounds back to the original matrix.
• gram = nothing: Hermitian positive definite Gram matrix G. When given, the certification is carried out in the weighted norm ‖x‖_G = √(x*Gx) by the exact change of variables à = L A L⁻¹ with G = L*L, so every returned resolvent bound is a bound in the G-norm. Unlike folding √κ₂(G) into Cbound, this costs no condition-number factor. See gram_transform.
• gram_factor = nothing, gram_factor_inv = nothing: precomputed (and, for the inverse, verified before use) factor data — see gram_transform.
• gram_kwargs = (;): extra keywords forwarded to gram_transform, e.g. (arithmetic = BigFloat, inverse_method = :verify).

The return value is a named tuple containing the computed Schur form, the accumulated certification log, and the resolvent bounds for both the Schur factor and the original matrix. The gram field holds the GramTransform that was applied, or nothing.

run_certification_ogita(A, circle; target_precision=256, kwargs...)

Optimized BigFloat certification using Ogita SVD refinement.

This function is specifically designed for BigFloat precision certification where computing fresh BigFloat SVDs at each sample point is expensive. Instead, it:

1. Computes Schur decomposition with BigFloat refinement
2. At each sample point z, computes Float64 SVD of (T - zI)
3. Refines the Float64 SVD to BigFloat using Ogita's algorithm (3-4 iterations)

Due to quadratic convergence, 4 Ogita iterations from Float64 (~10^-16 error) achieve ~10^-64 error, saturating 256-bit precision.

Performance

Typically 10-100x faster than computing fresh BigFloat SVDs at each point.

Arguments

• A: matrix to certify (will be converted to BigFloat BallMatrix)
• circle: CertificationCircle describing the contour
• target_precision::Int=256: precision in bits for BigFloat
• max_ogita_iterations::Int=3: Ogita iterations (3 is optimal for 256-bit precision)
• gram, gram_factor, gram_factor_inv, gram_kwargs: certify in the weighted norm ‖x‖_G = √(x*Gx) — see run_certification and gram_transform. The transform is applied after promotion to BigFloat, so the Gram factor is computed at target_precision.
• Other kwargs passed to standard certification

run_certification_parametric(A, circle; k=nothing, config=config_v2(), kwargs...)

Run the adaptive certification routine using the parametric Sylvester-based certifier.

This method uses block-diagonalization via Sylvester equation to certify the resolvent norm, which can be more efficient than full SVD for large matrices.

Arguments

• A: matrix to certify. Converted to BallMatrix when required.
• circle: CertificationCircle describing the contour.

Keyword Arguments

• k::Union{Nothing, Int}=nothing: Split index for the Schur block decomposition. If nothing, automatically selects k ≈ n/4.
• config::ResolventBoundConfig=config_v2(): Configuration specifying estimators. Options: config_v1(), config_v2(), config_v2p5(), config_v3().
• polynomial = nothing: optional coefficients for polynomial certification.
• η = 0.5: admissible threshold for adaptive refinement.
• check_interval = 100: number of arcs between progress reports.
• log_io = stdout: destination for log messages.
• Cbound = 1.0: constant for resolvent bound lifting.
• gram, gram_factor, gram_factor_inv, gram_kwargs: certify in the weighted norm ‖x‖_G = √(x*Gx) — see run_certification and gram_transform.

Example

using BallArithmetic
A = randn(ComplexF64, 50, 50)
circle = CertificationCircle(0.0, 1.5; samples=64)

# Use V2 configuration (default)
result = run_certification_parametric(A, circle)

# Use V3 configuration with Neumann bounds
result = run_certification_parametric(A, circle; config=config_v3())

# Specify custom split
result = run_certification_parametric(A, circle; k=10)

save_snapshot!(arcs, cache, log, pending, basepath, toggle)

Persist the current certification state to disk using alternating files.

schur_to_original_resolvent(resolvent_schur, ϵ; zmax, Cbound = 1.0)

The resolvent bound for the original matrix from the one for its approximate Schur factor, on a contour whose points satisfy |z| ≤ zmax. Returns Inf when the hypothesis below is not met.

Reference

A. Blumenthal, I. Nisoli and T. Taylor-Crush, A pseudospectral approach to rigorous numerical estimation of resonances of transfer operators, Found. Comput. Math. (accepted; arXiv:2507.09021), Lemma 5.1. With r_B(z) := ‖(z − B)⁻¹‖:

Let M = Z T Z* + E_M (26) and I = Z Z* + E_Z (27), and assume ‖E_M‖, ‖E_Z‖ ≤ ϵ and ‖Z‖, ‖Z⁻¹‖ ≤ 1 + ϵ. Suppose z ∈ ℂ ∖ σ(T) is such that

2ϵ(1 + ϵ)² r_T(z) · max{1, |z|} < 1/2.                                   (28)

Then

r_M(z) ≤ 2(1 + ϵ)² r_T(z) / (1 − 2ϵ(1 + ϵ)² r_T(z)).                      (29)

The factor max{1, |z|} is not decoration: Z T Z* is not similar to T unless Z is unitary, z − Z T Z* = Z(z − T)Z* + z(I − Z Z*), and the hypothesis is what absorbs the second term. With Z = √1.001·I, T = diag(1000, 0) and z = 1001 the matrix z − Z T Z* is singular while r_T(z) = 1; (28) fails there, and this function returns Inf.

resolvent_schur must bound r_T(z) on the whole contour, ϵ the four defects, and zmax the moduli of its points. (28) is checked with r_T and |z| replaced by these upper bounds, which is stronger than the hypothesis at each point, and (29), which increases with r_T, is evaluated there with every operation rounded upward. Cbound multiplies the result, for a caller converting to another norm.

schur_to_original_resolvent_defects(resolvent_schur, δ, e; zmax, Cbound = 1.0)

The resolvent bound for the original matrix A from the one for its approximate Schur factor T, in terms of the two defects separately: δ ≥ ‖I − Q*Q‖ and e ≥ ‖Q T Q* − A‖, on a contour with |z| ≤ zmax. With R ≥ ‖(zI − T)⁻¹‖ on the contour,

R₀ := R / (1 − δ − R·zmax·δ),        ‖(zI − A)⁻¹‖ ≤ R₀ / (1 − R₀ e),

provided both denominators are positive; Inf otherwise. It is sharper than schur_to_original_resolvent, which loses a factor 2: to first order this is R(1 + δ + R·zmax·δ + R e).

Where it comes from, and what is corrected

The route is that of Section 3.1 of

I. Nisoli, Certified spectral approximation of transfer operators and the Gauss map, arXiv:2602.19435 (2026): with the proxy S₀ := Q T Q* and E := A − S₀, Lemma 3.2 there (a Neumann series) gives ‖(zI − A)⁻¹‖ ≤ ‖(zI − S₀)⁻¹‖/(1 − ‖(zI − S₀)⁻¹‖‖E‖), and Theorem 3.4 combines it with a bound of ‖(zI − S₀)⁻¹‖ by ‖(zI − T)⁻¹‖.

That second bound is NOT taken from the paper. Its Lemma 3.3 states ‖(zI − S₀)⁻¹‖ ≤ κ(Q)²‖(zI − T)⁻¹‖ for every z, and the proof writes S₀ = U(HTH)U* for the polar decomposition Q = UH and then inverts zI − HTH as H(zI − T)H. These differ by z(I − H²): Q T Q* is similar to T only when Q is unitary. For Q = √1.001·I, T = diag(1000, 0) and z = 1001, κ(Q) = 1 and ‖(zI − T)⁻¹‖ = 1, while zI − S₀ is singular. The statement used here, with the term restored, is proved as follows.

Claim. Let δ = ‖I − Q*Q‖ < 1, P := Q*Q, and R ≥ ‖(zI − T)⁻¹‖. If R|z|δ < 1 − δ, then zI − S₀ is invertible and ‖(zI − S₀)⁻¹‖ ≤ R/(1 − δ − R|z|δ).

Proof. P is Hermitian with spectrum in [1 − δ, 1 + δ], so Q is invertible, Q⁻¹Q⁻* = P⁻¹, ‖Q⁻¹‖² = ‖P⁻¹‖ ≤ 1/(1 − δ), and ‖P⁻¹ − I‖ ≤ δ/(1 − δ). Then

zI − S₀ = Q (z Q⁻¹Q⁻* − T) Q* = Q ((zI − T) + z(P⁻¹ − I)) Q*,

and (zI − T) + z(P⁻¹ − I) = (zI − T)(I + (zI − T)⁻¹ z (P⁻¹ − I)), whose second factor is invertible by a Neumann series when R|z|δ/(1 − δ) < 1, with inverse of norm at most 1/(1 − R|z|δ/(1 − δ)). Hence

‖(zI − S₀)⁻¹‖ ≤ ‖Q⁻¹‖² · R / (1 − R|z|δ/(1 − δ)) ≤ R / (1 − δ − R|z|δ).  ∎

The corresponding result with the |z| dependence placed in a hypothesis, at the price of a factor 2, is Lemma 5.1 of Blumenthal, Nisoli and Taylor-Crush; see schur_to_original_resolvent.

The bound increases with R, δ, e and |z|, so it is evaluated at upper bounds of the four, every operation rounded outward. Cbound multiplies the result.

set_parametric_config!(precomp, config; k = 0)

Register in this process the data of the parametric certifier: precomp from sylvester_resolvent_precompute (it holds the similarity and the enclosure of its residual), the ResolventBoundConfig, and the split index k for reference.

set_schur_matrix!(T)

Store the Schur factor T used by dowork.