|
|
Ublas : |
Subject: Re: [ublas] boost::ublas lu_substitute runtime error
From: Ungermann, Jörn (j.ungermann_at_[hidden])
Date: 2011-08-17 02:59:44
Hi Andrey,
I guess the problem stems from the ill-conditionedness of the matrix you try
to invert.
The assertion complains about one step in the computation of the inverse not
being fully reversible to machine precision and is perhaps a bit
overzealous.
But the creators of ublas wanted foremost correct solutions, not fast ones,
so this is expected.
The following code replicates the problem on my end:
<<< snippet start ----
//#define BOOST_UBLAS_NDEBUG 1
#include <iostream>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include "storage_adaptors.hpp"
using namespace boost::numeric::ublas;
using namespace std;
template <class T>
static bool invertMatrix(matrix<T>& orig, matrix<T>& inverted) {
cout << "Orig Matrix size: " << orig.size1() << " " << orig.size2() <<
endl;
cout << "Inv Matrix size: " << inverted.size1() << " " <<
inverted.size2() << endl;
typedef permutation_matrix<std::size_t> pmatrix;
matrix<T> A(orig);
pmatrix pm(A.size1());
// perform LU-factorization
int res = lu_factorize(A,pm);
if( res != 0 ) return false;
cout << "calling assign, size a: " << A.size1() << " " <<
inverted.size1() << endl;;
inverted.assign(identity_matrix<T>(A.size1()));
cout << "A Matrix size: " << A.size1() << " " << A.size2() << endl;
cout << "PM Matrix size: " << pm.size() << endl;
lu_substitute(A, pm, inverted);
return true;
}
int main() {
double data[6][6] = {{15, 29700, 472042, 7.8021e+06, 1.32426e+08,
2.29091e+09},
{29700,1.32426e+08,2.29091e+09,4.01989e+10,7.13142e+11,1.27611e+13},
{472042,2.29091e+09,4.01989e+10,7.13142e+11,1.27611e+13,2.29941e+14},
{7.8021e+06,4.01989e+10,7.13142e+11,1.27611e+13,2.29941e+14,4.16694e+15},
{1.32426e+08,7.13142e+11,1.27611e+13,2.29941e+14,4.16694e+15,7.58705e+16},
{2.29091e+09,1.27611e+13,2.29941e+14,4.16694e+15,7.58705e+16,1.38694e+18}};
matrix<double> a = make_matrix_from_pointer(data);
matrix<double> b(a);
invertMatrix(a, b);
}
snippet stop ---->>>
using G. Winklers storage_adaptors.hpp
(http://www.guwi17.de/ublas/examples/storage_adaptors.hpp).
If you know what you are doing mathematically, you may simply comment in the
first line,
which disables this check and makes UBLAS run *much* faster, so it is a good
idea anyway.
You could also twiddle with the BOOST_UBLAS_TYPE_CHECK_EPS and
BOOST_UBLAS_TYPE_CHECK_MIN defines
To make the check more fore-giving, but choosing proper values here for
arbitrarily ill-conditioned
Matrixes sounds difficult.
The inverse produced by LU looks quite OK, in my eyes and is identical to
that produced by ATLAS
And python/scipy, so you should be fine. You might want to check from time
to time, that the product
of matrix and inverse does not look to different from the identity matrix,
though.
However, a condition number of 1.6e+18 for a 6x6 matrix is quite abysmal, so
you might want to look
into pseudo-inverses you can compute, e.g., via singular value
decomposition.
We use the ATLAS backend and numeric bindings to that purpose.
Regards,
Joern
> -----Original Message-----
> From: ublas-bounces_at_[hidden] [mailto:ublas-
> bounces_at_[hidden]] On Behalf Of mga-mga
> Sent: Dienstag, 16. August 2011 16:02
> To: ublas_at_[hidden]
> Subject: Re: [ublas] boost::ublas lu_substitute runtime error
>
> Hello!
>
> Thank you so much for your quick answer! Actually, there are printouts
> showing that both original matrix and the matrix holding the inverted
> version are both of the same size. Here is the procedure I used for
> inversion, with debug printouts (which are redundant, please do not pay
> attention) as appropriate:
>
> template <class T>
> static bool invertMatrix(matrix<T>& orig, matrix<T>& inverted)
> {
> cout << "Orig Matrix size: " << orig.size1() << " " << orig.size2()
> <<
> endl;
> cout << "Inv Matrix size: " << inverted.size1() << " " <<
> inverted.size2() << endl;
>
> typedef permutation_matrix<std::size_t> pmatrix;
> matrix<T> A(orig);
> pmatrix pm(A.size1());
>
> // perform LU-factorization
> int res = lu_factorize(A,pm);
> if( res != 0 ) return false;
>
> cout << "calling assign, size a: " << A.size1() << " " <<
> inverted.size1() << endl;;
> inverted.assign(identity_matrix<T>(A.size1()));
>
> cout << "A Matrix size: " << A.size1() << " " << A.size2() << endl;
> cout << "PM Matrix size: " << pm.size() << endl;
> lu_substitute(A, pm, inverted);
>
> return true;
> }
>
>
> And here are the printouts:
>
> Orig Matrix size: 7 7
> Inv Matrix size: 7 7
> calling assign, size a: 7 7
> A Matrix size: 7 7
> PM Matrix size: 7
> Check failed in file /usr/local/include/boost/numeric/ublas/lu.hpp at
> line
> 294:
> detail::expression_type_check (prod
> (triangular_adaptor<const_matrix_type, upper> (m), e), cm2)
>
> Sorry, I forgot to mention: I'm using RedHat + boost 1.42.0
>
> WBR, Andrey
>
> --
> View this message in context: http://boost.2283326.n4.nabble.com/boost-
> ublas-lu-substitute-runtime-error-tp3747071p3747294.html
> Sent from the Boost - uBLAS mailing list archive at Nabble.com.
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: j.ungermann_at_[hidden]