A Purely Algebraic Justification of the Kabsch-Umeyama Algorithm

TL;DR

This paper offers a purely algebraic proof of the Kabsch-Umeyama algorithm for the orthogonal Procrustes problem, avoiding the use of Lagrange multipliers and clarifying its mathematical foundations.

Contribution

It provides a transparent, algebraic justification of the Kabsch-Umeyama algorithm, enhancing understanding of its mathematical basis without relying on optimization techniques.

Findings

Algebraic proof of the Kabsch-Umeyama algorithm

Reduction of rigid motion problem to orthogonal Procrustes problem

Clarification of the algorithm's mathematical foundations

Abstract

The constrained orthogonal Procrustes problem is the least-squares problem that calls for a rotation matrix that optimally aligns two matrices of the same order. Over past decades, the algorithm of choice for solving this problem has been the Kabsch-Umeyama algorithm, which is effectively no more than the computation of the singular value decomposition of a particular matrix. Its justification, as presented separately by Kabsch and Umeyama, is not totally algebraic since it is based on solving the minimization problem via Lagrange multipliers. In order to provide a more transparent alternative, it is the main purpose of this paper to present a purely algebraic justification of the algorithm through the exclusive use of simple concepts from linear algebra. For the sake of completeness, a proof is also included of the well known and widely used fact that the orientation-preserving rigid motion problem, i.e., the least-squares problem that calls for an orientation-preserving rigid motion that optimally aligns two corresponding sets of points in d-dimensional Euclidean space, reduces to the constrained orthogonal Procrustes problem.

Equations

$$\displaystyle\sum_{i=1}^{n}||Uq_{i}-p_{i}||^{2}=\sum_{i=1}^{n}(Uq_{i}-p_{i})^{T}(Uq_{i}-p_{i})={\rm tr}\bigl{(}(UQ-P)^{T}(UQ-P)\bigr{)}$$

$$=$$

$$\bar{q}=\frac{1}{n}\sum _{i=1}^n{q}_i\mathrm{and}\bar{p}=\frac{1}{n}\sum _{i=1}^n{p}_i.$$

$$\Delta (\varphi )=\sum _{i=1}^n\Vert \varphi (q_i)-p_i{\Vert }^2,$$

$$||\varphi (q_i)-p_i|{|}^2-||\tau (q_i)-p_i|{|}^2=(U{q}_i+t-p_i{)}^T(U{q}_i+t-p_i)$$

$$\Delta (\varphi )-\Delta (\tau )$$

$$\sum _{i=1}^n\Vert U{q}_i+\bar{p}-U\bar{q}-p_i{\Vert }^2=\sum _{i=1}^n\Vert (U(q_i-\bar{q})-(p_i-\bar{p}){\Vert }^2$$

Full text

MSC: 15A18, 15A42, 65H17, 65K99, 93B60

Keywords: constrained, Frobenius, least squares, orientation preserving, orthogonal, Procrustes, rigid motion, rotation, SVD, trace

1 Introduction

In the orthogonal Procrustes problem [2, 8], given real matrices $P$ and $Q$ of size $d\times n$, the problem is that of finding a $d\times d$ orthogonal matrix $U$ that minimizes $\|UQ-P\|_{F}$, where $\|\cdot\|_{F}$ denotes the Frobenius norm of a matrix. On the other hand, in the constrained orthogonal Procrustes problem [4, 5, 10], the same function is minimized but $U$ is constrained to be a rotation matrix, i.e., an orthogonal matrix of determinant 1. By letting $p_{i}$, $q_{i}$, $i=1\ldots,n$, be the vectors in $\mathbb{R}^{d}$ that are the columns from left to right of $P$ and $Q$, respectively, since clearly $\|UQ-P\|_{F}^{2}=\sum_{i=1}^{n}\|Uq_{i}-p_{i}\|^{2}$, where $\|\cdot\|$ denotes the $d-$dimensional Euclidean norm, then an alternative formulation of the two problems above is that of finding an orthogonal matrix $U$ (of determinant 1 for the constrained problem) that minimizes $\sum_{i=1}^{n}\|Uq_{i}-p_{i}\|^{2}$. We note that minimizing matrices do exist for the two problems as the function being minimized is continuous and both the set of orthogonal matrices and the set of rotation matrices are compact (in some topology). Finally, in the same vein, another problem of interest is the orientation-preserving rigid motion problem which is that of finding an orientation-preserving rigid motion $\phi$ of $\mathbb{R}^{d}$ that minimizes $\sum_{i=1}^{n}\|\phi(q_{i})-p_{i}\|^{2}$. An affine linear function $\phi$, $\phi:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, is a rigid motion of $\mathbb{R}^{d}$ if it is of the form $\phi(q)=Uq+t$ for $q\in\mathbb{R}^{d}$, where $U$ is a $d\times d$ orthogonal matrix, and $t$ is a vector in $\mathbb{R}^{d}$. The rigid motion $\phi$ is orientation preserving if det$(U)=1$, i.e., the determinant of $U$ equals 1. With $\bar{p}$, $\bar{q}$ denoting the centroids of $\{p_{i}\}$, $\{q_{i}\}$, respectively, as will be shown in Section 3 of this paper, this problem can be reduced to the constrained orthogonal Procrustes problem by translating $\{p_{i}\}$, $\{q_{i}\}$ to become $\{p_{i}-\bar{p}\}$, $\{q_{i}-\bar{q}\}$, respectively, so that the centroid of each set becomes $0\in\mathbb{R}^{d}$.

With $P$, $Q$, $p_{i}$, $q_{i}$, $i=1,\ldots,n$, as above, in this paper we focus our attention mostly on the constrained orthogonal Procrustes problem, and therefore wish to find a $d\times d$ rotation matrix $U$ that minimizes $\sum_{i=1}^{n}\|Uq_{i}-p_{i}\|^{2}$.

With this purpose in mind, we rewrite $\sum_{i=1}^{n}\|Uq_{i}-p_{i}\|^{2}$ as follows, where given a square matrix $R$, ${\rm tr}(R)$ stands for the trace of $R$.

[TABLE]

[TABLE]

Since only the third term in the last line above depends on $U$, it suffices to find a $d\times d$ rotation matrix $U$ that maximizes ${\rm tr}(P^{T}UQ)$. Since ${\rm tr}(P^{T}UQ)={\rm tr}(UQP^{T})$ (note in general ${\rm tr}(AB)={\rm tr}(BA)$, $A$ an $n\times d$ matrix, $B$ a $d\times n$ matrix), denoting the $d\times d$ matrix $QP^{T}$ by $M$, this problem is equivalent to finding a $d\times d$ rotation matrix $U$ that maximizes ${\rm tr}(UM)$, and it is well known that one such $U$ can be computed from the singular value decomposition of $M$ [4, 5, 10]. This is done with the Kabsch-Umeyama algorithm [4, 5, 10] (see Algorithm Kabsch-Umeyama below, where $\mathrm{diag}\{s_{1},\ldots,s_{d}\}$ is the $d\times d$ diagonal matrix with numbers $s_{1},\ldots,s_{d}$ as the elements of the diagonal, in that order running from the upper left to the lower right of the matrix). A singular value decomposition (SVD) [6] of $M$ is a representation of the form $M=VSW^{T}$, where $V$ and $W$ are $d\times d$ orthogonal matrices and $S$ is a $d\times d$ diagonal matrix with the singular values of $M$, which are nonnegative real numbers, appearing in the diagonal of $S$ in descending order, from the upper left to the lower right of $S$. Finally, note that any matrix, not necessarily square, has a singular value decomposition, not necessarily unique [6].

   

  Algorithm Kabsch-Umeyama

   

  Compute $d\times d$ matrix $M=QP^{T}$.

  Compute SVD of $M$, i.e., identify $d\times d$ matrices $V$, $S$, $W$,

  so that $M=VSW^{T}$ in the SVD sense.

  Set $s_{1}=\ldots=s_{d-1}=1$.

  If $\det(VW)>0$, then set $s_{d}=1$, else set $s_{d}=-1$.

  Set $\tilde{S}=\mathrm{diag}\{s_{1},\ldots,s_{d}\}$.

  Return $d\times d$ rotation matrix $U=W\tilde{S}V^{T}$.

Algorithm Kabsch-Umeyama has existed for several decades [4, 5, 10], however the known justifications of the algorithm [4, 5, 10] are not totally algebraic as they are based on exploiting the optimization technique of Lagrange multipliers. It is the main purpose of this paper to justify the algorithm in a purely algebraic manner through the exclusive use of simple concepts from linear algebra. This is done in Section 2 of the paper. Finally, we note that applications of the algorithm can be found, notably in the field of functional and shape data analysis [1, 9], where, in particular, the shapes of two curves are compared, in part by optimally rotating one curve to match the other.

2 Algebraic Justification of the Kabsch-Umeyama Algorithm

We justify Algorithm Kabsch-Umeyama using exclusively simple concepts from linear algebra, mostly in the proof of the following useful proposition. We note that most of the proof of the proposition is concerned with proving (3) of the proposition. Thus, it seems reasonable to say that any justification of the algorithm that requires the conclusion in (3) but lacks a proof for it, is not exactly complete. See page 47 of the otherwise excellent thesis in [7] for an example of this situation. See [3] for an outline of this dissertation.

Proposition 1: If $D=\mathrm{diag}\{\sigma_{1},\ldots,\sigma_{d}\}$, $\sigma_{j}\geq 0$, $j=1,\ldots,d$, and $W$ is a $d\times d$ orthogonal matrix, then

1. tr$(WD)\leq\sum_{j=1}^{d}\sigma_{j}$.

2. If $B$ is a $d\times d$ orthogonal matrix, $S=B^{T}DB$, then $\mathrm{tr}(WS)\leq\mathrm{tr}(S)$.

3. If det$(W)=-1$, $\sigma_{d}\leq\sigma_{j}$, $j=1,\ldots,d-1$, then $\mathrm{tr}(WD)\leq\sum_{j=1}^{d-1}\sigma_{j}-\sigma_{d}$.

Proof: Since $W$ is orthogonal and if $W_{kj}$, $k,j=1,\ldots,d$, are the entries of $W$, then, in particular, $W_{jj}\leq 1$, $j=1,\ldots,d$, so that

$\mathrm{tr}(WD)=\sum_{j=1}^{d}W_{jj}\sigma_{j}\leq\sum_{j=1}^{d}\sigma_{j}$, and therefore statement (1) holds.

Accordingly, assuming $B$ is a $d\times d$ orthogonal matrix, since $BWB^{T}$ is also orthogonal, it follows from (1) that

$\mathrm{tr}(WS)=\mathrm{tr}(WB^{T}DB)=\mathrm{tr}(BWB^{T}D)\leq\sum_{j=1}^{d}\sigma_{j}=\mathrm{tr}(D)=\mathrm{tr}(S)$, and therefore (2) holds.

If det$(W)=-1$, we show next that a $d\times d$ orthogonal matrix $B$ can be identified so that with $\bar{W}=B^{T}WB$, then $\bar{W}=\left(\begin{smallmatrix}W_{0}&O\\ O^{T}&-1\\ \end{smallmatrix}\right)$, $W_{0}$ interpreted as the upper leftmost $d-1\times d-1$ entries of $\bar{W}$ and as a $d-1\times d-1$ matrix as well; $O$ interpreted as a vertical column or vector of $d-1$ zeroes.

With $I$ as the $d\times d$ identity matrix, then det$(W)=-1$ implies $\mathrm{det}(W+I)=-\mathrm{det}(W)\mathrm{det}(W+I)=-\mathrm{det}(W^{T})\mathrm{det}(W+I)=-\mathrm{det}(I+W^{T})=-\mathrm{det}(I+W)$ which implies det$(W+I)=0$ so that $x\not=0$ exists in $\mathbb{R}^{d}$ with $Wx=-x$. It also follows then that $W^{T}Wx=W^{T}(-x)$ which gives $x=-W^{T}x$ so that $W^{T}x=-x$ as well.

Letting $b_{d}=x$, vectors $b_{1},\ldots,b_{d-1}$ can be obtained so that $b_{1},\ldots,b_{d}$ form a basis of $\mathbb{R}^{d}$, and by the Gram-Schmidt process starting with $b_{d}$, we may assume $b_{1},\ldots,b_{d}$ form an orthonormal basis of $\mathbb{R}^{d}$ with $Wb_{d}=W^{T}b_{d}=-b_{d}$. Letting $B=(b_{1},\ldots,b_{d})$, interpreted as a $d\times d$ matrix with columns $b_{1},\ldots,b_{d}$, in that order, it then follows that $B$ is orthogonal, and with $\bar{W}=B^{T}WB$ and $W_{0}$, $O$ as previously described, noting $B^{T}Wb_{d}=B^{T}(-b_{d})=\left(\begin{smallmatrix}O\\ -1\\ \end{smallmatrix}\right)$ and $b_{d}^{T}WB=(W^{T}b_{d})^{T}B=(-b_{d})^{T}B=(O^{T}\,-1)$, then $\bar{W}=\left(\begin{smallmatrix}W_{0}&O\\ O^{T}&-1\\ \end{smallmatrix}\right)$. Note $\bar{W}$ is orthogonal and therefore so is the $d-1\times d-1$ matrix $W_{0}$.

Let $S=B^{T}DB$ and write $S=\left(\begin{smallmatrix}S_{0}&a\\ b^{T}&\gamma\\ \end{smallmatrix}\right)$, $S_{0}$ interpreted as the upper leftmost $d-1\times d-1$ entries of $S$ and as a $d-1\times d-1$ matrix as well; $a$ and $b$ interpreted as vertical columns or vectors of $d-1$ entries, and $\gamma$ as a scalar.

Note $\mathrm{tr}(WD)=\mathrm{tr}(B^{T}WDB)=\mathrm{tr}(B^{T}WBB^{T}DB)=\mathrm{tr}(\bar{W}S)$, so that $\bar{W}S=$ $\left(\begin{smallmatrix}W_{0}&O\\ O^{T}&-1\\ \end{smallmatrix}\right)$ $\left(\begin{smallmatrix}S_{0}&a\\ b^{T}&\gamma\\ \end{smallmatrix}\right)=$ $\left(\begin{smallmatrix}W_{0}S_{0}&W_{0}a\\ -b^{T}&-\gamma\\ \end{smallmatrix}\right)$ gives $\mathrm{tr}(WD)=\mathrm{tr}(W_{0}S_{0})-\gamma$.

We show $\mathrm{tr}(W_{0}S_{0})\leq\mathrm{tr}(S_{0})$. For this purpose let $\hat{W}=\left(\begin{smallmatrix}W_{0}&O\\ O^{T}&1\\ \end{smallmatrix}\right)$, $W_{0}$ and $O$ as above. Since $W_{0}$ is orthogonal, then clearly $\hat{W}$ is a $d\times d$ orthogonal matrix, and by (2), $\mathrm{tr}(\hat{W}S)\leq\mathrm{tr}(S)$ so that $\hat{W}S=$ $\left(\begin{smallmatrix}W_{0}&O\\ O^{T}&1\\ \end{smallmatrix}\right)$ $\left(\begin{smallmatrix}S_{0}&a\\ b^{T}&\gamma\\ \end{smallmatrix}\right)=$ $\left(\begin{smallmatrix}W_{0}S_{0}&W_{0}a\\ b^{T}&\gamma\\ \end{smallmatrix}\right)$ gives $\mathrm{tr}(W_{0}S_{0})+\gamma=\mathrm{tr}(\hat{W}S)\leq\mathrm{tr}(S)=\mathrm{tr}(S_{0})+\gamma$. Thus, $\mathrm{tr}(W_{0}S_{0})\leq\mathrm{tr}(S_{0})$.

Note $\mathrm{tr}(S_{0})+\gamma=\mathrm{tr}(S)=\mathrm{tr}(D)$, and if $B_{kj}$, $k,j=1,\ldots,d$ are the entries of $B$, then $\gamma=\sum_{k=1}^{d}B_{kd}^{2}\sigma_{k}$, a convex combination of the $\sigma_{k}$’s, so that $\gamma\geq\sigma_{d}$. It then follows that

$\mathrm{tr}(WD)=\mathrm{tr}(W_{0}S_{0})-\gamma\leq\mathrm{tr}(S_{0})-\gamma=\mathrm{tr}(D)-\gamma-\gamma\leq\sum_{j=1}^{d-1}\sigma_{j}-\sigma_{d}$, and therefore (3) holds. $\square$

Finally, the following theorem, a consequence of Proposition 1, justifies the Kabsch-Umeyama algorithm.

Theorem 1: Given a $d\times d$ matrix $M$, let $V$, $S$, $W$ be $d\times d$ matrices such that the singular value decomposition of $M$ gives $M=VSW^{T}$. If det$(VW)>0$, then $U=WV^{T}$ maximizes tr$(UM)$ over all $d\times d$ rotation matrices $U$. Otherwise, if det$(VW)<0$, with $\tilde{S}=\mathrm{diag}\{s_{1},\ldots,s_{d}\}$, $s_{1}=\ldots=s_{d-1}=1,\ s_{d}=-1$, then $U=W\tilde{S}V^{T}$ maximizes tr$(UM)$ over all $d\times d$ rotation matrices $U$.

Proof: Let $\sigma_{j}$, $j=1,\ldots,d$, $\sigma_{1}\geq\sigma_{2}\geq\ldots\geq\sigma_{d}\geq 0$, be the singular values of $M$, so that $S=\mathrm{diag}\{\sigma_{1},\ldots,\sigma_{d}\}$.

Assume det$(VW)>0$. If $U$ is any rotation matrix, then $U$ is orthogonal. From (1) of Proposition 1 since $W^{T}UV$ is orthogonal, then

$\mathrm{tr}(UM)=\mathrm{tr}(UVSW^{T})=\mathrm{tr}(W^{T}UVS)\leq\sum_{j=1}^{d}\sigma_{j}.$

On the other hand, if $U=WV^{T}$, then $U$ is clearly orthogonal, det$(U)=1$, and $\mathrm{tr}(UM)=\mathrm{tr}(WV^{T}VSW^{T})=\mathrm{tr}(WSW^{T})=\mathrm{tr}(S)=\sum_{j=1}^{d}\sigma_{j}.$

Thus, $U=WV^{T}$ maximizes tr$(UM)$ over all $d\times d$ rotation matrices $U$.

Finally, assume det$(VW)<0$. If $U$ is any rotation matrix, then $U$ is orthogonal and det$(U)=1$. From (3) of Proposition 1 since $W^{T}UV$ is orthogonal and det$(W^{T}UV)=-1$, then

$\mathrm{tr}(UM)=\mathrm{tr}(UVSW^{T})=\mathrm{tr}(W^{T}UVS)\leq\sum_{j=1}^{d-1}\sigma_{j}-\sigma_{d}.$

On the other hand, if $U=W\tilde{S}V^{T}$, then $U$ is clearly orthogonal, det$(U)=1$, and $\mathrm{tr}(UM)=\mathrm{tr}(W\tilde{S}V^{T}VSW^{T})=\mathrm{tr}(W\tilde{S}SW^{T})=\mathrm{tr}(\tilde{S}S)=\sum_{j=1}^{d-1}\sigma_{j}-\sigma_{d}.$

Thus, $U=W\tilde{S}V^{T}$ maximizes tr$(UM)$ over all $d\times d$ rotation matrices $U$. $\square$

3 Reduction of the Orientation-Preserving Rigid Motion Problem to the

Constrained Orthogonal Procrustes Problem

Although not exactly related to the main goal of this paper, for the sake of completeness, we show the orientation-preserving rigid motion problem reduces to the constrained orthogonal Procrustes problem. For this purpose, let $\bar{q}$ and $\bar{p}$ denote the centroids of the sets $\{q_{i}\}_{i=1}^{n}$ and $\{p_{i}\}_{i=1}^{n}$ in $\mathbb{R}^{d}$, respectively:

[TABLE]

First we prove a proposition that shows, in particular, that if $\hat{\phi}(\bar{q})\not=\bar{p}$, then $\phi=\hat{\phi}$ does not minimize

[TABLE]

the minimization occurring over either the set of all rigid motions $\phi$ of $\mathbb{R}^{d}$ or the smaller set of rigid motions $\phi$ of $\mathbb{R}^{d}$ that are orientation preserving.

Proposition 2: Let $\phi$ be a rigid motion of $\mathbb{R}^{d}$ with $\phi(\bar{q})\not=\bar{p}$ and define an affine linear function $\tau$, $\tau:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, by $\tau(q)=\phi(q)-\phi(\bar{q})+\bar{p}$ for $q\in\mathbb{R}^{d}$. Then $\tau$ is a rigid motion of $\mathbb{R}^{d}$, $\tau(\bar{q})=\bar{p}$, $\Delta(\tau)<\Delta(\phi)$, and if $\phi$ is orientation preserving, then so is $\tau$.

Proof: Clearly $\tau(\bar{q})=\bar{p}$. Let $U$ be a $d\times d$ orthogonal matrix and $t\in\mathbb{R}^{d}$ such that $\phi(q)=Uq+t$ for $q\in\mathbb{R}^{d}$. Then $\tau(q)=Uq-U\bar{q}+\bar{p}$ so that $\tau$ is a rigid motion of $\mathbb{R}^{d}$, $\tau$ is orientation preserving if $\phi$ is, and for $1\leq i\leq n$, we have

[TABLE]

It then follows that

[TABLE]

as $\phi(\bar{q})-\bar{p}$ is nonzero. Thus $\Delta(\tau)<\Delta(\phi)$. $\square$

Finally, the following corollary, a consequence of Proposition 2, shows that the problem of finding an orientation-preserving rigid motion $\phi$ of $\mathbb{R}^{d}$ that minimizes $\sum_{i=1}^{n}\|\phi(q_{i})-p_{i}\|^{2}$ can be reduced to a constrained orthogonal Procrustes problem which, of course, then can be solved with the Kabsch-Umeyama algorithm. Here $r_{i}=p_{i}-\bar{p}$, $s_{i}=q_{i}-\bar{q}$, for $i=1,\ldots,n$, and if $\bar{r}=\frac{1}{n}\sum_{i=1}^{n}r_{i}$, $\bar{s}=\frac{1}{n}\sum_{i=1}^{n}s_{i}$, then clearly $\bar{r}=\bar{s}=0$.

Corollary 1: Let $\hat{U}$ be such that $U=\hat{U}$ minimizes $\sum_{i=1}^{n}\|Us_{i}-r_{i}\|^{2}$ over all $d\times d$ rotation matrices $U$. Let $\hat{t}=\bar{p}-\hat{U}\bar{q}$, and let $\hat{\phi}$ be given by $\hat{\phi}(q)=\hat{U}q+\hat{t}$ for $q\in\mathbb{R}^{d}$. Then $\phi=\hat{\phi}$ minimizes $\sum_{i=1}^{n}\|\phi(q_{i})-p_{i}\|^{2}$ over all orientation-preserving rigid motions $\phi$ of $\mathbb{R}^{d}$.

Proof: One such $\hat{U}$ can be computed with the Kabsch-Umeyama algorithm.

By Proposition 2, in order to minimize $\sum_{i=1}^{n}\|\phi(q_{i})-p_{i}\|^{2}$ over all orientation-preserving rigid motions $\phi$ of $\mathbb{R}^{d}$, it suffices to do it over those for which $\phi(\bar{q})=\bar{p}$. Therefore, it suffices to minimize $\sum_{i=1}^{n}\|Uq_{i}+t-p_{i}\|^{2}$ with $t=\bar{p}-U\bar{q}$ over all $d\times d$ rotation matrices $U$, i.e., it suffices to minimize

[TABLE]

over all $d\times d$ rotation matrices $U$. But minimizing the last expression is equivalent to minimizing $\sum_{i=1}^{n}\|Us_{i}-r_{i}\|^{2}$ over all $d\times d$ rotation matrices $U$. Since $U=\hat{U}$ is a solution to this last problem, it then follows that $U=\hat{U}$ minimizes $\sum_{i=1}^{n}\|Uq_{i}+\bar{p}-U\bar{q}-p_{i}\|^{2}$ $=\sum_{i=1}^{n}\|Uq_{i}+t-p_{i}\|^{2}$ with $t=\bar{p}-U\bar{q}$ over all $d\times d$ rotation matrices $U$. Consequently, if $\hat{t}=\bar{p}-\hat{U}\bar{q}$, and $\hat{\phi}$ is given by $\hat{\phi}(q)=\hat{U}q+\hat{t}$ for $q\in\mathbb{R}^{d}$, then $\phi=\hat{\phi}$ clearly minimizes $\sum_{i=1}^{n}\|\phi(q_{i})-p_{i}\|^{2}$ over all orientation-preserving rigid motions $\phi$ of $\mathbb{R}^{d}$. $\square$