Skip to content

CommaInitializer and solve lead to rounding errors

Summary

Assignment using << and = do not yield the same results, when temporary expression objects (in this case invoking a solver) are involved.

While I don't know if this is indeed a bug, I find it very counter-intuitive. We've discovered this behaviour in a course with ~50 students at ETH Zurich, so I can confirm that it is not limited to my machine.

Environment

  • Operating System : Linux, Windows, and MacOS
  • Architecture : x64
  • Eigen Version : 3.4.0-1 (likely also earlier versions)
  • Compiler Version : GCC 11.1.0, Clang 12.0.1 (also earlier compilers)
  • Compile Flags : no specific ones
  • Vector Extension : SSE, SSe2, AVX, ... we tested it across many CPUs

Minimal Example

int main() {
  // system dimension
  const unsigned int n = 9;
  // random test vectors
  Eigen::MatrixXd R = Eigen::MatrixXd::Random(n, n);
  Eigen::VectorXd b = Eigen::VectorXd::Random(n);
  Eigen::VectorXd xa(n), xb(n);
  
  //----------------------------------------------------------------
  xa << R.triangularView<Eigen::Upper>().solve(b);
  xb = R.triangularView<Eigen::Upper>().solve(b);
  //----------------------------------------------------------------


  std::cout << "LSE solution A: " << xa.transpose() << std::endl;
  std::cout << "LSE solution B: " << xb.transpose() << std::endl;
  std::cout << "Error compared between each other: " << (xa - xb).norm() << std::endl;

}

(you can swap out the randomized matrix and vector for static ones; it doesn't change the result).

Steps to reproduce

  1. Compile (using plain gcc file.cc is enough) and run the example.

What is the current bug behavior?

(xa - xb).norm() != 0

What is the expected correct behavior?

We expect that both assignments yield the same result.

Anything else that might help

Including but not limited to the following: The values are equal (i.e. the norm of the difference is approx 0) if we either store the result in a variable:

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

or force the evaluation before <<-assigning it:

xa << R.triangularView<Eigen::Upper>.solve(b).evaluate();

Also, the issue only comes up in combination with solve.

If there is anything else I can provide or help with, just @ me :).

In case this is intended / known behaviour, it might be good to document it somewhere.

Edited by Stresspresso