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

Question

我从 Eigen 中的线性求解计算中得到了意外和无效的结果。我有一个 vector 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 求逆)，看看它会很有趣。