[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

factorization of singular matrices - avoiding negative eigenvalues

From: Mike Miller
Subject: factorization of singular matrices - avoiding negative eigenvalues
Date: Fri, 14 Oct 2005 18:11:20 -0500 (CDT)

We typically create a p-variate multivariate normal vector by Cholesky factorizing the var-covar matrix and multiplying the resultant by an iid normal(0,1) vector:

mvnorm = chol(cov)'*rand(p,1);

or for N mv-normal row vectors:

mvnorm = rand(N,p)*chol(cov);

That works great unless the variance covariance matrix is singular. For that case, I have the answer given in comments at the end of the page at this URL:

% But the Cholesky decomposition function doesn't always work...
% Consider Sigma=[1 1;1 1]. Now inv(Sigma) doesn't actually exist, but Matlab's
% mvnrnd provides samples with this covariance st x(1)~N(0,1) x(2)=x(1). The
% fast way to deal with this would do something similar to chol but be clever
% when the rows aren't linearly independent. However, I can't be bothered, so
% another way of doing the decomposition is by diagonalising Sigma (which is
% slower but works).
% if
%       [E,Lambda]=eig(Sigma)
% then
%       Sigma = E*Lambda*E'
% so
%       U = sqrt(Lambda)*E'
% If any Lambdas are negative then Sigma just isn't even positive semi-definite
% so we can give up.

That seems like the perfect solution until we realize that numerical imprecision causes the lambdas to be small negative numbers when they really should be zeros:

[e,lambda]=eig(ones(3)); sqrt(lambda(1))
ans = 0.00000000000000e+00 + 3.25207016166203e-09i

Thus, checking for negative lambdas will sometimes fail us but not checking for them will cause some complex-valued mv-normal data, which is very undesirable. What do you guys think is the best way of dealing with this problem? I would like to return an error if the matrix has a clearly negative eigenvalue, but not if the eigenvalue was an imperfectly estimated zero. When the zero eigenvalue is estimated as a small negative number, I could either use sqrt(abs(lambda)) instead of sqrt(lambda) or maybe real(sqrt(lambda)). Or maybe I could replace the near-zero values with zeros. I'm not concerned about speed for this step.

Any suggestions?

I would also love to know if any of you can suggest a way to produce a generalized Cholesky factorization for non-negative definite matrices that work when the matrix is singular. For singular non-ND matrices, I think there must be infinitely many possible upper triangular matrices, U, such that U'U returns the original non-ND matrix. I'd be happy with any one of those solutions so long as it contains only real values!


Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:
How to fund new projects:
Subscription information:

reply via email to

[Prev in Thread] Current Thread [Next in Thread]