
Hello Joerg,
Hi Patrick,
without a symmetric type I can just check for a diagonal matrix inside the pinv function and branch, see code below.
That in fact is not true. Defining the interface for a moore-penrose-inverse function could look more like this (answering myself again ;-) ): diagonal_matrix<T> pinv(const diagonal_matrix<T>&); matrix<T> pinv(const matrix<T>&); but banded_matrix<T> pinv(const banded_matrix<T>&) is not ok, because the pinv of a banded_matrix<T> is a rectangular matrix<T>. So branching inside the function ist not nice. At all I think it is a nice example, that in a mathematical kind of view, there could be a huge difference between banded and diagonal.
As far as I've understood you're looking for a type safe diagonal variant of banded_matrix. Would it be helpful, if we'd go to add something like
it would be very nice to have something like that. But I think it is only useful if there are more applications with a semantic mathematical meaning. In my case the difference in the code is not huge: 1.) without diagonal_matrix: boost::numeric::ublas::matrix<double> pinv(const boost::numeric::ublas::matrix<double>& m) { using namespace boost::numeric::ublas; matrix<double> U,V; banded_matrix<double> S; svd(m,U,S,V); // assuming we have one for (size_t i=0;i<S.size1();i++) if (S(i,i)!=0) S(i,i)=1/S(i,i); // pseudoinverse of S stored in S return prod(prod(V,S),trans(U)); } 2.) with diagonal_matrix: boost::numeric::ublas::diagonal_matrix<double> pinv(boost::numeric::ublas::diagonal_matrix<double>& S) { for (size_t i=0;i<S.size1();i++) if (S(i,i)!=0) S(i,i)=1/S(i,i); return S; } boost::numeric::ublas::matrix<double> pinv(const boost::numeric::ublas::matrix<double>& m) { using namespace boost::numeric::ublas; matrix<double> U,V; banded_matrix<double> S; svd(m,U,S,V); return prod(prod(V,pinv(S)),trans(U)); } while in general I like more the second one, but the first one is more efficient :-). So, in this case we talk about one line moved out and afterwards moved in again to be more effective. OK, to be more enthusiastic: I think a diagonal_matrix would be very nice, it would be useful for me anyway.
// Diagonal matrix class template<class T, class F, class A> class diagonal_matrix: public banded_matrix<T, F, A> { public: BOOST_UBLAS_USING banded_matrix<T, F, A>::operator =; typedef banded_matrix<T, F, A> matrix_type;
// Construction and destruction BOOST_UBLAS_INLINE diagonal_matrix (): matrix_type () {} BOOST_UBLAS_INLINE diagonal_matrix (std::size_t size1, std::size_t size2): matrix_type (size1, size2) {} BOOST_UBLAS_INLINE ~diagonal_matrix () {}
// Assignment BOOST_UBLAS_INLINE diagonal_matrix &operator = (const diagonal_matrix &m) { matrix_type::operator = (m); return *this; } };
and the corresponding diagonal_adaptor to the code base?
YES !!!! JUHU !!! ...but I think
BOOST_UBLAS_INLINE diagonal_matrix (std::size_t size1, std::size_t size2): matrix_type (size1, size2) {}
could be: BOOST_UBLAS_INLINE diagonal_matrix (std::size_t size): matrix_type (size,size) {} or both.
Are there some recommendations for namespaces for functions like "pinv" ? Inside the loop, shall I use size_type?
I currently haven't enough time to look into a SVD implementation. So I'll cc this reply to http://groups.yahoo.com/group/ublas-dev.
The SVD is not important (i do not know what is inside either). Just keep in mind that a SVD for a matrix<T>A produces 3 matrices: matrix<T>U,diagonal_matrix<T>S, matrix<T>V. I should have pointed that out more clearly. Sorry. Aeh, is this thread OT here? Shall I just send these topics to the mailing list uBlas-dev? I will make a reply to the forward as well. Thank you, greetings Patrick