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

\chapter{Eigenvalues and eigenvectors}

\section{Introduction}


A non-zero column-vector $\ve{v}$ is called an {\em
eigenvector}\index{eigenvector} of a matrix $A$ with an {\em
eigenvalue}\index{eigenvalue} $\lambda$, if
	\begin{equation}
A\ve{v}=\lambda \ve{v}\;.
	\end{equation}

If an $n\times n$ matrix $A$ is real and symmetric, $A^T=A$, then it has
$n$ real eigenvalues $\lambda_1,\dots,\lambda_n$, and its (orthogonalized)
eigenvectors $V=\{\ve{v}_{1},\dots,\ve{v}_{n}\}$ form a full
basis,
	\begin{equation}
VV^T=V^TV=\ve{1} \;,
	\end{equation}
in which the matrix is diagonal,
	\begin{equation}
V^{T}AV=
	\begin{bmatrix}
 \lambda_1 & 0         & \cdots & 0 \\
 0         & \lambda_2 &        & \vdots \\
 \vdots    &           & \ddots \\
 0         &  \cdots   &        & \lambda_n
	\end{bmatrix}
\;.
	\end{equation}
{\em Matrix diagonalization}\index{matrix diagonalization} means finding
all eigenvalues and (optionally) eigenvectors of a matrix.

Eigenvalues and eigenvectors enjoy a multitude of applications in
different branches of science and technology.

\section{Similarity transformations}
Orthogonal transformations\index{orthogonal transformation},
	\begin{equation}
A\rightarrow Q^{T}AQ\;,
	\end{equation}
where $Q^{T}Q=\ve{1}$, and, generally, similarity
transformations\index{similarity transformation},
	\begin{equation}
A\rightarrow S^{-1}AS\;,
	\end{equation}
preserve eigenvalues and eigenvectors. Therefore one of the strategies to
diagonalize a matrix is to apply a sequence of similarity transformations
(also called rotations) which (iteratively) turn the matrix into
diagonal form.

\subsection{Jacobi eigenvalue algorithm}
Jacobi eigenvalue algorithm is an iterative method to calculate the
eigenvalues and eigenvectors of a real symmetric matrix by a sequence
of Jacobi rotations\index{Jacobi rotation}.

Jacobi rotation is an orthogonal transformation which zeroes a pair of
the off-diagonal elements of a (real symmetric) matrix $A$,
	\begin{equation}
A \rightarrow A^\prime = J(p,q)^TAJ(p,q) \;:\;
A^\prime_{pq}=A^\prime_{qp}=0\;.
	\end{equation}
The orthogonal matrix $J(p,q)$ which eliminates the element $A_{pq}$
is called the Jacobi rotation matrix.  It is equal identity matrix except
for the four elements with indices $pp$, $pq$, $qp$, and $qq$,
	\begin{equation}
J(p,q)=
	\begin{bmatrix}
 1 &	    &	&   &	&   &  \\
   & \ddots &	 &   &	 & 0 &	\\ &	    & \cos\phi & \cdots &
   \sin\phi & & \\ &		   & \vdots & \ddots & \vdots &   &  \\
   &	      & -\sin\phi & \cdots & \cos\phi &  &   \\ &   0	 &   &
   &   & \ddots &   \\ &	&   &	&   &	& 1
	\end{bmatrix} \begin{array}{c} \\ \\ \leftarrow \mathrm{row}\;
	p \\ \\
\leftarrow \mathrm{row}\; q \\ \\ \\ \end{array} \;.
	\end{equation}
Or explicitly,
	\begin{eqnarray}
&&J(p,q)_{ij} = \delta_{ij}\;\forall\;ij \notin
 \{pq,qp,pp,qq\}\;;\nonumber\\
&&J(p,q)_{pp} = \cos\phi = J(p,q)_{qq} \;;\nonumber\\ &&J(p,q)_{pq} =
\sin\phi = -J(p,q)_{qp} \;.
	\end{eqnarray}
After a Jacobi rotation, $A\rightarrow A^\prime=J^TAJ$, the matrix
elements of $A'$ become
	\begin{eqnarray}
A^\prime_{ij}&=&A_{ij} \;\forall\; i\ne p,q\wedge j\ne p,q \nonumber\\
A^\prime_{pi}&=&A^\prime_{ip}=cA_{pi}-sA_{qi}\;\forall\;i\ne
p,q\;;\nonumber\\
A^\prime_{qi}&=&A^\prime_{iq}=sA_{pi}+cA_{qi}\;\forall\;i\ne
p,q\;;\nonumber\\
A^\prime_{pp}&=&c^2A_{pp}-2scA_{pq}+s^2A_{qq}\;;\nonumber\\
A^\prime_{qq}&=&s^2A_{pp}+2scA_{pq}+c^2A_{qq}\;;\nonumber\\
A^\prime_{pq}&=&A^\prime_{qp}=sc(A_{pp}-A_{qq})+(c^2-s^2)A_{pq}\;,
	\end{eqnarray}
where $c\equiv\cos\phi$, $s\equiv\sin\phi$.  The angle $\phi$ is chosen
such that after rotation the matrix element $A^\prime_{pq}$ is zeroed,
	\begin{equation}
\cot(2\phi)=\frac{A_{qq}-A_{pp}}{2A_{pq}}\;\Rightarrow\;
A^\prime_{pq}=0\;.
	\end{equation}

A side effect of zeroing a given off-diagonal element $A_{pq}$ by a
Jacobi rotation is that other off-diagonal elements are changed. Namely,
the elements of the rows and columns with indices equal to $p$ and $q$.
However, after the Jacobi rotation the sum of squares of all off-diagonal
elements is reduced. The algorithm repeatedly performs rotations until
the off-diagonal elements become sufficiently small.

The convergence of the Jacobi method can be proved for two strategies
for choosing the order in which the elements are zeroed:
\begin{enumerate}
\item {\em Classical method}: with each rotation the largest of
the remaining off-diagonal elements is zeroed.
\item {\em Cyclic method}: 
the off-diagonal elements are zeroed in strict order, e.g. row after row.
\end{enumerate}

Although the classical method allows the least number of rotations, it is
typically slower than the cyclic method since searching for the largest
element is an $O(n^2)$ operation.  The count can be reduced by keeping
an additional array with indexes of the largest elements in each row.
Updating this array after each rotation is only an $O(n)$ operation.

A {\em sweep} is a sequence of Jacobi rotations applied to all
non-diagonal elements.  Typically the method converges after a small
number of sweeps. The operation count is $O(n)$ for a Jacobi rotation
and $O(n^3)$ for a sweep.

The typical convergence criterion is that the sum of absolute values of
the off-diagonal elements is small, $\sum_{i<j}|A_{ij}|<\epsilon$,
where $\epsilon$ is the required accuracy. Other criteria
can also be used, like the largest off-diagonal element is small,
$\max{|A_{i<j}|}<\epsilon$, or the diagonal elements have not changed
after a sweep.

The eigenvectors can be calculated as $V=\ve{1}J_1J_2...$,
where $J_i$ are the successive Jacobi matrices. At each stage the
transformation is
	\begin{eqnarray}
V_{ij}&\rightarrow& V_{ij}\;,\; j\ne p,q\nonumber\\
V_{ip}&\rightarrow& cV_{ip}-sV_{iq}\\
V_{iq}&\rightarrow& sV_{ip}+cV_{iq}\nonumber
	\end{eqnarray}

Alternatively, if only one (or few) eigenvector $\ve{v}_{k}$
is needed, one can instead solve the (singular) system $(A-\lambda
_{k})\ve{v}=0$.

\begin{table}[h]\caption{Jacobi diagonalization in C. The header {\tt
matrix.h} is supposed to define the self-explanatory functions
{\tt vector\_get},
{\tt vector\_set},
{\tt matrix\_get},
{\tt matrix\_set}, and
{\tt matrix\_set\_identity}.}
\lstset{frame=single}
\lstinputlisting
%{/usr/users/fedorov/public_html/numeric/c++/jacobi.eigenvalue/jacobi.cc}
{./c/jacobi.eigenvalue/jacobi.c}
	\end{table}

\subsubsection{QR/QL algorithm}

An orthogonal transformation of a real symmetric matrix,
$\ma{A}\rightarrow \ma{Q}^{T}\ma{A}\ma{Q}=\ma{R}\ma{Q}$, where
$\ma{Q}$ is from the QR-decomposition of $\ma{A}$, partly turns
the matrix $\ma{A}$ into diagonal form. Successive iterations
eventually make it diagonal. If there are degenerate eigenvalues
there will be a corresponding block-diagonal sub-matrix. For
convergence properties it is of advantage to use {\em shifts}:
instead of QR[$\ma{A}$] we do QR[$\ma{A}-s\ve{1}$] and then
$\ma{A}\rightarrow \ma{R}\ma{Q}+s\ve{1}$. The shift $s$ can be
chosen as $\ma{A}_{nn}$. As soon as an eigenvalue is found the matrix
is deflated, that is the corresponding row and column are crossed out.

Accumulating the successive transformation matrices $\ma{Q}_{i}$
into the total matrix $\ma{Q}=\ma{Q}_{1}\dots \ma{Q}_{N}$, such that
$Q^{T}\ma{A}\ma{Q}=\Lambda$, gives the eigenvectors as columns of the
$\ma{Q}$ matrix.

If only one (or few) eigenvector $\ve{v}_{k}$ is needed one can instead
solve the (singular) system $(\ma{A}-\lambda _{k})\ve{v}=0$.

\subsubsection{Tridiagonalization.}
Each iteration of the QR/QL algorithm is an $O(n^{3})$ operation. On a
tridiagonal matrix it is only $O(n)$. Therefore the effective strategy is
first to make the matrix tridiagonal and then apply the QR/QL algorithm.
Tridiagonalization of a matrix is a non-iterative operation with a fixed
number of steps.


\section{Eigenvalues of updated matrix}

In practice it happens quite often that the (symmetric) matrix $A$ to
be diagonalized is given in the form of a diagonal matrix, $\ma{D}$,
plus an update matrix, $\ma{W}$,
	\begin{equation}
\ma{A}=\ma{D}+\ma{W} \:,
	\end{equation}
where the update $\ma{W}$ is a simpler --- in a certain sense --- matrix
which allows a more efficient calculation of the updated eigenvalues,
as compared to general diagonalization algorithms.

The most common updates are
	\begin{itemize}

	\item symmetric rank-1 update,
	\begin{equation}
\ma{W}=\ve{u}\ve{u}^T \:,
	\end{equation}
where $\ve{u}$ is a columnt-vector;

	\item symmetric rank-2 update,
	\begin{equation}
\ma{W}=\ve{u}\ve{v}^T + \ve{v}\ve{u}^T \:;
	\end{equation}

	\item symmetric row/column update -- a special case of rank-2 update,
	\begin{equation}
\ma{W}=
	\begin{bmatrix}
0 & \dots & u_1 & \dots & 0 \\
\vdots & \ddots & \vdots & \ddots & \vdots \\
u_1 & \dots & u_p & \dots & u_n \\
\vdots & \ddots & \vdots & \ddots & \vdots \\
0 & \dots & u_n & \dots & 0 \\
	\end{bmatrix}
\equiv \ve{e}(p)\ve{u}^T+\ve{u}\ve{e}(p)^T
\:,
	\end{equation}
where $\ve{e}(p)$ is the unit vector in the $p$-direction.

	\end{itemize}

\subsection{Rank-1 update}
We assume that the size-$n$ real symmetric matrix $A$ to be diagonalized
is given in the form of a diagonal matrix plus a rank-1 update,
	\begin{equation}
\ma{A}=\ma{D}+\sigma\ve{u}\ve{u}^T \:,
	\end{equation}
where $\ma{D}$ is a diagonal matrix with diagonal elements
$\{d_1,\dots,d_n\}$ and $\ve{u}$ is a given vector.  The diagonalization
of such matrix can be done in $O(m^2)$ operations, where $m\le n$ is the
number of non-zero elements in the update vector $\ve{u}$, as compared
to $O(n^3)$ operations for a general diagonalization~\cite{golub}.

The eigenvalue equation for the updated matrix reads
	\begin{equation}
\left(\ma{D}+\sigma\ve{u}\ve{u}^T\right)\ve{q}=\lambda \ve{q} \:,
	\end{equation}
where $\lambda$ is an eigenvalue and $\ve{q}$ is the corresponding
eigenvector.  The equation can be rewritten as
	\begin{equation}
\left(\ma{D}-\lambda\ma{1}\right)\ve{q} + \sigma\ve{u}\ve{u}^T\ve{q}=0 \:.
	\end{equation}
Multiplying from the left with
$\ve{u}^T\left(\ma{D}-\lambda\ma{1}\right)^{-1}$ gives
	\begin{equation}
\ve{u}^T\ve{q} +
\ve{u}^T\left(\ma{D}-\lambda\ma{1}\right)^{-1}\sigma\ve{u}\ve{u}^T\ve{q}=0 \:.
	\end{equation}
Finally, dividing by $\ve{u}^T\ve{q}$ leads to the (scalar) {\em secular
equation}\index{secular equation} (or {\em characteristic
equation}\index{characteristic equation}) in $\lambda$,
	\begin{equation}
1+\sum_{i=1}^{m}\frac{\sigma u_i^2}{d_i-\lambda} = 0 \:,
	\end{equation}
where the summation index counts the $m$ non-zero components of the
update vector $\ve{u}$.  The $m$~roots of this equation determine the
(updated) eigenvalues\footnote{Multiplying this equation
by $\prod_{i=1}^{m}(d_i-\lambda)$ leads to an equivalent polynomial
equation of the order $m$, which has exactly $m$ roots.}.

Finding a root of a rational function requires an iterative technique,
such as the Newton-Raphson method. Therefore diagonalization of an updated
matrix is still an iterative procedure.  However, each root can be found in
$O(1)$ iterations, each iteration requiring $O(m)$ operations. Therefore the
iterative part of this algorithm needs $O(m^2)$ operations.

Finding roots of this particular secular equation can be
simplified by utilizing the fact that its roots are bounded by
the eigenvalues $d_i$ of the matrix $\ma{D}$.  Indeed if we denote
the roots as $\lambda_1,\lambda_2,\dots,\lambda_n$ and assume that
$\lambda_i\le\lambda_{i+1}$ and $d_i\le d_{i+1}$, it can be shown that
\begin{enumerate} \item if $\sigma \ge 0$,
	\begin{eqnarray}
&d_i \le \lambda_i \le d_{i+1} \:,&\; i=1,\dots,n-1 \:, \\ &d_n \le
\lambda_n \le d_n + \sigma \ve{u}^T\ve{u} \:;&
	\end{eqnarray}
\item if $\sigma \le 0$,
	\begin{eqnarray}
&d_{i-1} \le \lambda_i \le d_i \:,&\; i=2,\dots,n \:, \\ &d_1 + \sigma
\ve{u}^T\ve{u} \le \lambda_1 \le d_1 \:.&
	\end{eqnarray}
\end{enumerate}

\subsection{Symmetric row/column update}

The matrix $\ma{A}$ to be diagonalized is given in the form
	\begin{equation}
\ma{A}=
\ma{D} + \ve{e}(p)\ve{u}^T+\ve{u}\ve{e}(p)^T =
	\begin{bmatrix}
d_1 & \dots & u_1 & \dots & 0 \\
\vdots & \ddots & \vdots & \ddots & \vdots \\
u_1 & \dots & d_p & \dots & u_n \\
\vdots & \ddots & \vdots & \ddots & \vdots \\
0 & \dots & u_n & \dots & d_n \\
	\end{bmatrix}
\:,
	\end{equation}
where $\ma{D}$ is a diagonal matrix with diagonal elements
$\{d_i|i=1,\dots,n\}$, $\ve{e}(p)$ is the unit vector in the $p$-direction,
and $\ve{u}$ is a given update vector where the $p$-th element
can be assumed to equal zero, $u_p=0$,
without loss of generality. Indeed, if the
element is not zero, one can simply redefine $d_p\rightarrow d_p+2u_p$,
$u_p\rightarrow 0$.

The eigenvalue equation for matrix $\ma{A}$ is given as
	\begin{equation}
(D-\lambda)\ve{x} + \ve{e}(p)\ve{u}^T\ve{x}+\ve{u}\ve{e}(p)^T\ve{x} = 0 \:,
	\end{equation}
where $\ve{x}$ is an eigenvector and $\lambda$ is the corresponding
eigenvalue.
The component number $p$ of this vector-equation reads
	\begin{equation}\label{eig:eq:rc1}
(d_p-\lambda)x_p + \ve{u}^T\ve{x} = 0 \:,
	\end{equation}
while the component number $k\ne p$ reads
	\begin{equation}
(d_k-\lambda)x_k + u_k x_p = 0 \:,
	\end{equation}
Dividing the last equation by $(d_k-\lambda)$, multiplying from the left
with $\sum_{k=1}^{n}u_k$, substituting $\ve{u}^T\ve{x}$ using
equation~(\ref{eig:eq:rc1}) and dividing by $x_p$ gives the secular
equation,
	\begin{equation}
-(d_p-\lambda)+\sum_{k\ne p}^{n}\frac{u_k^2}{d_k-\lambda} = 0 \:,
	\end{equation}
which determines the updated eigenvalues.

	\begin{equation}
	\end{equation}


\subsection{Symmetric rank-2 update}
A symmetric rank-2 update can be represented as two consecutive rank-1
updates,
	\begin{equation}
\ve{u}\ve{v}^T + \ve{v}\ve{u}^T = \ve{a}\ve{a}^T - \ve{b}\ve{b}^T \:,
	\end{equation}
where
	\begin{equation}
\ve{a}=\frac{1}{\sqrt{2}}(\ve{u}+\ve{v}) \:, \;
\ve{b}=\frac{1}{\sqrt{2}}(\ve{u}-\ve{v}) \:.
	\end{equation}
The eigenvalues can then be found by applying the rank-1 update method
twice.
