Different behaviours when initializing differently

Question

When trying to initialize a Vector using the result of some operation in Eigen, the result seems to be different depending on what syntax is used, i.e., on my machine, the assertion at the end of the following code fails:

const unsigned int n = 25;
Eigen::MatrixXd R = Eigen::MatrixXd::Random(n,n);
Eigen::VectorXd b = Eigen::VectorXd::Random(n);

Eigen::VectorXd x = Eigen::VectorXd::Zero(n);
Eigen::VectorXd y = Eigen::VectorXd::Zero(n);

y = R.triangularView<Eigen::Upper>().solve(b);
x << R.triangularView<Eigen::Upper>().solve(b);

assert((x-y).norm() < std::numeric_limits<double>::epsilon()*10E6);

I am aware of potential rounding errors, but in my understanding, the two evaluations of R.triangularView<Eigen::Upper>().solve(b) should have the exact same precision errors, and therefore the same result. This also only happens when initializing one variable with <<and the other withoperator=, but not if both variables are assigned to the same way.

When not using only backwards substitution on the upper triangular part but evaluating R.lu().solve(b)on both and comparing the result, the difference is far smaller, but still exists. Why are the two vectors different if assigned in nearly the same, deterministic way?

I tried this code on Arch Linux and Debian with a x86-64 architecture, using Eigen Version 3.4.0, with C++11, C++17, C++20, compiled with both clang and gcc.

Answer 1

The "official" answer by Eigen maintainer Antonio Sànchez is the following:

[...] In this case, the triangular solver itself is taking slightly different code paths:

• the comma-initializer version is using a path where the RHS could be a matrix
• the assignment version is using an optimized path where the RHS is known to be a compile-time vector

Some order of operations here is causing slight variations, though the end results should be within a reasonable tolerance. This also reproduces the same issue:

Eigen::VectorXd xa = R.triangularView<Eigen::Upper>().solve(b);
Eigen::MatrixXd xb = R.triangularView<Eigen::Upper>().solve(b);

https://godbolt.org/z/hf3nE8Prq

This is happening because these code-path selections are made at compile-time. The solver does an in-place solve, so ends up copying b to the output first, then solves. Because of this, the matrix/vector determination actually ends up using the LHS type. In the case of the comma initializer (<<), the expression fed into the << operator is treated as a general expression, which could be a matrix.
If you add .evaluate() to that, it first evaluates the solve into a temporary compile-time vector (since the vector b is a compile-time vector), so we get the vector path again. In the end, I don't think this is a bug, and I wouldn't necessarily call it "intended".[...]

It goes into the same direction as H.Rittich theorizes in their answer: operator<< and operator= simply lead to different code paths, which in turn allow for different optimizations.