Re: chol2inv for Octave?

[Mike Miller]

On Fri, 6 May 2005, John W. Eaton wrote:

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


That's quick service!!  I'm very impressed.

I'll make two recommendations regarding the description of chol2inv:

(1) I would use "U" instead of "R" in the description of chol2inv to 
emphasize that the input matrix must be an upper-triangular matrix.

(2) There is a question about the meaning of "a Cholesky factor." 
Cholesky's theorem states that for any symmetric positive-definite matrix 
M, there exists a unique upper-triangular matrix U with all diagonal 
elements greater than zero such that U'*U = M.  Thus, any upper-triangular 
matrix with positive diagonal elements can be considered a "Cholesky 
factor."  We can make that explicit in the description.

Thus, I recommend that the description for chol2inv read as follows:

-- Loadable Function:  chol2inv (U)
    Invert a symmetric, positive definite square matrix from its
    Cholesky decomposition, U.  Note that U should be an upper-triangular
    matrix with positive diagonal elements.  chol2inv(U) provides
    inv(U'*U) but it is much faster than using inv.  See also: chol,
    cholinv.

The only thing missing now is whether cholinv, chol2inv and chol will work 
on symmetric Hermitian matrices containing complex numbers.  I really only 
work with reals, so I don't know much about complex numbers.

One small point:  If "U" has negative diagonal elements, but none equal to 
zero, then U'*U is still symmetric positive definite, so the chol2inv 
might still run but I don't know what it requires (this must depend on 
LAPACK).  If any diagonal element of U equals 0, then U'*U will be 
singular and chol2inv should fail.

Thanks again!

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