[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 <=