# 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:

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];
A.z
]


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

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