Scientific field at the intersection of statistics, machine learning and applied mathematics
Probabilistic numerics is an active field of study at the intersection of applied mathematics, statistics, and machine learning centering on the concept of uncertainty in computation. In probabilistic numerics, tasks in numerical analysis such as finding numerical solutions for integration, linear algebra, optimization and simulation 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 (often, but not always, Bayesian 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:
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.
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 formula_1 at a number of points formula_2 are used to estimate the integral formula_3 of a function formula_4 against some measure formula_5. Bayesian quadrature consists of specifying a prior distribution over formula_6 and conditioning this prior on formula_1 to obtain a posterior distribution over formula_6, then computing the implied posterior distribution on formula_3. 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 formula_4 is expensive to evaluate and the dimension of the data is small to moderate.
Optimization.
Probabilistic numerics have also been studied for mathematical optimization, which consist of finding the minimum or maximum of some objective function formula_4 given (possibly noisy or indirect) evaluations of that function at a set of points.
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 formula_4 throughout the optimization procedure; this often takes the form of a Gaussian process prior conditioned on observations. This belief then guides the algorithm in obtaining observations that are likely to advance the optimization process. Bayesian optimization policies are usually realized by transforming the objective function posterior into an inexpensive, differentiable "acquisition function" that is maximized to select each successive observation location. 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. A welcome side effect from this approach is that uncertainty in the objective function, as measured by the underlying probabilistic belief, can guide an optimization policy in addressing the classic exploration vs. exploitation tradeoff.
Local optimization.
Probabilistic numerical methods have been developed in the context of stochastic optimization for deep learning, in particular to address main issues such as 
learning rate tuning and line searches, 
batch-size selection, early stopping, 
pruning, and first- and second-order search directions.
In this setting, the optimization objective is often an empirical risk of the form formula_13 defined by a dataset formula_14, and a loss formula_15 that quantifies how well a predictive model formula_16 parameterized by formula_17 performs on predicting the target formula_18 from its corresponding input formula_19 .
Epistemic uncertainty arises when the dataset size formula_20 is large and cannot be processed at once meaning that local quantities (given some formula_21) such as the loss function formula_22 itself or its gradient formula_23 cannot be computed in reasonable time.
Hence, generally mini-batching is used to construct estimators of these quantities on a random subset of the data. Probabilistic numerical methods model this uncertainty explicitly and allow for automated decisions and parameter tuning.
Linear algebra.
Probabilistic numerical methods for linear algebra
 have primarily focused on solving systems of linear equations of the form formula_24 and the computation of determinants formula_25.
A large class of methods are iterative in nature and collect information about the linear system to be solved via repeated matrix-vector multiplication formula_26 with the system matrix formula_27 with different vectors formula_28. 
Such methods can be roughly split into a solution- and a matrix-based perspective, depending on whether belief is expressed over the solution formula_29 of the linear system or the (pseudo-)inverse of the matrix formula_30. 
The belief update uses that the inferred object is linked to matrix multiplications formula_31 or formula_32 via formula_33 and formula_34.
Methods typically assume a Gaussian distribution, due to its closedness under linear observations of the problem. While conceptually different, these two views are computationally equivalent and inherently connected via the right-hand-side through formula_35.
Probabilistic numerical linear algebra routines have been successfully applied to scale Gaussian processes to large datasets. In particular, they enable "exact" propagation of the approximation error to a combined Gaussian process posterior, which quantifies the uncertainty arising from both the "finite number of data" observed and the "finite amount of computation" expended.
Ordinary differential equations.
Probabilistic numerical methods for ordinary differential equations formula_36, have been developed for initial and boundary value problems. Many different probabilistic numerical methods designed for ordinary differential equations have been proposed, and these can broadly be grouped into the two following categories:
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, latent force models, and to differential equations with a geometric structure such as symplecticity.
Partial differential equations.
A number of probabilistic numerical methods have also been proposed for partial differential equations. As with ordinary differential equations, the approaches can broadly be divided into those based on randomisation, generally of some underlying finite-element mesh and those based on Gaussian process regression.
Probabilistic numerical PDE solvers based on Gaussian process regression recover classical methods on linear PDEs for certain priors, in particular methods of mean weighted residuals, which include Galerkin methods, finite element methods, as well as spectral methods.
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 formula_44, expressed as a formal power series with random coefficients, and asked for "probable values" of formula_45 given this prior and formula_46 observations formula_47 for formula_48.
A later seminal contribution to the interplay of numerical analysis and probability was provided by Albert Suldin in the context of univariate quadrature. The statistical problem considered by Suldin was the approximation of the definite integral formula_49 of a function formula_50, under a Brownian motion prior on formula_51, given access to pointwise evaluation of formula_51 at nodes formula_53. Suldin 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.
Suldin's point of view was later extended by Mike Larkin.
Note that Suldin's Brownian motion prior on the integrand formula_51 is a Gaussian measure and that the operations of integration and of point wise evaluation of formula_51 are both linear maps.
Thus, the definite integral formula_49 is a real-valued Gaussian random variable.
In particular, after conditioning on the observed pointwise values of formula_51, it follows a normal distribution with mean equal to the trapezoidal rule and variance equal to formula_58.
This viewpoint is very close to that of Bayesian quadrature, seeing the output of a quadrature method not just as a point estimate but as a probability distribution in its own right.
As noted by Houman Owhadi and collaborators, 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 Packel observed, 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 leads to 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.
References.
 ;templatestyles src="Reflist/styles.css" / ;