|
|
Ublas : |
Subject: Re: [ublas] SVD decomposition and Moore penrose pseudo inverse
From: Ungermann, Jörn (j.ungermann_at_[hidden])
Date: 2014-05-27 07:22:58
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:ublas-bounces_at_[hidden]] On Behalf Of
> Dera_Augustin
> Sent: Dienstag, 27. Mai 2014 10:25
> To: ublas_at_[hidden]
> 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<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;
> }
>
>
>
> --
> 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
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: j.ungermann_at_[hidden]
------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------
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
------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------