Re: Re: uBlas - diagonal_matrix

28 Apr 2003

      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

Re: Re: uBlas - diagonal_matrix

Patrick Kowalzick