Inertia-Constrained Kreiss RNN

(1) Integration point: use this matrix as the transition operator of a vanilla RNN or linear state-space layer, h_(t+1)=A h_t+B u_t. Choose hidden width n and a small inertia index q, for example q=4, and fix J=diag(+1,...,+1,-1,...,-1). Store an unconstrained real matrix H0, symmetrize it as H=(H0+H0^T)/2, and set A=JH. For complex-valued models use conjugate transpose instead. (2) Pseudocode: construct A=J(H0+H0^T)/2; estimate the spectral radius rho_hat(A) with several Arnoldi or power iterations; set A_s=A*min(1,(1-epsilon)/max(rho_hat(A),epsilon)); sample z=rho_i exp(i theta_j), with rho_i in {1.01,1.05,1.2,1.5,2.0}; estimate sigma_max((zI-A_s)^(-1)) using a few power iterations, each requiring linear solves with zI-A_s; compute L_K as the mean squared value of (|z|-1) times this singular-value estimate; optimize task_loss+lambda_K L_K. Use A_s in the forward pass. (3) The J-selfadjoint parameterization and inertia factor come from the paper. The spectral radius, sampled resolvent maximum, and finite-horizon power growth max over k less than or equal to T of ||A_s^k||_2 are empirical estimates. For a generic dense H, use r at most 2q+1 rather than estimating d(A). (4) First experiment: compare a width-128 tanh RNN, a spectrally rescaled unconstrained RNN, and this q=4 model on sequential MNIST and the adding or copying task. Match initialization, optimizer, and parameter count as closely as possible. Sweep lambda_K and measure validation accuracy, gradient norm variance, loss after long sequences, sampled K(A), and max over k less than or equal to 500 of ||A_s^k||_2. A successful result is lower transient amplification and fewer exploding gradients with equal or better long-sequence accuracy at comparable training loss and FLOPs. Cache the resolvent penalty every 10 to 50 optimizer steps to reduce overhead.