Search code examples
c++matrixboostmatrix-multiplicationboost-ublas

Error in using Boost Ublas lu_Factorize

I am trying to make a fun project which makes a matrix power function using BOOST Ublas.It is very similar to this Numpy library power function for matrix It uses matrix exponentiation to compute the nth power of a matrix in logarithmic time. It has 3 cases to compute the nth power of a matrix-:

  1. If power > 0 use matrix exponentiation directly
  2. If power = 0 check if matrix has an inverse(check Using lu_factorize)and if yes,return identity matrix
  3. If Power < 0 find the inverse(if it exists) then use matrix exponentiation on it

I am good with algorithms and implementation but this is the first time i am using any open source library, i want to learn this,so that i could eventually contribute to boost.

This is my header file 

//  Distributed under the Boost Software License, Version 1.0. (See
//  accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt)
//
#ifndef BOOST_UBLAS_POW_HPP
#define BOOST_UBLAS_POW_HPP 
#include <iostream>
#include <vector>
#include <iomanip>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/mpl/set.hpp>
#include <boost/mpl/assert.hpp>
#include <boost/multiprecision/cpp_int.hpp>
using namespace boost::numeric::ublas;
namespace boost { namespace numeric { namespace ublas {

    typedef permutation_matrix<std::size_t> pmatrix;
    template< typename T,typename T2 >
        matrix<T> matrix_power(const matrix<T> input, T2 exponent)
        {
          matrix<T> resultant=input;
          BOOST_ASSERT_MSG(input.size1()==input.size2(),"Not a square matrix\n");
           if(exponent>0) 
            resultant=matrix_exponent(input,exponent);// this where you could directly use matrix exponentiation 

           else if(exponent==0)
            {
             pmatrix pm(input.size1());
             matrix<T> A(input);
             BOOST_ASSERT_MSG(lu_factorize(A, pm)==0,"Attempted to compute a  power of 0 in a matrix without inverse\n");
             resultant.assign(identity_matrix<T> (input.size1()));// if the matrix is invertible output identity matrix for the 0t hpower
            }
            else {
              matrix<T> A(input);
              pmatrix pm(A.size1());
              BOOST_ASSERT_MSG(lu_factorize(A, pm)==0,"Attempted to compute inverse in a singular matrix\n");
              resultant.assign(identity_matrix<T> (A.size1()));
              lu_substitute(A, pm, resultant);
              resultant=matrix_exponent(resultant,std::abs(exponent));
            }// in last case first we compute the inverse of the matrix and then we do matrix exponentiation 
          return resultant;
        }



    template< typename T,typename T2 >
        matrix<T> matrix_exponent(matrix<T> base,T2 exponent)    
       {
          matrix<T> resultant=base;
          exponent-=2;
          while(exponent>0)
            {
              if(exponent%2==1)resultant=prod(resultant,base);
              exponent=exponent >> 1;
              base=prod(base,base);
            }
          return resultant;  
       }
    }
  }
}
#endif

I am testing this header file using this 

#include "power.hpp"
using namespace boost::numeric::ublas;
using namespace std;
typedef permutation_matrix<std::size_t> pmatrix;
int main() 
{
matrix<double> m (3, 3);
for(int i=0;i<3;i++)for(int j=0;j<3;j++)m(i,j)=3*i+j+8;
m(0,0)=11;
matrix<double> c(3, 3);
int h=matrix_power(m,c,-1);// c stores -1 power of m
if(h)// h tells whether the power exists or not 
std:: cout << c << "\n\n\n";}

The functions work perfectly fine for power >0, because it use matrix exponentiation directly. The works much faster than repetitive multiplication, i have seen the for some 1000 iterations the running time of this and using a loop differ by a factor of 100-1000. You can observe that But for power <=0,i Get sometimes get incorrect answer. (I checked this using the idea that the product of a matrix and an inverse is identity matrix)

This probably has something to do with lu_factorize and lu_substitute which has certain checks to ensure variable types are correct.

Since their is no documentation available for lu_factorize, I am not sure how to use it. (I just read an example of it available here using Ublas lu_factorize and lu_substitute to compute matrix inverse ).I even read the source code but code not grasp much due to lack of experience.

I now have a few questions regarding the issue that i am experiencing. Since this is my first attempt at boost, if i ask something stupid please don't be hard on me. So these are my questions -

  1. The incorrect answer is probably due to incorrect conversion between data type,or something similar. How can i resolve this ? How do i ensure the usage of correct data types at every step ?
  2. How do i give errors when incorrect types are given as input by the user, I know that i could use boost assert but that is giving loads of incomprehensible compilation errors. what are the simplest ways to ensure that the input type is valid. for instance I want to give an error if the user gives a string for input. Can you provide an example for this?

  3. I have tried various methods to resolve the compilation errors one of which was to use #define BOOST_UBLAS_TYPE_CHECK 0 this helped bypass compilation errors for data types but then i got incorrect answers. can you explain how that works atleast in this case?

  4. Since this is my first attempt to make anything out of a boost, I can understand that this can certainly be made much better. What are the other things that must be present in this header file to ensure library standards like Error handling, Support for multiple compiler etc?

Solution

  • Your matrix_power() function does not return a value from the final else block. If I add return true at the end of that block then your code runs for me:

    $ clang++ --version
Apple LLVM version 7.0.0 (clang-700.1.76)
Target: x86_64-apple-darwin14.5.0
Thread model: posix

$ clang++ -std=c++11 -Wall -Wextra -g -I/opt/local/include lu.cpp -o lu
$ ./lu
[3,3]((-0.0644531,-0.00195312,0.861328),(1.09896,-0.0677083,-11.1406),(-0.9375,0.0625,9.4375))

    I was also able to remove that BOOST_UBLAS_TYPE_CHECK define without getting compile-time (or run-time) errors. I'm using Boost 1.59.