△ Mechanism confirmed, baseline not beaten 2026

Padé-Hermite Neural ODE Integrator

Implementation & benchmark of arXiv:2607.24592 — Arbitrary-Order Padé-Closed Anchored Two-Derivative Time Discretizations: $s$ Active Stages, Order $2s$, and $L$-Stability

Usefulness7/10

Difficulty7/10

Novelty7/10

Dynamics Speedup Stability Accuracy Diffusion sampling World model Inference speedup Numerical analysis Dynamical systems Linear algebra

Source paper: Arbitrary-Order Padé-Closed Anchored Two-Derivative Time Discretizations: $s$ Active Stages, Order $2s$, and $L$-Stability arXiv:2607.24592 ⓘ · analyzed Aug 31, 2026

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

Idea description

Replace explicit Runge-Kutta integration in a neural ODE or probability-flow ODE sampler with the anchored two-derivative method. Each stage uses both the neural vector field and its total time derivative, while the coupled implicit solve is designed so the accepted map has an L-stable Padé stability function. The method should allow larger steps on stiff trajectories without amplifying fast decaying modes.

Formulas

$$s\ \text{active unknown stages}\quad\Longrightarrow\quad\text{global order }2s,\qquad\text{with }L\text{-stability for every }s.$$

$$U_i=u_n+h\sum_{j=0}^{s}A_{ij}f_j+h^{2}\sum_{j=0}^{s}\widehat A_{ij}g_j,\qquad g_j=D_tf_j=\partial_t f(t_j,U_j)+J_u f(t_j,U_j).$$

$$\det(I-zA-z^{2}\widehat A)=Q_s(z),\qquad R_s(z)=\frac{P_s(z)}{Q_s(z)},\qquad R_s(z)\to0\ \text{as }z\to-\infty.$$

$$u_{n+1}=u_n+h\sum_{j=0}^{s}b_jf_j+h^{2}\sum_{j=0}^{s}\widehat b_jg_j,\qquad z=h\lambda\ \text{for }u'=\lambda u.$$

Mathematical statement

The construction has s active unknown stages and global order 2s, summarized by s active unknown stages implying global order 2s with L-stability for every s. For u'=f(t,u), let U_0=u_n be the known anchor, let U_i for i=1,...,s be active stages, let t_i=t_n+c_i h, let f_i=f(t_i,U_i), and let g_i=D_t f(t_i,U_i)=partial_t f(t_i,U_i)+J_u f(t_i,U_i), where J_u is the Jacobian of f with respect to u. The stage equations are U_i=u_n+h sum_{j=0}^s A_{ij}f_j+h^2 sum_{j=0}^s hatA_{ij}g_j. The accepted update is u_{n+1}=u_n+h sum_{j=0}^s b_j f_j+h^2 sum_{j=0}^s hatb_j g_j. The first 2s Hermite moment conditions produce order 2s. The two remaining coefficients in each stage row are selected using the [s-1/s+1] Padé approximant P_s(z)/Q_s(z) to exp(z), giving det(I-zA-z^2 hatA)=Q_s(z) and R_s(z)=P_s(z)/Q_s(z), with z=h lambda for the scalar test equation u'=lambda u. L-stability means A-stability together with R_s(z) tending to zero as z tends to negative infinity. In a neural network, g_i can be computed by a JVP of f in the direction f_i, plus an explicit-time derivative when f depends directly on t.

Implementation notes

Integrate this at the ODE solver layer rather than inside the neural network parameters. Begin with s=2, nominally fourth order, and represent the active stages as a tensor U with shape batch by s by state dimension. Precompute the paper's node locations c_i and coefficient arrays A, hatA, b, and hatb; if the tables are unavailable, generate them offline by solving the Hermite moment equations together with the Padé closure equations and verify the resulting scalar stability function. For each step, set U_0=u_n and initialize active stages with an explicit Euler or extrapolated predictor. Iterate a damped fixed-point or Newton solve. At every iteration compute f_j=f(t_n+c_j h,U_j), then compute g_j using a forward-mode JVP of f with respect to U_j in tangent direction f_j; add an explicit partial-time derivative if needed. Build residuals r_i=U_i-u_n-h sum_j A_ij f_j-h^2 sum_j hatA_ij g_j and iterate until the normalized residual is below a fixed tolerance. Accept with the output formula and reduce h after failed nonlinear convergence. The paper supplies the order and Padé stability target; nonlinear tolerances and adaptive step control must be tuned empirically. First compare against Dormand-Prince and Heun on a two-dimensional stiff test field and a small CIFAR-10 probability-flow model. Measure accepted step count, network evaluations, JVP cost, trajectory error, and sample quality at equal compute. Success means larger stable steps and lower error or better FID at equal function/JVP evaluations.

Verification

Mechanism confirmed, baseline not beaten

Stage 1 · Toy mechanism gate: Passed ✓

Stage 2 · Mechanism transferred to benchmark: Confirmed ✓

Stage 2 · Practical benchmark result: No measurable effect

Methodology: 8 paired seeds, permutation test (p<0.05) against a learning-rate-tuned baseline. How verification works

Stage 1 — Mechanism check agent confidence 8/10

Built an s=2 Padé-Hermite anchored two-derivative integrator with coefficients generated from the paper's polynomial-basis construction, plus numerical verification and stiff toy comparisons. The mechanism manifested: moment residual was 1.1e-16, observed order was 3.94/3.96/3.78 versus predicted order 4, the left-half-plane grid maximum was |R|=1.0 versus A-stability, and the negative-axis decay slope was -1.96 versus the predicted -2. Node-independent amplification matched to about 1e-15, while nonlinear fixed-point convergence failed at larger stiff steps, so the stability claim is supported but practical nonlinear robustness is not.

Agent confidence: 8/10

Baseline: Explicit Euler stability boundary observed at z=-2: |1+z| was 0.9 at -1.9, 1.0 at -2.0, and 1.1 at -2.1. RK4 on the nonlinear lambda=80 test failed at h=0.05, had error 6.77e-8 at h=0.025, and error 2.24e-9 at h=0.0125.
Idea: s=2 Padé-Hermite: observed global orders 3.94, 3.96, 3.78 versus predicted 4; |R(-100)|=5.42e-4 and |R(-1000)| approximately 5.42e-6, giving observed log-log slope -1.96 versus predicted -2; amplification was identical across internal nodes [0.5, 0.25, 0.75] within roughly 1e-15. On the nonlinear stiff test, h=0.05 and 0.025 failed the simple damped fixed-point solve, while h=0.0125 converged with error 4.59e-10 but required 1194 total iterations.

Limitations: Only a small CPU NumPy toy experiment was tested; no neural ODE, JVP implementation, CIFAR-10 probability-flow model, GPU timing, adaptive step controller, Newton/Krylov solve, FLOP-normalized comparison, or sample-quality metric was evaluated. The A-stability check was a finite complex grid, and nonlinear failures reflect the basic fixed-point solver rather than the exact linear stability function.

How to run: python3 experiment.py

Stage 2 — Benchmark vs. tuned baseline (arXiv:2607.24592) Sep 2, 2026

Ahead of the baseline (−0.0%) but not statistically significant (p=0.7814, wins 4 / 8 paired seeds).

Benchmark: Dynamics forecasting
Model: rnn_small
Paired seeds: 8
Baseline mean: 0
Idea mean: 0
Effect (Δ): -0 (−0.0%; negative = idea better)
Wins: 4 / 8 paired seeds
p-value: 0.7814 (permutation test, 20 000 shuffles)
Smallest detectable effect: ±0.2%
Mechanism: Confirmed ✓
Practical effect: No effect
Baseline tuning: swept over 9 configs

Limitations:

Tested only the registered built-in dynamics pendulum forecast, not CIFAR probability-flow sampling or a genuine multi-stage implicit nonlinear solve. The benchmark GRU lacks forward-mode AD, so the JVP used reverse-mode differentiation. Global order, A/L-stability, wall-clock speed, function/JVP counts, and sample quality were not tested.

How to run: python3 bench_ph.py

Verdict computed by deterministic test code from paired-seed statistics — not by the language model.

Stage 2 — Benchmark Aug 31, 2026 latest Failed ✗

Agent confidence: 9/10

Implemented and ran the Padé-Hermite rollout intervention on the registered dynamics/rnn_small benchmark. The idea produced a tiny nonsignificant improvement: baseline MSE 5.681601408014103e-06 versus idea MSE 5.6799913750182895e-06, paired delta -1.6100329958135262e-09, permutation p=0.78145; therefore the benchmark verdict is no measurable effect.

Baseline: Best tuned baseline lr=0.01,h=0.025; full 8-seed mean MSE 5.681601408014103e-06.
Idea: Best idea setting selected from the same lr/h union; full 8-seed mean MSE 5.6799913750182895e-06; 4/8 paired wins; p=0.78145.

Benchmark result