Re: Cholesky and non-definite

[John W. Eaton]

On 18-Nov-2004, E. Joshua Rigler <address@hidden> wrote:

| I still wouldn't mind knowing how chol.oct determines whether a matrix
| is positive definite or not.

In Octave's chol.cc, the code is

  int info;
  CHOL fact (m, info);
  if (nargout == 2 || info == 0)
    {
      retval(1) = static_cast<double> (info);
      retval(0) = fact.chol_matrix ();
    }
  else
    error ("chol: matrix not positive definite");

The statement

  CHOL fact (m, info);

creates a "CHOL" object.  If that is successful (info == 0) or you
requested the info output (nargout == 2), then the first return value
is the matrix factor R.

You can find the definition of CHOL in liboctave/dbleCHOL.cc.  This is
the code that initializes the CHOL object, computing the factorization
in the process:

  int
  CHOL::init (const Matrix& a)
  {
    int a_nr = a.rows ();
    int a_nc = a.cols ();

    if (a_nr != a_nc)
      {
        (*current_liboctave_error_handler) ("CHOL requires square matrix");
        return -1;
      }

    int n = a_nc;
    int info;

    chol_mat = a;
    double *h = chol_mat.fortran_vec ();

    F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1),
                               n, h, n, info
                               F77_CHAR_ARG_LEN (1)));

    if (f77_exception_encountered)
      (*current_liboctave_error_handler) ("unrecoverable error in dpotrf");
    else
      {
        // If someone thinks of a more graceful way of doing this (or
        // faster for that matter :-)), please let me know!

        if (n > 1)
          for (int j = 0; j < a_nc; j++)
            for (int i = j+1; i < a_nr; i++)
              chol_mat.elem (i, j) = 0.0;
      }

    return info;
  }

so info is just the value returned from the Lapack subroutine DPOTRF.

If you look at the file libcruft/lapack/dpotrf.f, you'll see

*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          > 0:  if INFO = i, the leading minor of order i is not
*                positive definite, and the factorization could not be
*                completed.

In Octave, INFO < 0 should never happen.

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