I have needed to implement some of my research code in C++ and I have been using the Eigen library for linear algebra. There are a lot of great things about Eigen, however, there are some hang-ups I have run into that do not seem to be documented anywhere. As I run into these issues I am going to document how to deal with them.

This post focuses on some peculiar broadcasting syntax, specifically, that the visitor needs to appear on the left hand side of the operation that is broadcasting. This feature conflicts with math in that scalar multiplication is commutative, i.e. $x * y = y * x$ and so is actually hard to track down from the compiler error.

As a concrete example, we consider a similar example that mimics the use-case I encountered in my research. The full code can be found at this gist (Note that the code in the gist uses ArrayXXd and ArrayXd rather than MatrixXd and VectorXd).

Suppose we have a $N \times p$ matrix $A$ and $p$-vector $a$ where we want to pointwise multiply the rows of $A$ by $a$. In Matlab this would be accomplished with the code: A .* repmat(a, 2, 1), and in numpy we would write A * a. The equivalent code in Eigen (with some setup code) looks like

    MatrixXd A(2,4);
    VectorXd a(4);

    A << 1, 2, 6, 9,
         3, 1, 7, 2;
    a << 0, 1, 2, 3;

    MatrixXd res(2,4);
    res = A.rowwise() * a.transpose();

    // res = 0 2 12 27
    //       0 1 14 6

This code will compile fine and produce the correct output. However, if we make a mathematically innocuous change so that the last line becomes

    res = a.transpose() * A.rowwise();

then the code no longer compiles with an error saying that operator* is not defined for the types of a.transpose() and A.rowwise(). To reiterate the main point: when broadcasting, the visitor must be on the left hand side of the operator being used. This is not the case with both Matlab or python where repmat(a, 2, 1) .* A and a * A will work as well.

Hopefully, this is just an oversight and there is an easy fix, and not a shortcoming of the expression generation that Eigen performs.