Boost C++ Libraries: Ticket #11295: Matrix Memory problem after using lu_factorize

2015csb1032@… — Fri, 31 Mar 2017 10:45:26 GMT

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.
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) =
y1 y2 y3
y4 y5 y6
y7 y8 y9

Then, L=
1 0 0
y4 1 0
y7 y7 1

and u =
y1 y2 y3
0 y5 y6
0 0 y9

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 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

anonymous — Fri, 16 Jun 2017 11:49:54 GMT

Thanks for the explanation. Didn't realize that both L and U are combined into a single matrix.

Thanks again, Michael

Replying to 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.
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) =
y1 y2 y3
y4 y5 y6
y7 y8 y9

Then, L=
1 0 0
y4 1 0
y7 y7 1

and u =
y1 y2 y3
0 y5 y6
0 0 y9

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 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