Question 47097 in Debian Question posted in
The oficial documentation can be found here.

Different behaviours when initializing differently – Debian

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.

Tags:

3

Answers


  1. 0

    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.

    Login or Signup to reply.
  2. 0

    "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
    Login or Signup to reply.
  3. 0

    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.

    Login or Signup to reply.
Please signup or login to give your own answer.

Back To Top
Search