不同的初始化方式会导致不同的行为。

Question

尝试使用Eigen中的某些操作结果来初始化一个向量时，似乎根据使用的语法不同而产生不同的结果。例如，在我的机器上，以下代码结束处的断言失败：

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=时才会出现这种情况，如果两个变量都以相同的方式分配，则不会出现这种情况。

当不仅在上三角部分使用向后代换，而是在R.lu().solve(b)两者上进行评估并比较结果时，差异要小得多，但仍然存在差异。为什么两个向量在几乎相同的确定性方式下被分配时会有所不同？

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

- pfisch

---

哦，<<行甚至无法编译。您使用的Eigen版本和头文件是什么？ - Swift - Friday Pie

我的头文件是#include <eigen3/Eigen/Dense> #include <cassert>，我使用的是Eigen 3.4.0。 - pfisch

当我运行你的示例时，x和y之间的差异在e-14的数量级上。范数大约为2e-12，断言通过是因为eps*10e6等于2e-9。更改c++版本不会影响结果。此外，执行x << y产生了预期的0.0差异。执行x = R.triangular...也产生了0.0，因此我可以确认插入与求解确实导致了差异。 - Matt

我对此感到相当困惑。@Matt 我之前也尝试过计算并使用<<存储值，但如果我将Expression模板传递给operator<<，问题似乎只会出现在这种情况下。我认为这应该等同于使用operator=，特别是因为Eigen在其文档中使用了类似的代码，也将表达式模板传递给该运算符。我不太明白结果如何可能在数值上有所不同。 - pfisch

请在以后的提问中学习如何创建一个适当的 [mre] 来展示给我们。同时，请花些时间阅读 帮助页面，参观 SO [导览], 阅读 [提问指南]，以及 这个问题清单。 - Some programmer dude

Answer 1

由 Eigen 维护者 Antonio Sànchez 给出的 "官方"答案 是:

[...] 在这种情况下，三角求解器本身正在采用略微不同的代码路径:

• 逗号初始化版本使用的是 RHS 可能是矩阵的路径
• 赋值版本使用的是已知为编译时向量的优化路径

某些运算顺序在这里会导致轻微变化，但最终结果应在合理的容差范围内。 这也会产生相同的问题：

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

https://godbolt.org/z/hf3nE8Prq

这是因为这些代码路径选择是在编译时进行的。求解器进行了原地求解，所以首先将b复制到输出中，然后进行求解。由于这个原因，矩阵/向量确定实际上使用了LHS类型。在逗号初始化程序(<<)的情况下，输入到<<运算符中的表达式被视为一般表达式，可能是一个矩阵。
如果你添加.evaluate()，它首先将求解评估为一个临时的编译时向量(因为向量b是一个编译时向量)，所以我们再次得到了向量路径。最后，我不认为这是一个错误，我也不会把它称为"有意的"[...]

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

- Teee

---

此外，他还写道这是一个错误，因为Block typedef应该将编译时已知的大小传播到Block模板，但它没有这样做。 - pfisch

Answer 2

您正在解决的线性系统所定义的矩阵条件数为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的值相同。

- H. Rittich

---

谢谢您的回答，但它并没有完全回答我的问题。我完全同意您的观点，在该范围内出现误差是可以预料的，并且可以加以缓解。让我困惑的是，operator<<和operator=产生不同的数值误差，因为我原本认为R.triangularView<Eigen::Upper>().solve(b)的两个计算过程都是对同一例程的调用，因此在给定相同输入时应具有相同的舍入误差。为什么这个例程会是非确定性的呢？ - pfisch

Eigen在很大程度上依赖于内联和编译器优化。函数不仅有一个调用，而是根据上下文为每个solve函数调用生成优化变体。因此，operator<<和operator=可能允许不同的优化，从而导致不同的机器代码。但两者都是确定性的。 - H. Rittich

Answer 3

“几乎相同”并不等同于“完全相同”。在这种情况下，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

- Matt