Everybody loves the Gaussian distribution (or if you prefer, the Normal distribution) to bits. It’s easy to see why: The Gaussian distribution is easy to handle, there are tons of closed-form results, and if you get a large enough sample size, everything turns Gaussian!

Because we love Gaussians so much, we love to sample from them. But how do you sample from a multivariate Gaussian? It doesn’t look like you can easily find its inverse function:

\(\mathcal{N}(\mu, \Sigma) = \frac{1}{\sqrt{(2\Pi)^d|\Sigma|}}e^{-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)}\)

In fact, it looks pretty difficult. So I asked myself, what really happens when you invoke rnorm(mu, sigma)?

Thankfully, because the Gaussian is so old and well-studied, there exist really good tutorial papers on this kind of stuff. Essentially, sampling from a multivariate Gaussian is done in three stages:

1. Generate \(n\) random variables \(x_i\). Then the central limit theorem allows us to construct samples from \(\mathcal{N}(0,1)\).
2. Repeat step 1. \(d\) times to get a \(d\)-dimensional vector \(x\). This is a sample from \(\mathcal{N}(0,I_d)\).
3. Now for the tricky part: Use the eigenvectors and eigenvalues of the covariance matrix \(\Sigma\) that we wish to sample from to transform the \(\mathcal{N}(0,I_d)\) sample into a \(\mathcal{N}(\mu,\Sigma)\) sample.

The last step bears some explanation: If \(\Phi\) is the matrix of normalised eigenvectors of \(\Sigma\), and \(Delta\) is a diagonal matrix with the eigenvalues along the diagonal, then:

\[y = \Delta^{\frac{1}{2}}\Phi x +\mu\]

is a sample from \(\mathcal{N}(\mu,\Sigma)\).

The above is taken from this excellent review by Istvan Hernadvolgyi (there are at least three accents in this name that I did not render properly, and for that, I apologise).