\documentclass[11pt]{article}
\usepackage{amsmath,amsfonts,amsthm,amssymb,html}
\usepackage[bottom,none]{draftcopy}

\newcommand{\x}{\boldsymbol{x}}
\newcommand{\E}{\boldsymbol{e}}
\newcommand{\y}{\boldsymbol{y}}
\newcommand{\z}{\boldsymbol{z}}
\newcommand{\0}{\boldsymbol{0}}
\renewcommand{\a}{\boldsymbol{a}}
\renewcommand{\b}{\boldsymbol{b}}
\renewcommand{\c}{\boldsymbol{c}}
\renewcommand{\v}{\boldsymbol{v}}
\renewcommand{\u}{\boldsymbol{u}}
\newcommand{\w}{\boldsymbol{w}}
\newcommand{\h}{\boldsymbol{h}}
\newcommand{\g}{\boldsymbol{g}}
\newcommand{\Eta}{\boldsymbol{\eta}}
\newcommand{\e}{\mathrm{e}}

\newcommand{\R}{\mathbb{R}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\V}{\mathcal{V}}
\newcommand{\W}{\mathcal{W}}
\newcommand{\U}{\mathcal{U}}

\renewcommand{\d}{\mathrm{d}}
\newcommand{\Ker}{\operatorname{Ker}}
\newcommand{\Range}{\operatorname{Range}}
\newcommand{\Span}{\operatorname{Span}}
\newcommand{\true}{{\operatorname{true}}}
\newcommand{\diag}{{\operatorname{diag}}}
\newcommand{\trunc}{{\operatorname{trunc}}}
\newcommand{\noise}{{\operatorname{noise}}}
\newcommand{\opt}{{\operatorname{opt}}}

\renewcommand{\theenumi}{\roman{enumi}}
\renewcommand{\labelenumi}{(\theenumi)}

\newtheorem{theorem}{Theorem}
\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{exercise}[theorem]{Exercise}
\newtheorem{example}[theorem]{Example}

\begin{document}
\title{Least squares, the singular value decomposition, and linear
inverse problems}
\author{Marcel Oliver}
\date{Perspectives of Mathematics II -- Spring 2010}
\maketitle

\section{Simple least squares}
\label{s.simple}

The \emph{simple least squares} is the following special case of a
\emph{linear regression} problem: Find the equation of a line which is
``closest'' to a given set of points in the plane.  More precisely,
given tuples of real numbers $(x_1, y_1), \dots, (x_n, y_n)$, find
numbers $a$ and $b$ such that
\begin{equation}
  f(a,b) \equiv \sum_{i=1}^n (a \, x_i + b - y_i)^2
\end{equation}
is minimal.

As $f$ is a smooth function defined for all $(a,b) \in \R^2$, calculus
tells us that it can only have a minimum provided its partial
derivatives vanish.  I.e.,
\begin{subequations}
\begin{gather}
  0 = \frac{\partial f}{\partial b}
    = 2 \sum_{i=1}^n a \, x_i + b - y_i 
  \label{e.dfdb}
\intertext{and}
  0 = \frac{\partial f}{\partial a}
    = 2 \sum_{i=1}^n x_i \, (a \, x_i + b - y_i) \,.
  \label{e.dfda}
\end{gather}
\end{subequations}
Defining the \emph{mean} of the $\{x_i\}$ and $\{y_i\}$ as
\begin{gather}
  \bar x = \frac1n \sum_{i=1}^n x_i
  \qquad \text{and} \qquad
  \bar y = \frac1n \sum_{i=1}^n y_i \,,
\end{gather}
equation \eqref{e.dfdb} can be written
\begin{subequations}
\begin{gather}
  b = \bar y - a \, \bar x \,.
  \label{e.b}
\end{gather}
Substituting into \eqref{e.dfda} and solving for $a$, we find
\begin{align}
  a 
  & = \frac{\displaystyle \sum_{i=1}^n x_i \, y_i 
            - n \, \bar x \, \bar y}%
           {\displaystyle \sum_{i=1}^n x_i^2 - n \, \bar x^2}
%      \notag \\
    = \frac{\displaystyle \sum_{i=1}^n (x_i - \bar x)(y_i - \bar y)}%
           {\displaystyle \sum_{i=1}^n (x_i - \bar x)^2} \,.
  \label{e.a}
\end{align}
\end{subequations}
Thus, we can first compute $a$ via \eqref{e.a} and then $b$ via
\eqref{e.b}. 

Two remarks are in order.  First, vanishing of the partial derivatives
is, in general, only a necessary condition for the existence of a
minimum.  Here, however, it is not difficult to show that $f(a,b) \to
\infty$ as $a,b \to \infty$ (in any way) provided there are at least
two distinct data points.  Thus, since the necessary condition yields
a unique solution, this solution must correspond to a minimum of $f$.  

Second, the minimization of the Euclidean norm of the error as opposed
to some other error measure is certainly convenient as it yields a
smooth and even linear relationship between measured $y$-values and
estimates for the coefficients.  However, it is often also preferred
for statistical reasons:  When the errors are assumed Gaussian, the
least-squares fit equals the so-called maximum-likelihood estimator.
When the errors are all from the same probability distribution,
least-squares is the so-called best unbiased linear estimator.   We
will not dwell on these issues for the time being.


\section{Calculus on normed vector spaces}

In the following, we give a simple framework for doing calculus on
arbitrary normed\footnote{A normed vector space $\V$ is a vector
space together with a map $\lVert \, \cdot \, \rVert \colon \V \to
[0,\infty)$ such that
\begin{enumerate}
\item $\lVert v \rVert = 0$ if and only if $v=0$,
\item $\lVert \lambda v \rVert = \lvert \lambda \rvert \, \lVert v
\rVert$ for all $\lambda \in \R$ and $v \in \V$, and
\item $\lVert v+w \rVert \leq \lVert v \rVert + \lVert w \rVert$ for
all $v,w \in \V$.
\end{enumerate}
The last property is called the \emph{triangle inequality}.  When $\V
\in \R^n$, we take without further discussion the \emph{Euclidean
Norm} defined via
\[
  \lVert \v \rVert^2
  = \sum_{i=1}^n v_i^2
  = \v^T \v
\]
for $\v = (v_1, \dots, v_n)^T \in \R^n$.}  vector spaces.\footnote{For
simplicity, we only consider vector spaces over $\R$.  In case that
the space is infinite dimensional, we would typically assume
completeness, i.e., that we have a Banach space.}  A typical first
year course on multivariate calculus typically considers the vector
space $\R^n$ only.  However, looking at more general vector spaces is
in some sense easier, and naturally includes differentiation with
respect to matrices as well as the calculus of variations, both of
which occur frequently in applications.  In this basic introduction,
we make no attempt to give a rigorous definition of differentiability
nor do we discuss any of the analytical subtleties of doing calculus
on infinite dimensional vector spaces such as function spaces.

Let $\V$ and $\W$ be vector spaces and $f \colon \V \to \W$ a
continuous function.  Central to the understanding of differentiation
in this setting is the notion of \emph{directional derivative}.
Namely, fix $v \in \V$ and choose an arbitrary $\delta v \in
\V$.\footnote{We must issue an important warning: In typical
generalizations, $v$ and $\delta v$ are taken from different sets.
For example, $\V$ may be an affine space $a + \U$.  Then the
``admissible variations'' $\delta v$ must be taken from $\U$.  Or $\V$
may be a manifold, then $\delta v$ is an element from its tangent
space $T\V$.}  We define the directional derivative as the
instantaneous rate of change of $f$ when changing $v$ by one unit of
$\delta v$.  We write
\begin{gather}
  \d f (v)[\delta v]
  = \lim_{h \to 0} \frac{f(v + h \, \delta v) - f(v)}h
  = \frac{\d f(v + h \, \delta v)}{\d h} \bigg|_{h = 0} \,.
  \label{e.df}
\end{gather}
If $\d f$ is continuous on $\V \times \V$, for each $v \in \V$, the
map $\d f(v)[\,\cdot\,] \colon \V \to \W$, the called the \emph{total
derivative}, is linear.\footnote{In the context of infinite
dimensional Banach spaces, the directional derivative is called the
G\^ateaux derivative and the total derivative leads to so-called
Fr\'echet derivative.  For G\^ateaux differentiability to imply
Fr\'echet differentiability, one needs, moreover, continuity of $\d f
\colon \V \to \mathcal L(\V,\W)$.}  I.e., for fixed $v \in \V$,
\begin{subequations}
\begin{gather}
  \d f(v)[\lambda \delta v] 
  = \lambda \, \d f(v)[\delta v] \,, \\
  \d f(v)[\delta v + \delta w] 
  = \d f(v)[\delta v] + \d f(v) [\delta w]
\end{gather}
\end{subequations}
for all $\delta v, \delta w \in \V$ and $\lambda \in \R$.

\paragraph{Practical computation}  

The following trick often simplifies the computation of the
derivative.  Writing
\begin{gather}
  \delta f \equiv \d f(v)[\delta v] \,,
  \label{e.deltaf}
\end{gather}
the $\delta$ symbol formally behaves like differentiation,\footnote{In
fact, with only a change of language we could call this a
differential.} so we can apply the product rule, chain rule, etc.\ as
usual.

\begin{example}
On $\V=\R^n$, fix an $A \in M(n \times n)$ and define
\begin{equation}
  f(\v) = \v^T A \v \,.
\end{equation}
Then, using the product rule, 
\begin{equation}
  \delta f = (\delta \v)^T A \v + \v^T A \delta \v
           = \v^T A^T \delta \v + \v^T A \delta \v
\end{equation}
so that $\d f(\v)[\delta \v] = \v^T(A+A^T) \delta \v$ or $\d f(\v) =
\v^T(A+A^T)$.  You probably know that the derivative of a real-valued
function on $\R^n$ can be computed via its Jacobian, which here is the
$1 \times n$ matrix
\begin{equation}
  \d f(\v) = (\partial_1 f, \dots, \partial_n f) \,.
\end{equation}
Sure enough, using this formula, you'll find the same answer, but with
a lot more pain.  Try it!
\end{example}

\begin{example}
Let 
\begin{equation}
  \V = \{ \phi \in C^2([0,1]) \colon \phi(0)=\phi(1)=0 \} 
\end{equation}
and fix $g \in C([0,1])$.  Define
\begin{equation}
  f(\phi) = \int_0^1 
  \bigl( \tfrac12 \, (\phi'(x)^2 + g(x) \, \phi(x) \bigr) \, \d x \,. 
\end{equation}
Then
\begin{equation}
  \delta f = \int_0^1
  \bigl( \phi'(x) \, \delta \phi'(x) + g(x) \, \delta \phi(x) \big) \,
  \d x 
  = \int_0^1 (-\phi''(x) + g(x)) \, \delta \phi(x) \, \d x 
  \label{e.functional-deriv}
\end{equation}
where, in the last equality, we have used integration by parts noting
that the boundary terms vanish due to the definition of $\V$.  Thus,
\eqref{e.functional-deriv} is an expression for $\d f(\phi)[\delta
\phi]$; we note that there is no elementary way to separate out
$\delta \phi$ in this expression.
\end{example}

\paragraph{Extrema} 

We now specifically consider maps $f \colon \V \to \R$, as in the
examples above, and ask for conditions under which such maps assume
local minima or maxima.  Throughout, we assume that $f \in
C^1(\V,\R)$, i.e., $f$ is at least once continuously differentiable.
In this case, $\d f(v)$ is a linear functional (i.e., a linear map
from $\V$ to $\R$) called the \emph{gradient} of $f$.  We say that $f$
has a local maximum at $v \in \V$ if there is $\varepsilon>0$ such
that $f(v) \geq f(w)$ for all $\lVert v-w \rVert < \varepsilon$; $f$
has a local minimum if $-f$ has a local maximum.  Then we have the
following necessary condition for the existence of a local extreme
point:

\begin{quote}
$f$ has a local minimum or maximum at $v \in V$ only if $V$ is a
\emph{stationary point} of $f$, i.e., if $\d f(v)=0$.
\end{quote}

The proof is a direct extension of the proof of the ``first
derivative test'' from single-variable calculus and shall not be given
here. 

Note that if $f$ is only defined on a domain $D \subsetneq V$, then it
may assume a local or global extremum on the boundary (with obvious
modification to the definition of the local extremum).  This case may
have to be considered separately.  The technique of Lagrange
multipliers below can be useful if the boundary can be described by a
constraint. 

\paragraph{Lagrange Multipliers}

Sometimes we need to find extreme values of $f$ where $v$ is subject
to additional constraints.  In the simplest case, assume that the
constraint is given by the equation $g(v)=0$ for some $g \in
C^1(\V,\R)$.  We have the following necessary condition.

\begin{quote}
$f$ has a local minimum or maximum on the set $\{v \in V \colon
g(v)=0\}$ only if there exists $\lambda \in \R$ such that 
\begin{equation}
  \d f(v) =\lambda \, \d g(v) \,.
  \label{e.lagrange}
\end{equation}
\end{quote}

We will not give a proof, but shall give the following geometric
motivation.  For the purposes of discussion, suppose that $\V =
\R^n$.  First, we notice that the gradient of a function is a vector
pointing in the direction of steepest ascent.\footnote{To see this,
recall that in Euclidean geometry,
\[
  \u^T \w = \lVert \u \rVert \, \lVert \w \rVert \, 
  \cos \angle(\u,\v) \,.
\]
Thus, if $\u = \d f(\v)$, the instantaneous rate of change in the
direction of $\delta \v$ can be written $\u^T \delta \v$ and is
therefore maximal among all vectors of length $\lVert \u \rVert$ if
$\delta \v = \u$.}  Thus, if we vary $\v$ in the constraint set $\{v
\in V \colon g(v)=0\}$ in the direction $\delta \v = \d g(\v)^T$ leave
the constraint set.  On the other hand, any variation $\delta \v \bot
\d g(\v)^T$ is tangent to the constraint set and therefore permitted.
Thus, when ``probing'' $f$ for a local extremum, we can allow $f$ to
ascent or decent in the ``forbidden'' direction $\d g(\v)^T$.  On the
other hand, when $f$'s steepest ascent direction $\d f(\v)$ aligns
with the ``forbidden'' direction $\d f(\v)$, then the instantaneous
rate of change of $f$ is zero in any direction tangent to the
constraint set.  This alignment is represented by the Lagrange
multiplier condition \eqref{e.lagrange}.

This argument also points out how to generalize to $K$ constraints:
Each constraint function $g_1, \dots, g_K$ induces a ``forbidden''
direction $\d g_i(\v)^T$ which must be excluded from the ``search''
when probing $f$ for an extremum.  Thus, the ascent direction of $f$
may lie in the ``forbidden'' space spanned by all the $\d g_i(\v)^T$.
Thus, a necessary condition for a local extremum under $K$ constraints
is the existence of a Lagrange multiplier $\boldsymbol \lambda \in
\R^K$ such that
\begin{equation}
  \d f(v) = \boldsymbol \lambda^T \d \g(v) \,.
  \label{e.lagrange-2}
\end{equation}


\section{General linear least squares}
\label{s.glls}

Recall that in Section~\ref{s.simple}, the model for our data is the
linear equation
\begin{equation}
  y = a \, x + b \,.
\end{equation}
We generalize as follows.  Assume $\x \in \R^m$ with $m \in \N$ while
we still take $y \in \R$.  (The situation discussed here easily
generalizes to vector-valued ``output data'' $y$ by considering each
of its components separately.)  So our data is again a set of tuples
$(\x_1, y_1), \dots, (\x_n, y_n)$.  A \emph{linear model} for our data
is a linear combination with unknown coefficients of a fixed finite
set of possibly nonlinear ``ansatz functions'' $h_i \colon \R^m \to
\R$ for $i=1, \dots, k$.  I.e., we seek to find $\v \in \R^k$ such
that 
\begin{equation}
  y = \h(\x)^T \v
\end{equation}
with $\h(\x) = (h_1(\x), \dots, h_k(\x))^T$ reproduces the data with
minimal mean-square error
\begin{align}
  f(\v) 
  & = \sum_{i=1}^n (\h(\x_i)^T \v - y_i)^2 \,.
\end{align}
Setting
\begin{equation}
  H = 
  \begin{pmatrix}
    h_1(\x_1) & \dots & h_k(\x_1) \\
    \vdots & & \vdots \\
    h_1(\x_n) & \dots & h_k(\x_n)
  \end{pmatrix} \qquad \text{and} \qquad
  \y = 
  \begin{pmatrix}
    y_1 \\ \vdots \\ y_n
  \end{pmatrix} \,,
\end{equation}
we can interpret this problem as the minimization of the so-called
residual $H\v - \y$ corresponding to the potentially inconsistent
linear system $H\v = \y$.  (In the following, nothing depends on the
fact that these equations come from a general linear fitting problem;
we are generally asking for the least-square solution of some
inconsistent system of linear equations.)  In matrix notation, we ask
for the minimum of
\begin{equation}
  f(\v) = \lVert H \v - \y \rVert^2 
        = (H \v - \y)^T (H \v - \y) \,.
\end{equation}
The necessary condition for a minimum of $f$ is
\begin{equation}
  0 = \delta f 
    = (H \delta \v)^T (H \v - \y) + (H \v - \y)^T H \delta \v
    = 2 \, (H \v - \y)^T H \delta \v
\end{equation}
for any $\delta \v \in \R^k$, which implies $(H \v)^T H = \y^T H$ or
\begin{equation}
  (\y - H \v)^T H = \boldsymbol 0 \,.
  \label{e.orth}
\end{equation}

\paragraph{Geometric interpretation}

In the language of Linear Algebra, equation \eqref{e.orth} can be seen
as an orthogonality condition, namely\footnote{We define, as usual,
\[
  \Range H = \{ H \u \colon \u \in \R^k \}
\]
and
\[
  \Ker H = \{ \u \in \R^k \colon H \u = \boldsymbol 0 \} \,.
\]}
\begin{equation}
  \y - H \v \bot \Range H \,,
\end{equation}
i.e., the process of minimizing $f(\v)$ yields a decomposition
\begin{equation}
  \y = \underbrace{\y - H \v}_{\bot \Range H}
     + \underbrace{H \v {\vphantom\y}}_{\in \, \Range H} \,.
  \label{e.range-decomposition}
\end{equation}
Thus, the point $H\v$ is the point on the hyperplane $\Range H$ of
closest Euclidean distance to $\y$.

\paragraph{Solvability}

In practical terms, we solve \eqref{e.orth} by solving the linear
system
\begin{equation}
  H^T H \v = H^T \y \,,
  \label{e.normal}
\end{equation}
called the \emph{normal equations}.  It is an easy exercise in Linear
Algebra to verify that $\Range H^T = \Range H^TH$, therefore
\eqref{e.normal} has a solution $\v$ for any $\y \in \R^n$.

In the typical case, the $k \times k$ matrix $H^T H$ is invertible, so
that we can write
\begin{equation}
  \v = (H^T H)^{-1} H^T \y \,.
  \label{e.v}
\end{equation}
This equation can be solved directly.  (For big problems, matrix
inversion is expensive and should be replaced by more efficient
computational methods, for details see \cite{B}.)

\begin{example}
We can recover the result of Section~\ref{s.simple} by setting $m=1$,
$k=2$, $h_1(x) = x$, $h_2(x) = 1$, and $\v=(a,b)^T$.  Equation
\eqref{e.v} is thus a compact way of writing \eqref{e.a} and
\eqref{e.b}.
\end{example}

\begin{example}
Suppose we would like to fit a polynomial of degree $k-1$ to a set of
data points taken at integer points $x_1 = 0, \dots, x_{n}=n-1$.  If
we choose as ansatz functions $h_1(x)=1, h_2(x)=x, \dots,
h_{k}(x)=x^{k-1}$, the matrix $H$ reads
\begin{equation}
  H = 
  \begin{pmatrix}
    1 & 0 & 0 & \cdots & 0  \\
    1 & 1 & 1 & \cdots & 1 \\
    1 & 2 & 2^2 & \cdots & 2^{k-1} \\
    \vdots & & & & \vdots \\
    1 & n-1 & (n-1)^2 & \cdots & (n-1)^{k-1}
  \end{pmatrix}
\end{equation}
This matrix is known as the \emph{Vandermonde matrix}.  It is known to
be poorly conditioned [to be discussed later in the course] for even
moderate values of $k$; as a result, polynomial least squares with a
degree higher than about $3$--$6$ are almost never a good idea.
\end{example}

When $H^T H$ is not invertible, i.e., when $\Ker H \neq \{ \boldsymbol
0 \}$, there is linear dependence in the data points or in the ansatz
functions.  In this case, the effective number of independent
observations is less than $k$.  In this case, we may wish to compute a
least-norm solution of \eqref{e.normal}, discussed in the next
section.


\section{Least norm solutions}

When a system of linear equations
\begin{equation}
  A \v = \b \,,
  \label{e.linear}
\end{equation}
where $A \in M(n \times k)$ is consistent, but does not have a unique
solution (we can express this situation in technical terms by saying
that $\b \in \Range A$ and $\Ker A \neq \{\boldsymbol 0\}$), we may
ask for a solution $\v$ which has the least norm amongst the family of
solutions.  (In the context of the general linear least squares
problem of Section~\ref{s.glls}, $A=H^H$ and $\b = H^T \y$.)

\paragraph{Geometric interpretation} 

If $\w$ is any solution of $A\w=\b$, then it can be decomposed as
\begin{equation}
  \w = \underbrace{\w - \v}_{\in \Ker A} + 
       \underbrace{\v}_{\bot \Ker A} \,,
\end{equation}
where $\v$ is the minimum norm solution.  Thus, $\lVert \v \rVert$ is
the distance of $\w$ to the hyperplane $\Ker A$.  This decomposition
is the domain-side analog to the range-side decomposition
\eqref{e.range-decomposition} that comes up when computing
least-square solutions. 

\paragraph{Formulation as constrained minimization}

The least norm problem can be formulated as the following constrained
minimization problem.
\begin{quote}
Find $\v$ minimizing $f(\v) = \v^T \v$ subject to $\g(\v) = A \v
- \b = 0$.
\end{quote}
Since $\d f(\v) = 2 \v^T$ and $\d \g(\v) = A$, the method of Lagrange
multipliers \eqref{e.lagrange-2} yields the necessary condition   $2
\v^T = \boldsymbol\lambda^T A$ or, setting $\boldsymbol \mu = -2
\boldsymbol \lambda$,
\begin{equation}
  I \v + A^T \boldsymbol \mu = \boldsymbol 0 \,.
\end{equation}
Together with the constraint $A\v = \b$, we can write the necessary
condition in block-matrix form,
\begin{equation}
  \begin{pmatrix}
    I & A^T \\ A & 0
  \end{pmatrix}
  \begin{pmatrix}
    \v \\ \boldsymbol \mu
  \end{pmatrix}
  = 
  \begin{pmatrix}
    \boldsymbol 0 \\ \b
  \end{pmatrix} \,.
  \label{e.least-norm}
\end{equation}
It is easy to check that this system is consistent as long as
\eqref{e.linear} is consistent.  Moreover, 
\begin{equation}
  \Ker 
  \begin{pmatrix}
    I & A^T \\ A & 0
  \end{pmatrix}
  = \biggl\{
  \begin{pmatrix}
    \boldsymbol 0 \\ \u
  \end{pmatrix}
  \in \R^{k+n} \colon \u \in \Ker A^T 
  \biggr\} \,.
\end{equation}
In other words, the solution of \eqref{e.least-norm} exists and, if
not unique, the non-uniqueness affects only the Lagrange multiplier
$\boldsymbol \mu$; the least-norm vector $\v$ is uniquely determined.

\section{Singular value decomposition (SVD)}

Recall from your linear algebra class that if $S \in M(k \times k)$ is
a symmetric matrix, it has an orthogonal diagonalization.  More
precisely, there exists an orthogonal matrix\footnote{A matrix  $V \in
M(k \times k)$ is orthogonal if $V^{-1}=V^T$.  This is equivalent to
$V^TV=I$ (the columns of $V$ are an orthogonal basis of the column
space) which is equivalent to $VV^T=I$ (the rows of $V$ are an
orthogonal basis of the row space of $V$).} $V \in M(k \times k)$ such
that 
\begin{equation}
  \label{e.orth-diag}
  D = V^T S V
\end{equation}
is diagonal.  The orthogonal diagonalizability of symmetric matrices
can be used to construct a decomposition of \emph{any} matrix into a
product of orthogonal and diagonal matrices.  We proceed as follows.

Let $A \in M(n \times k)$.  Without loss of generality, we may assume
that $n \geq k$.  Then $S = A^TA \in M(k \times k)$ is symmetric and
has an orthogonal diagonalization of the form \eqref{e.orth-diag}.
Moreover, $S$ is positive semidefinite.\footnote{A matrix $S \in
M(k\times k)$ is positive definite if $\v^T S \v>0$ for all nonzero
$\v \in \R^k$; it is positive semi-definite if $\v^T S \v \geq 0$ for
all $\v \in \R^k$.  It is easy to see that all eigenvalues of a
positive definite matrix must be positive; all eigenvalues of a
positive semidefinite matrix must be nonnegative.}  By permuting the
columns of $V$, we can always achieve that the diagonal entries
$d_{ii}$ of $D$ appear in decreasing order so that, in particular,
$d_{11}, \dots, d_{mm} > 0$ and $d_{m+1,m+1}, \dots, d_{kk} = 0$ for
some $m \leq k$.  We now define
\begin{equation}
  \Sigma_0 = \sqrt D \equiv 
  \begin{pmatrix}
    \sqrt{d_{11}} & 0 & & & \cdots & 0\\
    0 & \ddots & & & & \vdots \\
    & & \sqrt{d_{mm}} & \\
    & & & 0 \\
    \vdots & & & & \ddots \\
    0 & \cdots & & & & 0 
  \end{pmatrix}
\end{equation}
and\footnote{$\Sigma^\dagger_0$ is called the \emph{pseudo-inverse} of
$\Sigma_0$.  It is characterized by
\[
  \Sigma^\dagger_0 \Sigma_0 = 
  \begin{pmatrix}
    I_m & 0 \\ 0 & 0
  \end{pmatrix} \,,
\]
where $I_m$ denotes the $m\times m$ identity matrix.}
\begin{equation}
  \Sigma_0^\dagger = 
  \begin{pmatrix}
    d_{11}^{-1/2} & 0 & & & \cdots & 0\\
    0 & \ddots & & & & \vdots \\
    & & d_{mm}^{-1/2} & \\
    & & & 0 \\
    \vdots & & & & \ddots \\
    0 & \cdots & & & & 0 
  \end{pmatrix} \,.
\end{equation}
and
\begin{equation}
  U_0 = A V \Sigma_0^\dagger \,.
\end{equation}
First, we note that the first $m$ columns of $U_0$ are an orthonormal
set in the column space, as
\begin{equation}
  U_0^T U_0
  = \Sigma_0^\dagger V^T A^T A V \Sigma_0^\dagger
  = \Sigma_0^\dagger V^T S V \Sigma_0^\dagger
  = \Sigma_0^\dagger D \Sigma_0^\dagger
  =
  \begin{pmatrix}
  I_m & 0 \\ 0 & 0
  \end{pmatrix} \,.
  \label{e.orthonormal}
\end{equation}
Second,
\begin{equation}
  U_0 \Sigma_0 V^T
  = A V \Sigma_0^\dagger \Sigma_0 V^T
  = A \,.
  \label{e.restrictedsvd}
\end{equation}
The last identity is obvious if $A$ is nonsingular so that $m=k$ and
$\Sigma_0^\dagger \Sigma_0 = I$.  In general, note that $A^TA=VDV^T$
implies
\begin{subequations}
\label{e.subspaces}
\begin{equation}
  \Ker A^\bot = \Span \{ \v_{1}, \dots, \v_m \}
\end{equation}
and
\begin{equation}
  \Ker A = \Span \{ \v_{m+1}, \dots, \v_{k} \} \,,
\end{equation}
where $\v_1, \dots, \v_k$ denote the columns of $V$.  The last step
of \eqref{e.restrictedsvd} is then easily verified by restricting to
$\Ker A$ and $\Ker A^\bot$, respectively.  Equation
\eqref{e.restrictedsvd} is referred to as the \emph{restricted
singular value decomposition} of $A$.

Now let $\u_1, \dots, \u_m$ denote the first $m$ columns of $U_0$.
Note that \eqref{e.orthonormal} shows that 
\begin{equation}
  \Range A = \Span \{ \u_{1}, \dots, \u_m \} \,.
\end{equation}
It is often convenient to extend the orthonormal set $\u_1, \dots,
\u_m$ to an orthonormal basis $\u_1, \dots, \u_n$ of $\R^n$, so that
\begin{equation}
  \Range A^\bot = \Span \{ \u_{m+1}, \dots, \u_{n} \} \,.  
\end{equation}
\end{subequations}
Then, setting
\begin{equation}
  U = 
  \begin{pmatrix}
  | & & | \\
  \u_1 & \cdots & \u_n \\
  | & & |
  \end{pmatrix} 
  \in M(n \times n)
  \quad \text{and} \quad
  \Sigma = 
  \begin{pmatrix}
    \Sigma_0 \\ 0
  \end{pmatrix}
  \in M(n \times k) \,,
\end{equation}
we can write
\begin{equation}
  \label{e.svd}
  A = U \Sigma V^T
\end{equation}
This is the \emph{singular value decomposition} of $A$; the diagonal
entries of $\Sigma$, denoted $\sigma_1, \dots, \sigma_k$, are the
\emph{singular values} of $A$, the vectors $\v_1, \dots, \v_k$ are the
\emph{right singular vectors}, and the vectors $\u_1, \dots, \u_n$ are
the \emph{left singular vectors}.

\paragraph{Variational characterization}

Recall, from Homework~1, that the principal unit eigenvector of a
symmetric matrix $S \in M(k \times k)$ maximizes the function
\begin{equation}
  f(\v) = \v^T S \v
\end{equation}
among all unit vectors $\v \in \R^k$.  The singular value
decomposition of a matrix $A \in M(n \times k)$ has a similar
variational characterization.  Consider the function
\begin{equation}
  f(\u,\v) = \u^T A \v \,.
\end{equation}
Then the maximum of $f$ over the set of unit vectors $\u \in \R^n$ and
unit vectors $\v \in \R^k$ is the largest singular value of $A$;
vectors $\u$ and $\v$ at which the maximum is attained can be taken as
the first left singular vector $\u_1$ and the first right singular
vector $\v_1$, respectively, in the notation of \eqref{e.subspaces}.

To establish our claim, we proceed as follows.  Our task can be
formulated as an optimization problem for $f$ on $\R^{n+k}$ with the
two scalar constraints
\begin{gather}
  g_1(\u,\v) \equiv \u^T \u - 1 = 0 
  \quad \text{and} \quad
  g_2(\u,\v) \equiv \v^T \v - 1 = 0 \,.
\end{gather}
Clearly,
\begin{subequations}
\label{e.gradients}
\begin{gather}
  \d f(\u,\v)[\delta \u, \delta \v]
  = \v^T A^T \delta \u + \u^T A \delta \v \,, \\
  \d g_1(\u,\v)[\delta \u, \delta \v] 
  = 2 \, \u^T \delta \u \,,
\intertext{and}
  \d g_2(\u,\v)[\delta \u, \delta \v]
  = 2 \, \v^T \delta \v \,.
\end{gather}
\end{subequations}
Then, the necessary condition for a local maximum is given by the
Lagrange multiplier principle \eqref{e.lagrange-2} which reads, after
substituting in \eqref{e.gradients} and rearranging terms,
\begin{gather}
  \bigl(
    \v^T A^T \delta \u - 2 \, \lambda_1 \, \u^T \delta \u
  \bigr)
  + 
  \bigl(
    \u^T A \delta \v - 2 \, \lambda_2 \, \v^T \delta \v
  \bigr) = 0
\end{gather}
which, since $\delta \u$ and $\delta \v$ can be chosen independently,
implies
\begin{gather}
  \v^T A^T = 2 \, \lambda_1 \, \u^T
  \quad \text{and} \quad
  \u^T A = 2 \, \lambda_2 \, \v^T \,.
\end{gather}
Right-multiplying the first equality with $\u$ and the second equality
with $\v$, we conclude that
\begin{gather}
  \u^T A \v = 2 \, \lambda_1 = 2 \, \lambda_2 \equiv \sigma \,.
\end{gather}
The existence of a solution is clear as we are maximizing a smooth
function over a compact constraint set; hence the necessary condition
must be satisfied somewhere.  Denote some such solution triple $\u_1$,
$\v_1$, and $\sigma_1$.\footnote{We are not claiming uniqueness.  In
fact, it can be shown that the solution is unique up to a choice of
sign if the singular values of $A$ are nondegenerate; in the
degenerate case, it is clear that we cannot expect uniqueness.}  Now,
restricting $\u$ and $\v$ to the respective orthogonal complements of
the vectors already established, we can iterate this argument $k$
times to obtain a full set of ``singular values'' $\sigma_1, \dots,
\sigma_k$, ``right singular vectors'' $\v_1, \dots, \v_k$ and,
completing the orthogonal set to a basis, ``left singular vectors''
$\u_1, \dots, \u_n$.

Two claims remain to be verified.  First, that the result of the above
construction is indeed indeed a singular value decomposition.  Second,
that the largest singular value is indeed the global (not just some
local) maximum of $f$ as claimed.  The verification of both statements
is not difficult and shall be left as an exercise.


\paragraph{Matrix norms}

Given any norm $\lVert \, \cdot \, \rVert$ on $\R^n$, we can define
the norm of a matrix $A = M(n \times k)$ as the smallest number
$\lVert A \rVert$ such that
\begin{equation}
  \lVert A \x \rVert \leq \lVert A \rVert \, \lVert \x \rVert
  \label{e.prodnorm}
\end{equation}
holds true.  In other words,
\begin{equation}
  \lVert A \rVert = \max_{\x \neq \boldsymbol 0}
  \frac{\lVert A \x \rVert}{\lVert \x \rVert} \,.
\end{equation}
This definition generalizes to operators between Banach spaces and is
often referred to as the \emph{operator norm}.  Note that, if the
coefficients of $A$ are denoted $a_{ij}$,
\begin{equation}
  \lVert A \rVert_F^2
  = \sum_{\substack{i \in 1,\dots, n \\ j \in 1, \dots, k}}
    \lvert a_{ij} \rvert^2
\end{equation}
defines a norm, called the \emph{Frobenius norm}, which is is not an
operator norm.  

Consider, in particular, the operator norm induced by the Euclidean
norm on $\R^n$, which we generally assume in these notes.  Then
$\lVert \x \rVert^2 = \x^T\x$, so that, using the singular value
decomposition,
\begin{align}
  \lVert A \rVert^2
  & = \max_{\x \neq \boldsymbol 0}
      \frac{\x^T V \Sigma^T U^T U \Sigma V^T \x}{\x^T \x} 
      \notag \\
  & = \max_{\x \neq \boldsymbol 0}
      \frac{\x^T V \Sigma^T \Sigma V^T \x}{\x^T VV^T \x}
      \notag \\
  & = \max_{\y \neq \boldsymbol 0}
      \frac{\y^T \Sigma^T \Sigma \y}{\y^T \y} 
      \notag \\
  & = \sigma_1^2
  \label{e.svd-2norm}
\end{align}
In other words, the matrix norm induced by the Euclidean metric equals
the largest singular value of the matrix.

In contrast, it is not difficult to show that
\begin{equation}
  \lVert A \rVert_F^2
  = \sum_{i=1}^k \sigma_i^2 \,.
\end{equation}

\paragraph{Condition numbers}

Suppose, for simplicity, that $A \in M(n\times n)$ is a nonsingular
matrix.  We would like to characterize the stability of the solution
$\x$ to the system of linear equations
\begin{equation}
  A \x = \b 
\end{equation}
under perturbation of the right hand side.  By linearity, $\Delta \b =
A \Delta \x$, where $\Delta \b$ denotes an \emph{absolute} change of
the vector $\b$, and $\Delta \x$ the corresponding absolute change in
$\x$.  Thus, by \eqref{e.prodnorm}, we obtain $\lVert \Delta \x \rVert
\leq \lVert A^{-1} \rVert \, \lVert \Delta \b \rVert$ and $\lVert \b
\rVert \leq \rVert A \rVert \, \lVert \x \rVert$, so that
\begin{equation}
  \frac{\lVert \Delta \x \rVert}{\lVert \x \rVert}
  \leq \lVert A \rVert \, \lVert A^{-1} \rVert \,
  \frac{\lVert \Delta \b \rVert}{\lVert \b \rVert} \,.
\end{equation}
The worst-case amplification factor from a \emph{relative} change in
$\b$ to the relative change in $\x$,
\begin{equation}
  \kappa(A) = \lVert A \rVert \, \lVert A^{-1} \rVert \,,
\end{equation}
is called the \emph{condition number} of $A$.  Due to
\eqref{e.svd-2norm},
\begin{equation}
  \kappa(A) = \frac{\sigma_1}{\sigma_n} \,.
\end{equation}
This notion of condition number generalizes naturally to non-square
matrices $A$ of full rank.  When the condition number is large, the
solution $\x$ becomes very sensitive to error (e.g.\ measurement
error) in the data $\b$; we speak of an \emph{ill-conditioned}
problem.

\section{Approximate solution of ill-conditioned linear equations}

Let us consider the following prototypical version of an
ill-conditioned system of linear equations.  Our presentation closely
follows \cite{V}.  Let $\x \in \R^n$ represent the true pixel values
of a, for simplicity, one-dimensional image.  The image is taken by a
device which, due to lens imperfections, say, records a blurred set of
pixel values $\b \in \R^n$.  We assume a simple linear model for the
device imperfections, namely that
\begin{equation}
  b_i = \sum_{j=1}^n a_{i-j} \, x_j \,,
  \label{e.convolution}
\end{equation}
where $b_1, \dots, b_n$ and $x_1,\dots,x_n$ denote the components of
$\b$ and $\x$, respectively.  This type of expression is referred to
as a \emph{discrete convolution} with \emph{kernel} $a_k$; typical is
a bell-like shape of its graph, e.g.,
\begin{equation}
  a_k = c \, \e^{-\gamma k^2}
\end{equation}
for some given constants $c$ and $\gamma$.  Equation
\eqref{e.convolution} can be written as a system of linear equations
\begin{equation}
  \b = A \x
  \label{e.linear2}
\end{equation}
with $n \times n$ matrix
\begin{equation}
  A = 
  \begin{pmatrix}
    a_0  & a_{-1} & \cdots \\
    a_1 & \ddots & \ddots \\
    \vdots & \ddots & & \ddots  & \vdots\\
    & & \ddots & \ddots & a_{-1} \\
    & & \cdots & a_1 & a_0
  \end{pmatrix} \,.
\end{equation}
Matrices of this type are known as \emph{Toeplitz matrices}.  They are
typically very ill-conditioned\footnote{In fact, \eqref{e.convolution}
can be interpreted as a discrete approximation of the integral
equation
\[
  b(x) = \int_a^b k(x-y) \, f(y) \, \d y
\]
where, with a smooth kernel function $k$, the right hand side defines
a compact, hence non-invertible operator on a Hilbert space.  Such
problem are referred to as \emph{ill-posed}.  Direct discretization of
an ill-posed continuum problem always leads to an ill-conditioned
finite dimensional system.  In this sense, the ill-conditioning in our
example is a generic phenomenon.}, so that, even in the case that $A$
is non-singular, the direct solution of \eqref{e.linear2} may be
practically meaningless.

\paragraph{Regularization}

The ill conditioning of \eqref{e.linear2} can be understood in terms of
the singular values of $A$ as follows.  Suppose $\x_\true$ refers to
the true pixel values of the image; due to measurement error,
\begin{equation}
  A \x_\true = \b + \Eta \,,
\end{equation}
where $\Eta$ is the measurement error vector with norm $\lVert \Eta
\rVert > 0$.  For simplicity, we suppose that $A \in M(n\times n)$ is
nonsingular.  Thus, when attempting to compute the pixel values by
direct inversion of $AS$, we obtain
\begin{equation}
  \x_{\operatorname{naive}}
  = A^{-1} \b
  = \x_\true + V \Sigma^{-1} U^T \Eta \,.
  \label{e.naive}
\end{equation}
The ill-conditioning manifests itself in the occurrence of small
singular values for $A$, so that $\Sigma^{-1}$ in \eqref{e.naive} has
small denominators---some components of $\Eta$ with respect to the
orthonormal basis $\{ \v_i \}$ are being multiplied with large gains.
As a result, we cannot expect $\x_{\operatorname{naive}}$ and
$\x_\true$ to be close.

Since the problem is due to small denominators, let us define a
filtered solution
\begin{equation}
  \x_\alpha = V \Sigma_\alpha^{-1} U^T \b 
  \label{e.xalpha}
\end{equation}
where $\Sigma_\alpha^{-1}$ is the filtered inverse of $\Sigma$,
defined by
\begin{equation}
  \Sigma_\alpha^{-1} = \diag \{w_\alpha(\sigma_i^2) \sigma_i^{-1}
  \colon i=1,\dots,n \}
\end{equation}
with $w(\sigma^2)$ converging to zero at least linearly as $\sigma \to
0$ to avoid unbounded singular values of the filtered inverse in this
limit.  Among the possible choices for $w_\alpha$, the following are
particularly natural.

\paragraph{Truncation}

The simplest choice is the cut-off filter
\begin{equation}
  w_\alpha(\sigma^2) 
  = \begin{cases}
      1 & \text{if } \sigma^2 \geq \alpha \\
      0 & \text{if } \sigma^2 < \alpha \,.
    \end{cases}
\end{equation}

\paragraph{Tychonov regularization}

With the choice 
\begin{equation}
  w_\alpha(\sigma^2) 
  = \frac{\sigma^2}{\sigma^2 + \alpha} \,,
\end{equation}
we obtain
\begin{equation}
  w_\alpha(\sigma^2) \sigma^{-1} 
  = \frac{\sigma}{\sigma^2 + \alpha}
\end{equation}
or
\begin{equation}
  \Sigma_\alpha^{-1} 
  = (\Sigma^2 + \alpha I)^{-1} \Sigma \,,
\end{equation}
thereby providing a smooth transition between the filtered and
unfiltered regime.  Substituting this expression into
\eqref{e.xalpha}, we obtain
\begin{align}
  \x_\alpha 
  & = V (\Sigma^2 + \alpha I)^{-1} \Sigma U^T \b
      \notag \\
  & = V (\Sigma^2 + \alpha I)^{-1} V^T V \Sigma U^T \b
      \notag \\
  & = (V (\Sigma^2 + \alpha I) V^T)^{-1} V \Sigma U^T \b
      \notag \\
  & = (A^TA + \alpha I)^{-1} A^T \b \,.
\end{align}
Note that this expression formally resembles the normal equations
\eqref{e.v} of the linear least squares problem, suggesting a close
relationship between the two problems.  This expectation will be
substantiated later when we describe a variational formulation of the
Tychonov regularization.

\paragraph{Error analysis}

We provide a deterministic error analysis, i.e., we assume that $\Eta$
is a fixed measurement error vector of known magnitude.  We can then
define
\begin{align}
  \E_\alpha 
  & = \x_\true - \x_\alpha
      \notag \\
  & = \x_\true - V \Sigma_\alpha^{-1} U^T \b 
      \notag \\
  & = V I V^T \x_\true - V \Sigma_\alpha^{-1} U^T (A \x_\true - \Eta)
      \notag \\
  & = V (I - \Sigma_\alpha^{-1} \Sigma) V^T \x_\true - 
      V \Sigma_\alpha^{-1} U^T \Eta
      \notag \\
  & \equiv \E_\trunc + \E_\noise
  \label{e.ealpha}
\end{align}
Note that, in contrast to the noise term in \eqref{e.naive}, the
expression for $\E_\noise$ involves a filtered inverse, thus has
bounded gains.  We expect that $\lVert \E_\noise \rVert$ decreases
with increasing filter parameter $\alpha$, while the truncation error
$\lVert \E_\trunc \rVert$ increases; thus, there should be an optimal
filter parameter $\alpha$.  This is borne out by the following
estimation.  

First, due to \eqref{e.prodnorm} and the fact that the Euclidean
matrix norm of an orthogonal matrix equals $1$,
\begin{equation}
  \lVert \E_\noise \rVert
  \leq \lVert V \rVert \, 
       \lVert \Sigma_\alpha^{-1} \rVert \,
       \lVert U^T \rVert \, \lVert \Eta \rVert
  \leq \max_{\sigma > 0} \frac{w_\alpha(\sigma^2)}{\sigma} \, 
       \lVert \Eta \rVert
  \leq \frac1{\sqrt \alpha} \, \lVert \Eta \rVert \,.
\end{equation}
The last inequality above can be verified independently for the
truncation filter and for the Tychonov regularization.

Second, we seek to derive estimates on the truncation error.  Here we
are faced with a fundamental dilemma.  In order to conclude that
truncation error---the error committed by filtering---is small, we
need to ensure that the information which gets removed by the filter
is is some sense unimportant for the true solution.  In other words,
we need to provide additional \emph{a priori} information on the
problem.\footnote{In fact, in the absence of further information, we
can still conclude, directly from the definition of the $w_\alpha$
considered, that $\lVert \E_\trunc \rVert \to 0$ as $\alpha \to 0$.
This conclusion, however, is practically meaningless because we have
no control about the rate of convergence, so there is no way to ensure
that we are doing substantially better than the naive, direct
inversion \eqref{e.naive}.}  In a standard analysis of this problem,
one assumes the \emph{range condition}\footnote{In an infinite
dimensional Hilbert space setting, the range condition can be stated
concisely as
\[
  \x_\true \in \Range A^T \,.
\]
There, this condition is nontrivial.} 
\begin{equation}
  \x_\true = A^T \z
  \label{e.range}
\end{equation}
such that $\lVert \z \rVert$ is ``reasonable,'' meaning of the same
order of magnitude as $\lVert \x_\true \rVert$.  It is important to
realize that this condition, from the mathematical point of view, is
arbitrary.  Rather, it represents additional knowledge about the
specific problem providing additional restrictions on the the class of
vectors in which we seek a solution.  As such, it must be justified on
physical grounds which necessarily lie outside of the analysis
presented here.  

Proceeding with our estimation, we obtain, inserting the range
condition \eqref{e.range} into the expression for $\E_\trunc$,
\begin{align}
  \lVert \E_\trunc \rVert^2
  = \E_\trunc^T \, \E_\trunc^{\vphantom T}
  & = \z^T A V (I - \Sigma_\alpha^{-1} \Sigma) V^T
      V (I - \Sigma_\alpha^{-1} \Sigma) V^T A^T \z
      \notag \\
  & = \z^T U \, (I - \Sigma_\alpha^{-1} \Sigma)^2 \, \Sigma^2 \, 
      U^T \z
      \notag \\
  & \leq \max_{\sigma>0} (1-w_\alpha(\sigma^2))^2 \, \sigma^2 \, 
      \lVert \z \rVert^2
      \notag \\
  & \leq \alpha \, \lVert \z \rVert^2 \,.
\end{align}
Inserting the bounds for $\E_\trunc$ and $\E_\noise$ into
\eqref{e.ealpha}, we obtain
\begin{equation}
  \lVert \E_\alpha \rVert
  \leq \sqrt\alpha \, \lVert \z \rVert + 
       \frac1{\sqrt \alpha} \, \lVert \Eta \rVert \,.
  \label{e.ealpha-estimate}
\end{equation}
It is easy to show that the right hand side is minimized when the two
term balance, i.e., when
\begin{equation}
  \alpha = \frac{\lVert \Eta \rVert}{\lVert \z \rVert} 
\end{equation}
so that
\begin{equation}
  \lVert \E_\alpha \rVert 
  = 2 \, \sqrt{\lVert \Eta \rVert \, \lVert \z \rVert} \,.
  \label{e.apriori}
\end{equation}
We call this an \emph{a priori} estimate for the regularization
parameter since knowledge of the true solution $\x_\true$ is
required.  In the following, we discuss methods of how to choose the
regularization parameter in an optimal way knowing only the data.

\paragraph{Morozov discrepancy principle}

In contrast to what was done in \eqref{e.ealpha}, we can split the
error into two components as follows,\footnote{The definition
$\E_\alpha^* = A \E_\alpha$ seems arbitrary, but is motivated by the
need to avoid any inversion of the nonregularized $A$ as this
computation generally cannot be performed in a stable manner.}
\begin{equation}
  \E_\alpha^*
  \equiv A (\x_\true - \x_\alpha)
  = (\b - A \x_\alpha) + \Eta
  \equiv \E_\trunc^* + \E_\noise^* \,.
  \label{e.ealpha-alt}
\end{equation}
Since $\E_\trunc^*$ now only depends on the measured data $\b$, we can
compute an \emph{a posteriori} estimate of $\alpha$ by balancing the
norms of $\E_\trunc^*$ and $\E_\noise^*$ rather than those of
$\E_\trunc$ and $\E_\noise$ as in \eqref{e.ealpha-estimate}.  Notice
that 
\begin{equation}
  \E_\trunc^* = U (I - \Sigma \Sigma_\alpha^{-1}) U^T \b 
\end{equation}
so that
\begin{equation}
  \lVert \E_\trunc^* \lVert^2
  = \b^T U \, (I - \Sigma \Sigma_\alpha^{-1})^2 \, U^T \b
  \to
  \begin{cases}
    0 & \text{as } \alpha \to 0 \\
    \b^T \b & \text{as } \alpha \to \infty \,;
  \end{cases}
\end{equation}
the left hand expression is, moreover, smooth and monotonic.  Thus,
whenever $\lVert \b \rVert \geq \lVert \Eta \rVert$, the intermediate
value theorem guarantees the existence of an $\alpha_\opt$ such that
$\lVert \E_\trunc^*(\alpha_\opt) \rVert = \lVert \Eta \rVert$.  This
regularization parameter selection criterion is called the
\emph{Morozov discrepancy principle}.

\paragraph{Error analysis of the Morozov principle}

We now ask the question how the error $\E_\alpha$ behaves under the
Morozov discrepancy principle.  This question is only easy to answer
in the case of the Tychonov regularization where $\x_\alpha$ is known
to minimize the function 
\begin{equation}
  f_\alpha (\x) 
  = \lVert A \x - \b \rVert^2 + \alpha \, \lVert \x \rVert^2
  \label{e.variational}
\end{equation}
over all $\x \in \R^n$; see Homework~2.  Thus, in particular,
\begin{equation}
  f_\alpha(\x_\alpha) \leq f_\alpha(\x_\true)
\end{equation}
which we can write 
\begin{equation}
  \lVert A \x_\alpha - \b \rVert^2 
  + \alpha \, \lVert \x_\alpha \rVert^2
  \leq \lVert \Eta \rVert^2 + \alpha \, \lVert \x_\true \rVert^2 \,.
\end{equation}
When $\alpha$ is chosen according to the Morozov principle, $\lVert A
\x_\alpha - \b \rVert = \lVert \Eta \rVert$, so that the first terms
on left and right of this inequality drop out and we conclude that
\begin{equation}
  \lVert \x_\alpha \rVert \leq \lVert \x_\true \rVert \,.
\end{equation}
Then, 
\begin{align}
  \lVert \E_\alpha \rVert^2
    = \lVert \x_\true \rVert^2 - 2 \, \x_\alpha^T \x_\true
      + \lVert \x_\alpha \rVert^2 
      \notag \\
    \leq 2 \, \lVert \x_\true \rVert^2 - 2 \, \x_\alpha^T \x_\true 
    = 2 \, \E_\alpha^T \x_\true \,.
\end{align}
As before, we assume the range condition with $\x_\true = A^T \z$, so
that
\begin{align}
  \lVert \E_\alpha \rVert^2
  \leq 2 \, \E_\alpha^T A^T \z 
  = 2 \, (\E_\alpha^*)^T \z
  \leq 2 \, \lVert \E_\alpha^* \rVert \, \lVert \z \rVert \,.
  \label{e.error2}
\end{align}
When $\alpha$ is chosen according to the Morozov principle, $\lVert
\E_\alpha^* \rVert = 2 \, \lVert \Eta \rVert$.  Inserting this
estimate back into \eqref{e.error2} and taking the square root, we
find that
\begin{equation}
  \lVert \E_\alpha \rVert 
  = 2 \, \sqrt{\lVert \Eta \rVert \, \lVert \z \rVert}
\end{equation}
which coincides with the \emph{a priori} error estimate
\eqref{e.apriori}.

\paragraph{Final Remarks}

The range condition is crucial for getting good reconstructions.  In
the case when $A$ represents a discrete convolution or the inverse of
a discrete differential operator, the range condition qualitatively
expresses that $\x_\true$ consists of samples of a smooth function on
a grid.  A concise discussion is beyond the scope of these notes, but
the consequences of this statement are easily visible in simple
computational experiments.  In particular, the filtered inverse will
reconstruct or reproduce step functions only poorly---the edges will
get eroded.

In the variational characterization \eqref{e.variational}, the term
$\alpha \, \lVert \x \rVert^2$ can be seen as a ``penalty term'' which
might be replaced with very different types of penalties, not all of
which might have an easy description in terms of the singular value
decomposition.  One important example is to replace $\lVert \x
\rVert^2$ by the so-called \emph{total variation}.  The total
variation penalizes oscillatory behavior while allowing jumps.  The
book by Vogel \cite{V} gives more complete, yet relatively accessible
introduction to theory and computational methods for inverse problems.


\begin{thebibliography}{9}

% \bibitem{DT}
%   G.B. Dantzig, M.N. Thapa,
%   ``Linear Programming'',
%   Springer, New York, 1997.

% \bibitem{G}
%   D. Gale,
%   \textit{Linear Programming and the simplex method},
%   Notices of the AMS \textbf{54} (2007), 364--369. \\
%   \url{http://www.ams.org/notices/200703/fea-gale.pdf}

% \bibitem{P}
%   P. Pedregal, ``Introduction to optimization'', Springer, New York,
%   2004. 

\bibitem{B}
  \r{A}. Bj\"orck,
  ``Numerical Methods for Least Squares Problems,''
  SIAM, Philadelphia, 1996.

\bibitem{V}
  C. Vogel, 
  ``Computational Methods for Inverse Problems,''
   SIAM, Philadelphia, 2002.

\end{thebibliography}

\end{document}
