Checking singularity of a matrix

Here's my 2 cents. I would set up the following minimization problem

$$ pi(x) = frac{1}{2} (Ax)^T(Ax) $$

If $A$ has eigenvalues which are zero, there will exist a nonzero $x$ such that $pi(x)=0$. So, I would try computing the gradient of $pi$ wrt $x$ and use a gradient descent algorithm to drive $pi$ towards zero. If you get reasonably close to zero ($piapprox$ 1e-12), then the matrix is singular. The first variation of $pi$ can be computed to be

$$ deltapi = x^TA^TAdelta{x} = (Ax)^TAdelta{x} = g^Tdelta{x}, $$ where $g$ is the gradient. So $g$ is $$ g = A^TAx $$

You'd also need to avoid the $x=0$ case. Starting from a non zero random vector might help.

Here, you need to compute the product of $A^T$ with $Ax$. Given the constraint that you have ( you can only compute $Ax$), I'm not sure if this is possible. Maybe the experts can answer.

If you can compute products with $A$ and $A^T$, as you specify in a comment, you can run the classical sparse SVD algorithms such as scipy.sparse.linalg.svds, Matlab's svds, or Julia's Arpack.svds, which are based on Lanczos bidiagonalization. They are designed to compute singular values, and are likely to be more robust than a minimization routine coded by hand.

Then the distance (in Euclidean or Frobenius norm) from $A$ to the nearest singular matrix is precisely the smallest singular value $sigma_{min}$. You won't be able to return an exact yes/no answer using IEEE arithmetic, so this is the best you can do.

I also suggest looking into the condition number estimators, which will (with some degree of [un]reliability) predict how effectively numerically singular the matrix is.

In particular,

• "Spectral Condition-Number Estimation of Large Sparse Matrices" based on LSQR seems like an interesting choice.

attracted my attention one day. I would also suggest going over the references in this paper to get some perspective on the problem.

@Federico Poloni 's fine answer states the impossibility of getting an exact yes/no answer using IEEE arithmetic.

However, using interval arithmetic with outward rounding, it is possible to get a "not singular/don't know" answer. In particular, it may be possible to definitively conclude that the smallest singular value is strictly greater than zero.

See Verified bounds for singular values, in particular for the spectral norm of a matrix and its inverse, Siegfried M. Rump, BIT, 51(2):367384, 2011.

This interval arithmetic computation using outward rounding of the smallest singular value, can be done, for example, using INTLAB (developed by Siegfried M. Rump) under MATLAB.

The result will be either:

1. "not singular", i.e., the matrix is not singular, because its smallest singular value is definitely > 0

or

2. "don't know", i.e., not definitively determined whether or not the matrix is singular, because not definitively determined whether or not its smallest singular value = 0

"Standard" interval arithmetic with outward rounding (to include, INTLAB) is done using double precision floating point arithmetic. However, it is possible to perform interval arithmetic with outward rounding in higher precision. Doing so might allow a definitive "not singular" result for some matrices which "evaluate" to "don't know" using double precision interval arithmetic with outward rounding. In the (unimplementable) limit of infinite precision interval arithmetic with outward rounding, it can be definitively determined whether or not a given matrix is singular.

I solved the problem by computing the eigenvalue with smallest magnitude with ARPACK. Is this correct?