[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:
`
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
-------------------------------------------------------------

**factorization of singular matrices - avoiding negative eigenvalues**,
*Mike Miller* **<=**