CODEWITHSUNDEEP
c++performancematrixeigen3
Frank-Rene Schäfer
Frank-Rene Schäfer

Reputation: 3352

Eigen C++; In-place matrix multiplication

Using the Eigen C++ Matrix library, how can one efficiently multiply an n x n matrix A with a n x m matrix B and store the result in A? That is, how can I avoid generating a temporary n x m matrix and instead store the result directly in B?

For applications where m is very much larger (e.g. 100000) than n (e.g. 3), this definitely makes sense since it avoids the application of very a large array.

The following code, I cannot get to work:

     B.noalias() = A * B;

I would think, that what would have to happen internally is the following. Each column of B should be treated separately. The column under consideration column_i must be copied to a backup column column_tmp. Then,

     B(row_i, column_i) = A.row(row_i) * column_tmp; // dot-product

for all column_i = 0 to m. Is there a way in Eigen to do this efficiently and profit from its optimizations?

Upvotes: 2

Views: 4063

Answers (3)

ggael
ggael

Reputation: 29205

Yes with a 3x3 times 3xHuge your column-by-column evaluation would indeed make sense, but not in the general case. For instance, if n=m=1000, then the column-by-column evaluation would be orders of magnitude slower than the current Eigen logic.

If you write:

B.noalias() = A * B;

Eigen will follow a col-by-col evaluation (because A is small at compile-time), but the result will be wrong because B does alias, basically it will generate:

for j = 0..m-1
  B.col(j).noalias() = A * B.col(j);

To solve this issue elegantly we would need a way to say that only different columns do not alias... The proposed:

B.transpose() *= A.transpose();

is indeed an option let Eigen knows at compile time this kind of information, though having to transpose both sides is a bit cumbersome. And the proper evaluation logic still needs to be implemented on Eigen's side. Currently this information is not exploited.

Upvotes: 2

Massimiliano Janes
Massimiliano Janes

Reputation: 5624

The most explicit way of telling Eigen that you want the product to happen "in place" could be:

B.transpose() *= A.transpose();
// or B.applyOnTheLeft(A);

that said, there's no guarantee this won't incur in any temporary, you'll have to trust internal Eigen cost logic for that (Eigen designers probably know better :) ... or check it yourself via a debugger, after proper profiling suggested this to be a real problem and not just premature optimization).

On my Eigen copy (3.2.7), the above invokes MatrixBase::applyThisOnTheRight directly on the Transpose expression, that in turn sadly reduces to B=A*B internally; the same happens as of applyOnTheLeft, so you're out of luck in this case.

If you really need to avoid any nxm temporary, you can perform the product manually vector-wise, something like:

for(auto i=0;i<B.cols();++i)
    B.col(i) = A * B.col(i);

this will consume much less extra memory assuming B.rows()<<B.cols(), but you may miss some important optimizations here; indeed, I guess having a temporary could still give the best trade-off here.

Upvotes: 5

Caleth
Caleth

Reputation: 62576

Your example of B.noalias() = A * B; does not match "store the result in A". All you need for that is A *= B;. If you instead do intend to overwrite B, then you are lying with .noalias()

Upvotes: 1

Related Questions