LU Matrix Inversion

BOOST WIKI | RecentChanges | Preferences | Page List | Links List

Difference (from prior major revision) (minor diff, author diff)

Changed: 18c18
 void InvertMatrix? (const ublas::matrix& input, ublas::matrix& inverse) {
 bool InvertMatrix? (const ublas::matrix& input, ublas::matrix& inverse) {

Changed: 27,28c27,29
 lu_factorize(A,pm);
 int res = lu_factorize(A,pm); if( res != 0 ) return false;

 return true;

Changed: 40c43
 This code snipet is concise and effective, but when input matrix is singular, the function simply crashes. A more complex version is at the second section of Effective UBLAS/Matrix Inversion with the capability to detect and report singularity of the input matrix.
 This version is efficient as using LU. A more complex and probably inefficient version is at the second section of Effective UBLAS/Matrix Inversion.

Removed: 42d44
 edit: lu_factorize() returns 0 if it was successful. It returns (k+1) if it detects singularity after processing row k. So one should always check its return value.

The following code inverts the matrix input using LU-decomposition with backsubstitution of unit vectors. Reference: Numerical Recipies in C, 2nd ed., by Press, Teukolsky, Vetterling & Flannery.
``` #ifndef INVERT_MATRIX_HPP
#define INVERT_MATRIX_HPP
```

``` // REMEMBER to update "lu.hpp" header includes from boost-CVS
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>
```

``` namespace ublas = boost::numeric::ublas;
```

``` /* Matrix inversion routine.
Uses lu_factorize and lu_substitute in uBLAS to invert a matrix */
template<class T>
bool InvertMatrix? (const ublas::matrix<T>& input, ublas::matrix<T>& inverse) {
using namespace boost::numeric::ublas;
typedef permutation_matrix<std::size_t> pmatrix;
// create a working copy of the input
matrix<T> A(input);
// create a permutation matrix for the LU-factorization
pmatrix pm(A.size1());
```

``` 	// perform LU-factorization
int res = lu_factorize(A,pm);
if( res != 0 ) return false;
```

``` 	// create identity matrix of "inverse"
inverse.assign(ublas::identity_matrix<T>(A.size1()));
```

``` 	// backsubstitute to get the inverse
lu_substitute(A, pm, inverse);
```

``` 	return true;
}
```

``` #endif //INVERT_MATRIX_HPP
```
Hope someone finds this useful. Regards, Fredrik Orderud.

Matrix inversion with Singularity Detection

This version is efficient as using LU. A more complex and probably inefficient version is at the second section of Effective UBLAS/Matrix Inversion.

BOOST WIKI | RecentChanges | Preferences | Page List | Links List