Skip to main content


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_] :=

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}];

AllTrue[sigmas, PositiveDefiniteMatrixQ[#] && SymmetricMatrixQ[#]&]
=> True

Now we time this on sigmas:

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


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:

  1. Find a real matrix $A$ such that $AA^T = \Sigma$, typically the Cholesky decomposition of $\Sigma$.
  2. Sample $N$ independent standard normal variates and put them in a vector $z$.
  3. 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];

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

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


Comments powered by Disqus