GKWExperiments.jl

GKWExperiments

High-level interface for experimenting with the GKW transfer operator in ArbNumerics. The package re-exports helper routines for computing Hurwitz and Dirichlet zeta values as well as discretisations of the transfer operator that can be used in numerical experiments.

The certification infrastructure is provided by BallArithmetic.jl, which this package re-exports for convenience.

mid(x)

Return the midpoint of an ArbNumerics ball x.

This is a thin wrapper around ArbNumerics.midpoint that is convenient to re-export alongside the rest of the package API.

build_Ls_matrix_arb(s::ArbComplex{P}; K::Int=64, N::Int=1024, prec::Int=P)
    -> M::Matrix{ArbComplex{P}}

Return the matrix M of size (K+1)×(K+1) where the k-th column is the vector of Taylor coefficients b_0,…,b_K of (L_s (w-1)^k)(w) = ∑ b_j (w-1)^j, computed by:

1. precomputing a table Z[j+1,m] = ζ(2s + j, w_m + 1) on the unit circle,
2. evaluating (L_s (w-1)^k)(w_m) via the binomial/Hurwitz identity,
3. extracting coefficients by DFT.

Use N > K for reliable coefficient extraction.

coeffs_from_boundary(vals) -> c

Return the (unnormalized) DFT of boundary samples vals divided by N. On the unit circle |w-1|=1, this yields the Taylor coefficients at w=1: f(1+e^{iθ}) = ∑{k≥0} ak e^{ikθ} ⇒ a_k ≈ (DFT(vals)[k+1]) / N.

gkw_matrix_direct(s::ArbComplex{P}; K::Int=20, prec::Int=P) -> Matrix{ArbComplex{P}}

Assemble the (K+1)×(K+1) Galerkin matrix of the GKW transfer operator using its direct series expansion involving Hurwitz zeta values at the point a = 2.

The entry with zero-based indices (ℓ, k) is computed as

(-1)^{k+ℓ} ∑_{j=0}^k (-1)^j binomial(k, j) (2s + j)_ℓ / ℓ! * ζ(2s + j + ℓ, 2),

where (2s + j)_ℓ denotes the rising factorial. The keyword argument prec controls the Arb precision used for the intermediate Hurwitz zeta evaluations.

gkw_matrix_direct_fast(s::ArbComplex{P}; K::Int=20, prec::Int=P, threaded::Bool=false)
    -> Matrix{ArbComplex{P}}

Same result as gkw_matrix_direct, but ~K/2 times faster by precomputing a table of (2s+j)_ℓ / ℓ! via the recurrence

R[j, ℓ+1] = R[j, ℓ] × (2s + j + ℓ) / (ℓ + 1),   R[j, 0] = 1.

This reduces the inner-loop work from O(ℓ) to O(1) per term, changing the total complexity from O(K⁴) to O(K³) Arb multiplications.

All arithmetic is in ArbNumerics, so error balls are tracked rigorously. Set threaded=true to parallelize column assembly over Threads.nthreads().

values_Ls_fk_from_table!(vals, Z, k)

In-place fill of vals[m] = (L_s (w-1)^k)(w_m) using the precomputed table Z with the binomial identity: (Ls (w-1)^k)(w) = ∑{j=0}^k (-1)^{k-j} binom(k,j) ζ(2s+j, w+1).

• vals must have length N matching size(Z,2).
• Z is the table from zeta_shift_table_on_circle.
• k ≥ 0 is the monomial degree.

zeta_shift_table_on_circle(s::ArbComplex{P}; K::Int, N::Int=1024, prec::Int=P)
    -> Z::Matrix{ArbComplex{P}}, w::Vector{ArbComplex{P}}

Precompute Z[j+1,m] = ζ(2s + j, w_m + 1) for j = 0..K and grid points wm = 1 + e^{i θm} with θ_m = 2π m / N on the unit circle centered at 1.

This supports the identity (Ls (w-1)^k)(w) = ∑{j=0}^k (-1)^{k-j} binom(k,j) ζ(2s+j, w+1).

dirichlet_zeta(s::ArbComplex{P}; prec::Int = P) -> ArbComplex{P}

Compute the (ordinary) Dirichlet zeta function ζ(s) using libarb’s acb_dirichlet_zeta interface. The precision defaults to that of s if not provided.

hurwitz_zeta(s::ArbComplex{P}, a::ArbComplex{P}; prec::Int = P) -> ArbComplex{P}

Compute the Hurwitz zeta function ζ(s, a) using libarb’s acb_hurwitz_zeta interface. The precision defaults to that of s if not provided.

Constants

GKW-specific mathematical constants and bounds from the Hardy space framework. Implements the operator norm bounds C₂ and truncation error Δ from the reference, as well as H²(D_r) whitening transforms for polynomial perturbation analysis.

_arb_to_bigfloat_upper(x)

Convert an ArbReal ball to a rigorous BigFloat upper bound.

Returns midpoint(x) + radius(x) converted to BigFloat with upward rounding. At matching precision (BigFloat ≥ Arb), the midpoint conversion is exact and no conversion error needs to be tracked.

_arb_to_float64_upper(x)

Convert an ArbReal ball to a rigorous Float64 upper bound.

Returns midpoint(x) + radius(x) + conversion_error, converted to Float64 with upward rounding. The conversion error accounts for the precision lost when truncating the Arb midpoint to Float64.

compute_C2(N::Integer)

Compute the operator norm bound C₂ for the GKW transfer operator $L: H²(D₁) → H²(D_{3/2})$ using the splitting formula from Theorem 3.1:

\[C_2(N) = \sum_{n=1}^{N-1} \frac{\sqrt{2n+1}}{(n-1/2)^2} + \sqrt{2 + \frac{2}{N-1/2}} \cdot \zeta(3/2, N-1/2)\]

Larger N gives tighter bounds (see Table 1 in reference).

compute_Δ(K; N=1000)

Compute the truncation error bound Δ for the rank-K Galerkin approximation of the GKW operator (Corollary 4.1 in reference):

\[\|L_1 - (L_1)_K\|_{H^2(D_1)} \leq C_2 \cdot (2/3)^{K+1} =: \Delta(K)\]

h2_whiten(A::AbstractMatrix, r::Real)

Apply the H²(D_r) whitening transform to matrix A.

The H²(Dr) operator norm is ``\|M\|{(r)} = \|Cr M Cr^{-1}\|2$where$Cr = \text{diag}(r^n)_{n=0}^{N-1}``.

This function returns the whitened matrix $C_r A C_r^{-1}$, so that Euclidean norms of the result equal H²(D_r) norms of the original.

is_certified(cert, K; polynomial=nothing)

Check if an eigenvalue certification result satisfies the small-gain condition $α = \|resolvent\| \cdot Δ(K) < 1$ for the GKW operator.

Returns (is_certified::Bool, certified_resolvent_bound).

lr_power_bounds_from_Ak(alpha::AbstractVector, ε::Real)

Compute binomial bounds $β_m ≤ \|L_r^m\|$ for $m = 0, 1, \ldots, L$ using the telescoping identity:

\[\beta_m = \sum_{t=0}^m \binom{m}{t} \varepsilon^t \alpha_{m-t}\]

where alpha[ℓ+1] = ‖A_k^ℓ‖_{(r)} and ε = ‖L_r - A_k‖_{(r)}.

poly_bridge_constant_powers_from_coeffs(c, Ak; r::Real, εr::Real)

Compute the bridge constant $\mathcal{C}_r^{pow}$ for polynomial perturbation bounds from Section 9 of the reference:

\[\mathcal{C}_r^{pow} = \sum_{j \geq 1} |a_j| \sum_{\ell=0}^{j-1} \alpha_\ell \beta_{j-1-\ell}\]

where c = [a₀, a₁, ..., aₐ] are polynomial coefficients, Ak is the discretization matrix, r is the Hardy space radius, and εr = ‖L_r - A_k‖_{(r)} is the discretization error.

Returns (Cr, α, β, s) where:

• Cr is the bridge constant
• α[ℓ+1] = ‖A_k^ℓ‖_{(r)}
• β[m+1] bounds ‖L_r^m‖_{(r)}
• s[k+1] = (α ⋆ β)_k is the discrete convolution

poly_perturbation_bound_powers_from_coeffs(c, Ak; r::Real, εr::Real)

Compute the polynomial perturbation bound $\varepsilon_r \cdot \mathcal{C}_r^{pow}$.

Returns (bound, Cr, α, β, s) where bound = |εr| * Cr.

power_opnorms(B::BallMatrix, L::Integer)

Compute rigorous upper bounds on $\|B^ℓ\|_2$ for $ℓ = 0, 1, \ldots, L$.

Returns a vector norms where norms[ℓ+1] bounds $\|B^ℓ\|_2$. The matrix B should already be whitened if H²(D_r) norms are needed.

Polynomials

Polynomial manipulation utilities for eigenvalue certification. Provides functions for polynomial evaluation, convolution, scaling, powers, deflation, and change of basis.

coeffs_about_0_from_about_c(b, c)

Given b = [b₀, b₁, ..., bₐ] with p(z) = Σ bⱼ (z-c)ʲ, return a = [a₀, a₁, ..., aₐ] such that p(z) = Σ aₖ zᵏ.

Formula: aₖ = Σⱼ₌ₖᵈ bⱼ * binom(j, k) * (-c)^(j-k).

This is the inverse of coeffs_about_c_from_about_0.

coeffs_about_c_from_about_0(a, c)

Given a = [a₀, a₁, ..., aₐ] with p(z) = Σ aₖ zᵏ, return b = [b₀, b₁, ..., bₐ] such that p(z) = Σ bⱼ (z-c)ʲ.

Formula: bⱼ = Σₖ₌ⱼᵈ aₖ * binom(k, j) * c^(k-j).

This is useful for recentering a polynomial around a different expansion point, e.g., for local analysis near an eigenvalue.

deflation_polynomial(zeros, λ_tgt; q::Int=1)

Return coefficients c of p(x) = (α * P(x))^q where P(x) = ∏(1 - x/ζᵢ) zeros out the given zeros, and α is chosen so that p(λ_tgt) = 1.

This is used to deflate already-certified eigenvalues when certifying subsequent eigenvalues.

poly_scale(coeffs, α)

Scale a polynomial by a constant: return coefficients of α * p(x).

polypow(coeffs, n)

Compute the n-th power of a polynomial via repeated convolution. Returns coefficients of p(x)ⁿ.

polyval(coeffs, x)

Evaluate the polynomial with coefficients coeffs at point x. Uses Horner's method for numerical stability. coeffs = [a₀, a₁, ..., aₙ] represents p(x) = a₀ + a₁x + ... + aₙxⁿ.

polyval_derivative(coeffs, x)

Evaluate the polynomial p(x) and its derivative p'(x) simultaneously using a single Horner pass.

coeffs = [a₀, a₁, ..., aₙ] represents p(x) = a₀ + a₁x + ... + aₙxⁿ.

Returns (p_val, dp_val).

EigenspaceCertification

Certify GKW eigenvalues and project the constant function 1 onto eigenspaces.

This module provides the GKW-specific interface for eigenspace certification, leveraging BallArithmetic.jl's rigorous block Schur decomposition (VBD) methods.

GKWEigenCertificationResult

Result of certifying GKW eigenvalues and projecting the constant function onto eigenspaces.

Fields

• block_schur: The rigorous block Schur decomposition from BallArithmetic.jl
• projections_of_one: P_k * [1, 0, ..., 0] for each cluster
• projection_coefficients: Leading coefficient (first entry) of each projection
• gkw_matrix: The discretized L_s matrix with rigorous error bounds
• s_parameter: The s parameter used
• discretization_size: The matrix size K+1

Base.show(io::IO, result::GKWEigenCertificationResult)

Pretty-print certification results.

arb_to_ball_matrix(M_arb::Matrix{ArbComplex{P}}) where {P}

Convert an ArbNumerics complex matrix to a Float64 BallMatrix.

Warning: This truncates the Arb midpoint to Float64, losing ~450 bits of precision. The conversion error is tracked rigorously in the radius, but the center is only Float64-accurate. For full-precision conversion, use BallMatrix(BigFloat, M_arb) instead (provided by BallArithmetic's ArbNumericsExt).

Delegates to BallMatrix(M_arb) from BallArithmetic's ArbNumericsExt.

bigfloat_ball_to_float64_ball(M::BallMatrix) → BallMatrix{Float64}

Convert a BigFloat BallMatrix to Float64 with rigorous error enclosure.

The radius includes both the original BigFloat uncertainty and the conversion error from BigFloat to Float64 centers. Follows the same pattern as arb_to_ball_matrix.

certify_gkw_eigenspaces(s::ArbComplex{P};
                        K::Int=64,
                        num_clusters::Union{Int,Nothing}=nothing,
                        hermitian::Bool=false) where {P}

Certify eigenvalue clusters of the GKW transfer operator and compute rigorous projections of the constant function 1 onto each eigenspace.

Uses BallArithmetic.jl's rigorous block Schur decomposition (Miyajima VBD) to:

1. Automatically cluster eigenvalues based on Gershgorin disc overlap
2. Provide rigorous error bounds on the decomposition
3. Compute spectral projectors for each cluster

Arguments

• s::ArbComplex{P}: The GKW parameter (e.g., s=1 for classical Gauss map)
• K::Int=64: Discretization size (matrix is (K+1) × (K+1))
• num_clusters::Union{Int,Nothing}=nothing: If specified, only compute projections for the first num_clusters clusters (sorted by the VBD ordering)
• hermitian::Bool=false: Whether to treat the matrix as Hermitian

Returns

GKWEigenCertificationResult containing the block Schur decomposition, projectors, and projections of the constant function.

Example

using GKWExperiments, ArbNumerics

s = ArbComplex(1.0)  # Classical GKW operator
result = certify_gkw_eigenspaces(s; K=32)

# Access eigenvalue clusters
println("Number of clusters: ", num_clusters(result))
println("Cluster intervals: ", result.block_schur.vbd_result.cluster_intervals)

# Projection of constant function 1 onto first cluster
println("P₁(1) leading coeff = ", result.projection_coefficients[1])

# Verification bounds
println("Residual norm: ", result.block_schur.residual_norm)
println("Orthogonality defect: ", result.block_schur.orthogonality_defect)

compute_cluster_projector(block_schur::RigorousBlockSchurResult, cluster_idx::Int)

Compute the spectral projector for a single cluster from the block Schur decomposition.

The projector Pk projects onto the invariant subspace corresponding to cluster k. In the Schur basis, this is simply a block of the identity; we transform back to the original coordinates via Pk = Q * Pschurk * Q'.

Arguments

• block_schur: Result from rigorous_block_schur
• cluster_idx: Index of the cluster (1-indexed)

Returns

• BallMatrix: The spectral projector for the specified cluster

float64_ball_to_bigfloat_ball(M::BallMatrix) → BallMatrix{BigFloat}

Promote a Float64 BallMatrix to BigFloat precision.

Float64 → BigFloat conversion is exact, so centers and radii are preserved without any additional error.

project_constant_function(P::BallMatrix, K::Int)

Project the constant function 1 = [1, 0, 0, ..., 0] (in monomial basis) using projector P.

In the monomial basis {(w-1)^k} for k=0..K, the constant function 1 = (w-1)^0 is represented as the unit vector e_1 = [1, 0, 0, ..., 0].

Arguments

• P::BallMatrix: Spectral projector matrix
• K::Int: Discretization size (matrix is (K+1) × (K+1))

Returns

• BallVector: The projected vector P * e_1 with rigorous error bounds

InfiniteDimensionalLift

Rigorous finite-to-infinite dimensional certification for the GKW transfer operator.

This module implements the resolvent bridge and spectral stability framework from the reference paper, allowing certified bounds on the spectrum of the infinite-dimensional operator Lr : H²(Dr) → H²(D_r) from finite-dimensional Galerkin approximations.

All computations work in a single Hardy space H²(D_r) for r ∈ {1, 3/2}.

Key Results Implemented

1. Truncation Error (Corollary 4.1): ‖Lr - Ak‖_{(r)} ≤ C₂(2/3)^{K+1}

2. Resolvent Bridge (Lemma 4): If εK ‖R{Ak}(z)‖{(r)} < 1, then z ∈ ρ(L_r)

3. Projector Approximation: Bound on ‖P{Lr}(Γ) - P{Ak}(Γ)‖_{(r)}

4. Newton-Kantorovich Error Bounds: Rigorous eigenvalue/eigenvector error in H²(D_r)

DeflationCertificationResult

Result of certifying an eigenvalue via polynomial deflation.

Fields

• eigenvalue_center, eigenvalue_radius, eigenvalue_ball: enclosure of the true eigenvalue
• deflation_polynomial_coeffs: coefficients of the deflation polynomial p
• deflation_polynomial_degree: degree of p
• deflation_power: exponent q used in p(z) = (α∏(1-z/λ̂ᵢ))^q
• deflated_eigenvalues: eigenvalues zeroed out by the polynomial
• image_circle_radius: radius of the certification circle in p-space (around 1)
• image_certified_radius: certified inclusion radius in p-space
• poly_perturbation_bound: ‖p(Lr) - p(Ak)‖
• bridge_constant: 𝒞_r^{pow} from the bridge constant computation
• resolvent_Mr: resolvent bound ‖R{p(Ak)}‖ on the image circle
• small_gain_factor: α = epsp · Mr (must be < 1)
• p_derivative_at_target: |p'(λ_tgt)| used for back-mapping
• is_certified: whether the certification succeeded
• truncation_error: base ε_K = C₂(2/3)^{K+1}
• discretization_size: K+1
• hardy_space_radius: r for H²(D_r)
• certification_method: :direct, :parametric, or :ogita
• timing: elapsed seconds for the certification

InfiniteDimCertificationResult

Result of the finite-to-infinite dimensional eigenvalue certification.

All norms are in the H²(Dr) space for the specified `hardyspace_radius`.

Fields

• eigenvalue_center: Approximate eigenvalue from discretization
• eigenvalue_radius: Rigorous error bound |λtrue - λapprox|
• eigenvalue_ball: Ball containing the true eigenvalue
• eigenvector_error: Bound on ‖vtrue - vapprox‖_{(r)}
• truncation_error: ε_K = C₂(2/3)^{K+1}
• resolvent_bound: Certified bound on ‖(zI - Lr)⁻¹‖{(r)} on contour
• small_gain_factor: α = εK · ‖R{A_k}‖ (must be < 1 for certification)
• is_certified: Whether the eigenvalue is rigorously certified
• discretization_size: K+1 (matrix size)
• hardy_space_radius: r ∈ {1, 3/2} for H²(D_r)
• C2_bound: Operator norm bound C₂

OrdschurDirectResult

Result of direct resolvent certification using block Schur structure.

The block triangular inversion formula gives a rigorous resolvent bound:

(zI - T)⁻¹ = [(zI-T₁₁)⁻¹    (zI-T₁₁)⁻¹ T₁₂ (zI-T₂₂)⁻¹]
              [0               (zI-T₂₂)⁻¹                  ]

So ‖(zI-T)⁻¹‖₂ ≤ 1/σ₁₁ · (1 + ‖T₁₂‖/σ₂₂) + 1/σ₂₂

where σ₁₁ = σmin(zI-T₁₁), σ₂₂ = σmin(zI-T₂₂).

TwoStageCertificationResult

Result of the two-stage certification pipeline for a single eigenvalue.

Stage 1 (Klow): Resolvent certification on excluding circles proves simplicity and bounds ‖R{L_r}(z)‖ on the contour.

Stage 2 (K_high): NK certification gives tight eigenpair enclosures.

Transfer bridge: Stage 1 resolvent + Stage 2 truncation error give Riesz projector error bounds.

Fields

• eigenvalue_center, eigenvalue_index: identity of the eigenvalue
• stage1_*: Stage 1 resolvent certification at K_low
• stage2_*: Stage 2 NK certification at K_high
• transfer_*: reverse perturbation from Lr to A{K_high}
• riesz_*: Riesz projector approximation error
• hardy_space_radius, C2_bound: metadata

_ball_to_float64_upper(x) → Float64

Extract a rigorous Float64 upper bound from a Ball (any precision) or scalar.

_bigfloat_ordschur_block(T, Q, select_indices)

Reorder a Schur decomposition so that eigenvalues at select_indices are moved to the top-left block. Returns (T_ord, Q_ord, k) where k = length(select_indices) and T_ord[1:k,1:k] contains the selected eigenvalues.

Works for BigFloat (LAPACK ordschur is Float64 only).

_bigfloat_to_float64_upper(x::BigFloat) → Float64

Convert a BigFloat to Float64 with rigorous upper bound (nextfloat if needed).

_swap_schur_1x1!(T, Q, k)

Swap adjacent eigenvalues at positions k and k+1 in the Schur form using a Givens rotation. Modifies T and Q in place.

backmap_inclusion_radius(r_p, p_coeffs, lambda_tgt; order=1)

Back-map an inclusion radius from p-space to λ-space.

Given a certified radius r_p in p-space (i.e., |p(λ) - 1| ≤ r_p), compute the radius in λ-space using:

First-order (order=1): |λ - λtgt| ≤ rp / |p'(λ_tgt)|

This is valid when p is injective on B(λ_tgt, δ₁). The injectivity condition can be verified using the second-order bound.

Second-order (order=2): Rigorous bound using BallArithmetic.

Evaluates p'' on the disk B(λtgt, δ₁) using Ball interval arithmetic to obtain a rigorous upper bound M₂ = sup{z ∈ disk} |p''(z)|. Then solves the quadratic (M₂/2)δ² - |p'(λtgt)|δ + rp = 0. The smaller root δ₂ is a rigorous inclusion radius by the inverse function theorem: on |λ - λtgt| = δ₂ we have |p(λ) - 1| ≥ |p'|δ₂ - M₂δ₂²/2 = rp, so all solutions of |p(λ) - 1| ≤ r_p lie within the disk of radius δ₂.

Note: δ₂ ≥ δ₁ because the curvature correction widens the bound. The first-order bound is tighter but requires an injectivity assumption; the second-order bound is rigorous without additional assumptions.

Arguments

• r_p: certified radius in p-space
• p_coeffs: polynomial coefficients [a₀, a₁, ..., aₐ]
• lambda_tgt: target eigenvalue
• order: 1 for first-order, 2 for rigorous second-order

Returns

• (radius, dp_abs) where radius is the back-mapped inclusion radius and dpabs = |p'(λtgt)|

certified_resolvent_bound(resolvent_Ak::Real, truncation_error::Real)

Compute the certified resolvent bound for the infinite-dimensional operator.

If α = ε_K · ‖R_{A_k}(z)‖ < 1, then (Lemma 4):

\[\|R_{L_r}(z)\|_{(r)} \leq \frac{\|R_{A_k}(z)\|_{(r)}}{1 - \alpha}\]

Arguments

• resolvent_Ak: Upper bound on ‖(zI - Ak)⁻¹‖{(r)} from finite-dimensional computation
• truncation_error: ε_K = C₂(2/3)^{K+1}

Returns

• (resolvent_Lr::Float64, is_valid::Bool) where is_valid indicates if α < 1

certify_eigenvalue_deflation(A::BallMatrix, lambda_tgt::Number,
                              certified_eigenvalues::AbstractVector;
                              K::Int, r::Real=1.0, N::Int=5000,
                              q::Int=1,
                              image_circle_radius::Real=0.5,
                              image_circle_samples::Int=128,
                              method::Symbol=:direct,
                              backmap_order::Int=1,
                              use_tight_bridge::Bool=true)

Certify an eigenvalue via polynomial deflation.

Pipeline

1. Build deflation polynomial p that zeros out certified_eigenvalues and maps lambda_tgt to 1.
2. Run resolvent certification on a circle around 1 in p-space, using run_certification(A, circle; polynomial=poly_coeffs).
3. Compute polynomial perturbation bound ε_p = ε_K · C_r^{pow}.
4. Check small-gain: ε_p · M_r < 1.
5. Back-map via backmap_inclusion_radius to get eigenvalue enclosure.

Arguments

• A: GKW discretization matrix as BallMatrix
• lambda_tgt: target eigenvalue to certify
• certified_eigenvalues: previously certified eigenvalue centers to deflate
• K: discretization order (matrix is (K+1)×(K+1))
• r: Hardy space radius
• N: splitting parameter for C₂
• q: power of the deflation polynomial
• image_circle_radius: radius of circle around 1 in p-space
• image_circle_samples: number of samples on image circle
• method: :direct, :parametric, or :ogita
• backmap_order: 1 or 2 for back-mapping precision
• use_tight_bridge: if true, use poly_bridge_constant_powers_from_coeffs

Returns

DeflationCertificationResult

certify_eigenvalue_deflation_bigfloat(A_f64::BallMatrix, lambda_tgt::Number,
                                       certified_indices::AbstractVector{<:Integer};
                                       K::Int, r::Real=1.0, N::Int=5000,
                                       q::Int=1,
                                       image_circle_radius::Real=0.5,
                                       image_circle_samples::Int=256,
                                       backmap_order::Int=2,
                                       use_tight_bridge::Bool=true,
                                       use_ordschur::Bool=true,
                                       ordschur_indices=nothing,
                                       schur_data_bf=nothing)

Certify an eigenvalue via polynomial deflation using BigFloat Schur decomposition.

This function targets eigenvalues that are too small for direct resolvent certification (|λ| < 10⁻²⁰). It uses BigFloat precision for the Schur decomposition and deflation polynomial construction, then converts to Float64 for fast SVD certification.

Pipeline

1. Promote A_f64 to BigFloat BallMatrix.
2. Compute BigFloat Schur decomposition (or reuse schur_data_bf).
3. Build deflation polynomial p in BigFloat using Schur diagonal eigenvalues at certified_indices, normalized so p(lambda_tgt) = 1.
4. Evaluate p(T) in BigFloat via Horner on the upper triangular Schur factor.
5. Convert p(T) to Float64 BallMatrix (radii capture conversion error).
6. Certify σ_min(p(T) - zI) > 0 on the image circle around z = 1.
7. Apply Schur bridge to get M_r = ‖(zI - p(A_k))⁻¹‖.
8. Check small-gain: ε_{p,r} · M_r < 1.
9. Back-map via backmap_inclusion_radius to get eigenvalue enclosure.

Arguments

• A_f64: Float64 GKW discretization matrix as BallMatrix
• lambda_tgt: target eigenvalue to certify
• certified_indices: indices (in magnitude-sorted order) of already-certified eigenvalues to deflate — their BigFloat values are taken from diag(T)
• K: discretization order (matrix is (K+1)×(K+1))
• r: Hardy space radius (default 1.0)
• N: splitting parameter for C₂ computation
• q: power of the deflation polynomial
• image_circle_radius: radius of certification circle around 1 in p-space
• image_circle_samples: number of samples on image circle
• backmap_order: 1 (first-order) or 2 (rigorous second-order via p'')
• use_tight_bridge: use poly_bridge_constant_powers_from_coeffs for tighter bound
• use_ordschur: (default true) use ordered Schur decomposition to project away certified eigenvalues, certifying only the smaller p(T₂₂) block. This avoids ill-conditioning when the full p(T) has eigenvalues spanning many orders of magnitude. The block triangular bound (Weyl's inequality) is: σ_min(zI-p(T)) ≥ min(σ_min(zI-p(T₁₁)), σ_min(zI-p(T₂₂))) - ‖B_off‖.
• ordschur_indices: (default nothing → same as certified_indices) indices of eigenvalues to move to the top-left T₁₁ block via ordschur. To avoid ill-conditioning, pass ALL eigenvalues with |λ| > |λ_tgt| so that T₂₂ only contains small eigenvalues. The deflation polynomial zeros (from certified_indices) must be a subset.
• schur_data_bf: pre-computed BigFloat Schur data (5-tuple from compute_schur_and_error); if nothing, computed from A_f64 promoted to BigFloat

Returns

DeflationCertificationResult

certify_eigenvalue_lift(finite_dim_result::GKWEigenCertificationResult,
                        certification_data, cluster_idx::Int;
                        r::Real=1.0, N::Int=1000)

Certify that a finite-dimensional eigenvalue corresponds to a true eigenvalue of the infinite-dimensional GKW operator Lr : H²(Dr) → H²(D_r).

Combines the block Schur decomposition with the resolvent bridge to provide rigorous error bounds on eigenvalues and eigenvectors.

Arguments

• finite_dim_result: Result from certify_gkw_eigenspaces
• certification_data: Result from run_certification containing resolvent bounds
• cluster_idx: Which eigenvalue cluster to certify (1-indexed)
• r: Hardy space radius (1.0 or 1.5)
• N: Splitting parameter for C₂ computation

Returns

InfiniteDimCertificationResult with rigorous error bounds.

Example

using GKWExperiments, ArbNumerics

s = ArbComplex(1.0)  # Classical GKW
finite_result = certify_gkw_eigenspaces(s; K=32)

# Get eigenvalue from VBD
cluster = finite_result.block_schur.clusters[1]
λ_approx = finite_result.block_schur.vbd_result.eigenvalues[cluster[1]]

# Run resolvent certification around this eigenvalue
A = finite_result.gkw_matrix
circle = CertificationCircle(mid(λ_approx), 0.01; samples=128)
cert_data = run_certification(A, circle)

# Combine for infinite-dimensional certification
inf_result = certify_eigenvalue_lift(finite_result, cert_data, 1)

certify_eigenvalue_ordschur_direct(A_f64::BallMatrix, lambda_tgt::Number,
                                   ordschur_indices::AbstractVector{<:Integer};
                                   K::Int, r::Real=1.0, N::Int=5000,
                                   circle_radius::Real=0.0,
                                   circle_samples::Int=256,
                                   schur_data_bf=nothing)

Certify an eigenvalue via direct resolvent bound using block Schur structure.

Unlike polynomial deflation, this works in the ORIGINAL λ-space. No polynomial, no rescaling, no back-mapping. The circle is directly around lambda_tgt.

Why this works

After ordschur, T = [T₁₁ T₁₂; 0 T₂₂] where T₁₁ contains the large eigenvalues (already certified) and T₂₂ contains small eigenvalues including the target. Both blocks are well-conditioned for svdbox:

• T₁₁ entries O(1), σmin(zI-T₁₁) ≥ |λtgt| (distance from z≈λ_tgt to T₁₁ spectrum)
• T₂₂ entries O(|λ{k+1}|), σmin(zI-T₂₂) controlled by eigenvalue separation

The block triangular resolvent formula gives a rigorous bound WITHOUT needing to invert or factorize the full matrix.

Arguments

• A_f64: Float64 GKW discretization matrix as BallMatrix
• lambda_tgt: target eigenvalue to certify
• ordschur_indices: indices (in magnitude-sorted order) of eigenvalues to move to T₁₁ (should include all eigenvalues larger than lambda_tgt)
• K: discretization order (matrix is (K+1)×(K+1))
• r: Hardy space radius (default 1.0)
• N: splitting parameter for C₂ computation
• circle_radius: radius of certification circle around lambda_tgt; if 0 (default), auto-set to half the distance to nearest eigenvalue
• circle_samples: number of samples on circle
• schur_data_bf: pre-computed BigFloat Schur data (5-tuple from compute_schur_and_error); if nothing, computed from A_f64

Returns

OrdschurDirectResult

certify_eigenvalue_schur_direct(A_f64::BallMatrix, schur_position::Int;
                                K::Int, r::Real=1.0, N::Int=5000,
                                circle_radius::Real=0.0,
                                circle_samples::Int=256,
                                schur_data_bf=nothing)

Certify an eigenvalue via direct resolvent bound using the natural Schur order.

This is the rigorous version of certify_eigenvalue_ordschur_direct: it avoids ordschur entirely, using the fact that GenericSchur.jl sorts eigenvalues by magnitude on the diagonal (position 1 = largest). The block split at position p = schur_position gives:

• T₁₁ = T[1:p-1, 1:p-1] (large eigenvalues, well-separated from target)
• T₁₂ = T[1:p-1, p:end] (coupling block)
• T₂₂ = T[p:end, p:end] (target + smaller eigenvalues)

Since no ordschur is performed, the Schur error bounds (‖Q‖, ‖Q⁻¹‖, ‖E‖) from compute_schur_and_error apply directly — no untracked Givens rotation rounding.

The resolvent bound on the certification circle uses:

1. BigFloat SVD (GenericLinearAlgebra) for σ_min(z₀I-T₁₁) at circle center
2. Weyl propagation: σmin(zI-T₁₁) ≥ σmin(z₀I-T₁₁) - ρ for all z on circle
3. Float64 svdbox for σ_min(zI-T₂₂) at each sample point
4. Block triangular formula: ‖(zI-T)⁻¹‖ ≤ 1/σ₁₁·(1+‖T₁₂‖/σ₂₂) + 1/σ₂₂

Arguments

• A_f64: Float64 GKW discretization matrix as BallMatrix
• schur_position: position of target eigenvalue on Schur diagonal (1-indexed). The eigenvalue at diag(T)[schur_position] is the target; positions 1:p-1 form T₁₁.
• K: discretization order (matrix is (K+1)×(K+1))
• r: Hardy space radius (default 1.0)
• N: splitting parameter for ε_K computation
• circle_radius: certification circle radius around the target eigenvalue; if 0, auto-set to half the distance to nearest eigenvalue
• circle_samples: number of sample points on circle
• schur_data_bf: pre-computed BigFloat Schur data (5-tuple from compute_schur_and_error)

Returns

OrdschurDirectResult with certification_method = :schur_direct

certify_eigenvalue_simple(A::BallMatrix, λ_approx::ComplexF64, circle_radius::Real;
                           K::Int, r::Real=1.0, N::Int=1000, samples::Int=128)

Simplified interface for eigenvalue certification without pre-computed results.

Runs the resolvent certification and returns the infinite-dimensional bound.

Arguments

• A: GKW discretization matrix as BallMatrix
• λ_approx: Approximate eigenvalue to certify
• circle_radius: Radius of certification contour
• K: Discretization order (matrix size is K+1)
• r: Hardy space radius
• N: Splitting parameter for C₂
• samples: Number of samples on certification circle

Returns

InfiniteDimCertificationResult with rigorous bounds.

deflation_truncation_error(poly_coeffs::AbstractVector, Ak_norm::Real,
                            Lr_norm::Real, base_truncation_error::Real)

Compute the truncation error for p(Lr) - p(Ak) from Section 9.

\[\|p(L_r) - p(A_k)\|_{(r)} \leq \varepsilon_r \cdot \mathfrak{C}_r(p; A_k, L_r)\]

where the bridge constant is bounded by:

\[\mathfrak{C}_r \leq \sum_{j=1}^d |a_j| \cdot j \cdot C_2^{j-1}\]

Arguments

• poly_coeffs: Polynomial coefficients [a₀, a₁, ..., aₐ]
• Ak_norm: ‖Ak‖{(r)} (or use C₂ for crude bound)
• Lr_norm: ‖Lr‖{(r)} (or use C₂ for crude bound)
• base_truncation_error: εr = ‖Lr - Ak‖{(r)}

Returns

• Truncation error ‖p(Lr) - p(Ak)‖_{(r)}

eigenvalue_inclusion_radius(λ_approx::Number, resolvent_on_circle::Real,
                             circle_radius::Real, truncation_error::Real)

Compute the eigenvalue inclusion radius using the resolvent bridge.

If Γ is a circle of radius r around λapprox, and the small-gain condition `α = εK · sup{z∈Γ} ‖R{Ak}(z)‖ < 1` is satisfied, then Γ ⊂ ρ(Lr). By the maximum principle, L_r has exactly one eigenvalue inside Γ.

Arguments

• λ_approx: Approximate eigenvalue from discretization
• resolvent_on_circle: Maximum of ‖(zI - A_k)⁻¹‖ on Γ
• circle_radius: Radius r of the circle Γ
• truncation_error: ε_K

Returns

• (radius::Float64, is_certified::Bool)

newton_kantorovich_error(resolvent_Ak::Real, left_eigenvector_norm::Real,
                          right_eigenvector_norm::Real, truncation_error::Real;
                          kappa::Real=0.0)

Compute Newton-Kantorovich error bound for eigenvalue/eigenvector (Section 8).

All norms are in the same H²(Dr) space (same-space normalization ‖v‖{(r)} = 1).

The bound is:

\[\|(\lambda_A, v_A) - (\lambda, v)\| \leq \frac{\|DF_{A_k}^{-1}\|_{(r)}}{1-\kappa} \|A_k - L_r\|_{(r)}\]

where ‖DF{Ak}⁻¹‖ ≤ max{‖uA‖{(r)} + ‖SA‖{(r)}, ‖vA‖{(r)}} and ‖SA‖{(r)} ≤ sup{z∈Γ} ‖R{Ak}(z)‖{(r)}.

Arguments

• resolvent_Ak: Resolvent bound (gives ‖S_A‖)
• left_eigenvector_norm: ‖uA‖{(r)} in H²(D_r)
• right_eigenvector_norm: ‖vA‖{(r)} in H²(D_r)
• truncation_error: εK = ‖Ak - Lr‖{(r)}
• kappa: Contraction constant (default 0 for first-order bound)

Returns

• (eigenvalue_error, eigenvector_error, total_error)

projector_approximation_error(contour_length::Real, resolvent_on_contour::Real,
                               truncation_error::Real)

Bound the spectral projector approximation error (Equation proj-diff in reference).

\[\|P_{L_r}(\Gamma) - P_{A_k}(\Gamma)\|_{(r)} \leq \frac{|\Gamma|}{2\pi} \sup_{z\in\Gamma} \frac{\|R_{A_k}(z)\|_{(r)}^2 \varepsilon_K}{1 - \varepsilon_K \|R_{A_k}(z)\|_{(r)}}\]

Arguments

• contour_length: |Γ| = 2πr for a circle of radius r
• resolvent_on_contour: sup{z∈Γ} ‖R{Ak}(z)‖{(r)}
• truncation_error: ε_K

Returns

• (error_bound::Float64, is_valid::Bool)

projector_approximation_error_rigorous(contour_length::Float64,
    resolvent_Ak_high::Float64, eps_K_high::Float64)

Rigorous Riesz projector error bound with directed rounding.

\[\|P_{L_r}(\Gamma) - P_{A_{K_{high}}}(\Gamma)\| \leq \frac{|\Gamma|}{2\pi} \cdot \frac{\|R_{A_{K_{high}}}\|^2 \cdot \varepsilon_{K_{high}}}{1 - \varepsilon_{K_{high}} \cdot \|R_{A_{K_{high}}}\|}\]

Returns

• (error_bound, is_valid) tuple

resolvent_bridge_condition(resolvent_norm::Real, truncation_error::Real)

Check the small-gain condition for the resolvent bridge (Lemma 4 in reference).

The condition α = ε_K · ‖R_{A_k}(z)‖ < 1 ensures that z ∈ ρ(L_r).

Returns

• (is_satisfied::Bool, alpha::Float64) where alpha is the small-gain factor

reverse_transfer_resolvent_bound(M_inf::Float64, eps_K_high::Float64)

Reverse perturbation: given ‖R{Lr}(z)‖ ≤ Minf (from Stage 1), bound ‖R{A{Khigh}}(z)‖ at a higher truncation level.

Since Lr = A{Khigh} + (Lr - A{Khigh}) and ‖Lr - A{Khigh}‖ ≤ ε{K_high}, the standard Neumann perturbation gives:

‖R_{A_{K_high}}(z)‖ ≤ M_inf / (1 - M_inf · ε_{K_high})

All arithmetic uses directed rounding for rigorous bounds.

Returns

• (resolvent_Ak_high, alpha_high, is_valid) tuple

verify_spectral_gap(certification_result, contour_center::Number, contour_radius::Real;
                    N::Int=1000)

Verify that a contour lies in the resolvent set of the infinite-dimensional operator.

Uses the resolvent bridge: if ε_K · sup_{z∈Γ} ‖R_{A_k}(z)‖ < 1, then Γ ⊂ ρ(L_r).

Arguments

• certification_result: Result from run_certification containing resolvent bounds
• contour_center: Center of the certification contour (unused, for API consistency)
• contour_radius: Radius of the certification contour (unused, for API consistency)
• N: Splitting parameter for C₂ computation

Returns

• (is_gap::Bool, gap_margin::Float64) where gap_margin = 1 - α