Re: chol2inv for Octave?

[David Bateman]

Octave has the "chol" function to calculate the cholesky factorization 
of a symmetric positive definite matrix, however octave can't have use 
of the fact that the factorization is triangular in its solve function 
to just have to do a back substitution. There is a replacement for 
"chol" in octave-forge that defines a new type for triangular matrices, 
based on the matrix type that also defines a "\" operator that performs 
the backsubstitution, but doesn't use BLAS to do it and is so slower 
than it needs to be.

What really needs to be done is have octave recognize triangular, 
hermitian and positive definite matrices in its solve functions and 
treat them differently. This is what I did for the sparse solve 
functions (check the section  "Linear Algebra on Sparse Matrices" of the 
manual), where there are 9 different ways to treat the solve (though the 
cholesky still isn't implemented), and is definitely something I'd like 
to do for the matrix class as well. It is interesting to note that 
matlab itself does something along these lines. Check

http://www.mathworks.com/access/helpdesk/help/techdoc/ref/mldivide.html

to look at what matlab does to choose the algorithm to treat a 
particular matrix...

As for the solution proposed by Michael, no insult meant, but I don't 
think much of it. The problem with Michael's solution is that even if 
the version of chol from octave-forge is used, he doesn't profit from 
its reimplementation of the "\" and calls "Finv" instead. He might as 
well had just called "Finv" directly on the matrix a, and it would have 
been faster. Use something
like

r = chol (a);
x = r \ (r' \ b);

instead. So in short, use the version of the "chol" function from 
octave-forge and you should get speed similar to "R", though not as fast 
as it could be, and in a while all of this will be treated transparently 
within octave itself....

Cheers
David



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