[boost] Re: uBLAS vs MTL for sparse matrices

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.