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.

Answer 2

The condition number of the matrix that defines the linear system you are solving is at the order of 10⁷. Roughly speaking, this means that after solving this system numerically the last 7 digits will be incorrect. Thus, leaving you with roughly 9 correct digits or an error of around 10⁻⁹. It seems like

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

produce slightly different machine codes. Since your matrix is that illconditioned we expect an error of the order of 10⁻⁹. Or in other words, that the computed solutions differ by around 10⁻⁹.

You can verify the behavior using the code below. If you activate the line

R += 10*MatrixXd::Identity(n,n);

you decrease the condition number of the matrix, by adding a diagonal term, and hence the error is significantly reduced.

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/SVD>

using Eigen::MatrixXd;
using Eigen::VectorXd;
using Eigen::BDCSVD;

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

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

  // Uncomment to reduce the condition number
  // R += 10*MatrixXd::Identity(n,n);

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

  std::cout << "res(x): " << (b - R.triangularView<Eigen::Upper>() * x).norm() << std::endl;
  std::cout << "res(y): " << (b - R.triangularView<Eigen::Upper>() * y).norm() << std::endl;

  Eigen::BDCSVD<Eigen::MatrixXd> svd(R.triangularView<Eigen::Upper>());
  VectorXd sigma = svd.singularValues();
  std::cout << "cond: " << sigma.maxCoeff() / sigma.minCoeff() << std::endl;

  std::cout << "error norm: " << (x-y).norm() << std::endl;
  std::cout << "error max: " << (x-y).lpNorm<Eigen::Infinity>() << std::endl;

  return 0;
}

Note that Eigen heavily relies on function inlining and compiler optimization. For each call to the solve function the compiler generates an optimized solve function depending on the context. Hence, operator<< and operator= might allow for different optimizations and hence lead to different machine codes. At least with my compiler, if you compute

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

the values for x and y agree.

Answer 3

"Nearly the same" is not "the same". In this case y = R.triangularView<Eigen::Upper>().solve(b); uses the assign_op while x << R.triangularView<Eigen::Upper>().solve(b); uses CommaInitializer.

The CommaInitializer uses the block method to copy the data while the assign_op appears to do the alignment and other things to establish the Eigen::Matrix.

The type of R.triangularView<Eigen::Upper>().solve(b) is Eigen::Solve, not a VectorXd. You can explicitly copy the Solve into the VectorXd with VectorXd(R.triangularView<Eigen::Upper>().solve(b)) or more simply just use auto and let the compiler figure it out. For example:

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);
auto x = R.triangularView<Eigen::Upper>().solve(b); // x is an Eigen::Solve

// assert((x-y).norm() < std::numeric_limits<double>::epsilon()*10E6);
std::cout << (x-y).norm() << std::endl; // will print 0