Bug 157 - add *, /, *=, /= to VectorwiseOp
add *, /, *=, /= to VectorwiseOp
Status: RESOLVED FIXED
Product: Eigen
Classification: Unclassified
Component: Core - general
unspecified
All All
: --- enhancement
Assigned To: Nobody
:
: 151 374 (view as bug list)
Depends on:
Blocks:
  Show dependency treegraph
 
Reported: 2011-01-19 19:17 UTC by Evan Herbst
Modified: 2011-11-18 20:06 UTC (History)
5 users (show)



Attachments
Fixed VectorwiseOp.h (23.82 KB, text/x-chdr)
2011-11-11 16:56 UTC, Kibeom Kim
no flags Details
this patch implements some of the operators. (2.71 KB, patch)
2011-11-18 00:13 UTC, Kibeom Kim
no flags Details | Diff

Description Evan Herbst 2011-01-19 19:17:43 UTC
They make sense at least for Array.

What I want to do is divide each row of rho by its sum:

MatrixXf rho(n, m), r(n, m);
r = rho.array().colwise() / rho.rowwise().sum().array();

and I don't see how to do it without a loop.

Thanks.
Comment 1 Marton Danoczy 2011-01-25 09:55:12 UTC
You can do that with .replicate():
r = rho.array().colwise() / rho.rowwise().sum().array().replicate(n,1);

I think this should work as well:
r = rho.array().colwise() / rho.rowwise().sum().array().colwise().replicate(n);

But yes, me too I miss that feature a lot as it enables code to be a lot more concise.
Comment 2 Kibeom Kim 2011-11-11 00:02:42 UTC
*** Bug 374 has been marked as a duplicate of this bug. ***
Comment 3 Kibeom Kim 2011-11-11 00:03:23 UTC
Summary:

From the tutorial,
http://eigen.tuxfamily.org/dox/TutorialReductionsVisitorsBroadcasting.html

mat.colwise() += v; //works
mat.colwise() -= v; //works
mat.colwise() *= v; //not implemented
mat.colwise() /= v; //not implemented

Moreover,

mat2 = mat.colwise() + v; //works
mat2 = mat.colwise() - v; //works
mat2 = mat.colwise() * v; //not implemented
mat2 = mat.colwise() / v; //not implemented

mat2 = v + mat.colwise(); //not implemented
mat2 = v - mat.colwise(); //not implemented
mat2 = v * mat.colwise(); //not implemented
mat2 = v / mat.colwise(); //not implemented
Comment 4 Kibeom Kim 2011-11-11 12:34:15 UTC
In addition,

Eigen::ArrayXf  v(2);
v = 3/v; //not implemented
Comment 5 Kibeom Kim 2011-11-11 16:56:49 UTC
Created attachment 229 [details]
Fixed VectorwiseOp.h

Here is the fixed VectorwiseOp.h.

but need to implement operations for rhs VectorwiseOp. But I think that should be somewhere else, not this file.
Comment 6 Kibeom Kim 2011-11-11 19:31:05 UTC
*** Bug 151 has been marked as a duplicate of this bug. ***
Comment 7 Kibeom Kim 2011-11-18 00:13:05 UTC
Created attachment 230 [details]
this patch implements some of the operators.

this patch implements some of the operators.

Eigen::MatrixXf mat(2,4);
Eigen::ArrayXf  v(2);

mat.colwise() += v; //works
mat.colwise() -= v; //works
mat.colwise() *= v; //not implemented -> now implemented
mat.colwise() /= v; //not implemented -> now implemented

mat2 = mat.colwise() + v; //works
mat2 = mat.colwise() - v; //works
mat2 = mat.colwise() * v; //not implemented -> now implemented
mat2 = mat.colwise() / v; //not implemented -> now implemented

mat2 = v + mat.colwise(); //not implemented
mat2 = v - mat.colwise(); //not implemented
mat2 = v * mat.colwise(); //not implemented
mat2 = v / mat.colwise(); //not implemented

v = 3 / v; //not implemented -> now implemented
Comment 8 Gael Guennebaud 2011-11-18 00:48:26 UTC
Comment on attachment 230 [details]
this patch implements some of the operators.

Review of attachment 230 [details]:
-----------------------------------------------------------------

First, such / operators should be available only in array mode. Otherwise they make no sense to me. The correct way of doing this with linear algebra objects is, e.g.:

C = A * v.diagonal().inverse();

Second, the documentation should mention that these functions are actually equivalent to matrix * diagonal matrix products:
C = A * v.diagonal(); <=> rowwise vector multiply
C = A * v.diagonal().inverse(); <=> rowwise vector divison
Comment 9 Kibeom Kim 2011-11-18 01:10:01 UTC
Comment on attachment 230 [details]
this patch implements some of the operators.

Review of attachment 230 [details]:
-----------------------------------------------------------------

1.
I also think that / operator should be only available for array mode. (but can you teach me how it can be done in a correct way? I'm still learning the library structure and it's not so easy)

2.
C = A / v //newly implemented operator
C = A * v.diagonal().inverse(); //equivalent?
seems your example will work for a scalar A and "ArrayXf v;" but how about for "ArrayXXf v;" ?
Comment 10 Benoit Jacob 2011-11-18 05:13:35 UTC
The relevant part of the class hierarchy is:
 - DenseBase is the common base class for all dense matrix and array expressions
 - MatrixBase inherits DenseBase, is the common base class for Matrix expressions
 - ArrayBase inherits DenseBase, is the common base class for Array expressions
Comment 11 Gael Guennebaud 2011-11-18 08:09:32 UTC
(In reply to comment #9)
 
For your point 2, the equivalents I shown are for the matrix world only where v is a vector, and it is not equivalent to C = A/v but C = A.rowwise() / v where v also has to be a vector.
Comment 12 Benoit Jacob 2011-11-18 15:45:22 UTC
This whole area is really a mess. I propose that we:
  1) take Kibeom's patch adding multiplicative vectorwise ops
  2) add lots of static assertions making it very explicit when all these ops are supposed to work:
    2a) additive ops already require that the lhs and rhs are of the same kind (matrix vs vector)
    2b) multiplicative ops (as in Kibeom's patch) also pretty much have to be restricted to arrays (enabling them for matrices would require a complex dance with .array()/.matrix() which would be fine for *= and /= but would force * and / to return very deeply nested expressions.
  3) add unit tests. there just aren't enough unit tests here. currently, additive ops are only tested in array.cpp on arrays as far as I can see.
Comment 13 Benoit Jacob 2011-11-18 15:48:25 UTC
oh and:

 4) this stuff:

   for(Index j=0; j<subVectors(); ++j)
     subVector(j).array() *= other.derived().array();
is not the right way, is it? In operator + we return a sum with the extendedTo() expression. I suppose that there is no reason for a custom for loop here.
Comment 14 Benoit Jacob 2011-11-18 15:48:58 UTC
(I mean, this problem preexists, it's what += and -= do. I'm criticizing existing code, not Kibeom's patch)
Comment 15 Benoit Jacob 2011-11-18 17:12:31 UTC
Pushed Kibeom's patch: 13f72ac8283d

Pushed another patch implementing my proposal in comment 12 and comment 13, as well as adding a new unit test, vectorwiseop.cpp: d48a7c25e052
Comment 16 Kibeom Kim 2011-11-18 17:42:42 UTC
I suggest to define rowwise() and colwise() like this.

"something.rowwise() [operator] other"
as
"something.row(k) [operator] other for all k"

because this is a better way to define the function cleanly and consistently.
so for example,

1.
"array.rowwise() * v"
should equal to
"array.row(k) * v for all k"

2.
"mat.rowwise() * v"
should equal to
"mat.row(k) * v for all k"
(notice that this is just mat * v)

so for the implementation, VectorwiseOp will hands over all the actual operations to an underlying object, array or matrix.
Comment 17 Benoit Jacob 2011-11-18 19:46:08 UTC
The problem is that for a matrix,

   "mat.row(k) * v"

gives a 1x1 result, containing the real dot product of mat.row(k) with v.

But many people would expect rowwise to behave like it does for arrays. So I'm afraid that this would only add confusion.

Moreover, for a matrix, "mat.row(k) * v for all k" is already available: that's just the matrix-vector product,  mat*v.
Comment 18 Kibeom Kim 2011-11-18 20:06:36 UTC
actually, I think either approach is fine.
(although I prefer "mat.row(k) * v for all k" = mat*v slightly more)

But whatever, I think having a simple generalized statement of what rowwise() actually does is good.

How about this?


"something.rowwise() [operator] other"

is conceptually identical to

"something.array().row(k) [operator] other for all k"

and then we can pass all the VectorwiseOp's operator implementations to "array.row() [operator] other".

Note You need to log in before you can comment on or make changes to this bug.