Kabsch algorithm

Let {{mvar|P}} and {{mvar|Q}} be two sets, each containing {{mvar|N}} points in $\backslash mathbb\{R\}^n$. We want to find the transformation from {{mvar|Q}} to {{mvar|P}}. For simplicity, we will consider the three-dimensional case ($n\; =\; 3$).

The sets {{mvar|P}} and {{mvar|Q}} can each be represented by {{math|N × 3}} 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:

$$\backslash 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.

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.

The second step consists of calculating a matrix {{mvar|H}}. In matrix notation,

:$H\; =\; P^\backslash mathsf\{T\}Q\; \backslash ,$

or, using summation notation,

:$H\_\{ij\}\; =\; \backslash sum\_\{k\; =\; 1\}^N\; P\_\{ki\}\; Q\_\{kj\},$

which is a cross-covariance matrix when {{mvar|P}} and {{mvar|Q}} are seen as data matrices.

It is possible to calculate the optimal rotation {{mvar|R}} based on the matrix formula

:$R\; =\; \backslash left(H^\backslash mathsf\{T\}\; H\backslash right)^\backslash 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 {{mvar|H}} not having an inverse).

If singular value decomposition (SVD) routines are available the optimal rotation, {{mvar|R}}, can be calculated using the following algorithm.

First, calculate the SVD of the covariance matrix {{mvar|H}},

:$H\; =\; U\; \backslash Sigma\; V^\backslash mathsf\{T\}$

where {{mvar|U}} and {{mvar|V}} are orthogonal and $\backslash Sigma$ is diagonal. Next, record if the orthogonal matrices contain a reflection,

:$d\; =\; \backslash det\backslash left(U\; V^\backslash mathsf\{T\}\backslash right)\; =\; \backslash det(U)\; \backslash det(V).$

Finally, calculate our optimal rotation matrix {{mvar|R}} as

:$R\; =\; U\; \backslash begin\{pmatrix\}$

1 & 0 & 0 \\

0 & 1 & 0 \\

0 & 0 & d \end{pmatrix} V^\mathsf{T}.

This {{mvar|R}} minimizes $\backslash sum\_\{k\; =\; 1\}^N|R\; q\_k\; -\; p\_k|$, where $q\_k$ and $p\_k$ are rows in {{mvar|Q}} and {{mvar|P}} respectively.

Alternatively, optimal rotation matrix can also be directly evaluated as quaternion.{{Cite journal|last=Horn|first=Berthold K. P.|authorlink=Berthold K.P. Horn|date=1987-04-01|title=Closed-form solution of absolute orientation using unit quaternions|journal=Journal of the Optical Society of America A|language=EN|volume=4|issue=4|pages=629|doi=10.1364/josaa.4.000629|bibcode=1987JOSAA...4..629H|issn=1520-8532|citeseerx=10.1.1.68.7320|s2cid=11038004 }}{{Cite journal|last=Kneller|first=Gerald R.|date=1991-05-01|title=Superposition of Molecular Structures using Quaternions|journal=Molecular Simulation|volume=7|issue=1–2|pages=113–119|doi=10.1080/08927029108022453|issn=0892-7022}}{{cite journal |last1=Coutsias |first1=E. A. |last2=Seok |first2=C. |last3=Dill |first3=K. A. | title = Using quaternions to calculate RMSD | journal = J. Comput. Chem. | volume = 25 | issue = 15 | pages = 1849–1857 | year = 2004 | pmid = 15376254 | doi = 10.1002/jcc.20110|s2cid=18224579 }}{{cite journal | last = Petitjean | first = M. | title = On the root mean square quantitative chirality and quantitative symmetry measures | journal = J. Math. Phys. | volume = 40 | issue = 9 | pages = 4587–4595 | year = 1999 | doi = 10.1063/1.532988| bibcode = 1999JMP....40.4587P | url = https://hal.archives-ouvertes.fr/hal-02122820/file/PMP.JMP_1999.pdf }} 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.{{Cite journal|date=2011-08-24|title=Least constraint approach to the extraction of internal motions from molecular dynamics trajectories of flexible macromolecules|journal=J. Chem. Phys.|volume=135|issue=8|pages=084110|doi=10.1063/1.3626275|pmid=21895162|issn=0021-9606|last1=Chevrot|first1=Guillaume|last2=Calligari|first2=Paolo|last3=Hinsen|first3=Konrad|last4=Kneller|first4=Gerald R.|bibcode=2011JChPh.135h4110C}} In 2002 a generalization for the application to probability distributions (continuous or not) was also proposed.{{cite journal | last = Petitjean | first = M. | title = Chiral mixtures | journal = J. Math. Phys. | volume = 43 | issue = 8 | pages = 4147–4157 | year = 2002 | doi = 10.1063/1.1484559| bibcode = 2002JMP....43.4147P | s2cid = 85454709 | url = https://hal.archives-ouvertes.fr/hal-02122882/file/PMP.JMP_2002.pdf }}

The algorithm was described for points in a three-dimensional space. The generalization to {{mvar|D}} dimensions is immediate.

This SVD algorithm is described in more detail at https://web.archive.org/web/20140225050055/http://cnx.org/content/m11608/latest/

A Matlab function is available at http://www.mathworks.com/matlabcentral/fileexchange/25746-kabsch-algorithm

A [https://github.com/oleg-alexandrov/projects/blob/master/eigen/Kabsch.cpp C++ implementation] (and unit test) using Eigen

A Python script is available at https://github.com/charnley/rmsd. Another implementation can be found in

[https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.transform.Rotation.align_vectors.html SciPy].

A free PyMol plugin easily implementing Kabsch is [https://www.pymolwiki.org/index.php/Kabsch]. (This previously linked to CEalign [https://wiki.pymol.org/index.php/Cealign], but this uses the Combinatorial Extension (CE) algorithm.) VMD uses the Kabsch algorithm for its alignment.

The FoldX modeling toolsuite incorporates the Kabsch algorithm to measure RMSD between Wild Type and Mutated protein structures.

• Wahba's Problem
• Orthogonal Procrustes problem

{{reflist}}

• {{cite journal|last=Kabsch|first=Wolfgang|date=1976|title=A solution for the best rotation to relate two sets of vectors|journal=Acta Crystallographica|volume=A32|issue=5|page=922|doi=10.1107/S0567739476001873|bibcode=1976AcCrA..32..922K}}
• With a correction in {{cite journal|last=Kabsch|first=Wolfgang|date=1978|title=A discussion of the solution for the best rotation to relate two sets of vectors|journal=Acta Crystallographica|volume=A34|issue=5|pages=827–828|doi=10.1107/S0567739478001680|bibcode=1978AcCrA..34..827K|doi-access=free}}
• {{cite conference|last1=Lin|first1=Ying-Hung|last2=Chang|first2=Hsun-Chang|last3=Lin|first3=Yaw-Ling|date=December 15–17, 2004|title=A Study on Tools and Algorithms for 3-D Protein Structures Alignment and Comparison|conference=International Computer Symposium|location=Taipei, Taiwan}}
• {{cite journal|last=Umeyama|first=Shinji|date=1991|title=Least-Squares Estimation of Transformation Parameters Between Two Point Patterns|journal=IEEE Trans. Pattern Anal. Mach. Intell.|volume=13|issue=4|pages=376–380|doi=10.1109/34.88573}}

Category:Bioinformatics algorithms