[Math] Trying to use Cholesky decomposition of covariance matrix to sample error ellipsoid

eigenvalues-eigenvectorseuclidean-geometrylinear algebramatricessampling

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 (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.

Thanks in advance!

Best Answer

Your prescription sounds right to me.

Suppose $R=LL'$ is the desired covariance matrix with its Cholesky decomposition. Let e be a random (column) vector distributed as a zero-mean multivariate Gaussian with unit variance and no correlations. That is, $E[e]=0$ and $E[ee'] = I$. Then $L$ and $e$ are the results of your steps 1 and 2.

The matrix-vector product $x=Le$ is the result of your step (3). I'll assert that $x$ is also a multivariate Gaussian, since proof of that is not the issue here. We can see that the expected covariance of $x$ is

$$\text{Covar}(x) = E[xx']-E[x]E[x'] = E[ L e e' L'] - 0 = L E[ee'] L' = LIL' = LL' = R.$$

So $x$ should follow a multivariate Gaussian distribution with the correct covariance matrix, $R$. All I can think of is that if you want to error ellipse to be distributed about a point other than the origin, then (step (4)) you need to add that vector to $x$.

Oh, I know one potential problem: Perhaps you used the upper, rather than the lower Cholesky factor? If you let $x=L'e$, then your random vectors would have an expected covariance of $L'L$, which does not equal $R$ unless $R$ is diagonal (i.e., this only works if there's no correlation, which is surely not the case you have in mind).