{'fraction_non_alphanumeric': 0.03866968930588497, 'fraction_numerical': 0.019917969415433806, 'mean_word_length': 4.34374343533366, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Regularized Zero-Variance Control Variates

Zero-variance control variates (ZV-CV) are a post-processing method to reduce the variance of Monte Carlo estimators of expectations using the derivatives of the log target. Once the derivatives are available, the only additional computational effort lies in solving a linear regression problem. Significant variance reductions have been achieved with this method in low dimensional examples, but the number of covariates in the regression rapidly increases with the dimension of the target. In this paper, we present compelling empirical evidence that the use of penalized regression techniques in the selection of high-dimensional control variates provides performance gains over the classical least squares method. Another type of regularization based on using subsets of derivatives, or a priori regularization as we refer to it in this paper, is also proposed to reduce computational and storage requirements. Several examples showing the utility and limitations of regularized ZV-CV for Bayesian inference are given. The methods proposed in this paper are accessible through the R package ZVCV.


Introduction
Our focus in this paper is on calculating the expectation of a square integrable function ϕ(θ) with respect to a distribution with (Lebesgue) density p(θ), θ ∈ Θ ⊆ R d . Given independent and identically iid ∼ p(θ), the standard Monte Carlo estimator, is an unbiased estimator of E p [ϕ(θ)] = Θ ϕ(θ)p(θ)dθ and its variance is O(1/N ). Reducing the variance of this estimator by increasing N is often infeasible due to the cost of sampling from p(θ) and potentially the cost of evaluating ϕ(θ). If the samples are not iid then the functional form of the estimator is the same and the methods described in this work still apply.
(ZV-CV) (Assaraf and Caffarel, 1999;Mira et al., 2013) and control functionals (CF)  are two such methods. ZV-CV amounts to solving a linear regression problem and CF is a nonparametric alternative. These methods can be used as post-processing procedures after N samples, not necessarily iid, from p have been produced along with evaluations of ∇ θ log p(θ) and ϕ(θ) for each of the samples. Often ∇ θ log p(θ) is already available because derivative-based methods like Metropolis adjusted Langevin algorithm (MALA) (Roberts and Stramer, 2002;Girolami and Calderhead, 2011) or Hamiltonian Monte Carlo (HMC) (Duane et al., 1987;Girolami and Calderhead, 2011) have been used in the sampling algorithm.
The parametric approximation in ZV-CV is based on a polynomial in θ, so the number of coefficients to estimate rapidly increases both with the polynomial order and with d. As a result of restricting to a low polynomial order, ZV-CV tends to offer less substantial improvements than CF for challenging lowdimensional applications. This is not surprising given the good statistical properties of CF which have been described in Oates et al. (2019); Barp et al. (2021). However, CF has an O(N 3 ) computational cost, compared to ZV-CV with has computational cost of O(N ), and it also suffers from the curse of dimensionality with respect to d due to the use of non-parametric methods. Some results in Oates et al. (2017), shown mainly in the appendices, suggest that the performance of CF compared to ZV-CV may deteriorate in higher dimensions.
One aim of this work is to develop derivative-based control variate methods which are inexpensive, effective and capable of handling higher dimensions than existing derivative-based methods. The novel methods that we introduce are referred to as regularized ZV-CV and they are based on two types of regularization: penalization methods for linear regression and what we refer to as a priori regularization. Penalized ZV-CV allows higher order polynomials to be used than could be considered with ordinary least squares. This method is motivated by showing that L 2 penalized ZV-CV is equivalent to CF with a second-order differential operator and finite-dimensional polynomial kernel. A priori ZV-CV is most beneficial when N < d. Empirical results in Section 4 suggest that significant variance reductions can be achieved with a priori ZV-CV when N < d or with penalized ZV-CV when the polynomial order is pushed beyond the limits of what standard ZV-CV can handle. We have developed an R package, ZVCV (South, 2018), which implements several derivative-based variance reduction techniques including standard ZV-CV, CF and regularized ZV-CV.
An important application area for ZV-CV and regularized ZV-CV is Bayesian inference, where Monte Carlo integration is commonly used. The use of ZV-CV and CF to improve posterior expectations based on samples from Markov chain Monte Carlo (MCMC) (Metropolis et al., 1953) is well established (see e.g. Mira et al. (2013); Papamarkou et al. (2014); Oates et al. (2017)). ZV-CV and CF have also been applied to the power posterior (Friel and Pettitt, 2008) estimator of the normalising constant in an MCMC setting by Oates et al. (2016) and Oates et al. (2017), where they refer to this method as controlled thermodynamic integration (CTI). In this paper we go beyond existing literature and describe how regularized ZV-CV fits naturally into the context of sequential Monte Carlo (SMC) samplers (Del Moral et al., 2006;Chopin, 2002). In doing so, we provide a setting where adaptive methods can easily be applied to the CTI estimator. A novel reduced-variance normalising constant estimator using the standard SMC identity is also proposed.
An introduction to derivative-based Monte Carlo variance reduction methods is provided in Section 2. The main methodological contributions in terms of developing regularized ZV-CV methods can be found in Section 3. Section 4 contains a simulation study comparing methods and estimators on the novel application to SMC. A final discussion of limitations and possible future work is given in Section 5.

Control Variates based on Stein Operators
In this section, we recall previous work on control variate methods. The classical framework for control variates (Ripley, 1987;Hammersley and Handscomb, 1964) is to determine an auxiliary functioñ ϕ(θ) = ϕ(θ) + h(θ) such that E p [φ(θ)] = E p [ϕ(θ)] and V p [φ(θ)] < V p [ϕ(θ)], where V p denotes the variance with respect to p(θ). Estimator (1) can then be replaced with the unbiased, reduced variance estimator, A control variate which has been considered in (Assaraf and Caffarel, 1999;Mira et al., 2013;Barp et al., 2021) is where L is a second-order Langevin Stein operator (Stein, 1972;Gorham and Mackey, 2015) depending on p, ∆ θ is the Laplacian operator represented in coordinates as d j=1 ∇ 2 θ[j] on R d , · is the dot product operator such that a · b = a b and g : Θ → R is a twice continuously differentiable function to be specified.
Langevin Stein operators are helpful in generating control variates for two reasons. Firstly, they can be applied in Bayesian inference because they do not require the normalising constant of p. Furthermore, by definition a Stein operator L depending on p(θ) satisfies E p [Lg(θ)] = 0 for all functions g(θ) in a set called a Stein function class (see Section 2.2) and therefore E p [φ(θ)] = E p [ϕ(θ)] under mild conditions. Typically it is also a requirement that Stein operators L fully characterize p but this is not necessary for Stein-based control variates.
What remains is to choose g. The function g for whichφ(θ) is constant, and thus zero variance is achieved, is generally intractable. In practice, g is restricted to some function class G and is estimated based on samples targeting p.

Choice of Function g
Variance reduction is effected through judicious choice of g. Once a function class G has been selected, the function g ∈ G is estimated by performing a regression task. As described in Barp et al. (2021), a generalization of several existing regression methods for this problem is (ĉ,ĝ) ∈ arg min c∈R g∈G where PEN(g) is a penalty function to be specified and λ ≥ 0. This amounts to a penalized least squares approach to estimating ϕ(θ) using the functional form c − Lg(θ). This perspective on the optimization problem encompases ZV-CV, CF and neural control variates (Zhu et al., 2019) as special cases. Further details on how ZV-CV and CF fit into this framework are given below. The main developments in this paper are based on considering alternative penalty functions.
Two recent contributions in control variates have optimization functions which do not fit into this framework, though the developments in penalty functions that are proposed in this paper could still be considered in these alternative frameworks. Belomestny et al. (2017) consider empirical variance minimization since minimising a square error objective function may not be optimal. Brosse et al. (2019) consider an alternative optimization problem which is motivated by minimising the asymptotic variance of a Langevin diffusion, which may be more suitable when samples have been obtained using MCMC with multivariate normal random walk or MALA proposals.

Control Functionals
CF Barp et al., 2021) is based on choosing G ≡ H where H is a user-specified Hilbert space of twice differentiable functionals on Θ. The penalty term PEN(g) considered in CF is PEN(g) = g 2 H , where · H is the norm associated with the Hilbert space H. The existence of a solution pair (ĉ,ĝ) ∈ R × H, together with an explicit algorithm for its computation, was obtained in that work under the assumption that the Hilbert space H admits a reproducing kernel (see Berlinet and Thomas-Agnan (2011) for background). This method leads to estimators with super-root-N convergence under conditions described in Oates et al. (2019) and Barp et al. (2021). However, the cost associated with computation ofĝ is O(N 3 ), due to the need to invert a dense kernel matrix, and moreover this matrix is typically not well-conditioned. For applications that involve MCMC and SMC, typically N will be at least 10 3 and thus (in the absence of further approximations) the algorithm of Oates et al. (2017); Barp et al. (2021) can become impractical.

Zero-Variance Control Variates
ZV-CV (Assaraf and Caffarel, 1999;Mira et al., 2013) amounts to using G as the class of Qth order polynomial functions in θ, and λ = 0. The polynomials P (θ) that we consider have total degree Q ∈ Z ≥0 , meaning that the maximum sum of exponents is Q and the monomial basis is θ (3), where P j (θ) is the jth monomial in the polynomial and β ∈ R J is the vector of polynomial coefficients, gives The jth zerovariate (covariate in the regression), x j (θ) = LP j (θ) is a term containing θ and ∇ θ log p(θ). Its exact form is given in Appendix A of the Online Resources. For a Qth order polynomial when the dimension of θ is d, the constant J = d+Q d − 1 is the number of regression parameters, excluding the intercept which is in the null space of the Stein operator.
The standard approach in the literature for choosing β is to perform ordinary least squares (OLS) (Glasserman, 2003). This is equivalent to choosing λ = 0 in (4). The computational cost of ZV-CV is O(J 3 + N J 2 ), which scales better with N than CF which has computational cost O(N 3 ), where often J N . Unlike in CF, regularization methods have not previously been used in connection with ZV-CV.
Common practice is to default to Q = 2 in ZV-CV. Mira et al. (2013) consider Q = 1 to at most Q = 3 and find that Q = 2 is sufficient to achieve orders of magnitude variance reduction in most of their examples. Papamarkou et al. (2014) consider Q ≤ 2, pointing out that "first and second degree polynomials suffice to attain considerable variance reduction." Low polynomial orders are also typically used in applications, for example Baker et al. (2019) use Q = 1 and Oates et al. (2016) use Q ≤ 2. Oates et al. (2017) compare CF with ZV-CV using Q = 2 in most examples.
It has previously been proposed to increase the number of control variates as the sample size increases (see e.g. Portier and Segers (2019) and the appendices of Oates et al. (2017)). This approach can be motivated by the Stone-Weierstrass theorem (Stone, 1948), which states that polynomial functions can be used to uniformly approximate, to an arbitrary level of precision, continuous functions on closed intervals. However, the increased number of coefficients in higher order polynomials may not be feasible or efficient to estimate with OLS.
We demonstrate in Section 4 that the common practice of defaulting to Q = 2 with OLS is often suboptimal. Regularization approaches are proposed in Section 3 to enable higher degree polynomials to be employed whilst avoiding over-fitting of the regression model used to estimate the optimal coefficients of the polynomials.

Unbiasedness
Suppose that g(θ) and p(θ) are sufficiently regular so that log p(θ) has continuous first order derivatives and g(θ) has continuous first and second order derivatives. Also suppose that, if g is to be estimated, then the samples used in estimating g(θ) are independent of those used in evaluating (2). If Θ = R d , then we require that Θ is compact and has piecewise smooth boundary ∂Θ. Under these conditions, where ∂Θ is a surface integral over ∂Θ, n(θ) is the unit vector orthogonal to θ at the boundary ∂Θ and S(dθ) is the surface element at θ ∈ ∂Θ. When Θ = R d is unbounded, condition (5) becomes a tail condition which is satisfied if Γr p(θ)∇ θ g(θ) · n(θ)S(dθ) → 0 as r → ∞ where Γ r ∈ R d is a sphere centred at the origin with radius r and n is the unit vector orthogonal to θ at Γ r . This requirement is given in Equation 9 of Mira et al. (2013) and Assumption 2 of Oates et al. (2017) and it is a direct result of applying the divergence theorem to E p [Lg(θ)] = 0. In the ZV-CV context for Θ = R d , a sufficient condition for (5) is that the tails of p decay faster than polynomially (Appendix B of Oates et al., 2016).
The unbiased estimator which uses independent samples for estimation of g(θ) and evaluation of (2) is referred to as the "split" estimator. In practice, the so-called "combined" estimator which uses the full set of N samples for both estimation of g(θ) and evaluation of (2) can have lower mean square error than the split estimator but is no longer unbiased. If MCMC methods are employed then bias is unavoidable and the combined estimator is likely to be preferred.

Parameterization
An additional consideration when performing either ZV-CV or CF is the adopted parameterization. Any deterministic, invertible transformation of the random variables ψ = f (θ) can be used so one can estimate instead of (2), where p θ ≡ p is the probability density function for θ, p ψ is the probability density function for ψ obtained through a change of measure and {ψ i } N i=1 ∼ p ψ . For simplicity, the θ parameterization is used in notation throughout the paper. The best parametrization to adopt for any given application is an open problem. If the original parameterization does not satisfy boundary condition (5), one could consider a reparameterization such that the boundary condition is satisfied.

Regularized Zero-Variance Control Variates
The aim of this section is to develop methods which are computationally less demanding than CF and offer improved statistical efficiency over standard ZV-CV. We describe two types of regularization: regularization through penalized regression and a priori regularization. The latter is primarily for cases where not all derivatives of the log target are available or when N d. Combinations of the two regularization ideas are also possible. Methods to choose between control variates are described in Section 3.3.

Regularization Through Penalized Regression
As mentioned earlier, the number of regression parameters in ZV-CV grows rapidly with the order Q of the polynomial and with the dimension d of θ. Therefore, the polynomial order that could be considered is limited by the number of samples required to ensure existence of a unique solution to the OLS problem, eliminating the potential reduction that could be achieved using higher order polynomials. In this section, we propose to use penalized regression techniques to help overcome this problem.
In most contexts, using penalized regression reduces variance at the cost of introducing bias. Recall that the conditions for unbiasedness in Section 2.2 do not depend on the mechanism for estimating β, as long as the samples used in estimating β are independent of those used in evaluating (2). Thus, the use of penalized regression methods does not introduce bias into ZV-CV.
The regularization methods introduced in Sections 3.1.1 and 3.1.2 involve a penalty function on β so we use standardization for stability and to be able to employ a single λ. The regression problem becomes: where the subscript s is in reference to the response and predictors being standardized by their sample mean and standard deviation. Specifically, using the notationā = 1 )/σ xj for j = 1, . . . , J and β s represents the coefficients on this standardized scale. The estimated coefficients on the original scale areβ [j] The parameter λ is chosen to minimize the k-fold cross-validation mean square error. The first type of penalization that we consider is Tikhonov regularization (Tikhonov et al., 2013), or ridge regression as it is known when applied in regression (Hoerl and Kennard, 1970). This involves using a squared L 2 penalty, PEN(g) = β s 2 2 . Ridge regression mitigates overfitting and allows for estimation when the regression problem is ill-posed due to a small number of observations. Closed form solutions forĉ andβ are available, leading to the same computational cost as OLS of O(J 3 + N J 2 ). The use of L 2 penalized ZV-CV can also be motivated using the results of Belkin et al. (2019), who argued that using J ≥ N with an interpolation-based approach (i.e. CF or regularized ZV-CV) can lead to better mean square loss compared to restricting to J ≤ N (i.e. standard ZV-CV) in situations where there is no reason to pre-suppose the first J basis functions are also the most useful. The latter condition may be satisfied when ϕ is too complex to be well-approximated using control variates based on low order polynomials.
To motivate this particular form of penalization, we now consider interpreting this method as a computationally efficient variant of CF. To facilitate a comparison with the approach of Barp et al. (2021), we consider a particular instance of CF with a reproducing kernel Hilbert space H that is carefully selected to lead to an algorithm with lower computational cost. Namely, we select a polynomial kernel where P j (θ) denotes the jth of all J monomial terms in θ up to order Q. For such a kernel, a welldefined Hilbert space H = span{P j } j=1,...,J is reproduced and we have an explicit expression for the Hilbert norm , which reveals the method of Barp et al. (2021) as an L 2 -penalized regression method. As such, the optimization problem in ZV-CV with PEN(g) = β s 2 2 is equivalent to the optimization problem in CF (without the standardization of the response and predictors) and it can be solved as a least-squares problem with complexity O(J 3 + N J 2 ). The first main contribution of our work is to propose a more practical alternative to the method of Barp et al. (2021), which we recall has O(N 3 ) computational cost, by using such a finite-dimensional polynomial kernel. Our results in this direction are empirical (only) and we explore the properties of this method for various values of Q in Section 4. Tikhonov regularization has been applied implicitly in the context of CF but, to the best of our knowledge, this is the first time that general penalized regression methods have been proposed in the context of ZV-CV. Results in Section 4 demonstrate that the new estimators can offer substantial variance reduction in practice when the number of samples is small relative to the number of coefficients being estimated. The principal aim in the design of a control variate h is to accurately predict the value that the function ϕ takes at an input θ * not included in the training dataset {(θ i , ϕ(θ i ))} N i=1 . It is well-understood that L 1 -regularization can outperform L 2 -regularization in the predictive context when the function ϕ can be well-approximated by a relatively sparse linear combination of predictors. In our case, the unstandardized predictors are the functions in the set {1} ∪ {LP j } j=1,...,J . Given that low-order polynomial approximation can often work well for integrands ϕ of interest, it seems plausible that L 1 -regularization could offer an improvement over the L 2 -regularization used in Oates et al. (2017); Barp et al. (2021). Investigating this question is the second main contribution of our work.
In the context of ZV-CV, L 1 -penalization can be interpreted as using the least absolute shrinkage and selection operator (LASSO, Tibshirani (1996)). LASSO introduces an L 1 penalty PEN(g) = β s 1 where ||β s || 1 = j |β s [j]|. The effect of the penalty is that some coefficients are estimated to be exactly zero.

A priori Regularization
As an alternative to penalized regression methods, in this section we consider restricting the function g to vary only in a lower-dimensional subspace of the domain Θ ⊆ R d . More specifically, a subset of parameters S ⊆ {1, . . . , d} is selected prior to estimation and the function g is defined, in a slight abuse of notation, as g(θ) = P (θ[S]). The log target derivatives, ∇ θ log p(θ), only appear in the control variates (3) through the dot product ∇ θ g(θ) · ∇ θ log p(θ). Therefore if j / ∈ S then the derivative ∇ θ[j] log p(θ) is not required. We refer to this approach as a priori regularization.
A priori regularization makes ZV-CV feasible when some derivatives cannot be used, for example due to intractability, numerical instability, computational expense or storage constraints. An example of where some derivatives may be difficult to obtain is in Bayesian inference for ordinary differential equation (ODE) models. Evaluating ∇ θ log p(θ) requires the sensitivities of the ODE to be computed, which involves augmenting the system of ODEs with additional equations. If some additional equations render the system stiff, then more costly implicit numerical solvers need to be used and in such cases it would be useful to avoid including sensitivites corresponding to the difficult elements of θ. It may also be infeasible to use the O(d) storage required to run standard ZV-CV. Storing a subset of the parameters and derivatives for use in a priori regularization may, however, be achievable. Another benefit of a priori ZV-CV is that it reduces the number of coefficients to estimate, making estimation feasible when N d. Zhuo et al. (2018) consider similar ideas to a priori ZV-CV in the context of Stein variational gradient descent, where they use the conditional independence in p(θ) for probabilistic graphical models to separate high dimensional inference problems into a series of lower dimensional problems.
The downside of using a priori ZV-CV is that the potential for variance reduction is reduced, except for under both conditions (a) θ[S] is independent of θ[S] according to p(θ), whereS = {1, . . . , d} \ S, and (b) ϕ(θ) = ϕ(θ[S]). Outside of this situation, restricting the polynomial to g(θ) = P (θ[S]) will give varying levels of performance depending on the subset that is selected. Intuitively, one may wish to choose the subset of variables so that θ[S] and/or ∇ θ[S] log p(θ) have high correlations with ϕ(θ). In practice, this is easiest to do when there is a priori knowledge and therefore not all derivatives need to be calculated and stored. Given (b), it is suspected that this method will be more useful for individual parameter expectations than for expectations of functions of multiple parameters.
Estimators using this approach are unbiased under the same conditions as ZV-CV and penalized ZV-CV. This method is also applicable to CF, though nonlinear approximation may be more difficult in this non-parametric setting.

Automatic Selection of Control Variates
The performance of regularized ZV-CV depends upon the polynomial order, the penalization type and on S. We demonstrate in Section 4 that the common practice of defaulting to Q = 2 with OLS is often sub-optimal and also that the optimal control variate depends on a variety of factors including N and p(θ). It has previously been proposed to increase the number of control variates as the sample size increases (see e.g. Portier and Segers (2019) and the appendices of Oates et al. (2017)). However, in these existing works the mechanism whereby the complexity of the control variate was increased was not data-dependent.
To choose between control variates in this work, we use 2-fold cross-validation so that our selection is data-dependent. For each combination of penalization type and S, we start with polynomial order Q = 1 and we continue to increase the polynomial order until the average cross-validation error is larger for Q + 1 than for Q. The combination of regularization method and polynomial order which gives the minimum cross-validation error is selected and we perform estimation using that method on the full set of samples. The cross-validation error that we use here is the sums of square residuals in the hold-out set, averaged across the two folds.

Empirical Assessment
In this section, we perform comparisons of regularized ZV-CV to ZV-CV and CF on Bayesian inference examples. In Bayesian statistics, the posterior distribution of the parameters θ of a statistical model given observed data y is where the function (y|θ) is the likelihood function, p 0 (θ) encorporates prior information and Z is a normalising constant. Interest is in estimating posterior expectations Θ ϕ(θ)p(θ|y)dθ and the normalising constant or so-called "evidence" Z = Θ (y|θ)p 0 (θ)dθ for Bayesian model choice. Posterior expectations and Z are typically analytically intractable and challenging to estimate due to the potentially high dimensional integration required.  2017) have also applied ZV-CV and CF, respectively, to a thermodynamic integration (Gelman and Meng, 1998;Ogata, 1989;Friel and Pettitt, 2008) estimator for the evidence, calling the resulting method controlled thermodynamic integration (CTI). The thermodynamic integration estimator gives the log evidence as the sum of multiple expectations with respect to p t where p t = (y|θ) t p 0 (θ)/Z t and t is referred to as the inverse temperature. Oates et al. (2017) use population Monte Carlo (Jasra et al., 2007) to obtain the samples from p t for t = 0, . . . , T and they consider specifically t j = (j/T ) 5 . A total of 2(T + 1) expectations are involved, with ZV-CV applied to the estimator for each expectation.

ZV
We propose to use sequential Monte Carlo (SMC, Del Moral et al. (2006)) with the tuning method of Salomone et al. (2018) for sampling, rather than the standard choices of MCMC or population MCMC. The benefit of this approach is that the samples are roughly independent which can be preferable over the high autocorrelation that can be seen in MCMC samples. The standard SMC evidence estimator is the product of T expectations, so we consider improving this estimator using ZV-CV and CF. Further details about implementation in SMC and the advantages of this approach are given in Appendix B of the Online Resources. From the perspective of comparing variance reduction methods, the application of ZV-CV to posterior expectations and to multiple evidence estimators means that ZV-CV can be compared on a variety of functions ϕ(θ) and distributions p(θ).
We perform an empirical comparison of the following methods using examples of varying complexity: • vanilla: Monte Carlo integration without control variates.
• ZV Q : ZV-CV with OLS and order Q polynomial.
• l-ZV Q : ZV-CV with LASSO and order Q polynomial.
• r-ZV Q : ZV-CV with ridge regression and order Q polynomial.
• sub k -: This prefix indicates a priori ZV-CV with a subset of size k. Applications are limited to d > 1 dimensions. We only apply sub kideas to posterior expectations since ϕ(θ) is a function of a single parameter and a potentially reasonable subset may be known a priori.
• crossval: Control variate selection using 2-fold cross-validation. This method chooses between ZV Q , l-ZV Q , r-ZV Q and sub kwhere applicable.
Methods written with the prefix sub k -or the name l-ZV Q , r-ZV Q or crossval are novel for all Q and k. The main purposes of these comparisons are to investigate the performance of higher order polynomials, the utility of penalized regression and the ability to achieve variance reduction using a subset of derivatives. The purpose of the comparisons to CF is not necessarily to outperform CF, as CF can be infeasible to apply in its basic form for large N , but to benchmark the performance of these novel methods against CF. A variety of sample sizes, integrands ϕ(θ) and target distributions p are used for fair comparisons. We focus on sample sizes that are typical of SMC, ranging from N = 10 to N = 10, 000 but we note that larger sample sizes can be accommodated by the regularized ZV-CV methods which have a computational complexity of O(N ).
Estimators are compared on the basis of mean square error (MSE), where the gold standard of estimation is carefully chosen for each example. The main quantity of interest reported in this section is , the MSE of the vanilla Monte Carlo estimator estimated from 100 independent SMC runs divided by the estimated MSE for the method in question. This quantity is referred to as statistical efficiency and it is reported for each fixed N . Values above one are preferred.
Control variate methods are most valuable when the sampling algorithm is expensive, for example due to the cost of evaluating the likelihood, or when evaluation of the function ϕ(θ) is costly. The overall efficiency, as measured by , is also considered for these examples. Here time[·] is the average time across the 100 runs to compute the estimator in question, including the time spent running the SMC sampler. We note that the run time is subject to the efficiency of the code and here (penalized) ZV-CV is based on the R package glmnet, cross-validated ZV-CV is written as a loop in R and CF is implemented in C++. Nevertheless, our proposed methods offer improved overall efficiency in several of the applications considered. The computational benefits of our approach will improve with increasing model complexity in terms of likelihood calculations, since the overhead associated with penalized regression will become relatively negligible.
Two examples are described in detail in this section. Appendices E, F and G also include results for a 61-dimensional logistic regression example, a one-dimensional ODE example which motivates higher order polynomials and a challenging nine-dimensional ODE model, respectively.
In terms of bias, boundary condition (5) is satisfied using the specified parameterizations for all examples considered in this paper. This can be verified through the sufficient condition that the tails of p decay faster than polynomially and Θ = R d (Appendix B of Oates et al. (2016)). However, the estimators are generically biased due to the use of SMC, as they would be with MCMC. All results are based on combined estimators as opposed to split estimators, so all pairs {θ i , ϕ(θ i )} N i=1 are used to buildφ and also to estimate E p [φ(θ)].

Recapture Example
This 11-dimensional example demonstrates that reduced variance estimators can be obtained with the use of higher order polynomials and regularization. Marzolin (1988) collected data on the capture and recapture of the bird species Cinclus cinclus over six years. Like Brooks et al. (2000), Nott et al. (2018) and South et al. (2019b), we use a Bayesian approach to estimate the parameters of a Cormack-Jolly-Seber model (Lebreton et al., 1992) for the capture and recapture of this species. The parameters of the Cormack-Jolly-Seber model used here are the probability of survival from year i to i + 1, φ i , and the probability of being captured in year k, p k , where i = 1, . . . , 6 and k = 2, . . . , 7. Denote the number of birds released in year i as D i and the number of animals caught in year k out of the number released in year i as y ik . It is simple to show that the number released in year i that are never caught is d i = D i −  The likelihood is given by where θ = (φ 1 , . . . , φ 5 , p 2 , . . . , p 6 , φ 6 p 7 ). Following South et al. (2019b), the parameters φ 6 and p 7 are multiplied together due to a parameter identifiability issue.
The prior is θ[j] ∼ U(0, 1) for j = 1, . . . , 11. To satisfy the boundary conditon (5) and to improve the efficiency of MCMC proposals, all parameters are transformed to the real line using ψ The gold standard of evidence estimation for this example is the mean evidence estimate for l-ZV 1 at N = 5000. The posterior expectation gold standard is the average posterior mean for ZV 4 at N = 5000.

Posterior Expectations
The average statistical efficiency and overall efficiency across parameters is shown in Figure 1, excluding a priori regularization results for simplicity. Higher order polynomials become more efficient as N increases and the use of penalized regression means that higher order polynomials can be considered for smaller N . LASSO regression is preferable over ridge regression for this example.
Using a priori ZV-CV with S = j, where j is the index of the current parameter of interest, sub 1 -ZV 1 is on average roughly 10 times more efficient than vanilla Monte Carlo integration.
Cross-validation generally gives similar results to CF and to the best performing fixed method. More details about the selected control variates can be found in Appendix C of the Online Resources.

Evidence Estimation
Regularized ZV-CV and automatic control variates give improved statistical efficiency over ZV-CV and CF for the range of N that are considered here, as seen in Figure 2. However, there is less improvement in terms of overall efficiency due to the fact that multiple expectations are required for evidence estimation. This puts the more computationally intensive methods including higher order polynomials, cross-validation and CF at a significant disadvantage. We note that this example was selected to allow for extensive comparisons and the cost of post-processing would have less impact under more expensive likelihood functions.
The selected control variates for N = 50 and N = 1000 can be found in Appendix C of the Online Resources.

Log-Gaussian Cox Point Process Example
We now consider an example where the dimension can be adjusted. The log-Gaussian Cox point process example of Møller et al. (1998) consists of locations of 126 Scots pine saplings in a 10 × 10 m 2 plot. The plot can be discretised into n × n grid cells, so that the dimension d = n 2 of the problem can be varied. Here we consider n = 4, n = 8 and n = 16 so that we have Bayesian inference problems of size d = 16, d = 64 and d = 256.
The model specifications, including code for the log likelihood, log prior and their gradients, match that of Heng and Jacob (2019). After normalising the plot to fit onto a unit square, the number of points at grid cell (i, j) for i, j = 1, . . . , n is denoted y i,j . It is assumed that the y i,j are conditionally independent and Poisson distributed with mean where θ i,j has a Gaussian process prior with mean µ and covariance function Σ Our goal is to estimate the posterior means for the parameters θ i,j for i, j = 1, . . . , n, and we do so using SMC runs with N = 100 particles. We do not consider evidence estimation due to lack of a reliable gold standard. Due to memory and time constraints, the maximum polynomial order is constrained so that the maximum number of covariates in ZV-CV is 5000. The gold standards in this example are the average posterior expectations across many independent unbiased Riemann-manifold HMC runs (Heng and Jacob, 2019) with unbiased control variates as described in South et al. (2019a). Details of the gold standard are available in Appendix D.
Tables 1, 2 and 3 show the mean relative statistical, computational and overall efficiency for posterior expectations in dimensions d = 16, d = 64 and d = 256, respectively, when N = 100. In all three settings, the existing methods (vanilla MC, ZV-CV with OLS and CF) are outperformed by the novel approaches of a priori ZV-CV, cross-validation and LASSO with a higher order polynomial than OLS could handle. The best performing novel method has an overall efficiency which is better than the best performing existing method by a factor of over 1,600,000 for d = 16, over 1,800 for d = 64 and over 25 for d = 256. Results showing the competitive performance of the novel methods for N = 1, 000 and N = 10, 000 are given in Appendix D. Like the results for N = 100, a priori ZV-CV, cross-validation and LASSO outperform existing alternatives in the majority of settings.

Discussion
In this paper, we introduced two types of regularized ZV-CV: regularization through penalized regression and regularization by selecting a subset of parameters to include in the regression model. Higher order polynomial basis functions have the potential to outperform the commonly used polynomial with Q = 2 as N -the number of Monte Carlo, MCMC or SMC simulations -increases. Our penalized ZV-CV ensures that the resulting functional approximation problem remains well-defined when N is less than the number of control variate coefficients (J) while performing similarly to standard ZV-CV when N > J. For the examples considered here, we found that LASSO generally resulted in better performance than ridge regression. A priori ZV-CV led to significant improvements over vanilla Monte Carlo for posterior expectations, with little computational overhead.
One of the main applications of the proposed methods is in models where the dimension, d, is too high for standard variance reduction techniques to be efficient. Empirical evidence suggests that using ZV-CV and penalized ZV-CV, where Q is increased with N , offers better statistical performance than CF in high dimensions. However, the computational cost of (penalized) ZV-CV is O(N d+Q which may prohibit the application of these methods with large Q in high dimensions. This explosion in complexity for large Q and d is a disadvantage relative to CF when the sample size is comparable or less than the dimension, though the complexity is similar when N ≈ d and Q = 1 in (penalized) ZV-CV. Table 1: 16-dimensional Cox example: efficiency for marginal posterior expectations when N = 100, averaged over results for all 16 parameters. A "−"' indicates that the population size N = 100 is insufficient for standard ZV-CV. We refer the reader to the beginning of Section 4 for acronym definitions.
Approach Q Stat. Efficiency Comp. Efficiency Overall Efficiency Vanilla NA 1.0 × 10 0 1.0000 3.5 × 10 6 0.9996 3.5 × 10 6 sub 1 -ZV 4 4 6.9 × 10 7 0.9995 6.9 × 10 7 crossval NA 2.8 × 10 8 0.2125 5.5 × 10 7  One could consider speeding up these algorithms by using partial LASSO searches (e.g. Efron et al., 2004;Fan and Lv, 2008) or by using approximate solvers as proposed in Si et al. (2021). Alternatively, in very large dimensions, the a priori ZV-CV approach can be used to obtain variance reductions with a complexity that is O(N |S| 2 + |S| 3 ) where 1 ≤ |S| ≤ d. This a priori approach also offers benefits when not all derivatives are available, when N d, or when information about the relationships between the integrand and parameters is known (for example when p(θ) has a directed acyclic graph factorization). Leluc et al. (2019) provide additional theoretical support for LASSO-based control variate selection. The work of Leluc et al. (2019), which was publicly available after the pre-print of our paper , gives concentration inequalities for the integration error with LASSO-based control variates and also shows that the correct control variates are selected with high probability. The theoretical results are based on bounded control variates, which do not apply in ZV-CV and CF when p has unbounded support. Leluc et al. (2019) find empirically that a methodological adjustment of performing OLS for estimation once the control variates have been selected via LASSO is helpful in reducing the variance of the estimator. We point out that this modification is necessary to obtain the zero-variance property of ZV-CV. The optimal coefficients required to obtain zero-variance estimators cannot be obtained directly from penalized regression methods like LASSO and ridge regression with non-zero λ.
We have proposed the consideration of different penalty functions in the optimization problem for control variates, but we focus specifically on LASSO and ridge regression. Some other potentially useful regularization methods for the situation where N < d+Q d are elastic net (Zou and Hastie, 2005) and partial least squares (PLS, Wold (1975)). Elastic net is a compromise between LASSO and ridge regression which uses two tuning parameters. PLS is based on choosing the k < d+Q d − 1 independent linear combinations of covariates that explain the maximum variance in the response, where k is chosen through cross-validation. Active subspaces (Constantine, 2015) are a more recent dimension-reduction technique which use the derivatives of the function of interest to find the linear combinations of covariates that are best at predicting the function. It would be of interest in future research to compare our LASSO and ridge regression ZV-CV methods with these alternatives.
The concept of regularization by selecting a subset of parameters is referred to as nonlinear approximation in approximation theory and applied mathematics (DeVore, 1998), and there is some theoretical evidence to suggest that this can outperform linear approximation (e.g. penalized regression which is described in Section 3.1). Selecting a particular subset of monomials which are used in a polynomial interpolant is also the same idea as in sparse grid algorithms for numerical integration (Smolyak, 1963). These methods are known to work well in high dimensions and could be useful alternatives for selecting the subset of monomials in ZV-CV.
Stein-based control variates using neural networks have recently appeared in the literature (Zhu et al., 2019). Zhu et al. (2019) added details of penalization methods to their approach, where the control variates cannot be fitted exactly and stochastic optimization is required. Penalization methods are simpler and more stable in the linear regression context but in future research it would be of interest to compare to neural control variates with regularization. This alternative approach is likely to outperform ZV-CV in some applications, such as when ϕ(θ) is multi-modal.
Finding the optimal parameterization for a given application is a challenging open problem. Choosing the parameterization is a trade-off between making p φ simpler and making ϕ(f −1 (φ)) simpler. Another potential benefit of reparameterising for ZV-CV is that there is the potential to enforce more sparsity in the predictors for improved performance in L 1 penalization.
Derivatives are available in closed form or can be unbiasedly estimated for a large class of problems. ZV-CV has been applied in big data settings in the context of post-processing after stochastic gradient MCMC (Baker et al., 2019) and for models with intractable likelihoods (Friel et al., 2016). Regularized ZV-CV also applies in these settings. Regularized ZV-CV could also be used in exact approximate settings where a particle filtering estimate of the likelihood is used (see for example Dahlin et al. (2015) and Nemeth et al. (2016)). However, derivative-based methods are most appealing when the derivative of the log target can be obtained with little additional cost relative to the likelihood itself. An interesting avenue for future research may be to consider automatic differentiation.

A Covariates in the ZV-CV Regression
As described in Section 2.1 of the main paper, ZV-CV uses control variates of the form β x(θ) where β ∈ R J , J = d+Q d − 1, θ ⊆ R d and Q is the polynomial order. This appendix gives the general form of x(θ) for a fixed Q.
The jth element of x(θ) is: where j = 1, . . . , J. The matrix A ∈ Z J×d ≥0 has J rows where each row corresponds to a unique vector A j,· such that 1 ≤ d k=1 A j,k ≤ Q. In other words, A contains all permutations of powers to θ[1], . . . , θ[d] that lead to a sum of exponents between 1 and Q.
It is straightforward to verify that x = ∇ θ log p(θ) ∈ R d for the first order polynomial P

B Variance Reduction in Sequential Monte Carlo
In this appendix, the novel applications of Stein-based variance reduction to SMC are described.

B.1 Sequential Monte Carlo
SMC samplers are naturally adaptive and parallelizable alternatives to standard MCMC for sampling from the posterior of static Bayesian models (Del Moral et al., 2006). A set of N weighted samples, {W i j , θ i j } N i=1 , are moved through a sequence of distributions, p tj (θ|y), for j = 0, . . . , T . The distributions p tj , henceforth p j for brevity, are properly normalized and η j represents the unnormalized distributions, i.e. p j ∝ η j for j = 0, . . . , T . The samples, or particles, are moved through these distributions using importance sampling, resampling and move steps.
The importance sampling step reweights particles {W i j , θ i j } N i=1 from p j−1 to target p j using , for i = 1, . . . , N , where w i j is the unnormalized weight for particle θ i j = θ i j−1 and W i 0 = 1/N if independent and identically distributed samples are drawn from the initial distribution p 0 . The normalized weights are W i Resampling, most commonly multinomial resampling, is used to remove particles with negligible weights and replicate particles with high weights. After resampling, the weights are set to 1/N . Finally a move step, most commonly in the form of several iterations of a p j -invariant MCMC kernel, is used to diversify the particles. Derivative based proposals have recently been used for the MCMC kernel in SMC (Sim et al., 2012), which means that ZV-CV can easily be performed on expectations with respect to p j for j = 1, . . . , T .
It is straightforward to adapt SMC algorithms online. Recent work has proposed adapting the MCMC kernel parameters online using the population of particles. For example, Fearnhead and Taylor (2013) and Salomone et al. (2018) propose adaptation methods for generic MCMC kernels. Buchholz et al. (2018) propose methods for performing the notoriously challenging tuning of HMC kernel parameters in SMC. It is also possible to adaptively choose whether to perform the resampling and move steps based on some measure of the weight degeneracy (Del Moral et al., 2012). In this context, the weights are updated when resampling and move steps are not performed but the particle values remain the same.
Commonly used sequences for p j in the literature are data annealing and likelihood annealing. In data annealing SMC (Chopin, 2002), the data are introduced sequentially so the targets are p j (θ|y 1:j ) where y 1:j denotes the first j data points. Likelihood annealing smoothly introduces the effect of the likelihood to help explore complex targets (Neal, 2001) through the sequence p j (θ|y) = (y|θ) tj p 0 (θ)/Z j , the same sequence that is used in thermodynamic integration (TI, Gelman and Meng (1998); Ogata (1989)). Unlike in the MCMC setting, where t j for j = 0, . . . , T needs to be fixed for TI, the inverse temperatures in SMC can easily be adapted online, for example using the methods described in Jasra et al. (2011) and detailed in Appendix B.5. We use likelihood annealing SMC in all examples for these reasons.

B.2 Handling Weights in ZV-CV
A weighted least squares can be performed in ZV-CV to take into account these weights. Recall that the optimization problem for unweighted samples

B.3 Posterior Expectations
The most straightforward way to estimate posterior expectations in SMC is to perform Monte Carlo integration using the final set of particles obtained after all T iterations. When resampling is performed at the final iteration, as we have done here, the final set of samples is the unweighted set {θ i T } N i=1 . The vanilla Monte Carlo estimator for the jth marginal posterior mean, , where θ[j] denotes the jth marginal. We apply ZV-CV and CF to improve this estimator.
It is also possible to estimate the posterior mean by using weighted particles from previous likelihood annealing targets. This recycling would extend the work of Briol et al. (2017) for improving tail coverage with the split estimator and could also improve the performance of higher order polynomials due to higher degrees of freedom. The weighted least squares approach in (8) could be applied in this setting.

B.4 Evidence Estimation
An added benefit of SMC samplers over alternatives like standard MCMC is that an estimate of the normalising constant is produced as a by-product in SMC. Two evidence estimators that can be obtained as a by-product of SMC are described below. Both of these evidence estimators are based on a likelihood annealing schedule, but applications of Stein-based variance reduction techniques to the data annealing SMC evidence estimator for streaming data would also be possible.

B.4.1 Controlled Thermodynamic Integration
Oates et al. (2016) and Oates et al. (2017) have applied ZV-CV and CF, respectively, to the power posterior estimator for the evidence (Friel and Pettitt, 2008). The estimator, is based on TI and it gives the log evidence as an integral with respect to the inverse temperature t, where p t = (y|θ) t p(θ)/Z t . For points in a discrete set {t j } T j=0 of inverse temperatures where 0 = t 0 < . . . < t T = 1, this integral is estimated using quadrature methods.
The second order quadrature method of Friel et al. (2014) estimates (9) using A simpler, first order quadrature approximation which is equivalent to the first sum in (10) was described in Friel and Pettitt (2008).
It is straightforward to apply ZV-CV to (10) by noticing that the estimator is simply a sum of expectations where Oates et al. (2016) refer to the use of ZV-CV for this purpose as controlled TI (CTI) and they implement CTI in an MCMC framework. In the MCMC context, the inverse temperatures are fixed prior to the runs, making it difficult to balance low quadrature bias from using a large number of inverse temperatures with low computational effort from using a small number of inverse temperatures. Furthermore, sampling from the target distributions p tj for j = 1, . . . , T requires tuning, which is often done manually in MCMC. One contribution of our work is to use the CTI estimator in the SMC framework which allows for online choice of the inverse temperature schedule and online tuning of the MCMC proposal. Another contribution, which is described in Appendix B.5, is a method for adjusting inverse temperatures after the MCMC or SMC runs have completed. We note that (9) has been used in the SMC framework by Zhou et al. (2012), but CTI has not previously been used in SMC.

B.4.2 SMC Evidence Estimator
The standard SMC evidence estimator (Del Moral et al., 2006) is based on the telescoping product Assuming that independent and identically distributed samples are drawn from the prior p 0 , so Z 0 = 1, the normalising constant can be written as the product of expectations as follows: The standard estimator, Z = T j=1 N i=1 W i j−1 (y|θ i j−1 ) tj −tj−1 , is unbiased when adaptive methods are not used. We propose the use of ZV-CV on each of the expectations in (11) to obtain a lower variance estimator, i.e. we take p(θ) in Section 2 of the main paper to be the power posterior p j (θ|y). Although ZV-CV and CF can lead to unbiased estimators for each individual expectation, it is not clear whether the product of these estimators,Ẑ, remains unbiased. Nevertheless, we find in practice that the SMC estimator with ZV-CV has lower mean square error than the SMC estimator without control variates.

B.4.3 Comparisons
The CTI estimator appears to be more amenable to a control variates treatment, so we propose its use over the SMC estimator when practitioners wish to make the best use of derivatives for evidence estimation. The two evidence estimators both involve expectations of non-linear transformations to (y|θ), but the integrand is on the logarithmic scale for the CTI estimator which may make it a simpler function to estimate. However, the results for the challenging ordinary differential equation example in Appendix G indicate that the SMC estimator may be preferable when the prior is highly diffuse. In this context, extreme values in the draws from the prior can lead to high variance in the CTI estimator.
When adaptive resampling is performed, the particles do not have equal weights at every iteration. ZV-CV is able to take into account these weights by using weighted means and weighted linear regression to estimate the ZV-CV coefficients β. In the combined CF estimator, the information contained in the weights is lost due to its interpolating property. This may explain its poor performance for evidence estimation in Section 4.1 of the main paper.

B.5 Post-hoc Temperature Choice
Using an insufficient number of inverse temperatures can lead to significant bias in the power posterior log evidence estimator. However, it is difficult to know a priori how many inverse temperatures will be required to achieve reasonably small bias. This appendix briefly describes some existing methods for choosing the inverse temperature schedule before describing our post-hoc approach. The proposed method is useful for both CTI and SMC evidence estimation.
The simplest approach for choosing the inverse temperatures is to use a fixed schedule, for example t j = (j/T ) 5 for j = 0, . . . , T (Friel and Pettitt, 2008). If this schedule is conservative in that T is very high, then some costly resampling and move steps can be avoided by performing these steps only when an approximation to the effective sample size (ESS) becomes low. The ESS is the number of independent samples from the target that would be required to achieve the same variance of the estimator and the ESS at inverse temperature t j is approximated in SMC by 1/ Friel et al. (2014) describe a method for adaptively choosing the inverse temperature schedule with the goal of minimising the discretization error in the power posterior log evidence estimator. Their method is used to calibrate the inverse temperatures prior to implementing the full sampler and it may be useful when evidence estimation is the primary focus.
In the SMC context, the most popular method for adaptively choosing the inverse temperatures is based on fixing the approximated ESS at ρN using the bisection method, where 0 < ρ < N (Jasra et al., 2011). This approach maintains a fixed discrepancy between p j−1 and p j for j = 1, . . . , T when resampling and move steps are performed at each iteration. Zhou et al. (2015) use the conditional ESS (CESS), instead of the ESS when resampling is not performed at every inverse temperature, because this is a more accurate measure of the discrepancy between p j−1 and p j when p j−1 is approximated with a weighted sample.
Using the approaches above, it is difficult to be confident that the quadrature bias will be sufficiently low without being overly conservative. Performing ZV-CV on the power posterior log evidence estimator requires T regressions to be performed, that is one at each inverse temperature regardless of whether resample and move steps were performed. A conservative choice of inverse temperature schedule increases the post-processing time in ZV-CV.
We propose a post-hoc method to adjust the inverse temperatures when the original choice is either not conservative enough or too conservative. To start with, any approach can be used to provide T distinct sets of particles and inverse temperatures, Inverse temperatures at which the particles are not moved (for example due to adaptive resampling methods in SMC) are not included in these T inverse temperatures. If (12) is too conservative or not conservative enough, then a new set of inverse temperatures {t j }T j=0 is selected as follows.
Given an inverse temperaturet j−1 (starting att 0 = 0), the bisection method is used to selectt j such that CESSt j =ρN . This process continues until an inverse temperature oftT = 1 satisfies CESStT ≥ρN . Each of the inverse temperatures must be assigned a relevant particle population and this is done by selecting population Put simply, the method involves choosing inverse temperatures post-hoc so that the CESS is fixed at ρN where 0 <ρ < N . The inverse temperatures {t j } T j=0 need not appear in {t j }T j=0 , but they may for some choices of initial inverse temperatures andρ.

B.6 Implementation Details
The adaptive SMC methods used to select the tuning parameters and inverse temperatures are described in Salomone et al. (2018) and Jasra et al. (2011), respectively. The post-hoc method for adapting the inverse temperatures is described in Appendix B.5. Table 4 gives all tuning parameter specifications. The number of particles in the adaptive SMC run is N . Inverse temperatures in the adaptive SMC run are chosen to maintain an ESS of ρN and they are adjusted post-hoc to maintain a CESS ofρN . The MCMC moves targeting p j use MALA proposals of the form whereΣ j is the empirical covariance and h j is a tuning parameter. We specify a set of 20 values which are log-uniform on the range of h min to h max and, following Salomone et al. (2018), we select the value which maximizes the highest median estimated expected square jumping distance (ESJD, Pasarica and Gelman (2010)). Finally, we choose the number of MCMC repeats so that a given percentage of particles have a total absolute jumping distance greater than the mean Mahalanobis distance between particles before resampling (Salomone et al., 2018). To avoid confounding the effects of ZV-CV with the effects of kernel parameter and inverse temperature choice, we do a single adaptive SMC run for each example. The inverse temperatures and kernel parameters from this run are then used in 100 independent SMC runs for each value of N .  The full set of inverse temperatures are used for both the CTI and SMC estimators. Additional inverse temperatures without resample and move steps do not improve the vanilla SMC evidence estimator, but we find that they can lead to substantial reductions in the ZV-CV SMC evidence estimator.
To improve stability, the ZV-CV regression for the SMC evidence estimator is performed using the integrand divided by its maximum value and the results are adjusted to correct for this.

C Additional Results for the Recapture Example
The selected control variates for the 1100 posterior expectations at N = 50 and N = 1000 can be seen in Table 5. Control variate selection for evidence estimatation at N = 50 and N = 1000 is shown in Figure 3. It is clear that higher order polynomials are selected for larger N , and lasso is selected more often than ridge regression. These general trends are consistent with the best performing fixed methods.

D Additional Results for the Log-Gaussian Cox Process Example
A gold standard estimate is required for all expectations of interest in order to assess the statistical efficiency of the different methods. For the Log-Gaussian Cox Process example, the gold standard of approximation comes from unbiased Riemannian-manifold Hamiltonian Monte Carlo (RM-HMC) with coupling as per (Heng and Jacob, 2019). Control variates based on minimisation of an upper bound on the variance are also used, as described in South et al. (2019a). The resulting estimators are unbiased, meaning that many cores can be used in parallel to estimate the expectations.
For all three models, the estimates were averaged over 30 repeats where the following settings were used. The RM-HMC part of the algorithm had a step size of 0.11 with 10 steps. The probability of selecting the coupled random walk Metropolis Hastings (RWMH) algorithm rather than the RM-HMC sampler was 0.05 and the proposal standard deviation of the RWMH was 0.001. All algorithms had a burn-in of 70 iterations. The total number of MCMC iterations for d = 16, d = 64 and d = 256 were 100, 000, 15, 000 and 4000, respectively. Each of the 30 repeats involved running R = 32 independent sets of coupled Markov chains. Control variates were applied such that half of the R = 32 lots of samples were used to estimate control variate coefficients for the other half of the samples. This was done to maintain unbiasedness. The specific classes of control variates used were third order ZV-CV, second order ZV-CV and first order ZV-CV for d = 16, d = 64 and d = 256, respectively.
Tables 6-11 show the mean relative statistical, computational and overall efficiency for posterior expectations in dimensions d = 16, d = 64 and d = 256 when N = 1, 000 and N = 10, 000. For N = 1, 000, the overall efficiency of the best novel method compared to the best existing method (vanilla MC, ZV-CV with OLS or CF) is a factor of over 600 for d = 16, over 600 for d = 64 and over 1.35 for d = 256. For N = 10, 000, the overall efficiency of the best novel method compared to the best existing method is a factor of over 3 for d = 16, over 4 for d = 64 and over 0.98 (a slight decrease) for d = 256. Table 6: 16-dimensional Cox example: efficiency for marginal posterior expectations when N = 1, 000, averaged over results for all 16 parameters. A "−"' indicates that the population size N = 1, 000 is insufficient for standard ZV-CV. We refer the reader to the beginning of Section 4 for acronym definitions.

E Sonar Example
Using all derivative information to perform ZV-CV with higher order polynomials may simply be unrealistic for some examples, due to storage constraints or to the number of regression parameters required. The 61-dimensional example below falls into this class of problems. A polynomial with Q = 2 for this example requires a restrictive 1953 regression parameters while Q = 3 and Q = 4 polynomials require a potentially unrealistic number of samples with 41,664 and 677,040 regression parameters, respectively.
This binary regression problem is based on discriminating between sonar signals bouncing off a metal cylinder versus sonar signals bouncing off a roughly cylindrical rock. The 60 covariates represent the total energy within a given frequency band over a fixed period of time, with increasing aspect angle from covariate 1 to covariate 60. Observations y i for i = 1, . . . , 208 are coded as 0 for rock and 1 for metal. The corresponding log likelihood is where X ∈ R 208×61 is the matrix of covariates starting with a column of 1's for the intercept, y ∈ R 208 is the vector of indicators for the response and θ ∈ R 61 is the vector of coefficients. Using (14) with {0, 1} encoding of the response is equivalent to using log (y, X|θ) = − 208 i=1 log(1 + exp(−y i X i,· θ)) with {−1, 1} encoding (see e.g. Hastie et al. (2015)). Following Chopin and Ridgway (2017), we standardize the predictors (columns 2-61 of X) to have standard deviation 0.5 and we use N (0, 5 2 ) priors. The intercept has a N (0, 20 2 ) prior. This prior specification is chosen over the Cauchy priors of Gelman et al. (2008) so that the expectations exist and boundary condition (2.4) from the main paper is satisfied. The data used here is from the UCI machine learning repository (Dheeru and Karra Taniskidou, 2017) and was originally collected by Gorman and Sejnowski (1988). This example is more challenging than most standard logistic regression problems due to the the high number of covariates and high correlations between covariates (Chopin and Ridgway, 2017). The aspect angle is increasing from covariate 1 to covariate 60, so it is reasonable to assume that coefficients of nearby covariates will be correlated. This gives useful information in choosing the subset of parameters for a priori ZV-CV. We perform ZV-CV with subsets of 1 or 5 parameters in estimating marginal posterior expectations. The polynomials with a subset of five parameters are based on the four closest angles, for example when estimating the posterior expectation of θ[30], the polynomial is a function of a subset of the five parameters S = {28, 29, 30, 31, 32}. The gold standard of posterior expectation in this example is the average posterior expectation across 100 independent SMC runs using N = 10000 and ZV-CV with Q = 2. Due to memory and time constraints, a maximum polynomial order of Q = 2 was considered in the cross-validation. Table 12 shows the mean statistical efficiency for posterior expectations. For small N , one can obtain better results using a subset of parameters than with l-ZV 1 , whereas for large N polynomials with Q = 2 become more efficient. Out of the 6100 expectations for N = 50, 83.9% of the control variates selected based on cross-validation use the subset of five parameters, 3.1% use the full 61 parameters with penalized regression methods and 13.1% use the subset of one parameter. At N = 5000, 97.5% of selected control variates are l-ZV 2 , 2.5% are r-ZV 2 and less than 1% use Q = 1. As suspected, due to the high dimension in this example, ZV-CV methods outperform CF. Table 13 shows the overall efficiency for each posterior expectation. Like the recapture example, the likelihood function in this example is relatively inexpensive. The proposed methods offer up to two orders of magnitude in overall efficiency compared to vanilla MC, and a priori ZV-CV is the best performing method for N = 50. N ZV 1 l-ZV 1 r-ZV 1 ZV 2 l-ZV 2 r-ZV 2 sub 1 -ZV 1 sub 5 -ZV 1 crossval CF 50 1.3 1.  N ZV 1 l-ZV 1 r-ZV 1 ZV 2 l-ZV 2 r-ZV 2 sub 1 -ZV 1 sub 5 -ZV 1 crossval CF 50

F Van der Pol Example
This example based on the Van der Pol oscillatory differential equations (Van der Pol, 1926) demonstrates the potential to improve the performance of ZV-CV by using higher order polynomials. Regularization is not required in this one-dimensional example.
In the Van der Pol oscillatory differential equations, position x(t) is governed by the second order differential equation d 2 x dt 2 − θ(1 − x 2 ) dx dt + x = 0, which can also be written as a system of first order differential equations, where x 1 = x and x 2 = dx/dt. The parameter θ represents the non-linearity of the system and the damping strength. We use the same data as Oates et al. (2017) which is based on θ = 1, initial position x 1 = 0, initial velocity x 2 = 2 and noisy observations y t ∼ N (x 1 (t), 0.1 2 ) recorded at times t = 0, 1, . . . , 10. The prior is log(θ) ∼ N (0, 0.25 2 ). Derivatives are obtained by augmenting the system of differential equations with sensitivity equations (see e.g. Appendix 7 of Oates et al. (2017)). Oates et al. (2017) found that second order polynomials resulted in mediocre performance for evidence estimation in this example. Oates et al. (2017) used population MCMC with within-temperature and between-temperature proposals to obtain samples from the power posteriors. The sampling algorithm used here is SMC, which is simpler and easier to adapt. Given that the prior is such that θ > 0, we have performed a log transform so that the MCMC proposals and ZV-CV are both based on ψ = log θ.
Expectations are calculated using (2.5) of the main paper, for example the posterior expectation of θ can be expressed as E p ψ [exp(ψ)].
The gold standard evidence estimate 1 for this example is based on numerical integration. For the posterior mean, the gold standard of θ is based on the mean of 100 estimates using ZV 9 at N = 2000.

F.0.1 Posterior Expectations
With this simple integrand, ZV-CV is able to perform extremely well and outperforms CF. It is apparent from Figure 5 that larger Q is preferable for this problem, with the selected polynomial order from cross-validation increasing with N . Here the likelihood function is expensive so the majority of the time is taken up by sampling and the difference between the two measures of efficiency is negligible. Figure 6 shows the statistical and overall efficiency of the two evidence estimators for different N . It is clear from these figures that higher order polynomials can offer substantial improvements in statistical efficiency over low order polynomials given a sufficient number of samples. CF performs similarly to the best performing ZV-CV method for the CTI estimator and typically slightly better than the best performing ZV-CV method for the SMC evidence estimator, at the cost of relatively poor performance for some values of N due to outlying estimates. The overall trends are similar for overall efficiency, although the additional expectations lead to an increased cost in postprocessing which is especially noticeable for CF.

F.0.2 Evidence Estimation
The cross-validation method results in estimators with similar MSE to fixed polynomial orders of Q = 2 or Q = 3. Although the chosen control variates lead to efficiency improvements of one to three orders of   magnitude over the vanilla estimator, the improvements are not as substantial as for the optimal fixed Q ∈ {1, 2, . . . , 9}. The cross-validation method is stopping at a sub-optimal Q, potentially because there is little to no improvement from one (sub-optimal) polynomial order to the next (sub-optimal) polynomial order for many of the expectations involved here. This hypothesis is supported by the polynomial order selection results in Table 14. One potential solution in practice may be to force the cross-validation to compare up to at least some fixed polynomial order.
At N = 10, the evidence estimates for the CTI estimator appear to be approximately unbiased ( Figure  7a). The estimators also have low variance with the exception of the polynomial with Q = 9, for which there are not enough samples for a reasonable fit. The SMC combined estimator suffers from more bias (Figure 7b), but this bias disappears when the samples are split to facilitate independent estimation of g(θ) and evaluation of (2.1) of the main text (results not shown). q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q

G Ordinary Differential Equation Example
Potential limitations of ZV-CV and regularized ZV-CV are illustrated through this challenging example with non-linear posterior dependencies. For the mostpart, the performance with low sample sizes is poor and the performance with larger sample sizes is modest compared to the examples in Section 4.1 of the main paper and Appendices D and E of these Online Resources. The likelihood function is relatively expensive in this example, but we consider statistical efficiency only because this is sufficient to illustrate the potential limitations of ZV-CV and regularized ZV-CV.
Geyer (1991) describes the following system of coupled ordinary differential equations (ODEs) for modelling biochemical pathways, This example has been considered from a Bayesian context in Girolami (2008), South et al. (2019b) and Salomone et al. (2018). Following South et al. (2019b), observations y(t) of R pp (t) are observed with noise such that y(t) ∼ N (R pp (t), 0.02 2 ) for t = 0, 3, . . . , 57. The nine parameters of interest are θ = (k 1 , V 1 , K m1 , K m2 , V 2 , S 0 , D 0 , R 0 , R pp0 ), the last four being the initial values of S, D, R and R pp . We use the same simulated data as South et al. (2019b) which is based on θ = (0.05, 0.20, 0.1, 0.1, 0.1, 1.00, 0, 1.00, 0) and we use the priors from Girolami (2008), where G(a, b) represents the Gamma distribution with shape parameter a, scale parameter b and mean ab. The log-transformed parameters, ψ = log θ, are used for MCMC proposals and for ZV-CV through (2.5) of the main paper.
As explained in South et al. (2019b), this example has non-linear posterior dependencies (see Appendix G of South et al. (2019b) for an illustrative figure) and is challenging partially because of numerous practical and structural identifiability issues in the model. It is also interesting to note that no information is obtained about D 0 by observing R pp , so D 0 is practically non-identifiable. D 0 is therefore independent of other parameters and of the data, so the posterior marginal for D 0 is simply its prior, G(1, 0.1). This helps to explain the extremely good performance of ZV-CV and subset regularized ZV-CV in estimating its posterior mean. The parameter R pp0 is well identified based on the data and is therefore conditionally independent of other parameters given the data, making it another parameter which is potentially easier to estimate.
The gold standard of posterior expectation and evidence estimation is the mean from 100 fixed SMC runs with N = 10000, except for the posterior mean for D 0 for which we use the true value of 0.1. When using the SMC evidence estimator, it is possible to obtain a negative estimate for the evidence and this happened in some runs with low N for this example. To avoid this issue, we repeat the coefficient estimation with a fixed intercept ofĉ = 1 N N i=1 ϕ(θ i ) when the evidence estimate is negative. For CF, we instead replace the estimator with 1 N N i=1 ϕ(θ i ) in this case, which amounts to using vanilla MC.

G.0.1 Posterior Expectations
Tables 15-23 show the statistical efficiency for each posterior expectation. Given the complex target distribution and strong dependencies between parameters, it is not possible to achieve improvements on vanilla Monte Carlo integration using N < 100 for most parameters. CF outperforms ZV-CV for some marginals, but the improvements are mostly similar to ZV-CV. The worst performing marginal expectation improves on vanilla Monte Carlo integration by at factor of at most 5.4 at N = 1000.
There is generally no difference in efficiency between vanilla Monte Carlo integration and a priori ZV-CV, with the exceptions of D 0 and R pp0 . Even for a sample size as small as N = 10, using ZV-CV either with a full polynomial in ψ or with a polynomial only in log D 0 gives a posterior mean for D 0 that is correct to 15 significant figures. The improvement is 30 orders of magnitude smaller for the same expectation with CF. Using a polynomial in log R pp0 gives an estimate of R pp0 which is 1.5 times more efficient than vanilla Monte Carlo integration at N = 1000.
N ZV 1 l-ZV 1 r-ZV 1 ZV 2 l-ZV 2 r-ZV 2 ZV 3 l-ZV 3 r-ZV 3 ZV 4 l-ZV 4 r-ZV 4 sub 1 -ZV  N ZV 1 l-ZV 1 r-ZV 1 ZV 2 l-ZV 2 r-ZV 2 ZV 3 l-ZV 3 r-ZV 3 ZV 4 l-ZV 4 r-ZV 4 sub 1 -ZV  The vanilla Monte Carlo estimators of the evidence have very high variance in this example. From Figure 8, it can be seen that the distribution of SMC evidence estimates using vanilla Monte Carlo is positively skewed. The estimator underestimates with a high probability and overestimates by a large amount with small probability. The behaviour of the CTI estimator is similar. Figure 8b illustrates the performance of ZV-CV and CF for N = 1000 compared to vanilla Monte Carlo integration with the gold standard of N = 10000. ZV-CV and CF significantly reduce the variance, but they introduce a negative bias. As suspected, the MSE is reduced by using ZV-CV and CF but the amount by which it is reduced is not stable across N (Tables 24 and 25) owing to the introduction of bias. Further research may be required to investigate the performance of ZV-CV and CF when the vanilla estimator is highly positively skewed.
Despite the challenges in this example, automatic control variate selection performs reasonably well.   Table 24: ODE example: statistical efficiency of the CTI evidence estimates for a range of N . Ridge regression performs either similarly or worse than LASSO and has therefore been excluded from this table. The * indicates that cross-validation efficiency values for N = 100, N = 500 and N = 1000 are based on 98, 99 and 97 runs, respectively. The remaining 6 runs reached a polynomial order which required too much RAM.