SVD decomposition and Moore penrose pseudo inverse

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

SVD decomposition and Moore penrose pseudo inverse

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<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!!!


Reply | Threaded
Open this post in threaded view
|

Re: SVD decomposition and Moore penrose pseudo inverse

Karl Meerbergen-2
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
[hidden email]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [hidden email]


Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm for more information.

_______________________________________________
ublas mailing list
[hidden email]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [hidden email]
Reply | Threaded
Open this post in threaded view
|

Re: SVD decomposition and Moore penrose pseudo inverse

Nasos Iliopoulos
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&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
[hidden email]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [hidden email]


Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm for more information.


_______________________________________________
ublas mailing list
[hidden email]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [hidden email]


_______________________________________________
ublas mailing list
[hidden email]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [hidden email]
Reply | Threaded
Open this post in threaded view
|

Re: SVD decomposition and Moore penrose pseudo inverse

Dera_Augustin
In reply to this post by Karl Meerbergen-2
Thx Karl, you right that solved the problem but now i got this one:

obj\Debug\main.o:D:\Tomography\svd\..\..\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\lapack\gesvd.hpp|81|undefined reference to `dgesvd_'|
||=== Build finished: 1 errors, 0 warnings ===|

I've been checking and i can confirm i have the file: gesvd.hpp in my binding\lapack directory

Why is it asking for dgesvd???

Thx Nasos, I will have a look on it but wouldn't give up on my idea
Reply | Threaded
Open this post in threaded view
|

Re: SVD decomposition and Moore penrose pseudo inverse

Oswin Krause
Hi,

you need to link against lapack

Regards,
Oswin

On 23.05.2014 14:36, Dera_Augustin wrote:

> Thx Karl, you right that solved the problem but now i got this one:
>
> obj\Debug\main.o:D:\Tomography\svd\..\..\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\lapack\gesvd.hpp|81|undefined
> reference to `dgesvd_'|
> ||=== Build finished: 1 errors, 0 warnings ===|
>
> I've been checking and i can confirm i have the file: gesvd.hpp in my
> binding\lapack directory
>
> Why is it asking for dgesvd???
>
> Thx Nasos, I will have a look on it but wouldn't give up on my idea
>
>
>
>
> --
> View this message in context: http://boost.2283326.n4.nabble.com/SVD-decomposition-and-Moore-penrose-pseudo-inverse-tp4662776p4662796.html
> Sent from the Boost - uBLAS mailing list archive at Nabble.com.
> _______________________________________________
> ublas mailing list
> [hidden email]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: [hidden email]

_______________________________________________
ublas mailing list
[hidden email]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [hidden email]
Reply | Threaded
Open this post in threaded view
|

Re: SVD decomposition and Moore penrose pseudo inverse

Dera_Augustin
Link against Lapack? How? soory but i'm new and i've been searching on the web but can't get through...
Thx in advance
Reply | Threaded
Open this post in threaded view
|

Re: SVD decomposition and Moore penrose pseudo inverse

Thomas Klimpel-2
Dera_Augustin wrote:
> Link against Lapack? How? soory but i'm new and i've been searching on the web but can't get through...

You seem to be on windows, so I would suggest you take "clapack-3.2.1-CMAKE.tgz" from "http://www.netlib.org/clapack/", download a recent version of CMake from "http://www.cmake.org/cmake/resources/software.html", find the correct directory with a CMakeLists.txt file, and use CMake to create a Visual Studio Solution. Then build that solution and copy the blas and lapack libraries to a convenient place. Then you add them as "Additional Dependencies" under "Linker" to your project.

There are other options like using Intel-MKL or AMD-AMCL libraries, but as you claim to be a newbee...

Regards,
Thomas
_______________________________________________
ublas mailing list
[hidden email]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [hidden email]
Reply | Threaded
Open this post in threaded view
|

Re: SVD decomposition and Moore penrose pseudo inverse

Dera_Augustin
In reply to this post by Dera_Augustin
Hi,
Finally, I got svd decomposition, but i have problem on making S as Matrix with makeMatrix function: I've already used that makeMatrix function with other vector but here with S, it give the error
D:\Tomography\svd\main.cpp|66|warning: "/*" within comment|
D:\Tomography\svd\main.cpp||In function 'int main()':|
D:\Tomography\svd\main.cpp|60|error: no matching function for call to 'makeMatrix(int, int, boost::numeric::ublas::vector<double, boost::numeric::ublas::unbounded_array<double, std::allocator<double> > >&)'|
||=== Build finished: 1 errors, 1 warnings ===|

My makeMatrix function is like this:
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);
}

I guess it has something to do with S vector and it's type but i can't find how to fix it
I have to get S as matrix to finally calculate Moore penrose pseudo inverse.

Thx

My entire code:
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <stdexcept>
#include <algorithm>
#include <vector>
#include "lapack.h"
//#include <D:\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\lapack\gesvd.hpp>
#include <boost/numeric/bindings/lapack/gesvd.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);
}

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, ublas::column_major> 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);
    gesvd(CtC, S, U, V);

    ublas::matrix<double> St=makeMatrix(18,1,S);
    std::cout<< "St= " << St << std::endl;
    /*ublas::matrix<double> pinv1=prod(S,Vt);
    ublas::matrix<double> Ut=trans(U);
    ublas::matrix<double> pinv=prod(pinv1,Ut);
    std::cout<< "pinv= " << pinv << 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;
}
Reply | Threaded
Open this post in threaded view
|

Re: SVD decomposition and Moore penrose pseudo inverse

Jörn Ungermann
Hi,

you supply a ublas::vector to a routine that expects a std::vector.

Without having the opportunity to verify that on my side, I would suppose that copying your makeMatrix routine and changing "std::vector" type to "ublas::vector" in the copy should do it the quick way.
Are you sure, that you want to convert S to a 18x1 matrix though?

Further, you should use matrix-vector products whenever possible, which seems feasible for your apparent use-case.
Also, using the SVD for the Moore-Penrose-Inverse makes only sense, in case if the problem is numerically unstable (ill-posed). In this case, you want to remove the "small" singular values from S before proceeding. I wrote this routine to invert an LES with an available svd decomposition:

    #include <cassert>

    /// Solves u * diag(s) * vt * x = b; cutoff determines the minimum value for a singular value to be considered.
    template<typename Majority>
    void svd_solve(const ublas::matrix<double, Majority>& u,
                   const ublas::matrix<double, Majority>& vt,
                   const ublas::vector<double>& s, const ublas::vector<double>& b, ublas::vector<double>& x,
                   double cutoff)
    {
      assert(cutoff >= 0);
      ublas::vector<double> temp(s.size());
      gemv(CblasTrans, 1.0, u, b, 0.0, temp);
      for (ublas::vector<double>::size_type i = 0; i < temp.size(); ++ i) {
        if (s(i) > cutoff) {
          temp(i) /= s(i);
        } else {
          temp(i) = 0;
        }
      }
      gemv(CblasTrans, 1.0, vt, temp, 0.0, x);
    }

"gemv" calls the atlas gemv in this case, but any other (also ublas-supplied) matrix vector product will do here.

Cheers,
Jörn


> -----Original Message-----
> From: ublas [mailto:[hidden email]] On Behalf Of
> Dera_Augustin
> Sent: Dienstag, 27. Mai 2014 10:25
> To: [hidden email]
> Subject: Re: [ublas] SVD decomposition and Moore penrose pseudo inverse
>
> Hi,
> Finally, I got svd decomposition, but i have problem on making S as
> Matrix with makeMatrix function: I've already used that makeMatrix
> function with other vector but here with S, it give the error
> D:\Tomography\svd\main.cpp|66|warning: "/*" within comment|
> D:\Tomography\svd\main.cpp||In function 'int main()':|
> D:\Tomography\svd\main.cpp|60|error: no matching function for call to
> 'makeMatrix(int, int, boost::numeric::ublas::vector<double,
> boost::numeric::ublas::unbounded_array&lt;double,
> std::allocator&lt;double>
> > >&)'|
> ||=== Build finished: 1 errors, 1 warnings ===|
>
> My makeMatrix function is like this:
> 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); }
>
> I guess it has something to do with S vector and it's type but i can't
> find how to fix it I have to get S as matrix to finally calculate Moore
> penrose pseudo inverse.
>
> Thx
>
> My entire code:
> #include <iostream>
> #include <stdio.h>
> #include <stdlib.h>
> #include <stdexcept>
> #include <algorithm>
> #include <vector>
> #include "lapack.h"
> //#include
> <D:\boost_1_47_0\boost_1_47_0\boost\numeric\bindings\lapack\gesvd.hpp>
> #include <boost/numeric/bindings/lapack/gesvd.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); }
>
> 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, ublas::column_major> 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);
>     gesvd(CtC, S, U, V);
>
>     ublas::matrix<double> St=makeMatrix(18,1,S);
>     std::cout<< "St= " << St << std::endl;
>     /*ublas::matrix<double> pinv1=prod(S,Vt);
>     ublas::matrix<double> Ut=trans(U);
>     ublas::matrix<double> pinv=prod(pinv1,Ut);
>     std::cout<< "pinv= " << pinv << 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;
> }
>
>
>
> --
> View this message in context: http://boost.2283326.n4.nabble.com/SVD-
> decomposition-and-Moore-penrose-pseudo-inverse-tp4662776p4662921.html
> Sent from the Boost - uBLAS mailing list archive at Nabble.com.
> _______________________________________________
> ublas mailing list
> [hidden email]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: [hidden email]


------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------
Forschungszentrum Juelich GmbH
52425 Juelich
Sitz der Gesellschaft: Juelich
Eingetragen im Handelsregister des Amtsgerichts Dueren Nr. HR B 3498
Vorsitzender des Aufsichtsrats: MinDir Dr. Karl Eugen Huthmacher
Geschaeftsfuehrung: Prof. Dr. Achim Bachem (Vorsitzender),
Karsten Beneke (stellv. Vorsitzender), Prof. Dr.-Ing. Harald Bolt,
Prof. Dr. Sebastian M. Schmidt
------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------

_______________________________________________
ublas mailing list
[hidden email]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [hidden email]