# Sample from multivariate normal distribution with given positive-semidefinite covariance matrix

I want to draw a random vector from a multivariate normal distribution with given covariance matrix $Σ$.

I'm following this algorithm:

A widely used method for drawing a random vector $x$ from the $N$-dimensional multivariate normal distribution with mean vector $μ$ and covariance matrix $Σ$ works as follows:

1. Find any real matrix $A$ such that $AA^T = Σ$. When $Σ$ is positive-definite (...)

But my covariance matrix is not positive-definite. It's positive-semidefinite. So, I can't use Cholesky decomposition. Further there is written:

An alternative is to use the matrix $A = UΛ^{\frac{1}{2}}$ obtained from a spectral decomposition $Σ = UΛU^T$ of $Σ$.

So, I have to use spectral decomposition? There is no easier way to draw such a random vector?

For such $Σ$:

 0.666 -0.333
-0.333  0.666


Cholesky decompososition gives such $A$ (using this Online Matrix Calculator):

 0.816  0.000
-0.408  0.707


Spectral decomposition, eigenvectors real values $U$:

 0.707  0.707
-0.707  0.707


Eigenvalues $Λ$:

 0.999  0.000
0.000  0.333


Fortunately square root of diagonal matrix is easy :) $Λ^{\frac{1}{2}}$:

 0.999  0.000
0.000  0.577


$A=UΛ^{\frac{1}{2}}$:

 0.707  0.407
-0.707 0.407


So, $A$ returned by spectral decomposition is different than this returned by Cholesky. Why? Is this correct?

They both seem to work quite well. They give slightly different results for the estimated covariance matrices of the generated series, but I wouldn't be surprised if it's due to rounding errors somewhere in the computations.

Below is some R code which generates samples from $N(0, \boldsymbol{\Sigma})$.

n <- 10000000

X <- cbind(rnorm(n), rnorm(n))

sigma <- t(matrix(c(0.666, -0.333, -0.333, 0.666), nrow=2))
spectral <- eigen(sigma)

X.spectral <- t(spectral$vectors %*% sqrt(diag(spectral$values)) %*% t(X))
X.cholesky <- t(t(chol(sigma)) %*% t(X))
cov(X.spectral)
cov(X.cholesky)


So with my 10,000,000 samples, the covariance matrices are

> cov(X.spectral)
[,1]       [,2]
[1,]  0.6660626 -0.3331138
[2,] -0.3331138  0.6658130
> cov(X.cholesky)
[,1]       [,2]
[1,]  0.6660344 -0.3328923
[2,] -0.3328923  0.6656198