\documentstyle{article}
\setlength{\oddsidemargin}{0in}
\setlength{\topmargin}{-0.5in}
\setlength{\textwidth}{6.5in}
\setlength{\textheight}{9in}
\begin{document}

%For automatic commenting out of the postscript figures interchange
%the two definitions of comment below

\def\comment#1{#1}
%\def\comment#1{}

\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{prop}{Proposition}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{cor}{Corollary}[section]
\newtheorem{defn}{Definition}[section]
\newtheorem{hyp}{Hypothesis}[section]
\newcommand{\proof}{\noindent{\bf Proof }}
\newcommand{\la}{\lambda}
\newcommand{\ha}{\frac{1}{2}}
\def\squarebox#1{\hbox to #1{\hfill\vbox to #1{\vfill}}}
\newcommand{\qed}{\hfill\hfill\vbox{\hrule\hbox{\vrule\squarebox
   {.667em}\vrule}\hrule}\smallskip}
\newlength{\boxwidth}
\setlength{\boxwidth}{\textwidth}
\addtolength{\boxwidth}{-20 pt}
\newcommand{\tr}{\mbox{tr}}
\newcommand{\boxer}[1]{\begin{center} \fbox{\fbox{
      \begin{minipage}{\boxwidth}  #1 \end{minipage}}} \end{center}}



\vspace{0.1in}

\title{
How Many Eigenvalues of a Random Matrix are Real?
}
\author{
Alan Edelman\thanks{
Department of Mathematics
\& Lawrence Berkeley Laboratory,
University of California,
Berkeley, California 94720,
{\protect\verb+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.}\\
Eric Kostlan\thanks{
Arts and Sciences,
Kapi`olani Community College,
4303 Diamond Head Road,
Honolulu, HI 96816,
{\protect\verb+kostlan@kahuna.math.hawaii.edu+}}\\
Michael Shub\thanks{
TJ Watson Research Center, 32--2,
IBM,
Yorktown Heights, NY 10598-0218,
{\protect\verb+shub@watson.ibm.com+}}}
\date{July 3, 1993}
\maketitle
\setlength{\baselineskip}{.5cm}

\begin{abstract}
Let $A$ be an $n \times n$ matrix whose elements
are independent random variables with standard normal distributions.
As $n \rightarrow \infty$, the expected number of real
eigenvalues is asymptotic to $\sqrt{2n/\pi}$.
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
$\lambda/\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.
\end{abstract}

\section{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:
\boxer{
{\sc 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 \rightarrow \infty} \frac{E_n}{\sqrt{n}} = \sqrt{\frac{2}{\pi}}.$$
\\
{\sc Asymptotic Series.} As $n \rightarrow \infty,$
$$E_n=\sqrt{\frac{2n}{\pi}} \left( 1-\frac{3}{8n}-\frac{3}{128n^2}
+\frac{27}{1024n^3}+\frac{499}{32768n^4}+
O(\frac{1}{n^5}) \right) + \frac{1}{2}.$$
{\hfill \small Corollary \ref{series-cor}} }


\comment{
\begin{figure}
\input{psfig}
\begin{center}~
\psfig{figure=dot50.ps,height=3.5in}
\end{center}
\caption{2500 dots representing normalized eigenvalues of
fifty random matrices of size $n=50$. Clearly visible are the points
on the real axis.}
\label{howfdots}
\end{figure}
}

\comment{
\begin{figure}
\input{psfig}
\begin{center}~
\psfig{figure=eig2222.ps,height=3.5in}
\end{center}
{\caption {Histogram of normalized real eigenvalues for 2222 matrices
of size 50}}
\label{howf1}
\end{figure}
}

Let $A$ be a 50 $\times$ 50 random matrix.  Figure \ref{howfdots} above
plots normalized eigenvalues $\lambda/\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 \cite{Gi2}, which
states that as $n$ gets large, $\lambda/\sqrt{n}$ is uniformly
distributed 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 \rightarrow \infty$.  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 \cite{Ginibre} and is reported by Mehta \cite[p.300]{Me}
and Hwang \cite{hwang}.
In this case, the squares of the absolute values of
the eigenvalues are independent random variables with 
distributions that are $\chi^2$ with $2,4,\ldots,2n$
degrees of freedom \cite{ko0}.  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 \ref{density-cor}: If $\lambda_n$ denotes a
real eigenvalue of an $n$ by $n$ random matrix, then as $n \rightarrow
\infty$, $\lambda_n/\sqrt{n}$ is uniformly distributed on the
interval $[-1,1].$ 


This result  is the limit of the
probability density for $\lambda_n$
proven in Theorem \ref{density-theorem}: if $\lambda_n$ denotes a
real eigenvalue of an $n$ by $n$ random matrix, then its probability
density $f_n(\lambda)$ is given by $$ f_n(\lambda)=\frac{1}{E_n}
\left(
\frac{1}{\sqrt{2\pi}}
\left[\frac{\Gamma(n-1,\lambda^2)}{\Gamma(n-1)} \right]+
\frac{|\lambda^{n-1}|e^{-\lambda^2/2}}
{\Gamma(n/2)2^{n/2}}
\left[ \frac{\gamma( (n-1)/2,\lambda^2/2 )}{\Gamma( (n-1)/2)} \right]
\right) .
$$

A related function that we study closely in Section 3 is
the
unnormalized density for $\lambda_n$.
Given a fixed matrix $A$ we can define the empirical cumulative distribution
function of its real eigenvalues:
$$\#_{(-\infty,x]}(A) \equiv \{ \mbox{ number of real eigenvalues of } A
\le x
\}.$$
Let $$\rho_n(x)=\frac{d}{dx} {\rm \bf E}_A\#_{(-\infty,x]}(A),$$
where {\bf E}$_A$ denotes expectation for random $A$.
Then
$$f_n(\lambda)=\frac{1}{E_n}\rho_n(\lambda).$$
In fact, if $\Lambda$ is any Lebesgue measurable set of the real
line,
$$\int_{\lambda \in \Lambda} \rho_n(\la)d\la=
{\rm \bf E}_A\#_\Lambda(A) \equiv
\mbox{ The expected number of real eigenvalues in }\Lambda.$$
Most simply put, $\rho_n(\lambda)$ is a true density; it is
the ``expected number of eigenvalues per unit length'' in an infinitesimal
interval near $\lambda$.



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 \ref{table1} 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.




\begin{table}
\begin{center}
\begin{tabular}{|r|c|r|} \hline
$n$ & $E_n$ & \\ \hline
\rule[-2mm]{0cm}{6mm}
1&$1$ & 1.00000\\ \hline
\rule[-2mm]{0cm}{6mm}
2&$ \sqrt{2}$ & 1.41421\\ \hline
\rule[-2mm]{0cm}{6mm}
3&$1+{1\over 2} \sqrt{2}$ & 1.70711\\ \hline
\rule[-2mm]{0cm}{6mm}
4&${{11}\over 8} \sqrt{2}$ & 1.94454\\ \hline
\rule[-2mm]{0cm}{6mm}
5&$1+{{13}\over {16}} \sqrt{2}$ & 2.14905\\ \hline
\rule[-2mm]{0cm}{6mm}
6&${{211}\over {128}} \sqrt{2}$ & 2.33124\\ \hline
\rule[-2mm]{0cm}{6mm}
7&$1+{{271}\over {256}} \sqrt{2}$ & 2.49708\\ \hline
\rule[-2mm]{0cm}{6mm}
8&${{1919}\over {1024}} \sqrt{2}$ & 2.65027\\ \hline
\rule[-2mm]{0cm}{6mm}
9&$1+{{2597}\over {2048}} \sqrt{2}$ & 2.79332\\ \hline
\rule[-2mm]{0cm}{6mm}
10&${{67843}\over {32768}} \sqrt{2}$ &2.92799\\ \hline
\end{tabular}
\end{center}
\caption{Expected number of real eigenvalues}
\label{table1}
\end{table}
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.




\boxer{
{\sc Exact Formulas for $E_n$}. (Some notation is defined below the
box.)\\ If $n$ is even, $$E_n= \sqrt{2} \sum_{k=0}^{n/2-1}
\frac{(4k-1)!!}{(4k)!!},$$ while if $n$ is odd, $$E_n= 1+\sqrt{2}
\sum_{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}{\pi}} \
\frac{\Gamma(n+1/2)}{ \Gamma(n)} \ {}_2F_1(1,-1/2;n;1/2 )$$
$$=\frac{1}{2}+\sqrt{2}\ \frac{{}_2F_1(1,-1/2;n;1/2)}{B(n,1/2)}
{\hspace{0.55in}}$$ $$=\frac{1-(-1)^n}{2} + \sqrt{2}\
P_{n-2}^{(1-n,\frac{3}{2})}(3), \ \
\mbox{   if  $n>1$.}$$
Perhaps nicer yet, we have the generating function
$$\frac{z(1-z+z\sqrt{2-2z})}{(1-z)^2(1+z)}= \sum_{n=0}^\infty
E_nz^n.$$
\small \hfill Theorem \ref{Engeneratingfunction} 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!!=\left\{
\begin{array}{ll}
1 \times 3 \times 5 \times \ldots \times n & \mbox{if $n$ is odd} \\ 2
\times 4 \times 6 \times \ldots \times n & \mbox{if $n$ is even} .
\end{array}\right. $$
By convention, $0!!=(-1)!!=1$.

Mathematica users who may wish to compute $E_n$
may do so by typing
$$ \mbox{\small
\verb~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-\lambda 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{\pi}\Gamma((n+1)/2)}{\Gamma(n/2)}.$$
The  asymptotic number of real generalized eigenvalues is
$$\lim_{n \rightarrow \infty} \frac{E_n^{G}}{\sqrt{n}} = \sqrt{\frac{\pi}{2}}.$$
An asymptotic series for this expression as  $n \rightarrow \infty,$ is
$$E_n^G=\sqrt{\frac{\pi n}{2}} \left( 1-\frac{1}{4n}+\frac{1}{32n^2}
+\frac{5}{128n^3}-\frac{21}{2048n^4}+
O(\frac{1}{n^5}) \right) .$$


We also compute the
probability density for the real generalized eigenvalues in Theorem 6.2. If
$\lambda$ denotes a real generalized eigenvalue of a pair of
independent random matrices, then its probability density $f^G(\lambda)$
is given by $$ f^G(\lambda)=\frac{1}{\pi ( 1 + \lambda^2 )}, $$ that is,
$\lambda$ obeys the standard Cauchy distribution.  Equivalently,
$\mbox{atan}(\lambda)$ is uniformly distributed on $[-\frac{\pi}{2},
\frac{\pi}{2}]$.
Notice that the density function of a real generalized eigenvalue does
not depend on $n$.  We could also  define $\rho_n^G (\lambda)$ in
analogy to $\rho_n(\lambda)$, but this will not be of use to us.


\section{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
\cite{Edbib} 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 \cite{Ginibre}.  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
{\em all} of the eigenvalues are real.\footnote{This 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.} In
Section 3.5 of \cite{Gi3}, 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.

Research into the analogous question for polynomials has been much
more successful, as is well documented in \cite{BR}.
For example, in
the 1940s Kac \cite{kac1,kac2}, considered
an $n$th 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 \rightarrow \infty$,
$(2/\pi) \log(n)$ real roots.
Kostlan \cite{ko} was able to derive an integral formula for the
expected number of real roots of a polynomial with {\it
any} central normal distribution using the Poincar\'e formula of
integral geometry.  Furthermore, Kostlan \cite{ko}, and Shub and Smale
\cite{SS} 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) {\it 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
\cite{lui} 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 \cite{demck}
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
\cite{sommers1} that
the expected number of real eigenvalues is 
$O(\sqrt{n})$ 



\section{Eigenvalue Inflation}

We begin by defining a process that might be called eigenvalue {\em
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 \ge 0$, and $w=(w_1\ \ldots
\ 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




\be
\label{niceform}
A \equiv H(v)
\left(
\begin{array}{cccc}
&&&0\\ &A_0&&\vdots\\ & & &0\\ w_1 & \ldots & w_{n-1} & \lambda \\
\end{array}
\right)
H(v).
\ee
Here $H(v)$ is the linear operator that exchanges $v$ and $e_n=(0\ 0\
\ldots 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 \cite{GVL}.

If we make a change of variables from $A$ to $v$, $\lambda$,
$w$ and $A_0$, the following lemma tells us how to integrate.

\begin{lemma}
\label{counting-lemma}
Let $\Lambda$ be a Lebesgue measurable subset of the real line,
let
$\#_\Lambda(A)$ denote the number of real eigenvalues of $A$
in $\Lambda$,  and let
$J(\lambda,v,A_0,w)$ denote the Jacobian of the transformation
defined in (\ref{niceform}). Further let $\tau$ denote the density
function (Radon-Nikodym derivative) of any measure that is absolutely
continuous with respect to Lebesgue measure.  We then have
\be
\label{integral-d}
\int_A {\#}_\Lambda(A)\tau(A) dA  =
\int_{v,\lambda \in \Lambda,w,A_0} J(v,\la,w,A_0) \
\tau \! \left(A(\lambda,w,A_0,v)\right)\ dS(v)\ d\lambda\ dw\ dA_0 ,
\ee
where $dS(v)$ is the standard (Haar) measure on the unit sphere, and
where $dA$, $d\lambda$, $dw$ and $dA_0$ are standard Lebesgue
measures.  In particular, the expected number of real eigenvalues
is
\be
\label{integral-e}E_n \equiv
\int_A {\#}_\Re (A)\tau(A) dA  =
\int_{v,\lambda,w,A_0} J(v,\la,w,A_0) \
\tau \! \left(A(\lambda,w,A_0,v)\right)\ dS(v)\ d\lambda\ dw\ dA_0 ,
\ee
\end{lemma}

\proof
It is easy to see that as $A_0$ varies over all $n-1$ by $n-1$
matrices, $w$ varies over $\Re^{n-1}$, and $v$ varies over the unit
hemisphere in $\Re^n$, every matrix $A$  is
covered exactly $k$ times, where $k$ is the number of real
eigenvalues of $A$ in $\Lambda$,
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.


\begin{lemma}
\label{jacobian-lemma}
The Jacobian of the transformation defined in (\ref{niceform}) is
$$J(v,\la,w,A_0)=|\det(A_0-\lambda I)|.$$
\end{lemma}
\proof
The proof requires calculation of some differentials near fixed
$\lambda,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${\bf 	dv}$=0$ so that $H^T${\bf dv}, which is also the last
column of $H^T${\bf dH}, has the form $(dy_1\ \ldots dy_{n-1}\ 0)^T$.  The
element of surface area in this rotating coordinate system,
$dS=dy_1dy_2 \ldots dy_{n-1}$, is the natural element of surface area
on the unit sphere.  See Muirhead \cite[p.63]{Mu} for a slightly
similar treatment in a more general setting.

Let $M$ denote $HAH$.  Since $H^2=I$, we have $H${\bf dH}$=-${\bf dH}$H$.
Therefore $A=HMH$ and
{\bf dA}$=${\bf dH}$MH+H${\bf dM}$H+HM${\bf dH} or
$H${\bf dA}$H=${\bf dM}$+(H${\bf dH}$)M-M(H${\bf dH}$)$.
It follows that if we omit
the last component of the last column of $H${\bf dA}$H$ we obtain
$(A_0-\lambda I)\left(\begin{array}{c}dy_1\\ \vdots\\ dy_{n-1}
\end{array}\right).$
The other elements of $H${\bf dA}$H$ contain differential forms composed of
the corresponding element of {\bf dM} and the $dy_i$.  Taking
exterior products of the differential forms of the $n^2$ components
using standard techniques, we see that $$\bigwedge_{ij} dA_{ij}=
|\det(A_0-\la I)| \ dS(v)\ d\lambda\ 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 $\hat{A}$ in
directions orthogonal to $\hat{A}$.  A perturbation theory argument
would derive a relationship from
$$(\hat{A}+\epsilon
we_n^T)(e_n+\epsilon y)=(\lambda+\epsilon \lambda_1) (e_n+\epsilon
y),$$
with the assumption that $e_n^Tw=0$ and $e_n^Ty=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 $\lambda
I-A_0.$ This is more informative than saying the Jacobian is
$|\det(\lambda 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
$\tau(A)=(2\pi)^{-n^2/2}
\exp(-\ha\sum_{i,j=1}^n a_{ij}^2).$


\begin{theorem}
\label{integral-geometry-theorem}
$$
\int_A \#_\Lambda(A) \left[ (2\pi)^{-n^2/2}\exp(-\ha\sum_{i,j=1}^n a_{ij}^2) dA
\right]$$ $$ = \int_{\lambda \in \Lambda,A_0}
\frac{\pi^{1/2}}{2^{(n-1)/2}\Gamma(n/2)}
|\det(A_0-\lambda I)|
\left[
(2\pi)^{-1/2}\exp(-\ha\lambda^2) d\lambda \right]
\left[ (2\pi)^{-(n-1)^2/2}\exp(-\ha\sum_{i,j=1}^{n-1} (a_0)_{ij}^2)
dA_0 \right], $$ and where $dA$, $d\lambda$ and $dA_0$ are standard
Lebesgue measures.  For clarity, we have placed Gaussian measures in
brackets.
\end{theorem}

\proof
By Lemma
\ref{counting-lemma}
and Lemma \ref{jacobian-lemma}, it is clear that the variables $v$
and $w$ are independent of $\lambda$ and $A_0$ and also they are independent
of each other.
Thus we can readily integrate
out the $v$ and $w$ terms: $\int_v dS(v)=\ha \mbox{Vol($S^{n-1}$)}=\pi^{n/2}/\Gamma(n/2)$ (where $S^{n-1}$ denotes the unit
sphere in $\Re^n$), and $\int_w \exp(-\ha\sum
w_i^2)=(2\pi)^{(n-1)/2}.$ {From} these equations and the previous two
lemmas the theorem is immediate.

Taking $\Lambda$ to be $\Re$, we have  that
\begin{cor}
$$E_n=\frac{\pi^{1/2}}{2^{(n-1)/2}\Gamma(n/2)}
\mbox{\rm \bf E}_{A_0,\lambda}|\det(A_0-\lambda I)|,$$
where the {\rm \bf  E} denotes expectation over the variables
in the subscripts.
\end{cor}


\begin{defn}
Let  $D_{n-1}(\lambda) = {\rm \bf E}_{A_0}|\det(A_0-\lambda 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
$$\rho_n(\lambda)=\frac{e^{-\la^2/2}}{2^{n/2}\Gamma(n/2)}
D_{n-1}(\lambda),\ \ \ E_n=\int_{-\infty}^{\infty} \rho_n(\lambda)
d\lambda,\ \ \
\mbox{ and} \ \ f_n(\lambda) = \frac{1}{E_n} \rho_n(\lambda).$$
\end{defn}

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

$$\rho_n(\lambda)=\lim_{\Delta \rightarrow 0}
\frac{1}{\Delta}\left( \mbox{\rm expected
number of eigenvalues in }\left[\lambda-{\Delta}/{2},
\lambda+{\Delta}/{2} \right] \right).$$
Therefore,
\be
\label{rhoj}
\int_{\lambda \in \Lambda} \rho_n(\lambda) d\lambda
\ee
represents the expected number of eigenvalues in $\Lambda$;
$E_n$  is the expected number of
real eigenvalues (i.e. the expected number of eigenvalues in $\Re$);
$f_n(\la)$  is the derivative of the cumulative distribution function
of the real eigenvalues.
It  is sometimes called a {\it condensed
density function}, in contrast to join densities
\cite{BR}.  Since we consider all the real eigenvalues to be identical,
$f_n(\lambda)$ is nothing more than the {\it marginal (probability) density
function} of a single real eigenvalue.  In the next two sections we
obtain explicit closed-form expressions for $\rho_n(\lambda)$, $E_n$ and
$f_n(\lambda)$.




\section{Density Computation}

The computation of the density of a real eigenvalue of an $n$ by $n$
random matrix proceeds by evaluating $D_{n-1} =
{\rm \bf E}_{A_0}(|\det(A_0-\lambda I)|)$
where $A_0$ is an $n-1$ by $n-1$ random
matrix first in terms of objects known as
zonal polynomials\footnote{Zonal polynomials arise in
group representation theory and
the study of  related hyergeometric functions \cite{Mu}.},
and then in terms of more elementary functions.

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
\cite[p.447]{Mu}, we have that $$D_n(\la)={\rm \bf E}_A(|\det(A-\la I)|)
=\frac{2^{n/2} \Gamma((n+1)/2)}{\sqrt{\pi}}
{}_1F_1(-\ha;\frac{n}{2};-\frac{\la^2}{2}I_n).$$
The scalar valued
hypergeometric function of a  matrix argument
that appears in the formula above
arises in multivariate statistics \cite{Mu}, 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)=\mbox{constant} \times
 \int_{0<Y<I_n} \exp(\tr(XY))(\det Y)^{a-(n+1)/2}
\det(I_n-Y)^{c-a-(n+1)} dY,$$
where $\mbox{Re}(a),\mbox{Re}(c),\mbox{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 ${}_1F_1(a;c;0)=1$.


We introduce the
following abbreviation:
\begin{defn}
$$ F_n (\lambda) = {}_1F_1(-\ha;\frac{n}{2};-\frac{\la^2}{2}I_n).$$
\end{defn}

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  \cite{GuR} 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(\lambda)$ can be written in terms of incomplete gamma
functions.
\begin{theorem}
\label{FnFormula}
$$F_n(\lambda) = \frac{e^{\lambda^2/2} \Gamma(n,\lambda^2)
}{\Gamma(n)} + \frac{2^{n-1}}{\Gamma(n)} \left( \frac{\lambda^2}{2}
\right)^{n/2}
\gamma(\frac{n}{2},\frac{\lambda^2}{2}).$$
\end{theorem}
We postpone the proof of this theorem until the end of this section.

\begin{cor}
\label{Fngeneratingfunction}
The generating function of the $F_n(\lambda)$ is given by $$
\sum_{n=0}^\infty F_n(\lambda) z^n = 1 + \frac{z}{1-z} e^{(2z-1)\lambda^2/2}
+ z\sqrt{\pi \lambda^2/2}\ e^{z^2\lambda^2/2}
\{ \mbox{\rm erf} [z\sqrt{\lambda^2/2}]
+ \mbox{\rm erf} [(1-z)\sqrt{\lambda^2/2}] \} .  $$
\end{cor}
\proof
Rewrite the formula in Theorem
\ref{FnFormula} using
$$ \Gamma(n,\lambda^2)=\int_{\lambda^2}^\infty t^{n-1} e^{-t} dt
\ \ \mbox{\rm and}\ \ \gamma(\frac{n}{2},\frac{\lambda^2}{2})=
\int_0^{\lambda^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 $$\rho_n(\lambda)
=\frac{e^{-\lambda^2/2}}{2^{n/2}\Gamma(n/2)} D_{n-1}(\lambda).$$ Thus
\be
\label{rhon}
\rho_n(\lambda)=\frac{e^{-\lambda^2/2}}{\sqrt{2\pi}} F_{n-1}(\lambda).
\ee

\noindent
Using the duplication  formula \cite[6.1.18]{Abr}
to rewrite $\Gamma(n)$ in the second term
of the formula in Theorem \ref{FnFormula}  and then combining with
(\ref{rhoj}) and (\ref{rhon}) proves the following two corollaries.
\begin{cor}
\label{rhocorollary}
The expected number of eigenvalues on the interval $[a,b]$ is equal to
$$
\int_a^b \left(
\frac{1}{\sqrt{2\pi}}
\left[\frac{\Gamma(n-1,\lambda^2)}{\Gamma(n-1)} \right]+
\frac{|\lambda^{n-1}|e^{-\lambda^2/2}}
{\Gamma(n/2)2^{n/2}}
\left[ \frac{\gamma( (n-1)/2,\lambda^2/2 )}{\Gamma( (n-1)/2)} \right]
\right) d \lambda.
$$
\end{cor}

\begin{cor}
\label{density-theorem}
If $\lambda_n$ denotes a real eigenvalue of an $n$ by $n$ random
matrix, then its marginal probability density $f_n(\lambda)$ is given by $$
f_n(\lambda)=\frac{1}{E_n} \left(
\frac{1}{\sqrt{2\pi}}
\left[\frac{\Gamma(n-1,\lambda^2)}{\Gamma(n-1)} \right]+
\frac{|\lambda^{n-1}|e^{-\lambda^2/2}}
{\Gamma(n/2)2^{n/2}}
\left[ \frac{\gamma( (n-1)/2,\lambda^2/2 )}{\Gamma( (n-1)/2)} \right]
\right) .
$$
\end{cor}




The probability density for the normalized eigenvalue
$x=\lambda/\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.
\begin{cor} For all real values of $x$,
$$ \lim_{n \rightarrow \infty} g_n(x) = \left\{ \begin{array}{cc} 1/2 &
|x| < 1 \\ (2+\sqrt{2})/8 & |x| = 1 \\ 0 & |x| > 1
\end{array} \right.  $$
Furthermore, the functions $g_n(x)$ converge in the $L^p$ norms for all
$1 \le p < \infty$.
\end{cor}
\proof
First we analyze pointwise convergence.  We will show in Corollary
\ref{series-cor} that
\be
\label{asymptforEn}
\lim_{n \rightarrow \infty} \frac{E_n}{\sqrt{n}} = \sqrt{\frac{2}{\pi}}.
\ee
Furthermore \cite[6.5.34]{Abr},
\be
\label{gammalimit}
\lim_{m \rightarrow \infty} \frac{\Gamma(m+\alpha,m y)}{\Gamma(m+\alpha)}
= 1 - \lim_{m \rightarrow \infty} \frac{\gamma(m+\alpha,m
y)}{\Gamma(m+\alpha)} = \left\{ \begin{array}{cc} 1 & 0 \le y < 1 \\
1/2 & y = 1 \\ 0 & y > 1
\end{array} \right. ,
\ee
and, using Stirling's (asymptotic to equality) inequality
\be
\label{Stirling}
\Gamma(m+1) \ge m^m e^{-m} \sqrt{2\pi m},
\ee
we can easily establish that
\be
\label{chilimit}
\lim_{n \rightarrow \infty}
\frac{(\sqrt{n}x)^{n-1} e^{-nx^2/2}}{\Gamma(n/2)2^{n/2}}
= \left\{ \begin{array}{cc} 1/(2\sqrt{\pi}) & |x|=1 \\ 0 & |x| \ne 1
\end{array} \right. .
\ee
Combining (\ref{asymptforEn}), (\ref{gammalimit}) and
(\ref{chilimit}), we establish the desired pointwise convergence.
Using elementary calculus one can show that for all $y \ge 0 $ and $m
\ge 1/2$,
\be
\label{mylemma}
y^{2m-1} e^{m-my^2} \le e^{1-y}.
\ee
Furthermore, the Gaussian continued fraction for the incomplete gamma
function \cite[8.358]{GR} shows that for $y>\alpha-1$,
\be
\label{gammainequality}
\Gamma(\alpha,y) \le \frac{e^{-y} y^\alpha}{y-\alpha+1}.
\ee
Using (\ref{asymptforEn}), (\ref{Stirling}), (\ref{mylemma}) and
(\ref{gammainequality}), it is not hard to show that for all
sufficiently large $n$,
$$g_n(x) \le e^{1-|x|}.$$ Thus by the dominated convergence theorem,
the sequence $\{g_n\}$ converges in the $L^p$ norm for all $1 \le p <
\infty$.

Since $L^1$ convergence of densities implies convergence in
distribution, we have at once another corollary.

\begin{cor}
\label{density-cor}
If $\lambda_n$ denotes a real eigenvalue of an $n$ by $n$ random
matrix, then as $n \rightarrow \infty$, the
normalized eigenvalue  $\lambda_n/\sqrt{n}$
converges in distribution to a random variable uniformly distributed
on the interval $[-1,1].$ \\
\end{cor}


\comment{
\begin{figure}
\input{psfig}
\begin{center}~
\psfig{figure=manydense.ps,height=3in}
\end{center}
\caption{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]$.}
\label{howf3}
\end{figure}
}

Figure \ref{howf3} illustrates the convergence to the uniform density.




\noindent{\bf Proof of Theorem 4.1:}
Following Muirhead, we begin by considering ordered partitions
$\kappa$ of an integer $k$ $$
\kappa=(k_1,k_2,\ldots,k_n),\mbox{  where  }  k_1+\ldots+k_n=k
\mbox{ and } k_1 \ge k_2 \ge \ldots \ge k_n \ge 0.$$
The (confluent) hypergeometric function of a matrix argument is defined as
$$ {}_1F_1(a;b;xI)=\sum_{k=0}^\infty \frac{x^k}{k!}
\sum_\kappa \frac{(a)_\kappa}{(b)_\kappa} C_\kappa(I),$$
where $$(a)_\kappa=\prod_{i=1}^{n} (a-(i-1)/2)_{k_i},$$
$$(a)_k=a(a+1)\ldots (a+k-1)=\Gamma(a+k)/\Gamma(a)$$ and the zonal
polynomial \cite[p.237]{Mu} is $$
C_\kappa(I)=2^{2k}k!(n/2)_\kappa
\frac
{ \prod_{i<j}^n (2k_i-2k_j-i+j)} {\prod_{i=1}^n (2k_i+n-i)!}.  $$
Observe that
\be
\label{haeq}
(-\ha)_\kappa=\prod_{i=1}^n (-i/2)_{k_i} .
\ee
However, since $(-1)_k=0$ unless $k=0$ or $1$,
$(-\ha)_\kappa=0$ unless $1 \ge k_2 \ge k_3 \ge \ldots \ge k_n \ge 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(\lambda)-F_{n-1}(\lambda)$.
Since $C_\kappa(I_n)/(n/2)_\kappa$ 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
$$\kappa =(k-n+1,\underbrace{1,1,\ldots,1}_{n-1}), \ \
k \ge n .$$
We see from (\ref{haeq}) that
$$(-1/2)_{\kappa}=(-1/2)_{k-n+1}\frac{n!}{(-2)^{n-1}},$$ and that
\begin{eqnarray*}
\lefteqn {\frac{C_{\kappa}(I_n)}{(n/2)_{\kappa_k}}=
2^{2k}k! \frac{ \prod_{j=2}^n (2(k-n+1)-2-1+j) \prod_{j=3}^n (j-2)!}
{(2(k-n+1)+n-1)! \prod_{i=2}^n (2+n-i)!}}&& \\ &&= 2^{2k}k! \frac{
\{ (2k-n-1)!/(2k-2n)! \} \prod_{i=1}^{n-2} i!}{
(2k-n+1)! \prod_{i=1}^n i!} \\ &&= 2^{2k}k! \frac{
(2k-n-1)!}{(2k-n+1)!(2k-2n)!n!(n-1)!} .
\end{eqnarray*}
Therefore,
\be
\label{Fnfn-1}
F_n(\la)-F_{n-1}(\la)=
\frac{1}{(n-1)!(-2)^{n-1}}\sum_{k=n}^\infty
\frac{(-1/2)_{k-n+1}(2k-n-1)!2^{2k}}{(2k-n+1)!(2k-2n)!}(-\la^2/2)^k.
\ee
Letting $l=k-n$ and noting
$$\frac{(-1/2)_{l+1} 2^{2l}}{(2l)!} =\frac{-1}{2l!},$$
we can rewrite (\ref{Fnfn-1}) as
\begin{eqnarray*}
\lefteqn{
F_n(\la)-F_{n-1}(\la)=
\frac{2^n}{\Gamma(n)}
\left( \frac{\lambda^2}{2} \right)^n
\sum_{l=0}^\infty
\frac{(-\la^2/2)^l }{l!(2l+n+1)(2l+n)}} \\
&&=
\frac{2^n}{\Gamma(n)}
\left( \frac{\lambda^2}{2} \right)^n
\left[ \sum_{l=0}^\infty \frac{(-\la^2/2)^l }{l!(2l+n)}
- \sum_{l=0}^\infty \frac{(-\la^2/2)^l }{l!(2l+n+1)} \right] \\ &&=
\frac{2^{n-1}}{\Gamma(n)}
\left( \frac{\lambda^2}{2} \right)^{(n-1)/2}
\left[ \left( \frac{\lambda^2}{2} \right)^{1/2}
\sum_{l=0}^\infty \frac{(-1)^l (\lambda^2/2)^{l+n/2} }{l!(l+n/2)}
- \sum_{l=0}^\infty \frac{(-1)^l (\lambda^2/2)^{l+(n+1)/2}
}{l![l+(n+1)/2]} \right] \\ &&=
\frac{2^{n-1}}{\Gamma(n)}
\left( \frac{\lambda^2}{2} \right)^{(n-1)/2}
\left[ \left( \frac{\lambda^2}{2} \right)^{1/2}
\gamma(\frac{n}{2},\frac{\lambda^2}{2})
- \gamma(\frac{n+1}{2},\frac{\lambda^2}{2}) \right] \\ &&=
\frac{2^{n-1}}{\Gamma(n)} \left( \frac{\lambda^2}{2} \right)^{n/2}
\gamma(\frac{n}{2},\frac{\lambda^2}{2})
-
\frac{2^{n-2}}{\Gamma(n-1)} \left( \frac{\lambda^2}{2} \right)^{(n-1)/2}
\gamma(\frac{n-1}{2},\frac{\lambda^2}{2})
+ \left( \frac{\lambda^2}{2} \right)^{n-1}
\frac{2^{n-1} e^{-\lambda^2/2} }{\Gamma(n)} .
\end{eqnarray*}

To calculate $F_n(\lambda)$, 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 $$ \sum_{k=1}^n
\left( \frac{\lambda^2}{2} \right)^{k-1}
\frac{2^{k-1} e^{-\lambda^2/2} }{\Gamma(k)}
= e^{-\lambda^2/2} \sum_{k=1}^n
\frac{(\lambda^2)^{k-1} }{\Gamma(k)}
= \frac{ e^{+\lambda^2/2} \Gamma(n,\lambda^2)}{\Gamma(n)}.  $$ Thus we
see that $$F_n(\lambda) = \frac{e^{\lambda^2/2} \Gamma(n,\lambda^2)
}{\Gamma(n)} + \frac{2^{n-1}}{\Gamma(n)} \left( \frac{\lambda^2}{2}
\right)^{n/2}
\gamma(\frac{n}{2},\frac{\lambda^2}{2}),$$
as required.


\section{Expected Number of Real Eigenvalues}

To calculate the expected number of real eigenvalues, we need only
perform the integration indicated in Corollary \ref{rhocorollary},
taking the interval to be $[-\infty,\infty]$.  The integrals involved
may be found in classical references (e.g. \cite[6.455]{GR}).  This
produces a closed-form expression for $E_n$ in terms of the Gaussian
hypergeometric function: $$ E_n = \sqrt{\frac{2}{\pi}}
\frac{\Gamma(n-1/2)}{\Gamma(n)}
\left[ n-1  + \frac{1}{2}\  {}_2F_1 (1,n-1/2;(n+1)/2;1/2) \right].
$$ 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.

\begin{theorem}
\label{Engeneratingfunction}
The generating function of the $E_n$ is given by $$
\sum_{n=0}^\infty E_n z^n =
\frac{z(1-z+z\sqrt{2-2z})}{(1-z)^2(1+z)}.
$$
\end{theorem}
\proof
Using the generating function for $F_n$ (Corollary
\ref{Fngeneratingfunction}), we can easily produce the generating
function for the $\rho_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$.
\begin{cor}
\label{hyperEn}
In terms of Gaussian hypergeometric functions,
$$E_n=\frac{1}{2}+\sqrt{\frac{2}{\pi}} \ \frac{\Gamma(n+1/2)}{
\Gamma(n)} \ {}_2F_1(1,-1/2;n;1/2) =\frac{1}{2}+\sqrt{2}\
\frac{{}_2F_1(1,-1/2;n;1/2)}{B(n,1/2)}.$$
\end{cor}
\proof
Observe that \cite[9.111]{GR} $$\sqrt{2}\ {}_2F_1(1,-1/2;n;1/2) =
(n-1)\int_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 \ref{Engeneratingfunction}.

\begin{cor}
\label{series-cor}
We have the asymptotic series
$$E_n=\sqrt{\frac{2n}{\pi}} \left( 1-\frac{3}{8n}-\frac{3}{128n^2}
+\frac{27}{1024n^3}+\frac{499}{32768n^4}+
O(\frac{1}{n^5}) \right) + \frac{1}{2}$$
as $n \rightarrow \infty$.
\end{cor}

\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 $\Gamma(n+1/2)/\Gamma(n)$ can  be found in
\cite[43:6:10]{SO}.


\begin{cor}
\label{doublefactEn}
If $n$ is even, $$E_n= \sqrt{2} \sum_{k=0}^{n/2-1}
\frac{(4k-1)!!}{(4k)!!},$$ while if $n$ is odd, $$E_n= 1+\sqrt{2}
\sum_{k=1}^{(n-1)/2} \frac{(4k-3)!!}{(4k-2)!!}.$$
\end{cor}
\proof
Using Corollary \ref{hyperEn}, the Gauss recursion formulas for
Gaussian hypergeometric functions give 
$$\frac{1}{\sqrt{2}}(E_n- E_{n-2})=
\frac{\Gamma(n-3/2)}{\sqrt{\pi}\Gamma(n-1)} =
\frac{(2n-5)!!}{(2n-4)!!}.$$ Thus it is elementary to establish this
corollary by induction.
\begin{cor}
For $n>1$,
$$E_n=\frac{1-(-1)^n}{2} + \sqrt{2}\ P_{n-2}^{(1-n,\frac{3}{2})}(3).$$
\end{cor}
\proof
The Jacobi polynomials are Gaussian hypergeometric functions. To be
precise
\cite[8.962.1]{GR},
$$ P_{n-2}^{(1-n,\frac{3}{2})}(z) =(-1)^n \
\frac{4\Gamma(n+1/2)}{3\sqrt{\pi}\Gamma(n-1)}
\ {}_2F_1(2-n,3/2;5/2;(z+1)/2).$$
Rewrite the postulated
$E_n$ using this formula, and then proceed as in Corollary
\ref{hyperEn}, or as in Corollary \ref{doublefactEn}.

\section{Real Generalized Eigenvalues}

A ``generalized eigenvalue'' of the pair of matrices $(M_1,M_2)$ (or
of the {\em pencil} $M_1-\lambda M_2$), is defined to be a solution
$\lambda$ to the equation $$\det( M_1 - \lambda 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.

\begin{theorem}
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{\pi}\Gamma((n+1)/2)}{\Gamma(n/2)}.$$
\end{theorem}

Since the asymptotic series of the Euler Beta function is known
\cite[43:6:10]{SO} we have an immediate corollary.
\begin{cor}
We have the asymptotic series
$$E_n^G=\sqrt{\frac{\pi n}{2}} \left( 1-\frac{1}{4n}+\frac{1}{32n^2}
+\frac{5}{128 n^3}-\frac{21}{2048n^4}+
O(\frac{1}{n^5}) \right) $$
as $n \rightarrow \infty$.
\end{cor}

\begin{theorem}
If
$\lambda$ denotes a real generalized eigenvalue of a pair of
independent random matrices, then its probability density $f^G(\lambda)$
is given by $$ f^G (\lambda)=\frac{1}{\pi ( 1 + \lambda^2 )}, $$ that is,
$\lambda$ obeys the standard Cauchy distribution.  Equivalently,
$\mbox{atan}(\lambda)$ is uniformly distributed on $[-\frac{\pi}{2},
\frac{\pi}{2}]$.
\end{theorem}


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)$.
{\it
Thus
the expected number of real eigenvalues of an $n$ by $n$ Cauchy matrix
is equal to $$\frac{\sqrt{\pi}\Gamma((n+1)/2)}{\Gamma(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.

\begin{defn}
Let $\Delta_n$ denote the set of all
$n$ by $n$  singular matrices of Frobenius
norm one.
\end{defn}
Following standard notation, the Frobenius norm of a matrix $A$ is
defined as $\|A\|_F \equiv \sqrt{\sum a_{ij}^2}$.  In the language of
algebraic geometry, $\Delta_n$ is a real algebraic subvariety of
dimension $n^2-2$ of the unit sphere $S^{n^2-1}$ in $\Re^{n^2}$.
Now let $(M_1,M_2)$ be a pair of  matrices.
The intersection in $\Re^{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 $\Delta_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 $\Delta_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.  


\begin{lemma}
(Poincar\'e) 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}$.
\end{lemma}
This formula and its generalizations appear in integral geometry books
such as
\cite{Sant}.
Poincar\'e's formula reduces the problem of calculating
the expected number of real generalized eigenvalues to
finding the volume of $\Delta_n$.

The set $\Delta_n$ was studied by Demmel \cite{de1} and Edelman
\cite{Ed4} 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
$\kappa_D(M) \equiv \|M\|_F\|M^{-1}\|_2$.  Computing the volume of
$\Delta_n$ reduces to computing the asymptotics of the probability
that $\kappa_D> \alpha$ as $\alpha \rightarrow \infty$:

Let $V_\epsilon(\Delta_n)$ be the volume of an $\epsilon$ neighborhood
of $\Delta_n$ in $S^{n^2-1}$.  Clearly, $$ \mbox{Vol}(\Delta_n) =
\lim_{\epsilon \rightarrow 0} (2\epsilon)^{-1} V_\epsilon(\Delta_n).$$
By the definition of the Demmel condition number $\kappa_D$,
$$V_\epsilon(\Delta_n) = \mbox{Prob}[\kappa_D > 1/\epsilon]
\mbox{Vol}(S^{n^2-1}).$$
Edelman
\cite[Corollary 2.2]{Ed4} has shown that
$$\lim_{\epsilon \rightarrow 0} \epsilon^{-1}
\mbox{Prob}[\kappa_D > 1/\epsilon] = \frac{2\Gamma((n+1)/2) \Gamma(n^2/2)}
{\Gamma(n/2) \Gamma((n^2-1)/2)} .$$
We conclude that
$$\mbox{Vol}(\Delta_n) =
\frac{2\pi^{n^2/2}\Gamma((n+1)/2)}{\Gamma(n/2) \Gamma((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(\theta) M_1 -
\sin(\theta) M_2 ] = 0.$$
Since
each of the bivariate normals is invariant under rotation, we can
readily see that $(\cos(\theta),\sin(\theta))$ is uniformly
distributed on the unit circle.  Since $\lambda = \tan(\theta)$,
we have immediately that $\lambda$ is Cauchy.





\section{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/\pi}$.

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:

\begin{center}
{\bf Expected number of real eigenvalues: CM-5 Experiments using
LAPACK on 64 processors}
\end{center}
\begin{center}
\begin{tabular}{|c|c|r|r|r|} \hline
$n$ & trials & experimental $E_n$ & theoretical $E_n$ & minutes \\
\hline 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 \\ \hline
\end{tabular}
\end{center}

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} \mbox{ is} \left\{ \begin{array}{cc}
\mbox{normally distributed} & i\le j \\
\mbox{distributed like $\chi_{n-j}$} & i=j-1\\
0 & {\mbox{otherwise}}
\end{array} \right.
$$ 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 \cite{GVL}.  The subdiagonal is the length of the
column below it which is a $\chi$ distribution, the appropriate
elements are zeroed out creating Hessenberg form, and the remainder of
the matrix remains normally distributed because of the orthogonal
invariance.


\section{Extensions to Other Distributions}

Mehta \cite[Conjectures 1.2.1 and 1.2.2]{Me} 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 \rightarrow \infty$.  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.

\begin{center}
{\bf Expected number of real eigenvalues:
CM-5 Experiments using LAPACK on 128 processors }
\end{center}
\begin{center}
\begin{tabular}{|c|c|r|r||r|c|r|r|r|} \hline
&\multicolumn{3}{|c||}{uniform distribution [-1,1]} &
\multicolumn{3}{|c|}{discrete distribution \{-1,1\}} \\ \hline
$n$ & trials & experimental $E_n$    & minutes 
& trials & experimental $E_n$ & minutes \\ \hline
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 \\ \hline
\end{tabular}
\end{center}





\newpage

\noindent
{\Large \bf 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.


\def\\{\medskip\noindent\hangindent .3in{}}
\setlength{\baselineskip 11pt}
\frenchspacing
\begin{thebibliography}{99}

\bibitem{Abr}
M.~Abramowitz and I.A.~Stegun,
{\sl Handbook of Mathematical Functions},
Dover Publications, New York, 1965.

\bibitem{BR}
A.T.~Bharucha-Reid and M.~Sambandham,
{\sl Random Polynomials},
Academic Press, New York, 1986.

\bibitem{de1}
J.W.~Demmel,
The probability that a numerical analysis problem is
difficult,
{\sl Math. Comp. 50} (1988), 449--480.

\bibitem{demck}
J.W.~Demmel and A.~McKenney,
A test matrix generation suite,
Argonne National Lab, MCS-P69-0389 and LAPACK working note 9.
Available from \verb+netlib@na-net.ornl.gov+ or \verb+xnetlib+.

\bibitem{Ed}
A.~Edelman,
Eigenvalues and condition numbers of random matrices,
{\sl SIAM J. Matrix Anal. Appl. 9} (1988), 543--560.

\bibitem{Ed3}
A.~Edelman,
{\sl Eigenvalues and Condition Numbers of Random Matrices},
PhD thesis, Department of Mathematics, MIT, 1989.

\bibitem{Ed2}
A.~Edelman,
The distribution and moments of the smallest eigenvalue of a random
matrix of Wishart type,
{\sl Linear Alg. Appl. 159} (1991) 55--80.


\bibitem{Ed4}
A.~Edelman,
On the distribution of a scaled condition number,
{\sl Mathematics of Computation 58} (1992), 185--190.

\bibitem{Ed5}
A.~Edelman,
Random matrix eigenvalues meet numerical linear algebra,
{\sl SIAM News 24} (November 1991), 11.

\bibitem{Edbib}
A.~Edelman,
Bibliography of Random Eigenvalue Literature,
available electronically
by anonymous FTP from \verb+math.berkeley.edu+
in the directory  \verb+/pub/edelman+.


\bibitem{Ginibre}
J.~Ginibre,
Statistical ensembles of complex, quaternion and real matrices,
{\sl J. Math. Phys. 6} (1965), 440--449.

\bibitem{Gi2}
V.L.~Girko,
Circular law,
{\sl Theory Prob. Appl. 29} (1984), 694--706.



\bibitem{Gi3}
V.L.~Girko,
{\sl Theory of Random Determinants,}
Kluwer Academic Press, Boston, 1990.

\bibitem{GVL}
G.H.~Golub and C.F.~van Loan,
{\sl Matrix Computations, Second Edition},
Johns Hopkins University Press, Baltimore, 1989.

\bibitem{GR}
I.S.~Gradshteyn and I.M.~Ryzhik
{\sl Table of Integrals, Series and Products, Corrected and Enlarged Edition},
Academic Press, New York, 1980.


\bibitem{GuR}
R.D.~Gupta and D.S.P.~Richards,
Hypergeometric functions of scalar matrix argument are expressible
in terms of classical hypergeometric functions,
{\sl SIAM J. Math. Anal. 16} (1985), 852--858.



\bibitem{hwang}
C.R.~Hwang,
A brief survey on the spectral radius and the spectral
distribution of large random matrices with i.i.d. entries,
{\sl Random Matrices and their Applications},
Contemporary Mathematics {\bf50}, American Mathematical Society, Providence,
1986, pp. 145--152.\smallskip

\bibitem{kac1}
M.~Kac,
On the average number of real roots of a random algebraic equation,
{\sl Bull. Am. Math. Soc. 49} (1943), 314--320 and 938.

\bibitem{kac2}
M.~Kac,
On the average number of real roots of a random algebraic
equation (II),
{\sl Proc. London Math. Soc. 50} (1948), 390--408.

\bibitem{ko0}
E.~Kostlan,
Statistical complexity of dominant eigenvector calculation,
{\sl Linear Alg.  Appl. 162--164} (1992), 385--388.

\bibitem{ko}
E.~Kostlan,
On the distribution of the roots of random polynomials,
in {\sl
{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.


\bibitem{sommers1}
N.~Lehmann and H.-J.~Sommers,
Eigenvalue statistics of random real matrices,
{\sl Phys. Rev. Let. 67} (1991), 941--944.


\bibitem{lui}
S.H.~Lui,
private communication, 1992.

\bibitem{Me}
M.L.~Mehta,
{\sl Random Matrices,}
Academic Press, New York, 1991.

\bibitem{Mu}
R.J.~Muirhead,
{\sl Aspects of Multivariate Statistical Theory,}
John Wiley \& Sons, New York,
1982.

\bibitem{Sant}
L.A.~Santal\'o,
{\sl Integral Geometry and Geometric Probability,}
Vol ume 1 of
{\sl Encyclopedia of Mathematics and Its Applications,}
Addison-Wesley, Reading,
1976.

\bibitem{SS}
M.~Shub and S.~Smale,
Complexity  of Bezout's Theorem II: Volumes and
Probabilities,
{\sl  Proceedings of MEGA 92}, to appear.

\bibitem{sommers2}
H.-J.~Sommers, A. Crisanti, H. Sompolinsky, and Y. Stein,
Spectrum of large random asymmetric matrices,
{\sl Phys. Rev. Let., 60} (1988), 1895--1898.

\bibitem{SO}
J.~Spanier and K.B.~Oldham,
{\sl An Atlas of Functions,}
Hemisphere Publishing, Washington, 1987.

\end{thebibliography}


\end{document}
