Re: Inverse Matrix Function appears a bit wonky

[John W. Eaton]

On  5-Nov-2004, David Bateman <address@hidden> wrote:

| There is indeed a significant speed difference, which is why I submitted
| a patch to add the calc_cond option to the inverse functions. I noticed
| a significant speed penalty on this point in the sciview benchmarks.
| 
| What is always calculated by LAPACK is rcond, which is an estimate of
| the reciprocal of the condition number. If you request the condition
| number then "{z,d}gecon" is called that calculates not just an
| estimate, but the real thing.

I'm not sure what has happened, but the code in inv.cc never calls
Matrix::inverse with four arguments, so the calc_cond argument
always has the default value of 1 (which says to compute the condition
number).

If the condition number is not calculated, there seems to be no way to
detect at least some numerically singular matrices.  It seems that
DGETRF can return INFO == 0 even for exactly singular matrices like
[1,2,3;4,5,6;7,8,9].  For example, try the following program:

      program foo
      integer nr, nc
      parameter (nr = 3, nc = 3, lwork = 2*nc)
      double precision a(nr,nc), work(lwork)
      integer ipvt(nr), info
      data a /1.0d0, 4.0d0, 7.0d0,
     $        2.0d0, 5.0d0, 8.0d0,
     $        3.0d0, 6.0d0, 9.0d0/
      do i = 1, nr
        print *, (a(i,j), j = 1, nc)
      enddo
      call dgetrf (nc, nc, a, nr, ipvt, info)
      print *, 'dgetrf info: ', info
      call dgetri (nc, a, nr, ipvt, work, lwork, info)
      print *, 'dgetri info: ', info
      end

On my system, this prints:

    1.  2.  3.
    4.  5.  6.
    7.  8.  9.
   dgetrf info:  0
   dgetri info:  0

The LAPACK docs for DGETRF and DGETRI say:

  *  INFO    (output) INTEGER
  *          = 0:  successful exit
  *          < 0:  if INFO = -i, the i-th argument had an illegal value
  *          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
  *                has been completed, but the factor U is exactly
  *                singular, and division by zero will occur if it is used
  *                to solve a system of equations.
  *

  *  INFO    (output) INTEGER
  *          = 0:  successful exit
  *          < 0:  if INFO = -i, the i-th argument had an illegal value
  *          > 0:  if INFO = i, U(i,i) is exactly zero; the matrix is
  *                singular and its inverse could not be computed.

So if we don't compute the condition number, we have no way of warning
the user that the matrix that they just inverted is numerically
singular.  That seems bad to me (and is why this whole thread started
in the first place).

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