Monte Carlo Estimation and Variance Reduction

1. The Challenge: Intractable Expectations and the Curse of Dimensionality

In many areas of machine learning, from Bayesian inference to reinforcement learning, we encounter the need to compute an expectation of a function $f(\mathbf{x})$ with respect to a probability distribution $p(\mathbf{x})$:

A prime example, directly following our last lesson, is computing the Bayesian evidence (or marginal likelihood), $p(X)=\int p(X|\text{vtheta})p(\text{vtheta})d\text{vtheta}$, which is an expectation of the likelihood with respect to the prior.

For simple, one-dimensional problems, we might use classical numerical integration methods like the trapezoidal rule or Simpson's rule. These methods work by discretizing the integration domain into a grid. However, this approach collapses completely in high dimensions—a phenomenon often called the curse of dimensionality. If our data $\mathbf{x}$ lives in a $D$-dimensional space and we want to place just 10 grid points along each dimension, we would need to evaluate the function at $10^D$ locations. For even a modest $D=10$, this is $10^{10}$ evaluations, and for a high-resolution image, $D$ can be in the millions. This makes any form of grid-based approximation computationally impossible.

The following video provides a concise overview of why direct computation and grid approximations fail, setting the stage for a sampling-based approach.

2. The Monte Carlo Estimator: An Elegant Solution

The Monte Carlo method sidesteps the curse of dimensionality with a simple yet profound idea: instead of trying to evaluate the integrand everywhere, we can approximate the integral by taking a random sample of points and averaging the function's value at those points.

Formally, to estimate the expectation ${\mathbb{E}}_{\mathbf{x}\sim p(\mathbf{x})}[f(\mathbf{x})]$, we draw $N$ independent and identically distributed (i.i.d.) samples $\{{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_N\}$ from the distribution $p(\mathbf{x})$. The Monte Carlo estimator of the expectation is then simply the sample mean:

The power of this method comes from the law of large numbers, which guarantees that as $N\to \infty$, this sample mean converges to the true expectation. More importantly for our purposes, this estimator is unbiased. This means that for any finite $N$, the expected value of our estimator is equal to the true value we are trying to estimate.

Let's prove this. The samples ${\mathbf{x}}_i$ are i.i.d. from $p(\mathbf{x})$, so $\mathbb{E}[f({\mathbf{x}}_i)]$ is the same for all $i$ and is equal to the true expectation, which we'll call $I$. Using the linearity of expectation:

The following resource from the excellent "Physically Based Rendering" textbook formalizes this core estimator.

3. Quantifying the Error: The Variance of the Estimator

Knowing our estimator is unbiased is crucial, but it doesn't tell us how accurate a specific estimate from a finite number of samples $N$ is likely to be. For that, we need to analyze its variance. The variance tells us how much the estimator ${\hat{F}}_N$ "bounces around" its mean value (the true integral $I$) over different sets of random samples.

Let's denote the variance of the function $f(\mathbf{x})$ itself as ${\sigma }^2=\text{Var}(f(\mathbf{x}))$. Since the samples ${\mathbf{x}}_i$ are independent, the variance of the sum is the sum of the variances. We can derive the variance of our estimator ${\hat{F}}_N$ as follows:

This is a fundamental result. The variance of the Monte Carlo estimator decreases linearly with the number of samples $N$.

The standard error of the estimator, which is the standard deviation of its sampling distribution, is the square root of the variance:

This tells us that to halve the error, we must quadruple the number of samples. This $O(1/\sqrt{N})$ convergence rate might seem slow, but its most important feature is that it is completely independent of the dimension $D$ of the integration space. This is why Monte Carlo methods are the tool of choice for high-dimensional integration.

In practice, we don't know the true variance ${\sigma }^2$, but we can estimate it from our samples using the sample variance:

This allows us to estimate the error of our Monte Carlo result without knowing the true answer.

The following reading elaborates on this error analysis.

4. A Glimpse of Variance Reduction: Importance Sampling

While the $1/\sqrt{N}$ convergence rate is a great asset, the constant term ${\sigma }^2$ can be a major problem. If the variance of $f(\mathbf{x})$ is very large, we might need an astronomical number of samples to achieve acceptable accuracy. This often happens when the function $f(\mathbf{x})$ has large values in regions where the probability $p(\mathbf{x})$ is small. Plain Monte Carlo will rarely sample these "important" but rare regions, leading to a high-variance estimate.

Importance Sampling is a powerful variance reduction technique that addresses this. The core idea is to stop sampling from $p(\mathbf{x})$ and instead use a different, "proposal" distribution $q(\mathbf{x})$ that is chosen to over-sample the important regions. To correct for this change of distribution, we re-weight the samples.

The derivation is simple:

The new estimator is:

The terms $w({\mathbf{x}}_i)=p({\mathbf{x}}_i)/q({\mathbf{x}}_i)$ are called the importance weights. By choosing $q(\mathbf{x})$ to be large where $|f(\mathbf{x})p(\mathbf{x})|$ is large, we can dramatically reduce the variance of the estimator. This is a vast topic, but the following video gives a great conceptual introduction.

Conclusion

In this lesson, you've added a cornerstone of modern computational statistics and machine learning to your toolkit. You now have a principled method for tackling the intractable integrals that stymie purely analytical approaches.

Key Takeaways:

• Monte Carlo (MC) estimation approximates intractable expectations by averaging function evaluations over random samples drawn from the relevant probability distribution.
• The standard MC estimator, ${\hat{F}}_N=\frac{1}{N}\sum f({\mathbf{x}}_i)$, is unbiased, meaning its expected value is the true expectation.
• The estimator's error, quantified by the standard error, decreases at a rate of $O(1/\sqrt{N})$. Crucially, this convergence rate is independent of the problem's dimensionality, which is the method's primary advantage over grid-based techniques.
• The variance of the estimator is ${\sigma }^2/N$, where ${\sigma }^2$ is the variance of the function being evaluated. This variance can be estimated from the samples themselves, providing a measure of confidence in our result.
• Techniques like Importance Sampling can significantly reduce the variance of the estimator by intelligently choosing where to draw samples, leading to faster convergence.

We have established how to estimate expectations. But in deep learning, we don't just want to estimate values; we want to optimize them. What happens if an expectation is part of a loss function that we need to minimize via gradient descent? This requires differentiating through the Monte Carlo estimation process.

In our next lesson, we will derive the pathwise gradient estimator, more famously known as the reparameterization trick. This is the key technique that allows us to compute low-variance gradients for stochastic models, unlocking our ability to train a Variational Autoencoder.