Eigen::Ref calls copy constructor on matrix multiplication

Question

I have the following code:

#include <Eigen/Eigen>
#include <iostream>

int main(int argc, char **argv) {
    std::cout << EIGEN_WORLD_VERSION << "." << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION << std::endl;

    Eigen::Vector4f vector{
            0.0f,
            0.0f,
            1.0f,
            0.0f
    };

    const Eigen::Matrix4f M = Eigen::Matrix4f().setOnes();
    Eigen::Vector4f out_vec;
    const Eigen::Ref<const Eigen::MatrixXf> M_ref(M);
    out_vec = M * vector;
    out_vec = M_ref * vector;
    return 0;
}

For some reason M_ref * vector calls Matrix copy constructor and allocates memory (stops at __libc_malloc breakpoint). However it seems like unexpected behavior for Ref class. Did anyone encounter this?

If I replace const Eigen::Ref<const Eigen::MatrixXf> M_ref(M); with fixed size variant const Eigen::Ref<const Eigen::Matrix4f> M_ref(M); then no malloc but copy constructor is still called.

Eigen version is 3.4.0. GCC 11.4.0 -std=gnu++17

Answer 1

This is completely unrelated to Eigen::Ref. The problem is that when Eigen assigns a matrix-matrix or matrix-vector product to another matrix or vector, it does not know if arguments on the right side alias (i.e. share memory) with the target variable. It therefore evaluates into a temporary and copies (or moves) that into the target.

If you know at compile-time that no alias can occur, you can use .noalias():

out_vec.noalias() = M_ref * vector;

The reason out_vec = M * vector; did not produce memory allocation, but a copy, is that the result has known size of 4 elements and is thus stored only on the stack (that copy may stay in a register and get optimized away, depending on compiler options).

Generally, .noalias() is only required if the expression on the right involves a matrix-matrix or matrix-vector product. For element-wise expressions Eigen assumes that no alias occurs or that each result element only depends on the same element -- this can create unexpected results, when working with overlapping blocks, in which case you may need to add an (...).eval() around the expression, or store into a different variable first.