Strange behaviour of \ w.r.t. chol

[David Bateman]

Marco Caliari wrote:
> Dear David,
>
> matrix B is not singular, as you can see from its eigenvalues (and
> from Matlab). Moreover, the Cholesky factorization can be successfully
> computed, as you can see running chol(B). So, if the test, even if
> weak, detects B to be 'Positive Definite' and if chol(B) can be
> computed, why det(B)=0 and B\y is not solved using Cholesky?
>
> Best regards,
>

Ok, I didn't read your message completely. My remark still holds, just
not in your case.. In any case here is a patch for the matrix_type
documentation issue.

As for the issue you are see, the difference is the manner in which
Octave is handling near singular matrices relative to Matlab. The
typical lapack function returns two values that are useful to test for
singular matrix.. The "info" which is negative for an illegal input,
zero for a successful factorization and positive for a failed
factorization.  The second is the reciprocal condition number "rcond".
The dpotrf function does not return info > 0 values and we must rely on
the value of rcond to test singular. How Octave treats these values in
general is something like the following

          if (info != 0)
            info = -2;

          volatile double rcond_plus_one = rcond + 1.0;

          if (rcond_plus_one == 1.0 || xisnan (rcond))
            {
              info = -2;

              if (sing_handler)
            sing_handler (rcond);
              else
            (*current_liboctave_error_handler)
              ("matrix singular to machine precision, rcond = %g",
               rcond);
            }

Therefore, if rcond is less than eps, then the matrix is flagged as
singular and eventually xGELSD is called to factor it. How it appears
that matlab is treating this is more like

if (info != 0 || xisnan(rcond))
  {
    info = -2;

    if (sing_handler)
      sing_handler (rcond);
    else
      (*current_liboctave_error_handler)
    ("matrix singular to machine precision, rcond = %g",
     rcond);
  }
else
  {
    volatile double rcond_plus_one = rcond + 1.0;

    if (rcond_plus_one == 1.0)
      {
    volatile double rcond10_plus_one = rcond * 10.0 + 1.0;
    if (rcond10_plus_one == 1.0)
      {
        (*current_liboctave_warning_handler)
          ("matrix is near singular and the result may be inaccurate,
rcond = %g",
           rcond);
      }
    else
      {
        info = -2;

        if (sing_handler)
          sing_handler (rcond);
        else
          (*current_liboctave_error_handler)
        ("matrix singular to machine precision, rcond = %g",
         rcond);
      }
      }
  }

where for values of rcond less than eps but greater than some value the
factorization is accepted, but a warning issued. I have no idea what
value of rcond a warning is issued rather than falling back to the
xGELSD factorization. This seems to be a general difference for all near
singular matrices. For example


octave:18> a = eye(4,4)*5e16;a(4,4)=1; a(4,1) = 1; a(1,4) = 4;
octave:19> 1/cond(a)
ans =  2.0000e-17
octave:20> b = a \ ones(4,1)
warning: matrix singular to machine precision, rcond = 2e-17
warning: attempting to find minimum norm solution
warning: dgelsd: rank deficient 4x4 matrix, rank = 3
b =

   2.0000e-17
   2.0000e-17
   2.0000e-17
   3.0815e-33

whereas with matlabR2007b I get

>> a = eye(4,4)*5e16;a(4,4)=1; a(4,1) = 1; a(1,4) = 4;
>>  1/cond(a)

ans =

   2.0000e-17

>>  b = a \ ones(4,1)
Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 2.000000e-17.

b =

   -0.0000
    0.0000
    0.0000
    1.0000

as expected. Pushing this case to the extreme in matlab gives

>> a = eye(4,4)*5e256;a(4,4)=1; a(4,1) = 1; a(1,4) = 4;
>>  b = a \ ones(4,1)                                 
Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 2.000000e-257.

b =

   -0.0000
    0.0000
    0.0000
    1.0000

so even in that case the return of dpotrf is used. In fact it appears
that matlab only uses its equivalent of xGELSD for non-square matrices
and just accepts the output of the LU factorization.. Should we change
the behavior of Octave to match this? Historically Octave had an LU
factorization only and fell back to using xGELSS if the matrix was
singular. Check the documentation

www.mathworks.com/access/helpdesk/help/techdoc/ref/mldivide.htm

to see that in fact the fallback to treatment of the matrix with a QR
decomposition doesn't happen in matlab..

D.

-- 
David Bateman                                David.Bateman at motorola.com
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary