\documentclass[twoside]{article}

\usepackage{palatino}
\usepackage{camnum}

\def\ba {\MM{a}}  
\def\bb {\MM{b}}  
\def\bu {\MM{u}}  
\def\bv {\MM{v}}  
\def\bw {\MM{w}}  
\def\bx {\MM{x}}  
\def\by {\MM{y}}  

\def\la {\langle}
\def\ra {\rangle}

\def\gl {\lambda}

\def\T {\CC{T}}

\def\R {\BB{R}}

\def\A {{\cal A}}

\textheight=25cm
\textwidth=15cm
\topmargin=-1.5cm
\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
      7\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}

\section{Linear least squares}

\subsection{Statement of the problem}

Consider the problem of finding a vector $\bx\in\R^n$ 
such that $A\bx = \bb$ where a matrix $A\in\R^{m\times n}$ 
and a vector $\bb\in\R^m$ are given and $m>n$. 
When there are more equations than unknowns,
the system $A\bx=\bb$ is called {\em overdetermined}.
In general, an overdetermined system has no solution, 
but we would like to have $A\bx$ and $\bb$ close in a sense.
Problems of this form occur frequently when we collect 
$m$ observations (which, typically, are prone to measurement error) 
and wish to exploit them to form an $n$-variable linear model
(e.g., trying to put some planet observartions on an ellipse).
Choosing the Euclidean distance as such a proximity  measure,
one obtains

The {\em least squares problem\/}. 
Given $A\in\R^{m\times n}$ and $\bb\in\R^m$, 
find $\bx^*\in\R^n$ which minimizes $\|A\bx-\bb\|$.

Let $A = [\ba_1,\ldots,\ba_n]$ with $\ba_i \in \R^m$. 
Then the vectors $A\bx = \sum_{i=1}^n x_i \ba_i$, 
where $\bx$ runs through $\R^n$, form an $n$-dimensional plane in $\R^m$, 
say $\A$, spanned by $(\ba_i)_{i=1}^n$. 
Given $\bb\in\R^m$, the minimal distance $\|\ba^* - \bb\|$ 
between $\bb$ and vectors in $\A$ is attained
when $\ba^*$ is the foot of the perpendicular 
from $\bb$ onto the plane $\A$.

{\bf Theorem.} 
Given $\bb \in \R^m$ and $\A = {\rm span}(\ba_i)_{i=1}^n$, 
the vector $\ba^*$ minimizes $\|\ba - \bb\|$ over all $\ba\in\A$ 
iff $\la \ba^* - \bb,\ba\ra = 0$ for all $\ba\in \A$. \\
%
\hspace*{8pt}{\bf Proof.} 
Given $\ba_*,\ba,\bb$, consider the quadratic polynomial in $\gl$,
$$
   p(\gl)
 = \| (\ba^*-\gl\ba) - \bb\|^2
 = \| \ba^*-\bb\|^2 - 2\gl \la \ba^* - \bb,\ba \ra + \gl^2 \|\ba\|^2.
$$
Its least value occurs when $p'(\gl_*) = 0$, so that
$$
    \gl_* = \frac{ \la \ba_* - \bb,\ba \ra}{\|\ba\|^2}, \quad
    p(\gl_*) = \| \ba^*-\bb\|^2 - \gl^2_* \|\ba\|^2.
$$ 
If $\la \ba^* - \bb,\ba \ra$ is nonzero for some $\ba$, so is $\gl_*$, 
and then $p(\gl_*) < p(0) = \| \ba^*-\bb\|^2$, i.e., $\ba^*$ is not optimal.
Conversely, if $\la \ba^* - \bb,\ba \ra = 0$ for all $\ba\in\A$, 
then, for any $\by\in\A$, with $\ba := \ba^* - \by$, we find that
$$
   \| \by - \bb\|^2 = \| (\ba^*-\ba) - \bb\|^2
 = \| \ba^*-\bb\|^2 + \|\ba\|^2 > \| \ba^* - \bb\|^2,
$$
unless $\by=\ba^*$. \QED

{\bf Corollary.} 
$\bx^* \in \R^n$ minimizes $\|A\bx - \bb\|$ 
iff $\la A\bx^* - \bb,A\bx\ra = 0$ for all $\bx\in \R^n$.

{\bf Corollary.} 
$\bx^* \in \R^n$ minimizes $\|A\bx - \bb\|$ 
iff $A^{\T}A\bx^* = A^{\T}\bb$.

The equation $A^{\T}A\bx = A^{\T}\bb$ is called the {\em normal equation}.


\subsection{Normal equations}

One way of finding optimal $\MM{x}$ is by solving the $n\times n$
linear system $A^{\CC{T}}A\MM{x}=A^{\CC{T}}\MM{b}$ -- the method of {\em normal
equations.\/} This approach is popular in many applications. However,
there are three disadvantages. Firstly, $A^{\CC{T}}A$ might be singular, 
secondly sparse $A$ might be replaced by a dense $A^{\CC{T}}A$ and, finally,
forming $A^{\CC{T}}A$ might lead to loss of accuracy. Thus, suppose
that our computer works in the IEEE arithmetic standard ($\approx15$
significant digits) and let
$$A=\left[\begin{array}{cc} 10^8 & -10^8 \\ 1 & 1\end{array}\right]
\qquad \Longrightarrow \qquad A^{\CC{T}}A=\left[\begin{array}{rr} 10^{16}+1 &
-10^{16}+1 \\ -10^{16}+1 & 10^{16}+1 \end{array}\right]\approx 10^{16}
\left[\begin{array}{rr}1 & -1 \\ -1 & 1\end{array}\right].$$
Given $\MM{b}=[0, 2]^{\CC{T}}$ the solution of $A\MM{x}=\MM{b}$ is
$[1,1]^{\CC{T}}$, as can be easily found by Gaussian elimination. However,
our computer `believes' that $A^{\CC{T}}A$ is singular!

\subsection{QR and least squares}

{\bf Lemma} Let $A$ be any $m\times n$ matrix and let
$\MM{b}\in\BB{R}^m$. The vector $\MM{x}\in\BB{R}^n$ minimises
$\|A\MM{x}-\MM{b}\|$ iff it minimises $\|\Omega A\MM{x}-\Omega
\MM{b}\|$ for an arbitrary $m\times m$ orthogonal matrix $\Omega$.\\
\hspace*{8pt} {\bf Proof.} Given an arbitrary vector
$\MM{v}\in\BB{R}^m$, we have
\begin{displaymath}
  \|\Omega\MM{v}\|^2=\MM{v}^{\CC{T}}\Omega^{\CC{T}}\Omega\MM{v}
  =\MM{v}^{\CC{T}} \MM{v}=\|\MM{v}\|^2.
\end{displaymath}
In particular, $\|\Omega A\MM{x}-\Omega\MM{b}\|=\|A\MM{x}-\MM{b}\|$. \QED

{\bf An irrelevant, yet important remark} The property that orthogonal
matrices leave the Euclidean distance intact is called
\emph{isometry\/} and it has many important ramifications throughout
mathematics and mathematical physics.

{\bf Method of solution} Suppose that $A=QR$, with an orthogonal $Q$
and an upper triangular $R$. Because of the lemma, letting
$\Omega:=Q^{\CC{T}}$, 
$$ \|A\MM{x}-\MM{b}\|=\|Q^{\CC{T}}(A\MM{x}-\MM{b})\|=\|
R\MM{x}-Q^{\CC{T}}\MM{b}\|,
$$
therefore we seek $\MM{x}\in\BB{R}^n$ that minimises
$\|R\MM{x}-Q^{\CC{T}}\MM{b}\|$.\\
Suppose 
for simplicity that $\mbox{rank}\,R=\mbox{rank}\,A=n$. Then the bottom
$m-n$ rows of $R$ are zero. Therefore we find $\MM{x}$ by solving the
(nonsingular) linear system given by the first $n$ equations of
$R\MM{x}=Q^{\CC{T}}\MM{b}$. Similar (although more complicated) algorithm
applies when $\mbox{rank}\,R\leq n-1$. Note, recalling our former
remark, that we don't require $Q$ explicitly and need to evaluate only
$Q^{\CC{T}}\MM{b}$. 

\section{Polynomial interpolation}

\subsection{The interpolation problem}

Given $n+1$ distinct real points $x_0,x_1,\ldots,x_n$ and real numbers
$f_0,f_1,\ldots,f_n$, we seek a function $p:\BB{R}\rightarrow\BB{R}$ such
that $p(x_i)=f_i$, $i=0,1,\ldots,n$. Such a function is called an {\em
interpolant.\/}\\
We denote by $\BB{P}_n[x]$ the set of all real polynomials of degree at
most $n$ and observe that each $p\in\BB{P}_n[x]$ is uniquely defined by
its $n+1$ coefficients. In other words, we have $n+1$ degrees of
freedom, while interpolation at $x_0,x_1,\ldots,x_n$ constitutes $n+1$
conditions. This, intuitively, justifies seeking an interpolant from
$\BB{P}_n[x]$. 

\subsection{The Lagrange formula}

Although, in principle, we may solve a linear problem with $n+1$
unknowns to determine a polynomial interpolant, this can be
accomplished more easily by using the explicit {\em Lagrange
formula.\/} We claim that
$$p(x)=\sum_{k=0}^n f_k\prod_{\stackrel{\scriptstyle\ell=0}{\ell\neq k}}^n
\frac{x-x_\ell}{x_k-x_\ell},\qquad x\in\BB{R}.$$
Note that $p\in\BB{P}_n[x]$, as required. We wish to show that it
interpolates the data. Define
$$
    L_k(x)
 := \prod_{\stackrel{\scriptstyle\ell=0}{\ell\neq k}}^n
    \frac{x-x_\ell}{x_k-x_\ell},\qquad
    k=0,1,\ldots,n
$$
({\em Lagrange cardinal polynomials\/}). 
It is trivial to verify that $L_k(x_k)=1$ and $L_k(x_j)=0$ for $k\neq j$,
hence 
$$
   p(x_j) = \sum_{k=0}^n f_k L_k(x_j)=f_j,\qquad j=0,1,\ldots,n,
$$
and $p$ is an interpolant,

{\bf Uniqueness.} Suppose that both $p\in\BB{P}_n[x]$ and $q\in\BB{P}_n[x]$
interpolate to the same $n+1$ data. Then the $n$th degree 
polynomial $p-q$ vanishes at $n+1$ distinct points. But the only
$n$th-degree polynomial with $\geq n+1$ zeros is the zero
polynomial. Therefore $p-q\equiv 0$ and the interpolating polynomial is
unique.

\end{document}