factorization of singular matrices - avoiding negative eigenvalues

[Mike Miller]

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