Re: SVD decomposition and Moore penrose pseudo inverse

Karl Meerbergen <[email protected]>
Newsgroups gmane.comp.lib.boost.ublas
Message-ID <[email protected]> 
Hi Dera,

You should use a column_major matrix in order to use LAPACK routines. By 
default matrices are stored row_major.

Best regards,

Karl



On 05/23/2014 11:57 AM, Dera_Augustin wrote:
> 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<typename F::orientation_category,
> boost::numeric::ublas::column_major_tag>::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&lt;double,
> boost::numeric::ublas::basic_row_major&lt;unsigned int, int>,
> boost::numeric::ublas::unbounded_array<double, std::allocator&lt;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&lt;unsigned int, int>,
> boost::numeric::ublas::unbounded_array<double, std::allocator&lt;double> >
>> , VecS = boost::numeric::ublas::vector<double,
> boost::numeric::ublas::unbounded_array&lt;double, std::allocator&lt;double>
>>> , MatrU = boost::numeric::ublas::matrix<double,
> boost::numeric::ublas::basic_column_major&lt;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&amp;, VecS&amp;,
> MatrU&amp;, MatrV&amp;) [with MatrA =
> boost::numeric::ublas::matrix&lt;double,
> boost::numeric::ublas::basic_row_major&lt;unsigned int, int>,
> boost::numeric::ublas::unbounded_array<double, std::allocator&lt;double> >
>> , VecS = boost::numeric::ublas::vector<double,
> boost::numeric::ublas::unbounded_array&lt;double, std::allocator&lt;double>
>>> , MatrU = boost::numeric::ublas::matrix<double,
> boost::numeric::ublas::basic_column_major&lt;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<typename F::orientation_category,
> boost::numeric::ublas::column_major_tag>::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<const
> boost::numeric::ublas::matrix<double,
> boost::numeric::ublas::basic_row_major<unsigned int, int>,
> boost::numeric::ublas::unbounded_array<double, std::allocator&lt;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&lt;unsigned int, int>,
> boost::numeric::ublas::unbounded_array<double, std::allocator&lt;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&lt;unsigned int, int>,
> boost::numeric::ublas::unbounded_array<double, std::allocator&lt;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&lt;unsigned int, int>,
> boost::numeric::ublas::unbounded_array<double, std::allocator&lt;double> >
>> , VecS = boost::numeric::ublas::vector<double,
> boost::numeric::ublas::unbounded_array&lt;double, std::allocator&lt;double>
>>> , MatrU = boost::numeric::ublas::matrix<double,
> boost::numeric::ublas::basic_column_major&lt;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&amp;, VecS&amp;,
> MatrU&amp;, MatrV&amp;) [with MatrA =
> boost::numeric::ublas::matrix&lt;double,
> boost::numeric::ublas::basic_row_major&lt;unsigned int, int>,
> boost::numeric::ublas::unbounded_array<double, std::allocator&lt;double> >
>> , VecS = boost::numeric::ublas::vector<double,
> boost::numeric::ublas::unbounded_array&lt;double, std::allocator&lt;double>
>>> , MatrU = boost::numeric::ublas::matrix<double,
> boost::numeric::ublas::basic_column_major&lt;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<typename F::orientation_category,
> boost::numeric::ublas::column_major_tag>::value)"|
> ||=== Build finished: 2 errors, 0 warnings ===|
>
>
>
>
>
> My code is like this:
>
> #include <iostream>
> #include <stdio.h>
> #include <stdlib.h>
> #include <stdexcept>
> #include <algorithm>
> #include <vector>
> //#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 <typename T>
> ublas::matrix<T> makeMatrix(std::size_t m, std::size_t n, std::vector<T> &
> v)
> {
>      if(m*n!=v.size()) {
>          ; // Handle this case
>      }
>      ublas::unbounded_array<T> storage(m*n);
>      std::copy(v.begin(), v.end(), storage.begin());
>      return boost::numeric::ublas::matrix<T>(m, n, storage);
> }
> /*LU PSEUDO INVERSE CANT FIX THE PRECISION
>   template<class T>
> bool InvertMatrix(const ublas::matrix<T>& input, ublas::matrix<T>& inverse)
> {
> using namespace boost::numeric::ublas;
> typedef permutation_matrix<std::size_t> pmatrix;
> // create a working copy of the input
> matrix<T> 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<T>(A.size1()));
>
> // backsubstitute to get the inverse
> lu_substitute(A, pm, inverse);
>
> return true;
> }
> */
>
>
> int main()
> {
>      cout << "Hello world!" << endl;
>      std:: vector<double> 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<double> vecbi{3, 3, 8, 8, 5, 14, 10, 2, 1, 2, 8, 8, 10, 2,
> 4, 13, 2, 1};
>      std:: vector<double> 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<double> Ct=makeMatrix(18,12,vecCt);
>      //std::cout<< Ct << std::endl;
>      ublas::matrix<double> bi=makeMatrix(18,1,vecbi);
>      //std::cout<< bi << std::endl;
>      ublas::matrix<double> Dgk=makeMatrix(12,1,vecDgk);
>      //std::cout<< Dgk << std::endl;
>      ublas::matrix<double> C=trans(Ct);
>      ublas::matrix<double> CtC=prod(Ct,C);
>      /*ublas::matrix<double> 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<double> S(18);
>      gesdd(CtC, S, U, V);
>
>      std::cout << U << std::endl;
>      std::cout << S << std::endl;
>      std::cout << V << std::endl;
>      /*ublas::matrix<double> pinv (18,18);
>      InvertMatrix(CtC,pinv);
>      //std::cout<< pinv << std::endl;
>      ublas::matrix<double> Cpinv=prod(C,pinv);
>      ublas::matrix<double> correction=prod(Cpinv,bi);
>      std::cout<< "Correction: " << correction << std::endl;
>      ublas::matrix<double> CpinvCt=prod(Cpinv,Ct);
>      ublas::matrix<double> projection=prod(CpinvCt,Dgk);
>      std::cout<< "Projection: " << projection << std::endl;
>      ublas::matrix<double> 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-tp4662776.html
> Sent from the Boost - uBLAS mailing list archive at Nabble.com.
> _______________________________________________
> ublas mailing list
> [email protected]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: [email protected]


Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
lmpx.com only provides a reader for public news (NNTP) servers. It is not affiliated with the servers or forums shown here and is not responsible for the content of articles, which is written by their respective authors.