pith. sign in

arxiv: 2111.10947 · v3 · submitted 2021-11-22 · 🧮 math.NA · cs.NA· stat.CO

Comparison of Numerical Solvers for Differential Equations for Holonomic Gradient Method in Statistics

Nobuki Takayama , Takaharu Yaguchi , Yi Zhang This is my paper

Pith reviewed 2026-05-24 12:57 UTC · model grok-4.3

classification 🧮 math.NA cs.NAstat.CO
keywords holonomic gradient methodnumerical ODE solversnormalizing constantsdefinite integralslinear differential equationsstatisticsholonomic functionsparameter restriction
0
0 comments X

The pith

Numerical solvers for linear ODEs arising from holonomic systems are compared for accuracy in evaluating statistical normalizing constants.

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

The paper shows how definite integrals of holonomic functions, which satisfy linear partial differential equations, reduce to ordinary differential equations when parameters are restricted to a one-dimensional curve. Solving these ODEs numerically then yields the value of the original integral. This reduction is the basis of the holonomic gradient method, already used to compute normalizing constants that appear in several statistical models. The authors therefore compare different numerical integration techniques for the resulting ODEs, focusing on their performance when the goal is to recover those constants. A sympathetic reader would care because the choice of solver directly affects whether the method can be applied reliably to concrete statistical problems.

Core claim

Definite integrals with parameters of holonomic functions satisfy holonomic systems of linear partial differential equations. Restricting the parameters to a one-dimensional curve converts the system into a linear ordinary differential equation. The integral is recovered by numerically solving this ODE. The holonomic gradient method applies this procedure to normalizing constants in statistics, and the paper compares several numerical solvers for the ODE step.

What carries the argument

The holonomic gradient method, which evaluates a definite integral by numerically integrating the linear ODE obtained after restricting a holonomic PDE system to a curve in parameter space.

If this is right

  • The comparison identifies which ODE solvers maintain sufficient accuracy for practical computation of normalizing constants.
  • Efficiency differences among solvers determine which can be used for higher-dimensional statistical models.
  • Numerical stability of the chosen solver becomes a necessary condition for reliable use of the holonomic gradient method.
  • The method can be applied to additional normalizing constants once a suitable solver is selected.

Where Pith is reading between the lines

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

  • If one solver class consistently outperforms others on statistical examples, software implementations of the holonomic gradient method could standardize on that class.
  • The same comparison framework might be tested on integrals outside statistics where holonomic functions appear.
  • Error bounds derived from the ODE solvers could be combined with the holonomic system to give rigorous guarantees on the integral approximation.

Load-bearing premise

Restricting parameters to a one-dimensional curve yields a linear ODE whose numerical solution recovers the original definite integral without large discretization or truncation errors.

What would settle it

Apply the holonomic gradient method with each solver to a normalizing constant whose exact value is known independently, such as the integral defining the multivariate gamma function, and check whether the computed values agree with the known value to within expected numerical tolerance.

Figures

Figures reproduced from arXiv: 2111.10947 by Nobuki Takayama, Takaharu Yaguchi, Yi Zhang.

Figure 1
Figure 1. Figure 1: is a graph of Airy Ai function and Airy Bi function. The function F(t) = (Ai(t), Ai′ (t))T satisfies the condition 2 of the Situation 1 of the instability problem [PITH_FULL_IMAGE:figures/full_fig_p004_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: Solving initial value problem of Airy differential equation by de [PITH_FULL_IMAGE:figures/full_fig_p007_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Solving the Airy differential equation by chebfun [PITH_FULL_IMAGE:figures/full_fig_p011_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: The Runge-Kutta(yellow green), the implicit Runge-Kutta [PITH_FULL_IMAGE:figures/full_fig_p014_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: The implicit Runge-Kutta method starting at [PITH_FULL_IMAGE:figures/full_fig_p014_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Solving by the sparse interpolation A [PITH_FULL_IMAGE:figures/full_fig_p015_6.png] view at source ↗
Figure 8
Figure 8. Figure 8: Solving a boundary value problem by chebfun. Left: [PITH_FULL_IMAGE:figures/full_fig_p016_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: Robustness of chebfun. H10 1 (1, y) on [104 , 104 + 40] with randomly perturbed boundary values. The values should be magnified by 1082 . Example 17 (Boundary value problem for Hk n(x, y) for x = 1 and y ∈ [108 , 108 + 2 × 105 ].) We give the boundary values of H10 1 (1, y) and ∂H10 1 ∂y (1, y) at y = 108 and y = 108 + 2×105 . We apply the chebfun package for this boundary value problem. (defusing/Hkn/y202… view at source ↗
Figure 11
Figure 11. Figure 11: The relative error of H10 1 (1, y) of the defusing method. The relative error is defined as (Hd − H)/H where Hd is the value by the defusing method and H is the exact value. = 1 2 Z ∞ t dσ Z ∞ −∞ db Z 2π 0 dθ Z 2π 0 dφ(σ 2 − b 2 ) s1s2 (2π) 2 expn − 1 2 R o , (41) where R = s1 (b sin θ sin φ + σ cos θ cos φ − m11) 2 + s2 (σ sin θ cos φ − b cos θ sin φ − m21) 2 +s1 (σ cos θ sin φ − b sin θ cos φ) 2 + s2 (b… view at source ↗
Figure 12
Figure 12. Figure 12: The graph of F29(t) and simulation values in the left and relative errors in the right. The data points are marked with ’x’ [PITH_FULL_IMAGE:figures/full_fig_p019_12.png] view at source ↗
Figure 13
Figure 13. Figure 13: The graph of F29(t)’s with random errors for data points and simulation values and relative errors in the right. by a Monte-Carlo simulation. (ec/tryb6/yiex5b.r) (pi , qi)’s are data points for the sparse ex￾trapolation method B. We divide the interval [ts, te] = [3.8, 4.0] into 200 segments and use the trapezoidal formula for ℓ(f). The optimal solution [f0, f1, . . . , f29] = [0.060145405402867516, −1.20… view at source ↗
read the original abstract

Definite integrals with parameters of holonomic functions satisfy holonomic systems of linear partial differential equations. When we restrict parameters to a one dimensional curve, the system becomes a linear ordinary differential equation (ODE) with respect to a curve in the parameter space. We can evaluate the integral by solving the linear ODE numerically. This approach to evaluate numerically definite integrals is called the holonomic gradient method (HGM) and it is useful to evaluate several normalizing constants in statistics. We will discuss and compare methods to solve linear ODE's to evaluate normalizing constants.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

0 major / 2 minor

Summary. The manuscript describes how definite integrals of holonomic functions satisfy holonomic PDE systems that reduce to linear ODEs when parameters are restricted to a one-dimensional curve in parameter space. It proposes using numerical solution of these ODEs (the holonomic gradient method, HGM) to evaluate normalizing constants in statistics and states that it will discuss and compare methods for solving the resulting linear ODEs.

Significance. If the comparisons identify reliable, efficient solvers with controlled truncation and discretization error for HGM integrals, the work could supply practical guidance for statisticians computing normalizing constants of holonomic distributions. The approach leverages an existing theoretical framework rather than introducing new theory, so its value lies in the concrete numerical comparison and any accompanying code or error bounds.

minor comments (2)
  1. The abstract states the plan to compare solvers but supplies no concrete results, tables, or error metrics; the full manuscript should include at least one representative numerical example with reported accuracy and timing to substantiate the comparison claim.
  2. Notation for the curve restriction and the resulting ODE coefficients is not introduced in the provided abstract; the manuscript should define these quantities explicitly (e.g., the parametrization of the curve and the matrix of the linear system) before discussing solvers.

Simulated Author's Rebuttal

0 responses · 0 unresolved

We thank the referee for their summary and assessment of our manuscript. No specific major comments were provided in the report, so we have no point-by-point responses to address. We remain available to clarify any aspects of the work or provide additional details if the editor or referee requests them.

Circularity Check

0 steps flagged

No significant circularity

full rationale

The paper is a methods-comparison study of numerical solvers for linear ODEs that arise when parameters are restricted to a one-dimensional curve in the holonomic gradient method. No derivation chain, prediction, or fitted quantity is presented that reduces to its own inputs by construction, self-definition, or self-citation load-bearing. The abstract and content describe a comparison task without equations or results that are statistically forced or renamed known patterns. The work is self-contained against external benchmarks for solver performance.

Axiom & Free-Parameter Ledger

0 free parameters · 0 axioms · 0 invented entities

Abstract-only review provides no information on free parameters, axioms, or invented entities; ledger is empty by necessity.

pith-pipeline@v0.9.0 · 5619 in / 954 out tokens · 24913 ms · 2026-05-24T12:57:55.658132+00:00 · methodology

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

32 extracted references · 32 canonical work pages · 2 internal anchors

  1. [1]

    F. S. Acton, Numerical Methods that Work, 1990, Mathematica l Association of America

    work page 1990

  2. [2]

    ApproxFun.jl, Julia package for function approximation, https://github.com/JuliaApproximation/ApproxFun.jl, https://juliaapproximation.github.io/ApproxFun.jl/latest/

    work page

  3. [3]

    U.M.Ascher , R.M.M.Mattheij, R.D.Russel, Numerical Solution of Bound ary Value Problems for Ordinary Differential Equations, 1987, SIAM

    work page 1987

  4. [4]

    https://doi.org/10.1145/3208103

    F.Br´ ehard, N.Brisebarre, M.Joldes, Validated and numerically effic ient Chebyshev spectral methods for linear ordinary differential equations, ACM Transactio ns on Mathematical Soft- ware, 44 (2018), 1-42. https://doi.org/10.1145/3208103

    work page doi:10.1145/3208103 2018

  5. [5]

    M. Barkatou, An algorithm to Compute the Exponential Part of a Formal Fundamental Matrix Solution of a Linear Differential System, Applicable Algebra in En gineering, Commu- nication, and Computing 8 (1997), 1–23

    work page 1997

  6. [6]

    chebfun for Matlab, https://chebfun.org

    work page

  7. [7]

    Multiple precision evaluation of the Airy Ai function with reduced cancellation

    S. Chevillard, M. Mezzarobba, Multiple-precision evaluation of the Airy Ai function with reduced cancellation, arxiv:1212.4731

    work page internal anchor Pith review Pith/arXiv arXiv

  8. [8]

    Holonomic Gradient Method-Based CDF Evaluation for the Largest Eigenvalue of a Complex Noncentral Wishart Matrix

    F.H.Danufane, K.Ohara, N.Takayama, C.Siriteanu, Holonomic Gradie nt Method-Based CDF Evaluation for the Largest Eigenvalue of a Complex Noncentral Wish art Matrix. https://arxiv.org/abs/1707.02564

    work page internal anchor Pith review Pith/arXiv arXiv

  9. [9]

    DeSolve (R-package), Solvers for Initial Value Problems of Differ ential Equations, http://desolve.r-forge.r-project.org/

    work page

  10. [10]

    Dieci, M

    L. Dieci, M. Osborne, R. Russell, A Riccati Transformation Meth od for Solving Linear BFPs I, SIAM Journal on Numerical Analysis 25 (1988), 1055–1073

    work page 1988

  11. [11]

    T.A.Driscoll, N.Hale, Rectangular spectral collocation, IMA Journ al of Numerical Analysis, 36 (2016), 108–132 https://doi.org/10.1093/imanum/dru062(2016)

    work page doi:10.1093/imanum/dru062(2016 2016

  12. [12]

    GSL, GNU scientific library, https://www.gnu.org/software/gsl/

    work page

  13. [13]

    Springer, 1993

    E.Hailer, S.P.Norsett, G.Wanner, Solving Ordinary Differential Equ ations I, II. Springer, 1993

    work page 1993

  14. [14]

    H.Hashiguchi, Y.Numata, N.Takayama, A.Takemura, Holonomic gra dient method for the distribution function of the largest root of a Wishart matrix, Journ al of Multivariate Analysis, 117, (2013) 296-312

    work page 2013

  15. [15]

    References for HGM, http://www.math.kobe-u.ac.jp/OpenXM/Math/hgm/ref-hg m.html 20

    work page

  16. [16]

    Hibi and et al, Gr¨ obner Bases: Statistics and Software Sys tems, 2013, Springer

    T. Hibi and et al, Gr¨ obner Bases: Statistics and Software Sys tems, 2013, Springer

    work page 2013

  17. [17]

    van Hoeij, Formal Solutions and Factorization of Differential Operators with Power Series Coefficients, Journal of Symbolic Computation 24 (1997), 1–30

    M. van Hoeij, Formal Solutions and Factorization of Differential Operators with Power Series Coefficients, Journal of Symbolic Computation 24 (1997), 1–30

    work page 1997

  18. [18]

    https://doi.org/10.1145/3476446.3535505

    M.Joldes, Validated Numerics: Algorithms and Practical Applicatio ns in Aerospace, Proceed- ings of the 2022 International Symposium on Symbolic and Algebraic C omputation, 2022, 1–2. https://doi.org/10.1145/3476446.3535505

    work page doi:10.1145/3476446.3535505 2022

  19. [19]

    DEtools in Maple, https://www.maplesoft.com/support/help/maple/view.aspx?path=DEtools

    work page

  20. [20]

    NDSolve in Mathematica, https://reference.wolfram.com/language/tutorial/NDSolveImplicitRungeKutta. html

    work page

  21. [21]

    The GNU MPFR library, https://www.mpfr.org/

    work page

  22. [22]

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

    work page 2006

  23. [23]

    Olver, A

    S. Olver, A. Townsend, A fast and well-conditioned spectral me thod, SIAM Review 55 (2013), 462–489

    work page 2013

  24. [24]

    pySDC project, a Python implementation of the spectral difer red correction, https://parallel-in-time.org/pySDC/

    work page

  25. [25]

    Risa/Asir, a Computer algebra system, http://www.openxm.org

    work page

  26. [26]

    scipy.linalg.solve https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html#scipy.linalg.solve

    work page

  27. [27]

    scipy.integrate.solve ivp, https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

    work page

  28. [28]

    scipy.optimize.least squares, https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html

    work page

  29. [29]

    cvxopt.solver.qp, https://cvxopt.org/userguide/coneprog.html#quadratic-programming

    work page

  30. [30]

    N.Takayama, L.Jiu, S.Kuriki, Y.Zhang, Computations of the Expec ted Euler Characteristic for the Largest Eigenvalue of a Real Wishart Matrix, Journal of Mu ltivariate Analysis 179 (2020), 104642

    work page 2020

  31. [31]

    L. N. Trefethen, Approximation Theory and Approximation Pra ctice, 2020, SIAM

    work page 2020

  32. [32]

    M.Winkel, R.Speck, D.Ruprecht, A high-order Boris integrator, J ournal of Computational Physics 295 (2015), 456–474. 21

    work page 2015