How to detect singularity of a matrix before inverting it?

Question

I have an algorithm which involves generating a random square matrix and then inverting it. To avoid the program interrupted by numerically singular matrix, I wish to do something like the pseudo-code below;

# generate A
while(A is numerically singular){
# generate A again
}
B = inv(A)

However, I do not know how to write the condition A is singular.

One possible solution is calling the matrixcalc::is.singular.matrix() function, but matlib::inv determines singularity by pivoting and is.singular.matrix() determines it by determinant, so this would not be a valid solution as some matrices can have reasonable determinant but ill-conditioned pivot. For instance, in the following code:

library(matlib)
library(matrixcalc)
A <- matrix(c(10^9, 10^9-1, 1, 1), nrow=2, ncol=2)
check <- is.singular.matrix(A)
B <- matlib::inv(A)

While check will be FALSE, inv will return an error:

Error in Inverse(X, tol = sqrt(.Machine$double.eps), ...) : 
 X is numerically singular

Answer 1

Wrap the calculation in try and check it for class "try-error".

repeat {
 # generate A
 B <- try(solve(A), silent = TRUE)
 if (!inherits(B, "try-error")) break
}

Answer 2

How about Matrix::rcond (estimate the reciprocal condition number of a matrix) with some tolerance?

my_sing <- function(X, tol = 1e-10) {
 Matrix::rcond(X) < tol
}
my_sing(A) ## TRUE

Alternatively, you could pull out the code that matlib::Inverse() is using to do the check ...

my_sing2 <- function(X, tol = sqrt(.Machine$double.eps)) {
 n <- nrow(X)
 X <- matlib::gaussianElimination(X, diag(n), tol = tol)
 any(apply(abs(X[, 1:n]) <= sqrt(.Machine$double.eps), 1, all))
}

• Friede: I am not too sure this indeed needs an apply. Perhaps a rowSums comparison helps here?
• Ben Bolker: @Friede, I agree; I just lazily copied the code from the guts of matlib::Inverse().