Pith. sign in

REVIEW 3 major objections 4 minor 47 references

Nonlinear ODE integration parallelizes to O(log N) span per Newton step via associative scans.

Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →

T0 review

2026-08-04 00:19 UTC pith:PMR76POM

load-bearing objection A clean algorithmic extension with sound Newton-scan theory, but the runtime claims rest on fixed iteration counts that aren't verified at the tested step sizes. the 3 major comments →

arxiv 2511.01465 v2 pith:PMR76POM submitted 2025-11-03 math.NA cs.NA

A parallel-in-time Newton's method-based ODE solver

classification math.NA cs.NA MSC 34A3465L0568W1065Y05
keywords ordinary differential equationsparallel-in-time integrationNewton's methodparallel prefix sumassociative scaninitial value problemGPURunge-Kutta
verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

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

This paper claims that the sequential time-stepping of a nonlinear ODE solver can be replaced by a parallel root-finding problem: write the rollout as a system of nonlinear equations h(ξ)=0 and solve it with Newton's method. The Newton step turns out to follow an affine recursion (3.17/3.47), and affine recursions are exactly what parallel prefix sums (associative scans) can evaluate in O(log N) span. If correct, this breaks the square-root complexity barrier of Parareal, giving an asymptotically faster parallel-in-time integrator for nonlinear initial value problems. The authors prove quadratic convergence of the Newton iterations and demonstrate wall-clock speedups on a GPU for several benchmark ODEs, both explicit and implicit.

Core claim

The paper establishes that an ODE rollout can be flattened into a nonlinear system h(ξ)=0, and that the Newton step obeys an affine recursion (u_k = (I + ∂g_{k-1}/∂x_{k-1})u_{k-1} - h_k for explicit methods, and a similar diagonal-inverse form for implicit methods, Eqs. 3.17 and 3.47). Affine maps compose associatively, so a parallel prefix sum (associative scan) computes the whole Newton step in O(log N) time. The paper proves quadratic convergence of the Newton iterates under standard assumptions and reports GPU experiments on logistic, van der Pol, cart-pole, Dahlquist, and Robertson systems where the method runs faster than Parareal and sequential integration.

What carries the argument

The mechanism is the affine recursion identity for the Newton step: for explicit schemes u_k = (I + ∂g_{k-1}/∂x_{k-1})u_{k-1} - h_k, and for implicit schemes u_k = -(∂g_{k-1}/∂x_k - I)^{-1}((-I - ∂g_{k-1}/∂x_{k-1})u_{k-1} + h_k). These are affine maps whose composition is associative, so a parallel scan (a prefix-sum algorithm for associative operators) computes all u_k in O(log N) span rather than the O(N) sequential sweep.

Load-bearing premise

The runtime comparisons in Section 4 assume that the fixed Newton iteration counts (11 explicit, 5 or 21 implicit) converge to the discrete ODE solution for every δt used, but residuals are only reported for the largest δt of each scheme (1e-2 explicit, 1e-1 implicit).

What would settle it

Run the benchmarks with a tolerance-based stopping criterion (e.g., ||h(ξ)||∞ < 1e-10) at the smallest δt in Figures 1 and 3 (δt=1e-5 explicit, δt=1e-4 Dahlquist), and record the actual iteration count and wall-clock time. If the iteration count grows with N, or the total time no longer scales as O(log N) per iteration, the claimed logarithmic advantage over Parareal would be falsified for those regimes.

Watch this falsifier. Get emailed when new claim-graph text bears on it.

Share X LinkedIn Reddit HN

If this is right

  • Per Newton iteration, the parallel span drops from linear to logarithmic, which is asymptotically faster than Parareal's square-root complexity.
  • The method applies to both explicit and implicit ODE solvers (including Runge–Kutta and backward Euler/trapezoidal) without changing the local stepping scheme.
  • Quadratic convergence, proven in Theorems 3.2 and 3.4, implies only a few iterations are needed near a good initial guess, as seen in the residual plots.
  • The solver produces the entire discrete trajectory ξ in a single parallel solve, so it is well-suited to problems that need all time points at once.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • Replacing the fixed iteration counts with a tolerance-based stopping rule and measuring how iterations grow with N and δt would test whether the reported runtimes correspond to accurate solves; the paper only verifies residuals at the largest δt of each scheme.
  • The affine structure could be exploited with blocked or hierarchical scan implementations, or fused with automatic differentiation of the Jacobian, to further cut memory traffic and communication on GPUs.
  • For stiff problems, the implicit scheme's requirement to invert (∂g/∂x_k - I) may become ill-conditioned for some step sizes, suggesting the method's practical range of δt depends on this non-singularity condition.

Editorial analysis

A structured set of objections, weighed in public.

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

Referee Report

3 major / 4 minor

Summary. The paper proposes a parallel-in-time solver for nonlinear ODEs. The discrete integration scheme (explicit or implicit) is written as a root-finding problem h(ξ)=0 for all time states, and Newton's method is applied. The key observation is that each Newton step is an affine recursion (Eqs. (3.17) and (3.47)), which can be evaluated with a parallel associative scan in O(log N) span. Local quadratic convergence is proven (Theorems 3.2 and 3.4). Numerical experiments on a GPU compare runtimes with Parareal and sequential integration for several benchmark ODEs, reporting average runtimes for fixed Newton iteration counts.

Significance. The central algorithmic idea is clean and potentially significant: it extends the affine-scan parallelization technique (previously used for linear problems) to nonlinear ODEs in a principled way. The derivation of the Newton-step recursion is correct, the convergence proofs are standard, and the paper includes a reproducible code base. If the runtime advantages hold under accuracy-normalized conditions, this would be a useful contribution to parallel-in-time integration. However, the headline claims of 'improved runtime' and 'logarithmic complexity' are not fully supported by the experiments as presented, because the fixed-iteration Newton solves are not verified for the step sizes used in the runtime plots. The theoretical contribution alone is sufficiently interesting, but the numerical evidence needs strengthening to justify the performance claims.

major comments (3)
  1. [§4.1 and §4.2, Figures 1–4] The runtime plots (Figures 1 and 3) use fixed Newton iteration counts (11 for explicit; 5 and 21 for implicit) for every step size, but residual/error decay is reported only for δt=1e-2 (explicit, Figure 2) and δt=1e-1 (implicit, Figure 4). For the smaller δt values used in the runtime plots (down to 1e-5 and 5e-3), the paper does not show that the fixed iteration count yields an accurate discrete solution. Since Newton's convergence is only local (Theorems 3.2 and 3.4) and the initial guesses are fixed (ones or zeros), it is possible that more iterations are needed as δt decreases or N grows. If so, the timed iteration counts do not represent the cost of an accurate solve, and the claimed end-to-end O(log N) span is not established. The authors should report, for every δt in the runtime plots, the residual norm after the fixed iteration count, or switch to a tolerance-based stopping rul
  2. [§4.1 and Appendix A (Parareal comparison)] The Parareal comparison in Figure 1 is not accuracy-normalized. The paper specifies M=N^{1/2} coarse steps and 11 iterations, and the fine integrator performs N/M sequential steps with the same δt, but the accuracy of the Parareal solution after those 11 iterations is never reported. Without knowing whether Parareal has converged (or conversely, whether it converged in far fewer iterations), the runtime comparison is potentially misleading. The authors should report the error or residual for Parareal for each δt, or compare runtimes required to reach a target tolerance for both methods.
  3. [Section 5 and Abstract (complexity claim)] The statement that the method has 'logarithmic span complexity' is per Newton step, not end-to-end. The total parallel span is O(K log N), where K is the number of Newton iterations. The paper does not discuss how K grows with N or with the required accuracy. If K remains constant for all tested step sizes, that should be stated explicitly with evidence; otherwise the asymptotic claim is incomplete and potentially misleading. The authors should qualify the complexity statement (e.g., 'per Newton iteration' or 'for a fixed number of iterations') and comment on the dependence of K on N and δt, ideally with supporting numerical evidence.
minor comments (4)
  1. [Eq. (3.12)] The last component of h(ξ) should be x_N - x_{N-1} - g_{N-1}(x_{N-1}, δt), not g_{N-1}(x_N, δt), to be consistent with Eq. (3.11).
  2. [Eq. (3.43)] The last component should be g_{N-1}(x_{N-1}, x_N, δt) — a comma is missing before δt.
  3. [Proof of Theorem 3.2] Equation (3.21) invokes 'Taylor's theorem' for a vector-valued function. The statement is the mean value theorem; please make the justification precise or replace with the integral remainder form.
  4. [§4.1] The text says for the logistic equation the residual drops by 8 orders in 5 iterations and by 16 orders in 10, but Figure 2 appears to show a decay to ~1e-16 after 10 iterations; quantifying this in the caption or text would improve clarity.

Circularity Check

0 steps flagged

No significant circularity: the Newton-step scan is derived in-text from first principles; self-citations are contextual, not load-bearing.

full rationale

The derivation chain is self-contained. The paper first rewrites the discrete ODE recursion as a root-finding problem h(ξ)=0 (Eqs. 3.12 and 3.43). It then computes the Newton step from the lower-bidiagonal Jacobian (Eqs. 3.13 and 3.44), obtaining the affine recurrences (3.17) and (3.47). These recurrences are then parallelized by the associative-composition formulas (3.5)-(3.10), whose composition rule F_{j,i}=F_{k,i}F_{j,k}, c_{j,i}=F_{k,i}c_{j,k}+c_{k,i} is derived in the paper rather than assumed. The convergence theorems (3.2 and 3.4) are standard Newton proofs using the stated Lipschitz and invertibility assumptions; they do not presuppose the runtime or uniqueness conclusions. The citations to [37] and [24] are contextual references to prior work on affine recursion and Newton-based NMPC, but the paper reproduces the needed derivation internally, and the scan primitive itself is attributed to the external Blelloch results. No fitted parameter is renamed as a prediction, no quantity is defined in terms of the quantity it is used to predict, and no central claim reduces by construction to an input. The runtime comparisons with Parareal are empirical benchmarks rather than derived predictions. The potential concern that fixed Newton iteration counts are not verified at every step size in the runtime plots is a correctness/verification risk, not a circularity.

Axiom & Free-Parameter Ledger

4 free parameters · 5 axioms · 0 invented entities

The mathematical derivation introduces no free parameters; the listed free parameters are experimental settings. The proofs rest on standard Newton smoothness and local-convergence assumptions and on the affine-scan composition rules. No new physical or mathematical entities are postulated.

free parameters (4)
  • Newton iteration count (explicit) = 11
    Fixed for all N and delta-t in runtime comparisons (Section 4.1); no convergence criterion, so runtimes may correspond to different accuracy levels.
  • Newton iteration count (implicit) = 5 (Dahlquist), 21 (Robertson)
    Fixed for all N and delta-t (Section 4.2); residual convergence is shown only for delta-t = 1e-1.
  • Initial guess xi^(0) = ones for logistic/van der Pol; zeros for cart-pole/Dahlquist/Robertson
    Chosen by hand; Theorems 3.2 and 3.4 only guarantee local convergence from a sufficiently close initial guess, which is not quantified.
  • Parareal coarse-step count M = N^(1/2)
    Selected as the standard optimal Parareal choice; affects the comparison baseline but not the method's own derivation.
axioms (5)
  • domain assumption Assumption 3.1: Jacobian H of the discretized residual is Lipschitz continuous near the solution xi*.
    Standard Newton smoothness; needed for quadratic convergence in Theorems 3.2 and 3.4; not verified in experiments beyond residual plots.
  • ad hoc to paper Initial guess xi^(0) lies in the basin of attraction of xi*.
    Theorem 3.2/3.4 require a 'sufficiently close' initial guess; the paper uses fixed zeros or ones and demonstrates convergence only for selected delta-t (Figures 2 and 4), not for all runtime configurations.
  • domain assumption Assumption 3.3: det(dg_{i-1}/dx_i - I) != 0 for the implicit discretization.
    Ensures the diagonal blocks of the implicit Jacobian are invertible; standard for backward Euler and implicit Runge-Kutta under reasonable step sizes.
  • domain assumption Parallel scan composition of affine maps (Eqs. 3.5-3.7) has O(log N) span, treating each per-node matrix operation as O(1).
    The logarithmic span claim counts matrix multiplications and additions as constant cost, i.e., fixed state dimension dx; for large dx the per-node cost affects practical scaling.
  • standard math The discrete nonlinear system h(xi)=0 has a solution xi* for the chosen discretizations, and f is locally Lipschitz.
    Standard existence, uniqueness, and consistency assumptions for the ODE and numerical methods used in Section 2.

reviewed 2026-08-04 · how reviews work

0 comments
Cite this review

Pith. "Pith review of A parallel-in-time Newton's method-based ODE solver." pith.science (2026). https://pith.science/paper/PMR76POM

@misc{pith2026251101465,
  author       = {Pith},
  title        = {Pith review of: A parallel-in-time Newton's method-based ODE solver},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/PMR76POM}},
  note         = {Machine review of arXiv:2511.01465}
}
Share X LinkedIn Reddit HN
read the original abstract

In this article, we introduce a novel parallel-in-time solver for nonlinear ordinary differential equations (ODEs). We state the numerical solution of an ODE as a root-finding problem that we solve using Newton's method. The affine recursive operations arising in Newton's step are parallelized in time by using parallel prefix sums, that is, parallel scan operations, which leads to a logarithmic span complexity. This yields an improved runtime compared to the previously proposed Parareal method. We demonstrate the computational advantage through numerical simulations of various systems of ODEs.

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

47 extracted references · 2 linked inside Pith

  1. [1]

    D. S. Bernstein , Matrix Mathematics: Theory, Facts, and Formulas , Princeton University Press, 2nd ed., 2009

  2. [2]

    Bhatt, L

    R. Bhatt, L. Debreu, and A. Vidard , Introducing time parallelization within data assimilation , SIAM Journal on Scientific Computing, 47 (2025), pp. B533--B557

  3. [3]

    G. E. Blelloch , Scans as primitive parallel operations , IEEE Transactions on Computers, 38 (1989), pp. 1526--1538

  4. [4]

    G. E. Blelloch , Prefix sums and their applications , tech. report, School of Computer Science, Carnegie Mellon University Pittsburgh, PA, USA, 1990

  5. [5]

    Bosch, A

    N. Bosch, A. Corenflos, F. Yaghoobi, F. Tronarp, P. Hennig, and S. S \"a rkk \"a , Parallel-in-time probabilistic numerical ODE solvers , Journal of Machine Learning Research, 25 (2024), pp. 1--27

  6. [6]

    J. P. Boyd , Chebyshev and Fourier Spectral Methods , Courier Corporation, 2001

  7. [7]

    Bradbury, R

    J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. Vander P las, S. Wanderman- M ilne, and Q. Zhang , JAX : composable transformations of P ython+ N um P y programs , 2018, http://github.com/jax-ml/jax

  8. [8]

    Buonomo , The periodic solution of van der P ol's equation , SIAM Journal on Applied Mathematics, 59 (1998), pp

    A. Buonomo , The periodic solution of van der P ol's equation , SIAM Journal on Applied Mathematics, 59 (1998), pp. 156--171

  9. [9]

    J. C. Butcher , Numerical Methods for Ordinary Differential Equations , John Wiley & Sons, 2016

  10. [10]

    R. M. Corless, C. Y. Kaya, and R. H. Moir , Optimal residuals and the D ahlquist test problem , Numerical Algorithms, 81 (2019), pp. 1253--1274

  11. [11]

    V. A. Dobrev, T. Kolev, N. A. Petersson, and J. B. Schroder , Two-level convergence theory for multigrid reduction in time ( MGRIT ) , SIAM Journal on Scientific Computing, 39 (2017), pp. S501--S527

  12. [12]

    R. D. Falgout, S. Friedhoff, T. V. Kolev, S. P. MacLachlan, and J. B. Schroder , Parallel time integration with multigrid , SIAM Journal on Scientific Computing, 36 (2014), pp. C635--C661

  13. [13]

    L. Fang, S. Vandewalle, and J. Meyers , A parallel-in-time multiple shooting algorithm for large-scale PDE -constrained optimal control problems , Journal of Computational Physics, 452 (2022), p. 110926

  14. [14]

    M. J. Gander , 50 years of time parallel time integration , in Multiple Shooting and Time Domain Decomposition Methods: MuS-TDD, Heidelberg, May 6-8, 2013, Springer, 2015, pp. 69--113

  15. [15]

    M. J. Gander and S. Gu\"uttel , PARAEXP : A parallel integrator for linear initial-value problems , SIAM Journal on Scientific Computing, 35 (2013), pp. C123--C142

  16. [16]

    M. J. Gander, F. Kwok, and J. Salomon , PARAOPT : A parareal algorithm for optimality systems , SIAM Journal on Scientific Computing, 42 (2020), pp. A2773--A2802

  17. [17]

    M. J. Gander, M. Ohlberger, and S. Rave , A parareal algorithm with spectral coarse solver , arXiv preprint arXiv:2508.08873, (2025)

  18. [18]

    M. J. Gander and S. Vandewalle , Analysis of the parareal time-parallel time-integration method , SIAM Journal on Scientific Computing, 29 (2007), pp. 556--578

  19. [19]

    Gattiglio, L

    G. Gattiglio, L. Grigoryeva, and M. Tamborrino , Rand N et- P arareal: a time-parallel PDE solver using random neural networks , Advances in Neural Information Processing Systems, 37 (2024), pp. 94993--95025

  20. [20]

    Gattiglio, L

    G. Gattiglio, L. Grigoryeva, and M. Tamborrino , Prob- GP arareal: A probabilistic numerical parallel-in-time solver for differential equations , arXiv preprint arXiv:2509.03945, (2025)

  21. [21]

    Hairer and G

    E. Hairer and G. Wanner , Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems , Springer, 2nd ed., 1996

  22. [22]

    Hairer, G

    E. Hairer, G. Wanner, and S. P. N rsett , Solving Ordinary Differential Equations I: Nonstiff Problems , Springer, 2nd ed., 1993

  23. [23]

    Haut and B

    T. Haut and B. Wingate , An asymptotic parallel-in-time method for highly oscillatory PDE s , SIAM Journal on Scientific Computing, 36 (2014), pp. A693--A713

  24. [24]

    Iacob, H

    C. Iacob, H. Abdulsamad, and S. S \"a rkk \"a , A parallel-in-time N ewton’s method for nonlinear model predictive control , IEEE Transactions on Control Systems Technology, (2025)

  25. [25]

    Iqbal, H

    S. Iqbal, H. Abdulsamad, T. Cator, U. Braga-Neto, and S. S \"a rkk \"a , Parallel-in-time probabilistic solutions for time-dependent nonlinear partial differential equations , in 2024 IEEE 34th International Workshop on Machine Learning for Signal Processing (MLSP), 2024, pp. 1--6

  26. [26]

    B. Jin, Q. Lin, and Z. Zhou , Optimizing coarse propagators in parareal algorithms , SIAM Journal on Scientific Computing, 47 (2025), pp. A735--A761

  27. [27]

    Legoll, T

    F. Legoll, T. Leli \`e vre, K. Myerscough, and G. Samaey , Parareal computation of stochastic differential equations with time-scale separation: a numerical convergence study , Computing and Visualization in Science, 23 (2020), p. 9

  28. [28]

    F. L. Lewis, D. Vrabie, and V. L. Syrmos , Optimal Control , John Wiley & Sons, 3rd ed., 2012

  29. [29]

    Lions, Y

    J.-L. Lions, Y. Maday, and G. Turinici , R \'e solution d' EDP par un sch \'e ma en temps ``parar \'e el'' , Comptes Rendus de l'Acad \'e mie des Sciences-Series I-Mathematics, 332 (2001), pp. 661--668

  30. [30]

    E. N. Lorenz , Deterministic nonperiodic flow , in Universality in Chaos, 2nd edition, Routledge, 2017, pp. 367--378

  31. [31]

    Maday, M.-K

    Y. Maday, M.-K. Riahi, and J. Salomon , Parareal in time intermediate targets methods for optimal control problems , in Control and optimization with PDE constraints, Springer, 2013, pp. 79--92

  32. [32]

    T. P. Mathew, M. Sarkis, and C. E. Schaerer , Analysis of block parareal preconditioners for parabolic optimal control problems , SIAM Journal on Scientific Computing, 32 (2010), pp. 1180--1200

  33. [33]

    Nocedal and S

    J. Nocedal and S. J. Wright , Numerical Optimization , Springer, 2006

  34. [34]

    Pentland, M

    K. Pentland, M. Tamborrino, and T. J. Sullivan , Error bound analysis of the stochastic parareal algorithm , SIAM Journal on Scientific Computing, 45 (2023), pp. A2657--A2678

  35. [35]

    Rader, T

    J. Rader, T. Lyons, and P. Kidger , Optimistix: modular optimisation in JAX and E quinox , arXiv preprint arXiv:2402.09983, (2024)

  36. [36]

    H. H. Robertson , The solution of a set of reaction rate equations , in Numerical Analysis: An Introduction, J. Walsh, ed., Academic Press, London, England, 1966, pp. 178--182

  37. [37]

    S \"a rkk \"a and \'A

    S. S \"a rkk \"a and \'A . F. Garc \' a-Fern \'a ndez , Temporal parallelization of dynamic programming and linear quadratic control , IEEE Transactions on Automatic Control, 68 (2022), pp. 851--866

  38. [38]

    S \"a rkk \"a and \'A

    S. S \"a rkk \"a and \'A . F. Garc \' a-Fern \'a ndez , Temporal parallelisation of the HJB equation and continuous-time linear quadratic control , IEEE Transactions on Automatic Control, (2024)

  39. [39]

    S \"a rkk \"a and L

    S. S \"a rkk \"a and L. Svensson , Bayesian Filtering and Smoothing , vol. 17, Cambridge University Press, 2023

  40. [40]

    Satish, M

    N. Satish, M. Harris, and M. Garland , Designing efficient sorting algorithms for manycore gpus , in 2009 IEEE International Symposium on Parallel & Distributed Processing, IEEE, 2009, pp. 1--10

  41. [41]

    R. S. Sutton and A. G. Barto , Reinforcement Learning: An Introduction , vol. 1, MIT Press Cambridge, 1998

  42. [42]

    Tedrake , Underactuated robotics

    R. Tedrake , Underactuated robotics . Course Notes for MIT 6.832, 2023, https://underactuated.csail.mit.edu

  43. [43]

    Tronarp, H

    F. Tronarp, H. Kersting, S. S \"a rkk \"a , and P. Hennig , Probabilistic solutions to ordinary differential equations as nonlinear B ayesian filtering: a new perspective , Statistics and Computing, 29 (2019), pp. 1297--1315

  44. [44]

    Yaghoobi, A

    F. Yaghoobi, A. Corenflos, S. Hassan, and S. S \"a rkk \"a , Parallel square-root statistical linear regression for inference in nonlinear state space models , SIAM Journal on Scientific Computing, 47 (2025), pp. B454--B476

  45. [45]

    Yaghoobi and S

    F. Yaghoobi and S. S \"a rkk \"a , Parallel state estimation for systems with integrated measurements , IEEE Signal Processing Letters, (2024)

  46. [46]

    Y. Yang, Y. Wu, and J. Pan , Parallel dynamics computation using prefix sum operations , IEEE Robotics and Automation Letters, 2 (2017), pp. 1296--1303

  47. [47]

    write newline

    " write newline "" before.all 'output.state := FUNCTION fin.entry add.period write newline FUNCTION new.block output.state before.all = 'skip after.block 'output.state := if FUNCTION not #0 #1 if FUNCTION and 'skip pop #0 if FUNCTION or pop #1 'skip if FUNCTION new.block.checka empty 'skip 'new.block if FUNCTION field.or.null duplicate empty pop "" 'skip ...

This paper was first reviewed by deepseek-v4-flash on August 4, 2026.