\documentclass[11pt]{article}
\usepackage{latexsym}
%\input psfig.sty
%\input spec_frac
 \def\emptyline{\vspace{12pt}}
\usepackage{epsfig}
\usepackage{graphicx}
\newcounter{algo}
\renewcommand{\baselinestretch}{1.5}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% USER DEFINED COMMANDS *
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\bdm}{\begin{displaymath}}
\newcommand{\edm}{\end{displaymath}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINITIONS *
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\bbox#1{\mbox{\boldmath $#1$}}
\def\eq#1{Eq.~(\ref{eq#1})}  %label equations as \label{eq12} etc...
\def\Eq#1{Equation~(\ref{eq#1})} %label equations as \label{eq12} etc...
\def\fig#1{Fig.~\ref{fig#1}}  %label figures as \label{fig5} etc...
\def\Fig#1{Figure~\ref{fig#1}}
\def\dt{$\Delta t$}
\def\order#1{${\cal O}(\Delta t^{#1})$}

\begin{document}
\begin{flushleft}
        {\normalsize \it \bf The University of New Mexico} \\
        {\normalsize \it E.A. Coutsias} \\
        {\Large \bf 504 '07-Project}
        \end{flushleft}
{\bf Part A: "To build a better symmetric eigensolver!"}\\
We want to write a fast eigensolver for moderate-sized
symmetric matrices. It should be able to give efficiently either
a selected few or all eigenvalues and, optionally, eigenvectors.
All tricks you know must be incorporated, with key objective to
perform the calculation in the most efficient way possible:
\begin{itemize}
\item Conversion to (symmetric = tridiagonal) Hessenberg form.
\item Efficient/stable solution of triangular system (Rayleigh 
quotient iteration step) using either sparse LU with pivoting
or Givens QR.
\item for the case when all eigenvalues are required, follow the 
steps given in problem 29.1 of Trefethen and Bau. For the case
where only one eigenpair is sought, use Rayleigh quotient iteration.
\item if the largest eigenvalue is sought, use Gershgorin estimates
to get a good upper bound for its value before starting the iteration.
You must ensure that you did not converge to a lesser eigenvalue.
\end{itemize} 

\newpage

{\bf Part B: Application to the optimal superposition problem}\\
Demonstrate the power of your fast Rayleigh quotient
by solving the optimal superposition problem \cite{CSD}. 
The mathematical statement of the problem is:\\
\noindent {\em Given an ordered set of vectors ${\bf y}_k$ (target)
and a second set ${\bf x}_k$ (model), $1\le k\le N$, find an orthogonal
transformation $ {\cal U}$ and a translation {\bf r} such that the residual
$E$ 
\begin{equation}
E := \frac{1}{N}\sum_{k=1}^N |{\cal U} {\bf x}_k + {\bf r} - {\bf y}_k|^2
\label{res}
\label{eq1}
\end{equation}
is minimized.}

Shifting the two sets, $\{{\bf x}_k\}$ and $\{{\bf y}_k\}$, to their respective barycenters
$\bar{{\bf x}}$ and $\bar{{\bf y}}$, we introduce:
\[ {\tilde {\bf x}_k} := {\bf x}_k - \bar{{\bf x}}, \
{\tilde {\bf y}_k} := {\bf y}_k - \bar{{\bf y}}. \]
Below we drop the tildes (i.e., we will assume that both sets have
been shifted to bring their respective barycenters to the origin),
and then the residual becomes:
\be
E  = \frac{1}{N}\sum_{k=1}^N | {\bf x}'_k -
 {\bf y}_k |^2
\ee
where we set $ {\bf x}' := {\cal U} {\bf x} $.
                                                                                
In the sequel we will make use of the $3\times N$ matrices ${\cal X}$
and ${\cal Y}$ where the $k$-th column of ${\cal X}$ is the vector
${\bf x}_k$ and similarly for ${\cal Y}$.

We now promote ${\bf x}_k$ and ${\bf y}_k$ to {\em pure quaternions}
(see \ref{sec:quaternions}),
$x_k := (0, {\bf x}_k)$ with $x_k^c = - x_k$ and similarly for $y_k$.
The rotation ${\cal U}(q)$ on ${\bf x}_k$ is then written as
\[ \left( 0, {\cal U}(q) {\bf x}_k \right) = q x_k q^c. \]
The residual is written, in terms of quaternions, as
\begin{equation}
E_q = \frac{1}{N}\sum_{k=1}^N \left( q x_k q^c - y_k\right)
\left(q x_k q^c - y_k\right)^c.
\label{resq}
\end{equation}

 By employing the rules of quaternion algebra, as discussed in
\cite{CSD}, we find that the residual can be written as
\be
NE_q = \sum_{k=1}^N\left( |{\bf x}_k|^2 + |{\bf y}_k|^2 \right)
- 2{\cal Q}^T {\cal F} {\cal Q},
\ee
where
${\cal Q} = (q_0, q_1, q_2, q_3)^T$ is a 4-vector form of the
quaternion $q$ and 
\[ {\cal F} := -\sum_{k=1}^N {\cal A}_L(y_k) {\cal A}_R(x_k), \]
(see \ref{eq38}).
In terms of the matrix elements of the correlation matrix
${\cal R}$ (\ref{eq:Kabsch}),
\begin{equation}
\label{eq:Kabsch}
{\cal R} := {\cal X}{\cal Y}^T = \sum_{k=1}^N {\bf x}_k {\bf y}_k^T \rightarrow R_{ij}
= \sum_{k=1}^N x_{ik} y_{jk} \ ,\ i,j=1,2,3
%\left(
%\begin{array}{cccc}
%R_{11} & R_{12} & R_{13} \\
%R_{21} & R_{22} & R_{23} \\
%R_{31} & R_{32} & R_{33}
%\end{array}
%\right)
\ ,
\end{equation}
the explicit form of the matrix ${\cal F}$ is\\
${\cal F}=$
\begin{small}
\be
%{\cal F} =
\left(
\begin{array}{cccc}
R_{11} + R_{22} + R_{33} & R_{23} - R_{32}
& R_{31} - R_{13} & R_{12} - R_{21} \\
R_{23} - R_{32} & R_{11} - R_{22} - R_{33}
& R_{12} + R_{21} & R_{13} + R_{31} \\
R_{31} - R_{13} & R_{12} + R_{21}
& -R_{11} + R_{22} - R_{33} & R_{23} + R_{32} \\
R_{12} - R_{21} & R_{13} + R_{31}
& R_{23} + R_{32} &-R_{11} - R_{22} + R_{33}
\end{array}
\right) \ .
\label{eq:F_explicit}
\ee
\end{small}
The problem has in this way been reduced to that of finding the extrema
of a quadratic form ${\cal Q}^T {\cal F} {\cal Q}$ in the four variables
$q_i$, $i = 0,1,2,3$, subject to the constraint ${\cal Q}^T {\cal Q} =1$.
Note that here we are using the vector ${\cal Q}$, so that the squared norm
$q^cq$ is written equivalently as ${\cal Q}^T{\cal Q}$.
${\cal Q}^T {\cal F} {\cal Q}$ is the standard Rayleigh quotient
for a symmetric matrix ${\cal F}$,
and the maximum value achieved by ${\cal Q}^T {\cal F} {\cal Q}$ is
equal to its largest eigenvalue.
Thus, the desired minimization leads to the eigenproblem
\be
{\cal F} {\cal Q} = \lambda {\cal Q}.
\label{eig}
\ee
We see that the extremum $\lambda$ is equal to one of the
eigenvalues of a
$4\times 4$ symmetric, traceless matrix, and the corresponding eigenvector
gives one of the candidate rotations that extremize the residual, as
sought.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Quaternion notation and properties}
\label{sec:quaternions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                                                                
We review basic properties of quaternions below. We follow the notes
of Coutsias and Romero \cite{notes}, and refer the reader to these notes
for more
details. The book by Rappaport \cite{rap}, which discusses quaternions in the
context of molecular dynamics, is also useful as an introduction.
We define a quaternion $q$ as a 4-vector $q=(q_0,q_1,q_2,q_3)
\equiv(q_0,{\bf q})$, where ${\bf q} = (q_1,q_2,q_3)$,
with the obvious addition law and the fundamental multiplication law given by
\bea
a = (a_0, {\bf a}), &\quad& b = (b_0, {\bf b}), \nonumber \\
a+b &=& (a_0+b_0, {\bf a}+{\bf b}), \nonumber \\
\label{eq26}
a\,b &=& (a_0 b_0 - {\bf a}\cdot{\bf b},
a_0 {\bf b}+b_0 {\bf a}+ {\bf a}\times{\bf b}),
\eea
where the multiplication can be seen to be associative, $a(bc) = (ab)c$,
but not commutative, $ab \ne ba$, in general.
We also define the \textit{conjugate} quaternion $q^c$, the squared norm of a
quaternion $N(q)$, and the inverse quaternion $q^{-1}$ by
\bea
\label{eq27}
q^c &:=& (q_0,-{\bf q}), \qquad (a b)^c = b^c a^c, \\
\label{eq28}
N(q) &:=& q^cq = qq^c = q_0^2+{\bf q}\cdot{\bf q}, \\
\label{eq29}
q^{-1} &:=& \frac{q^c}{N(q)} = \frac{(q_0,-{\bf q})}{N(q)}.
\eea
                                                                                
An important class of quaternions is the \textit{pure} quaternions
which have the form $q = (0,{\bf q})$.
Pure quaternions can be considered as ordinary 3-vectors
${\bf q}=(q_1,q_2,q_3)$ which have been mapped to 4-vectors,
$(0,{\bf q}) = (0,q_1,q_2,q_3)$.
Note for pure quaternions $a,b$, \eq{26} and \eq{27} yield
\bea
\label{eq31}
a = (0,\bbox{a}), &\quad& b = (0,\bbox{b}), \nonumber \\
a\,b &=& (-\bbox{a}\cdot\bbox{b}, \, \bbox{a}\times\bbox{b}), \\
\label{eq32}
a^c &=& (0,-\bbox{a}) = -a.
\eea
                                                                                
The quaternion product can be formulated
as a matrix multiplication as follows:
\bea
\label{eq36}
pq & =: & {\cal A}_L(p) {\cal Q}, \qquad \textrm{ $p$ operates
on $q$ from the left,}\\
\label{eq37}
qp & =: & {\cal A}_R(p) {\cal Q}, \qquad \textrm{ $p$ operates on
$q$ from the right,}
\eea
where ${\cal A}_L(p)$ and ${\cal A}_R(p)$ are $4\times 4$ matrices,
and ${\cal Q}=(q_0,q_1,q_2,q_3)^T$ is the column 4-vector representation
of quaternion $q$.
By convention, the matrix multiplication form of the quaternion product
(the right-hand side of the above)
is defined to act from the left on a column 4-vector.
For $p=(p_0,p_1,p_2,p_3)$ the matrices ${\cal A}_L(p)$ and
${\cal A}_R(p)$ are given by
\be
{\cal A}_R(p) = \left(
\begin{array}{cccc}
p_0 & -p_1 & -p_2 & -p_3 \\
p_1 & p_0 & p_3 & -p_2 \\
p_2 & -p_3 & p_0 & p_1 \\
p_3 & p_2 & -p_1 & p_0
\end{array}
\right),
\quad
{\cal A}_L(p) = \left(
\begin{array}{cccc}
p_0 & -p_1 & -p_2 & -p_3 \\
p_1 & p_0 & -p_3 & p_2 \\
p_2 & p_3 & p_0 & -p_1 \\
p_3 & -p_2 & p_1 & p_0
\end{array}
\right).
\label{eq38}
\ee
                                                                                
Quaternions provide a natural coordinate system for ${\cal SU}(3)$, the
group of proper rotations of 3-space. Thus they can describe rotations without
the singularities of, say, the Euler angles rotation matrices. For the
3-vector $\bbox{r}'$ obtained by a rotation of the vector $\bbox{r}$
via the orthogonal rotation matrix ${\cal{U}}$ we have
\bea
\label{eq33}
r' = (0,{\bf r}'), &\quad& r = (0,{\bf r}), \nonumber \\
r' &=& q \, r \, q^c, \\
\label{eq34}
(0,\bbox{r'}) &=& (q_0,{\bf q}) \; (0,\bbox{r}) \; (q_0, -{\bf q})
\nonumber \\ &=& (0,\, {{\cal{U}}(q)}\bbox{r}).
\eea
Using the notation of Eq.\ (\ref{eq38}), we can express the rotation matrix
as a product of two matrices, each depending linearly on a unit quaternion $q$
($qq^c = 1$):
\[ \left(
\begin{array}{cc}
1 & {\bf 0}^T \\
{\bf 0} & {\cal U}(q)
\end{array}
\right)
= {\cal A}_L(q) {\cal A}_R(q^c). \]
%(One must keep in mind when mixing matrices and quaternions that
%\textit{quaternion operations are not associative with
%ordinary $4\times 4$ matrix multiplications}).
The $3 \times 3$ rotation matrix ${\cal{U}}$ is given by
\be
{\cal{U}}(q) =
\left(
\begin{array}{ccc}
q_0^2 + q_1^2-q_2^2-q_3^2 & 2(q_1 q_2-q_0 q_3)      & 2(q_1 q_3+q_0 q_2) \\
2(q_1 q_2+q_0 q_3)        & q_0^2-q_1^2+q_2^2-q_3^2 & 2(q_2 q_3-q_0 q_1) \\
2(q_1 q_3-q_0 q_2)        & 2(q_2 q_3+q_0 q_1)      & q_0^2-q_1^2-q_2^2+q_3^2
\end{array}
\right)\ .
\label{eq35}
\ee
Since ${\cal{U}}$ is orthogonal we have ${\cal{U}}^{-1} = {\cal{U}}^T$.
Note that by its construction, ${\cal{U}}(q)$ is a proper rotation whose
angle $\theta$ and axis ${\bf c}$ are seen in the form of $q$
\cite{notes, rap}:
\[ q = (q_0, {\bf q}) = (\cos{(\theta/2)}, \sin{(\theta/2)} {\bf c}). \]
It is easy to see that $\left\{ \theta, {\bf c} \right\}$ and
$\left\{ -\theta, -{\bf c} \right\}$
produce the same quaternion $q$ and hence the same rotation.
Also, $q$ and $-q$ produce the same rotation, a fact easily deduced from
the form of the rotation matrix ${\cal{U}}(q)$.

\begin{thebibliography}{99}
\bibitem{CSD}  E.A. Coutsias, C. Seok, K.A. Dill 
Using Quaternions to calculate RMSD  
(Journal of Computational Chemistry, 2004, vol 25, pp 1849-1857)
\bibitem{notes} Coutsias, E.\ A.; Romero, L.
The Quaternions with applications to Rigid Body Dynamics,
Sandia National Laboratories, 2004,
Technical Report SAND2004-0153.
\bibitem{rap} Rappaport, D.\ C. The Art of Molecular Dynamics
Simulation; Cambridge Univ.\ Press: New York, 1995, Chapter 8, pp. 191-221.
\end{thebibliography}

\end{document}
