Linear-Time Probabilistic Solution of Boundary Value Problems
Authors: Nicholas Krämer
Philipp Hennig
Venue: NeurIPS 2021 Poster (2021)
Keywords: Differential equations
Gaussian processes
probabilistic numerics

Abstract: We propose a fast algorithm for the probabilistic solution of boundary value problems (BVPs), which are ordinary differential equations subject to boundary conditions.  In contrast to previous work, we introduce a Gauss-Markov prior and tailor it specifically to BVPs, which allows computing a posterior distribution over the solution in linear time, at a quality and cost comparable to that of well-established, non-probabilistic methods.  Our model further delivers uncertainty quantification, mesh refinement, and hyperparameter adaptation. We demonstrate how these practical considerations positively impact the efficiency of the scheme. Altogether, this results in a practically usable probabilistic BVP solver that is (in contrast to non-probabilistic algorithms) natively compatible with other parts of the statistical modelling tool-chain.



Linear-Time Probabilistic Solutions of Boundary Value Problems
Nicholas Krämer University of Tübingen Tübingen, Germany
Philipp Hennig
University of Tübingen and
Max Planck Institute for Intelligent Systems
ni nicholas. kraemer@uni- wwwerrrnnnnteeeeennnn.emm
Tübingen, Germany
Abstract
We propose a fast algorithm for the probabilistic solution of boundary value problems (BVPs), which are ordinary differential equations subject to boundary conditions. In contrast to previous work, we introduce Gauss-Mans -Markov prior and tailor it specifically to BVPs, which allows computing posterior distribution over the solution in linear time, a quality and cost comparable to that well- established, non-probabillitic methods. Our model further delivers uncertainty quantification, mesh refinement, and hyperparameter adaptation. We demonstrate how these practical considerations positively impact the efficiency of the scheme. Altogether, this results a practically usable probabilistic BVP solver that is (in contrast non- non-pro pan--oobbbilittc algorithms) natively compatible with other parts of the statistical modelling tool-chain.
Boundary value problems in computational pipelines
This work develops class of algorithms for solving ODE boundary value problems; that is, ordinary differential equations (ODEs)
útt) f()(t),t)
(1)
subject to left- and right-hand side boundary conditions Ly(to) = and Ry(tmax) = 3max The vector field f Rd > Rdd, as well as E Rdz xd R E Rdxxd to E R. tmax E RR yo E Rd and Ymax E RR are given. It is loss of generality consider first-order boundary value problem because higher-order problems can be transformed into first-order problems [1]-
Loosely speaking. solving BVPs amounts to following the law of dynamical system when "connect- ing two points". This setting is relevant to several scientific applications of machine learning. As motivation, consider three examples, all of which depicted in Figure 1. First, recovering the trajectory of pendulum between two positions amounts to solving the ODE Ti() = 99.81 sin(y(t)) subject to the positions boundary conditions. If the positions were interpolated without the ODE knowledge, the output would be physically meaningless. Second, BVPs arise when inferring the evolution of the case counts of people that fall victim to an infectious disease. A lack of counts of (a specific subset of) non-infected people at the initial time-point can be made up for by available counts of infected people the final time-point of the integration domain. Third, efficient manifold learning necessitates repeated computation of (geodesic) distances between two points, which amounts to solving BVPs [2, 3]. Depending on application details. the ability to produce structured output uncertainty or to enhance the algorithm by including additional sources of information can be crucial, Probabilistic numerical algorithms respond these challenges by solving problems of numerical simulation with probabilistic inference. For initial value problems, probabilistic solvers share linear- time complexity, adaptive step-size selection, and high polynomial convergence rates with their
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

Solution Derivativve
100
2nd Principal comp.
Derivative
Case counts
1.57
Solution
-1.57
0.0
1.5
Time 50
Time
1st Principal comp.
Figure 1: Recovering the trajectory of pendulum between two positions is BVP (LEFT). Lack of initial values can be made up by boundary values in an SEIR model (MIDDLE). Straight lines on manifolds give distance measures and demand solving BVP (RIGHT; depicted are the mean and ten samples of the probabilistic solution; principal components of 1000 MNIST images of the digit "I"). In all three figures, the ball/diamond markers express boundary conditions.
non-probabilistic counterparts [4-7], and further provide functionality to quantify uncertainty within probabilistic programs [8, 9].
Probabilistic BVP solvers have not yet reached this level quality. Existing probabilistic treatments of BVPs [10-12] iteratively condition Gaussian process on approximately "solving the BVP". Each such iteration requires solving generic least-squares problem of size equal to the number of employed grid points. The resulting cubic complexity puts severe upper limits grid resolution. Traditional, non-pro con-probabiliiticc BVP solvers (for instance, those presented in [I]) are very efficient but do not provide probabilistic output. Thereby, they would have to serve black-boxes inside probabilistic programs. In this work, we close this gap. The main idea of this paper that computing probabilistic solution of BVPs is fast if the prior is Markovian (Section 2). Probabilistic modelling provides additional advantages. In particular, algorithmic parameters can be estimated automatically (including those that must be provided by the user in traditional methods; Section 3-5).
2 Boundary value problems as probabilistic inference tasks
2.1 Generative model
Let a > 0. We define the process Y = [Yo..., Yol"= 140. xc) + Rd(++1) as the solution of the stochastic differential equation
dy dt) (t) dt BdW Y(to) N(mo. o²Co).
driven by a d-dimensional Wiener process W [40, 00) + Rd with diffusion l = 0²1 E IR dxd and initial parameters mo E Rd((-1), Co E (d((++))x((++1)) [4]. For the moment, we set = 1, m (O....,O) 0), and Co = I, and will discuss parameters calibration later. Let A E [d((++)xd((++)) and E I) 1)xd such that the zeroth component Yo of is the integrated Wiener process. The ęth component Y, of Y is the ęth derivative of the integrated Wiener process. In this setup, Yo models the BVP solution y. and Yg models the şth derivative of the BVP solution y: Yt(t) x y( = 0,.,, v [6. Equation 2]. This is the prior for the probabilistic BVP solver. Other choices are possible, too [5, Section 2.1].
For ODE solvers, the likelihood best described terms an information operator [13], For BVPs, there two sources of information: first, the boundary conditions
(x) LYo(to) yo and RR(Y) RYo (fmax) Vmnxx
and second, the differential equation, encoded by the information operator
((Y) Yi f(Y((),,)
(4)
Similar models are used in the gradient matching literature [14, 15]. Different to the likelihoods from conventional Bayesian inference, information operators used in probabilistie numerics are frree-free and often map between (possibly infinite-dimeniin vector spaces of functions [16, Section 2].

---

Let T := (fo...., ty = mmax be a grid on [tv. f For now, we assume this grid is given; Section 4 introduces a strategy for iterative mesh-refinement based on error-control. We will abbreviate (V) (l,) and fu.n. (fo...., ), n 0, o,... N. Using T, well as the likelihoods in Equations (3) and (4), the approximate ODE solution is captured by the posterior distribution
(( = 0, (Y) 0, L (Y) = 0).
(5)
Unfortunately, the full posterior (5) is intractable because of the non-linearity of f (which implies non-linearity in all (x) We will thus approximate with Gaussian: the probabilistic BVP solution.
2.2 Approximate Gaussian posterior inference
While the full posterior in Equation (5) cannot be computed in closed form, maximum-a-ppsterne (MAP) estimate is obtained by finding the minimum
arg (T) {- logp(Y L,((Y) = 0, L (Y = 0, (o. 0}.
(6)
This constrained optimisation problem can be solved with the iterated extended Kalman smoother (IEKS). The IEKS is state-space implementation of Gauss-Newton algorithm [17]. As such, one step of the IEKS computes the closed-form minimum of Equation (6) with a Kalman smoother, where the non-linear No.N replaced by its first-order Taylor approximation around the previous iterate. Under mild assumptions on the non-linearity of f and the magnitude of the objective at the optimum, Gauss-Newton methods are locally convergent with linear rate [18].
Each iteration of the IEKS returns a mean and covariance function. Eventually, the scheme approaches variant of the Laplace approximation of the posterior (note the shorthand of Equation (5))
(l) ~ N (TTMPP()) (MA)(t)) R ((Y (r, (oNNeR).
(this is non-standard Laplace approximation in so far as it employs Gauss -Nwton approximation of the Hessian). more detailed explanation is in Appendix A. The mean m (t) is the MAP estimate, because it minimises the objective in Equation (6). The covariance CMAP(() () is the inverse (approximate) Hessian of the negative log-posterior distribution, evaluated at MMAP(().
The Gaussian posterior returned by the IEKS is a probabilistic BVP solution. Thus, this basic version of the algorithm is already valid BVP solver, But some degrees of freedom remain, whose efficient selection improves performance significantly. These will be the concern of the remainder of this work. Table presents an outline.
Table I: Configuration of the remaining degrees of freedom.
What?
How?
Where?
Initialisation of the IEKS
ODE filter with Gaussian bridge
Section 3
Mesh T
Error control
Section 4
Diffusion a
Quasi-maximum likelihood estimation
Section 5
Initial parameters mg. Co
Expectation-maxxiisatio
Section 5
3 An initial guess is not strictly necessary
Like every optimisation algorithm, the IEKS needs appropriate initialisation. Not only does the number of iterations depend on the proximity of the initial guess the optimum, but BVPs often allow multiple solutions, and the algorithm can find only one of them [19, p. 10]. probaabiiiiii solvers outsource this issue to the user by expecting that an initial guess is provided. While the same strategy is available for the probabilistic solver, there are natural alternatives in non non-iterative Gaussian smoothers (Section 3. I), which further benefit from combination with a bridge prior (Section 3.2).
For example, at the time of this writing, the BVP solvers in SciPy, Matlab, and DifferentialEquation.] require the user pass vector of initial guesses of the solution initial grid to the algorithm [20-22].
3

---

EN
(1o)
(tm)
((xx)
f
ER
Figure 2: Construct a bridge by considering boundary values first. Graphical depiction of the inference problem (LEFT). Samples from the Mouss-Markov (CENTRE) and bridge prior (RIGHT).
3.1 Initialisation with an extended Kalman smoother
If the target of Laplace approximation of the BVP posterior is relaxed only some Gaussian approximation, an initial guess can be computed with an extended Kalman smoother (EKS) [13, 5]. Like the IEKS, the EKS linearises the non-linear ODE measurements ta.N with first-order Taylor series. It differs from the IEKS in the position around which the approximation is constructed. The IEKS linearises all f..v at once after each completed forwards-backwant each b, on the fly during the forward pass, at the respective predictive mean [23]. In other words, the EKS does not need an initial guess, which why the tool of choice to construct one [5].
pass. The EKS linearises
the BVP is linear, the EKS computes the true posterior [23, 24]. If the BVP is non-linear, the EKS introduces significant linearisation error wherever the predictive distribution deviates strongly from the true posterior. Unfortunately, in its standard implementation, the EKS necessarily starts with incomplete information about the state (fo) and higher-order derivatives (initialisation of which is crucial to probabilistic initial value problem solvers as well [7]. Ensuring that the prior distribution satisfies the boundary conditions by construction solves this problem because the iteration can never drift too far away from the optimum. The following Section 3.2 explains more.
3.2 Changing the order of updates to build a bridge prior
Recall that there are three sources of information: the left-hand side boundary condition L,, the right-hand side boundary condition ln, and the ODE measurements fg:N- If the initial and terminal state of the prior distribution are forced to accommodate é, and Ê before conditioning on fo.N. samples from the resulting Gaussian bridge satisfy the boundary conditions by construction; see Figure 2. The linear-time complexity of Gaussian filtering/smoothingg is preserved through this change the order updates, because the Markov property of Y yields
N-1
(R) = (1) Cr, (R)
II p(Y (X))))) n=
(8)
The transition densities p( (t,),, kR) well as the initialisation |(,, are available in closed form (Appendix B). A reader familiar with the prediction-correctio nature of Gaussian filtering can think of the implementation as follows: Roughly speaking, each prediction step of the EKS with bridge prior involves extrapolating from the current state to the terminal state, conditioning on the boundary condition ln. and smoothing back to the current state. Therefore, the computational complexity of an EKS borward-backwer pass with the bridge prior is about twice as large compared to an EKS are in Appendix B.
pass with the conventional prior. Precise derivations
Figure 3 shows that this extra cost is made up for by the improved linearisation behaviour because encoding the boundary conditions into the prior improves the initialisation drastically. Following the neowards-backward
pass with the EKS, the IEKS requires only few more iterations to find fixed point similar to the truth. Not using either the bridge prior or the EKS results an initial guess that takes more iterations to find a fixed point of a lower approximation quality. Abandoning both options, which aligns with initialisation of traditional BVP solvers, is least efficient since converges to an inaccurate fixed-point.

---

Conventional & Guess
Bridge Guess
Conventional & EKS
Bridge EKS
Solution
1.5
0.5
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5 Time
1.0
0.0
0.5
1.0
Time
Time
Time
combination, EKS and bridge prior initialise well. Initialisation and five iterations the IEKS depicted from light to dark the 20th problem in [25] (truth in black). Without bridges, and with an initial guess of constant twos, the fixed-point of the Gauss-Newtor scheme is inaccurate on N = points (LEFT). Using either a bridge prior (CENTRE LEFT) the EKS (CENTRE RIGHT) lessens this problem. The bridge/EKS combination finds an accurate estimate almost immediately (RIGHT), because due the bridge prior, the EKS linearises around more accurate location than it would with a conventional prior during the first forward-pass.
Figure 3:
Linear BVPs undo the effect of the bridge prior because the full posterior is computed accurately with conventional Kalman smoother [23]. Likewise, an IEKS iteration linearises all lu:N at once, outside of the barudd-ba pass, which renders the bridge obsolete as well. In other words, the changed order of updates only relevant for the initialisation.
Of course, a fixed point of the IEKS not necessarily a reliable BVP solution: its accuracy depends on the number and distribution of mesh points. The following Section 4 develops a principled and probabilistic approach to error control and mesh refinement in the BVP solver.
Estimate the error and refine the mesh
So far, the mesh T was assumed as given. The larger the size of this mesh is, the more accurate the solution becomes; but computational cost grows linearly with the mesh size. Low error tolerances thus require smart meshing via error control. There estimators, all of which connect to the probabilistic formulation of solving BVPs.
(plus one) natural candidates for error
Standard deviation: The output of the IEKS is Gaussian process, which be evaluated at any point in the domain of the boundary value problem [24, Chapter 10]. Its associated standard deviation the alternatives explained below is that it comes
provides error estimator. The advantage (essentially) for free as part of the dense output of the posterior. A potential downside of this intrinsic estimator is its dependence the calibration of hyperparameter (more this in Section 5).
The inference problem (Equation (5) and Figure 2) is constructed by conditioning the prior Y attaining consistently small values in its residual ((Y) = Yt(t) f(Y(t), t). Recall that if Yo were the true ODE solution, and Y, were its derivative, {(Y) would be zero on the whole domain. Thus, the residual of the posterior mean of the approximate ODE solution estimates the error, which is a common approach in traditional, probabiliisicc algorithms as well (for instance [19] or [1, Section 9.5. On side note, considering the full posterior distribution implies that the residual would be a deterministic transformation of a random variable. Thus in principle a random variable might make
Residual:
appropriate model for the residual error than point estimate (see Remark 1). However, this quantity will reveal itself as inaccurate in the benchmarks below.
Remark 1. For a Gaussian process posterior Y, the law of E(Y') intractable general. Lineari- sation of (at the previous iterate, like in the IEKS) unlocks Gaussian approximation: denote the Gaussian random variable (t))= 2s (tY) An upper bound of the probability of |Z]| exceeding some tolerance,
p (1) 1ot2) <
[[(Z)(0]2) fIof,
is due to the Markov inequality and a third approach to error control. The numerator of the right-hand side will be treated as an error estimator in the benchmarks below. The main difference to the point estimate is that the probabilistic version punishes magnitude and uncertainty in the residual.
5

---

All three options (which we denote by generic e from now on) estimate the error at given t. For mesh refinement, however, it more instructive to consider the accumulated error on each interval
On -
An lee(t)] dt
1/2
n= 0,.,... N
(10)
If each 6n is sufficiently small, the BVP solution is adequately accurate and the mesh appropriately fine. On those intervals where ' is too large, we introduce new grid points as follows. Assuming that the integrated error is of order p > 0, En Ohh)), splitting the interval into two equally large parts reduces the error by a factor 2 and splitting into three equal parts by a factor 3 We use these threshold values to guide where to introduce one point and where to introduce two points. Like Kierzenka and Shampine [26]. we never introduce more than two at once. For the experiments herein, and y-tiimes integrated Wiener processes, we use p V + 1/2 (which has not been proved yet but seems like a reasonable conjecture in light of Theorem of Tronarp et al. [5] and our experiments).
The integral that underlies en can usually not be computed in closed form but needs to be approximated by a numerical integration scheme. We use Bayesian quadrature (BQ) [27]. Not only does it fit neatly into the probabilistic framework, but it also allows us to place quadrature nodes freely in each domain In, 1n+1). If viewed as an integral from to I, we choose quadrature nodes at O, 0.33, 0.5, 0.67, 1. These locations include the boundary points of the domains (0 and 1), as well the nodes that will be introduced in case the error is too large (either 0.5, or 0.33 and 0.67). This has the advantage that, reuse the evaluation of the posterior at the new mesh points. If the
at the start of the next iteration, residual estimates the error, there is another advantage. Since the IEKS approximation is a minimum of constrained optimisation task, the residual is zero at the boundaries of each interval. In this case, the integral is only computed on the three interior nodes. For the same reasons, non-probabilistic solvers with residual control usually employ Gauss-Lobatto schemes [19].
final motivation for BQ that we can tailor an integration kernel to e. For instance the following reproducing kernel Hilbert spaces (RKHSs) are known [5]: (i) the RKHS of y-times integrated Wiener process priors Y (.) is the Sobolev space of (v + 1)-times weakly differentiable functions; (ii) under some regularity assumptions on the ODE vector field, as well as the (assumed to be) unique solution of the ODE, the RKHS of the residual (-)) is the Sobolev space of y-times weakly differentiable functions. Therefore, we base the BQ scheme on (v - (v/1/2 order Matérn prior, which has the same native space as the residual [28] (we use an exponentiated quadratic kernel for 3 because the required kernel embeddings are casier to compute [29, Appendix JJ). The accuracy of the quadrature approximation improves with increasing v [29]. Matching the BQ kernel to the ODE prior parallels the choices of quadrature schemes in pro-probaiillitic solvers [19], and ensures that the accuracy of the numerical integration does not limit the validity of the error estimate.
Which one is the most reliable error estimate? As first testbed, we use the seventh in a collection of test problems for BVP solvers by Mazzia [25] (which will feature heavily in the remainder of this work). The derivative of the solution of this linear BVP approaches singularity if a specific parameter is chosen sufficiently small (we use 10 3). This poses challenges for error estimators and mesh-refinement strategies. The error estimates are visualised in Figure 4. They suggest that at high tolerances, the standard deviation is more accurate than the residual; at low tolerances, the situation is reversed. This trend is preserved when moving to more challenging setups (see Section 6).
To conclude, the probabilistic framework introduces three options for error estimation and comes with a natural algorithm to compute accumulated errors in BQ. With everything explained so far, we can solve BVPs with an algorithm that adaptively refines the mesh when the solution is not sufficiently accurate. After each mesh refinement, the iteration is restarted. While it may be clear that the initial guess for the new IEKS implementation should be the approximate posterior from the previous computation, beginning new Gauss-Newtr scheme offers the chance to update the choice of other hyperparameters and thus set up a more appropriate probabilistic model for free (Section 5).
5 Calibration of hyperparameters with maximum-likelihood and EM
Thus far, an approximate BVP solution has been computed with a, mgg and Co set to default values. Maximum- -ikimum-likelhoodd estimates of these hyperparameters can be computed by coordinate ascent,

---

N=
N 25
125
N 625
(Prob.) Residual Std. deviation
10.6
10a
10-5
www
10-14
Time
Time
Time
Time
Figure 4: Error estimation on the seventh testproblem in [25]. Evaluated N 5 (LEFT), 25 (CENTRE LEFT), N 125 (CENTRE RIGHT), and N 625 equidistant grid points (RIGHT). Standard deviation (TOP ROW) and residual (BOTTOM ROW) respectively the probabilistic residual (BOTTOM ROW). True error in black. A good estimate accurately measures the magnitude of the error well as the location of large deviation. On few points, the latter less important so the well-calibrated standard deviation provides good estimate. On many points, is underconfident. Since for large N, the location of the error becomes increasingly important, the residual should be used; the probabilistic residual consistently underconfident. The "winners" of each column have darker colour.
which repeats alternating updates
new == arg max, plg, foNN, l0, my new Cg
(1la)
m0 Crew arg max C, log p(L,, Lo.N LR new mo. Co).
(11b)
until some stopping criterion is satisfied [30]. A quasi-maximum likelihood update for new (Equation 1 la) is available closed form as a by-product of the forward-pass of each IEKS iteration. This is also true for the specific order of updates detailed previously in Section 3.2 (Proposition 2 below).
Proposition Assume that the initial covariance and the diffusion of the Wiener process depend multiplicatively the scalar o² (recall Equation (2)). L. L. and lu:N noise-free (which herein they always are), the covariance of the posterior process depends multiplicatively on o² and a quasi-maxinge maxismmmximmn likelihood estimate for is available in closed form.
The proof of this proposition similar to the proof of Proposition 4 of Tronarp et al. [13] yet requires few additional manipulations because of the boundary value information contained in the bridge. A derivation - and the precise formula for the quasi-MLE are in Appendix C.
While is tuned with quasi-maximum likelihood estimationn the parameters and Co of the initial distribution are separately calibrated with single step of the popectation-maximistt algorithm [31, 32] whenever the mesh needs to be refined, which implies restart of the IEKS. In other words, this "outer loop" around calls to the IEKS is already part of the computational budget; therefore, sensible updates to the initial distribution parameters mo and Co are free. The general idea of EM is to maximise lower bound of Equation (b instead of maximising directly, by computing alternating E- and M -tstps. For parameter estimates in state-space models, the E-step of the EM algorithm is the posterior distribution in Equation (5) (see e.g. [33]), a Gaussian approximation of which is available through the IEKS: recall (tAr(t) N (mmssr (MAr(()). The M Msste consists of [23, Theorem 12.5 and Algorithm 12.7]
(EM)
mew mMAr((t)
(12a)
o o²Cmsx((oo) (mp (ml))) m T
(12b)
EM steps always increase the likelihood, and for exponential families, convergence to a stationary point of the likelihood function is guaranteed [34, 32]. Thus, computing alternating and M -steps until convergence (which we do not do) would eventually yield good estimate of the parameters. But already the pre-asymptotic regime and for fixed total number IEKS iterations, making an EM update every few steps helps convergence of the IEKS in subsequent iterations (Figure 5).
7

---

No EM Updates
EM Updates every 5th Iteration
Updates every iteration
Solution
10a
A
V
8
8
Residual
10
A
10:14
0.0
0.5
1.0 0.0
0.5 Time
1.0
0.0
0.5 Time
1.0
Time
Figure 5: EM helps the IEKS overcome unknown initial conditions. Depicted are fixed total of the
first 25 IEKS iterations (light to dark in each respective colour) N 6 grid points, initialised with an EKS using a 7-times integrated Wiener process bridge prior on the 20th test problem [25]. Without any EM updates to the initial condition, the convergence of the IEKS is inhibited (LEFT). EM updates every fifth IEKS iteration lead to the residual converging to zero reliably (CENTRE). Too many EM updates
optimal either (RIGHT).
10
10
statistic
RMSE
10-6
RMSE
10-4
Time
10
10-1
10
10-9
10-2
10
10
10
SdPy
10-3.
10a
10a
10-2 10-1 10
10a
10º
10a
10*
10a
No. grid points
Time (s)
No. of grid points
No. of grid points
Figure 6: Results Bratu's problem. The higher-order solvers converge least fast the SeiPy reference (LEFT) and the wall-clock time grows linearly with the number of grid points (CENTRE RIGHT; linear complexity reference line in the background). Our research code slower by factor 100 than this highly optimised code base, which is a success compared to prior implementations of probabilistic BVP solvers (CENTRE LEFT). The -2-statistic remains within 95% confidence (RIGHT; intervals shaded in gray, mean (= 1) in black). To show mesh refinement, the initial grid consisted of only three points; the probabilistic solver initialises with EKS and bridge, and uses the standard deviation as an етгог estimate.
6
The solver converges quickly on test problems
place, evaluate the performance of the solver on a range of scenarios. All experiments use the CPU of a consumer-level laptop. An efficient probabilistic numerical method should provide both good point estimate (through its posterior mean) and error estimate (through its posterior covariance). First, the approximation error should decrease rapidly with the number of grid points; we report root-mean poot-mean-squar square errors the lower, the better. Second, the width of the posterior distribution should be representative of the numerical approximation error (which has, to some extent, been shown Section 4 already); posterior uncertainty is calibrated. A simulation of Bratu's 's problem [36] for varying tolerances and orders v suggests that the solver performs well in both metrics (Figure 6). Reassuringly, higher orders of the solver lead to faster convergence, which motivates the analysis of convergence rates akin to analysis of Tronarp et al. [5] for initial value problem solvers. The experiments also suggest that the uncertainties are calibrated but tend to be under-confident. This phenomenon is known from probabilistic initial value problem solvers (ef. 6]. been studied for general Gaussian process approximations by [37].
Now that all parts
the =-sstatistic [35]. it is close to the

---

10a
Tol 10
0.21
10-2
0.30.11
(S, 0.0011
(5,0.2)
RMSE
10-4
CLO, 0.006]
66, 0.31
16, 0.2)
(5, 0.002)
(5. 0.004)
Tol 10-6
(77, 2.01
10
1255 0.5)
CL 30,
112. 0.30
(31. 0.71
10-
1178, 0.01)
169, 0.03)
(S3, 0.008)
(292. 0.021
(11P, 0.011
10-10
7th Problem
23rd Problem
24th Problem
28th Problem
32nd Probiem
Figure 7: The solver efficiently computes (mostly) calibrated posteriors on many problems. Proba- bilistic solver (V; 6) versus SciPy's 's BVP solver (0). Markers are annotated with the number of grid points and runtime (in seconds). The tolerances are 10 (GRAY) and 10 6 (BLUB). The closer a coloured marker is to its reference line, the better. The fewer grid points and the less time required the better. Fill-color describes calibration: x² is within 80 % (GREEN), within 99% (oRANGE), or outside of these ranges (RED). SciPy does not allow notion of calibration.
Efficient mesh refinement and fast convergence are evident when considering wider range of test problems. Figure 7 depicts the results of simulating five BVPs (all from Mazzia [25]): the 7th problem approaches singularity in its derivative, the 23rd problem has a boundary layer at mmxx the 24th problem describes a fluid mechanical model of a shock wave, the 28th problem has a corner layer at min- and the 32nd problem involves fourth-order derivatives. On all problems, the probabílistic solver efficiently computes calibrated posteriors at specified tolerances.
7 Related work
How does the proposed algorithm fit into the context of state-of-the-ar probabilistic and non- probabilistic BVP solvers? Headway the probabilistic solution of BVPs has been made by Hennig and Hauberg [10], Arvanitidis et al. [1 I]. and John et al. [12]. Hennig and Hauberg [10] and Arvani- tidis et al. [11] focus the application of BVP solvers to Riemannian statistics. None of the three algorithms exploit the state space structure of the prior with its beneficial computational complexity, nor they concerned with error estimation, mesh refinement, and the other computational aspects to the extent that this work is. Other, conceptually different probabilistic algorithms for BVPs have been proposed by Skilling [38], Chkrebtii et al. [39], Conrad et al. [40], O'Leary and Harker [41].
In terms of accuracy and cost, the present approach should rather be compared to off-the-shell non- probabilistic BVP solvers: for instance, those im- plemented in Matlab [42, 19, 26], Python/SciPy [43], and Julia [44]. These toolboxes contain algorithms that implement collocation formulas and gain linear- time complexity from sparse system matrices. The Markov property makes our algorithm equally fast (in terms of the number of grid points N) (Table 2). The computational complexity of Algorithm 1 is (T JIEKS Nu3dd), where JIEKS is the number of IEKS iterations, and Mush is the number of mesh re- finements. In our experiments, we found IIEKS to be small, usually bounded by 10. The mesh refinement is designed to make Mesh as small as possible. Linear
Algorithm 1: BVP Solver
Input: BVP, mesh, order (v)) tolerances. Output: Probabilistic BVP Solution Initialise with bridge and ODE filter; while 3 Z interval with large error do
Run IEKS;
Update mi and Co (Equation (12));
Update o:
Compute error between gridpoints;
Refine mesh where necessary;
end
complexity in N stems from the state space implementation of the IEKS and could potentially be reduced to log N by temporal parallelisation [45]. The cubic complexity in and in d stems from the matrix-min matrix operations that are required in Kalman filter step [23]. Cubic complexity suggests that high-order BVPs should be solved directly, without transforming them into first-order [46]. This is not uncommon for BVP solvers [1, Section 5.6] and used herein (a version of e, that is suitable to high order ODEs is explained in Appendix D).

---

Table 2: Comparison of probabilistic and non-probabilistic BVP solvers.
Non-probabilisticc
Probabilistic (present work)
(NN) achieved by
Sparse matrices
Markov property
Error estimates
Residual (point estimate)
Many options, e.g. standard deviation
Initial guess
Mandatory
Optional
Uncertainty quantification
No
Yes
8
Conclusion
We have arguably provided the first practically usable probabilistic BVP solver. Our method achieves the same linear computational complexity as off-the--helff solvers, with high-quality point estimates and calibrated uncertainty. Algorithmic parameters can be set automatically by the method, including some that have to be set manually for non-probabilistic solvers. Our method thus closes methodological gap in the toolbox of probabilistic numerics.
Acknowledgements
The authors gratefully acknowledge financial support by the German Federal Ministry of Educa- tion and Research (BMBF) through Project ADIMEM (FKZ O1IS 8052B). They also gratefully acknowledge financial support by the European Research Council through ERC SIG Action 757275 / PANAMA; the DFG Cluster of Excellence "Machine Learning New Perspectives for Science", EXC 2064/1, project number 390727645; the German Federal Ministry of Education and Research (BMBF) through the Tübingen AI Center (FKZ: O1IS 18039A); and funds from the Ministry of Science, Re- search and Arts of the State of Baden-Württemben Württemberg. Moreover, the authors thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for supporting Nicholas Krämer.
The authors thank Nathanael Bosch, Filip Tronarp, and Jonathan Schmidt for valuable discussions. They are grateful to Georgios Arvanitidis for helping with the implementation of the manifold example in Figure 1.
References
[1] Uri M Ascher, Robert MM Mattheij, and Robert D Russell. Numerical Solution of Boundary Value Problems for Ordinary Differential Equations. SIAM, 1995.
[2] Manfredo Perdigao do Carmo. Riemannian Geometry. Birkhäuser, 1992.
[3] Georgios Arvanitidis, Lars Kai Hansen, and Søren Hauberg. Latent space oddity: On the curvature of deep generative models. In International Conference on Learning Representations, 2018.
[4] Hans Kersting, Tim Sullivan, and Philipp Hennig. Convergence rates of Gaussian ODE filters. Statistics and Computing, 3096)119 2020.
[5] Filip Tronarp, Simo Särkkä, and Philipp Hennig. Bayesian ODE solvers: The maximum posteriori estimate. Statisties and Computing, 31(3)), 2021.
[6] Nathanael Bosch, Philipp Hennig, and Filip Tronarp. Calibrated adaptive probabilistic ODE solvers. In International Conference on Artificial Intelligence and Statistics, 2021.
[7] Nicholas Krämer and Philipp Hennig. Stable implementation of probabilistic ODE solvers. arXiv:2012.10106, 2020.
[8] Hans Kersting, Nicholas Krämer, Martin Schiegg, Christian Daniel, Michael Tiemann, and Philipp Hennig. Differentiable likelihoods for fast inversion of 'likelihood-fres dynamical systems. In International Conference on Machine Learning. 2020.
[9] Jonathan Schmidt, Nicholas Krämer, and Philipp Hennig. A probabilistic state space model for joint inference from differential equations and data. arXiv:2103.10153 2021.
10

---

[10] Philipp Hennig and Søren Hauberg. Probabilistic solutions to differential equations and their application Riemannian statistics. In International Conference on Artificial Intelligence and Statistics, 2014.
[11] Georgios Arvanitidis, Soren Hauberg, Philipp Hennig, and Michael Schober. Fast and robust shortest paths on manifolds learned from data. In International Conference on Artificial Intelligence and Statistics, 2019.
[12] David John, Vincent Heuveline, and Michael Schober. GOODE: A Gaussian off-the-shelf ordinary differential equation solver. In International Conference on Machine Learning, 2019.
[13] Filip Tronarp, Hans Kersting, Simo Särkkä, and Philipp Hennig. Probabilistic solutions to ordinary differential equations as non-linear Bayesian filtering: new perspective. Statistics and Computing, 296)): 2019.
[14] Ben Calderhead, Mark Girolami, and Neil D Lawrence. Accelcrating Bayesian inference over nonlinear differential equations with Gaussian processes. In Advances in Neural Information Processing Systems, 2009.
[15] Philippe Wenk, Gabriele Abbati, Michael A Osborne, Bernhard Schölkopf, Andreas Krause, and Stefan Bauer. ODIN: ODE-informed regression for parameter and state inference in time-continuous dynamical systems. In AAAI Conference on Artificial Intelligence, 2020.
[16] Jon Cockayne, Chris J Oates, TJ Sullivan, and Mark Girolami. Bayesian probabilistic numerical methods. SIAM Review, ((()))56-789, 2019.
[17] Bradley M Bell. The iterated Kalman smoother as a Gauss-Newton method. SIAM Journal on Optimization, 4(3):226-636, 1994.
[18] Knoth. A globalization scheme for the generalized Gauss-Newton method. Numerische Mathematik, 56(6):591-607, 1989.
[19] Jacek Kierzenka and Lawrence F Shampine. A BVP solver based on residual control and the Matlab PSE. ACM Transactions on Mathematical Software, 27(3):299-316, 2001.
[20] SciPy's solve_bvp documentation. http://///cooccccccaoooooo/oo//s wttp/////scciip//eeeeeeeeeee cenerated/scipy
www.graaee.slvve........... html. Accessed: May 23, 2021.
[21] Matlab's bvpinit documentation. http:://ikk.aattooorrsssoooo/oop bypinit. html. Accessed: May 23, 2021.
[22] Julia's BoundaryValueerroblem documentation (in (iniffereenttillaaatatiooo addiffee..ciml..ate wttttaaleae/taaleeeeeeeeaaaaalll wttpt/taaaeeeeeeeeeeeeeeeeeeeeee http:ttyeeeeppppeeeeeeeeeeee wttp:tpetttttttttppepppppeeeeeee Accessed: May 23, 2021.
https:
[23] Simo Särkkä. Bayesian Filtering and Smoothing. Cambridge University Press, 2013.
[24] Simo Särkkä and Arno Solin. Applied Stochastic Differential Equations. Cambridge University Press, 2019,
[25] F Mazzia. Test set for boundary value problem solvers, release 0.5, 2014.
[26] Jacek Kierzenka and Lawrence F Shampine. A BVP solver that controls residual and error. Journal of Numerical Analysis, Industrial and Applied Mathematics, 3(1-2):27-41, 2008
(27] François-Xavier Briol, Chris J Oates, Mark Girolami, Michael A Osborne, and Dino Sejdinovic. Probabilistic integration: role in statistical computation? Statistical Science, 34(1):1-22, 2019.
[28] Holger Wendland. Scattered Data Approximation. Cambridge University Press, 2004.
[29] François-Xavier Briol, Chris J Oates, Mark Girolami, Michael A Osborne, and Dino Sejdinovic. Probabilistic integration. 152200/5 2015.
[30] Stephen Wright. Coordinate descent algorithms. Mathematical Programming, 151(1):3-34,
2015.

---

[31] Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 1-(1):1-22, 1977.
[32] Robert H Shumway and David S Stoffer. An approach to time series smoothing and forecasting using the EM algorithm. Journal of Time Series Analysis, 3(4):253- 264, 1982.
[33] Radford M Neal and Geoffrey Hinton. A view of the EM algorithm that justifies incremental, sparse, and other variants. In Learning in Graphical Models, pages 355-368. Springer, 1998.
[34] CF Jeff Wu. On the convergence properties of the EM algorithm. The Annals of Statistics, pages 95-103, 1983.
[35] Yaakov Bar-Shalom, X Rong Li, and Thiagalingam Kirubarajan. Estimation With Applications to Tracking and Navigation: Theory Algorithms and Software. John Wiley Sons, 2004.
[36] Gh Bratu, Sur les équations intégrales non linéaires. Bulletin de la Société Mathématique de France, 42:113-142. 142, 1914.
[37] Toni Karvonen, George Wynne, Filip Tronarp, Chris Oates, and Simo Särkkä. Maximum likelihood estimation and uncertainty quantification for gaussian process approximation of deterministic functions. SIAM/ASA Journal on Uncertainty Quantification, (3)-992-958, 2020.
[38] John Skilling. Bayesian solution of ordinary differential equations. In Maximun entropy and Bayesian methods, pages 23-37. Springer, 1992.
[39] Oksana A Chkrebtii, David A Campbell, Ben Calderhead, and Mark A Girolami. Bayesian solution uncertainty quantification for differential equations. Bayesian Analysis, 1239) 1267, 2016.
[40] Patrick R Conrad, Mark Girolami, Simo Särkkä, Andrew Stuart, and Konstantinos Zygalakis. Statistical analysis of differential equations: introducing probability measures on numerical solutions. Statistics and Computing, 27(4):1065-1000 1082, 2017.
[41] Paul O'Leary and Matthew Harker. Inverse boundary value problems with uncertain boundary values and their relevance to inclinometer measurements. In 2014 IEEE International Instru- mentation and Measurement Technology Conference (12MTC) Proceedings, pages 165-169.
IEEE, 2014.
[42] Lawrence F Shampine and Mark W Reichelt. The Matlab ODE suite. SIAM Journal on Scientific Computing, 18(1):1-22, 1997.
[43] Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy, David Courna- peau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, et al. SciPy 1.0: Fundamental algorithms for scientific computing in Python. Nature Methods, 1763):261-222, 2020.
[44] Christopher Rackauckas and Qing Nie. DifferentialEquaatinns performant and feature-rich ecosystem for solving differential equations in Julia. Journal of Open Research Software, 5(1), 2017.
[45] Fatemeh Yaghoobi, Adrien Corenflos, Sakira Hassan, and Simo Särkkä. Parallel iterated extended and sigma-point Kalman smoothers. In IEEE International Conference on Acoustics, Speech and Signal Processing, 2021.
[46] Nathanael Bosch, Filip Tronarp, and Philipp Hennig. Pick-and-mix information operators for probabilistic ODE solvers. arXiv:2110.10770, 2021.
12


Reviews:
[{"review_id": "qyKxpMLqwFD", "reviewer": "Reviewer_9rs3", "summary": "The manuscript entitled “Linear-Time Probabilistic Solution of Boundary Value Problems” introduces a Gauss–Markov prior and tailors it specifically to BVPs, which allows computing a posterior distribution over the solution in linear time, at a quality and cost comparable to that of well established, non-probabilistic methods. Existing probabilistic solvers has cubic complexity which puts severe upper limits on grid resolution. Traditional, non-probabilistic BVP solvers are very efficient but do not provide probabilistic output. Thereby, they would have to serve as black-boxes inside probabilistic programs. \nThe main contribution of this work is that this proposed solver can deliver uncertainty quantification, mesh refinement, and hyperparameter adaptation with linear time (although it is still much slower than the non-probabilistic BVP solvers). This results in a usable probabilistic BVP solver that is natively compatible with other parts of the statistical modelling tool-chain.\n\n", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "In general, this paper is in a good shape. I have some comments on this paper as follows:\n\n1.\tFigure 5 shows too many EM updates are not optimal either (RIGHT). The authors should give an explanation for this phenomenon.\n\n2.\tIn Figure 6, the legend sits on Subplot 3. However, this legend should belong to all these 4 subplots. Also, the legend shouldn’t cover the data points on the plots. \n\n3.\tAlthough the proposed solver has linear-time complexity in terms of the number of grids, it is about 100 times slower than the off-the-shelf probabilistic solver which also has linear-time complexity. Authors should comment on what causes such a slowness compared with the non-probabilistic solver.\n\n4.\tI don’t see the appendix for this paper although the paper says it has.\n\n\nThis paper tries to close the gap between the probabilistic and non-probabilistic BVP solvers such that the new proposed probabilistic BVP solver has a linear-time complexity which is much better than the exiting probabilistic BVP solvers, and it can deliver the algorithmic parameters automatically, including some that have to be set manually for non-probabilistic solvers. It beats non-probabilistic BVP solvers in a way that this solver is naturally compatible with other statistical modelling tool-chain. Although the experimental results show this solver is still much slower and gets poorer point estimates than the off-the-shelf non-probabilistic solvers. But it made some improvements for the probabilistic BVP solvers. This paper has its significance.\n\nIn summary, I recommend its publication in NeurIPS after those issues are addressed.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The author should address more of this solver’s limitations and discuss the potential improvement in the future. Limitations can be the Markov assumption, poorer point estimation than non-probabilistic solvers, the slowness of the computation compared with the off-the-shelf non-probabilistic solvers, ", "review_text": "", "extra_scores": {"ethical_concerns": "No", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5 hours"}, "extra_text": {}}, {"review_id": "cuYdOIJjPK5", "reviewer": "Reviewer_j7RB", "summary": "In this paper, a fast algorithm for the probabilistic solution of BVPs are proposed. Comparing with the classical deterministic methods, the proposed method allows computing a posterior distribution in linear time while further delivers uncertainty quantification, mesh refinement and hyperparameter adaptation.", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is generally written well. The mathematical derivations seem proper and the discussions on initial guess, mesh refinement and hyperparameter calibration are helpful.\n\nMy main concern about the paper is its lack of comparison with baseline methods. Bayesian approach for BVPs was discussed in many classical uncertainty quantification literatures, for example:\nChkrebtii, O. A., Campbell, D. A., Calderhead, B., & Girolami, M. A. (2016). Bayesian solution uncertainty quantification for differential equations. Bayesian Analysis, 11(4), 1239-1267.\nO'Leary, P., & Harker, M. (2014, May). Inverse boundary value problems with uncertain boundary values and their relevance to inclinometer measurements. In 2014 IEEE International Instrumentation and Measurement Technology Conference (I2MTC) Proceedings (pp. 165-169). IEEE.\nHowever, the authors didn't provide a thorough review on these previous works in the area nor situate their approach with respect to them.\n\nAnother possible pitfall is: BVP problems can be ill-posed. How does the proposed method perform on those problems?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "See main review.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "MQuEtZbOVS", "reviewer": "Reviewer_zNRW", "summary": "The paper discusses refinements to current probabilistic BVP solvers, namely the IEKS. The paper proposes new initialisation and meshing strategies for speed-up. This is an interesting paper with potentially high impact.\n", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Post-discussion update. The authors have adressed my concerns, and this work does improve over earlier work substantially, and there is sufficient discussion about this. Some issues still stand, but I'm nevertheless raising my score to 6, and would be happy to see this paper accepted. \n\n--------\n\nThe paper discusses refinements to current probabilistic BVP solvers, namely the IEKS. The paper proposes new initialisation and meshing strategies to speed-up the method. This is an interesting paper with potentially high impact.\n\nI’m having lots of trouble understanding the papers notation and nomenclature. The paper talks of posteriors without introducing any priors or likelihoods. The math is incomplete (see below). This makes following the paper unnecessarily difficult. In NIPS papers should be accessible to wide audience, and not limited by subdomain shortcuts. The sec 3 would be easier to follow if one would mathematically define at least the gist of bridge and EKS, currently the method is only in the appendix.\n\nThe paper is otherwise well written. All the methods are sensible and the individual improvements are demonstrated also empirically.\n\nThere is strangely no empirical comparisons between the final tuned IEKS and the vanilla IEKS. I don’t understand why this is missing, since the main contribution of the paper is to improve upon IEKS (or is it..). Maybe the authors feel that this would be too trivial experiment. However, the improvements gained from the individual “tricks” have varied effects, and its not clear which one improves what part of the performance. In fig5 one already gets very low residuals with vanilla method, and in fig 6 the reference method seems way faster. The authors claim that this work is the first practically usable probabilistic solver, but the impracticality of IEKS was not shown in any way. The results should also contain ablation studies of all parts or combinations of the new ideas.\n\nThe paper proposes well-thought refinements to probabilistic BVP, but the empirical evaluation of the refinements seems incomplete or lacking. \n\n\nOther comments:\n* The learning problem is unclear from page 1, what is it?\n* What do the ball/diamond mean in fig 1?\n* The paper’s notation is a bit complex. The start of sec 2.1. does not define \\nu. What is a single Y_i? The notation Y : t -> R is unnecessarily complex. \\Gamma does not connect anywhere. What is “integrated” wiener (reference)? How come Y has a derivative if its Wiener?  \n* Why is eq 5 a posterior? It looks just a conditional distribution since we haven't placed any priors to the losses \\ell. Similarly this does not look like MAP, but constrained global optimisation.\n* What is p(Y(t) | ..)? Please define\n* Why is \\ell_n nonlinear? Eq 4 looks linear. \n* eq 7 drops the zero-conditioning, how do we condition here if not 0?\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations could be clearer.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "3zghRhpLTkK", "reviewer": "Reviewer_iTr9", "summary": "This paper introduces a probabilistic boundary value problem solver that inherits much of its background from the similar  Gaussian filtering paradigm developed over several recent papers for IVP solvers. While the method initially seems like a minor extension to these, the paper discusses initialisation, calibration, mesh refinement to achieve pre-set tolerance level, things which are often neglected in these types of works. There are encouraging results (though mutliple approaches for eg. initialisation are presented and its not clear which is recommended in given cases) and the claim is made that the computational cost of the algorithm is competitive with state of the art classical solvers.", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This is a strong paper which attacks an important problem (BVPs) within the relatively new yet increasingly established probabilistic numerics paradigm. The work is rigorous and considers at length practical issues such as initialisation and computational cost (though in the latter case benchmarking is not provided, and neither is code) while also undertaking a detailed discussion of calibration, hyperparameter setting, and differences in prior choice. These issues are often neglected in similar papers. The detailed description of the model in Appendix A is informative and well written and it is a pity it does not fit into the main paper. I particularly appreciated the focus on the use of the method in practice (including what seems from the description like carefully designed code) though obviously the code has not been made available. I also appreciated the clear writing. I did wonder given the claims of 'practical usability' why no benchmarking was attempted.\n\nI do have some questions and clarifications for the authors:\n* Can this method easily extend to the situation where the boundary conditions are random (with known distribution, say, but not simply a fixed point). Would this also apply to the revised model in Sec 3.2 ?\n* Fig 3. Please provide some more detail on what happens when the number of points is increased. Is the poor EKS linearisation using the guess start made up for with more points, or is it long term catastrophic? Are there any issues with the EKS scheme (beyond the additional cost), say the failure to find alternative solutions where more than one exists etc? Some comment would be good.\n* It’s not discussed how to set \\nu. What is the intuition with the higher derivatives? I can see in Fig 6 that you have varied this parameter but how does the practioner decide? This seems like the major critical algorithm parameter not discussed.\n*Is BQ sufficiently accurate in practice with such a small number of nodes? How is the BQ uncertainty propagated to the main IEKS algorithm? That’s not entirely clear. How can you discriminate the source of error as coming from the ODE or from the BQ within the steps?\n* Calibration (right hand plot of Fig 6) is interesting. This seems to be not only quite poor (underconfident as stated in the main text) but also getting worse as the number of grid points increases. Are the strategies within the current algorithm (a different approach to setting certain parameters for example) that would mitigate this? How should a practitioner who presumably wants ‘ a well calibrated method’ interpret this plot?\n* Did you compare your method (in wall clock time rather than in computational complexity) to existing solvers including state of the art classical ones? Can you give anything quantitive?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "No plausible negative societal impacts from this work so not addressed.", "review_text": "", "extra_scores": {"ethical_concerns": "No plausible negative ethical concerns from this work.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2.5"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes and studies a probabilistic numerical method for the solution of boundary value problems (BVPs) specified by ODEs.  The contribution builds on earlier work that exploited the state-space formulation of Gaussian processes (GPs) to solve initial value problems specified by ODEs.  Considerable attention is paid to practical implementation and I believe the author(s)' claim that the method is useable out of the box (the author(s) promise that code will be released).  The NeurIPS community have historically been interested in the use of GPs to solve ODEs, and I expect that this work (which is the first practical method for BVPs) will again attract significant interest.
Decision: Accept (Poster)