\documentclass[12pt]{article}
\usepackage{multirow,amssymb,amsmath}

%\textwidth6.25in \textheight8.5in \oddsidemargin.25in \topmargin0in

%\renewcommand{\baselinestretch}{2.0}


\newtheorem{theorem}{Theorem}[section]
%\newtheorem{theorem}{Theorem}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{remark}[theorem]{Remark}

\makeatletter
\renewcommand{\section}{\@startsection {section}{1}{\z@}%
                                   {-3.5ex \@plus -1ex \@minus -.2ex}%
                                   {2.3ex \@plus.2ex}%
                                   {\normalfont\large\bfseries}}
\renewcommand{\subsection}[1]{\@startsection{subsection}{2}{\z@}%
    {-2.5ex\@plus -1ex \@minus -.2ex}%
    {-.0em}
    {\reset@font\bfseries}{\!\!#1}{$ $}}
\makeatother

\def\proof{\bigskip\noindent{\em Proof. }}
\def\qed{\hfill$\square$}
\def\ueq{\stackrel{\ast}{\sim}}
\def\oeq{\stackrel{\rm o}{\sim}}
\def\new{\bigskip\noindent}
\def\note{$\circ\ $}
\def\bull{$\bullet$}
\def\IR{{\mathbb R}}
\def\IC{{\mathbb C}}
\def\IF{{\mathbb F}}
\def\IP{{\mathbb P}}
\def\la{{\langle}}
\def\ra{{\rangle}}
\def\cO{{\mathcal O}}
\def\i{{\mathbf i}}

\def\beq{\begin{equation}}
\def\eeq{\end{equation}}

\begin{document}


\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}

\begin{center}
{\Large\bf Numerical Linear Algebra} \\
{Nikolai Chernov}
\end{center}

\section{Fundamentals}

\subsection{Matrices}

The set of $m\times n$ matrices ($m$ rows, $n$ columns) with
entries in a field $\IF$ is denoted by $\IF^{m\times n}$ (we will
only consider two fields: $\IF=\IC$ and $\IF=\IR$). For any matrix
$A\in\IF^{m\times n}$, we denote its entries by
$$
      A =
   \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\
                   a_{21} & a_{22} & \cdots & a_{2n} \\
                   \vdots & \vdots & \ddots & \vdots \\
                   a_{m1} & a_{m2} & \cdots & a_{mn}
   \end{bmatrix}
    \in \IF^{m\times n}.
$$
The vector space $\IF^n$ consists of column vectors with $n$
components:
$$
     x=\begin{bmatrix} x_1\\ \vdots \\ x_n \end{bmatrix}
     \in \IF^n.
$$
The product $y=Ax$ is a vector in $\IF^m$:
$$
      y= Ax =
   \left[\!\!\!
   \begin{array}{c} \phantom{\vdots} \\ a_1 \\ \phantom{\vdots} \end{array}
   \!\!\!\right | \left . \!\!\!\!
   \begin{array}{c} \phantom{\vdots} \\ a_2 \\ \phantom{\vdots} \end{array}
   \!\!\!\right | \left . \!\!\!\!
   \begin{array}{c} \phantom{\vdots} \\ \cdots \\ \phantom{\vdots} \end{array}
   \!\!\!\right | \left . \!\!\!\!
   \begin{array}{c} \phantom{\vdots} \\ a_n \\ \phantom{\vdots} \end{array}
   \!\!\!\right ]
   \begin{bmatrix} x_1\\ \vdots \\ x_n \end{bmatrix}
   =x_1\!\!\begin{bmatrix}
   \phantom{\vdots} \\ a_1 \\ \phantom{\vdots}
   \!\!\end{bmatrix}
   +x_2\!\!\begin{bmatrix}
   \phantom{\vdots} \\ a_2 \\ \phantom{\vdots}
   \!\!\end{bmatrix}
   +\cdots+x_n\!\!\begin{bmatrix}
   \phantom{\vdots} \\ a_n \\ \phantom{\vdots}
   \!\!\end{bmatrix}
$$
where $a_1,\ldots,a_n$ denote the columns of the matrix $A$ (note that
$y$ is a linear combination of the columns of $A$). Every matrix
$A\in\IF^{m\times n}$ defines a linear transformation
$$
   A\colon \IF^n\to\IF^m\qquad{\rm by}\qquad
   x\mapsto Ax.
$$
The {\em range} of $A$ is a subspace of $\IF^m$
$$
   {\rm Range}\, A = \{Ax\colon x\in\IF^n\} \subset \IF^m
$$
spanned by the columns of $A$. The {\em rank} of $A$ is defined by
$$
   {\rm rank}\, A = {\rm dim}({\rm Range}\, A).
$$
The {\em kernel} (also called the {\em nullspace}) of $A$ is a
subspace in $\IF^n$
$$
   {\rm Ker}\, A = \{x\colon Ax=0\} \subset \IF^n.
$$
Note that
$$
       {\rm dim}({\rm Range}\, A)+
       {\rm dim}({\rm Ker}\, A) = n.
$$
\note The transformation $A$ is surjective iff ${\rm Range}\,A=\IF^m$.\\
\note The transformation $A$ is injective iff ${\rm Ker}\,A=\{0\}$.\\
\note If $A$ is bijective, then $m=n$ and we call $A$ an {\em
isomorphism}.

\subsection{Square matrices} \label{sqmat}

Every square matrix $A\in\IF^{n\times n}$ defines a linear
transformation $\IF^n\to\IF^n$, called an {\em operator} on $\IF^n$.
The {\em inverse} of a matrix $A$ is uniquely defined by
$$
    A^{-1}A=AA^{-1}=I \qquad (\text{identity matrix}).
$$
A matrix $A$ is said to be {\em invertible} ({\em nonsingular}) iff
$A^{-1}$ exists, otherwise the matrix is {\em noninvertible} ({\em
singular}). The following are equivalent:\\

(a) $A$ is invertible

(b) rank$\,A=n$

(c) Range$\,A=\IF^n$

(d) Ker$\, A=\{0\}$

(e) $0$ is not an eigenvalue of $A$

(f) det$\,A\neq 0$\\

\noindent Note that
$$
      (AB)^{-1} = B^{-1}A^{-1}
$$
A matrix is $A$ {\em upper triangular} if $a_{ij}=0$ for all
$i>j$. A matrix $A$ is {\em lower triangular} if $a_{ij}=0$ for
all $i<j$. A matrix $D$ is {\em diagonal} if $d_{ij}=0$ for all
$i\neq j$, and in that case we write $D=\,{\rm diag} \{d_{11},
\ldots, d_{nn}\}$. Note: if $A$ and $B$ are upper (lower)
triangular, then so are $AB$ and $A^{-1}$ and $B^{-1}$ (if the
inverses exist).

We say that $\lambda\in\IF$ is an {\em eigenvalue} for a matrix
$A\in\IF^{n\times n}$ with an {\em eigenvector} $x\in\IF^n$ if
$$
      Ax=\lambda x\qquad\text{and}\qquad x\neq 0.
$$
Eigenvalues are the roots of the characteristic polynomial
$$
    C_A(\lambda)=\,{\rm det}(\lambda I-A)=0.
$$
By the fundamental theorem of algebra, every polynomial of degree
$n$ with complex coefficients has exactly $n$ complex roots
(counting multiplicities), i.e.
$$
   C_A(\lambda) = \prod_{i=1}^n (\lambda - \lambda_i).
$$
A matrix $A\in\IF^{n\times n}$ is {\em diagonalizable} (over $\IF$) iff
\beq
         A = X\Lambda X^{-1},   \tag{ED}
\eeq
where $\Lambda =\,\rm{diag}\{\lambda_1,\ldots,\lambda_n\}$ is a
diagonal matrix and $X\in\IF^{n\times n}$. In this case
$\lambda_1,\ldots,\lambda_n$ are the eigenvalues of $A$ and the
columns $x_1,\ldots,x_n$ of the matrix $X$ are the corresponding
eigenvectors. To see this, rewrite the above equation in the form
$AX=X\Lambda$:
$$
      \left[\!\!\!
   \begin{array}{c} \phantom{\lambda_1} \\ \phantom{\ddots} \\ \phantom{\lambda} \end{array}
   \begin{array}{c} \phantom{\lambda} \\ A \\ \phantom{\lambda} \end{array}
   \begin{array}{c} \phantom{\lambda} \\ \phantom{\ddots} \\ \phantom{\lambda_n} \end{array}
   \!\!\!\right ]
      \left[\!\!\!
   \begin{array}{c} \phantom{\vdots} \\ x_1 \\ \phantom{\vdots} \end{array}
   \!\!\!\right | \left . \!\!\!\!
   \begin{array}{c} \phantom{\vdots} \\ x_2 \\ \phantom{\vdots} \end{array}
   \!\!\!\right | \left . \!\!\!\!
   \begin{array}{c} \phantom{\vdots} \\ \cdots \\ \phantom{\vdots} \end{array}
   \!\!\!\right | \left . \!\!\!\!
   \begin{array}{c} \phantom{\vdots} \\ x_n \\ \phantom{\vdots} \end{array}
   \!\!\!\right ]
     =
      \left[\!\!\!
   \begin{array}{c} \phantom{\vdots} \\ x_1 \\ \phantom{\vdots} \end{array}
   \!\!\!\right | \left . \!\!\!\!
   \begin{array}{c} \phantom{\vdots} \\ x_2 \\ \phantom{\vdots} \end{array}
   \!\!\!\right | \left . \!\!\!\!
   \begin{array}{c} \phantom{\vdots} \\ \cdots \\ \phantom{\vdots} \end{array}
   \!\!\!\right | \left . \!\!\!\!
   \begin{array}{c} \phantom{\vdots} \\ x_n \\ \phantom{\vdots} \end{array}
   \!\!\!\right ]
      \left[\!\!\!
   \begin{array}{c} \lambda_1 \\ \phantom{\ddots} \\ \phantom{\lambda} \end{array}
   \begin{array}{c} \phantom{\lambda} \\ \ddots \\ \phantom{\lambda} \end{array}
   \begin{array}{c} \phantom{\lambda} \\ \phantom{\ddots} \\ \lambda_n \end{array}
   \!\!\!\right ].
$$
The formula (ED) is referred to as an \emph{eigendecomposition} of $A$.

Note that a real matrix $A\in\IR^{n\times n}$ may be diagonalizable over
$\IC$ but not over $\IR$. This happens, for example, when $A$ has distinct
complex eigenvalues.

The {\em trace} of a matrix $A\in\IF^{n\times n}$ is
$$
    {\rm tr}\,A = \sum_{i=1}^n a_{ii}
$$
Trace has the following properties:

(a) tr$\,AB=\,$tr$\,BA$;

(b) if $A=X^{-1}BX$, then tr$\,A=\,$tr$\,B$;

(c) tr$\,A=\lambda_1+\cdots+\lambda_n$ (the sum of all complex
eigenvalues).

\subsection{Transposed and adjoint matrices}

For any matrix $A=(a_{ij})\in\IF^{m\times n}$ we denote by
$A^T=(a_{ji})\in\IF^{n\times m}$ the {\em transpose} of $A$. Note
that
$$
      (AB)^T = B^TA^T, \qquad  \bigl(A^T\bigr)^T=A.
$$
If $A$ is a square matrix, then
$$
   {\rm det}\,A^T=\,{\rm det}\,A,
   \qquad \bigl(A^T\bigr)^{-1}=\bigl(A^{-1}\bigr)^T.
$$
For any matrix $A=(a_{ij})\in\IC^{m\times n}$ we denote by
$A^{\ast}=(\bar{a}_{ji})\in\IF^{n\times m}$ the {\em adjoint} of
$A$. Note that
$$
      (AB)^{\ast} = B^{\ast}A^{\ast},
      \qquad \bigl(A^{\ast}\bigr)^{\ast}=A.
$$
If $A$ is a square matrix, then
$$
   {\rm det}\,A^{\ast}=\,\overline{{\rm det}\,A},
   \qquad \bigl(A^{\ast}\bigr)^{-1}=\bigl(A^{-1}\bigr)^{\ast}.
$$
For $A\in\IR^{m\times n}$, we have $A^{\ast}=A^T$.

\subsection{Norms}
\label{8.1}

A {\em norm} on a vector space $V$ over $\IC$ is a real valued
function $\|\cdot \|$ on $V$ satisfying three axioms:
\begin{enumerate}
\item $\|v\|\geq 0$ for all $v\in V$ and $\|v\|=0$ if and only if
$v=0$.\item $\|cv\|=|c|\, \|v\|$ for all $c\in\IC$ and $v\in
V$.\item $\|u+v\|\leq \|u\|+\|v\|$ for all $u,v\in V$ ({\em
triangle inequality}).\end{enumerate} If $V$ is a vector space
over $\IR$, then an inner product is a function $V\to\IR$
satisfying the same axioms, except $c\in\IR$.

Some common norms on $\IC^n$ and $\IR^n$ are:
\begin{align}
   \|x\|_1 &=\sum_{i=1}^n |x_i|  \tag{1-norm}\\
   \|x\|_2 &=\biggl(\sum_{i=1}^n |x_i|^2\biggr)^{1/2}\tag{2-norm}\\
   \|x\|_p &=\biggl(\sum_{i=1}^n |x_i|^p\biggr)^{1/p} \tag{$p$-norm, $p\geq 1$}\\
   \|x\|_{\infty} &=\max_{1\leq i\leq n} |x_i| \tag{$\infty$-norm}
\end{align}
The 2-norm in $\IR^n$ corresponds to the Euclidean distance.

Some norms on $C[a,b]$, the space of continuous functions on
$[a,b]$:
\begin{align}
       \|f\|_1 &=\int_a^b|f(x)|\, dx   \tag{1-norm}\\
       \|f\|_2 &=\biggl(\int_a^b|f(x)|^2\, dx\biggr)^{1/2} \tag{2-norm}\\
       \|f\|_{\infty} &=\max_{a\leq x\leq b} |f(x)| \tag{$\infty$-norm}
\end{align}
We define
\beq
   {\mathbb S}_1=\{v\in V\colon\, \|v\|=1\}
    \tag{\bf Unit sphere}
\eeq
The vectors $v\in{\mathbb S}_1$ are called {\em unit vectors}. For any
$v\neq 0$, the vector $u=v/\|v\|$ is a unit vector.

The space $\IC^{m\times n}$ of matrices is isomorphic to
$\IC^{mn}$, hence we can define norms on it in a similar way. In
particular, an analogue of the 2-norm
$$
   \|A\|_{\rm F}=\biggl(\sum_i\sum_j |a_{ij}|^2\biggr)^{1/2}
$$
is known as {\em Frobenius norm}. Note that
$$
   \|A\|_F^2=\,{\rm tr}\,(A^{\ast}A)=\,{\rm tr}\,(AA^{\ast}).
$$
However, most important are matrix norms induced by vector norms:

\subsection{Induced matrix norms} \label{8.8}

Let $A\colon\IF^n\to\IF^m$, and let the spaces $\IF^n$ and $\IF^m$
be equipped with a certain norm $\|\cdot\|$. Then
$$
    \|A\|\colon
    =\sup_{\|x\|=1} \|Ax\|
    =\sup_{x\neq 0}\frac{\|Ax\|}{\|x\|}
$$
defines an {\em induced norm} (also called an {\em operator norm})
of $A$. Respectively, we obtain $\|A\|_2$, $\|A\|_1$, and
$\|A\|_{\infty}$, if the spaces $\IF^n$ and $\IF^m$ are equipped
with $\|\cdot \|_2$, $\|\cdot \|_1$, and $\|\cdot \|_{\infty}$.

The supremum here is always attained and can be replaced by
maximum. This follows from the compactness of ${\mathbb S}_1$ and
the continuity of $\|\cdot\|$. For the 2-norm, this can be also
proved by an algebraic argument, see \ref{12.15} below.

Remark: there are norms on $\IC^{n\times n}$ that are not induced
by any norm on $\IC^n$, for example $\|A\|=\max_{i,j}|a_{ij}|$
(see homework assignment).

We have simple rules for computing $\|A\|_1$ and $\|A\|_{\infty}$:
\begin{align}
  \|A\|_1&=\max_{1\leq j\leq n}\sum_{i=1}^m|a_{ij}|
     \tag{maximum column sum} \\
  \|A\|_{\infty}&=\max_{1\leq i\leq m}\sum_{j=1}^n|a_{ij}|
   \tag{maximum row sum}
\end{align}
but there is no explicit formulas for $\|A\|_2$ in terms of the
$a_{ij}$'s.

Any induced matrix norm satisfies
$$
   \|Ax\|\leq \|A\|\, \|x\| \quad {\rm and}\quad
   \|AB\|\leq \|A\|\, \|B\|.
$$

\subsection{Inner products}\label{8.11}

Let $V$ be a vector space over $\IC$. An {\em inner product} on $V$ is
a function on $V\times V\to\IC$, denoted by $\la\cdot,\cdot\ra$,
satisfying four axioms:
\begin{enumerate}\item $\la u,v\ra=\overline{\la v,u\ra}$ for
all $u,v\in V$. \item $\la cu,v\ra=c\la u,v\ra$ for all $c\in\IC$ and
$u,v\in V$. \item $\la u+v,w\ra=\la u,w\ra+\la v,w\ra$ for all
$u,v,w\in V$. \item $\la u,u\ra\geq 0$ for all $u\in V$, and $\la
u,u\ra=0$ iff $u=0$.
\end{enumerate}
Note that $\la u,v+w\ra=\la u,v\ra+\la u,w\ra$ and $\la u,cv\ra =
\bar{c}\la u,v\ra$.

If $V$ is a vector space over $\IR$, then an inner product is a
function $V\times V\to\IR$ satisfying the same axioms (except
$c\in\IR$, and there is no need to take a conjugate). A standard
inner product in $\IC^n$ is
$$
   \la x,y\ra=\sum_{i=1}^n x_i\bar{y}_i
   =x^T\bar{y}=y^{\ast}x.
$$
A standard inner product in $\IR^n$ is
$$
   \la x,y\ra=\sum_{i=1}^n x_i y_i
   =x^Ty=y^Tx.
$$
A standard inner product in $C([a,b])$, the space of continuous
functions, is:
$$
   \la f,g\ra = \int_a^b f(x)\,\bar{g}(x)\,dx
$$

\subsection{Cauchy-Schwarz inequality}
\label{8.16}

Let $V$ be an inner product space. Then $$ |\la u,v\ra |\leq \la
u,u\ra^{1/2}\la v,v\ra^{1/2}$$ for all $u,v\in V$.

The equality holds if and only if $u$ and $v$ are linearly
dependent.


\proof We do it in the complex case. Assume that $v\neq 0$. Consider
the function
\begin{align*}
   f(z)&= \la u-zv,u-zv\ra \\
   &= \la u,u\ra
   - z\la v,u\ra
   - \bar{z}\la u,v\ra
   + |z|^2\la v,v\ra
\end{align*}
of a complex variable $z$. Let $z=re^{i\theta}$ and $\la
u,v\ra=se^{i\varphi}$ be the polar forms of the numbers $z$ and $\la
u,v\ra$. Set $\theta=\varphi$ and assume that $r$ varies from $-\infty$
to $\infty$, then
$$
   0\leq f(z)= \la u,u\ra
   - 2sr + r^2\la v,v\ra
$$
Since this holds for all $r\in\IR$ (also for $r<0$, because the
coefficients are all nonnegative), the discriminant has to be $\leq 0$,
i.e. $s^2- \la u,u\ra\la v,v\ra\leq 0$. This completes the proof in the
complex case. It the real case it goes even easier, just assume
$z\in\IR$. The equality case in the theorem corresponds to the zero
discriminant, hence the above polynomial assumes a zero value, and
hence $u=zv$ for some $z\in\IC$. (We left out the case $v=0$, do it
yourself as an exercise.) \qed

\subsection{Induced norms} \label{8.17}

If $V$ is an inner product vector space, then
\beq
   \|v\|=\la v,v\ra^{1/2}
    \tag{induced norm}
\eeq
defines a norm on $V$ (to prove the triangle inequality, we need
Theorem~\ref{8.16}). In vector spaces over $\IR$, a norm is induced by
an inner product if and only if the function
\beq
  \la u,v\ra \colon = \tfrac 14
   \left( \|u+v\|^2-\|u-v\|^2\right)
      \tag{\bf Polarization identity}
\eeq
satisfies the axioms of inner products. (A similar but more
complicated polarization identity holds in vector spaces over
$\IC$.)

\subsection{Orthogonal vectors} \label{9.1}

Two vectors $u,v\in V$ are said to be {\em orthogonal} if $\la u,v\ra =0$.
In this case
\beq
      \|u+v\|^2=\|u\|^2+\|v\|^2
      \tag{\bf Pythagorean Theorem}
\eeq
Inductively, if $u_1,\ldots,u_k$ are mutually orthogonal, then
$$
  \|u_1+\cdots +u_k\|^2=\|u_1\|^2+\cdots +\|u_k\|^2
$$
If nonzero vectors $u_1,\ldots,u_k$ are mutually orthogonal, then
they are linearly independent.

\subsection{Orthonormal basis (ONB)}
\label{9.5}

The set $\{u_1,\ldots,u_n\}$ is an {\em orthonormal basis} (ONB)
in $V$, if it is a basis and all the vectors $u_i$ are mutually
orthogonal and have unit length (i.e., $\|u_i\|=1$ for all $i$).
Note: $\la u_i,u_j\ra = \delta_{ij}$ (this is the Kronecker delta
symbol: $\delta_{ij}=1$ if $i=j$ and $\delta_{ij}=0$ if $i\neq
j$).

If $\{u_1,\ldots,u_n\}$ is an ONB, then for any vector $v$
\beq
   v=\sum_{i=1}^n\la v,u_i\ra\, u_i
   \tag{\bf Fourier expansion}
\eeq
In other words, the numbers $\la v,u_i\ra$ ({\em Fourier coefficients})
are the coordinates of the vector $v$ in this basis.

\subsection{Orthogonal projection} \label{9.7}

Let $u,v\in V$, and $u\neq 0$. The {\em orthogonal projection} of
$v$ onto $u$ is
$$
    {\rm Pr}_uv=\frac{\la v,u\ra}{\|u\|^2}\,u
$$

\new Note that the vector $w\colon = v-{\rm Pr}_uv$ is orthogonal
to $u$. Therefore, $v$ is the sum of two vectors, Pr$_uv$
(parallel to $u$), and $w$ (orthogonal to $u$), see the diagram
below.

\begin{figure}[h]\caption{Orthogonal projection of $u$ onto $v$}\label{fig:1}
\setlength{\unitlength}{1.00mm} \special{em:linewidth 0.4pt}
\linethickness{0.4pt}
%\begin{picture}(125.00,87.50)
\begin{picture}(125.00,47.50)(25,0)
\put(45.00,10.00){\vector(2,1){55.00}}
\put(46.00,10.00){\vector(1,0){79.00}}
\put(100.00,37.00){\vector(0,-1){27.00}}
\put(61.00,14.00){\makebox(0,0)[cc]{$\theta$}}
\put(74.00,30.00){\makebox(0,0)[cc]{$v$}}
\put(123.00,3.00){\makebox(0,0)[cc]{$u$}}
\put(75.00,3.00){\makebox(0,0)[cc]{$\|v\|\cos\theta\cdot
(u/\|u\|)$}}
\end{picture}
\end{figure}

In the real case, for any nonzero vectors $u,v\in V$ let
\beq
   \cos\theta=\frac{\la v,u\ra}{\|u\|\, \|v\|}
   \tag{\bf Angle}
\eeq
By \ref{8.16}, we have $\cos\theta\in [-1,1]$. Hence, there is a unique
angle $\theta\in [0,\pi]$ with this value of cosine. It is called the
{\em angle between} $u$ and $v$.

Note that $\cos\theta=0$ if and only if $u$ and $v$ are
orthogonal. Also, $\cos\theta=\pm 1$ if and only if $u,v$ are
proportional, $v=cu$, then the sign of $c$ coincides with the sign
of $\cos\theta$.

\subsection{Gram-Schmidt orthogonalization} \label{9.11}

Let $u_1,\ldots,u_k$ be unit mutually orthogonal vectors. For
$v\in V$, set
$$
    w=v-\sum_{i=1}^k\la v,u_i\ra u_i.
$$
Then the vectors $u_1,\ldots,u_k,w$ are mutually orthogonal, and
$$
   {\rm span}\{u_1,\ldots,u_k,v\}=\,
   {\rm span}\{u_1,\ldots,u_k,w\}.
$$
In particular, $w=0$ if and only if $v\in\,$span$\{u_1,
\ldots,u_k\}$.

Next, let $\{v_1,\ldots,v_n\}$ be a basis in $V$. Define
$$
            w_1=v_1 \qquad \text{and} \qquad
            u_1=w_1/\|w_1\|,
$$
and then inductively, for $k\geq 1$,
\beq
   w_{k}=v_{k}-\sum_{i=1}^{k-1}\la v_{k},u_i\ra\, u_i,
   \quad \text{and} \quad
   u_k=w_k/\|w_k\|
    \tag{GS}
\eeq
This gives an orthogonal basis $\{u_1,\ldots,u_n\}$, which
`agrees' with the basis $\{v_1,\ldots,v_n\}$ in the following
sense:
$$
   {\rm span}\{v_1,\ldots,v_k\}=\,{\rm span}\{u_1,\ldots,u_k\}
   \quad \forall\ 1\leq k\leq n
$$
As a corollary, we conclude that every finite dimensional vector
space with an inner product has an ONB. Furthermore, every set of
orthonormal vectors $\{u_1,\ldots,u_k\}$ can be extended to an
ONB.

\subsection{Legendre polynomials} \label{9.12}

Let $V=\IP_n(\IR)$, the space of real polynomials of degree $\leq n$,
with the inner product given by $\la f,g\ra=\int_0^1 f(x)g(x)\, dx$.
Applying Gram-Schmidt orthogonalization to the basis
$\{1,x,\ldots,x^n\}$ gives the first $n+1$ of the so called {\em
Legendre polynomials} (see also homework assignment).

\subsection{Orthogonal complement} \label{9.14}

Let $S\subset V$ be a subset (not necessarily a subspace). Then
$$
 S^{\perp}\colon =
 \{v\in V\colon\, \la v,w\ra =0\ \
 {\rm for}\ {\rm all}\ \ w\in S\}
$$
is called the {\em orthogonal complement} to $S$. Note that
$S^{\perp}$ is a subspace of $V$. Also, $(W^\perp)^\perp\subset
W$, and if $V$ is finite dimensional, then $(W^\perp)^\perp = W$
(see homework assignment).

If $W$ is a finite dimensional subspace of $V$, then $V=W\oplus
W^{\perp}$. To prove this, choose an ONB $\{u_1,\ldots,u_n\}$ of
$W$ (by \ref{9.11}), and then for any $v\in V$
$$
  v-\sum_{i=1}^n\la v,u_i\ra\, u_i \in W^\perp.
$$
Note: the condition dim$\,W<\infty$ is essential. Let $V=C[a,b]$
with the standard inner product and $W\subset V$ be the set of
real polynomials restricted to the interval $[a,b]$. Then
$W^{\perp}=\{0\}$, and at the same time $V\neq W$.

Note: if $\{u_1,\ldots,u_n\}$ is an orthonormal subset of $V$,
then
\beq
  \|v\|^2\geq\sum_{i=1}^n|\la v,u_i\ra|^2
   \tag{\bf Bessel's inequality}
\eeq
If $\{u_1,\ldots,u_n\}$ is an ONB, then
$$
   \|v\|^2 = \sum_{i=1}^n |\la v,u_i\ra|^2
$$
and, more generally,
\beq
   \la v,w\ra = \sum_{i=1}^n
   \la v,u_i\ra \overline{\la w,u_i\ra}
    \tag{\bf Parceval's identity}
\eeq
This follows from the Fourier expansion, see \ref{9.5}.

\newpage

\section{Unitary matrices}
\label{10}

\subsection{Isometry} \label{9.20}

Let $V$ and $W$ be two inner product spaces (both real or both
complex). An isomorphism $T\colon V\to W$ is called an {\em isometry}
if it preserves the inner product, i.e. $$\la Tv,Tw\ra=\la v,w\ra$$ for
all $v,w\in V$. In this case $V$ and $W$ are said to be {\em
isometric}.

\new Polarization identity \ref{8.17} implies that
$T$ preserves inner product iff $T$ preserves the induced norm,
i.e. $\|Tv\|=\|v\|$ for all $v\in V$ (the latter need only be
checked for unit vectors).

\subsection{Characterization of isometries} \label{9.21}

Let dim$\, V<\infty$. A linear transformation $T\colon V\to W$ is an
isometry if and only if there exists an ONB $\{u_1,\ldots,u_n\}$ in $V$
such that $\{Tu_1,\ldots,Tu_n\}$ is an ONB in $W$.

\proof Use Parceval's identity \ref{9.14}. \qed \medskip

As a result, finite dimensional inner product spaces $V$ and $W$ (over
the same field) are isometric if and only if $\dim V=\dim W$. Thus, all
finite dimensional real inner product spaces can be identified with
$\IR^n$, $n\geq 1$, and complex -- with $\IC^n$, $n\geq 1$, with
standard inner products.

\subsection{Unitary matrices} \label{10.2}

A matrix $Q\in\IC^{n\times n}$ is said to be {\em unitary} if
$Q^{\ast}Q=I$ (in the real case, $Q$ is said to be {\em
orthogonal}).

\new Note: If $Q$ is unitary, then $QQ^{\ast}=I$, and so
$Q^{\ast}$ is a unitary matrix as well. Furthermore, $\bar{Q} Q^T
=I$ and $Q^T\bar{Q}=I$, hence $\bar{Q}$ and $Q^T$ are also unitary
matrices.

\subsection{Theorem} \label{10.3}

\begin{itemize}
\item[\textcopyright ] The linear
transformation of $\IC^n$ defined by a matrix $Q\in\IC^{n\times n}$
preserves the standard inner product iff $Q$ is
unitary.\item[\textregistered] The linear transformation of $\IR^n$
defined by a matrix $Q\in\IR^{n\times n}$ preserves the standard inner
product iff $Q$ is orthogonal.
\end{itemize}


\proof $\la Qx,Qy\ra =(Qy)^{\ast}Qx=y^\ast Q^\ast Qx=y^\ast x=\la
x,y\ra$. \qed

\subsection{Remarks}

Unitary $n\times n$ matrices make a group, denoted by $U(n)$.
Orthogonal $n\times n$ matrices make a group, denoted by $O(n)$.

For any matrix $A\in\IC^{m\times n}$ and a unitary matrix
$U\in\IC^{m\times m}$ and a unitary matrix $V\in\IC^{n\times n}$ we
have
$$
   \|A\|_2 = \|UA\|_2 = \|AV\|_2
$$
(see homework assignment).

\subsection{Examples of orthogonal matrices} \label{9.23}

$Q=\begin{bmatrix}0 & 1 \\ 1 & 0\end{bmatrix}$ defines a reflection
across the line $y=x$;\\

$Q=\left[\begin{array}{rr}
      \cos\theta & -\sin\theta \\
      \sin\theta & \cos\theta
   \end{array} \right]$
is a counterclockwise rotation by angle $\theta$.


\subsection{Characterization of unitary matrices}\label{10.6}

A matrix $Q\in \IC^{n\times n}$ is unitary iff its columns make an ONB.
The same holds if one replaces columns by rows.


\subsection{Theorem} \label{10.7}

If $Q$ is unitary, then $|{\rm det}\, Q|=1$.


\proof $1=\,$det$\, I=\,$det$\, Q^{\ast}Q=\,$det$\,
Q^\ast\cdot\,$det$\, Q= |\,$det$\, Q|^2$. \qed

\new Orthogonal matrices have determinant 1 or $-1$.
Orthogonal $n\times n$ matrices with determinant 1 make a subgroup of
$O(n)$, denoted by $SO(n)$.

\subsection{Theorem} \label{10.9}

If $\lambda$ is an eigenvalue of a unitary matrix, then
$|\lambda|=1$.

\proof If $Qx=\lambda x$ for some $x\neq 0$, then $\la x,x\ra=\la
Qx,Qx\ra = \la \lambda x, \lambda x\ra =|\lambda|^2\la x,x\ra$, so that
$|\lambda|^2=1$.

\subsection{Theorem} \label{10.10}

Let $T\colon V\to V$ be an isometry and dim$\,V<\infty$. If a subspace
$W\subset V$ is invariant under $T$, i.e.\ $TW\subset W$, then so is
its orthogonal complement $W^{\perp}$, i.e.\ $TW^{\perp}\subset
W^{\perp}$.

\proof Since Ker$\,T=\{0\}$, then $TW=W$, hence $\forall w\in W\
\exists w'\in W\colon\ Tw'=w$. Now, if $v\in W^{\perp}$, then $\la
w,Tv\ra=\la Tw',Tv\ra=\la w',v\ra=0$, hence $Tv\in W^{\perp}$.

\subsection{Corollary} \label{10.11}

For any isometry $T$ of a finite dimensional complex space $V$
there is an ONB of $V$ consisting of eigenvectors of $T$.

\proof Use induction on the dimension of the space and \ref{10.10}.

\subsection{Lemma} \label{10.11a}

Every operator $T\colon V\to V$ on a finite dimensional real space $V$
has either a one-dimensional or two-dimensional invariant subspace
$W\subset V$.

\proof Let $T$ be represented by a matrix $A\in\IR^{n\times n}$ in some
basis. If $A$ has a real eigenvalue, then $Ax=\lambda x$ and we get a
one-dimensional invariant subset span$\{x\}$. If $A$ has no real
eigenvalues, then the matrix $A$, {\em considered as a complex matrix},
has a complex eigenvalue $\lambda=a+\i b$, with $a,b\in\IR$ and $\i =
\sqrt{-1}$, and a complex eigenvector $x+\i y$, with $x,y\in\IR^n$. The
equation
$$
    A(x+\i y)=(a+\i b)(x+\i y)
$$
can be written as
\begin{align*}
         Ax &= a x-b y \\
         Ay &= b x+a y
\end{align*}
Thus, the space span$\{x,y\}$ is invariant. \qed

\new Note: one can easily verify that $x$ and $y$ are linearly
dependent iff $b=0$, i.e.\ $\lambda\in\IR$.


\subsection{Theorem} \label{10.12}

Let $T\colon V\to V$ be an orthogonal operator on a finite dimensional
real space $V$. Then $V=V_1\oplus\cdots\oplus V_m$, where $V_i$ are
mutually orthogonal subspaces, each $V_i$ is $T$-invariant, and either
dim$\,V_i=1$ or dim$\,V_i=2$.

\proof Use induction on the dimension of $V$ and apply \ref{10.10} and
\ref{10.11a}.

\new Note: the restriction of the operator $T$ to each
of the two-dimensional invariant subspaces is simply a rotation by some
angle (as in Example~\ref{9.23}), see homework assignment.

\subsection{Unitary and orthogonal equivalence}
\label{11.17}


\begin{itemize}\item[\textcopyright]
Two complex matrices $A,B\in \IC^{n\times n}$ are said to be {\em
unitary equivalent}
%($A\ueq B$)
if $B=P^{-1}AP$ for some unitary matrix
$P$. This can be also written as $B=P^{\ast}AP$.
\item[\textregistered] Two real matrices $A,B\in \IR^{n\times n}$ are
said to be {\em orthogonally equivalent}
%($A\oeq B$)
if $B=P^{-1}AP$ for some orthogonal matrix $P$. This can be also
written as $B=P^TAP$. \end{itemize}

By Theorem~\ref{10.11}, any unitary matrix $Q$ is unitary
equivalent to a diagonal matrix $D=\,{\rm
diag}\{d_1,\ldots,d_n\}$, whose diagonal entries belong to the
unit circle, i.e.\ $|d_i|=1$ for $1\leq i\leq n$.\medskip

By Theorem~\ref{10.12}, any orthogonal matrix $Q$ is orthogonally
equivalent to a block-diagonal matrix
$$
     B=\left[\begin{array}{cccc}
     R_{11} & 0 & \cdots & 0 \\
     0 & R_{22} & \cdots & 0 \\
     \vdots & \vdots & \ddots & \vdots \\
     0 & 0 & \cdots & R_{mm}
     \end{array}\right]
$$
where $R_{ii}$ are $1\times 1$ and $2\times 2$ diagonal blocks.
Furthermore, all $1\times 1$ blocks are either $R_{ii}=+1$ or
$R_{ii}=-1$, and the $2\times 2$ blocks are rotation matrices
$$
  R_{ii}=\left[\begin{array}{rr}\cos\theta & -\sin\theta \\
  \sin\theta & \cos\theta\end{array} \right]
$$

\newpage

\section{Adjoint and self-adjoint matrices}
\label{11}

In this chapter, $V$, $W$ denote finite dimensional inner product
spaces (unless stated otherwise).

\subsection{Riesz representation theorem} \label{11.1}

Let $f\in V^{\ast}$, i.e. $f$ is a linear functional on $V$. Then there
is a unique vector $u\in V$ such that
$$
   f(v)=\la v,u\ra\ \ \ \ \ \ \forall v\in V
$$

\proof Let $B=\{u_1,\ldots,u_n\}$ be an ONB in $V$. Then for any
$v=\sum c_iu_i$ we have $f(v)=\sum c_if(u_i)$ by linearity. Also, for
any $u=\sum d_iu_i$ we have $\la v,u\ra = \sum c_i\bar{d}_i$. Hence,
the vector $u=\sum \overline{f(u_i)}\,u_i$ will suffice. The uniqueness
of $u$ is obvious.

\subsection{Corollary} \label{11.2}

The identity $f\leftrightarrow u$ established in the previous
theorem is ``quasi-linear'' in the following sense:
$f_1+f_2\leftrightarrow u_1+u_2$ and $cf\leftrightarrow
\bar{c}u$. In the real case, it is perfectly linear, though, and hence
it is an isomorphism between $V^{\ast}$ and $V$.

\subsection{Remark}\label{11.2a}

If dim$\, V=\infty$, then Theorem~\ref{11.1} fails. Consider, for example,
$V=C[0,1]$ (real functions) with the inner product $\la
F,G\ra=\int_0^1F(x)G(x)\, dx$. Pick a point $t\in [0,1]$. Let $f\in
V^{\ast}$ be a linear functional defined by $f(F)=F(t)$. It does not
correspond to any $G\in V$ so that $f(F)=\la F,G\ra$. In fact, the lack of
such functions $G$ has led mathematicians to the concept of {\em
delta-functions}: a delta-function $\delta_t(x)$ is ``defined'' by three
requirements: $\delta_t(x)\equiv 0$ for all $x\neq t$,
$\delta_t(t)=\infty$ and $\int_0^1F(x)\delta_t(x)\, dx=F(t)$ for every
$F\in C[0,1]$. Clearly, no regular function satisfies these requirements,
so ``delta function'' is a somewhat abstract concept.

\subsection{Adjoint operator}
\label{11.4}

Let $T\colon V\to W$ be a linear transformation. Then there is a unique
linear transformation $T^{\ast}\colon W\to V$ such that $\forall v\in
V$ and $\forall w\in W$
$$
   \la Tv,w\ra = \la v,T^{\ast}w\ra.
$$
$T^\ast$ is called the {\em adjoint} of $T$.

\proof Let $w\in W$. Then $f(v)\colon =\la Tv,w\ra$ defines a linear
functional $f\in V^{\ast}$. By the Riesz representation theorem, there
is a unique $v'\in V$ such that $f(v)=\la v,v'\ra$. Then we define
$T^{\ast}$ by setting $T^{\ast}w=v'$. The linearity of $T^{\ast}$ is a
routine check. Note that in the complex case the conjugating bar
appears twice and thus cancels. The uniqueness of $T^{\ast}$ is
obvious. \qed

\subsection{Adjoint matrix} \label{11.6a}

Let $A\in \IC^{m\times n}$ be an $m\times n$ matrix that defines a
linear transformation $T_A\colon\IC^n\to\IC^m$ (with the standard inner
product). Then its adjoint $T_A^{\ast}\colon\IC^m\to\IC^n$ is defined
by the adjoint matrix $A^{\ast}$, i.e.\ $T_A^{\ast}=T_{A^{\ast}}$.

\proof $\la Ax,y\ra = y^\ast Ax= \bigl(A^\ast y\bigr)^\ast x=\la
x,A^{\ast}y\ra$.

\subsection{Theorem} \label{11.7}

Let $T\colon V\to W$ be a linear operator. Then
$$
   {\rm Ker}\,T^{\ast}=\, ({\rm Range}\, T)^{\perp}
$$
In particular, $W=\, {\rm Range}\, T\oplus\, {\rm Ker}\, T^{\ast}$.
(Proof is a routine check.)

% \proof It is a routine check. \qed

% In particular, for any $A\in\IC^{m\times n}$ $$
%    {\rm Ker}\,A^{\ast}=\, ({\rm Range}\, A)^{\perp} $$

\subsection{Selfadjoint operators and matrices} \label{11.8}

A linear operator $T\colon V\to V$ is said to be {\em selfadjoint} if
$T^{\ast}=T$.

A square matrix $A$ is said to be {\em selfadjoint} if $A^{\ast}=A$. In
the real case, this is equivalent to $A^T=A$, i.e.\ $A$ is a symmetric
matrix. In the complex case, selfadjoint matrices are often called {\em
Hermitian matrices}.

\subsection{Theorem} \label{11.9}

Let $A\in\IC^{n\times n}$ be a selfadjoint matrix. Then
\begin{itemize}\item[(a)] All eigenvalues of $A$ are real. \item[(b)]
If $\lambda_1\neq\lambda_2$ are two distinct eigenvalues of $A$ with
respective eigenvectors $x_1,x_2$, then $x_1$ and $x_2$ are
orthogonal.\end{itemize}

\proof If $Ax=\lambda x$, then $\lambda\la x,x\ra= \la Ax,x\ra=\la
x,Ax\ra=\bar{\lambda}\la x,x\ra$, hence $\lambda=\bar{\lambda}$ (since
$x\neq 0$). To prove (b), we have $\lambda_1\la x_1,x_2\ra=\la
Ax_1,x_2\ra=\la x_1,Ax_2\ra=\overline{\lambda_2}\la x_1,x_2\ra$, and in
addition $\overline{\lambda_2}=\lambda_2\neq\lambda_1$, hence $\la
x_1,x_2\ra=0$. \qed

\new{\em Remark}. Theorem~\ref{11.9} holds in the real
case, too. It implies that all the roots of the characteristic
polynomial $C_A(x)$ of a symmetric matrix $A$ are real, i.e.\
$C_A(x)=\prod_i (x-\lambda_i)$ where all $\lambda_i$'s are real
numbers. In particular, any selfadjoint operator, real or complex, has
at least one eigenvalue.

\subsection{Theorem} \label{11.11}

Let $T$ be a selfadjoint operator and a subspace $W$ be
$T$-invariant, i.e. $TW\subset W$. Then $W^{\perp}$ is also
$T$-invariant, i.e. $TW^{\perp}\subset W^{\perp}$.

\proof If $v\in W^{\perp}$, then for any $w\in W$ we have $\la Tv,w\ra
=\la v,Tw\ra=0$, so $Tv\in W^{\perp}$.

\subsection{Spectral theorem} \label{11.20}

Let $T\colon V\to V$ be a selfadjoint operator. Then there is an ONB
consisting of eigenvectors of $T$.

\proof Use induction on $n=\,$dim$\, V$ and apply \ref{11.9} and
\ref{11.11}.

\subsection{Corollary} \label{11.21}

\begin{itemize}\item[\textcopyright]
Any complex Hermitian matrix is unitary equivalent to a diagonal matrix
(with real diagonal entries). \item[\textregistered] Any real symmetric
matrix is orthogonally equivalent to a diagonal matrix. \end{itemize}

\subsection{Theorem (converse)} \label{11.22}

\begin{itemize}\item[\textcopyright]
If a complex matrix $A$ is unitary equivalent to a diagonal matrix with
real diagonal entries, then $A$ is Hermitian. \item[\textregistered] If
a real matrix $A$ is orthogonally equivalent to a diagonal matrix, then
$A$ is symmetric. \end{itemize}

\proof If $A=P^{-1}DP$ with $P^{-1}=P^{\ast}$, then
$A^{\ast}=P^{\ast}D^{\ast}(P^{-1})^{\ast}=A$, because $D^{\ast}=D$ (as
a real diagonal matrix).

\subsection{Remark} \label{11.22a}

Let $A=QDQ^{\ast}$, where $Q$ is a unitary matrix and $D$ is a diagonal
matrix. Denote by $q_i$ the $i$th column of $Q$ and by $d_i$ the $i$th
diagonal entry of $D$. Then
$$
         Aq_i = d_iq_i,\quad\quad 1\leq i\leq n
$$
i.e.\ the columns of $Q$ are the eigenvectors of $A$, whose
eigenvalues are the corresponding diagonal elements of $D$ (see
\ref{sqmat}).

\subsection{Theorem} \label{11.24}

If an operator $T$ is selfadjoint and invertible, then so is $T^{-1}$.
If a matrix $A$ is selfadjoint and nonsingular, then so is $A^{-1}$.

\proof By the Spectral Theorem~\ref{11.20}, there is an ONB $B$
consisting of eigenvectors of $T$. Now $T^{-1}$ has the same
eigenvectors, and its eigenvalues are the reciprocals of those of $T$,
hence they are real, too. Therefore, $T^{-1}$ is selfadjoint.

\subsection{Projections} \label{11.12}

Let $V$ be a vector space. A linear map $P\colon V\to V$ is called
a {\em projection} (or a {\em projector}) if $P^2=P$.

The following properties of projections can be easily verified
directly: Ker$\, P=\,{\rm Range}\, (I-P)$ and $V=\,{\rm Ker}\,
P\oplus\,{\rm Range}\, P$. For any vector $v\in V$ there is a unique
decomposition $v=v_1+v_2$ with $v_1\in\,{\rm Ker}\, P$ and
$v_2\in\,{\rm Range}\, P$ and $Pv=v_2$. We say that $P$ is the
projection on Range$\, P$ along Ker$\, P$.

Also, if $V=W_1\oplus W_2$, then there is a unique projection $P_1$ on
$W_1$ along $W_2$ and a unique projection $P_2$ on $W_2$ along $W_1$.
Moreover, $P_1+P_2=I$.

\subsection{Orthogonal projections} \label{11.14}

Let $V$ be an inner product vector space and $W\subset V$ a finite
dimensional subspace. Then the projection on $W$ along $W^{\perp}$
is called the {\em orthogonal projection} on $W$, denoted by
$P_W$.

Note: the assumption dim$\,W<\infty$ is made to ensure that $V=W\oplus
W^{\perp}$, recall \ref{9.14}.

\subsection{Theorem} \label{11.16}

Let $P$ be a projection. Then $P$ is an orthogonal projection if and
only if $P$ is selfadjoint.

\proof Let $P$ be a projection on $W_2$ along $W_1$, and $V=W_1\oplus
W_2$. For any vectors $v,w\in V$ we have $v=v_1+v_2$ and $w=w_1+w_2$
with some $v_i,w_i\in W_i$, $i=1,2$. Now, if $P$ is orthogonal, then
$\la Pv,w\ra=\la v_2,w\ra =\la v_2,w_2\ra=\la v,w_2\ra=\la v,Pw\ra$. If
$P$ is not orthogonal, then there are $v_1\in W_1$, $w_2\in W_2$ so
that $\la v_1,w_2\ra\neq 0$. Then $\la v_1,Pw_2\ra\neq 0=\la
Pv_1,w_2\ra$.

\newpage

\section{Positive definite matrices}
\label{12}

\subsection{Bilinear forms} \label{12.1}

A {\em bilinear form} on a complex vector space $V$ is a mapping
$f\colon V\times V\to\IC$ such that
\begin{align*}
    f(u_1+u_2,v) &=f(u_1,v)+f(u_2,v) \\
    f(cu,v) &=cf(u,v) \\
    f(u,v_1+v_2) &=f(u,v_1)+f(u,v_2) \\
    f(u,cv) &=\bar{c}f(u,v)
\end{align*}
for all vectors $u,v,u_i,v_i\in V$ and scalars $c\in\IC$. In other
words, $f$ is linear in the first argument and conjugate linear in the
second.

A bilinear form on a real vector space is a mapping $f\colon V\times
V\to\IR$ that satisfies the same properties, except $c$ is a real
scalar and so $\bar{c}=c$.

\subsection{Theorem} \label{12.3}

Let $V$ be a finite dimensional inner product space. Then for every
bilinear form $f$ on $V$, then there is a unique linear operator
$T\colon V\to V$ such that
$$
   f(u,v)=\la Tu,v\ra\ \ \ \ \ \ \forall u,v\in V
$$

\proof For every $v\in V$ the function $g(u)=f(u,v)$ is linear in $u$,
so by the Riesz representation theorem~\ref{11.1} there is a vector
$w\in V$ such that $f(u,v)=\la u,w\ra$. Define a map $S\colon V\to V$
by $Sv=w$. It is then a routine check that $S$ is linear. Setting
$T=S^{\ast}$ proves the existence. The uniqueness is obvious. \qed

\subsection{Corollary} \label{12.4}

\begin{itemize} \item[\textcopyright]
Every bilinear form on $\IC^n$ can be represented by $f(x,y)=\la
Ax,y\ra$ with $A\in\IC^{n\times n}$.
\item[\textregistered]
Every bilinear form on $\IR^n$ can be represented by $f(x,y)=\la
Ax,y\ra$ with $A\in\IR^{n\times n}$.
\end{itemize}

\new Bilinear forms generalize the notion of an inner product. In order
for a bilinear form to become an inner product, though, it needs two
additional properties: conjugate symmetry $f(x,y)=\overline{f(y,x)}$
and the positivity $f(x,x)>0$ for all $x\neq 0$. We will see next what
this means in terms of the matrix $A$ that defines the bilinear form
(in $\IC^n$ or $\IR^n$).

\subsection{Hermitian/symmetric forms} \label{12.5}

A bilinear form $f$ on a complex (real) vector space $V$ is called
Hermitian (resp., symmetric) if
$$
         f(u,v)=\overline{f(v,u)}\ \ \ \ \ \ \forall u,v\in V
$$
In the real case, the bar can be dropped.

\subsection{Quadratic forms}

For a Hermitian bilinear form $f$, the function $q\colon V\to \IR$
defined by $q(u)\colon =f(u,u)$ is called the {\em quadratic form}
associated with $f$. Note that $q(u)\in\IR$ even in the complex case,
because $f(u,u)=\overline{f(u,u)}$.

\subsection{Theorem} \label{12.6}

A linear operator $T\colon V\to V$ is selfadjoint if and only if the
bilinear form $f(u,v)=\la Tu,v\ra$ is Hermitian (in the real case,
symmetric).

\proof If $T$ is selfadjoint, then $f(u,v)= \la Tu,v\ra =\la u,Tv\ra=
\overline{\la Tv,u\ra} =\overline{f(v,u)}$. If $f$ is Hermitian, then
$\la u,Tv\ra =\overline{\la Tv,u\ra} =\overline{f(v,u)}=f(u,v)= \la
Tu,v\ra =\la u,T^{\ast}v\ra$, hence $T=T^{\ast}$. \qed

\new Therefore, Hermitian bilinear forms on $\IC^n$ are defined
by Hermitian matrices.

\subsection{Positive definite forms and matrices}
\label{12.9}

A Hermitian (symmetric) bilinear form $f$ on a vector space $V$ is said
to be {\em positive definite} if $f(u,u)>0$ for all $u\neq 0$.

A selfadjoint operator $T\colon V\to V$ is said to be {\em positive
definite} if $\la Tu,u\ra >0$ for all $u\neq 0$.

A selfadjoint matrix $A$ is said to be {\em positive definite} if
$\la Ax,x\ra >0$ for all $x\neq 0$ (equivalently, $x^{\ast}Ax
>0$ for all $x\neq 0$).

\new By replacing ``$>0$'' with ``$\geq 0$'', one gets {\em
positive semi-definite} bilinear forms, operators, and matrices.

\subsection{Theorem}

The following are equivalent: \begin{itemize}\item[(a)] a bilinear form
$f(u,v)$ is an inner product. \item[(b)] $f(u,v)=\la Tu,v\ra$, where
$T$ is a positive definite operator.
\item[(c)] in $\IC^n$ and $\IR^n$, $f(x,y)=\la Ax,y\ra$
with a positive definite matrix $A$.\end{itemize}

\subsection{Lemma} \label{12.10bb}

Let $A \in \IC^{n\times n}$ be a Hermitian matrix with eigenvalues
$\lambda_1, \ldots, \lambda_n$. Then
$$
    \|x\|_2^2 \min_{1\leq i\leq n} \lambda_i \leq
     \la Ax,x\ra \leq \|x\|_2^2 \max_{1\leq i\leq n} \lambda_i
$$
and each equality is attained if and only if $x$ is an eigenvector
corresponding to the smallest or largest eigenvalue, respectively, i.e. to
$$
  \lambda_{\max} = \max_{1\leq i\leq n} \lambda_i \qquad
  \lambda_{\min} = \min_{1\leq i\leq n} \lambda_i.
$$

\medskip\noindent{\em Proof}.
By Spectral Theorem~\ref{11.20}, there is an ONB consisting of
eigenvectors $u_1,\ldots,u_n$ of $A$. Then for any vector $x=\sum c_iu_i$
we have $Ax = \sum \lambda_ic_iu_i$ and
$$
     \la Ax,x\ra=\lambda_1|c_1|^2+\cdots +\lambda_n|c_n|^2
$$
which easily implies the lemma in view of $\|x\|_2^2 = \sum |c_i|^2$,
which follows from Parceval's identity \ref{9.14}. \qed

\subsection{Lemma} \label{12.10aa}

Let $A \in \IC^{n\times n}$ be a Hermitian matrix with eigenvalues
$\lambda_1, \ldots, \lambda_n$. Then
$$
   \|A\|_2 = \max_{1\leq i\leq n} |\lambda_i|.
$$

\medskip\noindent{\em Proof}. In the notation of the previous
proof,
$$
   \|Ax\|_2^2 = \la Ax,Ax\ra =\lambda_1^2|c_1|^2+\cdots +\lambda_n^2|c_n|^2
$$
(remember that $\lambda_i\in\IR$, but $c_i\in \IC$). Thus
$$
   \|Ax\|_2^2 \leq \max_{1\leq i\leq n} \lambda_i^2\sum |c_i|^2
   =\max_{1\leq i\leq n} \lambda_i^2 \|x\|_2.
$$
The equality is attained if $x$ is an eigenvector corresponding to the
largest (in absolute value) eigenvalue. This implies the lemma. \qed



\subsection{Theorem} \label{12.10}

A Hermitian $n\times n$ matrix $A \in \IC^{n\times n}$ is positive
definite iff its eigenvalues are positive.

\proof This follows from Lemma~\ref{12.10bb}. \qed

\new Similarly, $A$ is positive-semidefinite iff it is Hermitian
and all its eigenvalues are nonnegative.

\subsection{Corollary} \label{12.11}

If a matrix $A$ is positive definite, then so is $A^{-1}$. If an
operator $T$ is positive definite, then so is $T^{-1}$.

\new Interestingly, the condition $\la Ax,x\ra >0$
for all $x\in\IC^n$ implies that the matrix $A$ is Hermitian!

\subsection{Lemma} \label{12.7}

Let $A,B$ be complex matrices. If $\la Ax,x\ra=\la Bx,x\ra$ for all
$x\in \IC^n$, then $A=B$.

\proof See homework assignment.

\new Note: This lemma holds only in complex spaces, it
fails in real spaces (see again homework assignment).

\subsection{Corollary} \label{12.8}

If $A$ is a complex matrix such that $\la Ax,x\ra\in\IR$ for all
$x\in\IC^n$, then $A$ is Hermitian.

\proof $\la Ax,x\ra = \overline{\la Ax,x\ra}= \la x,Ax\ra=\la
A^{\ast}x,x\ra$, hence $A=A^{\ast}$ by \ref{12.7}. \qed

\subsection{Theorem} \label{12.13}

A matrix $A\in\IC^{n\times n}$ is positive definite iff there is a
nonsingular matrix $B$ such that $A=B^{\ast}B$.

\proof ``$\Leftarrow$'' If $A=B^{\ast}B$, then
$A^{\ast}=B^{\ast}(B^{\ast})^{\ast}=A$ and $\la Ax,x\ra=\la Bx,Bx\ra
>0$ for any $x\neq 0$, because $B$ is nonsingular.

\noindent ``$\Rightarrow$'' By \ref{11.21} and \ref{12.10},
$A=P^{-1}DP$, where $D$ is a diagonal matrix with positive diagonal
entries and $P$ a unitary matrix. If $D=\,{\rm
diag}\,\{d_1,\ldots,d_n\}$, then denote $D^{1/2}=\,{\rm
diag}\,\{\sqrt{d_1},\ldots,\sqrt{d_n}\}$. Now $A=B^2$ where
$B=P^{-1}D^{1/2}P$, and $B$ is selfadjoint by \ref{11.22}. \qed

\new{\em Remark}. A matrix $A\in\IC^{n\times n}$ is positive
semi-definite if and only if there is a matrix $B$ such that
$A=B^{\ast}B$.

\subsection{Definition} \label{17.14}

A matrix $A\in\IC^{m\times n}$ is said to have {\em full rank} if
$$
   {\rm rank}\, A=\min\{m,n\}
$$
Otherwise, $A$ is said to be {\em rank deficient}.

\subsection{Theorem} \label{16.8}

Let $A$ be a rectangular $m\times n$ matrix. Then the matrices
$A^{\ast}A$ and $AA^{\ast}$ are Hermitian and positive
semi-definite. If $m\neq n$ and $A$ has full rank, then the
smaller of the two matrices $A^{\ast}A$ and $AA^{\ast}$ is
positive definite.

\proof $\la A^{\ast}Ax,x\ra = \la Ax,Ax\ra \geq 0$. If $A$ has full
rank and $m>n$, then Ker$(A)=\{0\}$, hence $Ax\neq 0$ for any $x\neq
0$, so that $\la A^{\ast}Ax,x\ra>0$. The case $m<n$ is similar.

\subsection{Theorem} \label{12.15}

For any $A\in\IC^{m\times n}$ we have
$$
    \|A\|_2^2 = \|A^{\ast}\|_2^2
    = \|A^{\ast}A\|_2 = \|AA^{\ast}\|_2
    =\lambda_{\max}
$$
where $\lambda_{\max}$ is the largest eigenvalue of both $A^{\ast}A$
and $AA^{\ast}$.

\new {\em Remark}: this theorem gives a practical method to compute $\|A\|_2$ for
small matrices (when $m=2$ or $n=2$), since one of the two matrices
$A^{\ast}A$ and $AA^{\ast}$ is $2\times 2$, thus its eigenvalues are
easily computable.

\medskip\noindent \emph{Proof of Theorem \ref{12.15}} is long and will be
done step by step. In the proof, $\|\cdot\|$ will always denote the
2-norm.

\medskip\noindent\textbf{Lemma}. For every vector $z\in{\mathbb C}^n$ we
have $\|z\|=\max_{\|y\|=1}|\langle y,z\rangle |$.

\medskip\noindent{\em Proof}. Indeed, by the Cauchy-Schwarz inequality
$$
  |\la y,z\ra |\leq \la y,y\ra^{1/2}\la z,z\ra^{1/2} = \|z\|
$$
and the equality is attained whenever $y$ is parallel to $z$. So we can
set $y=\pm\frac{z}{\|z\|}$ and achieve the maximum. $\Box$ \medskip
\medskip

\noindent\textbf{Step 1}. To prove that $\|A\| = \|A^{\ast}\|$ we write
$$
  \|A\| = \sup_{\|x\|=1} \|Ax\| =
  \sup_{\|x\|=1}\sup_{\|y\|=1}|\langle y,Ax\rangle |=
  \sup_{\|x\|=1}\sup_{\|y\|=1}|\langle A^\ast y,x\rangle |
$$
$$
  =\sup_{\|y\|=1}\sup_{\|x\|=1}|\langle x, A^\ast y\rangle |
   = \sup_{\|x\|=1} \|A^\ast x\| =\|A^\ast\|
$$

\medskip\noindent\textbf{Step 2}. To prove that $\|A\|^2 =
\|A^{\ast}A\|$ we write
$$
  \|A^\ast A\| = \sup_{\|x\|=1} \|A^\ast Ax\| =
  \sup_{\|x\|=1}\sup_{\|y\|=1}|\langle y,A^\ast Ax\rangle |=
  \sup_{\|x\|=1}\sup_{\|y\|=1}|\langle Ay,Ax\rangle |
$$
Then again by the Cauchy-Schwarz inequality
$$
   |\langle Ay,Ax\rangle |\leq \|Ax\|\,\|Ay\| \leq \|A\|\,\|A\|=\|A\|^2
$$
hence $\|A^\ast A\| \leq \|A\|^2$. On the other hand,
$$
  \|A^\ast A\| = \sup_{\|x\|=1}\sup_{\|y\|=1}|\langle Ay,Ax\rangle |
  \geq \sup_{\|x\|=1}|\langle Ax,Ax\rangle |=\|A\|^2.
$$
Therefore,  $\|A^\ast A\| = \|A\|^2$.

\medskip\noindent\textbf{Step 3}. Using an obvious symmetry we conclude that
$\|A^\ast \|^2 = \|A A^\ast \|$


\medskip\noindent\textbf{Step 4}. By Lemma~\ref{12.10aa}, we have
$$
  \|A^\ast A\|_2 = \max |\lambda_i(A^\ast A)|
$$
Recall that $A^\ast A$ is a positive-semidefinite matrix, so its
eigenvalues $\lambda_i(A^\ast A)$ are real and $\geq 0$, hence $\max
|\lambda_i(A^\ast A)| = \lambda_{\max}(A^\ast A)$, the largest eigenvalue
of $A^\ast A$. The same argument applies to $AA^\ast$. In particular, we
see that
$$
 \lambda_{\max}(A^\ast A)=\lambda_{\max}(AA^\ast ).
$$
This completes the proof of Theorem \ref{12.15}.

\subsection{Corollary} \label{12.15aa}

If $\lambda_{\max}$ again denotes the largest eigenvalue of $A^\ast A$,
then
$$
      \|Ax\|_2=\|A\|_2\|x\|_2
      \qquad \Longleftrightarrow \qquad
      A^{\ast}Ax=\lambda_{\max}x.
$$
Hence, the supremum in \ref{8.8} is attained (on the eigenvectors of
$A^\ast A$ corresponding to $\lambda_{\max}$) and can be replaced by
maximum. Moreover, this implies that the 2-norm of a real matrix is the
same, whether it is computed in the complex space or in the real space.

\medskip\noindent \emph{Proof}. On the one hand
$$
      \|Ax\|_2^2 = \langle Ax,Ax\rangle = \langle A^\ast Ax,x\rangle
$$
and on the other hand
$$
         \|A\|^2 = \lambda_{\max},
$$
so for any vector $x$ with $\|x\|=1$ we have
$$
      \|Ax\|_2^2=\|A\|_2^2
      \qquad \Longleftrightarrow \qquad
      \langle A^{\ast}Ax,x\rangle=\lambda_{\max}.
$$
Then we use Lemma~\ref{12.10bb}. $\Box$


\newpage

\section{Singular value decomposition}

A matrix $A\in\IC^{m\times n}$ defines a linear transformation
$T\colon \IC^n\to\IC^m$. If $B$ is an ONB in $\IC^n$ and $B'$ an
ONB in $\IC^m$, then $T$ is represented in the bases $B$ and $B'$
by the matrix $U^\ast AV$, where $U\in\IC^{m\times m}$ and
$V\in\IC^{n\times n}$ are unitary matrices. The following theorem
shows that one can always find bases $B$ and $B'$ so that the
matrix $U^\ast AV$ will be diagonal.

\new Note: $D\in \IC^{m\times n}$ is said to be diagonal
if $D_{ij}=0$ for $i\neq j$. It has exactly $p=\min\{m,n\}$ diagonal
entries and can be denoted by $D=\,$diag$\{d_1,\ldots,d_p\}$.

\subsection{Singular value decomposition (SVD)}
\label{17.8}

Let $A\in\IC^{m\times n}$ have rank $r$ and let $p=\min\{m,n\}$. Then
there are unitary matrices $U\in\IC^{m\times m}$ and $V\in\IC^{n\times
n}$ and a real diagonal matrix $D=\, {\rm diag}
\{\sigma_1,\ldots,\sigma_p\}$ such that
\beq
         A=UDV^{\ast}
         \tag{\bf SVD}
\eeq
On the diagonal of $D$, exactly $r$ elements are nonzero. If we
require, additionally, that $\sigma_1\geq\cdots\geq \sigma_r>0$ and
$\sigma_{r+1}=\cdots =\sigma_p=0$, then the matrix $D$ is unique.

\proof Observe that (SVD) is equivalent to $A^{\ast}=VD^TU^{\ast}$,
hence $A$ has an SVD if and only if $A^{\ast}$ does. Without loss of
generality, we assume that $m\geq n$, then fix $l=m-n$ and use the
induction on $n$ (taking $m=n+l$).

Let $\sigma_1=\|A\|_2$. There is a unit vector $v_1\in\IC^n$ such that
$\|Av_1\|_2 = \|A\|_2$, see \ref{12.15aa}, and a unit vector $u_1\in\IC^m$
such that $Av_1=\sigma_1u_1$. Extend $v_1$ to an ONB $\{v_1,\ldots,v_n\}$
in $\IC^n$ and extend $u_1$ to an ONB $\{u_1,\ldots,u_m\}$ in $\IC^m$. Let
$V_1$ denote the matrix whose columns are $v_1,\ldots,v_n$ and $U_1$
denote the matrix whose columns are $u_1,\ldots,u_m$. Then we have
\beq
   U_1^{\ast}AV_1=S=\left[\begin{array}{cc}
   \sigma_1 & w^{\ast} \\ 0 & B \end{array}\right]
     \tag{$\ast$}
\eeq
If $n=1$, then $S$ is diagonal yielding an SVD. For $n>1$, observe
that
$$
   \|S\|_2 = \|S^\ast\|_2 \geq    \left\|\left[\begin{array}{cc}
   \sigma_1 & 0 \\ w & B^\ast \end{array}\right]
   \left[\begin{array}{c}
   1 \\ 0 \end{array}\right]\right\|_2
   = \left\|\left[\begin{array}{c}
   \sigma_1 \\ w \end{array}\right]\right\|_2
  =\sqrt{\sigma_1^2+w^\ast w}.
$$
On the other hand, the matrices $U_1$ and $V_1$ are unitary, hence
$\|S\|_2=\|A\|_2=\sigma_1$, thus $w=0$. By our inductive assumption
the matrix $B$ has an SVD $B=\hat{U}\hat{D}\hat{V}^\ast$. Now it is
easily verified that
$$
   A=U_1\begin{bmatrix} 1 & 0 \\ 0 & \hat{U} \end{bmatrix}
   \begin{bmatrix} \sigma_1 & 0 \\ 0 & \hat{D} \end{bmatrix}
   \begin{bmatrix} 1 & 0 \\ 0 & \hat{V}^\ast \end{bmatrix}V_1^\ast
$$
is an SVD of $A$.

To prove the uniqueness, observe that
$$
  A^\ast A= VD^TDV^\ast
  \qquad\text{and}\qquad
  AA^\ast=UDD^TU^\ast
$$
hence $\sigma_1^2,\ldots,\sigma_p^2$ are the eigenvalues of both
$A^\ast A$ and $AA^\ast$ (hence, these matrices have common
non-zero eigenvalues). Note also that the columns of $U$ are the
eigenvectors of $AA^\ast$ and the columns of $V$ are the
eigenvectors of $A^\ast A$, see \ref{11.22a}. \qed

\new Note: if $A\in\IR^{m\times n}$, then
\beq
            A=UDV^T
             \tag{\bf Real SVD}
\eeq
where $U\in\IR^{m\times m}$ and $V\in\IR^{n\times n}$ are
orthogonal matrices.


\subsection{Singular values and vectors}
\label{17.9}

The positive numbers $\sigma_1,\ldots,\sigma_r$ are called the {\em
singular values} of $A$. The columns $v_1,\ldots,v_n$ of the matrix $V$
(not those of $V^\ast$) are called the {\em right singular vectors} for
$A$, and the columns $u_1,\ldots,u_m$ of the matrix $U$ are called the
{\em left singular vectors} for $A$.

\subsection{Corollary} \label{17.10}

For $1\leq i\leq r$ we have
$$
       Av_i=\sigma_i u_i,\qquad
       A^{\ast}u_i=\sigma_i v_i
$$
We also have
\begin{align*}
   {\rm Ker}\, A&=\,{\rm span}\{v_{r+1},\ldots,v_n\},&
      \qquad
   {\rm Ker}\, A^{\ast}&=\,{\rm span}\{u_{r+1},\ldots,u_m\}\\
   {\rm Range}\, A&=\,{\rm span}\{u_1,\ldots,u_r\},&
      \qquad
   {\rm Range}\, A^{\ast}&=\,{\rm span}\{v_1,\ldots,v_r\}
\end{align*}
and
$$
   {\rm rank}\,A =\,{\rm rank}\, A^\ast = r.
$$

Here is an illustration:
$$
  \begin{array}{lclcl}
    & A & & A^{\ast} & \\ & & & & \\
   v_1 & \stackrel{\sigma_1}{\longrightarrow} & u_1 &
         \stackrel{\sigma_1}{\longrightarrow} & v_1 \\
   v_2 & \stackrel{\sigma_2}{\longrightarrow} & u_2 &
         \stackrel{\sigma_2}{\longrightarrow} & v_2 \\
    \vdots & \vdots & \vdots & \vdots & \vdots \\
   v_r & \stackrel{\sigma_r}{\longrightarrow} & u_r &
         \stackrel{\sigma_r}{\longrightarrow} & v_r \\
   v_{r+1} & \rightarrow 0 & u_{r+1} & \rightarrow 0 & \\
    \vdots & \vdots & \vdots & \vdots & \\
   v_n & \rightarrow 0 & u_m & \rightarrow 0 &
   \end{array}
$$

\subsection{Remarks}\label{17.10aa}

For any matrix $A\in\IC^{m\times n}$
$$
   \|A\|_2 = \|A^\ast\|_2 = \|D\|_2 = \sigma_1.
$$
If $A$ is a square invertible matrix, then
$$
   A^{-1} = VD^{-1}U^{\ast}
$$
and
$$
    \|A^{-1}\|_2 = \|D^{-1}\|_2 = \sigma_n^{-1}.
$$
If $A\in\IC^{n\times n}$ is Hermitian with eigenvalues
$\lambda_1,\ldots,\lambda_n$, then its singular values are
$|\lambda_1|,\ldots,|\lambda_n|$ (this follows from
Corollary~\ref{11.21}).

\subsection{Computing SVD}

In practice, to compute an SVD of a $m\times 2$ matrix, find the
eigenvalues $\lambda_1\geq\lambda_2$ and the corresponding unit
eigenvectors $v_1,v_2$ of the matrix $A^\ast A$, then compute
$\sigma_1=\sqrt{\lambda_1}$, $\sigma_2=\sqrt{\lambda_2}$ and
$u_1=\sigma_1^{-1} Av_1$, $u_2=\sigma_2^{-1} Av_2$, and lastly extend
$\{u_1,u_2\}$ to an ONB in $\IC^m$ arbitrarily.

\subsection{Reduced SVD} \label{SVDred}

Let $A\in\IC^{m\times n}$ with $m>n$ and rank$\,A=r$. Then there is a
matrix $\hat{U}\in\IC^{m\times n}$ with orthonormal columns, a unitary
matrix $V\in\IC^{n\times n}$ and a square diagonal matrix $\hat{D}=\,
{\rm diag} \{\sigma_1,\ldots,\sigma_n\}$ such that
\beq
     A=\hat{U}\hat{D}V^{\ast}
     \tag{reduced SVD}
\eeq
For $m<n$, the reduced SVD is similar.

\proof Use the (full) SVD $A=UDV^{\ast}$ given by \ref{17.8} and then
erase the last $m-n$ columns of $U$ and the bottom $m-n$ rows of $D$.

\subsection{Rank-one expansion} \label{17.11}

We have the following:
\beq
       A=\sum_{i=1}^r \sigma_iu_iv_i^{\ast}
        \tag{SVD expansion}
\eeq

\proof It is enough to observe that for every $v_j$
$$
       \biggl(\sum_{i=1}^r \sigma_iu_iv_i^{\ast}\biggr)v_j=\sigma_ju_j=Av_j
$$
because $v_i^{\ast}v_j=\delta_{ij}$.

\subsection{Remarks}

Recall the Frobenius norm of a matrix, cf.\ \ref{8.1}. We have
$$
   \|A\|_F^2 =\,\text{(see  homework)}\,=
   \|D\|_F^2=\sigma_1^2+\cdots+\sigma_r^2.
$$
In particular we see that $\|A\|_2 \leq \|A\|_F$, cf.\
\ref{17.10aa}.

The value of $\|A\|_F^2$ can be interpreted as the {\em energy} of
the matrix $A$. The energy is conserved under multiplication by
unitary matrices, and the SVD pulls all the energy of a matrix onto
its diagonal.

Now, for any unit vectors $u$ and $v$ we have
$$
       \|uv^{\ast}\|_F=1
$$
hence the SVD expansion presents $A$ as a sum of rank-one matrices so
that each partial sum captures as much energy of $A$ as possible.

\subsection{Low-rank approximation}

For any $1\leq k\leq r$, define
$$
       A_k=\sum_{i=1}^k \sigma_iu_iv_i^{\ast}
$$
Then
$$
   \sigma_{k+1} = \|A-A_k\|_2 =
   \inf_{\substack{B\in\IC^{m\times n}\\
  {\rm rank}\,B\leq k}} \|A-B\|_2
$$
(with the convention $\sigma_{r+1}=0$). Thus, $A_k$ gives the best
approximation to $A$ by rank $k$ matrices.

\proof Suppose there is some matrix $B$ with rank$\, B\leq k$ such
that $\|A-B\|_2<\sigma_{k+1}$. Then dim$({\rm Ker}\, B)\geq n-k$,
and for any nonzero vector $v\in\,{\rm Ker}\, B$
$$
   \|Av\|_2=\|(A-B)v\|_2\leq
   \|A-B\|_2\|v\|_2<\sigma_{k+1}\|v\|_2
$$
On the other hand, there is a $(k+1)$-dimensional subspace
span$\{v_1,\ldots,v_{k+1}\}$ on which
$\|Av\|_2\geq\sigma_{k+1}\|v\|_2$. Since $(n-k)+(k+1)>n$, these
two spaces have a common nonzero vector, a contradiction.

\subsection{Corollary}\label{17.20}

Let $A\in\IC^{n\times n}$ be a nonsingular matrix. Then
$$
   \min\left\{\|A-A_s\|_2\colon\
   A_s\ {\rm is}\ {\rm singular}\right \}=
    \sigma_n
$$

\proof Recall: $A_s$ is singular iff rank$A_s\leq n-1$, then use the
previous fact.

\subsection{Rank with tolerance $\varepsilon$} \label{17.17}

The {\em rank} of $A\in\IC^{m\times n}$ with {\em tolerance}
$\varepsilon>0$ (also called {\em numerical rank}) is defined by
$$
      {\rm rank}(A,\varepsilon)=
      \min_{\|E\|_2\leq\varepsilon}{\rm rank}(A+E)
$$

\medskip\noindent Note: rank$(A,\varepsilon)\leq\,$rank$\, A$. The
rank with tolerance $\varepsilon$ gives the minimum rank of $A$
under perturbations of 2-norm $\leq\varepsilon$. If $A$ has full
rank, but rank$(A,\varepsilon)<p=\min\{m,n\}$ for a small
$\varepsilon$, then $A$ is `nearly rank deficient'.

\subsection{Corollary} \label{17.18}

rank$(A,\varepsilon)$ equals the number of singular values of $A$
(counted with multiplicity) that are greater than $\varepsilon$.

\new Note: rank$(A^{\ast},\varepsilon)=\,$rank$(A,\varepsilon)$,
since $A$ and $A^{\ast}$ have the same singular values.

\subsection{Definition} \label{17.15}

The vector space $\IC^{m\times n}$ with the distance between
matrices defined by
$$
   \text{dist}(A,B) = \|A-B\|_2
$$
is a {\em metric space}. Then topological notions, like open sets,
dense sets, etc., apply. Note that if a matrix $E \in \IC^{m\times
n}$ has small components, say $|e_{ij}| < \varepsilon$ for all $1
\leq i \leq m$ and $1\leq j\leq n$, then its 2-norm is bounded by
$$
  \|E\|_2 \leq \|E\|_F\leq \varepsilon\sqrt{mn}.
$$
Thus if we perturb every component of a matrix $A\in \IC^{m\times
n}$ by less than $\varepsilon$, then the perturbed matrix will be
within distance $<\varepsilon \sqrt{mn}$ from $A$.

\subsection{Theorem} \label{17.16}

Full rank matrices make an open and dense subset of $\IC^{m\times
n}$.

\new Openness means that for any full rank matrix $A$ there
is an $\varepsilon>0$ such that $A+E$ has full rank whenever
$\|E\|_2\leq\varepsilon$. Denseness means that if $A$ is rank
deficient, then for any $\varepsilon>0$ there is $E$,
$\|E\|_2\leq\varepsilon$, such that $A+E$ has full rank.

\proof To prove denseness, let $A$ be a rank deficient matrix and
$A=UDV^{\ast}$ its SVD. For $\varepsilon>0$, put $D_{\varepsilon}
=\varepsilon I$ and $E=UD_{\varepsilon}V^{\ast}$. Then
$\|E\|_2=\|D_{\varepsilon}\|_2=\varepsilon$ and
$$
   {\rm rank}(A+E)=\,{\rm rank}(U(D+D_{\varepsilon})V^{\ast})
   =\,{\rm rank}(D+D_{\varepsilon})=\min\{m,n\}
$$
Openness follows from \ref{17.18}. Indeed, if $\sigma_{\min}>0$ is
the smallest singular value of a full rank matrix $A$, then every
matrix $B$ such that $\|A-B\|_2<\sigma_{\min}$ also has full rank.

\subsection{Theorem} \label{17.16}

Diagonalizable matrices make a dense subset of $\IC^{n\times n}$.

\proof See homework assignment.

\newpage

\section{Schur decomposition}

Recall that every complex matrix $A\in\IC^{n\times n}$ is equivalent to
a Jordan matrix. In other words, for any linear operator $T\colon
\IC^n\to\IC^n$ there exists a basis in which $T$ is represented by a
Jordan matrix. What if we restrict our interests to orthonormal bases
in $\IC^n$? In other words, to what extend one can simplify a complex
matrix by using unitary equivalence?

\subsection{Schur decomposition} \label{17.1}

Any matrix $A\in \IC^{n\times n}$ is unitary equivalent to an upper
triangular matrix $T$. Moreover, one can find $T$ is such a way that
the eigenvalues of $A$ appear in any given order on the diagonal of
$T$.

\new Note: the unitary equivalence
means that $Q^{\ast}AQ=T$ for some unitary matrix $Q$. The columns of
the matrix $Q$ are called {\em Schur vectors}.

\proof We use induction on $n$. The theorem is obvious for $n=1$.
Assume that it holds for matrices of order less than $n$. Let $\lambda$
be an eigenvalue of $A$ and let $x$ be a unit eigenvector for
$\lambda$. Let $Q_1$ be a unitary matrix whose first column is $x$
(note: such a matrix exists, because there is an ONB in $\IC^n$ whose
first vector is $x$, by \ref{9.11}, and then $Q_1$ can be constructed
so that the vectors of that ONB are columns of $Q_1$). Note that
$Q_1e_1=x$, and hence $Q_1^{\ast}x=e_1$, since $Q_1^{-1}=Q_1^{\ast}$.
Hence, we have,
$$
    Q_1^{\ast}AQ_1e_1=Q_1^{\ast}Ax=\lambda Q_1^{\ast}x=\lambda e_1
$$
so $e_1$ is an eigenvector of the matrix $Q_1^{\ast}AQ_1$ for the
eigenvalue $\lambda$. Thus,
$$
   Q_1^{\ast}AQ_1=\left[\begin{array}{cc}
   \lambda & w^\ast \\ 0 & B\end{array}\right]
$$
with some $w\in\IC^{n-1}$ and $B\in\IC^{(n-1)\times (n-1)}$. By
the inductive assumption, there is a unitary matrix
$\hat{Q}\in\IC^{(n-1)\times (n-1)}$ such that
$\hat{Q}^{\ast}B\hat{Q}=\hat{T}$, where $\hat{T}$ is upper
triangular. Let
$$
     Q=Q_1\,\left[\begin{array}{cc} 1 &
     0 \\ 0 & \hat{Q} \end{array}\right]
$$
which is a unitary matrix. Next,
\begin{align*}
      Q^{\ast}AQ &=\left[\begin{array}{cc}
        1 & 0 \\ 0 & \hat{Q}^{\ast}
          \end{array}\right]
         \left[\begin{array}{cc} \lambda & w^\ast \\ 0 & B \end{array}\right]
         \left[\begin{array}{cc} 1 & 0 \\ 0 & \hat{Q} \end{array}\right ]\\
          &=\left[\begin{array}{cc} \lambda & w^\ast \hat{Q} \\
          0 & \hat{Q}^{\ast}B\hat{Q} \end{array}\right]
          =\left[\begin{array}{cc} \lambda & w^\ast \hat{Q} \\
          0 & \hat{T} \end{array}\right]
\end{align*}
which is upper triangular, as required.

\subsection{Normal matrices} \label{17.2}

A matrix $A\in\IC^{n\times n}$ is said to be {\em normal} if
$AA^{\ast}=A^{\ast}A$.

\new Note: unitary and Hermitian matrices are normal.

\subsection{Lemma} \label{17.3}

If $A$ is normal and $Q$ unitary, then $B=Q^{\ast}AQ$ is also normal
(i.e., the class of normal matrices is closed under unitary
equivalence).

\proof Note that $B^{\ast}=Q^{\ast}A^{\ast}Q$ and
$BB^{\ast}=Q^{\ast}AA^{\ast}Q=Q^{\ast}A^{\ast}AQ=B^{\ast}B$.

\subsection{Lemma} \label{17.4}

If $A$ is normal and upper triangular, then $A$ is diagonal.

\proof We use induction on $n$. For $n=1$ the theorem is trivial.
Assume that it holds for matrices of order less than $n$. Compute the
top left element of the matrix $AA^{\ast}=A^{\ast}A$. On the one hand,
it is
$$
     \sum_{i=1}^n a_{1i}\bar{a}_{1i}=
     \sum_{i=1}^n |a_{1i}|^2
$$
On the other hand, it is just $|a_{11}|^2$. Hence,
$a_{12}=\cdots=a_{1n}=0$, and
$$
     A=\left[\begin{array}{cc} a_{11} & 0 \\ 0 & B\end{array}\right]
$$
One can easily check that $AA^{\ast}=A^{\ast}A$ implies
$BB^{\ast}=B^{\ast}B$. By the inductive assumption, $B$ is diagonal.
\qed

\new Note: Any diagonal matrix is obviously normal.

\subsection{Theorem} \label{17.5}

A matrix $A\in\IC^{n\times n}$ is normal if and only if it is unitary
equivalent to a diagonal matrix. In that case the Schur decomposition
takes form $$Q^{\ast}AQ=D$$ where $D$ is a diagonal matrix, and the
columns of $Q$ (Schur vectors) become eigenvectors of $A$


\proof It follows from \ref{17.1}--\ref{17.4} and  \ref{11.22a}. \qed


\subsection{Remark} \label{17.6}

Three classes of complex matrices have the same property: they are
unitary equivalent to a diagonal matrix (i.e., admit an ONB consisting
of eigenvectors). The difference between those classes lies in
restrictions on the eigenvalues: unitary matrices have eigenvalues on
the unit circle ($|\lambda|=1$), Hermitian matrices have real
eigenvalues ($\lambda\in\IR$), and now normal matrices have arbitrary
complex eigenvalues.

\subsection{Real Schur decomposition} \label{17.7}

If $A\in\IR^{n\times n}$, then there exists an orthogonal matrix $Q$
such that
$$
     Q^TAQ=\left[\begin{array}{cccc}
     R_{11} & R_{12} & \cdots & R_{1m} \\
     0 & R_{22} & \cdots & R_{2m} \\
     \vdots & \vdots & \ddots & \vdots \\
     0 & 0 & \cdots & R_{mm}
     \end{array}\right]
$$
where each diagonal block $R_{ii}$ is either a $1\times 1$ matrix or a
$2\times 2$ matrix.

\proof We use the induction on $n$. If the matrix $A$ has a real
eigenvalue, the we can reduce the dimension and use induction just like
in the proof of Schur Theorem~\ref{17.1}. If $A$ has no real
eigenvalues, then by Lemma~\ref{10.11a} there is a two-dimensional
subspace $W\subset \IR^2$ invariant under $A$. Let $\{x_1,x_2\}$ be an
ONB of $W$. The invariance under $A$ implies
\begin{align*}
         Ax_1 &= r_{11}\, x_1 + r_{21}\, x_2 \\
         Ax_2 &= r_{12}\, x_1 + r_{22}\, x_2
\end{align*}
with some $r_{ij}\in\IR$. We now extend $\{x_1,x_2\}$ to an ONB
$\{x_1,\ldots,x_n\}$ in $\IR^n$ and denote by $\tilde{Q}$ the
orthogonal matrix with columns $x_1,\ldots,x_n$. Observe that
$\tilde{Q}e_1=x_1$ and $\tilde{Q}e_2=x_2$, hence $\tilde{Q}^T x_1
=e_1$ and $\tilde{Q}^T x_2 =e_2$. Therefore,
$$
   \tilde{Q}^T A\tilde{Q}e_1
   =\tilde{Q}^T Ax_1
   =\tilde{Q}^T (r_{11}\, x_1 + r_{21}\, x_2)
   =r_{11}\, e_1 + r_{21}\, e_2
$$
and similarly
$$
   \tilde{Q}^T A\tilde{Q}e_2
   =\tilde{Q}^T Ax_2
   =\tilde{Q}^T (r_{12}\, x_1 + r_{22}\, x_2)
   =r_{12}\, e_1 + r_{22}\, e_2
$$
Thus,
$$
   \tilde{Q}^T A\tilde{Q}
      =\left[\begin{array}{cc}
      R_{11} & \tilde{R}_{12} \\ 0 & \tilde{R}_{22}
      \end{array}\right]
$$
where
$$
     R_{11}
      =\left[\begin{array}{cc}
      r_{11} & r_{12} \\ r_{21} & r_{22}
      \end{array}\right]
$$
$\tilde{R}_{12}$ is some $2\times (n-2)$ matrix and $\tilde{R}_{22}$ is
some $(n-2)\times (n-2)$ matrix. Now we can apply our inductive
assumption to the $(n-2)\times (n-2)$ matrix $\tilde{R}_{22}$.

\newpage

\section{Gaussian elimination and LU decomposition}

\subsection{Gaussian Elimination}

Let $A\in\IC^{n\times n}$ matrix with $a_{11} \ne 0$. Denote
$A^{(1)}=A$ and $a_{ij}^{(1)}=a_{ij}$. We define {\em multipliers}
$$
     m_{i1} = a^{(1)}_{i1}/a^{(1)}_{11}\qquad
     \text{for}\ \ i=2,\ldots,n
$$
and replace the $i$-th row $a_i'$ of the matrix $A$ with $a_i' -
m_{i1}a_1'$ for all $i=2,\ldots,n$. This creates zeros in the first
column of ${A}^{(1)}$, which then takes the form
$$
    {A}^{(2)} = \left[
     \begin{array}{cccc}
     a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)} \\
     0 & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} \\
    \vdots & \vdots& \ddots & \vdots \\
    0 & a_{n2}^{(2)} & \cdots & a_{nn}^{(2)}
    \end{array}\right]
$$
where
$$
    a_{ij}^{(2)} = a_{ij}^{(1)}-m_{i1} a_{1j}^{(1)}
    \qquad \text{for}\ \ 2\le i,j \le n
$$
Next, assume that $a_{22}^{(2)} \ne 0$. Then we can continue this
process and define multipliers
$$
   m_{i2} = a_{i2}^{(2)}/a_{22}^{(2)}\qquad \text{for}\ \ i=3,\ldots,n
$$
and replace the $i$-th row $a_i'$ of the matrix $A^{(2)}$ with $a_i' -
m_{i2}a_2'$ for all $i=3,\ldots,n$. This creates a matrix, ${A}^{(3)}$,
with zeros in the second column below the main diagonal, and so on. If
all the numbers $a_{ii}^{(i)}$, $1\le i \le n-1$, are different from
zero, then one ultimately obtains an upper triangular matrix
$$
    A^{(n)}=U = \left[
    \begin{array}{ccccc}
     a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \\
     0 & a_{22}^{(2)} & a_{23}^{(2)} & \cdots & a_{2n}^{(2)} \\
     0 & 0 & a_{33}^{(3)} & \cdots & a_{3n}^{(3)} \\
    \vdots & \vdots & \ddots & \ddots & \vdots \\
    0 & 0 & \cdots & 0 & a_{nn}^{(n)}
    \end{array}\right]
$$
The elements $a_{ii}^{(i)}$, $1\le i \le n$, are called {\em pivots};
the method fails if (and only if) one of the pivots $a_{ii}^{(i)}$,
$1\le i \le n-1$,  is zero.

\subsection{Principal minors}

Let $A\in\IC^{n\times n}$. For $1\le k \le n$, the $k$-th principal
minor of $A$ is the $k\times k$ matrix formed by the entries in the
first $k$ rows and the first $k$ columns of $A$ (i.e., the top left
$k\times k$ block of $A$). We denote the $k$-th principal minor of $A$
by $A_k$.

\subsection{Theorem (Criterion of failure)}

Gaussian elimination fails if and only if det$\, A_k=0$ for some
$k=1,\ldots,n-1$. This is because for each $k=1,\ldots, n$
$$
    \det A_k=a^{(1)}_{11}\cdots a^{(k)}_{kk}
$$

\subsection{Gauss matrices}

Assume that $A\in\IC^{n\times n}$ has non-singular principal minors up
to the order $n-1$, so that the Gaussian elimination works. For each
$j=1,\ldots,n-1$ the Gauss matrix $G_j$ is defined by
$$
    {G}_j = \left[\begin{array}{cccclcc}
     1 & 0 & \cdots & 0 & \cdots & 0 & 0 \\
     0 & 1 & \cdots & 0 & \cdots & 0 & 0 \\
     \vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \vdots \\
     0 & 0 & \ddots & 1 & \cdots & \vdots & \vdots \\
     \vdots & \vdots & \ddots & -m_{j+1,j} & \ddots& \vdots & \vdots\\
     \vdots & \vdots & \ddots & \vdots & \vdots & 1 & 0\\
     0 & 0 & \cdots & -m_{n,j} & \cdots & 0 & 1
      \end{array}\right]
$$
Note that $G_j = I - m^{(j)} e_j^\ast$ where
$$
    m^{(j)} = \left[\begin{array}{c}
      0\\ \vdots \\ 0 \\ m_{j+1,j} \\ \vdots \\ m_{n,j}
    \end{array}\right]
$$

\subsection{Lemma}

For each $j=1,\ldots,n-1$ we have ${G}_j {A}^{(j)} = {A}^{(j+1)}$, and
therefore
$$
    U=A^{(n)}={G}_{n-1} \cdots {G}_2 {G}_1 {A}
$$

\subsection{Lemma}

For each $j=1,\ldots,n-1$ we have
$$
     L_j:={G}^{-1}_j =I + m^{(j)} e_j^\ast
$$
so that
$$
    {L}_j = \left[\begin{array}{cccclcc}
     1 & 0 & \cdots & 0 & \cdots & 0 & 0 \\
     0 & 1 & \cdots & 0 & \cdots & 0 & 0 \\
     \vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \vdots \\
     0 & 0 & \ddots & 1 & \cdots & \vdots & \vdots \\
     \vdots & \vdots & \ddots & m_{j+1,j} & \ddots& \vdots & \vdots\\
     \vdots & \vdots & \ddots & \vdots & \vdots & 1 & 0\\
     0 & 0 & \cdots & m_{n,j} & \cdots & 0 & 1
      \end{array}\right]
$$
and
$$
    L: = (I+m^{(1)} e_1^\ast)(I+m^{(2)} e_2^\ast) \cdots
                (I+m^{(n-1)} e_{n-1}^\ast)
    = I + \sum_{k=1}^{n-1} m^{(k)} e_k^\ast
$$
so that
$$
   L = \left[\begin{array}{ccccl}
    1  & 0 & 0 & \cdots & 0\\
    m_{21} & 1 & 0 & \cdots & 0\\
    \vdots & \ddots & \ddots & \ddots & \vdots \\
    m_{n-1,1}& m_{n-1,2} & \cdots & 1 & 0\\
    m_{n1} & m_{n2} & \cdots & m_{n,n-1} & 1
     \end{array}\right]
$$
is a unit lower triangular matrix. \\

\subsection{Theorem (LU decomposition)}

Let $A\in\IC^{n\times n}$ have non-singular principal minors up to
the order $n-1$. Then there is a decomposition
$$
      A=LU
$$
where $L$ is a unit lower triangular matrix and $U$ is an upper
triangular matrix. In that case
$$
   \det A=\det U=u_{11}\cdots u_{nn}
$$
If in addition, $A$ is non-singular, then the LU decomposition is
unique.

\noindent {\em Proof}. The existence is the result of Gaussian
elimination. To prove uniqueness, we argue by way of contradiction. Let
$A=\tilde{L} \tilde{U} = L U$. As $A$ is non-singular, it follows that
$\tilde{U}$ is also non-singular, and hence ${L}^{-1} \tilde{L} =
U\tilde{U}^{-1}$. Now, ${L}^{-1}$ is  unit lower triangular, so that
${L}^{-1} \tilde{L}$ is also unit lower triangular. On the other hand,
$U\tilde{U}^{-1}$ is upper triangular. The only matrix with these
properties is the identity matrix $I$, and thus $\tilde{L} = L$ and
$\tilde{U}=U$.

\subsection{Forward and backward substitutions}

Assume that $A\in\IC^{n\times n}$ is nonsingular and is decomposed as
$A=LU$, where $L$ is lower triangular and $U$ upper triangular. To
solve a system of equations $A x = b $, one writes it as $LUx=b$ and
then solves it in two steps:\medskip

\noindent{\bf Step 1}. Denote $Ux=y$ and solve the lower triangular
system $Ly = b $ for $y$ via ``forward substitution'' (finding
$y_1,\ldots,y_n$ subsequently).\medskip

\noindent {\bf Step 2}. Solve the system $Ux=y$ for $x$ via ``back
substitution'' (finding $x_n,\ldots,x_1$ subsequently).

\subsection{Cost of computation}

The cost of computation is measured in ``flops'', where a flop ({\em
fl}oating point {\em op}eration) is  an arithmetic operation (addition,
subtraction, multiplication, division, or a root extraction). Let us
estimate the cost of the LU decomposition. The cost of computation of
$A^{(2)}$ is $n-1$ divisions to compute the multipliers and $2n(n-1)$
flops ($n-1$ rows with $2n$ flops per row) to make the zeros in the
first column, i.e. total of approximately $2n^2$ flops. The computation
of $A^{(3)}$ then takes $2(n-1)^2$ flops, and so on. Thus the total
computational cost for the LU factorization is
$$
2(n^2+(n-1)^2 + \cdots + 2^2) = \frac{2n(n+1)(2n+1)}{6} \approx
\frac{2n^3}{3}
$$
flops.

If one solves a system $Ax=b$, then the LU decomposition is followed by
solving two triangular systems. The cost to solve one triangular system
is about $n^2$ flops. Hence, the total cost is still $\approx 2n^3/3$.

Note that the LU decomposition takes most of the computations
required for solving a system $Ax=b$ Thus, this method is
particularly well suited to situations in which one is solving
systems $Ax = b$ for more than one vector $b$: each additional $b$
will require $\approx 2n^2$ flops.

\subsection{Computation of $A^{-1}$}

Assume that $A\in\IC^{n\times n}$ is non-singular and has non-singular
principal minors up to the order $n-1$, so that the Gaussian
elimination works. One can find the matrix $X=A^{-1}$ by solving the
system $AX=I$ for $X\in\IC^{n\times n}$. This amounts to solving $n$
systems of linear equations $Ax_k=e_k$, for $k=1,\ldots,n$, where $x_k$
stands for the $k$-th column of the matrix $X$. The computational
cost of this procedure is $2n^3/3+n\times 2n^2=8n^3/3$. \\

\subsection{Diagonally dominant matrices}

A matrix $A\in\IC^{n\times n}$ such that
$$
   |a_{ii}| > \sum_{j \ne i} |a_{ij}|
$$
for all $i$ is said to be {\em strictly row diagonally dominant}. If
$$
    |a_{jj}| > \sum_{i \ne j} |a_{ij}|
$$
for all $j$ the matrix is said to be {\em strictly column diagonally
dominant}.

\subsection{Theorem}

If a matrix $A$ is strictly row (or column) diagonally dominant, then
det$\, A_k\neq 0$ for all $1\leq k\leq n$. Hence, no zero pivots will
be encountered during Gaussian elimination.

\medskip\noindent Note that if $A$ is strictly column diagonally
dominant, then all the multipliers (i.e., the elements of $L_j$) have
absolute value less than one.

\subsection{Remarks}

It is easy to find matrices for which Gaussian elimination does not
work. For example, if $A=\left (\begin{array}{cc}0 & 1\\ 1 &
0\end{array}\right )$, then the method fails instantly. (Furthermore,
for this $A$ there is no decomposition $A=LU$ with $L$ lower triangular
and $U$ upper triangular.) In practice, when a pivot is close to zero,
then multipliers become very large and numerical computations tend to
become very inaccurate. One tries to avoid this problem by using an
appropriate pivoting strategy.

\subsection{Partial pivoting}

The idea is to avoid small (in absolute value) pivots by interchanging
rows, if necessary. At any step of Gaussian elimination, one looks for
the largest (in absolute value) element in the pivot column (at or
below the main diagonal). For a non-singular matrix it cannot happen
that all of these numbers are zero. Then the row containing that
largest element is interchanged with the current row. Now the largest
element is on the main diagonal. After that the usual elimination step
is performed.

\subsection{Remarks}

Partial pivoting ensures that all the multipliers (i.e., the elements
of $L_j$) have absolute value less than or equal to one.

If a matrix $A$ is strictly column diagonally dominant, then the
Gaussian elimination with no pivoting is equivalent to Gaussian
elimination with partial pivoting (i.e., no row interchanges
are necessary). \\

\subsection{Complete pivoting}

The method of {\em complete pivoting} involves both row and column
interchanges to make use of the largest pivot available. This method
provides additional insurance against a buildup of computational
errors.

\newpage

\section{Cholesky factorization}
\label{13}


\subsection{LDM$^{\ast}$ Decomposition} \label{13.1}

Assume that $A\in \IC^{n\times n}$ has non-singular principal
minors up to the order $n$. Then there are unique matrices $L,D,M$
such that $L,M$ are unit lower triangular and $D$ is diagonal, and
$$
   A=LDM^{\ast}
$$

\proof Let $A=LU$ be the LU decomposition of $A$ and
$u_{11},\ldots,u_{nn}$ denote the diagonal entries of $U$. Set $D=\,
{\rm diag}\{u_{11},\ldots,u_{nn}\}$. Then the matrix $M^{\ast}\colon
=D^{-1}U$ is unit upper triangular, and $A=LDM^{\ast}$.

To establish uniqueness, let $A=LDM^{\ast}=L_1D_1M^{\ast}_1$. By
the uniqueness of the LU decomposition, we have $L=L_1$. Hence,
$(D^{-1}_1D)M^{\ast}=M^{\ast}_1$. Since both $M^{\ast}$ and
$M^{\ast}_1$ are unit upper triangular, the diagonal matrix
$D^{-1}_1D$ must be the identity matrix. Hence, $D=D_1$, and then
$M=M_1$.

\subsection{Corollary} \label{13.2}

If, in addition, $A$ is Hermitian, then there exist unique unit
lower triangular matrix $L$ and diagonal $D$ with real diagonal
entries such that
$$
        A=LDL^{\ast}
$$

\proof By the previous theorem $A=LDM^{\ast}$. Then
$A=A^{\ast}=MD^{\ast}L^{\ast}$, and by the uniqueness of the
LDM$^{\ast}$ decomposition we have $L=M$ and $D=D^{\ast}$.

\subsection{Sylvester's Theorem} \label{13.3}

Let $A\in\IC^{n\times n}$ be Hermitian. Then $A$ is positive definite
if and only if $\det A_k>0$ for all $k=1,\ldots, n$.

\proof Let $A$ be positive definite. By \ref{12.13}, $\det A=\det
B\times\det B^{\ast}=|\det B|^2>0$. Any principal minor $A_k$ is also a
symmetric positive definite matrix, therefore by the same argument
$\det A_k>0$. Conversely, let $\det A_k>0$. By Corollary~\ref{13.2} we
have $A=LDL^{\ast}$. Denote by $L_k$ and $D_k$ the $k$-th principal
minors of $L$ and $D$, respectively. Then $A_k= L_kD_kL_k^{\ast}$. Note
that $\det D_k=\det A_k>0$ for all $k=1,\ldots,n$, therefore all the
diagonal entries of $D$ are real and positive. Lastly, $\la Ax,x\ra=\la
DL^{\ast} x,L^{\ast} x\ra= \la Dy,y\ra >0$ because $y=L^{\ast}x\neq 0$
whenever $x\neq 0$.


\subsection{Corollary} \label{13.4}

Let $A$ be a positive definite matrix. Then $A_{ii}>0$ for all
$i=1,\ldots,n$. Furthermore, let $1\leq i_1<i_2<\cdots <i_k\leq
n$, and let $A'$ be the $k\times k$ matrix formed by the
intersections of the rows and columns of $A$ with numbers
$i_1,\ldots,i_k$. Then $\det A'>0$.

\proof Just reorder the coordinates in $\IC^n$ so that $A'$ becomes a
principal minor.

\subsection{Cholesky Factorization} \label{13.5}

Let $A\in\IC^{n\times n}$ be Hermitian and positive definite. Then
there exists a unique lower triangular matrix $G$ with real positive
diagonal entries such that
$$
        A=GG^{\ast}
$$

\proof By Corollary~\ref{13.2} we have $A=LDL^{\ast}$. Let $D=\,
{\rm diag}\{d_1,\ldots,d_n\}$. As it was shown in the proof of
\ref{13.3}, all $d_i$ are real and positive. Let $D^{1/2}=\, {\rm
diag}\{\sqrt{d_1},\ldots,\sqrt{d_n}\}$. Then $D=D^{1/2}D^{1/2}$
and setting $G=LD^{1/2}$ gives $A=GG^{\ast}$. The diagonal entries
of $G$ are $\sqrt{d_1}, \ldots,\sqrt{d_n}$, so they are positive.
To establish uniqueness, let $A=GG^{\ast}={G}_1{G}_1^{\ast}$. Then
${G}_1^{-1}G=G_1^{\ast}(G^{\ast})^{-1}$. Since this is the
equality of a lower triangular matrix and an upper triangular one,
then both matrices are diagonal:
${G}_1^{-1}G=G_1^{\ast}(G^{\ast})^{-1}=D$. Hence, ${G}_1=GD^{-1}$
and ${G}_1^{\ast}=DG^{\ast}\Rightarrow {G}_1=GD$. Therefore,
$D=D^{-1}$, so the diagonal entries of $D$ are all $\pm 1$. Note
that the value $-1$ is not possible since the diagonal entries of
both $G$ and $G_1$ are real and positive. Hence, $D=I$, and so
$G_1=G$.

\subsection{Algorithm for Cholesky factorization} \label{13.6}

Here we outline the algorithm for computing the matrix $G=(g_{ij})$
from the matrix $A=(a_{ij})$, in the \textbf{real case},
$A\in\IR^{n\times n}$. Note that $G$ is lower triangular, so
$g_{ij}=0$ for $i<j$. Hence,
$$
   a_{ij}=\sum_{k=1}^{\min\{i,j\}}g_{ik}g_{jk}
$$
Setting $i=j=1$ gives $a_{11}=g_{11}^2$, so
$$
       g_{11}=\sqrt{a_{11}}
$$
(remember, $g_{ii}$ must be positive). Next, for
$2\leq i\leq n$ we have $a_{i1}=g_{i1}g_{11}$, hence
$$
      g_{i1}=a_{i1}/g_{11}\ \ \ \ \ \ i=2,\ldots,n
$$
This gives the first column of $G$. Now, inductively,
assume that we already have the first $j-1$ columns
of $G$. Then $a_{jj}=\sum_{k=1}^j g_{jk}^2$, hence
$$
   g_{jj}=\sqrt{a_{jj}-\sum_{k=1}^{j-1}g_{jk}^2}
$$
Next, for $j+1\leq i\leq n$ we have $a_{ij}=
\sum_{k=1}^jg_{ik}g_{jk}$, hence
$$
   g_{ij}=\frac{1}{g_{jj}}\left ( a_{ij}-\sum_{k=1}^{j-1}g_{ik}g_{jk}\right )
$$
\bigskip

\subsection{Cost of computation} \label{13.7}

The computation of $g_{ij}$ takes $\approx 2j$ flops for each
$i=j,\ldots,n$, so the total is
$$
   \sum_{j=1}^n 2j(n-j)\approx 2n\frac{n^2}{2}-
   2\,\frac{n^3}{3}=\frac{n^3}{3}
$$
Recall that the LU decomposition takes about $2n^3/3$ flops, so
the Cholesky factorization is nearly twice as efficient. It is
also more stable than the LU decomposition, see Chapter~\ref{15}.

\subsection{Remark} \label{13.8}

The above algorithm can be used to verify that a given symmetric
matrix, $A$, is positive definite. Whenever the square root extractions
in \ref{13.6} are all possible and non zero, i.e. whenever
$$
   a_{11}>0\ \ \ \ {\rm and}\ \ \ \
   a_{jj}-\sum_{k=1}^{j-1}g_{jk}^2>0 \ \ \ \forall j\geq 2
$$
the matrix $A$ is positive definite.

\newpage

\section{QR decomposition}
\label{QR}

\subsection{Gram-Schmidt orthogonalization (revisited)}
\label{QR1}

Let $v_1,\ldots,v_n\in\IC^n$ be a basis. The Gram-Schmidt
orthogonalization, see \ref{9.11}, gives an ONB $\{u_1,\ldots,u_n\}$ in
$\IC^n$ such that
\begin{align*}
     v_1 &= r_{11}u_1, \\
     v_2 &= r_{12}u_1 + r_{22}u_2, \\
     v_3 &= r_{13}u_1 + r_{23}u_2 + r_{33}u_3, \\
     &\quad \dots \\
     v_n &= r_{1n}u_1 + r_{2n}u_2 + \cdots + r_{n-1,n}u_n+r_{nn}u_n
\end{align*}
where $r_{ik}=\la v_k,u_i\ra$ for $i<k\leq n$ and
$r_{kk}=\|w_{k}\|$. Note that $r_{kk}>0$.

\subsection{QR decomposition}
\label{QR4}

For any $A\in\IC^{m\times n}$ with $m\geq n$ there exist a unitary
matrix $Q\in\IC^{m\times m}$ and an upper triangular matrix
$R\in\IC^{m\times n}$ (i.e., $r_{ij}=0$ for all $i>j$) such that
$$
         A=QR
$$
If, in addition, $A$ has full rank (i.e., rank$\,A=n$), then $Q$ can be
chosen so that the diagonal entries of $R$ are real and positive.

\proof First, assume that $m=n$ and $A$ is nonsingular. Let
$v_1,\ldots,v_n$ be the columns of $A$. Since they make a basis in
$\IC^n$, we can apply Gram-Schmidt orthogonalization and obtain a
system of equations \ref{QR1}. Let $Q$ be the matrix whose columns are
$u_1,\ldots,u_n$ and
$$
     R=\left[\begin{array}{cccc}
     r_{11} & r_{12} & \cdots & r_{1n} \\
     0 & r_{22} & \cdots & r_{2n} \\
     \vdots & \vdots & \ddots & \vdots \\
     0 & 0 & \cdots & r_{nn}
     \end{array}\right]
$$
Then the above equations can be rewritten in a matrix form
$$
            A=QR
$$
Since the columns of $Q$ make an ONB, it is a unitary matrix.

Next, assume that $m=n$ but the matrix $A$ is singular, so that its
columns $v_1,\ldots,v_n$ are linearly dependent. Then the Gram-Schmidt
algorithm \ref{9.11} can be adjusted accordingly. In that case, for
some $k$ we may have $v_k\in\,$span$\{v_1,\ldots,v_{k-1}\}$, and then
$w_k=0$, in the notation of \ref{9.11}. Now the vector $u_k=
w_k/\|w_k\|$ cannot be determined. Instead, we need to set $u_k$ to an
arbitrary unit vector orthogonal to $u_1,\ldots,u_{k-1}$, and then
continue the calculations as described in \ref{9.11}. In that case, in
the equations of \ref{9.11} we get $r_{k,k}=0$.

Now let $m>n$. We extend the matrix $A$ on the right by adding $m-n$
columns consisting of zeroes. Denote by $A'$ the resulting $m\times m$
matrix with columns
$$
     v_1,v_2,\ldots,v_{n-1},v_n,\ v_{n+1}=0,\ldots,v_m=0
$$
This is a square matrix, so it has a QR decomposition
$$
            A'=Q'R'
$$
The last $m-n$ columns of the matrix $A$ are zeroes, so we obtain
$$
      0=v_k=\sum_{i=1}^k r_{1k}'u_k'
$$
for all $k=n+1,\ldots,m$. Since $u_1',\ldots,u_m'$ are basis vectors,
$r_{ik}'=0$ for all $k>n$ and $1\leq i\leq m$, hence the last $m-n$
columns of $R'$ are zeroes.

Now we erase the last $m-n$ columns of the matrices $A'$ and $R'$ to
obtain the desired QR decomposition
$$
         A=QR
$$
Theorem \ref{QR4} is proved. \qed

\new Note: if $A$ is a real matrix, then the matrices $Q$
and $R$ are also real, and so $Q$ is orthogonal.

\subsection{Reduced (``skinny'') QR decomposition}
\label{QR5}

Let $A\in\IC^{m\times n}$ with $m\geq n$. Then there is a matrix
$\hat{Q}\in\IC^{m\times n}$ with orthonormal columns and an upper
triangular matrix $\hat{R}\in\IC^{n\times n}$ such that
$$
         A=\hat{Q}\hat{R}
$$

\proof By Theorem~\ref{QR4}, $A=QR$. Let $\hat{Q}$ be the left $m\times
n$ rectangular block of $Q$ (the first $n$ columns of $Q$). Let
$\hat{R}$ be the top $n\times n$ square block of $R$ (note that the
remainder of $R$ is zero). Then $A=\hat{Q}\hat{R}$.

\subsection{Corollary} \label{QR6}

If, in addition, $A$ has full rank (rank$\, A=n$), then\begin{itemize}
\item[(a)] The columns of $\hat{Q}$ make an ONB in the column space of
$A$ (this is the subspace in $\IC^n$ spanned by the columns of $A$).
\item[(b)] One can find $\hat{Q}$ so that the diagonal entries of
$\hat{R}$ will be real and positive ($r_{ii}>0$).\item[(c)] The
$\hat{Q}$ and $\hat{R}$ described in part (b) are unique. \end{itemize}

\proof (a) and (b) follow from \ref{QR4}. To prove (c), let
$A=\hat{Q}\hat{R}=\hat{Q}_1\hat{R}_1$. Then $A^{\ast}A=
\hat{R}^{\ast}\hat{R}=\hat{R}^{\ast}_1\hat{R}_1$. Since
$A^{\ast}A$ is positive definite, we can use the uniqueness of
Cholesky factorization and obtain $\hat{R}=\hat{R}_1$. Then also
$\hat{Q}_1=\hat{Q}\hat{R}\hat{R}_1^{-1}=\hat{Q}$. \qed

\subsection{Cost of computation} \label{QR7}

In order to compute the reduced QR decomposition \ref{QR5}, one needs
to apply Gram-Schmidt orthogonalization to the $n$ columns of the
matrix $A$ and compute the $n$ columns $u_1,\ldots,u_n$ of the matrix
$\hat{Q}$. The vector $u_k$ is found by
$$
   w_{k}=v_{k}-\sum_{i=1}^{k-1}\la v_{k},u_i\ra\, u_i,
   \quad \text{and} \quad
   u_k=\frac{w_k}{\|w_k\|},
$$
see \ref{9.11}. Here each scalar product $\la v_k,u_i\ra$ requires
$m$ multiplications and $m$ additions, and then subtracting every
term $\la v_k,u_i\ra\, u_i$ from $v_k$ requires $m$ multiplication
and $m$ subtractions, for each $i=1,\ldots,k$. The total is $4mk$
flops. The subsequent computation of $\|w_{k+1}\|$ and then
$u_{k+1}$ requires $3m$ flops, which is a relatively small number,
and we ignore it. The total flop count is
$$
   \sum_{k=1}^n 4mk \approx 2mn^2
$$

\subsection{Modified Gram-Schmidt orthogonalization}

The algorithm \ref{9.11} is often called {\em classical} Gram-Schmidt
orthogonalization, as opposed to the {\em modified} Gram-Schmidt
orthogonalization we present next. Given a basis $\{v_1,\ldots,v_n\}$
in $V$, we denote $v_i^{(1)}=v_i$ for $i=1,\ldots,n$, then compute
$$
           u_1=v_1^{(1)}/\|v_1^{(1)}\|,
$$
and then modify all the remaining vectors by the rule
$$
    v_i^{(2)} = v_i^{(1)} - \la v_i^{(1)},u_1 \ra\, u_1
    \qquad \text{for}\ \ 1<i\leq n.
$$
After that, inductively, for each $k\geq 2$ we compute
$$
   u_k=v_k^{(k)}/\|v_k^{(k)}\|,
$$
and then modify all the remaining vectors by the rule
$$
    v_i^{(k+1)} = v_i^{(k)} - \la v_i^{(k)},u_k \ra\, u_k
    \qquad \text{for}\ \ k<i\leq n.
$$
The modified and classical Gram-Schmidt methods produce the same
orthogonal basis $\{u_1,\ldots,u_n\}$ (i.e., they are mathematically
equivalent). They are based on a different logic, though.

The classical Gram-Schmidt computes $u_k$ by making the current vector
$v_k$ orthogonal to all the previously constructed vectors
$u_1,\ldots,u_{k-1}$, without touching the remaining vectors
$v_{k+1},\ldots,v_n$. The amount of work increases as $k$ grows from 1
to $n$. This is a ``lazy man schedule'' - do as little as possible and
leave the rest of the work ``for later''.

The modified Gram-Schmidt computes $u_k$ and makes all the remaining
vectors $v_{k+1},\ldots,v_n$ orthogonal to it. Once this is done, the
remaining vectors will be in the orthogonal complement to the subspace
span$\{u_1,\ldots,u_k\}$ and there is no need to involve the previously
constructed vectors anymore. The amount of work decreases as $k$ grows
from 1 to $n$. This is an ``industrious man schedule'' - do as much as
possible now and reduce the workload.

Overall, both methods require the same amount of flops. The modified
Gram-Scmidt has an important advantage that it gives more accurate
numerical results in computer calculations, see programming assignment.


\subsection{Computation of SVD} \label{SVDcom}

We note that there is no finite algorithms for the computation of SVD
decomposition \ref{17.8} or its reduced version \ref{SVDred} (except for
small matrices, see Remark~\ref{12.7}). The reason will be discussed in
Section~\ref{19}. In practice, SVD is computed by special iterative
algorithms. The computation of reduced SVD require approximately
$$
   2mn^2 + 11n^3\qquad {\rm flops}.
$$

\newpage


\section{Overdetermined linear systems}
\label{16}

\subsection{Definition}
\label{16.1}

A system of linear equations $Ax=b$ with $A\in\IC^{m\times n}$, $x\in
\IC^n$ and $b\in\IC^m$, is said to be {\em overdetermined} if $m>n$.
Since there are more equations than unknowns, the system usually has no
solutions.

\new Note that the matrix $A$ defines a linear
transformation $\IC^n\to\IC^m$. Then $Ax=b$ has a solution if and only
if $b\in\,{\rm Range}\,A$. In this case the solution is unique if and
only if Ker$\, A=\{0\}$, i.e. rank$\, A=n$.

\subsection{Least squares solution} \label{16.5}

Let $Ax=b$ be an overdetermined linear system. A vector $x\in\IC^n$
that minimizes the function
$$
        E(x)=\|b-Ax\|_2
$$
is called a {\em least squares solution} of $Ax=b$. The vector
$r=b-Ax$ is called the {\em residual vector} and $\|r\|_2$ the
{\em residual norm}.

\subsection{Normal equations} \label{16.7}

Let $Ax=b$ be an overdetermined linear system. Then the linear
system
$$
        A^{\ast}Ax=A^{\ast}b
$$
is called the {\em system of normal equations} associated with $Ax=b$.

\subsection{Theorem} \label{16.9}

Let $Ax=b$ be an overdetermined linear system. Then
\begin{itemize}
\item[(a)] A vector $x$ minimizes $E(x)=\|b-Ax\|_2$ if and only if
it is an exact solution of the system $Ax=\hat{b}$, where $\hat{b}$ is
the orthogonal projection of $b$ onto Range$\, A$.
\item[(b)] A vector $x$ minimizing $E(x)$ always exists. It is
unique if and only if $A$ has full rank.
\item[(c)] A vector $x$ minimizes $E(x)$ if and only if it is
a solution of the system of normal equations $A^{\ast}Ax
=A^{\ast}b$.
\end{itemize}

\proof Recall that $({\rm Range}\, A)^\perp =\,{\rm Ker}\, A^{\ast}$ by
\ref{11.7}. So, we have an orthogonal decomposition $\IC^m={\rm
Range}\,A\oplus \,{\rm Ker}\, A^{\ast}$, in particular $b=\hat{b}+r$,
where $\hat{b}\in{\rm Range}\, A$ and $r\in{\rm Ker}\, A^{\ast}$. Since
$\la \hat{b},r \ra=0$, Pythagorean Theorem implies
$$
   \|b-Ax\|_2^2=\|\hat{b}-Ax\|_2^2+\|r\|_2^2\geq \|r\|_2^2
$$
Hence, $\min_x E(x)=\|r\|_2$ is attained whenever $Ax=\hat{b}$. Since
$\hat{b}\in{\rm Range}\, A$, there is always an $x\in\IC^n$ such that
$Ax=\hat{b}$. The vector $x$ is unique whenever Ker$\, A=\{0\}$, i.e.\
whenever $A$ has full rank. This proves (a) and (b).

To prove (c), observe that if $x$ minimizes $E(x)$, then $r\in{\rm
Ker}\, A^{\ast}$, and thus $A^{\ast}(b-Ax)=0$. Conversely, if
$A^{\ast}(b-Ax)=0$, then $r\in{\rm Ker}\, A^{\ast}$, and by the
uniqueness of the decomposition $b=\hat{b}+r$ we have $Ax=\hat{b}$.

\new Next we give examples that lead to overdetermined systems and
least squares problems.

\subsection{Linear least squares fit} \label{16.2}

Let $(x_i,y_i)$, $1\leq i\leq m$, be points on the $xy$ plane. For
any straight line $y=a_0+a_1x$ one defines the ``combined
distance'' of that line from the given points by
$$
        E(a_0,a_1)=\biggl[\sum_{i=1}^m (a_0+a_1x_i-y_i)^2\biggr]^{1/2}
$$
The line $y=\hat{a}_0+\hat{a}_1x$ that minimizes the function
$E(a_0,a_1)$ is called the {\em least squares fit} to the points
$(x_i,y_i)$. This is a basic tool in statistics. Let
$$
    A=\left[\begin{array}{cc} 1 & x_1 \\ 1 & x_2 \\
     \vdots & \vdots \\ 1 & x_m \end{array}\right]
   \ \ \ \ {\bf x}=\left[\begin{array}{c} a_0 \\ a_1 \end{array}\right]
   \ \ \ \ {\bf b}=\left[\begin{array}{c} y_1 \\ y_2 \\ \vdots \\
    y_m\end{array}\right]
$$
Then
$$
        E(a_0,a_1)=\|{\bf b}-A{\bf x}\|_2
$$
Hence the least squares fit is equivalent to the least squares
problem $A{\bf x}={\bf b}$.

\subsection{Polynomial least squares fit} \label{16.27}\nopagebreak

Generalizing \ref{16.2}, one can fit a set of data points
$(x_i,y_i)$, $1\leq i\leq m$, by a polynomial
$y=p(x)=a_0+a_1x+\cdots +a_nx^n$ with $n+1\leq m$. The least
squares fit is based on minimizing the function
$$
   E(a_0,\ldots,a_n)=
   \biggl[\sum_{i=1}^m\Big (a_0+a_1x_i+
   \cdots +a_nx_i^n-y_i\Big )^2\biggl]^{1/2}
$$
This is equivalent to the least squares problem for an overdetermined
linear system
$$
   a_0+a_1x_i+\cdots +a_nx_i^n=y_i\ \ \ \ \ 1\leq i\leq m
$$
in which $a_0,\ldots,a_n$ are unknowns.


\subsection{Continuous least squares fit} \label{16.28}

Instead of fitting a discrete data set $(x_i,y_i)$ one can fit a
continuous function $y=f(x)$ on $[0,1]$ by a polynomial $y=p(x)\in
P_n(\IR)$. The least squares fit is based on minimization of
$$
   E(a_0,\ldots,a_n)=
   \biggl[\int_0^1|f(x)-p(x)|^2\, dx\biggr]^{1/2}
$$
The solution of this problem is the orthogonal projection of $f(x)$
onto $P_n(\IR)$.

To find the solution, consider a basis $\{1,x,\ldots,x^n\}$ in
$P_n(\IR)$. Then $a_0,\ldots,a_n$ can be found by solving the system of
equations (analogous to normal equations)
$$
     \sum_{j=0}^n a_j\la x^j,x^i\ra =
     \la f,x^i\ra \ \ \ \ \ \ \ 1\leq i\leq n
$$
The matrix of coefficients here is
$$
     \la x^j,x^i\ra =\int_0^1 x^{i+j}\, dx=\frac{1}{1+i+j}
$$
for $0\leq i,j\leq n$.

\new Next we present three methods of solving the least square
problem.

\subsection{Algorithm 1, based on normal equations} \label{A1}

This is the simplest one:

\new 1.  Form the matrix $A^{\ast}A$ and the vector
$A^{\ast}b$.\\ 2. Compute the Cholesky factorization
$A^{\ast}A=GG^{\ast}$.\\ 3. Solve the lower-triangular system
$Gz=A^{\ast}b$ for $z$.\\ 4. Solve the upper-triangular system
$G^{\ast}x=z$ for $x$.

\new The cost of this algorithm is dominated by steps 1 and 2.
Because of symmetry, the computation of $A^{\ast}A$ requires
$mn(n+1)\approx mn^2$ flops. The computation of $A^{\ast}b$ requires
only $2mn$ flops, a relatively small amount which we ignore. The
Cholesky factorization takes $n^3/3$ flops (see \ref{13.7}), a total of
$$
    \boxed{\phantom\int\approx mn^2+\tfrac 13
    n^3\qquad{\rm flops}\phantom\int}
$$

\subsection{Algorithm 2, based on QR decomposition} \label{A2}

Using the reduced QR decomposition \ref{QR5} allows us to rewrite
the system of normal equations as
$$
    \hat{R}^\ast\hat{Q}^{\ast}\hat{Q}\hat{R}x
    = \hat{R}^\ast\hat{Q}^\ast b
$$
If $A$ has full rank, the matrix $\hat{R}^\ast$ is nonsingular and
we cancel it out. Also, since the columns of $\hat{Q}$ are
orthonormal vectors, $\hat{Q}^{\ast}\hat{Q}=I$. Hence
$$
   \hat{R}x=\hat{Q}^{\ast}b
$$
This suggests the following algorithm:

\new 1. Compute the reduced QR decomposition $A=\hat{Q}\hat{R}$\\
2. Compute the vector $\hat{Q}^{\ast}b$.\\ 3. Solve the
upper-triangular system $\hat{R}x=\hat{Q}^{\ast}b$ for $x$.

\new The cost of this algorithm is dominated by step 1, the
reduced QR factorization, which requires
$$
    \boxed{\phantom\int \approx 2mn^2\qquad
     {\rm flops}\phantom\int}
$$
see \ref{QR7}. This is approximately twice as much as Algorithm 1
requires.

\subsection{Algorithm 3, based on SVD decomposition} \label{A3}

Using the reduced SVD decomposition \ref{SVDred} allows us to
rewrite the system of normal equations as
$$
    V\hat{D}\hat{U}^\ast\hat{U}\hat{D}V^{\ast}x
    = V\hat{D}\hat{U}^\ast b
$$
The matrix $V$ is unitary and we cancel it out. If $A$ has full
rank, the matrix $\hat{D}$ is nonsingular and we cancel it out.
Since the columns of $\hat{U}$ are orthonormal vectors,
$\hat{U}^{\ast}\hat{U}=I$. Hence
$$
   \hat{D}V^{\ast}x=\hat{U}^{\ast}b
$$
This suggests the following algorithm:

\new 1. Compute the reduced SVD decomposition $A=\hat{U}\hat{D}V^{\ast}$\\
2. Compute the vector $\hat{U}^{\ast}b$.\\ 3. Solve the diagonal
system $\hat{D}z=\hat{U}^{\ast}b$ for $z$. \\ 4. Set $x=Vz$.

\new The cost of this algorithm is dominated by step 1, the
reduced SVD decomposition, which requires
$$
    \boxed{\phantom\int \approx 2mn^2+11n^3
    \qquad {\rm flops}\phantom\int}
$$
see \ref{SVDcom}. This is approximately the same amount as in Algorithm
2 for $m\gg n$, but for $n\approx m$ this algorithm is much more
expensive.

\new The algorithm 1 is the simplest and the cheapest, but it often
gives inaccurate results in numerical calculations. The algorithms 2
and 3 are more complicated and expensive, but usually give more
accurate numerical results, for the reasons we learn in the next
sections.

\subsection{The case of a rank deficient $A$} \label{Adefic}

If rank$\,A<n$, then the least squares solution is not unique: the
set of solutions is the preimage of the vector $\hat{b}$ under the
map $A\colon \IC^n\to\IC^m$. In this case the ``best'' solution is
the one of minimal norm:
$$
      Ax_{\text{best}}=\hat{b}\quad\text{and}\quad
      \|x_{\text{best}}\|_2\leq\|x\|_2\quad
      \forall x\colon Ax=\hat{b}
$$
Algorithms 1 and 2 fail to find any solution for a rank deficient
matrix $A$. On the contrary, Algorithm 3 easily finds the minimal
norm solution.

\newpage

\section{Machine arithmetic}
\label{14}

\subsection{Decimal number system} \label{14.1}

In our decimal system, natural numbers are represented by a sequence of
digits
$$
    N= (a_n\cdots a_1a_0)_{10}=10^na_n+\cdots +10a_1+a_0
$$
where $a_i\in\{0, 1, 2, 3, 4, 5, 6, 7, 8, 9\}$ are digits.
Fractional numbers require an additional {\em fractional part}:
$$
     f=.d_1d_2\ldots d_t\ldots=
     10^{-1}d_1+10^{-2}d_2+\cdots+10^{-t}d_t+\cdots
$$
which may be finite or infinite.

\subsection{Floating point representation} \label{14.2}

Alternatively, any real number can be written as a product of a
fractional part with a sign and a power of ten:
$$
     r=\bigl(\pm .d_1d_2\ldots d_t\ldots\bigr)\times 10^e
$$
where $d_i$ are decimal digits and $e\in{\mathbb Z}$ is an integer. For
example, $18.2=0.182\times 10^2=0.0182\times 10^3$, etc. This is called
a {\em floating point} representation of decimal numbers. The part
$.d_1\ldots d_t$ is called the {\em mantissa} and $e$ is the {\em
exponent}. By changing the exponent $e$ with a fixed mantissa
$.d_1\ldots d_t$ we can move (``float'') the decimal point, for example
$0.182\times 10^2=18.2$ and $0.182\times 10^1=1.82$.

\subsection{Normalized floating point representation}
\label{14.3}

To avoid unnecessary multiple representations of the same number (as
that of 18.2 by $0.182\times 10^2$ and $0.0182\times 10^3$ above), we
require that $d_1\neq 0$. We say the floating point representation is
{\em normalized} if $d_1\neq 0$. Then $0.182\times 10^2$ is the only
normalized representation of the number 18.2.

For every positive real $r>0$ there is a unique integer $e\in{\mathbb
Z}$ such that $f\colon =10^{-e}r\in [0.1,1)$. Then $r=f\times 10^e$ is
the normalized representation of $r$. For most of real numbers, the
normalized representation is unique, however, there are exceptions,
such as
$$
     \bigl(.9999\ldots\bigr)\times 10^0 =
     \bigl(.1000\ldots\bigr)\times 10^1,
$$
which admit dual normalized representations. In such cases one of the
two representation has a finite fractional part, and the other --
infinite.


\subsection{Binary number system} \label{14.4}

In the binary number system, the base is 2 (instead of 10) and digits
take two values: 0 and 1. Any natural number $N$ can be written, in the
binary system, as a sequence of binary digits:
$$
    N=(a_n\cdots a_1a_0)_2=2^na_n+\cdots +2a_1+a_0
$$
For example, $5=101_2$, $11=1011_2$, $64=1000000_2$, etc. Binary
system, due to its simplicity, is used by all computers. In the modern
computer world, a {\em bit} means a {\em b}inary dig{\em it}.

\subsection{Other number systems} \label{14.5}

Now suppose we are working in a number system with base $\beta\geq 2$.
By analogy with \ref{14.1} and \ref{14.4}, any natural number $N$ can
be written, in that system, as a sequence of digits:
$$
    N=(a_n\cdots a_1a_0)_{\beta}=\beta^na_n+\cdots +\beta a_1+a_0
$$
where $a_i\in\{0, 1,\ldots,\beta-1\}$ are digits. Fractional
numbers require an additional {\em fractional part}:
$$
     f=.d_1d_2\ldots d_t\ldots=
     \beta^{-1}d_1+\beta^{-2}d_2+\cdots+\beta^{-t}d_t+\cdots
$$
which may be finite or infinite. The floating point representation of
real numbers in the system with base $\beta$ is given by
$$
     r=\bigl(\pm .d_1d_2\ldots d_t\ldots\bigr)\times \beta^e
$$
where $.d_1\ldots d_t$ is called the {\em mantissa} and $e\in{\mathbb
Z}$ is the {\em exponent}. Again, we say that the above representation
is normalized if $d_1\neq 0$, this ensures uniqueness for almost all
real numbers.

\subsection{Machine floating point numbers} \label{14.6}

Any computer can only handle finitely many numbers. Hence, the number
of digits $d_i$'s in the mantissa is necessarily bounded, and the
possible values of the exponent $e$ are limited to a finite interval.
Assume that the length of the mantissa, $t$, is fixed, and the exponent
is bounded by $L\leq e\leq U$. Then the four integers $\beta,t,L,U$
completely characterize the set of real numbers
$$
     r=\bigl(\pm .d_1d_2\ldots d_t\bigr)\times \beta^e,
     \quad L\leq e\leq U
$$
that a given machine system can handle. Note that, curiously, zero
cannot be represented in the above format, since $d_1\neq 0$. Thus, the
real machine systems usually include a few special numbers, like zero,
that require a different form of representation. We also note that any
given computer can usually adopt many possible machine systems, with
different values of $t,L,U$, see \ref{14.8}.

\subsection{Remark} \label{14.7}

The maximal (in absolute value) number that a machine system can handle
is $M=\beta^U(1-\beta^{-t})$. The minimal positive number is
$m=\beta^{L-1}$.

\subsection{Examples}\label{14.8}

Most modern computers conform to the IEEE floating-point standard
(ANSI/IEEE Standard 754-1985), which provides two machine
systems:\begin{itemize}\item[\textcircled{\scriptsize I}] {\em Single
precision} is characterized by $t=24$, $L=-125$ and
$U=128$.\item[\textcircled{\scriptsize II}] {\em Double precision} is
characterized by $t=53$, $L=-1021$ and $U=1024$.
\end{itemize}

\subsection{Relative errors} \label{14.9}

How can real numbers be represented in a machine system
characterized by $\beta,t,L,U$? Let $x\neq 0$ be a real number
with a normalized floating point representation
$$
   x=\bigl(\pm 0.d_1d_2\ldots\bigr)\times \beta^e
$$
where the number of digits may be finite or infinite. If $e>U$ or
$e<L$, then $x$ cannot be correctly represented in the machine
system (it is either ``too large'' or ``too small''). If $e\in
[L,U]$ is within the right range, then the mantissa has to be
reduced to $t$ digits (if it is longer or infinite). There are two
standard ways to do such a reduction:
\begin{itemize} \item[(a)] keep the first $t$ digits and chop off the
rest; \item[(b)] round off to the nearest available, i.e. use the
rules
$$
  \left \{ \begin{array}{ll}
  .d_1\ldots d_t & {\rm if}\ \ d_{t+1}<\beta/2  \\
  .d_1\ldots d_t+.0\ldots 01 \ \ \ & {\rm if}\ \ d_{t+1}\geq\beta/2
      \end{array}\right .
$$
\end{itemize}
Denote the obtained number by $x_c$ (the computer representation
of $x$). The relative error in this representation can be
estimated as
$$
     \frac{x_c-x}{x}=\varepsilon
       \ \ \ \ {\rm or} \ \ \ \
     x_c=x(1+\varepsilon)
$$
where the maximal possible value of $\varepsilon$ is
$$
      {\bf u}=\left\{\begin{array}{rl}
 \beta^{1-t}\ \ \ &
 {\rm for}\ {\rm chopped}\ {\rm arithmetic}\ {\rm (a)} \\
 \frac 12\beta^{1-t}\ \ \ &
 {\rm for}\ {\rm rounded}\ {\rm arithmetic}\ {\rm (b)}
      \end{array}\right .
$$
The number $\bf u$ is called the {\em unit round off} or the {\em
machine precision}.

\subsection{Examples} \label{14.10}

a) For the IEEE floating-point single precision standard with
chopped arithmetic ${\bf u}=2^{-23}\approx 1.2\times 10^{-7}$. In
other words, approximately 7 decimal digits are accurate.
\bigskip

b) For the IEEE floating-point double precision standard with chopped
arithmetic ${\bf u}=2^{-52}\approx 2.2\times 10^{-16}$. In other words,
approximately 16 decimal digits are accurate.

\subsection{Example} \label{14.11}

Consider the system of equations
$$
  \left[\begin{array}{cc} 0.01 & 2 \\ 1 & 3 \end{array}\right]
  \left[\begin{array}{l} x \\ y \end{array}\right]=
  \left[\begin{array}{l} 2 \\ 4 \end{array}\right]
$$
The exact solution is $x=\tfrac{200}{197}\approx 1.015$ and
$y=\tfrac{196}{197}\approx 0.995$.

Let us solve this system by using a chopped arithmetic with base
$\beta=10$ and $t=2$ (i.e. working with a two digit mantissa). If
we use the Gaussian elimination without pivoting, then the
computed solution will be $x_c=0.0$ and $y_c=1.0$. The value of
$x$ is 100\% off! Increasing the length of the mantissa to $t=3$
gives $x_c=2.0$ and $y_c=0.994$, not much of improvement, since
$x_c$ is still about 100\% off. We postpone the explanation until
\ref{15.11a}.

On the other hand, if we apply partial pivoting (interchanging the
rows), then the computed solution with $t=2$ will be $x_c=1.0$,
$y_c=1.0$, and with $t=3$ it will be $x_c=1.02$, $y_c=0.994$, which is
good. The table below shows that the relative error of the numerical
solutions is proportional to minimal round-off error $10^{-t}$, with a
factor of about 2 to 5.

\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline  & relative error & min.\ error & factor \\
\hline  $t=2$ & $1.5\times 10^{-2}$  & $10^{-2}$ & $1.5$ \\
\hline  $t=3$ & $4.8\times 10^{-3}$  & $10^{-3}$ & $4.8$ \\
\hline  $t=4$ & $2.2\times 10^{-4}$  & $10^{-4}$ & $2.2$ \\
\hline
\end{tabular}\vspace*{0.2cm}
\end{center}

\new Conclusion: Gaussian elimination without pivoting may lead to
catastrophic errors and unreliable numerical solutions. Pivoting
is more reliable... but see the next example:

\subsection{Example} \label{14.11a}

Consider another system of equations:
$$
  \left[\begin{array}{cc} 3 & 1 \\ 1 & 0.35 \end{array}\right]
  \left[\begin{array}{l} x \\ y \end{array}\right]=
  \left[\begin{array}{l} 5 \\ 1.7 \end{array}\right]
$$
The exact solution here is $x=1$ and $y=2$. The largest
coefficient is at the top left corner already, so there is no need
for pivoting.

Solving this system in a chopped arithmetic with $\beta=10$ and
$t=2$ gives $x_c=0$ and $y_c=5$, which is 150\% off. Increasing
the length of the mantissa to $t=3$ gives $x_c=0.883$ and
$y_c=2.35$, so the relative error is 17\%. With $t=4$, we obtain
$x_c=0.987$ and $y_c=2.039$, now the relative error is 2\%. The
table below shows that the relative error of the numerical
solutions is proportional to minimal round-off error $10^{-t}$,
with a factor of about 150 to 200.

\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline  & relative error & min.\ error & factor \\
\hline  $t=2$ & $1.5\times 10^{-0}$  & $10^{-2}$ & $150$ \\
\hline  $t=3$ & $1.7\times 10^{-1}$  & $10^{-3}$ & $170$ \\
\hline  $t=4$ & $2.0\times 10^{-2}$  & $10^{-4}$ & $200$ \\
\hline
\end{tabular}\vspace*{0.2cm}
\end{center}

\noindent We postpone a complete analysis of these two examples
until \ref{15.11b}.

\subsection{Computational errors} \label{14.12}

Let $x,y$ be two real numbers represented in a machine system by
$x_c,y_c$. An arithmetic operation $x\ast y$, where $\ast$ is one
of $+,-,\times,\div$, is performed by a computer in the following
way. The computer finds $x_c\ast y_c$ (first, exactly) and then
represents that number by the machine system. The result is
$z=(x_c\ast y_c)_c$. Note that, generally, $z$ is different from
$(x\ast y)_c$, which is the representation of the exact result
$x\ast y$. Hence, $z$ is not necessarily the best representation
for $x\ast y$. In other words, the computer makes additional round
off errors during computations. Assuming that
$x_c=x(1+\varepsilon_1)$ and $y_c=y(1+\varepsilon_2)$ we have
$$
    (x_c\ast y_c)_c=(x_c\ast y_c)\,(1+\varepsilon_3)=
     [\,x(1+\varepsilon_1)\,]\ast [\,y(1+\varepsilon_2)\,]
     \,(1+\varepsilon_3)
$$
where $|\varepsilon_1|,|\varepsilon_2|,|\varepsilon_3|\leq\bf u$.

\subsection{Multiplication and division} \label{14.13}

For multiplication, we have
$$
   z=xy(1+\varepsilon_1)(1+\varepsilon_2)(1+\varepsilon_3)\approx
       xy(1+\varepsilon_1+\varepsilon_2+\varepsilon_3)
$$
(here we ignore higher order terms), so the relative error is
(approximately) bounded by $3{\bf u}$. A similar estimate can be
made in the case of division:
$$
   z=\frac{x(1+\varepsilon_1)(1+\varepsilon_3)}{y(1+\varepsilon_2)}\approx
       \frac xy\,(1+\varepsilon_1-\varepsilon_2+\varepsilon_3)
$$
(here we use Taylor expansion for $(1+\varepsilon_2)^{-1}$ and
ignore higher order terms). Hence, machine multiplication and
machine division increase relative errors by a factor of three, at
most.


\subsection{Addition and subtraction} \label{14.14}

For addition, we have
$$
    z=(x+y+x\varepsilon_1+y\varepsilon_2)(1+\varepsilon_3)=
    (x+y)\left (1+\frac{x\varepsilon_1+y\varepsilon_2}{x+y}\right)
    (1+\varepsilon_3)
$$
The relative error is now small if $|x|$ and $|y|$ are not much bigger
than $|x+y|$. Again ignoring higher order terms, we can bound the
relative error of $z$ by
$$
      \frac{|x|+|y|}{|x+y|}\,{\bf u} + {\bf u}
$$
Thus, the addition increases relative errors by a factor of
$$
      \frac{|x|+|y|}{|x+y|}+1
$$
Similar estimates can be made in the case of subtraction $x-y$:
it increases relative errors by a factor
$$
      \frac{|x|+|y|}{|x-y|}+1
$$
We see that the addition and subtraction increase relative errors by a
variable factor which depends on $x$ and $y$. This factor may be
arbitrarily large if $x+y\approx 0$ in addition or $x-y\approx 0$ in
subtraction. This phenomenon is known as {\em catastrophic
cancellation}. It occurred in our Example~\ref{14.11}, when we
attempted to solve the system without pivoting.


\newpage

\section{Conditioning}
\label{15}

\subsection{Condition number of a function} \label{14.15}

Let $V$ and $W$ be two normed vector spaces, and $f\colon V\to W$ a
function. The {\em condition number} $\kappa$ of $f$ at a point $x\in
V$ is
$$
   \kappa = \kappa(f,x) = \lim_{\delta\to 0}
   \sup_{\|\Delta x\|\leq\delta}
   \left ( \frac{\|\Delta f\|}{\|f\|}
   \bigg/ \frac{\|\Delta x\|}{\|x\|} \right )
$$
where $\Delta f = f(x+\Delta x) - f(x)$. This is the maximal factor by
which relative errors are magnified by $f$ in the vicinity of the point
$x$. The condition number characterizes the {\em sensitivity} of $f(x)$
to small perturbations of $x$.

\subsection{Lemma} \label{14.15a}

Let $\IC^n\to\IC^n$ be a linear operator defined by a nonsingular
matrix $A\in\IC^{n\times n}$, and let $\|\cdot\|$ be a norm on
$\IC^n$. Then for every $x\in\IC^n$
$$
   \kappa(A,x)\leq \|A\|\,\|A^{-1}\|,
   \quad{\rm and}\quad
   \sup_{x}\kappa(A,x) = \|A\|\,\|A^{-1}\|
$$
where $\|A\|$ denotes the induced matrix norm.

\proof Let $y=Ax$. Since $\Delta y=A(\Delta x)$, then
$$
  \kappa(A,x)=\sup_{\Delta x\neq 0}
  \frac{\|A(\Delta x)\|}{\|\Delta x\|}
  \frac{\|x\|}{\|Ax\|}
  =\|A\|\frac{\|x\|}{\|Ax\|}
$$
and $\sup_x \|x\|/\|Ax\|=\sup_y \|A^{-1}y\|/\|y\| = \|A^{-1}\|$.


\subsection{Condition number of a matrix} \label{15.2}

For a nonsingular matrix $A\in\IC^{n\times n}$, the condition
number with respect to a given matrix norm $\|\cdot \|$ is
$$
    \kappa(A)=\|A\|\,\|A^{-1}\|
$$
We denote by $\kappa_1(A)$, $\kappa_2(A)$, $\kappa_{\infty}(A)$ the
condition numbers with respect to the 1-norm, 2-norm and $\infty$-norm,
respectively.

\subsection{Theorem} \label{15.3}

Let $A\in\IC^{n\times n}$ be a nonsingular matrix and
\begin{align*}
     Ax &= b \\
    (A+\Delta A)(x+\Delta x)&=b+\Delta b
\end{align*}
Assume that $\|\Delta A\|$ is small so that $\|\Delta A\|\,
\|A^{-1}\|<1$. Then
$$
   \frac{\|\Delta x\|}{\|x\|}\leq
   \frac{\kappa(A)}{1-\kappa(A)\frac{\|\Delta A\|}{\|A\|}}
   \left (\frac{\|\Delta A\|}{\|A\|}+\frac{\|\Delta b\|}{\|b\|}\right )
$$

\proof Expanding out the second equation, subtracting the first one and
multiplying by $A^{-1}$ gives
$$
   \Delta x=-A^{-1}\Delta A(x+\Delta x)+A^{-1}\Delta b
$$
Taking norms and using the triangular inequality gives
$$
   \|\Delta x\|\leq \|A^{-1}\|\, \|\Delta A\|\Big (\|x\|+\|\Delta x\|\Big )
         +\|A^{-1}\|\, \|\Delta b\|
$$
Using $\|b\|\leq \|A\|\, \|x\|$, this rearranges to
$$
   \Big (1-\|A^{-1}\|\, \|\Delta A\|\Big )\|\Delta x\|\leq
    \left (\|A^{-1}\|\, \|\Delta A\|+\|A^{-1}\|\, \|A\|
      \frac{\|\Delta b\|}{\|b\|}\right )\|x\|
$$
Recall that $\|\Delta A\|\, \|A^{-1}\|<1$, so the first factor above is
positive. The theorem follows immediately. \qed

\new {\em Interpretation}. Let $Ax=b$ be a given system of
linear equations to be solved numerically. A computer represents $A$ by
$A_c=A+\Delta A$ and $b$ by $b_c=b+\Delta b$ and then solves the system
$A_cx=b_c$. Assume that the computer finds the exact solution
$x_c=x+\Delta x$, i.e., $A_cx_c=b_c$. Then the relative error $\|\Delta
x\|/\|x\|$ can be estimated by \ref{15.3}. The smaller the condition
number $\kappa(A)$, the tighter (better) estimate on $\|\Delta
x\|/\|x\|$ we get. The value of $\kappa(A)$ thus characterizes the {\em
sensitivity} of the solution of the linear system $Ax=b$ to small
errors in $A$ and $b$.

\subsection{Corollary}

Consider the problem of solving a system of linear equations $Ax=b$
with a nonsingular matrix $A$. Then:
\begin{itemize} \item[(a)] If we fix $A$ and vary $b$, we get a map
$f_A\colon b\mapsto x$. The condition number of $f_A$ satisfies
$$
   \kappa(f_A,b)\leq \kappa(A)
   \quad{\rm and}\quad
   \sup_{b\in\IC^n}\kappa(A,b) = \kappa(A)
$$
\item[(b)] If we fix $b$ and vary $A$, we get a map
$f_b\colon A\mapsto x$. The condition number of $f_b$ satisfies
$$
   \kappa(f_b,A) \leq \kappa(A)
$$
\end{itemize}

\proof Both inequalities immediately follow from Theorem~\ref{15.3}. To
prove the equality in (a), note that $x=A^{-1}b$ and apply
Lemma~\ref{14.15a}. \qed

\new {\em Remark}. In the part (b), we actually have equality
$\kappa(f_b,A)=\kappa(A)$, but the proof is beyond the scope of
our course, see the textbook.


\subsection{Corollary} \label{15.4}

Assume that in Theorem~\ref{15.3} we have $\|\Delta A\|\leq {\bf
u}\|A\|$ and $\|\Delta b\|\leq {\bf u}\|b\|$, i.e. the matrix $A$ and
the vector $b$ are represented with the best possible machine accuracy.
Then
$$
   \frac{\|\Delta x\|}{\|x\|}\leq
   \frac{2{\bf u}\kappa(A)}{1-{\bf u}\kappa(A)}
$$

\subsection{Remark} \label{15.5}

Assume that ${\bf u}\approx 10^{-l}$, i.e. the machine system provides $l$
accurate digits. Then if $\kappa(A)\approx 10^k$ with $k<l$, then
$\|\Delta x\|/\|x\|\leq 10^{-(l-k)}$, i.e. even a numerical solution only
provides $l-k$ accurate digits.

\new Linear systems $Ax=b$ with small $\kappa(A)$ ($\sim 1$, $10$, $10^2$) are often
called {\em well-conditioned}. Those with large $\kappa(A)$ ($\sim
10^3$, $10^4$, etc.) are called {\em ill-conditioned}. Their numerical
solutions are unreliable and should be avoided.

\subsection{Proposition} \label{15.6}

\begin{enumerate}
\item We have
$$
       \kappa(A)=\frac{\sup_{\|x\|=1}\|Ax\|}{\inf_{\|x\|=1}\|Ax\|}
$$
\item If $a_j$ denotes the $j$-th column of $A$, then $\kappa(A)\geq \|a_j\|/\|a_i\|$
\item $\kappa(A)\geq 1$ and $\kappa(I)=1$
\item $\kappa_2(A)=1$ if and only if $A$ is a multiple of an unitary matrix.
\item For any unitary matrix $Q$,
   $$
        \kappa_2(QA)=\kappa_2(AQ)=\kappa_2(A)
   $$
\item If $D=\, {\rm diag}\{d_1,\dots,d_n\}$ then
   $$
       \kappa_2(D)=\kappa_1(D)=\kappa_{\infty}(D)=
       \frac{\max_{1\leq i\leq n}|d_i|}{\min_{1\leq i\leq n}|d_i|}
   $$
\item If $\sigma_{\max}$ is the largest singular value of $A$ and
   $\sigma_{\min}$ is its smallest singular value, then
   $$
       \kappa_2(A) = \sigma_{\max}/\sigma_{\min}
   $$
\item If $\lambda_{\max}(B)$ is the largest eigenvalue of $B=A^\ast A$ and
   $\lambda_{\min}(B)$ is its smallest eigenvalue, then
   $$
       \kappa_2(A^\ast A)=\kappa_2(AA^\ast)=
       [\kappa_2(A)]^2=\lambda_{\max}(B)/\lambda_{\min}(B)
   $$
   If $A$ is Hermitian with eigenvalues $\lambda_1,\ldots,\lambda_n$, then
   $$
      \kappa_2(A)=\frac{\max_i|\lambda_i|}{\min_i|\lambda_i|}
   $$
\item $\kappa_2(A)=\kappa_2(A^T)$.
\end{enumerate}


\subsection{Remark}\label{15.5a}

Another way to look at the condition number $\kappa(A)$ is the
following:
$$
   \min\left\{\frac{\|A-A_s\|_2}{\|A\|_2}\colon\
   A_s\ {\rm is}\ {\rm singular}\right \}=
    \frac{1}{\kappa_2(A)}
$$
hence $1/\kappa_2(A)$ is the relative distance from $A$ to the
nearest singular matrix. This follows from \ref{17.20}.

\subsection{Remark} \label{15.8}

Yet another way to look at the condition number $\kappa(A)$ is the
following. Since in practice the exact solution $x$ of the system
$Ax=b$ is rarely known, one can try to verify whether the so
called {\em residual vector} $r=b-Ax_c$ is small relative to $b$.
Since $Ax_c=b+r$, Theorem~\ref{15.3} with $\Delta A=0$ implies
that
$$
   \frac{\|x_c-x\|}{\|x\|}\leq \kappa(A)\frac{\|r\|}{\|b\|}
$$
If $A$ is well conditioned, the smallness of $\|r\|/\|b\|$ ensures
the smallness of the relative error $\|x_c-x\|/\|x\|$. If $A$ is
ill-conditioned, this does not work: the smallness of $r$ does not
guarantee the accuracy of $x_c$.

\subsection{Definition} \label{15.6a}

The condition number of a rectangular $m\times n$ matrix $A$ with $m>n$
is defined by
$$
    \kappa(A) \colon =
    \frac{\sup_{\|x\|=1}\|Ax\|}
    {\inf_{\|x\|=1}\|Ax\|},
$$
motivated by \ref{15.6} (1). Under this definition, the properties
2, 3, and 7 in \ref{15.6} will still hold, and in the property 8 we
only claim that
$$
       \kappa_2(A^T A)=
       [\kappa_2(A)]^2=\lambda_{\max}(B)/\lambda_{\min}(B)
$$
This definition gives a clear geometric interpretation of
$\kappa_2(A)$: it is the maximum distortion factor of the map
$A\colon \IC^n \to \IC^m$.

\newpage

\section{Stability}

\subsection{Round-off error analysis} \label{15.9}

When a system of linear equations $Ax=b$ is solved numerically,
the computer gets the machine representations $A_c$ and $b_c$ of
$A$ and $b$. Then, at best, it can find $x_c=A_c^{-1}b_c$, instead
of $x=A^{-1}b$, and we know that the resulting error will be
$$
     E_1\colon= \frac{\|x_c-x\|}{\|x\|}\lesssim 2{\bf u}\kappa(A)
$$
by \ref{15.4} (we ignored the small term ${\bf u}\kappa(A)$ in the
denominator). This may be bad enough already, when the problem is
ill-conditioned. However, in reality things appear to be even worse,
since the computer does not evaluate $x_c$ precisely, apart from
trivial cases. The computer executes a certain program that normally
involves many arithmetic operations. As the program runs, more and more
round-off errors are made at each step and the errors compound toward
the end. As a result, the computer finds a vector $\hat{x}_c\neq x_c$,
so that $A_c\hat{x}_c\neq b_c$. The actual output $\hat{x}_c$ depends
not only on the machine system but even more on the algorithm that is
used to evaluate it. See the diagram nearby.

Of course, we do not expect the final error
$$
      E_2\colon= \frac{\|\hat{x}_c-x\|}{\|x\|}
$$
to be smaller than $E_1$, but we hope that it will not be much
larger either. In other words, a good algorithm should not magnify
the errors caused already by conditioning. If this is the case,
the algorithm is said to be {\em stable}.

\subsection{Stable algorithms}

An algorithm for solving a system of linear equations $Ax=b$ is {\em
stable} if
$$
      \frac{\|\hat{x}_c-x\|}{\|x\|} \leq C{\bf u}\kappa(A)
$$
where $C>0$ is independent of $A,b$ and the machine system (but it
may depend on the size of the matrix, $n$).

\subsection{Backward error analysis} \label{15.10}

In order to estimate the final error $E_2$, a typical approach is
to ``trace the errors backwards'' and find another matrix,
$\hat{A}=A+\delta A$ such that $(A+\delta A)\hat{x}_c=b_c$. We
call $A+\delta A$ a {\em virtual} matrix, since it is neither
given nor computed numerically. Moreover, it is far from being
unique. One wants to find a virtual matrix as close to $A$ as
possible, to make $\delta A$ small, for the reasons made clear
below.

\setlength{\unitlength}{0.02in}

\begin{figure}
\begin{picture}(300,100)(0,0)

\put(80,80){\vector(1,0){120}} \put(80,50){\vector(1,0){120}}
\put(80,50){\vector(4,-1){120}} \put(80,20){\line(1,0){10}}
\put(100,20){\line(1,0){10}} \put(120,20){\line(1,0){10}}
\put(140,20){\line(1,0){10}} \put(160,20){\line(1,0){10}}
\put(180,20){\vector(1,0){10}}

\put(77,80){\circle*{2}} \put(77,50){\circle*{2}}
\put(77,20){\circle{2}} \put(203,80){\circle*{2}}
\put(203,50){\circle*{2}} \put(203,20){\circle*{2}}

\put(45,83){exact $A,b$} \put(37,53){machine $A_c,b_c$}
\put(37,23){`virtual' $\hat{A}_c,b_c$} \put(199,83){exact $x$}
\put(199,53){`idealistic' $x_c$} \put(199,23){computed $\hat{x}_c$}
\put(115,83){exact solution} \put(115,53){exact solution}
\put(115,13){exact solution} \put(136,36){computer}
\put(164,30){algorithm}

\end{picture}
%\caption{The scheme of round-off error analysis.}
\end{figure}

\subsection{Backward stable algorithms}

An algorithm for solving a system of linear equations $Ax=b$ is {\em
backward stable} if there exists a virtual matrix $A+\delta A$ such
that
$$
      \frac{\|\delta A\|}{\|A\|} \leq C{\bf u}
$$
where $C>0$ is independent of $A,b$ and the machine system (but it
may depend on the size of the matrix, $n$).

\subsection{Theorem}

A backward stable algorithm is stable.

\proof By Theorem~\ref{15.3},
$$
   \frac{\|\hat{x}_c-x\|}{\|x\|}\leq
   \frac{\kappa(A)}{1-\kappa(A)\frac{\|\delta A\|}{\|A\|}}
   \frac{\|\delta A\|}{\|A\|} \lesssim C{\bf u}\kappa(A).
   \qquad\square
$$

The proofs of stability (or instability) of algorithms of linear
algebra are quite involved. We only present relevant facts here,
without proofs.

\subsection{Theorem}
\label{15.11}

If one uses the LU decomposition $A=LU$ for solving a system
$Ax=b$, then there is a virtual matrix $A+\delta A$ such that
$$
  \|\delta A\| \leq C \|L\|\, \|U\| {\bf u},
$$
where $C>0$ is independent of $A$ and the machine system (but it
may depend on the size of the matrix, $n$). Thus, the LU algorithm
is unstable, its accuracy deteriorates when $\|L\|\,\|U\|$ is
large.

\subsection{Example} \label{15.11a}

In example \ref{14.11}, the LU decomposition (without pivoting) is
$$
  \left[\begin{array}{cc} 0.01 & 2 \\ 1 & 3 \end{array}\right]
 =\left[\begin{array}{cc} 1 & 0 \\ 100 & 1 \end{array}\right]\,
  \left[\begin{array}{cc} 0.01 & 2 \\ 0 & -197 \end{array}\right]
$$
hence $\|L\|\,\|U\|\sim 10^4$. This explains the huge errors of the
corresponding numerical solutions that we observed.

\subsection{Remarks}

\begin{itemize}\item[(a)]
Applying partial pivoting ensures that the entries of $L$ are uniformly
bounded: $|L_{ij}|\leq 1$. Also, it is observed in practice that in
most cases $\|U\|\leq C\|A\|$, hence the partial pivoting algorithm is
usually stable. \item[(b)] The LU decomposition with complete pivoting
is always stable, in this case one can prove that $\|U\|\leq C\|A\|$.
\item[(c)] The Cholesky factorization $A=GG^{\ast}$ of a positive
definite matrix $A$ is a particular form of the LU decomposition, so
the above analysis applies. In this case, we know that
$$
       a_{ii}=\sum_{j=1}^i g_{ij}^2
$$
see \ref{13.6}. Thus, one can easily prove that $\|G\|\leq
C\|A\|^{1/2}$, hence the Cholesky factorization is always stable.
\end{itemize}

\subsection{Example} \label{15.11b}

In Example \ref{14.11}, the matrix
$$
  \begin{bmatrix}0.01 & 2 \\ 1 & 3\end{bmatrix}
$$
has singular values $\sigma_1=3.7037$ and $\sigma_2=0.5319$, hence
its condition number is $\kappa(A)= \sigma_1/\sigma_2 = 6.96$.
This explains a moderate factor ($\leq 5$) by which relative
errors of the numerical solutions are related to the minimal error
$10^{-t}$ in Example~\ref{14.11}.

In Example \ref{14.11a}, the matrix
$$
  \begin{bmatrix}3 & 1 \\ 1 & 0.35\end{bmatrix}
$$
has singular values $\sigma_1=3.33$ and $\sigma_2=0.0150$, hence
its condition number is $\kappa(A)=222$. This explains a large
factor (up to 200) by which relative errors of the numerical
solutions are related to the minimal error $10^{-t}$ in
Example~\ref{14.11a}, even though we used a stable algorithm (the
LU with complete pivoting). Remember that a stable algorithm {\em
should not increase} errors already caused by conditioning, but it
cannot cancel them out.

\newpage

\section{Numerical solution of overdetermined systems}
\label{16a}

In Chapter~\ref{16} we discussed least squares solution of
overdetermined systems of equations
$$
        Ax = b
$$
where $A\in\IC^{m\times n}$ with $m\geq n$ and developed three
algorithms, based on normal equations, QR decomposition, and SVD
decomposition.

\subsection{Comparison of algorithms} \label{16.10}

\begin{itemize}
\item[(a)]
Algorithm \ref{A1} via normal equations  has many advantages, it is the
most compact and elegant one, and it is twice as cheap as the other
two. It should be used whenever the computations are precise. However,
if computations involve round-off errors, other considerations come
into play. If the matrix $A$ is ill-conditioned (i.e.\ its condition
number is large, $\kappa_2(A)\gg 1$), then by \ref{15.6a}
$$
      \kappa_2(A^{\ast}A)=[\kappa_2(A)]^2
$$
hence the condition of the matrix $A^\ast A$ will be much worse than
that of $A$, and solving the normal equations can be disasterous. For
example, let
$$
    A=\left[\begin{array}{cc} 1 & 1 \\ \varepsilon & 0 \\
     0 & \varepsilon \end{array}\right]
$$
then
$$
   A^\ast A=\left[\begin{array}{cc} 1+\varepsilon^2 & 1 \\
        1 & 1+\varepsilon^2 \end{array}\right]
$$
If $\varepsilon$ is so small that $\varepsilon^2<{\bf u}$ (for
example, $\varepsilon=10^{-4}$ in single precision), then the
matrix $A^\ast A$ will be stored in computer memory as a singular
matrix, and the algorithm via normal equations will fail. Still,
it is possible to find a good numerical solution to the original
system $Ax=b$, if one uses more elaborate methods. \item[(b)] Due
to \ref{15.6}, $\kappa_2(A)=\kappa_2(QA)$ for any unitary matrix
$Q\in\IC^{m\times m}$, hence $\kappa_2(A)=\kappa_2(R)$ in the
QR-based algorithm \ref{A2} and $\kappa_2(A)=\kappa_2(D)$ in the
SVD-based algorithm \ref{A3}. Hence, these methods are safe,
regarding the conditioning of the problem. The other aspect of
numerical algorithms is stability. It turns out that the
Gram-Schmidt orthogonalization we used in Chapter~\ref{QR} to
construct the QR decomposition is {\em unstable}, hence it leads
to unpredictable round-off errors in numerical computations. There
is a modified Gram-Schmidt algorithm (see programming assignment),
which is stable, but still far from optimal. Very good algorithms
for constructing the QR decomposition are based on reflection
matrices or rotation matrices, see below. \item[(c)] The SVD-based
algorithm \ref{A3} requires the computation of the SVD of the
matrix $A$, for which no simple algorithm exists, and such
algorithms are beyond the scope of this course. Standard software
packages (like MATLAB) use only stable algorithms for the SVD
computation. Practically, the SVD-based method \ref{A3} is as good
as the QR-based method \ref{A2} (see programming assignment). It
is observed that the SVD-based method is more reliable when the
matrix $A$ is nearly singular or just singular, cf.\ \ref{Adefic}.
\end{itemize}

\subsection{Hyperplanes and reflections} \label{16.11}

Let $V$ be a finite dimensional vector space. A subspace $W\subset V$
is called a {\em hyperplane} if dim$\, W=\,$dim$\,V-1$. Note that in
this case dim$\, W^{\perp}=1$.

Let $W\subset V$ be a hyperplane. For any vector $v\in V$ we have a
unique decomposition $v=w+w'$, where $w\in W$ and $w'\in W^{\perp}$.
The linear operator $P$ on $V$ defined by $Pv=w-w'$ is called a {\em
reflection} (or {\em reflector}) across the hyperplane $W$. It is
identical on $W$ and negates vectors orthogonal to $W$.

\subsection{Householder reflector matrices} \label{16.12}

Let $x\neq 0$ be a vector in $\IR^n$ or $\IC^n$. The matrix
$$
     P=I-2\,\frac{x x^\ast}{x^\ast x}
     =I-2\,\frac{x x^\ast}{\|x\|^2}
$$
is called the {\em Householder reflector matrix} corresponding to $x$.
Obviously, $P$ is unchanged if $x$ is replaced by $cx$ for any $c\neq
0$.

\subsection{Theorem} \label{16.13}

Let $P$ be a reflection matrix corresponding to a vector $x\neq 0$.
Then \begin{itemize}\item[(a)] $Px=-x$. \item[(b)] $Py=y$ whenever $\la
y,x\ra =0$. \item[(c)] $P$ is Hermitian (in the real case it is
symmetric). \item[(d)] $P$ is unitary (in the real case it is
orthogonal). \item[(e)] $P$ is involution, i.e. $P^2=I$.
\end{itemize}

\proof Direct calculation.

\subsection{Theorem} \label{16.14}

Let $y$ be a vector in $\IR^n$ or $\IC^n$. Choose a scalar $\sigma$ so
that $|\sigma|=\|y\|$ and $\sigma\cdot\la e_1,y\ra\in\IR$. Suppose that
$x=y+\sigma e_1\neq 0$. Let $P=I-2x x^\ast/\|x\|^2$ be the reflection
matrix defined in \ref{16.12}. Then $Py=-\sigma e_1$.

\proof Notice that $\la y-\sigma e_1,y+\sigma e_1\ra =\|y\|^2-\sigma\la
e_1,y\ra+\bar{\sigma}\la y,e_1\ra-|\sigma|^2=0$. Hence, $P(y-\sigma
e_1)=y-\sigma e_1$ by \ref{16.13} (b). Besides, $P(y+\sigma
e_1)=-y-\sigma e_1$ by \ref{16.13} (a). Adding these two equations
proves the theorem.

\subsection{Remark} \label{16.15}

\begin{itemize}\item[(a)]
To choose $\sigma$ in \ref{16.14}, write a polar representation for
$\la e_1,y\ra=\bar{y}_1=re^{i\theta}$ and then set $\sigma=\pm
\|y\|e^{-i\theta}$. \item[(b)] In the real case, we have $y_1\in\IR$,
and one can just set $\sigma=\pm \|y\|$.
\end{itemize}

\new Note: It is geometrically obvious that for any two
unit vectors $x,y\in\IR^n$ there is a reflection $P$ that takes $x$ to
$y$. In the complex space $\IC^n$, this is not true: for generic unit
vectors $x,y$ one can only find a reflection that takes $x$ to $cy$
with some scalar $c\in\IC$.

\subsection{Corollary} \label{16.16}

For any vector $y$ in $\IR^n$ or $\IC^n$ there is a scalar $\sigma$
(defined in \ref{16.14}) and a matrix $P$, which is either a reflection
or an identity, such that $Py=-\sigma e_1$.

\proof Apply Theorem~\ref{16.14} in the case $y+\sigma e_1\neq 0$ and
set $P=I$ otherwise.

\subsection{QR Decomposition via Householder reflectors} \label{16.17}

For any $A\in\IC^{m\times n}$ with $m\geq n$ there is a QR
decomposition with a unitary matrix $Q\in\IC^{m\times m}$ that is a
product of at most $n$ reflector matrices.

\proof We use induction on $n$. Let $n=1$, so that $A$ is a column
$m$-vector. By \ref{16.16} there is a matrix $P$ (a reflection or an
identity) such that $PA=-\sigma e_1$ for a scalar $\sigma$. Hence,
$A=PR$ where $R=-\sigma e_1$ is upper triangular. Now, let $n\geq 1$
and $a_1$ the first column of $A$. Again, by \ref{16.16} there is a
(reflection or identity) matrix $P$ such that $Pa_1=-\sigma e_1$.
Hence,
$$
    PA=\left[\begin{array}{rc} -\sigma &
     w^\ast \\ 0 & B\end{array}\right]
$$
where $w\in\IC^{n-1}$ and $B\in\IC^{(m-1)\times (n-1)}$. By the
inductive assumption, there is a unitary matrix $Q_1\in\IC^{(m-1)\times
(m-1)}$ and an upper triangular matrix $R_1\in\IC^{(m-1)\times (n-1)}$
such that $B=Q_1R_1$. Consider the unitary $m\times m$ matrix
$$
     Q_2=\left[\begin{array}{cc} 1 &
     0 \\ 0 & Q_1 \end{array}\right]
$$
By \ref{10.6}, the matrix $Q_2$ is unitary whenever $Q_1$ is.
Furthermore, if $Q_1$ is a product of $\leq n-1$ reflectors, then the
same is true for $Q_2$. Now one can easily check that $PA=Q_2R$ where
$$
     R=\left[\begin{array}{rc} -\sigma & w^\ast \\
      0 & R_1 \end{array}\right]
$$
is an upper triangular matrix. Hence, $A=QR$ with $Q=PQ_2$.

\subsection{Remark}

In the real case, there are two choices for the scalar $\sigma$, see
\ref{16.15}, that is $\sigma=\pm \|y\|$. The better one is
$$
     \sigma = \,{\rm sgn}(y_1)\,\|y\|
$$
i.e.\ the sign of $\sigma$ is determined by the sign of $y_1$. Then
computing the vector $x=y+\sigma e_1$ is always stable, there is no
danger of catastrophic cancellation.

\subsection{Givens rotation matrices} \label{16.24}

Let $1\leq p<q\leq m$ and $\theta\in [0,2\pi)$. The matrix
$G=G_{p,q,\theta}=(g_{ij})$ defined by $g_{pp}=\cos\theta$,
$g_{pq}=\sin\theta$, $g_{qp}=-\sin\theta$, $g_{qq}=\cos\theta$ and
$g_{ij}=\delta_{ij}$ otherwise is called a {\em Givens rotation matrix}
(or a {\em Givens rotator}). It defines a rotation through the angle
$\theta$ of the $x_px_q$ coordinate plane in $\IR^m$ with all the other
coordinates fixed. Obviously, $G$ is orthogonal.

\subsection{QR decomposition via Givens rotators} \label{16.25}

For any $A\in\IR^{m\times n}$ with $m\geq n$ there is a QR
decomposition with an orthogonal matrix $Q\in\IR^{m\times m}$ that is a
product of Givens rotator matrices.

\proof Let $a_j$ be the leftmost column of $A$ that contains a nonzero
entry below the main diagonal, $a_{ij}\neq 0$ with some $i>j$. Consider
the matrix $A'=GA$ where $G=G_{j,i,\theta}$ is a Givens rotation
matrix. One easily checks that\\ (a) the first $j-1$ columns of $A'$
are zero below the main diagonal;\\ (b) in the $j$-th column, only the
elements $a_{jj}'$ and $a_{ij}'$ will be different from the
corresponding elements of $A$, and moreover
$$
   a_{ij}'=-a_{jj}\sin\theta+a_{ij}\cos\theta
$$
Now we find $\sin\theta$ and $\cos\theta$
so that $a_{ij}'=0$. For example, let
$$
    \cos\theta=\frac{a_{jj}}{\sqrt{a_{jj}^2+a_{ij}^2}}
    \qquad\text{and}\qquad
    \sin\theta=\frac{a_{ij}}{\sqrt{a_{jj}^2+a_{ij}^2}}
$$
Note that one never actually evaluates the angle $\theta$, since
$G_{j,i,\theta}$ only contains $\cos\theta$ and $\sin\theta$.

In this way we eliminate one nonzero element $a_{ij}$ below the main
diagonal. Working from left to right, one can convert $A$ into an upper
triangular matrix $\tilde{G}A=R$ where $\tilde{G}$ is a product of
Givens rotators. Each nonzero element of $A$ below the main diagonal
requires one multiplication by a rotator. Then we get $A=QR$ with an
orthogonal matrix $Q=\tilde{G}^T$. \qed

\subsection{Cost of QR via Givens rotators} \label{16.25a}

The evaluation of $\cos\theta$ and $\sin\theta$ takes 6 flops (the
square root extraction is regarded as one flop), then the subsequent
multiplication of $A$ by $G_{j,i,\theta}$ takes $6n$ flops. Thus, if
$A$ originally had $p$ nonzero subdiagonal entries, then the QR
decomposition via Givens rotators takes $6pn$ flops.

When $p$ is close to its maximal value, $mn - n^2/2$, then the total
cost $\sim 6mn^2 - n^3/2$ greatly exceeds the cost of QR via
Householder reflectors or Gram-Schmidt decomposition. Hence Givens
rotators are very inefficient for generic matrices. But they work well
if the matrix $A$ is {\em sparse}, i.e. contains just a few nonzero
elements below the main diagonal. Then Givens rotators give the
quickest result. We will see such instances later.

\newpage

\section{Computation of eigenvalues: theory}
\label{18}

\subsection{Preface} \label{19.1}

Eigenvalues of a matrix $A\in\IC^{n\times n}$ are the roots of its
characteristic polynomial, $C_A(x)$. It is a consequence of the
famous Galois group theory (Abel's theorem) that there is no
formulas for the roots of a generic polynomial of degree $>4$.
Hence, there are no finite algorithms for computation of the roots
of polynomials (or eigenvalues, for that matter).

Thus, all the methods for computing eigenvalues of matrices of
size $n\geq 5$ are necessarily iterative, they provide successive
approximations to the eigenvalues, but never exact results.
Furthermore, even though for $n=3$ and $n=4$ exact formulas exist,
they are rather impractical and often numerically unstable, so
even in these cases iterative methods should be used instead.

Furthermore, matrix decompositions (such as Schur and SVD) that
involve eigenvalues cannot be implemented by finite algorithms.
Note that the SVD requires the eigenvalues of a Hermitian positive
semidefinite matrix $A^\ast A$, hence it is particularly important
to develop algorithms for this class of matrices.

If an eigenvalue $\lambda$ of a matrix $A$ is known, an
eigenvector $x$ can be found by solving the linear system
$(A-\lambda I)x=0$ (say, by LU decomposition). Conversely, if an
eigenvector $x$ is known, the eigenvalue $\lambda$ can be
immediately found by $(Ax)_i/x_i$ whenever $x_i\neq 0$. Hence,
eigenvalues and eigenvectors are often computed `in parallel'. In
this section, we develop a theoretical basis for computation of
eigenvalues and eigenvectors, while in the next section we
consider practical algorithms.


\subsection{Rayleigh quotient} \label{18.4}

Let $A\in\IC^{n\times n}$. We call
$$
      r(x)= \frac{x^{\ast}Ax}{x^{\ast}x}
      =\frac{\la Ax,x\ra}{\la x,x\ra},
      \qquad x\neq 0
$$
the {\em Rayleigh quotient} of $A$.

\new Note: $r(cx)=r(x)$ for $c\neq 0$,
hence $r(x)$ is constant on the line span$\{x\}$ (with the zero
vector removed). Since any nonzero vector is a scalar multiple of
a unit vector, $r(x)$ is completely defined by its values on the
unit sphere ${\mathbb S}_1$, on which
$$
          r(x)= x^{\ast}Ax=\la Ax,x\ra,
          \qquad \|x\|_2=1
$$
thus $r(x)$ is a quadratic function of the coordinates of $x$.

\new If $A$ is Hermitian, then $r(x)\in\IR$ for any nonzero
$x\in\IC^n$.

\new If $x$ is a unit eigenvector with an eigenvalue $\lambda$, then
$Ax=\lambda x$ and so
$$
       r(x)=\lambda
$$
If $x$ is an arbitrary unit vector, then $r(x)\,x=(x^{\ast}Ax)\,x$
is the orthogonal projection of the vector $Ax$ on the line
spanned by $x$. Hence
$$
   \|Ax-r(x)\,x\|_2=\min_{\mu\in\IC}\|Ax-\mu x\|_2
$$
If one regards $x$ as an `approximate' eigenvector, then the
Rayleigh quotient $r(x)$ is the best choice that one could make
for the associated `approximate' eigenvalue in the sense that this
value of $\mu$ comes closest (in the 2-norm) to achieving $Ax-\mu
x=0$.

\subsection{Theorem} \label{18.6}

Let $A\in\IC^{n\times n}$ and $x$ a unit eigenvector of $A$
corresponding to eigenvalue $\lambda$. Let $y$ be another unit
vector and $r=y^{\ast}Ay$. Then
$$
   |\lambda - r|\leq 2\, \|A\|_2\, \|x-y\|_2
$$
Moreover, if $A$ is a Hermitian matrix, then there is a constant
$C=C(A)>0$ such that
$$
   |\lambda - r|\leq C\, \|x-y\|_2^2
$$

\proof To prove the first part, put
$$
   \lambda-r=x^{\ast}A(x-y)+(x-y)^{\ast}Ay
$$
and then using the triangle inequality and Cauchy-Schwarz inequality
gives the result. Now, assume that $A$ is Hermitian. Then there is an
ONB of eigenvectors, and we can assume that $x$ is one of them. Denote
that ONB by $\{x,x_2,\ldots,x_n\}$ and the corresponding eigenvalues by
$\lambda,\lambda_2,\ldots, \lambda_n$. Let $y=cx+c_2x_2+\cdots
+c_nx_n$. Then
$$
     \|y-x\|^2=|c-1|^2+\sum_{i=2}^n |c_i|^2\geq \sum_{i=2}^n |c_i|^2
$$
On the other hand, $\|y\|=1$, so
$$
      \lambda=\lambda |c|^2+\sum_{i=2}^n\lambda |c_i|^2
$$
Now, $Ay=c\lambda x+\sum_{i=2}^nc_i\lambda_ix_i$, so
$$
     r=\la Ay,y\ra =\lambda |c|^2+\sum_{i=2}^n\lambda_i |c_i|^2
$$
Therefore,
$$
     \lambda-r=\sum_{i=2}^n(\lambda-\lambda_i)|c_i|^2
$$
The result now follows with
$$
      C=\max_{2\leq i\leq n}|\lambda-\lambda_i|
$$
The theorem is proved. \qed

\subsection{Lemma} \label{18.7}

Let $A\in\IC^{n\times n}$ be Hermitian with eigenvalues
$\lambda_1\leq\cdots\leq\lambda_n$. Then for any $\|x\|=1$
$$
       \lambda_1\leq x^{\ast}Ax\leq\lambda_n
$$

\proof Let $\{x_1,\ldots,x_n\}$ be an ONB of eigenvectors of $A$.
Let $x=c_1x_1+\cdots+c_nx_n$. Then
$x^{\ast}Ax=\lambda_1|c_1|^2+\cdots +\lambda_n|c_n|^2$. The result
now follows easily.

\subsection{Lemma} \label{18.8}

Let $L$ and $G$ be subspaces of $\IC^n$ and dim$\, G >\,$dim$\, L$.
Then there is a nonzero vector in $G$ orthogonal to $L$.

\proof By way of contradiction, if $G\cap L^{\perp}=\{0\}$, then
$G\oplus L^{\perp}$ is a subspace of $\IC^n$ with dimension  dim$\,
G+n-\,$dim$\, L$. But our assumptions imply that this value exceeds
$n$, a contradiction.


\subsection{Courant-Fisher Minimax Theorem}
\label{18.9}

Let $A\in\IC^{n\times n}$ be Hermitian with eigenvalues
$\lambda_1\leq\cdots\leq\lambda_n$. Then for every $i=1,\ldots,n$
$$
     \lambda_i=\min_{ {\rm dim}\, L=i}\ \max_{x\in L\setminus \{0\} }
       \frac{x^{\ast}Ax}{x^{\ast}x}
$$
where $L$ stands for a subspace of $\IC^n$.

\proof Let $\{u_1,\ldots,u_n\}$ be an ONB of eigenvectors of $A$
corresponding to the eigenvalues $\lambda_1,\ldots,\lambda_n$. If dim$\,
L=i$, then by Lemma~\ref{18.8} there is a nonzero vector $x\in L$
orthogonal to the space span$\{u_1,\ldots,u_{i-1}\}$. Hence, the first
$i-1$ coordinates of $x$ are zero and $x=\sum_{j=i}^nc_ju_j$, so
$$
       \frac{x^{\ast}Ax}{x^{\ast}x}=
       \frac{\sum_{j=i}^n|c_j|^2\lambda_j}{\sum_{j=i}^n|c_j|^2}\geq\lambda_i
$$
Therefore,
$$
     \max_{x\in L\setminus \{0\} } \frac{x^{\ast}Ax}{x^{\ast}x}\geq\lambda_i
$$
Now, take the subspace $L=\,$span$\{u_1,\ldots,u_i\}$. Obviously,
dim$\, L=i$ and for every nonzero vector $x\in L$ we have
$x=\sum_{j=1}^ic_ju_j$, so
$$
       \frac{x^{\ast}Ax}{x^{\ast}x}=
       \frac{\sum_{j=1}^i|c_j|^2\lambda_j}{\sum_{j=1}^i|c_j|^2}\leq\lambda_i
$$
The theorem is proved.

\subsection{Theorem} \label{18.10}

Let $A$ and $\Delta A$ be Hermitian matrices. Let $\alpha_1\leq
\cdots\leq\alpha_n$ be the eigenvalues of $A$, $\delta_{\min}$ and
$\delta_{\max}$ the smallest and the largest eigenvalues of $\Delta A$.
Denote the eigenvalues of the matrix $B=A+\Delta A$ by
$\beta_1\leq\cdots\leq\beta_n$. Then for each $i=1,\ldots,n$
$$
      \alpha_i+\delta_{\min}\leq\beta_i\leq\alpha_i+\delta_{\max}
$$

\proof Let $\{u_1,\ldots,u_n\}$ be an ONB of eigenvectors of $A$
corresponding to the eigenvalues $\alpha_1, \ldots,\alpha_n$. Let
$L=\,$span$\{u_1,\ldots,u_i\}$. Then, by \ref{18.9},
\begin{align*}
  \beta_i &\leq \max_{x\in L\setminus\{0\}}\frac{x^{\ast}Bx}{x^{\ast}x}\\
          &\leq \max_{x\in L\setminus\{0\}}\frac{x^{\ast}Ax}{x^{\ast}x}
           + \max_{x\in L\setminus\{0\}}\frac{x^{\ast}\Delta
           Ax}{x^{\ast}x}\\
     &\leq \alpha_i+\max_{x\in C^n\setminus\{0\}}\frac{x^{\ast}\Delta
     Ax}{x^{\ast}x}\\
         &= \alpha_i+\delta_{\max}
\end{align*}
which is the right inequality. Now apply this theorem to the matrices
$B$, $-\Delta A$ and $A=B+(-\Delta A)$. Then its right inequality, just
proved, will read $\beta_i\leq\alpha_i-\delta_{\min}$ (note that the
largest eigenvalue of $-\Delta A$ is $-\delta_{\min}$). The theorem is
completely proved. \qed

\subsection{Corollary} \label{18.10a}

Since $\|\Delta A\|_2 = \max\{|\delta_{\min}|,|\delta_{\max}|\}$,
we have
$$
      \alpha_i-\|\Delta A|_2\leq\beta_i\leq
      \alpha_i+\|\Delta A\|_2 \qquad \forall i=1,\ldots,n
$$

\subsection{Remark} \label{18.11}

Suppose one knows an approximate eigenvalue $\lambda$ and an
approximate eigenvector $x$ of a matrix $A$. We can assume that
$\|x\|=1$. To estimate the closeness of $\lambda$ to the actual
but unknown eigenvalue of $A$, one can compute the residual
$r=Ax-\lambda x$. Assume that $r$ is small and define the matrix
$\Delta A=-rx^{\ast}$. Then $\|\Delta A\|_2=\|r\|_2$ and
$$
     (A+\Delta A)x=Ax-rx^{\ast}x=\lambda x
$$
Therefore, $(\lambda,x)$ are an exact eigenpair of a perturbed
matrix $A+\Delta A$, and the norm $\|\Delta A\|_2$ is known. One
could then apply Corollary~\ref{18.10a} to estimate the closeness
of $\lambda$ to the actual eigenvalue of $A$, if the matrices $A$
and $\Delta A$ were Hermitian. Since this is not always the case,
we need to study how eigenvalues of a generic matrix change under
small perturbations of the matrix. This is the issue of {\em
eigenvalue sensitivity}.

\subsection{Bauer-Fike theorem} \label{18.12}

Let $A\in\IC^{n\times n}$ and suppose that
$$
      X^{-1}AX=D=\, {\rm diag}\{\lambda_1,\ldots,\lambda_n\}
$$
If $\mu$ is an eigenvalue of a perturbed matrix $A+\Delta A$,
then
$$
    \min_{1\leq i\leq n}|\lambda_i-\mu|
       \leq \kappa_p(X)\,\|\Delta A\|_p
$$
where $\|\cdot \|_p$ stands for any $p$-norm.

\proof If $\mu$ is an eigenvalue of $A$, the claim is trivial. If not,
the matrix $D-\mu I$ is invertible. Observe that
\begin{align*}
  X^{-1}(A+\Delta A-\mu I)X &= D+X^{-1}\Delta A X-\mu I \\
             &= (D-\mu I)\bigl [ I+(D-\mu I)^{-1}
             (X^{-1}\Delta A X) \bigr ]
\end{align*}
Since the matrix $A+\Delta A-\mu I$ is singular, so is the matrix
$I+(D-\mu I)^{-1}(X^{-1}\Delta A X)$. Then the Neumann lemma
implies
$$
 1\leq \|(D-\mu I)^{-1}(X^{-1}\Delta A X)\|_p
  \leq \|(D-\mu I)^{-1}\|_p\|X^{-1}\|_p\|\Delta A\|_p\|X\|_p
$$
Lastly, observe that $(D-\mu I)^{-1}$ is diagonal, so
$$
  \|(D-\mu I)^{-1}\|_p=\max_{1\leq i\leq n}\frac{1}{|\lambda_i-\mu|}
  =\frac{1}{\min_{1\leq i\leq n}|\lambda_i-\mu|}
$$
The theorem now follows. \qed

\subsection{Corollary}

If $A$ is a normal matrix, then in the above theorem
$$
    \min_{1\leq i\leq n}|\lambda_i-\mu|
       \leq \|\Delta A\|_2
$$
because $X$ is a unitary matrix and so $\kappa_2(X)=1$.

\new Theorem~\ref{18.2} answers the question raised in \ref{18.11}, it
gives an estimate on the error in the eigenvalue in terms of
$\|\Delta A\|$ and $\kappa(X)$. However, this answer is not good
enough -- it gives one estimate for all eigenvalues. In practice,
some eigenvalues can be estimated much better than others. It is
important then to develop finer estimates for individual
eigenvalues.

\subsection{Left eigenvectors} \label{18.1}

Let $A\in\IC^{n\times n}$. A nonzero vector $x\in\IC^n$ is called
a {\em left eigenvector} of $A$ corresponding to an eigenvalue
$\lambda$ if
$$
      x^{\ast}A=\lambda x^{\ast}
$$

\noindent Note that this is equivalent to
$A^{\ast}x=\bar{\lambda}x$, i.e. $x$ being an ordinary (right)
eigenvector of $A^{\ast}$ corresponding to the eigenvalue
$\bar{\lambda}$.

\subsection{Lemma} \label{18.2}

A matrix $A$ has a left eigenvector corresponding to $\lambda$ if
and only if $\lambda$ is an eigenvalue of $A$ (a root of the
characteristic polynomial of $A$).

\proof $x^{\ast}A=\lambda x^{\ast}$ for an $x\neq 0$ is equivalent
to $(A^{\ast}-\bar{\lambda}I)x=0$, which means that
det$(A^{\ast}-\bar{\lambda}I)=0$, or equivalently, det$(A-\lambda
I)=0$, i.e. $C_A(\lambda)=0$. \qed

\new This explains why we do not introduce a notion of a
{\em left} eigenvalue: the eigenvalues for left eigenvectors and
those for ordinary (right) eigenvectors are the same.

\subsection{Lemma} \label{18.3}

For any eigenvalue $\lambda$ of $A$ the dimension of the ordinary
(right) eigenspace equals the dimension of the left eigenspace
(i.e., the geometric multiplicity of $\lambda$ is the same).

\proof ${\rm dim}\,{\rm Ker}(A-\lambda I)=n-\,{\rm rank}(A-\lambda
I)= n-\, {\rm rank}(A^{\ast}-\bar{\lambda}I)=\, {\rm dim}\,{\rm
Ker}(A^{\ast} -\bar{\lambda}I)$.

\subsection{Lemma} \label{18.13}

Let $A\in\IC^{n\times n}$. Then we have:
\begin{itemize}\item[(a)] If $\lambda$ is an eigenvalue with a right
eigenvector $x$, and $\mu\neq\lambda$ is another eigenvalue with a left
eigenvector $y$, then $y^{\ast}x=0$. \item[(b)] If $\lambda$ is a
simple eigenvalue (of algebraic multiplicity one) with right and left
eigenvectors $x$ and $y$, respectively, then $y^{\ast}x\neq
0$.\end{itemize}

\proof To prove (a), observe that $\la Ax,y\ra=\lambda\la x,y\ra$ and,
by a remark after \ref{18.1}, $\la x,A^{\ast}y\ra= \la
x,\bar{\mu}y\ra=\mu\la x,y\ra$. Hence, $\lambda\la x,y\ra=\mu\la
x,y\ra$, which proves (a), since $\lambda\neq\mu$.

To prove (b), assume that $\|x\|=1$. By the Schur decomposition
theorem, there is a unitary matrix $R$ with first column $x$ such that
$$
       R^{\ast}AR=\left[\begin{array}{ll}
       \lambda & h^{\ast} \\ 0 & B
       \end{array}\right]
$$
with some $h\in\IC^{n-1}$ and $B\in\IC^{(n-1)\times(n-1)}$. Note also
that $Re_1=x$. Since $\lambda$ is a simple eigenvalue of $A$, it is not
an eigenvalue of $B$, so the matrix $\lambda I-B$ is invertible, and so
is $\bar{\lambda}I-B^{\ast}$. Let $z= (\bar{\lambda}I-B^{\ast})^{-1}h$.
Then $\bar{\lambda}z-B^{\ast}z=h$, hence
$$
     h^{\ast}+z^{\ast}B=\lambda z^{\ast}
$$
One can immediately verify now that
$$
     [\,1\,z^{\ast}]\,R^{\ast}AR=\lambda\, [\,1\,z^{\ast}]
$$
Put $w^{\ast}=[\,1\,z^{\ast}]\, R^{\ast}$. The above equation can be
rewritten as
$$
      w^{\ast}A=\lambda w^{\ast}
$$
Hence $w$ is a left eigenvector of $A$. By the simplicity of
$\lambda$, the vector $w$ is a nonzero multiple of $y$. However,
observe that
$$
       w^{\ast}x=[\,1\,z^{\ast}]\,R^{\ast}Re_1=1
$$
which proves the lemma. \qed

\subsection{Theorem} \label{18.14}

Let $A\in\IC^{n\times n}$ have a simple eigenvalue $\lambda$ with right
and left unit eigenvectors $x$ and $y$, respectively. Let
$E\in\IC^{n\times n}$ such that $\|E\|_2=1$. For small $\varepsilon$,
denote by $\lambda(\varepsilon)$ and $x(\varepsilon)$, $y(\varepsilon)$
the eigenvalue and unit eigenvectors of the matrix $A+\varepsilon E$
obtained from $\lambda$ and $x$, $y$. Then
$$
       |\lambda'(0)|\leq\frac{1}{|y^{\ast}x|}
$$

\proof It follows from the inverse function theorem that
$\lambda(\varepsilon)$ and $x(\varepsilon)$ are differentiable for
sufficiently small $\varepsilon$. Write the equation
$$
     (A+\varepsilon E)\, x(\varepsilon)=
     \lambda(\varepsilon)\, x(\varepsilon)
$$
and differentiate it in $\varepsilon$, set $\varepsilon=0$, and get
$$
    Ax'(0)+Ex=\lambda'(0)x+\lambda\, x'(0)
$$
Then multiply this equation through on the left by the vector
$y^{\ast}$, use the fact that $y^{\ast}A=\lambda y^{\ast}$ and get
$$
       y^{\ast}Ex=\lambda'(0)\, y^{\ast}x
$$
Now the result follows since $|y^{\ast}Ex|=|\la Ex,y \ra|\leq \|Ex\|\,
\|y\| \leq 1$. Note that $y^{\ast}x\neq 0$ by \ref{18.13}. \qed

\new Note: The matrix $A+\varepsilon E$ is the perturbation
of $A$ in the direction of $E$. If the perturbation matrix $E$ is
known, one has exactly
$$
       \lambda'(0) =\frac{y^{\ast}Ex}{y^{\ast}x}
$$
and so
$$
        \lambda(\varepsilon)=\lambda+\frac{y^{\ast}Ex}{y^{\ast}x}\varepsilon
           +O(\varepsilon^2)
$$
by Taylor expansion, a fairly precise estimate on
$\lambda(\varepsilon)$. In practice, however, the matrix $E$ is
absolutely unknown, so one has to use the bound in \ref{18.14} to
estimate the sensitivity of $\lambda$ to all small perturbations of
$A$.

\subsection{Condition Number of An Eigenvalue}
\label{18.15}

Let $\lambda$ be a simple eigenvalue (of algebraic multiplicity one) of
a matrix $A\in\IC^{n\times n}$ and $x,y$ the corresponding right and
left unit eigenvectors. The {\em condition number} of $\lambda$ is
$$
      K(\lambda)=\frac{1}{|y^{\ast}x|}
$$

\new The condition number $K(\lambda)$ describes the sensitivity
of a (simple) eigenvalue to small perturbations of the matrix. Large
$K(\lambda)$ signifies an {\em ill-conditioned} eigenvalue.

\subsection{Simple properties of $K(\lambda)$} \label{18.16}

\begin{itemize}\item[(a)] Obviously, $K(\lambda)\geq 1$. \item[(b)]
If a matrix $A$ is normal, then $K(\lambda)=1$ for all its eigenvalues.
\item[(c)] Conversely, if a matrix $A$ has all simple eigenvalues with
$K(\lambda)=1$, then it is normal.\end{itemize}

Proof. See homework assignment.

\new Normal matrices are characterized by the fact that the
Schur decomposition $Q^{\ast}AQ=T$ results in a diagonal matrix $T$.
One can expect that if the matrix $T$ is nearly diagonal (its
off-diagonal elements are small), then the eigenvalues of $A$ are
well-conditioned. On the contrary, if some off-diagonal elements of $T$
are substantial, then at least some eigenvalues of $A$ are
ill-conditioned.

\subsection{Remark} \label{18.17}

It remains to discuss the case of multiple eigenvalues (of algebraic
multiplicity $\geq 2$). If $\lambda$ is a multiple eigenvalue, the left
and right eigenvectors may be orthogonal even if the geometric
multiplicity of $\lambda$ equals one. Example: $A=\begin{bmatrix} 0 & 1
\\ 0 & 0\end{bmatrix}$, the right and left eigenvectors are
$x=e_1$ and $y=e_2$, respectively. Moreover, if the geometric
multiplicity is $\geq 2$, then for any right eigenvector $x$ there
is a left eigenvector $y$ such that $y^{\ast}x=0$ (see homework
assignment). Hence, the definition \ref{18.15} gives an infinite
value of $K(\lambda)$.

This does not necessarily mean that a multiple eigenvalue
is always ill-conditioned. It does mean, however, that an
ill-conditioned simple eigenvalue is `nearly multiple'.
Precisely, if $\lambda$ is a simple
eigenvalue of $A$ with $K(\lambda)>1$, then there is a matrix
$E$ such that
$$
      \frac{\|E\|_2}{\|A\|_2}\leq\frac{1}{\sqrt{K(\lambda)^2-1}}
$$
and $\lambda$ is a multiple eigenvalue of $A+E$. We leave out the
proof. We will not further discuss the sensitivity of multiple
eigenvalues.

\subsection{Theorem (1st Gershgorin)} \label{18.18}

Let $A\in\IC^{n\times n}$ be `almost diagonal'. Precisely, let $A=D+E$,
where $D=\,$diag$\{d_1,\ldots,d_n\}$ and $E=(e_{ij})$ is small. Then
every eigenvalue of $A$ lies in at least one of the circular disks
$$
      D_i=\{z\colon\, |z-d_i|\leq\sum_{j=1}^n|e_{ij}|\}
$$
Note: $D_i$ are called {\em Gershgorin disks}.

\proof Let $\lambda$ be an eigenvalue of $A$ with eigenvector $x$. Let
$|x_r|=\max_i\{|x_1|,\ldots,|x_n|\}$ be the maximal (in absolute value)
component of $x$. We can normalize $x$ so that $x_r=1$ and $|x_i|\leq
1$ for all $i$. On equating the $r$-th components in $Ax=\lambda x$ we
obtain
$$
  (Ax)_r=d_rx_r+\sum_{j=1}^n e_{rj}x_j
        =d_r+\sum_{j=1}^n e_{rj}x_j=\lambda x_r=\lambda
$$
Hence
$$
      |\lambda-d_r|\leq\sum_{j=1}^n|e_{rj}|\, |x_j|
                   \leq\sum_{j=1}^n|e_{rj}|
$$
The theorem is proved.

\subsection{Theorem (2nd Gershgorin)} \label{18.19}

If $k$ of the Gershgorin disks $D_i$ form a connected region which is
isolated from the other disks, then there are precisely $k$ eigenvalues
of $A$ (counting multiplicity) in this connected region.

\new We call connected components of the union $\cup D_i$
{\em clusters}.


\proof For brevity, denote
$$
         h_i=\sum_{j=1}^n|e_{ij}|
$$
Consider a family of matrices $A(s)=D+sE$ for $0\leq s\leq 1$. It is a
standard fact in complex analysis that the roots of a complex
polynomial change continuously with its coefficients. Hence the
eigenvalues of $A(s)$ depend continuously on $s$. The Gershgorin disks
$D_i(s)$ for the matrix $A(s)$ are centered at $d_i$ and have radii
$sh_i$. As $s$ increases, the disks $D_i(s)$ grow concentrically, until
they reach the size of the Gershgorin disks $D_i$ of
Theorem~\ref{18.18} at $s=1$. When $s=0$, each Gershgorin disk $D_i(0)$
is just a point, $d_i$, which is an eigenvalue of the matrix $A(0)=D$.
So, if $d_i$ is an eigenvalue of multiplicity $m\geq 1$, then exactly
$m$ disks $D_j(0)$ coincide with the point $d_i$, which make a cluster
of $m$ disks containing $m$ eigenvalues (more precisely, one eigenvalue
of multiplicity $m$). As the disks grow with $s$, the eigenvalues
cannot jump from one cluster to another (by continuity), unless two
cluster overlap and then make one cluster. When two clusters overlap
(merge) at some $s$, they will be in one cluster for all larger values
of $s$, including $s=1$. This proves the theorem. \qed

\new Note: If the Gershgorin disks are disjoint, then each
contains exactly one eigenvalue of $A$.

\newpage

\section{Computation of eigenvalues: power method}
\label{19}

To simplify the matter, we always assume that the matrix $A$ is
diagonalizable, i.e. it has a complete set of eigenvectors $x_1,
\ldots,x_n$ with eigenvalues $\lambda_1,\ldots,\lambda_n$, which
are ordered in absolute value:
$$
      |\lambda_1|\geq |\lambda_2|\geq\cdots\geq |\lambda_n|
$$

\subsection{Dominant eigenvalue and eigenvector}
\label{19.2}

Assume that $|\lambda_1|>|\lambda_2|$, i.e.\ the largest eigenvalue is
simple. We call $\lambda_1$ the {\em dominant eigenvalue} and $x_1$ a
{\em dominant eigenvector}.

\subsection{Power method: the idea} \label{19.3}

Let $\lambda_1$ be the dominant eigenvalue of $A$ and
$$
       q=c_1x_1+\cdots +c_nx_n
$$
an arbitrary vector such that $c_1\neq 0$. Then
\begin{align*}
       A^kq &= c_1\lambda_1^kx_1+\cdots +c_n\lambda_n^kx_n \\
            &= \lambda_1^k[c_1x_1+c_2(\lambda_2/\lambda_1)^kx_2+
                \cdots +c_n(\lambda_n/\lambda_1)^kx_n]
\end{align*}
Denote
$$
       q^{(k)}=A^kq/\lambda_1^k=c_1x_1+
     \underbrace{c_2(\lambda_2/\lambda_1)^kx_2+
                \cdots +c_n(\lambda_n/\lambda_1)^kx_n}_{\Delta_k}
$$

\subsection{Lemma} \label{19.4}

The vector $q^{(k)}$ converges to $c_1x_1$. Moreover,
$$
       \|\Delta_k\|=\|q^{(k)}-c_1x_1\|\leq\,{\rm const}\cdot r^k
$$
where $r=|\lambda_2/\lambda_1|<1$.

\new Therefore, the vectors $A^kq$ (obtained by the powers
of $A$) will align in the direction of the dominant eigenvector $x_1$
as $k\to\infty$. The number $r$ characterizes the speed of alignment,
i.e. the speed of convergence $\|\Delta_k\|\to 0$. Note that if
$c_2\neq 0$, then
$$
      \|q^{(k+1)}-c_1x_1\|/\|q^{(k)}-c_1x_1\|\to r
$$
The number $r$ is called the {\em convergence ratio} or the {\em
contraction number}.

\subsection{Linear, quadratic and cubic convergence}
\label{19.5}

We say that the convergence $x_k\to x$ is {\em linear} if
$$
     |x_{k+1}-x|\leq r|x_k-x|
$$
for some $r<1$ and all sufficiently large $k$. If
$$
     |x_{k+1}-x|\leq C|x_k-x|^a
$$
with some $C>0$ and $a>1$, then the convergence is said to be {\em
superlinear} (it is faster than linear). For $a=2$ the convergence is
{\em quadratic}, and for $a=3$ {\em cubic}.

In order to come close to within $\varepsilon$ of the limit, the linear
convergence takes $k\sim\cO(|\log\varepsilon|)$ iterations, while the
superlinear convergence takes
$k\sim\cO\bigl(|\log(|\log\varepsilon|)|\bigr)$ iterations.

\subsection{Remarks on convergence} \label{19.6}

If the convergence is linear, then by induction we have
$|a_k-a|\leq Cr^k$, where $C=|a_0-a|$. The sequence $Cr^k$
decreases to zero exponentially fast, which is very fast by
calculus standards. However, in numerical algebra standards are
different.

The linear convergence practically means that each iteration adds a fixed
number of accurate digits to the result. For example, if $r=0.1$,
then each iteration adds one correct decimal digit. This is not
superb, since it will take 6--7 iterations to reach the maximum accuracy
in single precision arithmetic and 15--16 iterations in double
precision. If $r=0.5$, then each iteration adds one binary
digit (a bit), hence one needs 22--23 iterations in single
precision and 52--53 in double precision. Now imagine how long it
might take when $r=0.9$ or $r=0.99$.

On the other hand, the quadratic convergence means that each
iteration doubles (!) the number of digits of accuracy. Starting
with just one accurate binary digit, one needs 4--5 iterations in
single precision arithmetic and only one more iteration in double
precision. The cubic convergence means that each iteration triples
(!!!) the number of digits of accuracy. See Example~\ref{19.16}
for an illustration.

\subsection{Scaling problem in the power method} \label{19.7}

In practice, the vector $q^{(k)}=A^kq/\lambda_1^k$ is inaccessible
because we do not know $\lambda_1$ in advance. But we cannot just
drop the factor $\lambda_1^k$, because then $\|A^kq\|\to\infty$ if
$|\lambda_1|>1$ and $\|A^kq\|\to 0$ if $|\lambda_1|<1$, possibly
causing overflow or underflow in numerical computations. Thus we
must somehow normalize, or scale, the vector $A^kq$.

\subsection{Power method: two scaling factors} \label{19.8}

Pick an initial vector $q_0$. For $k\geq 1$, define
$$
        q_k=Aq_{k-1}/\sigma_k
$$
where $\sigma_k$ is a properly chosen scaling factor. One common
choice is
$$
   \sigma_k=\|Aq_{k-1}\|, \qquad
   \text{then}\qquad \|q_{k}\|=1
$$
In this case one can approximate the eigenvalue $\lambda_1$ by the
Rayleigh quotient
$$
      \lambda_1^{(k)}=q_k^{\ast}Aq_{k}.
$$
Another popular choice for $\sigma_k$ is the largest (in absolute
value) component of the vector $Aq_{k-1}$. This ensures that the
largest (in absolute value) component of $q_k$ equals one, in
particular $\|q_{k}\|_{\infty}=1$. Assume that the vector $x$ has
one component with the largest absolute value. In that case
$\sigma_k$ itself is a good approximation for $\lambda_1$, and we
set
$$
       \lambda_1^{(k)} = \sigma_k.
$$

To estimate how close the unit vector $q_k$ is to the
one-dimensional eigenspace span$\{x_1\}$, denote by $p_{k}$ the
orthogonal projection of $q_{k}$ on span$\{x_1\}$ and by
$d_{k}=q_{k}-p_{k}$ the orthogonal component. Then $\|d_{k}\|$
measures the distance from $q_{k}$ to span$\{x_1\}$.

\subsection{Theorem (convergence of the power method)} \label{19.9}

Assume that $\lambda_1$ is the dominant eigenvalue, and $q_{0}=\sum
c_ix_i$ is chosen so that $c_1\neq 0$. Then the distance from
$q_{k}$ to the eigenspace ${\rm span}\{x_1\}$ converges to zero and
$\lambda_1^{(k)}$ converges to $\lambda_1$. Furthermore,
$$
  \|d_{k}\|\leq {\rm const}\cdot r^k\ \ \ \ \ \ \ \
  |\lambda_1^{(k)}-\lambda_1|\leq {\rm const}\cdot r^k
$$

\new Note: The sequence of vectors $q_{k}$ need not have a
limit, see examples.

\proof It is a direct calculation, based on the representation
$A^kq_{0}=\lambda_1^k(c_1x_1+\Delta_k)$ of \ref{19.3} and
Lemma~\ref{19.4}.

\subsection{Examples} \label{19.11}

\begin{itemize}
\item[(a)] Let $A=\left[\begin{array}{rr} 3 & 2 \\ 1 &
1\end{array}\right]$. Pick $q_0=(1,1)^T$ and use the second choice
of $\sigma_k$ in \ref{19.8}. Then $\sigma_1=5$ and $q_1=(1,0.4)^T$,
$\sigma_2=3.8$ and $q_2= (1,0.368)^T$, $\sigma_3=3.736$ etc. Here
$\sigma_k$ converges to the dominant eigenvalue
$\lambda_1=2+\sqrt{3}=3.732$ and $q_k$ converges to a dominant
eigenvector $(1,\sqrt{3}/2-1/2)^T= (1,0.366)^T$.
\item[(b)] Let $A=\left[\begin{array}{rr} -1 & 0 \\ 0 &
0\end{array}\right]$. Pick $q_0=(1,1)^T$ and use the first choice
of $\sigma_k$ in \ref{19.8}. Then $q_k=((-1)^k,0)$ does not have a
limit, it oscillates. With the second choice of $\sigma_k$ in
\ref{19.8}, we have $q_k=(1,0)$ and $\sigma_k=-1=\lambda_1$ for
all $k\geq 1$.
\end{itemize}

\subsection{Initial choice} \label{19.12}

The choice of the initial vector $q_0$ only has to fulfill the
requirement $c_1\neq 0$. Since the vectors with $c_1=0$ form a
hyperplane in $\IC^n$, one hopes that a vector $q_0$ picked ``at
random'' will not lie in that hyperplane ``with probability one''.
Furthermore, even if $c_1=0$, round-off errors will most likely pull
the numerical vectors $q_k$ away from that hyperplane. If that does not
seem to be enough, one can carry out the power method for $n$ different
initial vectors that make a basis, say $e_1,\ldots,e_n$. One of these
vectors surely lies away from that hyperplane.



\subsection{Inverse power method} \label{19.13}

Assume that $A$ is invertible. Then $\lambda_1^{-1},\ldots,\lambda_n^{-1}$
are the eigenvalues of $A^{-1}$, with the same eigenvectors $x_1,\ldots,
x_n$. Note that $|\lambda_1^{-1}|\leq\cdots\leq |\lambda_n^{-1}|$.

Assume that $|\lambda_n^{-1}|>|\lambda_{n-1}^{-1}|$. Then
$\lambda_n^{-1}$ is the dominant eigenvalue of $A^{-1}$ and $x_n$ a
dominant eigenvector. One can apply the power method to $A^{-1}$ and
find $\lambda_n^{-1}$ and $x_n$. The rate of convergence of iterations
will be characterized by the ratio $r=|\lambda_n/\lambda_{n-1}|$. This
is called the {\em inverse power method}.

Note: in practice, there is no need to compute the inverse matrix
$A^{-1}$ explicitly. To find $A^{-1}b$ for any given vector $b$,
one can just solve the system $Ax=b$ by the LU decomposition.

Now we know how to compute the largest and the smallest eigenvalues.
The following trick allows us to compute any simple eigenvalue.

\subsection{Power method with shift} \label{19.14}

Recall that if $\lambda$ is an eigenvalue of $A$
with eigenvector $x$, then
$\lambda-\rho$ is an eigenvalue of $A-\rho I$
with the same eigenvector $x$.

Assume that $\rho$ is a good approximation to a simple eigenvalue
$\lambda_i$ of $A$, so that $|\lambda_i-\rho| < |\lambda_j-\rho|$ for all
$j\neq i$. Then the matrix $A-\rho I$ will have the smallest
eigenvalue $\lambda_i-\rho$ with the eigenvector $x_i$.

The inverse power method can now be applied to $A-\rho I$ to
find $\lambda_i-\rho$ and $x_i$. The convergence of iterations
will be linear with ratio
$$
     r=\frac{|\lambda_i-\rho|}{\min_{j\neq i}|\lambda_j-\rho|}
$$
Hence, the better $\rho$ approximates $\lambda_i$, the faster
convergence is guaranteed.

By subtracting $\rho$ from all the eigenvalues of $A$ we shift the
entire spectrum of $A$ by $\rho$. The number $\rho$ is called the
{\em shift}. The above algorithm for computing $\lambda_i$ and
$x_i$ is called the {\em (inverse) power method with shift}.

The power method with shift allows us to compute all the simple
eigenvalues and eigenvectors of a matrix, but the convergence is
slow (just linear).

\subsection{Power method with Rayleigh quotient shift}
\label{19.15}

This is an improvement of the algorithm \ref{19.14}. Since at each
iteration of the inverse power method we obtain a better
approximation to the eigenvalue $\lambda_i$, we can use it as the
shift $\rho$ for the next iteration. So, the shift $\rho$ will be
updated at each iteration. This will ensure a faster convergence.

One chooses an initial vector $q_0$ and an initial approximation
$\rho_0$, and for $k\geq 1$ computes
$$
     q_{k}=\frac{(A-\rho_{k-1}I)^{-1}q_{k-1}}{\sigma_k}
$$
and
$$
     \rho_{k}=\frac{q_k^{\ast}Aq_k}{q_k^{\ast}q_k}
$$
where $\sigma_k$ a convenient scaling factor, for example,
$\sigma_k=\|(A-\rho_{k-1}I)^{-1}q_{k-1}\|$.

The convergence of the Rayleigh quotient iterations is, generally,
quadratic (better than linear). If the matrix $A$ is Hermitian, the
convergence is even faster -- it is cubic.


\subsection{Example} \label{19.16}

Consider the symmetric matrix
$$
   A=\left[\begin{array}{rrr} 2 & 1 & 1 \\
   1 & 3 & 1 \\ 1 & 1 & 4
   \end{array}\right]
$$
and let $q_0=(1,1,1)^T$ be the initial vector. When Rayleigh quotient
iteration is applied to $A$, the following values $\rho_k$ are computed
by the first three iterations:
$$
   \rho_0=5,\ \ \ \ \rho_1=5.2131,\ \ \ \
   \rho_2=5.214319743184
$$
The actual value is $\lambda=5.214319743377$. After only three
iterations, Rayleigh quotient method produced 10 accurate digits. The
next iteration would bring about 30 accurate digits -- more than enough
in double precision.

\subsection{Power method: pros and cons} \label{19.17}

\begin{itemize} \item[(a)]
The power method is classic. It is very simple and generally good.
\item[(b)]
An obvious concern is the numerical stability of the method.
The matrices used in the inverse power method with shift tend to
be exceedingly ill-conditioned. As a result, the numerical
solution of the system $(A-\rho I)x=b$, call it $x_c$, will
deviate significantly from its exact solution $x$. However, for
some peculiar reason (we do not elaborate) the difference $x_c-x$
tends to align with the vector $x$. Therefore, the normalized
vectors $x_c/\|x_c\|$ and $x/\|x\|$ are close to each other.
Hence, ill conditioning of the matrices does not cause a
serious problem here.
\item[(c)]
On the other hand, the power method is slow. Each iteration requires
solving a linear system of equations $(A-\rho_{k-1}I)x=b$, every time
with a new matrix, so the LU decomposition must be repeated, which
takes $2n^3/3$ fops. If we want to compute all $n$ eigenpairs, and make
$p$ iterations per eigenpair, then the total cost is $2pn^4/3$ flops.
\item[(d)]
Lastly, when the matrix $A$ is real, then in order to compute its
complex eigenvalues one has to deal with complex matrices
$(A-\rho_{k-1}I)$, which is inconvenient and expensive. It would be
nice to stick to real matrices as much as possible and obtain pairs of
complex conjugate eigenvalues only at final steps. The QR algorithm, to
be discussed below, provides such a luxury.
\end{itemize}

\newpage

\section{Computation of eigenvalues: QR algorithm}
\label{19a}

The QR algorithm (not to be confused with QR decomposition!) dates back
to the early 1960s, and in the recent decades it became the most widely
used method for calculating the complete set of eigenvalues and
eigenvectors.

\subsection{Pure QR algorithm} \label{19.18}

Let $A\in\IC^{n\times n}$. The algorithm starts with $A_0=A$ and
generates a sequence of matrices $A_k$ defined as follows:
$$
         A_{k-1}=Q_kR_k,\qquad A_k=R_kQ_k.
$$
That is, a QR decomposition  of $A_{k-1}$ is computed and then its
factors are recombined in reverse order to produce $A_k$. One iteration
of the QR algorithm is called a {\em QR step}.

\subsection{Lemma} \label{19.19}

\begin{itemize} \item[(a)]
All matrices $A_k$ in the QR algorithm are unitary equivalent, in
particular, they have the same eigenvalues. \item[(b)] If $A$ is a
Hermitian matrix, then all $A_k$'s are Hermitian matrices as well.
\end{itemize}

\proof To prove (a), we note that $A_{k}=Q_k^{\ast}A_{k-1}Q_k$. Now (b)
follows by induction from $A_{k}^{\ast} = Q_k^{\ast} A_{k-1}^{\ast}Q_k
= Q_k^{\ast}A_{k-1}Q_k =A_k$. \qed

\subsection{Theorem (convergence of the QR algorithm)}
\label{19.20}

Let $\lambda_1,\ldots,\lambda_n$ be the eigenvalues of $A$
satisfying
$$
       |\lambda_1|>|\lambda_2|>\cdots >|\lambda_n|>0
$$
Under one technical assumption, see below, the matrix
$A_k=(a_{ij}^{(k)})$ will converge to an upper triangular form, so that
\begin{itemize}\item[(a)] $a_{ij}^{(k)}\to 0$ as $k\to\infty$ for all
$i>j$. \item[(b)] $a_{ii}^{(k)}\to\lambda_i$ as $k\to\infty$ for all
$i$. \end{itemize} The convergence is linear, with the ratio
$$
     r=\max_k|\lambda_{k+1}/ \lambda_k|<1.
$$

This theorem is given without proof. The technical assumption here is
that the matrix $Y$ whose $i$-th row is a left eigenvector of $A$
corresponding to $\lambda_i$ for all $i$, must have an LU decomposition
(i.e.\ all its principal minors must be nonsingular).

\subsection{Remarks} \label{19.21}

\begin{itemize}\item[(a)]
If $A$ is Hermitian, then $A_k$ will converge to a diagonal matrix.
\item[(b)] All the matrices $A_k$ (and $R_k$) involved in the QR
algorithm have the same 2-condition number (by \ref{15.6}), thus the QR
algorithm is numerically superior to the power method. This is also
true for all the variations of the QR method discussed below.
\item[(c)]
On the other hand, the pure QR algorithm is quite expensive.
Especially, each iteration is very costly: the QR decomposition takes
$2n^3$ flops and the multiplication of two matrices $R_k$ and $Q_k$
(even if we take advantage of the triangular structure of $R_k$!) takes
$n^3$ flops, a total of $3n^3$ flops. Thus if we make $p$ iterations
then the total cost is $3pn^3$. This seems to be an improvement over
$2pn^4/3$ flops of the power method, but we also need to remember that
the convergence is slow (linear), thus it may require many more
iterations, according to \ref{19.5}. Fortunately, the computational
cost can be reduced with the help of Hessenberg matrices, see below.
\item[(d)] The pure QR algorithm fails on matrices with multiple eigenvalues
(see Example~\ref{19.27a} below) and on real matrices with complex
eigenvalues. Indeed, complex eigenvalues of a real matrix come in
conjugate pairs $a\pm ib$, which have equal absolute values
$|a+bi|=|a-bi|$, violating the main assumption of the theorem.
Furthermore, if $A$ is real, then all $Q_k$, $R_k$ and $A_k$ are real
matrices as well, and thus we cannot even expect the real diagonal
elements of $A_k$ to converge to the complex eigenvalues of $A$. In
this case $A_k$ may not converge to anything.
\end{itemize}

\subsection{Hessenberg matrix} \label{19.22}

$A\in \IC^{n\times n}$ is called an {\em (upper) Hessenberg matrix} if
$a_{ij}=0$ for all $i>j+1$, i.e.\ $A$ has the form
$$
      A =
    \left[\begin{array}{ccccc}
       \times & \times & \cdots & \cdots & \times \\
       \times & \times & \ddots & \ddots & \vdots \\
       0 & \times & \times & \ddots & \vdots \\
       \vdots & \ddots & \ddots & \ddots & \vdots \\
       0 & \cdots & 0 & \times & \times
     \end{array}\right]
$$

\subsection{Lemma} \label{19.24}

If an invertible matrix $A_0$ is Hessenberg, then all the matrices
$A_k$ generated by the QR algorithm are Hessenberg.

\proof By induction, let $A_{k-1}$ be Hessenberg. Then $A_{k-1}=Q_kR_k$
and so $Q_k=A_{k-1}R_k^{-1}$. Since this is a product of a Hessenberg
matrix and an upper triangular matrix, it is verified by direct
inspection that $Q_k$ is Hessenberg. Then, similarly, $A_k=R_kQ_k$ is a
product of an upper triangular and Hessenberg matrices, so it is
Hessenberg. \qed \medskip

For noninvertible Hessenberg matrices, one can show that $Q_k$
constructed by Gram-Schmidt orthogonalization, see Chapter~\ref{QR}, is
also Hessenberg, thus $A_k$ will again be a Hessenberg matrix.

\subsection{Cost of a QR step for Hessenberg matrices}

For Hessenberg matrices, the QR algorithm can be implemented with a
substantial reduction of computational cost. First, the QR
decomposition via Givens rotators takes $6n^2$ flops, according to
\ref{16.25a}. Second, $Q_k$ is a product of $n-1$ rotators, hence the
multiplication of $R_k$ by $Q_k$ takes $\leq 6n^2$ flops. The total for
the QR step is then $\leq 12n^2$ flops, a dramatic drop from $3n^3$
flops as required for a non-Hessenberg matrix $A$.


\subsection{Theorem} \label{19.23}

Every matrix $A\in\IC^{n\times n}$ is unitary equivalent to a
Hessenberg matrix, i.e.
$$
        A=Q^{\ast}HQ
$$
where $H$ is a Hessenberg matrix and $Q$ is a unitary matrix. There is
an exact finite algorithm for computing $H$ and $Q$.

\medskip Note: the existence of $H$ immediately follows from Schur
decomposition $A=Q^{\ast}TQ$, so we only need to show that $H$ and $Q$
here can be constructed by a finite algorithm. Remember that there is
no finite algorithm for Schur decomposition, since it involves the
eigenvalues of the matrix $A$, cf.\ \ref{19.1}.

\proof An explicit algorithm for computation of $Q$ and $H$ is known as
{\em Arnoldi algorithm}. The matrix equation $A=Q^{\ast}HQ$ can be
rewritten as $AQ^{\ast}=Q^{\ast}H$. Denote by $q_i$, $1\leq i\leq n$,
the columns of the unitary matrix $Q^{\ast}$ and by $h_{ij}$ the
entries of $H$. Equating the columns of the matrices $AQ^{\ast}$ and
$Q^{\ast}H$ (and remembering that $h_{ij}=0$ for $i>j+1$) we obtain a
system of equations
\begin{align*}
    Aq_1 &= h_{11}q_1+ h_{2,1}q_{2}\\
    Aq_2 &= h_{12}q_1+ h_{22}q_2+h_{3,2}q_{3}\\
       & \cdots \\
    Aq_i &= h_{1i}q_1+ \cdots +h_{ii}q_i+h_{i+1,i}q_{i+1}\\
       & \cdots \\
    Aq_n &= h_{1n}q_1+ \cdots +h_{n-1,n}q_{n-1}+h_{n,n}q_{n}\\
\end{align*}
(note that the last equation is slightly different from the others
since it terminates on a diagonal entry of $H$).

Now the Arnoldi algorithm goes along the lines similar to the classical
Gram-Schmidt orthogonalization. We pick an arbitrary unit vector $q_1$,
compute $v_1=Aq_1$, and represent
$$
   v_1 = \,{\rm Pr}_{q_1}v_1 + w_2 = \la v_1,q_1\ra\,
   q_1 + w_2
$$
where $w_2$ is orthogonal to $q_1$. Then we set $h_{11}=\la v_1,q_1\ra$
and $h_{21}=\|w_2\|$ and define $q_2=w_2/\|w_2\|$. This enforces the
first equation in the above system.

Generally, for every $i=1,\ldots,n-1$ we make four steps:

\new {\bf Step 1}:  $\ \ $compute $v_{i}=Aq_{i}$.\\ {\bf Step
2}: $\ \ $for all $j=1,\ldots,i$ compute $h_{ji}=\la v_{i}, q_j\ra$.\\
{\bf Step 3}: $\ \ w_{i} = v_{i}-\sum_{j=1}^i h_{ji}q_j$ (note: this
vector is orthogonal to $q_1,\ldots,q_i$).\\ {\bf Step 4}: $\ \
h_{i+1,i}=\|w_{i}\|$ and $q_{i+1}=w_{i}/h_{i+1,i}$, unless
$h_{i+1,i}=0$, see below.\\

If $h_{i+1,i}=0$ in Step 4, we pick an arbitrary unit vector $q_{i+1}$
orthogonal to $q_1,\ldots,q_i$. Finally, for $i=n$ we execute steps 1
and 2 only. \qed

\new Theorem \ref{19.23} shows that one can first transform $A$ to a
Hessenberg matrix $A_0=H$, which has the same eigenvalues as $A$ does
(by similarity) and then start the QR algorithm with $A_0$.

\subsection{Remarks on terminology}

The exceptional case $h_{i+1,i}=0$ in Step 4 is referred to as the {\em
breakdown} of the Arnoldi algorithm. This term is quite misleading,
since the method does not really break down. In fact, the resulting
Hessenberg matrix $H$ will have a simpler structure (an extra zero on
its subdiagonal), and then the matrix $H$ has form
$$
       H = \begin{bmatrix} H_1 & C \\ 0 & H_2 \end{bmatrix}
$$
where $H_1\in\IC^{i\times i}$ and $H_2\in\IC^{(n-i)\times (n-i)}$.
Clearly, the set of eigenvalues of $H$ is the union of the eigenvalues
of $H_1$ and $H_2$, and those can be found by the QR algorithm applied
to the two smaller matrices $H_1$ and $H_2$ separately. This leads to a
reduction of the problem.

The Arnoldi algorithm is often applied to very large (or even
infinite) matrices, where the complete construction of $H$ and $Q$
is out of the question. Then one can run the Arnoldi algorithm
partway to obtain approximations to $H$ and $Q$. Such a method is
referred to as {\em Arnoldi iterations}.

\subsection{Cost of Arnoldi algorithm}

The cost of Step 1 is $2n^2$ flops, the cost of Step 2 is $2ni$ flops,
the same for Step 3, and lastly Step 4 takes $3n$ flops. The total cost
is then
$$
   \sum_{i=1}^n (2n^2+4ni+3n) \sim 4n^3
$$
(actually, by choosing $q_1=e_1$ one can save some work and compute $H$
and $Q$ in $\tfrac{10}{3}n^3$ flops, see the textbook, but this is
little improvement). This cost is comparable to the cost of one QR step
for a generic (non-Hessenberg) matrix. The Arnoldi algorithm can be
regarded as a {\em pre-processing} of the matrix $A$; it only needs to
be done once, and then the QR algorithm is applied to the resulting
Hessenberg matrix.


\subsection{The case of Hermitian matrices}

If $A$ is a Hermitian matrix, then $H$ will be both Hermitian and
Hessenberg. Hence $H$ will be a {\em tridiagonal matrix}
($h_{ij}=0$ for all $|i-j| > 1$). Its construction by the Arnoldi
algorithm takes only $\tfrac 43 n^3$ flops (see the textbook).

%\new There is a further improvement of the QR algorithm that accelerates its convergence:

\subsection{Theorem} \label{19.25}

Assume that $A_0$, and hence $A_k$ for all $k\geq 1$, are
Hessenberg matrices. Then the convergence $a_{i,i-1}^{(k)}\to 0$
as $k\to\infty$ in Theorem~\ref{19.20} is linear with ratio
$r_i=|\lambda_{i}/\lambda_{i-1}|$. In addition, the convergence
$a_{nn}^{(k)} \to \lambda_n$ is linear with ratio
$r_n=|\lambda_{n}/\lambda_{n-1}|$.

\medskip This theorem is given without proof.

\medskip Note that $|\lambda_{i+1}/\lambda_i|<1$ for all $i$. The smaller
this ratio, the faster the convergence. It is also important to
note that each subdiagonal entry $a_{i,i-1}^{(k)}$ has its own
rate of convergence. This allows us to accelerate the convergence
of some selected entries, see below.

\subsection{QR algorithm with shift - 1} \label{19.26a}

One can modify the matrix $A$ to decrease the ratio
$|\lambda_{n}/\lambda_{n-1}|$ and thus make the convergence
\beq
  a_{n,n-1}^{(k)}\to 0\qquad\text{and}\qquad
  a_{nn}^{(k)}\to\lambda_n
     \tag{B}
\eeq
of the two bottom entries faster with the help of shifting, as in
\ref{19.14}. To achieve this, one applies the QR steps to the
matrix $A-\rho I$ where $\rho$ is an appropriate approximation to
$\lambda_n$. Then the convergence (B) will be linear with ratio
$r=|\lambda_n-\rho|/|\lambda_{n-1}-\rho|$. The better $\rho$
approximates $\lambda_n$ the faster the convergence.

\subsection{QR algorithm with shift - 2} \label{19.26b}

The approximation $\rho$ can be updated at every iteration, as in
\ref{19.15}, by using Rayleigh quotient
$$
        \rho=\rho_k=u_k^{\ast}A_ku_k
$$
where $u_k$ is an approximate unit eigenvector of the matrix $A_k$
corresponding to the smallest eigenvalue $\lambda_n$. In practice, a
simple and convenient choice is $\rho_k = a_{nn}^{(k)}$, which
corresponds to $u_k=e_n$. Then the {\em QR algorithm with shift} then
goes as follows:
$$
      A_{k-1}-\rho_{k-1}I=Q_kR_k\qquad
      R_kQ_k+\rho_{k-1}I=A_k.
$$
This is called the {\em Rayleigh quotient shift}. The convergence
(B) is quadratic for the same reasons as in \ref{19.15}.

\subsection{QR algorithm with shift - 3} \label{19.26c}

However, the other subdiagonal entries, $a_{i+1,i}^{(k)}$, $1\leq
i\leq n-2$, move to zero slowly (linearly, with variable ratios).
To speed them up, one uses the following trick. After making
$a_{n,n-1}^{(k)}$ practically zero, one ensures that
$a_{nn}^{(k)}$ is practically $\lambda_n$. Then one can partition
the matrix $A_k$ as
$$
       A_k= \left[\begin{array}{cc}
          \hat{A}_k & b_k \\
       0 & \lambda_n \end{array}\right]
$$
where $\hat{A}_k$ is an $(n-1)\times (n-1)$ Hessenberg matrix,
whose eigenvalues are (obviously)
$\lambda_1,\ldots,\lambda_{n-1}$. Then one can apply further steps
of the QR algorithm (with shift) to the matrix $\hat{A}_k$ instead
of $A_k$. This quickly produces its smallest eigenvalue,
$\lambda_{n-1}$, which can be split off in the same manner. This
procedure is called the {\em deflation} of the matrix $A$.

In practice, each eigenvalue of $A$ requires just 3-5 iterations
(QR steps), on the average. For Hermitian matrices, it takes just
2-3 iterations per eigenvalue. Thus the QR algorithm with shift
and deflation achieves a top speed.

\new It remains to discuss the problem of computing multiple real
eigenvalues and complex eigenvalues for real matrices described in
Remark~\ref{19.21} (d).

\subsection{Example} \label{19.27a}

Let $A=\begin{bmatrix} 0&1 \\ 1&0\end{bmatrix}$. Then the pure QR
algorithm gives
$$
   A_0 = Q_1R_1 =
   \begin{bmatrix} 0&1 \\ 1&0\end{bmatrix}
   \begin{bmatrix} 1&0 \\ 0&1\end{bmatrix}
$$
hence
$$
   A_1 = R_1Q_1 =
   \begin{bmatrix} 1&0 \\ 0&1\end{bmatrix}
   \begin{bmatrix} 0&1 \\ 1&0\end{bmatrix} = A
$$
so the process goes nowhere. The Rayleigh quotient shift $\rho=a_{22}$
has no effect either, since $a_{22}=0$. The reason of this failure is
that the eigenvalues of $A$, which are $+1$ and $-1$, have equal
absolute values -- this is a symmetry which confuses the QR algorithm,
it ``cannot decide'' which eigenvalue to approach. To break this
symmetry, one needs to choose the shift $\rho$ differently.

\subsection{Wilkinson iteration} \label{19.27}

At step $k$ of the QR algorithm with shift, consider the trailing
$2\times 2$ matrix at the bottom right of $A_k$:
$$
     B_k  = \left[\begin{array}{ll}
     a_{n-1,n-1}^{(k)} & a_{n-1,n}^{(k)} \\
     a_{n,n-1}^{(k)} & a_{n,n}^{(k)}
     \end{array} \right]
$$
Now set $\rho_k$ to the eigenvalue of $B_k$ that is closer to
$a_{n,n}^{(k)}$ (in case of a tie, either one can be taken). This is
called the {\em Wilkinson shift}.

The eigenvalues of a $2\times 2$ matrix can be easily (and precisely)
computed by the quadratic formula, whether they are real or complex. If
they are real, then the Wilkinson shift will help to break the symmetry
in the case of multiple eigenvalues of $A$.

\subsection{Example~\ref{19.27a} continued}

The Wilkinson shift here is either $\rho=1$ or $\rho=-1$. Let us choose
$\rho=-1$. Then
$$
      A_0-\rho I=
   \begin{bmatrix} 1&1 \\ 1&1\end{bmatrix}
      =Q_1R_1=
   \left[\begin{array}{rr} 1/\sqrt{2}&1/\sqrt{2} \\
   1/\sqrt{2}&-1/\sqrt{2}\end{array}\right]
   \begin{bmatrix} \sqrt{2}&\sqrt{2} \\ 0&0\end{bmatrix}
$$
and then
$$
      A_1=R_1Q_1+\rho I =
   \begin{bmatrix} 2&0 \\ 0&0\end{bmatrix}+
   \left[\begin{array}{rr} -1&0 \\
   0&-1\end{array}\right] =
   \left[\begin{array}{rr} 1&0 \\
   0&-1\end{array}\right]
$$
so the QR algorithm converges in one step.

\subsection{Wilkinson iteration for complex eigenvalues} \label{19.27b}

If the eigenvalues of $B_k$ are complex, then the shift $\rho_k$ will
be complex, and the matrix $A-\rho_k I$ will be complex, too. The QR
step will then produce a complex matrix $A_{k+1}$. In this case one can
use the following trick to avoid further working with complex matrices:
set $\rho_{k+1}=\bar{\rho}_k$ (which is the other complex eigenvalue of
$B_k$) for the next QR step. Then the resulting matrix $A_{k+2}$ will
be real again, see bellow.

Furthermore, the values $\rho_k$ and $\bar{\rho}_k$ will approximate
two complex eigenvalues of the matrix $A$. Therefore, the QR algorithm
with Wilkinson shift is able to compute the complex eigenvalues of a
real matrix $A$, thus resolving the concern raised in \ref{19.17} (d).

Actually, there is no need to compute the complex matrix $A_{k+1}$. One
can just combine the two QR steps together and construct $A_{k+2}$
directly from $A_k$ (bypassing $A_{k+1}$). This can be carried out
entirely in real arithmetic. The resulting all-real procedure is called
the {\em double-step QR algorithm with Wilkinson shift}.

\subsection{Lemma}

Let $A_0\in\IR^{n\times n}$ and $\rho$, $\bar{\rho}$ not eigenvalues of
$A_0$. Consider a pair of QR steps with complex conjugate shifts:
\begin{align*}
      A_0-\rho I &=Q_1R_1\qquad
      R_1Q_1+\rho I=A_1 \\
      A_1-\bar{\rho} I &=Q_2R_2\qquad
      R_2Q_2+\bar{\rho} I=A_2
\end{align*}
Since the matrices $A_0-\rho I$ and $A_1-\bar{\rho}I$ are nonsingular,
the above QR decompositions may be constructed so that $R_1$ and $R_2$
have positive real entries, cf.\ \ref{QR6}. In this case $A_2$ will be
real.

\proof First, note that $A_1=Q_1^\ast A_0 Q_1$ and  $A_2=Q_2^\ast A_1
Q_2$. Now
\begin{align*}
  (A_0-\bar{\rho} I)(A_0-\rho I)
   &= (A_0-\bar{\rho} I)Q_1R_1 \\
   &= Q_1Q_1^\ast(A_0-\bar{\rho} I)Q_1R_1 \\
   &= Q_1(A_1-\bar{\rho} I)R_1 \\
   &= Q_1Q_2R_2R_1,
\end{align*}
which is a QR decomposition of the matrix $(A_0-\bar{\rho} I)(A_0-\rho
I)$, in which the upper triangular matrix $R_2R_1$ has real positive
diagonal entries. Since the matrix
$$
  (A_0-\bar{\rho} I)(A_0-\rho I) =
  A_0^2 - (\rho+\bar{\rho})A_0 + \rho\bar{\rho}I
$$
is obviously real, the above QR decomposition is unique by
Corollary~\ref{QR6}, thus it must be real, hence the matrices $Q_1Q_2$ and
$R_2R_1$ are entirely real. Thus $A_2 = (Q_1Q_2)^\ast A_0 (Q_1Q_2)$ is
real. \qed

\end{document}
\end


The characteristic polynomial
$$
         C_A(z) = \,{\rm det}(A-zI)
$$
has real coefficients and no (real) roots, so its roots are pairs of
conjugate complex numbers (because $C_A(z)=0\, \Leftrightarrow\,
C_A(\bar{z})=0$). So, $C_A(z)$ is a product of quadratic polynomials
with no real roots: $C_A(z)=P_1(z)\cdots P_k(z)$, where deg$\,
P_i(z)=2$ and $k=n/2$. By Cayley-Hamilton theorem $C_A(A)=0$, hence
$$
    P_1(A)\cdots P_k(A)=0
$$
Thus, at least one matrix $P_i(A)$, $1\leq i\leq k$, is singular
(otherwise the product $C_A(A)$ would be invertible). Let
$P_i(z)=z^2+az+b$. Then the matrix $A^2+aA+bI$ has a nontrivial kernel,
i.e. $A^2x+aAx+bx=0$ for some vector $x\neq 0$. Let $y=Ax$ (note that
$y\neq 0$, otherwise $x$ would be an eigenvector with eigenvalue zero).
We get $Ay+ay+bx=0$, so that
$$
    Ax=y\qquad {\rm and} \qquad  Ay=-ay-bx
$$
Hence, the subspace span$\{x,y\}$ is invariant under $A$. It is
two-dimensional, since $x$ and $y$ are linearly independent (otherwise
$Ax=y=\lambda x$, and $\lambda$ would be an eigenvalue of $A$). \qed
