CODEWITHSUNDEEP
lapack
user402940
user402940

Reputation: 345

compute determinant from LU decomposition in lapack

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.

Upvotes: 8

Views: 7693

Answers (4)

siimurik
siimurik

Reputation: 1

C version:

/*
 Compile and excute with:
    $ gcc det.c -o det -llapacke
    $ ./det 
*/

#include <stdio.h>
#include <stdlib.h>
#include <lapacke.h>
#define dim 3

int main() {
    int n = dim; 
    int lda = n;

    printf("Initial matrix:\n");
    double A[dim*dim] = { 3.0, -1.0,  2.0, 
                         -3.0,  3.0, -1.0, 
                          6.0,  0.0,  4.0 };
    int ipiv[dim];
    int info;

    // Print matrix A
    for (int i = 0; i < n; i++) {
        printf("  [");
        for (int j = 0; j < n; j++) {
            printf("%f", A[i * n + j]);
            if (j < n - 1) printf(", ");
        }
        printf("]");
        if (i < n - 1) printf(",\n");
    }

    // Perform LU-decomposition to get ipiv values
    info = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, A, lda, ipiv);
    if (info < 0) {
        printf("\nError: Argument %d had an illegal value\n", -info);
        return -1;
    } else if (info > 0) {
        printf("\nError: Matrix is singular\n");
        return -1;
    }

    // Calculate the determinant
    double determinant = 1.0;
    for (int i = 0; i < n; i++) {
        determinant *= A[i * n + i];
    }
    int sign = 1;
    for (int i = 0; i < n; i++) {
        if (ipiv[i] != i + 1) {
            sign = -sign;
        }
    }
    determinant *= sign;

    // Print the result
    printf("\nDeterminant: %lf\n", determinant);

    return 0;
}

Upvotes: 0

PhiPsi.top
PhiPsi.top

Reputation: 11

Fortran code:

!LU decomposition
lda = n
call dgetrf( n, n, A,lda, ipiv, info)
    
!det_U
det_U = ONE
do i=1,n
  det_U =  det_U*A(i,i)
enddo

!det_L
det_L = ONE

!det_P
det_P = ONE
do i=1,n
  if (ipiv(i)/=i)then
      det_P = -det_P
  endif
enddo

!det
det = det_P*det_L*det_U

Upvotes: 0

francis
francis

Reputation: 9817

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...

Upvotes: 7

percusse
percusse

Reputation: 3106

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.

Upvotes: 2

Related Questions