当尝试使用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)进行评价并比较结果时，差异要小得多，但仍然存在。如果用几乎相同的确定性方式赋值，为什么这两个向量是不同的?
我在Arch Linux和Debian上用x86-64架构，使用Eigen版本3.4.0,c++ 11, c++ 17, c++ 20，用clang和gcc编译。
官方;Eigen维护者Antonio Sànchez的回答如下:
[…在这种情况下，三角解算器本身稍微有点不同的代码路径:
- 逗号初始化器版本使用了一个路径，其中RHS可以是矩阵
- 赋值版本使用的是优化路径，其中RHS已知是编译时向量
这里的一些操作顺序会引起轻微的变化最终结果应在合理的容忍范围内。这也复制相同的问题:
Eigen::VectorXd xa = R.triangularView<Eigen::Upper>().solve(b);
Eigen::MatrixXd xb = R.triangularView<Eigen::Upper>().solve(b);
这是因为这些代码路径选择是在编译时进行的。求解器做一个就地求解，因此首先将b复制到输出，然后求解。正因为如此，矩阵/向量的确定实际上结束了使用LHS类型。对于逗号初始化式(<<)表达式输入到<<运算符被视为一般的表达式，可以是一个矩阵。
如果你添加。evaluate()它首先将求解结果求值为一个临时编译时向量(因为vector b是编译时vector)，所以得到vector路径。最后，我不认为这是一个bug，我也不会必须称之为"有意的"[…]
它与h . rititch在他们的答案中理论的方向相同:operator<<和operator=只是导致不同的代码路径，从而允许不同的优化。
定义你正在求解的线性系统的矩阵的条件数是10⁷阶。粗略地说，这意味着在用数字解决这个系统后，最后7位数字将是不正确的。因此，给你留下大约9个正确的数字或10 - 9左右的错误。
y = R.triangularView<Eigen::Upper>().solve(b);
x << R.triangularView<Eigen::Upper>().solve(b);
产生的机器码略有不同。因为你的矩阵是无条件的，所以我们预计误差为10(9)。或者换句话说，计算出的解决方案相差10⁻(9)。
您可以使用下面的代码验证行为。如果激活
行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;
}
注意特征很大程度上依赖于函数内联和编译器的优化。对于每次调用solve函数，编译器都会根据上下文生成一个优化的solve函数。因此，operator<<和operator=可能允许不同的优化，从而导致不同的机器码。至少在我的编译器中，如果你计算
x << VectorXd(R.triangularView<Eigen::Upper>().solve(b));
的值x和y同意。
"近same"不是"相同"的。在本例中，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。您可以显式地将Solve与VectorXd(R.triangularView<Eigen::Upper>().solve(b))复制到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