How to find out if a matrix is singular? [duplicate]

Question

So how does one identify if a matrix is truly singular? (In MATLAB, without using paper and pencil or symbolic computations, or hand written row operations?) The textbooks sometimes tell students to use the determinant, so I'll start there.

In theory, one can simply test if the determinant of your matrix is zero. Thus

M = magic(4)
M =
    16     2     3    13
     5    11    10     8
     9     7     6    12
     4    14    15     1

det(M)
ans =
  -1.4495e-12

As it turns out, this matrix is indeed singular, so there is a way to write a row of M as a linear combination of the other rows (also true for the columns.) But we got a value that was not exactly zero from det. Is it really zero, and MATLAB just got confused? Why is that? The symbolic determinant is truly zero.

det(sym(M))
ans =
0

As it turns out, computation of the determinant is a terribly inefficient thing for larger arrays. So a nice alternative is to use the product of the diagonal elements of a specific matrix factorization of our square array. In fact, this is what MATLAB does inside det itself for non-symbolic inputs.

[L,U,P] = lu(M)
L =
            1            0            0            0
         0.25            1            0            0
       0.3125      0.76852            1            0
       0.5625      0.43519            1            1
U =
           16            2            3           13
            0         13.5        14.25        -2.25
            0            0      -1.8889       5.6667
            0            0            0   3.5527e-15
P =
     1     0     0     0
     0     0     0     1
     0     1     0     0
     0     0     1     0

See that the diagonal elements of L are ones, but U has non-zero diagonal elements. And there are nice properties about the determinant of a product of matrices, as well as the determinant of upper or lower triangular matrices.

prod(diag(U))
ans =
  -1.4495e-12

See that we got the same answer. Since a LU factorization is reasonably fast even for large arrays, it is a nice way to compute the determinant. The problem is, it uses floating point arithmetic. So those diagonal elements of U are real, floating point numbers. When we take the product, we get a result that is not exactly zero. That last element was just slightly non-zero.

There are other problems with det. For example, we know that the matrix eye(100) is extremely well conditioned. It is an identity matrix after all.

det(eye(100))
ans =
     1

Now, if we multiply a matrix by a constant, this does NOT change the status of the matrix as a singular one. But what happens with the determinant?

det(.1*eye(100))
ans =
   1e-100

So is this matrix singular? After all, if det(magic(4)) gives us 1e-12, then 1e-100 must correspond to a singular matrix! But it does not. Even worse,

det(.0001*eye(100))
ans =
     0

In fact, the determinant is scaled by 0.0001^100, which in matlab will be 1e-400. At least that would be true if matlab could represent a number that small using a double. It cannot do so. That number will underflow. Or as easily, we can make it overflow.

det(10000*eye(100))
ans =
   Inf

Clearly, all of these scaled identity matrices are equally non-singular, but det can be made to give us any answer we want to see! Therefore we must conclude that computing a determinant is a terrible thing to do to a matrix. I don't care what your textbook told you long ago, or what your boss or teacher told you. If somebody told you to compute the determinant for this purpose USING A COMPUTER, that was terrible advice. Period. Determinants simply have too many problems.

We can do other things to test for singularity. The best tool is to use rank. Thus, if the rank of an NxM matrix is less than min(N,M), then the matrix is singular. Here are a couple of tests:

rank(M)
ans =
     3

rank(.0001*eye(100))
ans =
   100

So rank is able to tell us that the 4x4 magic square is singular, but our scaled identity matrix is not singular. (A nice thing is that rank can test for singularity of a non-square matrix.)

We can also use cond to test for numerical singularity. The smallest possible condition number is 1.0, which corresponds to a very well behaved matrix. Large condition numbers are bad. In double precision, this means when the condition number is greater than about 1e15, your matrix is very problematic.

cond(M)
ans =
    8.148e+16

cond(.0001*eye(100))
ans =
     1

In fact, cond perceives that the scaled identity matrix is well conditioned after all. Large condition numbers are bad. For a double precision matrix, a condition number that is anywhere near 1e15 or so indicates a matrix that is probably numerically singular. So we see that M is clearly singular. Again, cond is able to work on non-square matrices.

Or, we could use rcond, a tool which estimates the reciprocal of the condition number. This is a nice tool for really large arrays. When rcond gives a number that is anywhere near eps, WATCH OUT!

rcond(M)
ans =
   1.3061e-17

rcond(.0001*eye(100))
ans =
     1

Finally, for the mathematical gear-heads (like me), we might pull out svd. After all, svd is the tool that cond and rank are based on. When one or more of the singular values of the matrix are tiny compared to the largest singular value, again we have singularity.

svd(M)
ans =
           34
       17.889
       4.4721
   4.1728e-16

Here we look at when a singular value is small compared to the largest singular value of the matrix. A nice thing is svd can tell us how close the matrix is to singularity, and if there are more than one small singular values, if gives us information about the rank of the matrix.

The nice thing is that none of the tools I've shown require the user to do elementary row operations or anything fancy.

BUT PLEASE DON'T USE DET! Yes, it shows up in textbooks. Yes, maybe your instructor or your boss told you to use it. They were just wrong, because tools like this fail when they are applied on a computer, which uses floating point arithmetic. And you simply don't want to compute a symbolic determinant, which will be terribly inefficient.

I'm sorry if this was a long read. I'll get off my soapbox now.

Calculate the rank and compare with the dimension. If the rank is lower than the dimension, then the matrix is singular.

The most reliable approach is to perform a singular value decomposition on the matrix. The ratio of the largest to the smallest singular values should be within some reasonable tolerance. This ratio is the condition number of the matrix. With double precision values, things are getting very dicy with double precision values when the condition number exceeds a million or more, and that's a rather high limit. Note that once you have the SVD, it is useful for a lot of other things than just computing the condition number.

Singular value decomposition is the swiss army chainsaw of numerical analysis; it can be a bit heavy handed of a tool if you know the matrix isn't singular / ill-conditioned. But if you don't know, it's a good tool to know. Particularly so with Matlab since it's a built-in tool.

I would use cond. This gives you a numerical estimate of how close to singular a matrix is (where Inf is a singular matrix).

For example:

m = randn(4);
cond(m)   %Well conditioned, usually in the 10's

m = diag([1e-6 1 2 1e6]);
cond(m)   %Less well conditioned, 1e12

m = diag([0 1 2 3]);
cond(m)   %Singular:  Inf