Re: [Boost-bugs] [Boost C++ Libraries] #11295: Matrix Memory problem after using lu_factorize

Subject: Re: [Boost-bugs] [Boost C++ Libraries] #11295: Matrix Memory problem after using lu_factorize
From: Boost C++ Libraries (noreply_at_[hidden])
Date: 2017-03-31 10:45:26

#11295: Matrix Memory problem after using lu_factorize
-------------------------------+---------------------
  Reporter:  michael.cortis@…  |      Owner:  guwi17
      Type:  Bugs              |     Status:  new
 Milestone:  To Be Determined  |  Component:  uBLAS
   Version:  Boost 1.55.0      |   Severity:  Problem
Resolution:                    |   Keywords:
-------------------------------+---------------------

Comment (by 2015csb1032@…):

 You interpreted the function in the wrong way, that's what the
 lu_factorize function does.
 The return value, 0 or non 0 value which you see, only specifies whether
 the matrix is singular or not.
 The return value is 0 if the matix is singular, and non 0 if the matrix is
 non-singular.

 Now intuitively, you would think that lu_factorize should give two
 matrices L and U on decomposition.
 That's what is done in a way.

 The matrix which you pass in lu_factorize(X), here X is passed by
 reference to lu_factorize function,which will '''decompose X and now, X
 will contain the information of the two factorized matrix L and U. Hence X
 is changed.'''[[BR]]
 If you don't want X to be changed. Then, create another matrix Y, do Y=X
 and then pass Y in lu_factorize : lu_factorize(Y)

 It's a typical approach in LU decomposition, that the diagonal elements of
 L contains 1. So, the two matrices
 L and U are fitted into 1 matrix Y, putting L in the lower part omitting
 the diagonal elements, and
 putting U in the upper part. L and U can be extracted from Y easily (shown
 in code).

 If, Y (after function ends) = [[BR]]
 y1 y2 y3[[BR]]
 y4 y5 y6[[BR]]
 y7 y8 y9[[BR]]

 Then, L=[[BR]]
 1 0 0[[BR]]
 y4 1 0[[BR]]
 y7 y7 1[[BR]]

 and u =[[BR]]
 y1 y2 y3[[BR]]
 0 y5 y6[[BR]]
 0 0 y9[[BR]]


 You can run the following code, it will show what I'm trying to say above.
 {{{
   ublas::matrix<double> X(3,3); X.clear();
   X(0,0) = 0.995182407377577; X(0,1) =-0.006473367705848; X(0,2)
 =-0.002032391957706;
   X(1,0) =-0.006473367705848; X(1,1) = 0.995182407377577; X(1,2)
 =-0.002032391957706;
   X(2,0) =-0.002032391957706; X(2,1) =-0.002032391957706; X(2,2) =
 0.936175146339137;

   ublas::matrix<double> Y(3,3);
   Y=X;
   ublas::lu_factorize(Y);

   ublas::matrix<double> L(3,3); L.clear();
   ublas::matrix<double> U(3,3); U.clear();

   for (int i=0; i<3; i++){
         for (int j=0; j<3; j++){
                 if (i<j){
                         U(i,j) = Y(i,j);
                         L(i,j) = 0;
                 }
                 else if (i>j){
                         L(i,j) = Y(i,j);
                         U(i,j) = 0;
                 }
                 else if (i==j){
                         L(i,j) = 1;
                         U(i,j) = Y(i,j);
                 }
         }
   }

   cout << "The original matrix X :"<<endl;
   cout<<setprecision(16)<<X<<endl<<endl;

   cout << "LU_decomposed matrix in 1 matrix :"<<endl;
   cout<<setprecision(16)<<Y<<endl<<endl;

   cout << "LU_decomposed matrices in 2 diff matrix, L :"<<endl;
   cout<<setprecision(16)<<L<<endl<<endl;

   cout << "LU_decomposed matrices in 2 diff matrix, U :"<<endl;
   cout<<setprecision(16)<<U<<endl<<endl;

   ublas::matrix<double> Z(3,3);
   cout << "The product Z=L*U, which should be equal to X (Z=L*U=X),
 Z:"<<endl;
   axpy_prod(L, U, Z, true);
   cout<<setprecision(16)<<Z<<endl<<endl;
 }}}

 Also, run the code with
 {{{
   X(0,0) = 1.0; X(0,1) = 2.0; X(0,2) = 3.0;
   X(1,0) = 4.0; X(1,1) = 5.0; X(1,2) = 6.0;
   X(2,0) = 7.0; X(2,1) = 8.0; X(2,2) = 9.0;
 }}}
  instead of X(0,0) = 0.995182407377577;... for better understanding.


 Replying to [ticket:11295 michael.cortis@…]:
 > Discovered that after using lu_factorize the values of the matrix have
 changed!
 >
 >
 > {{{
 > #include <iostream>
 >
 > #include <boost/date_time/posix_time/posix_time.hpp>
 > #include <boost/numeric/ublas/matrix.hpp>
 > #include <boost/numeric/ublas/matrix_proxy.hpp>
 > #include <boost/numeric/ublas/vector.hpp>
 > #include <boost/numeric/ublas/io.hpp>
 > #include <boost/numeric/ublas/lu.hpp>
 >
 > using namespace std;
 > using namespace boost::numeric;
 >
 > int main()
 > {
 >   ublas::matrix<double> X(3,3); X.clear();
 >   X(0,0) = 0.995182407377577; X(0,1) =-0.006473367705848; X(0,2)
 =-0.002032391957706;
 >   X(1,0) =-0.006473367705848; X(1,1) = 0.995182407377577; X(1,2)
 =-0.002032391957706;
 >   X(2,0) =-0.002032391957706; X(2,1) =-0.002032391957706; X(2,2) =
 0.936175146339137;
 >   cout<<setprecision(16)<<X<<endl;
 >   cout<<ublas::lu_factorize(X)<<endl;
 >   cout<<setprecision(16)<<X<<endl;
 >   cout<<ublas::lu_factorize(X)<<endl;
 >   cout<<setprecision(16)<<X<<endl;
 > }
 > }}}
 >
 > Output:
 >
 >
 > {{{
 >
 [3,3]((0.995182407377577,-0.006473367705848,-0.002032391957706),(-0.006473367705848,0.995182407377577,-0.002032391957706),(-0.002032391957706,-0.002032391957706,0.936175146339137))
 > 0
 >
 [3,3]((0.995182407377577,-0.006473367705848,-0.002032391957706),(-0.006504704723334175,0.9951403000320849,-0.002045612067272957),(-0.002042230592742885,-0.002055601674665375,0.9361667907625133))
 > 0
 >
 [3,3]((0.995182407377577,-0.006473367705848,-0.002032391957706),(-0.006536193440632496,0.9950979888485471,-0.002058896174255709),(-0.002052116855767581,-0.002079077442455779,0.9361583394521271))
 >
 > }}}
 >
 > Is this fixed in newer versions?
 >
 > Thanks
 >
 > Michael Cortis
 >
 > RA, Durham University

-- 
Ticket URL: <https://svn.boost.org/trac/boost/ticket/11295#comment:1>
Boost C++ Libraries <http://www.boost.org/>
Boost provides free peer-reviewed portable C++ source libraries.