\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{\ha}{\frac{1}{2}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\str}{\rule[-.45cm]{0cm}{0.9cm}}

\newtheorem{con}{Conjecture}[section]
\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{\per}{ {\mbox{\rm per}} }

\newcommand{\xh}{\hat{x}}
\newcommand{\yh}{\hat{y}}
\newcommand{\rh}{\hat{\rho}}
\newcommand{\dH}{ {\rm \bf dH}}
\newcommand{\E}{{\rm \bf E}}
\newcommand{\dZ}{ {\rm \bf dZ}}
\newcommand{\dM}{ {\rm \bf dM}}
\newcommand{\dR}{ {\rm \bf dR}} 
\newcommand{\dQ}{ {\rm \bf dQ}} 
\newcommand{\dA}{ {\rm \bf dA}} 
\newcommand{\su}{ \mbox{$\sum$}}
\newcommand{\erfc}{{\rm erfc}}
\newcommand{\pro}{\mbox{$\prod$}}
\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{
The  Circular Law and
the Probability that a Random 
Matrix Has $k$ Real Eigenvalues
}
\author{\sl{Alan Edelman} 
\\ Department of Mathematics \\
University of California \\
Berkeley, California 94720 \\
\verb+edelman@math.berkeley.edu+}
\date{March 23, 1993}

 
\maketitle
 
\begin{abstract}
Let $A$ be an $n$ by $n$ matrix  whose elements
are independent random variables with 
standard normal distributions.
Girko's controversial circular law states that the
distribution of appropriately normalized eigenvalues 
is asymptotically uniform in the unit disk in the complex
plane.  We derive the exact distribution
of the complex eigenvalues  for finite $n$, 
from which the 
circular law  for normally distributed matrices
may be derived.

Similar methodology allows us to 
derive a joint distribution formula for the real Schur
decomposition of $A$. Integration of this distribution
yields the probability that $A$ has exactly $k$ real eigenvalues.
For example, we show that the probability that
$A$ has all real eigenvalues is exactly $2^{-n(n-1)/4}$.
\end{abstract}

\section{Introduction}

This paper investigates the eigenvalues of a real random 
$n$ by $n$ matrix
of standard normals.
We have two questions to explore:

\begin{itemize}
\item What is the probability
$p_{n,k}$ that exactly  $k$ eigenvalues are real?
\item Why do the complex  eigenvalues when properly normalized roughly fall
uniformly 
in the unit disk as in Figure 1?  
\end{itemize}

Both questions can be answered by first factoring the matrix into 
some form of the real Schur decomposition, then interpreting this
decomposition as a change of variables, and finally performing a
wedge product derivation of the Jacobian of this
change of variables.  We demonstrate the power
of these techniques by obtaining  exact answers to these questions.
Previous investigations of thse issues have  lead to
inconclusive  results \cite{Ginibre} and very cumbersome
and probably  erroneous results
\cite{Gi3}.

This paper may be thought of as a sequel to \cite{Ed6},
where we answered the questions
\begin{itemize}
\item What is the expected number $E_n=\sum_k k p_{n,k}$
of real eigenvalues?
\item Why do the real eigenvalues when properly normalized
roughly fall uniformly in the interval $[-1,1]$?
\end{itemize}
In fact, our investigation into $p_{n,k}$ preceded that
of \cite{Gi3}, but when we saw that $\sum kp_{n,k}$, always
had a particularly simple form, 
we diverted our attention to understanding $E_n$ and
related issues. We felt that
the derivation of $E_n=\sum k p_{n,k}$  must somehow be simpler
than the derivation of the individual $p_{n,k}$.


Random eigenvalue experiments are irresistible given the availability of
modern interactive computing packages.
Numerical experiments beckon us for theoretical explanation
as surely as any experiment in mechanics did centuries ago.
If understanding is not sufficient motivation,
random matrix eigenvalues arise in models of nuclear physics
\cite{Me}, in multivariate statistics \cite{Mu}, and in other areas
of pure and applied mathematics including numerical analysis
\cite{Gi3,Edbib}.\footnote{Most
modern  papers on numerical linear algorithms, rightly or
wrongly, contain numerical experiments on random matrices.}


Perhaps we are still at the tip of the iceberg in our understanding
of eigenvalues of random matrices.
Few papers contain exact distributional results for non-symmetric
random matrices.  Given an exact distribution formula for a normally
distributed matrix, one expects such a formula to be asymptotically
valid for matrices of elements of mean $0$ and variance $1$.  This
is a central limit effect. 
At the present time, however, there is no
satisfying  theory that seems to allow
us to make the leap from normally distributed matrices to a wider
class.


\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}
}

The contrast between the real and complex eigenvalues
of a normally distributed matrix is illustrated in
Figure 1.  This figure plots the 2500 values of $\lambda/\sqrt{n}$
where $\lambda$ is an eigenvalue of the random matrix of dimension
$n=50$.   
Notice that the complex normalized eigenvalues may appear
to be roughly uniformly distributed in the unit disk.
This is 
Girko's controversial
circular law
\cite{Gi2},
which states that 
if the elements of a random matrix are independent with mean $0$ 
and variance $1$, then the distribution of the normalized eigenvalues is
asymptotically uniformly  distributed over the disk.
Nobody doubts the overwhelming experimental evidence;
the controversy is over whether this law has yet
been proved mathematically.

Many have tried
to understand Girko's ``proofs,'' without success.
This author has also not been able to verify Girko's
proof, though we are not sure if the errors are in the translation
from the original Russian.  We have contacted Girko on this
matter and eagerly await his response.

In this paper we rigorously prove a result that is 
at the same time  weaker 
and stronger than Girko's law.  We derive the exact distribution of
the complex eigenvalues of a real normally distributed matrix.
The limit of  the distribution of the normalized eigenvalues 
is indeed the uniform distribution on the circular disk.
This is weaker than Girko's law in the sense that our results
only apply for the normal distribution.   Intuitively, the normal
distribution models all distributions with elements of mean 0
and variance 1,  by central limit effects, but this sort of reasoning
alone is not rigorous mathematics.  Our result is stronger than
Girko's law in that it gives exact distributions for finite $n$.

A much simpler case occurs when the random matrix comes 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 reported by Mehta \cite[p.300]{Me}
and Hwang \cite{hwang}.



Two recent books \cite{Me,Gi3}
report on previous investigations concerning
the eigenvalues of a random Gaussian matrix with no symmetry conditions
imposed.
Our approach improves on the  work of Dyson 
\cite[Appendix 35]{Me},      
Ginibre  \cite{Ginibre}, and
Girko\cite[Theorem 3.5.1]{Gi3}.
Ginibre and Girko compute  expressions for the joint distribution
by diagonalizing $A=X\Lambda X^{-1}$.
Girko's form is in terms of the real eigenvalues and the real and
imaginary parts of the complex eigenvalues.
Dyson computed an
expression for the special case
when $A$ has all real eigenvalues using the real
Schur decomposition \cite{Murn} $A=QRQ^T$, where $R$ is upper triangular.
Girko further expresses the probability $c_k$ that $A$ 
has $k$ real eigenvalues as the solution to an $n$ by $n$ linear
system of equations, where the coefficients are multivariate integrals. 
Girko's approach is rather cumbersome, and does not seem to be correct,
at least in the English version.

After preparation of this manuscript, we learned that
Theorem 6.1 was independently discovered by Lehmann and
Sommers in \cite{ls} who consider the range of problems
from the completely symmetric case to the completely
antisymmetric case.
Our derivation follows a different viewpoint.

Our approach is simplified by directly considering the real Schur 
decomposition \cite{Murn} even when the matrix  does not have all real
eigenvalues. 
It is well known to researchers in matrix computations that
orthogonal decompositions are of great value \cite{GVL}.
It is perhaps not surprising that such decompositions are also 
appropriate for random matrix theory.


\section{Main results}


\boxer{
{\sc Probability of exactly $k$ real eigenvalues} \\
The probability $p_{n,k}$ that
a random  $A$ has $k$ real eigenvalues
has the form $r+s\sqrt{2}$, where $r$ and $s$ are rational.
In particular, the probability that a random matrix
has all real eigenvalues is
$$p_{n,n}=1/2^{n(n-1)/4}.$$
\hfill Cor \ref{cor71} and \ref{cor72}
}

In the tables  below, 
we list the probabilities for $n=1,\ldots,9$ both exactly and 
to five decimal places.


\renewcommand{\arraystretch}{1.3}
\begin{center}
\begin{tabular}{|cc|c|l|}
\hline
$n$ & $k$ & \multicolumn{2}{c|}{ $p_{n,k}$ }  \\ \hline
1
 &1& $ 1 $ &1\\
\hline
2
 &2& $\frac{1}{2} \sqrt{2} $ &0.70711\\
 &0& $ 1 -\frac{1}{2} \sqrt{2} $ &0.29289\\
\hline
3
 &3& $\frac{1}{4} \sqrt{2} $ &0.35355\\
 &1& $ 1 -\frac{1}{4} \sqrt{2} $ &0.64645\\
\hline
4
 &4& $ \frac{1}{8}  $ &0.125\\
 &2& $-\frac{1}{4} +\frac{11}{16} \sqrt{2} $ &0.72227\\
 &0& $ \frac{9}{8} -\frac{11}{16} \sqrt{2} $ &0.15273\\
\hline
5
 &5& $ \frac{1}{32} $ &0.03125\\
 &3& $-\frac{1}{16} +\frac{13}{32} \sqrt{2} $ &0.51202\\
 &1& $ \frac{33}{32} -\frac{13}{32} \sqrt{2} $ &0.45673\\
\hline
6
 &6& $ \frac{1}{256} \sqrt{2} $ &0.00552\\
 &4& $ \frac{271}{1024} -\frac{3}{256} \sqrt{2} $ &0.24808\\
 &2& $-\frac{271}{512} +\frac{107}{128} \sqrt{2} $ &0.65290\\
 &0& $ \frac{1295}{1024} -\frac{53}{64} \sqrt{2} $ &0.09350\\
\hline
\end{tabular}
\hspace{.5in}
\renewcommand{\arraystretch}{1.4}
\begin{tabular}{|cc|c|l|}
\hline
$n$ & $k$ & \multicolumn{2}{c|}{ $p_{n,k}$ }  \\ \hline
7
 &7& $ \frac{1}{2048} \sqrt{2} $ &0.00069\\
 &5& $ \frac{355}{4096} -\frac{3}{2048} \sqrt{2} $ &0.08460\\
 &3& $-\frac{355}{2048} +\frac{1087}{2048} \sqrt{2} $ &0.57727\\
 &1& $ \frac{4451}{4096} -\frac{1085}{2048} \sqrt{2} $ &0.33744\\
\hline
8
 &8& $ \frac{1}{16384}  $ &0.00006\\
 &6& $-\frac{1}{4096} +\frac{3851}{262144} \sqrt{2} $ &0.02053\\
 &4& $ \frac{53519}{131072} -\frac{11553}{262144} \sqrt{2} $ &0.34599\\
 &2& $-\frac{53487}{65536} +\frac{257185}{262144} \sqrt{2} $ &0.57131\\
 &0& $ \frac{184551}{131072} -\frac{249483}{262144} \sqrt{2} $ &0.06210\\
\hline
9
 &9& $ \frac{1}{262144} $ &0.00000\\
 &7& $-\frac{1}{65536} +\frac{5297}{2097152} \sqrt{2} $ &0.00356\\
 &5& $ \frac{82347}{524288} -\frac{15891}{2097152} \sqrt{2} $ &0.14635\\
 &3& $-\frac{82339}{262144} +\frac{1345555}{2097152} \sqrt{2} $ &0.59328\\
 &1& $ \frac{606625}{524288} -\frac{1334961}{2097152} \sqrt{2} $ &0.25681\\
\hline
\end{tabular}
\end{center}
\renewcommand{\arraystretch}{1.0}

\boxer{
{\sc Joint Eigenvalue Density given $k$ real eigenvalues} \\
The joint density  of the
ordered  real eigenvalues $\lambda_j$ and ordered 
complex eigenvalue pairs $x_j \pm iy_j$, $y_j>0$
given that $A$ has $k$ real eigenvalues is
$$
\frac{2^{l-n(n+1)/4}}{
\prod_{i=1}^n \Gamma(i/2)}
\Delta
 \exp \left(
\su(y_i^2-x_i^2)
- \su \lambda_i^2/2 \right)
\pro  \mbox{\rm  \erfc($y_i \sqrt{2}$)} ,
$$
where $\Delta$ is the 
magnitude of the  product of the differences of the eigenvalues
of $A$.
Integrating this formula over the $\lambda_j$, $x_j$ and $y_j>0$
gives $p_{n,k}$. \hfill Theorem \ref{theorem61}
}


\boxer{
{\sc Density of Non-real Eigenvalues} \\
The density of a random complex eigenvalue of a normally distributed
matrix is
$$\rho_n(x,y)=\sqrt{2/\pi} y 
e^{y^2-x^2}
\erfc(y \sqrt{2})
e_{n-2}(x^2+y^2),
$$
where $e_n(z)=\sum_{k=0}^n z^k/k!$.
Integrating this formula over the upper half plane
gives half the expected number of non-real eigenvalues.
\hfill Theorem \ref{theorem62}
}

Notice the factor $y$ in the density indicating a low 
density near the real axis.  Readers may notice 
the white space immediately above and below the
real axis in Figure 1.  We think of the real axis as
attracting  these nearly real eigenvalues.






\boxer{
{\sc The Circular Law }\\
Let $z$ denote a non-real eigenvalue normalized
by dividing by $\sqrt{n}$. As $n\rightarrow
\infty$, $z$ is uniformly distributed in the disk $|z|<1$.
Furthermore, as $n \rightarrow \infty$, each eigenvalue
is almost surely non-real.
\hfill Theorem \ref{theorem63}
}

It is tempting to believe that the eigenvalues are almost
surely not real, because the real line is a set of 
Lebesgue measure zero in the plane.  However this
reasoning is not correct. Indeed it is not true
for finite $n$. When $n=2$, we see that the probability
of real eigenvalues is greater than non-real eigenvalues.
The error arises because of the false intuition that
the density of the eigenvalues is absolutely continuous
with respect to Lebesgue measure.



\section{Random Determinants, Permanents, Derangements,
 and  Hypergeometric Functions
of Matrix Argument}

This section  may be read independently from the rest of the paper.
Here we relate random determinants with the theory
of permanents and hypergeometric functions of matrix argument.
We have not seen this connection in the literature.


We remind the reader that the permanent function of a matrix
is similar to the
usual definition of the  determinant except that the sign of the permutation
is replaced with a plus sign:
$$\per(A)=\sum_{\pi} a_{1\pi(1)}a_{2\pi(2)} \ldots a_{n\pi(n)}.$$
Generally the permanent is
considerably  more difficult to compute than the determinant.

Hypergeometric functions of a matrix argument are less familiar
than permanents.  They arise in multivariate statistics and more
recently in specialized fields of harmonic analysis.  Unlike
the matrix exponential, for example, these functions take matrix
arguments but yield scalar output.  They are 
much more complicated than
merely evaluating a scalar hypergeometric function at the eigenvalues
of the matrix.  Readers unfamiliar with the theory may safely 
browse or skip over allusions to this theory.
Those wanting to know more should
consult  \cite[Chapter 7]{Mu}.


\begin{theorem}
\label{theorem31}
If  $A$ is a random matrix with independent elements of
mean $0$ and variance $1$, and if $z$ is any scalar (even complex!)
then
$$\E \det(A^2+z^2I)=
\E \det(A+zI)^2=\per(J+z^2I)=n!e_n(z^2)
=n! \ _1F_1(-1;\frac{1}{2}n;-\frac{1}{2}z^2I)$$
where \E denotes
the expectation operator,  $J$ is the matrix of all ones, 
$$e_n(x)=\sum_{k=0}^n \frac{x^k}{k!}$$ is the truncated
Taylor series for $\exp(x)$, 
``per'' refers to the permanent function of a matrix,
and the hypergeometric function has a scalar multiple of the
identity as argument.
\end{theorem}

The alert reader might already detect our motivation
for Theorem \ref{theorem31}.  It is the source of
the term $e_{n-2}(x^2+y^2)$ in the density of the 
non-real eigenvalues that will appear again in
Theorem \ref{theorem62}.

Before we prove this theorem from a set of lemmas to follow,
we remark on two special elementary cases. 
Taking $z=0$ is particularly simple.
When $z=0$, the theorem states that the 
expected value of the determinant of $A^2$ is $n!$, and of course
the permanent of a matrix of all ones is $n!$.
Expanding and squaring the determinant of $A^2$, we see
cross terms have expected value $0$ and the other $n!$ squares
have expected value $1$.  


Now consider $z^2=-1$. The theorem states  that
$\E(\det(A^2-I))=n!\sum_{k=0}^n (-1)^k/k!$, an expression which we
recognize as $D_n$, the number of derangements on $n$ objects,
i.e. the number of permutations on $n$ objects that leave no object
fixed in place.  Many combinatorics books
point out that
$\per(J)=n!$ and $\per(J-I)=D_n$
(\cite[p. 28]{ryser}, \cite[p. 161]{aigner},
and \cite[p. 44]{minc}),
but we have not seen a general formula for  $\per(\lambda I -J)$,
the {\it permanental characteristic polynomial} of $J$,
in such books.

Theorem  3.1 is a synthesis of the following lemmas:

\begin{lemma}
\label{lemma31}
Let $A$,$B$, and $C$ be matrices whose elements are the 
indeterminants $a_{ij}$,$b_{ij}$ and $c_{ij}$, respectively.
Any term in the polynomial expansion of
$\det(A^2+B+C)$ that contains a $b_{ij}$ does not contain
an $a_{ij}.$
\end{lemma}

\proof
Let $X_{ij}$ denote the matrix obtained from $A^2+B+C$
by removing the $i$th row and the $j$th column.
Then
$$\det(A^2+B+C)=\pm b_{ij} \det(X_{ij}) +\mbox{ terms independent of
$b_{ij}$.}$$  Since $a_{ij}$ only appears in the $i$th row
or $j$th
column  of $A^2$, the lemma is proved.
\qed

\begin{cor}
Let $A$ be a matrix of independent random elements of mean $0$,
and let  $x$ and $y$ be  indeterminants.
Then
$$\E \det( (A+xI)^2+y^2 I))=\E \det(A^2+(x^2+y^2)I).$$
$$\mbox{\hspace*{1.5in}} = \E \det \left( A + \sqrt{x^2+y^2}I \right)^2.$$
\end{cor}

\proof
We use Lemma \ref{lemma31} by letting $B=2xA$ and $C=(x^2+y^2)I$.
Since,
$$
\E \det((A+xI)^2+y^2 I)=\E \det(A^2+2xA+(x^2+y^2)I),$$
it follows that any term multiplying $b_{ij}=2xa_{ij}$ in the expansion
of the determinant is
independent of $a_{ij}$ and so has expected value $0$.
Therefore the term $2xA$ makes no contribution to the expected
determinant, and may be deleted.
The second equality may be verified similarly.
\qed

\begin{lemma}
If $A$ is a random matrix of independent random elements of
mean $0$ and variance $1$, then
$$\E \det (A+zI)^2 = \per( J+z^2I).$$
\end{lemma}

\proof 
Expand the determinant of $A+zI$ into its $n!$ terms and
square.  The expected value of any cross terms must be $0$.
The expected value of the square of any term is 
$(1+z^2)^{s(\pi)}$, where $\pi$ is the permutation associated
with the term and $s(\pi)$ is the number of elements unmoved
by $\pi$.  The sum over all permutations
is clearly equal to $\per(J+z^2 I)$
because $s(\pi)$ counts the number of elements from the diagonal
in the expansion of the permanent.
\qed

\begin{lemma}
If $A$ is a random matrix with independent  normally distributed
elements then
$$\E \det(A+zI)^2 = n! {}_1F_1(-1;\ha n;-\ha z^2 I).$$
\end{lemma}

\proof
Let $M=A+zI$ and $z$ a real scalar.  
Since the quantity of interest is clearly a polynomial in $z$
the assumption that $z$ is real can be later relaxed.
The random matrix  $W=MM^T$ has
a non-central Wishart distribution \cite[p.\ 441]{Mu}
with $n$ degrees of freedom, covariance matrix $I$,
and noncentrality parameters $\Omega=z^2I$.  
The moments of the noncentral Wishart distribution 
were computed in 1955 and 1963 by Herz and Constantine
in terms of hypergeometric functions of matrix
argument. (See
\cite[Theorem 10.3.7, p.\ 447]{Mu}. )
In general, 
$$\E \det W^r = (\det \Sigma)^r 2^{mr} \frac{\Gamma_m(\frac{1}{2}n+r)}
{\Gamma_m(\frac{1}{2}n)}{}_1F_1(-r;\ha n;- \ha \Omega),$$
where $ \Gamma_m$ denotes the multivariate gamma function.
In our case, we care interested in $\E \det W$ which simplifies to
$$\E \det W= n! {}_1F_1(-1;\ha n;-\ha z^2 I).$$
\qed

We remark that the identity
$${}_1F_1(-1;\ha n;-\ha z^2 I)=e_n(z^2)$$
can be obtained directly from a zonal polynomial expansion
of the hypergeometric function of matrix argument using the
homogeneity of zonal polynomials and the value of certain
zonal polynomials at the identity matrix.
(The key expressions are Formula (18) on p.\ 237 and 
Formulas (1) and (2) on p.\ 258 of \cite{Mu}.)
The only partitions that play any role in the sum are those of the
form $(1,1,\ldots,1,0,\ldots,0)$ with at most $n$ ones.

We further remark that in light of Lemma 3.2, the assumption
of normality may be relaxed to any distribution with mean $0$
and variance $1$.  
Thus we see our first example where the assumption of normality
allows us to derive a formula which is correct in more general
circumstances.


 


\section{Orthogonal Matrix Decompositions}

This section may also be read independently from the remainder
of the paper.  Here we review the real Schur decomposition
in various forms that we will need to perform our calculations.
We begin by pointing out that standard numerical linear algebra
calculations may be best looked at geometrically with a change of  basis.


\subsection{Quaternions and Elementary Numerical Linear Algebra}

We wish to study orthogonal (unit determinant) similarity transformations
of $2$ by $2$ matrices:
\be
\label{sim2by2}
M^\prime=
\left(
\begin{array}{rr}
c & -s \\
s & c\\
\end{array}
\right)
M
\left(
\begin{array}{rr}
c & s \\
-s & c\\
\end{array}
\right) ,
\ee
where $c=\cos \theta$ and $s=\sin \theta$.

Researchers familiar with numerical linear algebra know that
$2$ by $2$ matrices are the foundation of many numerical
calculations.\footnote{Developers of LAPACK and its precursors
spent a disproportionate amount of time on the software for 
$2$ by $2$ and other small matrix calculations.}
However, we are not aware that
the elementary geometry hidden in (\ref{sim2by2}) is ever pointed out.

Consider the following 
orthogonal basis for $2$ by $2$ matrices:
$$I=
\left(
\begin{array}{rr}
 1 & 0  \\
 0 & 1  \\
\end{array}
 \right)
\ \ \
J=
\left(
\begin{array}{rr}
 0 & 1  \\
 -1 & 0  \\
\end{array}
 \right)
\ \ \
K=
\left(
\begin{array}{rr}
 1 & 0  \\
 0 & -1  \\
\end{array}
 \right)
\ \ \
L=
\left(
\begin{array}{rr}
 0 & 1  \\
 1 & 0  \\
\end{array}
 \right) .
$$
Any  $2$ by $2$
matrix $M$ can be written as $M=\alpha I + \beta J+
\gamma K + \delta L$.  This basis is a variation of the
usual quaternion basis.
(The usual quaternion basis is $I$, $J$, $iK$, and $iL$,
where $i=\sqrt{-1}$.  It is well known that real linear combinations
of these latter
four matrices form a noncommutative division algebra.)

In the $I,J,K,L$ basis, it is readily calculated that the similarity
in (\ref{sim2by2}) may be expressed
\be
\label{quat4}
\left(
\begin{array}{c}
\alpha^\prime \\
\beta^\prime \\
\gamma^\prime \\
\delta^\prime
\end{array}
\right)
=
\left(
\begin{array}{rrrr}
1 \\
& 1 \\
&& C & -S \\
&& S & C \\
\end{array}
\right)
\left(
\begin{array}{c}
\alpha \\
\beta  \\
\gamma \\
\delta \\
\end{array}
\right) ,
\ee
where $C=\cos 2\theta$ and $S=\sin 2\theta$.
In the language of numerical linear algebra,
\begin{quote}
a $2$ by $2$ orthogonal
similarity transformation is equivalent to a Givens rotation
of twice the angle applied to a vector in $\Re^4$.
\end{quote}

The set of matrices for which $\beta=0$ is exactly
the  three dimensional space of real $2$ by $2$ symmetric matrices.  
Choosing an angle $2\theta$ in the Givens rotation to zero out $\delta^\prime$
corresponds to a step of Jacobi's symmetric eigenvalue algorithm.
This is one easy way to see the familiar fact that only one angle
$\theta \in [0,\pi/2)$ will diagonalize a symmetric $2$ by $2$ matrix,
so long as the matrix is not a constant multiple of the identity.

A case of more concern to us is choosing an angle $2\theta$ 
in the Givens rotation
to set $\gamma^\prime=0$, in the general nonsymmetric case.
Doing so we conclude:

\begin{lemma}
Any $2$ by $2$ matrix $M$ is orthogonally similar to a matrix with
equal diagonals.  If the matrix 
is not equal to $\alpha I +\beta J$, then there is only one
angle
$\theta \in [0,\pi/2)$ that transforms  $M$ into this form.
\end{lemma}

Since the components in the $I$ and $J$ directions are invariant
under the Givens rotation (\ref{quat4}) we read right off that the
trace of the matrix and the difference of the off-diagonal terms
of the matrix is invariant under orthogonal similarities.
We readily conclude

\begin{lemma}
\label{unique-lemma}
Let $M$ be a $2$ by $2$ matrix with the non-real eigenvalues
$x \pm yi$. There is a unique matrix orthogonally similar to
$M$ of the form
$$Z=\left( \begin{array}{rr}
x & b \\
-c & x 
\end{array}
\right)
,
\
\
bc >0, \ \ b \ge c.
$$
\end{lemma}

Given the matrix
$M$ and its  eigenvalues
$x\pm yi$, $b$ and $c$ may be computed from the
equations $bc=y^2$ and $b+c=M_{12}-M_{21}$.
If $M$ is not of the form $\alpha I + \beta J$, then the orthogonal
matrix
with determinant $1$ is unique up to sign.  Notice that
adding $\pi$ to $\theta$ just changes the sign of the orthogonal
matrix.  Adding $\pi/2$ to $\theta$ has the effect of interchanging
the $b$ and the $c$.

\subsection{Real  Schur Decomposition}

The real Shur decomposition expresses $A$ as orthogonally
similar to an upper quasi-triangular matrix $R$ \cite[p.362]{GVL}.
To be precise $A=QRQ^T$ where 
\be
R=
\left(
\begin{array}{cccccc}
\lambda_1 & \cdots &  R_{1k} & R_{1,k+1} & \cdots & R_{1m} \\
          & \ddots & \vdots  & \cdots   & \cdots & \vdots \\
          &        & \lambda_k & R_{k,k+1} & \cdots & R_{km}\\ 
          &        &           & Z_{k+1}  & \cdots & R_{k+1,m} \\
          &        &           &    & \ddots & \vdots \\
          &        &           &    &        & Z_m \\
    
\end{array}
\right).
\label{schurform}
\ee
$R$ is an $n \times n$ matrix ($m=(n+k)/2$) with blocks
$$
R_{ij}\mbox{ of size }\left\{
\begin{array}{ll}
\mbox{1 by 1} & \mbox{if $i\le k,j \le k$} \\
\mbox{1 by 2} & \mbox{if $i\le k,j> k$} \\
\mbox{2 by 1} & \mbox{if $i>k,j\le k$} \\
\mbox{2 by 2} & \mbox{if $i>k,j> k$} 
\end{array}
\right.
$$
Here $R_{jj}$ is the real eigenvalue $\lambda_j$ of 
$A$ when $j\le k$ and $R_{jj}$ is a 
2 by 2 block
$$Z_j=
\left(
\begin{array}{rc}
x_j & b_j\\
-c_j & x_j
\end{array}
\right),
\ \ bc>0, \ \ b \ge c ,
$$
so that the complex eigenvalues of $A$ are
$x_j \pm y_j i$, where $y_j=
\sqrt{b_jc_j},$
for $j>k$.  
Finally, as indicated in (\ref{schurform}),
$R_{ij}$ is a zero block if $i>j$.


The block structure indicated in
(\ref{schurform}) is quite powerful in that it allows us to
simultaneously handle all the possibilities rather cleanly.
All that one must remember is that an index $i$ refers to a size
of one or two depending on whether $i \le k$ or $i>k$ respectively.


Since matrices with multiple eigenvalues form a set of
Lebesgue measure $0$, so long
as our probability distribution is absolutely continuous with
respect to Lebesgue measure (such as the normal distribution) we 
may disregard the possibility that $A$ may have  multiple eigenvalues.
We can then take the Schur form to be unique
if  
we make some further  (arbitrary) requirements on $R$ such as
$\lambda_1 > \cdots > \lambda_k$, $x_{k+1}>\cdots> x_m$.  
Similarly the orthogonal matrix $Q$ is
generically  unique if, for example,
we assume the first row of $Q$ is positive.
For $j \le k$, it is easy to see by induction that the $j$th column
of $Q$ can be chosen to be of arbitrary sign.
{From } Lemma 4.2, the next two columns of $Q$ are unique up
to sign, if the full matrix $Q$ is to have determinant $1$.
Allowing $Q$ to have determinant $\pm 1$ allows us to simply 
specify that the first row of $Q$ be generically positive.


\subsection{Incomplete Schur Decomposition}

Let a matrix $A$ have  a non-real eigenvalue pair $x\pm yi$.
A numerical algorithm to compute the eigenvalues of $A$ would
deflate out these eigenvalues with an appropriate orthogonal
transformation.  We wish to do the same by performing an incomplete
Schur decomposition.

\begin{defn}
We say that $A=QMQ^T$ is an {\bf incomplete Schur decomposition}
for the matrix $A$ with non-real eigenvalue pair $x \pm yi$
if
$$
A=Q\left(
\begin{array}{c|c}

\begin{array}{rr}
x & b \\
-c & x 
\end{array}
&
W  \\ \hline
\rule[-.15in]{0in}{.4in}0 &  A_1 \\
\end{array}
\right) Q^T, \ \ b \ge c, \ y=\sqrt{bc},
$$
where $A_1$ is an $n-2$ by $n-2$  matrix, and
$Q$ is an orthogonal matrix that is the product of two
Householder reflections as specified below.
\end{defn}

The form of $Q$ is important.  We know 
that generically there  is a unique 
$n$ by $2$ matrix
$H$ with positive first row for which  $$H^TAH=
\left(
\begin{array}{rr}
x & b \\
-c & x 
\end{array}
\right) \ \ b \ge c.
$$
Numerical linear algebra textbooks \cite{GVL} describe
how to construct Householder reflections $Q_1$ and $Q_2$ such
that  $Q_2Q_1H$ is  the first two columns of the identity matrix.
We take our $Q$ in the definition above to be the matrix $Q_2Q_1$.
Notice that $H$ is the first two columns of $Q^T$.

This decomposition will be used in the derivation of the circular
law and complex eigenvalue density of normally distributed matrices.




\section{Jacobian Computation}

The symmetries built into the orthogonal decompositions allow us to
compute the Jacobians of the various factorizations that are of interest
to us.

We begin with a lemma that is familiar to numerical analysts who
discuss the discrete Laplacian operator in more than one dimension:

\begin{lemma}
Let $X$ be an $m$ by $n$ matrix.  Define the linear operator
$$\Omega(X)= XA-BX,$$ 
where $A$ and $B$ are fixed square  matrices of dimension $m$ and 
$n$ respectively.
If $\lambda_A$ is an eigenvalue of $A$ and $\lambda_B$ is an
eigenvalue of $B$, then $\lambda_A-\lambda_B$ is an eigenvalue
of the operator $\Omega$.
\end{lemma}

\proof
We remark that the operator $\Omega$
can be  represented in terms of  the Kronecker product as
$$\Omega=A^T \otimes I-I\otimes B.$$
If $v_A$ is a left eigenvector of $A$, and $v_B$ is a right
eigenvector of $B$, then $v_Bv_A^T$ is an eigenvector of
$\Omega$.
So long as $A$ and $B$ have a full set of eigenvectors,
we have accounted for all of the eigenvalues of $\Omega$.
This restriction is not needed by continuity. \qed





We proceed to compute the Jacobians. 
Matrices and vectors of differential quantities
are in bold face Roman letters so as to distinguish them from the
notation for Lebesgue measure.  Exterior products of differential
quantities are expressed either in math italics or 
in wedge product notation.  Math italics denote the
wedge product over the independent elements of $dX$ without regard to
sign. 
Therefore if $X$ is arbitrary,
{\bf dX} denotes the matrix with elements $dx_{ij}$
and $dX$ denotes the natural element of integration in $\Re^{n^2}$
which could also be expressed as $\wedge_{ij} dx_{ij}$ 
or $\wedge{\rm \bf dX}$.  (A wedge without range also means the
exterior product over  the independent entries.)
If $X$ is diagonal, then $dX=dx_{11}\ldots dx_{nn}$ while
if $X$ is symmetric
or upper triangular  (antisymmetric), $dX$ is a wedge product
over the $n(n+1)/2$ ($n(n-1)/2$) independent elements.


If $Q$ is orthogonal, then the matrix $\dH=Q^T\dQ$
is antisymmetric, so 
$$dH=\bigwedge_{i>j} q_i^Tdq_j$$
is the natural element of integration
(for Haar measure) over the space of
orthogonal matrices.

\subsection{Real Schur Decomposition}

The Schur decomposition given in (\ref{schurform}),
indeed any matrix factorization, may
be thought of as a change of variables.   
Let $l=(n-k)/2$ denote the number of complex conjugate
eigenvalue  pairs, and
$R^U$ denotes the strictly upper triangular part of $R$.
(By this we mean  $R$ with the $\lambda_i$ and $Z_i$ replaced with zeros.)
The $n^2$ independent parameters of $A$ are expressed in terms
of $Q,$ $\Lambda=(\lambda_i)_{i=1}^k$, 
$Z=(Z_i)_{i= k+1}^m$ and $R^U$.  For a matrix with $k$ real
eigenvalues and $l=(n-k)/2$ complex conjugate pairs, the $n^2$
independent parameters are found in the new variables as follows:


\bigskip

\begin{center}
\begin{tabular}{|l|c|} \hline
& parameters \\ \hline
$Q$ & $n(n-1)/2$ \\
$\Lambda$ & $k$ \\
$Z$ & $3l$ \\
$R^U$ & $n(n-1)/2-l$ \\ \hline \hline
$A$ & $n^2$ \\ \hline
\end{tabular}
\end{center}

\bigskip

To obtain the Jacobian of this change of parameters, we express
$dA$ in terms of $dHdR^Ud\Lambda dZ$:


\begin{theorem}
\label{changevars}
Let $A$ be an $n \times n$ matrix written in real Schur form
$A=QRQ^T$. The Jacobian of the change of variables is 
$$dA= 2^l
\Delta_0 
\prod_{i>k} (b_i-c_i)
(dH \ dR^U \
d\Lambda  \ dZ )
,$$
where 
$$\Delta_0=
\prod_{i>j}
|\lambda(R_{ii})-\lambda(R_{jj})|.
$$
Here $\Delta_0$ denotes the absolute product of all the differences of an
eigenvalue of $R_{ii}$ and an eigenvalue of $R_{jj}$,
where $i>j$.   Every distinct pair of eigenvalues of $A$ appears 
as a term in $\Delta_0$ except for complex conjugate pairs.
For reference,
$$ dR^U=\bigwedge_{j>i} dR_{ij}$$ 
an exterior product over $n(n-1)/2-2l$
differentials of the strictly upper
triangular part of $R$; 
$$d\Lambda=d\lambda_1 \ldots d\lambda_k; $$
and 
$$dZ=\bigwedge_{j=k+1}^m dZ_j,$$
a product over $3l$  differentials.

\end{theorem}

\proof
Since $A=QRQ^T$, we know that
\begin{eqnarray*}
\lefteqn\dA  &&=Q\dR Q^T+\dQ RQ^T+QR\dQ^T
\\
  &&=Q(\dR+Q^T\dQ R+R\dQ^TQ)Q^T.
\end{eqnarray*}
Let $$\dH=Q^T\dQ$$ which is antisymmetric and let
$$\dM=\dR+\dH R-R\dH=
Q^T\dA   Q.$$
It is evident that $dA=dM$ since the orthogonal matrices
make no contribution to the Jacobian.
Our goal is to compute
$$dM=
\bigwedge_{i>j} dM_{ij} \bigwedge_{i=j} dM_{ij} 
\bigwedge_{i<j} dM_{ij}.$$

We use the same block decomposition for $H$ and $M$ as we did for $R$.
We begin with the most complicated case. If $i>j$, we have that 
\begin{eqnarray}
\label{ginlike}
\nonumber \lefteqn{
\dM_{ij} } \ &&=
\dH_{ij}R_{jj}-R_{ii}
{\rm  \bf dH}_{ij}+\\
&& \sum_{k<j} \dH_{ik}R_{kj}-\sum_{k>i}R_{ik}
\dH_{kj}.
\end{eqnarray}
The quantities in Equation (\ref{ginlike}) in  bold face are matrices
of differential quantities  with either one or two 
rows and one or two columns.
Following our notation convention, $dM_{ij}$ denotes the
exterior product of the one, two, or four independent differential
quantities in $\dM_{ij}$.
Notice that the  $\dH_{ik}$ and $\dH_{kj}$ inside
the summations have first index greater than $i$ or second index
smaller than $j$. Therefore,
if we order  the blocks by decreasing $i$ and then  increasing $j$,
 we see that
\be
\label{wedged}
\bigwedge_{i>j} dM_{ij} =
\bigwedge_{i>j} \wedge (\dH_{ij}R_{jj}-R_{ii}\dH_{ij});
\ee
the differentials in the summation play no role.


Lemma 5.1 states  that 
$$\wedge (\dH_{ij}R_{jj}-R_{ii}\dH_{ij})=
\prod |\lambda(R_{ii})-\lambda(R_{jj})| dH_{ij}.$$
There are either one, two,or four multiplicands in the product.
Since $R_{ii}$ and $R_{jj}$ are either 1 by 1 or 2 by 2 matrices,
we  can explicitly write the product as
\be
\label{del_prod}
\prod |\lambda(R_{ii})-\lambda(R_{jj})|=
\left\{
\begin{array}{ll}
|\lambda_j-\lambda_i|  & \mbox{if $1 \le j < i \le k$} \\
(\lambda_j-x_i)^2+y_i^2 & \mbox{if $1 \le j \le k < i$}\\
\left((x_j-x_i)^2+(y_j-y_i)^2 \right)
\left((x_j-x_i)^2+(y_j+y_i)^2 \right)
 & \mbox{if $k \le j < i$}
\end{array}
\right.
\ee
Putting this all together we have that 
$$\bigwedge_{i>j} dM_{ij}=\Delta_0 \bigwedge_{i>j} dH_{ij}.$$


We now consider $i=j$. In the following, we let an ellipsis
($\ldots$)
denote terms in $dH_{ij}$ in which $i \ne j$.  Such terms play no
further  role
in the Jacobian.  If  
$i \le k$, then
$$dM_{ii}=d\lambda_i + \ldots.$$ 
If $i>k$ then $dM_{ii}$ is a bit more complicated.
It is easy to see that in this case
$$dM_{ii}=\wedge (\dZ_i+\dH_{ii}Z_i-Z_i\dH_{ii})+\ldots.$$
Notice that since $\dH$ is
antisymmetric,  then if $i>k$, $\dH_{ii}$ has the form
$$\left( \begin{array}{cc}
0 & dh_i \\
-dh_i & 0\\
\end{array}
\right) .
$$
It follows that
\be
\label{e10}
(\dZ_i+\dH_{ii}Z_i-Z_i\dH_{ii})=
\left(
\begin{array}{cc}
dx_i+(b_i-c_i)dh_i & db_i \\
-dc_i & dx_i+(c_i-b_i)dh_i
\end{array}
\right)  .
\ee

With $dZ_i=db_idc_idx_i$, 
we have the exterior product of the elements in (\ref{e10})
is
$2(b_i-c_i)dZ_i dh_i .$
We therefore learn that
$$
\bigwedge_{i=j} dM_{ij} = \prod_{i>k} 2^l(b_i-c_i) d\Lambda dZ
\wedge_{i=k+1}^m dh_i + \ldots.
$$ 



Finally, it is easy to verify that if $i<j$,
$dM_{ij}=dR_{ij}+\ldots.$  Therefore
$$\bigwedge_{i<j} dM_{ij}=dR^U + \ldots  $$ 
completing the proof. \qed

Given the $2$ by $2$ matrix
$$
\left(
\begin{array}{rr}
x & b \\
-c & x
\end{array}
\right) , \ bc >0, \ b \ge c,
$$
let $\delta$ denote $b-c$.  Then
\begin{lemma}
\label{change-lemma}
The Jacobian of the change of variables from $b$ and $c$
to $\delta$ and $y$ is
$$db \ dc = \frac{2y}{\sqrt{\delta^2+4y^2}} \ dy \ d\delta . $$
\end{lemma}
\proof


Since
$$bc=y^2 \mbox{   and   }  b-c =\delta \ge 0,$$
it follows that
$$bdc + cdb=2ydy \mbox{   and   }  db-dc=d\delta.$$
Therefore
$(b+c)dbdc=2y d\delta dy$. The conclusion is derived from
the equation $(b+c)^2=\delta^2+4y^2.$ \qed

\subsection{Incomplete Schur Decomposition}

The incomplete Schur decomposition is also a change of parameters
with  counts indicated in the table below.


\bigskip

\begin{center}
\begin{tabular}{|l|c|} \hline
& parameters \\ \hline
$Q$ & $2n-3$ \\
$x,b,c$ & $3$ \\
$W$ & $2n-4$ \\
$A_1$ & $(n-2)^2$ \\ \hline \hline
$A$ & $n^2$ \\ \hline
\end{tabular}
\end{center}

\bigskip


Let  $H$ denote the first two columns of $Q^T$
as we did in Section 5.2.  The set of all possible $H$ is one quarter 
of the 
$2n-3$ dimensional Stiefel  submanifold of $\Re^{2n}$ 
\cite{Mu}. 
The one quarter arises from that fact that we are assuming
that the first row of $H$ contains positive elements.
The matrix of differentials $\dH$ is antisymmetric in 
its first two rows.  Thus there is one differential element in the
first two rows and $2n-4$ below the first two rows.

The natural element of integration on this
submanifold is $dS=\wedge ( Q^T\dH)$.  We use the notation
$dS$ to remind us that $dS$ is a (higher dimensional) \underline{\bf s}urface
element
of the \underline{\bf S}tiefel manifold.

\begin{theorem}
\label{theorem52}
Let $A$ be an $n \times n$ matrix written in incomplete Schur form
$A=QMQ^T$.  The Jacobian of the change of variables is
$$dA=     2 (b-c) \det
\left((A_1-xI)^2+y^2I \right) \ ( db \ dc \ dx \ dA_1 \ dW \ dS).$$
\end{theorem}

The proof of this theorem is very similar to that of Theorem
\ref{changevars}, though simpler. The determinant 
comes from Lemma 5.1  with matrices of dimension $2$ and $n-2$.
  
\section{Applications to the Normal Distribution}


The Jacobians computed in the previous section may be integrated
to compute densities and probabilities.  In the case of the
real Schur decomposition we obtain the joint density 
of the eigenvalues conditioned on $k$ eigenvalues being real.
We further obtain the probability that $k$ eigenvalues are real.
The incomplete Schur decomposition gives us the density of a complex
eigenvalue.  Its integral is the expected number of complex
eigenvalues.
 

\subsection{Applications of the Real Schur Decomposition}

In this section, we calculate the joint density of the eigenvalues
of $A$ and the probability that $A$ has exactly $k$ real eigenvalues
If the elements $a_{ij}$ of $A$ are independent standard normals.


\begin{theorem}
\label{theorem61}
Let  $\Theta_k$ denote the set of matrices $A$ with exactly $k$ real
eigenvalues.
Let $p_{n,k}$ denote the probability that $A\in \Theta_k$.
The ordered  real eigenvalues of $A$ are denoted
$\lambda_i, i=1,\ldots,k$, while the
$l=(n-k)/2$ ordered complex eigenvalue pairs are denoted
 $x_i \pm y_i\sqrt{-1}, i=k+1,\ldots,m$.
Let $$ c_{n,k}=\frac{2^{2l-n(n+1)/4}}{
\prod_{i=1}^n \Gamma(i/2)}. $$

The joint distribution of the real and complex eigenvalues
given that $A$ has $k$ real eigenvalues is
\be
\label{proden}
p_{n,k}^{-1}c_{n,k} \Delta_0
 \exp \left(
\su(y_i^2-x_i^2)
- \su \lambda_i^2/2 \right)
\pro \left[y_i \mbox{\rm  \erfc($y_i \sqrt{2}$)}\right] ,
\ee
where $\Delta_0$ is as in Theorem 5.1.
The probability that $A \in \Theta_k$ is
\be
\label{proitself}
p_{n,k}=\frac{c_{n,k}}{k!l!}\int
\Delta_0
 \exp \left(\su\lambda_i^2/2
+\su(y_i^2-x_i^2) \right)
\pro \left[y_i \mbox{\rm  \erfc($y_i \sqrt{2}$)}\right] 
d\lambda_1\ldots d\lambda_kdx_1\ldots dx_l dy_1 \ldots dy_l .
\ee


\vspace{-.15in}
\hspace*{0.8in}
$
{\tiny
\begin{array}
{l}x_i \in \Re \\ y_i \in \Re^+ \\ \lambda_i \in \Re
\end{array}}
$
\vspace{.2in}
\end{theorem}

\proof
If the elements of $A$ are independent standard normals then
the joint probability density for $A$ is 
$$
(2\pi)^{-n^2/2}\mbox{\rm etr} (-\frac{1}{2}AA^T) dA,
$$
where
$$\mbox{\rm etr}(X)=\exp(\mbox{\rm trace}(X)).$$

Theorem 5.1 states that 
$$
{(2\pi)^{-\frac{n^2}{2}}
\mbox{\rm etr} (-\frac{1}{2}AA^T) dA \ =} \
\hspace{4in}
$$
\be
\label{tellsus}
(2\pi)^{-\frac{n^2}{2}}
2^l 
\Delta_0 
\prod_{i>k} (b_i-c_i) 
\left(  dH   \right) 
\label{jdmess}
\left( e^{-\sum {\frac{r_{ij}^2}{2}}}dR^U \right)
\left( e^{-\sum \frac{\lambda_i^2}{2}} d\Lambda  \right)
\left( e^{-\sum \left( x_i^2+\frac{b_i^2}{2}+\frac{c_i^2}{2}
 \right) }dZ \right) .
\ee

The integral of  (\ref{jdmess}) over $Q,R^U,\Lambda$ and $Z$
(with the restrictions on these variables  specified in Section
4.2 to make the Schur factorization unique), counts every
matrix in $\Theta_k$ once generically. This
integral is
$$
p_{n,k}=(2\pi)^{-n^2/2}
\int_{A \in \Theta_k}
\mbox{\rm etr} (-\frac{1}{2}AA^T) dA,
$$

To obtain the joint density of the eigenvalues, 
we must integrate out all the variables other than the $\lambda_i$,
the $x_i$ and the $y_i$.  To obtain the probability $p_{n,k}$,
we must integrate over all the variables.

(i)The integral over the orthogonal matrices with positive elements
in the first row is
\be
\label{orthvol}
\int 
dH
=
\mbox{\rm Vol($O(n)$)}/2^n.
\ee
The derivation of the volume of the orthogonal group may be
found in \cite{Mu} to be
\be
\label{orthvol2}
\mbox{\rm Vol($O(n)$)}=\frac{2^n\pi^{n(n+1)/4}}{\prod_{i=1}^n \Gamma(i/2)}.
\ee
The $2^n$ in
(\ref{orthvol}) represents the volume of that part of the
orthogonal group which  has positive first row.

\be
\label{RU}
\mbox{
(ii) The integral over $R^U$ is $(2\pi)^{n(n-1)/4-l/2}.$
} \hspace{3.1in}
\ee

(iii) We make the change of variables from $b_i$ and $c_i$ to
$\delta_i$ and $y_i$ described in Lemma \ref{change-lemma}, and then
we integrate out $\delta_i$. Since $b_i^2+c_i^2=\delta^2+2y^2$,
we see that integrating out $d\delta_i$ amounts to computing
$$2ye^{-y_i^2}\int_{\delta_i=0}^\infty
\frac{ \delta_i e^{-\delta_i^2/2}}
{\left( \delta_i^2+4y^2 \right)^{1/2}}d\delta_i . $$
This integral can be obtained from  \cite[3.362.2, p.\ 315]{GR}
with a change of variables.  It equals
\be
\label{fromerf}
2\sqrt{2\pi}y_ie^{y_i^2}\erfc(y_i2^{1/2}).
\ee



We use horizontal braces to indicate the source of each
term in the combined expression for
the joint density of the real
and complex eigenvalues:
\be
\label{bigmess}
p_{n,k}^{-1} 
\overbrace{\underbrace{\str(2\pi)^{-n^2/2}2^l \Delta_0 
\exp\left(
-\su x_i^2
- \su \lambda_i^2/2
\right)}_{\mbox{(\ref{tellsus})}}}
\overbrace{
\underbrace{ \str \frac{\pi^{n(n+1)/4}}{\prod_{i=1}^n \Gamma(i/2)}}_{\mbox{
(\ref{orthvol}) 
and (\ref{orthvol2})}}
}
\overbrace{
\underbrace{\str (2\pi)^{n(n-1)/4-l/2}}_{\mbox{(\ref{RU})}}
}
\overbrace{
\underbrace{\str 2^l(2\pi)^{l/2} \pro [y_i e^{y_i^2} \erfc (y_i \sqrt{2})]}_{\mbox{
(\ref{fromerf})}}
}
\ee
The $p_{n,k}^{-1}$ is the normalization constant so that the joint
density integrates to $1$.  Therefore (\ref{bigmess}) is a probability
density which simplifies to (\ref{proden}).

Equation (\ref{proitself}) is an exact expression for $p_{n,k}$.
Notice that we removed the ordering of the variables, and compensated
for this by dividing by $k!$ and $l!$.  The integral in
(\ref{proitself}) would separate into univariate integrals if it
were not for the connecting term $\Delta_0$.\footnote{This reminds
us of the electrostatic interactions in the Hamiltonian
for a many electron system in quantum mechanics.  Without these
interactions Schr\"{o}dinger's equation could be solved exactly.}


If we wish to compute $p_{n,k}$ explicitly we must integrate out
the $x_i$, the $y_i$ and the $\lambda_i$ in
(\ref{proitself}).  We now describe the
approach we took to integrate out the $x_i$ and $y_i$.
We postpone discussion of the integration of the $\lambda_i$ to Section 7.


From (\ref{del_prod}) we see that
$\Delta_0$ is a polynomial in  the variables $x_i$ and $y_i$.
We can use the integration formula 
\cite[3.461.2]{GR}
\be
\label{intx}
\int_{-\infty}^\infty x^{n}e^{-x^2} dx=
\frac{(n-1)!!}{2^{n/2}} \sqrt{\pi}, \ \ \  \mbox{$n$ even}
\ee

and also

\begin{lemma}
\begin{eqnarray}
\label{eq1}
\int_0^\infty y^{2n+1} \erfc(y2^{1/2})e^{y^2} dy=
\frac{\Gamma(n+\frac{3}{2})}{\sqrt{\pi}(2n+2)2^{n+1}}{}_2F_1(n+1,n+\frac{3}{2};
n+2;\ha)\\
\label{eq2}
\hspace{1in}= \frac{n!(-1)^n}{2} \left( \sqrt{2}
\sum_{k=0}^n (-1)^k \frac{(1/2)_k}{k!}-1 \right),
\end{eqnarray}
where the Pochhammer symbol $(x)_k$ denotes the product
$(x)_k=x(x+1) \ldots(x+k-1)$.
\end{lemma}

\proof
The expression in (\ref{eq1})
for the integral as a hypergeometric function may
be found in
\cite[6.286.1]{GR}.  We did not find this form particularly
enlightening so we attempted to find another expression.

We outline our approach for obtaining the expression (\ref{eq2})
leaving the details to the reader.  
Replacing the erfc in the integrand with $\frac{2}{\sqrt{\pi}}
\int_{x=y\sqrt 2}^\infty e^{-x^2} dx $ allows us to interchange the order
of integration. A change of variables and
a careful application of \cite[2.321.2]{GR}
and standard Gamma integrals completes the task.

It is also possible, though tedious, to ignore the integral
and directly show equality 
of the hypergeometric function expression and (\ref{eq2}).
The function  $_2F_1(n,n+\ha;n+1;x)$ is obtained from
$_2F_1(1,3/2;2;x)$ by taking $n-1$ derivatives \cite[15.2.2]{Abr}.
We may also show that
$$_2F_1(1,3/2;2;x)=
\frac{2}{x}\left( (1-x)^{-1/2}-1 \right)$$
from
\cite[15.2.21 and 15.1.8]{Abr}. \qed

\subsection{Applications of the Incomplete Schur Factorization}

We now turn to the distribution of the complex eigenvalues
$x+yi$, $y > 0$ of the real matrix $A$.
Let $\rho_n(x,y)dxdy$ 
denote the expected number of complex
eigenvalues  of a random matrix $A$  in an infinitesimal area
$dx \ dy$ of the upper half plane.
More rigorously $\rho_n(x,y)$ is the Radon-Nickodym derivative
of the measure
$$\mu_n(\Omega)=\int_\Omega \rho_n(x,y) dxdy
= {\bf E}_A \{ \mbox{ Number of eigenvalues of $A$ contained
     in $\Omega$} \}
$$
defined on measurable subsets $\Omega$ of the
upper half  plane.
\begin{theorem}
\label{theorem62}
The complex eigenvalue density is
$$
\rho_n(x,y)=\sqrt{\frac{2}{\pi}}y e^{y^2-x^2} \erfc (y \sqrt{2})
e_{n-2}(x^2+y^2),
$$
where $e_n(z)=\sum_{k=0}^n z^k/k!$.
\end{theorem}

\proof
The techniques used in this proof are
very similar to those used in Theorem \ref{theorem61}.
Theorem \ref{theorem52}  states that 

\begin{eqnarray}
\label{incjac}
\lefteqn
{(2\pi)^{-\frac{n^2}{2}}
\mbox{\rm etr} (-\frac{1}{2}AA^T) dA \ =}  \\
\nonumber
&
(2\pi)^{-n^2/2} 2(b-c)
\det \left( (A_1-x)^2 + y^2 I \right)
e^{-(x^2+\frac{b^2}{2}+\frac{c^2}{2})}
\  db \ dc \ dx \
{\rm etr}(-\ha A_1A_1^T) \ dA_1 \
{\rm etr}(-\ha WW^T) dW \ dS.
\end{eqnarray}

(i) The volume of the Stiefel manifold may be found in
\cite{Mu}.  One quarter of the volume of the Stiefel
manifold (because of the sign choice on the first
row of $H$) is
\be
\label{t1}
\frac{(2\pi)^{n-1}}{2 \Gamma(n-1)}.
\ee

\be
\label{t2}
\mbox{(ii) The integral over $W$ is $(2\pi)^{n-2}$. \hspace{3.6in}}
\ee

(iii)
Exactly as before, the $b$ and $c$ components transform into
(\ref{fromerf}).

(iv)
We recognize the integral of
$
\det \left( (A_1-x)^2 + y^2 I \right)
{\rm etr}(-\ha A_1A_1^T) \ dA_1 $.
It is
$$(2\pi)^{(n-2)^2/2} {\bf E}_{A_1} 
\det \left( (A_1-x)^2 + y^2 I \right).$$
This explains our motivation for Section 3.
From Theorem \ref{theorem31}, we learn that  the value of the integral of the
terms containing $A_1$ is
\be
\label{t4}
(2\pi)^{(n-2)^2/2}\Gamma(n-1)e_{n-2}(x^2+y^2).
\ee


Combining terms we have  that $\rho_n(x,y)=$
\be
\label{smallmess}
\overbrace{\underbrace{
\str(2\pi)^{-n^2/2}2e^{-x^2}}_{\mbox{(\ref{incjac})}}}
\overbrace{
\underbrace{ \str \frac{(2\pi)^{n-1}}{2\Gamma(n-1)}}_{\mbox{
(\ref{t1})}}}
\overbrace{
\underbrace{\str (2\pi)^{n-2}
}_{\mbox{(\ref{t2})}}
}
\overbrace{
\underbrace{\str 2\sqrt{2\pi}y e^{y^2}
\erfc(y \sqrt{2})}_{\mbox{(\ref{fromerf})}}
}
\overbrace{
\underbrace{\str
(2\pi)^{(n-2)^2/2}\Gamma(n-1)e_{n-2}(x^2+y^2)
}_{\mbox{
(\ref{t4})}}
} .
\ee
This completes the proof. \qed

\begin{cor}
The expected number of non real eigenvalues of a random matrix
$A$ is
$$c_n= \sum_k (n-k)p_{n,k} = 2\int_{\mbox{\hspace*{-.2in}$\small 
\begin{array}{c}
 \\ x \in \Re  \hspace{.1in} \\ y \in \Re^+
\end{array}$ \hspace{-0.3in}}}
\rho_n(x,y)dx \ dy.
$$
\end{cor}
The factor of $2$ counts the complex conjugate pairs.  
We could proceed by integrating $\rho_n(x,y)$, to compute
$c_n$, but this is not necessary because in \cite{Ed6}
we computed
$$c_n=n-\ha-\sqrt{2} \frac{{}_2F_1(1,-1/2;n;\ha)}{B(n,1/2)},$$
expressing $c_n$ in terms of a hypergeometric function
and an Euler Beta function.  Other equivalent expressions
for $n-c_n$ may be found in \cite{Ed6}.
28~
We now turn to understanding this density in the
context of Girko's circular law (or Figure 1)
as $n \rightarrow \infty$.  It is interesting to normalize the
eigenvalues by dividing by $\sqrt{n}$.  Thus we introduce
$$\xh=x/\sqrt{n}, \ \ \
  \yh=y/\sqrt{n}, \ \ \
  \rh_n(\xh,\yh)=
n \sqrt{\frac{2n}{\pi }}
\yh
e^{n(\yh^2-\xh^2)}
\erfc (\yh \sqrt{2n})
e_{n-2}(n(\xh^2+\yh^2)).
$$
$\rh(\xh,\yh)$ is the density of the expected number 
of
normalized eigenvalues in an infinitesimal
area around $\xh+i\yh$.

\begin{theorem}
\label{theorem63}
The density of the normalized eigenvalues 
converges pointwise  to a very simple form as
$n \rightarrow \infty$:
\be
\label{circlelimit}
\lim_{n \rightarrow \infty}
n^{-1}\rh_n(\xh,\yh)=
\left\{
\begin{array}{ll}
\pi^{-1} & \xh^2+\yh^2 < 1 \\
0 &  \xh^2+\yh^2 > 1
\end{array}
\right.
\ee
Asymptotically, a random eigenvalue is
almost surely complex. When normalized by
dividing by $\sqrt{n}$, the eigenvalue 
converges in distribution to a random variable
that is
uniformly distributed in the unit disk
$\xh^2+\yh^2<1$.
\end{theorem}

\proof
The formulas
\cite[6.5.34]{Abr} and
\cite[7.1.13]{Abr} may be used to
verify (\ref{circlelimit}). The former states that
$$\lim_{n \rightarrow \infty}
\frac{e_n(n(\xh^2+\yh^2))}{e^{n(\xh^2+\yh^2)}}
=
\left\{
\begin{array}{l}
1 \mbox{ for } \xh^2+\yh^2 <1 \\
0 \mbox{ for } \xh^2+\yh^2 >1 \\
\end{array}
\right. 
,$$
while the latter states that 
$$
\lim_{n\rightarrow \infty}
\sqrt{n} \yh e^{2n \yh^2} \erfc( \yh \sqrt{2n})
= (2\pi)^{-1/2}.
$$

It is routine though tedious
to show that the convergence is 
dominated by an integrable function on the plane.
We need the inequality
$$e_n(nz)/\exp(nz)=\frac{1}{n!} \int_{nz}^\infty e^{-t}t^n dt
\le \frac{n^n}{n!} z^n e^{-nz} \frac{z}{z-1}$$
and Stirling's inequality \cite[6.1.38]{Abr}
$$n! \ge \sqrt{2\pi n}n^n e^{-n}.$$
We leave the details to the reader.  By the Lebesgue
dominated convergence theorem, we achieve convergence
in distribution. \qed


\section{Relationship with Symmetric Random Matrices}

We now turn to the task of integrating the $\lambda_i$
variables in  (\ref{proitself}).
The relevant part of the integrand containing $\lambda_i$ has
the form
\be
\label{lamform}
 \prod q(\lambda_i) \pro |\lambda_i-\lambda_j|
 e^{-\sum \lambda_i^2/2},
\ee
where $q(\lambda_i)$ is some polynomial function of  $\lambda_i$
originating in  $\Delta_0$.



This form is closely related to 
the joint density of the $k$ real eigenvalues 
of a random $k$ by $k$ symmetric matrix $S=(A+A^T)/2$, where
the elements of $A$ are normally distributed.  In the
physics literature the probability distribution for $S$ is
known as the Gaussian orthogonal ensemble. The joint probability
density for the eigenvalues
$\lambda_1 \ge \ldots \ge \lambda_k$ of $S$ is well known
\cite{Me,Ed3,Mu}:
$$
\frac{2^{-k/2}}{\prod_{i=1}^k \Gamma(i/2)} \prod |\lambda_i-\lambda_j|
e^{-\sum \lambda_i^2/2}.
$$

Therefore  the integral of (\ref{lamform}) is essentially an expectation
for the determinant of a polynomial in $S$. To be more precise:

\begin{lemma}
$$
\frac{1}{k!}
\int_{\lambda_i \in \Re}  \pro q(\lambda_i) \pro |\lambda_i-\lambda_j|
 e^{-\sum \lambda_i^2/2}
= \left\{ 2^{k/2}
 {\pro_{i=1}^k \Gamma(i/2)} \right\}
  {\bf   E}_S \det (q(S)), $$
where $E$ denotes expectation with respect to the symmetric random
matrices defined above.
\end{lemma}

Dividing by $k!$ allows us to symmetrize the integrand
so that we need not assume the $\lambda_i$ are ordered.

If $q$ is a polynomial with rational coefficients,
then  ${\bf E} \det (q(S))$ must be rational since all the moments
of a standard normal are rational.  We suspect that a better
understanding of ${\bf E}\det( q(S))$ is possible, but we are not
even aware of an exact formula for ${\bf E} \det{S^2}$.
Such an understanding would help us simplify our formula
for $p_{n,k}$. For now we are 
content  with this
application of Lemma 7.1.

\begin{theorem}
\label{theorem71}
$$
p_{n,k}=d_{n,k}^{-1}
\int_{\mbox{\hspace*{-.2in}$\small 
\begin{array}{c}
 \\ x_i \in \Re  \hspace{.1in} \\ y_i \in \Re^+
\end{array}$ \hspace{-0.3in}}}
{\bf E}_S
\left\{ \pro_{i=1}^k \det\left( (S-xI)^2+y^2I\right) \right\}
\Delta_{xy} e^{\su (y_i^2-x_i^2)}\pro \left[ y_i
\erfc (y_i \sqrt{2}) \right]
dx_1 \ldots dx_l dy_1 \ldots dy_l,
$$
where
$$d_{n,k}=
2^{\frac{n(n-3)}{4}+k/2} l! \pro_{i=k+1}^n \Gamma(i/2)
$$
and
$\Delta_{xy}$ denotes the multiplicands
in $\Delta_0$ that do not include a $\lambda$ term.
\end{theorem}



\begin{cor}
\label{cor71}
The probability of all real eigenvalues is
$p_{n,n}=2^{-n(n-1)/4}$.
\end{cor}

\proof When $k=n$ and $l=0$, $d_{n,n}=2^{n(n-1)/4}$
and the integrand is simply $1$. \qed

\begin{cor}
\label{cor72}
The probability $p_{n,k}$ has the form $r+s\sqrt{2}$
where $r$ and $s$ are rational. Furthermore $r$ and $s$ may be expressed
as an integer fractions with denominator equal to a power of $2$.
\end{cor}


\proof
The $\E_S$ term in the integrand is a polynomial
in $x$ and $y$ with integer coefficients.  Before taking
expectations, the determinant is a polynomial in $x$, $y$, and
the $S_{ij}$ with integer coefficients.  We recall that the
$n$th moment of the normal distribution is $(n-1)!!$ if
$n$ is even and $0$ if $n$ is odd.  Integrating out
the $x_i$ using Equation (\ref{intx}) gives a factor
of $pi^{l/2}$ which is canceled by the $\pi^{l/2}$
in $\pro_{k+1}^n \Gamma(i/2)$.  The integration of the
$y_i$ using (\ref{eq2})
leads only to rational numbers of the proper form.
\qed

Theorem  \ref{theorem71} is the basis of our earliest
 Mathematica program
(Apendix A) for computing $p_{n,k}$.  
We first evaluate 
${\bf E}_S \prod_{i=1}^k \det\left( (S-xI)^2+y^2I\right)$
by explicity computing the determinant symbollically in terms
of the elements of $S$ and the $x_j$ and $y_j$.
Then we replace the various powers of the elements of $S$ with
their expectation.  We then integrate out the $x_i$ and $y_i$
by replacing powers using formulas (\ref{intx}) and (\ref{eq2}).

We then modified the program to save some arithmetic.
It is readily shown that a certain tridiagonal matrix
with $\chi$ distributions on the three diagonals has
the same eigenvalue distribution as that of $S$.
Therefore, we may use $T$ in place of $S$ to compute
the expectations of the determinants of polynomials in $S$.
This is our current program in Appendix A.
The program may be summarized by saying that a polynomial
is computed in expanded form,  then powers of variaous variables are
symbolically replaced with the appropriate moments.




\section{Eigenvalue Roulette and
Numerical Simulations}

Fifty five years ago Littlewood and Offord \cite{LO}
proposed a game on random polynomials that they hardly had
the technology to play:

\begin{quote}
Let the reader place himself at the point of view of A in the
following situation. An equation of degree 30 with real coefficients
being selected by a referee, A lays B even odds that the equation
has not more than $r$ real roots. $\ldots$ What smallest possible
value of $r$ will make A safe [from being ruined]?
\end{quote}



 
With the benefit of easy to use interactive computer languages
best exemplified by Matlab, it is trivial to play their
game or our eigenvalue version of the game that we choose
to call ``eigenvalue roulette.''  As an example, when $n=8$, the
``wheel'' contains the numbers $0, 2,4,6$ and $8$ and the 
wheel may be spun simply by executing
\begin{verbatim}
              >> rand('normal'); a=rand(8); sum(abs(imag(eig(a)))<.0001
\end{verbatim}

The \verb+.0001+ is an arbitrary
numerical  tolerance that seems to be unimportant in practice.


As a 
more substantial check of our results,
we executed  LAPACK simultaneously on each of 128
processors of the  CM-5 supercomputer.
The experiments agreed perfectly with
the theory.

\section{Open Problems}

This section contains a number of conjectures that we strongly
suspect to be true.  

\begin{con}
Let $\{M_n\}$ be a sequence of $n$ by $n$ random matrices.
Let $t_n$ denote the number of real eigenvalues of $M_n$.
Then $\lim_{n \rightarrow \infty} t_n/n = \sqrt{2/\pi}$
almost surely.
\end{con}
Intuitively this states that if we plot $k/n$ against 
$np_{n,k}$ for $k=0, \ldots,n$, the points ``converge''
to the degenerate distribution $\delta(x-\sqrt{2/\pi})$.

\begin{con}
An explicit formula for $p_{n,k}$ that is easier to compute
from than \ref{theorem71}  may be obtained.
\end{con}

Our formula allowed us to compute $p_{n,k}$ explicitly when
either $k=n$ or $n<10$, but we have not yet succeeded in computing
other probabilities due to the computational complexity 
of our Mathematica program.  We suspect the program can
be rewritten so as to compute values for say $n=10$ or
$n=11$, but ultimately a better formulation will be needed.

\begin{con}
Girko's circular law
for arbitrary random matrices of iid elements of
mean $0$ and variance $1$ may be derived as a corollary
of Theorem \ref{theorem63} using some kind of central
limit theorem.
\end{con}


\vspace{0.5in}

\noindent
{\Large \bf Acknowledgements}

\vspace{0.1in}

The author would like to thank Eric Kostlan for many 
interesting discussions on these and related issues.

\newpage



\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{aigner}
M.~Aigner,
{\sl Combinatorial Theory},
Springer Verlag, New York, 1979.

\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}, to appear.

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


\bibitem{Ed6}
A.~Edelman, E.~Kostlan, and M.~Shub,
How many eigenvalues of a random matrix are real?,
{\sl J. Amer. Math. Soc.}, to appear, 1994.

\bibitem{Edbib}
A.~Edelman,
Bibliography of Random Eigenvalue Literature,
currently available electronically
by anonymous FTP from \verb+math.berkeley.edu+
in the directory  \verb+/pub/Alan.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{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{ls}
N.~Lehmann and H.-J.~Sommers,
Eigenvalue statistics of random real matrices,
{\sl Phys. Rev. Let. 67} (1991), 941--944.

\bibitem{LO}
J.E. Littlewood, A.C. Offord. On the number of real roots of
a random algebraic equation,
{\sl J. London Math. Soc. 13} (1938), 288--295.

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

\bibitem{minc}
H.~Minc,
{\sl Permanents},
Encyclopedia of Mathematics and its Applications, 6,
Addison-Wesley, 1978.



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

\bibitem{Murn}
F.D.~Murnaghan and A.~Wintner,
A canonical form for real matrices under orthogonal transformations,
{\sl Proc. Nat. Acad. Sci., 17} (1931), 417--420.


\bibitem{ryser}
H.J.~Ryser,
{\sl Combinatorial mathematics,}
Carus mathematical monographs; no. 14,
John Wiley \& Sons, New York,
1963.


\bibitem{sommers}
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{SP}
%R.~Schreiber and B.N.~Parlett,
%Block reflectors: theory and practice,
%{\sl SIAM J. Numer. Analy., 25} (1987), 189--205.

\end{thebibliography}

\newpage
\appendix
\section{Mathematica Program to compute $p_{n,k}$}

\small

\begin{verbatim}
(* Mathematica  program to compute the probability that an n x n random Gaussian matrix has k real eigenvalues.
              Alan Edelman -- UC Berkeley 12/92 *)

n=Input["n? "]; k=Input["k? "]; l=(n-k)/2;

(*A k x k Tridiagonal matrix with w on the diag and v off the diag *)
W=DiagonalMatrix[Table[w[i],{i,k}]]; Do[ (W[[i+1,i]]=v[k-i]; W[[i,i+1]]=v[k-i]), {i,k-1}];

(* Determinant involving x[i] and y[i] *)
II=IdentityMatrix[k];
z =If[ l==0 || k==0, 1, Det[(W-II*x).(W-II*x)+ II*y^2]];
z = Expand[Product[z/.x:>x[i]/.y:>y[i],{i,l}]];
z = z /. w[i_]^n_:>0 /; OddQ[n]; z = z /. w[i_]^n_:> (n-1)!!; z = z /. w[i_]->0;
z = z /. v[i_]^n_ :> 0 /; OddQ[n]; z = z /. v[i_]^n_ :> Gamma[(i+n)/2]/Gamma[i/2]; z = z /. v[i_]-> 0;

Do[(   (*Only introduce those x[ii] and y[ii] that we need so far*)
  z = z*y[ii]* Product[  ((x[j]-x[ii])^2+(y[j]-y[ii])^2)((x[j]-x[ii])^2+(y[j]+y[ii])^2), {j,ii+1,l}] ;
  z=Expand[z];
  (* Integrate x[ii] *)
  z = z /. x[ii]^n_ :> 0 /; OddQ[n];  z = z /. x[ii]^n_ :> (n-1)!!/(2^(n/2));  z = z /. x[ii]->0;
  (* Integrate y[ii] *)
  iy[n_]:=n!(-1)^n/2(Sqrt[2] Sum[ (-1)^k Pochhammer[1/2,k]/k!,{k,0,n}]-1);
  z = z /. y[ii]^n_ :> iy[(n-1)/2];  z = z /. y[ii]    :> iy[0];
),{ii,1,l}]
z=Simplify[ z Sqrt[Pi]^l /( 2^(n(n-3)/4+k/2) Product[Gamma[i/2],{i,k+1,n}] l!)];

\end{verbatim}



\end{document}
