What would be a reliable way to check if a square matrix is computationally singular? To diagnose computational singularity, one could, for example, check rank deficiency or compute singular-value decomposition. I have the following questions:
(1). Do the two methods listed above generally produce results that coincide with each other?
(2). If the answer to (1) is negative, which method is considered more reliable? 
ADDITIONAL INFO: I don't think I should've used the word "reliable" here. I am writing a C++ function using the Armadillo library. The function will be called in R. Part of the C++ function uses the Armadillo function solve(). Under singularity, the function would throw an exception. But I would like to catch the problem before solve() is executed and output a more meaningful error message to R. So I was thinking perhaps I could check rank deficiency or svd to accomplish that.
 A: Why not use a procedure that will both check the matrix for invertibility and find its inverse?  A QR decomposition will do that.  This is a numerically stable algorithm (employed by many least squares procedures, like lm in R) to express a matrix $A$ as a product
$$A = QR$$
where $Q$ is orthogonal and $R$ is triangular.  Triangular matrices are fast and easy to invert (using back substitution) and orthogonal matrices -- by definition -- are inverted by transposing them, whence
$$A^{-1} = R^{-1}Q^\prime.$$
Thus, even when $A$ is singular the right hand side can be computed.  Consider checking whether it is an inverse by comparing its product with $A$ to the identity matrix.  When agreement is satisfactory (to within a specified tolerance), you're good to go--and you already have a verified version of $A^{-1}.$
I find that typically this test fails in dramatic fashion when the range of singular values of $A$ exceeds $16$ orders of magnitude (roughly the precision of IEEE doubles): one or more off-diagonal elements of $$(R^{-1}Q^\prime)A$$ are $1$ or larger in size.

Here is an R prototype.  It uses the built-in qr function to compute the QR decomposition and applies qr.solve to find $R^{-1}Q^\prime.$  Statistical computing platforms with equivalent capabilities can be similarly programmed.  The principal idiosyncrasy is that the standard test of near-zero numbers, zapsmall, requires you to specify its tolerance as a number of (base 10) digits, while qr allows you to specify an actual tolerance.  I coded this example for approximate consistency of those values.
solve. <- function(A, tol=.Machine$double.eps) {
  B <- qr.solve(qr(A))
  Z <- zapsmall(B %*% A, digits=floor(log10(tol))) - diag(1, nrow(A))
  if(!all(Z==0)) message("System is effectively singular")
  B
}

The following lines generate a largish matrix $A$ that sometimes will be invertible and sometimes not.  Run it a few times (with different dimensions d) to get a feel for the results.
d <- 1000
A <- cbind(matrix(rnorm(d*(d-1)), d), rnorm(d, 0, 3e-16 * d^(2/3)))
system.time(solve.(A))

A: Partially answered in comments:

It depends on how you do it. For instance, if you threshold the
eigenvalues with an SVD to check rank, then obviously the two results
are the same! Also, if you use exact arithmetic the results will be
the same. Perhaps you could edit this question to stipulate the
algorithms you have in mind. When you do so, please also consider telling us a little more about what you mean by "reliable," given that
the singularity of a matrix depends partly on what computations you
intend to do with it.


*

*whuber

Since you will be using R, look at the R function (in base package) kappa which calculates or estimates the condition number, and the following posts:

*

*Kappa condition number in R


*Condition number calculation in R
