Quadrature

2020

1. A locally adaptive Bayesian cubature method

   Fisher, Matthew, Oates, Chris, Powell, Catherine, and Teckentrup, Aretha

   In International Conference on Artificial Intelligence and Statistics 2020

2018

1. Bayesian Quadrature for Multiple Related Integrals

   Xi, X., Briol, F.-X., and Girolami, M.

   ArXiv e-prints 2018

   Abs link

   Bayesian probabilistic numerical methods are a set of tools providing posterior distributions on the output of numerical methods. The use of these methods is usually motivated by the fact that they can represent our uncertainty due to incomplete/finite information about the continuous mathematical problem being approximated. In this paper, we demonstrate that this paradigm can provide additional advantages, such as the possibility of transferring information between several numerical methods. This allows users to represent uncertainty in a more faithfully manner and, as a by-product, provide increased numerical efficiency. We propose the first such numerical method by extending the well-known Bayesian quadrature algorithm to the case where we are interested in computing the integral of several related functions. We then demonstrate its efficiency in the context of multi-fidelity models for complex engineering systems, as well as a problem of global illumination in computer graphics.

2017

1. Convergence Analysis of Deterministic Kernel-Based Quadrature Rules in Misspecified Settings

   Kanagawa, Motonobu, Sriperumbudur, Bharath K., and Fukumizu, Kenji

   arXiv:1709.00147 [cs, math, stat] 2017

   Abs

   This paper presents convergence analysis of kernel-based quadrature rules in misspecified settings, focusing on deterministic quadrature in Sobolev spaces. In particular, we deal with misspecified settings where a test integrand is less smooth than a Sobolev RKHS based on which a quadrature rule is constructed. We provide convergence guarantees based on two different assumptions on a quadrature rule: one on quadrature weights, and the other on design points. More precisely, we show that convergence rates can be derived (i) if the sum of absolute weights remains constant (or does not increase quickly), or (ii) if the minimum distance between distance design points does not decrease very quickly. As a consequence of the latter result, we derive a rate of convergence for Bayesian quadrature in misspecified settings. We reveal a condition on design points to make Bayesian quadrature robust to misspecification, and show that, under this condition, it may adaptively achieve the optimal rate of convergence in the Sobolev space of a lesser order (i.e., of the unknown smoothness of a test integrand), under a slightly stronger regularity condition on the integrand.

2. Optimal Monte Carlo integration on closed manifolds

   Ehler, M., Graef, M., and Oates, C. J.

   ArXiv e-prints 2017

   Abs link

   The worst case integration error in reproducing kernel Hilbert spaces of standard Monte Carlo methods with n random points decays as 1/sqrt(n). However, re-weighting of random points can sometimes be used to improve the convergence order. This paper contributes general theoretical results for Sobolev spaces on closed Riemannian manifolds, where we verify that such re-weighting yields optimal approximation rates up to a logarithmic factor. We also provide numerical experiments matching the theoretical results for some Sobolev spaces on the unit sphere and on the Grassmannian manifold. Our theoretical findings also cover function spaces on more general sets such as the unit ball, the cube, and the simplex.

3. On the Sampling Problem for Kernel Quadrature

   Briol, François-Xavier, Oates, Chris J., Cockayne, Jon, Chen, Wilson Ye, and Girolami, Mark

   In Thirty-fourth International Conference on Machine Learning (ICML 2017) 2017

   Abs

   The standard Kernel Quadrature method for numerical integration with random point sets (also called Bayesian Monte Carlo) is known to converge in root mean square error at a rate determined by the ratio \s/d where \s and \d encode the smoothness and dimension of the integrand. However, an empirical investigation reveals that the rate constant \C is highly sensitive to the distribution of the random points. In contrast to standard Monte Carlo integration, for which optimal importance sampling is well-understood, the sampling distribution that minimises \C for Kernel Quadrature does not admit a closed form. This paper argues that the practical choice of sampling distribution is an important open problem. One solution is considered; a novel automatic approach based on adaptive tempering and sequential Monte Carlo. Empirical results demonstrate a dramatic reduction in integration error of up to 4 orders of magnitude can be achieved with the proposed method.

4. Fully symmetric kernel quadrature

   Karvonen, Toni, and Särkkä, Simo

   arXiv:1703.06359 [cs, math, stat] 2017

   Abs

   Kernel quadratures and other kernel-based approximation methods typically suffer from prohibitive cubic time and quadratic space complexity in the number of function evaluations. The problem arises because a system of linear equations needs to be solved. In this article we show that the weights of a kernel quadrature rule can be computed efficiently and exactly for up to tens of millions of nodes if the kernel, integration domain, and measure are fully symmetric and the node set is a union of fully symmetric sets. This is based on the observations that in such a setting there are only as many distinct weights as there are fully symmetric sets and that these weights can be solved from a linear system of equations constructed out of row sums of certain submatrices of the full kernel matrix. We present several numerical examples that show feasibility, both for a large number of nodes and in high dimensions, of the developed fully symmetric kernel quadrature rules. Most prominent of the fully symmetric kernel quadrature rules we propose are those that use sparse grids.

5. Classical quadrature rules via Gaussian processes

   Karvonen, Toni, and Särkkä, Simo

   In 2017 IEEE 27th International Workshop on Machine Learning for Signal Processing (MLSP) 2017

2016

1. Probabilistic Models for Integration Error in the Assessment of Functional Cardiac Models

   Oates, Chris J., Niederer, Steven, Lee, Angela, Briol, François-Xavier, and Girolami, Mark

   arXiv:1606.06841 [stat] 2016

   Abs

   This paper studies the numerical computation of integrals, representing estimates or predictions, over the output \f(x) of a computational model with respect to a distribution \p(}mathrm{d}x) over uncertain inputs \x to the model. For the functional cardiac models that motivate this work, neither \f nor \p possess a closed-form expression and evaluation of either requires {}approx 100 CPU hours, precluding standard numerical integration methods. Our proposal is to treat integration as an estimation problem, with a joint model for both the a priori unknown function \f and the a priori unknown distribution \p\. The result is a posterior distribution over the integral that explicitly accounts for dual sources of numerical approximation error due to a severely limited computational budget. This construction is applied to account, in a statistically principled manner, for the impact of numerical errors that (at present) are confounding factors in functional cardiac model assessment.

2. Bayesian Quadrature Variance in Sigma-Point Filtering

   Prüher, Jakub, and Šimandl, Miroslav

   2016

   Abs link Code

   Sigma-point filters are algorithms for recursive state estimation of the stochastic dynamic systems from noisy measurements, which rely on moment integral approximations by means of various numerical quadrature rules. In practice, however, it is hardly guaranteed that the system dynamics or measurement functions will meet the restrictive requirements of the classical quadratures, which inevitably results in approximation errors that are not accounted for in the current state-of-the-art sigma-point filters. We propose a method for incorporating information about the integral approximation error into the filtering algorithm by exploiting features of a Bayesian quadrature—an alternative to classical numerical integration. This is enabled by the fact that the Bayesian quadrature treats numerical integration as a statistical estimation problem, where the posterior distribution over the values of the integral serves as a model of numerical error. We demonstrate superior performance of the proposed filters on a simple univariate benchmarking example.

3. Convergence guarantees for kernel-based quadrature rules in misspecified settings

   Kanagawa, Motonobu, Sriperumbudur, Bharath K., and Fukumizu, Kenji

   2016

   Abs link

   Kernel-based quadrature rules are becoming important in machine learning and statistics, as they achieve super-\sqrtn convergence rates in numerical integration, and thus provide alternatives to Monte Carlo integration in challenging settings where integrands are expensive to evaluate or where integrands are high dimensional. These rules are based on the assumption that the integrand has a certain degree of smoothness, which is expressed as that the integrand belongs to a certain reproducing kernel Hilbert space (RKHS). However, this assumption can be violated in practice (e.g., when the integrand is a black box function), and no general theory has been established for the convergence of kernel quadratures in such misspecified settings. Our contribution is in proving that kernel quadratures can be consistent even when the integrand does not belong to the assumed RKHS, i.e., when the integrand is less smooth than assumed. Specifically, we derive convergence rates that depend on the (unknown) lesser smoothness of the integrand, where the degree of smoothness is expressed via powers of RKHSs or via Sobolev spaces.

2015

1. On the relation between Gaussian process quadratures and sigma-point methods

   Särkkä, S., Hartikainen, J., Svensson, L., and Sandblom, F.

   arXiv preprint stat.ME 1504.05994 2015

   Abs

   This article is concerned with Gaussian process quadratures, which are numerical integration methods based on Gaussian process regression methods, and sigma-point methods, which are used in advanced non-linear Kalman filtering and smoothing algorithms. We show that many sigma-point methods can be interpreted as Gaussian quadrature based methods with suitably selected covariance functions. We show that this interpretation also extends to more general multivariate Gauss–Hermite integration methods and related spherical cubature rules. Additionally, we discuss different criteria for selecting the sigma-point locations: exactness for multivariate polynomials up to a given order, minimum average error, and quasi-random point sets. The performance of the different methods is tested in numerical experiments.

2. Kernel-Based Just-In-Time Learning for Passing Expectation Propagation Messages

   Jitkrittum, W., Gretton, A., Heess, N., Eslami, S. M. A., Lakshminarayanan, B., Sejdinovic, D., and Szabó, Z.

   In Uncertainty in Artificial Intelligence (UAI) 31 2015

   Abs

   We propose an efficient nonparametric strategy for learning a message operator in expectation propagation (EP), which takes as input the set of incoming messages to a factor node, and produces an outgoing message as output. This learned operator replaces the multivariate integral required in classical EP, which may not have an analytic expression. We use kernel-based regression, which is trained on a set of probability distributions representing the incoming messages, and the associated outgoing messages. The kernel approach has two main advantages: first, it is fast, as it is implemented using a novel two-layer random feature representation of the input message distributions; second, it has principled uncertainty estimates, and can be cheaply updated online, meaning it can request and incorporate new training data when it encounters inputs on which it is uncertain. In experiments, our approach is able to solve learning problems where a single message operator is required for multiple, substantially different data sets (logistic regression for a variety of classification problems), where it is essential to accurately assess uncertainty and to efficiently and robustly update the message operator.

3. Frank-Wolfe Bayesian Quadrature: Probabilistic Integration with Theoretical Guarantees

   Briol, François-Xavier, Oates, Chris J., Girolami, Mark, and Osborne, Michael A.

   In Advances in Neural Information Processing Systems (NIPS) 2015

   Abs

   There is renewed interest in formulating integration as an 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 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 exponential and posterior contraction rates are proven to be superexponential. 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 model choice problem in cellular biology.

4. On the Equivalence between Quadrature Rules and Random Features

   Bach, Francis

   arXiv preprint arXiv:1502.06800 2015

   Abs

   We show that kernel-based quadrature rules for computing in tegrals can be seen as a special case of random feature expansions for positive definite kernels, for a particular decomposition that always exists for such kernels. We provide a theoretical analysis of the number of required samples for a given approximation error, leading to both upper and lower bounds that are based solely on the eigenvalues of the associated integral operator and match up to logarithmic terms. In particular, we show that the upper bound may be obtained from independent an d identically distributed samples from a specific non-uniform distribution, while the lower bo und if valid for any set of points. Applying our results to kernel-based quadrature, while our results are fairly general, we recover known upper and lower bounds for the special cases of Sobolev spaces. Moreover, our results extend to the more general problem of full function approxim ations (beyond simply computing an integral), with results in L2- and L∞-norm that match known results for special cases. Applying our results to random features, we show an improvement of the number of random features needed to preserve the generalization guarantees for learning with Lipshitz-continuous losses.

5. Probabilistic Integration: A Role for Statisticians in Numerical Analysis?

   Briol, François-Xavier, Oates, Chris J., Girolami, Mark, Osborne, Michael A., and Sejdinovic, Dino

   arXiv:1512.00933 [cs, math, stat] 2015

   Abs

   A research frontier has emerged in scientific computation, founded on the principle that numerical error entails epistemic uncertainty that ought to be subjected to statistical analysis. This viewpoint raises several interesting challenges, including the design of statistical methods that enable the coherent propagation of probabilities through a (possibly deterministic) computational pipeline. This paper examines thoroughly the case for probabilistic numerical methods in statistical computation and a specific case study is presented for Markov chain and Quasi Monte Carlo methods. A probabilistic integrator is equipped with a full distribution over its output, providing a measure of epistemic uncertainty that is shown to be statistically valid at finite computational levels, as well as in asymptotic regimes. The approach is motivated by expensive integration problems, where, as in krigging, one is willing to expend, at worst, cubic computational effort in order to gain uncertainty quantification. There, probabilistic integrators enjoy the "best of both worlds", leveraging the sampling efficiency of Monte Carlo methods whilst providing a principled route to assessment of the impact of numerical error on scientific conclusions. Several substantial applications are provided for illustration and critical evaluation, including examples from statistical modelling, computer graphics and uncertainty quantification in oil reservoir modelling.

2014

1. Gaussian Process Quadratures in Nonlinear Sigma-Point Filtering and Smoothing

   Särkkä, Simo, Hartikainen, Jouni, Svensson, Lennart, and Sandblom, Fredrik

   In FUSION 2014

   Abs

   This paper is concerned with the use of Gaussian process regression based quadrature rules in the context of sigma- point-based nonlinear Kalman filtering and smoothing. We show how Gaussian process (i.e., Bayesian or Bayes–Hermite) quadratures can be used for numerical solving of the Gaussian integrals arising in the filters and smoothers. An interesting additional result is that with suitable selections of Hermite polynomial covariance functions the Gaussian process quadratures can be reduced to unscented transforms, spherical cubature rules, and to Gauss- Hermite rules previously proposed for approximate nonlinear Kalman filter and smoothing. Finally, the performance of the Gaussian process quadratures in this context is evaluated with numerical simulations.

2. Sampling for Inference in Probabilistic Models with Fast Bayesian Quadrature

   Gunter, Tom, Osborne, Michael A., Garnett, Roman, Hennig, Philipp, and Roberts, Stephen

   In Advances in Neural Information Processing Systems (NIPS) 2014

   Abs Code

   We propose a novel sampling framework for inference in probabilistic models: an active learning approach that converges more quickly (in wall-clock time) than Markov chain Monte Carlo (MCMC) benchmarks. The central challenge in probabilistic inference is numerical integration, to average over ensembles of models or unknown (hyper-)parameters (for example to compute marginal likelihood or a partition function). MCMC has provided approaches to numerical integration that deliver state-of-the-art inference, but can suffer from sample inefficiency and poor convergence diagnostics. Bayesian quadrature techniques offer a model-based solution to such problems, but their uptake has been hindered by prohibitive computation costs. We introduce a warped model for probabilistic integrands (likelihoods) that are known to be non-negative, permitting a cheap active learning scheme to optimally select sample locations. Our algorithm is demonstrated to offer faster convergence (in seconds) relative to simple Monte Carlo and annealed importance sampling on both synthetic and real-world examples.

3. Just-In-Time Learning for Fast and Flexible Inference

   Eslami, S. M. Ali, Tarlow, Daniel, Kohli, Pushmeet, and Winn, John

   In Advances in Neural Information Processing Systems (NIPS) 27 2014

   Abs

   Much of research in machine learning has centered around the search for inference algorithms that are both general-purpose and efficient. The problem is extremely challenging and general inference remains computationally expensive. We seek to address this problem by observing that in most specific applications of a model, we typically only need to perform a small subset of all possible inference computations. Motivated by this, we introduce just-in-time learning, a framework for fast and flexible inference that learns to speed up inference at run-time. Through a series of experiments, we show how this framework can allow us to combine the flexibility of sampling with the efficiency of deterministic message-passing.

2012

1. Bayesian quadrature for ratios

   Osborne, M.A., Garnett, R., Roberts, S.J., Hart, C., Aigrain, S., and Gibson, N.

   In International Conference on Artificial Intelligence and Statistics 2012

   Abs

   We describe a novel approach to quadrature for ratios of probabilistic integrals, such as are used to compute posterior probabilities. This approach offers performance superior to Monte Carlo methods by exploiting a Bayesian quadrature framework. We improve upon previous Bayesian quadrature techniques by explicitly modelling the non-negativity of our integrands, and the correlations that exist between them. It offers most where the integrand is multi-modal and expensive to evaluate. We demonstrate the efficacy of our method on data from the Kepler space telescope.

2. Active Learning of Model Evidence Using Bayesian Quadrature.

   Osborne, M.A., Duvenaud, D.K., Garnett, R., Rasmussen, C.E., Roberts, S.J., and Ghahramani, Z.

   In Advances in Neural Information Processing Systems (NIPS) 2012

   Abs

   Numerical integration is a key component of many problems in scientific computing, statistical modelling, and machine learning. Bayesian Quadrature is a model-based method for numerical integration which, relative to standard Monte Carlo methods, offers increased sample efficiency and a more robust estimate of the uncertainty in the estimated integral. We propose a novel Bayesian Quadrature approach for numerical integration when the integrand is non-negative, such as the case of computing the marginal likelihood, predictive distribution, or normalising constant of a probabilistic model. Our approach approximately marginalises the quadrature model’s hyperparameters in closed form, and introduces an ac- tive learning scheme to optimally select function evaluations, as opposed to using Monte Carlo samples. We demonstrate our method on both a number of synthetic benchmarks and a real scientific problem from astronomy.

2007

1. Further Explorations of Likelihood Theory for Monte Carlo Integration

   Kong, Augustine, McCullagh, Peter, Meng, Xiao-Li, and Nicolae, Dan L.

   Advances in Statistical Modelling and Inference 2007

   Abs link

   Monte-Carlo estimation of an integral is usually based on the method of moments or an estimating equation. Recently, Kong, McCullagh, Meng, Nicolae and Tan (2003) proposed a likelihood based theory, which puts Monte-Carlo estimation of integrals on a firmer, less ad hoc, basis by formulating the problem as a likelihood inference problems for the baseline measure with simulated observations as data. In this paper, we provide further exploration and development of this theory. After an overview of the likelihood formulation, we first demonstrate the power of the likelihood-based method by presenting a universally improved importance sampling estimator. We then prove that the formal, infinite-dimensional Fisher-information based variance calculation given in Kong et al. (2003) is asymptotically the same as the sampling based "sandwich" variance estimator. Next, we explore the gain in Monte Carlo efficiency when the baseline measure can be parameterized. Furthermore, we show how the Monte Carlo integration problem can also be dealt with by the method of empirical likelihood, and how the baseline measure parameter can be properly profiled out to form a profile likelihood for the integrals only. As a byproduct, we obtain four equivalent conditions for the existence of unique maximum likelihood estimate for mixture models with known components. We also discuss an apparent paradox for Bayesian inference with Monte Carlo integration.

2004

1. On a Likelihood Approach for Monte Carlo Integration

   Tan, Zhiqiang

   Journal of the American Statistical Association 2004

   Abs link

   The use of estimating equations has been a common approach for constructing Monte Carlo estimators. Recently, Kong et al. proposed a formulation of Monte Carlo integration as a statistical model, making explicit what information is ignored and what is retained about the baseline measure. From simulated data, the baseline measure is estimated by maximum likelihood, and then integrals of interest are estimated by substituting the estimated measure. For two different situations in which independent observations are simulated from multiple distributions, we show that this likelihood approach achieves the lowest asymptotic variance possible by using estimating equations. In the first situation, the normalizing constants of the design distributions are estimated, and Meng andWong’s bridge sampling estimating equation is considered. In the second situation, the values of the normalizing constants are known, thereby imposing linear constraints on the baseline measure. Estimating equations including Hesterberg’s stratified importance sampling estimator, Veach and Guibas’s multiple importance sampling estimator, and Owen and Zhou’s method of control variates are considered.

2003

1. A theory of statistical models for Monte Carlo integration

   Kong, Augustine, McCullagh, Peter, Meng, Xiao-Li, Nicolae, Dan L., and Tan, Zhiquiang

   Journal of the Royal Statistical Society, Series B (Statistical Methodology) 2003

   Abs link

   The task of estimating an integral by Monte Carlo methods is formulated as a statistical model using simulated observations as data. The difficulty in this exercise is that we ordinarily have at our disposal all of the information required to compute integrals exactly by calculus or numerical integration, but we choose to ignore some of the information for simplicity or computational feasibility. Our proposal is to use a semiparametric statistical model that makes explicit what information is ignored and what information is retained. The parameter space in this model is a set of measures on the sample space, which is ordinarily an infinite dimensional object. None-the-less, from simulated data the base-line measure can be estimated by maximum likelihood, and the required integrals computed by a simple formula previously derived by Vardi and by Lindsay in a closely related model for biased sampling. The same formula was also suggested by Geyer and by Meng and Wong using entirely different arguments. By contrast with Geyer’s retrospective likelihood, a correct estimate of simulation error is available directly from the Fisher information. The principal advantage of the semiparametric model is that variance reduction techniques are associated with submodels in which the maximum likelihood estima-tor in the submodel may have substantially smaller variance than the traditional estimator. The method is applicable to Markov chain and more general Monte Carlo sampling schemes with multiple samplers.

2002

1. Bayesian Monte Carlo

   Ghahramani, Zoubin, and Rasmussen, Carl E

   In Advances in neural information processing systems 2002

   Abs

   We investigate Bayesian alternatives to classical Monte Carlo methods for evaluating integrals. Bayesian Monte Carlo (BMC) allows the in- corporation of prior knowledge, such as smoothness of the integrand, into the estimation. In a simple problem we show that this outperforms any classical importance sampling method. We also attempt more chal- lenging multidimensional integrals involved in computing marginal like- lihoods of statistical models (a.k.a. partition functions and model evi- dences). We find that Bayesian Monte Carlo outperformed Annealed Importance Sampling, although for very high dimensional problems or problems with massive multimodality BMC may be less adequate. One advantage of the Bayesian approach to Monte Carlo is that samples can be drawn from any distribution. This allows for the possibility of active design of sample points so as to maximise information gain.

2000

1. Deriving quadrature rules from Gaussian processes

   Minka, T.P.

   2000

   Abs link

   Quadrature rules are often designed to achieve zero error on a small set of functions, e.g. polynomials of specified degree. A more robust method is to minimize average error over a large class or distribution of functions. If functions are distributed according to a Gaussian process, then designing an average-case quadrature rule reduces to solving a system of 2n equations, where n is the number of nodes in the rule (O’Hagan, 1991). It is shown how this very general technique can be used to design customized quadrature rules, in the style of Yarvin & Rokhlin (1998), without the need for singular value decomposition and in any number of dimensions. It is also shown how classical Gaussian quadrature rules, trigonometric lattice rules, and spline rules can be extended to the average-case and to multiple dimensions by deriving them from Gaussian processes. In addition to being more robust, multidimensional quadrature rules designed for the average-case are found to be much less ambiguous than those designed for a given polynomial degree.

1998

1. Bayesian quadrature with non-normal approximating functions

   Kennedy, Marc

   Statistics and Computing 1998

   link

1996

1. Iterative rescaling for Bayesian quadrature

   Kennedy, MC, and O’Hagan, A

   Bayesian Statistics 1996

1991

1. Bayes–Hermite quadrature

   O’Hagan, A.

   Journal of statistical planning and inference 1991

   Abs

   Bayesian quadrature treats the problem of numerical integration as one of statistical inference. A prior Gaussian process distribution is assumed for the integrand, observations arise from evaluating the integrand at selected points, and a posterior distribution is derived for the integrand and the integral. Methods are developed for quadrature in p. A particular application is integrating the posterior density arising from some other Bayesian analysis. Simulation results are presented, to show that the resulting Bayes–Hermite quadrature rules may perform better than the conventional Gauss–Hermite rules for this application. A key result is derived for product designs, which makes Bayesian quadrature practically useful for integrating in several dimensions. Although the method does not at present provide a solution to the more difficult problem of quadrature in high dimensions, it does seem to offer real improvements over existing methods in relatively low dimensions.

1987

1. Monte Carlo is Fundamentally Unsound

   O’Hagan, A.

   Journal of the Royal Statistical Society. Series D (The Statistician) 1987

   Abs link

   We present some fundamental objections to the Monte Carlo method of numerical integration.