Jtdylktuy

How can the determinant of a 2x2 matrix
$$
begin{vmatrix}
a & b \
c & d
end{vmatrix}
=
a d - b c
$$
be computed in floating point without suffering unnecessary catastrophic cancellation?

For a challenging case, consider the determinant:

$$
begin{vmatrix}
10^9 & 10^9-1 \
10^9-1 & 10^9-2
end{vmatrix}
=
begin{vmatrix}
1 000 000 000 & 999 999 999 \
999 999 999 & 999 999 998
end{vmatrix}
= -1.
$$

All of the numbers above are exactly representable in IEEE double precision
arithmetic, including the determinant which is exactly $-1$.

But naively computing the determinant as $a b - c d$ in double precision produces $0$:
$$
begin{eqnarray}
a b - c d &=& 10^9 (10^9-2) - (10^9-1)^2 \
&=& (10^{18} - 2*10^9) - (10^{18} - 2*10^9 + 1) \
&=& (10^{18} - 2*10^9) - (10^{18} - 2*10^9) textrm{ (the 1 is lost in rounding)} \
&=& 0
end{eqnarray}
$$

A seemingly reasonable alternative approach might be to reduce
the matrix to triangular form using determinant-preserving row operations;
that is, gaussian elimination
(possibly with some form of row and/or column pivoting),
and then take the product of the diagonal entries in the reduced triangular matrix.

Trying that on the above example, we see that
the upper-left entry $a$, being the entry with maximum magnitude,
is already the optimal pivot, so no row or column permutation is done initially.
So gaussian elimination proceeds to zero out the lower-left entry
by subtracting $c/a=(10^9-1)/10^9$ times the first row from the second row, yielding:
$$
begin{eqnarray}
& &
begin{vmatrix}
10^9 & 10^9-1 \
0 & (10^9-2) - (10^9-1)*(10^9-1)/10^9
end{vmatrix} \
&=&
begin{vmatrix}
10^9 & 10^9-1 \
0 & (10^9-2) - (10^9 - 2 + 10^{-9})
end{vmatrix} \
&=&
begin{vmatrix}
10^9 & 10^9-1 \
0 & (10^9-2) - (10^9-2)
end{vmatrix}
textrm{ (the }10^{-9}textrm{ is lost in rounding)} \
&=&
begin{vmatrix}
10^9 & 10^9-1 \
0 & 0
end{vmatrix} \
&=& 0
end{eqnarray}
$$

So gaussian elimination, too, has failed on this example.

This example can be solved in double precision,
via the following ad-hoc sequence of determinant-preserving row and column operations.
Start with the original matrix:
$$
begin{pmatrix}
10^9 & 10^9-1 \
10^9-1 & 10^9-2
end{pmatrix}
$$
Subtract the first row from the second:
$$
begin{pmatrix}
10^9 & 10^9-1 \
-1 & -1
end{pmatrix}
$$
Subtract the second column from the first:
$$
begin{pmatrix}
1 & 10^9-1 \
0 & -1
end{pmatrix}
$$
Then the determinant can be read off as $(1*-1) - (10^9-1)*0 = -1$,
which is the correct answer.

Is there a general robust method?