Let and be two sets, each containing points in \mathbb{R}^n. We want to find the transformation from to . For simplicity, we will consider the three-dimensional case (n = 3). The sets and can each be represented by
matrices with the first row containing the coordinates of the first point, the second row containing the coordinates of the second point, and so on, as shown in this matrix: \begin{pmatrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ \vdots & \vdots & \vdots \\ x_N & y_N & z_N \end{pmatrix} The algorithm works in three steps: a
translation, the
computation of a covariance matrix, and
the computation of the optimal rotation matrix.
Translation Both sets of coordinates must be translated first, so that their
centroid coincides with the origin of the
coordinate system. This is done by subtracting the centroid coordinates from the point coordinates.
Computation of the covariance matrix The second step consists of calculating a matrix . In matrix notation, : H = P^\mathsf{T}Q \, or, using summation notation, : H_{ij} = \sum_{k = 1}^N P_{ki} Q_{kj}, which is a
cross-covariance matrix when and are seen as
data matrices.
Computation of the optimal rotation matrix It is possible to calculate the optimal rotation based on the matrix formula : R = \left(H^\mathsf{T} H\right)^\frac12 H^{-1}, but implementing a numerical solution to this formula becomes complicated when all special cases are accounted for (for example, the case of not having an inverse). If
singular value decomposition (SVD) routines are available the optimal rotation, , can be calculated using the following algorithm. First, calculate the SVD of the covariance matrix , : H = U \Sigma V^\mathsf{T} where and are orthogonal and \Sigma is diagonal. Next, record if the orthogonal matrices contain a reflection, : d = \det\left(U V^\mathsf{T}\right) = \det(U) \det(V). Finally, calculate our optimal rotation matrix as : R = U \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & d \end{pmatrix} V^\mathsf{T}. This minimizes \sum_{k = 1}^N\|R q_k - p_k\|_2^2, where q_k and p_k are rows in and respectively. Alternatively, optimal rotation matrix can also be directly evaluated as
quaternion. This alternative description has been used in the development of a rigorous method for removing rigid-body motions from
molecular dynamics trajectories of flexible molecules. In 2002 a generalization for the application to probability distributions (continuous or not) was also proposed.
Generalizations The algorithm was described for points in a three-dimensional space. The generalization to dimensions is immediate. == External links ==