SVD decomposition and moore penrose pseudo inverse of matrix

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

SVD decomposition and moore penrose pseudo inverse of matrix

Dera_Augustin
Hi all, First i'm newbie with boost::ublas and ublas binding I'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 from http://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)); #endif And 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 template bool InvertMatrix(const ublas::matrix& input, ublas::matrix& inverse) { using namespace boost::numeric::ublas; typedef permutation_matrix pmatrix; // create a working copy of the input matrix A(input); // create a permutation matrix for the LU-factorization pmatrix pm(A.size1()); // perform LU-factorization int 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 inverse lu_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!!!