Kabsch–Umeyama algorithm is a method for aligning and comparing the similarity between two sets of points. It finds the optimal translation, rotation and scaling by minimizing the root-mean-square deviation (RMSD) of the point pairs. The algorithm for finding the optimal rotation was first described by Kabsch (1976, 1978). In this article we're taking a look at a similar method proposed by Umeyama (1991) which, in addition to rotation, also supports scaling.

Kabsch–Umeyama algorithm is efficient and easy to implement because it only requires a few matrix operations. It has many applications in different domains including physics simulation. However, the algorithm requires a predefined mapping between the sets of points. If you don't have this mapping, your data is noisy or incomplete, or you need a non-rigid transformation, methods for point set registration may work better for you.

## Example

Before diving into the details, let's test out the algorithm in practise. Below we'll attempt to align the Little Dipper and Big Dipper constellations:

## Algorithm

Inputs of the algorithm are two paired sets of $$n$$ points represented as $$m$$-dimensional vectors. Let's say $$A = (\textbf{a}_1, \dots, \textbf{a}_n)$$ are the reference points which we want to align points $$B = (\textbf{b}_1, \dots, \textbf{b}_n)$$ with. We can write out the points as matrices which looks like the following in 2-dimensional case:

$A = \begin{bmatrix} x_1 & y_1 \\ x_2 & y_2 \\ \vdots & \vdots \\ x_n & y_n \end{bmatrix} \quad B = \begin{bmatrix} i_1 & j_1 \\ i_2 & j_2 \\ \vdots & \vdots \\ i_n & j_n \end{bmatrix}$

The algorithm starts by finding the centroids of the points which are their arithmetic means in this case:

$\pmb\mu_A = \frac{1}{n}\sum_{i=1}^n \textbf a_i \quad \pmb\mu_B = \frac{1}{n}\sum_{i=1}^n \textbf b_i$

We also need to the following variance to calculate scaling:

$\sigma_A^2 = \frac{1}{n}\sum_{i=1}^n ||\textbf a_i - \pmb\mu_A||^2$

Next, we calculate the following covariance matrix:

$H = \frac{1}{n}\sum_{i=1}^n (\textbf{a}_i - \pmb\mu_A)^T(\textbf{b}_i - \pmb\mu_B)$

and compute its singular value decomposition $$H = UDV^T$$.

Reflection is sometimes allowed in addition to translation, rotation and scaling. In order to detect and prevent reflection we need the following matrix:

\begin{aligned} d &= \text{sgn}(\text{det}(U)\,\text{det}(V^T)) \\ S &= \text{diag}(\underbrace{1, \dots, 1}_{m-1}, d) \end{aligned}

Finally, we can calculate the optimal rotation matrix $$R$$ and scale factor $$c$$:

\begin{aligned} R &= USV^T \\ c &= \frac{\sigma_A^2}{\text{tr}(DS)} \\ \end{aligned}

With these at hand, we can calculate the aligned points $$B' = (\textbf{b}'_1, \dots, \textbf{b}'_n)$$ where:

$\textbf{b}'_i = \pmb\mu_A + cR(\textbf{b}_i - \pmb\mu_B)$

Sometimes it's nicer to have a single translation vector $$\textbf{t}$$:

$\textbf{t} = \pmb\mu_A - cR\pmb\mu_B$

which can then be used to define the alignment as:

$\textbf{b}'_i = \textbf{t} + cR\textbf{b}_i$

## Implementation

Translating the previous steps into code should be straightforward with linear algebra libraries or languages with such things built-in. Here's an implementation in Python with NumPy:

def kabsch_umeyama(A, B):
assert A.shape == B.shape
n, m = A.shape

EA = np.mean(A, axis=0)
EB = np.mean(B, axis=0)
VarA = np.mean(np.linalg.norm(A - EA, axis=1) ** 2)

H = ((A - EA).T @ (B - EB)) / n
U, D, VT = np.linalg.svd(H)
d = np.sign(np.linalg.det(U) * np.linalg.det(VT))
S = np.diag( * (m - 1) + [d])

R = U @ S @ VT
c = VarA / np.trace(np.diag(D) @ S)
t = EA - c * R @ EB

return R, c, t


Here's how to align the constellations as shown above:

A = np.array([[ 23, 178],
[ 66, 173],
[ 88, 187],
[119, 202],
[122, 229],
[170, 232],
[179, 199]])
B = np.array([[232, 38],
[208, 32],
[181, 31],
[155, 45],
[142, 33],
[121, 59],
[139, 69]])

R, c, t = kabsch_umeyama(A, B)
# R = [[-0.81034281,  0.58595608]
#      [-0.58595608, -0.81034281]]
# c = 1.46166131
# t = [271.3345951, 396.07800317]

B = np.array([t + c * R @ b for b in B])
# B = [[ 29.08878779, 152.36814188]
#      [ 52.37669337, 180.03008629]
#      [ 83.50028582, 204.33920503]
#      [126.28647155, 210.02515345]
#      [131.40664707, 235.37261559]
#      [178.54823113, 222.56285654]
#      [165.79288328, 195.30194121]]

• Kabsch, W. (1976). "A solution for the best rotation to relate two sets of vectors". Acta Crystallographica. A32 (5): 922–923. doi:10.1107/S0567739476001873.
• Kabsch, W. (1978). "A discussion of the solution for the best rotation to relate two sets of vectors". Acta Crystallographica. A34 (5): 827–828. doi:10.1107/S0567739478001680
• Umeyama, S. (1991). "Least-squares estimation of transformation parameters between two point patterns". IEEE Transactions on Pattern Analysis and Machine Intelligence. 13 (4): 376–380. doi:10.1109/34.88573.