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": "B83n8FqNv7c", "reviewer": "Reviewer_Cw4U", "summary": "This paper extends filtering based probabilistic numerical methods for solving ODE initial value problems to the solution of boundary value ODE problems by introducing boundary value constraints through the likelihood.  ", "questions": "", "limitations": "", "rating": 4, "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 application is an important one, as pointed out in the introduction. However, the mathematical details provided in the paper seem insufficient to fully understand the proposed approach. \n\nAbstract\n\nI could not verify the claim in the abstract that \"In contrast to previous work, we introduce a Gauss–Markov prior and tailor it specifically to BVPs\", because the prior (solution to equation 2) is not directly constrained at the boundaries. In fact, it seems to me that the BVP constraints are incorporated through the likelihood (equation 3) rather than directly through the prior. However, the prior is often referred to as a \"bridge prior\" in the paper, which makes me wonder if I have misunderstood something. \n\nIntroduction \n\nThe authors should discuss the difference between their approach and other existing PNN methods for solving BVPs. Two that I am aware of are Bayesian Solution Uncertainty Quantification for Differential Equations by Chkrebti et al. and Bayesian Probabilistic Numerical Methods by Cockayne et al.\n\nSection 2\n\nThe prior should be defined outside of section 2.1, since \"Generative model\" typically refers to the likelihood rather than the prior.\n\nThe prior specification is incomplete. The prior is a diffusion process, specifically it is the solution to the SDE in equation (2). However, the form of the drift and diffusion terms are not specified. A question that naturally comes up is why a diffusion process is appropriate for modeling a smooth ODE solution. There is a vague suggestion that some of the components of the vector Y are actually integrated diffusion processes, which may at least guarantee that the prior over the solution is differentiable. The paper states that \"Other [prior] choices are possible\" without elaborating. It seems that other approaches in the literature may be more natural, for example jointly modeling the solution and derivatives through a Gaussian process prior allows better control of the smoothness. At a minimum I believe that a discussion of the prior and its features and drawbacks is required.  \n\nSection 3 \n\nLack of clarity in section 2 makes the rest of the paper difficult to understand. The details of the filtering algorithm seem likewise difficult to understand for a general mathematical audience. \n\nIn general, the paper is not self-contained. It references specific equations from various other papers that many readers will not have the time to track down. For example, the sentence starting on line 53 on page 2, \"The entries in A ... are specified by the integrated Wiener process and imply that $Y_q$ models the $q$th derivative of the BVP solution\" is incomplete. The reference to  \"[6, Equation 2]\", which presumably provides the form of A and an explanation of the derivative modeling, requires reading a different paper to be able to understand the proposed model. This happens several times throughout the paper.   ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Unsure", "review_text": "", "extra_scores": {"ethical_concerns": "No", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "6nVB6urUpy5", "reviewer": "Reviewer_JUUq", "summary": "The paper proposes a probabilistic method to solve ode with boundary conditions. The. method is Bayesian, where the posterior distribution is computed with a Gaussian process as the prior, conditioning on satisfying the ode and boundary conditions. The posterior distribution is approximated with a Gaussian distribution where the mean and variance are obtained from implementing an extended Kalman smoother. The paper discusses at length how to select the algorithm hyper parameters, such as initialization and mesh and illustrates the algorithm on several numerical examples.          ", "questions": "", "limitations": "", "rating": 4, "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": "I am not sure about significant novelty introduced compared to prior work on probabilistic methods to solve BVPs that are referenced. For example, the claims in [12] are quite similar to the claims of this paper: a probabilistic method based on a Gaussian posterior that performs on par with deterministic methods.  Both methods seem to involve learning Gaussian processes. What are the significant and new ideas presented in the paper to distinguish it from prior work? From the way the paper is written, it seems like another probabilistic method way to  solve the BVP, with more discussion on choosing the hyperparameters.  \n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"ethical_concerns": "None", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "3CBWJnQ0kdj", "reviewer": "Reviewer_odbV", "summary": "This paper details a probabilistic method for solving nonlinear first order differential equation systems (ODEs) subject to given boundary conditions. The authors title this a boundary value problem (BVP). The paper builds on existing methods for solving first order ODEs in reference [6] in the paper, and represents the solution of the discretized ODE as the  MAP estimate of a Bayesian posterior. The primary contribution of the paper is a practical algorithm for computing the MAP estimate using the so-called iterated extended Kalman smoother. The primary contributions of the paper is to describe methods to initialize the IEKS, adaptively select the discretization grid and calibrate hyperparameters.", "questions": "", "limitations": "", "rating": 4, "confidence": 3, "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": "Overall, I found the paper to be a little narrow, and the algorithmic contributions to be somewhat narrow.\nOriginality: The paper hews very closely to existing works. The main results on initialization, error estimation, and calibration are rather boilerplate. The utility of the proposed algorithms seems limited beyond the specific setting of this paper.\n\nQuality and Significance: The quality of the results and presentation are lacking. The results, as noted above, are pretty standard, and the proposed algorithm is really a mash-up of a many different algorithms. Overall I would say the results are not significant.\n\nClarity: I found the paper to be rather poorly written. In particular the description of the generative model leaves much to desire. I note that (2) is not the integrated Wiener process. Indeed, if one follows the description in [6], it is pretty clear that the multivariate process Y involves time integrals, and not SDEs. Indeed, I do not believe that the integrated Wiener process in [6] can be expressed as the solution of an SDE, like (2). Ultimately, however (2) doesn’t play a big role in the overall message of the paper, but it does show that the authors have a fundamental misunderstanding of the setting in [6]. Moving beyond this, the overall description of the posterior  is rather poor, and the reader is expected to interpolate the authors descriptions. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "It would be hard to ascribe negative societal impact of solving first order ODEs. I have no feedback to offer on this front.", "review_text": "", "extra_scores": {"ethical_concerns": "None.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2 hours"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper was uneasy to review. The AC asked 10+ ppl if they could review it, but still could only find 2 reviewers. With the help of SAC, we found a 3rd reviewer. Unfortunately, some reviewer didn't fully understand the technicality of the paper, and the last reviewer who knows the subject better did not appreciate the significance of the paper's results. Based on all reviewers' feedback, we can do nothing but reject this paper.

The fact that we had a very hard time finding good-quality reviewers for this paper hints that: the paper's topic is far from mainstream machine learning (which isn't necessarily an issue), but the paper is also not presented in a way to appeal to the machine learning community (which is a real issue). We recommend that the reviewers either submit the paper to venues more focused on differential equations, or substantially revise the paper to suit a machine learning conference.
Decision: Accept (Poster)