Warn about Ill-Conditioned Matrices

[gtamba]

I was able to reproduce #436 for every nxn square matrix whose elements are the natural numbers starting from 1 to n^n using this script. The output shows that each of these matrices has an almost zero determinant and returns an inverse when inverted.

These matrices cover a subset of the family of ill-conditioned matrices that are numerically unstable and sensitive to errors when precision is a requirement. (Like how singularity is defined EXACTLY at det = 0 and not the infinitesimally small neighborhood around it)

It is interesting to note that LAPACK does not supply a determinant routine. I googled rather hard to find out why but all I could find was this in the FAQ :-

Are there routines in LAPACK to compute determinants?

No. There are no routines in LAPACK to compute determinants. This is discussed in the "Accuracy and Stability" chapter in the LAPACK Users' Guide.

I decided to test out the same issues on Matlab and as expected an inverse was returned, however, Matlab was aware of it :-

>> a = [1.0, 2.0, 3.0; 4.0, 5.0, 6.0 ; 7.0 , 8.0 , 9.0]

a =

     1     2     3
     4     5     6
     7     8     9

>> inv(a)
Warning: Matrix is close to singular or badly scaled. 
Results may be inaccurate. RCOND =  1.541976e-18. 

ans =

   1.0e+16 *

   -0.4504    0.9007   -0.4504
    0.9007   -1.8014    0.9007
   -0.4504    0.9007   -0.4504

Determinant also fails to be zero.

>> det(a)

ans =

   6.6613e-16

It seems that Matlab uses a cond() and rcond() function to compute how numerically unstable the matrix is ; a higher condition number implies a more unstable matrix.

>> cond(a)

ans =

   3.8131e+16

I think the issue here is not with the decomposition algorithms, simply the hard limitations of machine floating point accuracies. The condition number in this case is incredibly high.

[TL ; DR;-]

I propose that we implement cond() and perhaps call it while calculating inverse to produce a warning. It might be useful as a public function for someone who is aware of this to check their workspace along the way.

[translunar]

Great research on this. Are you going to work on it?

[gtamba]

Yes, I'm still researching the best way to do this because it doesn't seem to be a trivial calculation, I'll send in a PR when I have something concrete.

[wlevine]

This is a good idea, but is there anyway to check programmatically if the warning appears? You wouldn't want to throw an exception because then you couldn't return the result of the calculation.

[wlevine]

It looks like to calculate condition numbers, you need to calculate matrix norms. #400 and #204 are unmerged pull requests to implement some matrix norms, they might be a good place to start.

[gtamba]

I think the best approach for this is using QR factorization mentioned in #440 , so I'm going to focus on exposing that from LAPACKE first since a lot of issues seem to want that as well.

The norms can be used in a more stable way on a QR factorization and the inverse won't have to be calculated every time. My idea is to use QR factorization for a general condition number calculation and the standard method (Norm of matrix times Norm of inverse) for when inverse() is called because we're making a call to inverse anyway, but I will benchmark both once QR becomes available.

[gtamba]

I see that there is a WIP PR already open using the matrix norms to calculate condition number using matrix norms. Ill work on getting QR to work in the meanwhile so that more options are available to benchmark against.

I've already made some progress but I need to tweak a few things to make it work properly. I've also found a condition number routine in LAPACK so I'll look at the feasibility of that as well.

@wlevine : I don't think we can catch a warning in Ruby unless we store it somewhere in a booleanish bit. I'll look into this though.

[translunar]

I believe norms were merged as part of #586, so this should be doable.

[translunar]

I have written some C++ code for computing rcond:

https://github.com/SciRuby/nmatrix/pull/594#issuecomment-300217643

This probably needs to be implemented in Java as well, for nmatrix-jruby (see #594).