Unit-circle-root Toeplitz mixer

Integrate this at the sequence-mixing operation, either as a standalone causal/noncausal FIR layer or as a replacement for a relative-position attention kernel. Let the hidden width be m and choose polynomial degree n-1. Parameterize n-1 unconstrained complex m x m matrices B_j, map each to a unitary U_j by a Cayley transform U_j=(I+iB_j)(I-iB_j)^{-1} when B_j is Hermitian, or by QR retraction after every optimizer step. Form Q(z)=A product_j(zI-U_j), expand it into coefficient matrices V_1,...,V_n, and compute S=sum_i V_i^*V_i+epsilon I. Normalize coefficients as V_i_tilde=V_i S^{-1/2}; optionally multiply by a scalar gain g clipped to a fixed range. For an input sequence x_1,...,x_L, compute the FIR output y_t=sum_{i=1}^n V_i_tilde x_{t-i} with zero, circular, or causal padding. For long sequences, evaluate the same operation with an FFT: H(omega_k)=sum_i V_i_tilde exp(i omega_k(n-i)) and y=IFFT(H(omega_k) FFT(x)). The paper-derived quantities are the polynomial, its unit-circle root construction, and the isometric coefficient normalization. Estimate no spectral quantities during training; verify the root constraint numerically by eigendecomposing each U_j and check coefficient conditioning through the largest singular value of H(omega_k). The first experiment should use a 2-layer Transformer or SSM on sequential MNIST, CIFAR-10 byte streams, or a 1B-token language-model slice, comparing an equal-parameter unconstrained FIR mixer and standard depthwise convolution. Measure training loss at equal FLOPs, gradient-norm spikes, validation accuracy/perplexity at lengths longer than training, and FFT versus direct runtime. A successful result is fewer instability spikes and better extrapolation to 2x context length without increasing parameter count; failure indicators are severe optimization slowdown from matrix products or worse loss at matched compute.