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*^{2} - 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.

## 8 Comments

There is a 2 in 2^64 (±0; 1 in 2^63) chance of drawing a double-precision zero from a random pattern of bits (there are invalid bit patterns, including ±infinity and a variety of forms of Not a Number patterns that were included in the "non-singular" 1x1 matrix, which is arguably a good way of putting it).

So drawing 4 random bit patterns would give you a space of (2^64)^4 = 2^256 possible double-precision numbers....

To complete the thought, the determinant is a single equation, which removes one variable from the set of equations. So (to first order) you have one variable (2^64 possible values) in 2^64^N possible values, or a 1 in 2^(64*(N-1)) chance of drawing a singular matrix (note that N is the number of matrix elements, which scales as V^2 with the number of variables. So the chance of drawing a random matrix scales inversely with the size of the matrix, by my reckoning.

inverse-exponentially (2^(-V^2)), I should say. It's been a long day. :)

Hi Rick,

I have a question but didn't know exactly where I should post it. It concerns finding the root of a square matrix (particularly a correlation matrix). I've been using the ROOT function which returns the Cholesky root. There seem to be several other ways to obtain a root matrix (e.g. the GSORTH call). Is there a method that is more computationally stable than others? I need to do this in some simulations and the ROOT function seems to balk once in a while. I would appreciate your comments.

Post any SAS-related questions to one of the discussion forums at SAS at http://communities.sas.com/index.jspa?view=overview

For SAS/IML questions, there is a SAS/IML Discussion Forum. The direct link is http://communities.sas.com/community/sas_iml_and_sas_iml_studio

To answer your question, if you have a symmetric positive semi-definite matrix, the Cholesky decomposition is your best choice. If ROOT is "balking" (I assume this means "giving an error") then your alleged correlation matrix is not actually positive semi-definite. This means it has a negative eigenvalue, which you can verify with the EIGENVAL function.

Pingback: Readers’ choice 2011: The DO Loop’s 10 most popular posts - The DO Loop

Pingback: Generating a random orthogonal matrix - The DO Loop

Pingback: Popular! Articles that strike a chord with SAS users - The DO Loop