###### tags: `vrmc` `monte carlo` `expository`

# Variance Reduction in Monte Carlo Methods: Bilinear Strategies

**Overview**: In this note, I will describe the challenge of variance reduction in the context of Monte Carlo methods, and describe two methods for accomplishing this goal, based around the notion of *bilinearity*.

## Estimating Integrals using Monte Carlo

The premise of the Monte Carlo method is that

1. Many quantities of interest can be written as integrals with respect to some probability measure.
2. Drawing samples from probability measures is sometimes simple.
3. Using samples from probability measures, one can estimate the value of integrals against the underlying measure.

More concretely, let $I$ be some quantity of interest, and suppose that we can write

\begin{align}
I = \int f(x) p(dx)
\end{align}

where $f$ is some function which we can evaluate, and $p$ is some probability measure from which we can draw samples. Sometimes, this representation will come directly from a probabilistic interpretation of the quantity at hand, as is common in Probability and Statistics. In other examples, the probabilistic interpretation will be imposed artificially in order to enable this methodology.

With this representation available to us, we can write that $I = \mathbf{E}_p [f]$, and then proceed to

1. Draw samples $x^1, \ldots, x^N \sim p$
2. Form the estimator

\begin{align}
\hat{I}_N = \frac{1}{N} \sum_{i = 1}^N f(x^i).
\end{align}

The law of large numbers guarantees that under reasonable conditions, $\lim_{N \to \infty} \hat{I}_N = I$, and other theorems can make this result quantitative. For example, if the function $f$ has finite variance under $p$, then one can use Chebyshev's inequality to say that

\begin{align}
\mathbf{P} ( |\hat{I}_N - I| \geqslant \varepsilon) \leqslant \frac{\text{Var}_p (f)}{n \varepsilon^2}.
\end{align}

This result suggests that one should typically expect that

\begin{align}
\hat{I}_N \approx I + O_p \left(\sqrt{\frac{\text{Var}_p (f)}{N}} \right),
\end{align}

and other results like the Central Limit Theorem reassure us that this is roughly the best rate of convergence that we can expect.

## Variance Reduction

On one hand, this is good news: as $N$ grows, we have that $\hat{I}_N \to I$, and so the problem of approximating $I$ can to some extent be solved through simulation. On the other hand, the $N^{-1/2}$ rate does not appear to be particularly fast, and moreover, it can be worsened by the presence of a large variance term $\text{Var}_p (f)$.

With this in mind, a natural strategy is to try to reduce the magnitude of this variance term. For fixed $(p, f)$, there is not much to be done: the variance is just a number. However, our real task is to estimate $I$, and we can happily replace $(p, f)$ with some alternative pair $(p', f')$ which also satisfies $I = \mathbf{E}_{p'} [f']$. So, how might we construct such alternative pairs?

## Bilinearity for Variance Reduction

A basic observation in the context of Monte Carlo is to note that the mapping

\begin{align}
(p, f) \mapsto \langle p, f \rangle := \int p(dx) f(x)
\end{align}

is bilinear in its inputs. By definition, this tells us that (with some minor caveats)

* For a probability measure $p$, and functions $f, g$, it holds that

\begin{align}
\langle p, f + g \rangle = \langle p, f \rangle + \langle p, g \rangle.
\end{align}

* For probability measures $p, q$, and a function $f$, it holds that

\begin{align}
\langle p + q, f \rangle = \langle p, f \rangle + \langle q, f \rangle.
\end{align}

In the following sections, I will give two examples of variance reduction strategies which can be viewed as algorithmic implementations of these observations.

## Control Variates

For our first example, fix a probability measure $p$ and a function $f$, and consider the problem of estimating $\langle p, f \rangle$. Often, when we write down a $p$ from which we can sample, we have some extra information available to us about the measure. If we know enough about a measure that we are able to draw samples from it, then it is reasonably likely that we know something about where the measure is supported, centered, or similar. A common form of information which we might know about a measure is its *moments*, that is, we know some function $\phi : \mathcal{X} \to \mathbf{R}^F$ such that

\begin{align}
\int p(dx) \phi(x) = 0.
\end{align}

For example, if $\phi(x) = x - \mu$, then we know that $x$ has mean $\mu$ under $p$.

Information of this form tells us that $\langle p, \phi \rangle = 0$. It is typical to refer to such $\phi$ as a *control variate* for $p$.

In light of our earlier observation, we now have something new to work with: by bilinearity, we know that for any $\beta \in \mathbf{R}^F$, it holds that

\begin{align}
\langle p, f - \beta^T\phi \rangle &= \langle p, f \rangle - \langle p, \beta^T \phi \rangle \\
&= \langle p, f \rangle - \beta^T\langle p, \phi \rangle \\
&= \langle p, f \rangle - \beta^T0_F \\
&= \langle p, f \rangle \\
&= I.
\end{align}

With this in mind, we now have a whole family of estimators which we can use to estimate $I$, indexed by $\beta \in \mathbf{R}^F$:

\begin{align}
\hat{I}_N ( \beta ) = \frac{1}{N} \sum_{i = 1}^N \left( f(x^i) - \beta^T \phi (x^i) \right).
\end{align}

We should thus choose $\beta$ such that the variance of the resulting estimator is minimal. Some routine calculations assure us that we should take $\beta$ as

\begin{align}
\beta_* = \mathbf{Cov}_p ( \phi (x))^{-1} \mathbf{Cov}_p ( \phi(x), f(x)),
\end{align}

and one can also prove that $\text{Var}_p \left( \hat{I}_N ( \beta_* ) \right) \leqslant \text{Var}_p \left( \hat{I}_N \right)$, i.e. that using control variates properly does not make matters worse. Depending on the application, control variates can make matters much better.

In practice, $\beta_*$ is usually not known analytically, and so one has to estimate it from samples as well. If the dimension of $\phi$ is large relative to the number of samples, some extra care is typically required.

Another challenge is the *ab-initio* identification of valid control variates $\phi$. This is highly dependent on the nature of the measure $p$, but there are a number of practically relevant scenarios in which such $\phi$ may be systematically constructed.

## Mixture Decompositions

Having used bilinearity on one side of the pairing to improve the variance of an estimator, we now turn to the other side.

Suppose now that we are working with a probability measure $p$ which admits a decomposition into a mixture of $K$ other probability measures, i.e. that we can write

\begin{align}
p(dx) = \sum_{k = 1}^K w_k p_k (dx),
\end{align}

where each $p_k$ is a probability measure, and the weights $w_k$ are nonnegative and sum to $1$. Sometimes, a probability measure is explicitly specified in this way; at other times, one can identify such a decomposition by inspection, considerations of symmetry, or partitioning of the underlying space.

In this case, bilinearity tells us that

\begin{align}
\langle p, f \rangle &= \left\langle \sum_{k = 1}^K w_k p_k, f \right\rangle \\
&= \sum_{k = 1}^K w_k \langle p_k, f \rangle.
\end{align}

As such, we might consider a divide-and-conquer strategy, i.e. draw samples from each of the $p_k$, use them to estimate $\langle p_k, f \rangle$, and then aggregate these estimates to estimate $\langle p, f \rangle$. We thus define the following estimator

1. For $k = 1, \ldots, K$
2. Draw samples $x_k^1, \ldots, x_k^{N_k} \sim p_k$
3. Form the preliminary estimator

\begin{align}
\hat{I}_{k, N_k} = \frac{1}{N_k} \sum_{i = 1}^N f(x_k^i).
\end{align}

3. Aggregate the preliminary estimators into the final estimator

\begin{align}
\hat{I}_{N_1,\ldots,N_K} = \sum_{k = 1}^K w_k \hat{I}_{k, N_k}.
\end{align}

This is known as *stratified sampling*. The interpretation is that the original measured is partitioned into $K$ distinct 'strata', within which we sample independently.

As before, it is in our interests to choose the parameters of our new estimator so as to minimise its variance. Here, the freedom is in how we allocate our sampling budget to each of the $p_k$.

Suppose that we simulate a total of $N$ samples, and we set $N_k = \alpha_k \cdot N$, with the $\alpha_k$ nonnegative and summing to $1$. Standard calculations allow us to calculate that

\begin{align}
\text{Var} \left( \hat{I}_{N_1,\ldots,N_K} \right) = N^{-1} \sum_{k = 1}^K w_k \alpha_k^{-1} \text{Var}_{p_k} (f).
\end{align}

Now, multiply this expression by $1 = \sum_{k = 1}^K \alpha_k$ and apply the Cauchy-Schwartz Inequality to note that

\begin{align}
\sum_{k = 1}^K w_k \alpha_k^{-1} \text{Var}_{p_k} (f) &= \left( \sum_{k = 1}^K w_k \alpha_k^{-1} \text{Var}_{p_k} (f) \right) \cdot \left( \sum_{k = 1}^K \alpha_k \right) \\
&\geqslant \left( \sum_{k = 1}^K \sqrt{w_k \alpha_k^{-1} \text{Var}_{p_k} (f)} \cdot \sqrt{\alpha_k}\right)^2 \\
&= \left( \sum_{k = 1}^K \sqrt{w_k \text{Var}_{p_k} (f)} \right)^2
\end{align}

with equality holding when $\alpha_k \propto \sqrt{w_k \cdot \text{Var}_{p_k} (f)}$. With some extra care, one can check that this is indeed the optimal choice of $\alpha$, and the resulting estimator is never worse than the original estimator. Depending on how imbalanced the values of $w_k \cdot \text{Var}_{p_k} (f)$ are, there is potential for great improvement.

One very simple setting to which stratified sampling applies is when $p$ admits symmetries. For example, suppose that $p$ is a probability measure on the real line, with a density which satisfies $p(x) = p(-x)$. Here, one can decompose $p = \frac{1}{2} p_1 + \frac{1}{2} p_2$, where

\begin{align}
p_1 (dx) &= 2 \cdot p(dx) \cdot \mathbf{I} [ x < 0 ] \\
p_2 (dx) &= 2 \cdot p(dx) \cdot \mathbf{I} [ x \geqslant 0 ].
\end{align}

Applying stratified sampling to this decomposition corresponds to recognising that $p$ places equal mass on both halves of the real line, and using that knowledge to reproduce this feature in our simulated sample. This is related to the technique of *antithetic* sampling, which goes a step further and generates the samples from $p_1$ and $p_2$ in a dependent fashion.

## Conclusion

In Monte Carlo, a key factor governing convergence rates is the variance of the estimator at hand. As such, a natural strategy for improving performance is to devise an estimator with lower variance. In this note, I have demonstrated that by simple noting the bilinearity of the key pairing which enables the Monte Carlo method, one can derive some simple strategies for variance reduction.

In the two strategies presented above, the dimension of the control variate ($F$) and the size of the decomposition ($K$) are treated as fixed. In some applications, it is appropriate to treat them as growing with $N$. In these scenarios, it is sometimes possible to achieve a faster rate of convergence, in the sense that $\text{Var}(\hat{I}^*_N) = o(N^{-1})$ as $N$ grows. The interested reader is encouraged to study some examples of this form.