# Trying to use Cholesky decomposition of covariance matrix to sample error ellipsoid

I'm trying to construct an error ellipsoid from a covariance matrix (which exists for a 3D point) and then sample consistent xyz points in this region. In a previous question when I asked about this (on math.stackexchange, here) a user suggested that I compute the Cholesky decomposition of the matrix and then:

sample Gaussian with identity matrix as the covariance (easy to do by sampling each coordinate as a Gaussian) and then apply the linear transformation from the Cholesky decomposition matrix to get your sampled point from the ellipsoid covariance Gaussian

I've had a go at following this but I assume I've misunderstood to at least some extent as the result isn't consistent with what I'd expect. The steps I took were:

1) Calculate the Cholesky decomposition of the covariance matrix.
2) Sample each initial vertex point as a Gaussian with width 1 to generate (x', y', z')
3) Multiply (x',y',z') by the Cholesky decomposition matrix for the newly generated point.

I'm certain this isn't correct, but don't have the experience to know exactly what is wrong.

To add some more details, the covariance matrix of the initial point is:

$\begin{pmatrix} 10.0115 & -10.6835 & 5.18024 \\ -10.6835 & 11.4009 & -5.52798 \\ 5.18024 & -5.52798 & 2.77646 \\ \end{pmatrix}$

I then calculate the Cholesky decomposition, which is:

$\begin{pmatrix} 3.164095& 0& 0\\ -3.376478& 0.017131& 0\\ 1.637195& -0.001619& 0.309921\\ \end{pmatrix}$

And the initial xyz point is at (35.5361, -37.2661, 22.521).

As a check before sampling the Gaussian that this multiplication would yield something sensible I multiply:

$\begin{pmatrix} 35.5361& 0& 0\\ 0& -37.2661& 0\\ 0& 0& 22.521\\ \end{pmatrix}$

by the Cholesky decomposition matrix (assuming that the Gaussian sampling would look something like this), and then I was expecting this to yield some new xyz point consistent with the original. The actual result is:

$\begin{pmatrix} 112.44& 0& 0\\ 0& -0.6384065& 0\\ 0& 0& 6.97973\\ \end{pmatrix}$

• You are missing 4) add the desired mean of the distribution (provided that in part 2 'width 1' = standard deviation 1 and mean is 0). If this does not fix the issue, can you add the results (e.g., some plot) and information about why they are not consistent with what you expect? The method you describe, if correctly implemented, samples from the multivariate Gaussian, but maybe you want something else than samples from the distribution. – Juho Kokkala Jun 9 '14 at 18:53
• Many thanks for your comment - I expanded on the original post to point out more precisely what I was doing, so perhaps you can see more clearly which step is incorrect (maybe my implementation isn't consistent with how I explained things). – anthr Jun 9 '14 at 19:10