# Sampling from Multivariate Gaussian distribution in Matlab

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. ```