不同初始化时的不同行为

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

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

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

Answer 1

Eigen 维护者 Antonio S\xc3\xa0nchez 的“官方”答案如下：

\n

\n

[...] 在这种情况下，三角求解器本身采用稍微\n不同的代码路径：

\n

\n
• 逗号初始值设定项版本使用的路径的 RHS 可以是\n矩阵
\n
• 分配版本使用优化路径，其中 RHS\nis 已知为编译时向量
\n

\n

尽管最终结果应该在合理的容差范围内，但此处的某些操作顺序会导致轻微的变化。这也重现了相同的问题：

\n

\n

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

\n

\n

https://godbolt.org/z/hf3nE8Prq

\n

\n

\n

发生这种情况是因为这些代码路径选择是在编译时进行的。求解器执行\n就地求解，因此最终首先将 b 复制到输出，然后求解。\n因此，矩阵/向量确定实际上\n使用 LHS 类型。对于逗号初始值设定项 (<<)，输入到 << 运算符的\n表达式被视为一般\n表达式，它可以是一个矩阵。
\n如果将 .evaluate() 添加到其中，\n首先将求解计算为临时编译时向量\n（因为向量 b 是编译时向量），因此我们再次获得向量\n路径。最后，我不认为这是一个错误，而且我不会\n必然将其称为“有意的”。[...]

\n

\n

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

\n