\documentclass[twoside]{article}

\usepackage{palatino}
\usepackage{camnum}

\textheight=23.5cm
\textwidth=15cm
\topmargin=-2cm
\oddsidemargin=0.5cm
\evensidemargin=0.5cm
\parindent=0pt
\parskip=6pt

\setcounter{section}{4}
\setcounter{subsection}{2}
\setcounter{equation}{0}

\begin{document}
\begin{center}
  {\large \textbf{Mathematical Tripos Part IB: Easter 2001}}

  {\Large \textbf{Numerical Analysis -- Lecture
      5\footnote{Corrections and suggestions to these notes should
        be emailed to \texttt{a.shadrin@damtp.cam.ac.uk}. All handouts
        are available on the WWW at the URL 
        \texttt{http://www.damtp.cam.ac.uk/user/na/PartIB/Handouts.html}.}}}
\end{center}

\subsection{The Gram--Schmidt algorithm}

Given a nonzero $m\times n$ matrix $A$ with the columns
$\MM{a}_1,\MM{a}_2,\ldots,\MM{a}_n\in\BB{R}^m$, we construct $Q$ \&
$R$ where $Q$ is orthogonal and $R$ upper-triangular and $A=QR$: in
other words,  
\begin{equation}
  \label{eq:4.1}
  \sum_{k=1}^\ell R_{k,\ell}\MM{q}_\ell=\MM{a}_k,\quad
  k=1,2,\ldots,n,\qquad \mbox{where}\qquad A=[\!
  \begin{array}{cccc}
    \MM{a}_1 & \MM{a}_2 & \cdots & \MM{a}_n
  \end{array}\!].
\end{equation}
Assuming $\MM{a}_1\neq\MM{0}$, we derive $\MM{q}_1$ and
$R_{1,1}$ from the equation \R{eq:4.1} for $k=1$. Since
$\|\MM{q}_1\|=1$, we let $\MM{q}_1=\MM{a}_1/\|\MM{a}_1\|$,
$R_{1,1}=\|\MM{a}_1\|$.\\
Next we form the vector $\MM{b}=\MM{a}_2-\langle\MM{q}_1,\MM{a}_2\rangle
\MM{q}_1$. It is orthogonal to $\MM{q}_1$, since
$$\langle\MM{q}_1, \MM{a}_2-\langle\MM{q}_1,\MM{a}_2\rangle
\MM{q}_1\rangle =\langle\MM{q}_1,\MM{a}_2\rangle
-\langle\MM{q}_1,\MM{a}_2\rangle \langle \MM{q}_1,\MM{q}_1\rangle=0.$$
If $\MM{b}\neq\MM{0}$, we set $\MM{q}_2=\MM{b}/\|\MM{b}\|$, hence
$\MM{q}_1$ and $\MM{q}_2$ are orthonormal. Moreover,
$$\langle\MM{q}_1,\MM{a}_2\rangle \MM{q}_1+\|\MM{b}\|\MM{q}_2
=\langle\MM{q}_1,\MM{a}_2\rangle \MM{q}_1+\MM{b}=\MM{a}_2,$$
hence, to obey \R{eq:4.1} for $k=2$, we let
$R_{1,2}=\langle \MM{q}_1,\MM{a}_2\rangle$, $R_{2,2}=\|\MM{b}\|$.

{\bf The Gram--Schmidt algorithm} The above idea can be extended to
all columns of $A$.

\vspace*{-3pt}
\textbf{\textit{Step 1\/}} Set $k:=0$, $j:=0$ ($k$ is the number of
columns of $Q$ that have been already formed and $j$ is the number of
columns of $A$ that have been already considered, clearly $k\leq j$);

\vspace*{-3pt}
\textbf{\textit{Step 2\/}} Increase $j$ by 1. If $k=0$ then set
$\MM{b}:=\MM{a}_j$, otherwise (i.e., when $k\geq1$) set
$R_{i,j}:=\langle\MM{q}_i,\MM{a}_j\rangle$, $i=1,2,\ldots,k$, and
$\MM{b}:=\MM{a}_j-\sum_{i=1}^k \langle\MM{q}_i,\MM{a}_j\rangle
\MM{q}_i$. {\em [Note: $\MM{b}$ is orthogonal to
  $\MM{q}_1,\MM{q}_2,\ldots,\MM{q}_k$.]\/}

\vspace*{-3pt}
\textbf{\textit{Step 3\/}} If $\MM{b}\neq\MM{0}$ increase $k$ by 1.
Subsequently, set $\MM{q}_k:=\MM{b}/\|\MM{b}\|$, $R_{k,j}:=\|\MM{b}\|$
and $R_{i,j}:=0$ for $i\geq k+1$. {\em [Note: Hence, each column 
of $Q$ has unit length, as required, $\MM{a}_j=\sum_{i=1}^k
R_{i,j}\MM{q}_j$ and $R$ is upper triangular, because $k\leq j$.]\/}

\vspace*{-3pt}
\textbf{\textit{Step 4\/}} Terminate if $j=n$, otherwise go to
\textbf{\textit{Step 2.\/}} 

Previous lecture $\Rightarrow$ Since the columns of $Q$ are
orthonormal, there are at most $m$ of them, i.e.\ the final value of
$k$ can't exceed $m$. If it is less then $m$ then a previous lemma
demonstrates that we can add columns so that $Q$ becomes $m\times m$
and orthogonal.

The disadvantage of Gram--Schmidt is its {\em ill-conditioning\/}.
Since we are using finite arithmetic, even small imprecisions in the
calculation of inner products rapidly lead to effective loss of
orthogonality. Thus, errors accumulate fast and even for
moderate values of $m$ it is no longer true that the computed off-diagonal
elements of $Q^{\CC{T}}Q$ are very small in magnitude.

On the other hand, orthogonality conditions are preserved well when one
generates a new orthogonal matrix by computing the product of two given
orthogonal matrices. Therefore algorithms that express $Q$ as a product of
simple orthogonal matrices are highly useful. This suggests an alternative
way forward.

\subsection{Orthogonal rotations}

Given real, $m\times n$ matrix $A_0=A$, we seek a sequence
$\Omega_1,\Omega_2,\ldots,\Omega_k$ of $m\times m$ orthogonal matrices
such that the matrix $A_i:=\Omega_i A_{i-1}$ has more zero elements below
the main diagonal than $A_{i-1}$ for $i=1,2,\ldots,k$ and so that the
manner of insertion of such zeros is such that $A_k$ is upper
triangular. We then let $R=A_k$, therefore
$$\Omega_k\Omega_{k-1}\cdots\Omega_2\Omega_1A=R$$
and $Q=(\Omega_k\Omega_{k-1}\cdots\Omega_1)^{-1}
=(\Omega_k\Omega_{k-1}\cdots\Omega_1)^{\CC{T}}=  
\Omega_1^{\CC{T}}\Omega_2^{\CC{T}}\cdots \Omega_k^{\CC{T}}$. Hence
$A=QR$, where $Q$ is 
orthogonal and $R$ upper triangular.

\subsection{Givens rotations}

We say that an $m\times m$ orthogonal matrix $\Omega_j$ is a {\em
Givens rotation\/} if it coincides with the unit matrix, except for
four elements. Specifically, we use the notation $\Omega^{[p,q]}$,
where $1\leq p<q\leq m$ for a matrix such that
$$\Omega^{[p,q]}_{p,p}=\Omega^{[p,q]}_{q,q}=\cos\theta,\qquad
\Omega_{p,q}^{[p,q]}=\sin\theta, \qquad
\Omega_{q,p}^{[p,q]}=-\sin\theta$$
for some $\theta\in[-\pi,\pi]$. The remaining elements of
$\Omega^{[p,q]}$ are those of a unit matrix. For example,
$$m=4\quad\Longrightarrow\quad \Omega^{[1,2]}=\left[
  \begin{array}{cccc}
    \cos\theta & \sin\theta &0&0\\
    -\sin\theta&\cos\theta & 0&0\\
    0&0&1&0\\
    0&0&0&1
  \end{array}\right],\quad 
\Omega^{[2,4]}=\left[
  \begin{array}{cccc}
    1&0&0&0\\
    0&\cos\theta&0&\sin\theta\\
    0&0&1&0\\
    0&-\sin\theta &0&\cos\theta 
  \end{array}\right].$$
Geometrically, such matrices correspond to the underlying coordinate
system being rigidly rotated along a two-dimensional plane. It is
trivial to confirm that they are orthogonal.

{\bf Theorem} Let $A$ be an $m\times n$ matrix. Then, for every $1\leq
p<q\leq m$, $i\in\{p,q\}$ and $1\leq j\leq n$, there exists
$\theta\in[-\pi,\pi]$ such that $(\Omega^{[p,q]}A)_{i,j}=0$. Moreover,
all the rows of $\Omega^{[p,q]}A$, except for the $p$th and the $q$th,
are the same as the corresponding rows of $A$, whereas the $p$th and
the $q$th rows are linear combinations of the `old' $p$th and $q$th
rows.\\
\hspace*{8pt}{\bf Proof.} Let $i=q$. If $A_{p,j}=A_{q,j}=0$ then any
$\theta$ will do, otherwise we let
$$\cos\theta:=A_{p,j}/\sqrt{A_{p,j}^2+A_{q,j}^2},\qquad
\sin\theta:=A_{q,j}/\sqrt{A_{p,j}^2+A_{q,j}^2}.$$
Hence 
$$(\Omega^{[p,q]}A)_{q,k}=-(\sin\theta)A_{p,k}+(\cos\theta)A_{q,k},\quad
k=1,2,\ldots,n \qquad\Rightarrow\qquad (\Omega^{[p,q]})_{q,j}=0.$$
Likewise, when $i=p$ we let  $\cos\theta:=A_{q,j}/\sqrt{A_{p,j}^2
+A_{q,j}^2}$, $\sin\theta:=-A_{p,j}/\sqrt{A_{p,j}^2+A_{q,j}^2}$.\\
The last two statements of the theorem are an immediate consequence of
the structure of $\Omega^{[p,q]}$. \QED

{\bf An example:} Suppose that $A$ is $3\times3$.
We can force zeros underneath the main diagonal as follows.\\[4pt]
$\MM{1}\hspace*{8pt}$ First pick $\Omega^{[1,2]}$ so that
$(\Omega^{[1,2]}A)_{2,1}=0$ $\quad\Rightarrow\quad$
$\Omega^{[1,2]}A=\left[ 
  \begin{array}{ccc}
    \times & \times & \times\\
    0 & \times & \times\\
    \times & \times & \times
  \end{array}\right]$.\\[4pt]
$\MM{2}\hspace*{8pt}$ Next pick $\Omega^{[1,3]}$ so that
$(\Omega^{[1,3]}\Omega^{[1,2]}A)_{3,1}=0$. Note that multiplication by
$\Omega^{[1,3]}$ doesn't alter the second row, hence
$(\Omega^{[1,3]}\Omega^{[1,2]}A)_{2,1}$ remains
zero $\quad\Rightarrow\quad$ $\Omega^{[1,3]}\Omega^{[1,2]}A =
\left[  
  \begin{array}{ccc}
    \times & \times & \times\\
    0 & \times & \times\\
    0 & \times & \times
  \end{array}\right]$.\\[4pt]
$\MM{3}\hspace*{8pt}$ Finally, pick $\Omega^{[2,3]}$ so that
$(\Omega^{[2,3]}\Omega^{[1,3]}\Omega^{[1,2]}A)_{3,2}=0$. Since both
second and third row of $\Omega^{[1,3]}\Omega^{[1,2]}A$ have a leading
zero, their linear combination preserves these zeros, hence also
$$(\Omega^{[2,3]}\Omega^{[1,3]}\Omega^{[1,2]}A)_{2,1}=
(\Omega^{[2,3]}\Omega^{[1,3]}\Omega^{[1,2]}A)_{3,1}=0.$$ 
It follows that $\Omega^{[2,3]}\Omega^{[1,3]}\Omega^{[1,2]}A$ is upper
triangular. Therefore
$$R=\Omega^{[2,3]}\Omega^{[1,3]}\Omega^{[1,2]}A =\left[ 
  \begin{array}{ccc}
    \times & \times & \times\\
    0 & \times & \times\\
    0 & 0 & \times
  \end{array}\right],\qquad
Q=(\Omega^{[2,3]}\Omega^{[1,3]}\Omega^{[1,2]})^{\CC{T}}.$$ 

\end{document}