\documentstyle[leqno,twoside,11pt,ltexproc,epsf]{article}

\begin{document}
\cleardoublepage
\pagestyle{myheadings}
\raggedbottom

% Floats
\setcounter{topnumber}{3}
\renewcommand{\topfraction}{1}

% Useful macros
\newlength{\boxwidth}
\setlength{\boxwidth}{\hsize}
\addtolength{\boxwidth}{-40 pt}
\def\squarebox#1{\hbox to #1{\hfill\vbox to #1{\vfill}}}
%\newcommand{\qed}{\hfill\hfill\vbox{\hrule\hbox{\vrule\squarebox
%   {.667em}\vrule}\hrule}\smallskip}
\newcommand{\dy}{\Delta}
\newcommand{\tr}{\mbox{tr}}
\newcommand{\boxer}[1]{\vspace{1pc}\hspace*{-1.2in} \begin{center} \fbox{\fbox{
     \begin{minipage}{\boxwidth} #1 \end{minipage}}} \end{center} 
     \vspace{1pc}}
\newcommand{\masbox}[1]{\begin{center} {\framebox[6.5in][c]{
     \begin{minipage}{\boxwidth} \vspace{0.1in} #1 \vspace{0.1in}
\end{minipage}}} \end{center}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\T}{^T\!}
\newcommand{\Yd}{{\dot Y}}
\newcommand{\Ydd}{{\ddot Y}}
\newcommand{\Yperp}{{Y_{\!\perp}}}
\newcommand{\Qd}{{\dot Q}}
\newcommand{\Qdd}{{\ddot Q}}
%\newtheorem{theorem}{Theorem}[section]
%\newcommand{\proof}{\noindent{\bf Proof }}
\def\myheading#1{\hline\multicolumn{3}{|l}{#1}}

% Derivatives
\def\pdx#1{{\partial #1\over\partial x}}
\def\pddx#1{{\partial^2 #1\over\partial x^2}}
\def\grad{\mathop{\rm grad}\nolimits}
\def\Hess{\mathop{\rm Hess}\nolimits}

% Real and Complex
\newcommand{\R}{{\bf R}}
\newcommand{\C}{{\bf C}}

% Boxed matrices
%\newdimen\sqbox \sqbox=48pt
%\newdimen\rcbox \rcbox=10pt

\def\boxit#1{\vbox{\hrule\hbox{\vrule\kern0pt
      \vbox{\kern0pt\hbox{#1}\kern0pt}\kern0pt\vrule}\hrule}}
\def\boxitop#1{\vtop{\hrule\hbox{\vrule\kern0pt
      \vbox{\kern0pt\hbox{#1}\kern0pt}\kern0pt\vrule}\hrule}}
\def\boxedQ{\vcenter{\boxit{\rule[-.4\sqbox]{0pt}{\sqbox}%
  \hbox to\sqbox{\hss$Q$\hss}}}}
\def\boxedY{\vcenter{\boxit{\rule[-.4\sqbox]{0pt}{\sqbox}%
  \hbox to\rcbox{\hss$Y$\hss}}}}

% Fix eqnarray's \def so that it uses the correct spacing (cf. latex.tex).
\makeatletter
\def\eqnarray{\stepcounter{equation}\let\@currentlabel=\theequation
  \global\@eqnswtrue
\global\@eqcnt\z@\tabskip\@centering\let\\=\@eqncr $$\halign to
\displaywidth\bgroup\@eqnsel\hskip\@centering
$\displaystyle\tabskip\z@{##}$&\global\@eqcnt\@ne
${{}##{}}$&\global\@eqcnt\tw@$\displaystyle\tabskip\z@{##}$\hfil
\tabskip\@centering&\llap{##}\tabskip\z@\cr}
% 12-point definitions
\setbox\strutbox=\hbox{\vrule height10pt depth4pt width\z@}%
\def\big#1{{\hbox{$\left#1\vbox to10\p@{}\right.\n@space$}}}%
\def\Big#1{{\hbox{$\left#1\vbox to13\p@{}\right.\n@space$}}}%
\def\bigg#1{{\hbox{$\left#1\vbox to17.5\p@{}\right.\n@space$}}}%
\def\Bigg#1{{\hbox{$\left#1\vbox to21\p@{}\right.\n@space$}}}%

% Add some handy plain TeX macros
\def\eqalign#1{\null\,\vcenter{\openup\jot\m@th
  \ialign{\strut\hfil$\displaystyle{##}$&$\displaystyle{{}##}$\hfil
\crcr#1\crcr}}\,}

\makeatother

\title{On Conjugate Gradient-Like Methods \\
 for Eigen-Like Problems}

\author{
Alan Edelman\thanks{Department of Mathematics Room 2-380,
Massachusetts Institute of Technology, Cambridge, MA 02139,
{\protect\verb+edelman@math.mit.edu+}, Supported by a fellowship from
the Alfred P. Sloan Foundation
and  NSF Grant 9404326-CCR.  }
\and
Steven T. Smith\thanks{M.I.T. Lincoln Laboratory, 244 Wood Street,
Lexington, MA 02173-9108, {\protect\verb+stsmith@ll.mit.edu+}.
This work was sponsored by ARPA under Air~Force contract
F19628-95-C-0002. Opinions, interpretations, conclusions, and
recommendations are those of the author and are not necessarily
endorsed by the United States Air~Force.
\newline 
\hspace*{0.2in} This paper may also appear in the 
journal {\it BIT}.
}
}
\date{}
\maketitle
\markboth{Edelman and Smith}{Conjugate Gradient-Like Methods for
Eigen-Like Problems}

\pagenumbering{arabic}


%\setlength{\baselineskip}{.5cm}



\begin{abstract}
Numerical analysts, physicists, and signal processing engineers have
proposed algorithms that might be called conjugate gradient for
problems associated with the computation of eigenvalues.  There are
many variations, mostly one eigenvalue at a time though sometimes
block algorithms are proposed.  Is there a correct ``conjugate
gradient'' algorithm for the eigenvalue problem? How are the
algorithms related amongst themselves and with other related
algorithms such as Lanczos, the Newton method, and the Rayleigh
quotient?
\end{abstract}


\section{Introduction}
If we minimize $y^T\!Ay$ on the unit sphere, perhaps with conjugate
gradient optimization, we compute the smallest eigenvalue of $A$
(assumed to be symmetric).  Our objective function is quadratic, so we
obtain the exact eigenvalue after $n$ steps of conjugate gradient.
The computation is mathematically equivalent to the Lanczos procedure.
At the $k$th step, we obtain the optimal answer in a $k$ dimensional
Krylov space.


The above paragraph may sound correct, but it is wrong!  The objective
function is quadratic, but only because of the constraint. The
Rayleigh quotient $r(y)=y^T\!Ay/y^T\!y$ can be minimized without
constraint, but it is not quadratic.  The conjugate gradient link with
the Lanczos process here is also dubious.  How tempting it is to hear
``eigenvalue'' juxtaposed with ``conjugate gradient'' and instantly,
almost by reflex, yell ``Lanczos.''  To add to the confusion,
conjugate gradient optimization of the Rayleigh quotient does compute
an answer in a Krylov space, but not the optimal answer.

This note, which complements and serves as a friendly introduction to
our recent paper \cite{us95}, sorts through this muddle.

\section{Conjugate Gradient: linear, nonlinear, and idealized}

Our first task is to eliminate any 
confusion among  the various conjugate
gradient algorithms with the Lanczos algorithm.
This confusion has appeared in the literature.
We  hope that a naming convention
will help distinguish the algorithms in the future:
\boxer{
{\sc Conjugate Gradient Algorithms}\\ \\
\begin{tabular}{ll}
LCG & {\bf Linear Conjugate Gradient:} \\
& Standard method for solving $Ax=b$, where $A$ is a symmetric \\
&positive definite matrix.  Found in all numerical linear algebra \\
& references. May be viewed as minimizing the objective function  \\
&   $f(x)=\frac{1}{2}x^T\!Ax - x^T\!b$. \\ \\
NCG & {\bf Nonlinear Conjugate Gradient:} \\
    & A well-known algorithm often used for unconstrained minimization. \\
    & Found in all numerical optimization texts.  May be used to \\
    & minimize very general objective functions.   Each step in \\
    & the algorithm requires the solution to a one dimensional \\
    & line minimization problem. \\ \\
ICG & {\bf Idealized Conjugate Gradient:} \\
    & A fictional algorithm introduced  here for purposes of exposition.  \\
    & We define ICG as  an algorithm whose $k$th iterate $x_k$    \\
    & satisfies $f(x_k)=\min_{ x \in {\cal K}_k}  f(x)$, where
     ${\cal K}_k$ denotes a  \\&  $k$ dimensional Krylov space. The 
     Krylov spaces are nested. \\
\end{tabular}
}


The three conjugate gradient algorithms are equivalent for the special
case of quadratic objective functions with positive definite Hessians.
To say this in a slightly different way, the LCG algorithm, already a
special case of the more general NCG algorithm, may be derived with
the ``one line minimization per iteration'' point of view, and
nevertheless, the algorithm has the ``global $k$-dimensional''
property of an ICG anyway.


To contrast, if the function is not quadratic with positive definite
Hessian, there is no expectation that the $k$th step of NCG will have
the ICG property of being the global minimum in a $k$-dimensional space.
The best that we hope for is that the ICG property holds approximately
as we approach the minimum.


Similar to ICG is $s$-step steepest descent. This algorithm minimizes
an objective function in the $s$-dimensional space spanned by the most
recent $s$ gradients.  If the gradients fall in a certain natural
Krylov space, then there is no distinction between ICG and $s$-step
steepest descent.

Block versions of NCG and ICG may be defined.  A block NCG replaces a
$1$-dimensional line search with a $p$-dimensional line search.  A
block ICG maximizes over block Krylov spaces of size $p$, $2p$, $3p$, etc.

Consider the case when the objective function is the Rayleigh
quotient, $f(y)=(y^T\!Ay)/(y^T\!y)$.
\begin{itemize}
\item  An NCG algorithm for Rayleigh quotient optimization
computes an answer in a Krylov space.  This answer is
not the  optimal choice from
that subspace. Therefore the NCG algorithm is not an  ICG
algorithm.  
The exact eigenvalue will not be obtained in $n$ steps.
Such an algorithm  is not equivalent to Lanczos.
\item  The Lanczos algorithm  (including the
computation of the tridiagonal's smallest eigenvalue) is an ICG
algorithm for the Rayleigh quotient.  Lanczos does compute the exact
eigenvalue in $n$ steps.  This algorithm is not an NCG algorithm
because it is not equivalent to an algorithm that performs line
searches.
\item
There is also the LCG link with Lanczos through
tridiagonalization. The residual vectors computed by conjugate
gradient are (up to normalization) Lanczos vectors that tridiagonalize
the matrix.  The link is through the tridiagonalization and not the
eigenvalue computation.  There is no link to the NCG algorithm for
eigenvalue computation.
\end{itemize}

Also we note that
\begin{itemize}
\item any algorithm that
requires that the symmetric matrix $A$ be positive definite is not a
correct NCG algorithm for the Rayleigh quotient.  Since all
derivatives of the objective function are invariant under shifts, so
should be the algorithm.
\end{itemize}
Understanding this last point clearly is the key sign that the reader
understands the distinction between NCG for the Rayleigh quotient and
LCG which is an optimization on quadratic functions.  LCG requires
positive definite quadratic optimization functions so that a minimum
exists.  The Rayleigh quotient always has a minimum for symmetric $A$;
no positive definite condition is needed.



The first choice in our design space for algorithms is to consider
whether to include the constraints or work with the Rayleigh quotient.
The second choice is whether to have a one eigenvalue algorithm or a
block algorithm.  We indicate these choices below:
\boxer{
{\sc Rayleigh Quotients} \\
\begin{tabular}{ll}
Unconstrained one eigenvalue Rayleigh quotient &
$r(y)=\frac{Y\T Ay}{y^Ty}$ \\
Constrained Rayleigh quotient &$ r(y)=Y\T Ay \ (y^Ty=1)$ \\ 
\\
Unconstrained block Rayleigh quotient & $R(Y)=\tr (Y\T AY)(Y^TY)^{-1}$ \\
Constrained block Rayleigh quotient & $R(Y)=\tr Y\T AY \quad (Y^TY=I_p)$
\end{tabular}}

In the above box $Y$ denotes an $n \times p$ matrix.
Of the four choices,
only the constrained $r(y)=Y\T Ay$ is
(generically)  non-degenerate in
that the global minima are points.




The third choice in our design space for 
an NCG
eigenvalue  algorithm  is how to 
pick the new search direction.  The various
options are
\boxer{
{\sc NCG search direction choices} \\
\begin{tabular}{ll}
FR & The Fletcher-Reeves approach \\
PR & The Polak-Ribi\`{e}re approach \\
CA & Conjugacy through the matrix $A$ \\
CSH & Conjugacy through the singular free space Hessian $H$ \\
CCH & Conjugacy through a non-singular constrained Hessian $H$ \\
\end{tabular}
}

All of these choices, except for CA, are reasonable.
There is no sound mathematical justification for CA, though
it has been proposed in the literature.

\section{Comparison between NCG and Lanczos}

Since Lanczos is an ICG algorithm for the
Rayleigh quotient, why would anyone consider
the construction of an NCG algorithm which
is guaranteed to give a worse answer for
the same number of matrix vector multiplies?
Ignoring the extra storage that may be
needed in the Lanczos algorithm, the
more practical point is that in many algorithms
it is not the eigenvalues of a matrix
that are desired at all.

In the applications of interest, a matrix
may be changing in  ``space'' or time.
What we want from an algorithm is
to be able to track the changes.
In many such applications, the problem
sizes are huge so computational
efficiency is of utmost importance.
We have two particular examples in mind.
In the Local Density Approximation (LDA)
to Schr\"{o}dinger's equation, the function
to be minimized is not the Rayleigh
quotient at all, but rather the LDA energy.
In many applications, the energy function
has the same degeneracy as the Rayleigh
quotient  in that it depends only on the
span of the columns of $Y$ rather than all of $Y$.
Lanczos would not directly apply unless one
wants to pretend that the function 
really is the Rayleigh quotient for a few
iterations.  Such pretense has lead to
instabilities in the past.

Another class of problems arise in signal
processing where the matrix is changing
in time.  Once again, Lanczos does not
directly apply, unless one wants
to pretend that the matrix is constant
for a few intervals.


\section{A History of CG for eigen-like problems}

The purpose of this section is to demonstrate
the need for a unifying theory of conjugate
gradient algorithms for the eigenproblem
by giving an annotated  chronology of the 
references.  We note that algorithms
have been proposed by researchers in a
number of disciplines, so it seems likely
that researchers in one field may well
be unaware of those from another.



In order to keep the chronology  focused, we only
included papers that linked conjugate
gradient to eigenvalue or eigenvalue-like
problems.  We found it convenient to divide the papers
into five broad categories:

\begin{itemize}
\item {\bf Single Eigenvalue (NCG) Algorithms} 
This refers to the
one eigenvalue at a time algorithm that  was originally
proposed by Bradbury and Fletcher \cite{bradbury66a}.
This algorithm minimizes $r(y)$.  Other eigenvalues
may be obtained by deflation.  The historical
development covers a number of points in the design
space.  

\item {\bf Block (NCG)  Algorithms}
These are attempts at computing multiple eigenvalues
simultaneously using a conjugate gradient style
algorithm.  The objective function is $R(Y)$.
It is our point of view that none of these algorithms
are a true conjugate gradient algorithm, though
many important features of a true algorithm may
be recognized in either explicit or rudimentary form
in the proposed algorithms. The block case exhibits
difficulties not found in the single eigenvalue case,
because not only is there the orthogonality constraint
$Y^TY=I$, but also there is the ``Grassmann equivalence''
$R(Y)=R(YQ)$ for orthogonal matrices $Q$.
This degeneracy may well be overcome without
differential geometry, but one effective and beautiful
theoretical approach towards working with
these equivalence classes is the total space/base space
viewpoint from  differential geometry.

\item  {\bf ICG: The Lanczos Link}
Here our goal was not to mention papers that were connected
to the Lanczos algorithm, but papers that looked at the
Lanczos algorithm as an optimization algorithm with
something of  a conjugate gradient flavor.   As mentioned
earlier, the famous link between Lanczos and LCG  is
of no relevance here and is not included.

\item {\bf  Differential Geometry Viewpoint}
These papers make the link between numerical linear
algebra, optimization, and differential geometry.
They give a new theoretical viewpoint on the
algorithms discussed in the other areas.
This viewpoints shows that constrained algorithms may 
be implemented without explicit constraints.
Ultimately, they answer the question of what
it {\it means} to do conjugate gradient minimization
for the eigenvalue problem.



\item {\bf Application I: LDA Schr\"{o}dinger's Equation}
Rightly or wrongly, conjugate gradient has become an extremely popular
new method for scientists working with the local density approximation
to Schr\"{o}dinger's Equation.  Here the problem is akin to minimizing
the block Rayleigh quotient, but it is in fact more complicated.  Many
researchers are currently trying variations on the same theme.  We
only need to list the few that we were most aware of.  One point of
view is that they are computing eigenvalues of a matrix that is in
some sense varying with space.

\item {\bf Application II: Adaptive Signal Processing}
Here the algorithm is an eigenvalue algorithm, but
the matrix may be thought of as changing with time.
The goal is to do subspace tracking.  All the papers in
this category have already been listed in a previous
category, but we thought it worthwhile to collect
these papers under one heading as well.
\end{itemize}


\def\tin#1{\parbox[t]{2.7in}{\renewcommand{\baselinestretch}{.1} \small {#1}
}\renewcommand{\baselinestretch}{.1} }

\newcommand{\hello}{ \protect \rule{0mm}{6mm}    }

\newpage

\begin{center}
{\large Chronology of Algorithms} \\
Key to search direction choices: FR=Fletcher-Reeves,
PR=Polak-Ribi\`{e}re, CA=Conjugacy through $A$,
CSH=Conjugacy through singular Hessian $H$, \\  CCH=
Conjugacy through constrained Hessian $H$.
\end{center}


\newcommand{\xxx}{  \begin{center}  \noindent
  \begin{tabular}{|l@{\hello}ll@{\hspace{-0in}}l|}
}

\newcommand{\yyy}{ \\
&  & 
\dimen0=\hsize\advance\dimen0by-3.7in\rule{\dimen0}{0in}
 &  \\
\hline \end{tabular} \end{center}}


\xxx
\myheading{
1. {\bf Single Eigenvalue Algorithms}} &
\\
& 1966 & Bradbury and Fletcher  \cite{bradbury66a} & 
\tin{ Notes  degeneracy  of  Rayleigh quotient.
  Proposes projection to  $\infty$- and 2-norm
unit spheres. FR.  }\\  
&1969 & Fox and Kapoor  \cite{foxkapoor69} & 
\tin{ Application in Structural Dynamics.
Downplays importance of the  constraint. FR.} \\
 &  1969   & Fried \cite{fried69} & \tin{ 
Application in Structural Dynamics. FR.} \\
& 1971 & Andersson \cite{andersson71} &
\tin{Compares two norms mentioned above \cite{bradbury66a}.} \\
&1971 & Geradin  \cite{geradin71a} & \tin{
Unconstrained CSH.
}\\
&  1972   & Fried \cite{fried72} &
\tin{Unconstrained FR.}
\\
& 1974 & Ruhe \cite{ruhe74} & 
\tin{Systematically compares above CG algorithms with
other non-CG  algorithms.  For CG, prefers
unconstrained approach in
\cite{foxkapoor69} and \cite{geradin71a}.
Observes spectral influence on convergence.
}\\
%& 1978 & Townsend and Johnson \cite{townsend78}. & \\ 
&1984 & Fuhrmann and Liu  \cite{fuhrmann84} &
\tin{Precursor to intrinsic geometry: Searches
along geodesics.  FR.}\\
& 1986 & Perdon and Gambolati \cite{perdon86} & 
\tin{CA with a proposed  preconditioner.} \\
&1986 & Chen, Sarkar, et al.   \cite{chen86} &
\tin{ Rederives known algorithms.   CA.
} \\
& 1987 & Haimi-Cohen and Cohen \cite{cohen87} \hspace*{3mm}  &
\tin{Rederives results such as those in \cite{bradbury66a}.} \\
&1989 & Yang, Sarkar, and Arvas \cite{yang89} &
\tin{
Considers effect of normalization.  Concludes
that their 1986 proposal is not competitive
with CSH or other choices.
FR or PR style.
\cite{townsend78}.} 
\yyy


\xxx
\myheading{
2. {\bf Block Algorithms:}} &
\\
&1971 & Als\'{e}n    \cite{alsen71a}
& \tin{  Assumes $A$ positive definite.
Exact line maximization (not minimization)
replaced by one step of orthogonal iteration.
Corresponding minimization
algorithm  would require  orthogonal inverse
iteration. FR.
} \\
%&1973 & Shavitt, et al. \cite{shavitt73} \\
&1982 & Sameh and Wisniewski  \cite{sameh82a} & 
\tin{
Not exactly a CG method
but very similar. 
Constrained objective function with unconstrained
block line searches.  Section 2.3 suggests
block line minimization
on unconstrained function as a least
squares problem to be solved with LCG.
Makes the Lanczos style link through simultaneous
iteration.
} \\ 
&1995 & Fu and Dowling \cite{fu95}&
\tin{Optimization
of unconstrained $R(Y)$.
Block search directions
arranged column by column.
Projects back
to constraint surface with Gram Schmidt. CA.} 
\yyy

\xxx
\myheading{
3. {\bf ICG: The Lanczos Link:}} & \\
&1951 & Karush \cite{karush51} & \tin{
An ICG algorithm (with restarts) for the Rayleigh quotient.
Not explicitly identified with Lanczos or CG.}
\\
&1974 & Cullum and Donath \cite{cullum74b,cullum74a} &
\tin{
Identifies block Lanczos as a block ICG.}
\\
& 1978 & Cullum \cite{cullum78}
& \tin{Shows that block NCG on the Rayleigh
quotient computes a (generally non-optimal)
 answer in the block
Krylov space over which the  ICG would finds the  minimum.} \\
&1985 & Cullum and Willoughby  \cite{cullum85}  &
\tin{Chapter 7 summarizes  optimization interpretation.
CG used at times  to mean NCG and other times ICG.
} 
%&1993 & omit veselic?  \cite{kovac93a} 
\yyy



\xxx
\myheading{
4. {\bf  Differential Geometry Approach:}} &
\\
& 1993 & Smith \cite{smith93a,smith94a} & 
\tin{Introduces Differential Geometry viewpoint.
Hessian is replaced with second covariant derivative.
Line searches replaced with geodesic searches.
} \\
&1995 & Edelman, Arias, Smith  \cite{us95}  &
\tin{ Works out   details of the Grassmann manifold
and Stiefel manifold approaches, including the
total space-base space formulation, parallel
transport, differing metrics, and the linear
algebra links.
} 
\yyy

\xxx
\myheading{
5. {\bf LDA Schr\"{o}dinger's Equation} }&
\\
& 1989 & Gillan \cite{gillan89a} & 
\tin{Projects onto constraint space, search
in tangent space, suggests NCG preconditioning, FR.} \\
& 1989 & \v{S}tich, Car, et al. \cite{stich89a} &
\tin{Minimizes Rayleigh quotient only (not LDA energy) for simplicity.} \\
& 1989 & Teter, Payne, and Allan  \cite{teter89a} & \\
& 1992 & Arias \cite{arias92a} & \tin{Forms
unconstrained analogue of Rayleigh quotient for LDA.}\\
& 1992 & Payne, Teter, Allan et al. \cite{ptaaj} &
\tin{Survey of minimization in LDA.} \\
& 1993 & Kresse and Hafner \cite{kresse93} &
\tin{Molecular dynamics of liquid metals.} \\
& 1994 & Sung, Kawai, and Weare \cite{sung94} &
\tin{Computes  structure of Lithium.}
\yyy


\xxx
\myheading{
6. {\bf Adaptive Signal Processing}} &
\\
&1986 & Chen, Sarkar, et al.   \cite{chen86} & \tin{See above.} \\
&1989 & Yang, Sarkar, and Arvas \cite{yang89} & \tin{See above.} \\
&1995 & Fu and Dowling \cite{fu95}& \tin{See above.}
\yyy


\vspace{0.2in}

\setlength{\textwidth}{6.5in}

There are many closely related
ideas such as steepest descents or coordinate
relaxation that are not discussed here.
It is instructive to consider the subtle variations
between similar algorithms, as well as the
key differences between algorithms that sound
the same, but are quite different.



Almost all
the algorithmic 
ideas carry over to the generalized eigenvalue
problem, as was recognized by most of the authors.
For purposes of exposition, we will only discuss
the ordinary eigenvalue problem in this paper.
We have chosen to confine our discussion to
conjugate gradient methods for the general eigenproblem.
The larger history of conjugate gradient methods
may be found in the survey by Golub and O'Leary
\cite{golub89}. The search direction choice(s)
appear in the annotations.




\section{The Differential Geometry Viewpoint for NCG}

In a recent paper \cite{us95}, we proposed an
algorithm to perform minimization on the
block Rayleigh quotient which has appealing
theoretical properties.
Our algorithm   takes into account
both the constraints $Y^TY=I_p$ 
and the degeneracy of the function ($R(Y)=R(YQ)$ for orthogonal
$p \times p$ matrices $Q$) in a clean manner.
Following the lead proposed by Smith \cite{smith93a,smith94a},
we derived an algorithm based on intrinsic
geometry that is  practical when expressed in  extrinsic
coordinates.


\begin{figure}[p]
\centerline{\epsfysize=2in\epsfbox{curved.ps}}
\caption{Conjugate gradient on Riemannian manifolds.}
\end{figure}

\begin{figure}[p]
\centerline{\epsfysize=4in\epsfbox{sfig.ps}}
\caption{Riemannian conjugate gradient, unconstrained 
block Polak-Ribi\`{e}re, and the $A$-conjugacy algorithms
compared.}
\end{figure}


Like French-English false cognates, there are
a few terms used  both in optimization and
differential geometry with somewhat different
meanings.  This may cause confusion to readers of
our paper \cite{us95} who are specialists in
optimization.


\def\descript#1{\vtop{\hsize=4in\raggedright\noindent#1}}

\boxer{
{\bf  Optimization vs.\  Differential Geometry } \\
\\
\begin{tabular}{lll}
{\bf  Metric} \\
Optimization:&  \descript{Sometimes the
inner product  defined by the  Hessian or its inverse.}
 \\
Diff.~Geom.:&
\descript{Any positive definite inner product defined on a tangent space
from which all geometrical information about a space may be derived.}
 \\  \\
{\bf  Curvature} & \\
Optimization:&  \descript{Hessian. Usage  refers to non-linearity 
of the graph of a function with
nonzero Hessian.} \\
Diff.~Geom.:& \descript{A rank four tensor that refers roughly to the
non-flatness of a space.}
\end{tabular}
}


The Grassmann  algorithm  \cite{us95}  does not suffer from any of the
deficiencies of the other block algorithms
for the eigenvalue problem.  There
is no requirement that the matrix be positive
definite, there is no  theoretical difficulty
related to the singular Hessian,  and the
algorithm converges quadratically.

The key components of the Grassmann algorithm are the
total space/base space point of view, the following
of geodesics, and the parallel transportation
of tangent vectors.  Below we plot a schematic
figure that describes conjugate gradient on
a Riemannian manifold:

In Figure~2, we show convergence curves in exact arithmetic for three
block algorithms.  The $A$-conjugacy algorithm is clearly inferior and
not worth considering.  The Riemannian algorithm and our own
unconstrained version of the block Polak-Ribi\`{e}re algorithm, which
approximates the Riemannian CG algorithm up to second order and is
less computationally intensive than the full-blown Riemannian
algorithm, are quadratically convergent.


\section{Beyond Conjugate Gradient}


Potential users of a conjugate gradient algorithm
may well consider the viewpoint  expressed
by Jorge Nocedal \cite{nocedal91}:
\begin{quotation}
The recent development of limited memory and
discrete Newton methods have narrowed the class
of problems for which conjugate gradient methods
are recommended.  Nevertheless, in my view,
conjugate gradient methods are still the
best choice for solving very large problems with
relatively inexpensive objective functions.
They can also be more suitable than limited memory
methods on several types of multiprocessor computers.
\end{quotation}

Though we believe that we have derived the
correct NCG algorithm for functions such
as the block Rayleigh quotient, it is
very possible that the best algorithms
for applications in the LDA community
and the signal processing community
may well be a Newton iteration rather than
an algorithm in the conjugate gradient
family. Rightly or wrongly,
the  CG algorithms may yet  remain popular 
for the largest problems because of the
simplicity of programming and limited 
memory requirements.


It is important to understand the relationship between
preconditioned conjugate gradient and Newton's method.
We sometimes consider the ideal preconditioner
to be one that is easily computed, and yet closely
resembles the inverse of the Hessian.  Multiplying
the gradient by the exact inverse of the Hessian
is exactly the Newton method. Therefore the Newton
method is equivalent to preconditioned conjugate
gradient without any use of a previous search direction.
One expects that
taking advantage of Hessian information, if convenient,
would lead to superior convergence.


We now explore Newton's method for invariant subspace computation more
closely.  Suppose (for simplicity only) that $A$ is symmetric.  We
would then wish to find matrices $Y$ and $B$ such that
\begin{equation} AY-YB=0. \end{equation} The above system is degenerate;
constraints of some form or another need to be imposed.  One approach
is the general \begin{equation} \mbox{affine constraint } Z^TY=I
\end{equation} introduced by Chatelin \cite{chatelin84,chatelin93} in
which case $B=Z^TAY$.  We observe that the algorithms proposed by
Dongarra, Moler, and Wilkinson \cite{dongarra83} and Stewart
\cite{stewart73}.  represent two special choices for the affine
constraint matrix $Z$.  In the Dongarra et al.\ case, $Z$ may be
obtained by inverting and transposing an arbitrary $p \times p$ minor
of the $n \times p$ matrix $Y$.  In Moler's Matlab notation,\ \
\verb+Z=zeros(n,p); Z(r,:)=inv(Y(r,:))'+,\ \ where \verb+r+ denotes a
$p$-vector of row indices.  For Stewart, $Z=Y(Y^TY)^{-1}$.  Demmel
\cite[Section3]{demmel87} originally observed that the three
algorithms were related.  Being linear, the affine constraint allows a
direct application of Newton's method without any difficulty.

An alternative choice for the constraint is the \begin{equation}
\mbox{orthogonality constraint } Y^TY=I. \end{equation} Newton's
method with this constraint is slightly more complicated because of
the nonlinearity.  The value of $B$ is now $Y\T AY$.  However the
equation $AY-Y(Y\T AY)=0$, even with the orthogonality constraints, is
degenerate.  This makes this problem even more difficult to handle.

The geometric approach to the Grassmann manifold gives an approach to
resolving this problem.  We do not derive the details here, but the
Grassmann point of view on the Newton method starts with the second
covariant derivative of the block Rayleigh quotient (see \cite{us95}):
\begin{equation} {\Hess } (\Delta_1,\Delta_2) =
\tr\bigl(\Delta_1^TA\Delta_2 - (\Delta_1^T\Delta_2)Y\T
AY\bigr). \end{equation} From there to pick the Newton search
direction, we must solve for $\Delta$ such that $Y\T\Delta=0$ in the
Sylvester equation \begin{equation} \Pi(A\Delta -\Delta(Y\T AY))=-G,
\label{eq:sylvester} \end{equation} where $\Pi=(I-YY^T)$ denotes the
projection onto the tangent space of the Grassmann manifold, and
$G=\Pi AY$ is the gradient.  We then may use this search direction to
follow a geodesic.

Equation~(\ref{eq:sylvester}) is really a column by column Rayleigh
quotient iteration in disguise.  To see this, write $Y\T AY=:Q\Theta
Q^T$, $\Theta$ diagonal, and project Equation (\ref{eq:sylvester})
onto Image($\Pi$).  The projected equation is \begin{equation}
\bar{A}\bar{\Delta}-\bar{\Delta}{\Theta}=-\bar{G}, \label{eq:ritzrqi}
\end{equation} where the barred quantities are $\bar{A}=\Pi A\Pi$,
$\bar{\Delta}=\Pi \Delta Q$ is a column matrix of Ritz vectors
\cite{parlett80}, and $\bar{G}=GQ$.  If $\bar{\Delta}$ is any solution
to the above equation, then $\Delta=\bar{\Delta}Q^T$ is a solution to
Equation~(\ref{eq:sylvester}).  One way to interpret
Equation~(\ref{eq:ritzrqi}) is that we turned the orthogonality
constraint $Y^TY=I$ into an affine constraint $\Delta^TY=0$ by
differentiating.

Solving for $\Delta$ and exponentiating amounts to performing the
Rayleigh quotient iteration on each of the Ritz vectors associated
with the subspace $Y$.  Therefore, Newton's method applied to the
function $\tr Y\T AY$ on the Grassmann manifold converges cubically.
This is a generalization of the identification between RQI and
Newton's method applied to Rayleigh's quotient on the sphere
\cite{smith93a,smith94a}.  This method also has very much in common
with Chatelin's method, yet it is the natural algorithm from the
differential geometry point of view given the orthogonality
constraints $Y^TY=I$.



\bibliography{smith}
%/usr/local/lib/texmf/tex/latex209/distribs/latex/sty
\bibliographystyle{siamproc}

\end{document}
