Boost logo

Ublas :

From: Albert Strasheim (fullung_at_[hidden])
Date: 2007-05-30 15:37:04

Hello all

On Wed, 30 May 2007, Gunter Winkler wrote:

> On Wednesday 30 May 2007 14:38, Albert Strasheim wrote:
> > I would like to write:
> >
> > datamatrix -= meanvec;
> > datamatrix /= varvec;
> >
> > and sometimes
> >
> > datarow -= meanvec;
> > datarow /= varvec;
> >
> > Unless I am mistaken, one can't do this with uBLAS as-is?
> You can do this by using outer products. You simply need a vector with all
> elements equal to 1.
> matrix A(5,2); vector b(5); scalar_vector e(2,1.0);
> A += outer_prod(b, e); // add b to each column of A

This seems like a viable approach for small A, but I'm running into
problems when dealing with large matrices. Code:

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>

using namespace boost::numeric::ublas;

int main(int argc, char* argv[])
  matrix<double> A(196608, 1024);
  vector<double> b(A.size2());
  scalar_vector<double> e(A.size1(), 1.0);
#if 0
  A += outer_prod(e, b);
  for (std::size_t i = 0; i < A.size1(); ++i)
    row(A, i) += b;
  return 0;

I think I adapted your code correctly to add b to each row of A. I did a
quick test that seemed to confirm this.

However, it seems the outer_prod doubles the memory requirements of
this program compared to using a for-loop. Matrix A is 1536 MB, which
can be allocated on my 32-bit Linux system without problems. However, I
get a bad_alloc exception when trying to add the outer product. Indeed,
for smaller sizes of A I can see the resident memory size of my program
doubles while it does the outer product.

Do you have any other suggestions for how to accomplish this without
needing extra memory?

> element_prod(A, outer_prod(b,e)); // a(i,j) *= b(i)
> > At some distant point in the
> > future one might want to support broadcasting [1].
> "Broadcasting" is simply computing outer products with all-one vectors.

So noted. Thanks.

> PS: please add you solution to the FAQ in

Will do once I get something that works.

Thanks again.