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