Multiplication of each matrix column by each vector element using Eigen C++ Library

Question

I need to multiply each matrix column by each vector element using Eigen C++ library. I tried colwise without success.

Sample data:

Eigen::Matrix3Xf A(3,2); //3x2
A << 1 2,
     2 2,
     3 5;

Eigen::Vector3f V = Eigen::Vector3f(2, 3);

//Expected result
C = A.colwise()*V;

//C
//2 6,
//4 6,
//6 15
//this means C 1st col by V first element and C 2nd col by V 2nd element.

Matrix A can have 3xN and V Nx1. Meaning (cols x rowls).

Focal

Answer 1

This is what I would do:

Code

Eigen::Matrix3Xf A(3, 2);  // 3x2
A << 1, 2, 2, 2, 3, 5;

Eigen::Vector3f V = Eigen::Vector3f(1, 2, 3);

const Eigen::Matrix3Xf C = A.array().colwise() * V.array();
std::cout << C << std::endl;

Example output:

 1  2
 4  4
 9 15

Explanation

You were close, the trick is to use .array() to do broadcasting multiplications.

colwiseReturnType doesn't have a .array() method, so we have to do our colwise shenanigans on the array view of A.

If you want to compute the element-wise product of two vectors (The coolest of cool cats call this the Hadamard Product), you can do

Eigen::Vector3f a = ...;
Eigen::Vector3f b = ...;
Eigen::Vector3f elementwise_product = a.array() * b.array();

Which is what the above code is doing, in a columnwise fashion.

Edit:

To address the row case, you can use .rowwise(), and you'll need an extra transpose() to make things fit

Eigen::Matrix<float, 3, 2> A;  // 3x2
A << 1, 2, 2, 2, 3, 5;

Eigen::Vector2f V = Eigen::Vector2f(2, 3);

// Expected result
Eigen::Matrix<float, 3, 2> C = A.array().rowwise() * V.transpose().array();
std::cout << C << std::endl;

Example output:

 2  6
 4  6
 6 15

Integument

Thank you Jacob, me and my colleague were trying to solve this. You know what is odd, the "unrolled" code runs a bit slower than using a for loop. – Perfection

Sorry, my falt, what I need is to multiply each column, there is an error in the question, vector V is 2x1. I'm going to edit it. Can you please add also this case to your solution? Thank you very much. – Focal

Happily! I wonder if it the loop is faster the compiler is doing a better job of generating SIMD instructions. – Integument

Thanks, working now. A lot of trouble to unroll the for loop but no speed up. Any advise? Thad is sad. – Focal

Hm, what is this for? Your other questions suggest some kind of raytracing application? If so, you'll get much more mileage out of an acceleration structure, like hierarchical bounding boxes, than speeding up code. I've also been not-so-impressed with Eigen's vectorization for this wacky array operations. Manually writing AVX/SSE intrinsics might get you further. – Integument

Wow, exactly ! We have no experience in AVX/SSE. – Perfection

This is quite confusing, I mean when to use mat.array(), since according to the official broadcasting code example, we don't have to use .array() at all, see here: eigen.tuxfamily.org/dox/… – Cereus

Thank you! I "think" I like Eigen, but this is another example of something that should be much easier made hard. – Lorettalorette

What would you expect? If you have Matrix/vector at hand then the solution I proposed is pretty concise and clear, and if you have Array at hand then you can omit the array() views, and the rowwise syntaxe becomes quite concise and clear too. – Ita

Answer 2

In other words, you want to scale each column by a different factor, that is, apply a non uniform scaling. Scaling are best represented as a diagonal matrix, thus:

C = A * V.asDiagonal();

Since Eigen is based on expression template, this does not create any temporary and amount to a code similar to Jacob's answer:

C = A.array().rowwise() * V.transpose().array();

Ita

I think the second code snippet doesn't work as-is, at least I get a Error C2338 YOU_PASSED_A_COLUMN_VECTOR_BUT_A_ROW_VECTOR_WAS_EXPECTED (which makes sense to me, looking at the code). – Pyrogenous

The simplicity (and mathematical elegance) of the first is much nicer than the second, but I think it is an O(N^3) operation for an NxN matrix? The second should be an O(N^2) if implemented well. – Lorettalorette

both are O(N^2) because V.asDiagonal() returns an expression of a diagonal matrix (it is a simple view), so operator * knows what to do! – Ita