9

I'm looking to generate correlated random variables. I have a symmetric, positive definite matrix. So I know that you can use the Cholesky decomposition, however I keep being told that this only works for Gaussian random variables?! Is that true?

Furthermore how does this compare to Eigen decomposition. For example using Cholesky decomposition we can write a random parameter as:

$x = \bar{x} + Lz$

where $L$ is the Cholesky decomposition (lower/upper triangular matrix) and $z$ is some vector of random variables. So one can sample the $z$'s and build up a pdf of x. Now we could also use Eigen decomposition and write x as:

$x = \bar{x} + U\lambda^{1\over2}z$

where $\lambda$ is a diagonal matrix of eigenvalues and $U$ is a matrix composed of the eigenvalues. So we could also build a pdf of this. But if we equate these $x$'s we find that $L = U\lambda^{1\over2}$ But this isn’t true as $L$ is triangular and $U\lambda^{1\over2}$ is not?! So I'm really, really confused. So to clarify the questions:

1) For Cholesky decomposition does the vector z have to be only Gaussian? 2) How does the eigenvalue compare with the Cholesky decomposition? They are clearly different factorisation techniques. So I don't see how the $x$'s above can be equivalent?

Thanks, as always, guys.

user2350366
  • 674
  • 5
  • 12
  • (1) No, any distribution. (2) Hmm, I'm not mathematician. As I know the two decompositions are indeed quite different. `Composed of the eigenvalues` Did you mean "eigenvectors"? BTW, it is possible to generate the variables you want also via principal components, with the help of eigenvalues and eigenvectors. But, sorry, I failed to understand what you imply by the two formulas you showed. – ttnphns Jan 30 '14 at 17:20
  • Take a look at [this blog post](http://sciencemeanderthal.wordpress.com/2012/06/28/cholesky-decomposition-of-variance-covariance-matrices-in-the-classic-twin-study). Hope this helps. – Aleksandr Blekh Dec 17 '14 at 02:23
  • It's a theorem that if $x$ has a multivariate normal distribution, and we multiply $y=Ax$, then $y$ will also have a multivariate normal distribution. However, if $x$ has some other distribution then the distribution of $y$ need not be of the same family as the distribution of $x$. However, the method that you've described will ensure that you get the desired covariance. – Brian Borchers Feb 22 '15 at 22:12
  • How about reading [my answer](http://stackoverflow.com/a/40316331/4891738) at Stack Overflow: [Obtain vertices of the ellipse on an ellipse covariance plot (created by `car::ellipse`)](http://stackoverflow.com/a/40316331/4891738). Although the question is asked in different application, the theory behind is the same. You will see nice figures for geometric explanation there. – Zheyuan Li Nov 04 '16 at 11:22

1 Answers1

7

1) Pretty much yes. The reason is that the $x_i$'s are going to end up being a linear combination of the $z_i$'s. That works out nicely for Gaussian deviates because any linear combination of Gaussian deviates is, itself, a Gaussian deviate. Unfortunately, this is not necessarily true of other distributions.

2) It's a little puzzling, I know, but they are equivalent. Let $\Sigma$ be your covariance matrix and suppose you have both the Cholesky factorization, $\Sigma=L L^T$ and the eigendecomposition, $\Sigma=U \lambda U^T$. The covariance of $L z$ is given by: $$ \begin{array}{} E[L z (L z)^T] & = & E[L z z^T L^T] \\ & = & L \ E[z z^T] \ L^T \\ & = & L \ I \ L^T \\ & = & L L^T \\ & = & \Sigma \end{array} $$ Similarly, the covariance of $U \lambda^\frac{1}{2} z$ is given by: $$ \begin{array}{} E[U \lambda^\frac{1}{2} z (U \lambda^\frac{1}{2} z)^T] & = & E[U \lambda^\frac{1}{2} z z^T \lambda^\frac{1}{2} U^T] \\ & = & U \lambda^\frac{1}{2} \ E[z z^T] \ \lambda^\frac{1}{2} U^T \\ & = & U \lambda^\frac{1}{2} \ I \ \lambda^\frac{1}{2} U^T \\ & = & U \lambda^\frac{1}{2} \lambda^\frac{1}{2} U^T \\ & = & U \lambda U^T \\ & = & \Sigma \end{array} $$ For purposes of computation, I suggest you stick with the Cholesky factorization unless your covariance matrix is ill-conditioned/nearly singular/has a high condition number. Then it's probably best to switch to the eigendecomposition.

Sycorax
  • 76,417
  • 20
  • 189
  • 313
Bill Woessner
  • 1,168
  • 11
  • 15
  • 2
    Oh, I should add that this technique works for ANY factorization of $\Sigma = F F^T$. And this factorization is not unique unless you force $F$ to be lower triangular (in which case, $F$ is the Cholesky factor). – Bill Woessner Oct 30 '15 at 14:54
  • Nice answer which is spot on together with the comment. – chichi Jul 21 '21 at 03:21