我从 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
给出了相同(错误)的答案。怎么回事?
最佳答案
以下是错误的(除了您在问题中遗漏的 ()
):
(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 线性系统求解,数值问题?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/51971562/