Boost logo

Boost :

From: Paul C. Leopardi (leopardi_at_[hidden])
Date: 2002-10-25 23:19:08


Michael,
I've tried a few of the things you suggested. Results below.
Best regards

> Date: Fri, 25 Oct 2002 12:13:52 +1000
> From: Michael Stevens <michael_at_[hidden]>
[...] I have reordered some of your comments below:
> > o Do you have any tips on where to tweak uBLAS to get better (ie. O(n) )
> performance?
> Try using std::map as the stroage for your sparse types. I would be
> interested in the performance. The proxy used in map_array has its problems.
> I'm trying to convince Joerg to use std::map as the default storage type
> for the sparse types!!

std::map is faster than map_array for this application. There seems to be an
asymptotic improvement, but I can't tell if the order has gone down to
O(n log n) without more extensive analysis and testing:

Size Time Time Slowdown Time Slowdown
       MTL std::map map_array
  32 10 40 4 30 3
  64 20 70 3.5 70 3.5
 128 80 640 8 830 10.4
 256 160 1440 9 1940 12.1
 512 700 20380 29.1 31770 45.4

The speed is still clearly unacceptable in comparison to MTL.

[...]
> > o Is there a more appropriate data type in uBLAS?
> I know Joerg (and a few others) are working on this. I particular a
> row/column compressed matrix type.

It looks like GluCat will need this. I understand this type to be the exact
equivalent of eg.
mtl::matrix< T, mtl::rectangle<>, mtl::compressed<>, mtl::row_major >

[...]
> > But with uBLAS, essentially the same operations seem to be O(n2):
>
> The O(n*n) suggests the matrices have become non sparse. Check non_zeros()!
> uBLAS (Boost version) will not itself loose sparseness when you use BLAS
> expressions. However it is possible, by careless use of indexing on non
> constant objects, to force the sparse types to create zero elements. If you
> do index over all elements make sure you do so via a constant reference so
> uBLAS knows not to create zeros. Alternatively use iterators, which only
> access the non zeros.

>From your comments and from examining the uBLAS documentation, it looks like
the following two routines are likely to be much slower with uBLAS
sparse_matrix than with MTL. Both of these access the matrix via row and
column indices.

The first, inner(), accesses the constant matix y:

  /// Inner product: sum(x(i,j)*y(i,j))/x.nrows()
  template< class Matrix_T, class Scalar_T >
  Scalar_T
  inner(const Matrix_T& x, const Matrix_T& y)
  {
    Scalar_T result = 0;
    for (typename Matrix_T::const_iterator1 i = x.begin1();
         i != x.end1(); ++i)
      for (typename Matrix_T::const_iterator2 j = i.begin();
         j != i.end(); ++j)
      {
        Scalar_T yij = y(j.index1(),j.index2());
        if (yij != Scalar_T(0))
          result += (*j) * yij;
      }
    return result / x.size1();
  }

The second, kron(), builds the matrix z, element-by-element:

  /// Kronecker tensor product of matrices - as per Matlab kron
  template< class Matrix_T >
  void
  kron(const Matrix_T& x, const Matrix_T& y, Matrix_T& z)
  {
    const int rx = x.size1();
    const int cx = x.size2();
    const int ry = y.size1();
    const int cy = y.size2();
    z.resize (rx*ry, cx*cy);
    z.clear ();
    for ( typename Matrix_T::const_iterator1 i = x.begin1();
               i != x.end1(); ++i)
      for ( typename Matrix_T::const_iterator2 j = i.begin();
               j != i.end(); ++j)
        for ( typename Matrix_T::const_iterator1 k = y.begin1();
               k != y.end1(); ++k)
          for (typename Matrix_T::const_iterator2 l = k.begin();
               l != k.end(); ++l)
              z(j.index1()*ry + l.index1(), j.index2()*cy + l.index2()) =
                      (*j) * (*l);
  }

I tried using row proxies speed up either or both of these. Proxies worked for
kron(), but had essentially no effect on speed.

The following code for inner():

  /// Inner product: sum(x(i,j)*y(i,j))/x.nrows()
  template< class Matrix_T, class Scalar_T >
  Scalar_T
  inner(const Matrix_T& x, const Matrix_T& y)
  {
    Scalar_T result = 0;
    for (typename Matrix_T::const_iterator1 i = x.begin1();
         i != x.end1(); ++i)
    {
      boost::numeric::ublas::matrix_row<Matrix_T> yr(y, i.index1());
      for (typename Matrix_T::const_iterator2 j = i.begin();
           j != i.end(); ++j)
      {
        Scalar_T yij = yr(j.index2());
        if (yij != Scalar_T(0))
          result += (*j) * yij;
      }
    }
    return result / x.size1();
  }

gave the following error message (with BNU substituted for
boost::numeric::ublas to reduce clutter):

../glucat/matrix_imp.h:152: no matching function for call to `
 BNU::matrix_row<BNU::sparse_matrix<long double, BNU::row_major,
 std::map<size_t, long double, std::less<size_t>,
 std::allocator<std::pair<const size_t, long double> > > > >::
 matrix_row(const BNU::sparse_matrix<long double, BNU::row_major,
 std::map<size_t, long double, std::less<size_t>,
 std::allocator<std::pair<const size_t, long double> > > >&, size_t)'

/usr/local/include/boost/numeric/ublas/matrix_proxy.hpp:31: candidates are:
 BNU::matrix_row<BNU::sparse_matrix<long double, BNU::row_major,
 std::map<size_t, long double, std::less<size_t>,
 std::allocator<std::pair<const size_t, long double> > > > >::
 matrix_row(const BNU::matrix_row<BNU::sparse_matrix<long double,
 BNU::row_major, std::map<size_t, long double, std::less<size_t>,
 std::allocator<std::pair<const size_t, long double> > > > >&)

/usr/local/include/boost/numeric/ublas/matrix_proxy.hpp:63:
 BNU::matrix_row<E>::matrix_row(M&, M::size_type) [with M =
 BNU::sparse_matrix<long double, BNU::row_major,
 std::map<size_t, long double, std::less<size_t>,
 std::allocator<std::pair<const size_t, long double> > > >]

/usr/local/include/boost/numeric/ublas/matrix_proxy.hpp:60:
 BNU::matrix_row<E>::matrix_row() [with M =
 BNU::sparse_matrix<long double, BNU::row_major,
 std::map<size_t, long double,std::less<size_t>,
 std::allocator<std::pair<const size_t, long double> > > >]

It looks like there is no such thing in uBLAS as a const matrix_row, ie. a
matrix_row is assumed to be non const.


Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk