Boost logo

Ublas :

Subject: Re: [ublas] Patch/proposal for overloading operator * in ublas
From: Jesse Perla (jesseperla_at_[hidden])
Date: 2009-08-18 21:58:03

On Tue, Aug 18, 2009 at 11:12 AM, Rutger ter Borg <rutger_at_[hidden]>wrote:

> Karl Meerbergen wrote:
> after all these discussions... what about if I would start stopping
> reinventing the wheel, and use GLAS as a high-level "engine" for the
> bindings? I see BLAS/LAPACK support is there but not as complete as it
> could
> be. Could you bring me up to speed w.r.t. the expression engine of GLAS?

(The following is a plea to extend the existing ublas. I am trying to focus
on economics research and keep finding myself implementing basics... I
imagine many people doing scientific computing are in the same boat).

GLAS would be great, and is a possible front-end to UBLAS in the future.
But it isn't deployed or in boost yet.

The discussion on boost developer and here make it pretty clear that an
immediate, temporary fix would be extremely useful (and hopefully help focus
the efforts of other people on extending instead of replacing). I really
think that a little bit of effort on Matlab'ing UBLAS could go a long way.
Then the longer term approach can be evaluated.

As much as I love the bindings, some of these might be of higher importance
to pound through. It seems that the ublas interface hasn't changed much
since 2002ish, so it probably is time.

I have spent a lot of time using ublas in the last year and trying to get
Matlab programmers to use it. From my experience, the following is a list
of the things that would really help people out in order of importance (I
have tried to remove all features dependent on external libraries).

*1. Operator **

   - Operator * semantics relatively compatible with MATLAB. Including
   inner/outer product and nested products.

*2. Data Adaptors for C/Fortran*

   - There is code floating around for adaptors, but we really need to have
   a formal way of adapting a C or Fortran vectors for ublas operations
   - Also, we need to have a better way to be able to pass the underlying
   UBLAS data to C or Fortran read/write routines. I know all of the
   difficulties with this, but it is essential for practical use. There is
   just too much Fortran code out there. I keep passing in: to just get a double* and hope for the best. I know
   other people are in this boat.
   - For data updating operations, it is perfectly reasonable to force
   people to create an adaptor to get a pointer to the underlying data (which
   may copy the data to a temporary if the storage ordering isn't what people
   thought it would be).
   - I would guess that these exact same problems keep creeping up in the
   bindings, especially for the higher level ones. Maybe the code to do this
   is already built somewhere.

*3. Some temporary solution to the assignment of fixed matrices*.

   - This has been a huge stumbling block for people to use the library...
   much more than I would have guessed. Perhaps something like boost::assign
   could be used until C++0X uniform initialization
   - We should assume that everything is added in row major (and that the
   assignment would take care of the proper storage orderings)

using namespace boost::assign;
matrix<double> A(2,2);
A += 1, 2,
        2, 1;

   - I believe that the += and , overloads for this could work as the ,
   might build a vector, then the += could do assignment to the underlying
   data. It might be reasonable to only allow the assignment for certain types
   of matrices and all of them preallocated. I assume this is possible since
   Blitz does it.

*4. Matrix/Vector Stacking

   - One of the things matlab is very good at is stacking matrices and
   - e.g. Y = [A B]; or Y = [A; B]; in matlab to stack horizontally or
   - This pattern comes out all over the place in algorithms, and would make
   an enormous difference in legibility of code if it were implemented
   - I imagine that the results would be a new matrix and vector expression
   that keeps track of stacking. It would need to implement the () and [] to
   go to the appropriate one it is storing... but with the expression
   templates, I can see this being very efficient.
   - What should the semantics be? It could overload an operator like , |
   and && but probably better to just have functions.

matrix<double> A, B;
vector<double> a, b;

matrix<double> C = stack_vertical(A, B);
auto Cp = stack_vertical(A, B); //This is lazy. The Cp expression has done
no copying.
matrix<double> C = stack_horizontal(a, b);

//or with overloaded operators, something like:
matrix<double C = A & B; //Horizontal
matrix<double C = A & B |
                            B & A; //Stacks like [A B; B A] in matlab.

//or could have some kinds of tags to denote newlines:
matrix<double> C = stack(A, B);
auto Cp = stack(A, matrix::newline, B); //This is lazy. The Cp expression
has done no copying.
matrix<double> C= stack(A, matrix::newline, B); //This is lazy. The Cp
expression has done no copying.

*Interoperability with boost::multi_array*

   - I haven't heard many people talk about this, but it is crucial for a
   huge number of applications.
   - Matlab allows you to get 1 or 2D slices from their multidimensional
   array and then use normal matrix operations on it.
   - We need the ability to adapt boost::multi_array slices/references/etc.
   as ublas matrix and vector types with read/write operations (It is critical
   that this allows modifiable references).
   - This shouldn't be too hard since the multi_array manages all of the
   strides/etc. and has very well thought through iterators.
   - It would be amazing if there was a cast from the multi_array<double, 2>
   reference type to whatever matrix_expression of doubles we build, and
   multi_array<double, 1> to vector_expression.
   - This would be a HUGE value for a lot of work and showcase
   interoperability of boost, and suggest to people they shouldn't try to
   reinvent the wheel here
   - What would I like to be able to do:

boost::multi_array<double, 3> covariances(boost::extents[2[[2][10]); //10
covariance matrices over time.

//Fill in these covariances over time...
boost::matrix<double> A, B;
covariances[boost::indices[range(0,2)][range(0,2)][0] ] = A * B; //Matrix
multiplication and assign to 0'th time period...

boost::matrix<double> C;
cout << 2 * C * covariances[1]; //This gets the 1'th slice along the

*6. A couple of simple linear algebra operations*

   - While in general it is a bindings issue, many people want a couple of
   very basic linear algebra operations out of the box. It is perfectly
   reasonable to provide a few of them with decent semantics, even if they
   aren't in BLAS>
   - For the most part, I believe that there are existing ublas only
   implementations of these. In particular:
      - inverse A = inv(B). Various implementations are out there using the
      LU decomposition, etc. in UBLAS
      - determinant. existing using lu
      - cholesky().
      - kron
      - trace
      - mldivide (see the Kian Ming implementation)
      - A simple solve() using the LU decomposition.
      - I have collected a bunch of these, but Kian Ming's are probably much

*7. Move Semantics*

   - To enable matlab like access, we would really like to be able to write
   our own functions and do things like:

matrix<double> A_inv = inv(A);

   - For some inv() function. Of course, this is only efficient if done
   with rvalue references and move semantics, but my understanding is that
   boost::move will work with older compilers if we enabled some of the ublas
   constructors to be move aware. Am I wrong? If not, this could be extremely
   useful and enable much cleaner semantics on

It appears that Kian Ming has already done much of this work. Are there
people who are competent enough to merge these sorts of changes into a
future ublas version?