2
$\begingroup$

How to sample from a multivariate normal given the $P^TLDL^TP$ decomposition of $\Sigma$, and $\mu$?

Here, $P$ is a permutation matrix, $L$ is lower-triangular, and $D$ is diagonal

Given the Cholesky decomposition $LL^T$ of $\Sigma$, and $\mu$, we can sample by doing:

$sample = \mu + L * z$

... where $z$ is a vector of univariate standard normals.

I imagine that there must be some way of tweaking this formula to work with the $P^TLDL^TP$ decomposition, but my maths isn't quite good enough to see how to do that?

$\endgroup$
  • 3
    $\begingroup$ Perhaps your notation is confusing you, so let's write the decomposition as $P'\Lambda D\Lambda'P = LL'$. Now isn't it obvious that $L = P'\Lambda \sqrt{D}$ where the square root is applied to the diagonal elements of $D$? $\endgroup$ – whuber Jan 29 '13 at 4:00
  • $\begingroup$ I thought that too, until I tried it. Unfortunately, it is not true that $L = P^T\Lambda \sqrt(D)$, because $L$ is lower-triangular, and $P^T\Lambda \sqrt(D)$ is not, because it's been permutated. $\endgroup$ – Hugh Perkins Jan 29 '13 at 4:11
  • 2
    $\begingroup$ But it's not necessary for $L$ to be lower triangular for this method to work! All that is required is that $LL'=\Sigma$. $\endgroup$ – whuber Jan 29 '13 at 4:13
  • $\begingroup$ So you knew the answer all along :-). $\endgroup$ – whuber Jan 29 '13 at 4:18
  • $\begingroup$ Well, I guessed at two plausible answers, but I lacked a way of proving/disproving either of them. $\endgroup$ – Hugh Perkins Jan 29 '13 at 4:21
1
$\begingroup$

With credit to whuber, I'm just typing up the exchange in the comments so that this question is marked as having an answer.

We can write $$ P^\prime \Lambda D \Lambda^\prime P=LL^\prime $$

and $ L=P^\prime\Lambda\sqrt{D} $ where we understand $\sqrt{D}$ to be the square root of the diagonal of $D$. The only requirement for this to work is that $LL^\prime=\Sigma,$ so it's not important in this instance that $L=P^\prime\Lambda\sqrt{D}$ may be permuted.

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