## control functionals for Monte Carlo integration

**A** paper on control variates by Chris Oates, Mark Girolami (Warwick) and Nicolas Chopin (CREST) appeared in a recent issue of Series B. I had read and discussed the paper with them previously and the following is a set of comments I wrote at some stage, to be taken with enough gains of salt since Chris, Mark and Nicolas answered them either orally or in the paper. Note also that I already discussed an earlier version, with comments that are not necessarily coherent with the following ones! *[Thanks to the busy softshop this week, I resorted to publish some older drafts, so mileage can vary in the coming days.]*

First, it took me quite a while to get over the paper, mostly because I have never worked with reproducible kernel Hilbert spaces (RKHS) before. I looked at some proofs in the appendix and at the whole paper but could not spot anything amiss. It is obviously a major step to uncover a manageable method with a rate that is lower than √n. When I set my PhD student Anne Philippe on the approach via Riemann sums, we were quickly hindered by the dimension issue and could not find a way out. In the first versions of the nested sampling approach, John Skilling had also thought he could get higher convergence rates before realising the Monte Carlo error had not disappeared and hence was keeping the rate at the same √n speed.

The core proof in the paper leading to the 7/12 convergence rate relies on a mathematical result of Sun and Wu (2009) that a certain rate of regularisation of the function of interest leads to an average variance of order 1/6. I have no reason to mistrust the result (and anyway did not check the original paper), but I am still puzzled by the fact that it almost immediately leads to the control variate estimator having a smaller order variance (or at least variability). On average or in probability. (I am also uncertain on the possibility to interpret the boxplot figures as establishing super-√n speed.)

Another thing I cannot truly grasp is how the control functional estimator of (7) can be both a mere linear recombination of individual unbiased estimators of the target expectation and an improvement in the variance rate. I acknowledge that the coefficients of the matrices are functions of the sample simulated from the target density but still…

Another source of inner puzzlement is the choice of the kernel in the paper, which seems too simple to be able to cover all problems despite being used in every illustration there. I see the kernel as centred at zero, which means a central location must be know, decreasing to zero away from this centre, so possibly missing aspects of the integrand that are too far away, and isotonic in the reference norm, which also seems to preclude some settings where the integrand is not that compatible with the geometry.

I am equally nonplussed by the existence of a deterministic bound on the error, although it is not completely deterministic, depending on the values of the reproducible kernel at the points of the sample. Does it imply anything restrictive on the function to be integrated?

A side remark about the use of intractable in the paper is that, given the development of a whole new branch of computational statistics handling likelihoods that cannot be computed at all, intractable should possibly be reserved for such higher complexity models.

June 28, 2016 at 1:45 am

Thanks again for your interest and helpful feedback Christian – I hope we did manage to cover at least some of your points in the revision!

Perhaps I can add a footnote on our more recent developments? Decent progress has been made on better understanding the convergence rates of these estimators, and there is a pre-print of this new work available [1]. To spoil the plot of the new paper, we can show that convergence rates are governed by the (minimum of the) number of derivatives of the target density and the number of derivatives of the integrand. This means that problems which are smooth (in both of these senses) admit more rapid estimators – as you might perhaps expect from analogy with quasi Monte Carlo methods.

The more I think about these methods the closer they become to “kernel quadrature”; albeit the function space is rather non-standard. Nevertheless, these methods should be interpretable in terms of integral operators and eigen-decompositions, just like kernel quadrature (an example of that being Bayesian quadrature). That is the sort of direction I would like to move in with this work, going forward!

[1] Oates CJ, Cockayne J, Briol F-X, Girolami M. Convergence Rates for a Class of Estimators Based on Stein’s Identity. http://arxiv.org/abs/1603.03220