Loading ...
Sorry, an error occurred while loading the content.
 

Re: [ublas-dev] Fw: [Boost-Users] Re: uBlas - diagonal_matrix

Expand Messages
  • jhr.walter@t-online.de
    Hello Patrick, ... Yep. ... But the difference in speed could be ;-) BTW: prod(prod()) is dangerous. Recommended usage is either prod(prod ())
    Message 1 of 3 , Apr 30, 2003
      Hello Patrick,

      you wrote:

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

      Yep.

      > > 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));
      > }

      But the difference in speed could be ;-) BTW: prod(prod()) is dangerous.
      Recommended usage is either prod(prod<some_matrix_type>()) or separating the
      products.

      > while in general I like more the second one, but the first one is
      > more efficient :-).

      Yes.

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

      OK.

      [snip]

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

      Oops. Fixed, thanks.

      Best,

      Joerg
    Your message has been successfully submitted and would be delivered to recipients shortly.