Re: [ublas] Matrix-vector and vector-vector operations

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);
#else
  for (std::size_t i = 0; i < A.size1(); ++i)
  {
    row(A, i) += b;
  }
#endif
  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 
> 
> http://www.crystalclearsoftware.com/cgi-bin/boost_wiki/wiki.pl?Frequently_Asked_Questions_Using_UBLAS

Will do once I get something that works.

Thanks again.

Cheers,

Albert