compute determinant from LU decomposition in lapack

Question

Lapack, most probably, doesn't have any routine for computing determinant. However we can compute it using either LU, QR or SVD decomposition. I prefer to use LU decomposition. Now lapack uses some dgetrf subroutine to factorize a matrix A into PLU format with some IPIV array. I don't have much idea how to deal with this information. To compute the determinant, I just multiply diagonal elements of U matrix. But what is L and U in PLU format and how to extract them. I am programming in C.

Answer 1

Lapack's dgetrf() computes a A=P*L*U decomposition for a general M-by-N matrix A. Assuming an invertible square matrix A, its determinant can be computed as a product:

• U is an upper triangular matrix. Hence, its determinant is the product of the diagonal elements, which happens to be the diagonal elements of the output A. Indeed, see 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 is a lower triangular matrix featuring unit diagonal elements which are not stored. Hence, its determinant is always 1.

• P is a permutation matrix coded as a product of transpositions( i.e. 2-cycles or swap) . Indeed, see dgetri() to understand how it is used. Hence, its determinant is either 1 or -1, depending on whether the number of transpositions is even or odd. As a result, the determinant of P can be computed as:

  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;
      }
  }
  

The complexity of this method is dominated by the cost of the Gaussian elimination using partial pivoting, that is O(2/3n^3).

You will likely turn to full pivoting using dgetc2() or QR decomposition to improve the accuracy. As signaled in Algebraic and Numerical Techniques for the Computation of Matrix Determinants by Pan et. al., combining equation 4.8, 4.9 and proposition 4.1, the final error on the determinant likely scales like ed=(a+eps*a*n^4)^{n-1}*eps*an^5=a^n*(1+eps*n^4)^{n-1}*n^5*eps where eps is the precision of double (about 1e-13), a is the largest magnitude of all elements in matrix A and n is the size of the matrix. It means that the computed determinant is not very significant for "large" matrices: see tables, it particular the relative error as PLU decomposition is used! The article also provides an algorithm to track the propagation of error and produce a better estimate of the error.

You can also try the Faddeev–Le Verrier algorithm...

Answer 2

L is a unit diagonal matrix so its determinant is always unity. For P it is either 1 or -1 depending on the odd/even number of permutations but ipiv is a swap array not the permutation array so you can convert the following simple Python loop to C

det = 1

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

But note that multiplication can both overflow or underflow if the conditioning of the matrix is problematic in terms of conditioning. Instead just use SVD.