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 →
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 →
A parallel-in-time Newton's method-based ODE solver
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
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.
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
- 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.
Referee Report
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)
- [§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
- [§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.
- [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)
- [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).
- [Eq. (3.43)] The last component should be g_{N-1}(x_{N-1}, x_N, δt) — a comma is missing before δt.
- [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.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
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
free parameters (4)
- Newton iteration count (explicit) =
11
- Newton iteration count (implicit) =
5 (Dahlquist), 21 (Robertson)
- Initial guess xi^(0) =
ones for logistic/van der Pol; zeros for cart-pole/Dahlquist/Robertson
- Parareal coarse-step count M =
N^(1/2)
axioms (5)
- domain assumption Assumption 3.1: Jacobian H of the discretized residual is Lipschitz continuous near the solution xi*.
- ad hoc to paper Initial guess xi^(0) lies in the basin of attraction of xi*.
- domain assumption Assumption 3.3: det(dg_{i-1}/dx_i - I) != 0 for the implicit discretization.
- 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).
- standard math The discrete nonlinear system h(xi)=0 has a solution xi* for the chosen discretizations, and f is locally Lipschitz.
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}
}
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.
Reference graph
Works this paper leans on
-
[1]
D. S. Bernstein , Matrix Mathematics: Theory, Facts, and Formulas , Princeton University Press, 2nd ed., 2009
2009
-
[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
2025
-
[3]
G. E. Blelloch , Scans as primitive parallel operations , IEEE Transactions on Computers, 38 (1989), pp. 1526--1538
1989
-
[4]
G. E. Blelloch , Prefix sums and their applications , tech. report, School of Computer Science, Carnegie Mellon University Pittsburgh, PA, USA, 1990
1990
-
[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
2024
-
[6]
J. P. Boyd , Chebyshev and Fourier Spectral Methods , Courier Corporation, 2001
2001
-
[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
2018
-
[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
1998
-
[9]
J. C. Butcher , Numerical Methods for Ordinary Differential Equations , John Wiley & Sons, 2016
2016
-
[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
2019
-
[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
2017
-
[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
2014
-
[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
2022
-
[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
2013
-
[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
2013
-
[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
2020
-
[17]
M. J. Gander, M. Ohlberger, and S. Rave , A parareal algorithm with spectral coarse solver , arXiv preprint arXiv:2508.08873, (2025)
Pith/arXiv arXiv 2025
-
[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
2007
-
[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
2024
-
[20]
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)
arXiv 2025
-
[21]
Hairer and G
E. Hairer and G. Wanner , Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems , Springer, 2nd ed., 1996
1996
-
[22]
Hairer, G
E. Hairer, G. Wanner, and S. P. N rsett , Solving Ordinary Differential Equations I: Nonstiff Problems , Springer, 2nd ed., 1993
1993
-
[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
2014
-
[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)
2025
-
[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
2024
-
[26]
B. Jin, Q. Lin, and Z. Zhou , Optimizing coarse propagators in parareal algorithms , SIAM Journal on Scientific Computing, 47 (2025), pp. A735--A761
2025
-
[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
2020
-
[28]
F. L. Lewis, D. Vrabie, and V. L. Syrmos , Optimal Control , John Wiley & Sons, 3rd ed., 2012
2012
-
[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
2001
-
[30]
E. N. Lorenz , Deterministic nonperiodic flow , in Universality in Chaos, 2nd edition, Routledge, 2017, pp. 367--378
2017
-
[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
2013
-
[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
2010
-
[33]
Nocedal and S
J. Nocedal and S. J. Wright , Numerical Optimization , Springer, 2006
2006
-
[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
2023
-
[35]
J. Rader, T. Lyons, and P. Kidger , Optimistix: modular optimisation in JAX and E quinox , arXiv preprint arXiv:2402.09983, (2024)
Pith/arXiv arXiv 2024
-
[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
1966
-
[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
2022
-
[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)
2024
-
[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
2023
-
[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
2009
-
[41]
R. S. Sutton and A. G. Barto , Reinforcement Learning: An Introduction , vol. 1, MIT Press Cambridge, 1998
1998
-
[42]
Tedrake , Underactuated robotics
R. Tedrake , Underactuated robotics . Course Notes for MIT 6.832, 2023, https://underactuated.csail.mit.edu
2023
-
[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
2019
-
[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
2025
-
[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)
2024
-
[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
2017
-
[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.
discussion (0)
Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.