Fast Bayesian Inference with Batch Bayesian Quadrature via Kernel Recombination
Authors: Masaki Adachi
Satoshi Hayakawa
Martin Jørgensen
Harald Oberhauser
Michael A Osborne
Venue: NeurIPS 2022 Accept (2022)
Keywords: Bayesian Quadrature
Kernel Quadrature
Gaussian Process
Active Learning
Model Evidence
Approximate Bayesian Computation

Abstract: Calculation of Bayesian posteriors and model evidences typically requires numerical integration. 
Bayesian quadrature (BQ), a surrogate-model-based approach to numerical integration, is capable of superb sample efficiency, but its lack of parallelisation has hindered its practical applications. 
In this work, we propose a parallelised (batch) BQ method, employing techniques from kernel quadrature, that possesses an empirically exponential convergence rate.
Additionally, just as with Nested Sampling, our method permits simultaneous inference of both posteriors and model evidence.
Samples from our BQ surrogate model are re-selected to give a sparse set of samples, via a kernel recombination algorithm, requiring negligible additional time to increase the batch size.
Empirically, we find that our approach significantly outperforms the sampling efficiency of both state-of-the-art BQ techniques and Nested Sampling in various real-world datasets, including lithium-ion battery analytics.

Fast Bayesian Inference with Batch Bayesian Quadrature via Kernel Recombination
Masaki Adachi"
Machine Learning Research Group, University of Oxford Toyota Motor Corporation
masaki @robots o..ac. uk
Satoshi Hayakawa", Harald Oberhauser Mathematical Institute, University of Oxford (hayakawa, (berhauser)@maths ox ac uk
Martin Jørgensen, Michael A. Osborne
Machine Learning Research Group, University of Oxford (martinj, mosb)@robots ox.ac.ut
Abstract
Calculation of Bayesian posteriors and model evidences typically requires numer- ical integration. Bayesian quadrature (BQ), surrogate-modell-aate approach to numerical integration, is capable of superb sample efficiency, but its lack of paralleli- sation has hindered its practical applications. In this work, we propose parallelised (batch) BQ method, employing techniques from kernel quadrature, that possesses empirically exponential convergence rate. Additionally, just as with Nested Sam- pling, our method permits simultaneous inference of both posteriors and model evi- dence. Samples from our BQ surrogate model are re-selected to give sparse set of samples, via a kernel recombination algorithm, requiring negligible additional time to increase the batch size. Empirically, we find that our approach significantly out- performs the sampling efficiency both state-of-the-art BQ techniques and Nested Sampling in various real-world datasets, including lithium-ion battery analytics
1 Introduction
Many applications in science, engineering, and economics involve complex simulations explain the structure and dynamics of the process. Such models are derived from knowledge of the mechanisms and principles underlying the data-generating process, and are critical for scientific hypothesis- building and testing. However, dozens plausible simulators describing the same phenomena often exist, owing to differing assumptions levels of approximation. Similar situations can be found in selection of simulator- datasets, and in many data assimilation applications (28).
control models, selection of machine learning models on large-scale
In such settings, with multiple competing models, choosing the best model for the dataset crucial. Bayesian model evidence gives clear criteria for such model selection. However, computing model evidence requires integration over the likelihood, which is challenging, particularly when the likelihood is mon-closed-foo and/or expensive. The ascertained model is often applied to
"qqual contribution
Code:
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

produce posteriors for prediction and parameter estimation afterwards. There are many algorithms specialised for the calculation of model evidences or posteriors, although only limited number of Bayesian inference solvers estimate both model evidence and posteriors in one go. As such, costly computations are often repeated (at least) twice. Addressing this concern, nested sampling (NS) [66 was developed to estimate both model evidence and posteriors simultaneously, and has been broadly applied, especially amongst astrophysicists for cosmological model selection [63]. However, NS is based on a MMonte Carlo (MC)sampler sampler, and its slow convergence rate is practical hindrance.
To aid NS, and other approaches, parallel computing widely applied to improve the speed of wall-clock computation. Modern computer clusters and graphical processing units enable scientists to query the likelihood in large batches. However, parallelisation can, at best, linearly accelerate NS, doing little to counter NS's inherently slow convergence rate as MC sampler.
This paper investigates batch (Bayesian quadrature (BQ) 665 for fast Bayesian inference. BQ solves the integral as an inference problem, modelling the likelihood function with a probabilistic model (typically a GGusssian process (GP)). Gunter et al. [37] proposed |Warped sequential active Bayesian integration (WSABD) which adopts active learning to select samples upon uncertainty over the integrand. WSABI showed that BQ with expensive GP calculations could achieve faster convergence in wall time than cheap MC samplers. Wagstaff et al. (78) introduced batch WSABI, achieving even faster calculation via parallel computing and became the fastest BQ model to date. We improve upon these existing works for a large-scale batch case.
2 Background
Vanilla Bayesian quadrature While BQ in general is the method for the integration, the functional approximation nature permits solving the following integral Z and obtaining the surrogate function of posterior p(z) simultaneously in the Bayesian inference context:
p(x)= Luve (x)
Luve (2)
J Lave (z (x)
where both (x)) (e.g. a likelihood) and π (ar) (e.g. a prior) are non-negative, and x E Rd is a sample, and is sampled from prior x ~ ((x). BQ solves the above integral an inference problem,
modelling a likelihood function {(z) by a GP in order to construct surrogate model of the expensive true likelihood L (z). The surrogate likelihood function f(x) is modelled:
Ely GP(4:my; Cp),
(2a)
my(x) K(x, K(x,X) (XX,
(2h)
Cy K(x,
(z, X))
(2c)
where X Raxd the matrix of observed samples, y E R is the observed true likelihood values, K is the kernel. Due to linearity, the mean and variance of the integrals are simply
E [y] / (y(z))(z))
(3a)
Varlz [y] = // Cy(l, dx dz'.
(3b)
In particular, (3) becomes analytic when (Ix) is Gaussian and K is squared exponential kernel, ((X,x)) = u /2 WWW(X; W), where kernel variance and W is the diagonal covariance matrix whose diagonal elements are the lengthscales of each dimension. Since both the mean and variance of the integrals can be calculated analytically, posterior and model evidence can be obtained simultaneously. Note that non-Gaussian prior and kernel are also possible to be chosen for modelling via kernel recombination (see Supplementary). Still, we use this combination throughout this paper for simplicity.
"In GP modelling, the GP likelihood function is modelled as GP K). and 2 is the resulting posterior GP. Throughout the paper, we refer to symmetric positive semi definite kernel just as a kernel. The notations and refer to being sampled from and being conditioned, respectively.

---

Warped Bayesian quadrature (WSABI) WSABI with linearisation approximation (WSABI D adoots the square-root warping GP for non con-negatvvity with linearisation approximation of the transform م a + _ The square-root GP defined as l ~ GP(L; my, Cy), and we have the following linear approximation:
y GP( ,
(4a)
m³ (x) + m
(4b)
(x,,') my(rr) (1,, (my()))
(4c)
Gaussianity implies the model evidence and posterior p()) remain analytical (see Supplementary).
Kernel quadrature in general kernel quadrature (KQ) the group of numerical integration rules for calculating the integral of function classes that form the reproducing kernel Hilbert space (RKHS)
With KQ rule Qw.x given by weights (w,)== and points (xx))==> approximate the integral by the weighted sum
n Qu = Σ w,)((;)) ~ { MMo)m(aade
(5)
where h is function of RKHS associated with the kernel K. We define its worst-case error by wee(Q == su xh) h(x)r(z)) dr). Surprisingly, it is shown in [48] that we have
ZarlZ] yy] inf wce(Quxx22
Thus, the point configuration in KQ with small worst-case emror gives good way select points to reduce the integral variance in Bayesian quadrature.
Random Convex Hull Quadrature (RCHQ) Recall from and (6) that we wish approximate the integral of function h in the current RKHS. First, we prepare n test functions 22...... Pa-1 based on M sample points using the Nyström approximation of the kernel: (p((x) == u, K (Xays: x). where u, RR- the i-th eigenvector of K (Xnyss Xnys). If we let A be the i-th eigenvalue of the same matrix, the following gives a practical approximation [55
Ko(2, v) = Σ (())))))))))
Next, we consider extracting a weighted set of points (Wquad- Xquad) from a set of N points X
with positive weights Wree. We do it by the so-called kernel recombination algorithm [59] 75), so that the measure induced by (Wquad, X quua)) exactly integrates the above test functions f1, ....... Pp-1 with respect to the measure given by (Wre, xzcc
In the actual implementation of multidimensional case, we execute the kernel recombination not by the algorithm [75] with the best known computational complexity OCC, n log( / /n) (where C. is the cost of evaluating (p=))=] at a point), but the one of [59] using an LP solver (Gurobi (38) for this time) with empirically faster computational time. We also adopt the randomized
SVD 40 for the Nyström approximation, so we have a computational time empirically faster than O(NM + M log + Mn log( N/n))
44 in practice.
3 Related works
Bayesian inference for intractable likelihood Inference with intractable likelihoods a long- standing problem, and plethora of methods have been proposed. Most infer posterior and evidence separately, and hence are not our fair competitors, as solving both is more challenging. For posterior inference, Markov Chain Monte Carlo [62 43), particularly Hamilton Monte Carlo [47), is the gold standard. In a Likelihood-free inference context, kernel density estimation (KDE) with Bayesian
= 0.8 x min(y). See Supplementary for the details

---

optimisation 39 and neural networks [36] surrogates are proposed for simulation- based inference
25 In a Bayesian coreset context, scalable Bayesian inference [60], sparse variational inference
20, active learning [68] have been proposed for large-scale dataset inference. However, all of the above only calculate posteriors, not model evidences. For evidence inference, Annealed Importance Sampling [64], and Bridge Sampling 99. Sequential Monte Carlo (SMC)
19
27 are popular, but only
estimate evidence, not the posterior.
Bayesian quadrature The early works on BQ, which directly replaced the likelihood function with a GP [65] 66, 69], did not explicitly handle non non-negative integrand constraints. Osborne et al. [67] introduced logarithmic warped GP to handle the non-negativity of the integrand, and introduced active learning for BQ, a method that selects samples based on where the variance of the integral will be minimised. Gunter et al. [37] introduced square-root GP to make the integrand function closed- form and to speed up the computation. Furthermore, they accelerated active learning by changing the objective from the variance of the integral Varl (Zy to simply the variance of the integrand Var((())T(<)]] Wagstaff et al. (78) introduced the first batch BQ. Chai et al. [22] generalised warped GP for BQ, and proposed rroii-transformation. BQ has been extended to machine learning tasks (model selection [2I], manifold learning (30], kernel learning 4I) with new aaquisition luunction (AF) designed for each purpose. For posterior inference, VBMC M has pioneered that BQ can infer posterior and evidence in one go via variational inference, and 2 has extended it to the noisy likelihood case. This approach is an order of magnitude more expensive than the WSABI approach because it essentially involves BQ inside of an optimisation loop for variational inference. This paper employs WSABI-L and its AF for simplicity and faster computation. Still, our approach a natural extension of BQ methods and compatible with the above advances. (e.g. changing the RCHQ kernel into the prospective uncertainty sampling AF for VMBC)
Kernel quadrature There are number of KQ algorithms from herding/optimisstion 22 6 48
to random sampling 8. In the context of BQ, Trank-Wolfe Bayesian quadrature (FWBQ) 13
using kernel herding has been proposed. This method proposes to do BQ with the points given by herding, but the guaranteed exponential convergence ((exp(<en) is limited to finie-dimensiona. RKHS, which is not the case in setting. For general kernels, the convergence rate drops to O(//Vn) [6. Recently, a random sampling method with good convergence rate was proposed for dimeite-dimensional RKHS (4). Based Carathéodory's theorem for the convex hull, it efficiently calculates a reduced measure for a larger empirical measure. In this paper, we call it random convex
hull quadrature (RCHQ) and use it in combination with the BQ methods. (See Supplementary)
4 Proposed Method: BASQ
4.1
General idea
We now introduce our algorithm, named Bayesian alternately subsampled quadrature (BASQ)
Kernel recombination for batch selection Batch WSABI (78) selects batch samples based on the AF, taking samples having the maximum AF values greedily via gradient- based optimisation with multiple starting points. The computational cost of this sampling scheme increases exponentially with the batch size and/or the problem dimension. Moreover, this method is often only locally optimal (79 We adopt scalable, gradient-free optimisation via KQ algorithm based on kernel recombination 44
Surprisingly, Huszár et al. [48] pointed out the equivalence between BQ and KQ with optimised
weights. KQ can select n samples from the N candidate samples X to efficiently reduce the worst-case error. The problem in batch BQ is selecting n samples from the probability measure (2) that minimises integral variance. When subsampling N candidate samples X
7 T(z), we can regard this samples X as an empirical measure Temp(z) (x) approximating the true probability measure (x) n < N. Therefore, applying KQ to select n samples that can minimise Var (ar(((()))))))) Temp (4:) is equivalent to selecting n batch samples for batch BQ. As more observations make surrogate model t(x) more accurate, the empirical integrand model ((mp))))(( approaches to the true integrand model frme (x)T(()) This subsampling scheme allows us to apply any KQ methods for batch BQ. However, such a dual quadrature scheme tends to be computationally demanding. Hayakawa et al. 44 proposed an efficient algorithm based on kernel recombination, RCHQ, which automatically

---

(c) Uncertainty Sampling
(1) Recombination
Log number of sssmples
Log wall time (S)
(a)
8-
(b)
(owe(s)
likelibood
6
5
prior (ss)
3
2
posterior ((l)
1
3
5
2500
5000
1
Log batch size
wall time (8)
Log MAE of integral
Log MAE of inteerrll
GP-modelled likelihood #(+))
0
(c)
(d)
-3
-3
oequsstion function ApA
sampling funcrion
subsaapples
updated GP ((+)
6
1000
2000
500 wall time (s)
1000
number samples
densé
sparse
BASQ (Ours)
batch WSABI
Log MAE -1.82
Log MAE -5.42
Figure I: Performance comparison of our algorithm BASQ against batch WSABI 78 All evaluation
was performed with the likelihood of a mixture of N-dimensional Gaussians. (a), (b), (c), (d) 10-dimensional Gaussians, (e), (f) univariate Gaussians.
returns a sparse set of n samples based on Carathéodory's theorem. The computational cost of batch size for our algorithm, BASQ, is lower than O M M22 log Mn² log(N/m) /n) [44).
Alternately subsampling The performance of RCHQ relies upon the quality of a predefined kemel. Thus, we add BQ elements to KQ return; making RCHQ an online algorithm. Throughout the sequential batch update, we pass the BQ-updated kernels to RCHQ This enables RCHQ to exploit the function shape information This corresponds to the batch maximisation of the model evidence for ((z). Then, BQ optimises the hyperparameters based on samples from true likelihood (x), which corresponds to an optimised kernel preset for RCHQ in the next round. These alternate updates characterise our algorithm, BASQ (See Supplementary)
This
select the best batch samples for minimising Var(E(x)) Temp (+)).
Importance sampling for uncertainty We added one more BQ element RCHQ; uncertainty sampling. RCHQ relies upon the quality of subsamples from the prior distribution. However, sharp, multimodal, or high-dimensional likelihoods make it challenging to find prior subsamples that overlap over the likelihood distribution. This typically happens in Bayesian inference when the likelihood expressed as the product of Gaussians, which gives rise very sharp likelihood peak. The likelihood of big data tends to become multimodal [70]. Therefore, we adopt importance sampling, gravitating subsamples toward the meaningful region, and correct the bias via its weights. We propose a mixture of prior and AF guiding proposal distribution. The prior encourages global (rather than just local) optimisation, and the AF encourages sampling from uncertain areas for faster convergence. However, sampling from AF is expensive. We derive efficient sparse Gaussian mixture sampler. Moreover, introducing square-root warping 37 enables the sampling distribution to factorise, yielding faster sampling.
Summary of contribution
We summarised the key differences between batch WSABI and BASQ (a) shows that BASQ is more scalable in batch size. (b) clarifies that BASQ can sample
in Figure 6
'For updating kernel, kernel to updated is Cy. not the kernel K. The kernel K just corresponds to the prior belief the distribution of f, so once we have observed the samples X (and y), the variance to be minimised becomes Cy.
'We randomly varied number components between 10 and 15, setting their variance uniformly random between and 4, and setting their means uniformly at random within the box bounded by [-3.3] in all dimensions. The weights of Gaussians were randomly generated from uniform distribution, but set to be one after integration, mean absolute error (MAE) was adopted for the evaluation of integral estimation.
5

---

Table 1: BASQ algorithm
Algorithm 1: Bayesian Alternately Subsampled Quadrature (BASQ)
Notation:
Viit initial guess, k: convergence criterion, n: batch size. :M: Nyström sample size, N: recombination sample size, f: GP-modelled likelihood, m true likelihood,
[æ)) prior, p(z): posterior. A(r): AF, K: kerne,,
Snyss Sree: samplers for Nyström and recombination, respectively
Input: prior Tæ)), true likelihood function lrue Output: posterior p(r)), the mean and variance of model evidence IE|Z yy], varlZ ly]
1: Cy. ; A(xx) + InitialiseGPs(rnnnt))
# Initialise GPs with initial guess Set samplers
2: Snys. Sree f- ([0Sampler(A(x)
3: while varlZ y) > k:
4: Xnys Says (M)) (M)
Samples M points for test function yp # Samples N points for recombination
5:
(Wree: X See (N)
6: PI,... pn-1 + Nyström (Nyys Ci²y)
# Define n 1 test functions
7:
solve a kernel recombination problem
8: Find an n-point subset Xyudd X,e and Wqual
9: s.t. w Pi = Weo (Xre), w "quad 1
w rec
10:
return Xquad
The sparse set of n samples
11:
y Parallel mne quaa)))
# Parallel computing of likelihood
12:
K + Update( (Xquad. y)
# Train GPs
13:
C},((z) + Opte)yeessThee)pdde((((()
Type II MLE optimisation
14:
Snys, Sree f ResetSampler( (()))) T(x))
# Reset samplers with the updated A(x)
15:
E[Z y). var Z y) +- BayesQuad( Cj. (x)) # Calculate via Eqs.
and 4
16:
return p(r), E Z y var [ y
10 to 100 times many samples in the same time budget WSABI, supported by the efficient sampler and RCHQ. (c) states the convergence rate of BASQ is faster than WSABI, regardless of computation time. (d) demonstrates the combined acceleration in wall time. While the batch WSABI reached 10-1 after 1,000 seconds passed, BASQ was within seconds. (e) and (f), visualised how RCHQ selects sparser samples than batch WSABI. This clearly explains that gradient-free kernel recombination is better in finding the global optimal than multi-star start optimisation. These results demonstrate that we were able to combine the merits of BQ and KQ (see Supplementary). We further tested with various synthetic datasets and real-world tasks in the fields of lithium ion batteries and material science. Moreover, we mathematically analyse the convergence rate with proof in Section 6
4.2 Algorithm
Table Dillustrates the pseudo-code for our algorithm. Rows 4 10 correspond to RCHQ, and rows 11 15 correspond BQ. We can use the variance of the integral Varl / y convergence criterion. For hyperparameter optimisation, we adopt the type-Il maximum likelihood estimation (MLE) to
optimise hyperparameters via L-BFGS [8] for speed.
Importance sampling for uncertainty Lemma in the supplementary proves the optimal upper bound of the proposal distribution g(x) R / ({(x,))(xx) (a, where f(xx) == (σ). However, sampling from square-root variance intractable, so we linearised to g(x)
0.5(1 + Of (z). To correct the lincarisation error, the coefficient 0.5 was changed into the hyperparameter r, which is defined as follows:
((x) - ())r)f())++)))))) ( 05r51
(8)
A(x) = (2,)))(())
99)
wis(x) = f(z)/g()),
(10)

---

where Wis is the weight of the importance sampling. While = 1 becomes the pure uncertainty sampling, r 0 is the vanilla MC sampler.
Efficient sampler Sampling from A(z), a mixture of Gaussians, is expensive, as some mixture weights are negative, preventing the usual practice of weighted sampling from each Gaussian. As the warped kernel C (z. x) is also computationally expensive, we adopt a factorisation trick:
1 Z = ((x) (()) = a+ / [(()) do Rs a + 2 / ((x))L(()))) (z)
(II)
where we have changed the distribution of interest to f(x) (z) This is doubly beneficial. Firstly, the distribution of interest fxx will be updated over iterations. The previous f(x) π (a) means the subsample distribution eventually obeys prior, which is disadvantageous if the prior does not overlap the likelihood sufficiently. On the contrary, the new f(r) narrows its region via |nū (z)]- Secondly, the likelihood function changed to [,(x)), thus the kernels shared with RCHQ changed into cheap warped kernel Cy (z. æt. This reduces the computational cost of RCHQ, and the sampling cost of A(x). Now A(x) (x)))((z)) which is also Gaussian mixture, but the number of components is significantly lower than the original AF ( As Cy(x) is positive, the positive weights of the Gaussian mixture should cover the negative components. Interestingly. in many cases, the positive weights vary exponentially, which means that limited number of components dominate the functional shape. Thus, we can ignore the trivial components for sampling [ Then we adopt SMC [53] to sample A(x). We have highly-calibrated proposal distribution of sparse Gaussian mixture, leading to efficient resampling from real A(r) (see Supplementary).
Variants of proposal distribution Although (8) has mathematical motivation, sometimes we wish to incorporate prior information not included in the above procedure. We propose two additional "biased"' proposal distributions. The first case is where we know both the maximum likelihood points and the likelihood's unimodality. This is typical in Bayesian inference because we can obtain (aab-)mxximum points via maximum posteriori probability (MAP) In this case, we know exploring around the perfect initial guess is optimal rather than unnecessarily exploring an uncertain region. Thus, we introduce the finitial guess believer (IGB) proposal distribution, JIGB(z) This is written as gIGB (x) = (1 ( + +r Ei= WiLIGB.N (x;; X;, W), where WG, = {0 if Y S 0, else 1), X; X. This means exploring only the vicinity of the observed data X. The second case is where we know the likelihood is multimodal. In this case, determining all peak positions is most beneficial. Thus more explorative distribution is preferred. As such, we introduce the (nccrtaintyy believer (UB) proposal distribution, 9UB (2)( This is written as gUB (I) A(z), meaning pure uncertainty sampling. To contrast the above two, we term the proposal distribution in Eq. (8) as integral variance reduction (IVR) [IVR)9IvR(() (IvR(z). (z).
5
Experiments
Given our new model BASQ, with three variants of the proposal distribution, IVR, IGB, and UB, we now test for speed against MC samplers and batch WSABI. We compared with three NS methods
17 According the review (16). MLFriends is the state-of-the- art NS sampler to date. The code is implemented based on 777 34 42 36) [O) (7), and code around kernel recombination 24 4 with additional modification. All experiments on synthetic datasets were averaged over 10 repeats, computed in parallel with multicore CPUs, without GPU for fair comparison X The posterior distribution of NS estimated via KDE with weighted samples 33 For maximum speed performance, batch size was optimised for each method in each dataset, in fairness to the competitors. Batch WSABI needs to optimise batch size to balance the likelihood query cost and sampling cost, because sampling cost increases rapidly with batch size, as shown in Figure a).. Therefore, it has an optimal batch size for faster convergence. By wall time cost, we exclude the cost of integrand evaluation; that is, the wall time cost is the overhead cost of batch evaluation. Details can be found in the Supplementary.
77
29
46 14 [5), coded with [72
Negattive elements in the matrices only exist in K (X,X) which can be drawn from the memory of the GP regression model without additional calculation. The number of positive components is half of the matrix on average. resulting in Onn /2). Then, taking the threshold via the inverse the recombination sample size N. the number components becomes Neump & n²? resulting in sampling complexity Onn /2 + neomp N).
Performed on MacBook Pro 2019, 2.4 GHz 8-Core Intel Core 19, 64 GB 2667 MHz DDR4
7

---

5.1 Synthetic problems
We evaluate all methods on three synthetic problems. The goal is to estimate the integral and posterior of the likelihood modelled with the highly multimodal functions. Prior was set to -oodimeeniooall multivariate normal distribution, with a zero mean vector, and covariance whose diagonal elements are 2. The optimised batch sizes for each methods are BASQ: 100, batch WSABI: 16. The synthetic likelihood functions are cheap (0.5 ms on average). This is advantageous setting for NS: Within 10 seconds, the batch WSABI, BASQ, and NS collected 32, 600, 23225 samples, respectively. As for the metrics, posterior estimation was tested with with|KKllbacke KattlKKlbbck-Leibblee (KL) juoon random 10,000 samples from true posterior. Evidence was evaluated with MAE, and ground truth was derived analytically.
Likelihood functions Branin-Hoo 49 is 8 modal function in two-dimensional space. Ackley (73
is a highly multimodal function with point symmetric periodical peaks in two-dimensional space. Oscillatory function (32) is highly multimodal function with reflection symmetric periodical peaks
of hoghly-eorrelated ellipsoids in two-dimensional space.
5.2 Real-world dataset
We consider three real-world applications with expensive likelihoods, which are simulator-based and hierarchical GP. We adopted the empirical metric due to no ground truth. For the posterior, we calculate the true conditional posterior distribution along the line passing through ground truth parameter points. Then, evaluate the posterior with root mean squared error (RMSE) against 50 test samples for each dimension. For integral, likelihoods makes the sample size per wall time amongst the methods no significant difference, whereas rejection sampling based NS dismiss more than 50% of queried samples. The batch sizes BASQ: 32, batch WSABI: (see Supplementary)
compare the model evidence itself. Expensive
Parameter estimation of the lithium-ion battery simulator The simulator is the SPMe [61]. estimating 3 simulation parameters given time-series voltage-current signal (the diffusivity of lithium-ion on the anode and cathode, and the experimental noise variance). Prior is modified to log multivariate normal distribution from Each query takes 1.2 seconds on average.
Parameter estimation of the phase-field model : The simulator is the phase-fiel field model [6] 10 esstmaatnn 4 simulation parameters at given time-series two-dimensional morphological image (temperature, interaction parameter. Bohr magneton coefficient, and gradient energy coefficient). Prior is log multivariate normal distribution. Each query takes 7.4 seconds on average.
Hyperparameter marginalisation of hierarchical GP model The hierarchical GP model was designed for analysing the large-scale battery time-series dataset from solar off-grid system field data For fast estimation of parameters in each GP, the recursive technique 74] is adopted. The task to marginalise 5 GP hyperparameters at given hyperprior, which is modified to log multivariate normal distribution from (3]. Each query takes 1.1 seconds average.
3
5.3
Results
We find BASQ consistently delivers strong empirical performance, as shown in Figure 2 On all benchmark problems, BASQ-IVR, IGB, or UB outperform baseline methods except in the battery simulator evidence estimation. The very low-dimensional and sharp unimodal nature this likelihood could be advantageous for biased greedy batch WSABI, as IGB superiority supports this viewpoint. This suggests that BASQ could be generally fast Bayesian solver as far as we investigated. In the multimodal setting of the synthetic dataset, BASQ-UB outperforms, whereas IVR does in simulator- based likelihoods. When comparing each proposal distribution, BASQ-IVR was the performant. Our results support the general use of IVR, or UB if the likelihood is known to be highly multimodal.
SPMe code used translated into Python from MATLAB II] 12 This open- code published under the BSD 3-clause Licence. See more information on [I]
"Coode used was from [54] 3 All rights of the code are reserved by the authors. Thus, we do not redistribute
the original code.
8

---

Hoo Fusction
Oscillatery
simulator
fiel model
Hiczashical
Lnn Kt. Enverggee
of
150
300
150
300
150
300
2250
4500
2250
4500
2250
4500
Log MAE of tneesni
30-
0.6
15
150
300
150 wall time (6)
300
150 time
300
2250 wall time (s)
2250 wall time (4)
4503
2250
4500
time (4)
wall time
BASQ-IVR
BASQJIGB
BASQ-UB
batch WSABI
(MulaNess)
(dyrumiic)
(MLFricmk)
Figure 2: Time in seconds vs. KL divergence for posterior, and MAE for evidence in the synthetic datasets. Time in seconds vs. RMSE for posterior and evidence itself in the real-world dataset.
Convergence analysis
We analysed the convergence over single iteration simplified version of BASQ, which assumes the BQ is modelled with vanilla BQ. without batch and hyperparameter updates. Note that the kernel K" on Rdd in this section refers the given covariance kernel of GP at each step. We diseuss the convergence of BASQ in one iteration. We consider the importance sampling: Let f be a probability density on Rd and be another density such that f = Ag with nonnegative function X. Let call such a pair, (f.g). density pair with weight A.
We approximate the kernel with Kn = Ei= (x))())) In general, we can apply the ker- nel recombination algorithm [59] 75 with the weighted sample (Wresy X obtain weighted
point set (w (Wqund: Xquad) of size satisfying w Pi (Xquad)
Pi (i 1,...,, 1) wToc. By modifying the kernel recombination algorithm, we can require "quaa Ky (Xqund) S. Wi
and w quad
wT
///2 (X.c), where k//2²(x (x) == vK(I,z)- Ko(x, x) [44]. We call such Wguad gquad) proper kernel recombination of (Wree: X.co) with Ko T We have the following guarantee (proved in Supplementary):
1/2
Theorem 1. Suppose J / K(K,, x)f dr < c0, GP(m, K), and we are given an (n - 1)- dimensional kernel Ko such that K₁ == = Ko also a kernel. Let (f.g) be density pair with weight X. Let Xiee be an N-point independent sample from and Wiee == X(Xre). Then, if (w yunds Xyuaa) a proper kernel recombination of (Wrerr X for Ko, satisfies
22 |2 ((Ki(e,))(())))
1/2
+ CK.
Ex,
V
var(Zz
Zas
Xquna]
+V
+
(12)
N
where 2 == = (xx) dr and CK.fa == ((() (x) If K(x, v) da
The above approximation has one source of randomness which stems from sampling N points X from g. One can also apply this estimate with random kernel and thereby introduce another source of randomness. particular, when we use the Nyström approximation for K, (that ensures Ki is kernel (4)). then one can show that K1 (I., x)) (x)) dæ can be bounded by
/ (K(xx dr < no, +
(n.K Ammxxx
Σ m=
+ +Op
(13)
VM
where on is the n-th eigenvalue of the integral operator L2( 3 h H (K(,))h())f()) dz, Kmax sup, K(x, x). However, note that unlike Eq. (12. this inequality only applies with high probability due to the randomness Ko; see Supplementary for details.
Note that the inequality constraint on the diagonal value here is only needed for theoretical guarantee, and skipping it does reduce the empirical performance 44

---

If, for example, K is Gaussian kernel on R and f is a Gaussian distribution, we have ((ex((-cn-/a) for some constant с > O (see Supplementary), So in (12) we also achieve an empirically exponential rate when N CK. RCHQ works well with moderate M in practise. Note that unlike the previous analysis (00 do not have to assume that the space is compact. 12
7 Discussion
We introduced a batch BQ approach, BASQ, capable of simultaneous calculation of both model evidence and posteriors. BASQ demonstrated faster convergence (in wall-clock time) on both synthetic and real-world datasets, when compared against existing BQ approaches and state-of-the-an NS. Further, mathematical analysis shows the possibility to converge exponentially-fas under natural assumptions. As the BASQ framework is general purpose. this can be applied to other active learning GP-based applications. such as Bayesian optimisation [52]. dynamic optimisation like control [26]. and probabilistic numerics like ODE solvers 45) Although it scales to the number of data seen in large-scale GP experiments, practical BASQ usage is limited to fewer than 16 dimensions (similar to many GGP-based algorithms). However, RCHQ is agnostic to the input space, allowing quadrature in manifold space, An appropriate latent variable warped GP modelling, such as GPLVM (58), could pave the way to high dimensional quadrature in future work. In addition, while WSABI modelling limits the kernel to a squared exponential kernel, RCHQ allows to adopt other kernels or priors without bespoke modelling BQ models. (See Supplementary). As for the mathematical proof, do not incorporate batch and hyperparameter updates, which should be addressed in future work. The generality of our theoretical guarantee with respect to kernel and distribution should be useful for extending the analysis to the whole algorithm.
Acknowledgments and Disclosure of Funding
We thank Saad Hamid and Xingchen Wan for the insightful discussion of Bayesian quadrature, and Antti Aitio and David Howey for fruitful discussion of Bayesian inference for battery analytics and for sharing his codes the single particle model with electrolyte dynamics, and hierarchical GPs. We would like to thank Binxin Ru, Michael Cohen, Samuel Daulton, Ondrej Bajgar, and anonymous reviewers for their helpful comments about improving the paper. Masaki Adachi was supported by the Clarendon Fund, the Oxford Kobe Scholarship, the Watanabe Foundation, the British Council Japan Association, and Toyota Motor Corporation. Satoshi Hayakawa was supported by the Clarendon Fund, the Oxford Kobe Scholarship, and the Toyota Riken Overseas Scholarship. Harald Oberhauser was supported by the DataSig Program [EP/S026347/1), the Hong Kong Innovation and Technology Commission (InnоHK Project CIMDA), and the Oxford-Man Institute. Martin Jørgensen was supported by the Carlsherg Foundation.
References
[1] Luigi Acerbi. Variational bayesian monte carlo. Advances in Neural Information Processing Systems, 31, 2018.
[2] Luigi Acerbi. Variational bayesian monte carlo with noisy likelihoods. Advances in Neural Information Processing Systems, 33:821 1-8222, 2020.
[3] A. Aitio and D. A. Howey. Predicting battery end of life from solar off-grid system field data using machine learning. Joule, 5(12):320432200 2021.
[4] A. Aitio, S. G. Marquis, P. Ascencio, and D. A. Howey. Bayesian parameter estimation applied to the Li-ion battery single particle model with electrolyte dynamics. IFAC-PPpperSOmLLine, 53(2):12477-12504 2020.
[5] F. Bach. On the equivalence between kernel quadrature rules and random feature expansions. Journal of Machine Learning Research, 18:714, 2017.
The rate N 1/2 for the expected integral variance on non-compact sets can also be found 35 [Corollary 2.8]. In the case of Gaussian integration distribution and the Gaussian kernel exponential rates of convergence of
the BQ integral variance can be found [8] [Theorem l]. [57 [Theorem 4.1), [56) [Theorem 3.2), [1 Theorems
2.5 and 2.10]
10

---

[6] F. Bach, S. Lacoste-Julien, and G. Obozinski. On the equivalence between herding and conditional gradient algorithms. In International Conference Machine Learning (ICML). 2012.
[7] M. Balandat, B. Karrer, D. Jiang. S. Daulton, B. Letham, A. G. Wilson, and E. Bakshy. BoTocch: framework for efficient Monte-Carlo Bayesian optimization. International Conference Neural Information Processing Systems (NeurIPS), 33:21524 21538, 2020.
[8] A. Belhadji. R. Bardenet, and P. Chainais. Kernel quadrature with DPPs. In International Conference on Neural Information Processing Systems (NeurIPS), 2019. URL: URL::httss: hpppss//proceedngg... www.imuurric.o cc/ Paper.
[9] C. H. Bennett. Efficient estimation of free energy differences from Monte Carlo data, Journal of Computational Physicx, 22:245 268, 1976.
[10] Berahas, J. Nocedal, and M. Takác. multi-batch I-hfgs method for machine learning. International Conference on Neural Information Processing Systems (NeurIPS), 29, 2016.
[11] A. M. Bizeray, Reniers, and D. A. Howey. Spectral Li-ion SPM. zenodo. 212178" 2016.
[12] A. M. Bizeray, S. Zhao, S. R. Duncan, and A. Howey. Lithium-ion battery thermal-elecmmhhmiical model-based state estimation using orthogonal collocation and a modified extended Kalman filter. Journal of Power Sources, 296:400-413 4120000-11.22 2015.
[13] F.-X Briol, C. J. Oates, M. Girolami, and M. A. Osborne. Frank-Wolfe Bayesian quadrature: Prob- abilistic integration with theoretical guarantees. In International Conference Neural Information Processing Systems (NewrIPS), 2015. URL: http:///pooeeedsss.eeeeeee hecour/pp.cc/ppppp/////////11111 011266660065400166711166666///62 Pap.aaaaaaaaaaaraaeeaaaaaaaa
[14] J. Buchner. A statistical test for nested sampling algorithms. Statistics and Computing, 26(1):383-392, 2016.
[15] J. Buchner. Collaborative nested sampling: Big data complex physical models. Publications of the Astronomical Society of the Pucific. 131(1004):10800) 2019.
[16] J. Buchner. Nested sampling methods. arXiv preprint arXiv:2101.09633 2021.
[17] J. Buchner. UltraNest-a robust, general purpose Bayesian inference engine. arXiv preprint arXiv:2101.090 2021.
[18] R. H. Byrd, P. La, J. Nocedal, and C. Zhu. A limited memory algorithm for bound constrained optimization. SIAM Journal on scientific computing, 16(5):1190-120) 1995.
[19] T. Campbell and B. Beronov. Sparse variational inferenee: Bayesian coresets from scratch. International Conference on Neural Information Processing Systems (NeurIPS), 32, 2019.
[20] T. Campbell and T. Broderick. Automated scalable Bayesian inference via Hilbert coresets. The Journal of Machine Learning Research, 20(1):551 -5551-888 2019.
[21] H. Chai, F. Ton, M. A. Osboree, and R. Gamett. Automated model selection with Bayesian quadrature. In International Conference on Machine Learning (ICML), volume 97, pages 931 940, 2019. URL: https::/prooeediing mlr httpp//eeaatiaaa..ererreee/////a
html
[22] H. Chai and R. Garnett. Improving quadrature for constrained integrands. In In IEEE 16th International Conference on Data Mining (ICDM), page 1107, 2016. URL: http:: http:ppp:///oo..r///...o 1109/ICDM, 2016 0144
[23] Y. Chen, M. Welling, and A. Smola. Super-samples from kernel herding. In International Conference on Uncertainty in Artificial Intelligence (UAI), 2010.
[24] F. Cosentino, H. Oberhauser, and A. Abate. A randomized algorithm to reduce the support of discrete measures. International Conference on Neural Information Processing Systems (NeurIPS), 33:15100- 10, 2020.
[25] K. Cranmer, J. Brehmer, and G. Louppe. The frontier of simulation-based inference. Proceedings of the National Academy of Sciences, 117449)33005-33066 2020.
[26] M. P. Deisenroth, C. E. Rasmussen, and J. Peters. Gaussian process dynamic programming, Neurocomput- ing. 12(7-9):1508-11 2009.

---

[27] X. Didelot, R. G. Everitt, A. M. Johansen, and D. J. Lawson. Likelihood- free estimation of model evidence. Bayesian analysis, 6(1):49-76, 2011.
[28] G. Evensen. Sequential data assimilation with nonlinear quasi geostrohhc model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research: Oceans, 10111)12043311016,, 1994.
[29] F. Feroz, P. Hobson, and M. Bridges. MultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics. Monthly Notices of the Royal Astronomical Society, 398(4):1601-1614 2009.
[30] C. Fröhlich, A. Gessmer, P. Hennig, B. Schölkopf, and G. Arvanitidis. Bayesian quadrature on Riemannian data manifolds. In International Conference on Machine Learning (ICMI., 2021. URI.: htIl::/ttps:/// http:://cc:
/t/vviruaa/2021//ssttrr//// 9655
[31] J. Gardner, G. Pleiss, K. Q. Weinberger, D. Bindel, and A. G. Wilson. GPyTorch: Blackbox matrix matrix Gaussian process inference with GPU acceleration. International Conference on Neural Information Processing Systems (NeurlPS), 31, 2018.
[32] A. Genz. Testing multidimensional integration routines. In In Proc. of international conference on Tools, methods and languages for scientific and engineering computation, pages 81-94, 1984.
[33] F. J. G. Gisbert. Weighted samples, kemel density estimators and convergence. Empirical Economics, 28:335-351, 2003. URL: https: //do::///do...oo org/ 10. 1107/10001100000111111
[34] GPy. GPy: A gaussian process framework in python. httt://ggttbb..com//hffffeeeeed 2012.
since
[35] D.M. Gräf. Efficient algorithums for the computation of optimal quadrature points on Riemannian manifolds. PhD thesis, Chemnitz University of Technology. 2013.
[36] D. Greenberg, M. Nonnenmacher, and J. Macke. Automatic posterior transformation for likelihood-frie inference. In International Conference on Machine Learning (ICML), pages 2404-2414, 2019.
[37] T. Gunter, M. A. Osbore, R. Gamett, P. Hennig, and S. J. Roberts. Sampling for inference in prob- abilistic models with fast Bayesian quadrature. In International Conference on Neural Information Processing Systems (NeuriPS), 2014, URL: ttattp::/pproceediggg Paper.
[38] Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2022. URL: http:tttt////.ggrrot coml
[39] M. U. Gutmann, J. Corander, et al. Bayesian optimization for likelihood-free inference of simulator-based statistical models. Journal of Machine Learning Research, 2016.
[40] N. Halko, P. G. Martinsson, and J. Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM review, 53(2):217 2011.
[41] S. Hamid, S. Schulze, M. A. Osborne, and S. J. Roberts. Marginalising over stationary kercels with Bayesian quadrature. International Conference on Artificial Intelligence and Statisties (AISTATS), 2021. URL: htt.htttt///proceediiiii
mlr. pwtpps//ee//////eeeeeeppepeateee http1////aaamaaa///aaaa htt.dd/2aaaaamm22aapap.mmm
[42] C. R. Harris, J. Millman. S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus. S. Hoyer, M. H. Kerkwijk, M. Brett, A. Haldane, J. F. Río, M. Wiebe, P. Peterson, P. Gérard- Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. F. Oliphant. Array programming with NumPy. Nature, 585(7825):357-366 2020.
101::001038/6000 /41///111866-6 1w111110//2110000-----11 2649-2
[43] W. K. Hastings. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1):97- 109, 1970. URL: https:///o1:1000 wtt.//199//bomeeeee bwo.et/57. 1.97
[44] S. Hayakawa, H. Oberhauser, and T. Lyons. Positively weighted kernel quadrature via subsampling. In International Conference on Neural Information Processing Systems (NeurIPS), 2022. [oi: 48550/ arXiv 2107 095597
[45] P. Hennig, M. A. Osbome, and M. Girolami. Probabilistic numerics and uncertainty in computations. Pro- ceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2119)220151444 2015.
12

---

[46] E. Higson, W. Handley, M. Hobson, and A. Lasenby. Dynamic nested sampling: an improved algorithm for parameter estimation and evidence calculation. Stat. Comput., 29:891-913, 2019. URL: https //doi. ottg//001007///1222-0033
[47] M. D. Hoffman and A. Gelman. The No-U-Tum sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15()):1993-6623 2014.
[48] F. Huszár and D. Duvenaud. Opiinally-weighted herding is Bayesian quadrature. In International Conference on Uncertainty in Artificial Intelligence (UAI). 2012.
[49] D. R. Jones. A taxonomy of global optimization methods based on response surfaces. Journal of global optimization, 21: 21:345-388 383, 2001.
[50] M. Kanagawa and P. Hennig. Convergence guarantees for adaptive Bayesian quadrature methods. Advances in Neural Information Processing Systems, 32, 2019.
[51] T. Karvonen, C. Oates, and M. Girolami. Integration in reproducing kernel hilbert spaces of gaussian kernels. Mathematics of Computation, 906331):2209 2233, 2021.
[52] T. Kathuria, A. Deshpande, and P. Kohli. Batched Gaussian process bandit optimization via determinantal point processes. Advances in Neural Information Processing Systems, 29, 2016.
[53] G. Kitagawa. A Monte Carlo filtering and smoothing method for non- Gaussian nonlinear state space models. In Proceedings of the 2nd U.S.JJppn Joint Seminar on Statistical Time Series Analysis, page 110, 1993. (tt::httss:////a /ww. ism. ac. jp/ kitagawa/ Batpaa//933__U_ 11093__S__00 S__JUS-Jaa Japan, pdf
[54] T. Koyama. Computational Engineering for Materials Design, Computational Microstructures, New and Expanded Edition: Microstructure Formation Analysis by the Phase-Field Method. Uchida Rokakuho, Tokyo, Japan, 2019.
[55] S. Kumar, M. Mohri, and A. Talwalkar. Sampling methods for the Nyström method. The Journal of Machine Learning Research, (3(1)991111000 2012.
[56] F. Kuo, Sloan, and H. Wožniakowski. Multivariate integration for analytic functions with gaussian kernels. Mathematics of Computation, 86(304):829-855 2017.
[57] F. Y. Kuo and H. Wožniakowski. Gauss-Hermite quadratures for functions from Hilbert spaces with Gaussian reproducing kernels. BIT Numerical Mathematics, 52(2):425 436, 2012.
[58] N. Lawrence. Gaussian process latent variable models for visualisation of high dimensional data. Advances in neural information processing systems, 16, 2003.
[59] C. Litterer and T. Lyons. High order recombination and application cubature on Wiener space. The Annals of Applied Probability, 1301-1337,, 2012.
[60] D. Manousakas. Z. Xu, C. Mascolo. and T. Campbell. Bayesian pseudocoresets. International Conference on Neural Information Processing Systems (NeurIPS). 33:14950-144 14960,,1100011100 2020.
[61] S. G. Marquis, V. Sulzer, R. Timms, C. P. Please, and S. J. Chapman. An asymptotic derivation of single particle model with electrolyte. Journal the Electrochemical Society, 166(15))436933437766 2019. URL: htttt:////:/d ieoooi.orgeo org/10. 1119//.1144111 0341915jes
[62] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equation of state calculations by fast computing machines. The journal of chemical physics, 21(6):108--1092 1953.
[63] P. Mukherjee, D. Parkinson, and A. R. Liddle. A nested sampling algorithm for cosmological model selection. The Astrophysical Journal, 638(2)LL51, 2006.
[64] R.M. Neal. Annealed importance sampling. Statistics and Computing, 11(2):125-139, 2001.
[65] A. O'Hagan. Bayes-Heemmite quadrature. Journal of Statistical Planning and Inference, 29:245, 1991. URL: httpss://oi. orto:io..rgg1111 10. 1016/0378- 3758 (91)9002-10
[66] A. O'Hagan. Bayesian quadrature with non normal approximating functions. Statistics and Computing, 8:365, 1998.
[67] A. Osborne, D. Duvenaud, R. Garnett, Rasmussen,
Roberts, and Ghahramani. Active learning of model evidence using Bayesian quadrature. In International Conference on Neural Information Processing Sysiems (New/IPS), 2012. URL: https: http:://pprccedinng: neur neurips. www..pp..c/pppprr/////// 6854430004955aaadddaaaaa36666666 Paper. pdf
13

---

[68] R. Pinsler. J. Gordon, E. Nalisnick, and J. M. Hernández-Lobato. Bayesian batch active learning as sparse subset approximation. International Conference on Neural Information Processing Systems (NeurIPS), 32, 2019.
[69] C. E. Rasmussen and Z. Chahramani, Bayesian Monte Carlo. In International Conferencee Neural Infor- mation Processing Systems (NeurIPS), 2003, URL: httpt://tpmmmgg... eng. cam cam.cc.................... uaaaaaaaccuuuoouuuaappppsssss RasCha03, pdf
[70] S. Roberts, Osborne, Ebden, Reece, N. Gibson, Aigrain. Gaussian processes time-series modelling. Phil. Trans. R. Soc. A., 301:201112110 10550, 2013.
[71] Skilling. Nested sampling for general Bayesian computation. Bayesian Analysis, 1(4):833 859, 2006. URL: httpp:://oi.orgg// org/10. 1214/06- 2214/00-BA177
[72] S. Speagle. dynesty: dynamic nested sampling package estimating Bayesian posteriors and
evidences. Monthly Notices of the Royal Astronomical Society, 493(3):3132-31333 -158,18 2020. URL: https
1tw1/000m3///////// maras/ staa278
[73] S. Surjanovic and D. Bingham. Virtual library of simulation experiments: Test functions and datasets. Retrieved May 12, 2022, from hommhttt:///ww/ sfu. ca/" ano
[74] S. Särkkä and J. Hartikainen. Infinite dimensionll Kalman filtering approach to spatio-temporal Gaussian process regression. In Artificial Intelligence and Statistics, pages 993-1001, 2012.
[75] M. Tchernychova. Carathéodory cubature measures. PhD thesis, University of Oxford, 2015.
[76] P. Virtanen, R. Gommers,
Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, Peterson. W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, 1. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perk- told, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-211 2020, doi: 000.:001000/////00000000 019--6666--0011-----2--
[77] E. WagstafT. Bayesian
since 2019.
[78] E. Wagstaff, S. Hamid, and M. A. Osborne. Batch selection for parallelisation of Bayesian quadrature. arXiv preprint, 2018. URL:
[79] J. Wilson, F. Hutter, and M. Deisenroth. Maximizing acquisition functions for Bayesian optimization. International Conference on Neural Information Processing Systems (NeurIPS), 31, 2018
14

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contribu- tions and scope? [Yes]
(b) Did you describe the limitations your work? [Yes] section
(c) Did you discuss any potential negative societal impacts of your work? [N/A]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] section and Supplementary
(b) Did you include complete proofs of all theoretical results? [Yes] see section Blnd Supplementary
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as URL)? [Yes] supplemental material
(h) Did you specify all the training details (e.g., data splits, hyperparamcters, how they were chosen)?
[Yes] see section S and Supplementary
(c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [Yes] see section [
(d) Did you include total amount of compute internal cluster, or cloud provider)? [Yes] see section [
type resources used (e.g.. type GPUs,
4. If you using existing assets (e.g.. code, data, models) curating/releasing assets...
(a) your work uses existing assets, did you cite the creators? [Yes] see section S
(b) Did you mention license of the assets? [Yes] see section S
(c) Did you include any new assets either in the supplemental material or as URL? [Yes]
(d) Did you discuss whether and how consent was obtained from people whose data you're us-
ing/curating? [Yes] see section [
(c) Did you discuss whether the data you using/curating contains personally identifiable informa- tion or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable?
[N/A]
(b) Did you deseribe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
(C) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
15


Reviews:
[{"review_id": "eZlSovWCTlZ", "reviewer": "Reviewer_CnRy", "summary": "This paper combines the batch-version of WSABI Bayesian quadrature method [30], which uses a square-root transform to encode non-negativity of the integrand function (in this case a likelihood function), with the kernel recombination method for sampling of quadrature points from [37] and the Nyström approximation of the kernel.", "questions": "ON THEORY:\n\n1. Theorem 1, if correct (see the next two points), guarantees exponential convergence over a single iteration in the sense that the BQ integral variance decays exponentially fast. This is a consequence of the use of the extremely strong squared exponential GP prior. But no discussion or elaboration on the practical meaning of this statement can be found in the paper. Does the integral estimate produced by the method converge exponentially to the integral of the true likelihood function $\\ell_\\text{true}$? If so, what assumptions need to be imposed on $\\ell_\\text{true}$? To include an an unconditional statement that the proposed method \"possesses a provably-exponential convergence rate\" in the abstract is quite misleading.\n2. It is claimed after Theorem 1 that the convergence of the variance is exponential if the kernel has exponentially decaying Mercer eigenvalues. But what about the second term on the right-hand side of Equation (12)? I do not see how this term is to decay exponentially.\n3. Equation (12) in the supplement uses the property in Equation (6) of the main paper which says that BQ variance equals the minimal squared worst-case error in the RKHS of $K$. But this characterization only applies to the variance of the vanilla BQ, not the variance of WSABI. Therefore it seems to me that Theorem 1 does not apply to BASQ as proposed in Section 4 but to a variant which uses vanilla BQ instead of WSABI. Is that correct?\n\nMINOR COMMENTS:\n\n- p2, lines 38-39: \"Choosing an appropriate kernel and prior makes the integrals analytical.\" I suppose this refers to the kernel mean embeddings etc. that need to be computed in BQ. However, how the sentence in inserted in the introduction makes it sound like the likelihood integrals that are being computed with BQ are available analytically. There is a similar issue on p2 line 54, \"In particular, the integral becomes analytic...\".\n- Section 2: In this section, BQ is presented as a method specifically for computation of quantities related to posteriors (\"The objective of BQ is to compute the following integrals and obtain surrogate function of posterior:\") while in reality the integrals computed by BQ do not have to have anything to do with posteriors or Bayesian inference. It would be good to clarify this, e.g. saying that this is just what BQ is used for in this particular paper.\n- p2, line 45: \"following integrals\" There is only one integral in Equation (1).\n- p2, lines 60-61: \"WSABI ... sets ...\" It is unclear what this sentence is trying to say. What is being \"set\" to something?\n- Section 2, in particular \"Vanilla Bayesian quadrature\": One unacquianted with GPs would probably have a very hard time trying to understand what happens in this section. For example, it is not said that the likelihood function is being modelled as GP(0, K) and that Equation (2) is the resulting posterior GP. The section also contains a lot of undefined notation.\n- p2, line 62: Why is the coefficient 0.8 used in the definition of alpha?\n- Section 2: I do not see why different notation for point sets is used in \"Vanilla Bayesian quadrature\" and \"Kernel quadrature in general\".\n- p5, lines 160-1: A word must be missing (\"selects sparse samples than batch WSABI\"). Or \"sparser\"?\n- Eq (9): dx in the denominator looks different than elsewhere in the paper.\n- p6, line 178: Isn't $C_y^L$ the warped kernel?\n- Lines 243-4: \"that can earn per wall time\" (?)\n- Lines 246-7: \"Still, the time becomes much longer when it increases for accuracy.\" (?)\n- Line 248: \"Yet, the original data size is roughly 1,000 times.\" (?)\n- Line 268: \"explorative UB proposal distribution exceeds\" (?)\n- Line 299: \"It also works well with a moderate M in practise.\" What works?\n- Ref [9]: \"F. X. Briol\" should be \"F.-X. Briol\"\n- Ref [58]: \"Sarkka\" should be \"Särkkä\"\n- Supplement, line 15 has \\propto while Line 21 has \\approx.\n", "limitations": "One of the main limitations of the proposed method (which is shared by all methods that utilize WSABI-BQ) is that it really only applies to $\\pi$ being a mixture of Gaussian densities and $K$ the squared exponential kernel because in other cases the integrals of the functions in Equation (4) have not been computed analytically. This limitation should be stated more openly as at the moment the only allusion to it in the abstract and introduction is the sentence \"Choosing an appropriate kernel and prior makes the integrals analytical.\" on p2, lines 38-9.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "I do not find this paper particularly exciting: the proposed BASQ method is a conceptually rather straightforward combination of a few existing methods. As far as I can see, the paper contains no interesting new ideas of any generality. The main contribution is that the proposed method outperforms some existing approaches in a number of examples. While the examples do seem impressive, I am not sure if this is sufficient and find the paper borderline, though tend to recommend rejection. However, I have no serious objections to acceptance if the consensus is that a contribution like this is sufficient for this conference.\n\nI have certain concerns with the convergence theory in Section 6 and how it is advertised, which I detail in the next section of this review. Whether or not these concerns are valid has little impact on my assessment of this paper because I do not expect a paper like this to have any useful convergence theory.\n\nThe paper is generally reasonably well written and effortless to read. However, there are many unclear statements and undefined notation, examples of which I give in the next section. In general, I find that it would be helpful to replace many of the verbal descriptions with mathematics.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "U6xXDbXe6eF", "reviewer": "Reviewer_8yqu", "summary": "This paper makes a contribution to speeding up Bayesian Quadrature (BQ), a popular method for computing the posterior and evidence. The proposed method builds on the weighted kernel quadrature method of Hayakawa et al. and uses the (proven) advantage of kernel recombination to improve sample selection for worst case error reduction and hence improve the accuracy of the quantities estimated by standard BQ methodology. It is shown that the method performs favourably compared to the state of the art WSABI.\n", "questions": "To me, the weakest point of the paper is a somewhat unclear treatment of the relationship between KQ and BQ, w.r.t reference [40] and the explicit contribution made. It is said that \"the point configuration in KQ with a small worst-case error gives a good way to select points to reduce the integral variance in Bayesian quadrature.\" Now, to the best of my knowledge, (naive) kernel herding is limited by the computation of certain integrals required when sequentially selecting points. As far as I understand, the reference [37] is a practical way of circumventing this limitation? It would be nice to expand a bit more on this crucial point upon which the paper is ultimately based. It would also be good to mention some potential areas, if any, where the proposed method performs worse than existing ones. \n", "limitations": "Yes.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "In my view, this paper is an obvious improvement on the WSABI methodology with few drawbacks (other than those intrinsic to BQ utilizing Gaussian processes in general, such an upper bound of applicability to 10-15 dimensions) and a very natural alternative to it that I can see being used in practice. It is well-written, complete and easy to read.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "MXFSYjTO5ys", "reviewer": "Reviewer_Aavt", "summary": "This paper proposes a novel quadrature procedure that combines elements of warped Bayesian quadrature and random convex hull quadrature. The resulting method enjoys benefits from both component aspects: the computational tractability of RCHQ and the integrand-modeling/active learning techniques of square-root warped BQ (WSABI). The authors empirically evaluate their method against a batch version of WSABI and various implements of nested sampling; they also provide some theoretical analysis on the convergence rate of their method. ", "questions": "- Did the authors consider alternative low-rank approximation methods for the kernel matrix? Or is there something fundamental about the Nyström approximation that would invalidate this method if an alternative were used?\n\n- Similarly, is there something about the WSABI-L approximation that makes it particularly well suited for the method proposed by the authors? Did they also experiment with the WSABI-M approximation or other previously-proposed transforms?\n\n- In the empirical evaluation, are confidence intervals omitted for the real-word experiments or are they just imperceptibly small? In general, are there any meaningful takeaways from the fact that batch WSABI performs identically to the BASQ variants on the real-world benchmarks in terms of log MAE? ", "limitations": "The authors' discussion of limitations in Section 7 is sufficiently detailed in my opinion.", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall, I think this paper presents some interesting ideas and represents a promising line of inquiry. However, in its current state, I believe there are some issues (detailed below) which keep this version from being a strong candidate for acceptance. \n\n1. This work makes a litany of smaller, practical contributions which are certainly novel, although the core approach is a combination of existing ideas. However, I do think the methods described in Section 4.2 are sufficient in terms of getting this work above the bar in terms of originality. \n\n2. My primary concern with this work is the level of analysis performed: given that this work combines many complex components, hyperparameters, design choices, heuristics, etc... an ablation study/sensitivity analysis of some sort seems warranted to assess which aspects are providing what benefit or if the sum is greater than the whole of its parts. Some notable components that I feel warrant a deeper dive include the hyperparameters M, N and r; the GP model on the integrand; sampling methodology (some of these are expanded on in my questions below). \n\nRelatedly the subsample proposal distributions are somewhat poorly motivated heuristics: it is not clear to me why the uncertainty sampling acquisition function need be mixed with the prior to ensure global optimization. Calling this mixing the \"integral variance reduction\" distribution is also a bit strange: there is a well-known connection between DPPs and kernel quadrature and the correctly defined DPP would in some sense be the actual integral variance reducing distribution (although I suspect that sampling from this distribution would likely be computationally intractable for this purpose). \n\n3. Overall, I think this manuscript is reasonably clear/well-written: there are a few minor grammatical errors but they did not hinder my ability to comprehend the material being presented. As a minor nitpick, I would have liked to see more main paper real estate devoted to the figures, both in terms of actual size and interpretation: Figure 1(e) and (f) are a bit difficult to parse and the text doesn't provide much insight but ultimately, they did add to my appreciation of the proposed method. Similarly, the bottom-right quadrant of Figure 2 provides little compelling evidence in favor of their method. \n\n4. My impression of this work is that it is of low to middling impact, given the somewhat niche use case of low-dimensional, non-negative batch quadrature; this is somewhat validated by the limited existing work/baselines to compare against in this space (of course, the lack of existing work is not being held against the authors, I bring it up merely as a measure of potential significance). Furthermore, related to my notes about quality, the myriad design choices represent a nontrivial barrier for other practitioners/researchers looking to experiment with/improve upon the work presented here. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept