初始化不同时的不同行为

Question

当尝试使用 Eigen 中某些操作的结果初始化 Vector 时，结果似乎因使用的语法而异，即在我的机器上，以下代码末尾的断言失败：

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);

我知道可能存在舍入错误，但据我所知，两个评估 R.triangularView<Eigen::Upper>().solve(b) 应该具有完全相同的精度误差，因此结果相同。这也只会在用 << 初始化一个变量而用 operator= 初始化另一个变量时发生，但如果两个变量都以相同的方式分配则不会发生。

当不只对上三角部分使用反向替换而是对两者都进行评估并比较结果时，差异要小得多，但仍然存在。如果以几乎相同的确定性方式分配，为什么这两个向量不同？

我在 Arch Linux 和具有 x86-64 架构的 Debian 上尝试了这段代码，使用 Eigen 版本 3.4.0，C++11、C++17、C++20，编译为clang 和 gcc。

Answer 1

定义您正在求解的线性系统的矩阵的条件数约为 10⁷。粗略地说，这意味着在对这个系统进行数值求解后，最后 7 位数字将不正确。因此，给你留下大约 9 个正确的数字或大约 10⁻⁹ 的错误。好像

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

生成略有不同的机器代码。由于您的矩阵是病态的，我们预计误差约为 10⁻⁹。或者换句话说，计算出的解决方案相差大约 10⁻⁹。

您可以使用以下代码验证行为。如果您激活行

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

您通过添加对角线项减少了矩阵的条件数，因此误差显着减少。

#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;
}

请注意，Eigen 严重依赖函数内联和编译器优化。对于每次调用 solve 函数，编译器都会根据上下文生成优化的 solve 函数。因此，operator<< 和 operator= 可能允许不同的优化，从而导致不同的机器代码。至少对于我的编译器，如果你计算

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

x 和 y 的值一致。

Answer 2

“差不多”不是“一样”。在这种情况下 y = R.triangularView<Eigen::Upper>().solve(b); 使用 assign_op 而 x << R.triangularView<Eigen::Upper>().solve(b); 使用 CommaInitializer.

CommaInitializer 使用 block 方法复制数据，而 assign_op 似乎进行对齐和其他操作以建立 Eigen::Matrix.

R.triangularView<Eigen::Upper>().solve(b) 的类型是 Eigen::Solve，而不是 VectorXd。您可以使用 VectorXd(R.triangularView<Eigen::Upper>().solve(b)) 显式地将 Solve 复制到 VectorXd 中，或者更简单地使用 auto 并让编译器计算出来。例如：

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

Answer 3

Eigen 维护者 Antonio Sànchez 的 "official" answer 如下：

[...] 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".[...]

它与 H.Rittich 在他们的回答中的理论化方向相同：operator<< 和 operator= 只是导致不同的代码路径，这反过来又允许不同的优化。

初始化不同时的不同行为

Different behaviours when initializing differently

c++

eigen

eigen3