Kabsch Algorithm

The Kabsch algorithm finds the rotation, and with it the translation, that best aligns two sets of 3D points whose correspondences are known, in the least-squares sense and in closed form. Given the same points measured in two coordinate frames, such as features seen by a depth camera from two positions, a model’s keypoints and their measured positions, or an estimated trajectory and its ground truth, it returns the rigid motion between the frames without iteration and without an initial guess. The problem is old and has several names: absolute orientation in photogrammetry (Horn, 1987), the constrained orthogonal Procrustes problem in linear algebra (Lawrence, Bernal and Witzgall, 2019), and, for the rotation-only alignment of vectors, Wahba’s problem (SciPy documentation). It has several closed-form solutions that give the same answer: an SVD of a $3 \times 3$ cross-covariance matrix (Kabsch, 1976, and Arun, Huang and Blostein, 1987, as Lawrence et al. and Lorusso et al. describe them), the eigenvector of a $4 \times 4$ symmetric matrix as a unit quaternion (Horn, 1987), and a polar decomposition (Horn, Hilden and Negahdaripour, 1988). Umeyama (1991) added a scale factor, giving the similarity transform that monocular trajectories and reconstructions need; Lawrence et al. call the SVD form the Kabsch–Umeyama algorithm. It is the inner step of point-to-point iterative closest point, which supplies the correspondences by iteration; this article covers the closed-form step itself.

Problem

Let $\mathbf{p}_1, \dots, \mathbf{p}_n$ be points measured in frame A and $\mathbf{q}_1, \dots, \mathbf{q}_n$ the same points measured in frame B, with $\mathbf{p}_i$ and $\mathbf{q}_i$ the same physical point. Find the rotation $R \in SO(3)$, translation $\mathbf{t}$ and, optionally, scale $s > 0$ that minimize the weighted sum of squared residuals

$$E(R, \mathbf{t}, s) = \sum_{i=1}^{n} w_i \,\lVert s R\mathbf{p}_i + \mathbf{t} - \mathbf{q}_i \rVert^2 ,$$

with weights $w_i > 0$, all equal in the basic algorithm. With $s = 1$ this is the rigid, or Euclidean, fit; with $s$ free it is the similarity fit of the geometric transformation hierarchy. The constraint $\det R = +1$ matters: without it the minimizer can be a reflection, an orthogonal matrix with determinant $-1$, which no physical motion produces. In the notation of camera pose, the result maps frame-A coordinates into frame B, $\mathbf{x}_B = R\mathbf{x}_A + \mathbf{t}$; swapping the roles of the two sets gives the inverse.

Horn (1987) counts the constraints: a rigid motion has six degrees of freedom and a similarity seven, three points known in both frames give nine equations, and two points are not enough, since the rotation about the line through them stays free. With more than three points the equations are inconsistent under noise, and the least-squares fit uses all of them rather than discarding information, the defect Horn found in the earlier photogrammetric methods built on three selected points.

Inputs and Outputs

• Inputs: two $n \times 3$ arrays of corresponding points, $n \ge 3$ and not all on one line; optional weights; whether to estimate a scale.
• Outputs: $R$ (with $\det R = +1$), $\mathbf{t}$, optionally $s$, and the residual error. The singular values, or Horn’s eigenvalues, are worth returning too, because they show how well the rotation is determined.

Intuition

The translation is the easy part: at the optimum the two centroids must coincide, so the translation moves the rotated centroid of the first set onto the centroid of the second. Once both sets are centred, only a rotation (and a scale) about the common centroid remains. The best rotation is the one that makes the rotated vectors $R\mathbf{p}'_i$ point as nearly as possible along their partners $\mathbf{q}'_i$: it maximizes the sum of their dot products. That sum depends on the data only through the nine sums of products collected in a $3 \times 3$ cross-covariance matrix, so the whole problem, however many points there are, reduces to a decision about one $3 \times 3$ matrix. The SVD reads off the rotation that best matches that matrix’s principal directions; the quaternion method turns the same maximization into a $4 \times 4$ eigenvalue problem.

Algorithm

With weights normalized to sum to one:

1. Centroids: $\bar{\mathbf{p}} = \sum_i w_i \mathbf{p}_i$, $\bar{\mathbf{q}} = \sum_i w_i \mathbf{q}_i$; centre the points, $\mathbf{p}'_i = \mathbf{p}_i - \bar{\mathbf{p}}$ and $\mathbf{q}'_i = \mathbf{q}_i - \bar{\mathbf{q}}$.
2. Cross-covariance: $H = \sum_i w_i \mathbf{p}'_i \mathbf{q}'^\top_i$.
3. SVD: $H = U S V^\top$, $S = \operatorname{diag}(\sigma_1, \sigma_2, \sigma_3)$ with $\sigma_1 \ge \sigma_2 \ge \sigma_3 \ge 0$.
4. Reflection test: $d = \operatorname{sign}\det(VU^\top)$, that is $+1$ or $-1$.
5. Rotation: $R = V \operatorname{diag}(1, 1, d)\, U^\top$.
6. Scale (optional, Umeyama): $s = (\sigma_1 + \sigma_2 + d\,\sigma_3) \big/ \sum_i w_i \lVert \mathbf{p}'_i \rVert^2$; otherwise $s = 1$.
7. Translation: $\mathbf{t} = \bar{\mathbf{q}} - sR\bar{\mathbf{p}}$.

Horn’s method replaces steps 3–5 by forming a symmetric $4 \times 4$ matrix from the entries of $H$ and taking its eigenvector for the largest eigenvalue as a unit quaternion.

Mathematical Formulation

Translation

For fixed $R$ and $s$, setting the derivative of $E$ with respect to $\mathbf{t}$ to zero gives $\mathbf{t} = \bar{\mathbf{q}} - sR\bar{\mathbf{p}}$ with weighted centroids. Horn (1987, §2.C) derives it without calculus: in centred coordinates the error is $\sum_i \lVert \mathbf{q}'_i - sR\mathbf{p}'_i - \mathbf{t}' \rVert^2$ with $\mathbf{t}' = \mathbf{t} - \bar{\mathbf{q}} + sR\bar{\mathbf{p}}$, the cross term vanishes because centred coordinates sum to zero, and what remains is smallest at $\mathbf{t}' = \mathbf{0}$. The weighted version, with weighted centroids and a weighted cross-covariance, is his Appendix A2.

Rotation: the SVD solution

With the translation eliminated and $s = 1$, expanding the squares gives

$$\sum_i w_i \lVert R\mathbf{p}'_i - \mathbf{q}'_i \rVert^2 = \sum_i w_i \lVert \mathbf{p}'_i \rVert^2 + \sum_i w_i \lVert \mathbf{q}'_i \rVert^2 - 2 \operatorname{tr}(RH),$$

because rotation preserves lengths and $\sum_i w_i \mathbf{q}'^\top_i R\mathbf{p}'_i = \operatorname{tr}(R H)$. Minimizing the error means maximizing $\operatorname{tr}(RH)$. With $H = USV^\top$, $\operatorname{tr}(RH) = \operatorname{tr}(ZS) = \sum_j z_{jj}\sigma_j$ for the orthogonal matrix $Z = V^\top R U$. The diagonal entries of an orthogonal matrix are at most 1, so the trace is at most $\sigma_1 + \sigma_2 + \sigma_3$, reached at $Z = I$, that is $R = VU^\top$. This is the solution of the orthogonal Procrustes problem, which the singular value decomposition article introduces as a projection onto orthogonal matrices; Lorusso, Eggert and Fisher (1995) note that it had been known for some time, citing Schönemann (1966).

The reflection correction

$VU^\top$ is orthogonal but not always a rotation: $\det(VU^\top) = \det V \det U$ is $-1$ in some cases. Then every rotation $R$ gives a $Z = V^\top R U$ with $\det Z = -1$, and for such $Z$ the trace $\operatorname{tr}(ZS)$ is at most $\sigma_1 + \sigma_2 - \sigma_3$, reached at $Z = \operatorname{diag}(1, 1, -1)$, because $\sigma_3$ is the smallest singular value. Lawrence, Bernal and Witzgall (2019) prove this bound with linear algebra alone (their Proposition 1 and Theorem 1), noting that the published justifications by Kabsch and by Umeyama rest on Lagrange multipliers. The best rotation is therefore

$$R = V \operatorname{diag}\!\left(1,\, 1,\, d\right) U^\top, \qquad d = \operatorname{sign}\det(VU^\top),$$

and the minimum of the error is $\sum_i w_i \lVert \mathbf{p}'_i \rVert^2 + \sum_i w_i \lVert \mathbf{q}'_i \rVert^2 - 2(\sigma_1 + \sigma_2 + d\,\sigma_3)$. Lawrence et al. state the algorithm in this form and cite Kabsch’s two papers and Umeyama for it.

When does the unconstrained solution come out as a reflection? Lorusso et al. (1995), describing Arun et al.’s method, name two cases: planar point sets, and large amounts of noise. For the planar case Arun et al. negate the third column of $V$, the one belonging to the zero singular value. For coplanar data that fix costs nothing, since $\sigma_3 = 0$ and the trace is unchanged, and the result is the least-squares rotation. For non-planar data the reflection wins only when the data are severely corrupted, as Horn et al. (1988) put it, and the corrected solution is then the best rotation, but a poor fit. Horn et al. note that, in the non-singular case, their solution is a rotation when the determinant of the cross-covariance matrix is positive and a reflection when it is negative; since $\det H = \det U \det S \det V^\top$ and $\det S > 0$, this is the same test as long as $\sigma_3 > 0$. For a matrix of rank 2, $\det H$ is zero and its computed sign is rounding noise, so they handle that case separately, choosing a sign that makes the determinant of the result positive. scikit-image’s implementation of Umeyama’s method, whose comments cite Umeyama’s equations (38)–(43), does the same: it uses the sign of $\det H$, but switches to the sign of $\det U \det V$ when the matrix has rank one less than the dimension, rank 2 in 3D (0.25.2 source). Zhang and Scaramuzza (2018, Algorithm 1) and Lawrence et al. use the $\det U \det V$ test throughout, which covers both cases. The example below shows what goes wrong with the $\det H$ test alone.

Scale: Umeyama’s similarity

With a scale, the centred error is $s^2 \sum_i w_i \lVert \mathbf{p}'_i \rVert^2 - 2s \operatorname{tr}(RH) + \sum_i w_i \lVert \mathbf{q}'_i \rVert^2$. For any $s > 0$ the best $R$ is still the one that maximizes $\operatorname{tr}(RH)$, so the rotation does not depend on the scale, and minimizing over $s$ gives

$$s = \frac{\sigma_1 + \sigma_2 + d\,\sigma_3}{\sum_i w_i \lVert \mathbf{p}'_i \rVert^2}, \qquad \mathbf{t} = \bar{\mathbf{q}} - sR\bar{\mathbf{p}} .$$

This is Umeyama’s (1991) result as Zhang and Scaramuzza (2018, Algorithm 1) give it, there with $\Sigma = \frac{1}{n}\sum_i \mathbf{q}'_i\mathbf{p}'^\top_i$, the transpose of $H$, so their $U$ and $V$ are our $V$ and $U$. Horn (1987, §2.D) derived the same scale, $s = D/S_l$ in his notation, with $D$ the sum of dot products of the right points and the rotated left points and $S_l$ the sum of squared lengths of the centred left points. He also noted its drawback (§2.E): fitting the inverse transformation and inverting does not give the same scale. He proposed a symmetric scale, $s = (\sum_i \lVert \mathbf{q}'_i \rVert^2 / \sum_i \lVert \mathbf{p}'_i \rVert^2)^{1/2}$, the ratio of the RMS spreads, which needs no rotation, and recommended it when both sets are measured with similar errors; an asymmetric scale suits the case where one set is far more precise than the other. The choice does not change the optimal rotation. Absolute trajectory error explains which of these scales the trajectory metric uses.

Horn’s unit-quaternion method

Horn (1987) represents the rotation by a unit quaternion $\mathring{q} = (q_0, q_x, q_y, q_z)$, scalar first, as in rotation representations, and writes the quantity to maximize, $\sum_i \mathbf{q}'_i \cdot R\mathbf{p}'_i$, as the quadratic form $\mathring{q}^\top N \mathring{q}$. With the entries of the cross-covariance $M = H$ written $S_{xy} = \sum_i x'_{p,i}\, y'_{q,i}$ and so on (Horn’s left set is our $\mathbf{p}$, his right set our $\mathbf{q}$),

$$N = \begin{bmatrix} S_{xx} + S_{yy} + S_{zz} & S_{yz} - S_{zy} & S_{zx} - S_{xz} & S_{xy} - S_{yx} \\ S_{yz} - S_{zy} & S_{xx} - S_{yy} - S_{zz} & S_{xy} + S_{yx} & S_{zx} + S_{xz} \\ S_{zx} - S_{xz} & S_{xy} + S_{yx} & -S_{xx} + S_{yy} - S_{zz} & S_{yz} + S_{zy} \\ S_{xy} - S_{yx} & S_{zx} + S_{xz} & S_{yz} + S_{zy} & -S_{xx} - S_{yy} + S_{zz} \end{bmatrix} .$$

$N$ is symmetric with zero trace. The unit vector that maximizes $\mathring{q}^\top N\mathring{q}$ is the eigenvector of the most positive eigenvalue (his Appendix A3), and the solution is unique if that eigenvalue is distinct. Because the composite product $\mathring{q}\,\mathring{r}\,\mathring{q}^*$ preserves cross products, a unit quaternion always represents a rotation, never a reflection (§3.E and Appendix A6), so the method needs no determinant test. Horn computed the eigenvalues from the characteristic quartic: its cubic coefficient is zero (the trace), its quadratic coefficient is $-2\operatorname{tr}(M^\top M)$, and its linear coefficient is $-8\det M$, which vanishes when either set is coplanar, leaving a quadratic in $\lambda^2$ (§§4.D, 5). Today a library eigensolver for symmetric matrices does the same job. The quaternion converts to a matrix with the formula of his §3.E, and $\mathring{q}$ and $-\mathring{q}$ give the same rotation.

The two methods agree exactly in exact arithmetic. It can be checked that the eigenvalues of $N$ are $\sigma_1 + \sigma_2 + d\sigma_3$, $\sigma_1 - \sigma_2 - d\sigma_3$, $-\sigma_1 + \sigma_2 - d\sigma_3$ and $-\sigma_1 - \sigma_2 + d\sigma_3$; the largest equals the maximal trace of the SVD solution, and the gap to the next one is $2(\sigma_2 + d\sigma_3)$. This is our derivation, not a statement of Horn’s; it shows that the rotation is unique exactly when $\sigma_2 + d\sigma_3 > 0$, which fails for collinear points ($\sigma_2 = \sigma_3 = 0$).

Orthonormal matrices and the polar decomposition

Horn, Hilden and Negahdaripour (1988) solved the same problem directly with $3 \times 3$ matrices: the rotation is the orthonormal factor of the polar decomposition of the cross-covariance matrix, which they compute from the eigenvalues and eigenvectors of a symmetric $3 \times 3$ product of that matrix with its transpose, through its inverse square root. Lorusso et al. (1995) write it in the notation above as $R = H^\top \left( \sum_k \lambda_k^{-1/2} \mathbf{u}_k \mathbf{u}_k^\top \right)$, with $\lambda_k$ and $\mathbf{u}_k$ the eigenvalues and eigenvectors of $HH^\top$; it equals $VU^\top$ when $H$ is non-singular. Horn et al. note that the product is a rotation when the determinant of the cross-covariance is positive and a reflection when it is negative, give a separate formula for rank 2, the coplanar case, and observe that when the points are collinear the problem has no unique solution. They also anticipate numerical problems when one of the eigenvalues is nearly zero, as happens when one set lies almost in a plane. Lorusso et al. found the proposed rank-2 correction not valid in all cases and give an alternative. In the same paper Horn et al. note that Arun et al. had independently solved an equivalent problem with an SVD, and they compare their solution with the older approach of fitting an unconstrained linear map and taking the nearest orthonormal matrix, which does not solve the original least-squares problem; the SVD article’s example measures that difference.

Weights and noise

Least squares is the maximum-likelihood estimate when the residuals are independent, isotropic and Gaussian with equal variances; weights $w_i = 1/\sigma_i^2$ extend it to points with different but still isotropic noise. It follows, by our reasoning, that the rotation stays the maximum-likelihood one when both sets carry isotropic noise, since a rotation leaves an isotropic covariance unchanged and the residual $R\mathbf{p}_i + \mathbf{t} - \mathbf{q}_i$ then has variance $\sigma_{p,i}^2 + \sigma_{q,i}^2$ in every direction. Real 3D measurements are rarely isotropic: points from stereo or a depth camera are far less certain along the viewing ray, with depth error growing with the square of depth, as disparity explains. Scalar weights capture the change of uncertainty with depth but not its direction; a full covariance per point leaves the closed form behind and is usually handled with nonlinear least squares. The visual odometry article reports the related finding of Nistér et al. (2004), as Scaramuzza and Fraundorfer describe it, that 3D-to-2D motion estimation was greatly superior to 3D–3D alignment for stereo, because triangulated points are much more uncertain in depth.

SciPy’s Rotation.align_vectors solves the weighted rotation-only problem, which it calls the Kabsch algorithm and Wahba’s problem. It aligns vectors, so point sets must be centred first, and it can return a sensitivity matrix which, multiplied by the harmonic mean of the observation variances, approximates the covariance of the rotation error when the weights are inversely proportional to those variances (SciPy 1.16 documentation).

Parameters

• Weights. Uniform by default. Inverse variances when the measurement quality of the points differs, for instance with depth.
• Scale. Fixed at 1 for sensors with metric scale; estimated for monocular reconstructions and trajectories. Which of Horn’s three scales to use depends on which set carries the noise (see above, and the example).
• Reflection test. $\operatorname{sign}\det(VU^\top)$, or equivalently $\det U\det V$, is safe for all non-collinear data; the sign of $\det H$ needs a separate rank-deficient branch.
• Direction. Fitting $\mathbf{p} \to \mathbf{q}$ or $\mathbf{q} \to \mathbf{p}$ gives exactly inverse rigid transforms, but not inverse scales unless the symmetric scale is used (Horn, 1987, §2.E).

Complexity

Forming the centroids and the nine sums of $H$ takes $O(n)$ time and $O(1)$ extra memory; the SVD of a $3 \times 3$ matrix or the eigen-decomposition of the $4 \times 4$ matrix $N$ costs a constant. Lorusso et al. (1995) note that all four methods they compared are $O(n)$, and in their 1995 timings on a SPARCstation 10 the polar-decomposition method was the quickest for fewer than about 100 points while the differences between the SVD, quaternion and polar methods became negligible, about 1%, for large sets; the cost of the decomposition, not the sums, sets the constant, so the choice among the methods is not a matter of speed.

Implementation

The example builds every point set from simulated RGB-D image measurements. Points of a known scene are projected into two depth cameras ($f = 500$ px, $640 \times 480$), the second about 0.9 m from the first and turned toward the scene; each point’s pixel is rounded to the pixel grid, its depth receives Gaussian noise with standard deviation $1.5\,\text{mm} \times (z / 1\,\text{m})^2$, a modelling choice that imitates depth from disparity, and it is back-projected through the rounded pixel. The task is the relative pose $\mathbf{x}_2 = R\mathbf{x}_1 + \mathbf{t}$ of the two cameras from 40 known correspondences. The script compares Kabsch (with Umeyama’s scale) and Horn’s quaternion method on a general scene (1); counts reflections on a planar poster, tests the two determinant rules on exactly planar data and simulates a handedness bug (2); shrinks a pole toward a line (3); compares unweighted and inverse-variance weighted fits on points 1–6 m deep (4); and compares Horn’s three scales on a small object whose second copy is in units half as large (5). It uses NumPy 2.3.5 and SciPy 1.16.3 and runs in under a second. After the first run, two printed lines were changed without changing any computation: part 2b’s counts now print as plain integers rather than in NumPy’s np.int64(...) form, and part 5, whose label promised the standard deviation of each scale but printed only the forward one’s, now prints all three.

import numpy as np
from scipy.spatial.transform import Rotation

rng = np.random.default_rng(0)
f, cx, cy, W, H = 500.0, 319.5, 239.5, 640, 480         # RGB-D camera (pixels)

def look_at(C, target, roll_deg=0.0):
    """World-to-camera (R, t), x right, y down, z forward."""
    z = (target - C) / np.linalg.norm(target - C)
    x = np.cross([0.0, 1.0, 0.0], z); x /= np.linalg.norm(x)
    R = Rotation.from_euler("z", roll_deg, degrees=True).as_matrix() @ np.stack([x, np.cross(z, x), z])
    return R, -R @ C

# Camera 1 defines the world frame; camera 2 stands 0.9 m from camera 1 and looks at a point 2.5 m ahead.
R2, t2 = look_at(np.array([0.8, -0.2, 0.3]), np.array([0.0, 0.0, 2.5]), roll_deg=5.0)
R_true, t_true = R2, t2                                   # X_c2 = R_true X_c1 + t_true

def pixels(X):
    return f * X[:, 0] / X[:, 2] + cx, f * X[:, 1] / X[:, 2] + cy

def in_both(X):
    ok = np.ones(len(X), bool)
    for Xc in (X, X @ R2.T + t2):
        u, v = pixels(Xc)
        ok &= (Xc[:, 2] > 0) & (u >= 0) & (u < W) & (v >= 0) & (v < H)
    return ok

def draw_visible(sampler, n=40):
    X = np.empty((0, 3))
    while len(X) < n:
        Y = sampler(n); X = np.vstack([X, Y[in_both(Y)]])
    return X[:n]

def sigma_z(z):
    return 1.5e-3 * z ** 2                                # depth noise (m), grows with z^2

def measure(Xc):
    """RGB-D measurement: pixel rounded to the grid, noisy depth, back-projected."""
    u, v = np.round(pixels(Xc)[0]), np.round(pixels(Xc)[1])
    z = Xc[:, 2] + rng.normal(0, sigma_z(Xc[:, 2]))
    return np.column_stack([(u - cx) * z / f, (v - cy) * z / f, z])

def observe(X):
    return measure(X), measure(X @ R2.T + t2)             # same points seen in both cameras

def kabsch(P, Q, w=None, rule="det(VU^T)"):
    """R, t minimizing sum w_i ||R p_i + t - q_i||^2; also Umeyama's scale and its translation."""
    w = np.ones(len(P)) / len(P) if w is None else w / w.sum()
    p0, q0 = w @ P, w @ Q
    Hm = (P - p0).T @ ((Q - q0) * w[:, None])             # H = sum w p' q'^T = U S V^T
    U, S, Vt = np.linalg.svd(Hm)
    if rule == "none":
        d = 1.0                                           # plain orthogonal Procrustes
    elif rule == "det(H)":
        d = -1.0 if np.linalg.det(Hm) < 0 else 1.0
    else:
        d = -1.0 if np.linalg.det(Vt.T @ U.T) < 0 else 1.0
    R = Vt.T @ np.diag([1.0, 1.0, d]) @ U.T
    s = (S[0] + S[1] + d * S[2]) / (w @ np.sum((P - p0) ** 2, axis=1))
    return R, q0 - R @ p0, s

def horn(P, Q, w=None):
    """Horn (1987): unit quaternion (w, x, y, z) = top eigenvector of the 4x4 matrix N."""
    w = np.ones(len(P)) / len(P) if w is None else w / w.sum()
    p0, q0 = w @ P, w @ Q
    (Sxx, Sxy, Sxz), (Syx, Syy, Syz), (Szx, Szy, Szz) = (P - p0).T @ ((Q - q0) * w[:, None])
    N = np.array([[Sxx + Syy + Szz, Syz - Szy, Szx - Sxz, Sxy - Syx],
                  [Syz - Szy, Sxx - Syy - Szz, Sxy + Syx, Szx + Sxz],
                  [Szx - Sxz, Sxy + Syx, -Sxx + Syy - Szz, Syz + Szy],
                  [Sxy - Syx, Szx + Sxz, Syz + Szy, -Sxx - Syy + Szz]])
    lam, E = np.linalg.eigh(N)                            # ascending eigenvalues
    a, b, c, d = E[:, -1]
    R = np.array([[a*a + b*b - c*c - d*d, 2 * (b*c - a*d), 2 * (b*d + a*c)],
                  [2 * (c*b + a*d), a*a - b*b + c*c - d*d, 2 * (c*d - a*b)],
                  [2 * (d*b - a*c), 2 * (d*c + a*b), a*a - b*b - c*c + d*d]])
    return R, q0 - R @ p0, lam

def ang(Ra, Rb):                                          # geodesic angle in degrees
    return np.degrees(Rotation.from_matrix(Ra @ Rb.T).magnitude())

# 1. A general 3D scene: Kabsch, Horn and Umeyama on the same noisy data (200 draws).
X = draw_visible(lambda n: np.column_stack([rng.uniform(-0.8, 0.8, n), rng.uniform(-0.6, 0.6, n),
                                            rng.uniform(1.5, 3.5, n)]))
res = []
for _ in range(200):
    P, Q = observe(X)
    Rk, tk, s = kabsch(P, Q)
    Rh, th, _ = horn(P, Q)
    res.append([ang(Rk, Rh), np.linalg.norm(tk - th), ang(Rk, R_true), np.linalg.norm(tk - t_true), s])
res = np.array(res)
print("1. volume scene, 200 draws")
print("   max |Kabsch - Horn|: %.1e deg, %.1e m" % (res[:, 0].max(), res[:, 1].max()))
print("   median error: rotation %.3f deg, translation %.1f mm; Umeyama scale %.4f to %.4f"
      % (np.median(res[:, 2]), 1000 * np.median(res[:, 3]), res[:, 4].min(), res[:, 4].max()))

# 2a. A planar poster: how often is the uncorrected Procrustes solution a reflection?
xy = np.column_stack([rng.uniform(-0.6, 0.6, 40), rng.uniform(-0.45, 0.45, 40), np.zeros(40)])
poster = xy @ Rotation.from_euler("y", 20, degrees=True).as_matrix().T + [0.1, 0.0, 2.5]
refl, err_c = 0, []
for _ in range(1000):
    P, Q = observe(poster)
    refl += np.linalg.det(kabsch(P, Q, rule="none")[0]) < 0
    err_c.append(ang(kabsch(P, Q)[0], R_true))
print("2a. planar poster, 1000 draws: uncorrected V U^T is a reflection in %d; "
      "corrected median error %.3f deg" % (refl, np.median(err_c)))

# 2b. Exactly planar, noise-free data: which determinant test is safe?
bad = {"det(H)": 0, "det(VU^T)": 0}
for k in range(1000):
    Pp = rng.uniform(-1, 1, (40, 3)) * [1, 1, 0] @ Rotation.random(random_state=k).as_matrix().T
    Qp = Pp @ Rotation.random(random_state=1000 + k).as_matrix().T + rng.uniform(-1, 1, 3)
    for rule in bad:
        bad[rule] += int(np.linalg.det(kabsch(Pp, Qp, rule=rule)[0]) < 0)
print("2b. exactly planar, 1000 configurations: reflections returned with test", bad)

# 2c. A handedness bug: camera 2's points written with y up instead of y down.
P, Q = observe(X)
Qf = Q * [1, -1, 1]
for rule in ["det(VU^T)", "none"]:
    R, t, _ = kabsch(P, Qf, rule=rule)
    rms = np.sqrt(np.mean(np.sum((P @ R.T + t - Qf) ** 2, axis=1)))
    print("2c. y flipped, test %-9s det(R) = %+.0f, RMS residual %.1f mm" % (rule, np.linalg.det(R), 1000 * rms))
R, t, _ = kabsch(P, Q)
print("    correct data:            RMS residual %.1f mm" % (1000 * np.sqrt(np.mean(np.sum((P @ R.T + t - Q) ** 2, axis=1)))))

# 3. Near-collinear points: a pole 1.55 m long whose cross-section radius r shrinks.
a0, a1 = np.array([-0.5, -0.6, 2.3]), np.array([0.4, 0.6, 2.7])
axis = (a1 - a0) / np.linalg.norm(a1 - a0)
e1 = np.cross(axis, [1.0, 0, 0]); e1 /= np.linalg.norm(e1); e2 = np.cross(axis, e1)
axis2 = R_true @ axis                                     # pole direction in camera 2
print("3. pole, 200 draws per radius: median errors (deg) and Horn's eigenvalue gap")
for r in [0.3, 0.1, 0.03, 0.01, 0.003, 0.0]:
    rad, phi = r * np.sqrt(rng.uniform(0, 1, 40)), rng.uniform(0, 2 * np.pi, 40)
    pole = (a0 + np.outer(rng.uniform(0, 1, 40), a1 - a0)
            + np.outer(rad * np.cos(phi), e1) + np.outer(rad * np.sin(phi), e2))
    out = []
    for _ in range(200):
        P, Q = observe(pole)
        Rk, Rh = kabsch(P, Q)[0], horn(P, Q)[0]
        lam = horn(P, Q)[2]
        wvec = Rotation.from_matrix(Rk @ R_true.T).as_rotvec()
        about = abs(wvec @ axis2)
        out.append([ang(Rk, R_true), ang(Rh, R_true), np.degrees(about),
                    np.degrees(np.linalg.norm(wvec - (wvec @ axis2) * axis2)),
                    (lam[-1] - lam[-2]) / lam[-1], ang(Rk, Rh)])
    m = np.median(out, axis=0)
    print("   r = %5.3f m: Kabsch %7.3f  Horn %7.3f  (about axis %7.3f, across %.3f)  gap %.1e  max |K-H| %.0e"
          % (r, m[0], m[1], m[2], m[3], m[4], np.max(np.array(out)[:, 5])))

# 4. Weighting: points from 1 to 6 m deep, where depth noise grows 36-fold.
def frustum(n):
    z, u, v = rng.uniform(1, 6, n), rng.uniform(0, W, n), rng.uniform(0, H, n)
    return np.column_stack([(u - cx) * z / f, (v - cy) * z / f, z])
Xd = draw_visible(frustum)
wts = 1.0 / (sigma_z(Xd[:, 2]) ** 2 + sigma_z((Xd @ R2.T + t2)[:, 2]) ** 2)
out = []
for _ in range(200):
    P, Q = observe(Xd)
    for w in (None, wts):
        R, t, _ = kabsch(P, Q, w)
        out.append([ang(R, R_true), 1000 * np.linalg.norm(t - t_true)])
out = np.array(out).reshape(200, 2, 2)
med = np.median(out, axis=0)
print("4. depths 1-6 m, 200 draws: median error unweighted %.3f deg / %.1f mm, "
      "weighted %.3f deg / %.1f mm; weighted better in %d of 200 (rotation)"
      % (med[0, 0], med[0, 1], med[1, 0], med[1, 1], np.sum(out[:, 1, 0] < out[:, 0, 0])))

# 5. Scale: a 20 cm object 3 m away; the second set is in units half as large (s = 0.5).
cube = rng.uniform(-0.1, 0.1, (40, 3)) + [0.1, 0.05, 3.0]
sc = []
for _ in range(500):
    P, Q = observe(cube)
    Q = 0.5 * Q
    s_fwd = kabsch(P, Q)[2]                               # Umeyama: scale P onto Q
    s_inv = 1.0 / kabsch(Q, P)[2]                         # reverse fit, inverted
    s_sym = np.sqrt(np.sum((Q - Q.mean(0)) ** 2) / np.sum((P - P.mean(0)) ** 2))   # Horn's symmetric
    sc.append([s_fwd, s_inv, s_sym])
sc = 100 * (np.array(sc) / 0.5 - 1)
print("5. scale, true 0.5, 500 draws: mean error forward %+.2f%%, inverse of reverse %+.2f%%, "
      "symmetric %+.2f%% (s.d. %.2f%%, %.2f%%, %.2f%%)" % (*sc.mean(0), *sc.std(0)))

Output:

1. volume scene, 200 draws
   max |Kabsch - Horn|: 3.1e-14 deg, 4.5e-15 m
   median error: rotation 0.128 deg, translation 6.6 mm; Umeyama scale 0.9943 to 1.0059
2a. planar poster, 1000 draws: uncorrected V U^T is a reflection in 516; corrected median error 0.433 deg
2b. exactly planar, 1000 configurations: reflections returned with test {'det(H)': 278, 'det(VU^T)': 0}
2c. y flipped, test det(VU^T) det(R) = +1, RMS residual 554.4 mm
2c. y flipped, test none      det(R) = -1, RMS residual 18.5 mm
    correct data:            RMS residual 18.5 mm
3. pole, 200 draws per radius: median errors (deg) and Horn's eigenvalue gap
   r = 0.300 m: Kabsch   0.363  Horn   0.363  (about axis   0.257, across 0.176)  gap 3.5e-01  max |K-H| 5e-14
   r = 0.100 m: Kabsch   0.771  Horn   0.771  (about axis   0.704, across 0.184)  gap 6.5e-02  max |K-H| 6e-13
   r = 0.030 m: Kabsch   2.219  Horn   2.219  (about axis   2.217, across 0.170)  gap 4.2e-03  max |K-H| 7e-12
   r = 0.010 m: Kabsch   8.485  Horn   8.485  (about axis   8.478, across 0.212)  gap 6.5e-04  max |K-H| 7e-11
   r = 0.003 m: Kabsch  33.923  Horn  33.923  (about axis  33.920, across 0.220)  gap 1.0e-04  max |K-H| 3e-09
   r = 0.000 m: Kabsch 122.592  Horn 122.592  (about axis 122.591, across 0.227)  gap 8.0e-05  max |K-H| 7e-09
4. depths 1-6 m, 200 draws: median error unweighted 0.131 deg / 12.1 mm, weighted 0.078 deg / 3.1 mm; weighted better in 157 of 200 (rotation)
5. scale, true 0.5, 500 draws: mean error forward -1.80%, inverse of reverse +1.35%, symmetric -0.24% (s.d. 1.77%, 1.83%, 1.77%)

Agreement. Kabsch’s SVD and Horn’s eigenvector give the same rotation to within $10^{-13}$ degrees in all 200 draws, as they must: they maximize the same trace. Lorusso et al. (1995) reached the same conclusion experimentally, finding differences in accuracy many orders of magnitude below normal noise levels. The errors that remain, about a tenth of a degree and a few millimetres, come from the measurements, not from the solvers; the Umeyama scale, estimated although the true scale is 1, stays within 0.6% of it.

Reflections. On the poster the third singular value of $H$ is made of noise alone, so noise decides whether $VU^\top$ is a reflection: about half the draws (516 of 1000) give $\det = -1$, the planar case of Arun et al. After the correction every draw yields a rotation, with a median error of 0.43°, larger than on the volume scene; in our reading because the rotation about the poster’s in-plane axes is pinned only by the spread across a 1.2 m by 0.9 m plane, with no depth extent. On exactly planar, noise-free data the $\det H$ rule returns a reflection in 278 of 1000 configurations: $\det H$ is a rounding residue whose sign is meaningless, which is why Horn et al. (1988) and scikit-image treat rank 2 separately, while the $\det(VU^\top)$ rule never fails. The handedness bug shows the other side: when camera 2’s points are written with $y$ up instead of $y$ down, the corrected algorithm still returns a rotation, the best one, but its RMS residual of 554 mm is thirty times the residual on correct data. The uncorrected orthogonal fit matches the flipped data as well as the correct fit matches correct data (18.5 mm) and has $\det = -1$. In our reading, a correction that silently turns a reflection into a rotation should therefore be paired with a check of the residual, or of the uncorrected determinant, whenever the frames come from different software.

Near-collinear points. As the pole’s radius shrinks from 30 cm to 1 cm, the median rotation error grows from 0.36° to 8.5°, and almost all of it is a rotation about the pole’s axis: the error across the axis stays near 0.2° throughout, in our reading set by the 1.55 m length. At $r = 3$ mm, below the depth noise of the construction at 2.5 m (9.4 mm), the angle about the axis is essentially guessed, and for exactly collinear points the median error is 123°. Horn’s relative eigenvalue gap shrinks with the radius, and the two methods drift apart (to about $10^{-8}$ degrees) only as the problem becomes ill-conditioned, still negligibly. Lorusso et al. (1995) observed the same breakdown in their 3-D-to-1-D experiments: under isotropic noise the error rose as the line’s diameter fell and, in the best cases, broke down completely when the diameter reached the noise level. The small eigenvalue gap or the ratio $\sigma_2/\sigma_1$ is the warning sign to check before trusting a rotation.

Weighting. With points from 1 to 6 m deep the depth noise of the construction varies 36-fold in standard deviation, from 1.5 to 54 mm, and the unweighted fit lets the far, noisy points pull the estimate. Inverse-variance weights cut the median rotation error from 0.131° to 0.078° and the translation error from 12.1 to 3.1 mm, and give the smaller rotation error in 157 of 200 draws. The weights use only the depth variance; the remaining anisotropy is ignored.

Scale. All three scales are noisy, with a spread of about 1.8% per draw, but their means differ systematically. Umeyama’s forward scale is low by 1.8% on average, and the inverse of the reverse fit is high by 1.35%. In our reading this is the familiar bias of fitting with noise in the regressor: noise inflates the spread $\sum \lVert \mathbf{p}'_i \rVert^2$ in the denominator of the forward scale, while the reverse fit’s denominator is the spread of the second set. The second camera is closer to the object than the first, so its depth noise is smaller and the reverse bias is the smaller one. Horn’s symmetric scale, the ratio of the two spreads, is biased by roughly half the difference of the two (−0.24%), and would be unbiased if the two sets’ noise were in proportion to their scales. None of the three is right in general. By the same reasoning, the forward scale is unbiased when the set being scaled is the precise one, and the inverse of the reverse fit when the target is; Horn’s advice for asymmetric errors leaves open which one to choose. For a metric such as absolute trajectory error, which measures residuals after scaling the estimate onto the ground truth, the forward scale is the one that minimizes those residuals.

Properties and Behavior

• Exactness and global optimality. The solution is the global minimizer of the least-squares error over all rotations (and scales), not a local one; it needs no starting point, which is what Horn (1987) presents as the advantage over the iterative photogrammetric methods of the time.
• Equivalence of the methods. SVD, quaternion and polar-decomposition solutions are the same estimator. In the experiments of Lorusso et al. (1995), on random point sets from 4 to 10,000 points and noise levels up to $\sigma = 1$, the four methods they compared (the three above and Walker, Shao and Volz’s dual quaternions) differed in accuracy by at most $10^{-13}$ in the noise-free case. On degenerate data the SVD and quaternion methods were usually the most stable, the polar-decomposition method was less stable for planar data, and the dual-quaternion method usually broke down first; they concluded that the SVD method gave the best overall accuracy and stability. The journal version of this comparison (Eggert, Lorusso and Fisher, 1997) could not be read for this article; the statements here come from the 1995 conference paper.
• Symmetry. Rigid fits are exactly symmetric: aligning $\mathbf{q}$ to $\mathbf{p}$ gives the inverse of aligning $\mathbf{p}$ to $\mathbf{q}$. Scaled fits are symmetric only with Horn’s symmetric scale (Horn, 1987, §2.E, Appendix A1).
• Degenerate configurations. Collinear points (and any set of fewer than three distinct points) leave the rotation about the line undetermined (Horn et al., 1988). Coplanar points determine the rotation uniquely, but the cross-covariance has rank 2, the reflection ambiguity is decided by noise, and implementations must use a test that survives rank deficiency. Nearly collinear sets give large errors about the line, as in the example. The quantities $\sigma_2 + d\sigma_3$, or Horn’s eigenvalue gap, measure how far a configuration is from degeneracy.
• Exact orthogonality. An estimated $R$ is always exactly orthogonal up to rounding, with no squaring-up step, unlike the unconstrained linear fit that Horn (1987) criticizes; its geodesic error from the true rotation is measured as in rotation representations.

Limitations

• Correspondences must be known. The method aligns given pairs; it does not find them. Without correspondences, iterative closest point guesses them from nearest neighbours and calls this solution at every iteration.
• Not robust to outliers. One wrong correspondence enters the squared error with full weight. A common remedy is to wrap the closed form in RANSAC, fitting it to small random samples of correspondences (three is the minimum); the Tanks and Temples evaluation, for example, fits a similarity to random samples of six camera centres with Umeyama’s method.
• Isotropic noise model. Anisotropic or correlated errors, such as depth errors along viewing rays, make the least-squares solution suboptimal; weights only partly help.
• Ill-conditioning near degeneracy. Long thin or nearly planar configurations give a well-defined but poorly determined answer, which the method returns without warning unless the singular values are inspected.
• Rigid or similarity only. Deforming objects, independently moving parts and lens or calibration errors that warp one set are outside the model.

Variants

• Similarity (Umeyama, 1991). The scaled fit above, used for monocular trajectories and reconstructions; absolute trajectory error and camera pose cover its use in evaluation, including the yaw-only variant of Zhang and Scaramuzza (2018) for visual-inertial estimates.
• Horn’s three-point solutions. For exactly three points Horn (1987, §2.A) gives a triad construction that discards information, and in §5 a least-squares solution for coplanar sets, which first rotates one plane onto the other and then solves a one-parameter rotation within the plane.
• Dual quaternions. Walker, Shao and Volz represent rotation and translation together as a dual quaternion and can include matched normal vectors with their own weights; Lorusso et al. (1995) describe the method and include it in their comparison.
• Weighted vector alignment. Wahba’s problem, the attitude-determination form, aligns weighted direction vectors without translation; SciPy’s align_vectors solves it with the Kabsch algorithm and also allows one infinitely weighted vector that must be matched exactly (SciPy 1.16 documentation).
• Any dimension. The SVD form works unchanged in $d$ dimensions with $\operatorname{diag}(1, \dots, 1, d)$, as Lawrence et al. (2019) state it; scikit-image’s EuclideanTransform and SimilarityTransform estimate 2D and 3D transforms this way.
• Inside other estimators. The closed form is the point-to-point step of ICP; it aligns trajectories for absolute trajectory error on TUM RGB-D and EuRoC MAV; it aligns predicted and true skeletons in the Procrustes-aligned error of 3D human pose estimation; and 3D–3D visual odometry, which aligns points triangulated at consecutive times, poses the same problem. Horn (1987) names further uses: relating a photogrammetric stereo model to geodetic coordinates, and the hand–eye transform between a camera and a manipulator in robotics, where it connects to robotic manipulation. Point cloud registration, hand–eye calibration and photogrammetry build on it in turn.

History

The problem arose independently in several fields. Psychometricians solved the orthogonal Procrustes problem with the SVD in the 1950s and 1960s, as the SVD article recounts. Kabsch (1976) published his solution for the best rotation relating two sets of vectors in a crystallography journal, followed by a discussion of the solution in 1978; neither paper could be read for this article, and its name is attached to the algorithm by later users, such as SciPy and Lawrence et al. (2019). Photogrammetry knew the task as absolute orientation, solved by empirical, graphical and iterative methods or from three selected points (Horn, 1987, §1.A). In 1987 Horn published his unit-quaternion closed form, and Arun, Huang and Blostein their SVD solution; Horn, Hilden and Negahdaripour added the orthonormal-matrix form in 1988, having learned of Arun et al.’s work in the meantime. Umeyama (1991) gave the rotation-guaranteed solution with scale, and Lorusso, Eggert and Fisher (1995) compared four of the methods experimentally.