Is there some method to compute the inverse of a matrix and obtain information if the matrix is invertible, like ei_compute_inverse_in_size2_case_with_check ?
As most of the time computing if the matrix is invertible is as costly as inverting the matrix, doing the check and then inverting the matrix seems not optimal.

If the method is not available:
How complex it is to have this feature ?

I propose the function:
bool compute_inverse_with_check( const MatrixType& matrix, MatrixType* result )

and the associated method of matrix_base:
bool computeInverseWithCheck( MatrixType* result ) const

I have modified the file Inverse.hpp and MatrixBase.hpp (just added the declaration of one method there) to include those features, it compiles and seems to execute fine on some random examples.

You can find the files there:
http://rapidshare.com/files/246412035/m ... s.tgz.html

What do you think of that feature ? Of the way I modified your code ?

I am not very happy with the fact that there is a lot of copy/paste from the original code of compute_inverse.
But the only way I see, for the moment, to factorize the code is to use macros (do you see a better choice ?)
for instance:

template
bool ei_compute_inverse_in_size3_case_with_check(const XprType& matrix, MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
const Scalar det_minor00 = matrix.minor(0,0).determinant();
const Scalar det_minor10 = matrix.minor(1,0).determinant();
const Scalar det_minor20 = matrix.minor(2,0).determinant();
const Scalar det = ( det_minor00 * matrix.coeff(0,0)
- det_minor10 * matrix.coeff(1,0)
+ det_minor20 * matrix.coeff(2,0) );
if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false;
const Scalar invdet = Scalar(1) / det;
result->coeffRef(0, 0) = det_minor00 * invdet;
result->coeffRef(0, 1) = -det_minor10 * invdet;
result->coeffRef(0, 2) = det_minor20 * invdet;
result->coeffRef(1, 0) = -matrix.minor(0,1).determinant() * invdet;
result->coeffRef(1, 1) = matrix.minor(1,1).determinant() * invdet;
result->coeffRef(1, 2) = -matrix.minor(2,1).determinant() * invdet;
result->coeffRef(2, 0) = matrix.minor(0,2).determinant() * invdet;
result->coeffRef(2, 1) = -matrix.minor(1,2).determinant() * invdet;
result->coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet;
return true;
}

can be replaced by:

#define MINOR_AND_DETERMINANT_COMPUTATION_3D
typedef typename MatrixType::Scalar Scalar;
const Scalar det_minor00 = matrix.minor(0,0).determinant();
const Scalar det_minor10 = matrix.minor(1,0).determinant();
const Scalar det_minor20 = matrix.minor(2,0).determinant();
const Scalar det = ( det_minor00 * matrix.coeff(0,0)
- det_minor10 * matrix.coeff(1,0)
+ det_minor20 * matrix.coeff(2,0) );

#define INVERSE_FROM_COMATRIX_3D
...

template
bool ei_compute_inverse_in_size3_case_with_check(const XprType& matrix, MatrixType* result)
{
MINOR_AND_DETERMINANT_COMPUTATION_3D
if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false;
INVERSE_FROM_COMATRIX_3D
return true;
}

#undef MINOR_AND_DETERMINANT_3D
#undef INVERSE_FROM_COMATRIX_3D

and use this macro where it is usefull in ei_compute_inverse_in_size3_case for instance.
However, I don't think that way of coding is very nice especially for two instances of the code only.

I just had a look to main page and have seen that the current versioning tool is mercurial.
I did a checkout with svn, so perhaps the files I modified were not up-to-date.
Is the svn repository updated too ?

I will get the mercurial repository and see it tomorrow.
Sorry for the post.

Of course Eigen allows you to check invertibility and to combine this computation with the computation of the inverse so that there is no redundancy. The Eigen way to do that is:

Code: Select all

  MatrixXd m;

  ...

  LU lu(m);
  if(lu.isInvertible())
  {
     the matrix is invertible, its inverse is given by lu.inverse();
  }
  else
  {
     the matrix is singular
  }


This method is fast for large enough matrices. But it's true that it's inherently slow for very small matrices: up to size 4x4.

Are you specifically interested in <= 4x4 matrices and do you specifically need inverse-with-check for them? It's true that we overlooked this case, but that's because we didn't believe that there would be much of a use case for that. What's your use case?

It sure would be possible to add a inverseWithCheck() method, that uses the LU method in the generic case, and specialize to special code in small size.

Don't worry, we won't have to use macros, templates and inline functions allow to do a lot of things

myguel wrote:I just had a look to main page and have seen that the current versioning tool is mercurial.
I did a checkout with svn, so perhaps the files I modified were not up-to-date.
Is the svn repository updated too ?

The SVN repo isn't updated anymore, it will be removed soon.

OK I looked at your files but it would be 10x better for me to have a diff:

if you used svn, you can do:

cd eigen2/
svn diff > my_diff

it creates a file my_diff that you can send us or paste here as code.

Ok, I need the 3x3 case.
Here is the diff:

Code: Select all

diff -r 2dcba9122a77 Eigen/src/Core/MatrixBase.h
--- a/Eigen/src/Core/MatrixBase.h   Fri Jun 19 20:46:55 2009 +0200
+++ b/Eigen/src/Core/MatrixBase.h   Sun Jun 21 23:12:49 2009 +0200
@@ -643,6 +643,7 @@
     const PartialLU partialLu() const;
     const PlainMatrixType inverse() const;
     void computeInverse(PlainMatrixType *result) const;
+   bool computeInverseWithCheck( PlainMatrixType *result ) const;
     Scalar determinant() const;
 
 /////////// Cholesky module ///////////
diff -r 2dcba9122a77 Eigen/src/LU/Inverse.h
--- a/Eigen/src/LU/Inverse.h   Fri Jun 19 20:46:55 2009 +0200
+++ b/Eigen/src/LU/Inverse.h   Sun Jun 21 23:12:49 2009 +0200
@@ -73,6 +73,30 @@
   result->coeffRef(2, 0) = matrix.minor(0,2).determinant() * invdet;
   result->coeffRef(2, 1) = -matrix.minor(1,2).determinant() * invdet;
   result->coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet;
+}
+
+template
+bool ei_compute_inverse_in_size3_case_with_check(const XprType& matrix, MatrixType* result)
+{
+  typedef typename MatrixType::Scalar Scalar;
+  const Scalar det_minor00 = matrix.minor(0,0).determinant();
+  const Scalar det_minor10 = matrix.minor(1,0).determinant();
+  const Scalar det_minor20 = matrix.minor(2,0).determinant();
+  const Scalar det = ( det_minor00 * matrix.coeff(0,0)
+        - det_minor10 * matrix.coeff(1,0)
+        + det_minor20 * matrix.coeff(2,0) );
+  if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false;
+  const Scalar invdet = Scalar(1) / det;
+  result->coeffRef(0, 0) = det_minor00 * invdet;
+  result->coeffRef(0, 1) = -det_minor10 * invdet;
+  result->coeffRef(0, 2) = det_minor20 * invdet;
+  result->coeffRef(1, 0) = -matrix.minor(0,1).determinant() * invdet;
+  result->coeffRef(1, 1) = matrix.minor(1,1).determinant() * invdet;
+  result->coeffRef(1, 2) = -matrix.minor(2,1).determinant() * invdet;
+  result->coeffRef(2, 0) = matrix.minor(0,2).determinant() * invdet;
+  result->coeffRef(2, 1) = -matrix.minor(1,2).determinant() * invdet;
+  result->coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet;
+  return true;
 }
 
 template
@@ -160,6 +184,59 @@
     }
   }
 }
+
+
+template
+bool ei_compute_inverse_in_size4_case_with_check(const XprType& matrix, MatrixType* result)
+{
+  if(ei_compute_inverse_in_size4_case_helper(matrix, result))
+  {
+    // good ! The topleft 2x2 block was invertible, so the 2x2 blocks approach is successful.
+    return true;
+  }
+  else
+  {
+    // rare case: the topleft 2x2 block is not invertible (but the matrix itself is assumed to be).
+    // since this is a rare case, we don't need to optimize it. We just want to handle it with little
+    // additional code.
+    MatrixType m(matrix);
+    m.row(0).swap(m.row(2));
+    m.row(1).swap(m.row(3));
+    if(ei_compute_inverse_in_size4_case_helper(m, result))
+    {
+      // good, the topleft 2x2 block of m is invertible. Since m is different from matrix in that some
+      // rows were permuted, the actual inverse of matrix is derived from the inverse of m by permuting
+      // the corresponding columns.
+      result->col(0).swap(result->col(2));
+      result->col(1).swap(result->col(3));
+     return true;
+    }
+    else
+    {
+      // last possible case. Since matrix is assumed to be invertible, this last case has to work.
+      // first, undo the swaps previously made
+      m.row(0).swap(m.row(2));
+      m.row(1).swap(m.row(3));
+      // swap row 0 with the the row among 0 and 1 that has the biggest 2 first coeffs
+      int swap0with = ei_abs(m.coeff(0,0))+ei_abs(m.coeff(0,1))>ei_abs(m.coeff(1,0))+ei_abs(m.coeff(1,1)) ? 0 : 1;
+      m.row(0).swap(m.row(swap0with));
+      // swap row 1 with the the row among 2 and 3 that has the biggest 2 first coeffs
+      int swap1with = ei_abs(m.coeff(2,0))+ei_abs(m.coeff(2,1))>ei_abs(m.coeff(3,0))+ei_abs(m.coeff(3,1)) ? 2 : 3;
+      m.row(1).swap(m.row(swap1with));
+     if( ei_compute_inverse_in_size4_case_helper(m, result) )
+     {
+        result->col(1).swap(result->col(swap1with));
+        result->col(0).swap(result->col(swap0with));
+        return true;
+     }
+     else{
+        return false; }
+    }
+  }
+
+}
+
+
 
 /***********************************************
 *** Part 2 : selector and MatrixBase methods ***
@@ -254,4 +331,84 @@
   return result;
 }
 
+
+/*****************************************
+ * Compute inverse with invertibility check
+*********************************************/
+
+template
+struct ei_compute_inverse_with_check
+{
+  static inline bool run(const MatrixType& matrix, MatrixType* result)
+  {
+     typedef typename MatrixType::Scalar Scalar;
+     LU lu( matrix );
+     if( !lu.isInvertible() ) return false;
+     lu.computeInverse(result);
+     return true;
+  }
+};
+
+template
+struct ei_compute_inverse_with_check
+{
+  static inline bool run(const MatrixType& matrix, MatrixType* result)
+  {
+     if( 0 == result->coeffRef(0,0) ) return false;
+   
+     typedef typename MatrixType::Scalar Scalar;
+     result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
+     return true;
+  }
+};
+
+template
+struct ei_compute_inverse_with_check
+{
+  static inline bool run(const MatrixType& matrix, MatrixType* result)
+  {
+    return ei_compute_inverse_in_size2_case_with_check(matrix, result);
+  }
+};
+
+template
+struct ei_compute_inverse_with_check
+{
+  static inline bool run(const MatrixType& matrix, MatrixType* result)
+  {
+    return ei_compute_inverse_in_size3_case_with_check(matrix, result);
+  }
+};
+
+template
+struct ei_compute_inverse_with_check
+{
+  static inline bool run(const MatrixType& matrix, MatrixType* result)
+  {
+    return ei_compute_inverse_in_size4_case_with_check(matrix, result);
+  }
+};
+
+/** lu_module
+  *
+  * If the matrix is invertible, computes the matrix inverse of this matrix
+  * and returns true otherwise return false.
+  *
+  * note This matrix must be invertible, otherwise the result is undefined.
+  *
+  * param result Pointer to the matrix in which to store the result. Undefined
+  *  if the matrix is not invertible.
+  * return true if the matrix is invertible false otherwise.
+  *
+  * sa inverse()
+  */
+template
+inline bool MatrixBase::computeInverseWithCheck(PlainMatrixType *result) const
+{
+  ei_assert(rows() == cols());
+  EIGEN_STATIC_ASSERT(NumTraits::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT)
+  return ei_compute_inverse_with_check::run(eval(), result);
+}
+
+
 #endif // EIGEN_INVERSE_H


OK, thanks for pinging me on IRC, I was forgetting about this.

I looked at your patch but it's looking like some <...> were eaten away ? like after "template"