|
|
Ublas : |
Subject: Re: [ublas] SVD decomposition and Moore penrose pseudo inverse
From: Nasos Iliopoulos (nasos_i_at_[hidden])
Date: 2014-05-23 07:28:52
Or you can use Marco's implementation from here:
https://github.com/sguazt/boost-ublasx/blob/master/boost/numeric/ublasx/operation/svd.hpp
-Nasos
On 05/23/2014 06:00 AM, Karl Meerbergen wrote:
> 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<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<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<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<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 for more
> information.
>
>
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: athanasios.iliopoulos.ctr.gr_at_[hidden]