Beurling–Kato Polynomial Stability Regularizer

Integrate the regularizer into the recurrent transition or SSM state-update matrix, rather than the input or output projections. For an RNN with h_{t+1}=phi(T h_t+B x_t), apply it to T. For a residual or deep-equilibrium block, apply it to the linear transition reused during iteration. Begin with three fixed real probes: w_1(z)=1-z, w_2(z)=1+z, and w_3(z)=1-z^2. Estimate each ||w_j||_infty by evaluating the polynomial on a dense unit-circle grid z_m=exp(2 pi i m/M), using max_m |w_j(z_m)|. During each optimizer step, construct P_0=I and P_k=P_{k-1}T, then form q_j(T)=sum_k a_{jk}P_k. Estimate ||q_j(T)||_2 with one or two differentiable power iterations: initialize a random vector v, compute u=Mv/(||Mv||_2+epsilon), then v=M^T u/(||M^T u||_2+epsilon), and use ||Mv||_2 as the estimate. Add L_BK to the task loss with a small coefficient, using the hinge-square penalty in the displayed formula. Every 100 steps, separately measure ||T^n||_2 for n in {1,2,4,8,16}; these measurements may initially be detached from backpropagation. The polynomial bank, margin delta, horizon set, and bound B are empirical choices; the mathematical ingredients are the polynomial operator evaluation and norm-gap target. First test a one-layer tanh RNN and a two-layer linear SSM on sequential MNIST and the Adding Problem. Compare no regularizer, spectral-radius regularization, and this method at equal hidden size and FLOPs. Record validation accuracy, training loss, gradient norms, maximum sampled power norm, and exploding or vanishing-gradient frequency. A successful result is reduced transient amplification and fewer long-sequence failures at equal or better accuracy, with less than 10 percent training overhead.