AverageZero{B<:Basis}

Yield a basis of the space of average zero vectors

Replacement of DualComposedWithDynamic.

aux_normalized_projection_error(B::Basis)

Return a constant Eh (typically scales as h ~ 1/n) such that

$|||P_h f|||\leq |||f|||+ Eh * ||f||_s$

Must be rounded up correctly!

aux_weak_bound(B::Basis)

Return a constant $M₂$ such that for a vector $v ∈ Uₕ$

$|||v|||\leq M_2||v||$

Must be rounded up correctly!

bound_linalg_norm_L1_from_weak(B::Basis)

Return a constant $A$ such that for a vector $v ∈ Uₕ$

$||v||_1\leq A||v||$

Must be rounded up correctly!

bound_linalg_norm_L∞_from_weak(B::Basis)

Return a constant $A$ such that for a vector $v ∈ Uₕ$ ||v||_\infty \leq A||v|| Must be rounded up correctly!

bound_weak_norm_abstract(B::Basis, D=nothing; dfly_coefficients=nothing)

Returns an a priori bound on the weak norm of the abstract operator $L$

bound_weak_norm_from_linalg_norm(B::Basis)

Return constants W₁, W₂ such that for a vector $v ∈ Uₕ$

$||v||\leq W_1||v||_1+W_2||v||_{\infty}$

Must be rounded up correctly!

evaluate(B::Basis, i, x)

Evaluate the i-th basis element at x

evaluate_integral(B::Basis, i; T = Float64)

Value of the integral on [0,1] of the i-th basis element

integral(B::Basis, v; T = Float64)

Return the integral of the function with coefficients v in the basis B

integral_covector(B::Basis)

Return a covector that represents the integral in the basis B

invariant_measure_strong_norm_bound(B::Basis, D::Dynamic)

Bounds $||u||_s$, where $u$ is the invariant measure normalized with $i(u)=1$.

is_integral_preserving(B::Basis)

Integral-preserving discretizations may specialize this to "true"

is_refinement(Bfine::Basis, Bcoarse::Basis)

Check if Bfine is a refinement of Bcoarse

Example

julia> using RigorousInvariantMeasures

julia> B = Ulam(1024)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 1025))

julia> Bfine = Ulam(2048)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 2049))

julia> is_refinement(Bfine, B)
true

julia> Bfine = Ulam(2049)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 2050))

julia> is_refinement(Bfine, B)
false

one_vector(B::Basis)

Vector that represents the function 1 in the basis B

projection(B::Basis, f::Function; kwargs...)

Project a function f onto the basis B. The result type and available keyword arguments are basis-dependent. Currently implemented via the TaylorModels extension for the Ulam basis.

strong_norm(B::Basis)

Return the type of the strong norm of the basis

Example

julia> using RigorousInvariantMeasures

julia> B = Ulam(1024)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 1025))

julia> strong_norm(B)
TotalVariation

strong_weak_bound(B::Basis)

Return a constant $M₁n$ such that for a vector $v ∈ Uₕ$

$||v||_s\leq M1n*||v||$

Must be rounded up correctly!

Return constants $S₁, S₂$ such that for a vector $v ∈ Uₕ$

$||v||\leq S_1||v||_s+S_2|||v|||$

Must be rounded up correctly!

weak_norm(B::Basis)

Return the type of the weak norm of the basis

Example

julia> using RigorousInvariantMeasures

julia> B = Ulam(1024)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 1025))

julia> weak_norm(B)
L1

weak_projection_error(B::Basis)

Return a constant Kh (typically scales as h ~ 1/n) such that

$||P_h f-f||\leq Kh ||f||_s$

Must be rounded up correctly! This function is not exported explictly but is used in all the estimates.

Example

julia> using RigorousInvariantMeasures;

julia> B = Ulam(1024)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 1025))

julia> RigorousInvariantMeasures.weak_projection_error(B)
0.00048828125

Ulam

Ulam basis on [0,1] associated to the partition $p = \{x_0 = 0, x_1, \ldots, x_n=1\}$

Ulam(n::Integer)

Equispaced Ulam basis on [0,1] of size n

p1 * p2  for two `ProjectedFunction{<:Ulam}` on the same basis

Componentwise multiplication of the cell-value vectors. For Ulam each $p_i.v$ is the cell-value vector of the piecewise-constant reconstruction; their pointwise product is again piecewise-constant on the same partition with cell value $p_1.v[i] \cdot p_2.v[i]$.

Bounds combine via Hölder:

\[\|fg - φ_{p_1}φ_{p_2}\|_{L^1} \;\leq\; \|f\|_{L^∞}\,\|g - φ_{p_2}\|_{L^1} + \|φ_{p_2}\|_{L^∞}\,\|f - φ_{p_1}\|_{L^1}\]

so

weak_dual_bound = p1.weak_dual_bound * p2.weak_dual_bound
proj_error =
    p1.weak_dual_bound * p2.proj_error
    + p2.weak_dual_bound * p1.proj_error

Base.getindex(B::Ulam, i::Int)

Returns the i-th element of the Ulam basis as a function.

Example

julia> using RigorousInvariantMeasures

julia> B = Ulam(16)
Ulam{LinRange{Float64, Int64}}(LinRange{Float64}(0.0, 1.0, 17))

julia> B[1](1/32)
1

julia> B[2](1/32)
0

iterate(S::AverageZero{Ulam{T}}, state = 1) where{T}

Return the elements of the basis of average $0$ functions in the Ulam Basis

Given a preimage of an interval I_i, this iterator returns its relative intersection with all the elements of the Ulam basis that have nonzero intersection with it

Base.length(B::Ulam)

Returns the size of the Ulam basis (the size of the underlying vector -1)

Base.length(S::AverageZero{Ulam})

Return the size of the Average Zero space

integral_pairing(ϕ::Observable{<:Ulam}, ρ::AbstractVector, ρ_w_error;
                 ρ_dual_weak_bound = …)

For Ulam the pairing simplifies: $ρ_N$ is piecewise constant and $ϕ.v[i] = N\,\int_{I_i} ϕ\,dx$ is N times the exact cell integral, so

\[\frac{1}{N}\sum_i ϕ.v[i]\,ρ_N[i] \;=\; \sum_i ρ_N[i]\,\int_{I_i} ϕ\,dx \;=\; \int_0^1 ϕ\,ρ_N\,dx.\]

There is no projection-error term $\|ϕ - ϕ_N\|_w\,\|ρ_N\|_{w^*}$ to add here — the discrete pairing already equals $\int ϕ\,ρ_N$ exactly (modulo floating-point rounding, which interval lifting of ρ covers). The only error contribution is from $ρ \neq ρ_N$: $|\int ϕ\,(ρ - ρ_N)\,dx| \leq \|ϕ\|_{L^∞}\,\|ρ - ρ_N\|_{L^1}$.

The ρ_dual_weak_bound kwarg is accepted for API symmetry with the Fourier method but is unused for Ulam.

nonzero_on(B::Ulam, (a, b))

Returns the indices of the elements of the Ulam basis that intersect with the interval y We do not assume an order of a and b; this should not matter unless the preimages are computed with very low accuracy. We assume, though, that y comes from the (possibly inexact) numerical approximation of an interval in $[0,1]$, i.e., we restrict to $y \cap [0,1]$

normbound(B::Ulam, ::Type{TotalVariation}, v)

Rigorous upper bound for the total variation of the piecewise-constant function on $[0,1]$ with Ulam coefficients v (its variation is exactly $\sum_j |v_{j+1}-v_j|$; the sum is evaluated in interval arithmetic).

projection_defect_coefficients(B::Ulam, D::Dynamic; dfly_coefficients)

Sharp Ulam discretization-defect constants $(c_s, c_w) = (1 + A,\; B)$, i.e.

\[\|(L - Q_h)\,v\|_{L^1} \le K_h\big[(1+A)\,\mathrm{Var}(v) + B\|v\|_{L^1}\big], \qquad K_h = 1/(2n).\]

Proof. For Ulam the projection $\Pi$ is the conditional expectation on the cells, so $\|\Pi\|_{L^1} \le 1$, $\|(I-\Pi)g\|_{L^1} \le K_h \mathrm{Var}(g)$, the discretized operator is integral preserving ($Q_h = \Pi L \Pi$, no correction term), and $\|L\|_{L^1} \le 1$. Decomposing $L - Q_h = (I-\Pi)L + \Pi L(I-\Pi)$ (on the range of $\Pi$ the second term even vanishes):

\[\|(I-\Pi)Lv\|_{L^1} \le K_h \mathrm{Var}(Lv) \le K_h(A\,\mathrm{Var}(v) + B\|v\|_{L^1}), \qquad \|\Pi L(I-\Pi)v\|_{L^1} \le \|(I-\Pi)v\|_{L^1} \le K_h \mathrm{Var}(v). \qquad\square\]

This improves the generic compatible-discretization constants $(2(\|L\|_w+A), 2B)$ of Lemma 3.5 in [Galatolo–Monge–Nisoli–Poloni, Chaos Solitons & Fractals 170 (2023) 113329] by roughly a factor $4$.

relative_measure((a,b)::Tuple{<:Interval,<:Interval}, (c,d)::Tuple{<:Interval,<:Interval})

Relative measure of the intersection of (a,b) wrt the whole interval (c,d) Assumes that a,b and c,d are sorted correctly

Hat

This type represents a Hat basis on $S^1$. It contains a vector with the midpoints of the hats.

Hat(n::Integer)

This constructs a Hat basis on $S^1$ on equispaced points

HatFunctionOnTorus

Hat function (on the torus)

This is a piecewise linear function such that: f(x) = 0 if x <= lo f(mi) = 1 f(x) = 0 if x >= ho

Evaluate a HatFunctionOnTorus correctly on an IntervalOnTorus

Assumption: this is only called on functions defined on our partition, so either mi==0, hi==0, or the three values are in increasing order

IntervalOnTorus

A separate type for intervals on the torus (mod 1) to "remind" us of the quotient

The interval is normalized in the constructor: the caller may assume that

• 0 <= inf(i) < 1
• sup(i) < inf(i) + 1 OR i==interval(0,1)

Base.getindex(B::Hat, i::Int)

Make so that B[j] returns a HatFunctionOnTorus with the j-th basis element

Given a preimage $y$ of a point $x$, this iterator returns $\phi_j(y)/T'(y)$

Base.length(B::Hat{T})

Return the size of the Hat basis

nonzero_on(B::Hat, dual_element)

Return the range of indices of the elements of the basis whose support intersects with the given dual element (i.e., a pair (y, absT')). The range may end with length(B)+1; this must be interpreted "mod length(B)": it means that it intersects with the hat function peaked in 0 as well (think for instance y = 0.9999).

normbound(B::Hat, ::Type{L1}, v)

Rigorous upper bound for the $L^1$ norm (the auxiliary norm of the Hat basis) of the continuous piecewise-linear function on the torus with Hat coefficients v. On each piece $[p_i,p_{i+1}]$ the trapezoidal value $(p_{i+1}-p_i)\,(|v_i|+|v_{i+1}|)/2$ is an upper bound for $\int|f|$ (exact when $f$ does not change sign, e.g. for a density); summing over the $n$ pieces, wrap-around included, bounds $\|f\|_{L^1}$. Evaluated in interval arithmetic.

normbound(B::Hat, ::Type{Lipschitz}, v)

Rigorous upper bound for the Lipschitz seminorm (the strong norm of the Hat basis) of the continuous piecewise-linear function on the torus $S^1$ with Hat coefficients v (its values at the nodes $p_1,\dots,p_n$). On a linear piece the slope is constant, so the seminorm is exactly the largest absolute slope over the $n$ pieces, including the wrap-around piece joining the last node to the first:

\[\mathrm{Lip}(f) = \max_i \frac{|v_{i+1}-v_i|}{p_{i+1}-p_i},\]

with $v_{n+1} := v_1$ and $p_{n+1}-p_n$ the length of the wrap-around piece. The quotients are evaluated in interval arithmetic.

Hat function (on the reals)

This is a piecewise linear function such that: f(x) = 0 if x <= lo f(mi) = 1 f(hi)

Evaluate a HatFunction (on the real line).

Must work correctly when S is an interval.

Base.getindex(B::HatNP, i::Int)

makes so that B[j] returns a HatFunction with the j-th basis element

Given a preimage y of a point x, this iterator returns \phi_j(y)/T'(y)

Base.length(B::HatNP)

Return the size of the HatNP basis

Return the range of indices of the elements of the basis whose support intersects with the given dual element (i.e., a pair (y, absT')). The range may end with length(B)+1; this must be interpreted "mod length(B)": it means that it intersects with the hat function peaked in 0 as well (think for instance y = 0.9999).

Equispaced partition of size n of [0,1]

Return (in an iterator) the pairs (i, (x, |T'(x)|)) where x is a preimage of p[i], which describe the "dual" L* evaluation(p[i])

Given a preimage y of a point x, this iterator returns $\phi_j(y)/T'(y)$

Return the size of the C2 basisBase.length(S::AverageZero) = length(S.basis)-1

Return the range of indices of the elements of the basis whose support intersects with the given dual element (i.e., a pair (y, absT', derder)).

Chebyshev{S<:NormKind, WK<:NormKind, T<:AbstractVector} <: Basis

Chebyshev basis $φ_i(x) = T_{i-1}(2x-1)$ on $[0,1]$, carrying its strong and weak norms as type parameters in the same style as FourierAnalytic.

Two weak norms are available:

• C1 (the default, and what Chebyshev(n, k) builds) — the setting of Nisoli & Taylor-Crush, Rigorous Computation of Linear Response for Intermittent Maps, J. Stat. Phys. 190 (2023) 192, whose Theorems 3.13/3.14 supply the projection errors used below.
• L2, selected by Chebyshev(n, k, L2). The L2 here is $L^2(μ)$ for the arcsine measure $dμ = dx/(π\sqrt{x(1-x)})$, the one that makes the $T_m$ orthogonal, so Parseval holds on the coefficients and the norm interface collapses to the Fourier-like one-liners. Use l2_measure_conversion_bounds and gram_restrict_to_average_zero to move between this and $L^2(dx)$.

The strong norm is $W^{k,1}$ in both cases, so dfly(W{k,1}, L1, D) applies unchanged — the Lasota–Yorke inequality is a statement about function spaces and does not see the basis.

Base.getindex(B::Chebyshev, i::Int)

Make so that B[j] returns a HatFunctionOnTorus with the j-th basis element

Return the size of the Chebyshev basis

weak_projection_error(B::Chebyshev{Eρ, L2})

$\|f - p_n\|_{L^2(dx)} \le 4ρ^{-n}/(ρ-1)$ on the unit ball of the Eρ norm, from Trefethen, ATAP, Thm 8.2 (8.3) — the interpolation form, since this basis interpolates at the Chebyshev points. dx is a probability measure on $[0,1]$, so the sup-norm bound carries to $L^2$ with constant 1.

_cheb_gram_factor(n, T; precision_bits = _default_chol_bits(T))

Cholesky factor U of the Lebesgue Gram matrix and its verified inverse, cached by (n, T, precision_bits).

_eval_T

Eval the Chebyshev polynomial up to degree n on an array of points in [-1, 1].

Not satisfactory, the intervals explode

_verified_inverse(G::Matrix{Interval{T}}) where {T}

Verified inverse of a general interval matrix, by one Krawczyk-verified linear solve per column. The Lebesgue Gram matrix is only mildly ill-conditioned (cond ~ 1.3n measured), so this converges comfortably.

For a triangular matrix prefer _verified_inverse_upper_triangular, which is far cheaper and just as tight.

_verified_inverse_upper_triangular(U)

Rigorous enclosure of $U^{-1}$ for an upper-triangular interval matrix, by back-substitution: column j solves $Ux = e_j$, whose entries above j vanish, so the work is $O(n^3/6)$ and every operation is a single interval divide or fused sum — no linear-system verification is needed, because back-substitution is exact as a formula and interval arithmetic carries the rounding.

Use this instead of _verified_inverse whenever the matrix is a Cholesky factor, which is always the case here. The general Krawczyk route costs one verified solve per column and is enormously more expensive at high precision: at n = 129 and 333 bits it takes 51.12 s against 0.37 s for this routine — a 138× difference — while returning the identical enclosure (max radius 1.3959819931627471e-75 both ways, $\|UV - I\| = 3.63\cdot10^{-74}$ both ways). It was the dominant cost of the whole BigFloat pipeline.

bernstein_expansion(f, ρ; n = 1024) -> ρ_image

Rigorous lower bound on

\[\min_{z \in ∂E_ρ} \; \texttt{bernstein\_parameter}(f(z)) ,\]

obtained by covering the boundary parameter with n intervals, enclosing the image of each arc under f, and taking the smallest enclosure. f must accept a Complex{Interval}.

f maps $E_ρ$ strictly outside itself — it expands the ellipse — exactly when the returned value exceeds ρ; see expands_bernstein_ellipse. This is the Chebyshev analogue of enclosing the image of an annulus and reading off its inner and outer radii, as done for the Fourier/analytic setting; an ellipse needs only the one number because $∂E_ρ$ is a single curve.

Example

The Chebyshev polynomials are exactly the maps $T_m(E_ρ) = E_{ρ^m}$, so bernstein_expansion(z -> 2z^2 - 1, ρ) returns $ρ^2$.

bernstein_l1_l2_equivalence(B::Chebyshev{Eρ}) -> √n

The constant in $\|f\|_{E_ρ} \le \sqrt n\,\|f\|_{E_ρ^{(2)}}$ on the span of the basis (Cauchy–Schwarz, n coefficients). The reverse inequality holds with constant 1.

bernstein_l2_resolvent_bound(B::Chebyshev, block::BallMatrix, z = 1.0; T = Float64)

Rigorous bound on $\|(zI - Q_N)^{-1}\|$ in the $ℓ^2$ Bernstein norm, for block the mean-zero restriction returned by gram_restrict_to_average_zero.

block lives in the Gram coordinates $y = Uc$, so the Bernstein weight has to be pulled back: a mean-zero vector is $c = U^{-1}[0; y']$ and its norm is $\|W y'\|_2$ with $W = (D U^{-1})[:, 2:\mathrm{end}]$. Writing $S = W^{*}W = R^{*}R$ for the Cholesky factor R, the induced norm is $\|R\,M\,R^{-1}\|_2$, so the bound is the verified SVD of $R (zI - Q_N) R^{-1}$ inverted.

bernstein_parameter(z)

The parameter $ρ ≥ 1$ of the Bernstein ellipse through the point z.

$E_ρ$ has foci $\pm 1$ and major axis $ρ + ρ^{-1}$, so the focal-distance characterisation of the ellipse gives

\[s := |z-1| + |z+1| = ρ + ρ^{-1}, \qquad ρ = \frac{s + \sqrt{s^2-4}}{2}.\]

Working through s keeps everything real: no complex square root, and hence no branch cut to worry about (the two Joukowski preimages $w$ and $w^{-1}$ of z give the same answer by construction). z may be a complex interval, and the result then encloses the parameters of all points it contains.

bernstein_point(ρ, θ)

The point of the Bernstein ellipse $∂E_ρ$ at parameter θ ∈ [0,1], namely $z = (w + w^{-1})/2$ with $w = ρ\,e^{2πiθ}$, which expands to

\[z = \frac{ρ + ρ^{-1}}{2}\cos 2πθ \;+\; i\,\frac{ρ - ρ^{-1}}{2}\sin 2πθ .\]

ρ and θ may be intervals, in which case the result encloses the corresponding arc.

bound_linalg_norm_L1_from_weak(B::Chebyshev)

Returns A such that ||ĉ||_{ℓ¹} ≤ A·||f||_{C1}.

Chebyshev coefficients satisfy |cⱼ| ≤ 2·||f||∞ for j ≥ 1 and |c₀| ≤ ||f||∞, so ||ĉ||{ℓ¹} ≤ (2n-1)·||f||∞ ≤ (2n-1)·||f||_{C1} where n = length(B).

bound_linalg_norm_L∞_from_weak(B::Chebyshev)

Returns A such that ||ĉ||_{ℓ∞} ≤ A·||f||_{C1}.

maxj |cⱼ| ≤ 2·||f||∞ ≤ 2·||f||_{C1}.

bound_weak_norm_abstract(B::Chebyshev, D; dfly_coefficients)

A priori bound on $\|L\|_{L^2 \to L^2}$.

Cauchy–Schwarz against the measure $\sum_k |g_k'|\,δ_{g_k(x)}$ gives $|Lf|^2 \le (L\mathbf 1)(L|f|^2)$; integrating and using $\int L|f|^2 = \int|f|^2$ yields

\[\|L\|_{L^2\to L^2} \;\le\; \Big(\sup_{[0,1]} L\mathbf 1\Big)^{1/2} \;\le\; \Big(\sum_k \frac{1}{\min_{[X_1,X_2]}|T_k'|}\Big)^{1/2},\]

which is what we compute when the dynamic is available. For the Lanford map this gives 1.0732 against a measured $\|Q\|_{L^2} = 1.0716$.

The package-wide convention elsewhere (Ulam, Hat, Fourier) is B + 1 from the Lasota–Yorke coefficients. That is fine when B > 0, but it degenerates to exactly 1 for the analytic strong norms, where B = 0 by construction — and 1 is not an upper bound for $\|L\|_{L^2}$ unless $\sup L\mathbf 1 \le 1$. Hence the direct computation here. We fall back to B + 1 only when no dynamic is supplied or a branch derivative fails to be bounded away from zero.

bound_weak_norm_from_linalg_norm(B::Chebyshev)

Returns (W₁, W₂) such that ||f||_{C1} ≤ W₁·||ĉ||_{ℓ¹} + W₂·||ĉ||_{ℓ∞}.

Since |Tⱼ(x)| ≤ 1: ||f||∞ ≤ ||ĉ||{ℓ¹}. For the derivative: ||f'||∞ ≤ Σ|cⱼ|·2j² ≤ 2(n-1)²·||ĉ||{ℓ¹} where n = degree.

certify_spectral_gap(B::Chebyshev, Q::DiscretizedOperator;
                      samples=256, radius_factor=0.5)

Run the full spectral gap certification pipeline for a Chebyshev discretized operator.

Uses BallArithmetic's CertifScripts to:

1. Compute Schur decomposition with rigorous error bounds
2. Find the spectral gap (distance from eigenvalue 1 to next largest eigenvalue)
3. Certify the resolvent on a circle separating eigenvalue 1 from the rest
4. Compute projector error bound δP for `restricttoaveragezero`

Returns a named tuple with fields:

• spectral_gap, second_largest, ρ (contour radius)
• M_inf (certified resolvent bound on contour)
• projector_error (δ_P for projector inflation)
• certification, schur_data, eigenvalue_index

eval_Clenshaw_BackwardFirst

Eval a polynomial in Chebyshev basis, ClenshawBackward, using ball arithmetic Following Viviane Ledoux, Guillaume Moroz "Evaluation of Chebyshev polynomials on intervals andapplication to root finding"

expands_bernstein_ellipse(f, ρ; n = 1024) -> (expands, ρ_image)

Whether f maps $∂E_ρ$ strictly outside $E_ρ$, together with the certified image parameter from bernstein_expansion.

For an interval map this is the analytic expansion condition behind an Eρ Lasota–Yorke inequality: if the forward map expands the ellipse then its inverse branches contract into it, so the transfer operator preserves analyticity on $E_ρ$.

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

gram_matrix(B::Chebyshev; T = Float64)
inv_gram_matrix(B::Chebyshev; T = Float64)

Gram matrix $G_{ij} = \langle φ_i, φ_j \rangle$ of the Chebyshev basis $φ_i(x) = T_{i-1}(2x-1)$ on $[0,1]$, and its inverse.

The inner product is taken against the arcsine probability measure

\[dμ(x) = \frac{dx}{π\sqrt{x(1-x)}}, \qquad \int_0^1 dμ = 1,\]

the measure that makes the Chebyshev polynomials orthogonal. Substituting $t = 2x-1$ turns it into the familiar $dt/(π\sqrt{1-t^2})$, so

\[G = \operatorname{diag}(1, \tfrac12, \tfrac12, \dots, \tfrac12), \qquad G^{-1} = \operatorname{diag}(1, 2, 2, \dots, 2).\]

(Against the unnormalized weight $1/\sqrt{1-t^2}$ every entry is $π$ times these.) Both are returned as Diagonal of intervals; every entry is exactly representable, so the enclosures are thin.

Transporting an $\ell^2$ bound to the function space

For $f = \sum_i c_i φ_i$ we have $\|f\|_{L^2(μ)}^2 = c^* G c$, so an operator with coefficient matrix Q satisfies

\[\|Q\|_{L^2(μ) \to L^2(μ)} = \|G^{1/2}\,Q\,G^{-1/2}\|_{\ell^2}.\]

This is the point of using this weight rather than Lebesgue: G is diagonal, so the similarity is a cheap rescaling and BallArithmetic's $\ell^2$ estimators (upper_bound_L2_opnorm, svd_bound_L2_opnorm_inverse, …) apply to gram_sqrt(B) * Q * inv_gram_sqrt(B) directly. Under Lebesgue the Gram matrix is dense and ill-conditioned and no such shortcut exists.

See also gram_sqrt, inv_gram_sqrt.

gram_restrict_to_average_zero(B::Chebyshev, BM::BallMatrix; T = Float64)

Restrict BM to the average-zero subspace by the Gram change of variables, returning (block, chol) where block is an (n-1) × (n-1) BallMatrix.

Unlike restrict_to_average_zero, which builds a certified Riesz projector via a Schur decomposition, this uses the fact that in the Gram-transformed coordinates the restriction is a plain submatrix.

Writing $G$ for the Lebesgue Gram matrix (gram_matrix with measure = :lebesgue) and $G = U^{*}U$ for its Cholesky factor:

• $T_0 = 1$, so the constant function is $e_1$ and the integral covector is exactly $v = G e_1$ — the first column of $G$;
• $U$ is upper triangular, hence $U e_1 \parallel e_1$ and $U^{-*} v = U_{11} e_1$, so average-zero $\{v \cdot c = 0\}$ becomes $\{y_1 = 0\}$ in the coordinates $y = Uc$;
• since the transfer operator preserves the integral, $v^{*}Q = v^{*}$, and therefore $\tilde A = U Q U^{-1}$ has first row exactly $e_1^{*}$.

The restriction is then $\tilde A[2:\mathrm{end}, 2:\mathrm{end}]$ — the same shape as restrict_to_average_zero(B::Fourier, …), of which this is the special case $G = I$. Because $\|M\|_G = \|U M U^{-1}\|_2$ exactly, an $\ell^2$ bound on the returned block is the $L^2(dx)$ bound on the average-zero subspace, with no condition-number penalty.

gram_sqrt(B::Chebyshev; T = Float64)
inv_gram_sqrt(B::Chebyshev; T = Float64)

The symmetric factors $G^{1/2}$ and $G^{-1/2}$ of the Chebyshev gram_matrix — i.e. $\operatorname{diag}(1, 1/\sqrt2, \dots)$ and $\operatorname{diag}(1, \sqrt2, \dots)$.

These, not G itself, are what an $\ell^2$ operator-norm estimator needs: $\|Q\|_{L^2(μ)} = \|G^{1/2} Q G^{-1/2}\|_{\ell^2}$. Unlike G, the entries are irrational, so the returned intervals are genuine (thin but not exact) enclosures.

gram_sqrt(B::Chebyshev; T = Float64)
inv_gram_sqrt(B::Chebyshev; T = Float64)

The symmetric factors $G^{1/2}$ and $G^{-1/2}$ of the Chebyshev gram_matrix — i.e. $\operatorname{diag}(1, 1/\sqrt2, \dots)$ and $\operatorname{diag}(1, \sqrt2, \dots)$.

These, not G itself, are what an $\ell^2$ operator-norm estimator needs: $\|Q\|_{L^2(μ)} = \|G^{1/2} Q G^{-1/2}\|_{\ell^2}$. Unlike G, the entries are irrational, so the returned intervals are genuine (thin but not exact) enclosures.

invariant_measure_strong_norm_bound(B::Chebyshev{W{k,l}}, D; dfly_coefficients)

The classical DFLY bound $\|h\|_s \le B/(1-A)$ on the invariant density, as for Ulam, Hat and Fourier. The analytic (Eρ) strong norm is degenerate, B = 0, and is handled separately in AnalyticDFLY.jl, where the bound is A.

Without this method distance_from_invariant fails with a MethodError on any W^{k,1} Chebyshev basis.

is_refinement(Bf::Chebyshev, Bc::Chebyshev)

Whether Bf refines Bc, i.e. whether the Chebyshev points of Bc are a subset of those of Bf. With degrees $n_f$ and $n_c$ the points are $\cos(jπ/n)$, so nesting holds exactly when $n_c \mid n_f$.

Note the argument order: the contract of is_refinement is (fine, coarse), as for Ulam and HatNP. This method used to be written (Bc, Bf) with the body length(Bc) < length(Bf), so a correct is_refinement(fine, coarse) call returned false and norms_of_powers_from_coarser_grid logged "The fine basis is not a refinement of the coarse basis" on every coarse–fine run.

l2_measure_conversion_bounds(B::Chebyshev; T = Float64) -> (c_leb, C_n)

The two constants relating the Lebesgue and arcsine $L^2$ norms on the span of this basis:

\[\|f\|_{L^2(dx)} \le c_{\mathrm{leb}}\,\|f\|_{L^2(μ)}, \qquad \|f\|_{L^2(μ)} \le C_n\,\|f\|_{L^2(dx)}.\]

The first is uniform and needs no linear algebra: $dx/dμ = π\sqrt{x(1-x)}$ is bounded by $π/2$, so $c_{\mathrm{leb}} = \sqrt{π/2}$ on all of $L^2(μ)$.

The second cannot hold uniformly — $dμ/dx$ blows up at the endpoints — and is genuinely finite-dimensional: $C_n^2 = λ_{\max}(G_μ^{1/2} G_L^{-1} G_μ^{1/2})$. It is evaluated sharply here, from the verified inverse and the diagonal $G_μ^{1/2}$. Bounding it instead by $\|G_L^{-1}\|_2$ — valid since $\|G_μ^{1/2}\| = 1$, and cheaper in that it needs no inverse — costs a factor tending to $\sqrt2$, the extremal direction lying in the $\tfrac12$-eigenspace of $G_μ$. Measured growth of the sharp constant is $C_n \approx 0.75\sqrt{n}$.

Use these to carry a resolvent or spectral bound proved in $L^2(μ)$ — where Parseval makes it an $\ell^2$ statement about the coefficient matrix — back to the $L^2(dx)$ setting the rest of the package works in.

normbound(B::Chebyshev, ::Type{L2}, v)

$\|v\|_{L^2(dx)} = \sqrt{c^{*}G_L c}$ for the coefficient vector c, with G_L the Lebesgue Gram matrix — the basis is not orthonormal, so this is not ‖c‖_{ℓ²}.

opnormbound(B::Chebyshev, ::Type{L2}, M)

$\|M\|_{L^2(dx)} = \|U M U^{-1}\|_2$ with G_L = U^{*}U; exact, with no condition-number penalty. See gram_restrict_to_average_zero.

The conjugated matrix is measured with _l2_opnorm_ball, which prefers the verified-SVD enclosure over the cheap min(Collatz, √(‖·‖₁‖·‖_∞)) bound.

opnormbound(B::Chebyshev, ::Type{L2}, v::AbstractVector)

A column vector is the operator $\mathbb R \to U_h$, $t \mapsto t\,v$, whose operator norm is just $\|v\|_{L^2(dx)}$. This is the Q.e of a NonIntegralPreservingDiscretizedOperator.

opnormbound(B::Chebyshev, ::Type{L2}, w::Adjoint)

A covector is the operator $U_h \to \mathbb R$, $c \mapsto w^{*}c$, whose operator norm is the dual norm

\[\|w\|_{*} = \sup_{c\neq 0}\frac{|w^{*}c|}{\|c\|_{L^2(dx)}} = \sqrt{w^{*}G_L^{-1}w} = \|U^{-*}w\|_2 ,\]

with $G_L = U^{*}U$. Note the $G_L^{-1}$: measuring a covector with G_L, as if its entries were the coefficients of a function, is a different (and for an ill-conditioned Gram matrix a very different) quantity. This is the Q.w of a NonIntegralPreservingDiscretizedOperator — for Chebyshev that is the integral covector, which is genuinely not e₁.

restrict_to_average_zero(B::Chebyshev, BM::BallMatrix, f; certification=nothing)

Restrict BM to the average-zero subspace U⁰ using a certified spectral projector.

For Chebyshev, integral_covector(B) ≠ [1,0,...,0], so simple submatrix extraction does not work. Instead, we compute the Riesz projector P₁ onto the eigenspace of eigenvalue 1 (the invariant measure direction) via Schur decomposition, then return (I - P₁) * BM.

If certification is provided (from certify_spectral_gap), the projector radii are inflated by the certified projector error δ_P to account for the Schur factorization error ‖ZTZ* - A‖₂.

strong_weak_bound(B::Chebyshev{Eρ, L2})

$\|v\|_{E_ρ} \le M\,\|v\|_{L^2(dx)}$ on the span of the basis. On $E_ρ$ one has $|T_k| \le (ρ^k + ρ^{-k})/2 \le ρ^k$, so $\|v\|_{E_ρ} \le ρ^{n}\|\hat c\|_{ℓ^1}$, and bound_linalg_norm_L1_from_weak supplies $\|\hat c\|_{ℓ^1}$ in terms of the weak norm.

strong_weak_bound(B::Chebyshev{S, L2})

$\|v\|_{W^{k,1}} \le M \|v\|_{L^2}$ on the span of the basis, obtained by composing the existing C1 estimate with $\|v\|_{C^1} \le M' \|v\|_{L^2}$, where $M'$ comes from $\|v\|_∞ \le \|\hat c\|_{ℓ^1}$ and the Markov bound $\|v'\|_∞ \le 2(n-1)^2 \|\hat c\|_{ℓ^1}$ already used by bound_weak_norm_from_linalg_norm.

to_symmetric_interval(f)

Conjugate a map of $[0,1]$ into a map of $[-1,1]$, t ↦ 2f((t+1)/2) - 1.

Bernstein ellipses live around $[-1,1]$ while the dynamics in this package live on $[0,1]$, so this is what to wrap a branch in before handing it to bernstein_expansion.

weak_by_strong_and_aux_bound(B::Chebyshev)

Returns (S₁, S₂) such that ||f||_{C1} ≤ S₁·||f||_{W^{k,1}} + S₂·||f||_{L1}.

For k ≥ 2: Sobolev embedding gives ||f||∞ ≤ ||f||{L1} + ||f'||{L1} and ||f'||∞ ≤ ||f'||{L1} + ||f''||{L1}, so ||f||{C1} ≤ 2·||f||{W^{k,1}}.

For k = 1: uses Markov inequality ||p'||∞ ≤ 2n²·||p||∞ for polynomials of degree n on [0,1], giving ||f||{C1} ≤ (1 + 2n²)·||f||{W^{1,1}}.

weak_projection_error(B::Chebyshev{S, L2})

$L^2$ projection error for the $W^{ν,1}$ unit ball, $ν$ = B.k.

This follows from the $W^{k,1}$ decay of the Chebyshev coefficients, in two steps, both from Nisoli & Taylor-Crush:

• Theorem 3.12 (their statement of Trefethen, ATAP, Thm 7.1): if $f^{(ν)}$ has bounded variation $V$ then $|\hat b_m| \le 2V/(π\,m(m-1)\cdots(m-ν))$;
• Theorem 3.13: hence $\|f - π_n f\|_∞ \le 2V/(π ν\, n(n-1)\cdots(n+1-ν))$.

Both $dx$ and $dμ$ are probability measures on $[0,1]$, so $\|\cdot\|_{L^2} \le \|\cdot\|_∞$ with constant 1 and the same bound serves either measure with no conversion factor. It is the quantity already computed by aux_normalized_projection_error.

Theorem 3.13 is stated for the interpolant $π_n$, which is what this basis uses (the coefficients come from an FFT at the Chebyshev points), so the aliasing is already accounted for.

p1 * p2  for two `ProjectedFunction{<:Fourier}` on the same basis

Multiply two Fourier projections by treating each operand as a trigonometric polynomial (the discrete coefficient vector) plus an opaque L²-error bound. All bounds are derived from the finite coefficient vectors alone — provenance (the original W^{k,1} seminorm, function class, …) is not used. This keeps multiplication composable: the result is again a trigonometric polynomial with an L² error bound, ready to be passed to * or + again.

The arithmetic:

1. Convolve the two N-mode coefficient vectors → 2N−1 modes.
2. Truncate the central N modes back into the basis layout.
3. Parseval gives $\|\text{discarded modes}\|_{\ell^2}$ directly from the convolution output.

The bounds:

• $\|fg\|_{L^2} \leq \|\hat f\|_{\ell^1} \cdot \|g\|_{L^2}$ (Young's inequality, with $\|\hat f\|_{\ell^1}$ the surrogate for $\|f\|_{L^\infty}$).
• ``\|fg - πN(φf φg)\|{L^2} \leq \|\hat f\|{\ell^1}\,p2.\text{proj_error}
  • \|\hat g\|{\ell^1}\,p1.\text{proj_error} + \|\text{discarded modes}\|_{\ell^2}``.

Each input's proj_error is the L² distance from the original function to its trigonometric-polynomial approximant; the multiplication treats that distance as the only thing it knows.

Make so that B[j] returns the basis function of coordinate j

dual_nodes(B::Fourier, computed_dual) -> (x, labels, weights)

Flat description of a computed dual, used by the fast assembler.

eval_on_dual(B, computed_dual, ϕ_m) for the basis function $ϕ_m(x) = e^{2πimx}$ is, for every Fourier basis, of the form

\[w^{(m)}[\ell] = \sum_{j\,:\,\texttt{labels}[j] = \ell} \texttt{weights}[j]\; e^{2πi m x_j},\]

with real weights. Exposing the three vectors lets assemble_common sweep all frequencies at once instead of re-evaluating a rigorous complex exponential for every (basis function, node) pair; see the comment there.

Weights being real is what makes the $m \mapsto -m$ conjugate symmetry valid, so any new Fourier dual must preserve that.

integral_pairing(ϕ::Observable{<:Fourier}, ρ, ρ_w_error;
                 ρ_dual_weak_bound = weak_dual_norm_bound(ϕ.B, ρ))

For Fourier bases (weak L²), the pairing $\int_0^1 ϕ_N(x)\,ρ_N(x)\,dx$ equals $\sum_n \hat ϕ_n\,\overline{\hat ρ_n}$. The result is real-valued for real signals; we return the real part of the sum.

is_refinement(Bf::Fourier, Bc::Fourier)

Whether Bf refines Bc. The Fourier spaces are nested frequency truncations, $\mathcal U_{N_c} \subseteq \mathcal U_{N_f}$ iff $N_c \le N_f$.

Note the argument order: the contract of is_refinement is (fine, coarse), as for Ulam and HatNP. This method used to be written (Bc, Bf), so a correct is_refinement(fine, coarse) call returned false and norms_of_powers_from_coarser_grid logged "The fine basis is not a refinement of the coarse basis" on every coarse–fine run. It is also now ≤, matching Ulam, so a basis refines itself.

Return the size of the Fourier basis

FourierAnalytic(k::Integer, n::Integer; η = 0.1, T = Float64)
FourierAnalytic(p::AbstractVector, k::Integer)
FourierAnalytic(k_freq::Integer, n::Integer, ::Type{W{j,l}}; T = Float64)

Fourier basis of $2k+1$ frequencies, collocated at n points, for densities that are analytic on a strip. Coefficients are stored in FFT order, $[0:k; -k:-1]$.

The strong norm defaults to Aη(η) — the weighted $\ell^1$ norm on a strip of half-width $η$ — with L2 as the weak norm; this is the setting in which the analytic Lasota-Yorke inequality is proved. The third constructor swaps the strong norm for the Sobolev W{j,l} seminorm instead.

Assembly requires the FFTWExt extension: load FFTW alongside this package, or assemble(B, D) will raise a MethodError.

Return the weak-strong norm bound when restricted on the finite dimensional subspace. For Aη, by Cauchy-Schwarz: $||v||_{Aη} = \sum |ĉ_k| e^{2πη|k|} \leq (\sum e^{4πη|k|})^{1/2} \cdot ||ĉ||_{ℓ²}$ Geometric series: $\sum_{|k|\leq N} e^{4πη|k|} = 1 + 2(e^{4πη} - e^{4πη(N+1)})/(1 - e^{4πη})$

For W^{k,1}: Bernstein + Cauchy-Schwarz on [0,1]: $||v^{(j)}||_{L¹} ≤ ||v^{(j)}||_{L²} ≤ (2πN)^j ||v||_{L²}$ Sum: $M₁n = \sum_{j=0}^{k} (2πN)^j$