# How to treat numerical errors in determinants of singular matrices when using LU decomposition

I want to calculate the determinant of a matrix. Currently I'm using LU decomposition.

To check my algorithm I wrote a unit test with random matrices. In one part I set one row to be equal to another row to test if the resulting determinant would be 0. This is the matrix:

$\begin{pmatrix} 0 & 0 & 33 & 17 & 42 & 44 & 20 & 28 & 43 \\ 38 & 46 & 11 & 39 & 21 & 24 & 24 & 29 & 4 \\ 52 & 13 & 43 & 3 & 54 & 31 & 32 & 50 & 32 \\ 2 & 10 & 41 & 12 & 10 & 28 & 32 & 14 & 16 \\ 22 & 34 & 44 & 52 & 4 & 36 & 8 & 44 & 44 \\ 20 & 0 & 19 & 11 & 52 & 19 & 54 & 43 & 18 \\ 17 & 20 & 13 & 49 & 9 & 23 & 22 & 21 & 20 \\ 52 & 13 & 43 & 3 & 54 & 31 & 32 & 50 & 32 \\ 11 & 14 & 37 & 43 & 21 & 24 & 49 & 32 & 8 \end{pmatrix}$

As can be seen the third and 8th row are equal. By reusing code from a Gauß-Jordan inversion I basically can get something alike to an LU-decomposition, i.e. I rewrite $A\vec{b}=I$ to $U\vec{b}=L$, so that $A=L^{-1}U$. The determinant can then be written as a product of the diagonal elements of L. It basically boils down to:

$|A| = 38 \cdot (-49.95) \cdot 33 \cdot (-20.7) \cdot (-58.6) \cdot (-14.61) \cdot 44.61 \cdot \left(-4.8 \cdot 10^{-15}\right) \cdot (-11.33) = 0.0027$

when using double precision. I expected a relative error of roughly DBL_EPSILON * 9, that's why my unit test failed. For single precision floating point I would get approx. 10000 instead of 0!

So the problem being here is that the error rises not additively, but multiplicatively for values near zero, so that even though the each diagonal element has roughly 10*DBL_EPSILON (-4.8e-15 instead of 0) relative error, the error on the calculated determinant becomes 1e12*DBL_EPSILON.

How can this be treated? Are there more stable algorithms for which I can get better precision on $|A|$, or should I just hard-code manually check whether one of the factors is near 0? Is such a check safe, or would it dirty other cases?

• The LUP decomposition is better if I recall correctly. – xavierm02 Feb 9 '17 at 13:17
• @xavierm02 If you mean using row permutations to avoid 1/0, then I'm already using that, didn't notice it was called differently. The problem still is the same. – mxmlnkn Feb 9 '17 at 13:50
• You don't just want to avoid dividing by $0$: "But we also want to avoid dividing by a small value - even if $A$ is nonsingular - because numerical instabilities can result. We therefore try to pivot on a large value." [Cormen, 3rd ed., p. 823] – xavierm02 Feb 9 '17 at 14:09
• @xavierm02 So I now do a maximum search for the pivot element and I still have this problem. I hope I'm doing it right. I could post the code, but it's a bit complex and lengthy, because it is in-place. But I still think it's a problem not dependent on my implementation. – mxmlnkn Feb 10 '17 at 0:43
• Is it out of question to compute it exactly using rationnal numbers or lists of floats? Also, since you're on integers, you could try computing it modulo several numbers, and then guess the value using the Chinese remainder theorem. I know this works well when trying to find something that can then easily be checked (like the solution to a linear system) but I'm not sure how to check the determinant... – xavierm02 Feb 10 '17 at 11:00