\documentclass[11pt]{article}
\usepackage{fullpage,amsmath,amsfonts}

\newcommand{\x}{\boldsymbol{x}}
\newcommand{\y}{\boldsymbol{y}}
\newcommand{\f}{\boldsymbol{f}}
\renewcommand{\v}{\boldsymbol{v}}
\newcommand{\w}{\boldsymbol{w}}
\renewcommand{\r}{\boldsymbol{r}}
\renewcommand{\a}{\boldsymbol{a}}
\renewcommand{\b}{\boldsymbol{b}}
\renewcommand{\c}{\boldsymbol{c}}
\renewcommand{\k}{\boldsymbol{k}}
\newcommand{\z}{\boldsymbol{z}}
\newcommand{\e}{\mathrm{e}}
\renewcommand{\i}{\mathrm{i}}
\renewcommand{\d}{\boldsymbol{d}}
\renewcommand{\Re}{\operatorname{Re}}
\renewcommand{\Im}{\operatorname{Im}}
\newcommand{\Range}{\operatorname{Range}}
\newcommand{\Ker}{\operatorname{Ker}}
\newcommand{\Span}{\operatorname{Span}}
\renewcommand{\dim}{\operatorname{dim}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}

\def\nl{\vphantom\vdots\\}

\begin{document}

\title{Numerical Methods II}
\author{Problem Set 6}
\date{due in class, April 28, 2004}
\maketitle

\begin{enumerate}

\item 
\begin{enumerate}
\item Let $T_k$ denote the Chebychev polynomial of order $k$ on the
interval $[-1,1]$,
\begin{equation}
  T_k(x) = \cos(k \, \arccos x) \,.
\end{equation}
Show that
\begin{equation*}
  T_k(x) = \tfrac12 \,
    \Bigl(
      x + \i \sqrt{1-x^2}
    \Bigr)^k + \tfrac12 \,
    \Bigl(
      x - \i \sqrt{1-x^2}
    \Bigr)^k \,.
\end{equation*}

\item Show that, with $\kappa = \lambda_{\operatorname{max}} /
\lambda_{\operatorname{min}}$,
\begin{equation*}
  T_k \biggl(
      \frac{\lambda_{\operatorname{max}}+\lambda_{\operatorname{min}}}% 
           {\lambda_{\operatorname{max}}-\lambda_{\operatorname{min}}}
    \biggr)^{-1}  
  \leq 2 \,
       \biggl(
         \frac{\sqrt \kappa - 1}{\sqrt \kappa + 1}
       \biggr)^k \,.
\end{equation*}

\end{enumerate}

\item (From \emph{An Introduction to the Conjugate Gradient Method
Without the Agonizing Pain} by J.R. Shewchuk.)
  
Suppose you wish to solve $A\x=\b$ for a symmetric, positive-definite
$n \times n$ matrix $A$.  Unfortunately, the trauma of your linear
algebra course has caused you to repress all memories of the Conjugate
Gradient algorithm.  Seeing you in distress, the Good Eigenfairy
materializes and grants you a list of $d$ distinct eigenvalues (but
not the eigenvectors) of $A$.  However, you do \emph{not} know how
many times each eigenvalue is repeated.

Clever person that you are, you mumbled the following algorithm in
your sleep this morning:
\begin{verbatim}
  Choose an arbitrary starting point x(1);
  for i = 1:d
    r(i) = b - A*x(i);
    Remove an arbitrary eigenvalue from the list and call it l(i);
    x(i+1) = x(i) + r(i)/l(i);
  end
\end{verbatim}
No eigenvalue is used twice; on termination, the list is empty.

\begin{enumerate}
\item Show that upon termination of this algorithm, $\x_{d+1}$ is the
solution to $A\x=\b$.

\item Although this routine finds the exact solution after $d$
  iterations, you would like each intermediate iterate $x_i$ to be as
  close to the solution as possible.  Give a \emph{crude} rule of
  thumb for how you should choose an eigenvalue from the list on each
  iteration.  (In other words, in what order should the eigenvalues be
  used?)

\item What could go terribly wrong with this algorithm if floating
  point roundoff occurs?

\item Give a rule of thumb for how you should choose an eigenvalue
  from the list on each iteration to prevent floating point roundoff
  error from escalating.

\emph{Hint:} The answer is not the same as that of question (b).

\end{enumerate}

\item Recall that the linear Conjugate Gradient method for minimizing
  the quadratic function $f(\x) = \tfrac12 \, \x^T A \x - \x^T \b$
can be written
\begin{align*}
  \r_k &= -\boldsymbol{\nabla} f(\x_k) \,, \\
  \d_k & = \begin{cases}
             \r_1 & \text{for } k=1 \\
             \r_k + \beta_k \, \d_{k-1} & \text{for } k \geq 2 \,,
             \text{ where } \displaystyle
             \beta_k = -\frac{\r_k^T A \d_{k-1}}{\d_{k-1}^T A
             \d_{k-1}} \,,
           \end{cases} \\
  \x_{k+1} & = \x_k + \alpha_k \, \d_k \,, \text{ where }
  \alpha_k = \frac{\d_k^T\r_k}{\d_k^TA\d_k} \,.
\end{align*}
Show that the following two expressions for $\beta_k$ are exact in the
linear case.
\begin{enumerate}
\item $\displaystyle \beta_k^{\text{FR}}
  = \frac{\r_{k}^T \r_k}{\r_{k-1}^T \r_{k-1}} \quad $
  (Fletcher--Reeves)

\item $\displaystyle \beta_k^{\text{PR}}
  = \frac{\r_k^T(\r_k-\r_{k-1})}{\r_{k-1}^T \r_{k-1}} \quad$
  (Polak--Ribi\`ere)
\end{enumerate}
\emph{Remark:} These two formulae are useful for nonlinear problems,
as they do not involve the Hessian matrix of $f$ explicitly.  In this
case, the two expressions may differ.  Fletcher--Reeves converges
provided the starting point is sufficiently close to the extremum.
Polak--Ribi\`ere often converges more quickly, but convergence is only
guaranteed with the choice $\beta_k = \max \{ 0, \beta_k^{\text{PR}}
\}$.

\item \textbf{Project:} Implement a nonlinear Conjugate Gradient
  method.  Note that the parameter $\alpha_k$ cannot be computed
  explicitly as in the linear case.  Instead, you need to implement a
  one-dimensional ``line search'' routine!

\item \textbf{Project:} Compare the convergence of Fletcher--Reeves
  vs.\ Polak--Ribi\`ere when using your nonlinear Conjugate Gradient
  algorithm to find the minimum of the function
\[
  f(x_1,x_2) = -\frac1{1+x_1^2 + (1 + 10 \, \sin^2 x_1) \, x_2^2} \,.
\]
with starting point $(10,10)$.

\emph{Extra credit:} Does periodic restart of the CG iteration improve
convergence? Can you get rid of the requirement to explicitly compute
the gradient of $f$ by using a Broyden-type argument?

\end{enumerate}
  



\end{document}
