How/why does the ICP algorithm using quaternions work? I've come across the Iterative Closest Point algorithm using quaternions (as described in "A Method for Registration of 3-D Shapes" by Besl and McKay) and I'm wondering, why it works. 
To me it seems like some magical algorithms, because I have no idea why the eigenvector corresponding to the maximum of the eigenvalue of the matrix Q turns out to be the optimal rotation. Also, I don't know why the matrix Q has to look like it does. 
I'd appreciate any explanation as well as any links to sources that derive or proof the algorithm.   
Thanks in advance! 
 A: That paper uses the Horn method for computing the best rotation (unit quaternion) aligning two sets or corresponding 3D vectors. 
Horn method is described in detail in the paper: "Closed-form solution of the absolute orientation using unit quaternions".
In his paper Horn is looking for the unit quaternion maximizing the following functional:
$arg \max_q E(q) = \sum_i E_i(q)$
S.t. $q q^* = 1$
Where $q$ is the unit quaternion to find and $E_i(q)$ is:
$E_i(q) = (q x_i q^*) \cdot y_i$
Given the sets of corresponding pure quaternion vectors $\{x_i\}$ and $\{y_i\}$.
The first thing that Horn does is to switch from quaternion algebra to linear algebra by expressing the quaternion products in terms of matrix products:
$arg \max_q E(q) = \sum_i E_i(q)$
$E_i(q) = q^T N_i q$
S.t $q^T q = 1$
This time $q$ is a  $4 \times 1$ column vector and $N_i$ is a $4 \times 4$ matrix formed using only $x_i$ and $y_i$.
Departing from quaternion expression:
$ (q x_i q^*) \cdot y_i = (q x_i) \cdot (y_i q)$
That can be transformed to linear algebra language as:
$(R_x q)^T (R_y q) =  q^T ( R_x^T R_y) q = q^T N_i q$
The maximization can be expressed as:
$arg \max_q E(q) = q^T (\sum_i N_i) q$
S.t $q^T q = 1$
Since $N = \sum_i N_i$ is symmetric (because is a sum of symmetric matrices), according to the spectral theorem, the quadratic form $q^T N q$ subject to constraint $q^T q = 1$ has a max value when $q$ is one of the eigenvectors of $N$, specifically the one with the largest eigenvalue.
Proof:
Since $N$ is real symmetric matrix the spectral theorem guarantee that $N$ can be decomposed as $N = V^T A V$. Where $V$ is orthonormal matrix which columns are eigenvectors of $N$ and $A$ is diagonal matrix with eigenvalues of $N$.
$q^T (V^T A V) q = p^T A p$
$p^T A p = \sum_i {a_i p_i^2}$
It turns out that:
$\sum_i {a_i p_i^2} < a_0 \sum_i { p_i^2}$ being $a_0$ the largest eigenvalue, since $\sum_i { p_i^2} = 1$. 
$\sum_i {a_i p_i^2} < a_0$ 
So we can pick $q$ to be the eigenvector (column of $V$)  corresponding to the largest eigenvalue $a_0$ to achieve the largest posible value of the quadratic form.
That concludes the proof.
You can find it more detailed at: Horn paper Appendix 3 or here: 
http://sfb649.wiwi.hu-berlin.de/fedc_homepage/xplore/tutorials/mvahtmlnode16.html
