Optimality Criteria for Probabilistic Numerical Methods

TL;DR

This paper explores an optimality criterion from Bayesian experimental design for probabilistic numerical methods, highlighting differences from average-case approaches and proposing regimes for developing optimal methods.

Contribution

It introduces a novel optimality criterion for probabilistic numerical methods based on Bayesian experimental design, contrasting it with traditional average-case analysis.

Findings

Optimal information differs from average-case optimal information.

Different regimes for developing optimal probabilistic numerical methods.

Connection established between Bayesian experimental design and numerical analysis.

Abstract

It is well understood that Bayesian decision theory and average case analysis are essentially identical. However, if one is interested in performing uncertainty quantification for a numerical task, it can be argued that standard approaches from the decision-theoretic framework are neither appropriate nor sufficient. Instead, we consider a particular optimality criterion from Bayesian experimental design and study its implied optimal information in the numerical context. This information is demonstrated to differ, in general, from the information that would be used in an average-case-optimal numerical method. The explicit connection to Bayesian experimental design suggests several distinct regimes in which optimal probabilistic numerical methods can be developed.

Equations

$${\text{ACE}}_p(e,d_e)$$

$$d_e^{\ast }$$

$$e^{\ast }$$

$${\pi }_{Y|e}(S)=\iint 1_S(y)\mathrm{d}{\pi }_{Y|x,e}(y)\mathrm{d}{\pi }_X(x),$$

$$\text{BR}(e,d_e)$$

$$d_e^{\ast }$$

$$e^{\ast }$$

$$d_e^{\ast }=\frac{1}{2}\sum _{i=1}^n(x(t_{i-1})+x(t_i))(t_i-t_{i-1}),$$

$$\int \ell (x,d_e^{\prime }(y))\mathrm{d}{\pi }_{Y|x,e}(y)$$

$${\mathcal{A}}_e^{\ast }(y_e)$$

$${\left[\frac{\mathrm{d}\varphi }{\mathrm{d}a}\right]}_{i,j}$$

$$\varphi (a)$$

$$\text{BED}(e)$$

$${\mathcal{E}}_{\text{BED}}^{\ast }$$

$$u(e,y)$$

$$\text{BED}(e)$$

$$=\iint \ell (x^{\prime },d_e^{\ast }(y))\mathrm{d}{\pi }_{X|y,e}(x^{\prime })\mathrm{d}{\pi }_{Y|e}(y)$$

$$=\iint \ell (x^{\prime },d_e^{\ast }(y))\mathrm{d}{\pi }_{Y|x^{\prime },e}(y)\mathrm{d}{\pi }_X(x^{\prime })=\text{BR}(e,d_e^{\ast })$$

$$\text{BR}(e,d_e^{\ast })$$

$$u(e,y)$$

$$\ell (x,a)$$

$$u(e,y)$$

$$=\int \log \left(\frac{\mathrm{d}{\pi }_{X|y,e}}{\mathrm{d}{\pi }_X}(x)\right)\mathrm{d}{\pi }_{X|y,e}(x)=:D_{\text{KL}}({\pi }_{X|y,e}\Vert {\pi }_X).$$

$$d_e(y)=N\left(\sum _{i=1}^n\frac{(x(t_{i-1})+x(t_i))}{2}(t_i-t_{i-1}),\sum _{i=1}^n\frac{(t_i-t_{i-1}{)}^3}{12}\right).$$

$$\underset{{\pi }_X}{\min }\underset{{\tilde{\pi }}_X}{\max }\iint \ell (x,d_e^{\ast }(y))\mathrm{d}{\pi }_{Y|x,e}(y)\mathrm{d}{\tilde{\pi }}_X(x).$$

$$u(e,x,y)$$

$$\text{BPN}(e)$$

$$\int \ell (x,x^{\prime })\mathrm{d}{\pi }_{X|y,e}(x^{\prime })$$

$$\text{BPN}(e)$$

$${\mathcal{E}}_{\text{BPN}}^{\ast }$$

$$-\Delta x(t)$$

$$x(t)$$

Full text

1 Introduction

To fix notation, consider the task of approximating a quantity of interest $\phi\colon\mathcal{X}\rightarrow\Phi$ based on finite information provided through a map $y_{e}\colon\mathcal{X}\rightarrow\mathcal{Y}_{e}$. There may be several such maps with equivalent computational cost, and these are indexed by $e\in\mathcal{E}$. This abstract formulation covers most basic numerical tasks, including numerical integration, numerical optimisation and the numerical solution of a differential equation [64]. Thus a numerical method is considered as a map $d_{e}\colon\mathcal{Y}_{e}\rightarrow\Phi$ and the set of all possible numerical methods is denoted $\mathcal{D}_{e}$.

The performance of a numerical method can be quantified in several ways, but a common scenario is for $\mathcal{X}$ and $\Phi$ to be normed spaces, in which case the average-case error of a numerical method $d_{e}$ can be defined as

[TABLE]

where $p\in[1,\infty)$ and $\pi_{X}$ is a Borel probability distribution on $\mathcal{X}$, to be specified. Average-case analysis can be motivated by considering the state $x$ to be the realisation of a random variable $X\sim\pi_{X}$ and then assessing the average performance of a numerical method [51].

In information-based complexity there is interest in optimal numerical methods, their associated optimal information and the tractability of the numerical task itself. For instance, an average-case-optimal numerical method is defined as

[TABLE]

and average-case-optimal information is defined as

[TABLE]

See [51, 57, 58, 64]. Here and in the sequel we simply assume all minima are realised, so (1) and (2) are presented as $\operatorname*{arg\,min}$ rather than $\operatorname*{arg\,inf}$. In [31] and later [21] it was observed that an average-case-optimal numerical method is mathematically identical to a Bayes rule in Bayesian decision theory when $\pi_{X}$ is interpreted as the prior distribution (see Section 2).

There has been recent interest in the development of probabilistic numerical methods, which treat the approximation of $\phi(x)$ as a statistical estimation task [33, 29]:

Definition 1** ([18]; Defn. 2.2).**

A probabilistic numerical method is a map $d_{e}\colon\mathcal{Y}_{e}\rightarrow\mathcal{P}_{\Phi}$, where $\mathcal{P}_{\Phi}$ is the set of all Borel probability distributions on $\Phi$.

The motivation for probabilistic numerical methods is to provide formal uncertainty quantification for the quantity of interest. Recall that the push-forward $\phi_{\#}\pi$ of a distribution $\pi$ is defined as $(\phi_{\#}\pi)(S)=\pi(\phi^{-1}(S))$ for all measurable $S$. Let $\pi_{Y|e}\coloneqq(y_{e})_{\#}\pi_{X}$ and denote the disintegration111If all random variables admit densities with respect to a reference measure then $\pi_{X|y,e}$ is the usual conditional distribution of $X$ given $y_{e}(X)$. A formal treatment of disintegration of measure is presented in [13]. of $\pi_{X}$ along the map $y_{e}$ by $\{\pi_{X|y,e}\}_{y\in\mathcal{Y}_{e}}$.

Definition 2** ([18]; Defn. 2.5).**

A probabilistic numerical method $d_{e}\colon\mathcal{Y}_{e}\rightarrow\mathcal{P}_{\Phi}$ is said to be Bayesian with prior $\pi_{X}$ if $d_{e}(y)=\phi_{\#}\pi_{X|y,e}$ for $\pi_{Y|e}$-almost all $y\in\mathcal{Y}_{e}$.

Thus the output of a Bayesian probabilistic numerical method is a posterior distribution providing uncertainty quantification for the quantity of interest. The foundations of Bayesian probabilistic numerical methods were established in [18] and Bayesian methods (in the strict sense of Defn. 2) have been studied for numerical integration [6, 32], global optimisation [28, 37, 38], ordinary [62] and partial [19, 49] differential equations and linear algebra [3, 17, 27]. Moreover, these methods are starting to see practical application [14, 40, 50].

The distributional output of a Bayesian probabilistic numerical method could of course be reduced to a point estimator for the quantity of interest, i.e. a (traditional) numerical method. Standard calculations from Bayesian decision theory, for example following [31], imply that the mean of $\phi_{\#}\pi_{X|y,e}$ is an average-case-optimal numerical method with $p=2$ (see Section 2.2). Conversely, [21] illustrated that the trapezoidal rule for numerical integration is a Bayes rule, which is in turn an average-case-optimal numerical method, when $\pi_{X}$ is the standard Weiner measure on $C(0,1)$. Despite their elegance, these connections rather overlook the main objective of a probabilistic numerical method, which is to provide formal uncertainty quantification for the quantity of interest. Indeed, to draw an analogy, it is well-understood that different strategies are required for the contrasting goals of parameter estimation and prediction in the statistical context [42]. In particular, the notion of optimal information put forward in average-case analysis (equivalently, Bayesian decision theory) is not necessarily an appropriate criteria on which to base the design of a Bayesian probabilistic numerical method.

The aim of this chapter is therefore twofold: First, we review connections between average-case-optimal numerical methods, average-case-optimal information and approaches to Bayesian experimental design from the statistical literature. Second, we discuss how optimal information could be defined for a Bayesian probabilistic numerical method, where the goal is to perform uncertainty quantification as opposed to direct estimation of a quantity of interest. In particular, we explore a particular criterion proposed in the recent work of [18] that has the advantage of being straightforward to approximate.

2 Bayesian Decision Theory

In this section we present a succinct overview of statistical decision theory [4], recalling that this is mathematically identical to average-case analysis, as explained in [31]. Throughout we use calligraphic font to represent (topological) state spaces, e.g. $\mathcal{X}$. Upper-case letters are used to denote (Borel) random variables, e.g. $X$. For convenience, we identify $X$ with a random variable on $\mathcal{X}$ whose distribution is denoted as $\pi_{X}$. The elements of a state space are denoted with lower-case letters, e.g. $x\in\mathcal{X}$.

Consider again an index set $\mathcal{E}$. In what follows we adopt the terminology of Bayesian decision theory and, in particular, refer to $e\in\mathcal{E}$ as an experiment. Our presentation will be more general than in Section 1 in three respects, which are all consequences of the broader context in which Bayesian decision-theoretic methods have been developed:

First, we associate to each experiment $e\in\mathcal{E}$ and each state $x\in\mathcal{X}$ a random variable $Y_{e}(x)$ taking values in $\mathcal{Y}_{e}$. The conditional distribution of this random variable is denoted $\pi_{Y|x,e}$. Often $Y_{e}(x)=y_{e}(x)$ with probability one, where $y_{e}(x)$ is a deterministic function of the state $x$, as in Section 1, in which case $\pi_{Y|x,e}=\delta(y_{e}(x))$ is an atomic measure centred on $y_{e}(x)\in\mathcal{Y}_{e}$. The more general formulation allows for the possibility that information about the state $x$ is corrupted by noise. Let $\pi_{Y|e}$ be the distribution defined through marginalisation as

[TABLE]

where $1_{S}$ is the indicator function for the measurable set $S\subset\mathcal{Y}_{e}$. The computational cost of each experiment is considered to be identical, but experiments differ in what information about $x$ is provided in $Y_{e}(x)$.

Second, we formulate the goal of Bayesian decision theory as the selection of an action $a$ from a specified set $\mathcal{A}$. Often, as in Section 1, the aim is to approximate $\phi(x)$ and the set of actions $\mathcal{A}$ is identical to $\Phi$. The more general formulation allows for more compex policy and control strategies to be considered in the statistical context. Correspondingly, a decision rule $d_{e}\in\mathcal{D}_{e}$ is considered to be a function $d_{e}\colon\mathcal{Y}_{e}\rightarrow\mathcal{A}$.

Third, we consider an arbitrary loss function $\ell\colon\mathcal{X}\times\mathcal{A}\rightarrow\mathbb{R}$. This includes the case $\ell(x,a)=\|\phi(x)-a\|_{\Phi}^{p}$ from Section 1. The choice of the loss function should be informed by the reason why the experiment is being conducted and, in an estimation context, should reflect the quantities of interest. See Chapter 2 of [4] for further discussion of this point.

Example 1** (Numerical integration).**

To illustrate the notation, consider the numerical task of integrating a continuous function $x\in\mathcal{X}=C(0,1)$. The aim here is to select an element $a\in\mathcal{A}=\mathbb{R}$ which represents an approximation to the quantity of interest $\phi(x)=\int_{0}^{1}x(t)\,\mathrm{d}t$. Here the experiment set $\mathcal{E}$ consists of vectors $e=[t_{1},\dots,t_{n-1}]$, along with fixed endpoints $t_{0}=0$, $t_{n}=1$, and $t_{0}\leq\dots\leq t_{n}$, corresponding to maps $y_{e}(x)=[x(t_{0}),\dots,x(t_{n})]$. One natural loss function is $\ell(x,a)=(\phi(x)-a)^{2}$. This problem has been well studied and will be used in an illustrative capacity in the sequel.

2.1 Bayes Risk

The framework of statistical decision theory allows comparison of different decision rules based on their Bayes risk:

[TABLE]

The Bayes risk quantifies the prior expected loss associated to a decision rule $d_{e}$ and therefore forms a natural criterion for the selection of a decision rule in the Bayesian context. We therefore define a Bayes rule for an experiment $e\in\mathcal{E}$ to be

[TABLE]

In a restricted context222Indeed, as highlighted in [31], Eq. (3) is identical to $\text{ACE}_{p}(e,d_{e})^{p}$ when $\pi_{Y|x,e}=\delta(y_{e}(x))$, $\mathcal{A}=\Phi$ and $\ell(x,a)=\|\phi(x)-a\|_{\Phi}^{p}$. the definition of a Bayes rule is mathematically identical to the definition of an average-case-optimal numerical method in Eq. (1). Similarly, we can define an optimal experiment

[TABLE]

which is analogous to average case optimal information in Eq. (2).

Example 2** (Numerical integration, continued).**

[54, 55]** considered the standard Wiener measure $\pi_{X}$, a Gaussian measure on $C(0,1)$ characterised by the property that for all $t,t^{\prime}\in[0,1]$ we have $\int x(t)\,\mathrm{d}\pi_{X}(x)=0$ and $\int x(t)x(t^{\prime})\,\mathrm{d}\pi_{X}(x)=\min(t,t^{\prime})$. In that work it was shown that for each experiment $e=[t_{1},\dots,t_{n-1}]$ there exists a unique Bayes rule

[TABLE]

which we recognise as the trapezoidal rule. Moreover, it was shown that there is a unique optimal experiment with $t_{i}=\frac{i}{n}$. Further contributions in this direction include [21, 33, 34, 41, 43, 44]. See [51] for a book-length treatment.

Remark 1** (Admissibility).**

It is important to note that other notions of optimality for decision rules, such as admissibility, need not coincide with the Bayesian notion of optimality. A decision rule $d_{e}\in\mathcal{D}_{e}$ is called admissible if there exists no $d_{e}^{\prime}\in\mathcal{D}_{e}$ such that

[TABLE]

for all $x\in\mathcal{X}$, with strict inequality for some $x\in\mathcal{X}$. The simplest illustration is estimation of $x\in\mathbb{R}$ based on $Y|x\sim N(x,1)$ and with $\ell(x,x^{\prime})=(x-x^{\prime})^{2}$, where an admissible decision rule is $d(y)=y$, but this is not a Bayes rule for any proper prior on $\mathbb{R}$. (Throughout this contribution, $N(\mu,\Sigma)$ denotes a Gaussian distribution with mean $\mu$ and covariance $\Sigma$.) Some results on when Bayes rules are, and are not, admissible can be found in the references provided in Section 8.4 of [4]. The famous result of [61] demonstrated that, under certain conditions, all admissible decision rules are so-called generalised Bayes rules, a definition in which improper priors are permitted. For the purposes of uncertainty quantification, however, we wish to remain in the Bayesian framework and consider only Bayes rules that arise from, and can be understood in terms of, a probabilistic model.

2.2 Bayes Acts

Following the definition of a Bayes rule it is reasonable to ask what actions a Bayes rule would select. To this end, denote the set of Bayes acts as

[TABLE]

The proof of the following result is provided in Appendix A.1:

Proposition 1**.**

A decision rule d_{e}^{\ast}\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{D}_{e}} is a Bayes rule if and only if $d_{e}^{\ast}(y)\in\mathcal{A}_{e}^{\ast}(y)$ for $\pi_{Y|e}$-almost all $y\in\mathcal{Y}_{e}$.

Often the set $\mathcal{D}_{e}$ will be infinite-dimensional, so that the search for a Bayes rule involves an optimisation problem which is infinite-dimensional. However, the set $\mathcal{A}_{e}$ is often finite-dimensional. Hence, the action of a Bayes rule is something that can often be computed.

The next result is technical and will be needed later, to deduce that the mean of $\pi_{X|y,e}$ is a Bayes act for a certain family of loss functions. Recall that a function $\phi\colon\mathcal{X}\rightarrow\mathbb{R}^{m}$ is called coercive if for all $c>0$ there exists a compact set $K\subset\mathcal{X}$ such that $\|\phi(x)\|_{2}\geq c$ for all $x\in\mathcal{X}\setminus K$. The proof of the following is provided in Appendix A.2:

Proposition 2**.**

Consider $\mathcal{A}=\mathcal{X}=\mathbb{R}^{d}$. Let $\ell(x,a)=\|\phi(x)-\phi(a)\|_{2}^{2}$ where $\phi\colon\mathcal{X}\rightarrow\mathbb{R}^{m}$, $m\in\mathbb{N}$. Assume that $\phi$ is twice continuously differentiable, that $\int\|\phi(x)\|_{2}\,\mathrm{d}\pi_{X|y,e}(x)$ is finite, and that the matrix

[TABLE]

has full row rank at all $a\in\mathcal{A}$. Then any Bayes act $a\in\mathcal{A}_{e}^{\ast}(y_{e})$ satisfies

[TABLE]

Moreover, if there exists a unique solution to Eq. (6) and the function $\phi$ is coercive, then this solution is a Bayes act.

For simplicity we have presented this result in finite dimensions and for the standard Euclidean norm, but it can be naturally extended to infinite dimensions and to an arbitrary Hilbert space $\mathcal{X}=\mathcal{A}=\Phi$.

Example 3** (Linear regression).**

Let $\pi_{X}=N(\mu_{0},\Sigma_{0})$, $X\in\mathbb{R}^{d}$, and $\pi_{Y|x,e}=N(A_{e}x,\Sigma)$, $Y\in\mathbb{R}^{n}$, where the matrix $\Sigma_{0}$ is positive definite and the matrix $A_{e}\in\mathbb{R}^{n\times d}$ is determined by the choice of experiment $e\in\mathcal{E}$. Consider a loss $\ell(x,x^{\prime})=(x-x^{\prime})^{\top}\Lambda(x-x^{\prime})$ where $\Lambda$ is a positive semi-definite matrix with a square root $\Lambda^{\frac{1}{2}}$. Then a Bayes decision rule $d_{e}^{\ast}$ is defined through the Bayes act(s) $a\in\mathbb{R}^{d}$ which, from Proposition 2, satisfy $\Lambda^{\frac{1}{2}}a=\Lambda^{\frac{1}{2}}\mu_{y,e}$ where $\pi_{X|y,e}=N(\mu_{y,e},\Sigma_{e})$ and $\Sigma_{e}=(A_{e}^{\top}A_{e}+\Sigma_{0}^{-1})^{-1}$, $\mu_{y,e}=\Sigma_{e}(A_{e}^{\top}y+\Sigma_{0}^{-1}\mu_{0})$. If $\Lambda$ is positive definite then it also follows from Proposition 2 that $\mu_{y,e}$ is the unique Bayes act.

The explicit connection between Bayesian decision theory and average-case analysis [31] can be exploited to obtain information-based complexity results for Bayes rules and optimal experiments in the decision-theoretic context. However, we recall that the differing goals of parameter estimation and uncertainty quantification need not lead to the same notions of an optimal experiment. Indeed, two posteriors $\pi_{X|y,e}$ can lead to the same Bayes act but provide very different uncertainty quantification for a quantity of interest; see Figure 1. To understand this point in detail, we turn to the more general framework of Bayesian experimental design in Section 3.

3 Bayesian Experimental Design

The term “Bayesian experimental design” has various historical usages, but in this work we are broadly consistent with the presentation of [12]. In this usage, Bayesian experimental design (BED) requires a utility function $u$ to be specified, such that $u(e,x,y)\in\mathbb{R}$ represents the “usefulness” of the information $y\in\mathcal{Y}_{e}$ when the true state of nature is $x\in\mathcal{X}$. Armed with a utility function, we can construct a design criterion

[TABLE]

and the following notion of an optimal experiment is considered:

[TABLE]

See [16, 36, 39]. In general optimal experimental designs will be analytically intractable, and in practice sophisticated numerical approaches to approximate the integral in Eq. (7) and to perform the global multivariate optimisation may be needed, e.g. [2, 39, 45, 46]. In an applied statistical context, the utility is something that itself may be elicited [63] and robustness to perturbation of the utility can be investigated [52].

The Bayesian experimental design framework is not concerned with selecting an action; it concerns only selection of an experiment. However, in terms of characterising an optimal experiment, the Bayesian experimental design framework is strictly more general than the Bayesian decision-theoretic approach, as we demonstrate next.

3.1 Design Criteria via Decision Theory

First, we explain how the decision-theoretic approach can be recovered in the experimental design framework [36, 20]. This is achieved by considering utilities $u(e,x,y)$ which are independent of the state $x$, in particular

[TABLE]

Under this choice, Eq. (7) becomes

[TABLE]

and so $\mathcal{E}_{\text{BED}}^{\ast}=\mathcal{E}_{\text{BDT}}^{\ast}$.

To overcome the fact that $\text{BR}(e,d_{e}^{\ast})$ is analytically intractable in general, a suite of approximations specific to the decision-theoretic context have been developed. Most approximations in the literature can be obtained by fixing a specific loss function $\ell$ and then making a Gaussian approximation to the posterior; see [20] for a more theoretical treatment. In this section we briefly present the approach and, for concreteness, we suppose $\mathcal{X}=\mathbb{R}^{m}$ so that $x$ is a column vector. Therefore, we consider Gaussian approximations $\pi_{X|y,e}\approx N(\mu_{y,e},\Sigma_{e})$. Note that the covariance $\Sigma_{e}$ is considered to be independent of the information $y_{e}$; a motivating example with this property was illustrated in Example 3. Then different choices of loss function $\ell$ lead to different approximations, the co-called alphabet criteria:

[TABLE]

for some positive semi-definite matrix $\Lambda$. Armed with an explicit approximation to $\text{BR}(e,d_{e}^{\ast})$, global multivariate optimisation over $e\in\mathcal{E}$ can be attempted; see the survey in [53].

Example 4** (Bayes A-, $c$- and E-optimality).**

Consider a loss $\ell(x,x^{\prime})=\|x-x^{\prime}\|_{\Lambda}^{2}$, where $\|x\|_{\Lambda}\coloneqq\|\Lambda^{\frac{1}{2}}x\|_{2}$ and $\Lambda$ is a positive semi-definite matrix with a square root $\Lambda^{\frac{1}{2}}$. Suppose that the posterior $\pi_{X|y,e}$ is a Gaussian $N(\mu_{y,e},\Sigma_{e})$. Then a Bayes decision rule $d_{e}^{\ast}$ is defined through the Bayes act(s) $a$, which satisfy $\Lambda^{\frac{1}{2}}a=\Lambda^{\frac{1}{2}}\mu_{y,e}$ due to Proposition 2. Now observe that $\int\ell(x,a)\,\mathrm{d}\pi_{X|y,e}(x)=\int(x-\mu_{y,e})^{\top}\Lambda(x-\mu_{y,e})\,\mathrm{d}\pi_{X|y,e}(x)=\text{\emph{tr}}(\Lambda\Sigma_{e})$, which is independent of $y_{e}$. It follows that $\text{\emph{BR}}(e,d_{e}^{\ast})=\text{\emph{tr}}(\Lambda\Sigma_{e})$. Selecting $e$ to minimise $\text{\emph{tr}}(\Lambda\Sigma_{e})$, or $\text{\emph{tr}}(\Sigma_{e})$ in the common case where $\Lambda=I$, is called Bayes A-optimal [47, 7, 8, 9, 22, 10]. In the special case where $\Lambda=cc^{\top}$ is a rank-1 matrix, the optimality criterion is called Bayes $c$-optimality [24]. Relatedly, minimisation of $\sup_{\|c\|=1}\text{\emph{tr}}(cc^{\top}\Sigma_{e})$ is called Bayes E-optimality [11].

Example 5** (Bayes D-optimality).**

Consider a loss $\ell(x,x^{\prime})=0$ if $\|x-x^{\prime}\|_{\Lambda}\leq\epsilon$ and, otherwise, $\ell(x,x^{\prime})=1$. Suppose that the posterior $\pi_{X|y,e}$ is a Gaussian $N(\mu_{y,e},\Sigma_{e})$. A Bayes decision rule $d_{e}^{\ast}$ is defined through the Bayes act(s), which one can verify include $\mu_{y,e}$. Now observe that $\int\ell(x,\mu_{y,e})\,\mathrm{d}\pi_{X|y,e}(x)=\int 1_{\|x-\mu_{y,e}\|_{\Lambda}>\epsilon}\,\mathrm{d}\pi_{X|y,e}(x)$, which is equal to one minus the probability that $\|Z\|_{2}\leq\epsilon$ where $Z\sim N(0,\Lambda^{1/2}\Sigma_{e}\Lambda^{1/2})$. For small $\epsilon$, this is $1-O(c_{d}^{-1}\text{\emph{det}}(\Lambda^{1/2}\Sigma_{e}\Lambda^{1/2})^{-d/2}\epsilon^{d})$ where $c_{d}\coloneqq 2^{\frac{d}{2}}\Gamma(\frac{d}{2}+1)$. Note in particular that this is independent of $y_{e}$. It follows that $\text{\emph{BR}}(e,d_{e}^{\ast})$ is minimised when $\text{\emph{det}}(\Lambda^{1/2}\Sigma_{e}\Lambda^{1/2})$ is minimised. This criterion is called Bayes D-optimality [56].

The notions above can be naturally extended to infinite-dimensional state spaces; see e.g. [1]. Given that the alphabet criteria are based on a Gaussian approximation of the posterior, it is incumbent on the analyst to attempt to verify the appropriateness of such approximations, for which a variety of diagnostics are available [15]. The alphabet criteria can also be derived as prior expectations of classical experimental design criteria, usually based on the Fisher information, an approach that was called pseudo-Bayesian in [53].

3.2 Design Criteria via Information Theory

To perform Bayesian experimental design a decision-theoretic basis is not essential. Indeed, other utilities — often based on information theory — have been proposed without reference to a decision-theoretic framework. For instance, [35] proposed

[TABLE]

where $D_{\text{KL}}$ denotes the Kullback-Leibler divergence (or information gain) between the posterior and the prior.

This criterion is not explicitly motivated by estimation of a quantity of interest and may be better suited to tasks such as uncertainty quantification for the state $x\in\mathcal{X}$. However, it is interesting to note that the proposal of [35] can actually also be conceived from the decision-theoretic framework, albeit in a way that is non-standard [5]. Indeed, consider an action set $\mathcal{A}=\mathcal{P}_{\mathcal{X}}$, whose elements $a\in\mathcal{A}$ are probability distributions on $\mathcal{X}$ and, for simplicity, existence of their Radon–Nikodym derivatives with respect to $\pi_{X}$ is assumed. Consider also a loss function of the form

[TABLE]

for arbitrary $f\colon\mathcal{X}\rightarrow\mathbb{R}$. Such a loss function is known as a proper scoring rule, since it can be shown that the unique Bayes act is to honestly report posterior belief; $a=\pi_{X|y,e}$, see Theorem 2 of [5]. Eq. (8) implies, in the particular case of $f\equiv 0$, that

[TABLE]

The theory of proper scoring rules can be developed further and the reader is referred to [23]. For other connections to information theory, we refer the reader to [26].

Next we explicitly address the challenge of designing a utility function that is commensurate with the aims of a Bayesian probabilistic numerical method.

4 Optimality Criteria for Probabilistic Numerical Methods

The purpose of a Bayesian probabilistic numerical method is to provide formal uncertainty quantification for the state $x\in\mathcal{X}$ or a derived quantity of interest $\phi(x)\in\Phi$. To build intuition, we immediately provide an example of such distributional output:

Example 6** (Numerical integration, continued).**

Consider the numerical integration task, set out in Example 1. The output of a Bayesian probabilistic numerical method, with prior $\pi_{X}$ equal to standard Weiner measure defined in Example 2 and information $y_{e}(x)=[x(t_{0}),\dots,x(t_{n})]$, is

[TABLE]

An overview of Bayesian probabilistic numerical methods for integration is provided in [6].

The distributional output can be naturally reduced to a point estimator; for instance the mean is recognised as the trapezoidal rule, a classical numerical method. In this spirit of estimation, one can apply the alphabet criteria from Section 3.1; the univariate Gaussian nature of the output implies that A-, $c$-, E- and D- optimal experiments $e$ all seek to minimise the variance $\frac{1}{12}\sum_{i=1}^{n}(t_{i}-t_{i-1})^{3}$. Thus a decision-theoretic perspective leads to an optimal experiment with uniformly spaced nodes $t_{i}=\frac{i}{n}$ .

However, the motivation for the development of probabilistic numerical methods is that the full distributional output provides formal uncertainty quantification for the quantity of interest. This objective is related but not identical to achieving strong performance at the estimation task. Indeed, for such uncertainty quantification to be useful, we put forward three desiderata to guide the design of a probabilistic numerical method $d_{e}$:

The output $d_{e}$ should “put mass close” to the true quantity of interest $\phi$, in a sense to be specified. 2. 2.

The output $d_{e}$ should be “calibrated”, in the sense that the true value $\phi(x)$ is well modelled as a random sample from the distributional output. 3. 3.

The search for an optimal experiment $e\in\mathcal{E}$ should, practically, not be too difficult.

Two important remarks are in order: First, notice the tension between (1) and (2), since (1) prefers an atomic mass placed on a “best” estimate, whilst (2) encourages probability density to be spread out so that the support at least contains the true value of the quantity of interest. Second, observe that for a Bayesian probabilistic numerical method it is not possible in general to ensure (2) through suitable choice of experiment $e\in\mathcal{E}$, since the output $d_{e}$ depends crucially on the prior $\pi_{X}$. Hence, our focus in the remainder is on satisfying (1) and (3) for a Bayesian probabilistic numerical method, with (2) being assumed through suitable elicitation of $\pi_{X}$. A detailed discussion of prior elicitation was provided in the supplement of [18] and sophisticated approaches to data-dependent elicitation have recently been developed, e.g. [30].

Remark 2** (Dependence on $\pi_{X}$).**

The dependence of a Bayesian decision rule $d_{e}$ on the prior $\pi_{X}$ is a fundamental aspect of Bayesian decision theory, thus the difficulties just described with regard to (2) are also fundamental. The alternative approach to decision theory pursued by Wald [60] circumvents the elicitation of a prior $\pi_{X}$ by casting the decision problem in an adversarial context, such that an adversary selects $\tilde{\pi}_{X}$, which need not coincide with $\pi_{X}$, from which $X\sim\tilde{\pi}_{X}$ is generated. The decision-maker then adopts a Bayes rule $d_{e}^{*}$, as in (4), based on a prior $\pi_{X}$ that solves the minimax optimisation problem

[TABLE]

*Wald’s approach is considered in detail in the context of probabilistic numerical methods in [48, 49], and forms an interesting point of contrast to the present discussion. *

4.1 The BPN Criterion

In what follows we consider the action space $\mathcal{A}$ to be equal to the state space $\mathcal{X}$. In that context, [18] proposed the utility

[TABLE]

which differs from Eq. (8) in that it depends explicitly on the true state $x$ and measures the posterior expected loss with respect to $x$. Thus two posteriors that give rise to the same Bayes act(s) need not give rise to the same utility function in Eq. (9); see Figure 2. The associated experimental design criterion, denoted BPN, is a special case of Eq. (7) with the above utility:

[TABLE]

To build intuition, for $\mathcal{X}$ a metric space with metric $d$ and $\ell(x,x^{\prime})=d(x,x^{\prime})^{p}$, the inner integral

[TABLE]

is the $p$th Wasserstein distance between the posterior $\pi_{X|y,e}$ and an atom on the true state $\delta(x)$. Thus it is seen that desideratum (1) is satisfied. Note also that that desideratum (3) is satisfied. Indeed, in contrast to the decision-theoretic approach to experimental design, where Gaussian approximations such as the alphabet criteria are practically necessary for an optimal experiment to be computed (to circumvent the repeated calculation of a Bayes act), no Gaussian approximation is needed in order to numerically compute $\text{BPN}(e)$. Indeed, we can re-write the BPN criterion in Eq. (10) as

[TABLE]

From this perspective it is clear that $\text{BPN}(e)$ can be directly computed using a Monte Carlo method to sample from $\pi_{X|y,e}$ in order to approximate the inner integral $(\ast)$. Although not closed form, we argue that the sophistication of modern Monte Carlo methods mean that desideratum (3) is essentially satisfied. The intergal $(\ast)$ has been termed a distance-based information function in [26] for the case $\ell(x,x^{\prime})=\|\phi(x)-\phi(x^{\prime})\|_{\Phi}^{2}$.

Example 7** (Numerical integration, continued).**

Unpacking Example 6, we note that $\pi_{X|y,e}$ is a collection of independent Brownian bridges $X_{i}\colon[t_{i-1},t_{i}]\rightarrow\mathbb{R}$ with $X_{i}(t_{i-1})=x(t_{i-1})$, $X_{i}(t_{i})=x(t_{i})$ and covariance function $\text{\emph{cov}}(X_{i}(t),X_{i}(t^{\prime}))=\frac{(t_{i}-t^{\prime})(t-t_{i-1})}{(t_{i}-t_{i-1})}$, $t_{i-1}\leq t\leq t^{\prime}\leq t_{i}$. Through direct calculation we obtain the closed form $\text{\emph{BPN}}(e)=\frac{1}{6}\sum_{i=1}^{n}(t_{i}-t_{i-1})^{3}$, which is easily minimised to obtain an optimal experiment by taking $t_{i}=\frac{i}{n}$. That this should coincide with the BDT-optimal experiment is explained later in Proposition 3.

In general the set of optimal experiments according to this criterion will be denoted

[TABLE]

For the remainder, our focus is on the properties of the BPN criterion and in particular we will ask, for this choice, whether $\mathcal{E}_{\text{BPN}}^{\ast}\stackrel{{\scriptstyle?}}{{=}}\mathcal{E}_{\text{BDT}}^{\ast}$.

4.2 Numerical Illustration

It was explained in Section 4.1 that the BPN criterion can be computed in circumstances that are quite general, as opposed to the BDT criteria for which Gaussian approximations are typically needed. The aim of this section is to illustrate the use of the BPN criteria on a problem for which BDT cannot easily be applied.

Consider the canonical linear elliptic partial differential equation

[TABLE]

on $D\coloneqq(0,1)^{2}$. Let $\ell(x,x^{\prime})=\left(\int(x(t)-x^{\prime}(t))^{p}\,\mathrm{d}t\right)^{1/p}$ and let the prior $\pi_{X}$ be Gaussian with $\int x(t)\,\mathrm{d}\pi_{X}(x)=0$, $\int x(t)x(t^{\prime})\,\mathrm{d}\pi_{X}(x)=\exp(-\|t-t^{\prime}\|_{2}^{2})$. An experiment $e\in\mathcal{E}$ determines the set of locations $\{t_{i}\}_{i=1}^{m}\subset D$, $\{t_{i}\}_{i=m+1}^{n}\subset\partial D$ at which the functions $f$ and $g$, respectively, are evaluated. To simplify the presentation we suppose that $m\ll n$ and focus on the information bottleneck, which is the placement of the points $\{t_{i}\}_{i=1}^{m}$ in $D$. To minimise $\text{BPN}(e)$ we recognised the triple integral in Eq. (12) as a single joint Gaussian integral and, to approximate this integral, we employed the standard Monte Carlo method.

For $p=2$ it can be shown that the optimal experiment according to BPN coincides with the optimal experiment according to BDT (see Proposition 3). The points $\{t_{i}\}_{i=1}^{m}$ were optimised sequentially in Figure 3, for convenience, and we note that a space-filling design was obtained. For $p=\infty$, in contrast, it is clear from Figure 4 that the optimal design from BPN is no longer space-filling. However, since it is not easily computed, it is unclear whether or not the optimal experiment according to BPN resembles the optimal experiment according to BDT.

4.3 Analysis and Open Questions

Our final aim is to address the question of whether $\mathcal{E}_{\text{BPN}}^{\ast}\stackrel{{\scriptstyle?}}{{=}}\mathcal{E}_{\text{BDT}}^{\ast}$. To this end, we first present a positive result in Proposition 3 and then present a negative result in Proposition 4.

4.3.1 A Positive Result

Under an assumption on the form of the loss $\ell$, the following equivalence was established in [18]:

Proposition 3**.**

Consider a loss function of the form $\ell(x,x^{\prime})=\|\phi(x)-\phi(x^{\prime})\|_{\Phi}^{2}$ where $\phi\colon\mathcal{X}\rightarrow\Phi$ takes values in an inner product space $\Phi$, with inner product $\langle\cdot,\cdot\rangle_{\Phi}$ and induced norm $\|\varphi\|_{\Phi}=\langle\varphi,\varphi\rangle_{\Phi}^{1/2}$. Suppose that the first conclusion of Proposition 2 holds (i.e. Eq. (6)). Then $\mathcal{E}_{\text{\emph{BPN}}}^{\ast}=\mathcal{E}_{\text{\emph{BDT}}}^{\ast}$.

For completeness, a slightly more concise proof of this result is included in Appendix A.3. This form of loss is commonly encountered and includes the Wasserstein distance with $p=2$ as a particular instance (see Eq. (11)). A consequence of Proposition 3 is that one can minimise the BPN criterion as an alternative to the BDT criteria in circumstances where the hypotheses of the Proposition hold.

4.3.2 A Negative Result

A converse result to Proposition 3 can also be constructed. The following is inspired by, but slightly more elegant than, a corresponding result in [18]:

Proposition 4**.**

Suppose that the state space $\mathcal{X}$ can be partitioned into three disjoint subsets, each with positive probability under $\pi_{X}$. Then there exists a loss function $\ell$ and a set of candidate experiments $\mathcal{E}$ such that $\mathcal{E}_{\text{\emph{BPN}}}^{\ast}\neq\mathcal{E}_{\text{\emph{BDT}}}^{\ast}$.

The proof of this result is constructive and is included in Appendix A.4. It is clear that the assumptions on $\mathcal{X}$ and $\pi_{X}$ are weak, to the point of being trivial.

4.3.3 Open Questions

An explicit characterisation of the loss functions $\ell$ for which $\mathcal{E}_{\text{BPN}}^{\ast}=\mathcal{E}_{\text{BDT}}^{\ast}$ is not, at least to our knowledge, available at present. In particular, the analytic intractability of optimal experiments in all but the simplest of numerical tasks leaves it unclear whether there exist a numerical task of practical importance for which $\mathcal{E}_{\text{BPN}}^{\ast}\neq\mathcal{E}_{\text{BDT}}^{\ast}$.

In general, any utility function $u$ from the Bayesian experimental design literature provides a criterion $\text{BED}(e)$ that can be studied from an information-based complexity standpoint. In particular, the issue of how $\text{BED}(e^{\ast})$ scales with $\text{dim}(\mathcal{X})$ when $e^{\ast}\in\mathcal{E}_{\text{BED}}^{\ast}$, so-called tractability, can be studied.

5 Discussion

The aim of this article was to build on the earlier work of [31], drawing attention to a wider range of optimality criteria that are used in the Bayesian experimental design literature and considering their use in the numerical context. In particular, we explained how the competing objectives of parameter estimation and uncertainty quantification lead, in general, to different notions of optimal information for a probabilistic numerical method.

One criterion, called BPN, was explored in detail. This formalised the idea that probability mass ought to be located close to the true quantity of interest. However, several other factors are also relevant in the design of a probabilistic numerical method and were not discussed. Indeed, the BPN criterion does not encode the notion that a posterior ought to be well calibrated, nor the notion that inferences should be robust to prior mis-specification, nor does it attempt more nuanced control (e.g. at the wall clock level) of computational cost [59]. In addition, no attempt was made to extend the notions introduced in this work to the adaptive context, where notions from sequential Bayesian experimental design are needed [25]. These are active areas of research and we look forward to seeing how these ideas are developed.

Appendix A Proofs

This appendix collects together proofs for the results quoted in the main text.

A.1 Proof of Proposition 1

Let $d_{e}^{\circ}\in\mathcal{D}_{e}$ be any decision rule such that $d_{e}^{\circ}(y)\in\mathcal{A}_{e}^{\ast}(y)$ for all $y\in\mathcal{Y}$. Then, for any other $d_{e}\in\mathcal{D}_{e}$,

[TABLE]

Thus $\text{BR}(e,d_{e}^{\circ})\leq\text{BR}(e,d_{e})$ and $d_{e}^{\circ}$ is a Bayes rule. Conversely, let $d_{e}^{\ast}\in\mathcal{D}_{e}$ be a Bayes rule. Then

[TABLE]

The definition of $d_{e}^{\circ}$ implies that $(\ast)\geq 0$. Thus $\ell(x,d_{e}^{\ast}(y))=\ell(x,d_{e}^{\circ}(y))$ holds $\pi_{Y|e}$-almost everywhere and so $d_{e}^{\ast}(y)\in\mathcal{A}_{e}^{\ast}(y)$ for $\pi_{Y|e}$-almost all $y\in\mathcal{Y}_{e}$, as required.

A.2 Proof of Proposition 2

Fix $e\in\mathcal{E}$ and $y_{e}\in\mathcal{Y}_{e}$. Consider the function

[TABLE]

and recall that the set $\mathcal{A}_{e}^{\ast}(y_{e})$ of Bayes acts is defined as the elements $a\in\mathcal{A}$ for which $f(a)$ is minimised. The assumptions on $\ell$ and $\phi$ imply that $f$ can be twice differentiated. The first derivative is easily computed: the integrability assumption on $\phi$ permits differentiation under the integral sign, yielding

[TABLE]

Since the matrix $\frac{\mathrm{d}\phi}{\mathrm{d}a}$ has full row rank, and since the space $\mathcal{A}$ is without boundary, if $a$ locally extremises $f$, then the term $(\ast)$ must vanish. This last condition is precisely Eq. (6).

For the converse result, assume that $\pi_{X|y,e}$ is a non-atomic distribution on $\mathcal{X}$, since otherwise the result is trivial. Let $\frac{\partial^{2}\phi}{\partial a_{i}\partial a_{j}}$ be the $m$-dimensional column vector with entries $\frac{\partial^{2}\phi_{k}}{\partial a_{i}\partial a_{j}}$. Then $f$ is seen to have continuous second derivative

[TABLE]

since $\phi$ is assumed to be twice continuously differentiable. The result will follow if $f$ is coercive, since a unique local extremum of a coercive function with continuous second derivatives must be a global minimum of that function. To this end, consider the reverse triangle inequality

[TABLE]

where $g$ is the quadratic $g(z)=c_{1}-2c_{2}z+z^{2}$. From Jensen’s inequality, $c_{2}^{2}\leq c_{1}$ with strict inequality if $\pi_{X|y,e}$ is not an atomic distribution on $\mathcal{X}$. Thus, the coercivity of $\phi$ implies the coercivity of $f$, and the proof is complete.

A.3 Proof of Proposition 3

To start, we indicate how the assumptions on the loss function will be used. Observe that

[TABLE]

where $\bar{\phi}_{y,e}=\int\phi(x)\,\mathrm{d}\pi_{X|y,e}(x)$. Indeed,

[TABLE]

and

[TABLE]

where the first and third terms have identical integrals under $\mathrm{d}\pi_{X|y,e}(x)\,\mathrm{d}\pi_{X|y,e}(x^{\prime})$ and

[TABLE]

Next, we use the fact that $\mathrm{d}\pi_{Y|x,e}(y)\,\mathrm{d}\pi_{X}(x)=\mathrm{d}\pi_{X|y,e}(x)\,\mathrm{d}\pi_{Y|e}(y)$ to observe that, for the loss $\ell(x,x^{\prime})=\|\phi(x)-\phi(x^{\prime})\|_{\Phi}^{2}$,

[TABLE]

Third, we recall that, for the loss $\ell(x,x^{\prime})=\|\phi(x)-\phi(x^{\prime})\|_{\Phi}^{2}$, a Bayes act $a\in\mathcal{A}$ must satisfy $\phi(a)=\bar{\phi}_{y,e}$ under the first conclusion of Proposition 2 (i.e. Eq. (6)). Therefore,

[TABLE]

where the final equality follows the same argument used to obtain Eq. (14). This final expression, in Eq. (16), is observed to be exactly half of the expression for $\text{BPN}(e)$ obtained in Eq. (15). Therefore, minimisation of Eq. (16) is equivalent to minimisation of Eq. (15) and it follows that $\mathcal{E}_{\text{BDT}}^{\ast}=\mathcal{E}_{\text{BPN}}^{\ast}$, as claimed.

A.4 Proof of Proposition 4

To start, suppose we are presented with a state space $\mathcal{X}$, that can be partitioned into three measurable subsets $①$, $②$, and $③$, with respective strictly positive probabilities $\pi_{\scriptsize①}$, $\pi_{\scriptsize②}$, and $\pi_{\scriptsize③}$. Without loss of generality we assume that $0<\pi_{\scriptsize①}\leq\pi_{\scriptsize②}\leq\pi_{\scriptsize③}<1$, noting that $\pi_{\scriptsize①}+\pi_{\scriptsize②}+\pi_{\scriptsize③}=1$. In what follows we exhibit a collection of candidate experiments $\mathcal{E}$ that, in effect, reveal some information about the state $x$ at the coarse level of the partition $\{①,②,③\}$. Let $s\colon\mathcal{X}\rightarrow\{①,②,③\}$ be the map that assigns each $x\in\mathcal{X}$ to its partition element. Consider a set of just two candidate experiments, in each case reporting a deterministic observation of the latent state:

[TABLE]

The aim is to decide whether experiment $e_{1}$ or $e_{2}$ should be performed. The informativeness of an experiment will be measured as described in the main text, based on the so-called 0-1 loss

[TABLE]

The 0-1 loss depends on the state $x$ only through the indicator function $1_{s(x)={\scriptsize①}}$, so that our task is equivalent to guessing whether the true $s(x)$ is $①$ or not, based on information obtained in the experiment. Thus, in a small abuse of notation, we may re-define the action space to be $\mathcal{A}=\{①,\neg①\}$.

Bayesian Decision Theory First we derive the BDT-optimal experiment(s) $\mathcal{E}_{\text{BDT}}^{\ast}$. Without loss of generality we can restrict attention, in the search for a Bayes rule, to the set of non-random decision rules of the form $d_{e}(y)=1_{y=0}a_{0}+1_{y=1}a_{1}$ for some $a_{1},a_{2}\in\mathcal{A}$ to be specified. The actions $a_{0},a_{1}$ are selected to minimise the risk

[TABLE]

Inputting the form of each experiment, we obtain

[TABLE]

The expression for $\text{BR}(e_{1},d_{e_{1}})$ is minimised by $a_{1}=\neg①$ and $a_{0}=\neg①$, since $\pi_{\scriptsize①},\pi_{\scriptsize②},\pi_{\scriptsize③}>0$. Similarly, the expression for $\text{BR}(e_{2},d_{e_{2}})$ is minimised by $a_{1}=\neg①$ and $a_{0}=\neg①$, since $\pi_{\scriptsize①}\leq\pi_{\scriptsize②}$ and $\pi_{\scriptsize③}>0$. It follows that $\text{BR}(e_{2},d_{e_{2}}^{\ast})=\text{BR}(e_{1},d_{e_{1}}^{\ast})=\pi_{\scriptsize①}$ and both experiments are optimal. Therefore, $\mathcal{E}_{\text{BDT}}^{\ast}=\{e_{1},e_{2}\}$.

Bayesian Experimental Design Now we derive the BPN-optimal experiment(s) $\mathcal{E}_{\text{BPN}}^{\ast}$. Through direct computation we see that

[TABLE]

Since $\pi_{\scriptsize②}\geq\pi_{\scriptsize①}>0$, it follows that experiment $e_{1}$ is the only optimal experiment. Therefore, $\mathcal{E}_{\text{BPN}}^{\ast}=\{e_{1}\}\neq\{e_{1},e_{2}\}=\mathcal{E}_{\text{BDT}}^{\ast}$, as claimed.

Acknowledgements:

The authors are grateful to an anonymous reviewer for their thoughtful comments, to Dave Woods for discussion of this work, as well as to the organisers and participants of the RICAM workshop on Multivariate Algorithms and Information-Based Complexity, for which this paper was prepared.

CJO and MG were supported by the Lloyd’s Register Foundation programme on Data-Centric engineering at the Alan Turing Institute, UK. TJS was supported by the Excellence Initiative of the German Research Foundation (DFG) through the Freie Universität Berlin. MG was supported by the EPSRC grants [EEP/P020720/1, EP/R018413/1, EP/R034710/1, EP/R004889/1] and a Royal Academy of Engineering Research Chair.