Matrix Multiplication - Evaluation Direction

[Keith O'Hara]

Hi Klaus and team,

I'm new to Blaze, but I'm a long-time user of another template math library. When playing around with a few tests I noticed something unexpected (on my part) in the expression template handling of matrix-matrix multiplication.

In the example below I am multiplying 3 matrices: A * B * C. With the dimensions I set, the evaluation should occur from right to left, but Blaze doesn't seem to recognize this without some assistance. (The difference in computation time is significant.)

Is this expected behavior? This works in simpler cases, such as matrix * matrix * vector; are these optimizations available only for Level 2-type operations?

blaze::DynamicMatrix<double,blaze::columnMajor> A(1200UL,900UL), B(900UL,600UL), C(600UL,300UL), Z(1200UL,300UL);

blaze::randomize(A);
blaze::randomize(B);
blaze::randomize(C);

Z = A*B*C;

for (size_t j=1; j<=1000; j++)
{
    Z = A*B*C;
//  Z = A*(B*C); compare to this
}

I'm compiling on MacOS with Clang 6 and flags:

-std=c++14 -Ofast -march=native -DNDEBUG

Best, Keith

[Klaus Iglberger]

Hi Keith!

Thanks a lot for pointing this out. You are absolutely correct, in the triple matrix product in your example needs to be evaluated from right to left and Blaze does not perform the necessary optimization automatically. The reason is that this is only one side of the medallion. Consider this counter example:

blaze::DynamicMatrix<double,blaze::columnMajor> A(300UL,600UL), B(600UL,900UL), C(900UL,1200UL), Z(300UL,1200UL);

blaze::randomize(A);
blaze::randomize(B);
blaze::randomize(C);

Z = A*B*C;  // Needs to be evaluated from left to right

In this example the triple matrix product needs to be evaluated from left to right. As this example proves, both evaluation directions are possible. However, the information from which side the expression needs to be evaluated depends on runtime parameters. The generation of the expression templates, however, is a compile time process, which uses the natural evaluation order (which is left to right for multiplication).

In case of the matrix-matrix-vector multiplication, the optimization can be performed at compile time, since the number of columns of the vector is known at compile time. This is why this optimization is implemented in Blaze. In fact, Blaze was the first library to implement this optimization (see the Meeting C++ 2014 talk Expression Templates Revisited).

If you prefer the right to left evaluation order, it is easy in Blaze to implement the necessary optimization manually. The following code snippet contains an overload of operator*(), which you can copy and paste into your code to enforce the right to left evaluation order:

namespace blaze {

template< typename MT1, typename MT2, bool SF, bool HF, bool LF, bool UF, typename MT3, bool SO >
decltype(auto) operator*( const TDMatTDMatMultExpr<MT1,MT2,SF,HF,LF,UF>& lhs, const DenseMatrix<MT3,SO>& rhs )
{
   return lhs.leftOperand() * ( lhs.rightOperand() * (~rhs) );
}

} // namespace blaze

In summary: We will not change the evaluation order at this point in time. However, this doesn't mean we will not do anything. In Blaze 3.4 or 3.5 we plan to provide an optimized triple matrix product, which should go beyond the performance of two consecutive matrix multiplications. In addition, thanks to this optimization, the evaluation order will become irrelevant. Till then, it will unfortunately be necessary to think about the correct evaluation order.

Thanks again for creating this issue,

Best regards,

Klaus!

[Keith O'Hara]

Great! I look forward to those future developments.

Thanks for your reply, Klaus!

Best, Keith