CODEWITHSUNDEEP
c++boost
frazman
frazman

Reputation: 33273

matrix inversion in boost

I am trying to do a simple matrix inversion operation using boost. But I am getting an error. Basically what I am trying to find is inversted_matrix = inverse(trans(matrix) * matrix) But I am getting an error 

Check failed in file boost_1_53_0/boost/numeric/ublas/lu.hpp at line 299: 
detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, 
upper> (m), e), cm2) 
terminate called after throwing an instance of 
'boost::numeric::ublas::internal_logic' 
  what(): internal logic 
Aborted (core dumped)

My attempt: 

#include <boost/numeric/ublas/matrix.hpp> 
#include <boost/numeric/ublas/vector.hpp> 
#include <boost/numeric/ublas/io.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> 

namespace ublas = boost::numeric::ublas; 
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; 
}

int main(){ 
    using namespace boost::numeric::ublas; 
    matrix<double> m(4,5); 
    vector<double> v(4); 
    vector<double> thetas; 
    m(0,0) = 1; m(0,1) = 2104; m(0,2) = 5; m(0,3) = 1;m(0,4) = 45; 
    m(1,0) = 1; m(1,1) = 1416; m(1,2) = 3; m(1,3) = 2;m(1,4) = 40; 
    m(2,0) = 1; m(2,1) = 1534; m(2,2) = 3; m(2,3) = 2;m(2,4) = 30; 
    m(3,0) = 1; m(3,1) = 852; m(3,2) = 2; m(3,3) = 1;m(3,4) = 36; 
    std::cout<<m<<std::endl; 
    matrix<double> product = prod(trans(m), m); 
    std::cout<<product<<std::endl; 
    matrix<double> inversion(5,5); 
    bool inverted; 
    inverted = InvertMatrix(product, inversion); 
    std::cout << inversion << std::endl; 
}

Upvotes: 6

Views: 12115

Answers (3)

Cloud Cho
Cloud Cho

Reputation: 1774

I think your matrix dimension, 4 by 5, caused the error. Like what Maarten Hilferink mentioned, you may try with a square matrix like 5 by 5. Here are requirement to have an inverse:

  1. The matrix must be square (same number of rows and columns).
  2. The determinant of the matrix must not be zero (determinants are covered in section 6.4). This is instead of the real number not being zero to have an inverse, the determinant must not be zero to have an inverse.

Upvotes: 0

Maarten Hilferink
Maarten Hilferink

Reputation: 723

As m has only 4 rows, prod(trans(m), m) cannot have a rank higher than 4, and as the product is a 5x5 matrix, it must be singular (i.e. it has determinant 0) and calculating the inverse of a singular matrix is like division by 0. Add independent rows to m to solve this singularity problem.

Upvotes: 5

Anycorn
Anycorn

Reputation: 51505

Boost Ublas has runtime checks to ensure among other thing numerical stability. If you look at source of the error, you can see that it tries to make sure that U*X = B, X = U^-1*B, U*X = B (or smth like that) are coorect to within some epsilon. If you have too much deviation numerically this will likely not hold.

You can disable checks via -DBOOST_UBLAS_NDEBUG or twiddle with BOOST_UBLAS_TYPE_CHECK_EPSILON, BOOST_UBLAS_TYPE_CHECK_MIN.

Upvotes: 7

Related Questions