# Determinant of an ill conditioned matrix

I have the following ill conditioned matrix. I want to find its determinant. How is it possible to calculate it without much error $$\left[\begin{array}{cccccc} 4.6169e90&1.0383e69&0&0&0&0\\ 1.5357e82&1.6153e65&1.2641e-100&4.4193e-109&2.5510e-128&1.5120e-131\\ 8.9950e76&1.3138e60&1.3720e-060&6.6491e-063&1.0853e-066&4.1555e-067\\ 1.7734e75&2.5766e58&2.7499e-063&1.0104e-065&8.0964e-070&2.7419e-070\\ 3.4969e73&5.0551e56&5.9760e-066&1.6975e-068&7.0870e-073&2.1437e-073\\ 6.8969e71&9.9210e54&1.3969e-068&3.1164e-071&7.1027e-076&1.9340e-076 \end{array}\right]$$ Let us call this matrix $A$. I want to find the determinant of $3.8233e17\times A$. MATLAB gives it as $4.5836e-013$. But is there a better way to do it in more accurate way. One information that may be helpful in answering this question is that I know the accurate $\log$ values of each element of the matrix. For e.g., I know the accurate log value of $A(1,1) = 4.6169e90$ and similarly for each element.

• Maybe this is a stupid comment because I don't know very much about this. But, anyway, have you tried to expand the determinant on the first row?
– mfl
Jul 7, 2014 at 23:54
• The calculation of the determinant in my problem is part of a bigger problem for which I am running a MATLAB code. I have to calculate the determinant several times at each iteration. So, it is not possible to do it with hand. Jul 8, 2014 at 0:01
• I refer to expand it using MATLAB. What is the result if you compute the two determinants instead of only one?
– mfl
Jul 8, 2014 at 0:09

This problem is difficult for numerical rather than computational reasons. Part of the problem is that you really need to be confident that the matrix is full rank, because if it is not, then a single error can make a determinant very large when it should actually be zero. For illustration, suppose we were trying to compute the determinant of

$$A =\begin{bmatrix} M & M \\ M & M \end{bmatrix}$$

where $M$ is very large. This determinant is of course $0$. Let's say some roundoff gave us

$$B = \begin{bmatrix} M & M \\ M & (M+\varepsilon) \end{bmatrix}$$

for some small $\varepsilon$. Now $A$ has determinant $0$ but $B$ has determinant $M\varepsilon$, which may be quite large. The coefficient gets much larger in higher dimensions (though it is always first order, as the determinant is a polynomial in the entries).

If you are confident that the matrix is full rank, then my best suggestion would be to perform an SVD, check to see that all the singular values are nonzero, then if they are not, do it again in higher precision.

Edit: there is one more thing you can do. Because $\text{det}(A)=\text{det}(A^T)$, you can perform column operations. In this case what you should do is rescale the matrix so that the largest entry in each column is $1$. You will multiply column $x_i$ by a number $c_i$, which will also multiply the determinant by $c_i$. Accordingly you will want to divide the final result by $c_i$, so that you get the determinant of the matrix you started with. You will still find that the matrix is extremely ill-conditioned afterward; for example you will still have a column with one entry of order $1$ and another of order $10^{-40}$. But the conflict between the first two columns and the rest will be gone.

I think the way to deal with such a problem is either using arbitrary-precision real field implementations, e.g. Multiple Precision Toolbox for MATLAB, RealField in Sage, mpmath for Python,

or an interval-arithemtic library, such as INTLAB for Matlab, see for the INTLAB program examples Computer-assisted Proofs and Self-validating Methods and THE DETERMINANT OF AN INTERVAL MATRIX USING GAUSSIAN ELIMINATION METHOD.

I know pari/GP calculator is not matlab or octave, but this might be better for what you need:

{
M = [4.6169e90, 1.0383e69, 0, 0, 0, 0;
1.5357e82, 1.6153e65, 1.2641e-100, 4.4193e-109, 2.5510e-128, 1.5120e-131;
8.9950e76, 1.3138e60, 1.3720e-060, 6.6491e-063, 1.0853e-066, 4.1555e-067;
1.7734e75, 2.5766e58, 2.7499e-063, 1.0104e-065, 8.0964e-070, 2.7419e-070;
3.4969e73, 5.0551e56, 5.9760e-066, 1.6975e-068, 7.0870e-073, 2.1437e-073;
6.8969e71, 9.9210e54, 1.3969e-068, 3.1164e-071, 7.1027e-076, 1.9340e-076];

M = M * 3.8233e17;
detM = matdet(M);
printf("%4d %20.15g\n", default(realprecision), detM);
}

pari/GP will use 38 digits by default. When I run this program using gp, this is what I get:

38 4.58945521135206 e-13

You can change the default precision by using "default(realprecision,100)" for example, to make it work with at least 100 digits. Multiprecision is built into the calculator from the start, and you don't need to hope that a "toolbox" is installed or not, or will work with your program.