--- title: Zero-Variance Control Variates url: https://www.emergentmind.com/topics/zero-variance-control-variates-zvcv type: topic --- # Zero-Variance Control Variates Zero-Variance Control Variates (ZVCV) are a class of control-variate methods for Monte Carlo estimation that seek to reduce, and in ideal cases eliminate, estimator variance by adding functions with zero expectation under the target distribution. In the regression formulation of Monte Carlo integration, control variates appear as regressors with known mean, and the resulting estimator is the ordinary least squares estimator for the intercept in a multiple linear regression model; in the Stein-based formulation, zero-mean functions are generated by applying a Stein operator to a suitable test function. The “zero-variance” designation refers to the exact-fit case in which the integrand, or an appropriate transformed integrand such as a pathwise gradient, lies in the span of the chosen zero-mean features plus a constant, so that the residual is constant and the Monte Carlo estimator becomes exact [1801.01797]. ## 1. Formal definition and regression structure In the classical Monte Carlo setting, let \(X \sim \pi\), let \(f \in L^2(P)\), and let \[ I = E[f(X)] = P(f). \] The naive Monte Carlo estimator is \[ \hat{I}_n = P_n(f) = \frac{1}{n}\sum_{i=1}^n f(X_i), \qquad \operatorname{Var}(\hat{I}_n) = \frac{1}{n}\sigma^2(f), \qquad \sigma^2(f) = P\big[(f-P(f))^2\big]. \] If \(h=(h_1,\dots,h_m)^\top\) is a collection of control variates satisfying \(P(h_j)=0\), then for any \(\beta \in \mathbb{R}^m\), \[ I = P(f-\beta^\top h), \qquad \hat{I}_n(\beta)=P_n(f-\beta^\top h), \qquad \operatorname{Var}\{\hat{I}_n(\beta)\}=\frac{1}{n}\sigma^2(f-\beta^\top h). \] The population-optimal coefficients are \[ \beta_{\mathrm{opt}} = P(hh^\top)^{-1}P(hf), \] and the variance is minimized at \(\beta_{\mathrm{opt}}\) through the orthogonality relations \(P(\varepsilon)=0\) and \(P(h\varepsilon)=0\), where \(\varepsilon=f-I-\beta_{\mathrm{opt}}^\top h\) [1801.01797]. This representation is equivalent to a multiple linear regression model with intercept \(\mu=I\): \[ f(X_i)=\mu+\beta_{\mathrm{opt}}^\top h(X_i)+\varepsilon(X_i). \] Accordingly, Monte Carlo integration with multiple control variates can be implemented by ordinary least squares. In the uncentered empirical form, the empirical covariance matrix is \[ G_n = P_n(hh^\top)-P_n(h)P_n(h)^\top, \] and the OLS slope and intercept are \[ \hat{\beta}_n^{\mathrm{OLS}}=G_n^{-1}\{P_n(hf)-P_n(h)P_n(f)\}, \] \[ \hat{I}_n^{\mathrm{OLS}} = P_n(f)-\{P_n(fh^\top)-P_n(f)P_n(h^\top)\}G_n^{-1}P_n(h). \] In the centered formulation, with \(g_j=h_j-P(h_j)\), the normal equations become \[ \hat{\beta}=(G^\top G)^{-1}G^\top(f-\bar f\,\mathbf 1), \qquad \hat I_{\mathrm{CV}}=\bar f-\hat{\beta}^\top \bar g. \] The zero-variance condition is the exact-fit case. If \(f-I\) lies exactly in \(\operatorname{span}\{h_1,\dots,h_m\}\), meaning that there exists \(\beta\) with \(f=I+\beta^\top h\) almost surely, then \(\varepsilon \equiv 0\) and \[ \operatorname{Var}\big(\hat I_n^{\mathrm{OLS}}\big)=\frac{1}{n}P(\varepsilon^2)=0. \] In that case the OLS Monte Carlo rule integrates constants and the control variates exactly, and the estimator is exact for every sample set [1801.01797]. ## 2. Stein operators and the zero-mean construction A major route to ZVCV constructs zero-mean functions through Stein identities. For a differentiable target density \(\pi\), the Langevin–Stein operator acting on a sufficiently smooth vector field \(\phi:\Theta\to\mathbb R^d\) is \[ T_\pi \phi(x)=\nabla \log \pi(x)\cdot \phi(x)+\nabla \cdot \phi(x), \] and under suitable boundary and regularity conditions, \[ E_\pi[T_\pi \phi(X)]=0. \] A common scalar-function specialization sets \(\phi=\nabla g\), giving the second-order Langevin Stein operator \[ \mathcal L g(x)=\Delta g(x)+\nabla g(x)\cdot \nabla \log \pi(x). \] The resulting function \(h_g(x)=\mathcal L g(x)\) has zero mean under \(\pi\), so it can be added to the integrand without changing the target expectation [1811.05073]. This operator-theoretic construction makes the “zero-variance” principle explicit. If \[ f(x)=\mu+\mathcal A\varphi(x) \] for a Stein operator \(\mathcal A\) and a test function \(\varphi\) in the Stein class, then \(g(x):=\mathcal A\varphi(x)\) is a zero-mean control variate and the adjusted estimator equals the constant \(\mu\) for every sample set. In this sense, classical ZVCV is the exact representability of \(f-\mu\) in the image of the Stein operator [2109.08944]. In parametric ZVCV, \(g\) is typically restricted to a polynomial family. If \(P_j\) are monomials of total degree \(1,\dots,Q\), then the control-variate features are \[ \psi_j(x)=\mathcal L P_j(x)=\Delta P_j(x)+\nabla P_j(x)\cdot \nabla \log \pi(x). \] For \(Q=1\), the polynomial is linear and \(\Delta P_j=0\), so the features reduce to score components. The control-variate estimator then becomes \[ \hat\mu_{\mathrm{ZV}}=\frac{1}{n}\sum_{i=1}^n\big[f(x_i)+\beta^\top \psi(x_i)\big], \] with \(\beta\) fitted by least squares or weighted least squares when importance or sequential Monte Carlo weights are present [1811.05073]. Stein-based ZVCV also admits nonparametric generalizations. In control functional methods and their vector-valued extensions, one works in a reproducing kernel Hilbert space and applies the Stein operator to kernel sections. In the scalar case this recovers kernel-based zero-variance control variates; in the vector-valued case, a matrix-valued Stein reproducing kernel couples multiple integrals and yields control variates with componentwise zero mean under multiple targets. This suggests a continuum from finite-dimensional polynomial ZVCV to kernel-based control functionals, distinguished chiefly by the choice of function class in which the Stein image is approximated [2109.08944]. ## 3. Asymptotics with many control variates A distinctive theoretical development studies Monte Carlo integration when the number of control variates grows with the sample size. In the triangular-array formulation, \[ f_n=\mu_n+\beta_n^\top h_n+\varepsilon_n, \qquad P(\varepsilon_n)=0, \qquad P(\varepsilon_n h_{n,j})=0, \qquad \sigma_n^2=P(\varepsilon_n^2), \] with \(m_n \to \infty\). The relevant geometric quantity is the leverage function \[ q_n(x)=h_n(x)^\top P(h_n h_n^\top)^{-1}h_n(x), \qquad P(q_n)=m_n, \] together with the leverage condition \[ \sup_{x\in S} q_n(x)=o\!\left(\frac{n}{m_n}\right). \] This condition implies \(m_n=o(n^{1/2})\) and controls invertibility and stability of empirical Gram matrices [1801.01797]. Under the leverage condition, \[ \frac{\sqrt n}{\sigma_n}\big(\hat I_n^{\mathrm{OLS}}-\mu_n\big) = \frac{\sqrt n}{\sigma_n}P_n(\varepsilon_n)+o_P(1) = O_P(1), \] so in particular, \[ \hat I_n^{\mathrm{OLS}}-\mu_n=O_P(\sigma_n/\sqrt n). \] The central limit theorem is governed by a necessary and sufficient Lindeberg condition: \[ \text{for all }\delta>0:\quad P\Big[\big(\varepsilon_n/\sigma_n\big)^2\mathbf 1\{|\varepsilon_n/\sigma_n|>\delta \sqrt n\}\Big]=o(1). \] When it holds, \[ \frac{\sqrt n}{\sigma_n}\big(\hat I_n^{\mathrm{OLS}}-\mu_n\big)\xrightarrow{d}\mathcal N(0,1), \qquad \hat\sigma_{n,\mathrm{OLS}}^2/\sigma_n^2 \to 1 \quad\text{in probability.} \] Hence practical confidence intervals take the form \[ \hat I_n^{\mathrm{OLS}}\pm z_{1-\alpha/2}\,\hat\sigma_{n,\mathrm{OLS}}/\sqrt n. \] The nonstandard feature is the scaling by \(\sigma_n\), the regression residual standard deviation. If the linear span of \(\{1,h_1,\dots,h_{m_n}\}\) is dense in a function space containing \(f\), then \[ f=\mu+\beta_{m_n}^\top h_{m_n}+r_{m_n}, \qquad \|r_{m_n}\|_{L^2(P)}\to 0, \] so \(\sigma_n^2=P(r_{m_n}^2)\to 0\). The integration error then shrinks at rate \(\sigma_n/\sqrt n\), faster than the classical \(n^{-1/2}\) Monte Carlo rate. Explicit examples include post-stratification on \([0,1]\), univariate Legendre polynomials on \([-1,1]\), and multivariate tensor-product Legendre systems on \([-1,1]^d\), with approximation rates depending on the regularity of \(f\) and on dimensionality [1801.01797]. The exact ZVCV case sits at the edge of this asymptotic theory. If \(\sigma_n=0\) along a subsequence, the estimator is exact, but the CLT normalization becomes degenerate and the Lindeberg condition fails. ZVCV therefore represents exactness, not Gaussian fluctuation. ## 4. MCMC, Poisson equations, and Bayesian post-processing In Markov chain Monte Carlo, ZVCV is closely connected to the Poisson equation for the chain. For a Markov kernel \(P\) with stationary distribution \(\pi\), the Poisson equation \[ g-Pg=f-\pi(f) \] implies that \[ U(x)=g(x)-(Pg)(x) \] has zero mean under \(\pi\). If the exact solution \(g\) were available, then \(f+U\equiv \pi(f)\), yielding literal zero variance for ergodic averages. Practical MCMC ZVCV replaces the unknown Poisson solution by a finite span \(G=(G_1,\dots,G_K)^\top\) and defines control variates \[ U_a(x)=\sum_{k=1}^K a_k\big(\psi_k(x)-(P\psi_k)(x)\big). \] For reversible chains, the optimal coefficients admit an explicit finite-dimensional representation, and consistent estimators can be obtained from the same MCMC run [1008.1355]. This construction is particularly effective for conjugate random-scan Gibbs samplers, where \(P\psi\) is analytically tractable through conditional expectations. In the Gaussian example of the reversible-MCMC framework, the solution to the Poisson equation lies in the span of coordinate functions, and coordinate-based control variates can achieve dramatic reductions in asymptotic variance. The same methodology extends to certain Metropolis–Hastings samplers and hybrid Metropolis-within-Gibbs algorithms when \(P\psi\) is tractable or can be approximated [1008.1355]. A second major Bayesian post-processing line uses Stein operators rather than the Markov operator. In this setting, one fits \(f\) by a finite-dimensional space \[ F=\operatorname{span}\{1\}\oplus L\Phi, \] where \(L\) is the Stein operator and \(\Phi\) is typically a polynomial basis, then augments this with a reproducing-kernel correction. The resulting semi-exact control functionals (SECF) are exact on a specified finite-dimensional class and retain a kernel-based nonparametric component. In the Gaussian Bernstein–von–Mises regime, with \(\Phi=P^r\), the SECF estimator is exact on the polynomial class \(P_0^r=\operatorname{span}\{1\}\oplus P^r\). The same work establishes a bias-correction property when the Markov chain is not invariant for the posterior, proving \(O_P(n^{-1/2})\) convergence under the stated assumptions [2002.00033]. For stochastic gradient MCMC, ZVCV can be applied as a post-processing step using the gradient of the log posterior. With a linear polynomial \(Q(\theta)=a^\top\theta\), the Stein feature reduces to a linear combination of score components, and the adjusted output becomes \[ \hat g(\theta_k)=g(\theta_k)+a^\top z_k, \qquad z_k=\frac12\,\widehat{\nabla f}(\theta_k). \] The method is described as available “for free” because the gradients are already computed during stochastic-gradient sampling. When noisy minibatch gradients replace exact scores, the estimator remains unbiased under unbiased gradient estimation, but the attainable variance reduction depends on the gradient-noise level; the paper analyzes this effect and shows that variance-reduced gradient estimators preserve stronger ZVCV gains at large dataset size [1706.05439]. ## 5. High-dimensional, vector-valued, and algorithmic extensions A central practical difficulty is dimensionality. For polynomial ZVCV of total degree \(Q\) in dimension \(d\), the number of monomials is \[ \binom{d+Q}{Q}, \] so the number of nonconstant regression covariates is \[ p=\binom{d+Q}{Q}-1. \] This combinatorial growth rapidly makes ordinary least squares unstable or infeasible. Regularized ZV-CV addresses the problem with penalized regression, \[ \arg\min_\beta \left[ \sum_{i=1}^n \big(f(x_i)+\beta^\top \psi(x_i)\big)^2 + \lambda\{\alpha\|\beta\|_1+(1-\alpha)\tfrac12\|\beta\|_2^2\} \right], \] interpolating between LASSO, ridge, and elastic net. A second device, called a priori regularization, restricts the polynomial to a subset of coordinates \(S\subset\{1,\dots,d\}\), reducing the number of features to \[ p_S=\binom{|S|+Q}{Q}-1. \] These methods were introduced precisely to stabilize high-dimensional fits and to reduce computational and storage costs [1811.05073]. A related recent development replaces explicit penalization by ensembles of ordinary least squares fits on random subsets of polynomial Stein features. In ensemble ZVCV, one constructs \(k\) component regressions using only \(J^*