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; } } }
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; }
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??
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.