Re: [eigen] Problem inverting a Matrix4f with Eigen 2.0.0

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Your matrix is indeed not very well conditioned for floats. Are you
using fixed size matrices (Matrix4f) or dynamically sized ones
(MatrixXf) ? I'm asking because there is a special path to invert 4x4
fixed size matrices, and as far as I remember there was bug for some
special configurations in the early versions of Eigen 2.0, so I would
recommend you to upgrade it to the latest. Then if you use Eigen2.0 LU
(the default for MatrixXf), then I confirm it detects your matrix is
not well conditioned enough for float and stops the decomposition.
Here are the results of A^1 * A using Eigen 2.0.15: (LU = full
pivoting):

4x4 float:
1 0.001343 0.002823 -0.0004272
-1.268e-11 1 -1.332e-07 1.49e-08
1.609e-11 7.916e-08 1 -1.49e-08
-9.095e-13 -1.49e-08 -2.98e-08 1

4x4 double:
1 3.411e-13 1.407e-12 -1.137e-13
-1.133e-20 1 -2.047e-16 2.776e-17
-2.393e-20 -2.03e-16 1 2.776e-17
1.694e-21 4.163e-17 5.551e-17 1

LU float:
-3.681e-32 -2.173e-28 -7.705e-28 8.745e-29
4.337e-42 2.561e-38 9.078e-38 -1.03e-38
-1.039e+11 -1.176e+15 -1.503e+15 2.467e+14
7.456e-42 8.063e-38 1.056e-37 8.356e-38

LU double:
1 3.268e-13 1.393e-12 2.274e-13
7.623e-21 1 1.908e-17 0
1.673e-20 1.301e-16 1 -2.776e-17
2.541e-21 1.388e-17 5.551e-17 1

LLT float:
1 0.0006409 0.001884 -0.0003662
6.821e-13 1 -6.519e-08 2.98e-08
8.583e-12 1.583e-08 1 2.235e-08
4.547e-13 0 2.235e-08 1

For the record here is with Eigen 3 (LU = Partial pivoting)

4x4 float:
           1 -7.62939e-05 -5.34058e-05  -0.00012207
-9.77707e-12            1 -3.91155e-08  1.49012e-08
  1.3074e-12 -4.37722e-08            1 -7.45058e-09
 4.54747e-13  7.45058e-09  1.49012e-08            1

4x4 double:
           1 -1.93268e-12 -4.05009e-12  3.41061e-13
-1.40819e-20            1 -2.61943e-16  2.77556e-17
-1.03762e-20 -3.64292e-17            1  1.38778e-17
           0  1.38778e-17  2.77556e-17            1

LU float:
           1 -0.000488281   0.00012207   -0.0211182
-1.59162e-12            1 -5.40167e-08   8.9407e-07
 5.00222e-12 -1.67638e-08            1  7.52509e-07
-4.54747e-13            0            0            1

LU double:
           1 -1.26334e-11 -1.99378e-11   6.9349e-11
 4.36222e-20            1  7.70217e-16 -2.96985e-15
 2.83756e-20  3.33067e-16            1 -2.44249e-15
 8.47033e-22  1.38778e-17  2.77556e-17            1

LLT float:
           1 -0.000335693  0.000915527 -0.000427246
 7.95808e-12            1 -5.58794e-09            0
  1.3074e-12 -1.39698e-08            1  2.98023e-08
 4.54747e-13            0  2.23517e-08            1


cheers,
gael.

On Sat, Jul 17, 2010 at 10:01 AM, Christoph Hertzberg
<chtz@xxxxxxxxxxxxxxxxxxxxxxxx> wrote:
> Robert Lupton the Good wrote:
>>
>> I'm using Eigen 2.0.0 to invert some small (4x4) matrices --- I need
>> the covariance matrix, so there's no need to tell me to
>> solve-not-invert!
>>
>> Anyway, with floats, we fail to invert although the determinant is
>> decent, ~ 3.6e4.  Is there a better way to check for a failed
>> inversion?  I didn't see anything in the manual.
>
> determinant != condition:
> http://en.wikipedia.org/wiki/Condition_number
>
> Matlab says the condition of your matrix is 3.8e8, which is pretty bad.
> (btw: is there a way to calculate condition with Eigen?)
>
> I think invert() uses LU-Decomposition. Have you tried using slightly more
> stable Cholesky instead? I guess something like this could work:
>
> Matrix4f inverse;
> inverse.setIdentity();
> fisher.ldlt().solveInPlace(inverse);
>
> hth, Christoph
>
>
> --
> ----------------------------------------------
> Dipl.-Inf. Christoph Hertzberg
> Cartesium 0.051
> Universität Bremen
> Enrique-Schmidt-Straße 5
> 28359 Bremen
>
> Tel: (+49) 421-218-64252
> ----------------------------------------------
>
>
>

#include <iostream>
#include <Eigen/LU>
#include <Eigen/Cholesky>
using namespace Eigen;

int main()
{
  Matrix4d m;
  m << 1.353e-05, 0.07932, 0.09445, -0.01072,
      0.07932, 1394, 557.7, -188.4,
      0.09445, 557.7, 1977, -224.4,
      -0.01072, -188.4, -224.4, 2231;

  std::cout << "4x4 float:\n" << m.cast<float>().inverse() << "\n\n";
  std::cout << "4x4 double:\n" << m.inverse() << "\n\n";
  
  std::cout << "LU float:\n" << m.cast<float>().lu().inverse() << "\n\n";
  std::cout << "LU double:\n" << m.lu().inverse() << "\n\n";
  Matrix4f invf = Matrix4f::Identity();
  m.cast<float>().llt().solveInPlace(invf);
  std::cout << "LLT float:\n" << invf << "\n\n";

  std::cout << "4x4 float:\n" << m.cast<float>().inverse()*m.cast<float>() << "\n\n";
  std::cout << "4x4 double:\n" << m.cast<double>().inverse()*m << "\n\n";
  std::cout << "LU float:\n" << m.cast<float>().lu().inverse()*m.cast<float>() << "\n\n";
  std::cout << "LU double:\n" << m.lu().inverse()*m << "\n\n";
  std::cout << "LLT float:\n" << invf*m.cast<float>() << "\n\n";
  
  return 0;
}