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

Question

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

Answer 1

拉派克 dgetrf() 计算 A=P*L*U一般 M×N 矩阵 A 的分解。假设可逆方阵 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-cycles 或 swap)。确实，见 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)。

您可能会使用 dgetc2() 转向完全旋转或 QR 分解以提高准确性。如 Algebraic and Numerical Techniques for the Computation of Matrix Determinants 中所述通过潘等。 al.，结合方程 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 是矩阵的大小。这意味着计算出的行列式对于“大”矩阵不是很重要:参见表格，它特别是使用 PLU 分解时的相对误差!文章还提供了一种算法来跟踪错误的传播并产生更好的错误估计。

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