---
abstract: 'There is renewed interest in formulating integration as a statistical inference problem, motivated by obtaining a full distribution over numerical error that can be propagated through subsequent computation. Current methods, such as Bayesian Quadrature, demonstrate impressive empirical performance but lack theoretical analysis. An important challenge is therefore to reconcile these probabilistic integrators with rigorous convergence guarantees. In this paper, we present the first probabilistic integrator that admits such theoretical treatment, called Frank-Wolfe Bayesian Quadrature (FWBQ). Under FWBQ, convergence to the true value of the integral is shown to be up to exponential and posterior contraction rates are proven to be up to super-exponential. In simulations, FWBQ is competitive with state-of-the-art methods and out-performs alternatives based on Frank-Wolfe optimisation. Our approach is applied to successfully quantify numerical error in the solution to a challenging Bayesian model choice problem in cellular biology.'
author:
- |
    François-Xavier Briol\
    Department of Statistics\
    University of Warwick\
    `f-x.briol@warwick.ac.uk`\
    Chris J. Oates\
    School of Mathematical and Physical Sciences\
    University of Technology, Sydney\
    `christopher.oates@uts.edu.au`\
    Mark Girolami\
    Department of Statistics\
    University of Warwick\
    & The Alan Turing Institute for Data Science\
    `m.girolami@warwick.ac.uk`\
    Michael A. Osborne\
    Department of Engineering Science\
    University of Oxford\
    `mosb@robots.ox.ac.uk`\
bibliography:
- 'NIPS\_bib.bib'
title: 'Frank-Wolfe Bayesian Quadrature: Probabilistic Integration with Theoretical Guarantees'
---

Introduction
============

Computing integrals is a core challenge in machine learning and numerical methods play a central role in this area. This can be problematic when a numerical integration routine is repeatedly called, maybe millions of times, within a larger computational pipeline. In such situations, the cumulative impact of numerical errors can be unclear, especially in cases where the error has a non-trivial structural component. One solution is to model the numerical error statistically and to propagate this source of uncertainty through subsequent computations. Conversely, an understanding of how errors arise and propagate can enable the efficient focusing of computational resources upon the most challenging numerical integrals in a pipeline.

Classical numerical integration schemes do not account for prior information on the integrand and, as a consequence, can require an excessive number of function evaluations to obtain a prescribed level of accuracy [@OHagan1984]. Alternatives such as Quasi-Monte Carlo (QMC) can exploit knowledge on the smoothness of the integrand to obtain optimal convergence rates [@Dick2010]. However these optimal rates can only hold on sub-sequences of sample sizes $n$, a consequence of the fact that all function evaluations are weighted equally in the estimator [@Owen2014]. A modern approach that avoids this problem is to consider arbitrarily weighted combinations of function values; the so-called *quadrature rules* (also called cubature rules). Whilst quadrature rules with non-equal weights have received comparatively little theoretical attention, it is known that the extra flexibility given by arbitrary weights can lead to extremely accurate approximations in many settings (see applications to image de-noising [@Chen2015] and mental simulation in psychology [@Hamrick2013mental]).

Probabilistic numerics, introduced in the seminal paper of [@Diaconis1988], aims at re-interpreting numerical tasks as inference tasks that are amenable to statistical analysis.[^1] Recent developments include probabilistic solvers for linear systems [@Hennig2015solvers] and differential equations [@Conrad2015; @Schober2014]. For the task of computing integrals, Bayesian Quadrature (BQ) [@OHagan1991] and more recent work by [@Oates2015] provide probabilistic numerics methods that produce a full posterior distribution on the output of numerical schemes. One advantage of this approach is that we can propagate uncertainty through all subsequent computations to explicitly model the impact of numerical error [@Hennig2015ProbNum]. Contrast this with chaining together classical error bounds; the result in such cases will typically be a weak bound that provides no insight into the error structure. At present, a significant shortcoming of these methods is the absence of theoretical results relating to rates of posterior contraction. This is unsatisfying and has likely hindered the adoption of probabilistic approaches to integration, since it is not clear that the induced posteriors represent a sensible quantification of the numerical error (by classical, frequentist standards).

This paper establishes convergence rates for a new probabilistic approach to integration. Our results thus overcome a key perceived weakness associated with probabilistic numerics in the quadrature setting. Our starting point is recent work by [@Bach2012], who cast the design of quadrature rules as a problem in convex optimisation that can be solved using the Frank-Wolfe (FW) algorithm. We propose a hybrid approach of [@Bach2012] with BQ, taking the form of a quadrature rule, that (i) carries a full probabilistic interpretation, (ii) is amenable to rigorous theoretical analysis, and (iii) converges orders-of-magnitude faster, empirically, compared with the original approaches in [@Bach2012]. In particular, we prove that super-exponential rates hold for posterior contraction (concentration of the posterior probability mass on the true value of the integral), showing that the posterior distribution provides a sensible and effective quantification of the uncertainty arising from numerical error. The methodology is explored in simulations and also applied to a challenging model selection problem from cellular biology, where numerical error could lead to mis-allocation of expensive resources.

Background {#section:sigmapoint}
==========

Quadrature and Cubature Methods
-------------------------------

Let $\mathcal{X} \subseteq \mathbb{R}^d$ be a measurable space such that $d \in \mathbb{N}_{+}$ and consider a probability density $p(x)$ defined with respect to the Lebesgue measure on $\mathcal{X}$. This paper focuses on computing integrals of the form $\int f(x) p(x) \mathrm{d}x$ for a test function $f:\mathcal{X} \rightarrow \mathbb{R}$ where, for simplicity, we assume $f$ is square-integrable with respect to $p(x)$. A *quadrature rule* approximates such integrals as a weighted sum of function values at some design points $\{x_i\}_{i=1}^n \subset \mathcal{X}$: $$\int_\mathcal{X} f(x) p(x) \mathrm{d}x \approx \sum_{i=1}^n w_i f(x_i).$$ Viewing integrals as projections, we write $p[f]$ for the left-hand side and $\hat{p}[f]$ for the right-hand side, where $\hat{p} = \sum_{i=1}^n w_i \delta(x_i)$ and $\delta(x_i)$ is a Dirac measure at $x_i$. Note that $\hat{p}$ may not be a probability distribution; in fact, weights $\{w_i\}_{i=1}^n$ do not have to sum to one or be non-negative. Quadrature rules can be extended to multivariate functions $f:\mathcal{X} \rightarrow \mathbb{R}^d$ by taking each component in turn.

There are many ways of choosing combinations $\{x_i,w_i\}_{i=1}^n$ in the literature. For example, taking weights to be $w_i = 1/n$ with points $\{x_i\}_{i=1}^n$ drawn independently from the probability distribution $p(x)$ recovers basic Monte Carlo integration. The case with weights $w_i= 1/n$, but with points chosen with respect to some specific (possibly deterministic) schemes includes kernel herding [@Chen2010] and Quasi-Monte Carlo (QMC) [@Dick2010]. In Bayesian Quadrature, the points $\{x_i\}_{i=1}^n$ are chosen to minimise a posterior variance, with weights $\{w_i\}_{i=1}^n$ arising from a posterior probability distribution.

Classical error analysis for quadrature rules is naturally couched in terms of minimising the worst-case estimation error. Let $\mathcal{H}$ be a Hilbert space of functions $f: \mathcal{X}\rightarrow \mathbb{R}$, equipped with the inner product $\langle \cdot,\cdot \rangle_{\mathcal{H}}$ and associated norm $\|\cdot\|_{\mathcal{H}}$. We define the *maximum mean discrepancy* (MMD) as: $$\text{MMD}\big(\{x_i,w_i\}_{i=1}^n \big) {\coloneqq}\sup_{\substack{f \in \mathcal{H}: \|f\|_{\mathcal{H}}=1}} \big|p[f] - \hat{p}[f] \big|.$$ The reader can refer to [@Sriperumbudur2009] for conditions on $\mathcal{H}$ that are needed for the existence of the MMD. The rate at which the MMD decreases with the number of samples $n$ is referred to as the ‘convergence rate’ of the quadrature rule. For Monte Carlo, the MMD decreases with the slow rate of $\mathcal{O}_P(n^{-1/2})$ (where the subscript $P$ specifies that the convergence is in probability). Let $\mathcal{H}$ be a RKHS with reproducing kernel $k: \mathcal{X}\times \mathcal{X} \rightarrow \mathbb{R}$ and denote the corresponding canonical feature map by $\Phi(x) = k(\cdot,x)$, so that the mean element is given by $\mu_p(x) = p[\Phi(x)] \in \mathcal{H}$. Then, following [@Sriperumbudur2009] $$\text{MMD}\big(\{x_i,w_i\}_{i=1}^n \big) =  \| \mu_p - \mu_{\hat{p}} \|_{\mathcal{H}}.$$ This shows that to obtain low integration error in the RKHS $\mathcal{H}$, one only needs to obtain a good approximation of its mean element $\mu_p$ (as $\forall f \in \mathcal{H}$: $p[f] = \langle f , \mu_p \rangle_{\mathcal{H}}$). Establishing theoretical results for such quadrature rules is an active area of research [@Bach2015].

Bayesian Quadrature {#subsec:BQ}
-------------------

Bayesian Quadrature (BQ) was originally introduced in [@OHagan1991] and later revisited by [@Rasmussen2003; @Gunter2014] and [@Osborne2012]. The main idea is to place a functional prior on the integrand $f$, then update this prior through Bayes’ theorem by conditioning on both samples $\{x_i\}_{i=1}^n$ and function evaluations at those sample points $\{\mathrm{f}_i\}_{i=1}^n$ where $\mathrm{f}_i = f(x_i)$. This induces a full posterior distribution over functions $f$ and hence over the value of the integral $p[f]$. The most common implementation assumes a Gaussian Process (GP) prior $f \sim \mathcal{GP}(0,k)$. A useful property motivating the use of GPs is that linear projection preserves normality, so that the posterior distribution for the integral $p[f]$ is also a Gaussian, characterised by its mean and covariance. A natural estimate of the integral $p[f]$ is given by the mean of this posterior distribution, which can be compactly written as $$\hat{p}_{\text{BQ}}[f] = \mathrm{z}^T K^{-1} \mathrm{f}.
\label{BQmeaneq}$$ where $\mathrm{z}_i = \mu_p(x_i)$ and $K_{ij} = k(x_i,x_j)$. Notice that this estimator takes the form of a quadrature rule with weights $\mathrm{w}^{\text{BQ}} =\mathrm{z}^T K^{-1}$. Recently, [@Sarkka2015] showed how specific choices of kernel and design points for BQ can recover classical quadrature rules. This begs the question of how to select design points $\{x_i\}_{i=1}^n$. A particularly natural approach aims to minimise the posterior uncertainty over the integral $p[f]$, which was shown in [@Huszar2012 Prop. 1] to equal: $$v_{\text{BQ}}\big(\{x_i\}_{i=1}^n \big) \; = \; p [\mu_p] - \mathrm{z}^T K^{-1} \mathrm{z} \; = \; \text{MMD}^2\big(\{x_i,w_i^{\text{BQ}}\}_{i=1}^n \big).
\label{eq:variance}$$ Thus, in the RKHS setting, minimising the posterior variance corresponds to minimising the worst case error of the quadrature rule. Below we refer to Optimal BQ (OBQ) as BQ coupled with design points $\{x_i^\text{OBQ}\}_{i=1}^n$ chosen to globally minimise . We also call Sequential BQ (SBQ) the algorithm that greedily selects design points to give the greatest decrease in posterior variance at each iteration. OBQ will give improved results over SBQ, but cannot be implemented in general, whereas SBQ is comparatively straight-forward to implement. There are currently no theoretical results establishing the convergence of either BQ, OBQ or SBQ.

[*Remark:*]{} is independent of observed function values $\mathrm{f}$. As such, no active learning is possible in SBQ (i.e. surprising function values never cause a revision of a planned sampling schedule). This is not always the case: For example [@Gunter2014] approximately encodes non-negativity of $f$ into BQ which leads to a dependence on $\mathrm{f}$ in the posterior variance. In this case sequential selection becomes an [*active*]{} strategy that outperforms batch selection in general.

Deriving Quadrature Rules via the Frank-Wolfe Algorithm {#section:introFW}
-------------------------------------------------------

Despite the elegance of BQ, its convergence rates have not yet been rigorously established. In brief, this is because $\hat{p}_{\text{BQ}}[f]$ is an orthogonal projection of $f$ onto the [*affine*]{} hull of $\{\Phi(x_i)\}_{i=1}^n$, rather than e.g. the [*convex*]{} hull. Standard results from the optimisation literature apply to bounded domains, but the affine hull is not bounded (i.e. the BQ weights can be arbitrarily large and possibly negative). Below we describe a solution to the problem of computing integrals recently proposed by [@Bach2012], based on the FW algorithm, that restricts attention to the (bounded) convex hull of $\{\Phi(x_i)\}_{i=1}^n$.

The Frank-Wolfe (FW) algorithm (Alg. \[alg:FWalgorithm\]), also called the conditional gradient algorithm, is a convex optimization method introduced in [@Frank1956]. It considers problems of the form $\min_{g \in \mathcal{G}} J(g)$ where the function $J:\mathcal{G}\rightarrow\mathbb{R}$ is convex and continuously differentiable. A particular case of interest in this paper will be when the domain $\mathcal{G}$ is a compact and convex space of functions, as recently investigated in [@Jaggi2013]. These assumptions imply the existence of a solution to the optimization problem.

At each iteration $i$, the FW algorithm computes a linearisation of the objective function $J$ at the previous state $g_{i-1} \in \mathcal{G}$ along its gradient $(DJ)(g_{i-1})$ and selects an ‘atom’ $\bar{g}_i \in \mathcal{G}$ that minimises the inner product a state $g$ and $(DJ)(g_{i-1})$. The new state $g_{i} \in \mathcal{G}$ is then a convex combination of the previous state $g_{i-1}$ and of the atom $\bar{g}_{i}$. This convex combination depends on a step-size $\rho_i$ which is pre-determined and different versions of the algorithm may have different step-size sequences.

function $J$, initial state $g_1=\bar{g}_1 \in \mathcal{G}$ (and, for FW only: step-size sequence $\{\rho_i\}_{i=1}^n$). Compute $\bar{g}_i = \text{argmin}_{g \in \mathcal{G}} \big\langle g , (DJ)(g_{i-1}) \big\rangle_{\times} $ Update $g_{i} = (1 - \rho_i) g_{i-1} + \rho_i \bar{g}_i$

Our goal in quadrature is to approximate the mean element $\mu_p$. Recently [@Bach2012] proposed to frame integration as a FW optimisation problem. Here, the domain $\mathcal{G} \subseteq \mathcal{H}$ is a space of functions and taking the objective function to be: $$J(g) = \frac{1}{2}\big\| g - \mu_p \big\|^2_{\mathcal{H}}.$$ This gives an approximation of the mean element and $J$ takes the form of half the posterior variance (or the MMD$^2$). In this functional approximation setting, minimisation of $J$ is carried out over $\mathcal{G} = \mathcal{M}$, the marginal polytope of the RKHS $\mathcal{H}$. The marginal polytope $\mathcal{M}$ is defined as the closure of the convex hull of $\Phi(\mathcal{X})$, so that in particular $\mu_p \in \mathcal{M}$. Assuming as in [@Lacoste-Julien2015] that $\Phi(x)$ is uniformly bounded in feature space (i.e. $\exists R>0: \forall x \in \mathcal{X}$, $\|\Phi(x)\|_{\mathcal{H}} \leq R$), then $\mathcal{M}$ is a closed and bounded set and can be optimised over.

In order to define the algorithm rigorously in this case, we introduce the Fréchet derivative of $J$, denoted $DJ$, such that for $\mathcal{H}^*$ being the dual space of $\mathcal{H}$, we have the unique map $DJ:\mathcal{H} \rightarrow \mathcal{H}^*$ such that for each $g \in \mathcal{H}$, $(DJ)(g)$ is the function mapping $h \in \mathcal{H}$ to $(DJ)(g)(h) = \big\langle g - \mu, h \big\rangle_\mathcal{H}$. We also introduce the bilinear map $\langle \cdot, \cdot \rangle_{\times}: \mathcal{H} \times \mathcal{H}^* \rightarrow \mathbb{R}$ which, for $F \in \mathcal{H}^*$ given by $F(g) = \langle g, f \rangle_\mathcal{H}$, is the rule giving $\langle h, F \rangle_{\times} = \langle h, f \rangle_{\mathcal{H}}$.

A particular advantage of this method is that it leads to ‘sparse’ solutions which are linear combinations of the atoms $\{\bar{g}_i\}_{i=1}^n$ [@Bach2012]. In particular this provides a weighted estimate for the mean element: $$\label{eq:FWsparse}
\hat{\mu}_{\text{FW}} {\coloneqq}g_n = \sum_{i=1}^n \Big( \prod_{j=i+1}^{n} \big( 1 - \rho_{j-1} \big) \rho_{i-1} \Big) \bar{g}_i  {\coloneqq}\sum_{i=1}^n w_i^{\text{FW}} \bar{g}_i ,$$ where by default $\rho_0 = 1$ which leads to all $w_i^{\text{FW}} \in [0,1]$ when $\rho_i = 1/(i+1)$. A typical sequence of approximations to the mean element is shown in Fig. \[fig:designpoints\] (left), demonstrating that the approximation quickly converges to the ground truth (in black). Since minimisation of a linear function can be restricted to extreme points of the domain, the atoms will be of the form $\bar{g}_i = \Phi(x_i^{\text{FW}}) = k(\cdot ,x_i^{\text{FW}})$ for some $x_i^{\text{FW}} \in \mathcal{X}$. The minimisation in $g$ over $\mathcal{G}$ from step 2 in Algorithm \[alg:FWalgorithm\] therefore becomes a minimisation in $x$ over $\mathcal{X}$ and this algorithm therefore provides us design points. In practice, at each iteration $i$, the FW algorithm hence selects a design point $x_i^{\text{FW}} \in \mathcal{X}$ which induces an atom $\bar{g}_i$ and gives us an approximation of the mean element $\mu_p$. We denote by $\hat{\mu}_{\text{FW}}$ this approximation after $n$ iterations. Using the reproducing property, we can show that the FW estimate is a quadrature rule: $$\hat{p}_{\text{FW}}[f]
{\coloneqq}\big\langle f,\hat{\mu}_{\text{FW}} \big\rangle_{\mathcal{H}}
= \Big\langle f , \sum_{i=1}^n w_i^{\text{FW}} \bar{g}_i \Big\rangle_{\mathcal{H}}
= \sum_{i=1}^n  w_i^{\text{FW}} \big\langle f , k(\cdot,x_i^{\text{FW}}) \big\rangle_{\mathcal{H}}
= \sum_{i=1}^n w_i^{\text{FW}} f(x_i^{\text{FW}}).$$ The total computational cost for FW is $\mathcal{O}(n^2)$. An extension known as FW with Line Search (FWLS) uses a line-search method to find the optimal step size $\rho_i$ at each iteration (see Alg. \[alg:FWalgorithm\]). Once again, the approximation obtained by FWLS has a sparse expression as a convex combination of all the previously visited states and we obtain an associated quadrature rule. FWLS has theoretical convergence rates that can be stronger than standard versions of FW but has computational cost in $\mathcal{O}(n^3)$. The authors in [@Garber2015] provide a survey of FW-based algorithms and their convergence rates under different regularity conditions on the objective function and domain of optimisation.

[*Remark:*]{} The FW design points $\{x_i^{\text{FW}}\}_{i=1}^n$ are generally not available in closed-form. We follow mainstream literature by selecting, at each iteration, the point that minimises the MMD over a finite collection of $M$ points, drawn i.i.d from $p(x)$. The authors in [@Lacoste-Julien2015] proved that this approximation adds a $\mathcal{O}(M^{-1/4})$ term to the MMD, so that theoretical results on FW convergence continue to apply provided that $M(n) \rightarrow \infty$ sufficiently quickly. Appendix A provides full details. In practice, one may also make use of a numerical optimisation scheme in order to select the points.

A Hybrid Approach: Frank-Wolfe Bayesian Quadrature {#section:BayesianFW}
==================================================

To combine the advantages of a probabilistic integrator with a formal convergence theory, we propose Frank-Wolfe Bayesian Quadrature (FWBQ). In FWBQ, we first select design points $\{x_i^{\text{FW}}\}_{i=1}^n$ using the FW algorithm. However, when computing the quadrature approximation, instead of using the usual FW weights $\{w_i^{\text{FW}}\}_{i=1}^n$ we use instead the weights $\{w_i^{\text{BQ}}\}_{i=1}^n$ provided by BQ. We denote this quadrature rule by $\hat{p}_{\text{FWBQ}}$ and also consider $\hat{p}_{\text{FWLSBQ}}$, which uses FWLS in place of FW. As we show below, these hybrid estimators (i) carry the Bayesian interpretation of Sec. \[subsec:BQ\], (ii) permit a rigorous theoretical analysis, and (iii) out-perform existing FW quadrature rules by orders of magnitude in simulations. FWBQ is hence ideally suited to probabilistic numerics applications.

For these theoretical results we assume that $f$ belongs to a finite-dimensional RKHS $\mathcal{H}$, in line with recent literature [@Bach2012; @Garber2015; @Jaggi2013; @Lacoste-Julien2015]. We further assume that $\mathcal{X}$ is a compact subset of $\mathbb{R}^d$, that $p(x) > 0$ $\forall x \in \mathcal{X}$ and that $k$ is continuous on $\mathcal{X} \times \mathcal{X}$. Under these hypotheses, Theorem \[theorem:posteriormean\] establishes consistency of the posterior mean, while Theorem \[theo1\] establishes contraction for the posterior distribution.

\[theorem:posteriormean\] The posterior mean $\hat{p}_{\emph{FWBQ}}[f]$ converges to the true integral $p[f]$ at the following rates: $$\Big|p[f] - \hat{p}_{\emph{FWBQ}}[f]\Big| \leq \text{MMD}\big(\{x_i,w_i\}_{i=1}^n\big) \leq \left\{ \begin{array}{cl} \frac{2 D^2}{R}  n^{-1} & \text{for FWBQ}\\
 \sqrt{2}D \exp(-\frac{R^2}{2 D^2} n) & \text{for FWLSBQ}  \end{array} \right.$$ where the FWBQ uses step-size $\rho_i = 1/(i+1)$, $D \in (0,\infty)$ is the diameter of the marginal polytope $\mathcal{M}$ and $R \in (0,\infty)$ gives the radius of the smallest ball of center $\mu_p$ included in $\mathcal{M}$.

Note that all the proofs of this paper can be found in Appendix B. An immediate corollary of Theorem \[theorem:posteriormean\] is that FWLSBQ has an asymptotic error which is exponential in $n$ and is therefore superior to that of any QMC estimator [@Dick2010]. This is not a contradiction - recall that QMC restricts attention to uniform weights, while FWLSBQ is able to propose arbitrary weightings. In addition we highlight a robustness property: Even when the assumptions of this section do not hold, one still obtains atleast a rate $\mathcal{O}_P (n^{-1/2})$ for the posterior mean using either FWBQ or FWLSBQ [@Dunn1980].

[*Remark*]{}: The choice of kernel affects the convergence of the FWBQ method [@Hennig2015ProbNum]. Clearly, we expect faster convergence if the function we are integrating is ‘close’ to the space of functions induced by our kernel. Indeed, the kernel specifies the geometry of the marginal polytope $\mathcal{M}$, that in turn directly influences the rate constant $R$ and $D$ associated with FW convex optimisation.

Consistency is only a stepping stone towards our main contribution which establishes posterior contraction rates for FWBQ. Posterior contraction is important as these results justify, for the first time, the probabilistic numerics approach to integration; that is, we show that the [*full*]{} posterior distribution is a sensible quantification (at least asymptotically) of numerical error in the integration routine:

\[theo1\] Let $S \subseteq \mathbb{R}$ be an open neighbourhood of the true integral $p[f]$ and let $\gamma = \inf_{r \in S^C} | r - p[f]| >0$. Then the posterior probability mass on $S^c = \mathbb{R} \setminus S$ vanishes at a rate: $$\emph{prob}(S^c) \leq \left\{ \begin{array}{cl} \frac{2\sqrt{2}D^2}{\sqrt{\pi}R \gamma} n^{-1} \exp \Big(- \frac{\gamma^2 R^2}{8 D^4} n^2 \Big) 
& \text{for FWBQ} \\ 
\frac{2 D}{\sqrt{\pi} \gamma} \exp\Big( - \frac{R^2}{2 D^2}  n - \frac{\gamma^2}{2\sqrt{2}D} \exp\big( \frac{R^2}{2D^2} n\big)\Big)
& \text{for FWLSBQ} \end{array} \right.$$ where the FWBQ uses step-size $\rho_i = 1/(i+1)$, $D \in (0,\infty)$ is the diameter of the marginal polytope $\mathcal{M}$ and $R \in (0,\infty)$ gives the radius of the smallest ball of center $\mu_p$ included in $\mathcal{M}$.

The contraction rates are exponential for FWBQ and super-exponential for FWLBQ, and thus the two algorithms enjoy both a probabilistic interpretation and rigorous theoretical guarantees. A notable corollary is that OBQ enjoys the same rates as FWLSBQ, resolving a conjecture by Tony O’Hagan that OBQ converges exponentially \[personal communication\]:

\[theorem:OBQ\_rates\] The consistency and contraction rates obtained for FWLSBQ apply also to OBQ.

Experimental Results {#section:Experimental_Results}
====================

Simulation Study
----------------

To facilitate the experiments in this paper we followed [@Bach2015; @Bach2012; @Rasmussen2003; @Lacoste-Julien2015] and employed an exponentiated-quadratic (EQ) kernel $k(x, x') {\coloneqq}\lambda^{2} \exp ( -\nicefrac{1}{2 \sigma^2} \|x-x'\|^2_2)$. This corresponds to an infinite-dimensional RKHS, not covered by our theory; nevertheless, we note that all simulations are practically finite-dimensional due to rounding at machine precision. See Appendix E for a finite-dimensional approximation using random Fourier features. EQ kernels are popular in the BQ literature as, when $p$ is a mixture of Gaussians, the mean element $\mu_p$ is analytically tractable (see Appendix C). Some other $(p,k)$ pairs that produce analytic mean elements are discussed in [@Bach2015].

For this simulation study, we took $p(x)$ to be a 20-component mixture of 2D-Gaussian distributions. Monte Carlo (MC) is often used for such distributions but has a slow convergence rate in $\mathcal{O}_P(n^{-1/2})$. FW and FWLS are known to converge more quickly and are in this sense preferable to MC [@Bach2012]. In our simulations (Fig. \[fig:sim study\], left), both our novel methods FWBQ and FWLSBQ decreased the MMD much faster than the FW/FWLS methods of [@Bach2012]. Here, the same kernel hyper-parameters $(\lambda,\sigma) = (1,0.8)$ were employed for all methods to have a fair comparison. This suggests that the best quadrature rules correspond to elements [*outside*]{} the convex hull of $\{\Phi(x_i)\}_{i=1}^n$. Examples of those, including BQ, often assign negative weights to features (Fig. S1 right, Appendix D).

The principle advantage of our proposed methods is that they reconcile theoretical tractability with a fully probabilistic interpretation. For illustration, Fig. \[fig:sim study\] (right) plots the posterior uncertainty due to numerical error for a typical integration problem based on this $p(x)$. In-depth empirical studies of such posteriors exist already in the literature and the reader is referred to [@Chen2015; @Hamrick2013mental; @OHagan1991] for details.

Beyond these theoretically tractable integrators, SBQ seems to give even better performance as $n$ increases. An intuitive explanation is that SBQ picks $\{x_i\}_{i=1}^n$ to minimise the MMD whereas FWBQ and FWLSBQ only minimise an approximation of the MMD (its linearisation along $DJ$). In addition, the SBQ weights are optimal at each iteration, which is not true for FWBQ and FWLSBQ. We conjecture that Theorem \[theorem:posteriormean\] and \[theo1\] provide upper bounds on the rates of SBQ. This conjecture is partly supported by Fig. \[fig:designpoints\] (right), which shows that SBQ selects similar design points to FW/FWLS (but weights them optimally). Note also that both FWBQ and FWLSBQ give very similar result. This is not surprising as FWLS has no guarantees over FW in infinite-dimensional RKHS [@Jaggi2013].

Quantifying Numerical Error in a Proteomic Model Selection Problem
------------------------------------------------------------------

A topical bioinformatics application that extends recent work by [@Oates2014] is presented. The objective is to select among a set of candidate models $\{M_i\}_{i=1}^m$ for protein regulation. This choice is based on a dataset $\mathcal{D}$ of protein expression levels, in order to determine a ‘most plausible’ biological hypothesis for further experimental investigation. Each $M_i$ is specified by a vector of kinetic parameters $\theta_i$ (full details in Appendix D). Bayesian model selection requires that these parameters are integrated out against a prior $p(\theta_i)$ to obtain marginal likelihood terms $L(M_i) = \int p(\mathcal{D}|\theta_i) p(\theta_i) \mathrm{d}\theta_i$. Our focus here is on obtaining the [*maximum a posteriori*]{} (MAP) model $M_j$, defined as the maximiser of the posterior model probability $L(M_j) / \sum_{i=1}^m L(M_i)$ (where we have assumed a uniform prior over model space). Numerical error in the computation of each term $L(M_i)$, if unaccounted for, could cause us to return a model $M_k$ that is different from the true MAP estimate $M_j$ and lead to the mis-allocation of valuable experimental resources.

The problem is quickly exaggerated when the number $m$ of models increases, as there are more opportunities for one of the $L(M_i)$ terms to be ‘too large’ due to numerical error. In [@Oates2014], the number $m$ of models was combinatorial in the number of protein kinases measured in a high-throughput assay (currently $\sim 10^2$ but in principle up to $\sim 10^4$). This led [@Oates2014] to deploy substantial computing resources to ensure that numerical error in each estimate of $L(M_i)$ was individually controlled. Probabilistic numerics provides a more elegant and efficient solution: At any given stage, we have a fully probabilistic quantification of our uncertainty in each of the integrals $L(M_i)$, shown to be sensible both theoretically and empirically. This induces a full posterior distribution over numerical uncertainty in the location of the MAP estimate (i.e. ‘Bayes all the way down’). As such we can determine, on-line, the precise point in the computational pipeline when numerical uncertainty near the MAP estimate becomes acceptably small, and cease further computation.

The FWBQ methodology was applied to one of the model selection tasks in [@Oates2014]. In Fig. \[fig:model posteriors\] (left) we display posterior model probabilities for each of $m = 352$ candidates models, where a low number ($n = 10$) of samples were used for each integral. (For display clarity only the first 50 models are shown.) In this low-$n$ regime, numerical error introduces a second level of uncertainty that we quantify by combining the FWBQ error models for all integrals in the computational pipeline; this is summarised by a box plot (rather than a single point) for each of the models (obtained by sampling - details in Appendix D). These box plots reveal that our estimated posterior model probabilities are completely dominated by numerical error. In contrast, when $n$ is increased through 50, 100 and 200 (Fig. \[fig:model posteriors\], right and Fig. S2), the uncertainty due to numerical error becomes negligible. At $n = 200$ we can conclude that model $26$ is the true MAP estimate and further computations can be halted. Correctness of this result was confirmed using the more computationally intensive methods in [@Oates2014].

In Appendix D we compared the relative performance of FWBQ, FWLSBQ and SBQ on this problem. Fig. S1 shows that the BQ weights reduced the MMD by orders of magnitude relative to FW and FWLS and that SBQ converged more quickly than both FWBQ and FWLSBQ.

Conclusions
===========

This paper provides the first theoretical results for probabilistic integration, in the form of posterior contraction rates for FWBQ and FWLSBQ. This is an important step in the probabilistic numerics research programme [@Hennig2015ProbNum] as it establishes a theoretical justification for using the posterior distribution as a model for the numerical integration error (which was previously assumed [@Rasmussen2003; @Gunter2014; @Oates2015; @Osborne2012; @Sarkka2015 e.g.]). The practical advantages conferred by a fully probabilistic error model were demonstrated on a model selection problem from proteomics, where sensitivity of an evaluation of the MAP estimate was modelled in terms of the error arising from repeated numerical integration.

The strengths and weaknesses of BQ (notably, including scalability in the dimension $d$ of $\mathcal{X}$) are well-known and are inherited by our FWBQ methodology. We do not review these here but refer the reader to [@OHagan1991] for an extended discussion. Convergence, in the classical sense, was proven here to occur exponentially quickly for FWLSBQ, which partially explains the excellent performance of BQ and related methods seen in applications [@Gunter2014; @Osborne2012], as well as resolving an open conjecture. As a bonus, the hybrid quadrature rules that we developed turned out to converge much faster in simulations than those in [@Bach2012], which originally motivated our work.

A key open problem for kernel methods in probabilistic numerics is to establish protocols for the practical elicitation of kernel hyper-parameters. This is important as hyper-parameters directly affect the scale of the posterior over numerical error that we ultimately aim to interpret. Note that this problem applies equally to BQ, as well as related quadrature methods [@Bach2012; @Rasmussen2003; @Gunter2014; @Oates2015] and more generally in probabilistic numerics [@Schober2014]. Previous work, such as [@Hamrick2013mental], optimised hyper-parameters on a per-application basis. Our ongoing research seeks automatic and general methods for hyper-parameter elicitation that provide good frequentist coverage properties for posterior credible intervals, but we reserve the details for a future publication.

### Acknowledgments {#acknowledgments .unnumbered}

The authors are grateful for discussions with Simon Lacoste-Julien, Simo S[ä]{}rkk[ä]{}, Arno Solin, Dino Sejdinovic, Tom Gunter and Mathias Cronj[ä]{}ger. FXB was supported by EPSRC \[EP/L016710/1\]. CJO was supported by EPSRC \[EP/D002060/1\]. MG was supported by EPSRC \[EP/J016934/1\], an EPSRC Established Career Fellowship, the EU grant \[EU/259348\] and a Royal Society Wolfson Research Merit Award.

Supplementary Material {#supplementary-material .unnumbered}
======================

Appendix A: Details for the FWBQ and FWLSBQ Algorithms {#appendix-a-details-for-the-fwbq-and-fwlsbq-algorithms .unnumbered}
------------------------------------------------------

A high-level pseudo-code description for the Frank-Wolfe Bayesian Quadrature (FWBQ) algorithm is provided below.

function $f$, reproducing kernel $k$, initial point $x_0 \in \mathcal{X}$. Compute design points $\big\{x_i^{\text{FW}} \big\}_{i=1}^n$ using the FW algorithm (Alg. 1). Compute associated weights $\big\{w_i^{\text{BQ}}\big\}_{i=1}^n$ using BQ (Eqn. 4). Compute the posterior mean $\hat{p}_{\text{FWBQ}}[f]$, i.e. the quadrature rule with $\big\{ x_i^{\text{FW}}, w_i^{\text{BQ}} \big\}_{i=1}^n$. Compute the posterior variance $v_{\text{BQ}}\big(\{x_i^{\text{FW}}\}_{i=1}^n \big)$ using BQ (Eqn. 5). Return the full posterior $\mathcal{N} \big(\hat{p}_{\text{FWBQ}},v_{\text{BQ}}(\{x_i^{\text{FW}}\}_{i=1}^n)\big)$ for the integral $p[f]$.

Frank-Wolfe Line-Search Bayesian Quadrature (FWLSBQ) is simply obtained by substituting the Frank-Wolfe algorithm with the Frank-Wolfe Line-Search algorithm. In this appendix, we derive all of the expressions necessary to implement both the FW and FWLS algorithms (for quadrature) in practice. All of the other steps can be derived from the relevant equations as highlighted in Algorithm \[alg:FWBQalgorithm\] above.

The FW/FWLS are both initialised by the user choosing a design point $x_1^{\text{FW}}$. This can be done either at random or by choosing a location which is known to have high probability mass under $p(x)$. The first approximation to $\mu_p$ is therefore given by $g_1 = k(\cdot,x_1^{\text{FW}})$. The algorithm then loops over the next three steps to obtain new design points $\{x_i^{\text{FW}}\}_{i=2}^n$:

[*Step 1) Obtaining the new Frank-Wolfe design points $x_{i+1}^{\text{FW}}$.*]{}

At iteration $i$, the step consists of choosing the point $\bar{x}_i^{\text{FW}}$. Let $\{w_l^{(i)}\}_{l=1}^{i-1}$ denote the Frank-Wolfe weights assigned to each of the previous design points $\{x_l^{\text{FW}}\}_{l=1}^{i-1}$ at this new iteration, given that we choose $x$ as our new design point. The choice of new design point is done by computing the derivative of the objective function $J(g_{i-1})$ and finding the point $x^*$ which minimises the inner product:

$${\arg\min}_{g \in \mathcal{G}} \big\langle g , (DJ)(g_{i-1}) \big\rangle_{\times}$$

To do so, we need to obtain an equivalent expression of the minimisation of the linearisation of $J$ (denoted $DJ$) in terms of kernel values and evaluations of the mean element $\mu_p$. Since minimisation of a linear function can be restricted to extreme points of the domain, we have that $${\arg\min}_{g \in \mathcal{G}} \big\langle g , (DJ)(g_{i-1}) \big\rangle_{\times} 
\; = \; {\arg\min}_{x \in \mathcal{X}} \big\langle \Phi(x) , (DJ)(g_{i-1}) \big\rangle_{\times}.$$ Then using the definition of $J$ we have: $${\arg\min}_{x \in \mathcal{X}} \big\langle \Phi(x) , (DJ)(g_{i-1}) \big\rangle_{\times} 
\; = \; {\arg\min}_{x \in \mathcal{X}} \big\langle \Phi(x) , g_{i-1} - \mu_p \big\rangle_{\mathcal{H}} ,$$ where $$\begin{split}
\big\langle \Phi(x) , g_{i-1} - \mu_p \big\rangle_{\mathcal{H}} 
\quad & = \quad
\Big\langle \Phi(x), \sum_{l=1}^{i-1} w_l^{(i-1)} \Phi(x_l) - \mu_p \Big\rangle_{\mathcal{H}} \\
& = \quad
\sum_{i=1}^{i-1} w_l^{(i-1)} \big\langle   \Phi(x) , \Phi(x_l) \big\rangle_{\mathcal{H}} - \big\langle \Phi(x) ,  \mu_p \big\rangle_{\mathcal{H}} \\
& = \quad \sum_{l=1}^{i-1} w_l^{(i-1)} k( x,x_l) - \mu_p(x).
\end{split}$$

Our new design point $x_i^{\text{FW}}$ is therefore the point $x^*$ which minimises this expression. Note that this equation may not be convex and may require us to make use of approximate methods to find the minimum $x^*$. To do so, we sample $M$ points (where $M$ is large) independently from the distribution $p$ and pick the sample which minimises the expression above. From [@Lacoste-Julien2015] this introduces an additive error term of size $\mathcal{O}(M^{-1/4})$, which does not impact our convergence analysis provided that M(n) vanishes sufficiently quickly. In all experiments we took $M$ between $10,000$ and $50,000$ so that this error will be negligible.

It is important to note that sampling from $p(x)$ is likely to not be the best solution to optimising this expression. One may, for example, be better off using any other optimisation method which does not require convexity (for example, Bayesian Optimization). However, we have used sampling as the result from [@Lacoste-Julien2015] discussed above allows us to have a theoretical upper bound on the error introduced.

[*Step 2) Computing the Step-Sizes and Weights for the Frank-Wolfe and Frank-Wolfe Line-Search Algorithms.*]{}

Computing the weights $\{w_l^{(i)}\}_{l=1}^n$ assigned by the FW/FWLS algorithms to each of the design points is obtained using the equation: $$w_l^{(i)}=   \prod_{j=l+1}^{i} \big( 1 - \rho_{j-1} \big) \rho_{l-1}$$ Clearly, this expression depends on the choice of step-sizes $\{\rho_l\}_{l=1}^i$. In the case of the standard Frank-Wolfe algorithm, this step-size sequence is a an input from the algorithm and so computing the weights is straightforward. However, in the case of the Frank-Wolfe Line-Search algorithm, the choice of step-size is optimized at each iteration so that $g_i$ minimises $J$ the most.

In the case of computing integrals, this optimization step can actually be obtained analytically. This analytic expression will be given in terms of values of the kernel values and evaluations of the mean element.

First, from the definition of $J$ $$\begin{split}
J \big( (1-\rho) g_{i-1} + \rho \Phi(x_{i}) \big) 
 \quad & = \quad \frac{1}{2} \big\langle (1-\rho) g_{i-1} + \rho \Phi(x_{i}) - \mu_p , (1-\rho) g_{i-1} + \rho \Phi(x_{i}) - \mu_p \big\rangle_{\mathcal{H}} \\
& = \quad \frac{1}{2} \Big[ (1- \rho)^2 \big\langle g_{i-1},g_{i-1} \big\rangle_{\mathcal{H}} + 2 (1 - \rho) \rho \big\langle g_{i-1}, \Phi(x_{i}) \big\rangle_{\mathcal{H}} \\ 
&  \qquad + 2 \rho^2 \big\langle \Phi(x_{i}), \Phi(x_{i}) \big\rangle_{\mathcal{H}}  - 2 (1 - \rho) \big\langle g_{i-1}, \mu_p \big\rangle_{\mathcal{H}} \\
&  \qquad  - 2\rho \big\langle \Phi(x_{i}), \mu_p \big\rangle_{\mathcal{H}} + \big\langle \mu_p,\mu_p \big\rangle_{\mathcal{H}} \Big]. 
\end{split}$$ Taking the derivative of this expression with respect to $\rho$, we get: $$\begin{split}
\frac{\partial J \big( (1-\rho) g_{i-1} + \rho \Phi(x_{i}) \big)  }{\partial \rho}  
\quad & = \quad \frac{1}{2} \Big[ - 2(1- \rho) \big\langle g_{i-1} , g_{i-1} \big\rangle_{\mathcal{H}} + 2 (1 - 2 \rho) \big\langle g_{i-1}, \Phi(x_{i})\big\rangle_{\mathcal{H}}    \\
&  \qquad + 2 \rho \big\langle \Phi(x_{i}), \Phi(x_{i}) \big\rangle_{\mathcal{H}} + 2  \big\langle g_{i-1} , \mu_p \big\rangle_{\mathcal{H}} - 2
 \big\langle \Phi(x_{i}), \mu_p \big\rangle_{\mathcal{H}}  \Big] \\ 
&  = \quad \rho \Big[ \big\langle g_{i-1}, g_{i-1} \big\rangle_{\mathcal{H}} - 2 \big\langle g_{i-1} , \Phi(x_{i}) \big\rangle_{\mathcal{H}} + \big\langle \Phi(x_{i}) , \Phi(x_{i}) \big\rangle_{\mathcal{H}}  \\
 & = \quad \rho \big\| g_{i-1} - \Phi(x_{i}) \big\|^2_{\mathcal{H}} - \big\langle g_{i-1} - \Phi(x_{i}) , g_{i-1} - \mu_p \big\rangle_{\mathcal{H}} .\\
 \end{split}$$

Setting this derivative to zero gives us the following optimum: $$\rho^* \quad = \quad \frac{\Big\langle  g_{i-1} - \mu_p  , g_{i-1} - \Phi(x_{i}) \Big\rangle_{\mathcal{H}}}{\Big\| g_{i-1} - \Phi(x_{i}) \Big\|^2_{\mathcal{H}} }.$$

Clearly, differentiating a second time with respect to $\rho$ gives $\| g_{i-1} - \Phi(x_{i}) \|_{\mathcal{H}}^2$, which is non-negative and so $\rho^*$ is a minimum. One can show using geometrical arguments about the marginal polytope $\mathcal{M}$ that $\rho^*$ will be in $[0,1]$ [@Jaggi2013].

The numerator of this line-search expression is $$\begin{split}
\Big\langle g_{i-1} - \mu_p, g_{i-1} - \Phi(x_{i})\Big\rangle_{\mathcal{H}} \quad 
& = \quad \big\langle g_{i-1}, g_{i-1} \big\rangle_{\mathcal{H}} 
- \big\langle \mu_p, g_{i-1} \big\rangle_{\mathcal{H}} \\
& \qquad -  \sum_{l=1}^{i-1} w_l^{(i-1)} k(x_{l},x_{i})  + \mu_p(x_{i})  \\
& = \quad \sum_{l=1}^{i-1} \sum_{m=1}^{i-1} w_l^{(i-1)} w_m^{(i-1)} k(x_l,x_m) 
 \\
& \qquad - \sum_{l=1}^{i-1} w_l^{(i-1)} \Big[ k(x_l,x_{i})  + \mu_p(x_l) \Big] + \mu_p(x_{i}). \\
\end{split}$$ Similarly the denominator is $$\begin{split}
\big\| g_{i-1} - \Phi(x_{i})\big\|^2_{\mathcal{H}} \quad 
& = \quad \big\langle g_{i-1} - \Phi(x_{i}), g_{i-1} - \Phi(x_{i}) \big\rangle_{\mathcal{H}} \\
& = \quad \big\langle g_{i-1}, g_{i-1} \big\rangle_{\mathcal{H}}  
- 2 \big\langle g_{i-1} , \Phi(x_{i}) \big\rangle_{\mathcal{H}} 
+ \big\langle \Phi(x_{i}), \Phi(x_{i}) \big\rangle_{\mathcal{H}} \\
& = \quad \sum_{l=1}^{i-1} \sum_{m=1}^{i-1} w_l^{(i-1)} w_m^{(i-1)} k(x_l,x_m)  
- 2 \sum_{l=1}^{i-1} w_l^{(i-1)} k(x_l,x_{i}) 
+  k(x_{i},x_{i}).\\
\end{split}$$

Clearly all expressions provided here can be vectorised for efficient computational implementation.

[*Step 3) Computing a new approximation of the mean element.*]{}

The final step consists of updating the approximation of the mean element, which can be done directly by setting: $$g_{i} = (1 - \rho_i) g_{i-1} + \rho_i \bar{g}_i$$

Appendix B: Proofs of Theorems and Corollaries {#appendix-b-proofs-of-theorems-and-corollaries .unnumbered}
----------------------------------------------

The posterior mean $\hat{p}_{\emph{FWBQ}}[f]$ converges to the true integral $p[f]$ at the following rates: $$\Big|p[f] - \hat{p}_{\emph{FWBQ}}[f]\Big| \leq \text{MMD}\big(\{x_i,w_i\}_{i=1}^n\big) \leq \left\{ \begin{array}{cl} \frac{2 D^2}{R}  n^{-1} & \text{for FWBQ}\\
 \sqrt{2}D \exp(-\frac{R^2}{2 D^2} n) & \text{for FWLSBQ}  \end{array} \right.$$ where the FWBQ uses step-size $\rho_i = 1/(i+1)$, $D \in (0,\infty)$ is the diameter of the marginal polytope $\mathcal{M}$ and $R \in (0,\infty)$ gives the radius of the smallest ball of center $\mu_p$ included in $\mathcal{M}$.

The posterior mean in BQ is a Bayes estimator and so the MMD takes a minimax form [@Huszar2012]. In particular, the BQ weights perform no worse than the FW weights: $$\text{MMD} \Big( \big\{x_i^{\text{FW}} ,w_i^{\text{BQ}} \big\}_{i=1}^n \Big)
 \; = \; \inf_{\textrm{w} \in \mathbb{R}^n} \text{MMD} \Big( \big\{x_i^{\text{FW}},w_i \big\}_{i=1}^n \Big)  
 \; \leq \; \text{MMD} \Big( \big\{x_i^{\text{FW}},w_i^{\text{FW}} \big\}_{i=1}^n \Big). \label{eq:minimax}$$ Now, the values attained by the objective function $J$ along the path $\{g_i\}_{i=1}^n$ determined by the FW(/FWLS) algorithm can be expressed in terms of the MMD as follows: $$J(g_n) = \frac{1}{2} \big\|\hat{\mu}_{\text{FW}} - \mu_p \big\|^2_{\mathcal{H}} = \frac{1}{2} \text{MMD}^2\Big( \big\{x_i^{\text{FW}},w_i^{\text{FW}} \big\}_{i=1}^n \Big).
\label{FWobjvals}$$ Combining and gives $$\Big|p[f] - \hat{p}_{\text{FWBQ}}[f] \Big| \; \leq \; \text{MMD}\Big( \big\{x_i^{\text{FW}},w_i^{\text{BQ}} \big\}_{i=1}^n \Big) \big\|f \big\|_{\mathcal{H}} \; \leq \; 2^{1/2} J^{1/2}(g_n),$$ since $\|f\|_{\mathcal{H}} \leq 1$. To complete the proof we leverage recent analysis of the FW algorithm with steps $\rho_i = 1/(n+1)$ and the FWLS algorithm. Specifically, from [@Bach2012 Prop. 1] we have that: $$J(g_n) \leq \left\{ \begin{array}{cl} \frac{2D^4}{R^2} n^{-2} & \text{for FW with step size } \rho_i = 1/(i+1) \\
 D^2 \exp(-R^2 n/D^2 ) & \text{for FWLS} \end{array} \right.$$ where $D$ is the diameter of the marginal polytope $\mathcal{M}$ and $R$ is the radius of the smallest ball centered at $\mu_p$ included in $\mathcal{M}$.

Let $S \subseteq \mathbb{R}$ be an open neighbourhood of the true integral $p[f]$ and let $\gamma = \inf_{r \in S^C} | r - p[f]| >0$. Then the posterior probability mass on $S^c = \mathbb{R} \setminus S$ vanishes at a rate: $$\emph{prob}(S^c) \leq \left\{ \begin{array}{cl} \frac{2\sqrt{2}D^2}{\sqrt{\pi}R \gamma} n^{-1} \exp \Big(- \frac{\gamma^2 R^2}{8 D^4} n^2 \Big) 
& \text{for FWBQ, } \rho_i = 1/(i+1) \\ 
\frac{2 D}{\sqrt{\pi} \gamma} \exp\Big( - \frac{R^2}{2 D^2}  n - \frac{\gamma^2}{2\sqrt{2}D} \exp\big( \frac{R^2}{2D^2} n\big)\Big)
& \text{for FWLSBQ} \end{array} \right.$$ where $D \in (0,\infty)$ is the diameter of the marginal polytope $\mathcal{M}$ and $R \in (0,\infty)$ gives the radius of the smallest ball of center $\mu_p$ included in $\mathcal{M}$.

We will obtain the posterior contraction rates of interest using the bounds on the MMD provided in the proof of Theorem 1. Given an open neighbourhood $S \subseteq \mathbb{R}$ of $p[f]$, we have that the complement $S^c = \mathbb{R} \setminus S$ is closed in $\mathbb{R}$. We assume without loss of generality that $S^c \neq \emptyset$, since the posterior mass on $S^c$ is trivially zero when $S^c = \emptyset$. Since $S^c$ is closed, the distance $\gamma = \inf_{r \in S^c} \bigl|r - p[f]\bigr| > 0$ is strictly positive. Denote the posterior distribution by $\mathcal{N}(m_n,\sigma_n^2)$ where we have that $m_n {\coloneqq}\hat{p}_{\text{FWBQ}}[f]$ where $\hat{p}_{\text{FWBQ}} = \sum_{i=1}^n w_i^{\text{BQ}} \delta(x_i^{\text{FW}})$ and $\sigma_n {\coloneqq}\text{MMD}(\{x_i^{\text{FW}},w_i^{\text{BQ}}\}_{i=1}^n)$. Directly from the supremum definition of the MMD we have: $$\Big|p\big[f\big] - m_n \Big| \leq \sigma_n \big\|f\big\|_{\mathcal{H}}. \label{eq:CS}$$ Now the posterior probability mass on $S^c$ is given by $$M_n = \int_{S^c} \phi(r|m_n,\sigma_n) \mathrm{d}r,$$ where $\phi(r|m_n,\sigma_n)$ is the p.d.f. of the posterior normal distribution. By the definition of $\gamma$ we get the upper bound: $$\begin{aligned}
M_n & \leq & \int_{-\infty}^{p[f] - \gamma} \phi(r|m_n,\sigma_n) \mathrm{d}r + \int_{p[f] + \gamma}^\infty \phi(r|m_n,\sigma_n) \mathrm{d}r \\
& = & 1 + \Phi\Big(\underbrace{\frac{p[f] - m_n}{\sigma_n}}_{(*)} - \frac{\gamma}{\sigma_n}\Big) - \Phi\Big(\underbrace{\frac{p[f] - m_n}{\sigma_n}}_{(*)} + \frac{\gamma}{\sigma_n}\Big).\end{aligned}$$ From we have that the terms $(*)$ are bounded by $\|f\|_{\mathcal{H}} \leq 1 <\infty$ as $\sigma_n \rightarrow 0$, so that asymptotically we have: $$\begin{aligned}
M_n & \lesssim & 1 + \Phi\big(- \gamma / \sigma_n\big) - \Phi\big(\gamma / \sigma_n \big) \\
& = & \text{erfc}\big(\gamma/\sqrt{2}\sigma_n\big) \sim \big(\sqrt{2}\sigma_n / \sqrt{\pi} \gamma \big) \exp\big(- \gamma^2 / 2 \sigma_n^2 \big). \label{eq:final}\end{aligned}$$ Finally we may substitute the asymptotic results derived in the proof of Theorem 1 for the MMD $\sigma_n$ into to complete the proof.

The consistency and contraction rates obtained for FWLSBQ apply also to OBQ.

By definition, OBQ chooses samples that globally minimise the MMD and we can hence bound this quantity from above by the MMD of FWLSBQ: $$\text{MMD}\Big(\big\{x_i^{\text{OBQ}},w_i^{\text{BQ}}\big\}_{i=1}^n \Big) 
= \inf_{\{x_i\}_{i=1}^n \in \mathcal{X}} \text{MMD}\Big(\big\{x_i,w_i^{\text{BQ}}\big\}_{i=1}^n \Big)  \leq  \text{MMD}\Big( \big\{x_i^{\text{FW}},w_i^{\text{BQ}} \big\}_{i=1}^n \Big).$$ Consistency and contraction follow from inserting this inequality into the above proofs.

Appendix C: Computing the Mean Element for the Simulation Study {#appendix-c-computing-the-mean-element-for-the-simulation-study .unnumbered}
---------------------------------------------------------------

We compute an expression for $\mu_p(x) = \int_{- \infty}^{\infty} k(x,x') p(x') \mathrm{d}x'$ in the case where $k$ is an exponentiated-quadratic kernel with length scale hyper-parameter $\sigma$: $$k \big(x,x' \big)
 \; := \; \lambda^2 \exp \Big( \frac{-  \sum_{i=1}^d (x_{i} - x_{i}')^2  }{2 \sigma^2 }  \Big) 
 \; = \;  \lambda^2 (\sqrt{2 \pi} \sigma)^d  \phi \big(x \big| x',  \Sigma_{\sigma} \big),$$ where $\Sigma_\sigma$ is a d-dimensional diagonal matrix with entries $\sigma^2$, and where $p(x)$ is a mixture of d-dimensional Gaussian distributions: $$p(x) \quad = \quad \sum_{l=1}^L \rho_l \hspace{1mm} \phi \big(x\big|\mu_l,\Sigma_l \big).$$ (Note that, in this section only, $x_i$ denotes the $i$th component of the vector $x$.) Using properties of Gaussian distributions (see Appendix A.2 of [@Rasmussen2006]) we obtain $$\begin{split}
\mu_p(x) \quad 
& = \quad \int_{- \infty}^{\infty} k(x,x') p(x') \mathrm{d}x' \\
& = \quad \int_{- \infty}^{\infty} \lambda^2 (\sqrt{2 \pi} \sigma)^d \phi\big(x' \big| x,  \Sigma_{\sigma} \big) \times \Big( \sum_{l=1}^L \rho_l \hspace{1mm} \phi\big(x'\big|\mu_l,\Sigma_l \big)\Big) \mathrm{d}x' \\
& = \quad \lambda^2 (\sqrt{2 \pi} \sigma)^d \sum_{l=1}^L \rho_l \int_{- \infty}^{\infty} \phi\big(x' \big| x,  \Sigma_{\sigma} \big) \times  \phi\big(x'\big|\mu_l,\Sigma_l \big) \mathrm{d}x' \\
& = \quad \lambda^2 (\sqrt{2 \pi} \sigma)^d  \sum_{l=1}^L \rho_l \int_{- \infty}^{\infty} a_l^{-1} \phi\big(x'\big|c_l,C_l \big) \mathrm{d}x' \\
& = \quad \lambda^2  (\sqrt{2 \pi} \sigma)^d \sum_{l=1}^L \rho_l a_l^{-1} .\\
\end{split}$$ where we have: $$a_l^{-1} \; = \; (2 \pi)^{-\frac{d}{2}} \big| \Sigma_{\sigma} + \Sigma_{l} \big|^{-\frac{1}{2}} \exp \big( -\frac{1}{2} \big(x - \mu_l \big)^T \big( \Sigma_{\sigma} + \Sigma_{l} \big)^{-1} \big(x - \mu_l \big) \big).$$ This last expression is in fact itself a Gaussian distribution with probability density function $\phi(x|\mu_l, \Sigma_l + \Sigma_\sigma)$ and we hence obtain: $$\mu_p(x) \quad := \quad \lambda^2 \big(\sqrt{2 \pi} \sigma \big)^d \sum_{l=1}^L \rho_l \text{ } \phi\big(x|\mu_l, \Sigma_l + \Sigma_\sigma\big).$$ Finally, we once again use properties of Gaussians to obtain $$\begin{split}
 \int_{- \infty}^{\infty} \mu_p(x) p(x) \mathrm{d}x \quad  
& = \quad \int_{- \infty}^{\infty} \Big[ \lambda^2 \big(\sqrt{2 \pi} \sigma \big)^d \sum_{l=1}^L \rho_l \text{ } \phi\big(x|\mu_l, \Sigma_l + \Sigma_\sigma\big) \Big] \\
& \quad \times \Big[ \sum_{m=1}^L \rho_m \hspace{1mm} \phi\big(x\big|\mu_m,\Sigma_m \big) \Big] \mathrm{d}x \\
& = \quad \lambda^2 \big(\sqrt{2 \pi} \sigma \big)^d \sum_{l=1}^L \sum_{m=1}^L \rho_l \rho_m \int_{- \infty}^{\infty} \phi\big(x|\mu_l, \Sigma_l + \Sigma_\sigma\big) \phi\big(x\big|\mu_m,\Sigma_m \big) \mathrm{d}x \\
& = \quad \lambda^2 \big(\sqrt{2 \pi} \sigma \big)^d \sum_{l=1}^L \sum_{m=1}^L \rho_l \rho_m a_{lm}^{-1} \\
& = \quad \lambda^2 \big(\sqrt{2 \pi} \sigma \big)^d \sum_{l=1}^L \sum_{m=1}^L \rho_l \rho_m \phi\big(\mu_l|\mu_m,\Sigma_l+\Sigma_m+\Sigma_{\sigma} \big).
\end{split}$$ Other combinations of kernel $k$ and density $p$ that give rise to an analytic mean element can be found in the references of [@Bach2015].

Appendix D: Details of the Application to Proteomics Data {#appendix-d-details-of-the-application-to-proteomics-data .unnumbered}
---------------------------------------------------------

[*Description of the Model Choice Problem*]{}

The ‘CheMA’ methodology described in [@Oates2014] contains several elements that we do not attempt to reproduce in full here; in particular we do not attempt to provide a detailed motivation for the mathematical forms presented below, as this requires elements from molecular chemistry. For our present purposes it will be sufficient to define the statistical models $\{M_i\}_{i=1}^m$ and to clearly specify the integration problems that are to be solved. We refer the reader to [@Oates2014] and the accompanying supplementary materials for a full biological background.

Denote by $\mathcal{D}$ the dataset containing normalised measured expression levels $y_S(t_j)$ and $y_S^*(t_j)$ for, respectively, the unphosphorylated and phosphorylated forms of a protein of interest (‘substrate’) in a longitudinal experiment at time $t_j$. In addition $\mathcal{D}$ contains normalised measured expression levels $y_{E_i}^*(t_j)$ for a set of possible regulator kinases (‘enzymes’, here phosphorylated proteins) that we denote by $\{E_i\}$.

An important scientific goal is to identify the roles of enzymes (or ‘kinases’) in protein signaling; in this case the problem takes the form of variable selection and we are interested to discover which enzymes must be included in a model for regulation of the substrate $S$. Specifically, a candidate model $M_i$ specifies which enzymes in the set $\{E_i\}$ are regulators of the substrate $S$, for example $M_3 = \{E_2,E_4\}$. Following [@Oates2014] we consider models containing at most two enzymes, as well as the model containing no enzymes.

Given a dataset $\mathcal{D}$ and model $M_i$, we can write down a likelihood function as follows: $$\begin{aligned}
L(\theta_i,M_i) = \prod_{n=1}^N \phi\left( \frac{y_S^*(t_{n+1}) - y_S^*(t_n)}{t_{n+1} - t_n} \left|  \frac{-V_0 y_S^*(t_n)}{y_S^*(t_n) + K_0} + \sum_{E_j \in M_i} \frac{V_j y_{E_j}^*(t_n) y_S(t_n)}{y_S(t_n) + K_j} , \sigma_{\text{err}}^2  \right. \right). \label{chemalike}\end{aligned}$$ Here the model parameters are $\theta_i = \{\mathrm{K} , \mathrm{V}, \sigma_{\text{err}} \}$, where $(\mathrm{K})_j = K_j$, $(\mathrm{V})_j = V_j$, $\phi$ is the normal p.d.f. and the mathematical forms arise from the Michaelis-Menten theory of enzyme kinetics. The $V_j$ are known as ‘maximum reaction rates’ and the $K_j$ are known as ‘Michaelis-Menten parameters’. This is classical chemical notation, not to be confused with the kernel matrix from the main text. The final parameter $\sigma_{\text{err}}$ defines the error magnitude for this ‘approximate gradient-matching’ statistical model.

The prior specification proposed in [@Oates2014] and followed here is $$\begin{aligned}
\mathrm{K} & \sim & \phi_T \big(K \big| 1, 2^{-1} \mathrm{I} \big), \\
\sigma_{\text{err}} | \mathrm{K} & \sim & p(\sigma_{\text{err}}) \propto 1/\sigma_{\text{err}}, \\
\mathrm{V} | \mathrm{K},\sigma & \sim & \phi_T \big(V \big| 1, N \sigma_{\text{err}}^2 \big(\mathrm{X}(\mathrm{K})^T\mathrm{X}(\mathrm{K})\big)^{-1} \big),\end{aligned}$$ where $\phi_T$ denotes a Gaussian distribution, truncated so that its support is $[0,\infty)$ (since kinetic parameters cannot be non-negative). Here $\mathrm{X}(\mathrm{K})$ is the design matrix associated with the linear regression that is obtained by treating the $\mathrm{K}$ as known constants; we refer to [@Oates2014] for further details.

Due to its careful design, the likelihood in Eqn. \[chemalike\] is partially conjugate, so the following integral can be evaluated in closed form: $$L(\mathrm{K},M_i) = \int_{0}^{\infty} \int_{0}^{\infty} L(\theta_i,M_i) p(\mathrm{V},\sigma_{\text{err}} | \mathrm{K}) \mathrm{d}\mathrm{V} \mathrm{d}\sigma_{\text{err}}.$$ The numerical challenge is then to compute the integral $$L(M_i) = \int_{0}^{\infty} L(\mathrm{K},M_i)  p(\mathrm{K}) \mathrm{d}\mathrm{K},$$ for each candidate model $M_i$. Depending on the number of enzymes in model $M_i$, this will either be a 1-, 2- or 3-dimensional numerical integral. Whilst such integrals are not challenging to compute on a per-individual basis, the nature of the application means that the values $L(M_i)$ will be similar for many candidate models and, when the number of models is large, this demands either a very precise calculation per model or a careful quantification of the impact of numerical error on the subsequent inferences (i.e. determining the MAP estimate). It is this particular issue that motivates the use of probabilistic numerical methods.

[*Description of the Computational Problem*]{}

We need to compute integrals of functions with domain $\mathcal{X} = [0,\infty)^d$ where $d \in \{1,2,3\}$ and the sampling distribution $p(x)$ takes the form $\phi_T(x | 1,2^{-1} \mathrm{I})$. The test function $f(x)$ corresponds to $L(\mathrm{K},M_i)$ with $x = \mathrm{K}$. This is given explicitly by the $g$-prior formulae as: $$\begin{aligned}
L(\mathrm{K},M_i) & = & \frac{1}{(2 \pi)^{N/2}} \frac{1}{(N+1)^{d/2}} \Gamma\left(\frac{N}{2}\right) b_N^{-\frac{N}{2}}, \\
b_N & = & \frac{1}{2} \left( \mathrm{Y}^T\mathrm{Y} + \frac{1}{N} 1^T \mathrm{X}^T\mathrm{X} 1 - \mathrm{V}_N^T \mathrm{\Omega}_N \mathrm{V}_N \right), \\
\mathrm{V}_N & = & \mathrm{\Omega}_N^{-1} \left( \frac{1}{N} \mathrm{X}^T\mathrm{X} 1 + \mathrm{X}^T\mathrm{Y} \right), \\
\mathrm{\Omega}_N & = & \left(1 + \frac{1}{N} \right) \mathrm{X}^T \mathrm{X}, \\
(\mathrm{Y})_n & = &  \frac{y_S^*(t_{n+1}) - y_S^*(t_n)}{t_{n+1} - t_n}, \\\end{aligned}$$ where for clarity we have suppressed the dependence of $\mathrm{X}$ on $\mathrm{K}$. For the Frank-Wolfe Bayesian Quadrature algorithm, we require that the mean element $\mu_p$ is analytically tractable and for this reason we employed the exponentiated-quadratic kernel with length scale $\lambda$ and width scale $\sigma$ parameters: $$k(x,x') = \lambda^2 \exp\left(- \frac{\sum_{i=1}^d (x_i - x_i')^2}{2 \sigma^2}\right).$$ For simplicity we focussed on the single hyper-parameter pair $\lambda = \sigma = 1$, which produces: $$\begin{aligned}
\mu_p(x) & = & \int_{0}^{\infty} k(x,x') p(x') \mathrm{d}x' \\
& = & \int_{0}^{\infty} \exp\left(-\sum_{i=1}^d (x_i - x_i')^2\right) \phi_T \big(x' \big|1,2^{-1}\mathrm{I} \big) \mathrm{d}x' \\
& = & 2^{-d/2} \big(1 + \text{erf}(1) \big)^{-d} \prod_{i=1}^d \exp\left(-\frac{(x_i-1)^2}{2}\right) \left(1 + \text{erf}\left(\frac{x_i + 1}{\sqrt{2}}\right)\right),\end{aligned}$$ where $\phi_T$ is the p.d.f. of the truncated Gaussian distribution introduced above and $\text{erf}$ is the error function. To compute the posterior variance of the numerical error we also require the quantity: $$\int_{0}^{\infty} \int_{0}^{\infty} k(x,x') p(x) p(x') \mathrm{d}x \mathrm{d}x' = \int_{0}^{\infty} \mu_p(x) p(x) \mathrm{d}x = \left\{ \begin{array}{cl} 0.629907... & \text{for } d = 1 \\ 0.396783... & \text{for } d = 2 \\ 0.249937... & \text{for } d = 3 \end{array} \right. ,$$ which we have simply evaluated numerically. We emphasise that principled approaches to hyper-parameter elicitation are an important open research problem that we aim to address in a future publication (see discussion in the main text). The values used here are scientifically reasonable and serve to illustrate key aspects of our methodology.

FWBQ provides posterior distributions over the numerical uncertainty in each of our estimates for the marginal likelihoods $L(M_i)$. In order to propagate this uncertainty forward into a posterior distribution over posterior model probabilities (see Figs. 3 in the main text and \[fig:model posteriors2\] below), we simply sampled values $\hat{L}(M_i)$ from each of the posterior distributions for $L(M_i)$ and used these samples values to construct posterior model probabilities $\hat{L}(M_i) / \sum_j \hat{L}(M_j)$. Repeating this procedure many times enables us to sample from the posterior distribution over the posterior model probabilities (i.e. two levels of Bayes’ theorem). This provides a principled quantification of the uncertainty due to numerical error in the output of our primary Bayesian analysis.

[*Description of the Data*]{}

The proteomic dataset $\mathcal{D}$ that we considered here was a subset of the larger dataset provided in [@Oates2014]. Specifically, the substrate $S$ was the well-studied 4E-binding protein 1 (4EBP1) and the enzymes $E_j$ consisted of a collection of key proteins that are thought to be connected with 4EBP1 regulation, or at least involved in similar regulatory processes within cellular signalling. Full details, including experimental protocols, data normalisation and the specific choice of measurement time points are provided in the supplementary materials associated with [@Oates2014].

For this particular problem, biological interest arises because the data-generating system was provided by breast cancer cell lines. As such, the textbook description of 4EBP1 regulation may not be valid and indeed it is thought that 4EBP1 dis-regulation is a major contributing factor to these complex diseases (see [@Weinberg2006]). We do not elaborate further on the scientific rationale for model-based proteomics in this work.

![ Comparison of quadrature methods on the proteomics dataset. *Left:* Value of the MMD$^2$ for FW (black), FWLS (red), FWBQ (green), FWLSBQ (orange) and SBQ (blue). Once again, we see the clear improvement of using Bayesian Quadrature weights and we see that Sequential Bayesian Quadrature improves on Frank-Wolfe Bayesian Quadrature and Frank-Wolfe Line-Search Bayesian Quadrature. *Right:* Empirical distribution of weights. The dotted line represent the weights of the Frank-Wolfe algorithm with line search, which has all weights $w_i=1/n$. Note that the distribution of Bayesian Quadrature weights ranges from $-17.39$ to $13.75$ whereas all versions of Frank-Wolfe have weights limited to $[0,1]$ and have to sum to $1$.[]{data-label="fig:proteinsignalling"}](MMD2protein4.pdf){width="49.00000%"}

Appendix E: FWBQ algorithms with Random Fourier Features {#appendix-e-fwbq-algorithms-with-random-fourier-features .unnumbered}
--------------------------------------------------------

In this section, we will investigate the use of random Fourier features (introduced in [@Rahimi2007]) for the FWLS and FWLSBQ algorithms. An advantage of using this type of approximation is that the cost of manipulating the Gram matrix, and in particular of inverting it, goes down from $\mathcal{O}(n^3)$ to $\mathcal{O}(nD^2)$ for some user-defined constant $D$ which controls the quality of approximation. This could make Bayesian Quadrature more competitive against other integration methods such as MCMC or QMC. Furthermore, the kernels obtained using this method lead to finite-dimensional RKHS, which therefore satisfy the assumptions required for the theory in this paper to hold. This will be the aspect that we will focus on. In particular, we will show empirically that exponential convergence may be possible even when the RKHS is infinite-dimensional.

We will re-use the $20$-component mixture of Gaussians example with $d=2$ from our simulation studies, but using instead a random Fourier approximation of the exponentiated-quadratic (EQ) kernel $k(x,x'):= \lambda^2\exp(-1/2\sigma^2\|x-x'\|_2^2)$ with $(\lambda,\sigma)=(1,0.8)$ and $M=10000$.

Following Bochner’s theorem, we can always express translation invariant kernels in Fourier space: $$k(x,x') = \int_\mathcal{W} g(w) \exp\big(jw(x-x')\big)\mathrm{d}w 
= \mathbb{E}\Big[\exp\big(jw^Tx \big)\exp\big(jw^Tx' \big)\Big]$$ where $w \sim g(w) $ for $g(w)$ being the Fourier transform of the kernel. One can then use a Monte Carlo approximation of the kernel’s Fourier expression with $D$ samples whenever $g$ is a p.d.f.. Our approximated kernel will then lead to a $D$-dimensional RKHS and will be given by: $$k(x,x') \approx \frac{1}{D} \sum_{j=1}^D z_{w_j,b_j}(x)z_{w_j,b_j}(x') = \hat{k}_D(x,x')$$ where $z_{w_j,b_j}(x)=\sqrt{2}\cos(w_j^Tx+b_j)$ and $b_j \sim [0,2\pi]$ uniformly. Random Fourier features approximations are unbiased and, in the specific case of a $d$-dimensional EQ kernel with $\lambda=1$, we have to samples from the following Fourier transform: $$g(w) = \Big(\frac{2\pi}{\sigma^2}\Big)^{-\frac{d}{2}} \exp \Big(-\frac{\sigma^2\|w\|_2^2}{2} \Big)$$ which is a $d$-dimensional Gaussian distribution with zero mean and covariance matrix with all diagonal elements equal to $(1/\sigma^2)$.

The impact on the MMD from the use of random Fourier features to approximate the kernel for both the FWLS and FWLSBQ algorithms is demonstrated in Figure \[fig:RFF\_MMD\]. In this example, the quadrature rule uses the kernel with random features but the MMD is calculated using the original $\mathcal{H}$-norm. The reason for using this $\mathcal{H}$-norm is to have a unique measure of distance between points which can be compared.

Clearly, we once again have that the rate of convergence of the FWLSBQ is much faster than FWLS when using the exact kernel. The same phenomena is observed for the method with high number of random features ($D=5000$). This suggests that both the choice of design points and the calculation of the BQ weights is not strongly influenced by the approximation. It is also interesting to notice that the rates of convergence is very close for the exact and $D=5000$ methods (atleast when $n$ is small), potentially suggesting that exponential convergence is possible for the exact method. This is not so surprising in itself since using a Gaussian kernel represents a prior belief that the integrand of interest is very smooth, and we can therefore expect fast convergence of the method.

However, in the case with a smaller number of random features is used ($D=1000$), we actually observe a very poor performance of the method, which is mainly due to the fact that the weights are not well approximated anymore.

In summary, the experiments in this section suggest that the use of random features is a potential alternative for scaling up Bayesian Quadrature, but that one needs to be careful to use a high enough number of features. The experiments also give hope of having very similar convergence for infinite-dimensional and finite-dimensional spaces.

[^1]: A detailed discussion on probabilistic numerics and an extensive up-to-date bibliography can be found at <http://www.probabilistic-numerics.org>.
