Singular matrices & determinants

 High Precision Lens by Kevin Dooley

High Precision Lens by Kevin Dooley

There are a great many problems in wave propagation and scattering that can be reduced to solving equations of the form  \[\det\mathbf{M} = 0,\]where \(\mathbf{M}\) is some, usually, large full matrix. Ocassionally the matrix \(\mathbf{M}\) will be sufficiently small or have a structure that allows the determinant to be evaluated symbolically. In most cases, however, the matrix is simply too large or complicated to handle analytically and numerics become necessary.

computing the determinant of large matrices is almost always a bad idea

The question then becomes, how to determine if a matrix is singular? The obvious answer would be to compute the determinant. However, computing the determinant of large matrices is almost always a bad idea. Using the numerical determinant to, ahem, determine if a matrix is singular is also often inaccurate.

For instance, MATLAB employs LU decomposition to compute the determinant as does NumPy. And whilst LU factorisation is reasonably efficient for large arrays, it does not behave particularly well for ill-conditioned matrices. It is trivially easy to obtain an under or overflow, try typing 

det(1e-3*eye(100))

into MATLAB, for example. Indeed, it is obvious that the determinant of a matrix can be made arbitrarily small or large regardless of its condition number

the determinant of a matrix can be made arbitrarily small or large

These issues can be mitigated, to some extent, by partial pivoting, but given that we are looking for cases where the matrix is singular, we are always going to have problems with LU decomposition.

Numerically, the best option is to perform a Singular Value Decomposition (SVD). The SVD of a matrix \(\mathbf{M}\) is simply a factorisation involving and diagonal matrix \(\mathbf{\Sigma}\) and two unitary matrices \(\mathbf{U}\) and \(\mathbf{V}\) such that \[ \mathbf{M} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^\dagger.\] The determinant of \(\mathbf{M}\) is therefore the product of the diagonal entries of \(\mathbf{\Sigma}\). Hence, the matrix \(\mathbf{M}\) is singular if one of the diagonal entries of \(\mathbf{\Sigma}\) is zero.

A small value (near machine precision) of the ratio of the smallest to largest diagonal element of \(\mathbf{\Sigma}\) indicates that \(\mathbf{M}\) is singular.

Most numerical packages make use of algorithm implemented in the LAPACK library (as well as the GSL), which is was presented by Demmel and Kahan for bidiagonal matrices. Without going into too much detail the algorithm first employs Householder transformations to reduce the matrix to the bidiagonal form. The calculation of the singular values of the bidiagonal matrix is done iteratively via a variant of the QR algorithm.

MATLAB also provides other commands such are cond and rank that can be used to determine if a matrix is singular, but they all make use of SVD. Beware of rcond as this used LU decomposition rather than SVD.

Moving away from ill-condition matrices and singular values, direct computation of the determinant is almost always a bad idea - it is usually extremely inefficient and, as we have seen, can often be misleading.

the determinant, though a convenient notion theoretically, rarely finds a useful role in numerical algorithms
— L. N. Trefethen and D. Bau, Numerical Linear Algebra (SIAM, 1997)