From bug-octave-request at bevo dot che dot wisc dot edu Mon Nov  4 14:52:19 2002
Subject: Re: chol problem
From: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
To: Rolf Fabian <fabian at TU-Cottbus dot De>,    "'bug-octave UWISC'" <bug-octave@bevo.che.wisc.edu>
Date: Mon, 4 Nov 2002 14:52:17 -0600

On  3-Nov-2002, Paul Kienzle <pkienzle at jazz dot ncnr dot nist dot gov> wrote:

| On Sun, Nov 03, 2002 at 01:34:33PM +0100, Rolf Fabian wrote:
| > There seems to be a bug in the Cholesky decomposition
| > chol(A) for input of complex diagonal matrices A.
| > 
| > In my opinion, a positive definite DIAGONAL matrix
| > must always be real (with strictly positive elements).
| > As a consequence, if diagonal A has any nonvanishing
| > imag parts, 'chol' should always return an error.
| > 
| > (1) This call correctly invokes an error:
| > :> C = chol( A=diag( [-1+i,1+i] ) )
| >    |- error: chol: matrix not positive definite
| > 
| > (2) But a bug occurs e.g. for:
| > :> C = chol( A=diag( [+1+i,1+i] ) )	# No error!
| > C =
| >   1  0
| >   0  1
| > 
| > :> RES = C'*C - A	# residuals, approx. ZERO expected
| > 			# if C is a Cholesky factor
| > RES =
| >   0 - 1i  0 + 0i
| >   0 + 0i  0 - 1i
| > 
| > If ALL real parts of A are (strictly) positive, Octave's
| > 'chol' doesn't detect, that A is actually NOT positive
| > definite. Hence the resulting C is garbage as can be
| > easily seen in our example by the nonvanishing RES matrix.
| > 
| > I use V2.1.14, ported to OS/2 by K. Gebhardt. This is a
| > bit exotic. Probably, the problem does not occur on other
| > platforms, or has been corrected in later versions.

2.1.14 is pretty ancient (3-1/2 years old now) and no one seems to be
maintaining the OS/2 port now (I don't think I've heard from Klaus in
several years).

| This behaviour still exists in octave 2.1.38 compiled for
| RedHat 7.1 with gcc-2.96 and lapack-3.0-9.rpm.

I think this is a bug in Lapack, because it seems to me that it should
diagnose invalid input.  The documentation for Lapack that I have says
that ZPOTRF performs the factorization for positive definite matrices.
It doesn't say precisely what it will do if the input matrix is not
positive definite.  However, it does say that it will return INFO > 0
if the leading minor of order i is not positive definite and the
factorization could not be completed.

Should Octave try to decide whether the matrix is positive definite
before calling ZPOTRF?  What's the most efficient and reliable way to
do that?

jwe



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