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?