Lyapunov Canonical-Angle Regularizer

Integrate this into the recurrent or state-space block, rather than into the output head. Let the hidden dimension be n, with transition A, input projection B, and readout C. First stabilize A, for example with A=rho A_raw/(max(1,||A_raw||_2)+1e-4), where rho is a learned or fixed scalar constrained below one. Every K optimizer steps, construct finite-horizon Gramians using T unrolled powers: initialize P=0, Q=0, Xc=B, Xo=C^T; for k=0,...,T-1, update P=P+Xc Xc^T, Q=Q+Xo Xo^T, then Xc=A Xc and Xo=A^T Xo. Replace P and Q by P+epsilon I and Q+epsilon I, with epsilon=1e-4 trace(P)/n. Compute Q^{1/2} by eigendecomposition, form S=Q^{1/2} P Q^{1/2}, and compute its eigenvalues s_i. Set lambda_i=1/(1+s_i), then add alpha times L_angle to the task loss. Use lambda_target=0.5 initially, or instead penalize spectral variance. Detach the eigensystem every K steps if memory is limiting. The theorem supplies the canonical-angle interpretation and the PQ to (I+PQ)^{-1} conversion; finite horizon, ridge, and target spectrum are empirical choices. First test a 256-unit Elman RNN and a small diagonal-plus-low-rank SSM on sequential MNIST, comparing against weight decay alone. Track validation accuracy, gradient norms over sequence time, Gramian condition numbers, and failures from exploding gradients. A positive result is equal or higher accuracy with fewer unstable runs and improved extrapolation from sequence length 100 to 500.