\chapter{History and applications}
\label{chapter-ap}
\section{Numerical analysis}
\label{section-na}
\index{Numerical analysis ff}
\index{Random matrices in numerical analysis ff}

The limitations on solving large systems of equations are 
computer memory and speed.  The speed
of computation, however, is not only measured by clocking hardware;
it also depends on numerical stability, 
and for iterative methods,
on convergence rates.
If we had a machine that performed Gaussian elimination\index{Gaussian
elimination}
on a random $n$ by $n$ matrix
for $n=$ 1,000,000 in 
16-digit
precision with a computation time of one year,
after that year, would we have any correct digits\index{Loss
of precision} at all? 
Would  doubling the precision and a wait of two years suffice?
What about iterative methods?  How long do we need to wait
to obtain convergence?

Such questions are not new. 
In 1946, von Neumann\index{von Neumann
etc ff} and his associates saw $n=100$ 
on the horizon, and hoped to understand the
numerical stability question
mathematically.  Today, perhaps the equivalent number is $n=$ 50,000,
but our understanding, though far better, is still imperfect. 

We review here the 
work of von Neumann and his colleagues 
Bargmann, Goldstine, and Montgomery \cite{NG},
for three reasons:
to set the stage for random condition
numbers, to illustrate historical lessons,   
and to fill a gap suggested by
Wilkinson\index{Wilkinson, J.H.}\cite{Wilkinson-2}: 
\begin{quote}
\ldots although their paper deservedly received widespread
recognition very few have been able to give a concise sketch of
the basis of their analysis.
\end{quote}
Given the difficulty of extracting the various ideas from
their work,\footnote{Wilkinson referred to the paper's ``indigestibility."}
we will recast their ideas in modern notation, 
supplying equation and page references
from \cite{NG}  for the reader wishing to go back
to the original source.

Goldstine and von Neumann's
analysis \cite{NGc,NGd} 
has three components: the choice of algorithm,
the rounding errors, and the implications.  
For their algorithm, they made the unfortunate	
choice of
solving $Ax=b$ by forming 
$\ata$ and solving the positive definite system
$\ata x=A^Tb$ (the ``normal equations").
This would be achieved by
inverting $\ata$
and then forming $M=(\ata)^{-1}A^T$ numerically. We 
know now that this 
algorithm unnecessarily squares the condition\index{Condition Number}
number $\kappa,$ 
and thus is unstable.
In fact, Goldstine and von Neumann recognized the importance
of the condition number\footnote{``\ldots the actual estimate
of loss of precision\index{Loss of precision}
\ldots depends not on $n$ only, but also on the
ratio $l$ of the upper and lower absolute bounds of the matrix. \ldots
It appears to be the `figure of merit' expressing the difficulties caused
by \ldots solving simultaneous equation systems \ldots"
(\cite{NGa}, p.~14).}
(though they did not use this term), 
yet failed to observe 
that squaring it was unnecessary.\footnote{To their credit,
Goldstine and von Neumann did remark, however,
that their reasons for forming $\ata$ ``may not 
be absolutely and permanently valid" (\cite{NGc}, p.~514).
They were merely attempting to work with a model with which
they could carry out the analysis. With this model, they avoided
the issue of pivoting.}
 
Through a meticulous account
of the numerical errors in Gaussian elimination, Goldstine and von Neumann 
\index{Error analysis of Gaussian elimination --von Neumann estimate} 
\index{Gaussian elimination - von Neumann error estimate}
bounded the matrix residual computed in fixed point arithmetic 
by\footnote
{This is a rough average of Equations
(7.5$^\prime$) and (7.5$_I^\prime$) (\cite{NGc}, p.~551).
See the comment after Equation (7.16) for
the justification for taking an average.}
\be
\label{GVN}
\|AM-I\|_2 \le 20 \kappa^2n^2\beta^{-s},
\ee
where 
$\beta$ is the base and $s$ is the number of digits.
The factor $n^2$ arose from a worst case analysis of the accumulation
of errors.
In their later 
paper, they reduced the factor $20n^2$ to  $200n$ by assuming
that the errors accumulate statistically 
(\cite{NGd}, Equation (9.19)).
Thus, they obtained the estimate
\be
\label{GVN-2}
\|AM-I\|_2 \approxlt 200 \kappa^2 n \beta^{-s}.
\ee 
In our summary so far, we have 
deliberately not introduced any random matrices;
the only randomness in \r{GVN-2} is
not in the matrix $A$, but
in the assumption of statistical distributions of errors.

How is $\kappa$ estimated  in 
\r{GVN-2}? This is where
von Neumann and his colleagues introduced the
assumption of 
random matrices\index{Random matrices--in numerical analysis models}
distributed with elements from
independent normals.
Various ``rules of thumb"\footnote{
Sometimes these estimates are referred to as an expectation and sometimes 
as a bound that holds with high, though unspecified, probability.} 
are given 
for $\kappa$ when the matrices are so distributed: 
$n$ (\cite{NGa}, p.14), $\sqrt{10}n$ (\cite{NGb}, p.477), and
$10n$ (\cite{NGc}, p.555).
We will show  in Section 
\ref{section-lim-con}
that the probabilities that
$\kappa$ is less than $n$, $\sqrt{10}n$, and $10n$
converge to (approximately)  $0.02$, $0.44$, and $0.80$ respectively.
Thus the estimates have the right order of magnitude, but
even the final estimate of $10n$
will fail for $20\%$ of random  matrices.

To ``reduce the probabilistic uncertainties to reasonably
safe levels," 
in their published papers \cite{NGc,NGd}, 
Goldstine and von Neumann made the assumptions that 
\be
\label{GVN-3}
\begin{array}{c}
\lmax \le 100n, \\
\lmin \ge 1/100n, \\
\kappa \le 10n,
\end{array}
\ee
where $\lmax$ and $\lmin$ are the 
extreme eigenvalues of $\ata.$
(\cite{NGc}, Equations (7.14$_I^\prime$.a, b, and c), written
in terms of the eigenvalues rather than the singular values.)
The result for $\lmax$ is rigorously proved to hold 
with probability greater than $0.99$ (\cite{NGd}, Equation (8.8)),  
but for the estimate of $\lmin$ (which we will show fails for
20\% of random matrices), 
they refer to a work of V. Bargmann\index{Bargmann see also von Neumann} 
(``Statistical distribution of proper values")
which can not be located.\footnote
{In telephone communications during  April
of 1989, Bargmann and Goldstine\index{Goldstine, see also von Neumann}
informed me 
that there probably is no such document.}

Substituting \r{GVN-3} into \r{GVN-2},
and taking ten digit precision,
we obtain that the right hand side of
\r{GVN-2} is less than unity if $n<79$.
Actually, Goldstine and von Neumann made a further 
approximation (\cite{NGd}, Equation (9.18$^{\prime\prime}$)) 
and obtained as one of their ``bottom lines" the value
$n<86$ (\cite{NGd}, Equation (9.19.b)). 

The principal direct methods today for solving linear
equations are Gaussian elimination\index{Gaussian elimination}
with partial pivoting
and the QR  factorization.
Thanks to the work of Wilkinson\index{Wilkinson}, we no longer use \r{GVN-2}
to estimate the errors, but rather we
bound
the relative error of the computed solution $\bar{x}$ 
\index{Error analysis of Gaussian elimination, Wilkinson bound}
by
\be
\label{Wil-1}
\frac{\|x-\bar{x}\|_2}{\|x\|_2}\le \kappa f(n) \epsilon.
\ee
Here $\epsilon$ is floating-point machine epsilon,
$\kappa$ is the condition number, and $f(n)$  
includes 
the possibly combined effects of
a polynomial and, for Gaussian elimination, the
{\it growth factor.}\index{Growth factor}\index{Gaussian elimination--
growth factor}
Both factors are 
often dismissed as pessimistic and ignorable in practice. 



In 1979, Birkhoff and Gulati\index{Birkhoff, G. and Gulati, S.}
\cite{BG} resurrected 
Goldstine and von Neumann's\index{von Neumann, J.}
model, but incorporated Wilkinson's\index{Wilkinson} bound.
They observed the orthogonal invariance\index{Orthogonal invariance}
of random matrices
with independent normally distributed elements
and performed several numerical experiments with
random number generators and EISPACK.\index{EISPACK}
They further suggested that random matrices with normal
elements are not ill-conditioned enough for testing
Gaussian elimination\index{Gaussian elimination}
and proposed a method that would
be.

Once again drawing from Goldstine and von Neumann's model,
Smale\index{Smale}, in a 1984 seminar, asked 
for
$$L_n=E(\log \kappa),$$
the expected logarithm of the condition number of
$n$ by $n$
matrices with normally distributed
elements. This  question is posed in \cite{Sm} 
alongside 
other random algorithm questions
(such as arise in linear programming)
where the average tends to be far less costly 
than the worst case.

To understand the significance of $\log \kappa,$
remember that
the condition number\index{condition number--definition}
can also be defined as the smallest number
$\kappa$ (depending on $A$ only)
such that for all $x$ and $\delta x$,
if $Ax=b$ and $A(x+\delta x)=b+\delta b$, then
\be
\label{Sm-1}
\frac{\N{\delta x}}{\N{x}}\le \kappa \frac{\N{\delta b}}{\N{b}}.
\ee
As a consequence,
$$(\log{\N{\delta x}}-\log{\N{x}}) - 
  (\log{\N{\delta b}}-\log{\N{b}}) \le \log \kappa.  $$
When taking the logarithm to the base 10, the equation above
has the interpretation that the number of correct digits in $x$
can differ from the number of correct digits in $b$ by at most
$\log \kappa.$ Smale thus calls $\log \kappa$ the
{\em loss of precision}\index{loss of precision}.
For another interpretation of $\log \kappa$, consider \r{Wil-1} 
with $f(n)\approx 1$.  The logarithm of the condition
number
then bounds the exponent in the relative error, giving
a kind of ``Richter scale" 
where a relative  error of size $\epsilon$ gives a loss of precision of $0$,
and every tenfold increase in the relative error increments this measure
by $1$.

In response to Smale's question,
bounds for $L_n$ were obtained by
Kostlan\index{Kostlan, E.} \cite{Ko} and Ocneanu\index{Ocneanu}(see \cite{Sm}),
who found
that the exponential of the expected logarithm 
to the base $e$ (i.e., the geometric mean\index{Geometric mean
--of condition numbers})
satisfies  the following inequalities:
\be
n^{2/3+o(1)} \le \exp(L_n) \le en^{5/2}.
\ee
In Section \ref{section-lim-con},
we will replace their bounds 
with $\exp(L_n)=4.65n+o(n)$. 

The gap between $2/3$ and $5/2$ in the exponents above represents
a world of difference in terms of comparing the average behavior
of the conjugate gradient\index{Conjugate gradient iteration} 
iteration with that of Gaussian elimination\index{Gaussian elimination}.
To see this,
assume the condition number of a random matrix $A$ is order $n^\alpha$,
and we apply conjugate gradient towards the linear system $A^TAx=A^Tb$.
The number of operations per iteration is given by the two matrix vector
multiplications, hence it is order $n^2$, while the number of iterations is
order of the condition number of $A$ which is order $n^\alpha$.
Thus the total number of operations is order $n^{2+\alpha}$. 
If $\alpha$ were $2/3$, then we would know that conjugate gradient
is better than Gaussian elimination for random matrices; 
the reverse would hold if $\alpha$ were $5/2$.  Our conclusion 
that $\alpha=1$ shows that both are $n^3$ operations, and the choice
of method  is more subtle; it
depends on the various constants involved.


Though most often one hears about the condition number\index{Condition number-
of a problem}
in the context
of the sensitivity of solving linear systems, 
in fact every problem (as explained further in numerical analysis
texts such as \cite{Henrici})
has a condition number that expresses the magnification of the
relative output error given the relative input error.
Problems are {\em ill-posed} when the condition number is infinite;
they are {\em ill-conditioned} when the condition number is large.
Demmel\index{Demmel, J. ,ff}
\cite{De} explains that the ill-posed
problems in a given class often form an {\em algebraic variety},
while the ill-conditioned problems form a {\em tubular neighborhood}.
Though Demmel's model is very general, we will explain it for
the case of real and complex $n$ by $n$ matrices,
which can be identified with the Euclidean spaces
$\real^{n^2}$ and $\real^{2n^2}$.
The ill-posed matrices lie on a submanifold 
defined by $\det(A)=0.$
Since $\det(A)$ is a polynomial, algebraic geometers refer to this
manifold as an algebraic variety.
$A$ is ill-conditioned when it is, in 
a sense that can be made precise, close to this manifold.
It naturally follows that ill-conditioned matrices lie in a tubular
neighborhood around the manifold,  i.e.,  the set of points within
a certain distance of the manifold. 
Demmel found it convenient to
define the condition number\index{Condition number--Demmel's}
$\kappa_D=\|A\|_F\N{A^{-1}}$.

With his model, Demmel concluded that for complex matrices
with normally distributed elements 
\index{Condition number distribution--Demmel's}
\be
\label{Dem-1}
\frac{(1-x^{-1})^{2n^2-2}}{2n^4x^2} \le
\mbox{Prob}(\kappa_D \ge x) \le 
\frac{e^2n^5(1+n^2/x)^{2n^2-2}}{x^2}.
\ee
He also computed the large $x$ behavior for a fixed $n$:
\be
\label{Dem-2}
\mbox{Prob}(\kappa_D \ge x)=\frac{n(n^2-1)}{x^2}+o(\frac{1}{x^2}).
\ee
The bounds given in \r{Dem-1} are not very tight.  For example,
for $x \approx n$ and $n$ large we get a lower bound
growing like $e^{-2n}$ and an upper bound growing like
$e^{2n}.$
In Section \ref{section--Demmel}, we
will derive the exact result 
$$P(\kappa_D \ge x)=1-(1-n/x^2)^{n^2-1},$$
from which we can readily derive Demmel's asymptotic result
\r{Dem-2} as a corollary.




When solving linear systems by
Gaussian elimination
with partial (or complete) pivoting,
the term
$f(n)$ in \r{Wil-1}
is always factored as $g(n)p(n)$ where $g(n)$ is the
{\em  growth factor}\index{Gaussian elimination--growth factor}
\index{Growth factor}
and $p(n)$ is a polynomial.
The worst case bound for $g(n)$ 
is $2^{n-1}$ for partial pivoting and is conjectured to be $n$
for complete pivoting, but what about the average $g(n)$? 
For partial pivoting, the growth factors that appear
in practice are much smaller than $2^{n-1}$, but the
this phenomenon represents an important
gap in our understanding of partial pivoting.
Recently, Trefethen\index{Trefethen, L.N.} and Schreiber\index{Schreiber, R.S.}
\cite{ST}
have
studied this question and have given evidence that
the growth factor is another example 
of a quantity whose average is far better than its worst case.
Through careful experiments and the beginnings of a theory,
they found
that the average growth
factors
behave 
approximately
like $n^{2/3}$ for partial pivoting  and $n^{1/2}$
for complete pivoting.
