Implicit Control Variate
11 min read
We want the average of a large population but can afford to inspect only a small sample. One way to choose that sample is to flip a coin for each item. Now both the selected values and the number selected fluctuate, and an estimate normalized by the expected sample count pays for both sources of noise.
Subtracting a common baseline before sampling can reduce this noise. The best baseline is the population average itself, which is exactly the number we do not know. Choosing a fixed number of distinct items removes the need to know it: the baseline correction cancels on every sample, including for the best baseline. The resulting variance depends on differences between the values rather than their common level. This cancellation is a special case of the classical sampling variance formulas ([Sen, ’53], [YG, JRSS B ’53], [HT, JASA ’52]).
Naive Estimation
How noisy is the simplest design, and what is it paying for? First consider independent Bernoulli sampling. 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 inclusion indicator of 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 variance contains the squares of the raw values \(x_i^2\), rather than their squared deviations from the average. When the values cluster around a large common offset, that offset dominates every \(x_i^2\) and inflates the variance without increasing their spread. Subtracting a baseline removes the offset, and a design that fixes its sample size turns out to do the subtraction on its own.
Constant Adjustment
What does subtracting a constant do to that variance? A common fix picks a constant \(c\), estimates the average of the shifted values \(x_i - c\), and adds \(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. Tuning \(c\) to \(\bar{x}\) turns each \(x_i^2\) into \((x_i - \bar{x})^2\), removing the common offset from the variance. The minimizing baseline is the unknown average itself. 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 [Williams, ML ’92], does the latter, subtracting a running average of past returns from each new one.
Collecting the \(c\) terms shows what the baseline does:
\[\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 multiplies \(|S| - m\), the deviation of the sample size from its mean. It is therefore a control variate on the sample size: tuning \(c\) fits a regression of the estimate on how many values the coins kept. A design with a fixed sample size makes this correction identically zero.
Fixing the Size
What happens to the correction when the sample size cannot fluctuate? 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 estimate, and hence its variance, is constant in \(c\).
The Variance Depends Only on Differences
If the baseline cancels, what is the fixed-size variance? The fixed-size variance depends only on the pairwise differences between values.
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 . \]
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) . \]
The bracket can be rewritten as a sum of pairwise 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 converts the raw squares into pairwise differences. Independent coin flips have no such covariance, so the common offset remains in their variance.
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. To assess that variance, we compare it with sampling with replacement below.
The Finite Population Correction
How does that variance compare with the textbook estimator that draws with replacement? 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 sampling without replacement fixes the number of distinct values and makes the baseline cancel. Compared with coin-flip sampling using its optimal baseline \(c = \bar{x}\), its variance differs by a single factor:
\[ \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, and sampling without replacement requires no tuned baseline.
In the first plot, the parabola is the coin-flip variance as \(c\) sweeps across the values, and the dashed line for \(\hat{x}^{\mathrm{wor}}\) equals its minimum. 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 use no baseline, so \(\hat{x}^{\mathrm{coin}}\) sits far above the other two because its variance contains the squares of the raw values. The two fixed-draw designs start together at \(m = 1\) and separate as \(m\) grows; the variance of \(\hat{x}^{\mathrm{wor}}\) falls to zero as the sample approaches the whole population. For a nonconstant population, with-replacement sampling continues to draw duplicates, so its variance remains positive at \(m=n\).
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
Where else does fixing a count cancel an offset? 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 common because every model’s masking game includes 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.
When a design can fix a nuisance statistic exactly, such as the sample size or per-stratum counts, do so. The baseline you would otherwise tune is subtracted for you, at exactly the value you could not have computed.