Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Does Eigen assume aliasing?

Tags:

c++

eigen

Matrix multiplication is the only operation in Eigen that assumes aliasing by default.

MatrixXf mat1(2,2); 
mat1 << 1, 2,  4, 7;
MatrixXf mat2 = mat1;
auto result = mat1 * mat2;

Eigen evaluates the product mat1 * mat2 in a temporary matrix which is then used to initialize result after the computation. As result does not appear in the right hand side, we do not need aliasing:

MatrixXf result;
result.noalias() = mat1 * mat2;

Now, the product mat1 * mat2 is directly evaluated into result.

So far, so good. But what happens in this case?

template <typename T1, typename T2>
auto multiplication(const T1& A, const T2& B) // I'm using C++17, decltype not needed
{
    return A*B;
}

int main()
{
    auto result = multiplication(mat1, mat2); // say mat1 and mat2 are the same as above
    // or even this
    mat1 = multiplication(mat1, mat2);
    return 0;
}

I would say no aliasing occurs as the multiplication(m1,m2) is an rvalue and constructed directly in result thanks to RVO. And I would say the same for the line mat1 = multiplication(mat1, mat2). We could then say that there is a way to multiply mat1 with another matrix and store the result in the same matrix mat1 without using a temporary matrix (so, avoiding aliasing).

Question:

Does Eigen assume aliasing here or is my assumption correct?

like image 249
mfnx Avatar asked Aug 09 '19 10:08

mfnx


1 Answers

You should also read the Common Pitfall regarding using the auto keyword.

If you write

MatrixXf mat1, mat2;
auto result = mat1 * mat2;

or

template <typename T1, typename T2>
auto multiplication(const T1& A, const T2& B) { return A*B; }

then the type of auto actually is just something like Product<MatrixXf, MatrixXf> or Product<T1,T2>, i.e., no computation at all happens at that point.

Therefore

MatrixXf mat1 = MatrixXf::Random(2,2), mat2 = MatrixXf::Random(2,2);

auto result = multiplication(mat1, mat2); // no computation happens here

// this is safe (Eigen assumes aliasing can happen):
mat1 = result; // equivalent to directly assign mat1 = mat1 * mat2;
// Pitfall: "result" now refers to a modified `mat1` object!

// this will give undefined results (you may need bigger matrices to actually see this):
mat1.noalias() = mat1*mat2; // tell Eigen this does not alias, but actually it does.

Addendum: In the comments the difference between assignment and initialization has been pointed out. In fact, during initialization Eigen assumes no aliasing happens, e.g., the following directly assigns to result (without temporaries):

MatrixXf result = mat1 * mat2; // initialization, not assignment!

Addendum 2: If you wrote (assuming the return type of foo is Object):

Object A;
A = foo(A);

there must be some kind of implicit assignment happening (with C++11 likely a move-assignment, if Object allows that). This is different to

Object A;
Object B = foo(A); // RVO possible (depending on foo).
like image 104
chtz Avatar answered Nov 18 '22 10:11

chtz