|
|
Ublas : |
Subject: [ublas] SVD decomposition and moore penrose pseudo inverse of matrix
From: Dera_Augustin (deraaugustin_at_[hidden])
Date: 2014-05-23 05:51:11
Hi all,First i'm newbie with boost::ublas and ublas bindingI'm trying to
solve an orthogonal projection problem.I've been trying to get the pseudo
inverse of my matrix with numerical recipe of lu decomposition and back
substitute that i found online but the fact that i have no clou about how to
fix the tolerance or precision i mean the digits after the decimal point so
i leave it to focus on Moore penrose pseudo inverse.If someone there have an
idea?My problem now also is to get svd decompostion, i've written the code
below and i get error of static assertion failed with the binding
:#include<boost/numeric/bindings/traits/ublas_matrix.hpp> that i've got
fromhttp://boost.cvs.sourceforge.net/viewvc/boost-sandbox/boost-sandbox/boost/numeric/bindings/lapack/The
exact error is on line 49 of the file ublas_matrix.hpp#ifdef
BOOST_NUMERIC_BINDINGS_FORTRAN
BOOST_STATIC_ASSERT((boost::is_same::value));#endifAnd the build message
error i got
is..\..\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\traits\matrix_traits.hpp|51|instantiated
from
'boost::numeric::bindings::traits::matrix_traits<boost::numeric::ublas::matrix<double,
boost::numeric::ublas::basic_row_major<unsigned int, int>,
boost::numeric::ublas::unbounded_array<double, std::allocator<double> > >
>'|..\..\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\lapack\gesdd.hpp|632|instantiated
from 'int boost::numeric::bindings::lapack::gesdd(char, char, MatrA&, VecS&,
MatrU&, MatrV&) [with MatrA = boost::numeric::ublas::matrix<double,
boost::numeric::ublas::basic_row_major<unsigned int, int>,
boost::numeric::ublas::unbounded_array<double, std::allocator<double> >
>, VecS = boost::numeric::ublas::vector<double,
boost::numeric::ublas::unbounded_array<double, std::allocator<double>
> >, MatrU = boost::numeric::ublas::matrix<double,
boost::numeric::ublas::basic_column_major<unsigned
in|..\..\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\lapack\gesdd.hpp|647|instantiated
from 'int boost::numeric::bindings::lapack::gesdd(MatrA&, VecS&,
MatrU&, MatrV&) [with MatrA =
boost::numeric::ublas::matrix<double,
boost::numeric::ublas::basic_row_major<unsigned int, int>,
boost::numeric::ublas::unbounded_array<double, std::allocator<double> >
>, VecS = boost::numeric::ublas::vector<double,
boost::numeric::ublas::unbounded_array<double, std::allocator<double>
> >, MatrU = boost::numeric::ublas::matrix<double,
boost::numeric::ublas::basic_column_major<unsigned int, int>,
boo|D:\Tomography\svd\main.cpp|81|instantiated from
here|..\..\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\traits\ublas_matrix.hpp|49|error:
static assertion failed:
"(boost::is_same::value)"|..\..\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\traits\matrix_traits.hpp|51|instantiated
from 'boost::numeric::bindings::traits::matrix_traits,
boost::numeric::ublas::unbounded_array<double, std::allocator<double> > >
>'|..\..\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\traits\matrix_traits.hpp|102|instantiated
from 'int boost::numeric::bindings::traits::matrix_size1(M&) [with M = const
boost::numeric::ublas::matrix<double,
boost::numeric::ublas::basic_row_major<unsigned int, int>,
boost::numeric::ublas::unbounded_array<double, std::allocator<double> >
>]'|..\..\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\lapack\gesdd.hpp|278|instantiated
from 'int boost::numeric::bindings::lapack::gesdd_work(char, char, const
MatrA&) [with MatrA = boost::numeric::ublas::matrix<double,
boost::numeric::ublas::basic_row_major<unsigned int, int>,
boost::numeric::ublas::unbounded_array<double, std::allocator<double> >
>]'|..\..\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\lapack\gesdd.hpp|606|instantiated
from 'int boost::numeric::bindings::lapack::gesdd(char, char, MatrA&, VecS&,
MatrU&, MatrV&) [with MatrA = boost::numeric::ublas::matrix<double,
boost::numeric::ublas::basic_row_major<unsigned int, int>,
boost::numeric::ublas::unbounded_array<double, std::allocator<double> >
>, VecS = boost::numeric::ublas::vector<double,
boost::numeric::ublas::unbounded_array<double, std::allocator<double>
> >, MatrU = boost::numeric::ublas::matrix<double,
boost::numeric::ublas::basic_column_major<unsigned
in|..\..\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\lapack\gesdd.hpp|647|instantiated
from 'int boost::numeric::bindings::lapack::gesdd(MatrA&, VecS&,
MatrU&, MatrV&) [with MatrA =
boost::numeric::ublas::matrix<double,
boost::numeric::ublas::basic_row_major<unsigned int, int>,
boost::numeric::ublas::unbounded_array<double, std::allocator<double> >
>, VecS = boost::numeric::ublas::vector<double,
boost::numeric::ublas::unbounded_array<double, std::allocator<double>
> >, MatrU = boost::numeric::ublas::matrix<double,
boost::numeric::ublas::basic_column_major<unsigned int, int>,
boo|D:\Tomography\svd\main.cpp|81|instantiated from
here|..\..\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\traits\ublas_matrix.hpp|49|error:
static assertion failed: "(boost::is_same::value)"|||=== Build finished: 2
errors, 0 warnings ===|My code is like this:#include #include #include
#include #include #include //#include
<D:\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\lapack\gesvd.hpp>#include
<boost/numeric/bindings/lapack/gesdd.hpp>#include
<D:\boost_1_47_0\boost_1_47_0\boost\numeric\ublas\storage.hpp>#include
<D:\boost_1_47_0\boost_1_47_0\boost\cast.hpp>#include
<D:\boost_1_47_0\boost_1_47_0\boost\numeric\ublas\vector.hpp>#include
<D:\boost_1_47_0\boost_1_47_0\boost\numeric\ublas\vector_proxy.hpp>#include
<D:\boost_1_47_0\boost_1_47_0\boost\numeric\ublas\matrix.hpp>#include
<D:\boost_1_47_0\boost_1_47_0\boost\numeric\ublas\triangular.hpp>#include
<D:\boost_1_47_0\boost_1_47_0\boost\numeric\ublas\lu.hpp>#include
<D:\boost_1_47_0\boost_1_47_0\boost\numeric\ublas\io.hpp>#include<boost/numeric/bindings/traits/ublas_matrix.hpp>#include<boost/numeric/bindings/traits/ublas_vector.hpp>using
namespace std;namespace ublas=boost::numeric::ublas;using namespace
boost::numeric::bindings::lapack;template ublas::matrix
makeMatrix(std::size_t m, std::size_t n, std::vector & v){
if(m*n!=v.size()) { ; // Handle this case }
ublas::unbounded_array storage(m*n); std::copy(v.begin(), v.end(),
storage.begin()); return boost::numeric::ublas::matrix(m, n,
storage);}/*LU PSEUDO INVERSE CANT FIX THE PRECISION templatebool
InvertMatrix(const ublas::matrix& input, ublas::matrix& inverse) {using
namespace boost::numeric::ublas;typedef permutation_matrix pmatrix;// create
a working copy of the inputmatrix A(input);// create a permutation matrix
for the LU-factorizationpmatrix pm(A.size1());// perform LU-factorizationint
res = lu_factorize(A,pm);if( res != 0 ) return false;// create identity
matrix of "inverse"inverse.assign(ublas::identity_matrix(A.size1()));//
backsubstitute to get the inverselu_substitute(A, pm, inverse);return
true;}*/int main(){ cout << "Hello world!" << endl; std:: vector
vecDgk{1, 6.5, 0.5, 3.636, 5.91, 0.454, 2.91, 0.727, 0.363, 0.696, 0.174,
1.13}; std:: vector vecbi{3, 3, 8, 8, 5, 14, 10, 2, 1, 2, 8, 8, 10, 2, 4,
13, 2, 1}; std:: vector vecCt
{0,0.5,0,0,0,0,0,0,0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,0,0,0,1,0,0,1,1,1,1,0,0,1,0.5,1,0,0,0,0,1,1,0,0,0,0,0.5,0,1,1,0,0,0,0,1,0,1,0,0,0,1,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,1,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,1,1,1,0,0,0,0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,1,1,0,0,1,0,0,1,0,0,1,0,0,0};
//{{0,0.5,0,0,0,0,0,0,0,0,0,0}, {0,0.5,0,0,0,0,0,0,0,0,0,0},
{0,0,0,1,0,0,1,1,1,1,0,0},
{0,0,0,1,0,0,1,1,1,1,0,0},{1,0.5,1,0,0,0,0,1,1,0,0,0},
{0,0.5,0,1,1,0,0,0,0,1,0,1}, {0,0,0,1,1,1,0,0,0,0,0,0},
{1,0,0,0,0,0,0,1,0,0,1,0}, {0,0,1,0,0,1,0,0,1,0,0,0},
{0,0,0,0,0,0,0,0,0,1,1,1}, {1,1,1,0,0,0,0,0,0,0,0,0},
{0,0,0,1,0,0,1,0,0,1,0,0}, {0,0,0,1,1,1,0,0,0,0,0,0},
{1,0,0,0,0,0,0,1,0,0,1,0}, {0,0,0,0,0,0,1,1,1,0,0,0},
{0,1,0,0,1,0,0,0,0,0,0,1}, {0,0,0,0,0,0,0,0,0,1,1,1},
{0,0,1,0,0,1,0,0,1,0,0,0}}; ublas::matrix Ct=makeMatrix(18,12,vecCt);
//std::cout<< Ct << std::endl; ublas::matrix bi=makeMatrix(18,1,vecbi);
//std::cout<< bi << std::endl; ublas::matrix Dgk=makeMatrix(12,1,vecDgk);
//std::cout<< Dgk << std::endl; ublas::matrix C=trans(Ct);
ublas::matrix CtC=prod(Ct,C); /*ublas::matrix invCtC;
svd_moore_penrose_pseudoinverse(matrix_type const & CtC, matrix_type&
invCtC);*/ std::cout<< "CtC= " << CtC << std::endl;
ublas::matrix<double, ublas::column_major> U(18,18);
ublas::matrix<double, ublas::column_major> V(18,18); ublas::vector S(18);
gesdd(CtC, S, U, V); std::cout << U << std::endl; std::cout << S <<
std::endl; std::cout << V << std::endl; /*ublas::matrix pinv (18,18);
InvertMatrix(CtC,pinv); //std::cout<< pinv << std::endl; ublas::matrix
Cpinv=prod(C,pinv); ublas::matrix correction=prod(Cpinv,bi);
std::cout<< "Correction: " << correction << std::endl; ublas::matrix
CpinvCt=prod(Cpinv,Ct); ublas::matrix projection=prod(CpinvCt,Dgk);
std::cout<< "Projection: " << projection << std::endl; ublas::matrix
dEstimees; dEstimees=correction-projection + Dgk; std::cout<<
"dEstimees: " << dEstimees << std::endl;*/ return 0;}Thx for any help you
can provide!!!
-- View this message in context: http://boost.2283326.n4.nabble.com/SVD-decomposition-and-moore-penrose-pseudo-inverse-of-matrix-tp4662775.html Sent from the Boost - uBLAS mailing list archive at Nabble.com.