[Home]Effective UBLAS/Matrix Inversion

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

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
 bool readDLMmatrix(const char * filename, 
                    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;
     readDLMmatrix("A.dlm", 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;
 }

Warning!!!!! Compiling your code with

gcc version 3.4.6 20060404 (Red Hat 3.4.6-3)

Gives me these warning:

./mysource.cpp:232: warning: name lookup of `k' changed ./mysource.cpp:191: warning: matches this `k' under ISO standard rules ./mysource.cpp:197: warning: matches this `k' under old rules

Line 232: for (int j = k+1; j < 2*size; j++)

Line 190: // Doing partial pivot Line 191: for (int k = 0; k < size; k++)

Line 196: // Swap rows to eliminate zero diagonal elements. Line 197: for (int k = 0; k < size; k++)

So, what k are you using at line 232? the one at 191 or the one at 197? Can you change your variables??

A Matrix Inverse implementation using QR decomposition

I (Eustache Diemert, eustache_at_diemert_dot_fr), gives this example code to compute square matrix inverse using QR decomposition (more precisely Householder reflections). Here it is:

 template<class T>
 void TransposeMultiply? (const ublas::vector<T>& vector, 
		         ublas::matrix<T>& result,
		         size_t size)
 {
  result.resize (size,size);
  result.clear ();
  for(unsigned int row=0; row< vector.size(); ++row)
    {
      for(unsigned int col=0; col < vector.size(); ++col)
	result(row,col) = vector(col) * vector(row);

    }
 }

 template<class T>
 void HouseholderCornerSubstraction? (ublas::matrix<T>& LeftLarge?, 
				   const ublas::matrix<T>& RightSmall?)
 {
  using namespace boost::numeric::ublas;
  using namespace std; 
  if( 
     !( 
          (LeftLarge?.size1() >= RightSmall?.size1())
       && (LeftLarge?.size2() >= RightSmall?.size2())
       ) 
     )
    {
      cerr << "invalid matrix dimensions" << endl;
      return;
    }  

  size_t row_offset = LeftLarge?.size2() - RightSmall?.size2();
  size_t col_offset = LeftLarge?.size1() - RightSmall?.size1();

  for(unsigned int row = 0; row < RightSmall?.size2(); ++row )
    for(unsigned int col = 0; col < RightSmall?.size1(); ++col )
      LeftLarge?(col_offset+col,row_offset+row) -= RightSmall?(col,row);
 }

 template<class T>
 void HouseholderQR? (const ublas::matrix<T>& M, 
		   ublas::matrix<T>& Q, 
		   ublas::matrix<T>& R)
 {
  using namespace boost::numeric::ublas;
  using namespace std;  

  if( 
     !( 
           (M.size1() == M.size2())
      ) 
    )
    {
      cerr << "invalid matrix dimensions" << endl;
      return;
    }
  size_t size = M.size1();

  // init Matrices
  matrix<T> H, HTemp;
  HTemp = identity_matrix<T>(size);
  Q = identity_matrix<T>(size);
  R = M;

  // find Householder reflection matrices
  for(unsigned int col = 0; col < size-1; ++col)
    {
      // create X vector
      ublas::vector<T> RRowView? = column(R,col);      
      vector_range< ublas::vector<T> > X2 (RRowView?, range (col, size));
      ublas::vector<T> X = X2;

      // X -> U~
      if(X(0) >= 0)
	X(0) += norm_2(X);
      else
	X(0) += -1*norm_2(X);      

      HTemp.resize(X.size(),X.size(),true);

      TransposeMultiply?(X, HTemp, X.size());

      // HTemp = the 2UUt part of H 
      HTemp *= ( 2 / inner_prod(X,X) );

      // H = I - 2UUt
      H = identity_matrix<T>(size);
      HouseholderCornerSubstraction?(H,HTemp);

      // add H to Q and R
      Q = prod(Q,H);
      R = prod(H,R);
    }
 }

And a dummy main() to test it :

 int main()
 {
  using namespace boost::numeric::ublas;
  using namespace std;
  matrix<double> A (3,3);
  A(0,0) = 1;
  A(0,1) = 1;
  A(0,2) = 0;
  A(1,1) = 1;
  A(1,0) = 0;
  A(1,2) = 0;
  A(2,2) = 1;
  A(2,0) = 1;
  A(2,1) = 0; 
  cout << "A=" << A << endl;

  cout << "QR decomposition using Householder" << endl;
  matrix<double> Q(3,3), R(3,3);
  HouseholderQR? (A,Q,R);
  matrix<double> Z = prod(Q,R) - A;
  float f = norm_1 (Z);
  cout << "Q=" << Q <<endl;
  cout << "R=" << R << endl;
  cout << "|Q*R - A|=" << f << endl;

  return 0;
 }

f value should be approximately 4e-16.


BOOST WIKI | Effective UBLAS | RecentChanges | Preferences | Page List | Links List
Edit text of this page | View other revisions
Last edited May 2, 2007 5:25 pm (diff)
Search:
Disclaimer: This site not officially maintained by Boost Developers