|
|
Ublas : |
Subject: [ublas] SVD decomposition and Moore penrose pseudo inverse
From: Dera_Augustin (deraaugustin_at_[hidden])
Date: 2014-05-23 05:57:32
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.