Pith. sign in

REVIEW 3 major objections 6 minor 31 references

Quantifying uncertainty in the numerical integration of evolution equations based on Bayesian isotonic regression

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

Pith's one-line read The paper claims that discretization errors in ODE solvers can be quantified as monotone error variances using Bayesian isotonic regression with a horseshoe-type prior, and that approximating the log-chi-square distribution by a…

desk verdict The Gibbs sampler as printed does not target the posterior it claims; the paper's core algorithm is wrong as written, but the modeling idea is solid and worth a serious referee. read the letter →

arxiv 2411.08338 v1 pith:PGN33AOH submitted 2024-11-13 math.NA cs.NAphysics.data-anstat.ME

classification math.NAcs.NAphysics.data-anstat.ME MSC 65L0562F1562G08
keywords discretizationerrorquantificationBayesianisotonicregressionhorseshoepriorGibbssamplerGaussianmixturemodelordinarydifferentialequationsprobabilisticnumericsuncertainty
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

This paper develops a Bayesian way to put error bars on the numerical solution of an ordinary differential equation when only noisy observations of the true solution are available. The key idea is to treat the discretization error at each observation time as a zero-mean Gaussian whose variance grows monotonically as integration proceeds, and to estimate those variances from the residuals between observations and numerical approximations. The authors show that a horseshoe-type shrinkage prior on the log-increments of the variances, together with a 10-component Gaussian mixture approximation of the log-chi-square distribution, turns the posterior into a Gibbs sampler whose full conditionals are all standard distributions. In numerical tests on the FitzHugh–Nagumo and Kepler systems, the resulting credible intervals for the absolute residual and error capture the observed scale of the error. If the framework is right, it gives practitioners a principled way to quantify and propagate discretization uncertainty in Bayesian inverse problems for ODEs.

What carries the argument

The engine is the variable transformation $\eta_1 = \log\sigma_1^2$ and $\eta_j = \log\sigma_j^2 - \log\sigma_{j-1}^2$, which turns the isotonic constraint $\gamma^2 \le \sigma_1^2 \le \cdots \le \sigma_n^2$ into lower-tail constraints on real-valued $\eta_j$'s. On these increments the paper puts a horseshoe-type scale-mixture prior: $\eta_1 \mid \tau_1 \sim N_{\ge\log\gamma^2}(\log\gamma^2,\tau_1)$ and $\eta_j \mid \tau_j,\lambda \sim N_{\ge 0}(0,\lambda\tau_j)$, with half-Cauchy marginals for $\sqrt{\tau_j}$ and $\sqrt{\lambda}$ giving strong shrinkage with occasional jumps. The second mechanism is the log-square trick: writing $z_i = \log r_i^2$ makes the cumulative $\eta$'s appear as location parameters in a regression of $z_i$ on log-$\chi^2_1$ noise. Replacing that noise density by a 10-component Gaussian mixture with tabulated weights, means, and variances lets the sampler treat the problem as conditionally Gaussian, so all Gibbs full conditionals become truncated normal, gamma, or generalized inverse Gaussian draws. This yields an $O(n^2)$-per-iteration sampling algorithm for the posterior of the discretization error variances.

What would settle it

Compute the actual discretization error sequence for the simple harmonic oscillator integrated with the Störmer–Verlet method over many periods, where the error variance oscillates rather than increasing monotonically. If the proposed method's 90% credible intervals cover the true error in substantially fewer than 90% of time points, the monotone-variance model is contradicted.

Watch

Extended reading notes

Core claim

The central claim is that discretization error variances $\sigma_i^2$ of a numerical ODE integrator can be estimated from observable residuals $r_i = v_i - O(x_i)$ by a Bayesian isotonic regression model, provided the variances are assumed nondecreasing and the observation noise variance $\gamma^2$ is known. Writing $\eta_1 = \log\sigma_1^2$ and $\eta_j = \log\sigma_j^2 - \log\sigma_{j-1}^2$ for $j\ge 2$ converts monotonicity into the simple constraints $\eta_1 \ge \log\gamma^2$ and $\eta_j \ge 0$. Placing a half-shrinkage (horseshoe) prior on the $\eta_j$'s induces strong shrinkage of adjacent variances toward equality while still allowing occasional sharp jumps. The log-squared residuals satisfy $z_i = \log r_i^2 = \eta_1 + \cdots + \eta_i + \log\varepsilon_i^2$, where $\log\varepsilon_i^2$ has the log-$\chi^2_1$ distribution; approximating that density by a 10-component Gaussian mixture with fixed weights, means, and variances makes every full conditional in the Gibbs sampler a standard truncated normal, gamma, or generalized inverse Gaussian distribution. The paper demonstrates on the FitzHugh–Nagumo model and a Kepler-orbit velocity observation that the resulting posterior means and 95% credible intervals for $\sigma_i$ and for $\sqrt{\sigma_i^2-\gamma^2}$ track the actual residual and error scales.

Load-bearing premise

The load-bearing premise is that discretization errors are independent, zero-mean Gaussian random variables whose variances never decrease; the paper itself flags this in Remark 2, and if the true error is correlated, biased, or has a decreasing phase, the posterior credible intervals will be miscalibrated.

Editorial extensions

If this is right

  • Posterior means and 95% credible intervals for each $\sigma_i$ (and for the pure discretization part $\sqrt{\sigma_i^2-\gamma^2}$) provide an uncertainty band around the numerical solution at every observation time.
  • Because all full conditionals are standard distributions, the Gibbs sampler is straightforward to implement, making error quantification practical when the ODE state dimension is large relative to the number of observations.
  • The construction extends Bayesian isotonic regression from estimating normal means to estimating normal variances, so it can be reused for other problems with a monotone variance structure.
  • In the FitzHugh–Nagumo and Kepler experiments, a large proportion of observed residuals and true errors fall inside the 90% credible intervals, indicating the estimated variances capture the error scale.

Reading between the lines

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

  • A natural next step is to couple this error quantification with Bayesian parameter estimation for the ODE itself, so that discretization uncertainty is propagated into posterior distributions of unknown parameters rather than being ignored.
  • The diagonal-covariance restriction could be lifted using a Wishart-type prior on full matrices, an extension the paper explicitly leaves for future work.
  • The fixed 10-component mixture approximation could be audited by comparing posterior intervals under more components or alternative approximations, and the posterior output could serve as a diagnostic for whether the monotone-variance assumption actually holds.
  • Because the method yields a full posterior over error variances, it could be used to choose step sizes adaptively by shrinking steps where the estimated variance is large, complementing traditional local error control.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

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

Referee Report

3 major / 6 minor

Summary. The manuscript proposes a Bayesian framework for quantifying discretization errors in numerical solutions of ordinary differential equations. The residuals between noisy observations and numerical approximations are modeled as independent zero-mean Gaussians with monotonically nondecreasing variances, and a horseshoe-type shrinkage prior is placed on the log-increments of these variances. To make posterior sampling tractable, the log-chi-square distribution arising from a log-squared transformation is approximated by a 10-component Gaussian mixture, and a Gibbs sampler is constructed. The method is demonstrated on the FitzHugh-Nagumo and Kepler models, where posterior credible intervals are compared with actual residual and error values.

Significance. The paper addresses an important problem in probabilistic numerics: giving uncertainty quantification for deterministic ODE discretization errors. Extending Bayesian isotonic regression from normal means to variances is a useful conceptual step, and the Gaussian-mixture approximation of the log-chi-square distribution is a sensible device for enabling conjugate-style updates. If the sampling algorithm is corrected, the framework could be practically valuable for inverse problems in ODE models. The authors are also explicit about the strong modeling assumptions (Gaussianity, independence, monotone variances), which is commendable. However, the central algorithmic claim is undermined by an incorrect full-conditional update for the mixture indicators, so the numerical results as presented cannot be used to validate the stated posterior.

major comments (3)
  1. [Section 2.3.2 and Algorithm 1] The update for the auxiliary mixture indicators s_i is not a valid Gibbs step. The paper states that s_i can be sampled from the prior probability mass (5), and Algorithm 1 step 4 draws s(t) ~ P(s). In the augmented model, the full conditional is p(s_i = k | z_i, eta_i) proportional to w_k N(z_i - eta_i; m_k, v_k^2), which depends on the data through z_i - eta_i. Because the subsequent update of eta_j in step 5-6 uses the current values of tau and lambda, the overall transition kernel depends on the current state in a way that does not have the intended posterior as its stationary distribution. Thus the sampler as printed does not target pi(Sigma, Upsilon | r), and the coverage results in Section 3 cannot be attributed to the stated posterior.
  2. [Section 2.3.2, formula for w^2_j] The variance in the conditional posterior of eta_j is printed as w^2_j = (tilde{w}^2_j lambda tau_1) / (tilde{w}^2_j + lambda tau_j), and the same expression is repeated in Algorithm 1 line 10. Since the prior for eta_j is N_{>=0}(0, lambda tau_j), the prior precision is 1/(lambda tau_j), so the correct posterior precision is 1/tilde{w}^2_j + 1/(lambda tau_j), giving w^2_j = (tilde{w}^2_j lambda tau_j) / (tilde{w}^2_j + lambda tau_j). The numerator should contain tau_j, not tau_1. This is a concrete formula error that would bias the sampling of eta_j if the algorithm is implemented as written.
  3. [Section 3, coverage assessment] The coverage checks in Figures 5-8 are in-sample: the same residuals are used to estimate the variance parameters and then to check whether the 90% credible intervals contain the residuals or errors. The text reports only that 'a large proportion of data points fall within the 90% credible intervals' for single trajectories, without a quantitative coverage rate. Because the model is fit to the same data being evaluated, the intervals are expected to show high coverage by construction. To support the central claim that the method accurately quantifies discretization error, the authors should report repeated-simulation coverage (e.g., the fraction of held-out errors covered across many trajectories) or otherwise assess out-of-sample behavior.
minor comments (6)
  1. [Section 2.3.2] The sentence 'Since the latent parameters s are independent of other parameters' is only true in the prior; in the posterior, s is correlated with eta and the data. This wording facilitates the incorrect update described in the first major comment.
  2. [Algorithm 1, line 10] The line reads 'Sample tau^{(t)}_1 ~ P(tau_j | eta^{(t)}_j, nu^{(t)}_j)'; the subscript on the sampled variable should be j, not 1.
  3. [Figure 1 caption] The caption contains a typo: 'the Gaussian mixture model g(epsilon) give in (4)' should read 'given in (4)'.
  4. [Remark 4] The phrase 'In such a case case' contains a duplicated word and should be 'In such a case'.
  5. [Section 2.3.1] The GIG parameter in the update for lambda is typeset as '-n + 2/2', which is ambiguous; it should be '-(n-2)/2'.
  6. [Section 3] No convergence diagnostics (trace plots, R-hat) are reported for the Gibbs sampler, and no code is released for the proposed method. Given the algorithmic error in the printed sampler, reproducible code and convergence checks would be essential for validation.

Circularity Check

0 steps flagged · score 2.0 of 10

No significant circularity: the Bayesian model is built from explicit priors and standard mixture approximations, and the in-sample coverage checks are model calibration, not disguised predictions.

full rationale

The central derivation chain is a Bayesian hierarchical model: residuals r_i = v_i - O(x_i) are modeled as N(0, sigma_i^2) with sigma_i^2 = gamma^2 + V_i (Section 2.1); the monotone variances are parameterized by log-increments eta_j with a horseshoe-type shrinkage prior adapted from Okano et al. [16]; and the log-chi-square density is replaced by the published 10-component Gaussian mixture of Omori et al. [29] (Section 2.3.2). These are explicit statistical inputs (prior and approximation), not outputs of the inference, so the posterior is defined by Bayes theorem rather than being equivalent to its inputs by construction. The numerical checks in Section 3 plot posterior credible intervals for the absolute residual |r_i| and for the discretization error sqrt(sigma_i^2 - gamma^2), then observe that a large proportion of the data points fall inside. These are in-sample posterior predictive checks, not out-of-sample predictions, so the intervals are not 'predictions' in the sense of being statistically forced by construction; moreover, the discretization error itself (x(t_i)-x_i) is not directly fed into the likelihood, giving the error check partial diagnostic content. Self-citations [12-14] and [16] (with co-author Irie) are used for background, motivation, and prior choice, but they are not load-bearing justifications for the central claim; Remark 2 explicitly acknowledges the Gaussian zero-mean assumption as a simplification. Two non-circular correctness issues should be noted: Section 2.3.2 and Algorithm 1 sample s_i from the prior P(s) rather than the conditional posterior p(s_i|z,eta), and the printed variance formula for eta_j has tau_1 in the numerator where tau_j is required. These may affect the numerical results but are not circularity. Overall, no step reduces by construction to its input, so the circularity score is low.

Assumptions & free parameters 2 free parameters · 5 assumptions · 0 invented entities

The model relies on a small number of structural assumptions: independent Gaussian discretization errors, known diagonal observation covariance, monotone variances, and the accuracy of the 10-component Gaussian mixture approximation. The variance sequence sigma_i^2 is the inferential target and is fitted to residuals; no additional physically invented entities are introduced.

free parameters (2)
  • sigma_i^2 (i=1,...,n) = posterior means displayed in Figures 3-8
    Discretization error variances at each observation time are the target parameters, estimated from the residuals r_i.
  • Prior hyperparameters (shapes and rates of Gamma priors on nu_j and xi) = fixed at 1/2 or 1
    Chosen by hand to implement the half shrinkage (horseshoe) prior; not fitted to data but influence the posterior shrinkage level.
assumptions (5)
  • domain assumption Discretization errors O(x(t_i))-O(x_i) are independent zero-mean Gaussian with variance V_i.
    This is the statistical model for deterministic discretization error; the authors note in Remark 2 it is adopted for simplicity and may not capture true error dynamics.
  • domain assumption The variances satisfy 0 < gamma^2 <= sigma_1^2 <= sigma_2^2 <= ... <= sigma_n^2.
    Monotonicity of the total variances is imposed through the prior support; if the true discretization error variance decreases at some time, the estimator is restricted.
  • domain assumption Observation noise covariance Gamma is known and diagonal, and the observation operator is handled directly.
    The likelihood treats Gamma as fixed and Sigma_i as diagonal; the univariate case is developed in full and multivariate extension is left for future work.
  • domain assumption The 10-component Gaussian mixture g(epsilon_tilde) accurately approximates the log-chi-square_1 density f(epsilon_tilde).
    The Gibbs sampler is based on this approximation; its accuracy is illustrated in Figure 1 but not formally bounded.
  • standard math The Gibbs sampler converges to the target posterior, and GIG variate generation is reliable in the edge cases discussed in Remark 4.
    No convergence diagnostics or formal ergodicity results are provided; the validity of the algorithm rests on standard MCMC theory.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Quantifying uncertainty in the numerical integration of evolution equations based on Bayesian isotonic regression." pith.science (2026). https://pith.science/paper/PGN33AOH

@misc{pith2026241108338,
  author       = {Pith},
  title        = {Pith review of: Quantifying uncertainty in the numerical integration of evolution equations based on Bayesian isotonic regression},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/PGN33AOH}},
  note         = {Machine review of arXiv:2411.08338}
}
abstract

This paper presents a new Bayesian framework for quantifying discretization errors in numerical solutions of ordinary differential equations. By modelling the errors as random variables, we impose a monotonicity constraint on the variances, referred to as discretization error variances. The key to our approach is the use of a shrinkage prior for the variances coupled with variable transformations. This methodology extends existing Bayesian isotonic regression techniques to tackle the challenge of estimating the variances of a normal distribution. An additional key feature is the use of a Gaussian mixture model for the $\log$-$\chi^2_1$ distribution, enabling the development of an efficient Gibbs sampling algorithm for the corresponding posterior.

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

31 extracted references · 21 canonical work pages

  1. [1]

    II: Stiff and Differential-Algebraic Problems

    Hairer, E., Wanner, G.: Solving Ordinary Differential Equations. II: Stiff and Differential-Algebraic Problems. Springer Series in Computational Mathematics, vol. 14. Springer, Berlin (1991). https://doi.org/10.1007/978-3-662-09947-6

  2. [2]

    I: 16 Nonstiff Problems, 2nd edn

    Hairer, E., Nørsett, S.P., Wanner, G.: Solving Ordinary Differential Equations. I: 16 Nonstiff Problems, 2nd edn. Springer Series in Computational Mathematics, vol

  3. [3]

    Conrad, P.R., Girolami, M., S¨ arkk¨ a, S., Stuart, A., Zygalakis, K.: Statistical analysis of differential equations: introducing probability measures on numeri- cal solutions. Stat. Comput. 27(4), 1065–1082 (2017) https://doi.org/10.1007/ s11222-016-9671-0

  4. [4]

    Kersting, H., Sullivan, T.J., Hennig, P.: Convergence rates of Gaussian ODE filters. Stat. Comput. 30(6), 1791–1816 (2020) https://doi.org/10.1007/ s11222-020-09972-4

  5. [5]

    Schober, M., S¨ arkk¨ a, S., Hennig, P.: A probabilistic model for the numerical solution of initial value problems. Stat. Comput. 29(1), 99–122 (2019) https: //doi.org/10.1007/s11222-017-9798-7

  6. [6]

    Tronarp, F., Kersting, H., S¨ arkk¨ a, S., Hennig, P.: Probabilistic solutions to ordi- nary differential equations as nonlinear Bayesian filtering: a new perspective. Stat. Comput. 29(6), 1297–1315 (2019) https://doi.org/10.1007/s11222-019-09900-1

  7. [7]

    Tronarp, F., S¨ arkk¨ a, S., Hennig, P.: Bayesian ODE solvers: the maximum a posteriori estimate. Stat. Comput. 31(3), 23–18 (2021) https://doi.org/10.1007/ s11222-021-09993-7

  8. [8]

    Springer, Berlin (1993)

Show all 31 references
  1. [9]

    Abdulle, A., Garegnani, G.: Random time step probabilistic methods for uncer- tainty quantification in chaotic and geometric numerical integration. Stat. Comput. 30(4), 907–932 (2020) https://doi.org/10.1007/s11222-020-09926-w

  2. [10]

    Lie, H.C., Stuart, A.M., Sullivan, T.J.: Strong convergence rates of probabilistic integrators for ordinary differential equations. Stat. Comput. 29(6), 1265–1283 (2019) https://doi.org/10.1007/s11222-019-09898-6

  3. [11]

    Calcolo 59(1), 13–33 (2022) https://doi.org/10.1007/s10092-022-00457-6

    Lie, H.C., Stahn, M., Sullivan, T.J.: Randomised one-step time integration meth- ods for deterministic operator differential equations. Calcolo 59(1), 13–33 (2022) https://doi.org/10.1007/s10092-022-00457-6

  4. [12]

    Cambridge University Press, Cambridge (2022)

    Hennig, P., Osborne, M.A., Kersting, H.P.: Probabilistic Numerics—Computation as Machine Learning. Cambridge University Press, Cambridge (2022). https:// doi.org/10.1017/9781316681411

  5. [13]

    Marumo, N., Matsuda, T., Miyatake, Y.: Modelling the discretization error of initial value problems using the Wishart distribution. Appl. Math. Lett. 147, 108833–9 (2024) https://doi.org/10.1016/j.aml.2023.108833

  6. [14]

    SIAM/ASA J

    Matsuda, T., Miyatake, Y.: Estimation of ordinary differential equation models with discretization error quantification. SIAM/ASA J. Uncertain. Quantif. 9(1), 302–331 (2021) https://doi.org/10.1137/19M1278405 17

  7. [15]

    Matsuda, T., Miyatake, Y.: Piecewise monotone estimation in one-parameter exponential family (2023)

  8. [16]

    Oates, C.J., Cockayne, J., Aykroyd, R.G., Girolami, M.: Bayesian probabilistic numerical methods in time-dependent state estimation for industrial hydrocyclone equipment. J. Amer. Statist. Assoc. 114(528), 1518–1531 (2019) https://doi.org/ 10.1080/01621459.2019.1574583

  9. [17]

    Okano, R., Hamura, Y., Irie, K., Sugasawa, S.: Locally adaptive Bayesian isotonic regression using half shrinkage priors. Scand. J. Statist. 51(1), 109–141 (2024) https://doi.org/10.1111/sjos.12676

  10. [18]

    Cai, B., Dunson, D.B.: Bayesian multivariate isotonic regression splines: Appli- cations to carcinogenicity studies. J. Amer. Statist. Assoc. 102(480), 1158–1171 (2007) https://doi.org/10.1198/016214506000000942

  11. [19]

    Curtis, S.M., Ghosh, S.K.: A variable selection approach to monotonic regression with Bernstein polynomials. J. Appl. Stat. 38(5), 961–976 (2011) https://doi. org/10.1080/02664761003692423

  12. [20]

    Bio- metrics 60(2), 398–406 (2004) https://doi.org/10.1111/j.0006-341X.2004.00184

    Neelon, B., Dunson, D.B.: Bayesian isotonic regression and trend analysis. Bio- metrics 60(2), 398–406 (2004) https://doi.org/10.1111/j.0006-341X.2004.00184. x

  13. [21]

    Shively, T.S., Sager, T.W., Walker, S.G.: A Bayesian approach to non-parametric monotone function estimation. J. R. Stat. Soc. Ser. B Stat. Methodol. 71(1), 159–175 (2009) https://doi.org/10.1111/j.1467-9868.2008.00677.x

  14. [22]

    In: Artificial Intelligence and Statistics, pp

    Carvalho, C.M., Polson, N.G., Scott, J.G.: Handling sparsity via the horseshoe. In: Artificial Intelligence and Statistics, pp. 73–80 (2009). PMLR

  15. [23]

    Biometrika 97(2), 465–480 (2010) https://doi.org/10.1093/biomet/asq017

    Carvalho, C.M., Polson, N.G., Scott, J.G.: The horseshoe estimator for sparse sig- nals. Biometrika 97(2), 465–480 (2010) https://doi.org/10.1093/biomet/asq017

  16. [24]

    Dagpunar, J.S.: An easily implemented generalised inverse Gaussian generator. Comm. Statist. Simulation Comput. 18(2), 703–710 (1989) https://doi.org/10. 1080/03610918908812785

  17. [25]

    Diploma Thesis, Technical University Graz, Austria (1989)

    Lehner, K.: Erzeugung von Zufallszahlen aus zwei exotischen Verteilungen. Diploma Thesis, Technical University Graz, Austria (1989)

  18. [26]

    H¨ ormann, W., Leydold, J.: Generating generalized inverse Gaussian ran- dom variates. Stat. Comput. 24(4), 547–557 (2014) https://doi.org/10.1007/ s11222-013-9387-3

  19. [27]

    PhD thesis, Massachusetts Institute of Technology (1988) 18

    Nelson, D.B.: The time series behavior of stock market volatility and returns. PhD thesis, Massachusetts Institute of Technology (1988) 18

  20. [28]

    Harvey, A., Ruiz, E., Shephard, N.: Multivariate stochastic variance models. Rev. Econ. Stud. 61(2), 247–264 (1994) https://doi.org/10.2307/2297980

  21. [29]

    Kim, S., Shephard, N., Chib, S.: Stochastic volatility: likelihood inference and comparison with ARCH Models. Rev. Econ. Stud. 65(3), 361–393 (1998) https: //doi.org/10.1111/1467-937X.00050

  22. [30]

    Omori, Y., Chib, S., Shephard, N., Nakajima, J.: Stochastic volatility with lever- age: Fast and efficient likelihood inference. J. Econom. 140(2), 425–449 (2007) https://doi.org/10.1016/j.jeconom.2006.07.008

  23. [31]

    arXiv preprint arXiv:2404.13986 (2024) 19

    Hiraki, D., Chib, S., Omori, Y.: Stochastic volatility in mean: Efficient analysis by a generalized mixture sampler. arXiv preprint arXiv:2404.13986 (2024) 19

Pith tools

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