From help-request at octave dot org Fri May 6 04:12:26 2005 Subject: Re: chol2inv for Octave? From: David Bateman To: Help-Octave Date: Fri, 06 May 2005 11:08:13 +0200 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 -------------------------------------------------------------