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

Question

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

Answer 1

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

• U是上三角矩阵。因此，它的行列式是对角元素的乘积，恰好是输出的对角元素A。事实上，看看输出 A 是如何定义的：

  退出时，因子 L 和 U 来自因式分解A = P*L*U；不存储 L 的单位对角线元素。

• 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)。

您可能会转向使用完全旋转dgetc2()或 QR 分解来提高准确性。正如Pan 等人在用于计算矩阵行列式的代数和数值技术中所指出的那样。人，结合方程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是的双（约1E-13）的精度，一个是在基质A中的所有元素的最大幅值和n是大小的矩阵。这意味着计算出的行列式对于“大”矩阵不是很重要：参见表格，它特别是使用 PLU 分解时的相对误差！文章还提供了一种算法来跟踪错误的传播并产生更好的错误估计。

您还可以尝试Faddeev–Le Verrier 算法...