calvin mccarter

tl;dr: Don’t use mvnrnd in Matlab for large problems; do it manually instead.

The first improvement uses the Cholesky decomposition, allowing us to sample from a univariate normal distribution. The second improvement uses the Cholesky decomposition of the sparse inverse covariance matrix, not the dense covariance matrix. The third improvement avoids computing the inverse, instead solving a (sparse) system of equations.

n = 10000;
Lambda = gallery('tridiag',n,-0.3,1,-0.3); % sparse
tic;
x_mvnrnd = mvnrnd(zeros(n,1),inv(Lambda));
toc;

tic;
z = randn(n,1); % univariate random
Sigma = inv(Lambda);
A = chol(Sigma,'lower'); % sparse
x_fromSigma = A*z; % using cholesky of Sigma
toc;

tic;
z = randn(n,1); % univariate random
L_Lambda = chol(Lambda,'lower'); % sparse
A_fromLambda = (inv(L_Lambda))'; % sparse
x_fromLambda = A_fromLambda*z;
toc;

tic;
z = randn(n,1); % univariate random
L_Lambda = chol(Lambda,'lower'); % sparse
x_fromLambda = L_Lambda'\z;
toc;

Results:

Elapsed time is 4.514641 seconds.
Elapsed time is 2.734001 seconds.
Elapsed time is 1.740317 seconds.
Elapsed time is 0.012431 seconds.

---

This was originally published here: https://calvinmccarter.wordpress.com/2015/01/06/multivariate-normal-random-number-generation-in-matlab/