Inverse matrix output for zero determinant

[shahsaurabh0605]

invert function gives inverse dense matrix output even if it has zero determinant.

[gtamba]

Could you share the exact matrix that you're trying to invert? When I try it with a singular matrix it gives me what I think should be the correct error :

[14] pry(main)> n = NMatrix.new(3,[1,2,3,4,5,6,2,4,6])
=> 
[
  [1, 2, 3]   [4, 5, 6]   [2, 4, 6] ]
[15] pry(main)> k = n.inverse
ZeroDivisionError: Expected Non-Singular Matrix.
from ~/.rvm/gems/ruby-2.2.2/gems/nmatrix-0.2.0/lib/nmatrix/math.rb:83:in `__inverse__'

I get a similar error for the invert function as well.

[shahsaurabh0605]

I tried the above matrix and it gives me same ZeroDivisionError but try this matrix: n= NMatrix.new(3,[1,2,3,4,5,6,7,8,9]) It has determinant zero but still it gives an inverse matrix!

[translunar]

Most likely this is because NMatrix#det is not equal to NMatrix#det_exact. In other words, the approximation being used to compute the determinant is returning a nonzero value, and so it's unaware the function is nonsingular. I'm curious: if you plug in this matrix and try to compute the inverse using dgetri, does it give the same answer?

[gtamba]

I tried out shahsaurabh0605's array input and the inverse function does indeed misbehave however the determinants seem to be correct.

n = NMatrix.new(3, [1,2,3,4,5,6,7,8,9])
=> 
[
  [1, 2, 3]   [4, 5, 6]   [7, 8, 9] ]
pry> n.inverse
=> 
[
  [  643371375338640.9, -1286742750677283.8, 
  643371375338642.2]   [-1286742750677283.8, 
 2573485501354568.5, -1286742750677284.5]  
 [  643371375338642.6, -1286742750677284.5,   643371375338642.2] ]

So inverse misbehaves with this matrix

 pry> n.det
=> 0
 pry> n.det_exact
=> 0

Could it be that the loss of precision in decimal during the factorization/elimination process prevents it from being unable to compute the inverse of the matrix?

[translunar]

That's not what I'm getting:

[1] pry(main)> NMatrix.new(3,[1.0,2,3,4,5,6,7,8,9]).det_exact
=> 0.0
[2] pry(main)> NMatrix.new(3,[1.0,2,3,4,5,6,7,8,9]).det
=> 6.661338147750939e-16

You probably need to make sure you're using dtype: :float64. By passing 1 as the first element in the array instead of 1.0, it may be treating it as an int.

[gtamba]

I see, that makes sense. Apologies for missing that.

Just a few doubts I had :-

1. Is the inverse calculated using Adjoint and Co-Factors or is it done by elimination?

2. Is it too expensive in the general case to just compute NMatrix#det for the check?

[shahsaurabh0605]

I see that the matrix inverse is calculated using Gauss-Jordan elimination method. So is it using NMatrix#det or NMatrix#det_exact to check matrix is non-singular? Because if not, then there must be some other floating point error.

[gtamba]

This line checks for singularity while inverting :

        if (matrix[k * (M + 1)] == (DType)(0)) {
          rb_raise(rb_eZeroDivError, "Expected Non-Singular Matrix.");
        }

So the determinant does not come into picture here. What I think happens is if the diagonal element on that row cannot be established as a pivot, then a singular matrix is detected. (This happens because of more subtle reasons I guess, but if you can't multiply a subsequent diagonal element to eliminate elements lower in that column, then the whole Gauss-Jordan method is upset , mainly because such matrices are singular and an inverse doesn't exist anyway.)

Still unsure about the exact cause, but it could be a floating point error as shahsaurabh0605 already said.

2 simples fixes I can think of right now might be :

a) Check the determinant before computing inverse. (it's as costly as inversion though, O(n^3) if I'm correct)

b) Verify that the output matrix multiplied by the input matrix gives us an identity matrix, else raise an exception. This may seem like a wasteful option considering that we'd have already inverted a matrix but consider that matrix multiplication is faster than calculating determinants and you're more likely to get a non-singular matrix than a non-invertible one as input.

[gtamba]

I've written a small fix that uses NMatrix#det_exact to reject singular matrices of size 3 or smaller. It solves the issue for this particular matrix and seems like the best solution (<= 3x3 Matrix Determination calculation isn't very costly) until more singular matrices that return inverses show up.

Would it be okay to send a pull request?

[translunar]

I'm actually not sure about this one. This is a particularly common problem with matrix decompositions. Usually it's solved by adding to the diagonal until the matrix is no longer singular. I'm not convinced it makes sense to address one edge case when it's not really possible to address all of them.

[gtamba]

I'll try to find more matrices that cause this problem in the meanwhile then. I suggested the above fix because in det_exact's C++ implementation, the computations are hardcoded multiplications for square matrices of shape 3 or less, (det_exact does not work on larger matrices because of the load) therefore it would run in about constant time.

This particular array returns an inverse in the lapacke and atlas plugins as well even though lapacke uses LU decomposition to calculate the inverse.

[translunar]

I'd much rather we have an exposed #inverse_exact method which relies on #det_exact. How about that?

[shahsaurabh0605]

Yes, I sent mail indicating the same. Problem is #inverse_exact method is not correctly implemented in the gem. If assigned, I would like to execute it correctly.

[gtamba]

#inverse_exact is implemented in ext/nmatrix/math.cpp already and only works for matrices of size 3 or 2. (Ref Version: 0.1.0.rc2 / 2014-03-12)

(uppercase comments denote code blocks edited out)

/*
     * Calculate the exact inverse for a dense matrix (A [elements]) of size 2 or 3. Places the result in B_elements.
     */
    template <typename DType>
    void inverse_exact(const int M, const void* A_elements, const int lda, void* B_elements, const int ldb) {

       //INITIALIZE
      //SOLVE FOR 1X1
     // SOLVE FOR 2X2

      else if (M == 3) {
        // Calculate the exact determinant.
        DType det;
        det_exact<DType>(M, A_elements, lda, reinterpret_cast<void*>(&det));
        if (det == 0) {
          rb_raise(nm_eNotInvertibleError, 
              "matrix must have non-zero determinant to be invertible (not getting this error does not mean matrix is invertible if you're dealing with floating points)");
        }

       //CALCULATE INVERSE FOR 3X3 MATRIX

      } else if (M == 1) {
        B[0] = 1 / A[0];
      } else {
        rb_raise(rb_eNotImpError, "exact inverse calculation needed for matrices larger than 3x3");
      }
    }

I'm not sure if this function is already exposed or not but it doesn't work for matrices larger than 3 since we don't have a fast algorithm for calculating the exact determinant of large matrices yet.

Is the proposal to implement _exact methods for a general square matrix? Sorry for the silly questions and thanks!

[v0dro]

@mohawkjohn why even keep this function if #inverse can very well calculate inverses for 3x3 matrices?

[gtamba]

Here's another singular square matrix that returns an inverse :-

[22] pry(main)> NMatrix.new([3,3], [1.0,2.0,1.0,-2.0,-3.0,1.0,3.0,5.0,0]).det
=> 1.3322676295501884e-15

[23] pry(main)> NMatrix.new([3,3], [1.0,2.0,1.0,-2.0,-3.0,1.0,3.0,5.0,0]).det_exact
=> 0.0

[24] pry(main)> NMatrix.new([3,3], [1.0,2.0,1.0,-2.0,-3.0,1.0,3.0,5.0,0]).inverse
=> 
[
  [-3.752999689475412e+15,  3.752999689475412e+15, 3.752999689475412e+15]
  [ 2.251799813685247e+15, -2.251799813685247e+15,   -2251799813685246.8]
  [    -750599937895081.6,      750599937895082.6,     750599937895082.2]
]

I had 'nmatrix/lapacke' plugged in for this one.

[lprimeroo]

Any progress on this issue? As far as I can derive from the issue thread above, the problem is caused by matrices of sizes < 3. Can't we just hard-convert the elements to dtype :float64 for n < 3? Or deprecate inverse_exact altogether?

[translunar]

@v0dro We keep it because the computation of exact inverse for 3x3 matrices (e.g., position, velocity) is a common one, and you want it to be exact. #inverse is an approximation.

@saru95 This comment gives the fix: https://github.com/SciRuby/nmatrix/issues/436#issuecomment-175048806