Procrustes Analysis

1. Background Knowledge

Procrustes analysis (in this rigid case) is about finding the best-fit rigid transform (rotation + translation, no scaling or shearing) that aligns one set of points to another in a least squares sense. You are given corresponding points $\{p_i\}$ and $\{q_i\}$, and you want the rotation matrix $R$ (orthonormal, $R^TR = I$, $\det(R)=1$) and translation vector $t$ that minimize

$$\sum_i \|R p_i + t - q_i\|^2.$$

This is a classic problem in computer vision and shape analysis, used in tasks like shape matching, pose estimation, and image alignment.

The key idea is to decouple translation from rotation: if you subtract the centroid from each point set, you remove translation, leaving only rotation to be estimated. The optimal rotation between two centered point sets can be found via the Singular Value Decomposition (SVD) of a covariance (cross-correlation) matrix built from the two sets. This gives a closed-form, globally optimal solution (for the least squares criterion).

---

2. Algorithm / Approach

The general pattern:

• Step 1: Center each point set by subtracting its centroid.
• Step 2: Compute covariance between centered source and target points.
• Step 3: Use SVD of this covariance to extract the optimal rotation.
• Step 4: Recover translation from the centroids and the found rotation.
• Optionally, enforce proper rotation ($\det(R)=1$) to avoid reflections.

This is a one-shot, non-iterative algorithm: no gradient descent needed, just linear algebra (means + SVD).

---

3. Step-by-Step Strategy

Assume you have $N$ corresponding points in $d$ dimensions:

• Source: P \in \mathbb{R}^{N \times d}$(rows are$p_i^T$)
• Target: $Q \in \mathbb{R}^{N \times d}$ (rows are $q_i^T$)

1. Compute centroids

mu_P = P.mean(axis=0) # shape (d,)
mu_Q = Q.mean(axis=0) # shape (d,)

2. Center the point sets

P_centered = P - mu_P
Q_centered = Q - mu_Q

3. Compute covariance matrix

• Typically: H = P_{\text{centered}}^T Q_{\text{centered}} (shape $d \times d$).

H = P_centered.T @ Q_centered

4. SVD of covariance

• Compute $H = U S V^T$.

U, S, Vt = np.linalg.svd(H)

5. Compute rotation

• Naively: $R = V U^T$ (or $R = Vt.T @ U.T$ in code).
• Then optionally fix for improper rotation (reflection) if $\det(R) < 0$.

R = Vt.T @ U.T
if np.linalg.det(R) < 0:
Vt[-1, :] *= -1
R = Vt.T @ U.T

6. Compute translation

• Using centroids: $t = \mu_Q - R \mu_P$.

t = mu_Q - R @ mu_P

7. Transform source points (if needed)

P_aligned = (R @ P.T).T + t # shape (N, d)

---

4. Common Pitfalls

• Mismatched correspondence: The method assumes point $p_i$ corresponds to $q_i$. Wrong ordering = wrong alignment.
• Shape and orientation of matrices:
• Be consistent whether points are rows or columns.
• Make sure $H$ is $d \times d$, not $N \times N$.
• Reflections:
• SVD can produce a rotation with $\det(R) = -1$ (a reflection). For a pure rigid transform, enforce $\det(R)=1$ by flipping the sign of the last singular vector.
• Degenerate cases:
• If all points are collinear or identical, $H$ becomes rank-deficient and the rotation may not be unique.
• Numerical precision:
• For very large or very small coordinate values, centering and SVD can suffer from numeric issues; using double precision and standard libraries usually suffices.

---

5. Time & Space Complexity

Let $N$ be the number of points and $d$ the dimension (usually $d=2$ or $3$).

• Time complexity

• Computing centroids: $O(Nd)$

• Centering: $O(Nd)$

• Covariance matrix $H$: $O(Nd^2)$

• SVD on $d \times d$ matrix: $O(d^3)$

• Overall: $O(Nd^2)$, which is linear in $N$ for fixed $d$.

• Space complexity

• Storing $P, Q$: $O(Nd)$

• Centered copies and $H$: also $O(Nd)$ + $O(d^2)$

• Overall: $O(Nd)$ extra space (or in-place with some care).