How Many Eigenvalues of a Random Matrix are Real?

Alan Edelman^{∗} Eric Kostlan^{†} Michael Shub^{‡}

July 3, 1993

Abstract Let A be an n × n matrix whose elements are independent random variables with standard normal distributions. As n → ∞, the expected number of real eigenvalues is asymptotic to \sqrt{2n/π}. We obtain a closed form expression for the expected number of real eigenvalues for finite n, and a formula for the density of a real eigenvalue for finite n. Asymptotically, a real normalized eigenvalue λ/\sqrt{n} of such a random matrix is uniformly distributed on the interval [−1, 1]. Analogous, but strikingly different, results are presented for the real generalized eigenvalues. We report on numerical experiments confirming these results and suggesting that the assumption of normality is not important for the asymptotic results.

1 Statements of Results

Consider a random matrix whose elements are independent random variables from a standard (mean zero, variance one) normal distribution. Unless otherwise stated, we omit the distribution and simply use the term “random matrix” to denote a matrix with independent standard normally distributed elements. Other distributions are considered in Section 8.

Here is one of our main results:

Asymptotic Number of Real Eigenvalues. If E_{n} denotes the expected number of real eigenvalues of an n by n random matrix, then

lim_{n→∞} \frac{E_{n}}{\sqrt{n}} = \sqrt{\frac{2}{π}} .

Asymptotic Series. As n → ∞,

E_{n} = \sqrt{\frac{2n}{π}} (1 − \frac{3}{8n} − \frac{3}{128n^{2}} + \frac{27}{1024n^{3}} + \frac{499}{32768n^{4}} + O(\frac{1}{n^{5}} )) + \frac{1}{2}.

Corollary 5.2

Let A be a 50 × 50 random matrix. Figure 1 above plots normalized eigenvalues λ/\sqrt{50} in the complex plane for fifty matrices. Thus there are 2500 dots in the figure. There are a number of striking features in the diagram. Firstly, nearly all the normalized eigenvalues appear to fall in the interior of the unit disk. This is Girko’s as yet unverified circular law [12], which states that as n gets large, λ/\sqrt{n} is uniformly distributed Figure 1: 2500 dots representing normalized eigenvalues of fifty random matrices of size n = 50. Clearly visible are the points on the real axis.

∗Department of Mathematics & Lawrence Berkeley Laboratory, University of California, Berkeley, California 94720, edelman@math.berkeley.edu, Supported by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S. Department of Energy under Contract DE-AC03-76SF00098.

†Arts and Sciences, Kapi‘olani Community College, 4303 Diamond Head Road, Honolulu, HI 96816, kostlan@kahuna.math.hawaii.edu

‡TJ Watson Research Center, 32–2, IBM, Yorktown Heights, NY 10598-0218, shub@watson.ibm.com

Figure 2: Histogram of normalized real eigenvalues for 2222 matrices of size 50 in the unit disk. It follows that the proportion of eigenvalues on the real line (also strikingly visible to the eye) must tend to 0 as n → ∞. Our results show how fast this converges.

A simpler version of this circular law occurs when the random matrix has elements from a complex normal distribution, i.e. the real and imaginary parts of each element are independent standard normals. In this case the exact distribution for the eigenvalue distribution and radius can be found in Ginibre [11] and is reported by Mehta [24, p.300] and Hwang [17]. In this case, the squares of the absolute values of the eigenvalues are independent random variables with distributions that are χ^{2} with 2, 4, …, 2n degrees of freedom [20]. The spectral radius is then the maximum of such random variables.

Figure 2 takes a closer look at the real eigenvalues again taking n = 50, but this time we took over 2200 random matrices, and histogrammed the real normalized eigenvalues. Notice the data suggests that the density is nearly uniform on [−1, 1]. The plotted curve is the exact density for n = 50. This suggests the form of the asymptotic density of real eigenvalues that we prove in Corollary 4.5: If λ_{n} denotes a real eigenvalue of an n by n random matrix, then as n → ∞, λ_{n}/\sqrt{n} is uniformly distributed on the interval [−1, 1].

This result is the limit of the probability density for λ_{n} proven in Theorem 4.3: if λ_{n} denotes a real eigenvalue of an n by n random matrix, then its probability density f_{n}(λ) is given by

f_{n}(λ) = \frac{1}{E_{n}} (\frac{1}{\sqrt{2π}} [\frac{Γ(n − 1, λ^{2})}{Γ(n − 1)}] + \frac{|λ^{n−1}|e^{−λ^{2}/2}}{Γ(n/2)2^{n/2}} [\frac{γ((n − 1)/2, λ^{2}/2)}{Γ((n − 1)/2)}]) .

A related function that we study closely in Section 3 is the unnormalized density for λ_{n}. Given a fixed matrix A we can define the empirical cumulative distribution function of its real eigenvalues:

#_{(−∞,x]}(A) ≡ { number of real eigenvalues of A ≤ x}.

Let

ρ_{n}(x) = \frac{d}{dx}E_{A}#_{(−∞,x]}(A),

where E_{A} denotes expectation for random A. Then

f_{n}(λ) = \frac{1}{E_{n}} ρ_{n}(λ).

In fact, if Λ is any Lebesgue measurable set of the real line,

∫_{λ∈Λ} ρ_{n}(λ)dλ = E_{A}#_{Λ}(A) ≡ The expected number of real eigenvalues in Λ.

Most simply put, ρ_{n}(λ) is a true density; it is the “expected number of eigenvalues per unit length” in an infinitesimal interval near λ.

We provide a Mathematica expression for E_{n} below and list enough values of E_{n} to suggest a conjecture which turns out to be true. Table 1 tabulates E_{n} for n from 1 to 10 and suggests a difference in the algebraic form of E_{n} for n even or odd.

We see that a 10 by 10 random matrix can be expected to have fewer than 3 real eigenvalues. More striking is the observation that if n is even, E_{n} is a rational multiple of \sqrt{2}, while if n is odd, E_{n} is one more than a rational multiple of \sqrt{2}. We like to think of this as the “extra” real eigenvalue guaranteed to exist since n is odd. Also notice that the denominators in the ratios are always powers of 2. The observations above and many others may be derived from the exact formulas below.

n  E_{n}
1  1  1.00000
2  \sqrt{2}  1.41421
3 1 + \frac{1}{2} \sqrt{2} 1.70711
4 \sqrt[\frac{11}{8}]{2} 1.94454
5 1 + \frac{13}{16} \sqrt{2} 2.14905
6 \sqrt[\frac{211}{128}]{2} 2.33124
7 1 + \frac{271}{256} \sqrt{2} 2.49708
8 \sqrt[\frac{1919}{1024}]{2} 2.65027
9 1 + \frac{2597}{2048} \sqrt{2} 2.79332
10  \sqrt[\frac{67843}{32768}]{2} 2.92799

Table 1: Expected number of real eigenvalues

Exact Formulas for E_{n}. (Some notation is defined below the box.) If n is even,

E_{n} = \sqrt{2} ∑_{k=0}^{n/2−1} \frac{(4k − 1)!!}{(4k)!!},

while if n is odd,

E_{n} = 1 + \sqrt{2} ∑_{k=1}^{(n−1)/2} \frac{(4k − 3)!!}{(4k − 2)!!}.

Alternatively, for both even and odd n,

E_{n} = \frac{1}{2} + \sqrt{\frac{2}{π}} \frac{Γ(n + 1/2)}{Γ(n)} _{2}F_{1}(1, −1/2; n; 1/2)

= \frac{1}{2} + \sqrt{2} \frac{_{2}F_{1}(1, −1/2; n; 1/2)}{B(n, 1/2)}

= \frac{1 − (−1)^{n}}{2} + \sqrt{2} P_{n−2}^{(1−n, \frac{3}{2} )}(3), if n > 1.

Perhaps nicer yet, we have the generating function

\frac{z(1 − z + z\sqrt{2 − 2z})}{(1 − z)^{2}(1 + z)} = ∑_{n=0}^{∞} E_{n}z^{n}.

Theorem 5.1 and corollaries.

In the formulas above we use the Euler Beta function, a Jacobi polynomial evaluated at three, and also the familiar double factorial (also known as the semifactorial) notation defined by

n!! = \{ \begin{matrix}1 × 3 × 5 × … × n & if n is odd \\ 2 × 4 × 6 × … × n & if n is even. \end{matrix}

By convention, 0!! = (−1)!! = 1.

Mathematica users who may wish to compute E_{n} may do so by typing

                      e[n ]:=(1-(-1)^n)/2 + Sqrt[2] JacobiP[n-2, 1-n, 3/2, 3]

As an example, the above Mathematica expression effortlessly computed the expected number of real eigenvalues of a 100 by 100 random matrix:

E_{100} = 75002314698289190681410505950979137956286758500731773968829\sqrt{2}/2^{193}.

In Section 6, we consider the generalized eigenvalue problem

det(M_{1} − λM_{2}) = 0,

where M_{1} and M_{2} are independent and random. One might guess that questions about generalized eigenvalues would be more difficult than corresponding questions about eigenvalues, but in fact they are simpler.

If E_{n}^{G} denotes the expected number of real generalized eigenvalues of a pair of independent n by n random matrices, then

E_{n}^{G} = \frac{\sqrt{π}Γ((n + 1)/2)}{Γ(n/2)}.

The asymptotic number of real generalized eigenvalues is

lim_{n→∞} \frac{E_{n}^{G}}{\sqrt{n}} = \sqrt{\frac{π}{2}} .

An asymptotic series for this expression as n → ∞, is

E_{n}^{G} = \sqrt{\frac{πn}{2}} (1 − \frac{1}{4n} + \frac{1}{32n^{2}} + \frac{5}{128n^{3}} − \frac{21}{2048n^{4}} + O(\frac{1}{n^{5}} )) .

We also compute the probability density for the real generalized eigenvalues in Theorem 6.2. If λ denotes a real generalized eigenvalue of a pair of independent random matrices, then its probability density f^{G}(λ) is given by

f^{G}(λ) = \frac{1}{π(1 + λ^{2})},

that is, λ obeys the standard Cauchy distribution. Equivalently, atan(λ) is uniformly distributed on [−\frac{π}{2} , \frac{π}{2} ]. Notice that the density function of a real generalized eigenvalue does not depend on n. We could also define ρ_{n}^{G}(λ) in analogy to ρ_{n}(λ), but this will not be of use to us.

2 Motivation, History, Background

Eigenvalues of random matrices arise in many applications areas, perhaps the most well-known areas are nuclear physics, multivariate statistics, and as test matrices for numerical algorithms. See [10] for references to some of these numerous applications. We strongly suspect that random eigenvalue models have been considered in any area where eigenvalues have been considered. The subject is also a favorite for pure mathematicians because it touches on harmonic analysis, combinatorics, and integral geometry.

The first investigation of the eigenvalues of real non-symmetric matrices with normally distributed entries began with Ginibre [11]. He attempted to calculate the probability distribution of the real eigenvalues under the assumption that some fixed number k of them are real, but only succeeded in the case when all of the eigenvalues are real.^{1} In Section 3.5 of [13], Girko derives formulas for the distribution of the eigenvalues under the same assumption that a fixed number are real. Unfortunately, derivations are tedious and the text, at least in translation, contains sufficiently many typographical errors as to make the derivations difficult to check.

1This is an extremely rare event for n not too small. It occurs with probability 2^{−n(n−1)/4}, a fact that will be derived in an upcoming paper.

Research into the analogous question for polynomials has been much more successful, as is well documented in [2]. For example, in the 1940s Kac [18, 19], considered an nth degree polynomial whose coefficients are independent standard normals. He derived an integral formula for the expected number of real roots and was able to show that there are, asymptotically as n → ∞, (2/π) log(n) real roots. Kostlan [21] was able to derive an integral formula for the expected number of real roots of a polynomial with any central normal distribution using the Poincaré formula of integral geometry. Furthermore, Kostlan [21], and Shub and Smale [27] were able to apply geometric methods to show that if the coefficients have independent central normal distributions with variances equal to the binomial coefficients, then the expected number of real roots is exactly the square root of the degree. That these geometric methods, unlike the purely analytic methods of Kac and others, give results for (even underdetermined) systems of equations, demonstrates the power of integral geometry.

Thus from the pure mathematics side, the problem of computing the expected number of real eigenvalues grew out of an attempt to apply integral geometry to linear algebra. The ease with which integral geometry gives the expected number of real generalized eigenvalues (Section 6) gave added hope that the problem of the expected number of real eigenvalues could be solved.

From the applied mathematics side, we wished to respond to a question by Shiu-Hong Lui [23] who was testing homotopy methods to find the eigenvalues and eigenvectors of a general real matrix using random test matrices. Random matrices are often used to test algorithms because of the small effort involved in producing them. As an example, the eigenvalues of random matrices are computed in the LAPACK test suite [4] though LAPACK makes no effort to count the number of eigenvalues that are real.

The physics community has also addressed this problem. Experimental evidence is presented in [22] that the expected number of real eigenvalues is O(\sqrt{n})

3 Eigenvalue Inflation

We begin by defining a process that might be called eigenvalue inflation because it inverts the usual numerical process known as eigenvalue deflation. Let A_{0} be any real n − 1 by n − 1 matrix, v be any unit n-vector such that v_{n} ≥ 0, and w = (w_{1} … w_{n}) be any n − 1 dimensional row vector. We can “inflate” the set of eigenvalues of A_{0} by building the n by n matrix

A ≡ H(v) ( \begin{matrix}& & & 0 \\ & A_{0} & & ⋮ \\ & & & 0 \\ w_{1} & … & w_{n−1} & λ \end{matrix} ) H(v).  (1)

Here H(v) is the linear operator that exchanges v and e_{n} = (0 0 … 0 1)^{T}. For definiteness, let H(v) denote reflection across the angle bisector of v and e_{n}. In numerical linear algebra, reflections of the sort that exchange an arbitrary vector v with e_{n} are usually called Householder reflections, they are orthogonal and symmetric, see [14].

If we make a change of variables from A to v, λ, w and A_{0}, the following lemma tells us how to integrate.

Lemma 3.1 Let Λ be a Lebesgue measurable subset of the real line, let #_{Λ}(A) denote the number of real eigenvalues of A in Λ, and let J(λ, v, A_{0}, w) denote the Jacobian of the transformation defined in (1). Further let τ denote the density function (Radon-Nikodym derivative) of any measure that is absolutely continuous with respect to Lebesgue measure. We then have

∫_{A} #_{Λ}(A)τ(A)dA = ∫_{v,λ∈Λ,w,A_{0}} J(v, λ, w, A_{0}) τ(A(λ, w, A_{0}, v)) dS(v) dλ dw dA_{0},  (2)

where dS(v) is the standard (Haar) measure on the unit sphere, and where dA, dλ, dw and dA_{0} are standard Lebesgue measures. In particular, the expected number of real eigenvalues is

E_{n} ≡ ∫_{A} #_{ℜ}(A)τ(A)dA = ∫_{v,λ,w,A_{0}} J(v, λ, w, A_{0}) τ(A(λ, w, A_{0}, v)) dS(v) dλ dw dA_{0},  (3)

Proof It is easy to see that as A_{0} varies over all n − 1 by n − 1 matrices, w varies over ℜ^{n−1}, and v varies over the unit hemisphere in ℜ^{n}, every matrix A is covered exactly k times, where k is the number of real eigenvalues of A in Λ, unless A falls on the set (of measure zero) of matrices with an eigenvector v where v_{n} = 0 or the set (of measure zero) of matrices with multiple eigenvalues.

Lemma 3.2 The Jacobian of the transformation defined in (1) is

J(v, λ, w, A_{0}) = | det(A_{0} − λI)|.

Proof The proof requires calculation of some differentials near fixed λ, v, A_{0}, and w so that we omit the dependence of H on v etc. In the following, matrices and vectors of differential quantities are in bold face Roman letters so as to distinguish them from the notation for Lebesgue measure.

Notice that v^{T} dv= 0 so that H^{T} dv, which is also the last column of H^{T} dH, has the form (dy_{1} … dy_{n−1} 0)^{T}. The element of surface area in this rotating coordinate system, dS = dy_{1}dy_{2} … dy_{n−1}, is the natural element of surface area on the unit sphere. See Muirhead [25, p.63] for a slightly similar treatment in a more general setting.

Let M denote HAH. Since H^{2} = I, we have HdH= −dHH. Therefore A = HMH and dA=dHMH + HdMH + HMdH or HdAH =dM+(HdH)M − M(HdH). It follows that if we omit the last component

of the last column of HdAH we obtain (A_{0} − λI) ( dy_{1} ; ⋮ ; dy_{n−1} ) . The other elements of HdAH contain

differential forms composed of the corresponding element of dM and the dy_{i}. Taking exterior products of the differential forms of the n^{2} components using standard techniques, we see that

⋀_{ij} dA_{ij} = | det(A_{0} − λI)| dS(v) dλ dw dA_{0},

completing the derivation.

This derivation in terms of differentials almost hides the action on the tangent spaces. Consider the tangent space at e_{n} and ask how does that map to the tangent space at Â in directions orthogonal to Â. A perturbation theory argument would derive a relationship from

(Â + ϵwe_{n}^{T})(e_{n} + ϵy) = (λ + ϵλ_{1})(e_{n} + ϵy),

with the assumption that e_{n}^{T} w = 0 and e_{n}^{T} y = 0. A quick calculation shows that the relationship between the last n − 1 components of w as a function of those of y is given by λI − A_{0}. This is more informative than saying the Jacobian is | det(λI − A_{0})|, because it gives a clear interpretation to the matrix itself.

We now specialize to the case when the matrix A has independent standard normally distributed elements, or, in other words, where τ(A) = (2π)^{−n^{2}/2} exp(−\frac{1}{2} ∑_{i,j=1}^{n} a_{ij}^{2}).

Theorem 3.1

∫_{A} #_{Λ}(A) [(2π)^{−n^{2}/2} exp(−\frac{1}{2} ∑_{i,j=1}^{n} a_{ij}^{2})dA]

= ∫_{λ∈Λ,A_{0}} \frac{π^{1/2}}{2^{(n−1)/2}Γ(n/2)}| det(A_{0} − λI)| [(2π)^{−1/2} exp(−\frac{1}{2}λ^{2})dλ] [(2π)^{−(n−1)^{2}/2} exp(−\frac{1}{2} ∑_{i,j=1}^{n−1}(a_{0})_{ij}^{2})dA_{0}] ,

and where dA, dλ and dA_{0} are standard Lebesgue measures. For clarity, we have placed Gaussian measures in brackets.

Proof By Lemma 3.1 and Lemma 3.2, it is clear that the variables v and w are independent of λ and A_{0} and also they are independent of each other. Thus we can readily integrate out the v and w terms: ∫_{v} dS(v) = \frac{1}{2}Vol(S^{n−1}) = π^{n/2}/Γ(n/2) (where S^{n−1} denotes the unit sphere in ℜ^{n}), and ∫_{w} exp(−\frac{1}{2} ∑ w_{i}^{2}) = (2π)^{(n−1)/2}. From these equations and the previous two lemmas the theorem is immediate.

Taking Λ to be ℜ, we have that

Corollary 3.1

E_{n} = \frac{π^{1/2}}{2^{(n−1)/2}Γ(n/2)}E_{A_{0},λ}| det(A_{0} − λI)|,

where the E denotes expectation over the variables in the subscripts.

Definition 3.1 Let D_{n−1}(λ) = E_{A_{0}}| det(A_{0} − λI)|, where the expected value is taken over all n − 1 by n − 1 matrices A_{0} with independent standard normally distributed elements. Also define

ρ_{n}(λ) = \frac{e^{−λ^{2}/2}}{2^{n/2}Γ(n/2)}D_{n−1}(λ), E_{n} = ∫_{−∞}^{∞} ρ_{n}(λ)dλ, and f_{n}(λ) = \frac{1}{E_{n}} ρ_{n}(λ).

From the discussion above, all of these quantities are related statistically to expectations concerning the real eigenvalues of a random matrix:

ρ_{n}(λ) = lim_{Δ→0} \frac{1}{Δ} (expected number of eigenvalues in [λ − Δ/2, λ + Δ/2]) .

Therefore,

∫_{λ∈Λ} ρ_{n}(λ)dλ  (4)

represents the expected number of eigenvalues in Λ; E_{n} is the expected number of real eigenvalues (i.e. the expected number of eigenvalues in ℜ); f_{n}(λ) is the derivative of the cumulative distribution function of the real eigenvalues. It is sometimes called a condensed density function, in contrast to join densities [2]. Since we consider all the real eigenvalues to be identical, f_{n}(λ) is nothing more than the marginal (probability) density function of a single real eigenvalue. In the next two sections we obtain explicit closed-form expressions for ρ_{n}(λ), E_{n} and f_{n}(λ).

4 Density Computation

The computation of the density of a real eigenvalue of an n by n random matrix proceeds by evaluating D_{n−1} = E_{A_{0}}(| det(A_{0} − λI)|) where A_{0} is an n − 1 by n − 1 random matrix first in terms of objects known as zonal polynomials^{2}, and then in terms of more elementary functions.

2Zonal polynomials arise in group representation theory and the study of related hyergeometric functions [25].

For simplicity we calculate D_{n} instead of D_{n−1}. Let A be an n by n random matrix. From Theorem 10.3.7 of [25, p.447], we have that

D_{n}(λ) = E_{A}(| det(A − λI)|) = \frac{2^{n/2}Γ((n + 1)/2)}{\sqrt{π}}_{1}F_{1}(−\frac{1}{2}; \frac{n}{2}; −\frac{λ^{2}}{2} I_{n}).

The scalar valued hypergeometric function of a matrix argument that appears in the formula above arises in multivariate statistics [25], and should not be confused with the matrix valued function obtained by applying a Taylor series to the matrix. A useful expansion for this hypergeometric function which may be taken as its definition may be found in the proof of Theorem 4.1. An alternative definition for a symmetric n by n matrix X is that

1F_{1}(a; c; X) = constant × ∫_{0<Y <I_{n}} exp(tr(XY))(det Y)^{a−(n+1)/2} det(I_{n} − Y)^{c−a−(n+1)}dY,

where Re(a), Re(c), Re(c−a) > (n−1)/2. The integration is over the symmetric matrices Y for which Y and I_{n} − Y are positive definite matrices. The measure is Lebesgue measure on the upper triangular elements of Y. Lastly, the constant is chosen so that _{1}F_{1}(a; c; 0) = 1.

We introduce the following abbreviation:

Definition 4.1

F_{n}(λ) = _{1}F_{1}(−\frac{1}{2}; \frac{n}{2}; −\frac{λ^{2}}{2} I_{n}).

It is not generally known when hypergeometric functions with a scalar matrix argument can be written as a finite sum of hypergeometric functions with scalar arguments. Gupta and Richards [16] have explored when certain hypergeometric functions of a scalar matrix argument can be written as infinite sums of simpler expressions. In our case, F_{n}(λ) can be written in terms of incomplete gamma functions.

Theorem 4.1

F_{n}(λ) = \frac{e^{λ^{2}/2}Γ(n, λ^{2})}{Γ(n)} + \frac{2^{n−1}}{Γ(n)} (\frac{λ^{2}}{2} )^{n/2} γ(\frac{n}{2}, \frac{λ^{2}}{2} ).

We postpone the proof of this theorem until the end of this section.

Corollary 4.1 The generating function of the F_{n}(λ) is given by

∑_{n=0}^{∞} F_{n}(λ)z^{n} = 1 + \frac{z}{1 − z} e^{(2z−1)λ^{2}/2} + z\sqrt{πλ^{2}/2} e^{z^{2}λ^{2}/2}\{erf[z\sqrt{λ^{2}/2}] + erf[(1 − z)\sqrt{λ^{2}/2}]\}.

Proof Rewrite the formula in Theorem 4.1 using

Γ(n, λ^{2}) = ∫_{λ^{2}}^{∞} t^{n−1}e^{−t}dt and γ(\frac{n}{2}, \frac{λ^{2}}{2} ) = ∫_{0}^{λ^{2}/2} t^{n/2−1}e^{−t}dt.

Switching the order of summation and integration, the generating function can be written as a sum of two integrals. These integrals are easy to evaluate.

In the previous section we established that

ρ_{n}(λ) = \frac{e^{−λ^{2}/2}}{2^{n/2}Γ(n/2)}D_{n−1}(λ).

Thus

ρ_{n}(λ) = \frac{e^{−λ^{2}/2}}{\sqrt{2π}} F_{n−1}(λ).  (5)

Using the duplication formula [1, 6.1.18] to rewrite Γ(n) in the second term of the formula in Theorem 4.1 and then combining with (4) and (5) proves the following two corollaries.

Corollary 4.2 The expected number of eigenvalues on the interval [a, b] is equal to

∫_{a}^{b} (\frac{1}{\sqrt{2π}} [\frac{Γ(n − 1, λ^{2})}{Γ(n − 1)}] + \frac{|λ^{n−1}|e^{−λ^{2}/2}}{Γ(n/2)2^{n/2}} [\frac{γ((n − 1)/2, λ^{2}/2)}{Γ((n − 1)/2)}]) dλ.

Corollary 4.3 If λ_{n} denotes a real eigenvalue of an n by n random matrix, then its marginal probability density f_{n}(λ) is given by

f_{n}(λ) = \frac{1}{E_{n}} (\frac{1}{\sqrt{2π}} [\frac{Γ(n − 1, λ^{2})}{Γ(n − 1)}] + \frac{|λ^{n−1}|e^{−λ^{2}/2}}{Γ(n/2)2^{n/2}} [\frac{γ((n − 1)/2, λ^{2}/2)}{Γ((n − 1)/2)}]) .

The probability density for the normalized eigenvalue x = λ/\sqrt{n} is g_{n}(x) = \sqrt{n}f_{n}(x\sqrt{n}). We wish to understand the limiting behavior of this function as n becomes large.

Corollary 4.4 For all real values of x,

lim_{n→∞} g_{n}(x) = \{ \begin{matrix}1/2 & |x| < 1 \\ (2 + \sqrt{2})/8 & |x| = 1 \\ 0 & |x| > 1 \end{matrix}

Furthermore, the functions g_{n}(x) converge in the L^{p} norms for all 1 ≤ p < ∞.

Proof First we analyze pointwise convergence. We will show in Corollary 5.2 that

lim_{n→∞} \frac{E_{n}}{\sqrt{n}} = \sqrt{\frac{2}{π}} .  (6)

Furthermore [1, 6.5.34],

lim_{m→∞} \frac{Γ(m + α, my)}{Γ(m + α)} = 1 − lim_{m→∞} \frac{γ(m + α, my)}{Γ(m + α)} = \{ \begin{matrix}1 & 0 ≤ y < 1 \\ 1/2 & y = 1 \\ 0 & y > 1 \end{matrix} ,  (7)

and, using Stirling’s (asymptotic to equality) inequality

Γ(m + 1) ≥ m^{m}e^{−m}\sqrt{2πm},  (8)

we can easily establish that

lim_{n→∞} \frac{(\sqrt{n}x)^{n−1}e^{−nx^{2}/2}}{Γ(n/2)2^{n/2}} = \{ \begin{matrix}1/(2\sqrt{π}) & |x| = 1 \\ 0 & |x| ≠ 1 \end{matrix} .  (9)

Combining (6), (7) and (9), we establish the desired pointwise convergence. Using elementary calculus one can show that for all y ≥ 0 and m ≥ 1/2,

y^{2m−1}e^{m−my^{2}} ≤ e^{1−y}.  (10)

Furthermore, the Gaussian continued fraction for the incomplete gamma function [15, 8.358] shows that for y > α − 1,

Γ(α, y) ≤ \frac{e^{−y}y^{α}}{y − α + 1}.  (11)

Using (6), (8), (10) and (11), it is not hard to show that for all sufficiently large n,

g_{n}(x) ≤ e^{1−|x|}.

Thus by the dominated convergence theorem, the sequence {g_{n}} converges in the L^{p} norm for all 1 ≤ p < ∞. Since L^{1} convergence of densities implies convergence in distribution, we have at once another corollary.

Figure 3: Density of normalized eigenvalues for n = 5, 10, 20, 50, 100, 200. The bigger n is the closer it resembles the uniform density on [−1, 1].

Corollary 4.5 If λ_{n} denotes a real eigenvalue of an n by n random matrix, then as n → ∞, the normalized eigenvalue λ_{n}/\sqrt{n} converges in distribution to a random variable uniformly distributed on the interval [−1, 1].

Figure 3 illustrates the convergence to the uniform density. Proof of Theorem 4.1: Following Muirhead, we begin by considering ordered partitions κ of an integer k

κ = (k_{1}, k_{2}, …, k_{n}), where k_{1} + … + k_{n} = k and k_{1} ≥ k_{2} ≥ … ≥ k_{n} ≥ 0.

The (confluent) hypergeometric function of a matrix argument is defined as

1F_{1}(a; b; xI) = ∑_{k=0}^{∞} \frac{x^{k}}{k!} ∑_{κ} \frac{(a)_{κ}}{(b)_{κ}} C_{κ}(I),

where

(a)_{κ} = ∏_{i=1}^{n}(a − (i − 1)/2)_{k_{i}},

(a)_{k} = a(a + 1) … (a + k − 1) = Γ(a + k)/Γ(a)

and the zonal polynomial [25, p.237] is

C_{κ}(I) = 2^{2k}k!(n/2)_{κ} \frac{∏_{i<j}^{n}(2k_{i} − 2k_{j} − i + j)}{∏_{i=1}^{n}(2k_{i} + n − i)!}.

Observe that

(−\frac{1}{2})_{κ} = ∏_{i=1}^{n}(−i/2)_{k_{i}}.  (12)

However, since (−1)_{k} = 0 unless k = 0 or 1, (−\frac{1}{2})_{κ} = 0 unless 1 ≥ k_{2} ≥ k_{3} ≥ … ≥ k_{n} ≥ 0. In other words, we are only interested in partitions where possibly only the first component is not 0 or 1.

We now focus on F_{n}(λ) − F_{n−1}(λ). Since C_{κ}(I_{n})/(n/2)_{κ} is independent of n, the only difference between the expansion for F_{n} and F_{n−1} is the summation over partitions with exactly n non-zero components. To be precise, we may restrict our attention to partitions of the form

κ = (k − n + 1, \underbrace{1, 1, …, 1}_{n−1}), k ≥ n.

We see from (12) that

(−1/2)_{κ} = (−1/2)_{k−n+1} \frac{n!}{(−2)^{n−1}} ,

and that

\frac{C_{κ}(I_{n})}{(n/2)_{κ_{k}}} = 2^{2k}k!\frac{∏_{j=2}^{n}(2(k − n + 1) − 2 − 1 + j) ∏_{j=3}^{n}(j − 2)!}{(2(k − n + 1) + n − 1)! ∏_{i=2}^{n}(2 + n − i)!}

= 2^{2k}k!\frac{\{(2k − n − 1)!/(2k − 2n)!\} ∏_{i=1}^{n−2} i!}{(2k − n + 1)! ∏_{i=1}^{n} i!}

= 2^{2k}k!\frac{(2k − n − 1)!}{(2k − n + 1)!(2k − 2n)!n!(n − 1)!}.

Therefore,

F_{n}(λ) − F_{n−1}(λ) = \frac{1}{(n − 1)!(−2)^{n−1}} ∑_{k=n}^{∞} \frac{(−1/2)_{k−n+1}(2k − n − 1)!2^{2k}}{(2k − n + 1)!(2k − 2n)!} (−λ^{2}/2)^{k}.  (13)

Letting l = k − n and noting

\frac{(−1/2)_{l+1}2^{2l}}{(2l)!} = \frac{−1}{2l!},

we can rewrite (13) as

F_{n}(λ) − F_{n−1}(λ) = \frac{2^{n}}{Γ(n)} (\frac{λ^{2}}{2} )^{n} ∑_{l=0}^{∞} \frac{(−λ^{2}/2)^{l}}{l!(2l + n + 1)(2l + n)}

= \frac{2^{n}}{Γ(n)} (\frac{λ^{2}}{2} )^{n} [∑_{l=0}^{∞} \frac{(−λ^{2}/2)^{l}}{l!(2l + n)} − ∑_{l=0}^{∞} \frac{(−λ^{2}/2)^{l}}{l!(2l + n + 1)}]

= \frac{2^{n−1}}{Γ(n)} (\frac{λ^{2}}{2} )^{(n−1)/2} [(\frac{λ^{2}}{2} )^{1/2} ∑_{l=0}^{∞} \frac{(−1)^{l}(λ^{2}/2)^{l+n/2}}{l!(l + n/2)} − ∑_{l=0}^{∞} \frac{(−1)^{l}(λ^{2}/2)^{l+(n+1)/2}}{l![l + (n + 1)/2]} ]

= \frac{2^{n−1}}{Γ(n)} (\frac{λ^{2}}{2} )^{(n−1)/2} [(\frac{λ^{2}}{2} )^{1/2} γ(\frac{n}{2}, \frac{λ^{2}}{2} ) − γ(\frac{n + 1}{2}, \frac{λ^{2}}{2} )]

= \frac{2^{n−1}}{Γ(n)} (\frac{λ^{2}}{2} )^{n/2} γ(\frac{n}{2}, \frac{λ^{2}}{2} ) − \frac{2^{n−2}}{Γ(n − 1)} (\frac{λ^{2}}{2} )^{(n−1)/2} γ(\frac{n − 1}{2}, \frac{λ^{2}}{2} ) + (\frac{λ^{2}}{2} )^{n−1} \frac{2^{n−1}e^{−λ^{2}/2}}{Γ(n)} .

To calculate F_{n}(λ), we sum the preceding formula over n. The first two terms of the formula telescope and it is only the last term that must be summed. However

∑_{k=1}^{n} (\frac{λ^{2}}{2} )^{k−1} \frac{2^{k−1}e^{−λ^{2}/2}}{Γ(k)} = e^{−λ^{2}/2} ∑_{k=1}^{n} \frac{(λ^{2})^{k−1}}{Γ(k)} = \frac{e^{+λ^{2}/2}Γ(n, λ^{2})}{Γ(n)}.

Thus we see that

F_{n}(λ) = \frac{e^{λ^{2}/2}Γ(n, λ^{2})}{Γ(n)} + \frac{2^{n−1}}{Γ(n)} (\frac{λ^{2}}{2} )^{n/2} γ(\frac{n}{2}, \frac{λ^{2}}{2} ),

as required.

5 Expected Number of Real Eigenvalues

To calculate the expected number of real eigenvalues, we need only perform the integration indicated in Corollary 4.2, taking the interval to be [−∞, ∞]. The integrals involved may be found in classical references (e.g. [15, 6.455]). This produces a closed-form expression for E_{n} in terms of the Gaussian hypergeometric function:

E_{n} = \sqrt{\frac{2}{π}} \frac{Γ(n − 1/2)}{Γ(n)} [n − 1 + \frac{1}{2} _{2}F_{1}(1, n − 1/2; (n + 1)/2; 1/2)] .

We wish to rewrite E_{n} in various forms, each form having its own advantages. The above form was not included in the first section of this paper, because we found it unenlightening. In principle, manipulations of Gaussian hypergeometric functions should be able to prove the equality of any two formulas for E_{n}. However, it is easier to check formulas for E_{n} by computing their generating functions and then comparing them to the result in the following theorem.

Theorem 5.1 The generating function of the E_{n} is given by

∑_{n=0}^{∞} E_{n}z^{n} = \frac{z(1 − z + z\sqrt{2 − 2z})}{(1 − z)^{2}(1 + z)}.

Proof Using the generating function for F_{n} (Corollary 4.1), we can easily produce the generating function for the ρ_{n} and integrate it to produce the generating function appearing in this theorem.

The following corollary will be convenient for computing the asymptotic character of E_{n} for large n.

Corollary 5.1 In terms of Gaussian hypergeometric functions,

E_{n} = \frac{1}{2} + \sqrt{\frac{2}{π}} \frac{Γ(n + 1/2)}{Γ(n)} _{2}F_{1}(1, −1/2; n; 1/2) = \frac{1}{2} + \sqrt{2} \frac{_{2}F_{1}(1, −1/2; n; 1/2)}{B(n, 1/2)}.

Proof Observe that [15, 9.111]

\sqrt{2} _{2}F_{1}(1, −1/2; n; 1/2) = (n − 1) ∫_{0}^{1}(1 − t)^{n−2}\sqrt{2 − t} dt.

Interchanging summation and integration, we can therefore write the generating function for the postulated E_{n} as a single integral. This integral will evaluate to an algebraic function. We then compare this with the generating function in Theorem 5.1.

Corollary 5.2 We have the asymptotic series

E_{n} = \sqrt{\frac{2n}{π}} (1 − \frac{3}{8n} − \frac{3}{128n^{2}} + \frac{27}{1024n^{3}} + \frac{499}{32768n^{4}} + O(\frac{1}{n^{5}} )) + \frac{1}{2}

as n → ∞.

Proof The standard series for the hypergeometric function serves as an asymptotic formula for large n since n appears in the denominator. An asymptotic formula for Γ(n + 1/2)/Γ(n) can be found in [29, 43:6:10].

Corollary 5.3 If n is even,

E_{n} = \sqrt{2} ∑_{k=0}^{n/2−1} \frac{(4k − 1)!!}{(4k)!!},

while if n is odd,

E_{n} = 1 + \sqrt{2} ∑_{k=1}^{(n−1)/2} \frac{(4k − 3)!!}{(4k − 2)!!}.

Proof Using Corollary 5.1, the Gauss recursion formulas for Gaussian hypergeometric functions give

\frac{1}{\sqrt{2}}(E_{n} − E_{n−2}) = \frac{Γ(n − 3/2)}{\sqrt{π}Γ(n − 1)} = \frac{(2n − 5)!!}{(2n − 4)!!}.

Thus it is elementary to establish this corollary by induction.

Corollary 5.4 For n > 1,

E_{n} = \frac{1 − (−1)^{n}}{2} + \sqrt{2} P_{n−2}^{(1−n, \frac{3}{2} )}(3).

Proof The Jacobi polynomials are Gaussian hypergeometric functions. To be precise [15, 8.962.1],

P_{n−2}^{(1−n, \frac{3}{2} )}(z) = (−1)^{n} \frac{4Γ(n + 1/2)}{3\sqrt{π}Γ(n − 1)} _{2}F_{1}(2 − n, 3/2; 5/2; (z + 1)/2).

Rewrite the postulated E_{n} using this formula, and then proceed as in Corollary 5.1, or as in Corollary 5.3.

6 Real Generalized Eigenvalues

A “generalized eigenvalue” of the pair of matrices (M_{1}, M_{2}) (or of the pencil M_{1} − λM_{2}), is defined to be a solution λ to the equation

det(M_{1} − λM_{2}) = 0.

In this section we show how symmetry can be used to obtain the expected number of real generalized eigenvalues and their density.

Theorem 6.1 If E_{n}^{G} denotes the expected number of real generalized eigenvalues of a pair of independent n by n random matrices, then

E_{n}^{G} = \frac{\sqrt{π}Γ((n + 1)/2)}{Γ(n/2)}.

Since the asymptotic series of the Euler Beta function is known [29, 43:6:10] we have an immediate corollary.

Corollary 6.1 We have the asymptotic series

E_{n}^{G} = \sqrt{\frac{πn}{2}} (1 − \frac{1}{4n} + \frac{1}{32n^{2}} + \frac{5}{128n^{3}} − \frac{21}{2048n^{4}} + O(\frac{1}{n^{5}} ))

as n → ∞.

Theorem 6.2 If λ denotes a real generalized eigenvalue of a pair of independent random matrices, then its probability density f^{G}(λ) is given by

f^{G}(λ) = \frac{1}{π(1 + λ^{2})},

that is, λ obeys the standard Cauchy distribution. Equivalently, atan(λ) is uniformly distributed on [−\frac{π}{2} , \frac{π}{2} ].

Since a standard Cauchy random variable can be defined as the ratio of two independent standard normals, it seems appropriate to call the random matrix M = M_{2}^{−1}M_{1} a “(standard) Cauchy matrix.” Clearly the eigenvalues of M are just the generalized eigenvalues of the pair (M_{1}, M_{2}). Thus the expected number of real eigenvalues of an n by n Cauchy matrix is equal to

\frac{\sqrt{π}Γ((n + 1)/2)}{Γ(n/2)},

and a real eigenvalue of a Cauchy matrix is Cauchy.

We now prove these results. A straight forward calculation using Jacobians would be possible here, but we prefer to use the more elegant tools of integral geometry.

Definition 6.1 Let Δ_{n} denote the set of all n by n singular matrices of Frobenius norm one.

Following standard notation, the Frobenius norm of a matrix A is defined as ∥A∥_{F} ≡ \sqrt{∑ a_{ij}^{2}}. In the language of algebraic geometry, Δ_{n} is a real algebraic subvariety of dimension n^{2} − 2 of the unit sphere S^{n^{2}−1} in ℜ^{n^{2}}. Now let (M_{1}, M_{2}) be a pair of matrices. The intersection in ℜ^{n^{2}} of the plane spanned by M_{1} and M_{2} and the sphere S^{n^{2}−1} is a great circle. Real generalized eigenvalues correspond to (pairs of antipodal) intersections of Δ_{n} with this great circle.

Thus when we consider real generalized eigenvalues of the random pair (M_{1}, M_{2}), we are considering intersections of Δ_{n} with random great circles in S^{n^{2}−1}. From the choice of probability measure for the pair (M_{1}, M_{2}), it is not hard to show that the random circles have the standard (Haar) measure. This is a classical set up for integral geometry. We wish to know the expected number of intersections of a fixed variety and a random variety.

Lemma 6.1 (Poincaré) Let V be a variety in S^{m} of dimension m−1. The expected number of intersections of V and a random great circle (with the normalized Haar measure) is equal to twice the volume of V divided by the volume of S^{m−1}.

This formula and its generalizations appear in integral geometry books such as [26]. Poincaré’s formula reduces the problem of calculating the expected number of real generalized eigenvalues to finding the volume of Δ_{n}.

The set Δ_{n} was studied by Demmel [3] and Edelman [8] in the context of studying the probability that a numerical analysis problem is difficult. In particular, they investigated the probability distribution of Demmel’s scaled condition number κ_{D}(M) ≡ ∥M∥_{F} ∥M^{−1}∥_{2}. Computing the volume of Δ_{n} reduces to computing the asymptotics of the probability that κ_{D} > α as α → ∞:

Let V_{ϵ}(Δ_{n}) be the volume of an ϵ neighborhood of Δ_{n} in S^{n^{2}−1}. Clearly,

Vol(Δ_{n}) = lim_{ϵ→0}(2ϵ)^{−1}V_{ϵ}(Δ_{n}).

By the definition of the Demmel condition number κ_{D},

V_{ϵ}(Δ_{n}) = Prob[κ_{D} > 1/ϵ]Vol(S^{n^{2}−1}).

Edelman [8, Corollary 2.2] has shown that

lim_{ϵ→0} ϵ^{−1}Prob[κ_{D} > 1/ϵ] = \frac{2Γ((n + 1)/2)Γ(n^{2}/2)}{Γ(n/2)Γ((n^{2} − 1)/2)}.

We conclude that

Vol(Δ_{n}) = \frac{2π^{n^{2}/2}Γ((n + 1)/2)}{Γ(n/2)Γ((n^{2} − 1)/2)}.

Dividing this by the volume of S^{n^{2}−2} gives the expected number of real generalized eigenvalues.

We now deduce the density function of a real generalized eigenvalue. Consider the pair (M_{1}, M_{2}) to be a collection of n^{2} bivariate normals. The generalized eigenvalue equation may be rewritten

det[cos(θ)M_{1} − sin(θ)M_{2}] = 0.

Since each of the bivariate normals is invariant under rotation, we can readily see that (cos(θ), sin(θ)) is uniformly distributed on the unit circle. Since λ = tan(θ), we have immediately that λ is Cauchy.

7 Numerical Experiments

Fairly early into our investigation, we had some notion that the expected number of real eigenvalues must be roughly 0.8\sqrt{n} from numerical simulations. We were later pleased to learn that this 0.8 was the number \sqrt{2/π}.

With the investigation completed, we can now provide the numerical experiments alongside the exact theoretical results. The numerical experiments were performed using the newly released LAPACK eigenvalue algorithms which we ran on 64 processors of the CM-5 parallel supercomputer. We are pleased to report that the LAPACK algorithm on the CM-5 computed results consistent with our theorems:

Expected number of real eigenvalues: CM-5 Experiments using LAPACK on 64 processors

n trials experimental E_{n} theoretical E_{n} minutes
80 640  7.6  7.603  1
160 640  10.7  10.569  7
320 640  14.9  14.756  51
640 128  20.8  20.673  82
900 64  24.5  24.427  107

We used the CM-5 in what is sometimes called “embarrassingly parallel mode” because each individual matrix never crossed any processor boundaries. Indeed, a 900 by 900 double precision real matrix is about the largest that can fit on any one processor. The results of the computations were sent to the CM-5’s host using the CM-5’s message passing language CMMD.

In order to save some computing time, rather than working with a dense matrix with normally distributed elements, we defined random upper Hessenberg matrices A with exactly the same eigenvalue distribution as a matrix with normally distributed elements. This upper Hessenberg matrix is defined by

a_{ij} is \{ \begin{matrix}normally distributed & i ≤ j \\ distributed like χ_{n−j} & i = j − 1 \\ 0 & otherwise \end{matrix}

To prove that this random matrix does indeed have the same eigenvalue distribution, merely consider the standard reduction to upper Hessenberg form using Householder matrices as described in books such as [14]. The subdiagonal is the length of the column below it which is a χ distribution, the appropriate elements are zeroed out creating Hessenberg form, and the remainder of the matrix remains normally distributed because of the orthogonal invariance.

8 Extensions to Other Distributions

Mehta [24, Conjectures 1.2.1 and 1.2.2] conjectures from extensive numerical experience that the statistical properties of matrices with independent identically distributed entries behave as if they were normally distributed as n → ∞. Mehta focuses on the symmetric or Hermitian cases, but surely the idea is quite general.

Through our own numerical experience, we believe that any eigenvalue property of most any well-behaved distribution can be modeled by the normal distribution. Below are some numerical experiments performed on matrices whose entries came from the uniform distribution on [−1, 1] and also the discrete distribution {−1, 1}. Notice that both of these measures have mean zero and finite variance. Though we have not tested this, we suspect that these are the crucial hypotheses. As indicated in the caption, our CM-5 was upgraded to 128 processors before running these experiments.

Expected number of real eigenvalues: CM-5 Experiments using LAPACK on 128 processors

uniform distribution [-1,1]  discrete distribution {-1,1}
n trials experimental E_{n} minutes trials experimental E_{n} minutes
80 3200  7.6  3.5 3200  7.5  3.3
160 3200  10.6  24.5 3200  10.5  24.1
320 3200  14.9  191 3200  14.8  188
640 896  21.1  412  640  20.8  308
900 384  24.6  499  384  24.7  500

Acknowledgements

We would like to thank Shiu Hong Lui for piquing our interest in this problem as well as Zhimin Yan who encouraged us to carry through a zonal polynomial solution to our problem. We further wish to thank the authors and architects of LAPACK, Matlab, and the CM-5 supercomputer without which such extensive experimentation would not have been possible.

References

[1] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, Dover Publications, New York, 1965.

[2] A.T. Bharucha-Reid and M. Sambandham, Random Polynomials, Academic Press, New York, 1986.

[3] J.W. Demmel, The probability that a numerical analysis problem is difficult, Math. Comp. 50 (1988), 449–480.

[4] J.W. Demmel and A. McKenney, A test matrix generation suite, Argonne National Lab, MCS-P69-0389 and LAPACK working note 9. Available from netlib@na-net.ornl.gov or xnetlib.

[5] A. Edelman, Eigenvalues and condition numbers of random matrices, SIAM J. Matrix Anal. Appl. 9 (1988), 543–560.

[6] A. Edelman, Eigenvalues and Condition Numbers of Random Matrices, PhD thesis, Department of Mathematics, MIT, 1989.

[7] A. Edelman, The distribution and moments of the smallest eigenvalue of a random matrix of Wishart type, Linear Alg. Appl. 159 (1991) 55–80.

[8] A. Edelman, On the distribution of a scaled condition number, Mathematics of Computation 58 (1992), 185–190.

[9] A. Edelman, Random matrix eigenvalues meet numerical linear algebra, SIAM News 24 (November 1991), 11.

[10] A. Edelman, Bibliography of Random Eigenvalue Literature, available electronically by anonymous FTP

from math.berkeley.edu in the directory /pub/edelman.

[11] J. Ginibre, Statistical ensembles of complex, quaternion and real matrices, J. Math. Phys. 6 (1965), 440–449.

[12] V.L. Girko, Circular law, Theory Prob. Appl. 29 (1984), 694–706.

[13] V.L. Girko, Theory of Random Determinants, Kluwer Academic Press, Boston, 1990.

[14] G.H. Golub and C.F. van Loan, Matrix Computations, Second Edition, Johns Hopkins University Press, Baltimore, 1989.

[15] I.S. Gradshteyn and I.M. Ryzhik Table of Integrals, Series and Products, Corrected and Enlarged Edition, Academic Press, New York, 1980.

[16] R.D. Gupta and D.S.P. Richards, Hypergeometric functions of scalar matrix argument are expressible in terms of classical hypergeometric functions, SIAM J. Math. Anal. 16 (1985), 852–858.

[17] C.R. Hwang, A brief survey on the spectral radius and the spectral distribution of large random matrices with i.i.d. entries, Random Matrices and their Applications, Contemporary Mathematics 50, American Mathematical Society, Providence, 1986, pp. 145–152.

[18] M. Kac, On the average number of real roots of a random algebraic equation, Bull. Am. Math. Soc. 49 (1943), 314–320 and 938.

[19] M. Kac, On the average number of real roots of a random algebraic equation (II), Proc. London Math. Soc. 50 (1948), 390–408.

[20] E. Kostlan, Statistical complexity of dominant eigenvector calculation, Linear Alg. Appl. 162–164 (1992), 385–388.

[21] E. Kostlan, On the distribution of the roots of random polynomials, in From Topology to Computation: Proceedings of the Smalefest edited by M.W. Hirsch, J. Marsden and M.Shub, Springer Verlag, New York, 1993, Ch. 38, 419–431.

[22] N. Lehmann and H.-J. Sommers, Eigenvalue statistics of random real matrices, Phys. Rev. Let. 67 (1991), 941–944.

[23] S.H. Lui, private communication, 1992.

[24] M.L. Mehta, Random Matrices, Academic Press, New York, 1991.

[25] R.J. Muirhead, Aspects of Multivariate Statistical Theory, John Wiley & Sons, New York, 1982.

[26] L.A. Santaló, Integral Geometry and Geometric Probability, Vol ume 1 of Encyclopedia of Mathematics and Its Applications, Addison-Wesley, Reading, 1976.

[27] M. Shub and S. Smale, Complexity of Bezout’s Theorem II: Volumes and Probabilities, Proceedings of MEGA 92, to appear.

[28] H.-J. Sommers, A. Crisanti, H. Sompolinsky, and Y. Stein, Spectrum of large random asymmetric matrices, Phys. Rev. Let., 60 (1988), 1895–1898.

[29] J. Spanier and K.B. Oldham, An Atlas of Functions, Hemisphere Publishing, Washington, 1987.