help-octave
[Top][All Lists]

## Re: factorization of singular matrices - avoiding negative eigenvalues

 From: Brendan Drew Subject: Re: factorization of singular matrices - avoiding negative eigenvalues Date: Fri, 14 Oct 2005 18:02:07 -0600 User-agent: Mozilla Thunderbird 1.0 (X11/20041206)

My suggestion would be to be very careful with the Cholesky factorization. Numerical Recipes in C has the Cholesky factorization code and in there is a test to see if one of the values computed is negative. This test can easily be modified: if something is of sufficiently small magnitude and is negative, print a warning saying that you are going to recondition the matrix and then change that value to an equally small positive value. Negative values of sufficiently large magnitude are flagged as errors as one would hope. This way semi-positive definite matrices remain such, those matrices which due to numerical issues are nearly not semi-positive definite are handled specially and all other cases bail out as one would expect Just my \$currency_of_choice 0.02. I'm sure that the particularly clever amongst us could prove error bounds as a result of applying this technique --- this is just one of many tricks I've been taught to deal with almost-misbehaved covariance matrices. Of course, this probably has lots of nasty implications that I haven't thought about, so take everything that I've said with a suitably large grain (mine?) of salt ;-)

Mike Miller wrote:

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:

http://www.gatsby.ucl.ac.uk/~iam23/code/mvnrnd.m

% 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!

Mike

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

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

--
Brendan Drew
Research Scientist
PercepTek, Inc.
Littleton, CO 80125
Tel: 720-344-1037 x 126
Fax: 720-344-2360

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

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------