Does Eigen assume aliasing?

Question

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?

Correspondence

You seem to treat something like auto result = mat1 * mat2; as an assignment. – Joiejoin

auto result = mat1 * mat2; is (in C++17) an initialization, not an assignment. It does not use operator=. – Pridgen

@MaxLanghof Ok, I'll edit that sentence (I copy pasted it from Eigen docs). – Correspondence

@MaxLanghof Could you please point me to the sentence that is wrongly expressed? I understand assignment vs initialization, yet fail to see which sentence I expressed bad. – Correspondence

Yes that was a mistake and I already deleted the comment. You were correctly referring to the code above (and not below) the sentence in question, my bad. – Pridgen

Can you multiply two square matrices and put the result in the first one without using any additional storage? Try doing this with exactly two postit notes, prlencil and eraser. – Thieve

Answer 1

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).

Michigan

I did not think about the auto pitfall. I was able to confirm the type of result (Product<T1,T2>). If I rewrote my template like template <typename T> T multiplication(const T& A, const T& B) { return A*B; }, then auto result = multiplication(A,B);, as I understand from your answer, result would be initialized and no aliasing would occur. Correct? – Correspondence

Correct and it should be optimized by RVO -- as long as it still compiles (it won't compile, e.g., for Matrix4f A, Vector4f B, in contrast to your template<T1, T2> solution). – Michigan

Ok, but as you said, auto deduces a type Product<T1,T2>. Maybe, I could rewrite the template as: template<typename T1, typename T2> multiplication(const T1& A, const T2& B) -> decltype((A*B).eval()). If Matrix4f A and Vector4f B, then auto result = multiplication(A,B) would then work: the type of result would be Vector4f and the initialization of result would not imply temporaries. – Correspondence

Yes, auto result = multiplication(A,B); would not have temporaries, but B = multiplication(A,B); would have. And MatrixXf result; result = multiplication(A,B); would also create a temporary (but likely gets move-assigned). Overall, I do not see what you want to achieve in contrast to simply writing result = A*B; – Michigan

Nothing really. The function multiplication was just an attempt to provide a minimal example. The function I use in my code is not a simple product. I could (and should) have left the multiplication part out. – Correspondence