# Mathematica: Faster Multinormal Sampling with Different Covariance Matrices

Let's say we are given a set of different covariance matrices. We wish to draw a single sample of the multivariate normal distribution from each of these covariance matrices.

Assuming the mean vector is zero, we would code this
for a single draw from a single covariance matrix
`sigma` in Mathematica:

MultinormalDraw[sigma_] := RandomVariate[MultinormalDistribution[sigma]]

To demonstrate why this is slow, we will use a hundred randomly generated $200\times 200$ covariance matrices (that is, matrices that are both symmetric and positive definite):

sigmas = Table[ M = RandomVariate[NormalDistribution[], {200, 200}]; M.Transpose[M], 100 ]; AllTrue[sigmas, PositiveDefiniteMatrixQ[#] && SymmetricMatrixQ[#]&] (* => True *)

Now we time this on `sigmas`:

Map[MultinormalDraw, sigmas] // AbsoluteTiming // First (* => 5.6585 *)

## Solution

In this case it is actually faster to implement a multinormal sampler from scratch. The procedure is simple. Let $\Sigma$ be a $N\times N$ covariance matrix. Then, according to Wikipedia:

- Find a real matrix $A$ such that $AA^T = \Sigma$, typically the Cholesky decomposition of $\Sigma$.
- Sample $N$ independent standard normal variates and put them in a vector $z$.
- Calculate $\mu + Az$, this vector is multivariate normally distributed with mean $\mu$ and covariance $\Sigma$.

Mathematica has already provided us all the ingredients we need. We write (assuming $\mu = 0$):

FasterDraw[sigma_] := Block[{A, z}, A = Transpose[CholeskyDecomposition[sigma]]; z = RandomVariate[NormalDistribution[], Length@sigma]; A.z ]

If we map this function over the list `sigmas`, we
obtain a significant speedup:

Map[FasterDraw, sigmas] // AbsoluteTiming // First (* => 0.109244 *)

## Comments

Comments powered by Disqus