Why is the determinant of this matrix unstable?

Question

A colleague has encountered an instability when computing the determinant of a matrix.

For this matrix : $$\begin{bmatrix} 123.585 &-1.33458 &78.9871 &43.2638\\ -1.33458 &190.581 &88.7218 &100.524\\ 78.9871 &88.7218 &167.465 &0.244369\\ 43.2638 &100.524 &0.244369 &143.544 \end{bmatrix}$$ the determinant is: $2077.34$

While, for this matrix (rounding of the digits) $$\begin{bmatrix} 123.6 &-1.335 &78.99 &43.26\\ -1.335 &190.6 &88.72 &100.5\\ 78.99 &88.72 &167.5 &0.2444\\ 43.26 &100.5 &0.2444 &143.5 \end{bmatrix}$$

This gives the determinant : $138949$

Why/Where is the instability coming from ?

Bonus : is there a way to have a stability ?

Answer 1

Let us call $M_1$ and $M_2$ the two matrices.

They have the resp. spectra :

$$ S_1=\{0.000294498, 121.945486280, 177.683446086, 325.545773134\}$$

$$S_2= \{0.019697264, 121.950561034, 177.688725537, 325.541016164\}$$

(I am fully conscious that there are too many decimals, considering the roundoff or entries).

Remark : $M_1$ and $M_2$ are symmetric definite positive (all their eigenvalues are $>0$) : a category of matrices that can not be ranked among the worst a priori.

There are different levels of explanations. I will only consider two.

Here is a first one : the "condition numbers" (ratios $\sigma_1/\sigma_4$ of singular values) of these matrices are very high, especially for $M_1$ : $1.1 \times 10^6$ which is huge, indicating that intrinsically $M_1$ is very unstable.

Having a closer look at the spectra, we see that, set apart the first pair of eigenvalues, the rsp. ratios of the other ones are very close to $1$.

Therefore, it will not come as a surprize (the determinant of a matrix being the product of its eigenvalues) that the ratio of the first two eigenvalues, $\approx 66.8840$ is very close to the ratio of determinants which is $\approx 66.8879$.

This shows that this issue of sensibility of determinants can be transfered into an issue of sensibility of the smallest eigenvalue.

I will only provide an insight about the situation by using the pencil of matrices

$$A_x=M_2+x(M_1-M_2)$$

and the attached following functions :

$f(x):=\det(A_x)$ and $g(x) := $ smallest eigenvalue of $A_x$,

for $x=0$ (case of matrix $M_2$) to $x=1$ (case of matrix $M_1$).

What do we see ? That when $x=1$, we are very close to a "catastrophic" situation where, at approximately $x \approx 1.0152$, i.e., at a very small distance of $M_1$, we have a zero determinant... As we have said, this trend can be as well rendered by function $g$ (that has been represented on the same graphic, but with ordinate axis unit multiplied by $10^7$...

Fig. 1 : Graphical representation of functions $f$ (blue) and $g$ ($\times 10^7$) (red) displaying a very similar behavior.

How to go further ? There are many sophisticated ways to analyse sensitivity of eigenvalues. I will not enter into this highly specialized domain, using very specific tools like pseudospectra, etc.

Answer 2

If you take the product along the main diagonal in the second matrix, you get

$$ 123.6 \times 190.6 \times 167.5 \times 143.5 \approx 5.662 \times 10^8. $$

This is roughly $4000$ times as large as the determinant you found.

The product of the main diagonal of the first matrix is also roughly that large.

The only reason the determinants in either case are as small as they are is because the product along the main diagonal, and other positive products that go into the determinant was nearly canceled (to within one part in $4000$ or better) by negative products.

If you do this calculation with numbers that are not precise to much better than $1$ part in $4000$, for example by rounding them as in the second matrix, you get catastrophic cancellation.

Answer 3

Let $A$ be the original matrix and $B$ be the rounded one. Both are symmetric and therefore have very stable eigen-elements. For example, if $Av=\lambda v$, where $v$ is unitary, then $\Delta \lambda=v^T(\Delta A) v$. In particular $||\Delta \lambda||\leq ||\Delta A||_2=||B-A||_2\approx 0.052$.

Yet, $\det(A)$ is stable only when $A$ has no eigenvalues close to $0$. Unfortunately, $A$ admits the couple of eigen-element $\lambda\approx 0.00029,v\approx 1/2[-1,-1,1,1]^T$. Then

$\Delta\lambda\approx 0.0194$ and $\dfrac{\Delta(\det(A))}{\det(A)}\approx\dfrac{\Delta\lambda}{\lambda}\approx 66.9$.

The key of the problem is the bad condition number of $A$. More precisely, the true key is the formula

$\det(B)-\det(A)\approx tr(adj(A)(B-A))\approx 136883$. (cf., above, the @Ian comment), where $adj(A)$ is the adjugate of $A$.

This large value does not depend on the choice of $B$; indeed $B-A$ is quasi random and has (with probability $1$) a non zero component on the eigenvector associated to the maximum eigenvalue of the operator $X\mapsto adj(A)X$; but it depends on the fact that $||adj(A)||_2\approx 7.10^6$ is large (compare with $||A||_2\approx 326$).

What the OP didn't tell us is that

which explains many things.