Numerical-Range Stabilization for Nonnormal State Dynamics

Integrate the regularizer into the transition operator of a linear state-space model or Elman-style RNN, specifically the matrix A in h_{t+1}=A h_t+B x_t. Choose center c and radius r; for contractive dynamics use c=0 and r<1, while for oscillatory systems choose c near the desired eigenvalue cluster. At each training step compute: loss = task_loss; nr = 0; for each theta in linspace(0, 2*pi, m, endpoint=False), form H = 0.5*(exp(-1j*theta)*(A-c*I) + exp(1j*theta)*(A.conj().T-conj(c)*I)); compute h = largest_eigenvalue(H); accumulate nr += relu(h-r)**2; then optimize loss + lambda_nr*nr/m. For hidden dimensions up to roughly 256, use a dense Hermitian eigensolver. For larger dimensions, estimate the largest eigenvalue with 3-10 Lanczos or power iterations, warm-starting each angular direction from the previous step. The exact mathematical quantity is the support function h_A(theta); the finite angular grid approximates the supremum over all directions, so evaluate a denser grid after training to detect missed violations. The first cheap experiment should use 64- or 128-dimensional linear SSMs and tanh RNNs on sequential or permuted MNIST, comparing this method against eigenvalue clipping and spectral normalization at equal parameter count and optimizer budget. Measure training-loss descent, gradient-norm spikes, activation growth, rollout norms ||A^k|| for k up to 1000, and validation accuracy. Use deliberately nonnormal upper-triangular initializations to expose the difference from spectral-radius control. Success means fewer transient explosions, lower worst-case rollout amplification, and improved long-horizon accuracy without reducing short-horizon performance. Ablate lambda_nr, angular-grid size m, center c, and radius r.