# Padé-contracting residual dynamics

- ID: 2038
- Canonical URL: https://synthcore.org/idea/2038/pad-contracting-residual-dynamics
- API JSON: https://synthcore.org/api/idea/2038.json
- API Markdown: https://synthcore.org/api/idea/2038.md
- Verification status: unverified
- Source: [arXiv:2608.09239](https://arxiv.org/abs/2608.09239)
- Category: dynamics
- Solves: stability, accuracy
- ML areas: training-dynamics, ssm
- Math tags: dynamical-systems, linear-algebra, numerical-analysis
- Ratings: usefulness 6/10; difficulty 6/10; novelty 6/10

## Idea description

Replace an explicit residual layer x_{k+1}=x_k+hLx_k with a first-subdiagonal Padé rational layer. For the lowest nontrivial approximant, use R_{1,2}(z)=(1+z/3)/(1-2z/3+z^2/6), so x_{k+1}=R_{1,2}(hL)x_k; parameterize L to have a negative-semidefinite symmetric part, preventing exploding activations even for large learned step sizes.

## Mathematical statement

The paper studies the linear system u'(t)=Lu(t) on a real Hilbert space, where L is seminegative: <Lv,v> <= 0 for every vector v. It defines the symmetric positive-semidefinite energy form [w,v]_L=-<Lw,v>-<w,Lv> and seminorm |v|_L^2=[v,v]_L. The continuous energy decreases because the symmetric part of L is nonpositive. The paper's new ingredient is an explicit Cholesky-type factorization of the coefficient matrix in the discrete energy identity for first-subdiagonal Padé approximants, yielding unconditional contractivity of the rational time step. The implementable lowest-order first-subdiagonal exponential Padé map is R_{1,2}(z)=(1+z/3)/(1-2z/3+z^2/6), with z=hL. Here h>0 is the layer step, L is the learned generator, and the inverse denominator is applied to the numerator output. Enforce seminegativity with L=A-A^T-B^T B, since <Lv,v>=-||Bv||^2 <= 0.

## Key formulas

- $$\frac{\mathrm{d}}{\mathrm{d}t}u(t)=Lu(t),\qquad \left\langle Lv,v\right\rangle\leq 0\quad\forall v.$$
- $$[w,v]_L:=-\left\langle Lw,v\right\rangle-\left\langle w,Lv\right\rangle,\qquad |v|_L^2=[v,v]_L\geq 0.$$
- $$R_{1,2}(z)=\frac{1+z/3}{1-2z/3+z^2/6},\qquad x_{k+1}=\left(I-\frac{2h}{3}L+\frac{h^2}{6}L^2\right)^{-1}\left(I+\frac{h}{3}L\right)x_k.$$
- $$L=A-A^\top-B^\top B\quad\Longrightarrow\quad \langle Lv,v\rangle=-\|Bv\|^2\leq0,\qquad [v,v]_L=2\|Bv\|^2.$$

## Implementation notes

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.

## Disclaimer

AI-generated research hypothesis, automatically tested. Not peer-reviewed.