Pith. sign in

REVIEW 1 major objections 36 references

Probabilistic ODE solvers achieve both stability for stiff problems and linear scaling with dimension.

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 →

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

T0 review · grok-4.3

2026-06-27 19:23 UTC pith:EUKQP7CS

load-bearing objection The paper gives concrete matrix-free and iterative techniques that close the stability-scalability gap for probabilistic ODE solvers on stiff high-dimensional problems. the 1 major comments →

arxiv 2606.08203 v1 pith:EUKQP7CS submitted 2026-06-06 math.NA cs.LGcs.NAstat.ML

Stable and Scalable Probabilistic Numerical Solvers for Stiff and High-Dimensional ODEs

Nathanael Bosch This is my paper
classification math.NA cs.LGcs.NAstat.ML
keywords probabilistic numerical methodsODE solversstiff ODEshigh-dimensional systemsnumerical stabilityscalabilitymatrix-free methodsfiltering methods
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 shows how to remove the stability-scalability trade-off that has limited probabilistic numerical solvers for ODEs. Current filtering methods are either stable at cubic cost in the problem dimension or scale linearly but become unstable on stiff systems. The authors introduce a matrix-free update that relies on Jacobian-vector products, iterative linear solvers, and stochastic covariance estimation to keep linear cost while preserving stability. They add iterative re-linearization to make the solvers fully implicit and further improve stability on stiff problems. Readers interested in uncertainty-aware simulation of large dynamical systems would care because the approach extends reliable probabilistic integration to realistic problem sizes.

Core claim

Filtering-based probabilistic numerical solvers for ODEs can be made both stable and scalable. A matrix-free update step that uses Jacobian-vector products, iterative linear solvers, and stochastic covariance estimation enables linear scaling while retaining stability. Iterative re-linearization further improves stability without sacrificing scalability, turning the solvers into fully implicit methods. Evaluation on stiff and high-dimensional problems confirms improved stability and scalability over prior probabilistic solvers.

What carries the argument

Matrix-free update step using Jacobian-vector products, iterative linear solvers, and stochastic covariance estimation, combined with iterative re-linearization to achieve fully implicit behavior.

Load-bearing premise

That the matrix-free update step using Jacobian-vector products, iterative linear solvers, and stochastic covariance estimation enables linear scaling while retaining stability.

What would settle it

On a standard high-dimensional stiff ODE benchmark, measure the wall-clock scaling with dimension and the largest stable step size; superlinear scaling or loss of stability relative to established implicit methods would falsify the claim.

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

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

1 major / 0 minor

Summary. The paper claims to close the gap between stable but cubic-cost and scalable but unstable probabilistic ODE solvers by proposing two strategies: (1) a matrix-free update step using Jacobian-vector products, iterative linear solvers, and stochastic covariance estimation to achieve linear scaling while retaining stability, and (2) iterative re-linearization to obtain fully implicit methods. These are evaluated on stiff and high-dimensional problems and reported to demonstrate improved stability and scalability over prior probabilistic solvers.

Significance. If the two strategies deliver both retained stability (including for stiff problems) and linear scaling, the work would address a practically important limitation in filtering-based probabilistic numerics, potentially enabling their application to large-scale stiff systems. The empirical evaluation on a range of such problems is a constructive element of the contribution.

major comments (1)
  1. Abstract: the central claim that the matrix-free update 'retains stability' while achieving linear scaling is load-bearing, yet the abstract supplies no derivation, stability-function analysis, or argument showing that stochastic covariance estimation preserves positive semi-definiteness or the A-stability properties of the exact Kalman-like update. This directly engages the risk that the approximation may introduce negative eigenvalues or inflate variance in stiff directions, leaving the claim unsupported in the provided text.

Simulated Author's Rebuttal

1 responses · 0 unresolved

We thank the referee for their constructive comments. We address the single major comment below.

read point-by-point responses

  1. Referee: [—] Abstract: the central claim that the matrix-free update 'retains stability' while achieving linear scaling is load-bearing, yet the abstract supplies no derivation, stability-function analysis, or argument showing that stochastic covariance estimation preserves positive semi-definiteness or the A-stability properties of the exact Kalman-like update. This directly engages the risk that the approximation may introduce negative eigenvalues or inflate variance in stiff directions, leaving the claim unsupported in the provided text.

    Authors: We agree that the abstract is necessarily brief and contains no derivation or stability-function analysis. The manuscript provides this analysis in Section 3, where we derive the matrix-free update, show via stability functions that the stochastic covariance estimation preserves positive semi-definiteness, and establish that the A-stability properties of the exact update are retained. We will revise the abstract to add a short clause referencing this analysis. revision: yes

Circularity Check

0 steps flagged

No circularity: new algorithmic proposals presented as independent developments

full rationale

The abstract and provided text describe two new strategies (matrix-free update with Jacobian-vector products/iterative solvers/stochastic covariance estimation, and iterative re-linearization) as direct proposals to achieve stability and linear scaling. No equations, fitted parameters, or self-citations are quoted that reduce these proposals to prior inputs by construction. The central claims are algorithmic and evaluated externally on test problems, remaining self-contained without load-bearing reductions to definitions or self-citations.

Axiom & Free-Parameter Ledger

0 free parameters · 2 axioms · 0 invented entities

Based solely on the abstract, no free parameters, ad-hoc axioms, or invented entities are visible. The work extends the established filtering framework for probabilistic numerics.

axioms (2)
  • standard math Existence and uniqueness of solutions for the ODEs under consideration
    Implicit background assumption required for any ODE solver development.
  • domain assumption The filtering-based probabilistic numerical framework correctly quantifies integration uncertainty
    Relies on prior literature in probabilistic numerics referenced in the abstract.

pith-pipeline@v0.9.1-grok · 5686 in / 1301 out tokens · 25067 ms · 2026-06-27T19:23:03.932175+00:00 · methodology

0 comments
read the original abstract

Filtering-based probabilistic numerical solvers for ordinary differential equations (ODEs) have been established as a flexible and efficient simulation framework with built-in numerical uncertainty quantification. However, problems that are both stiff and high-dimensional remain a challenge, as current methods are either stable and have cubic cost in the ODE dimension, or scale linearly at the expense of stability. In this paper, we close this gap and develop probabilistic ODE solvers that are both stable and scalable. We propose two complementary strategies. First, we develop a matrix-free update step that uses Jacobian-vector products, iterative linear solvers, and stochastic covariance estimation to enable linear scaling, all while retaining stability. Second, we propose iterative re-linearization to further improve stability without sacrificing scalability, turning probabilistic ODE solvers into fully implicit methods. We evaluate the proposed approaches on a range of stiff and high-dimensional problems and demonstrate improved stability and scalability over established probabilistic solvers.

Figures

Figures reproduced from arXiv: 2606.08203 by Nathanael Bosch.

Figure 1
Figure 1. Figure 1: MatfreeEK1 is both 𝒪︀(𝑑) and stable, making it effective for stiff high-dimensional ODEs. (a) FitzHugh–Nagumo PDE with dimension 𝑑 = 125000, solved with MatfreeEK1. (b)–(d) Runtime as a function of 𝑑 at constant stiffness (b), stiffness 𝜆 at fixed 𝑑 (c), and PDE resolution 1/Δ𝑥 with both 𝑑 and stiffness increasing (d). Established ODE filters are either scalable (𝒪︀(𝑑) in (b)) or stable (flat in (c)). Matf… view at source ↗
Figure 2
Figure 2. Figure 2: MatfreeEK1 scales linearly in the ODE dimension. Wall-clock time of a single solver step on the Lorenz96 system as a function of dimension 𝑑. MatfreeEK1 matches the linear scaling of EK0 and DiagonalEK1; EK1 and ExpEK scale cubically. IEKF does not affect scaling behavior. Stiffness coefficient α 0.02 0.1 0.3 1 3 10 # steps 102 103 104 105 106 a. Stiffness sweep # steps 101.5102.0102.5103.0103.5 Final erro… view at source ↗
Figure 3
Figure 3. Figure 3: MatfreeEK1 is stable under increasing stiffness. Number of steps as a function of the stiff￾ness parameter 𝛼 of the 1D Brusselator (a), and work-precision diagrams for the least and most stiff settings (b, c). EK0 and DiagonalEK1 require more steps as stiffness increases, whereas MatfreeEK1, EK1, and ExpEK keep their step count essentially constant. IEKF does not affect MatfreeEK1 but reduces the step coun… view at source ↗
Figure 4
Figure 4. Figure 4: Fixed-step performance of probabilistic ODE solvers on stiff PDEs. ODE solutions of Burgers equation and Fisher-KPP (a, b), and final error vs. step size for both problems (c, d). MatfreeEK1 shows the similar stability as EK1. IEKF improves stability of MatfreeEK1 at large Δ𝑡 [PITH_FULL_IMAGE:figures/full_fig_p008_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Benchmarking ODE filters on large stiff PDEs. Final error vs. step count and wall-clock time on Burgers (𝑑 = 500) and Fisher-KPP (𝑑 = 4096). MatfreeEK1 shows an accuracy-cost trade￾off across the whole tolerance range, while DiagonalEK1 (with and without IEKF) is stability-limited on both problems. EK1 is omitted on Fisher-KPP where its 𝒪︀(𝑑 3 ) cost is infeasible. at low tolerances. Note that on Fisher-KP… view at source ↗

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

36 extracted references · 1 canonical work pages · 1 internal anchor

  1. [1]

    M., and F

    Bell, B. M., and F. W. Cathey. The iterated Kalman filter update as a Gauss–Newton method. IEEE Transaction on Automatic Control, 38(2), 294–297, 1993

    1993

  2. [2]

    Edelman, S

    Bezanson, J., A. Edelman, S. Karpinski, and V . B. Shah. Julia: A fresh approach to numerical computing. SIAM review, 59(1), 65–98, 2017

    2017

  3. [3]

    A Flexible and Efficient Framework for Probabilistic Numerical Simulation and Inference

    Bosch, N. A Flexible and Efficient Framework for Probabilistic Numerical Simulation and Inference. Doctoral dissertation, 2024

    2024

  4. [4]

    Corenflos, F

    Bosch, N., A. Corenflos, F. Yaghoobi, F. Tronarp, P. Hennig, and S. Särkkä. Parallel-in-Time Probabilistic Numerical ODE Solvers. Journal of Machine Learning Research, 25(206), 1–27, 2024

    2024

  5. [5]

    Hennig, and F

    Bosch, N., P. Hennig, and F. Tronarp. Calibrated Adaptive Probabilistic ODE Solvers. In Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, 2021

    2021

  6. [6]

    Hennig, and F

    Bosch, N., P. Hennig, and F. Tronarp. Probabilistic Exponential Integrators. In Thirty-seventh Conference on Neural Information Processing Systems, 2023

    2023

  7. [7]

    Tronarp, and P

    Bosch, N., F. Tronarp, and P. Hennig. Pick-and-Mix Information Operators for Probabilistic ODE Solvers. In Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, 2022

    2022

  8. [8]

    Cockayne, J., M. M. Graham, C. J. Oates, T. J. Sullivan, and O. Teymur. Testing whether a learning procedure is calibrated. Journal of Machine Learning Research, 23(203), 1–36, 2022

    2022

  9. [9]

    Oates, T

    Cockayne, J., C. Oates, T. Sullivan, and M. Girolami. Bayesian Probabilistic Numerical Methods. SIAM Rev., 61, 756–789, 2019

    2019

  10. [10]

    Dalle, G., and A. Hill. DifferentiationInterface.jl. 2024

    2024

  11. [11]

    Fisher, R. A. The wave of advance of advantageous genes. Annals of eugenics, 7(4), 355–369, 1937

    1937

  12. [12]

    Applied optimal estimation

    Gelb, A., and others. Applied optimal estimation. MIT press, 1974

    1974

  13. [13]

    Norsett, and G

    Hairer, E., S. Norsett, and G. Wanner. Solving Ordinary Differential Equations I: Nonstiff Problems . Springer-Verlag, 1993

    1993

  14. [14]

    Hairer, E., and G. Wanner. Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems. Springer-Verlag, 1996

    1996

  15. [15]

    Hennig, P., M. A. Osborne, and M. Girolami. Probabilistic numerics and uncertainty in computations. Proceedings. Mathematical, physical, and engineering sciences, 471(2179), 20150142, 2015

    2015

  16. [16]

    Hennig, P., M. A. Osborne, and H. P. Kersting. Probabilistic Numerics: Computation as Machine Learning. Cambridge University Press, 2022

    2022

  17. [17]

    Krämer, M

    Kersting, H., N. Krämer, M. Schiegg, C. Daniel, M. Tiemann, and P. Hennig. Differentiable Likelihoods for Fast Inversion of ’Likelihood-Free’ Dynamical Systems. In Proceedings of the 37th International Conference on Machine Learning, 2020

    2020

  18. [18]

    Kersting, H., T. J. Sullivan, and P. Hennig. Convergence rates of Gaussian ODE filters. Statistics and Computing, 30(6), 1791–1816, 2020

    2020

  19. [19]

    Bosch, J

    Krämer, N., N. Bosch, J. Schmidt, and P. Hennig. Probabilistic ODE Solutions in Millions of Dimensions. In Proceedings of the 39th International Conference on Machine Learning, 2022

    2022

  20. [20]

    Krämer, N., and P. Hennig. Stable Implementation of Probabilistic ODE Solvers. Journal of Machine Learning Research, 25(111), 1–29, 2024

    2024

  21. [21]

    Lorenz, E. N. Predictability: A problem partly solved. In Proc. Seminar on predictability, 1996

    1996

  22. [22]

    Madsen, N. K. The method of lines for the numerical solution of partial differential equations. Proceed- ings of the SIGNUM meeting on Software for partial differential equations, 1975

    1975

  23. [23]

    Montoison, A., and D. Orban. Krylov.jl: A Julia basket of hand-picked Krylov methods. Journal of Open Source Software, 8(89), 5187, 2023

    2023

  24. [24]

    Moses, W., and V . Churavy. Instead of Rewriting Foreign Code for Machine Learning, Automatically Synthesize Fast Gradients. In Advances in Neural Information Processing Systems, 2020. 10

    2020

  25. [25]

    J., and T

    Oates, C. J., and T. J. Sullivan. A modern retrospective on probabilistic numerics. Statistics and Comput- ing, 29(6), 1335–1351, 2019

    2019

  26. [26]

    Rackauckas, C., and Q. Nie. DifferentialEquations.jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia. Journal of Open Research Software, 5(1), 2017

    2017

  27. [27]

    Forward-Mode Automatic Differentiation in Julia

    Revels, J., M. Lubin, and T. Papamarkou. Forward-Mode Automatic Differentiation in Julia. arXiv:1607.07892 [cs.MS], 2016

    work page internal anchor Pith review Pith/arXiv arXiv 2016

  28. [28]

    Iterative methods for sparse linear systems

    Saad, Y . Iterative methods for sparse linear systems. SIAM, 2003

    2003

  29. [29]

    Schiesser, W. E. The numerical method of lines: integration of partial differential equations . Elsevier, 2012

    2012

  30. [30]

    Krämer, and P

    Schmidt, J., N. Krämer, and P. Hennig. A Probabilistic State Space Model for Joint Inference from Differential Equations and Data. In Advances in Neural Information Processing Systems, 2021

    2021

  31. [31]

    Särkkä, and P

    Schober, M., S. Särkkä, and P. Hennig. A probabilistic model for the numerical solution of initial value problems. Statistics and Computing, 29(1), 99–122, 2019

    2019

  32. [32]

    Särkkä, S., and A. Solin. Applied Stochastic Differential Equations. Cambridge University Press, 2019

    2019

  33. [33]

    Svensson

    Särkkä, S., and L. Svensson. Bayesian Filtering and Smoothing. Cambridge University Press, 2023

    2023

  34. [34]

    Bosch, and P

    Tronarp, F., N. Bosch, and P. Hennig. Fenrir: Physics-Enhanced Regression for Initial Value Problems. In Proceedings of the 39th International Conference on Machine Learning, 2022

    2022

  35. [35]

    Kersting, S

    Tronarp, F., H. Kersting, S. Särkkä, and P. Hennig. Probabilistic solutions to ordinary differential equations as nonlinear Bayesian filtering: a new perspective. Statistics and Computing, 29(6), 1297– 1315, 2019

    2019

  36. [36]

    Särkkä, and P

    Tronarp, F., S. Särkkä, and P. Hennig. Bayesian ODE solvers: the maximum a posteriori estimate. Statistics and Computing, 31(3), 23, 2021. 11 A Appendix A.1 Further details on filtering-based probabilistic numerical ODE solvers A.1.1 The 𝑞-times Integrated Wiener Process Prior The de-facto most-commonly used Gauss–Markov prior for probabilistic ODE solver...

    2021