Superposition

3-D Rotation Problem: Given S = s0, s1, ..., sn-1, and T = t0, t1, ... tn-1, where each si and each ti is a 3-D vector, ⟨?x, ?y, ?z⟩, find a rotation, R, by some angle θ, about some axis, u, through the origin, ⟨0, 0, 0⟩, so as to minimise the sum of squared errors between S' = map R S and T, that is  ∑i (norm(s'i - ti))2.
 
Note, it is given that si ~ ti,  0 ≤ i < n.
To find the general rigid-body transformation of S that gives the least-squares superposition of S on T (a version of the Total Least-Squares (TLS) problem), first make the centre of gravity (CoG) of T the origin, then translate S's CoG (carrying S with it) to that origin, and finally solve the rotation problem above.
 
This problem appears in the superposition of two protein structures in Bioinformatics.
 
Kearsley's Method: This uses the quaternion representation of 3-D rotations, and an Eigen-value finding algorithm, such as the Jacobi algorithm, say. A pure rotation is represented by a constrained quaternion -- a general quaternion having four degrees of freedom and a 3-D rotation only three -- so a Lagrange multiplier also gets a look in. All this leads to an Eigen-value problem on a 4×4 matrix. The smallest Eigen-value, λ, is the least-square error, and the corresponding Eigen-vector yields the required rotation (quaternion).
S:  
T:
output
— L.A. © 

xm2 = ∑i xm*xm; xp2 = ∑i xp*xp;
  ym2 = ∑i ym*ym; yp2 = ∑i yp*yp;
    zm2 = ∑i zm*zm; zp2 = ∑i zp*zp;
xmym = ∑i xm*ym;
  xmyp = ∑i xm*yp;
    xmzm = ∑i xm*zm;
      xmzp = ∑i xm*zp;
xpym = ∑i xp*ym;
  xpyp = ∑i xp*yp;
    xpzm = ∑i xp*zm;
      xpzp = ∑i xp*zp;
ymzm = ∑i ym*zm;
  ymzp = ∑i ym*zp;
ypzm = ∑i yp*zm;
  ypzp = ∑i yp*zp;
where
   xm = si[0]-ti[0], xp = si[0]+ti[0],
   ym = si[1]-ti[1], yp = si[1]+ti[1],
   zm = si[2]-ti[2], zp = si[2]+ti[2];

M = [ [xm2+ym2+zm2,   ypzm-ymzp,   xmzp-xpzm,   xpym-xmyp],
      [  ypzm-ymzp, xm2+yp2+zp2,   xmym-xpyp,   xmzm-xpzp],
      [  xmzp-xpzm,   xmym-xpyp, xp2+ym2+zp2,   ymzm-ypzp],
      [  xpym-xmyp,   xmzm-xpzp,   ymzm-ypzp, xp2+yp2+zm2] ];

val_vec = Jacobi(M);  // Eigen-values and -vectors
val = val_vec[0], vec = val_vec[1];
small = 0;
for(i = 1; i < val.length; i ++)
  if(val[i] < val[small]) small=i;
return rotationQ( vec[small] );

Thanks to Arun K for bringing this method to my attention.
Search for [Kearsley superposition] in the [Bib].