从 lapack 中的 LU 分解计算行列式

Question

Lapack 很可能没有任何计算行列式的例程。然而，我们可以使用 LU、QR 或 SVD 分解来计算它。我更喜欢使用 LU 分解。现在 lapack 使用一些 dgetrf 子程序将矩阵 A 分解为带有一些 IPIV 数组的 PLU 格式。我不知道如何处理这些信息。为了计算行列式，我只需乘以 U 矩阵的对角线元素。但是PLU格式的L和U是什么以及如何提取。我正在用 C 编程。

Answer 1

L 是一个单位对角矩阵，所以它的行列式总是单位。对于 P，它是 1 或 -1，具体取决于 odd/even 排列数，但 ipiv 是交换数组而不是排列数组，因此您可以将按照简单的 Python 循环到 C

det = 1

for i in range(len(ipiv)):
    det *= u(i,i)
    if ipiv[i] != i:
        det *= -1

但请注意，如果矩阵的条件在条件方面存在问题，则乘法可能会溢出或下溢。相反，只需使用 SVD。

Answer 2

Lapack 的dgetrf() 计算一般 M×N 矩阵 A 的 A=P*L*U 分解。假设一个可逆方阵 A，其行列式可以计算为乘积：

U是上三角矩阵。因此，它的行列式是对角元素的乘积，恰好是输出的对角元素A。确实，参见 how the output A is defined:

On exit, the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored.
L 是一个下三角矩阵，具有未存储的单位对角线元素。因此，它的行列式总是 1.
P 是编码为转置（即 2 循环或交换）乘积的置换矩阵。实际上，请参阅 dgetri() 以了解它的使用方式。因此，它的行列式是 1 或 -1，这取决于转置数是偶数还是奇数。因此，P 的行列式可以计算为：
```
int j;
double detp=1.;
for( j=0;j<n;j++){
    if(j+1!=ipiv[j]){
        // j+1 : following feedback of ead : ipiv is from Fortran, hence starts at 1.
        // hey ! This is a transpose !
        detp=-detp;
    }
}
```

该方法的复杂度主要由使用部分旋转的高斯消元的成本决定，即 O(2/3n^3)。

您可能会使用 Pan 等人的 dgetc2() or QR decomposition to improve the accuracy. As signaled in Algebraic and Numerical Techniques for the Computation of Matrix Determinants 转向完全旋转。等，结合等式 4.8、4.9 和命题 4.1，行列式可能尺度上的最终误差如 ed=(a+eps*a*n^4)^{n-1}*eps*an^5=a^n*(1+eps*n^4)^{n-1}*n^5*eps 其中 eps 是 double 的精度（大约 1e-13），a 是最大幅度的矩阵 A 中的所有元素，n 是矩阵的大小。这意味着计算的行列式对于 "large" 矩阵不是很重要：见表，它特别是使用 PLU 分解时的相对误差！本文还提供了一种算法来跟踪误差的传播并产生更好的误差估计。

您也可以尝试 Faddeev–Le Verrier algorithm...

从 lapack 中的 LU 分解计算行列式

compute determinant from LU decomposition in lapack

lapack