Calculating determinant from LUDecomposition

As far as i know LU decomposition allow you to calculate matrix determinant following easy and cheap formula: Det[A] = Det[L] Det[U] = Det[U]

Trying this out in Mathematica 7 gives me correct result by absolute value, i.e. it ignores negative determinants and transforms them to positive ones.

Sample code:

matrixA = {{-1, 2, 3, 5}, {-7, -4, 5, 4}, {-89, 7, 8, -6}, {8, 6, -1, 4}};


Det[matrixA] gives out -2067

but

{lu, p, c} = LUDecomposition[matrixA]

u = lu SparseArray[{i_, j_} /; j >= i -> 1, {4, 4}]


Det[u] is 2067

Well, the question is obvious - how to get correct determinant in Mathematica using LU decomposition?

Well, this is because you forgot to take into account the permutation matrix that is output by the LU decomposition routine. You have to remember two facts: 1.) Gaussian elimination generally performs row interchanges for numerical stability reasons, and 2.) the determinant of a matrix changes sign if you interchange rows.

In Mathematica at least, the function that will rescue you is Signature[], which gives the signature of the permutation required to turn the permutation matrix (which Mathematica outputs as a scrambled list of the numbers from 1 to n, where n is the size of your matrix) into the identity matrix.

So, to use your example in Mathematica:

m={{-1, 2, 3, 5}, {-7, -4, 5, 4}, {-89, 7, 8, -6}, {8, 6, -1, 4}};


we compare

Det[m]
-2067


with the following:

{lu, p, c} = LUDecomposition[m]
{{{-1, 2, 3, 5}, {7, -18, -16, -31}, {-8, -11/9, 31/9, 55/9}, {89, 19/2, -963/31, 2067/62}}, {1, 2, 4, 3}, 1}


Now, LUDecomposition[] here outputs three things: the merged $\mathbf{L}$ and $\mathbf{U}$ matrices, the permutation, and the condition number. We can get the tentative determinant by multiplying together the diagonal elements of lu, thus:

Tr[lu, Times]
2067


Here is where Signature comes in:

Signature[{1, 2, 4, 3}]
-1


Note that to turn {1, 2, 4, 3} into {1, 2, 3, 4}, one needs one swap; namely, to swap the third and fourth elements. Since 1 is odd, the signature is -1.

Thus,

Signature[p]*Tr[lu, Times]
-2067


From the documentation: "Therefore multiplying L and U recreates the original matrix permuted with the row pivot vector." So you need to multiply by Signature[p] to get the correct sign.