Zero-Variance Control Variates

• Zero-Variance Control Variates is a variance reduction technique that adds zero-mean functions to the integrand to improve Monte Carlo estimation accuracy.
• It employs regression and Stein-based formulations to achieve an exact fit when the integrand lies in the span of chosen zero-mean features, resulting in zero estimator variance in ideal cases.
• The method extends to high-dimensional, Bayesian, and MCMC settings, offering practical improvements through regularized, ensemble, and vector-valued extensions.

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 (Portier et al., 2018).

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 $f \in L^2(P)$0 and $f \in L^2(P)$1, where $f \in L^2(P)$2 (Portier et al., 2018).

This representation is equivalent to a multiple linear regression model with intercept $f \in L^2(P)$3: $f \in L^2(P)$4 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

$f \in L^2(P)$5

and the OLS slope and intercept are

$f \in L^2(P)$6

$f \in L^2(P)$7

In the centered formulation, with $f \in L^2(P)$8, the normal equations become

$f \in L^2(P)$9

The zero-variance condition is the exact-fit case. If $I = E[f(X)] = P(f).$0 lies exactly in $I = E[f(X)] = P(f).$1, meaning that there exists $I = E[f(X)] = P(f).$2 with $I = E[f(X)] = P(f).$3 almost surely, then $I = E[f(X)] = P(f).$4 and

$I = E[f(X)] = P(f).$5

In that case the OLS Monte Carlo rule integrates constants and the control variates exactly, and the estimator is exact for every sample set (Portier et al., 2018).

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 $I = E[f(X)] = P(f).$6, the Langevin–Stein operator acting on a sufficiently smooth vector field $I = E[f(X)] = P(f).$7 is

$I = E[f(X)] = P(f).$8

and under suitable boundary and regularity conditions,

$I = E[f(X)] = P(f).$9

A common scalar-function specialization sets $$\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].$$0, giving the second-order Langevin Stein operator

$$\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].$$1

The resulting function $$\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].$$2 has zero mean under $$\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].$$3, so it can be added to the integrand without changing the target expectation (South et al., 2018).

This operator-theoretic construction makes the “zero-variance” principle explicit. If

$$\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].$$4

for a Stein operator $$\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].$$5 and a test function $$\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].$$6 in the Stein class, then $$\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].$$7 is a zero-mean control variate and the adjusted estimator equals the constant $$\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].$$8 for every sample set. In this sense, classical ZVCV is the exact representability of $$\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].$$9 in the image of the Stein operator (Sun et al., 2021).

In parametric ZVCV, $h=(h_1,\dots,h_m)^\top$0 is typically restricted to a polynomial family. If $h=(h_1,\dots,h_m)^\top$1 are monomials of total degree $h=(h_1,\dots,h_m)^\top$2, then the control-variate features are

$h=(h_1,\dots,h_m)^\top$3

For $h=(h_1,\dots,h_m)^\top$4, the polynomial is linear and $h=(h_1,\dots,h_m)^\top$5, so the features reduce to score components. The control-variate estimator then becomes

$h=(h_1,\dots,h_m)^\top$6

with $h=(h_1,\dots,h_m)^\top$7 fitted by least squares or weighted least squares when importance or sequential Monte Carlo weights are present (South et al., 2018).

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 (Sun et al., 2021).

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,

$h=(h_1,\dots,h_m)^\top$8

with $h=(h_1,\dots,h_m)^\top$9. The relevant geometric quantity is the leverage function

$P(h_j)=0$0

together with the leverage condition

$P(h_j)=0$1

This condition implies $P(h_j)=0$2 and controls invertibility and stability of empirical Gram matrices (Portier et al., 2018).

Under the leverage condition,

$P(h_j)=0$3

so in particular,

$P(h_j)=0$4

The central limit theorem is governed by a necessary and sufficient Lindeberg condition: $P(h_j)=0$5 When it holds,

$P(h_j)=0$6

Hence practical confidence intervals take the form

$P(h_j)=0$7

The nonstandard feature is the scaling by $P(h_j)=0$8, the regression residual standard deviation. If the linear span of $P(h_j)=0$9 is dense in a function space containing $\beta \in \mathbb{R}^m$0, then

$\beta \in \mathbb{R}^m$1

so $\beta \in \mathbb{R}^m$2. The integration error then shrinks at rate $\beta \in \mathbb{R}^m$3, faster than the classical $\beta \in \mathbb{R}^m$4 Monte Carlo rate. Explicit examples include post-stratification on $\beta \in \mathbb{R}^m$5, univariate Legendre polynomials on $\beta \in \mathbb{R}^m$6, and multivariate tensor-product Legendre systems on $\beta \in \mathbb{R}^m$7, with approximation rates depending on the regularity of $\beta \in \mathbb{R}^m$8 and on dimensionality (Portier et al., 2018).

The exact ZVCV case sits at the edge of this asymptotic theory. If $\beta \in \mathbb{R}^m$9 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 $$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).$$0 with stationary distribution $$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).$$1, the Poisson equation

$$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).$$2

implies that

$$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).$$3

has zero mean under $$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).$$4. If the exact solution $$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).$$5 were available, then $$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).$$6, yielding literal zero variance for ergodic averages. Practical MCMC ZVCV replaces the unknown Poisson solution by a finite span $$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).$$7 and defines control variates

$$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).$$8

For reversible chains, the optimal coefficients admit an explicit finite-dimensional representation, and consistent estimators can be obtained from the same MCMC run (Dellaportas et al., 2010).

This construction is particularly effective for conjugate random-scan Gibbs samplers, where $$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).$$9 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 $\beta_{\mathrm{opt}} = P(hh^\top)^{-1}P(hf),$0 is tractable or can be approximated (Dellaportas et al., 2010).

A second major Bayesian post-processing line uses Stein operators rather than the Markov operator. In this setting, one fits $\beta_{\mathrm{opt}} = P(hh^\top)^{-1}P(hf),$1 by a finite-dimensional space

$\beta_{\mathrm{opt}} = P(hh^\top)^{-1}P(hf),$2

where $\beta_{\mathrm{opt}} = P(hh^\top)^{-1}P(hf),$3 is the Stein operator and $\beta_{\mathrm{opt}} = P(hh^\top)^{-1}P(hf),$4 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 $\beta_{\mathrm{opt}} = P(hh^\top)^{-1}P(hf),$5, the SECF estimator is exact on the polynomial class $\beta_{\mathrm{opt}} = P(hh^\top)^{-1}P(hf),$6. The same work establishes a bias-correction property when the Markov chain is not invariant for the posterior, proving $\beta_{\mathrm{opt}} = P(hh^\top)^{-1}P(hf),$7 convergence under the stated assumptions (South et al., 2020).

For stochastic gradient MCMC, ZVCV can be applied as a post-processing step using the gradient of the log posterior. With a linear polynomial $\beta_{\mathrm{opt}} = P(hh^\top)^{-1}P(hf),$8, the Stein feature reduces to a linear combination of score components, and the adjusted output becomes

$\beta_{\mathrm{opt}} = P(hh^\top)^{-1}P(hf),$9

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 (Baker et al., 2017).

5. High-dimensional, vector-valued, and algorithmic extensions

A central practical difficulty is dimensionality. For polynomial ZVCV of total degree $\beta_{\mathrm{opt}}$0 in dimension $\beta_{\mathrm{opt}}$1, the number of monomials is

$\beta_{\mathrm{opt}}$2

so the number of nonconstant regression covariates is

$\beta_{\mathrm{opt}}$3

This combinatorial growth rapidly makes ordinary least squares unstable or infeasible. Regularized ZV-CV addresses the problem with penalized regression,

$\beta_{\mathrm{opt}}$4

interpolating between LASSO, ridge, and elastic net. A second device, called a priori regularization, restricts the polynomial to a subset of coordinates $\beta_{\mathrm{opt}}$5, reducing the number of features to

$\beta_{\mathrm{opt}}$6

These methods were introduced precisely to stabilize high-dimensional fits and to reduce computational and storage costs (South et al., 2018).

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 $\beta_{\mathrm{opt}}$7 component regressions using only $\beta_{\mathrm{opt}}$8 selected features per component, then aggregates the fitted intercepts. With the second-order Langevin–Stein feature map $\beta_{\mathrm{opt}}$9, each component defines

$f \in L^2(P)$00

and the final estimator is

$f \in L^2(P)$01

with weights given by simple averaging, a second OLS step, or a constrained Markowitz optimization. The paper emphasizes semi-exact selection, under which all base-order monomials are always included; for Gaussian $f \in L^2(P)$02 and polynomial $f \in L^2(P)$03, each component is then zero-variance when the base degree is high enough. Empirically, ensemble ZVCV is reported as competitive with regularized ZVCV in statistical efficiency and substantially faster in runtime (Nguyen et al., 1 Sep 2025).

ZVCV has also been generalized to multiple related integrals. In vector-valued control variates, the integrands $f \in L^2(P)$04 are treated jointly, and the Stein construction is lifted to a matrix-valued reproducing kernel Hilbert space. If

$f \in L^2(P)$05

then the vector control variate

$f \in L^2(P)$06

is componentwise zero mean, and the estimator collapses to $f \in L^2(P)$07 with zero variance jointly. This framework recovers classical ZVCV when $f \in L^2(P)$08 and extends it to multifidelity modelling, thermodynamic integration, and related multi-task settings (Sun et al., 2021).

6. Variational inference, multimodality, and practical limitations

ZVCV has recently been adapted to pathwise gradient estimators in variational inference. If $f \in L^2(P)$09 with $f \in L^2(P)$10 and

$f \in L^2(P)$11

then the pathwise Monte Carlo gradient estimator is

$f \in L^2(P)$12

ZVCV augments this with zero-mean Stein features $f \in L^2(P)$13: $f \in L^2(P)$14 For a Gaussian base $f \in L^2(P)$15, the first-order $f \in L^2(P)$16-space Stein features reduce to

$f \in L^2(P)$17

The optimal coefficients satisfy

$f \in L^2(P)$18

and the zero-variance condition is again exact linear representability of $f \in L^2(P)$19 in the span of the zero-mean features plus a constant. The method is attractive because it requires only sampling from $f \in L^2(P)$20 and knowledge of $f \in L^2(P)$21, but the reported empirical gains are modest in high dimensions, and simply increasing the number of gradient samples $f \in L^2(P)$22 is often more effective in ELBO-versus-wall-clock terms (Ng et al., 2024).

Multimodality exposes a further limitation of standard Stein features. For multimodal targets, oscillatory or mode-local Stein features can fail to capture cross-mode offsets in the integrand. A recent construction therefore introduces density-ratio-based zero-mean features. If $f \in L^2(P)$23 is a reference distribution approximating the modes of $f \in L^2(P)$24, if $f \in L^2(P)$25 satisfies $f \in L^2(P)$26, and if

$f \in L^2(P)$27

then $f \in L^2(P)$28. The proposed features use a mixture $f \in L^2(P)$29 and basis functions of the form

$f \in L^2(P)$30

The paper reports that combining these density-ratio features with Stein-based features reduces variance more effectively than either family alone on a bimodal example. This suggests that the adequacy of the zero-mean feature class, rather than the control-variate principle itself, is the crucial issue in strongly multimodal regimes (Yamashita et al., 4 Jun 2026).

Across these variants, the recurring limitations are consistent. ZVCV requires zero-mean feature constructions whose expectations are known or guaranteed by Stein identities, regularity and boundary conditions sufficient for the zero-mean property, and numerically stable regression or approximation. High leverage, ill-conditioned Gram matrices, weak correlation between the integrand and the feature class, or biased reuse of the same samples for fitting and estimation can degrade performance. Computational trade-offs also matter: increasing the number of control variates can or cannot be more efficient than increasing the Monte Carlo sample size, depending on approximation quality, leverage, and the cost of building and solving the regression problem (Portier et al., 2018).

7. Conceptual synthesis and relation to adjacent methods

ZVCV can be understood as a unifying principle rather than a single algorithm. In the OLS Monte Carlo framework, it is the exact-fit case of regression control variates. In the Stein framework, it is exact representability of the centered integrand in the image of a Stein operator. In the MCMC Poisson-equation framework, it is exact solution of the Poisson equation. In kernel, vector-valued, variational, regularized, ensemble, and multimodal extensions, the same principle persists: construct a rich class of zero-mean functions, fit the integrand or gradient against that class, and use the residual as the variance-reduced quantity.

Several neighboring methods can be viewed as structured relaxations of classical ZVCV. Control functionals replace finite polynomial spans by RKHS-based Stein images; semi-exact control functionals combine finite-dimensional exactness with kernel interpolation; vector-valued control variates share information across related tasks; regularized and ensemble ZVCV modify the regression stage to remain usable when the polynomial basis is over-parameterized; and pathwise-gradient ZVCV moves the construction from expectation estimation to stochastic-gradient estimation. A plausible implication is that contemporary work on ZVCV is less concerned with the zero-variance ideal in isolation than with designing feature spaces, operators, and solvers that make the ideal numerically approachable in realistic Bayesian and Monte Carlo pipelines (Nguyen et al., 1 Sep 2025).

Under exact representability, ZVCV yields zero variance. Under dense approximation and suitable leverage or regularity conditions, it yields faster-than-$f \in L^2(P)$31 Monte Carlo error. When neither exactness nor strong approximation is available, it remains a regression-based variance-reduction method whose effectiveness is governed by approximation error, conditioning, and computation. This combination of exactness theory, approximation theory, and regression geometry is the defining feature of the ZVCV literature.