Pulse Inversion

Question

I am trying to perform a simple matrix inverse operation using boost. But I get an error message. Basically what I'm trying to find is inversted_matrix = inverse matrix (trans (matrix) *) But I get 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; } 

Answer 1

Boost Ublas has runtime checks to ensure, among other things, numerical stability. If you look at the source of the error, you will see that it is trying to make sure that U*X = B, X = U^-1*B, U*X = B (or something like that) are corrected up to some epsilon. If you have too many deviations, it most likely will not.

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

Answer 2

Since m has only 4 rows, prod(trans(m), m) cannot have a rank higher than 4, and since the product is a 5x5 matrix, it must be singular (i.e., has determinant 0) and calculates the inverse singularity. The matrix is ​​similar to division by 0. Add independent lines to m to solve this singularity problem.