Suppose you have two n × n matrices, A and B, and you would like to find a rotation matrix Ω that lines up B with A. That is, you’d like to find Ω such that
A = ΩB.
This is asking too much, except in the trivial case of A and B being 1 × 1 matrices. You could view the matrix equation above as a set of n² equations in real numbers, but the space of orthogonal matrices only has n(n − 1) degrees of freedom [1].
When an equation doesn’t have an exact solution, the next best thing is to get as close as possible to a solution, typically in a least squares sense. The orthogonal Procrustes problem is to find an orthogonal matrix Ω minimizing the distance between A and ΩB That is, we want to minimize
|| A − ΩB ||
subject to the constraint that Ω is orthogonal. The matrix norm used in this problem is the Frobenius norm, the sum of the squares of the matrix entries. The Frobenius norm is the 2-norm if we straighten the matrices into vectors of dimension n².
Peter Schönemann found a solution to the orthogonal Procrustes problem in 1964. His solution involves singular value decomposition (SVD). This shouldn’t be surprising since SVD solves the problem of finding the closest thing to an inverse of an non-invertible matrix. (More on that here.)
Schönemann’s solution is to set M = ABT and find its singular value decomposition
M = UΣVT.
Then
Ω = UVT.
Python code
The following code illustrates solving the orthogonal Procrustes problem for random matrices.
import numpy as np
n = 3
# Generate random n x n matrices A and B
rng = np.random.default_rng(seed=20260211)
A = rng.standard_normal((n, n))
B = rng.standard_normal((n, n))
# Compute M = A * B^T
M = A @ B.T
# SVD: M = U * Sigma * V^T
U, s, Vt = np.linalg.svd(M, full_matrices=False)
# R = U * V^T
R = U @ Vt
# Verify that R * R^T is very nearly the identity matrix
print("||R^T R - I||_F =", np.linalg.norm(R.T @ R - np.eye(n), ord="fro"))
In this example the Frobenius norm between RRT and I is 4 × 10−16, so essentially RRT = I to machine precision.
Related posts
[1] Every column of an orthogonal matrix Ω must have length 1, so that gives n constraints. Furthermore, each pair of columns must be orthogonal, which gives n choose 2 more constraints. We start with Ω having n² degrees of freedom, but then remove n and then n(n − 1)/2 degrees of freedom.
n² − n − n(n − 1)/2 = n(n − 1)/2