Patch/proposal for overloading operator * in ublas

[Jesse Perla]

As discussed in previous posts, people really can't get past the lack of * overloading in ublas.  I know this was partially a philosophical/design decision, but it is extremely important to "joe user".  Can this be easily remedied?

Here is a simple patch for discussion that implements:
    matrix<double> A;
    vector<double> y;

    cout << 2.0 * y << endl;
    cout << y * 2 << endl;
    cout << A * 2 << endl;
    cout << 2 * A  << endl;
    cout << A * y << endl;
    cout << y * A  << endl;
    cout << A * A << endl;

The patch can be applied in the ublas directory.  All it does is:

• For the scalar * vector and scalar * matrix operator overloads, it uses boost::enable_if to remove anything from those overloads that doesn't fit boost::is_arithmetic.  This is to prevent a clash with the new operator *
  
• Then, new * operators are added for matrix*vector, vector*matrix, and matrix*matrix, which forward directly to the appropriate prod() dispatcher.
  

What is missing:

• Is this as simple as I think it is?  Is there any chance of overhead here?
• What about y * y for two vectors?  I think we should hold off on this because of the inner_prod vs. outer_prod confusion.  Perhaps with the extended numeric_traits created for the bindings this could be used to emulate some kind of matlab approach with explicitly stated orientation (row vs. column vectors), but it may be more confusing than not.  I personally hate the conformance bugs that come out of * meaning both inner and outer product in matlab.
  
• What about b *= A; and A *= B?   Am I right to think this kind of aliased multiplication should be using axpy_prod(b, A, b, false); and axpy_prod(A, B, A, false); ?  Should we just call this in the overload?  This operation comes up all over the place in the statistical work that i use it for.
  

The big one, of course, is nested matrix multiplications..  I don't entirely know how to do it, but here is my proposal:

• Operator * will do a (relatively) inefficient nested matrix multiplication using a heap allocated, row-major matrix<T> as the temporary object.  What type to use for T?  Use some kind of type promotion on the value_type inside of the two argument types....
• The docs tell people what it does, and how to have finer control with the prod(A, B, temp); approach.
  
• In line 4834, etc. of matrix_expression it tests for the complexity at compile time.  Instead, in the operator* overload, we could detect the complexity at compile time and tag-dispatch to a function that would create the appropriate temporary object based on the rules I mentioned above.  It would then call the prod() dispatcher.  We probably would have to do this for the vector*matrix operator overloads as well.
• Would this work?  Any overhead or gotchas?  I don't think I am capable of implementing it myself, but can help experts out with testing, etc?
  

-Jesse


_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Rutger ter Borg-2]

Jesse Perla wrote:

> As discussed in previous posts, people really can't get past the lack of *
> overloading in ublas.  I know this was partially a philosophical/design
> decision, but it is extremely important to "joe user".  Can this be easily
> remedied?
>

I would second this, I don't think there are technical limitations for it?

Perhaps

x = trans(a) * b     // inner prod
x = prec(trans(a)*b) // prec inner prod
x = a * trans(b)     // outer prod

might be considered too? If it will work like that work matrices, why not
for vectors? Fewer keywords, same expressiveness.

Cheers,

Rutger



_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Karl Meerbergen-2]

Rutger ter Borg wrote:

> Jesse Perla wrote:
>
>  
>> As discussed in previous posts, people really can't get past the lack of *
>> overloading in ublas.  I know this was partially a philosophical/design
>> decision, but it is extremely important to "joe user".  Can this be easily
>> remedied?
>>
>>    
>
> I would second this, I don't think there are technical limitations for it?
>
> Perhaps
>
> x = trans(a) * b     // inner prod
> x = a * trans(b)     // outer prod
>  

That is what glas uses, after a long discussion in the mailing list
several years ago.

Karl

> might be considered too? If it will work like that work matrices, why not
> for vectors? Fewer keywords, same expressiveness.
>
> Cheers,
>
> Rutger
>
>
>
> _______________________________________________
> ublas mailing list
> ublas@...
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: karl.meerbergen@...
>  

Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm

[Karl_Meerbergen.vcf]

---

begin:vcard
fn:Karl Meerbergen
n:Meerbergen;Karl
org:K.U. Leuven;Department of Computer Science
adr:Room 02.014;;Celestijnenlaan 200A;Heverlee (Leuven);;B-3001;Belgium
email;internet:Karl.Meerbergen@...
tel;work:[phone]
tel;fax:[phone]
tel;cell:[phone]
note:http://www.kuleuven.be/cwis/email_disclaimer.htm
url:http://www.cs.kuleuven.be/~karlm
version:2.1
end:vcard

---

_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Jesse Perla]

On Tue, Aug 18, 2009 at 7:37 AM, Rutger ter Borg <rutger@...> wrote:

x = trans(a) * b     // inner prod

x = prec(trans(a)*b) // prec inner prod
x = a * trans(b)     // outer prod

Sounds like the winner.  And matlab people will find it intuitive.  And it is consistent with the traits discussed for the new numeric bindings to LAPACK if I recall.   Any ideas on the nested prods?  Could my suggestion to dispatch on complexity to introduce temporaries work?

-Jesse

_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Rutger ter Borg-2]

Karl Meerbergen wrote:

> That is what glas uses, after a long discussion in the mailing list
> several years ago.
>
> Karl

Karl,

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?

Are you considering for submission / review at some point?

If yes :-), do you have an overview of currently supported syntax /
expressions for LAPACK / generic BLAS related stuff?

Have you considered using proto?

Cheers,

Rutger




_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[K.M.A.Chai]

Hi,

    Perhaps this is a nice time to introduce a matlab-like extension to
ublas that I've been working on:

http://homepages.inf.ed.ac.uk/s9810791/umatlab.tbz

----
This package provides "syntactic sugar" for doing matrix programming in
C++. Among others, it is built on top of

- ublas in boost
- numeric-bindings in boost
- multi_array in boost
- matio
- lbfgs
- toms865.

MATLAB is used as a reference when creating this package. Some
features of this package are

- A matrix may be represented in terms of its QR or Cholesky
factorisations. This is to speed up subsequent processing.

- Loosen the semantics of * in ublas to allow matrix multiplication

- For a given matrix operation that has a BLAS2 or BLAS3 equivalent,
it is possible to easily choose between using native ublas
implementation, or using the bindings. For example,
    Y = usebindings<lower>( X * trans(alias(X)) + b);
may be computed either using the C++ implementation in ublas, or a
blas3 bindings. Selection is via the macro UMATLAB_NO_USEBINDINGS.

- Matrix adaptors that behaves like the pass-by-reference version of
reshape in MATLAB.

- Loading and saving of matlab files through matio, with cell-support
via multi-array.
---

   I would, however, like to say that I've currently stopped working on
this. Anyone is welcomed to use and extend this, though.

   FYI, I'm using this for learning Gaussian processes models.


Kian Ming (Adam)


--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.

_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Hongyu Miao]

It looks like GLAS has done a lot of work to make the Matrix/Vector objects easy to use, which is what UBLAS really needs. I've been using UBLAS for two years and I am also using MATLAB in the mean time; anybody who has experience with both will know what I want to say, :). Really looking forward to the MATLAB-like manipulation of matrices and vectors in UBLAS.

Best,
------------------
Jacky Miao
2009-08-18

-------------------------------------------------------------
From：Rutger ter Borg
Date：2009-08-18 11:12:30
To：ublas
CC：
Subject：Re: [ublas] Patch/proposal for overloading operator * in ublas

Karl Meerbergen wrote:

> That is what glas uses, after a long discussion in the mailing list
> several years ago.
>
> Karl

Karl,

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?

Are you considering for submission / review at some point?

If yes :-), do you have an overview of currently supported syntax /
expressions for LAPACK / generic BLAS related stuff?

Have you considered using proto?

Cheers,

Rutger




_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: jackymiao@...
_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Jesse Perla]

On Tue, Aug 18, 2009 at 11:12 AM, Rutger ter Borg <rutger@...> 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 (read/write).
• 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: mymat.data().begin() 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 vectors.
• e.g. Y = [A B]; or Y = [A; B]; in matlab to stack horizontally or vertically.
• This pattern comes out all over the place in algorithms, and would make an enormous difference in legibility of code if it were implemented efficiently.
• 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.
  

e.g.
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...

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


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 better.
  

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?

_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Rutger ter Borg-2]

Jesse Perla wrote:

[snip]

> 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?

Dear Jesse,

thanks for the valuable feedback, I'm curious to see the response from the
uBLAS developers.

Cheers,

Rutger



_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Jesse Perla]

On Tue, Aug 18, 2009 at 12:17 PM, <K.M.A.Chai@...> wrote:

Hi,
This package provides "syntactic sugar" for doing matrix programming in
C++. Among others, it is built on top of

This looks excellent.  Do you have any regression tests/examples/docs on what these patches do (besides the obvious ones that replace matlab functions of course)?

Also, were you able to deal with the conversion between multi_array slices to ublas?

-Jesse

_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[K.M.A.Chai]

On Wed, 19 Aug 2009, Jesse Perla wrote:

>
> This looks excellent.  Do you have any regression tests/examples/docs on what these patches do
> (besides the obvious ones that replace matlab functions of course)?
>
> Also, were you able to deal with the conversion between multi_array slices to ublas?
>

Hi Jesse,

    Thanks for the interest. As I said in the previous post, I'm using it
for learning GP models. I've a package (also in the state of flux) for
this at http://homepages.inf.ed.ac.uk/s9810791/gpml-c++.tbz.
Unfortunately, there's also no documentation associated with this.

    I'll like to clarify that the only patches are under the directory
umatlab-patches. Much of the package should be seen as "add-on" instead of
patches.

    Unfortunately, it doesn't do the conversion between multi_array slices
to ublas. What it can do is load a cell variable from a mat file into a
multi_arrray of ublas vectors/matrices.

    I can, however, provide the following code fragment for doing
regularised linear regression, and returning the MSE on a test set.

-----
    typedef ublas::vector<double> dvector;
    typedef ublas::symmetric_matrix<double, lower, column_major> dsymmat;
    typedef umatlab::chol_representation<double, lower> dcholrep;

    try {
1    dsymmat XX(usebindings<lower>(trans(alias(X)) * X));
      noalias(XX) += regulariser * identity_matrix<double>(XX.size2());
2    dvector alpha(mldivide( dcholrep(XX), trans(X) * y));
      dvector ypredict( usebindings( Xtest * alpha ) );
      MSE = mean( pow(ypredict - ytest, 2) );
    }
    catch( bad_rcond& rcond ) {
      cout << "#W Bad rcondition: " << rcond.rcond << endl;
      MSE = NAN;
    }
-------

Annotation:
1. trans(alias(Xtr)) * Xtr says that we are doing trans(X)*X, where the
two X's are the same. "alias" is a new container introduced into the
package to annotate such instances for the compiler. Cf. the "noalias" of
ublas.

1. usebindings will parse the expression given as its arguments and try to
match a BLAS2/BLAS3 operation to it. I believe in this case "syrk" is
invoked. If UMATLAB_NO_USEBINDINGS is defined, then usebindings will be a
an identity operation, and the blas implementation will be used.

2. Notice that trans(Xtr)*Tautr is a valid operation, i.e., transpose
matrix times vector

2. dcholrep re-represent the matrix as its Cholesky factors. TOMS-865
algorithm (patched) for cholesky decomposition is currently used.

2. mldivide does a matrix left divide operation (see matlab mldivide).
Currently, it simply depatches the operation depending on whether the
matrix representation is QRP, cholesky, symmetrix matrix, triangular
matrix, or simply a square matrix. A vector or a matrix is returned
depending on its second argument.


Kian Ming (Adam)
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.

_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Rutger ter Borg-2]

K.M.A.Chai@... wrote:
>
> Hi Jesse,
>
>     Thanks for the interest. As I said in the previous post, I'm using it
> for learning GP models. I've a package (also in the state of flux) for
> this at http://homepages.inf.ed.ac.uk/s9810791/gpml-c++.tbz.
> Unfortunately, there's also no documentation associated with this.
>

Hello Kian Ming,

nice to see someone with a similar research interests -- I'll come back to
this within a couple of weeks.

Cheers,

Rutger



_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Gunter Winkler]

Am Tuesday 18 August 2009 schrieb K.M.A.Chai@...:

> Hi,
>
>     Perhaps this is a nice time to introduce a matlab-like extension
> to ublas that I've been working on:
>
> http://homepages.inf.ed.ac.uk/s9810791/umatlab.tbz
>
> ----
> This package provides "syntactic sugar" for doing matrix programming
> in C++. Among others, it is built on top of
>
> - ublas in boost
> - numeric-bindings in boost

This looks very interesting. Do you have contributions which should be
included in the official uBLAS release?

Since the development of uBLAs is quite frozen for a long time it would
be nice to have some new features.

mfg
Gunter




_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Gunter Winkler]

Am Wednesday 19 August 2009 schrieb Rutger ter Borg:

> Jesse Perla wrote:
>
> [snip]
>
> > 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?
>
> Dear Jesse,
>
> thanks for the valuable feedback, I'm curious to see the response
> from the uBLAS developers.

Let's start working on it:

https://svn.boost.org/trac/boost/ticket/3397

I can spent some time reviewing patches, but I won't have enough time to
do the merge myself.

mfg
Gunter


_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Gunter Winkler]

Hello Jesse,

after silently following the discussion for a while I'd like to comment
the findings.

Am Monday 17 August 2009 schrieb Jesse Perla:

> As discussed in previous posts, people really can't get past the lack
> of * overloading in ublas.  I know this was partially a
> philosophical/design decision, but it is extremely important to "joe
> user".  Can this be easily remedied?
>
> Here is a simple patch for discussion that implements:
>     matrix<double> A;
>     vector<double> y;
>
>     cout << 2.0 * y << endl;
>     cout << y * 2 << endl;
>     cout << A * 2 << endl;
>     cout << 2 * A  << endl;

so far I agree

>     cout << A * y << endl;

this is the common matrix-vector multiplication and should definitely
provided.

>     cout << y * A  << endl;

do we really need this? I'd prefer to fix vectors as column vectors and
write (trans(y) * A) or (trans(A) * y)
Otherwise we would have to tag vectors with an orientation.

>     cout << A * A << endl;

here I also agree.

>
> The patch can be applied in the ublas directory.  All it does is:
>
>    - For the scalar * vector and scalar * matrix operator overloads,
> it uses boost::enable_if to remove anything from those overloads that
> doesn't fit boost::is_arithmetic.  This is to prevent a clash with
> the new operator * - Then, new * operators are added for
> matrix*vector, vector*matrix, and matrix*matrix, which forward
> directly to the appropriate prod() dispatcher.

at first glance the patch looks ok. I'll spent some time testing it.

>
> What is missing:
>
>    - Is this as simple as I think it is?  Is there any chance of
> overhead here?

There is a big problem: How to compute the type of the result? If both
arguments are dense, then the answer is simply "dense". But what is the
type of "banded*dense", "sparse*dense" (or even worse "packed vector *
sparse vector", "compressed matrix * sparse vector"). Here we need a
big effort to invent a useful traits mechanism or we simply state to
use the "most dense" or "most sparse" representation.

>    - What about y * y for two vectors?  I think we should hold off on
> this because of the inner_prod vs. outer_prod confusion.  Perhaps
> with the extended numeric_traits created for the bindings this could
> be used to emulate some kind of matlab approach with explicitly
> stated orientation (row vs. column vectors), but it may be more
> confusing than not.  I personally hate the conformance bugs that come
> out of * meaning both inner and outer product in matlab.

here I'd suggest to use the same solution as GLAS. We define that all
vectors are column vectors and use

x = trans(a) * b     // inner prod
x = a * trans(b)     // outer prod

x = prec(trans(a)*b) // prec inner prod


>    - What about b *= A; and A *= B?   Am I right to think this kind
> of aliased multiplication should be using axpy_prod(b, A, b, false);
> and axpy_prod(A, B, A, false); ?  Should we just call this in the
> overload? This operation comes up all over the place in the
> statistical work that i use it for.

Yes, I think this is a good idea. However we have to assure that
axpy_prod correctly handles the aliasing. Currently it just works but
this is currently not required to work.

> The big one, of course, is nested matrix multiplications..  I don't
> entirely know how to do it, but here is my proposal:
>
>    - Operator * will do a (relatively) inefficient nested matrix
>    multiplication using a heap allocated, row-major matrix<T> as the
> temporary object.  What type to use for T?  Use some kind of type
> promotion on the value_type inside of the two argument types....

I think it could be helpful to have a default mechanism to create
temporary types and provide means of fine tuning. However I expect that
this requires weeks of developement. (Do you have some students to
assign the task to?)

>    - The docs tell people what it does, and how to have finer control
> with the prod(A, B, temp); approach.

this looks quite nice.

>    - In line 4834, etc. of matrix_expression it tests for the
> complexity at compile time.  Instead, in the operator* overload, we
> could detect the complexity at compile time and tag-dispatch to a
> function that would create the appropriate temporary object based on
> the rules I mentioned above.  It would then call the prod()
> dispatcher.  We probably would have to do this for the vector*matrix
> operator overloads as well.

Does the product of three general matrices occur often? IMO it most time
appears as symmetric product like trans(X)*A*X, or a scaled product
like invD * A * D. Here I'd prefer to provide specialized functions.

>    - Would this work?  Any overhead or gotchas?  I don't think I am
> capable of implementing it myself, but can help experts out with
> testing, etc?

Don't be shy to start developing uBLAS. Right now there is no active
developer. (note to myself: Spending a few hours per week does not
count as "active".)

mfg
Gunter


_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[K.M.A.Chai]

Hi Gunter,

    Thanks for the interest. My feeling is that ublas and umatlab
really have different aims, so most of parts of umatlab may be out of
place in ublas (though, lu.hpp also seems out of place in ublas).

    That said, it will be nice if ublas can be easily extendable. A lot of
people has been asking about overloading * in ublas. However, currently
ublas is pretty greedy: it tries to handle all of operator*, so that
extending operator* in a naive way (by the end user) will lead to the
compiler complaining. Attached are two patches to make ublas less greedy
in this aspect.

    There are some other patches for ublas: I'll submit them once I recall
why they're there at all :)

   Another possible addition would be the parser that breaks down an
expression into the consitutents for blas calls.

Kian Ming (Adam)



On Thu, 3 Sep 2009, Gunter Winkler wrote:

> Am Tuesday 18 August 2009 schrieb K.M.A.Chai@...:
>> Hi,
>>
>>     Perhaps this is a nice time to introduce a matlab-like extension
>> to ublas that I've been working on:
>>
>> http://homepages.inf.ed.ac.uk/s9810791/umatlab.tbz
>>
>> ----
>> This package provides "syntactic sugar" for doing matrix programming
>> in C++. Among others, it is built on top of
>>
>> - ublas in boost
>> - numeric-bindings in boost
>
> This looks very interesting. Do you have contributions which should be
> included in the official uBLAS release?
>
> Since the development of uBLAs is quite frozen for a long time it would
> be nice to have some new features.
>
> mfg
> Gunter
>
>
>

--

The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.


--- vector_expression.hpp.orig 2008-05-10 15:12:53.000000000 +0100
+++ vector_expression.hpp 2008-05-10 15:14:52.000000000 +0100
@@ -1170,7 +1170,9 @@ namespace boost { namespace numeric { na
     // (t * v) [i] = t * v [i]
     template<class T1, class E2>
     BOOST_UBLAS_INLINE
+    typename enable_if< is_convertible<T1, typename E2::value_type >,    
     typename vector_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::result_type
+    >::type
     operator * (const T1 &e1,
                 const vector_expression<E2> &e2) {
         typedef typename vector_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::expression_type expression_type;
@@ -1395,7 +1397,9 @@ namespace boost { namespace numeric { na
     // (v * t) [i] = v [i] * t
     template<class E1, class T2>
     BOOST_UBLAS_INLINE
+    typename enable_if< is_convertible<T2, typename E1::value_type >,    
     typename vector_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::result_type
+    >::type
     operator * (const vector_expression<E1> &e1,
                 const T2 &e2) {
         typedef typename vector_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::expression_type expression_type;

--- matrix_expression.hpp.orig 2008-05-01 11:11:57.690870000 +0100
+++ matrix_expression.hpp 2008-05-07 13:58:59.998389000 +0100
@@ -2940,7 +2940,9 @@ namespace boost { namespace numeric { na
     // (t * m) [i] [j] = t * m [i] [j]
     template<class T1, class E2>
     BOOST_UBLAS_INLINE
+    typename enable_if< is_convertible<T1, typename E2::value_type >,
     typename matrix_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::result_type
+    >::type
     operator * (const T1 &e1,
                 const matrix_expression<E2> &e2) {
         typedef typename matrix_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::expression_type expression_type;
@@ -3373,7 +3375,9 @@ namespace boost { namespace numeric { na
     // (m * t) [i] [j] = m [i] [j] * t
     template<class E1, class T2>
     BOOST_UBLAS_INLINE
+    typename enable_if< is_convertible<T2, typename E1::value_type>,
     typename matrix_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::result_type
+    >::type
     operator * (const matrix_expression<E1> &e1,
                 const T2 &e2) {
         typedef typename matrix_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::expression_type expression_type;

_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Jesse Perla]

Thanks for the response Gunter.  Everyone definetely appreciates the time you are able to put into the project.  See below:

On Thu, Sep 3, 2009 at 6:39 PM, Gunter Winkler <guwi17@...> wrote:

>     cout << y * A  << endl;

do we really need this? I'd prefer to fix vectors as column vectors and
write (trans(y) * A) or (trans(A) * y)
Otherwise we would have to tag vectors with an orientation.

I think you are probably right.  Having a standard orientation makes life a lot easier, and eliminates all of the conformability errors that plague writing matlab code (where a vector is an N by 1 or 1 by N matrix).

at first glance the patch looks ok. I'll spent some time testing it.

My guess is that K.M.A.Chai's are a much better place to start.  I don't trust my own code very well, and I think that implementation went a lot further.

 

>
> What is missing:

There is a big problem: How to compute the type of the result? If both
arguments are dense, then the answer is simply "dense". But what is the
type of "banded*dense", "sparse*dense" (or even worse "packed vector *
sparse vector", "compressed matrix * sparse vector"). Here we need a
big effort to invent a useful traits mechanism or we simply state to
use the "most dense" or "most sparse" representation.

Is this a problem with just the typical matrix-matrix operations or only with nested ones?  For the existing ones, since the operator overload just forwards to the existing function I think that nothing really changes and it just uses the current traits with result_type?  Or are you saying that the current result type is not especially efficient.
 
But for the nested products, it certainly is an issue.  My proposal is the following:  If it generates a temporary on its own, it chooses the most dense, most dynamic one possible for now (i.e. ublas::matrix<double>), etc.  And to find the data type if does some kind of type promotion on the value_type's of the arguments to the biggest type (i.e. matrix<float> * compressed_matrix<double> -> matrix<double>).

The good thing about this approach is that people will still have the old method to control the type of the temporary, and they should never count on the type of the temporary remaining the same between versions(which I don't think would bug anyone since the temporary is between 2 ublas calls)...so if we later put in a great system for traits of the multiplication, it would happen seamlessly.
 

>    - What about y * y for two vectors?  I think we should hold off on
here I'd suggest to use the same solution as GLAS. We define that all
vectors are column vectors and use

x = trans(a) * b     // inner prod

x = a * trans(b)     // outer prod

x = prec(trans(a)*b) // prec inner prod

Yup, GLAS compatibility is perfect.  The sooner the GLAS, the happier I am.
 

Does the product of three general matrices occur often? IMO it most time
appears as symmetric product like trans(X)*A*X, or a scaled product
like invD * A * D. Here I'd prefer to provide specialized functions.

Yes, there is a lot of that, so it probably makes sense to have that specialization.  And if there are more efficient algorithms then so much the better.  If we go down that route there are a lot of great specializations to implement (for exmaple, I would love a specific function to do a quadratic form v' A v)

But 3 general comes up frequently enough that we would need to have it work.  For example, anyone from the matlab world will likely type in "invD * A * D" directly or invert to solve regressoins... we can try to untrain bad matlab form, but still need to allow these sorts of things.  Again, I think that having a slow and safe default temporary as discussed above is perfectly fine for this.

 

>    - Would this work?  Any overhead or gotchas?  I don't think I am

> capable of implementing it myself, but can help experts out with
> testing, etc?

Don't be shy to start developing uBLAS. Right now there is no active
developer. (note to myself: Spending a few hours per week does not
count as "active".)

I will do my best, but I am an economist, and this kind of library is extremely difficult to do right.  ublas is such an essential library that it is shameful there isn't more direct help on it, especially considering all of the time that people spend playing around with linear algebra libraries out there.  I think I will make a plea to the boost community...

-Jesse


_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[K.M.A.Chai]

Hi Gunter,

   The following inclusions into matrix_expression.hpp will allow
"matrix * matrix" and "matrix * vector".

<patch>
   template<class E1, class E2>
   BOOST_UBLAS_INLINE
   typename matrix_matrix_binary_traits<typename E1::value_type, E1,
                                        typename E2::value_type, E2>::result_type
   operator * (const matrix_expression<E1> &e1, const matrix_expression<E2> &e2) {
     return prod (e1, e2);
   }

   template<class E1, class E2>
   BOOST_UBLAS_INLINE
   typename matrix_vector_binary1_traits<typename E1::value_type, E1,
                                         typename E2::value_type, E2>::result_type
   operator * (const matrix_expression<E1> &e1, const vector_expression<E2> &e2) {
     return prod (e1, e2);
   }
</patch>

   A relook at this again makes me wonder why "matrix * vector" is in
matrix_expression.hpp, since this is really a vector expression.


Kian Ming (Adam)



--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.

_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Rutger ter Borg-2]

K.M.A.Chai@... wrote:

>     That said, it will be nice if ublas can be easily extendable. A lot of
> people has been asking about overloading * in ublas. However, currently
> ublas is pretty greedy: it tries to handle all of operator*, so that
> extending operator* in a naive way (by the end user) will lead to the
> compiler complaining. Attached are two patches to make ublas less greedy
> in this aspect.

I agree that uBLAS should be extensible -- e.g., a bit like a toolbox in
Matlab and GLAS.

For this to work, one would need to easily export and extend the domain
specific language (DSL). Perhaps it is a good time to let Boost.Proto handle
uBLAS' expression machinery -- is has the facilities we're looking for.

What would be nice is to be able to work with expression trees. In this case
I'm referring to Proto's semantic actions, or expression transformation.
E.g., rewrite an expression to make use of BLAS calls (or, some other
"backend"). Pattern matching on the expression, rewrite, and transparently
an other expression is used. All in compile time.

Cheers,

Rutger



_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Jesse Perla]

On Sep 8, 2009 2:41am, Rutger ter Borg <rutger@...> wrote:
>
> What would be nice is to be able to work with expression trees. In this case
> I'm referring to Proto's semantic actions, or expression transformation.
> E.g., rewrite an expression to make use of BLAS calls (or, some other
> "backend"). Pattern matching on the expression, rewrite, and transparently
> an other expression is used. All in compile time.

It definitely would be nice. I have found that as soon as you start down the path of generic and functional programming in C++, you start wishing you could mix/match this stuff. For example, lambdas binding to specific locations in vectors, auto-differentiation of vectors/matrices, etc. While I am not a good enough programmer to understand the details (and hope I never do...), it sure smells like having a consistent expression template library would make this stuff fit together seamlessly.

That said, I fear that if we worry too much about the long-term goal, we won't make the necessary patches to make ublas a good enough solution to hold us over.

If the GLAS guys are watching this list, I would love to hear their thoughts on the longterm vision and when they think we might be in a position to evolve boost's matrix library to the next level. Based on my experiences, it sure seems to me that GLAS is the right approach, but it is a hell of a lot of work.

-Jesse
_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...