Fast method to check if a Matrix is singular? (non-invertible, det = 0)

Question

What is the fastest algorithm (a link to C or C++ example would be cool) to check if a small square matrix (<16*16 elements) is singular (non-invertible, det = 0) ?

Answer 1

Best way is to compute the condition number via SVD and check if it is greater than 1 / epsilon, where epsilon is the machine precision.

If you allow false negatives (ie. a matrix is defective, but your algorithm may not detect it), you can use the max(a_ii) / min(a_ii) formula from the Wikipedia article as a proxy for the condition number, but you have to compute the QR decomposition first (the formula applies to triangular matrices): A = QR with R orthogonal, then cond(A) = cond(Q). There are also techniques to compute the condition number of Q with O(N) operations, but there are more complex.

Answer 2

I agree with Gaussian elimination. http://math.nist.gov/javanumerics/jama/doc/Jama/LUDecomposition.html documents LU decomposition - after constructing the LU decomposition from a matrix you can call a method on it to get the determinant. My guess is that it is at least worth timing this to compare it with any more specialised scheme.

Answer 3

In R the difference is close, but the qr decomposition wins out. Try the following:

    Sig <- matrix(rnorm(16),4,4)
    
    a <- Sys.time()
    for(j in 1:1000){
      out <- svd(Sig)
      cond <- any(out$d < 1e-10)
    }
    Sys.time()-a
    
    a <- Sys.time()
    for(j in 1:1000){
      out <- qr(Sig)
      cond <- out$rank < NROW(Sig)
    }
    Sys.time()-a

Answer 4

The fastest way is probably to hard code a determinant function for each size matrix you expect to deal with.

Here is some psuedo-code for N=3, but if you check out The Leibniz formula for determinants the pattern should be clear for all N.

function is_singular3(matrix) {
    det = matrix[1][1]*matrix[2][2]*matrix[3][3]
        - matrix[1][1]*matrix[2][3]*matrix[3][2]
        - matrix[1][2]*matrix[2][1]*matrix[3][3]
        + matrix[1][2]*matrix[2][3]*matrix[3][1]
        + matrix[1][3]*matrix[2][1]*matrix[3][2]
        - matrix[1][3]*matrix[2][2]*matrix[3][1];

     if(det==0) return true
     else return false
}