(Pseudo)-Inverse of N by N matrix with zero determinant

Toon picture Toon · Dec 17, 2012 · Viewed 11.7k times · Source

I would like to take the inverse of a nxn matrix to use in my GraphSlam.

The issues that I encountered:

  • .inverse() Eigen-library (3.1.2) doesn't allow zero values, returns NaN
  • The LAPACK (3.4.2) library doesn't allow to use a zero determinant, but allows zero values (used example code from Computing the inverse of a matrix using lapack in C)
  • Seldon library (5.1.2) wouldn't compile for some reason

Did anyone successfully implemented an n x n matrix inversion code that allows negative, zero-values and a determinant of zero? Any good library (C++) recommendations?

I try to calculate the omega in the following for GraphSlam: http://www.acastano.com/others/udacity/cs_373_autonomous_car.html


Simple example:

[ 1 -1  0 0 ]
[ -1 2 -1 0 ]
[ 0 -1  1 0 ]
[ 0  0  0 0 ]

Real example would be 170x170 and contain 0's, negative values, bigger positive values. Given simple example is used to debug the code.


I can calculate this in matlab (Moore-Penrose pseudoinverse) but for some reason I'm not able to program this in C++.

A = [1 -1 0 0; -1 2 -1 0; 0 -1 1 0; 0 0 0 0]
B = pinv(A)
B=
[0.56   -0.12  -0.44  0]
[-0.12  0.22   -0.11  0]
[-0.44  -0.11   0.56  0]
[0  0  0   0]

For my application I can (temporarily) remove the dimension with zero's.
So I am going to remove the 4th column and the 4th row.
I can also do that for my 170x170 matrix, the 4x4 was just an example.

A:

[ 1 -1  0 ]
[ -1 2 -1 ]
[ 0 -1  1 ]

So removing the 4th column and the 4th row wouldn’t bring a zero determinant. But I can still have a zero determinant if my matrix is as above. This when the sum of each row or each column is zero. (Which I will have all the time in GraphSlam)

The LAPACK-solution (Moore-Penrose Inverse based) worked if the determinant was not zero (used example code from Computing the inverse of a matrix using lapack in C).
But failed as a "pseudoinverse" with a determinant of zero.


SOLUTION: (all credits to Frank Reininghaus), using SVD(singular value decomposition)
http://sourceware.org/ml/gsl-discuss/2008-q2/msg00013.html

Works with:

  • Zero values (even full 0 rows and full 0 columns)
  • Negative values
  • Determinant of zero

A^-1:

[0.56   -0.12  -0.44]
[-0.12  0.22   -0.11]
[-0.44  -0.11   0.56]

Answer

ggael picture ggael · Dec 18, 2012

If all you want is to solve problem of the form Ax=B (or equivalently compute products of the form A^-1 * b), then I recommend you not to compute the inverse or pseudo-inverse of A, but directly solve for Ax=b using an appropriate rank-revealing solver. For instance, using Eigen:

x = A.colPivHouseholderQr().solve(b);
x = A.jacobiSvd(ComputeThinU|ComputeThinV).solve(b);