Numerical problems using eigenvectors of large data with mvnpdf

Hello,

I am computing the eigenvectors of some large data with something like this:

XN is a matrix of size m x n containing some normalized data per row, where n is much larger than m. I would not like to use svd directly, but rather try an approach that computes the eigenvalues of just the m dimensions, as below:

[UU, S, ~] = svd((1/nr_groups) * (XN * XN'));

SIQ = pinv(sqrt(S));

U = (1/sqrt(nr_groups)) * (XN' * UU) * SIQ;

W = U(:,1:K); end

Then I compute the covariance of some data: xsigma = cov(data);

and then I project this covariance using the projection matrix W computed above:

psigma = W' * xsigma * W;

Then I want to use this psigma, along with some mean data with mvnpdf. This generally works fine but on some occasion I got this error:

Error using mvnpdf SIGMA must be a square, symmetric, positive definite matrix.

Basically for some value of psigma if I do the following:

mvnpdf(data, miu, psigma)

I get an error.

But if I copy the upper part into the lower part with something like this:

ts = triu(psigma) + triu(psigma, 1)';

and then use it:

mvnpdf(data, miu, ts)

I don't get the error anymore.

It sounds like a numerical problem, but doing the triu trick seems a bit odd to me.

Also, here is the problematic covariance matrix:

   1.0e-03 *

   0.031810217117700   0.002944646130627   0.062672109562341   0.027222257033394
   0.002944646130627   0.006499918688582   0.009385092678967   0.003494166116038
   0.062672109562340   0.009385092678967   0.153523119650464   0.066650626341073
   0.027222257033394   0.003494166116038   0.066650626341073   0.036189655028314

If you try that with mvnpdf you get the error, while if you "correct" the symmetry by the trick mentioned above it seems to work.

Any help is welcome. Thanks.

[Omar Choudary]

Update:

I've tried using the SVD directly and I still have the same problem.

That is, I compute [U, S, V] = SVD(X) where X is the covariance matrix of some data set.

Then I take W = U(:,1:K) as the projection matrix for the first K principal directions.

Next I compute the covariance matrix xsigma for some other data (using cov), and then I project this matrix using W as follows:

psigma = W' * xsigma * W;

Well now for some datasets the projection of xsigma into psigma leads to a matrix that is not perfectly symmetric, as shown in the earlier example. I think that the multiplication above introduces some small numerical errors which cause the matrix not to be perfectly symmetric anymore and that causes an algorithm such as mvnpdf to reject psigma.

Any thoughts?