Black Box Probabilistic Numerics
Authors: Onur Teymur
Christopher Neal Foley
Philip Gavin Breen
Toni Karvonen
Chris J. Oates
Venue: NeurIPS 2021 Poster (2021)
Keywords: Bayesian Statistics
Probabilistic Numerical Methods
Numerical Analysis
Algorithms
Gaussian Processes

Abstract: Probabilistic numerics casts numerical tasks, such the numerical solution of differential equations, as inference problems to be solved. One approach is to model the unknown quantity of interest as a random variable, and to constrain this variable using data generated during the course of a traditional numerical method. However, data may be nonlinearly related to the quantity of interest, rendering the proper conditioning of random variables difficult and limiting the range of numerical tasks that can be addressed. Instead, this paper proposes to construct probabilistic numerical methods based only on the final output from a traditional method. A convergent sequence of approximations to the quantity of interest constitute a dataset, from which the limiting quantity of interest can be extrapolated, in a probabilistic analogue of Richardson’s deferred approach to the limit. This black box approach (1) massively expands the range of tasks to which probabilistic numerics can be applied, (2) inherits the features and performance of state-of-the-art numerical methods, and (3) enables provably higher orders of convergence to be achieved. Applications are presented for nonlinear ordinary and partial differential equations, as well as for eigenvalue problems—a setting for which no probabilistic numerical methods have yet been developed.

Black Box Probabilistic Numerics
Onur Teymur University Kent
Christopher N. Foley University of Cambridge Optima Partners
Philip G. Breen
Roar Al
Alan Turing Institute
Toni Karvonen
Chris. Oates
University Helsinki Alan Turing Institute
Newcastle University
Alan Turing Institute
Abstract
Probabilistic numerics casts numerical tasks, such the numerical solution of differ- ential equations, as inference problems be solved. One approach is to model the unknown quantity of interest as random variable, and to constrain this variable using data generated during the course traditional numerical method. However, data may be nonlinearly related to the quantity of interest, rendering the proper conditioning of random variables difficult and limiting the range of numerical tasks that can be addressed. Instead, this paper proposes to construct probabilistic numerical methods based only on the final output from a traditional method. A con- vergent sequence of approximations to the quantity of interest constitute a dataset, from which the limiting quantity of interest can be extrapolated, in a probabilistic analogue of Richardson's 's deferred approach to the limit. This black box approach (1) massively expands the range of tasks to which probabilistic numerics can be applied, (2) inherits the features and performance of ffte-offthhear numerical methods, and (3) enables provably higher orders of convergence to be achieved. Applications are presented for nonlinear ordinary and partial differential equations, as well as for eigenvalue problems-a setting for which no probabilistic numerical methods have yet been developed.
Introduction
Probabilistic numerics (PN) has attracted significant recent interest from researchers in machine learn- ing, motivated by the possibility of incorporating probabilistic descriptions of numerical uncertainty into applications probabilistic inference and decision support [1 2122 PN treats the intermediate
calculations performed in running traditional (i.e. non-probabilistic) numerical procedure as data,
3 Conjugate which can be used to constrain a random variable model for the quantity of interest
Gaussian inference has been widely exploited, with an arsenal of PN methods developed for linear algebra II cubature 12 30 optimisation
[ k and differential equations (5) However, nonlinear tasks pose major technical challenge to this approach, as well to computational statistics general, due to the absence of explicit conditioning formulae. Compared to traditional numerical methods, which have benefited from century more of sustained research effort, the current scope of PN is limited. The performance gap is broadly characterised by the absence of certain important functionalitie--addee adaptivity, numerical well-conditioning, efficient use of computational resource all of which contribute to limited applicability in real-world settings.
This article proposes pragmatic solution that enables state-of- oftte-of-theart numerical algorithms to be immediately exploited the context of PN. The idea, which we term black box probabilistic numerics (BBPN), is statistical perspective on Richardson's 's deferred approach the limit (RDAL) [56]. The starting point for BBPN a sequence of increasingly accurate approximations produced by
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

traditional numerical method as its computational budget is increased. Extrapolation of this sequence (to the unattainable limit of 'infinite computational budget') is formulated as a prediction task, to which statistical techniques can be applied. For concreteness, we perform this prediction using Gaussian processes (GPs) [57], but other models could be used. Note that we do not aim to remove the numerical analyst from the loop; the performance of BBPN is limited by that of the numerical method on which is based.
There are three main advantages of BBPN compared to existing methods in PN: (1) BBPN is applicable to any numerical task for which there exists a traditional numerical method; (2) state- of-the-ar art performance and functionality are automatically inherited from the underlying numerical method: (3) BBPN achieves a provably higher order of convergence relative to a single application of the numerical method on which it is based, in an analogous manner to RDAL. The main limitations of BBPN, compared to existing methods in PN, are: (1) multiple realisations of a traditional numerical method are required (i.e. one datum is not sufficient in an extrapolation task), and (2) joint statistical model has to be built for not just the quantity of interest (as in standard PN), but also for the error associated with the output of a traditional numerical method. The capacity of generic statistical models, such as GPs, to learn salient aspects of this structure from data and to produce meaningful predictions over range of real-world numerical problems, demands to be investigated.
The article is organised as follows: In Section we recall classical RDAL. In Section we lift RDAL to the space of probability distributions, exploiting GPs to instantiate BBPN and providing a theoretical guarantee that higher order convergence is achieved by using GPs within BBPN. In Section we present a detailed empirical investigation into BBPN, demonstrating its effectiveness on challenging tasks that go beyond the capability of current methods in PN, while also highlighting potential pitfalls. As part of this we perform comparison of the uncertainty quantification properties of BBPN against earlier approaches. A closing discussion is contained in Section S
4
2 Turning Lead into Gold
Our starting point is the celebrated observation of Richardson [56). that multiple numerical approx- imations can be combined to produce an approximation more accurate than any of the individual approximations. To see this, consider an intractable scalar quantity of interest q" E IR, and suppose that q" can be approximated by a numerical method q that depends on a parameter h > 0, such that
{h) = q" + Ch + 1(h++))
(1)
for some C E I (which may be unknown) and a 0 (which is assumed known, and called the order of the method). Clearly ((h) converges to q" as h + 0, but we also suppose that the cost of computing (h) increases in the same limit, with exact evaluation of q()) requiring a hypothetically infinite computational budget. Proposition Tis ii the cornerstone of RDAL. It demonstrates that two evaluations of a numerical method of order a can be combined to obtain a numerical method of order a + 1. An elementary proof of this foundational result is provided in Appendix A
Let q be a numerical method of order a, as in ①.. Fix E (0.1) 1) and let q- (h) denote the height at which a straight line drawn through the points (h" q()) and ((yh)" (((h)) intersects the vertical axis in R.2 Then by is a numerical method of order a + 1.
Proposition 1.
Now consider the natural generalisation of Proposition I in which we compute approximations (h+,) along decreasing sequence of values (ha))= One can then fit a smooth interpolant the points {(hg, (0:,((hh))))) give an estimate for the quantity of interest. This simple idea is widely used in numerical analysis; its potential to radically improve solution accuracy, given only a sequence of simple calculations as input, prompted Press al. (58] p. 922] describe it as "turning lead into gold". The practical success of RDAL depends on the choice of interpolant, with polynomial interpolation being most commonly used. Unqualified, RDAL usually understood to refer to an order n polynomial fitted to points, which produces numerical method of order a see Theorem 9.1 of [59]. Higher-order polynomial extrapolation is known perform poorly unless the values (h,),,, are able be chosen specifically to mitigate Runge's 's phenomenon (60]. motivating the Bulirsch-Stur algorithm [6].. which instead fits a rational function interpolant. This allows both greater expressiveness and robustness than polynomial interpolation (though not necessarily as efficiently (58)). These methods are all situated within the broad category of extrapolation methods in numerical analysis; a comprehensive historical survey can be found in 62
(generalising the straight line through two points), then extrapolate this to / = 0, to

---

100: 5,00
1801 4,00
100: 5,00
0.08
3.04
h- 0.02
0.01
0.5
0.5
0.5
1.9
1855
1.05
18
1.8
1.25
1.75
1.7
1.85
1.5
Unear
Quadratto
Postericr mosn
Cubic
Peataricr main al-te
Data
1.45
Posterior moan sh-2e
1.4
1.4
0
0.01
0.00
0.00
0.04
0.08
0.05
0.07
0.08
0.08
0.01
0.02
0.00
0.04
0.05
0.05
0.07
0.08
0.08
0.1
h
Figure 1: Richardson's 's deferred approach to the limit (RDAL), applied to the method of Riemann sums with integration bandwidth h and oscillatory integrand fxx) displayed in the four top panes. The bottom left pane shows linear, quadratic and cubic interpolants converging on the true value of the integral, denoted by blue star. The strength of RDAL is seen in the fact that the cubic interpolant gives a better estimate than that given by the finest- -rii Riemann sum with h = 0.01. (This example is closely related to Romberg's 's method.) The bottom right pane illustrates black box probabilistic numerics (BBPN), in which a Gaussian process is fitted to the same data. The GP specification is crucial the performance of BBPN, and described in Section 3
Figure Q presents a simple visual demonstration of RDAL, applied to the method of Riemann sums for oscillatory ID integrand. While RDAL gives improved approximations, estimation uncertainty is provided. The only attempt of which we are aware to provide uncertainty quantification for RDAL is due to [63, who focused on the Navier-Stokes equation and a specific scalar quantity of interest. Here we go further, proposing the general framework of BBPN and introducing novel methodology that goes beyond scalar quantities of interest. The right-hand pane of Figure [ displays the outcome of the method we are about to introduce, applied to the same lask-observe that the true value of the integral falls within the +:2g credible produced using BBPN. Details of the simulations this figure are contained in Appendix C.T
quantification of
3
Methodology
The core idea of BBPN to model (h)) as a stochastic process Q(h) rather than fit a deterministic interpolant in RDAL. The distribution of the marginal random variable ((0) is then interpreted as a representation of the epistemic uncertainty in the quantity of interest (0). Conjugate Gaussian inference be performed in BBPN, since one needs only to construct interpolant. This means the challenge of nonlinear conditioning encountered in PN is avoided, massively extending the applicability of PN. In addition being able leverage state-of-the-an numerical methods, the BBPN approach enjoys provably higher orders of convergence relative to a single application of the numerical method on which it is based; see Section 3.3
3.1 Notation and Setup
Our starting point is to generalise ① to encompass essentially all numerical tasks, following the abstract perspective [4 To do so, we observe that any quantity of interest g" can characterised by sufficiently large collection of real values q" (t), with t ranging over appropriate index set T.
3

---

Definition 2. A traditional numerical method is defined as a map q ho) X T - R. for some ho > 0 such that, for all E T, the function h + q(,, t) is continuous at with limit q(0, (t).
For example, (univariate) initial value problem (IVP), in which is interpreted as time, can be solved using a traditional numerical method whose time step size h trades off approximation error against computational cost. The output q(h, t) of such method represents an approximation to the true solution q (t) of the IVP, at each time for which the solution is defined. In general, depending on the numerical task, the index could be tempora,,,,,,,, discrete, even unstructured set, while the meaning of the index h will depend on the numerical method.
Definition 2 thus encompasses, among other things: (1) adaptive integrators for time-evolving partial differential equations (PDEs), where > O represents a user-specified error tolerance, and the spatio-temporal domain of the solution is indexed by T; (2) iterative methods for approximating the singular values of a d x d matrix, where for example h := w with w the number of iterations performed, and the ordered singular values are indexed by T { .l.... d]: and (3) the simultancous approximation of multiple related quantities of interest, where T indexes not only the domain(s) on which the individual quantities of interest are defined, but also the multiple quantities of interest themselves (this situation arises, for example, in inverse problems that are -DE-ttrstrreeed (65).
The perspective in Definition 2 is abstract but, as these examples make clear, it will typically only be possible to compute q at certain input values (such as h - for w E N. or for just finite collection of inputs if the index set 7' is infinite), and furthermore each evaluation is likely to be associated with a computational cost. Thus complete information regarding the map q : [0, ho) x T - R will not be available in general, and there will therefore remain epistemic uncertainty in its complete description. Our aim in Section 3.2 in line with the central philosophical perspective of PN, is to characterise this uncertainty using a statistical model.
3.2 Black Box Probabilistic Numerics
The proposed BBPN approach begins with a prior stochastic process Q and constrains this prior using data D. Concretely, we assume that the real values (1,z) are provided at a finite set of resolutions h1 > ..... >> hn > 0 and distinct ordinates t,,,..... fi,m, E T. Note that the number of t-ordinates m, can depend on h4. Our dataset therefore contains the following information on q:
D= ((he, g(ee, 44))): i- i.,,...,,,
., m;}
The stochastic process obtained by conditioning Q on the dataset D, denoted D, implies a marginal distribution for Q(0, -). which we interpret as a statistical prediction for the unknown quantity of interest q" (-). In order for uncertainty quantification in this model to be meaningful, one either requires expert knowledge about the numerical method that generated D. or one must employ a stochastic process that is able adapt the data, so that its predictions can be calibrated.
Our goal is to specify a stochastic process model Q(h, t) that behaves in a desirable way under extrapolation to h 0. To this end, we decompose
Q(h,t) Q1) + E(h,
where Q² (t) is a prior model for the unknown quantity of interest q" (t), and F(h,,t) is a prior model for the error of the numerical method. It will be assumed that Q" and E are independent (denoted Q' IL E). meaning that prior belief about the quantity of interest is independent of prior belief regarding the performance of the numerical method. (This assumption is made only to simplify the model specification, but if detailed insight into the error structure of numerical method is available then this can be exploited.) Compared to the existing PN methods cited in Section I a prior model for the error E is an additional requirement in BBPN.
The error E(h, 1) assumed to vanist as - O. meaning that a stationary stochastic process model for Ehh, f), and hence for Q(h, t), is inappropriate, and can result in predictions that are both severely biased as well as under-confident; see Appendix C.2 In the next section, we propose a parsimonious non-stationary GP model for Qhh, t) of the form which combines knowledge of
1 'This statement covers several potentially subtle notions from numerical analysis such as well posednesss of the problem and numerical stability of the algorithm: these are studied in detail in their own right in the literature, and for purposes it suffices to assume that the error behaves well in the limit.

---

the order of the numerical method (only) with setting is practically relevant- the order of a numerical method is typically one of the first theoretical properties that researchers aim to establish while, conversely, for more complex numerical methods the order may actually be the only salient high-level error-characterising property that is known, and thus represent the limit of mathematical insight into the method.
estimation of GP hyperparameters. This
3.3 Gaussian Process BBPN
Gaussian processes provide a convenient model for Q* and E, since they produce an explicit form for the conditional Q) D. The details of conjugate Gaussian inference are standard (see e.g. (57] and so relegated to Appendix B.T our foeus here is the specification of GP priors for Q* and E.
The notation Q ~ GP (1aq kQ) will be used to denote that Q is GP with mean function diQ and covariance function kq. With no loss of generality, in what follows we consider centred processes (i.e. μ2 = 0). It will be assumed that T = Ti X x Tp, with each Z; either a discrete or a continuous subset of Euclidean space, with the Euclidean distance between elements 44, E Ti being denoted ||t, - l- (Typical applications involve small p: for example, the domain of spatio-tempora PDE is typically decomposed as T = T2 T2 where Ti indexes time and T2 indexes all spatial dimensions, so that 2.)
Prior for Q": In the absence of detailed prior belief about q² we consider the following default prior model. Let G ~ paki),, 2 = (Z1,..., .., Z,)) ~ N(0,02T) and let Z 1L. G. Let bi,...,b by be finite collection of basis functions and set b(t) = ((()),..,,,,(())) Then set
('(t)= Z b(t) + G(t)
where 0², PG > are parameters be estimated. The basis will be groblem-specinee and could be a polynomial basis, Fourier basis, or any number of other bases depending on context. The case v = 1 with a constant intercept is closely related to ordinary kriging and the case e >> 1 closely related to universul kriging [66] p. 8). The apparent redundancy in the parameterisation due to the product a² PG will be explained later. Using the notation = (4,,... (p), we consider tensor product covariance model ka(t, t') IT 1 KG. (t,,,t),) kg. (t.,
(llt; (e,), for some radial basis functions 6., scaled to satisfy Φ; (0) = 1, and length parameters E,,5 > 0 to be estimated.
Prior for E: The process E(h, t) is model for the numerical error g(h, t) qt), t > O, which may be highly structured. flexible prior model is therefore required. Moreover the error will. by definition. depend on the order of the numerical method: for successful extrapolation we must therefore encode this order into the model for E. It was earlier argued that a stationary GP is inappropriate, since the error is assumed to be (hh)). However, we observe that h=" (g(h, () 10)-gga q" (()) ((), suggesting that this quantity can be modelled using a stationary GP. We therefore take E ~ GP a² pgkg), where PE > 0 is a parameter to be estimated, and k t), (l',t')) = (hh')" & (// - h1l/()) (t,
(4) for a radial basis function vi, scaled to satisfy (0)) = 1, and a length-scall parameter é, > 0 to be estimated. Note how 4 separates the h and dependence of E' in the prior, and adopts the same covariance model kG that was used to model the t dependence of G. This can be motivated by the alternative perspective that follows from observing that is GP with covariance function
ko((h, f), (K, ()) (K)) o² { b(t)-b(t')
k f), (h', (H',())) (1,00) b(t') + (tk (1 PE kt(t,t')
}
(5)
where km/ke, is a kernel only depending on h. Written this way, the model is seen perform universal kriging over T with a covariance adjusted by a multiplicative error arising from non-zero values of h.
Higher-order convergence:
The GP specification just described is not arbitrary; it ensures that
the higher-order convergence property of RDAL is realised in BBPN. Consider again the setting in
Proposition I Suppose that there exist L, Eu > 0 and 3 E (0, 1] such that [1 (==) < 8 for all
8 E (0,))) Then the posterior mean E [x] satisfies 14* [[Q(O)) = O(h + as + 0. Thus if vé is Lipschitz (i.e. B = 1). BBPN achieves the same higher-order convergence, + 1, as RDAL. In this context we recall that any Matérn covariance function of smoothness at least 1/2 is
Lipschitz. The proof is provided in Appendix B.2
5

---

Model parameters: The free parameters of our prior model, o², PG, PE. En. and the 4,, for i I....... .. are collectively denoted 0. For all experiments in this article, was set using the maximum likelihoodd estimator 0M. Our choice of parameterisation ensures that the maximum likelihood estimate for the overall scale o², denoted La has a closed form expression terms of the remaining parameters. This is analytically derived in Appendix B33 and can be plugged straight into the likelihood. Gradients with respect to the remaining 3 + p parameters are derived in Appendix [B.3 and gradient-based optimisation of the log- oggliike was implemented for the remaining parameters.
Remark GP interpolation, as with classical RDAL, is not parameterisation invariant. Thus some care is required to employ a parameterisation of h that is amenable to the construction of a GP interpolant. The effect of differing parameterisations explored in Appendix C.2
Remark 4. The classical definition of RDAL presupposes that, in order to employ the method, the order a must be known a priori 167]. However if a not known, the probabilistic perspective affords us the opportunity to learn as an additional parameter in the statistical model-a procedure with no classical analogue. The feasibility of learning explored Section
Code: Software for BBPN, including code to reproduce the experiments in Section can be downloaded from gi thub.ccom wwt.c//ooo///t///tcc/ee /ot////t///bbbbb// bbbn..
4 Experimental Assessment
This section reports a rigorous experimental assessment of BBPN. Firstly, Section demonstrates that BBPN is competitive with existing PN methods in the context of ordinary differential equations (ODEs). This result is somewhat surprising, given the black box nature of BBPN compared to the bespoke nature of existing PN methods for ODEs. Secondly, in Section 4.2
demonstrate the versatility of BBPN by applying it to the nonlinear problem of eigenvalue computation, for which no PN methods currently exist. Finally, in Section 4 we use BBPN to provide uncertainty quantification for state-of-the- the-art numerical methods that aim to approximate the solution of nonlinear PDEs.
Default Settings: We use Matérn(1/2) kernels for Φ; and vù, i.e. (t,t,, f)
exp(-IIt; /
and similarly mutatis mutandis for gb. These kernels impose a minimal continuity assumption on without additional levels of smoothness being assumed. Sensitivity of results to the choice of kernel investigated in Appendix C.2
Performance Metrics: PN is distinguished from traditional numerical analysis by its aim to provide probabilistic uncertainty quantification, but nevertheless approximation accuracy remains important. To perform an assessment on these terms, we considered two orthogonal metrics. Firstly, we compute the error of the point estimate (mean), denoted W == |IE[Q -) D] q-)(l where the norm is taken over E T" where T' is either T itself or set of representative elements from T. Secondly, and most importantly from the point of view of PN, we consider the surprise S = [ICQ((0,)) - (E(Q(O, [)] - q" (-))l, where CIQ(0,-)) D] denotes the posterior covariance matrix. If the true quantity of interest g" was genuinely sample from DD then S" would follow x2 distribution with ' degrees of freedom. This observation enables the calibration of PN method to be assessed 68. Both metrics naturally require an accurate approximation to q* to act as the ground truth, which is available using brute force computation in Sections 4.1 and 4.2 but not in Section 4.3 The role of Section 4.3 is limited to demonstrating BBPN on problem class that is challenging even for state-of-the-ue methods.
4.1 Ordinary Differential Equations
The numerical solution of ODEs has received considerable attention in PN, with several sophisticated methods available to serve benchmarks. Here we consider numerical solution of the following Lotka-Voltere IVP, a popular test case in the PN literature: dy di = f(t.t)
0.5y/ - 0.05 //
20
y()) = 20
0..5g2 + 0.05y1y2
2 AAltemattve approaches, such as cross-validatio could also be used; see Chapter 5 of [57 Our choice of maximum likelihood was motivated by the absence of any degrees of freedom (such as the number of folds of cross-validation), which permits more objective empirical assessment.

---

Chiowbiti at at (2016)
Teymur al. (016))
Schober at at. (P018)
RAPN
15
20
10
20
10
15
20
10
15
30
25
20
-2
OVERCONFIDENT
lng. W
-6
100
Olkr.
-10
Teym.
Can: 01
Noach
88PN 01
8BRN Oz
-10
UNDERCONFOENT
-18
-15 -1.5 -2
-1.5
-2.5
-3.5
-4.5 -5.5
-2.5
-3.5
-6.
-55
1og
log_h
Figure 2: Ordinary differential equations. Top: Output from three existing PN algorithms and BBPN, applied to the
33
IVP. Bottom left: The error log2 W at the final time
point fend 20, a function of the time step size h. Bottom right: The surprise log, at
lend 20, with the central 95% probability band of Xå random variable shaded. Methods shown with (where applicable) their order: Chkr. [3 Conr. O1 [38]: Teym. 02 [1 Scho. O1
40 Tron. 02 47
Bosch 02 55
BBPN 01 & O2; and (traditional) Euler.
The aim in what follows approximate the quantity of interest q" ((end) for fend = 20. The top row of Figure displays output from three distinct PN methods due to [9,411 as well as BBPN. (For these plots the coarse step-size h. 0.5 was used, so the probabilistic output can be easily visualised.) In each case, these methods treat a sequence of evaluations of the gradient field f as data which are used to constrain a random variable model for the unknown solution of the ODE. Their fundamentally differing character makes direct comparisons challenging, particularly if we are to account for computational cost. However, each algorithm has a recognisable discretisation parameter h, it remains instructive study their h 0 limit. (In most cases h represents time step size, but the method of [55) is step-size adaptive; in this case h is an error tolerance that is user-specifiedd) The methods of 38 and thus have significant computational cost. The methods of [40] Gaussian filtering and are less computationally demanding, though in disregarding nonlinearities the description of uncertainty they provide not as rich. Interestingly, the output from 39 becomes overconfident h 0, with S being incompatible with a xá random variable, while the output from 40 becomes somewhat pessimistic in the same limit. Aside from these two outputs, the other PN methods considered appear be reasonably calibrated.
30 and [1] require parallel simulations to produce empirical credible sets, and [55] are based on
To illustrate BBPN, our data consist of the final states produced by either an Euler (order 1) or an Adams-Bashfoon (order 2) algorithm, which performed at different resolutions (h, 2 " 1....6). The datase is augmented cumulatively, so that for 1, 1...., 1' are used. The finest resolution in each case, h, is simply denoted h. For this experiment we use a prior with constant intercept, i.e. v= and by (t) = 1. The BBPN output, shown in the bottom row of Figure 2 is observed to be calibrated, and the (order 2) output provides the most accurate approximation among all calibrated PN methods considered. Note in particular how BBPN accelerates the convergence of the Euler method from first order second order, akin to RDAL.
i', all data generated by runs
'In this experiment the two components of were treated as priori independent, but this is not specific requirement of BBPN and dependence between outputs in principle also be encoded into the GP model.
7

---

10
Data =
6.5
7.5
2006-0
5.5
QR 1m3,m,mm.
m
OR: Im 10, mm
(20O
Shifted Power Methed
3.5
0.2
0.4
0.6
0.8
0.2
0.4
0.8
0.8
0.05
0.1
0.1
0.2
0.3
h
Figure Eigenvalue problems. Left and centre: QR algorithm. Right: Shifted power method. Details of each simulation are given in the main text. All plots show red shaded 12 credible intervals, numerical data black circles, and eigenvalues as blue stars.
In terms of computational cost, BBPN requires running a traditional numerical method at different resolutions, as well the fitting of GP. In this experiment, the computational cost of BBPN was intermediate between the filtering approach of [40 and the sampling approaches of 39 and Further details, including the sources of all these codes, are given Appendix C.3
4.2 Eigenvalue Problems
The calculation of eigenvalues is an important numerical task that has not yet received attention in
PN. In this section we apply BBPN to (1) the QR algorithm for matrix eigenvalue discovery, and (2)
a cutting-edge adaptive algorithm for the tensor eigenproblem, called the shifted power method [69
In these examples, the order is unknown and append it to estimated using maximum likelihood.
additional parameter to be
QR Algorithm: Take T f1..... n). Let w E N and define h w-1 Given matrix A E IR xn its Schur form Aoo is approximated by matrices A. == Ry 1Q Qw- where Rw and result of performing QR decomposition An = A. Then qhh, is the vector diag( A1-1). and o² = diag( (iag(AA)) ) is the vector of eigenvalues available in closed form (see Appendix [C.4. consider the following family of sparse matrices that arise as the discrete Laplace operator in the solution of the Poisson equation by finite difference method with five-point stencil. Let the / matrix B and the ml x ml matrix A be defined by
Q w arise
A Qw- Rw and where
of A. As test problem, whose eigenvalues
-1
-1
-
4
-1
- I
B
-1
B =
A
- B
BBPN output for this problem displayed Figure the left-hand pane, we take / 5, m and perform 5 QR iterations, displaying all 10 eigenvalues. In the centre pane, we take = 10, 10 and perform
iterations. For clarity, this pane only displays the largest few eigenvalues of this 100 x 100 matrix, and also show zoomed-in crop better demonstrate the extrapolation quality. Both examples show the convergence of Q(,,)) -) to q" (.) as w + oo. Recall that a is inferred in these simulations -the maximum likelihood values were, respectively, 1.0186 and 1.0167.
The extrapolation performed by our GP model is seen visually be effective and almost all true eigenvalues are contained within the #20 credible intervals plotted. For comparison, in Appendix C2 we contrast the result using stationary GP model (i.e. a = 0). The extrapolating properties of that GP are immediately seen to be unsatisfactory, and we support this observation by examining
calibration the two approaches, similar manner to Section 4.1 This analysis strongly supports proposed GP specification in Section

---

This iterative algorithm, due to 69 and implemented in 70 finds (random) eigenpairs of higher-order tensor systems, and we include demonstrate BBPN on challenging problem in linear algebra. For A a symmetric mth-order n-dimensional real tensor, and x an n-dimensional vector, define
Shifted Power Method:
(Ax m 1
Q11 Rin Pim i1 I,...,
and say that À R is an eigenvalue of if there exists E R" such that Ax"'! Ax x 1. Here we take n = m = 6 and produce a random symmetric tensor using the create prstteeprollem function of [70]. Two parameterisations of were considered; h w and h == w 2 where w denotes the number of iterations performed, with results based on the latter parameterisation presented in the right-hand pane of Figure (The maximum likelihood value for in this example was 1.3318.) It can be seen that, after iterations, BBPN is more accurate than each of the individual approximations on which it was trained. The choice of parameterisation affects the performance of BBPN, an issue we explore further in Appendix
In this example there is additional computational cost to BBPN in the data collection stage, since the dataset is generated during a single run of an iterative numerical method. Therefore the only overhead is due to fitting the GP; though for this example this cost is itself negligible. All details, including a systematic assessment error W and surprise S as is varied for both the QR algorithm and the shifted power method, are given in Appendix
Remark 5. In this section we have implicitly modelled eigenvalues as priori independent, for simplicity of exposition. Heuristics from random matrix theory suggest that, when treated probabilis- tically, eigenvalues may be better modelled with some non-triviee non-trivial dependence structure. We note that this additional structure can be easily incorporated into the prior GP model for BBPN.
4.3 Partial Differential Equations
To demonstrate the potential of BBPN on challenging problem for which state-of-the-art numerical methods are required, we consider numerical solution of the Auramoto-Sivashing equation [7] 72
d.u. 8% + 8z + ud_u 0.
This equation is used to study a variety of seaction-diffusiio systems, producing complex spatio- temporal dynamics which exhibit temporal chaos- the characteristics of which depend strongly on the amplitude of the initial data and the domain length. We consider flat initial condition u(x; 0) with periodic boundary condition on the domain 0 S x S 1, and numerically solve to t 200. Our quantity of interest is therefore u(x,20 200) for the domain
[0, 1].
To obtain numerical solutions to O we transfer the problem into Fourier space and apply the popular fourth-order diffeeiffernniigg ETD RK4 numerical scheme; see Appendix C.5 and [73]. ETD RK4 was designed to furnish fourth-orda order numerical solutions to time-dependent PDES which combine low- order nonlinear terms with higher-order linear terms, as in (6). Computing accurate approximations to the solution of chaotic, stiff PDEs challenging problem for existing PN methods because computationally demanding high-order approximations across both spatial and temporal domains are required. Here, we assess BBPN applied to three sequences of five runs of ETD RK4, with minimum temporal step size h 5l and, for simplicity, fixed spatial step size dx = 0.001 throughout f A reference solution was generated by taking h = 0.0005, but this cannot of course be guarantced to be an accurate approximation to the true solution of (
Results shown in Figure are encouraging; not only can accurate approximations be produced, but the associated uncertainty regions appear to be reasonably well calibrated, insofar as the magnitude of the uncertainty consistent with the magnitude of the discrepancy between the posterior mean and the reference solution. Full details of these simulations are contained in Appendix C
*For the h = 0.002 simulation in Figure 4 we have h; E (0.002, 0.005, 0.01, 0.02, 0.05). for the ht 0.005 simulation we have h. E (0.005, 0.01, 0.02, 0.05, 0.1}. and for the h = 0.01 simulation we have h, E {0.01, 0.02, 0.05, 0.2)

---

P. C.0O2
Porto nin
h 0.002
0.1
0.2
03
0.4
0.5
0.8
0.7
0.8
0.006
u9.200
Plalkrm
103
120
163
0.01
0002000
180
0.01
0
0.2
014
0.6
0.8
011
01.2
0.31
17.4 0.5
0.55
11.7
0.81
0.8
Figure 4: Partial differential equations. Left: Solution to the Saaamott--Svvshinskee equation. Right: Approximation of the solution the final time point 200) using BBPN, based on minimum time step sizes f E {0.002, 0.005, 0.01}. Posterior mean (blue) and credible regions (shaded) are displayed. reference solution (dashed black) is obtained by taking h 0.0005.
5
Discussion
This paper presented black box probabilistic numerics, a simple yet powerful framework that bridges the gap between existing PN methods and the numerical state-of-the-ard, Positive results were presented on the important problems numerically approximating ODEs, eigenvalues, and PDEs. Our main technical contribution is probabilistic generalisation of Richardson's 's deferred approach to the limit, which may be of independent interest.
The main drawbacks, compared to existing PN, are possibly increased computational cost and the additional requirement to model the error traditional numerical method. Compared to existing PN, in which detailed modelling of the inner workings of numerical algorithms are exploited, only the order of the numerical method used in BBPN (and we can even dispense with that, as in Section which may reduce its expressiveness in some settings. However, despite the black box approach, BBPN was no less accurate than existing PN in our experiments, and in fact the higher-order convergence property may enable BBPN to out-perform existing PN.
Some avenues for further research (that we did not consider in the present article) include the use of more flexible and/or computationally cheaper alternatives to GPs, the adoption of principles from experimental design to sequentially select resolutions h; given an overall computational budget, and the simultaneous use of different traditional (or even probabilistic) numerical methods within BBPN.
Acknowledgments and Disclosure of Funding
This work was supported by the Lloyd's Register Foundation programme data-centric engineering the Alan Turing Institute, UK. The authors wish to thank Jon Cockayne and lsse Ipsen for feedback on earlier versions of the manuscript, and Nicholas Krämer for assistance with the probnum package.
10

---

References
[1] Hennig, M. A. Osborne, and M. Girolami, "Probabilistic numerics and uncertainty in compu-
tations," Proc. R. Soc. A, vol. 471, no. 2179, p. 20150142, 2015.
[2] C. J. Oates and J. Sullivan, "A modern retrospective on probabilistic numerics," Stat. Comput., vol. 29, no. 6, pp. 1335-1351, 2019.
[3] Cockayne. C. J. Oates,
Sullivan, and M. Girolami, "Bayesian probabilistic numerical methods," SIAM Rev., vol. 61, no. 4, pp. 756-789, 2019.
[4] J. Cockayne, J. Oates, I. C. Ipsen, M. Girolami, et al., "A Bayesian conjugate gradient method (with discussion)," Bayesian Anal., vol. 14, no. 3, pp. 937-1012, 2019.
[5] S. Bartels, J. Cockayne, I. C. Ipsen, and P. Hennig, "Probabilistic linear solvers: unifying view," Stat. Comput., vol. 29, no. 6, pp. 1249-1263, 2019.
[6] Wenger and P. Hennig. "Probabilistic linear solvers for machine learning." NewrIPS 34, 2020.
[7] Hennig, "Probabilistic interpretation of linear solvers," SIAM J. Optim., vol. 25, no. pp. 234-260, 2015.
[8] T. W. Reid, I. C. Ipsen, J. Cockayne, and C. J. Oates, "A probabilistic numerical extension of the conjugate gradient method," arXiv:2008. 03225, 2020.
[9] F. Schäfer, T. J. Sullivan, and H. Owhadi, "Compression, inversion, and approximate PCA of dense kernel matrices at near-linear computational complexity," Multiscale Model. Simul., vol. 19, no. 2, pp. 688-7366 2021.
[10] S. Bartels and P. Hennig, "Probabilistic approximate least-squares," Proc. 19th Int. Conf. Artificial Intelligence and Statistics, vol. 51, pp. 676-684, 2016.
[11] J. Cockayne, C. Ipsen, C. J. Oates, and T. W. Reid, "Probabilistic iterative methods for linear systems," arXiv:20/2. /2615, 2020.
[12] Diaconis, "Bayesian numerical analysis," Statistical Decision Theory and Related Topics IV, vol. 1, pp. 163-175, 1988.
[13] A. O' O'Hagan, "Some Bayesian numerical analysis," Bayesian Statistics 4: Proc. 4th Valencia International Meeting, pp. 345-363, 1992.
[14] M. Fisher, C. Oates, C. Powell, and A. Teckentrup, "A locally adaptive Bayesian cubature method," PMLR, vol. 108, pp. 1275-1275, 2020.
[15] Prüher and Särkkä, "On the use of gradient information in Gaussian process quadratures." in IEEE 26th Int. Workshop Machine Learning for Signal Processing, 2016.
[16] A. Gessner, J. Gonzalez, and M. Mahsereci, "Active multi-information source Bayesian quadra- ture," PMLR, vol. 115, pp- 712-721, 2020.
[17] T. Karvonen, S. Särkkää and C. J. Oates, "Symmetry exploits for Bayesian cubature methods," Stat. Comput., vol. 29, no. 6, pp. 1231-1248, 2019.
[18] H. R. Chai and R. Garnett, "Improving quadrature for constrained integrands," PMLR, vol. 89, pp. 2751-2759, 2019.
[19] R. Jagadeeswaran and F. Hickernell, "Fast automatic Bayesian cubature using lattice sampling." Stat. Comput., vol. 29, no. 6, pp. 1215-1229, 2019.
[20] T. Karvonen and S. Särkkä, "Classical quadrature rules via Gaussian processes," IEEE 27th Int. Workshop on Machine Learning for Signal Processing, 2017.
[21] T. Karvonen, C. Oates, and S. Särkkä, "A Bayes-Sane Sard cubature method," NeurIPS 31, 2018.
[22] M. Osborne, R. Garnett, Roberts, C. Hart, S. Aigrain, and Gibson, "Bayesian quadrature
for ratios," PMIR, vol. 22, pp. 332-840, 2012.
[23] X. Xi, F.-X. Briol, and M. Girolami, "Bayesian quadrature for multiple related integrals," PMLR, vol. 80, pp. 5373-5382, 2018.
[24] F-X. Briol, C. J. Oates, M. Girolami, and M. A. Osborne, "Frank-War Wolfe Bayesian quadrature: probabilistic integration with theoretical guarantees," NeurlPS 28, 2015.
[25] Gunter, M. A. Osborne, R. Garnett, P. Hennig, and Roberts, "Sampling for inference in probabilistic models with fast Bayesian quadrature," " NeurIPS 27, 2014.

---

[26] M. Kennedy, "Bayesian quadrature with non non-normal approximating functions," Stat. Comput., vol. 8, no. 4, pp. 365-375, 1998.
[27] A. O'Hagan, "Bayes Hermite quadrature." J. Stat. Plan. Infer, vol. 29, no. 3, pp. 245-260,
1991.
[28] F. Larkin, "Gaussian measure in Hilbert space and applications in numerical analysis." Rocky Mountain Math., vol. 2, no. 3, pp. 379-422, 1972.
[29] C. E. Rasmussen and Z. Ghahramani. "Bayesian Monte Carlo," NeurIPS 15, pp. 505-512, 2003.
[30] F.-X. Briol, C. J. Oates, M. Girolami, M. Osbore,, D. Sejdinovic, et al., "Probabilistic integration: a role in statistical computation??" Stat. Sci., vol. 34, no. pp- 1 1-22, 2019.
[31] J. Močkus, "On Bayesian methods for seeking the extremum and their application.." " IFIP
Congress, pp. 195-200, 1977.
[32] Močkus, V. Tiesis, and Zilinskas, "The application of Bayesian methods for seeking the extremum," Towards Global Optimization, vol. 2, no. 2, pp. 117-129, 1978.
[33] J. Močkus, Bayesian Approach to Global Optinization. Springer, 1989.
[34] J. Snoek, H. Larochelle, and R. P. Adams, "Practical Bayesian optimization of machine learning algorithms," NeurlPS 25, 2012.
[35] Hennig and M. Kiefel, QQuasi-Newton methods: new direction," JMLR, vol. 14, no. 1, 843-865. 2013.
[36] M. Mahsereci and P. Hennig, "Probabilistic line searches for stochastic optimization," " NeurlPS
28, 2015.
[37] Skilling, "Bayesian solution of ordinary differential equations." Maximum Entropy and Bayesian Methods, pp. 23-37, 1992.
[38] P. R. Conrad, M. Girolami, S. Särkkä, A. Stuart, and K. Zygalakis, "Statistical analysis of differential equations: Introducing probability measures on numerical solutions," Stat. Comput., vol. 27, no. 4, pp. 1065-1082, 2016.
[39] O. A. Chkrebtii, D. A. Campbell, B. Calderhead, and M. A. Girolami, "Bayesian solution uncertainty quantification for differential equations," Bayesian Anal., vol. 11, no. 4, pp. 1239- 1 1267, 2016.
[40] M. Schober, S. Särkkä, and P. Hennig, "A probabilistic model for the numerical solution of initial value problems," Stat. Comput., vol. 29, no. 1, pp. 99-122, 2019.
[41] O. Teymur, K. Zygalakis, and B. Calderhead, "Probabilistic linear multistep methods," " NeurlPS 29, 2016.
[42] O. Teymur, B. Calderhead, H. C. Lie, and T. J. Sullivan, "Implicit probabilistic integrators for ODEs," NeurlPS 31, 2018.
[43] Hennig and S. Hauberg. "Probabilistic solutions to differential equations and their application to Riemannian statistics," PMLR, vol. 33, pp. 347-355, 2013.
[44] H. Kersting and P. Hennig, "Active uncertainty calibration in Bayesian ODE solvers," Proc. 32nd Conf. Uncertainty Artificial Intelligence, pp. 309-318, 2016.
[45] M. Schober, D. Duvenaud, and P. Hennig, "Probabilistic ODE solvers with Runge- Kutta means." NeurlPS 27, 2014.
[46] H. Owhadi, "Bayesian numerical homogenization." Multiscale Model. Sim., vol. 13, no. 3, pp. 812-828, 2015.
[47] Tronarp, H. Kersting, S. Särkkä, and Hennig, "Probabilistic solutions to ordinary differential
equations as nonlinear Bayesian filtering: new perspective," Stat. Comput., vol. 29, no. 6, pp. 1297-1315, 2019.
[48] Wang. J. Cockayne, and C. Oates, "On the Bayesian solution of differential equations," arXiv:1805.077 1805.1800009 2018.
[49] J. Cockayne, C. Oates, T. Sullivan, and M. Girolami, "Probabilistic numerical methods for PDE- constrained Bayesian inverse problems," AIP Conference Proceedings, vol. 1853, p. 060001, 2017.
12

---

[50] O. A. Chkrebtii and D. A. Campbell, "Adaptive step-size selection for state- space probabilistic differential equation solvers," Stat. Comput, vol. 29, no. 6, pp. 1285-1295, 2019.
[51] H. Owhadi and C. Scovel, Operator-Adapte Adapted Wavelets, Fast Solvers, and Numerical Homoge- nization: From a Game Theoretic Approach to Numerical Approximation and Algorithm Design. Cambridge University Press, 2019.
[52] A. Abdulle and G. Garegnani. "Random time step probabilistic methods for uncertainty quan- tification in chaotic and geometric numerical integration," Stat. Comput., vol. 30, no. 4, pp. 907-
932, 2020.
[53] H. Kersting, Sullivan, and Hennig, "Convergence rates of Gaussian ODE filters," Stat. Comput., vol. 30, 6, pp. 1791-1816, 2020.
[54] J. Wang, J. Cockayne, O. Chkrebtii, T. Sullivan, and C. J. Oates, "Bayesian numerical methods for nonlinear partial differential equations," arXiv:2 /2587, 2021.
[55] N. Bosch, P. Hennig, and F. Tronarp, "Calibrated Adaptive Probabilistic ODE Solvers," PMLR, vol. 130, pp. 3466-3474, 2021.
[56] L. F. Richardson and J. Gaunt, "VIII. The deferred approach the limit," Philos. Trans.
Soc. A, vol. 226, no. 336-646, pp. 299-361, 1927.
[57] C. Rasmussen and C. Williams, Gaussian Processes for Machine Learning. MIT Press, 2006.
[58] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and Flannery, Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge University Press, 2007.
[59] E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems. Springer, 2008
[60] C. Runge, "Über empirische Funktionen und die Interpolation zwischen äquidistanten Ordi- naten," Zeitschrift für Mathematik und Physik, vol. 46, no. 224-243, p. 20, 1901.
[61] R. Bulirsch and J. Stoer, 'Fehlerabschääzunnge und Extrapolation mit rationalen Funktionen bei Verfahren vom Richardson- Typus," Numer. Math., vol. 6, no. 1, pp. 413 413-427, 1964.
[62] D. C. Joyce. "Survey of extrapolation processes in numerical analysis." SIAM Rev. vol. 13,
no. 4, pp. 435-490, 1971.
[63] T. A. Oliver, N. Malaya, R. Ulerich, and R. D. Moser, "Estimating uncertainties in statistics computed from direct numerical simulation," Phys. Fluids, vol. 26, no. 3, p. 035101, 2014.
[64] B. Chartres and R. Stepleman, "A general theory of convergence for numerical methods," SIAM
J. Numer: Anal., vol. no. 3. pp. 476-492, 1972.
[65] J. C. de los Reyes, Numerical PDE- -DE-Constrai Optimization. Springer, 2015.
[66] M. L. Stein, Interpolation of Spatial Data: Some Theory for Kriging. Springer, 2012.
[67] C. Burg and Erwin, "Application of Richardson extrapolation to the numerical solution of partial differential equations," Numerical Methods for PDEs, vol. 25, no. 4, pp. 810-832, 2009.
[68] J. Cockayne, M. M. Graham, J. Oates, T. J. Sullivan, and O. Teymur, "Testing whether learning procedure calibrated," arXiv:2012.12277 2021.
[69] T. G. Kolda and J. R. Mayo, "Shifted power method for computing tensor eigenpairs," SIAM J.
Matrix Anal. A., vol. 32, no. 4, pp. 1095-1124, 2011.
[70] B. W. Bader, T. Kolda, al., "Tensor toolbox for MATLAB, version 3.2.1." 2021.
[71] Y. Kuramoto, "Diffusion-induced chaos in reaction systems," Prog. Theor. Phys. Supp., vol. 64,
pp. 346-367, 1978.
[72] G. Sivashinsky, "Nonlinear analysis of hydrodynamic instability in laminar flames," Acta Astronaut., vol. 4, no. 11, pp. 1177-1206, 1977.
[73] A. Kassam and L. Trefethen, "Fourth-order time stepping for stiff PDEs," SIAM J. Sci. Comput., vol. 26, pp. 2214-1233. 2005.
[74] T. Karvonen, G. Wynne, F. Tronarp. C. Oates, and S. Särkkä, "Maximum likelihood estimation and uncertainty quantification for Gaussian process approximation of deterministic functions," SIAM/ASA J. Uncertainty Quantification, vol. 8, no. 3, pp. 926-958, 2020.
[75] H. Kersting and M. Mahsereci, "A Fourier state space model for Bayesian ODE filters." arXiv:2007.09118, 2020.
13


Reviews:
[{"review_id": "m8gVhD1HV1i", "reviewer": "Reviewer_fzBQ", "summary": "This paper proposes a new \"black box\" approach to probabilistic numerics (PN)\nthat leverages existing (non-probabilistic) numerical methods to estimate,\nprobabilistically, the true, intractable quantity of interest (e.g. the solution\nof an ODE or of an eigenvalue problem). The method is demonstrated in three\nexperiments and is (i) compared to other PN approaches for solving ODEs, (ii)\napplied to eigenvalue problems (for which no PN method has yet been developed),\nand (iii) used for uncertainty quantification on a challenging nonlinear PDE\nproblem. Due to its generality, the proposed method expands probabilistic\nnumerical approaches to a larger range of problems.\n", "questions": "", "limitations": "", "rating": 7, "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": "__Originality__\nThe presented BBPN method presents, to the best of my knowledge, a completely\nnovel approach to the field of PN.\nThe relation to the well-established Richardson extrapolation is made very\nclear, and relevant PN approaches for the considered problems seem adequately\ncited.\n\n\n__Quality__\nThe submission is technically sound and the paper provides all necessary\nderivations.\nThe experimental results cover an appropriate range of experiments, though a\nwider range of comparisons to a wider range of PN methods that go beyond ODE\nsolvers, e.g. for quadrature or solving linear systems, would provide\nadditional, useful evidence to best assess the empirical qualities of the\nproposed approach.\nSimilarly, the flexibility could be further demonstrated by more experiments\nthat do not have such a clear notion of \"step size\" as in the ODE and PDE\nexample.\nBut again, while such additions would be very welcome I want to emphasize that\nthe current evaluation is already very useful and covers a good range of\nexperiments.\n\n\n__Clarity__\nThe presentation is very clear. The paper provides relevant background on the\nRichardson extrapolation and on PN and presents the method in an understandable\nand concise manner. The appendix provides many helpful details.\n\n\n__Significance__\nI believe that the paper provides a significant contribution to the field of\nprobabilistic numerics. Due to the universal applicability of the presented\nmethod, it potentially provides a new baseline for any future PN method that is\ndeveloped for novel problem classes (e.g. for eigenvalue problems as mentioned\nin the paper).\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations are discussed in Section 5 (Discussion) and Section 1\n(Introduction): \n- BBPN has possibly increased computational cost (due to the multiple calls to\n  the numerical method that BBPN requires); \n- BBPN requires a statistical model of the numerical error;\n- BBPN only makes use of the order of the method - in comparison, other non-BB\n  PN methods could potentially benefit from considering additional information\n  about the numerical method. \n\nIn addition, the important dependence on the parameterisation of h is discussed\nin more detail in the supplementary C.2.2. An additional experiment where h is\nnot equal to a step size, e.g. for ODE solvers with adaptive step-size\nselection, might provide additional relevant experimental evidence.  \nIn a similar light, comparisons to other PN methods on problems without such a\nclear notion of step-size, such as quadrature or linear systems, could provide\nfurther supporting evidence.\n\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "U0KAvLfJLoj", "reviewer": "Reviewer_7xVs", "summary": "This paper introduces a black-box approach to construct a probabilistic numerical method by treating a sequence of increasingly accurate approximations generated by a classic numerical method as observations. This probabilistic analogue of Richardson's deferred approach to the limit allows higher orders of convergence to be achieved and extends the PN approach to previously unexplored tasks.", "questions": "", "limitations": "", "rating": 7, "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": "\nThe idea of extending Richardson's deferred approach to the limit (RDAL) by conditioning a stochastic process on the output of a sequence of numerical methods seems very natural once one accepts the central premise of PN that epistemic uncertainty is to be quantified. I therefore find the approach presented in this paper to be quite attractive. Overall the paper is well-executed both in terms of its coherence, mathematical presentation, as well as experimental evaluation. There are some inherent limitations that come with the choice of method in my opinion and some important questions that are left somewhat unexplored. However, overall the positives by far outweigh the limitations, which is why I would recommend acceptance. See below for details.\n\n## Novelty\nFrom a classical perspective the added functionality of a certain class of interpolating functions given by the posterior mean of the GP that satisfy Proposition 1, may appear somewhat limited, since they directly extend the existing idea behind RDAL. However, from a PN perspective the added functionality of uncertainty quantification for the quantity of interest is useful. In particular probabilistic numerical methods for linear algebra have been almost exclusively focussed on the solution of linear systems, while this approach naturally extends to the estimation of eigenvalues (Section 4.2). Also inferring the order of a numerical method via hyperparameter optimization is an interesting and novel functionality to the best of my knowledge, albeit in practice possibly of limited importance (see limitations below).\n\n## Theoretical Contribution\nWhile the underlying result of obtaining a higher-order of convergence is clearly not new, the authors are very upfront about this. The kernel construction which preserves this property however is elegant and presents a significant contribution.\n\n## Experiments\nThe experimental evaluation of the presented BBPN approach is well-presented and also includes calibration. Often in practice PN methods suffer from miscalibration due to the problem of non-linear dependence of the quantity of interest on the data. The experimental demonstrations of how BBPN applies to different problems in numerics is commendable. However, clearly these experiments are illustrative in nature and may mask some of the limitations of the approach, in particular for high-dimensional quantities of interest (i.e. eigenvalue estimation for matrices of n >> 100) or expensive numerical methods, where computational budget does not allow running the method more than a few times.\n\n## Clarity and Reproducibility\nThe paper is generally well-written with a clear line of exposition, supported by explanatory figures. The paper comes with code reproducing some of the figures and has sufficiently detailed experimental descriptions in the appendix.\n\n### Typos\n- l30: computational resources", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "## Limitations\nTo my mind, the presented approach has the following limitations, which lead to some questions, which as a reader I felt were not entirely answered:\n- As the authors state themselves BBPN requires an underlying numerical method to be run to convergence for a sequence of parameters h. This makes BBPN (and RDAL) roughly a constant factor more expensive than simply running a classical method. When is this tradeoff worth it for the gain in accuracy? Also, likely there are settings (e.g. linear dependence on quantity of interest on data), where using PN methods which condition directly on the observations of the problem are much cheaper. I realize that wall-clock times are difficult to compare across implementations. However, running the same numerical method multiple times plus the added cost of GP inference in comparison to the cost of a single run to the highest accuracy is valid to my mind and would give some insight here.\n- Performing regression on a dataset with few observations (outputs of numerical methods) is a setting where Gaussian processes are well-established. However, in high-dimensional settings and when performing extrapolation, their performance may be quite poor. To my understanding this corresponds to the case where one cannot run the baseline numerical method to high-precision or the quantity of interest is high-dimensional (e.g. eigenvalue problems for realistically sized matrices n >> 100). Does this limit BBPN to problems where the QoI is low-dimensional?\n- In the same vein, inferring the order of the numerical method (\\alpha), while an attractive feature, may be computationally quite costly relative the benefit it provides. Also from the paper i was not able to judge how well the convergence order was inferred, since for the presented eigenvalue estimation problem the true convergence order is not known.\n\n## Societal Impact\nNo direct negative impact due to the methodological nature of this paper.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "Tk_uj82hExU", "reviewer": "Reviewer_TEPo", "summary": "The authors present a probabilistic numeric algorithm (BBPN) which operates analogously to Richardson's deferred approach to the limit (RDAL) in traditional numerical analysis, and has an advantage over other PN approaches in that it is applicable in any setting where an existing numerical method is available. The authors provide a detailed description and analysis of BBPN, including theoretical guarantees on the order of convergence. They finish with an experimental assessment of BBPN, by comparison to both probabilistic and traditional numerical methods.", "questions": "", "limitations": "", "rating": 7, "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": "This paper is a valuable, well-presented and very well-written contribution to the field of probabilistic numerics. The paper is well structured, striking a good balance between background, motivation, theory and experiment. BBPN is well motivated and clearly described, and the experimental results are convincing. I especially appreciated the presence of the uncertainty quantification metric in the experimental comparison against other probabilistic numeric methods for ODEs. I would welcome more discussion of the uncertainty quantification properties of BBPN earlier in the paper, since uncertainty quantification is of critical importance for probabilistic numeric methods. I understand that there may be no good way to expand this beyond what is already present in the experimental section, given the tight nine-page constraint.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I do not see any potential negative societal impact, and the authors adequately address the limitations of their work.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3.5"}, "extra_text": {}}, {"review_id": "8hLApg1zz25", "reviewer": "Reviewer_reNR", "summary": "This paper presents a probabilistic version of Richardson's \"deferred approach to the limit\" (RDAL) method.  RDAL is an established numerical analysis method (developed in 1927) with many applications. In RDAL several low-accuracy approximations are extrapolated to produce a more accurate approximation (that is costly to come with directly). \nIn the framework proposed in this draft, namely \"Black box probabilistic numerics\" (BBPN), the aforementioned extrapolation/curve fitting is carried out via Gaussian Processes and therefore comes with all the benefits for Bayesian approaches such as quantification of the  estimation uncertainty. This approach is simple and it can potentially be combined with Bayesian methods other than GPs as well. \nIn the experiment section the following applications are presented: nonlinear ODE/PDE i.e. ordinary/partial differential equations and eigenvalue problem.   \n", "questions": "", "limitations": "", "rating": 7, "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 main idea of this paper boils down to probabilistic curve fitting via Gaussian Processes. This idea is simple, up to my knowledge original and potentially significant.  Nevertheless, I think the contribution and impact  of this paper could be enriched if the authors would conduct more thorough experiments that quantitatively compare BBPN against state-of-the-art. In the current draft, only in the first experiment such a quantitatively comparison is provided.\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes,they mention that the limitations of their proposed algorithm are: increased computational costs as well as additional assumptions required for modelling uncertainty. \nSince this work is theoretical it does not have a direct societal impact.\n\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper suggests a "black box" approach to probabilistic numerics. The idea is to treat the output of existing (non-probabilistic) numerical methods as observations, after which the true/unknown quantities can be treated via probabilistic inference. Examples are given applying this idea to ODEs, to eigenvalue problems, and to uncertainty with nonlinear PDEs. All reviewers felt the paper was novel, technically sound, useful, and clearly written.
Decision: Accept (Poster)