% Except otherwise noted,
% this work is © 2009 Dmitri Fedorov, under CC-BY-SA license, version 3+.

\chapter{Power iteration methods and Krylov subspaces}

\section{Introduction}
Power method is an iterative method to calculate an eigenvalue and the
corresponding eigenvector of a matrix $A$ using the {\em power iteration}
	\begin{equation}
\ve{x}_{i+1}=\ma{A}\ve{x}_{i}\;.
	\end{equation}
The iteration converges to the eigenvector with the largest eigenvalue. The
eigenvalue can be estimated using the {\em Rayleigh quotient}
	\begin{equation}
\lambda[\ve{x}_i]=
\frac{\ve{x}_{i}^{T}\ma{A}\ve{x}_{i}}{\ve{x}_i^T\ve{x}_i}=
\frac{\ve{x}_{i+1}^{T}\ve{x}_{i}}{\ve{x}_i^T\ve{x}_i}.
	\end{equation}

The {\em inverse power iteration}\index{inverse power iteration}  with
the inverse matrix,
	\begin{equation}
\ve{x}_{i+1}=\ma{A}^{-1}\ve{x}_{i} \:,
	\end{equation}
converges to the smallest eigenvalue of matrix $A$.

Alternatively, the {\em shifted inverse iteration},
	\begin{equation}\label{eq:ipms}
\ve{x}_{i+1}=(\ma{A}-s\ma{1})^{-1}\ve{x}_{i} \:,
	\end{equation}
where $\ma{1}$ signifies a unity matrix of the same size as $\ma{A}$,
converges to the eigenvalue closest to the given number $s$.

The {\em inverse iteration method} is a refinement of the
inverse power method where the trick is not to invert the matrix
in~(\ref{eq:ipms}) but rather solve the linear system
	\begin{equation}
(\ma{A}-s\ma{1})\ve{x}_{i+1}=\ve{x}_{i}
	\end{equation}
using e.g. QR-decomposition.

The better approximation $s$ to the sought eigenvalue is chosen, the
faster convergence one gets. However, incorrect choice of $s$ can lead to
slow convergence or to the convergence of a different eigenvector. Usually
in practice the method is used when good approximation for the eigenvalue
is known, and hence one needs only few (quite often just one) iteration.

One can update the estimate for the eigenvalue using the Rayleigh
quotient\index{Rayleigh quotient} $\lambda[\ve{x}_i]$ after
each iteration and get faster convergence for the price of $O(n^3)$
operations per QR-decomposition; or one can instead make more
iterations (with $O(n^2)$ operations per iteration) using the same
matrix $(\ma{A}-s\ma{1})$. The optimal strategy is probably an
update after several iterations.

\section{Krylov subspaces}
When calculating an eigenvalue of a matrix $\ma{A}$ using the power
method, one starts with an initial random vector $\ve{b}$
and then computes iteratively the sequence $A\ve{b},
\ma{A}^2\ve{b},\dots,\ma{A}^{n-1}\ve{b}$ normalising and storing the
result in $\ve{b}$ on each iteration.  The sequence converges to
the eigenvector of the largest eigenvalue of $\ma{A}$.

The set of vectors
	\begin{equation}
{\cal K}_n=
\left\{\ve{b},\ma{A}\ve{b},\ma{A}^2\ve{b},\dots,\ma{A}^{n-1}\ve{b}\right\}\:,
	\end{equation}
where $n<\mathrm{rank}(\ma{A})$, is called the order-$n$ {\em Krylov
matrix}\index{Krylov matrix}, and the subspace spanned by these vectors
is called the order-$n$ {\em Krylov subspace}\index{Krylov subspace}.
The vectors are not orthogonal but can be made so e.g. by Gram-Schmidt
orthogonalisation.

For the same reason that $\ma{A}^{n-1}\ve{b}$ approximates the dominant
eigenvector one can expect that the other orthogonalised vectors
approximate the eigenvectors of the $n$ largest eigenvalues.

Krylov subspaces are the basis of several successful iterative methods in
numerical linear algebra, in particular: Arnoldi and Lanczos methods for
finding one (or a few) eigenvalues of a matrix; and  GMRES (Generalised
Minimum RESidual) method for solving systems of linear equations.

These methods are particularly suitable for large sparse matrices as they
avoid matrix-matrix operations but rather multiply vectors by matrices
and work with the resulting vectors and matrices in Krylov subspaces of
modest sizes.

\section{Arnoldi iteration}
Arnoldi iteration\index{Arnoldi iteration} is an algorithm where the
order-$n$ orthogonalised Krylov matrix $\ma{Q}_n$ for a given matrix
$\ma{A}$ is built using stabilised Gram-Schmidt process:
	\begin{code}
start with a set $\ma{Q}=\{\ve{q}^1\}$ consisting of one random
normalised vector $\ve{q}^1$; \\
repeat for $k=2$ to $n$~: \\
\hspace*{1em} make a new vector $\ve{q}^k=\ma{A}\ve{q}^{k-1}$ \\
\hspace*{1em} orthogonalise $\ve{q}^k$ to all vectors $\ve{q}^i\in \ma{Q}$
storing ${\ve{q}^i}^\dagger \ve{q}^k\rightarrow h_{i,k-1}$ \\
\hspace*{1em} normalise $\ve{q}^k$  storing $\|\ve{q}^k\|\rightarrow
h_{k,k-1}$ \\
\hspace*{1em} add $\ve{q}^k$ to the set $\ma{Q}$ \\
	\end{code}
By construction the matrix $\ma{H}_n$ made of the elements $h_{jk}$ is an
upper Hessenberg matrix,
	\begin{equation}
\ma{H}_n = \begin{bmatrix}
   h_{1,1} & h_{1,2} & h_{1,3} & \cdots  & h_{1,n} \\
   h_{2,1} & h_{2,2} & h_{2,3} & \cdots  & h_{2,n} \\
   0       & h_{3,2} & h_{3,3} & \cdots  & h_{3,n} \\
   \vdots  & \ddots  & \ddots  & \ddots  & \vdots  \\
   0       & \cdots  & 0     & h_{n,n-1} & h_{n,n} 
\end{bmatrix}\;,
	\end{equation}
which is a partial orthogonal reduction of $A$ into Hessenberg form,
	\begin{equation}
\ma{H}_n=\ma{Q}_n^\dagger \ma{A} \ma{Q}_n\:.
	\end{equation}

The matrix $\ma{H}_n$ can be viewed as a representation of $\ma{A}$ in
the Krylov subspace ${\cal K}_n$. The eigenvalues and eigenvectors of the
matrix $\ma{H}_n$ approximate the largest eigenvalues of matrix $\ma{A}$.

Since $\ma{H}_n$ is a Hessenberg matrix of modest size its eigenvalues can
be relatively easily computed with standard algorithms.

In practice if the size $n$ of the Krylov subspace becomes too large the
method is restarted.

\section{Lanczos iteration}
Lanczos iteration\index{Lanczos iteration} is Arnoldi iteration for
Hermitian matrices, in which case the Hessenberg matrix $\ma{H}_n$
of Arnoldi method becomes a tridiagonal matrix $\ma{T}_n$.

The Lanczos algorithm thus reduces the original hermitian $N\times N$
matrix $\ma{A}$ into a smaller $n\times n$ tridiagonal matrix $\ma{T}_n$
by an orthogonal projection onto the order-$n$ Krylov subspace. The
eigenvalues and eigenvectors of a tridiagonal matrix of a modest size
can be easily found by e.g. the QR-diagonalisation method.

In practice the Lanczos method is not very stable due to round-off errors
leading to quick loss of orthogonality. The eigenvalues of the resulting
tridiagonal matrix may then not be a good approximation to the original
matrix. Library implementations fight the stability issues by trying to
prevent the loss of orthogonality and/or to recover the orthogonality
after the basis is generated.

\section{Generalised minimum residual (GMRES)}

GMRES\index{GMRES} is an iterative method for the numerical solution of
a system of linear equations,
	\begin{equation}\label{eq-gmres-axb}
\ma{A}\ve{x}=\ve{b}\:,
	\end{equation}
where the exact solution $\ve{x}$ is approximated by the
vector $\ve{x}_n\in{\cal K}_n$ that minimises the residual
$\ma{A}\ve{x}_n-\ve{b}$
in the Krylov subspace ${\cal K}_n$ of matrix $\ma{A}$,
	\begin{equation}
\ve{x} \approx \ve{x}_n \leftarrow
\min_{\ve{x}\in{\cal K}_n}\|\ma{A}\ve{x}-\ve{b}\|\:.
	\end{equation}

The vector $\ve{x}_n\in{\cal K}_n$ can be represented as
$\ve{x}_n=\ma{Q}_n\ve{y_n}$ where $\ma{Q}_n$ is the projector on
the space ${\cal K}_n$ and $\ve{y_n}$ is an $n$-dimensional vector.
Substituting this into~(\ref{eq-gmres-axb}) gives
	\begin{equation}
\ma{A}\ma{Q}_n\ve{y_n}=\ve{b} \;,
	\end{equation}
which can be solved by least-squares method.

