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.

Dish

Hmm, << line doesn't even compile. What version of Eigen do you use and what headers? – Jarodjarosite

My headers are #include <eigen3/Eigen/Dense> #include <cassert>, I use Eigen 3.4.0. – Dish

When I run you're example the differences between x and y are on the order of e-14. The norm is around 2e-12 and the assert passes because eps*10e6 is 2e-9. Changing versions of c++ didn't affect the results. Also, doing x << y produce a difference of 0.0 as expected. Doing x = R.triangular... also produced 0.0 so I can confirm the insertion with solve does cause a difference. – Tad

Hmm, looks like implementation issue, on Eigen 3.2 the use of temporary from solve() and << doesn't compile by design. Apparently insertion acts upon result element by element. – Jarodjarosite

I'm quite confused by this. @Tad I also tried calculating it before and storing the values using << afterwards, and the problem seems to only occur if I pass the Expression template to operator<<. I assumed this would be equivalent to using operator=, especially since Eigen uses similiar code on their documentation, that also passes an expression template to that operator. I don't quite understand how the result could potentially differ numerically. – Dish

For future questions please learn how to create a proper minimal reproducible example to show us. Also please take some time to read the help pages, take the SO tour, read How to Ask, as well as this question checklist. – Acrogen

Not familiar with Eigen, but do note epsilon could potentially be defined differently on different implementations – Peso

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.

Eleonoreleonora

Adding to that, he also writes that this is a bug, since the Block typedef should propagate sizes known at compile time to the Block template, which it doesn't. – Dish

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.

Navvy

Thank you for your answer, but it doesn't quite answer my question. I fully agree with you that errors from the exact result in that range are to be expected, and that they can be mitigated. What confuses me is that operator<< and operator= produce different numerical errors, as I would have thought that the two evaluations of R.triangularView<Eigen::Upper>().solve(b) are calls to the same routine and should therefore have the same rounding errors when given the same input. Why would this routine be non-deterministic? – Dish

@Dish Eigen heavily relies on inlining and compiler optimization. There is not just one call to a function but the compiler generates optimized variants for each call to the solve function depending on the context. Hence, operator<< and operator= might allow for different optimizations and hence lead to different machine codes. Both of them deterministic, though. – Navvy

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

Tad