Boost logo

Ublas :

Subject: Re: [ublas] SVD decomposition and Moore penrose pseudo inverse
From: Karl Meerbergen (karl.meerbergen_at_[hidden])
Date: 2014-05-23 06:00:36


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
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: karl.meerbergen_at_[hidden]

Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm