Non-normal covariance equalizer

Integrate the method into the hidden-state update of a small tanh RNN, linear RNN, or diagonal-plus-low-rank SSM. The exact integration point is the recurrent state covariance, computed from the batch of hidden states after the recurrent update and before the output projection. Maintain an exponential moving average $\widehat{Q}_t=\alpha\widehat{Q}_{t-1}+(1-\alpha)C_t$, where $C_t=(1/B)\sum_{b=1}^{B}(h_{b,t}-\bar h_t)(h_{b,t}-\bar h_t)^\top+\delta I$; $B$ is batch size, $h_{b,t}$ is the hidden state, $\bar h_t$ is the batch mean, and $\delta$ prevents singularity. Every $K$ optimizer steps, symmetrize $\widehat Q$, compute its eigenvalues with torch.linalg.eigvalsh, normalize them into $p_i$, and add $\beta[\log r_{target}-H(\widehat Q)]_+^2$ to the task loss. Detach the EMA update if needed, but retain gradients through the current covariance estimate for the regularizer. The Lyapunov series is the mathematical interpretation; initially estimate $Q$ from activations rather than explicitly solving the equation. For large hidden dimensions, use randomized Lanczos estimates of the leading eigenvalues and stochastic trace estimates. Estimate $\rho(J)$ with 5--10 power iterations and optionally add the spectral-radius barrier only when it activates. First test on permuted-MNIST and the adding problem with a 128-unit tanh RNN. Compare plain training, spectral-radius regularization, and covariance regularization at matched FLOPs. Record accuracy, hidden-state effective rank, top-eigenvalue share, gradient norms, and failure rate from exploding or vanishing gradients. Success means reduced variance concentration and fewer unstable runs at equal accuracy, or better long-sequence accuracy at equal parameter count. Sweep $\beta$, $\alpha\in\{0.95,0.99\}$, and target effective-rank fractions of 0.25, 0.5, and 0.75 of the hidden dimension.