C++ Eigen LU 分解 - 行列式的符号

Posted 2023-03-12

【中文标题】C++ Eigen LU 分解 - 行列式的符号

【英文标题】：C++ Eigen LU decomposition - sign of determinant

【发布时间】：2020-11-09 12:55:56

【问题描述】：

我使用 C++ 14 和 Eigen，我想使用 LU 分解来计算方阵的行列式。 (link) 但我有一些问题：A 是主矩阵（n 大小），rez 是 A 的 LU 形式。

PartialPivLU<MatrixXd> rez = PartialPivLU<MatrixXd>(A);
MatrixXd r1 = rez.matrixLU().triangularView<UpLoType::Upper>();
double det_r1 = 1;
for(int i=0;i<n;i++)
    det_r1 = det_r1 * r1(i,i);
cout<<det_r1<<endl;
cout<<A.determinant();

r1 的行列式是主对角线元素的乘积。问题是 det(A) 不等于 det(r1)。

例如 det(A) = 500 和 det(r1) = -500。问题是关于 r1 的符号，我怎样才能得到符号？

【问题讨论】：

rez.determinant() 呢？

我需要使用主对角线元素的乘积来计算行列式

由于您可以访问排列的索引，因此您可以计算此排列是奇数还是偶数

@dacian 这就是 Eigen 在内部做的事情。但它也会跟踪排列的数量，因此也能正确计算符号（只需检查源代码以了解发生了什么）。

计算排列的奇偶性需要一点努力，例如查看循环的长度。你可以在这里找到一些信息：en.wikipedia.org/wiki/Parity_of_a_permutation

【参考方案1】：

您链接到的文档说这是 LU 分解部分旋转。这意味着您将 A 分解为 PLU，其中 P 是置换矩阵。它的行列式是排列的符号。

从您链接到的文档中：

该类表示平方可逆的 LU 分解 矩阵，部分旋转：矩阵 A 分解为 A = PLU 其中 L 是单位下三角形，U 是上三角形，P 是 置换矩阵。

在您的情况下，置换矩阵必须表示奇数置换，因此行列式为 -1。

【讨论】：

如何计算 P 的行列式？因为 P 是排列（类），而不是矩阵。或者，我怎样才能得到排列的符号？ @gspr

如果您真的想要行列式，为什么不直接使用PartialPivLU 中的determinant 函数？