How can the determinant of a $2\times2$ 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{align} a b - c d &= 10^9 (10^9-2) - (10^9-1)^2 \\ &= (10^{18} - 2\times10^9) - (10^{18} - 2\times10^9 + 1) \\ &= (10^{18} - 2\times10^9) - (10^{18} - 2\times10^9) \quad \text{(the 1 is lost in rounding)} \\ &= 0. \end{align}
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{align} & \begin{vmatrix} 10^9 & 10^9-1 \\ 0 & (10^9-2) - (10^9-1)\times(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} \quad \text{(the $10^{-9}$ is lost in rounding)} \\ &= \begin{vmatrix} 10^9 & 10^9-1 \\ 0 & 0 \end{vmatrix} \\ &= 0. \end{align}
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\times-1) - (10^9-1)\times0 = -1$, which is the correct answer.
Is there a general robust method?
Depending on the problem you’re dealing with, it might be useful to notice that $ad - bc = (a-b)d-b(c-d)$. For the example you provide, it pretty much solves the issue. Just note that it will only make things better if you can expect $|(a-b)d|$ and $|b(c-d)|$ to be much smaller than $|ad|$ and $|bc|$.
Got the idea from here:
https://www.johndcook.com/blog/2018/09/26/polygon-area/