chc.toeplitz

Stability tier: evolving.

Structured-operator primitives for LTI dynamics: Toeplitz matvec, Levinson-Durbin, deconvolution.

A Toeplitz operator is a linear time-invariant system -- its action is convolution with the impulse response. These primitives back the structured route to the dynamic causal effect (chc.irf): the fast forward operator, the Yule-Walker AR fit (reflection coefficients and innovation power), deconvolution (recover the excitation from the response), and the Gohberg-Semencul fast inverse.

Two ways to apply T^{-1}: solve_toeplitz is a one-shot dense solve; Gohberg-Semencul (gohberg_semencul_generators + gohberg_semencul_apply) compresses the whole inverse into two generator vectors (the first and last columns of T^{-1}), then applies T^{-1} v in O(L log L) per vector by FFT. The win is amortised repeated application (deconvolve many signals through one operator) plus the compression itself -- real at any scale, decisive at the L ~ 1e3-1e4 long sequences the idea targets. NumPy float64 (a numerical routine, off jit paths).

toeplitz_matvec(first_col, first_row, x)

Toeplitz matrix-vector product T x in O(L log L) via circulant embedding + FFT.

T has first column first_col and first row first_row (shared corner entry). The 2L-circulant with first column [first_col, 0, reverse(first_row[1:])] contains T in its leading block, so one FFT convolution recovers T x. Differentiable (JAX).

circulant_symbol(kernel)

Eigenvalues of the circulant whose first column is kernel, as a half spectrum.

A circulant is diagonalised EXACTLY by the DFT: C = F* diag(lambda) F with lambda the DFT of its first column (validation/spectral_circulant.mac STEP 1, verified symbolically at n = 4 and n = 6). For a real kernel the spectrum is Hermitian, so the n//2 + 1 bins returned here carry all of it -- the missing bins are conjugates and have the same modulus.

circulant_matvec(kernel, x)

Circulant matvec C x in O(n log n), with C's first column given by kernel.

This is toeplitz_matvec without the embedding: a Toeplitz matrix needs a 2n-circulant to hold it, but a circulant is already one, so the n-point transform is exact and half the work. irfft also makes the output real by construction rather than by discarding an imaginary part that ought to be zero. Tested against the dense product and against the embedded route.

circulant_operator_norm(kernel)

Exact spectral norm ||C||_2 = max_k |lambda_k| -- attained, not bounded.

Circulants are normal (they commute with their adjoint, being polynomials in the cyclic shift), so their singular values are the moduli of their eigenvalues and the maximum is achieved on the corresponding Fourier mode. Contrast chc.residual.LipschitzResidual.lipschitz_constant, which rests on the Schur bound sigma_max(W) <= sqrt(||W||_1 ||W||_inf): an inequality with no witness. Proved in proofs/spectral_circulant.v (two_mode_bounded for <=, two_mode_norm_attained for the equality).

sample_autocorrelation(x, max_lag)

Biased sample autocorrelation r[0..max_lag] -- always PSD, so Levinson stays safe.

Divides every lag by n (not n - k): the biased estimator is guaranteed PSD, so its Toeplitz matrix is valid and Levinson-Durbin stays stable even from few samples, where the unbiased estimator can be non-PD with reflection coefficients outside (-1, 1). The robust distillation of Gohberg-Semencul covariance estimation (arXiv:2311.14995).

levinson_durbin(autocorrelation)

Solve the Yule-Walker system for an AR(p) model from its autocorrelation r[0..p].

Returns (ar_coeffs, reflection_coeffs, prediction_error): the AR coefficients a with x_t ~ sum_i a_i x_{t-i}, the reflection (PARCOR) coefficients, and the prediction-error variance -- the innovation power left after the AR structure is removed. O(p^2).

solve_toeplitz(first_col, first_row, rhs)

Deconvolution: solve T u = rhs -- recover an LTI system's excitation from its response.

One-shot dense solve. For many right-hand sides through one operator, use Gohberg-Semencul (compute generators once, then apply T^{-1} in O(L log L) per vector).

gohberg_semencul_generators(first_col, first_row)

The two Gohberg-Semencul generators: the first and last columns of T^{-1}.

Solve T x = e_1 and T y = e_n; these vectors encode the whole inverse (displacement rank 2). Pass them to gohberg_semencul_apply. Requires x[0] != 0 (generic for a nonsingular Toeplitz). Setup is a dense solve; the point is that application is then fast.

gohberg_semencul_covariance(snapshots, order, size=None)

Few-sample Toeplitz covariance estimate via AR / Gohberg-Semencul (arXiv:2311.14995).

Averages the biased autocovariance over snapshots (rows), fits an AR(order) by Levinson (the paper's closed-form projected-least-squares estimator a = R_w^{-1} r), and extends it to the maximum-entropy Toeplitz covariance. The estimate is positive definite and full rank even when the sample covariance is singular (N < size). Its inverse is the GS precision with generator alpha = (1/sigma^2) [1, -a] (the AR whitening filter, Eq 40), recoverable via gohberg_semencul_generators. Distils the paper's PLS estimator.

It is a regularised estimator: low variance (a win at small N) bought with model bias. Use it for few samples of a stationary, roughly low-order-AR process; on a non-Toeplitz / non-stationary covariance it plateaus at a bias floor and the raw sample covariance wins once N is large. Empirically ~30x lower NMSE than the SCM at N << size on matched AR data, fading to a tie as N grows.

gohberg_semencul_apply(x, y, v)

Apply T^{-1} v in O(L log L) from the generators x, y (Gohberg-Semencul).

T^{-1} = (1/x_0) [ L(x) U(J y) - L(S y) U(S J x) ] with L/U lower/upper-triangular Toeplitz, J reverse and S down-shift; each factor is applied by FFT. Amortises over many right-hand sides that share the operator.