Effective Eigen: Using the Eigen::Ref Class

Imagine you’re reviewing a junior developer’s code, and they have written this function to find the intersection of two lines, each defined by a pair of points. The points are in 2-D projective space, so are specified by three homogenous coordinates:

Eigen::Vector3d intersection(Eigen::Vector3d pa_1, Eigen::Vector3d pa_2,
                             Eigen::Vector3d pb_1, Eigen::Vector3d pb_2) {
    Eigen::Vector3d line_a = pa_1.cross(pa_2);
    Eigen::Vector3d line_b = pb_1.cross(pb_2);
    Eigen::Vector3d i = line_a.cross(line_b);
    if (std::abs(i(2)) > 1e-12) {
        i = i/i(2);
    }
    return i;
}

One common suggestion would be to make the function arguments const refs, to avoid copying them unnecessarily:

Eigen::Vector3d intersection(const Eigen::Vector3d& pa_1, const Eigen::Vector3d& pa_2,
                             const Eigen::Vector3d& pb_1, const Eigen::Vector3d& pb_2) {
    Eigen::Vector3d line_a = pa_1.cross(pa_2);
    Eigen::Vector3d line_b = pb_1.cross(pb_2);
    Eigen::Vector3d i = line_a.cross(line_b);
    if (std::abs(i(2)) > 1e-12) {
        i = i/i(2);
    }
    return i;
}

While this is generally good advice in C++, unfortunately it’s not so simple in Eigen. It may provide improved performance when used like this,

Eigen::Vector3d a1, a2, b1, b2;
...
Eigen::Vector3d i = intersection(a1, a2, b1, b2);

However, it provides no performance benefit when the function arguments aren’t Eigen::Vector3d objects to begin with, for instance, when the inputs are columns from a matrix:

Eigen::Matrix3Xd m;
...
Eigen::Vector3d i = intersection(m.col(0), m.col(1), m.col(2), m.col(3));

Enter Eigen::Ref

Because it often makes sense to use the same function on a whole matrix and a block of a matrix, Eigen provides the Eigen::Ref class to represent either. In the above example, the const Eigen::Vector3d& can be replaced by const Eigen::Ref<const Eigen::Vector3d>. The function signature becomes:

Eigen::Vector3d intersection(const Eigen::Ref<const Eigen::Vector3d> pa_1,
                             const Eigen::Ref<const Eigen::Vector3d> pa_2,
                             const Eigen::Ref<const Eigen::Vector3d> pb_1,
                             const Eigen::Ref<const Eigen::Vector3d> pb_2) {

The C++ ref (&) is dropped (although it doesn’t have to be), and the type is wrapped by const Eigen::Ref<>.

Besides the performance benefit of avoiding copies and unnecessary constructions, an Eigen::Ref can be used for output (and in/out) parameters as well, so you can avoid an extra line or two of code to use your own temporary matrix object:

void foo_in_out(Eigen::Ref<Eigen::Matrix3d> m) {
  // rotate the corners
  double t = m(0, 0);
  m(0, 0) = m(2, 0);
  m(2, 0) = m(2, 2);
  m(2, 2) = m(0, 2);
  m(0, 2) = t;
}

...

foo_in_out(m.topLeftCorner<3, 3>());

Benchmarking results

I benchmarked the three implementations of intersection, when called directly on Eigen::Vector3d objects and when called on columns of a 3x4 matrix. The code is here, and here are the results:

2023-11-21T08:50:37-05:00
Running ./eigen_ref
Run on (16 X 5000 MHz CPU s)
CPU Caches:
  L1 Data 48 KiB (x8)
  L1 Instruction 32 KiB (x8)
  L2 Unified 1280 KiB (x8)
  L3 Unified 18432 KiB (x1)
Load Average: 0.20, 0.25, 0.27
***WARNING*** CPU scaling is enabled, the benchmark real time measurements
may be noisy and will incur extra overhead.
----------------------------------------------------------------------------
Benchmark                                  Time             CPU   Iterations
----------------------------------------------------------------------------
BM_vector_intersection_naive           0.602 ns        0.602 ns   1152199989
BM_vector_intersection_const_ref       0.602 ns        0.602 ns   1159881265
BM_vector_intersection_eigen_ref       0.602 ns        0.602 ns   1158563404
BM_mat_col_intersection_naive           6.62 ns         6.62 ns    105187801
BM_mat_col_intersection_const_ref       6.62 ns         6.62 ns    105236839
BM_mat_col_intersection_eigen_ref       3.20 ns         3.20 ns    218533472

Interestingly, all three methods perform the same when called with plain Eigen::Vector3d arguments. Probably the compiler is smart enough to optimize all three down to the same machine code. Running on a matrix column shows the benefit of the Eigen::Ref type, where a C++ const ref argument performs no better than the naive implementation, but the Eigen::Ref case is twice as fast.