γₙ constants for floating point error estimation, as in [Higham, Accuracy and Stability of Numerical Algorithms]

Estimate ||I + Q + Q^2 + … || (infinite sum) using a list of computed norm bounds norms[k] ≥ ||Q^k||.

Returns the maximum number of (structural) nonzeros in a row of A

norms_of_powers(BM_U0::BallMatrix, ::Type{L2}, m)

Direct L^2-operator-norm-of-iterates loop on a BallMatrix that the caller has already restricted to U^0 (the average-zero subspace). Returns norms[k] ≥ ‖BM_U0^k‖_{L²→L²} for k = 1, …, m, computed rigorously via interval matrix multiplication and _l2_opnorm_ball, which takes the best of the verified SVD enclosure, the Collatz/Perron–Frobenius bound, and the Hölder bound √(‖A‖_1‖A‖_∞) — the last two are what upper_bound_L2_opnorm already minimises over, using the radius-inflated |c| + r.

This is the basis-free entry point used by callers (e.g. SelfConsistentExperiments.certify_contraction_iterate) that assemble DT_N(u_N)|_{U_0} outside the DiscretizedOperator abstraction — for example the self-consistent linearisation with a rank-one feedback correction, which is not a transfer operator itself. No Basis / restrict_to_average_zero step is performed: the caller is responsible for ensuring BM_U0 already represents the operator on the average-zero subspace.

The coarse-fine refinement (refine_norms_of_powers) is not applied here, since uniform DFLY constants on the discretisation family are not assumed in the Fourier case.

Estimates the norms ||Q||, ||Q^2||, ... ||Q^m|| on U^0.

U is the matrix [ones(1,n-1); -I_(n-1,n-1)]. It is currently assumed that f*U==0 (i.e., all elements of f are equal).

f must be an interval vector.

The following constants may be specified as keyword arguments:

normQ, normE, normv0, normEF, normIEF, normN

otherwise they are computed (which may be slower).

e and f must be specified in case isintegralpreserving==false In case isintegralpreserving is true, they may be specified but they are then ignored.

Implementation note: currently we perform this computation one column at a time, to be able to scale (slowly) to cases with large size; for moderate sizes, it would indeed be better to do the computation all columns at the same time, in BLAS level 3.

Estimates ||Q||, ||Q^2||, ... ||Q^m|| on U^0 using L2 norm via BallArithmetic.

Computes matrix powers as BallMatrices, which automatically track rounding errors. Uses upperboundL2_opnorm for rigorous spectral norm bounds.

Arrays of bounds to $||Q^k||_{w → s} = \sup_{||f||_w=1} ||Q^k f||_s$ and to $||Q^k||_{w}$ coming theoretically from iterated DFLY inequalities (the "small matrix method").

Returns two arrays (strongs, norms) of length m: strongs[k] bounds $||Q^k f||_s$, norms[k] bounds $||Q^k f||$)

Estimate norms of powers from those on a coarser grid (see paper for details)

norms_of_powers_resolvent(m, ρ, M_inf)

Compute power norm bounds via certified resolvent on a contour. Returns norms[k] ≥ ||Q^k|_{complement of eigenvalue 1}||₂.

Uses the Cauchy integral bound: ||Q^k(I-P₁)|| ≤ ρ^{k+1} · M_∞

where ρ is the radius of the certified contour (must satisfy max|λⱼ| < ρ < 1 for j≥2) and M_∞ = sup_{z∈Γ} ||(zI-Q)⁻¹|| is the certified resolvent bound from certify_spectral_gap.

Array of "trivial" bounds for the powers of a DiscretizedOperator (on the whole space) coming from from ||Q^k|| ≤ ||Q||^k

Compute better and/or more estimates of power norms using the fact that $||Q^{k+h}|| ≤ ||Q^k|| \cdot ||Q^h||$. This uses multiplicativity, so it will not work for mixed norms, e.g., $||Q^k||_{s → w}$, or $||M^k|_{U^0}||$ (unless M preserves $U^0$, which is the case for $Q|_{U^0}$).

Restrict a matrix to the average-zero subspace. Default implementation: requires specialization per basis type.

Bounds rigorously the distance of w from the fixed point of Q (normalized with integral = 1), using a vector of bounds norms[k] ≥ ||Qh^k|{U_h^0}||.

Uses power norm bounds already computed for a coarse operator to estimate the same norms for a finer operator

Return a numerical approximation to the (hopefully unique) invariant vector of the dynamic with discretized operator Q.

The vector is normalized so that integral_covector(B)*w ≈ 1

Arrays of bounds to ||Q^k||{w → s} = sup{||f||w=1} ||Q^k f||s and to ||Q^k||_{w} coming theoretically from iterated DFLY inequalities (the "small matrix method").

Returns two arrays (strongs, norms) of length m: strongs[k] bounds ||Q^k f||_s, norms[k] bounds ||Q^k f||)

Estimate norms of powers from those on a coarser grid (see paper for details) TODO: Check if it works for other basis types

Estimates the norms ||Q||, ||Q^2||, ... ||Q^m|| on U^0.

U is the matrix [ones(1,n-1); -I_(n-1,n-1)]. It is currently assumed that f*U==0 (i.e., all elements of f are equal).

f must be an interval vector.

The following constants may be specified as keyword arguments:

normQ, normE, normv0, normEF, normIEF, normN

otherwise they are computed (which may be slower).

e and f must be specified in case isintegralpreserving==false In case isintegralpreserving is true, they may be specified but they are then ignored.

norms_of_powers_sequence_noise(B, N, Q, MK, f, shifts)

Like norms_of_powers_noise, but for a nonautonomous (sequential) cocycle in which a different forcing acts at each step. The k-th step is

v  ←  Shift(shifts[k]) ∘ MK ∘ Q  v ,

where Shift(s) is the (fractional) circular shift by s Ulam bins on the torus. A fractional shift Shift(s) = (1−f)·circshift(·,n₀) + f·circshift(·,n₀+1) (n₀ = ⌊s⌋, f = s − n₀) is a convex combination of two basis permutations: mass-preserving, with L¹ operator norm 1. An integer shift (f = 0) is an exact permutation and introduces no floating-point error; a fractional shift costs one extra rounded mul-add per entry, tracked in the accumulated radius ϵ exactly like the noise step (operator norm 1). Either way the rigorous error propagation of norms_of_powers_noise carries over.

Using fractional shifts lets the forcing track the real kick s_n = A·k·sin(2π θ_n) smoothly in the driving phase θ_n = θ₀ + nω (the irrational rotation), with no rounding of the kick.

shifts has length m (one shift per step); the returned vector norms[k] bounds ‖ L_{shifts[k]} ⋯ L_{shifts[1]} |_{U⁰} ‖, the norm of the length-k forced cocycle on the zero-average subspace. With shifts .= 0 it reproduces norms_of_powers_noise exactly.

Array of "trivial" bounds for the powers of a DiscretizedOperator (on the whole space) coming from from ||Q^k|| ≤ ||Q||^k

Return an upper bound to Q_h*w - w in the given norm

Aη(η) <: NormKind

Analytic norm on the strip of half-width $η$ around the real axis, i.e. $\|f\|_{A_η} = \sum_k |\hat f_k| e^{η|k|}$ on the Fourier coefficients.

Since $|e^{2\pi i k z}| \le e^{2\pi |k| \operatorname{Im} z|}$, this weighted $\ell^1$ norm dominates $\sup |f|$ on the strip, so bounds stated for the sup norm remain valid. It is the strong norm of FourierAnalytic.

Eρ is the Chebyshev counterpart: a strip is the natural neighbourhood of the circle, a Bernstein ellipse the natural neighbourhood of $[-1,1]$.

Eρ(ρ) <: NormKind

Analytic norm on the Bernstein ellipse $E_ρ$ — the image of the circle $|w| = ρ$ ($ρ > 1$) under the Joukowski map $z = (w + w^{-1})/2$, i.e. the ellipse with foci $\pm 1$ and semi-axes $(ρ \pm ρ^{-1})/2$.

$\|f\|_{E_ρ} = \sum_k |\hat b_k| ρ^k$ on the Chebyshev coefficients, matching the convention of Aη, which is likewise a weighted $\ell^1$ norm on the Fourier coefficients. Since $|T_k| \le (ρ^k + ρ^{-k})/2 \le ρ^k$ on $E_ρ$, this dominates $\sup_{E_ρ}|f|$, so bounds stated for the sup norm (such as the Trefethen projection estimates) remain valid. This is the Chebyshev counterpart of Aη, which measures analyticity on a strip for the Fourier bases: a strip is the natural neighbourhood of the circle, an ellipse the natural neighbourhood of $[-1,1]$.

See bernstein_parameter and bernstein_expansion for the rigorous test that a map expands such an ellipse.

L1 <: NormKind

The $L^1$ norm against Lebesgue measure $dx$ on the domain. This is the weak norm of the Ulam and hat bases. See L1μ for the counterpart against the arcsine measure, which is what the Chebyshev bases use.

L2μ <: NormKind
L1μ <: NormKind

$L^2$ and $L^1$ against the arcsine measure $dμ = dx/(π\sqrt{x(1-x)})$ on $[0,1]$, as opposed to L2 and L1, which are against $dx$.

These are the natural weak and auxiliary norms for the Chebyshev bases, because the two facts the analysis rests on both live against $μ$:

• the $T_m$ are orthogonal in $L^2(μ)$, so Parseval holds on the coefficients and the norm interface reduces to the same one-liners as the Fourier bases;
• the coefficients obey $|\hat b_k| = |2\int f T_k\,dμ| \le 2\|f\|_{L^1(μ)}$, uniformly in $k$. There is no such bound against $L^1(dx)$: a width-$ε$ sliver at an endpoint has $\|f\|_{L^1(dx)} \sim ε$ but $\|f\|_{L^1(μ)} \sim \sqrt{2ε}/π$.

Keeping both against $μ$ also makes $\|v\|_{L^1(μ)} \le \|v\|_{L^2(μ)}$ a plain Cauchy–Schwarz with constant 1, since $μ$ is a probability measure; pairing $L^1(dx)$ with $L^2(dμ)$ instead costs a factor $π/(2\sqrt2)$.

Use l2_measure_conversion_bounds to return to $L^2(dx)$.

L2 <: NormKind

The $L^2$ norm against Lebesgue measure $dx$ on the domain, and the induced $\ell^2$ operator norm on matrices. See L2μ for the counterpart against the arcsine measure.

L2μ <: NormKind
L1μ <: NormKind

$L^2$ and $L^1$ against the arcsine measure $dμ = dx/(π\sqrt{x(1-x)})$ on $[0,1]$, as opposed to L2 and L1, which are against $dx$.

These are the natural weak and auxiliary norms for the Chebyshev bases, because the two facts the analysis rests on both live against $μ$:

• the $T_m$ are orthogonal in $L^2(μ)$, so Parseval holds on the coefficients and the norm interface reduces to the same one-liners as the Fourier bases;
• the coefficients obey $|\hat b_k| = |2\int f T_k\,dμ| \le 2\|f\|_{L^1(μ)}$, uniformly in $k$. There is no such bound against $L^1(dx)$: a width-$ε$ sliver at an endpoint has $\|f\|_{L^1(dx)} \sim ε$ but $\|f\|_{L^1(μ)} \sim \sqrt{2ε}/π$.

Keeping both against $μ$ also makes $\|v\|_{L^1(μ)} \le \|v\|_{L^2(μ)}$ a plain Cauchy–Schwarz with constant 1, since $μ$ is a probability measure; pairing $L^1(dx)$ with $L^2(dμ)$ instead costs a factor $π/(2\sqrt2)$.

Use l2_measure_conversion_bounds to return to $L^2(dx)$.

'Absolute value' definition that returns mag(I) for an interval and abs(x) for a real

Computes a rigorous upper bound for z*z'

_l2_opnorm_ball(BM::BallMatrix)

Rigorous upper bound on $\|BM\|_2$, taking the best of the two estimators available: the cheap Collatz/Hölder bound upper_bound_L2_opnorm, and the sharp bound obtained from BallArithmetic's verified SVD, svd_bound_L2_opnorm.

The cheap bound is $\sqrt{\|A\|_1 \|A\|_\infty}$, which overestimates badly for the matrices arising here; the SVD bound is essentially tight but costs an $O(n^3)$ decomposition and needs GenericSchur loaded to work in extended precision. Should the SVD be unavailable or fail, the cheap bound is returned on its own. Because the result is the min of the two, it is a valid upper bound either way.

Rigorous upper bound on a vector norm. Note that Linf, L1 are the "analyst's" norms

Rigorous upper bound on the L2 norm of a vector, using Parseval identity.

Certified upper bound to ||A|| (of specified NormKind)

Certified upper bound to the L2 operator norm of a matrix via BallArithmetic.

upper_bound_L2_opnorm is min(Collatz, √(‖·‖₁‖·‖_∞)): cheap, but on the matrices that arise here it overestimates the spectral norm by 1.4–1.7×. The verified-SVD enclosure is sharp to ~1e-8 relative and costs a few ms at the sizes we use, and the weak norm feeds straight into γ_N, so we prefer it and keep the min with the cheap bound — the result can only improve.

Midpoints may be real or complex, and the element type may be BigFloat — svdbox handles it through BallArithmetic's GenericSchurExt, so load GenericSchur alongside this package to get the sharp bound in extended precision. Should the SVD be unavailable or fail, we fall back to the cheap bound; since we return the min of the two, the result is a valid upper bound either way.

These functions compute a rigorous upper bound for the 2-norm of a vector; we have a specialized version for complex numbers to avoid taking the sqrt root and squaring again

Types to compute norms iteratively by "adding a column at a time".

Create a new NormCacher to compute the normbound of the empty matrix with n rows

Update a NormCacher to add one column to the matrix it is computing a norm of. This column may be affected by an error ε (in the same norm).

Return the norm of the matrix the NormCacher is working on.

convergencerateabstract(Bas::Ulam, D::Dynamic, norms)

Estimate the strong norm of $||L^n|_{U_0}||_s$ from norms, the bounds on the weak norm of the discretized operator

$||L_{h}^n|_{U_0}||_w$

This method was developed in Stefano Galatolo, Isaia Nisoli, Benoît Saussol. An elementary way to rigorously estimate convergence to equilibrium and escape rates. Journal of Computational Dynamics, 2015, 2 (1) : 51-64. doi: 10.3934/jcd.2015.2.51

eig_costants_small_matrix(A)

Return the dominant eigenvalue and associated left eigenvector of the small matrix

abstract_weak_norm_bound(B::Basis; dfly_coefficients)

Float64 upper bound on $\|L_k u\|_w / \|u\|_w$ for $u \in U_h$, derived entirely from the basis interface via the vector inequality $\|v\|_w \le S_1\|v\|_s + S_2 |||v|||$ (weak_by_strong_and_aux_bound) combined with the DFLY inequality and the norm-equivalence constants on $U_h$:

\[\|L_k u\|_w \;\le\; S_1\bigl(a\,M_{1n} + b\,M_2\bigr)\|u\|_w + S_2\,M_2\,\|u\|_w, \qquad u \in U_h.\]

Assumes the aux norm is contracted by $L$ ($|||L u||| \le |||u|||$), which holds for transfer operators of integral-preserving maps with the typical aux = L¹ convention. For bases where the standard bound_weak_norm_abstract gives a sharper estimate, prefer that or pass a computed bound from powernormbounds directly.

bN_constant(a, b, M, N)

Float64 upper bound on the Lasota–Yorke iteration constant

\[b_N := b\sum_{j=0}^{N-1} a^{N-1-j} M^j,\]

from Lemma A.4 of [Nisoli, "Certified spectral approximation of transfer operators and the Gauss map"]. With a, b, M ≥ 0, this is the same $b_N$ appearing in

\[\|L^N u\|_s \;\le\; a^N \|u\|_s + b_N \|u\|_w,\]

and inside the perturbation factor of Proposition A.14. N ≥ 1; for N = 1 the sum has the single term $j = 0$, giving $b_1 = b$.

Caller note (package convention)

This helper is paper-faithful, with a single $M$. In the package the DFLY is on (strong, aux), so the correct value of M here is the aux-norm contractivity of $L$ — typically 1 for L¹-preserving transfer operators with aux_norm = L¹. The b passed in should be $b^{\text{pkg}} \cdot M_2$ where $M_2 = \mathrm{aux\_weak\_bound}(B)$.

In particular, do not pass bound_weak_norm_abstract(B, D) as M: that quantity (which can exceed 1) leads to b_N = O(M^N) blowing up and to Prop A.14 never closing. See coarse_fine_weak_resolvent for the basis-aware wiring that handles this conversion automatically.

coarse_fine_weak_resolvent(B::Basis, D::Dynamic, abs_z, K_z;
                           N, M_k_weak=..., kwargs...)

Float64 upper bound on the fine weak resolvent $\mathcal{R}_w(z, L_k)$ of the finite-rank discretization, from Proposition A.14 (DFLY appendix) of [Nisoli, "Certified spectral approximation of transfer operators and the Gauss map"]:

\[\mathcal{R}_w(z, L_k) \;\le\; \frac{S_N^{(k)}(z) + |z|^{-N}\,E_{s\to w}\,K(z)\,\bigl(a^N E_{k,w\to s} + b_N\bigr)} {1 - \widetilde{\beta}_N(z)},\]

\[\widetilde{\beta}_N(z) \;:=\; |z|^{-N}\,\Delta_k\,K(z)\, \bigl(a^N E_{k,w\to s} + b_N\bigr), \qquad S_N^{(k)}(z) := \frac{1}{|z|}\sum_{\ell=0}^{N-1}\frac{\|L_k^{\ell}\|_w}{|z|^{\ell}}.\]

Caller-supplied inputs

• abs_z — magnitude $|z|$; must satisfy $|z| > a$ for a useful bound.
• K_z — upper bound on the strong resolvent $\mathcal{R}_s(z, L)$ of the infinite-dimensional $L$ at this $z$ (typically from strong_resolvent_lift at a coarse level).
• N — Laurent-truncation parameter, ≥ 1.

Basis-derived defaults (override via keyword)

Internally b is multiplied by M_2 = aux_weak_bound(B) before being fed into b_N, so the paper formula

\[b_N = b \sum_{j=0}^{N-1} a^{N-1-j} M^j\]

is applied with b ← b · M_2 and M ← M_aux. This is the iterated (strong, weak) DFLY constant derived from the package's (strong, aux) DFLY plus the aux-to-weak conversion — and stays bounded by b · M_2 / (1 - a) when M_aux = 1.

Returns Inf when the perturbation condition $\widetilde{\beta}_N(z) < 1$ fails (typically because $N$ is too small or the discretization is too coarse).

Example

B = Ulam(1024)
D = mod1_dynamic(x -> 2 * x)
norms = powernormbounds(B, D)              # ‖Q^ℓ|_U‖_w sequence
M_k_seq = max.(1.0, norms)                 # whole-space upper bounds

# K_z certified at a coarse level via strong_resolvent_lift / a Schur
# resolvent bound on a contour.
R_w_fine = coarse_fine_weak_resolvent(B, D, abs_z, K_z;
                                       N = 8, M_k_weak = M_k_seq)

coarse_fine_weak_resolvent_auto_N(B::Basis, D::Dynamic, abs_z, K_z;
                                   μ=nothing, kwargs...)

Convenience wrapper applying Corollary A.16: choose

\[N_k := \left\lceil\frac{\log E_{k,w\to s}}{|\log(a/\mu)|}\right\rceil,\]

where $\mu \in (a, M)$ is a target spectral radius (defaulting to the geometric mean $\sqrt{aM}$ when not supplied). This is the truncation depth that balances $a^{N_k}$ against $(\mu/M)^{N_k}$ in the bound.

Returns (R_w_fine, N_k). Remaining keywords are forwarded to coarse_fine_weak_resolvent and share its basis-interface defaults.

projector_distance_bound(B_coarse, B_fine, D, Q_coarse, Q_fine, radius;
                          eigenvalue = 1.0, samples = 256, N, kwargs...)

Float64 upper bound on the Riesz projector distance $\|P_L - P_{L_{\rm fine}}\|_{s \to w}$ for the spectral projector onto the eigenspace of $L$ and $L_{\rm fine}$ enclosed by the small circle

\[\Gamma = \{ z \in \mathbb{C} : |z - \mathrm{eigenvalue}| = r \}.\]

Combines Lemma A.11 (strong-to-weak resolvent perturbation) with the contour-integral representation of the Riesz projector:

\[\| P_L - P_{L_{\rm fine}} \|_{s \to w} \;\le\; \frac{\ell(\Gamma)}{2\pi} \sup_{z \in \Gamma} \mathcal{R}_w(z, L_{\rm fine}) \,\delta_{\rm fine}\, \mathcal{R}_s(z, L) \;=\; r \cdot \sup_{z \in \Gamma} \mathcal{R}_w(z, L_{\rm fine}) \,\delta_{\rm fine}\, K(z).\]

How the pieces are obtained

The strategy is the coarse-fine one: validated certification is done on the coarse matrix only (small enough for cheap Schur), and both fine-level quantities are then propagated.

• CertifScripts.run_certification on the coarse Q_coarse.L with a CertifScripts.CertificationCircle of radius r centred at eigenvalue gives $\sup_{z \in \Gamma} \mathcal{R}_w(z, L_{\rm coarse})$.
• $K(z) := \mathcal{R}_s(z, L)$ is obtained from the coarse weak resolvent via strong_resolvent_lift (Prop A.7), evaluated at the worst contour magnitude $|z|_{\min} = |\text{eigenvalue}| - r$.
• $\sup_{z \in \Gamma} \mathcal{R}_w(z, L_{\rm fine})$ is propagated to the fine level via coarse_fine_weak_resolvent (Prop A.14), also at the worst contour magnitude.
• $\delta_{\rm fine}$ is weak_projection_error(B_fine).

This way CertifScripts is run only on the coarse matrix; the fine level is reached entirely by the propagation pipeline.

Inputs

• B_coarse, B_fine — coarse and fine Basis (same family).
• D::Dynamic — the dynamic.
• Q_coarse, Q_fine — DiscretizedOperator(B, D) for each level.
• radius — circle radius r; must satisfy 0 < r < |eigenvalue| and in practice r < 1 − λ₂ (else Γ encloses extra spectrum).
• eigenvalue::Number — circle centre (default 1.0).
• samples::Integer — number of contour samples for CertifScripts (default 256).
• N::Integer — Laurent-truncation depth for Prop A.14.
• M_k_weak — bound on $\|L_{\rm fine}\|_w$ (default upper_bound_L2_opnorm(BallMatrix(Q_fine.L))).
• Remaining kwargs (M_aux, M_2, Δ_k, E_sw, E_k_ws, dfly_coefficients, schur_data) are forwarded — see coarse_fine_weak_resolvent.

Connection to the eigenvector distance

For an L¹-preserving transfer operator with simple eigenvalue 1, both Riesz projectors are rank-1 of the form $P f = v\,\langle\ell, f\rangle$ where $\ell$ is the integral functional. Hence

\[\|P_L - P_{L_{\rm fine}}\|_{s \to w} \;=\; \|v - v_{\rm fine}\|_w \cdot \|\ell\|_{s \to \mathbb{R}},\]

so the projector-distance bound translates directly into a bound on the invariant-measure distance up to the factor $\|\ell\|_{s\to\mathbb R}$ (typically ≤ 1 for the integral functional on regular spaces).

Returns a named tuple with projector_distance, R_w_coarse, K_z, R_w_fine, δ_k, radius, and the underlying cert/schur_data. If Prop A.7 or A.14 fails to close, projector_distance is Inf.

strong_resolvent_lift(B::Basis, D::Dynamic, abs_z, R_w_coarse; kwargs...)

Float64 upper bound on the strong resolvent of the infinite-dimensional operator,

\[\mathcal{R}_s(z, L) \;\le\; \frac{b\,\mathcal{R}_w(z, L_{k_0})\,E_{s\to w} + 1} {|z| - a - b\,\mathcal{R}_w(z, L_{k_0})\,\delta_{k_0}},\]

derived from Proposition A.7 of the paper. The numerator bounds $\|u\|_s$ for the unique solution of $(zI - L)u = f$ with $\|f\|_s = 1$; the denominator is the perturbation gap

\[|z| > a + b\,\mathcal{R}_w(z, L_{k_0})\,\delta_{k_0}.\]

Use this to lift a coarse-level weak resolvent certificate $R_{w,\text{coarse}} \ge \mathcal{R}_w(z, L_{k_0})$ to a strong resolvent bound for the infinite-dimensional $L$. Returns Inf when the perturbation gap fails to close.

Standing constants are pulled from the basis interface; override via keyword for sharper specific-case bounds: