Re: Re: uBlas - diagonal_matrix

From: Patrick Kowalzick (yg-boost-users_at_[hidden])
Date: 2003-04-28 04:07:37

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