A few sharp-eyed readers questioned the validity of a technique that I used to demonstrate two ways to solve linear systems of equations. I generated a random n x n matrix and then proceeded to invert it, seemingly without worrying about whether the matrix even has an inverse!

I responded to the comments by saying that "the set of singular (n x n) matrices has dimension n² - 1 within [the set of all n x n matrices]. Therefore if you choose a matrix at random, you are going to choose a singular matrix with probability zero."

This article provides additional details. I assume throughout that the elements of the matrices are real numbers.

Multiplicative inverses

What is a singular matrix? It is a matrix that does not have a multiplicative inverse. That is, a square matrix A is singular if there is no square matrix B for which AB = BA = I.

Scalar matrices (n=1) are just real numbers. All real numbers have a multiplicative inverse except for zero. For n > 1, a matrix is invertible provided that the determinant of the matrix is not zero.

The dimension of the space of singular matrices

How many singular matrices are there, relative to the nonsingular (invertible) matrices? For n=1, the solution is easy to visualize. There is one singular matrix (0) and the rest are nonsingular. The set of singular matrices is a point, which has dimension zero. The set of nonsingular matrices is one-dimensional, because it is the union of two intervals {x | x<0} and {x | x>0}. Therefore the dimension of the set of singular matrices is one less than the dimension of the set of nonsingular matrices.

For n = 2, the set of all 2 x 2 matrices is equivalent to four-dimensional Euclidean space. The determinant vanishes for a singular matrix, therefore the set of singular matrices is equivalent to the set of points {(a,b,c,d} | ac-bd=0}. This is a three-dimensional surface in four-dimensional space. The standard visualization technique for a 3D surface is to slice the surface along the first dimension and looking at the two-dimensional surfaces that result. For each fixed value of a ≠ 0, the slice is a simple quadratic surface known as a hyperbolic paraboloid. For a = 0, the surface is the union of the two planes {c=0} and {d=0}.

In general, the dimension of the set of singular matrices is one less than the dimension of the set of invertible matrices. The reason is that the determinant of a matrix is a multivariable polynomial of its elements. It is shown in courses in multivariable calculus that the zero set of a (nonzero) polynomial is a set in the domain of the function that has measure zero. (For example, see Fleming (1977), Functions of Several Variables (Second Edition), p. 155.) This is obvious in the case of univariate polynomials: the zeros of a polynomial are points, which are zero dimensional subsets of the real line. For a polynomial function of two variables, the zero sets are (generically) the union of curves. For quadratic polynomials, these curves have familiar names: circles, ellipses, hyperbolas, and so forth.

Why is this important? Because if you choose a matrix at random, it will be singular with probability zero.

Moving from real to finite-precision arithmetic

"But, Rick," someone might object, "the computer doesn't compute with real numbers. How do things change when you move to finite precision arithmetic?"

The answer is that the probability of generating a set of singular matrices is no longer zero. It is now a finite—but very small!—probability.

Rather than trudge through finite-precision calculations, let's do an experiment. Let's generate random numbers in SAS and find out how many of them are exactly zero. I could sample from a uniform distribution such as U[-1,1], but instead I'll sample from the standard normal distribution, which generates values near zero more often than values far from zero. The following DATA step generates a billion random numbers:

data a; 
call streaminit(12345);
do i = 1 to 1e9; /* one billion numbers */
   z = rand("normal"); 
   if z=0 then output;
end;
run;

NOTE: The data set WORK.A has 0 observations and 2 variables.

How many times does the program generate an exact zero in a billion attempts? None. Skeptics who think this result depends on the random seed 12345 are welcome to run the experiment with other seeds.

What happens for random 2 x 2 matrices in which each element is drawn from the standard normal distribution? (This is the example from my previous post.) The following DATA step, which might require several minutes to run, also returns no observations:

data c; 
drop j;
call streaminit(54321);
array c[4];
do i = 1 to 1e9; /* one billion matrices */
   do j = 1 to 4;
      c[j] = rand("normal"); 
   end;
   if c[1]*c[4] - c[3]*c[2]=0 then output;
end;
run;

NOTE: The data set WORK.C has 0 observations and 5 variables.

My conclusion? Go ahead and generate a random matrix from the uniform or normal distributions. Chances are very good that the matrix is nonsingular.