Noise kernel and noise-specialized estimates

DiscretizedNoiseKernelUlam{T,S}

Ulam discretization of a noise kernel whose bounds are of type T, one of the types Interval{T} is built on.

T is a parameter rather than fixed to Float64 because w and z are the scratch buffers the multiplication writes through: they are touched once per output entry, so at the abstract Vector they cost a boxed load every iteration, 2.2 MB per application at k = 16384; and fixed to Vector{Float64} they would silently narrow a Vector{Interval{BigFloat}} argument to double precision, which the error term added afterwards does not account for.

UniformNoiseUlam([T,] ξ, B::Ulam, boundarycondition = :periodic)

Ulam discretization of the uniform noise kernel of half-width ξ, with bounds of type T (Float64 by default). Apply it to a Vector{Interval{T}} of the matching type: a vector of a wider type would be narrowed to T on the way into the scratch buffers, and the error term is not written to cover that.

Accumulate the reflecting-boundary convolution of the float vector src into dest, which must have length k and is overwritten.

This is definition 8 of Galatolo-Monge-Nisoli read literally: extend src by zero outside 1:k, convolve on ℤ, then push the mass that fell outside back in with the mirror map π of reflect_outward_idx. Every step preserves mass, so the result is Markov whatever the parity of the window.

Reflecting the input instead and reading off 1:k is the adjoint of this, not the same operator; with the even window this kernel builds it is not even mass-preserving, losing O(1/k) of the mass on a density that is not symmetric.

Wrap src periodically into the workspace w, so that w[j] = src[j-l mod k].

struct UniformKernelUlam{BC,T,S} <: NoiseKernel

Ulam discretization of the uniform noise kernel with half-width l on a partition B::Ulam. The type parameter BC specifies the boundary condition:

• :periodic → periodic wrap-around
• :reflecting → mirror reflection (projection π)

T is the type of the bounds, Float64 by default. The scratch buffers carry it, so applying a kernel to a Vector{Interval{T′}} with T′ wider than T would narrow the midpoints on the way in, which the error term added afterwards does not account for; build the kernel at the type you mean to use.

S is the summation scheme of the window sums, :sliding by default:

• :sliding → one running sum, updated by adding the entering entry and subtracting the leaving one, with Kahan compensation. The error of an entry is that of a compensated sum of every term processed so far, so it is bounded in terms of the whole vector and not of the window.
• :block → the block scheme of van Herk (Pattern Recognit. Lett. 13 (1992) 517–521) and Gil and Werman (IEEE TPAMI 15 (1993) 504–507), written there for running maxima and valid for any associative operation. The extended vector is cut into blocks of length n = 2l+1, sums are accumulated inside each block from the right and from the left, and every window, which meets at most two blocks, is the sum of one right-accumulated and one left-accumulated partial sum. No entry is subtracted, and each window sum is a floating-point sum of its own n terms, so the error of an entry is at most γₙ times the kernel applied to |v| at that entry. The cost is three additions per entry, as for :sliding.

Fields

• B::Ulam : the Ulam partition
• l::Int : half-width of the averaging window
• scratch_ext::Vector{T} : workspace of length k+2l for building extended vector
• scratch_sum::Vector{T} : workspace of length k for storing window sums
• scratch_suf::Vector{T} : workspace of length k+2l for the right-accumulated block sums of the :block scheme (empty for :sliding)

UniformKernelUlam(::Val{BC}, [T,] B::Ulam, l::Int; summation = :sliding)

Constructor for a UniformKernelUlam{BC,T,S} kernel on the partition B with half-width l, bounds of type T (Float64 by default) and summation scheme S = summation, either :sliding or :block. Allocates the necessary scratch buffers of sizes k+2l and k (k = length(B)).

*(K::UniformKernelUlam, v::AbstractVector)

Non-mutating application of the kernel to vector v. This makes a copy of v internally and dispatches to the appropriate mul! method.

mul!(K::UniformKernelUlam{BC,T,:block}, v::Vector{Interval{T}})

In-place application of the kernel to a vector of intervals with the block scheme. With m and r the midpoints and radii of v, the result at entry i is centred at the computed (Km)_i with radius

γₙ (K|m|)_i + (K r)_i ,

the first term the rounding of the centre, the second the propagated radii (K is entrywise nonnegative). Both (K|m|)_i and (Kr)_i are sums of nonnegative terms, which the block scheme computes with relative error at most γₙ, so each is bounded above by its computed value divided by 1 - γₙ; the radius is formed with upward rounding from these. The radius is local to the window, unlike the uniform one of the :sliding scheme.

mul!(K::UniformKernelUlam{BC,T,:block}, v::Vector{T})

In-place application of the kernel with the block scheme. Every entry of the result is a window sum of n = 2l+1 terms, formed with n-1 additions and no subtraction, divided by n; by the bound for summation in any order [Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., (4.3)], together with the final division,

|fl((Kv)_i) - (Kv)_i| ≤ γₙ (K|v|)_i ,

so that the error is local and, in any monotone norm, at most γₙ ‖K‖ ‖v‖.

mul!(K::UniformKernelUlam, v::Vector{Interval})

In-place application of the uniform kernel to a vector of intervals. The operation is performed on midpoints with sliding window sums, and then a uniform interval error bound is added to account for radii.

mul!(K::UniformKernelUlam{BC,T}, v::Vector{T})

In-place application of the uniform kernel to a real vector v. Uses a preallocated scratch extension vector and sliding-window sum with Kahan summation for numerical stability.

UniformKernelUlamPeriodic(B::Ulam, l::Int)

Construct a periodic uniform Ulam kernel on the partition B with half-width l (window size = 2l+1).

This operator acts as a stochastic convolution with uniform noise, where indices outside [1,k] wrap around periodically.

UniformKernelUlamReflecting(B::Ulam, l::Int)

Construct a reflecting uniform Ulam kernel on the partition B with half-width l (window size = 2l+1).

This operator acts as a stochastic convolution with uniform noise, where indices outside [1,k] are mapped back into [1,k] by the reflecting projection π (period-2 mirror).

_block_window_sums!(sums, ext, suf, k, n)

Write into sums[i], for i = 1:k, the floating-point sum of ext[i:i+n-1], computed without subtraction. The extended vector ext (length k+n-1) is cut into blocks 1:n, n+1:2n, …; suf[t] is the sum from t to the end of its block, accumulated from the right. A window starting at a block start is that block, suf[i]; any other window starting in block b ends in block b+1, and is suf[i] plus the sum from the start of block b+1 to the end of the window, accumulated from the left in p. Each window sum therefore uses n-1 additions of its own n terms and nothing else.

dfly(::Type{TotalVariation}, ::Type{L1}, N::UniformKernelUlam)

Return the Doeblin–Fortet–Lasota–Yorke inequality coefficients for the operator acting from L¹ to bounded variation (Total Variation). For the uniform kernel on k bins with half-width l, the effective noise size is ξ = (2l+1)/k, and the inequality is bounded by (0, 1/(2ξ)).

get_idx(::Val{:periodic}, i, k)

Return the periodic index corresponding to i in 1:k.

get_idx(::Val{:reflecting}, i, k)

Return the reflecting index corresponding to i in 1:k, according to the projection π.

nonzero_per_row(M::UniformKernelUlam)

Return the number of nonzeros per row of the transition matrix associated with the discretized uniform noise kernel. This equals the window size 2l+1.

opnormbound(B::Ulam, ::Type{L1}, M::UniformKernelUlam)

Return an upper bound on the operator norm in L¹. For the uniform noise kernel, the L¹ operator norm is exactly 1.

opradius(::Type{L1}, M::UniformKernelUlam)

Return the "radius" term in the L¹ Lasota–Yorke inequality. For the uniform kernel this is zero, but can be extended for interval-aware variants.

reflect_outward_idx(i, k)

Reflecting index mapping (projection π): maps an index i on ℤ into 1:k by period-2k reflection symmetry.

wrap_idx(i, k)

Periodic index mapping: wraps i into the range 1:k.

Information about a single preimage of an interval I under a branch Tₖ.

_reflect(i, size)

Reflecting boundary condition index mapping (0-based internally). Maps index i (1-based) into 1:size using reflection at boundaries.

estimate_L1_NL_discretization_error(var_NLf, derivs, δ, varRho)

Estimate ||NL(1-π)NLf||₁ — the L1 norm of the second-order discretization error.

For each interval: min(δ*varRho/4*deriv[i], 1) * var_NLf[i], summed and scaled by δ/2.

estimate_Wnorm_of_discretization_error(prep, density_vec, B_density::Ulam)

Estimate ||(1-π)Lf||_W where π is the Ulam projection and f is the density represented by density_vec on the basis B_density.

For each interval I and each preimage, computes two alternative bounds and takes the minimum:

• Variation-based: (δ/4) * estimate_var
• Measure-based: meas_in_ab

Returns (est, meas_Lf, var_Lf):

• est: the W-norm estimate (scalar, rounded up)
• meas_Lf: vector of measure of Lf per estimation interval
• var_Lf: vector of variation of Lf per estimation interval

estimate_measure_after_noise(meas_Lf, K::UniformKernelUlam)

Estimate the measure of NLf per interval by averaging meas_Lf over the noise window.

Returns meas_NLf vector of length Kest.

estimate_variation_after_noise(meas_Lf, var_Lf, K::UniformKernelUlam)

Estimate Var(NLf) — the variation of Lf after applying the noise operator N.

For each interval i, takes the minimum of:

• Measure-based: boundary values divided by noise width
• Variation-based: averaged variation over noise window

Returns var_NLf vector of length Kest.

noise_error_aposteriori(prep1, prep2, meas_vec, numeric_error,
                        B_density::Ulam, K::UniformKernelUlam,
                        α, sumCi)

Compute the a posteriori L1 error bound for a noisy dynamical system using the actual computed density meas_vec.

Arguments

• prep1: output of prepare_Wnorm_estimate
• prep2: output of prepare_derivative_bounds
• meas_vec: the computed density vector (approximation of the invariant density)
• numeric_error: the numerical error from the eigenvector computation (e.g., residualboundnoise)
• B_density: the Ulam basis on which meas_vec lives
• K: the UniformKernelUlam noise kernel
• α: spectral contraction rate (should be < 1)
• sumCi: Σ Cᵢ from the norms of powers

Returns

(L1Error, apriori_error, Linf_estimate):

• L1Error: the a posteriori L1 error bound (should be tighter than generic)
• apriori_error: the a priori error bound for comparison
• Linf_estimate: per-interval L∞ error estimates (vector)

noise_error_apriori(K_size, noise_size_rel, α, sumCi)

Simple a priori error bound for the noisy system.

error ≤ (1/(1-α)) * (2*sumCi + 1) * (δ/2) * varRho

where δ = 1/K_size and varRho = 2/noise_size_rel.

prepare_Wnorm_estimate(D::PwMap, B_basis::Ulam, B_est::Ulam)

For every interval Iᵢ in the estimation partition B_est, compute preimage data under each branch of D.

Returns Vector{Vector{PreimageInfo}} of length length(B_est):

• prep[i] contains one PreimageInfo per branch of D that has a preimage overlapping with Iᵢ.

This is the Julia equivalent of the Python preparation_for_estimate1.

prepare_derivative_bounds(D::PwMap, B::Ulam)

For each interval in the partition of B, compute an upper bound for |T'|. This is used to estimate the W → W norm.

Returns a Vector{Float64} of length length(B).

total_variation(v::AbstractVector)

Compute a rigorous upper bound for the total variation of a step-function represented by the vector v:

Var(v) = Σᵢ |v[i] − v[i+1]|

Uses directed rounding (round-up) to keep the result rigorous.

plot_noisy_system(D::PwMap, K::UniformKernelUlam, w, L; n_samples=1000, n_points=500)

Plot a 3-panel visualization of a dynamical system with uniform noise. Requires using Plots to load the extension that implements this function.

See the PlotsExt extension for full documentation.