Pith. sign in

REVIEW 4 major objections 4 minor 15 references

Training Stiff Neural Ordinary Differential Equations with Explicit Exponential Integration Methods

T0 review · 4 major / 4 minor · reviewed 2026-08-12 · deepseek-v4-flash

Pith's one-line read This paper establishes that the explicit integrating factor Euler (IF Euler) method can train stiff neural ODEs stably and cost-effectively, succeeding on the stiff Van der Pol oscillator where all tested implicit single-step methods fail.

desk verdict Classical IF Euler applied to stiff neural ODEs: the linear results are clean and exact, the Van der Pol claim is interesting but under-documented, and the higher-order method discussion is outsourced to a companion paper. read the letter →

arxiv 2412.01181 v1 pith:BIFYLZIP submitted 2024-12-02 math.NA cs.AIcs.LGcs.NAcs.SC

classification math.NAcs.AIcs.LGcs.NAcs.SC MSC 65L0465L2065F6068T07
keywords neuralordinarydifferentialequationsstiffODEsexponentialintegrationintegratingfactorEulerVanderPoloscillatorpolynomialnetworksdiscretize-optimizematrix
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper argues that an explicit exponential integration method, the integrating factor Euler (IF Euler) scheme, can train stiff neural ODEs stably and cheaply, a task where standard implicit single-step solvers fail on the stiff Van der Pol oscillator. Because IF Euler avoids solving nonlinear systems at each step, it is much cheaper per step than backward Euler, trapezoid, or Radau methods. The method recovers stiff linear systems almost exactly even from very few training points, and it scales to a 10-dimensional stiff linear model. Its first-order accuracy, however, means that for nonlinear stiff problems it trades accuracy for stability, which the paper presents as an open limitation.

What carries the argument

The engine of the method is the integrating factor transformation $w(t) = e^{-Lt}y(t)$, which converts the ODE into $\dot{w} = e^{-Lt}N(t, e^{Lt}w)$ and lets forward Euler handle only the slow nonlinear remainder while the stiff linear part is solved exactly by the matrix exponential. The Jacobian $L$ is evaluated once per integration interval at the initial condition, so the matrix exponential $e^{Lh}$ acts as a fixed stabilizing preconditioner for that interval. This is combined with a $\pi$-net polynomial neural network for the right-hand side, whose polynomial outputs allow exact equation recovery, and the discretize-optimize training paradigm, which makes gradients independent of the forward solver's accuracy.

What would settle it

Take a stiff nonlinear system whose Jacobian's dominant eigenvalue changes by orders of magnitude inside a single integration interval (for example, a Van der Pol oscillator with a time-dependent $\mu(t)$ varying between 1 and 1000), and train it with IF Euler using the paper's large-step protocol; if training diverges or the recovered system is not stiff, the frozen-Jacobian assumption is the limiting factor and the method's robustness does not extend beyond systems with locally stable fast manifolds.

Watch

Extended reading notes

Core claim

The central claim is that the explicit integrating factor Euler method is the most robust exponential integration method for training stiff neural ODEs. The scheme decomposes the right-hand side into a linear part frozen at the interval's initial Jacobian $L = df/dy(y_0)$ and a nonlinear remainder $N(t,y) = f(t,y) - Ly$, then exactly integrates the linear part with the matrix exponential while applying forward Euler to the transformed remainder: $y_{n+1} = e^{Lh}(y_n + hN_n)$. On the stiff Van der Pol oscillator with $\mu = 1000$, every implicit single-step method the authors tried diverged during training, whereas IF Euler remained stable even with large step sizes and recovered equations whose dominant terms approximate the true system as training data increase. The method also reproduces stiff linear systems to near machine precision, matching or beating Radau5 with scarce data.

Load-bearing premise

The method's stability rests on assuming that the Jacobian evaluated once at the start of each integration interval continues to represent the stiff dynamics for the entire step; if the stiffness changes direction or magnitude within a step, the matrix exponential may stabilize the wrong linearization and the training can fail or learn incorrect dynamics.

Editorial extensions

If this is right

  • Stiff neural ODE training can be done with explicit steps, removing the per-step nonlinear solves that make implicit methods costly and hard to differentiate.
  • Linear or nearly linear stiff dynamics can be identified from very few data points, as the 10D linear model is recovered almost exactly from 17 points, which matters for scientific system identification.
  • The stiff Van der Pol oscillator becomes trainable with large step sizes, a benchmark that previously defeated implicit single-step training in the authors' tests.
  • The first-order accuracy ceiling means that the practical benefit for nonlinear stiff problems is stability at coarse resolution, not high accuracy; a stable higher-order exponential method would be the next step.
  • The same explicit-exponential strategy could be applied to other differentiable simulations, such as neural PDEs, where repeated low-cost integration is required.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • If the frozen-Jacobian assumption is the real source of stability, IF Euler should extend most naturally to systems whose stiff subspace is roughly constant along the trajectory; systems with rapidly rotating or sign-changing Jacobians would likely need a per-step Jacobian update.
  • The Van der Pol success may owe to its slow-fast structure, where a locally fixed linearization captures the fast stable manifold; a stiff system without such a clean timescale separation should be a harder test.
  • One could test this by running IF Euler on a stiff system with a time-dependent stiffness parameter (e.g., $\mu(t)$) and comparing against a variant that recomputes $L$ at each step.
  • A higher-order extension might come from Rosenbrock-W methods, which share the same one-Jacobian-per-step structure but add order conditions without reintroducing nonlinear solves; the paper's negative result on other exponential methods suggests this route is not obvious.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 4 minor

Summary. The paper proposes applying the explicit exponential integrating factor Euler (IF Euler) method to train neural ODEs on stiff systems. The authors derive the scheme yn+1 = e^{Lh}(yn + hNn) with L frozen at the initial Jacobian of each integration interval, test it on a one-dimensional linear stiff equation, a ten-dimensional linear stiff system, a three-dimensional nonlinear stiff system, and the stiff Van der Pol oscillator (μ=1000), and compare against backward Euler, trapezoid, Radau3, and Radau5. The main claims are that IF Euler recovers stiff linear systems to numerical precision, trains the Van der Pol model where the implicit methods are asserted to fail, and is more efficient than implicit methods because it avoids nonlinear solves. The paper also states that higher-order explicit exponential methods were not found to be viable, referring to a companion paper.

Significance. If the nonlinear results were fully supported, the paper would provide a simple, low-cost alternative to implicit solvers for stiff neural ODEs, potentially widening the use of neural ODEs in scientific applications. The derivation of IF Euler is clear and the linear-system experiments are detailed and internally consistent; however, the exact recovery of linear systems is a by-construction property of integrating the linear part exactly, so it does not by itself validate the method for general stiff problems. The paper's broader significance hinges on Examples 3 and 4 and on the comparative claims against implicit and higher-order exponential methods, all of which are currently under-documented. The manuscript is clearly written but does not yet provide enough evidence to establish the central claims.

major comments (4)
  1. [Section III D, Example 4] The statement that 'every implicit integration technique we tried failed to properly train the stiff neural ODE' is made without any supporting data. Table 11 and Figure 11 report results only for IF Euler; no loss curves, parameter errors, or descriptions of the convergence failures are provided for backward Euler, trapezoid, Radau3, or Radau5. Since the paper's central contrast between IF Euler and implicit methods rests on this claim, the authors should supply quantitative evidence of the implicit-method failures, such as best-attained parameter errors and a statement of the failure mode (e.g., NaN loss, non-convergence of Newton iterations).
  2. [Section II E, Eqs. (12) and (17)] The method freezes the Jacobian L = df/dy(y0) at the initial condition of each integration interval. For nonlinear systems (Examples 3 and 4), no error bound or quantitative criterion is given for when this fixed linearization remains representative over the integration step. The matrix exponential may stabilize a linearization that is far from the local dynamics, and the explicit treatment of N(t,y) in Eq. (17) is not analyzed. The favorable Van der Pol result could be specific to that problem. The authors should provide a local error estimate in terms of the Jacobian variation over the interval, or at least report the Jacobian variation along the training trajectories and test the sensitivity to using the midpoint or endpoint linearization instead of the initial-condition linearization.
  3. [Sections I, II E, and Conclusion, Ref. 75] The claim that no higher-order explicit exponential integration method is viable for stiff neural ODEs is deferred to a companion paper (Ref. 75) by the same authors, which is not included. This assertion is load-bearing for the paper's framing as leaving higher-order methods as an open problem and for the exclusive focus on IF Euler. Without presenting at least a summary of the accuracy and stability results for representative higher-order exponential methods (e.g., ETD2RK or a higher-order Lawson method) on a stiff benchmark, the claim should be weakened to a statement about the methods the authors tried, and the universality removed.
  4. [Section III C and D, Examples 3 and 4] The experiments report a single trial per configuration with no seeds, error bars, or key hyperparameters (optimizer, learning rate, number of epochs, network width, initialization). The quantitative interpretation is further obscured because 'success' is not defined: in Example 4 the recovered equations deviate substantially from the true system (Table 11, n=100: the y-coefficient is 604.86 instead of 1000 and the x^2 y coefficient is -573.02 instead of -1000), and the error decays slowly with n. The authors should state a success criterion (e.g., relative parameter error below a threshold or trajectory error on a hold-out set) and report mean and variance over multiple random seeds.
minor comments (4)
  1. [Section II E, Eq. (17)] The text calls IF Euler an 'A-stable first-order' scheme, but the stability statement applies only to the exact integration of the linear part; the explicit Euler treatment of the nonlinear remainder can still impose step-size restrictions. This should be clarified to avoid overstating the stability of the full scheme.
  2. [Figure 5 caption] The caption contains a duplicated word: 'explicit exponential integrating factor Euler method method' should read 'explicit exponential integrating factor Euler method'.
  3. [Tables 7-10] The IF Euler rows for Example 3 contain many spurious higher-degree terms (e.g., y1^2 and y1 y2 terms with large coefficients) that are not present in the true model; this is a relevant observation that could be reported in the main text as a consequence of the first-order error and the frozen Jacobian, rather than appearing only in the tables.
  4. [Section III D, Eq. (21)] The text says the Van der Pol stiffness is parameterized by mu and then sets mu=1000, but Eq. (21) is written in the standard two-dimensional form with the coefficient 1000 appearing explicitly; stating the relation between Eq. (21) and the usual x'' + mu(x^2-1)x' + x = 0 form would help readers.

Circularity Check

1 steps flagged · score 4.0 of 10

Self-contained training experiments; the only circular element is the claim that no higher-order explicit exponential method is viable, which is deferred entirely to a companion paper by the same authors (Ref. 75).

  1. self citation load bearing [Section II.E, paragraph after Eq. (17); see also Section I and the Conclusion]
    "Despite an extensive search for a higher-order exponential integration method, our recent analysis (see Ref. 75) indicates that exponential integration methods fail to improve upon the first-order accuracy of IF Euler while remaining stable, revealing the IF Euler method as the only reliable choice for repeated, inexpensive integration in applications such as neural ODEs."

    The paper's selection of IF Euler and its negative claim about all higher-order explicit exponential methods rest entirely on Ref. 75, which is authored by the same two authors and is not included or reproduced. This negative result is load-bearing for the paper's framing ('the only reliable choice') and for the Conclusion's statement that a higher-order method could not be identified; the present paper provides no experiment, table, or theorem that excludes higher-order exponential methods. The justification therefore reduces to a self-citation chain rather than to evidence in the manuscript. The core training demonstrations (Examples 1-4) are independently executed in this paper, so the circularity is partial rather than total.

full rationale

The derivation of the IF Euler update (Eq. 17) from the integrating factor transformation (Eqs. 14-16) is straightforward and not circular: applying forward Euler to w' = exp(-Lt)N(t, exp(Lt)y) yields y_{n+1} = exp(Lh)(y_n + h N_n) by definition. The exactness on the stiff linear systems (Examples 1 and 2) is a known mathematical property of exponential integrators when N = 0 and L is the true Jacobian; presenting it as a 'notable finding' is an overstatement but not a circular reduction, since the network still has to fit the parameters from data. Examples 3 and 4 are empirical benchmarks with the method's own recovered equations reported; the frozen-Jacobian assumption (Eq. 12) is an unvalidated modeling choice and a correctness risk, but it is not circular because no claim is derived from it. The one genuine circular element is the unsupported, load-bearing self-citation to Ref. 75 for the claim that IF Euler is the only viable exponential integration method; because the central stiff-training results are self-contained, the appropriate score is 4, not higher.

Assumptions & free parameters 0 free parameters · 4 assumptions · 0 invented entities

The paper adds no fitted constants; its central method rests on standard exponential-integrator theory plus two empirical premises that are asserted rather than demonstrated: that implicit methods fail on the Van der Pol example, and that no higher-order explicit exponential method is viable. The latter is supported only by a self-cited companion paper.

assumptions (4)
  • standard math The matrix exponential exactly solves the linear part dy/dt = Ly.
    Used in Section II.E to derive IF Euler and to explain exactness on linear systems.
  • domain assumption The Jacobian L evaluated at the initial condition of each subinterval adequately represents the stiff linear dynamics over that interval.
    Section II.E Eq. 12; the method freezes L = df/dy(y0) over the step; no error analysis for strongly varying Jacobians is provided.
  • ad hoc to paper Higher-order explicit exponential integrators fail to outperform IF Euler for stiff neural ODEs.
    Asserted in Section II.E citing only the authors' companion paper Ref 75; no independent verification or derivation in this paper.
  • ad hoc to paper Implicit single-step methods fail to train the stiff Van der Pol neural ODE.
    Asserted in Section III.D without presenting failure data; used as the key comparison showing IF Euler's advantage.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Training Stiff Neural Ordinary Differential Equations with Explicit Exponential Integration Methods." pith.science (2026). https://pith.science/paper/BIFYLZIP

@misc{pith2026241201181,
  author       = {Pith},
  title        = {Pith review of: Training Stiff Neural Ordinary Differential Equations with Explicit Exponential Integration Methods},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/BIFYLZIP}},
  note         = {Machine review of arXiv:2412.01181}
}
read the original abstract

Stiff ordinary differential equations (ODEs) are common in many science and engineering fields, but standard neural ODE approaches struggle to accurately learn these stiff systems, posing a significant barrier to widespread adoption of neural ODEs. In our earlier work, we addressed this challenge by utilizing single-step implicit methods for solving stiff neural ODEs. While effective, these implicit methods are computationally costly and can be complex to implement. This paper expands on our earlier work by exploring explicit exponential integration methods as a more efficient alternative. We evaluate the potential of these explicit methods to handle stiff dynamics in neural ODEs, aiming to enhance their applicability to a broader range of scientific and engineering problems. We found the integrating factor Euler (IF Euler) method to excel in stability and efficiency. While implicit schemes failed to train the stiff Van der Pol oscillator, the IF Euler method succeeded, even with large step sizes. However, IF Euler's first-order accuracy limits its use, leaving the development of higher-order methods for stiff neural ODEs an open research problem.

Figures

Figures reproduced from arXiv: 2412.01181 by the authors.

Figure 1
Figure 1. FIG. 1. The neural network architecture of [PITH_FULL_IMAGE:figures/full_fig_p003_1.png] view at source ↗
Figure 2
Figure 2. FIG. 2. Comparison of the integration of the deterministic stiff van der Pol oscillator with [PITH_FULL_IMAGE:figures/full_fig_p004_2.png] view at source ↗
Figure 3
Figure 3. FIG. 3. Illustration of (a) Discretize-Optimize and (b) Optimize-Discretize methods. For Discretize-Optimize, black and red lines denote the [PITH_FULL_IMAGE:figures/full_fig_p005_3.png] view at source ↗
Figures from the paper (8 more)
Figure 4
Figure 4. Figure 4: FIG. 4. For the equation [PITH_FULL_IMAGE:figures/full_fig_p009_4.png]
Figure 5
Figure 5. Figure 5: FIG. 5. The fractional parameter relative error (non-percentage) is plotted against the number of training data points for the equation [PITH_FULL_IMAGE:figures/full_fig_p009_5.png]
Figure 6
Figure 6. Figure 6: FIG. 6. The training data for the 10-dimensional linear stiff system of ODEs (Example 2). The stiffness can be seen by the fact that each of [PITH_FULL_IMAGE:figures/full_fig_p012_6.png]
Figure 7
Figure 7. Figure 7: FIG. 7. For Example 2, we plot the fractional parameter relative error (not a percentage) against the number of training data points for [PITH_FULL_IMAGE:figures/full_fig_p018_7.png]
Figure 8
Figure 8. Figure 8: FIG. 8. Training data for Example 3, a 3D nonlinear stiff ODE system. This is the training data for the model corresponding to [PITH_FULL_IMAGE:figures/full_fig_p020_8.png]
Figure 9
Figure 9. Figure 9: FIG. 9. For Example 3, we plot the fractional parameter relative error (not a percentage) against the number of training data points for the [PITH_FULL_IMAGE:figures/full_fig_p021_9.png]
Figure 10
Figure 10. Figure 10: FIG. 10. The training region for the stiff Van der Pol model (see Eqn.21) [PITH_FULL_IMAGE:figures/full_fig_p027_10.png]
Figure 11
Figure 11. Figure 11: FIG. 11. For the stiff Van der Pol model, we plot the fractional parameter relative error (not a percentage) against the number of training data [PITH_FULL_IMAGE:figures/full_fig_p027_11.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

15 extracted references · 11 canonical work pages

  1. [4]

    Rows correspond to equations for y′ 0 to y′ 9, and columns represent coefficients for y0 to y9 and b. Neural ODEs for Stiff Systems: Explicit Exponential Integration Methods 15 y0 y1 y2 y3 y4 y5 y6 y7 y8 y9 b y′ 0 -12.2841 7.5659 -1.3450 0.8123 -0.7411 0.8725 -1.3970 1.6440 -0.0788 -0.6361 8.9738 y′ 1 10.0457 -28.0456 10.3687 -2.8062 1.6789 -1.8580 2.4636...

  2. [5]

    For comparison, we show the explicit exponential integrating factor Euler method method alongside a few implicit schemes

    The fractional parameter relative error (non-percentage) is plotted against the number of training data points for the equation y′ = −10000y. For comparison, we show the explicit exponential integrating factor Euler method method alongside a few implicit schemes. Neural ODEs for Stiff Systems: Explicit Exponential Integration Methods 10 y′ = −10000y n Equ...

  3. [6]

    The stiffness can be seen by the fact that each of the variables approach steady-state at a different timescale

    The training data for the 10-dimensional linear stiff system of ODEs (Example 2). The stiffness can be seen by the fact that each of the variables approach steady-state at a different timescale. Neural ODEs for Stiff Systems: Explicit Exponential Integration Methods 13 y0 y1 y2 y3 y4 y5 y6 y7 y8 y9 b y′ 0 -14.7836 14.5781 -8.4924 6.1276 -7.1492 9.2093 -22...

  4. [7]

    Rows correspond to equations for y′ 0 to y′ 9, and columns represent coefficients for y0 to y9 and b. Neural ODEs for Stiff Systems: Explicit Exponential Integration Methods 16 y0 y1 y2 y3 y4 y5 y6 y7 y8 y9 b y′ 0 -10.6184 5.8565 -0.6182 0.4069 -0.3637 0.4038 -0.3397 0.0666 0.1203 0.1399 1.4865 y′ 1 5.6450 -21.2199 6.0343 -0.6120 0.3982 -0.4303 0.5109 -0....

  5. [8]

    Rows correspond to equations for y′ 0 to y′ 9, and columns represent coefficients for y0 to y9 and b. Neural ODEs for Stiff Systems: Explicit Exponential Integration Methods 17 y0 y1 y2 y3 y4 y5 y6 y7 y8 y9 b y′ 0 -9.9886 4.9797 1.9113e-02 -1.3712e-02 1.4039e-02 -1.7570e-02 3.9613e-02 -7.1463e-02 5.3269e-02 -1.4328e-02 -1.1274e-03 y′ 1 5.0368 -20.0400 5.0...

  6. [12]

    Locally Regularized Neural Differential Equations: Some Black Boxes Were Meant to Remain Closed!

    pp. 8325–8335. 60A. Pal, A. Edelman, and C. Rackauckas, “Locally regularized neural dif- ferential equations: Some black boxes were meant to remain closed!” (2023), arXiv:2303.02262 [cs.LG]. 61T. Kumara, A. Kumara, and P. Pala, “A physics-constrained neural ordi- nary differential equations approach for robust learning of stiff chemical kinetics,” arXiv p...

  7. [17]

    Rows correspond to equations for y′ 0 to y′ 9, and columns represent coefficients for y0 to y9 and b. Neural ODEs for Stiff Systems: Explicit Exponential Integration Methods 14 y0 y1 y2 y3 y4 y5 y6 y7 y8 y9 b y′ 0 -10.00000012 5.00000255 -4.9835e-06 8.1128e-06 -5.9082e-05 1.6135e-04 -1.9551e-03 7.2901e-03 -1.2028e-02 9.9142e-03 6.3076e-07 y′ 1 4.99999616 ...

  8. [63]

    Discovering governing equa- tions from data by sparse identification of nonlinear dynamical systems,

    pioneered the success- ful training of stiff neural ODEs, showcasing robust training and accurate recovery of stiff dynamics using single-step im- plicit methods such as backward Euler, trapezoidal method, Radau3, and Radau5. These implicit approaches handle stiff- ness effectively but involve solving a nonlinear system at each time step, which can be com...

Show all 15 references
  1. [226]

    Stabilized neural ordinary differential equations for long-time forecasting of dynamical systems,

    44A. J. Linot, J. W. Burby, Q. Tang, P. Balaprakash, M. D. Graham, and R. Maulik, “Stabilized neural ordinary differential equations for long-time forecasting of dynamical systems,” Journal of Computational Physics474, 111838 (2023). 45S. I. Holt, Z. Qian, and M. van der Schaa...

  2. [1960]

    Performance evaluation of single-step explicit exponential integration methods on stiff ordinary differential equations,

    pp. 128–132. 75C. Fronk and L. Petzold, “Performance evaluation of single-step explicit exponential integration methods on stiff ordinary differential equations,” (2024), arXiv:2411.19374 [math.NA]. 76K. Hornik, M. B. Stinchcombe, and H. L. White, “Multilayer feedfor- ward net...

  3. [2003]

    High performance computing of the matrix exponential,

    pp. 111–120. 103P. Ruiz, J. Sastre, J. Ibáñez, and E. Defez, “High performance computing of the matrix exponential,” Journal of Computational and Applied Mathe- matics 291, 370–379 (2016), mathematical Modeling and Computational Methods. 104M. Arioli, B. Codenotti, and C. Fass...

  4. [2020]

    Neural controlled differen- tial equations for irregular time series,

    pp. 3870–3882. 18P. Kidger, J. Morrill, J. Foster, and T. Lyons, “Neural controlled differen- tial equations for irregular time series,” Advances in Neural Information Processing Systems 33, 6696–6707 (2020). 19P. Kidger, “On neural differential equations,” arXiv preprint arXi...

  5. [2021]

    Neural jump stochastic differential equations,

    pp. 7829–7838. 21J. Jia and A. R. Benson, “Neural jump stochastic differential equations,” Advances in Neural Information Processing Systems 32 (2019). 22R. T. Chen, B. Amos, and M. Nickel, “Learning neural event functions for ordinary differential equations,” arXiv preprint a...

  6. [2022]

    Proximal implicit ode solvers for accelerating learning neural odes,

    pp. 8811–8832. 46J. Baker, H. Xia, Y . Wang, E. Cherkaev, A. Narayan, L. Chen, J. Xin, A. L. Bertozzi, S. J. Osher, and B. Wang, “Proximal implicit ode solvers for accelerating learning neural odes,” arXiv preprint arXiv:2204.08621 (2022). 47R. Malpica Galassi, P. P. Ciottoli,...

  7. [2023]

    Generalized integrating factor methods for stiff pdes,

    pp. 372–388. 69S. Krogstad, “Generalized integrating factor methods for stiff pdes,” Jour- nal of Computational Physics 203, 72–88 (2005). 70B. Owren and A. Marthinsen, “Runge-kutta methods adapted to mani- folds and based on rigid frames,” BIT Numerical Mathematics 39, 116– 1...

Pith tools

Reviewed August 12, 2026 · model on record in the stance chip above.