Probabilistic numerics

Probabilistic numerics is a scientific field at the intersection of statistics, machine learning and applied mathematics, where tasks in numerical analysis including finding numerical solutions for integration, linear algebra, optimisation and differential equations are seen as problems of statistical, probabilistic, or Bayesian inference.

Introduction

A numerical method is an algorithm that approximates the solution to a mathematical problem (examples below include the solution to a linear system of equations, the value of an integral, the solution of a differential equation, the minimum of a multivariate function). In a probabilistic numerical algorithm, this process of approximation is thought of as a problem of estimation, inference or learning and realised in the framework of probabilistic inference.

Formally, this means casting the setup of the computational problem in terms of a prior distribution, formulating the relationship between numbers computed by the computer (e.g. matrix-vector multiplications in linear algebra, gradients in optimization, values of the integrand or the vector field defining a differential equation) and the quantity in question (the solution of the linear problem, the minimum, the integral, the solution curve) in a likelihood function, and returning a posterior distribution as the output. In most cases, numerical algorithms also take internal adaptive decisions about which numbers to compute, which form an active learning problem.

Many of the most popular classic numerical algorithms can be re-interpreted in the probabilistic framework. This includes the method of conjugate gradients, Nordsieck methods, Gaussian quadrature rules, and Quasi-Newton methods. In all these cases, the classic method is based on a regularized least-squares estimate that can be associated with the posterior mean arising from a Gaussian prior and likelihood. In such cases, the variance of the Gaussian posterior is then associated with a worst-case estimate for the squared error.

Probabilistic numerical methods promise several conceptual advantages over classic, point-estimate based approximation techniques:

• They return structured error estimates (in particular, the ability to return joint posterior samples, i.e. multiple realistic hypotheses for the true unknown solution of the problem)
• Hierarchical Bayesian inference can be used to set and control internal hyperparameters in such methods in a generic fashion, rather than having to re-invent novel methods for each parameter
• Since they use and allow for an explicit likelihood describing the relationship between computed numbers and target quantity, probabilistic numerical methods can use the results of even highly imprecise, biased and stochastic computations. Conversely, probabilistic numerical methods can also provide a likelihood in computations often considered "likelihood-free" elsewhere
• Because all probabilistic numerical methods use essentially the same data type – probability measures – to quantify uncertainty over both inputs and outputs they can be chained together to propagate uncertainty across large-scale, composite computations
• Sources from multiple sources of information (e.g. algebraic, mechanistic knowledge about the form of a differential equation, and observations of the trajectory of the system collected in the physical world) can be combined naturally and inside the inner loop of the algorithm, removing otherwise necessary nested loops in computation, e.g. in inverse problems.

These advantages are essentially the equivalent of similar functional advantages that Bayesian methods enjoy over point-estimates in machine learning, applied or transferred to the computational domain.

Numerical Tasks

Integration

Main article: Bayesian quadrature

Bayesian quadrature with a Gaussian process conditioned on $n=0,3,\ {\text{and}}\ 8$ evaluations of the integrand (shown in black). Shaded areas in the left column illustrate the marginal variance. The right figure shows the prior ($n=0$) and posterior ($n=3,8$) Gaussian distribution over the value of the integral, as well as the true solution.

Probabilistic numerical methods have been developed for the problem of numerical integration, with the most popular method called Bayesian quadrature .

In numerical integration, function evaluations $f(x_{1}),\ldots ,f(x_{n})$ at a number of points $x_{1},\ldots ,x_{n}$ are used to estimate the integral $\textstyle \int f(x)\nu (dx)$ of a function $f$ against some measure $\nu$. Bayesian quadrature consists of specifying a prior distribution over $f$ and conditioning this prior on $f(x_{1}),\ldots ,f(x_{n})$ to obtain a posterior distribution over $f$, then computing the implied posterior distribution on $\textstyle \int f(x)\nu (dx)$. The most common choice of prior is a Gaussian process as this allows us to obtain a closed-form posterior distribution on the integral which is a univariate Gaussian distribution. Bayesian quadrature is particularly useful when the function $f$ is expensive to evaluate and the dimension of the data is small to moderate.

Optimization

Probabilistic approaches have also been studied for mathematical optimization. Perhaps the most notable effort in this direction is Bayesian optimization, a general approach to optimization grounded in Bayesian inference. Bayesian optimization algorithms operate by maintaining a probabilistic belief about an objective function throughout optimization; this often takes the form of a Gaussian process prior conditioned on observations obtained during optimization. This belief then guides the algorithm in obtaining observations that are likely to advance the optimization process. One prominent approach is to model optimization via Bayesian sequential experimental design, seeking to obtain a sequence of observations yielding the most optimization progress as evaluated by an appropriate utility function.

Linear Algebra

Probabilistic numerical methods for linear algebra have primarily focused on the solution of linear systems of the form $Ax=b$.

,

Ordinary Differential Equations

Probabilistic numerical methods for ordinary differential equations ${\dot {f}}(t)=w(t,f(t))$, have been developed for initial and boundary value problems. Many different probabilistic numerical methods designed for ordinary differential equations have been proposed. These can broadly be grouped into the two following categories:

• Randomisation-based methods: These algorithms are based on perturbing standard deterministic numerical methods in a random fashion. For example, this has been achieved by adding Gaussian perturbations on the solution of one-step integrators or by perturbing randomly their time-step. This defines a probability model for the solution of the differential equation.
• Gaussian process regression. These algorithms are based on posing the problem as a Gaussian process regression problem, in a similar fashion to Bayesian cubature - But with different and often non-linear observation models. In its infancy, this was based on naive Gaussian process regression. This was later abandoned in favour of Gaussian Markov models, $\mathrm {d} x(t)=Ax(t)\mathrm {d} t+B\mathrm {d} v(t),$ where $x(t)$ is a vector containing the, say, first $\nu$ derivatives of $y(t)$. Inference can thus be implemented efficiently with Kalman filtering based methods.

The boundary between these two categories is not sharp, indeed a Gaussian process regression approach based on randomised data was developed as well. These methods have been applied to problems in computational Riemannian geometry, inverse problems, and latent force models.

Partial Differential Equations

,

History and Related Fields

The interplay between numerical analysis and probability is touched upon by a number of other areas of mathematics, including average-case analysis of numerical methods, information-based complexity, game theory, and statistical decision theory. Precursors to what is now being called “probabilistic numerics” can be found as early as the late 19th and early 20th century.

The origins of probabilistic numerics can be traced to a discussion of probabilistic approaches to polynomial interpolation by Henri Poincaré in his Calcul des Probabilités: in modern terminology, Poincaré considered a Gaussian prior distribution on a function $f\colon \mathbb {R} \to \mathbb {R}$, expressed as a formal power series with random coefficients, and asked for “probable values” of $f(x)$ given this prior and $n\in \mathbb {N}$ observations $f(a_{i})=B_{i}$ for $i=1,\dots ,n$.

A later seminal contribution to the interplay of numerical analysis and probability was provided by Al′bert Sul′din in the context of univariate quadrature. The statistical problem considered by Sul′din was the approximation of the definite integral $\textstyle \int _{a}^{b}u(t)\,\mathrm {d} t$ of a function $u\colon [a,b]\to \mathbb {R}$, under a Brownian motion prior on $u$, given access to pointwise evaluation of $u$ at nodes $t_{1},\dots ,t_{n}\in [a,b]$. Sul′din showed that, for given quadrature nodes, the quadrature rule with minimal mean squared error is the trapezoidal rule; furthermore, this minimal error is proportional to the sum of cubes of the inter-node spacings. As a result, one can see the trapezoidal rule with equally-spaced nodes as statistically optimal in some sense — an early example of the average-case analysis of a numerical method.

As discussed in and interplays between numerical approximation  and statistical inference can also be traced back to Palasti and Renyi, Sard,  Kimeldorf and Wahba (on the correspondence between Bayesian estimation and spline smoothing/interpolation) and Larkin (on the correspondence between Gaussian process regression and numerical approximation). Although the approach of modelling a perfectly known function as a sample from a random process may seem counterintuitive, a natural framework for understanding it can be found in information-based complexity (IBC), the branch of computational complexity founded on the observation that numerical implementation requires computation with partial information and limited resources. In IBC, the performance of an algorithm operating on incomplete information can be analyzed  in the  worst case  or the average case (randomized) setting  with respect to the missing information. Moreover,  as observed by Packel, the average case setting could be interpreted as a mixed strategy in an adversarial game obtained by lifting a (worst case) minmax problem  to a minmax problem over mixed (randomized) strategies. This observation   initiates a natural connection between numerical approximation  and  Wald's decision theory, evidently influenced by von Neumann's  theory of games. To describe this connection consider the optimal recovery setting of Micchelli and Rivlin in which one tries to approximate an unknown function from a finite number of linear measurements on that function. Interpreting this optimal recovery problem as a zero-sum game where Player I selects the unknown function and Player II selects its approximation, and using relative errors in a quadratic norm to define losses, Gaussian priors emerge as optimal mixed strategies for such games, and the covariance operator of the optimal Gaussian prior is determined by the quadratic norm used to define the relative error of the recovery.

Software

• ProbNum: Probabilistic Numerics in Python.
• ProbNumDiffEq.jl: Probabilistic numerical ODE solvers based on filtering implemented in Julia.

See also

• Average-case analysis
• Information-based complexity