Incorrect result for inverse in Atlas plugin

[arafatkatze]

Atlas uses clapack_getrf and clapack_getri Lapacke uses lapacke_getrf and lapacke_getri methods to compute the inverse of a matrix.

The results seem to almost equal for non - singular matrices.

Lapacke

[1] pry(main)> require 'nmatrix/lapacke'=> true
[2] pry(main)> b = NMatrix.new([3, 3], [1.0,2.0,1.0,-2.0,-3.0,1.0,3.0,5.0,214], dtype: :float64)
=> 
[
  [ 1.0,  2.0,   1.0]   [-2.0, -3.0,   1.0]   [ 3.0,  5.0, 214.0] ]
[3] pry(main)> b.det
=> 213.99999999999997
[4] pry(main)> b.invert
=> 
[
  [  -3.0233644859813085,   -1.9766355140186918,  0.023364485981308247]
  [    2.014018691588785,     0.985981308411215, -0.014018691588784993]
  [-0.004672897196261678, 0.0046728971962616845,  0.004672897196261682]
]

Atlas

[1] pry(main)> require 'nmatrix/atlas'=> true
[2] pry(main)> b = NMatrix.new([3, 3], [1.0,2.0,1.0,-2.0,-3.0,1.0,3.0,5.0,214], dtype: :float64)
=> 
[
  [ 1.0,  2.0,   1.0]   [-2.0, -3.0,   1.0]   [ 3.0,  5.0, 214.0] ]
[3] pry(main)> b.det=> 213.99999999999997
[4] pry(main)> b.invert
=> 
[
  [   -3.023364485981309,  -1.9766355140186918,   0.02336448598130841]
  [   2.0140186915887854,    0.985981308411215, -0.014018691588785047]
  [-0.004672897196261738, 0.004672897196261627,  0.004672897196261682]
]

But for Singular matrices their is a huge difference.

Lapacke

[1] pry(main)> require 'nmatrix/lapacke'
=> true
[2] pry(main)> b = NMatrix.new([3, 3], [1.0,2.0,1.0,-2.0,-3.0,1.0,3.0,5.0,0], dtype: :float64)
=> 
[
  [ 1.0,  2.0, 1.0]   [-2.0, -3.0, 1.0]   [ 3.0,  5.0, 0.0] ]
[3] pry(main)> b.invert
=> 
[
  [-3.752999689475412e+15,  3.752999689475412e+15, 3.752999689475412e+15]
  [ 2.251799813685247e+15, -2.251799813685247e+15,   -2251799813685246.8]
  [    -750599937895081.6,      750599937895082.6,     750599937895082.2]
]

Atlas

[1] pry(main)> require 'nmatrix/atlas'
=> true
[2] pry(main)> b = NMatrix.new([3, 3], [1.0,2.0,1.0,-2.0,-3.0,1.0,3.0,5.0,0], dtype: :float64)
=> 
[
  [ 1.0,  2.0, 1.0]   [-2.0, -3.0, 1.0]   [ 3.0,  5.0, 0.0] ]
[3] pry(main)> b.invert
=> 
[
  [ 2.0,  0.5,  0.5]   [-3.0,  2.5, -0.2]   [ 5.0, -2.5, -0.0] ]

Note that the although both the results for inverse are actually incorrect, result by lapacke method classify the matrix as ill conditioned matrix but the same is not true for atlas. A more elaborate description can be found in the pr #474 and in issue #470 .

[translunar]

From Wikipedia: "A square matrix that is not invertible is called singular or degenerate."

I would not expect two different implementations to give the same result for an invalid input. I don't think we can guard against this situation without adding a great deal of overhead. I support having functions the user can call, like #nonsingular? to determine if the matrix is invertible, but calling them needs to be optional (because such functions are expensive).

Does this answer your question?

[wlevine]

It's true that the result isn't well defined, but it still looks like something weird is going on in the ATLAS case. It looks like it isn't even trying to compute the inverse, in which case it should be telling us and we should raise an exception.

[translunar]

Oh, okay, misunderstood. I see it now.

[arafatkatze]

@wlevine That's exactly what i wanted to say. In Atlas inverse,after LU decomposition by clapack_getrf.The inverse is not calculated by clapack_getri. This only happens in case of ill conditioned matrices.

[translunar]

@Arafatk Do you get the same behavior when you do it with dtype :object instead? Have you looked at the ATLAS API linkages to see if you can track down where it's failing to give an error message?

[Shailesh-Tripathi]

Hi, Can't this be resolved by just putting a singularity check at the top of the function? And if singular just raise an exception and exit.

[arafatkatze]

@shailesh210 No, there are precision issues with using doubles in inversion and other types of matrix operations.

I recommend that you search through similar issues and read a bit about ill conditioned matrices.

[Shailesh-Tripathi]

Hi @Arafatk Ok I read about the issues. But then what can be done in case of ill-conditioned matrix? Raise an exception or work with them?

[translunar]

@shailesh210 What do we do elsewhere in the code?

[Shailesh-Tripathi]

Hi @mohawkjohn Sir that is what I am asking too :) How to deal with them? What do you suggest.

[translunar]

It's not a rhetorical question. Part of participating in an open source project is doing research and studying — especially in science. Questions you need to answer to address this problem:

1. How does NMatrix handle similar situations elsewhere in the code? (Separately, is that approach the correct one?)
2. How do other similar libraries handle this type of situation?

You don't have to call me sir. We're not very formal. :)

[Shailesh-Tripathi]

Hi @mohawkjohn Sorry for the confusion. I thought I was not clear that's why I asked it again. I will read about ill-contioned problems and discuss with you when I get some solution. :)

[Shailesh-Tripathi]

Hi @mohawkjohn Since I don't have much experience with this, kindly suggest me some sources from where I can learn. Also, what about Nmatrix behaviour with these matrices?

@Arafatk

[translunar]

@shailesh210 Matlab and Numpy are both good places to start.

[arafatkatze]

@shailesh210 Yes, Numpy and Matlab are pretty cool, and if you get confused about what might be the best solution to a problem you can look up in the numpy documentation and even the implementation links are provided in the documentation.
Ultimately there is no "perfect" way to do things and you would have to make a few decisions for yourself.

[Shailesh-Tripathi]

Hi @mohawkjohn and @Arafatk I was going through the issue but I am not able to reproduce the error. Yes I am getting the two outputs different but not the kind stated above.

### require 'nmatrix/atlas'
=> true
[2] pry(main)> b = NMatrix.new([3, 3], [1.0,2.0,1.0,-2.0,-3.0,1.0,3.0,5.0,0], dtype: :float64)
=> 
[
  [ 1.0,  2.0, 1.0]
  [-2.0, -3.0, 1.0]
  [ 3.0,  5.0, 0.0]
]
[3] pry(main)> b.invert
=> 
[
  [  3.602879701896396e+16, -3.602879701896397e+16, -3.602879701896397e+16]

 [-2.1617278211378376e+16,  2.161727821137838e+16,  2.161727821137838e+16]

 [  7.205759403792793e+15, -7.205759403792794e+15, -7.205759403792794e+15]
]

require 'nmatrix/lapacke'
=> false
[7] pry(main)> b = NMatrix.new([3, 3], [1.0,2.0,1.0,-2.0,-3.0,1.0,3.0,5.0,0], dtype: :float64)
=> 
[
  [ 1.0,  2.0, 1.0]
  [-2.0, -3.0, 1.0]
  [ 3.0,  5.0, 0.0]
]
[8] pry(main)> b.invert
=> 
[
  [-6.004799503160659e+15, 6.004799503160659e+15, 6.004799503160659e+15]

 [    3602879701896395.5,   -3602879701896395.5,   -3602879701896395.5]

 [   -1200959900632131.2,    1200959900632132.2,    1200959900632131.8]
]

Due to the slight different implementation of clapack_getrf() and lapacke_getrf() the outputs vary but they seem legible as both are using LU decomposition.

[translunar]

Okay. Having given it some more thought, I'm going to go back on my earlier comment. Inversion of a singular matrix should probably raise an exception. The function should not be giving different outputs with different libraries; it should refuse to give output.

[micw523]

I don't know if this issue has been solved, but from my experience dealing with numpy and MATLAB is that the code generally ignores if the matrix is singular and just returns garbage. LAPACK uses the LU decomposition for inversion and it is (almost) impossible that matrix will be purely singular when decomposed.

I have no experience with the ATLAS library and do not know the behavior of that library for inverting a matrix.

Numpy owners deliberately ignored this issue as said in the following post. They decided that the user should know when inverting a near-singular matrix. https://github.com/numpy/numpy/issues/2074

Deciding if a matrix is singular is not a trivial problem for computation. I suggest that maybe a rank function could be added using SVD so that the user can use it to check the rank of the matrix, if such function is not available.