C++ Eigen 线性系统求解，数值问题？

Question

我从 Eigen 中的线性求解计算中得到了意想不到的无效结果。我有一个向量 x 和矩阵 P.

在 MATLAB 符号中，我这样做：

x'/P*x

（类似于x'*inv(P)*x，但没有直接反转）

这是二次形式，应该会产生正结果。矩阵 P 是正定的并且不是病态的。在 MATLAB 中，结果是正的，但很大。

在 C++ Eigen 中，我的斜杠运算符是这样的：

x/P 实现为：

(P.transpose()).ColPivHouseholderQr().solve(x.transpose).transpose()

大体上给出了正确答案，但似乎比 MATLAB 更脆弱。在我现在看到的情况下，它为 x'/P*x 提供了一个非常大的负数，就好像有溢出和环绕或类似的东西。

有没有办法对其进行碎片整理？

编辑：做一些实验我发现 Eigen 也无法反转这个矩阵，而 MATLAB 没有问题。矩阵条件数为3e7，还不错。 Eigen 通过线性求解和简单的 x.transpose() * P.inverse() * x 给出了相同（错误）的答案。怎么回事？

Answer 1

以下是错误的（除了你在问题中遗漏的()）：

(P.transpose()).ColPivHouseholderQr().solve(x.transpose()).transpose();

因为

x'*inv(P) = ((x'    *inv(P))')'
          = (inv(P)'* (x')'  )'
          = (inv(P) *  x     )'  % Note: P=P'

在没有 -DNDEBUG 的情况下构建时，您的表达式实际上应该在 Eigen 中提出断言，因为 x.transpose() 只有一行，但 P 有（通常）更多行。

要计算对称正定 P 的 x'*inv(P)*x，我建议使用 Cholesky 分解 L*L'=P，它给出 x'*inv(P)*x = (L\x)'*(L\x)，在 Eigen 中是：

typedef Eigen::LLT<Eigen::MatrixXd> Chol;
Chol chol(P); // Can be reused if x changes but P stays the same
// Handle non positive definite covariance somehow:
if(chol.info()!=Eigen::Success) throw "decomposition failed!";
const Chol::Traits::MatrixL& L = chol.matrixL();
double quadform = (L.solve(x)).squaredNorm();

对于矩阵 Eigen 求逆失败（Matlab 确实求逆），看看它会很有趣。

C++ Eigen 线性系统求解，数值问题？

C++ Eigen linear system solve, numerical issues?

c++

matlab

eigen