chol2inv for Octave?

[John W. Eaton]

On  5-May-2005, Mike Miller wrote:

| A friend who uses R a lot was just showing me that to compute the inverse 
| of a symmetric positive-definite matrix, R has the function "chol2inv" 
| which takes the Cholesky upper triangle of the symmetric PD matrix as 
| input and returns the inverse of the symmetric PD matrix as output. 
| Thus, one may get the inverse of a symmetric PD matrix 'M' by doing this:
| 
| chol2inv(chol(M))
| 
| In R, this runs twice as fast as solve(M), which is R's command for 
| computing the inverse.  Comparing R to Octave on the same machine and 
| using a 1000 x 1000 matrix M, I find that R's solve(M) is slower by about 
| 25% than Octave's inv(M), but R's chol2inv(chol(M)) is much faster.  In 
| fact, chol2inv(chol(M) in R is about as fast as inv(chol(M)) in Octave, 
| but R is computing the inverse of M and Octave is computing the inverse of 
| chol(M).
| 
| So I'm wondering why Octave doesn't have a chol2inv function!  ;-)
| Maybe it has one and I just can't find it.

How about the following?  It provides the following functions

 -- Loadable Function:  cholinv (A)
     Use the Cholesky factorization to compute the inverse of of the
     symmetric positive definite matrix A.  See also: chol, chol2inv.


 -- Loadable Function:  chol2inv (R)
     Invert a symmetric, positive definite square matrix from its
     Cholesky decomposition, R.  Note that no check is performed to
     ensure that R is actually a Cholesky factor.  See also: chol,
     cholinv.

and also versions callable from C++.

The actual changes to the code are checked in to the main branch of
the CVS archive.  These changes should also apply to the 2.1.x
sources, but I have no plans to include them in any future 2.1.x
snapshot.

jwe


libcruft/ChangeLog:

2005-05-06  John W. Eaton  <address@hidden>

        * lapack/dpotri.f, lapack/dlauum.f, lapack/dlauu2.f,
        lapack/zpotri.f, lapack/zlauum.f, lapack/zlauu2.f: New files.


liboctave/ChangeLog:

2005-05-06  John W. Eaton  <address@hidden>

        * dbleCHOL.cc (CHOL::init): Use xelem instead of elem for indexing
        chol_mat.
        (chol2mat_internal, chol2mat, CHOL::inverse): New functions.
        * dbleCHOL.h (chol2mat_internal, chol2mat, CHOL::inverse):
        Provide decls.

        * CmplxChol.cc (ComplexCHOL::init): Use xelem instead of elem for
        indexing chol_mat.
        (chol2mat_internal, chol2mat, ComplexCHOL::inverse): New functions.
        * CmplxCHOL.h (chol2mat_internal, chol2mat, CmplxCHOL::inverse):
        Provide decls.

src/ChangeLog:

2005-05-06  John W. Eaton  <address@hidden>

        * DLD-FUNCTIONS/chol.cc (Fcholinv): New function.
        (Fchol2inv): New function.



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