Implicit Control Variate
10 min read
The average of a finite population is often too expensive to compute exactly, so we take a small sample and use the sample average to estimate it. A small sample is noisy, and the standard fix is to subtract a baseline constant from the values before averaging and add it back at the end. The best baseline is the population average itself, the one number we do not have. Fix how many values you sample, though, and that best baseline subtracts itself. The fact already sits inside the classical Sen–Yates–Grundy variance formula for the Horvitz–Thompson estimator ([Sen, ’53], [YG, JRSS B ’53], [HT, JASA ’52]); we pull out the special case where it is easiest to see.
Naive Estimation
To see what the baseline does, we need the design it corrects and a formula for its variance. A population holds \(n\) real values \(x_1, \dots, x_n\), and the target is their average:
\[ \bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i. \]
The simplest random design flips an independent coin for each value: keep value \(i\) with probability \(m/n\), and call the kept set \(S\). The sample then has size about \(m\): because the coins are independent, the size \(|S|\) is random and equals \(m\) only on average. The estimate divides by the target \(m\) regardless, averaging the survivors as if the size had come out right. Write \(\hat{x}^{\mathrm{coin}}\), tagging the estimator by its design, since two more designs will follow:
\[ \hat{x}^{\mathrm{coin}} = \frac{1}{m} \sum_{i \in S} x_i. \]
Write \(\mathbf{1}[i \in S]\) for the coin that keeps value \(i\), so the estimate is a weighted sum of these indicators. We can check that it is unbiased:
\[\begin{align} \mathbb{E}[\hat{x}^{\mathrm{coin}}] &= \frac1{m} \sum_{i=1}^n x_i \, \mathbb{E}\big[\mathbf{1}[i \in S]\big] \\ &= \frac1{m} \sum_{i=1}^n x_i \cdot \frac{m}{n} \\ &= \frac1{n} \sum_{i=1}^n x_i \\ &= \bar{x}, \end{align}\]
where the first line is linearity of expectation, the second uses that each coin keeps its value with probability \(m/n\), and the third cancels the \(m\) to leave the population average.
Since the coins are independent, the variance of the estimate is the sum of the per-value variances, and each coin is Bernoulli with probability \(m/n\):
\[ \mathrm{Var}(\mathbf{1}[i \in S]) = \frac{m}{n}\Big(1 - \frac{m}{n}\Big). \]
Substituting this gives the variance of the estimate:
\[\begin{align} \mathrm{Var}(\hat{x}^{\mathrm{coin}}) &= \frac1{m^2} \sum_{i=1}^{n} x_i^2 \, \mathrm{Var}(\mathbf{1}[i \in S]) \\ &= \frac1{m^2} \sum_{i=1}^{n} x_i^2 \cdot \frac{m}{n}\Big(1 - \frac{m}{n}\Big) \\ &= \frac{1 - m/n}{m} \cdot \frac{1}{n} \sum_{i=1}^{n} x_i^2, \end{align}\]
where the first line uses independence (no cross terms between coins), the second substitutes the Bernoulli variance, and the third collects the constants.
The formula charges for the squares of the raw values \(x_i^2\), not for how far they spread around the average. When the values cluster around some big number, that shared offset dominates every \(x_i^2\) and inflates the variance, even though it tells us nothing about the average.
Constant Adjustment
Here is a common fix: pick a constant \(c\), estimate the average of the shifted values \(x_i - c\), and add \(c\) back at the end:
\[ \hat{x}^{\mathrm{coin}}_c = c + \frac{1}{m} \sum_{i \in S} (x_i - c). \]
The same argument gives unbiasedness for every \(c\), and the same variance computation goes through with \(x_i^2\) replaced by \((x_i - c)^2\):
\[ \mathrm{Var}(\hat{x}^{\mathrm{coin}}_c) = \frac{1 - m/n}{m} \cdot \frac{1}{n} \sum_{i=1}^{n} (x_i - c)^2 . \]
This is a parabola in \(c\), and a sum of squares \(\sum_i (x_i - c)^2\) is smallest when \(c\) is the average. So tuning \(c\) to \(\bar{x}\) turns each \(x_i^2\) into \((x_i - \bar{x})^2\): the offset drops out, and only the genuine spread around the average is left to pay for. The catch is the minimizer: the best baseline is the answer you are estimating. In practice you spend some of your queries on a pilot sample for a rough \(\bar{x}\), or tune the baseline online. REINFORCE, the classic policy-gradient method in reinforcement learning, does the latter, subtracting a running average of past returns from each new one.
What the baseline actually does is easier to see once the \(c\) terms sit in one place. We split the shifted sum into its two pieces and collect them:
\[\begin{align} \hat{x}^{\mathrm{coin}}_c &= c + \frac1m \sum_{i \in S} x_i - \frac1m \sum_{i \in S} c \\ &= c + \hat{x}^{\mathrm{coin}} - \frac{c \, |S|}{m} \\ &= \hat{x}^{\mathrm{coin}} - \frac{c}{m} \big( |S| - m \big), \end{align}\]
where the first line distributes the sum, the second recognizes the sample average and counts the \(|S|\) copies of \(c\), and the third collects the two \(c\) terms. The baseline never touches the values at all. It multiplies \(|S| - m\), how far the coins overshot a target size they hit only on average. The baseline is therefore a control variate on the sample size, centered at its mean: tuning \(c\) fits a regression of the estimate on how many values the coins kept. So a design whose sample size cannot fluctuate should not need the baseline at all.
Fixing the Size
Fix the size directly: draw the sample uniformly from the \(\binom{n}{m}\) subsets of exactly \(m\) distinct values. This is simple random sampling, the basic form of sampling without replacement. The estimate is the same sample average as before; only the design feeding it \(S\) has changed. Write it \(\hat{x}^{\mathrm{wor}}\). Now \(|S| = m\) on every draw, and the baseline correction from the last display vanishes term by term:
\[\begin{align} \hat{x}^{\mathrm{wor}}_c &= c + \frac1m \sum_{i \in S} (x_i - c) \\ &= c + \hat{x}^{\mathrm{wor}} - \frac{c \, |S|}{m} \\ &= c + \hat{x}^{\mathrm{wor}} - c \\ &= \hat{x}^{\mathrm{wor}} , \end{align}\]
where the second line splits off the sample average and the \(|S|\) copies of \(c\), and the third sets \(|S| = m\) so the two \(c\) terms cancel. The corrected estimate is identically the plain one, for every \(c\). That includes the optimal \(c = \bar{x}\), which you no longer need to know. The parabola has collapsed to a single point.
The Variance Only Sees Differences
What variance does the collapsed parabola sit at? It cannot depend on \(c\), and the formula shows why: the data enters only through how far apart the values are.
Claim:
\[ \mathrm{Var}(\hat{x}^{\mathrm{wor}}) = \frac{1 - m/n}{m} \cdot \frac{1}{2n(n-1)} \sum_{i \ne j} (x_i - x_j)^2 . \]
In words, the variance is half the average squared gap between two distinct population values, scaled by a factor \(\frac{1 - m/n}{m}\) that hits zero once the sample is the whole population. Adding the same \(c\) to every value leaves every gap alone, so the variance does not move, matching the pathwise cancellation above.
We can check the smallest case by hand. With \(n = 2\) values and a sample of \(m = 1\), the estimate is \(x_1\) or \(x_2\) with equal probability, sitting half the gap from its mean on either side, so its variance is the squared half-gap:
\[ \mathrm{Var}(\hat{x}^{\mathrm{wor}}) = \Big( \frac{x_1 - x_2}{2} \Big)^2 = \frac{(x_1 - x_2)^2}{4} . \]
The formula agrees: the leading factor is \(\frac{1 - 1/2}{1} = \frac{1}{2}\), the sum over the two ordered pairs is \(2 (x_1 - x_2)^2\), and \(\frac{1}{2} \cdot \frac{2 (x_1 - x_2)^2}{2 \cdot 2 \cdot 1} = \frac{(x_1 - x_2)^2}{4}\).
This is the simplest case of the Sen–Yates–Grundy formula ([Sen, ’53], [YG, JRSS B ’53]), the variance companion to the Horvitz–Thompson estimator ([HT, JASA ’52]), which covers any fixed-size design, equal probabilities or not. What the claim does not say is whether that variance is any good.
Proof
The only inputs are how often a single value, and a pair of values, land in the sample. Write the estimate as a weighted sum of indicators over the entire population:
\[ \hat{x}^{\mathrm{wor}} = \frac{1}{m} \sum_{i=1}^{n} x_i \, \mathbf{1}[i \in S] . \]
The variance of a weighted sum collects a variance term for each value and a covariance term for each ordered pair:
\[ \mathrm{Var}(\hat{x}^{\mathrm{wor}}) = \frac{1}{m^2} \left( \sum_{i=1}^{n} x_i^2 \, \mathrm{Var}(\mathbf{1}[i \in S]) + \sum_{i \ne j} x_i x_j \, \mathrm{Cov}(\mathbf{1}[i \in S], \mathbf{1}[j \in S]) \right) . \]
We need the two indicator moments. Each indicator is Bernoulli with probability \(m/n\), so its variance is:
\[ \mathrm{Var}(\mathbf{1}[i \in S]) = \frac{m}{n}\Big(1 - \frac{m}{n}\Big) = \frac{m(n-m)}{n^2} . \]
For a pair, once \(i\) takes a seat, \(j\) has to land one of the \(m-1\) remaining seats among the \(n-1\) remaining values, which gives the joint inclusion probability:
\[ \mathbb{E}[\mathbf{1}[i \in S]\,\mathbf{1}[j \in S]] = \Pr(i, j \in S) = \frac{m}{n} \cdot \frac{m-1}{n-1} . \]
Subtracting the product of the means gives the covariance:
\[\begin{align} \mathrm{Cov}(\mathbf{1}[i \in S], \mathbf{1}[j \in S]) &= \frac{m(m-1)}{n(n-1)} - \frac{m^2}{n^2} \\ &= \frac{m\big[(m-1)n - m(n-1)\big]}{n^2(n-1)} \\ &= -\frac{m(n-m)}{n^2(n-1)} , \end{align}\]
where the middle line puts both terms over the common denominator \(n^2(n-1)\) and the last uses \((m-1)n - m(n-1) = -(n-m)\). The covariance is negative (one seat taken is one seat fewer), and it is exactly \(-\frac{1}{n-1}\) times the variance. That shared factor \(\frac{m(n-m)}{n^2}\) pulls out of both sums:
\[ \mathrm{Var}(\hat{x}^{\mathrm{wor}}) = \frac{1}{m^2} \cdot \frac{m(n-m)}{n^2} \left( \sum_{i=1}^{n} x_i^2 - \frac{1}{n-1} \sum_{i \ne j} x_i x_j \right) . \]
Now turn the bracket into differences. Expand \((x_i - x_j)^2 = x_i^2 + x_j^2 - 2 x_i x_j\) and sum over all ordered pairs \(i \ne j\):
\[ \sum_{i \ne j} (x_i - x_j)^2 = \sum_{i \ne j} x_i^2 + \sum_{i \ne j} x_j^2 - 2 \sum_{i \ne j} x_i x_j . \]
In the first sum each \(x_i^2\) appears \(n-1\) times, once for every partner \(j \ne i\), so \(\sum_{i \ne j} x_i^2 = (n-1) \sum_i x_i^2\); the second sum is the same by symmetry. Hence
\[ \sum_{i \ne j} (x_i - x_j)^2 = 2(n-1) \sum_{i=1}^{n} x_i^2 - 2 \sum_{i \ne j} x_i x_j , \]
and dividing by \(2(n-1)\) recovers the bracket exactly:
\[ \sum_{i=1}^{n} x_i^2 - \frac{1}{n-1} \sum_{i \ne j} x_i x_j = \frac{1}{2(n-1)} \sum_{i \ne j} (x_i - x_j)^2 . \]
Substituting this back and simplifying the constants gives the claim:
\[ \mathrm{Var}(\hat{x}^{\mathrm{wor}}) = \frac{1}{m^2} \cdot \frac{m(n-m)}{n^2} \cdot \frac{1}{2(n-1)} \sum_{i \ne j} (x_i - x_j)^2 = \frac{1 - m/n}{m} \cdot \frac{1}{2n(n-1)} \sum_{i \ne j} (x_i - x_j)^2 . \]
The negative covariance is the engine: independent coin flips have none, so their squares never pair up into differences, and the common offset never cancels.
The Finite Population Correction
To see how good that variance is, compare it against the textbook estimator: draw \(m\) values uniformly with replacement and average them, and call it \(\hat{x}^{\mathrm{wr}}\). The number of draws is fixed at \(m\), but the draws are independent, so a value can be drawn twice.
A single uniform draw equals \(x_i\) with probability \(1/n\), so it has mean \(\bar{x}\) and variance \(\sigma^2\), the population variance. We can write \(\sigma^2\) two ways, around the mean and over pairs:
\[ \sigma^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^2 = \frac{1}{2 n^2} \sum_{i \ne j} (x_i - x_j)^2 . \]
Expanding either set of squares gives \(\frac{1}{n} \sum_i x_i^2 - \bar{x}^2\), so the two forms agree.
Averaging \(m\) independent draws divides the variance by \(m\):
\[ \mathrm{Var}(\hat{x}^{\mathrm{wr}}) = \frac{\sigma^2}{m} = \frac{1}{2 n^2 m} \sum_{i \ne j} (x_i - x_j)^2 . \]
Both variances are the same sum over pairwise differences, so dividing one by the other cancels the sum and leaves a single ratio:
\[ \frac{\mathrm{Var}(\hat{x}^{\mathrm{wor}})}{\mathrm{Var}(\hat{x}^{\mathrm{wr}})} = \frac{n - m}{n - 1} \le 1 . \]
That factor is the finite population correction, and it is why \(\hat{x}^{\mathrm{wr}}\) is the weakest of the three: with replacement is the only design that can draw the same value twice, so some of its \(m\) draws are spent re-measuring values it already holds. Sample half the population and the factor is about \(\frac{1}{2}\): with replacement would need twice as many draws to reach the same variance.
The three designs differ in one respect that matters here, whether the sample size can fluctuate:
- Coin flips, \(\hat{x}^{\mathrm{coin}}\): the size \(|S|\) is random, only \(m\) on average.
- With replacement, \(\hat{x}^{\mathrm{wr}}\): a fixed number of draws, but a random number of distinct values.
- Without replacement, \(\hat{x}^{\mathrm{wor}}\): a fixed size of exactly \(m\) distinct values.
Only the last pins the size completely, and only it makes the baseline cancel. What it gives up is a single factor, measured against the coin-flip variance tuned to the bottom of the parabola at \(c = \bar{x}\):
\[ \mathrm{Var}(\hat{x}^{\mathrm{wor}}) = \frac{n}{n-1} \, \mathrm{Var}(\hat{x}^{\mathrm{coin}}_{\bar{x}}) . \]
The gap is one part in \(n-1\), a tenth of a percent for a thousand values. Without replacement sits at the floor of the parabola with nothing tuned.
In the first plot, the parabola is the coin-flip variance as \(c\) sweeps across the values, and the dashed line is \(\hat{x}^{\mathrm{wor}}\), sitting at the bottom of a curve it never had to search. The solid line above it is \(\hat{x}^{\mathrm{wr}}\), at roughly twice the variance, since the plot samples half the population. The comparison changes with the sample size.
In the second plot the coin flips carry no baseline, so \(\hat{x}^{\mathrm{coin}}\) sits far above the other two, still charged for the squares of the raw values. The two fixed-draw designs start together at \(m = 1\) and separate as \(m\) grows, with \(\hat{x}^{\mathrm{wor}}\) diving as the sample approaches the whole population, where it is exact. With replacement keeps re-drawing duplicates and never gets there.
The cancellation never used the fact that the target was an average, only that the design fixes how many \(c\) terms the sample contains.
Applications
The same move helps wherever a fixed linear combination is estimated by sampling. The example I keep running into is the Shapley value: a signed combination of the \(2^d\) coalition values of a game \(f\) on \(d\) players, with coefficients summing to zero, so the target ignores any constant offset. An i.i.d. coalition sampler does not ignore it: the offset adds variance on every draw, and no importance distribution can absorb it, because sampling probabilities cannot be proportional to signed coefficients. Sampling without replacement can: stratify the coalitions so the coefficient is constant within each stratum (by size, and by whether they contain the player), fix the count drawn from each stratum, and the offset cancels on every draw, exactly as \(c\) did above. Games with a large common offset are the norm, since every model’s masking game carries its base prediction.
Paired sampling is the same kind of move, cancelling the entire even part of the game pathwise; see the even–odd decomposition. The tuned single constant appears in the Shapley taxonomy as an estimator of its own.
The general lesson: when a design can pin a nuisance statistic exactly (the sample size, the per-stratum counts), pin it. The baseline you would otherwise tune is subtracted for you, at exactly the value you could not have computed.