# Effective UBLAS/Matrix Inversion

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

Showing revision 13
A more simple and efficient version is here LU Matrix Inversion

### A Matrix Inverse implementation

I wrote a simple matrix inversion algorithm that uses uBlas. I don't know wether it's usefull to you? Any comments? I'm new to uBlas and looking for a usefull inversion algorithm using uBlas that works good on dense 100x100 matrices. (I'malso trying to use Wiki for the first time, so don't bother all the edits.)

The result of the following function is that E2 becomes E1^(-1)*E2 and E1 becomes a unit matrix.

``` template<class E1, class E2>
void inverse (matrix_expression<E1> &e1, matrix_expression<E2> &e2) {
typedef BOOST_UBLAS_TYPENAME E2::size_type size_type;
typedef BOOST_UBLAS_TYPENAME E2::difference_type difference_type;
typedef BOOST_UBLAS_TYPENAME E2::value_type value_type;
```

```    BOOST_UBLAS_CHECK (e1 ().size1 () == e2 ().size1 (), bad_size ());
BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size2 (), bad_size ());
size_type size = e1 ().size1 ();
for (size_type n = 0; n < size; ++ n) {
// processing column n
// find the row that has the largest number at this column (in absolute value)
size_type best_row = index_norm_inf(row(e1(), n));
// check wether this number is'nt zero
BOOST_UBLAS_CHECK (e1 () (best_row, n) != value_type (), singular ());
```

```       { // swap this row with the n-th row
vector<value_type> temp = row(e1(), best_row);
row(e1(), n) = row(e1(), best_row);
row(e1(), best_row) = temp;
}
// do the same on the destination matrix
{ // swap this row with the n-th row
vector<value_type> temp = row(e2(), best_row);
row(e2(), n) = row(e2(), best_row);
row(e2(), best_row) = temp;
}
```

```       // now eliminate all elements below and above this row
for (size_type i = 0; i < size; ++ i)
if (i!=n) {
value_type t = -e1 () (i, n)/ e1 () (n, n);
row(e1(), i) += t*row(e1(), n);
row(e2(), i) += t*row(e2(), n);
} else {
value_type t = 1 / e1 () (i, n);
row(e1(), i) *= t;
row(e2(), i) *= t;
}
}
}
```

### Another Matrix Inverse implementation

I (Yi Wang, wangy01@mails.tsinghua.edu.cn) used the above algorithm, but cannot always get the right answer. I had not debug the above code, but wrote a new one. This code uses Gauss-Jordan decomposition with PARTIAL pivot. The code is tested by (1) inversing random permutations of rows of an identity matrix, and (2) inversing randomly generated large matrices (by MATLAB), and comparing the result with that of MATLAB inv(...) command.

``` #include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix_expression.hpp>
```

``` #include <iostream>
#include <fstream>
#include <vector>
```

``` /**
* Invert a matrix via gauss-jordan algorithm (PARTIAL PIVOT)
*
* @param m The matrix to invert. Must be square.
* @param singular If the matrix was found to be singular, then this
*        is set to true, else set to false.
* @return If singular is false, then the inverted matrix is returned.
*         Otherwise it contains random values.
*/
template<class T>
//#define T double /// for debug
boost::numeric::ublas::matrix<T>
gjinverse(const boost::numeric::ublas::matrix<T> &m,
bool &singular)
{
using namespace boost::numeric::ublas;
```

```     const int size = m.size1();
```

```     // Cannot invert if non-square matrix or 0x0 matrix.
// Report it as singular in these cases, and return
// a 0x0 matrix.
if (size != m.size2() || size == 0)
{
singular = true;
matrix<T> A(0,0);
return A;
}
```

```     // Handle 1x1 matrix edge case as general purpose
// inverter below requires 2x2 to function properly.
if (size == 1)
{
matrix<T> A(1, 1);
if (m(0,0) == 0.0)
{
singular = true;
return A;
}
singular = false;
A(0,0) = 1/m(0,0);
return A;
}
```

```     // Create an augmented matrix A to invert. Assign the
// matrix to be inverted to the left hand side and an
// identity matrix to the right hand side.
matrix<T> A(size, 2*size);
matrix_range<matrix<T> > Aleft(A,
range(0, size),
range(0, size));
Aleft = m;
matrix_range<matrix<T> > Aright(A,
range(0, size),
range(size, 2*size));
Aright = identity_matrix<T>(size);
```

```     // Swap rows to eliminate zero diagonal elements.
for (int k = 0; k < size; k++)
{
if ( A(k,k) == 0 ) // XXX: test for "small" instead
{
// Find a row(l) to swap with row(k)
int l = -1;
for (int i = k+1; i < size; i++)
{
if ( A(i,k) != 0 )
{
l = i;
break;
}
}
```

```             // Swap the rows if found
if ( l < 0 )
{
std::cerr << "Error:" <<  __FUNCTION__ << ":"
<< "Input matrix is singular, because cannot find"
<< " a row to swap while eliminating zero-diagonal.";
singular = true;
return Aleft;
}
else
{
matrix_row<matrix<T> > rowk(A, k);
matrix_row<matrix<T> > rowl(A, l);
rowk.swap(rowl);
```

``` #if defined(DEBUG) || !defined(NDEBUG)
std::cerr << __FUNCTION__ << ":"
<< "Swapped row " << k << " with row " << l
<< ":" << A << "\n";
#endif
}
}
}
```

```     // Doing partial pivot
for (int k = 0; k < size; k++)
{
// normalize the current row
for (int j = k+1; j < 2*size; j++)
A(k,j) /= A(k,k);
A(k,k) = 1;
```

```         // normalize other rows
for (int i = 0; i < size; i++)
{
if ( i != k )  // other rows  // FIX: PROBLEM HERE
{
if ( A(i,k) != 0 )
{
for (int j = k+1; j < 2*size; j++)
A(i,j) -= A(k,j) * A(i,k);
A(i,k) = 0;
}
}
}
```

``` #if defined(DEBUG) || !defined(NDEBUG)
std::cerr << __FUNCTION__ << ":"
<< "GJ row " << k << " : " << A << "\n";
#endif
}
```

```     singular = false;
return Aright;
}
```

``` /**
* Load a matrix from a ASCII text file generated by the MATLAB
* command:
*   A = rand(500);
*   disp(det(A));
*   dlmwrite('A.dlm', A, 'delimiter', ' ');
*
* @param m The matrix to read in.
* @param filename The text filename, e.g., "A.dlm"
* @return Whether the read operation success.
*/
template <class T>
//#define T double /// for debug
boost::numeric::ublas::matrix<T> &m )
{
std::vector<double> v;
```

```     std::ifstream ifs(filename);
```

```     if ( ifs.bad() )
{
std::cerr << "Error:" << __FUNCTION__ << ":"
<< "Cannot open file : " << filename << "\n";
return false;
}
```

```     while ( ifs.good() )
{
double val;
ifs >> val;
```

```         if ( ! ifs.eof() )
v.push_back(val);
else
break;
}
```

```     if ( ifs.eof() )
{
double size = sqrt( (double)v.size() );
if ( size * size == v.size() )
{
m.resize( (int)size, (int)size );
int k = 0;
for (int i = 0; i < size; i++)
for (int j = 0; j < size; j++)
m(i,j) = v[k++];
}
}
```

```     return true;
}
```

``` int main()
{
using namespace boost::numeric::ublas;
```

``` //     matrix<double> e1(2,2);
//     e1(0,0)=0;
//     e1(0,1)=2;
//     e1(1,0)=4;
//     e1(1,1)=4;
```

``` //     matrix<double> e1(3,3);
//     int k = 0;
//     static double mat_val[] = { 0,0,1, 1,0,0, 0,1,0 };
//     for (int i = 0; i < e1.size1(); i++)
//         for (int j = 0; j < e1.size2(); j++)
//             e1(i,j) = mat_val[k++];
```

```     matrix<double> e1;
```

```     std::cout << "e1:" << e1 << std::endl;
```

```     bool singular;
matrix<double> e2 = gjinverse(e1, singular);
```

```     if ( singular )
std::cout << "singluar matrix cannot be inverted.\n";
else
{
std::cout << "e2:" << e2 << "\n";
```

```         matrix<double> I = prod(e1, e2);
std::cout << "prod(e1, e2):\n";
for (int i = 0; i < I.size1(); i++)
{
for (int j = 0; j < I.size2(); j++ )
{
if ( I(i,j) < 1e-10 )
I(i,j) = 0;
std::cout << I(i,j) << ",";
}
std::cout << "\n";
}
}
```

```     return 0;
}
```

### Matrix inversion using LU decomposition

Check LU Matrix Inversion for an example on how to invert matrices using the LU-decomposition functions found in uBLAS.

### Correct Another Matrix Inverse implementation

If you try to use both above methods for Sparse Matrix you get a message about singular, but current matrix is not singular!

For example,

``` Q =
```

``` 1	0	0	0	0	0	0	0
```

``` 1	1	1	1	0	0	0	0
```

``` 0	1	2	3	0	-1	0	0
```

``` 0	0	1	3	0	0	-1	0
```

``` 0	0	0	0	1	0	0	0
```

``` 0	0	0	0	1	1	1	1
```

``` 0	0	1	0	0	0	0	0
```

``` 0	0	0	0	0	0	1	3
```

``` Q^{-1} =
```

``` 1	0	0	0	 0	 0	 0	  0
```

``` -1.25	1.25	-0.25	-0.167 0.25	 -0.25 -0.583 0.083
```

``` 0	0	0	0	 0	 0	 1	  0
```

``` 0.25	-0.25	0.25	0.167	 -0.25 0.25	 -0.417 -0.083
```

``` 0	0	0	0	 1	 0	 0	  0
```

``` -0.5	0.5	-0.5	0.333	 -0.5	 0.5	 0.167  -0.167
```

``` 0.75	-0.75	0.75	-0.5	 -0.75 0.75	 -0.25  -0.25
```

``` -0.25	0.25	-0.25	0.167	 0.25	 -0.25 0.083  0.417
```

Code of swap rows is incorrect placement and is result from this error.

I (Andrey Zakharov,anzakhar@narod.ru) correct the "Another Matrix Inverse implementation" by Yi Wang so:

``` #include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix_expression.hpp>
```

``` #include <iostream>
#include <fstream>
#include <vector>
```

``` /**
* Invert a matrix via gauss-jordan algorithm (PARTIAL PIVOT)
*
* @param m The matrix to invert. Must be square.
* @param singular If the matrix was found to be singular, then this
*        is set to true, else set to false.
* @return If singular is false, then the inverted matrix is returned.
*         Otherwise it contains random values.
*/
template<class T>
//#define T double /// for debug
boost::numeric::ublas::matrix<T>
gjinverse(const boost::numeric::ublas::matrix<T> &m,
bool &singular)
{
using namespace boost::numeric::ublas;
```

```     const int size = m.size1();
```

```     // Cannot invert if non-square matrix or 0x0 matrix.
// Report it as singular in these cases, and return
// a 0x0 matrix.
if (size != m.size2() || size == 0)
{
singular = true;
matrix<T> A(0,0);
return A;
}
```

```     // Handle 1x1 matrix edge case as general purpose
// inverter below requires 2x2 to function properly.
if (size == 1)
{
matrix<T> A(1, 1);
if (m(0,0) == 0.0)
{
singular = true;
return A;
}
singular = false;
A(0,0) = 1/m(0,0);
return A;
}
```

```     // Create an augmented matrix A to invert. Assign the
// matrix to be inverted to the left hand side and an
// identity matrix to the right hand side.
matrix<T> A(size, 2*size);
matrix_range<matrix<T> > Aleft(A,
range(0, size),
range(0, size));
Aleft = m;
matrix_range<matrix<T> > Aright(A,
range(0, size),
range(size, 2*size));
Aright = identity_matrix<T>(size);
```

``` /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// INCORRECT PLACEMENT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//
//     // Swap rows to eliminate zero diagonal elements.
//     for (int k = 0; k < size; k++)
//     {
//         if ( A(k,k) == 0 ) // XXX: test for "small" instead
//         {
//             // Find a row(l) to swap with row(k)
//             int l = -1;
//             for (int i = k+1; i < size; i++)
//             {
//                 if ( A(i,k) != 0 )
//                 {
//                     l = i;
//                     break;
//                 }
//             }
//
//             // Swap the rows if found
//             if ( l < 0 )
//             {
//                 std::cerr << "Error:" <<  __FUNCTION__ << ":"
//                           << "Input matrix is singular, because cannot find"
//                           << " a row to swap while eliminating zero-diagonal.";
//                 singular = true;
//                 return Aleft;
//             }
//             else
//             {
//                 matrix_row<matrix<T> > rowk(A, k);
//                 matrix_row<matrix<T> > rowl(A, l);
//                 rowk.swap(rowl);
//
// #if defined(DEBUG) || !defined(NDEBUG)
//                 std::cerr << __FUNCTION__ << ":"
//                           << "Swapped row " << k << " with row " << l
//                           << ":" << A << "\n";
// #endif
//             }
//         }
//     }
//
/////////////////////////////////////////////////////////////////////////////////////////////////////////
```

```     // Doing partial pivot
for (int k = 0; k < size; k++)
{
///////////////////////////////////////////////////////////////////////////////////////////////////////////
// RIGHT PLACEMENT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
```

```     // Swap rows to eliminate zero diagonal elements.
for (int k = 0; k < size; k++)
{
if ( A(k,k) == 0 ) // XXX: test for "small" instead
{
// Find a row(l) to swap with row(k)
int l = -1;
for (int i = k+1; i < size; i++)
{
if ( A(i,k) != 0 )
{
l = i;
break;
}
}
```

```             // Swap the rows if found
if ( l < 0 )
{
std::cerr << "Error:" <<  __FUNCTION__ << ":"
<< "Input matrix is singular, because cannot find"
<< " a row to swap while eliminating zero-diagonal.";
singular = true;
return Aleft;
}
else
{
matrix_row<matrix<T> > rowk(A, k);
matrix_row<matrix<T> > rowl(A, l);
rowk.swap(rowl);
```

``` #if defined(DEBUG) || !defined(NDEBUG)
std::cerr << __FUNCTION__ << ":"
<< "Swapped row " << k << " with row " << l
<< ":" << A << "\n";
#endif
}
}
}
```

``` /////////////////////////////////////////////////////////////////////////////////////////////////////////
// normalize the current row
for (int j = k+1; j < 2*size; j++)
A(k,j) /= A(k,k);
A(k,k) = 1;
```

```         // normalize other rows
for (int i = 0; i < size; i++)
{
if ( i != k )  // other rows  // FIX: PROBLEM HERE
{
if ( A(i,k) != 0 )
{
for (int j = k+1; j < 2*size; j++)
A(i,j) -= A(k,j) * A(i,k);
A(i,k) = 0;
}
}
}
```

``` #if defined(DEBUG) || !defined(NDEBUG)
std::cerr << __FUNCTION__ << ":"
<< "GJ row " << k << " : " << A << "\n";
#endif
}
```

```     singular = false;
return Aright;
}
```

BOOST WIKI | Effective UBLAS | RecentChanges | Preferences | Page List | Links List
Edit revision 13 of this page | View other revisions | View current revision
Edited March 1, 2007 8:04 am (diff)
Search:
Disclaimer: This site not officially maintained by Boost Developers