help-octave
[Top][All Lists]
Advanced

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

```

reply via email to

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