Padé-contracting residual dynamics

Integrate this at the residual or SSM transition, before the pointwise MLP or readout. Start with a channel dimension d small enough that a dense solve is feasible. Store trainable matrices A and B, construct L=A-A.T-B.T@B, and use a positive learned step h=softplus(rho), optionally clipped to a chosen maximum. For each batch vector x, compute rhs=(I+h*L/3)@x. Form D=I-2*h*L/3+(h*h/6)*(L@L), then solve D@y=rhs with a batched LU solve. The core pseudocode is: L=A-A.T-B.T@B; h=softplus(rho); rhs=x+(h/3)*L@x; D=I-(2*h/3)*L+(h*h/6)*(L@L); y=solve(D,rhs); y=y+mlp(y). For the first test, omit the extra MLP residual so the effect of the rational transition is isolated. The mathematical part is the Padé rational map and seminegative parameterization. Estimate the actual discrete contraction empirically with max_batch ||y||/||x||, and estimate the energy seminorm using 2||B@x||^2; do not assume floating-point solves preserve the exact identity. First test on CIFAR-10 with a 4-layer channel-mixing residual MLP and on a synthetic long-horizon system-identification task with a non-normal ground-truth L. Compare explicit Euler, unconstrained residual, skew-symmetric Cayley, and the Padé block at matched matrix multiplications. Success means fewer exploding or NaN runs, stable training at larger h, lower long-horizon rollout error, or faster loss descent at equal depth. A follow-up can replace dense solves with iterative GMRES or low-rank B if the dense inverse dominates FLOPs.