# Kabsch algorithm

> Mediated Wiki article. Canonical URL: https://mediated.wiki/source/Kabsch_algorithm
> Markdown URL: https://mediated.wiki/source/Kabsch_algorithm.md
> Source: https://en.wikipedia.org/wiki/Kabsch_algorithm
> Source revision: 1334377094
> License: Creative Commons Attribution-ShareAlike 4.0 International (https://creativecommons.org/licenses/by-sa/4.0/)

{{Short description|Type of algorithm}}
The '''Kabsch algorithm''', also known as the '''Kabsch-Umeyama algorithm''',<ref>{{Cite journal |last1=Lawrence |first1=Jim |last2=Bernal |first2=Javier |last3=Witzgall |first3=Christoph |date=2019-10-09 |title=A Purely Algebraic Justification of the Kabsch-Umeyama Algorithm |url=https://nvlpubs.nist.gov/nistpubs/jres/124/jres.124.028.pdf |journal=Journal of Research of the National Institute of Standards and Technology |language=en |volume=124 |pages=124028 |doi=10.6028/jres.124.028 |issn=2165-7254 |pmc=7340555 |pmid=34877177}}</ref> named after Wolfgang Kabsch and Shinji Umeyama, is a method for calculating the optimal [rotation matrix](/source/rotation_matrix) that minimizes the [RMSD](/source/RMSD) ([root mean square](/source/root_mean_square)d deviation) between two paired sets of points. It is useful for [point-set registration](/source/point-set_registration) in [computer graphics](/source/computer_graphics), and in [cheminformatics](/source/cheminformatics) and [bioinformatics](/source/bioinformatics) to compare molecular and [protein](/source/protein) structures (in particular, see [root-mean-square deviation (bioinformatics)](/source/root-mean-square_deviation_(bioinformatics))).

The algorithm only computes the rotation matrix, but it also requires the computation of a translation vector. When both the translation and rotation are actually performed, the algorithm is sometimes called partial [Procrustes superimposition](/source/Procrustes_superimposition) (see also [orthogonal Procrustes problem](/source/orthogonal_Procrustes_problem)).

== Description ==

Let {{mvar|P}} and {{mvar|Q}} be two sets, each containing {{mvar|N}} points in <math>\mathbb{R}^n</math>. We want to find the transformation from {{mvar|Q}} to {{mvar|P}}.  For simplicity, we will consider the three-dimensional case (<math>n = 3</math>).
The sets {{mvar|P}} and {{mvar|Q}} can each be represented by  {{math|''N'' × 3}} [matrices](/source/matrix_(mathematics)) 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: 

<math display="block">\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}</math>

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](/source/centroid) coincides with the origin of the [coordinate system](/source/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 {{mvar|H}}. In matrix notation,

:<math> H = P^\mathsf{T}Q \, </math>

or, using summation notation,

:<math> H_{ij} = \sum_{k = 1}^N P_{ki} Q_{kj}, </math>

which is a [cross-covariance matrix](/source/cross-covariance) when {{mvar|P}} and {{mvar|Q}} are seen as [data matrices](/source/design_matrix).

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

:<math> R = \left(H^\mathsf{T} H\right)^\frac12 H^{-1}, </math>

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](/source/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}},

:<math> H = U \Sigma V^\mathsf{T} </math>

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

:<math> d = \det\left(U V^\mathsf{T}\right) = \det(U) \det(V).</math>

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

:<math> R = U \begin{pmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & d \end{pmatrix} V^\mathsf{T}. </math>

This {{mvar|R}} minimizes <math>\sum_{k = 1}^N\|R q_k - p_k\|_2^2</math>, where <math>q_k</math> and <math>p_k</math> are rows in {{mvar|Q}} and {{mvar|P}} respectively.

Alternatively, optimal rotation matrix can also be directly evaluated as [quaternion](/source/quaternion).<ref>{{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 }}</ref><ref>{{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}}</ref><ref name="Coutsias2004">{{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 }}</ref><ref name="Petitjean1999">{{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 }}</ref> This alternative description has been used in the development of a rigorous method for removing rigid-body motions from [molecular dynamics](/source/molecular_dynamics) trajectories of flexible molecules.<ref>{{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}}</ref> In 2002 a generalization for the application to probability distributions (continuous or not) was also proposed.<ref name="Petitjean2002">{{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 }}</ref>

=== Generalizations ===

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

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

A [Matlab](/source/Matlab) function is available at http://www.mathworks.com/matlabcentral/fileexchange/25746-kabsch-algorithm

A [Python](/source/Python_(programming_language)) 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](/source/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](/source/Visual_Molecular_Dynamics) uses the Kabsch algorithm for its alignment.

The [FoldX](/source/FoldX) modeling toolsuite incorporates the Kabsch algorithm to measure RMSD between Wild Type and Mutated protein structures.

== See also ==

* [Wahba's Problem](/source/Wahba's_problem)
* [Orthogonal Procrustes problem](/source/Orthogonal_Procrustes_problem)

== References ==
{{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

---
Adapted from the Wikipedia article [Kabsch algorithm](https://en.wikipedia.org/wiki/Kabsch_algorithm) by Wikipedia contributors ([contributor history](https://en.wikipedia.org/wiki/Kabsch_algorithm?action=history)). Available under [Creative Commons Attribution-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-sa/4.0/). Changes may have been made.