2
$\begingroup$

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?

$\endgroup$
3
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.