6
$\begingroup$

It is well known that the solution to the orthogonal Procrustes problem $$ \textrm{arg min}_{\Omega \in \text{SO}(n)} ||Y - \Omega X||^2_2, $$ can be expressed in terms of an SVD of the covariance matrix $M = X Y^T = U\Sigma V^T$ and, ignoring $\det(M) = -1$, is $$ \Omega = V U^T. $$ However, the solution is non-unique if $\text{rank}(M) < n - 1$

I am trying to figure out how to get the minimal rotation in the case where $\text{rank}(M) < n-1$. This is important for me because interpolating the rotation produces an odd (i.e., non-direct) path (grey) between $X$ (red) and $Y$ (black).

Visualization of the problem

Formally (I think), I can express this as a regularized Procrustes problem $$ \textrm{arg min}_\Omega ||Y - \Omega X||^2_2 + \varepsilon ||\text{Log}(\Omega)||^2, $$ where $\varepsilon > 0$ and $\text{Log}$ is the matrix logarithm.

I have searched this website and the literature, but I can only find solutions that assume that $M$ is full rank:

  • Popular answer on the derivation of the solution to Procrustes problem
  • Golub and van Loan, Matrix Computations (3rd edition, page 601)
  • Dryden and Mardia, Statistical Shape Analysis (2nd edition, page 70, lemma 4.2) mentions that the solution might not be unique, but does not offer any solution

I think I could try to turn this into an optimization problem on the Rotation manifold and try to solve it using gradient descent. However, a direct solution based on the SVD would be much more appealing.


Below are two examples coded in R. The SVD approach succeeds for the first example but fails for the second.

# Rendering helper
library(rgl)
make_3d_plot <- function(X, Y, rot){
  clear3d()
  space <- 8
  decorate3d(xlim = c(-space, space), ylim  = c(-space, space), zlim  = c(-space, space))
  spheres3d(t(X), radius = 0.4, col = "black")
  spheres3d(t(Y), radius = 0.4, col = "red")
  
  rot_vec <- expm::logm(rot)
  
  for(t in seq(0, 1, length.out = 30)){
    spheres3d(t(expm::expm(t * rot_vec) %*% X), radius = 0.1, col = "grey")
  }
}

# SVD produces the optimal rotation
X <- rbind(0:4, 0, 0)
Y <- rbind(0, 0:4, 0)
crossprod <- X %*% t(Y)

svd <- svd(crossprod)
rot <- svd$v %*% t(svd$u)
make_3d_plot(X, Y, rot)
# SVD chooses a suboptimal rotation
set.seed(1)
proj <- qr.Q(qr(matrix(rnorm(3 * 3), nrow = 3, ncol = 3)))
X <- proj %*% rbind(0:4, 0, 0)
Y <- proj %*% rbind(0, 0:4, 0)
crossprod <- X %*% t(Y)

svd <- svd(crossprod)
rot <- svd$v %*% t(svd$u)
if(Matrix::det(rot) == -1){
  rot <- svd$v %*% diag(c(1, 1, -1)) %*% t(svd$u)  
}
make_3d_plot(X, Y, rot)

R function based on Mike Hawk's answer. It works exactly as advertised :)

procrustes <- function(Y, X){
  stopifnot(nrow(Y) == nrow(X))
  stopifnot(ncol(Y) == ncol(X))
  
  if(nrow(Y) == 0){
    matrix(nrow = 0, ncol = 0) 
  }else if(nrow(Y) == 1){
    matrix(1)
  }else if(nrow(Y) == 2){
    a <- sum(X * Y)
    b <- sum(Y[1,] * X[2,] - Y[2,] * X[1,])
    theta <- atan2(a, b) - pi/2
    matrix(c(cos(theta), sin(theta), -sin(theta), cos(theta)), ncol = 2)
  }else{
    M <- X %*% t(Y)
    svd <- svd(M)
    sel <- abs(svd$d) > 1e-12
if(any(! sel)){
  # Cross product is not full rank, recurse on subproblem
  U1 <- svd$u[,sel,drop=FALSE]
  V1 <- svd$v[,sel,drop=FALSE]
  U2 <- svd$u[,!sel,drop=FALSE]
  V2 <- svd$v[,!sel,drop=FALSE]
  fact <- Matrix::det(svd$v %*% t(svd$u))
  U2[,ncol(U2)] <- U2[,ncol(U2)] * fact
  sub_rot <- procrustes(t(V2), t(U2))
  Ucompl <- cbind(U1, U2 %*% t(sub_rot))
  Vcompl <- cbind(V1, V2)
  Vcompl %*% t(Ucompl)
}else{
  U <- svd$u
  V <- svd$v
      fact <- Matrix::det(V %*% t(U))
      U[,ncol(U)] <- U[,ncol(U)] * fact
      V %*% t(U)
    }
  }
}
$\endgroup$
3
  • $\begingroup$ Also check ICP literature. How rigorous does the formulation need to be? Your figure doesn't look like it is small angle, but if you can do a small angle approximation in Rodrigues' formula it will make things a lot easier. Aside from that, you can probably form a dual problem. In general, for ill-defined problems you must have either a prior or a regularization. In the past, everyone was all about problems and formulations needing to be convex, but with plug-and-play priors and machine learning that are not convex and producing great results that constraint has been loosened a lot. $\endgroup$ Dec 12, 2022 at 16:04
  • $\begingroup$ not sure if it matter yet, but do we always expect $\textrm{rank}(\mathbf{Y})=\textrm{rank}(\mathbf{X})$? $\endgroup$ Dec 13, 2022 at 14:51
  • 1
    $\begingroup$ @JohnMadden, no, not necessarily. I have usually focused on $rank(M)$ because that can again be different from both $rank(Y)$ and $rank(X)$. $\endgroup$
    – const-ae
    Dec 13, 2022 at 15:42

1 Answer 1

5
+50
$\begingroup$

As you probably know, the locus of solutions will be all orthogonal matrices such that $\Omega(u_i)=v_i$, where $U\Sigma V^T$ is the svd of $XY^t$ ($\Sigma$ is $r\times r$, where $r$ is the rank of $XY^t$).

Denote $u_1, \dots, u_r$ the columns of $U$ and similarly for $V$. Define the subspace $H_U=span\{u_1, \dots, u_r\}$ and let $H_U^{\perp}$ be its orthogonal complement. Finally, let $u_{r+1},\dots, u_{n}$ be any orthogonal basis of $H_U^{\perp}$, and similarly for $v_{r+1},\dots, v_n$.

It is clear that any $\Omega$ must satisfy $\Omega(H_U^{\perp})=H_V^{\perp}$. So there exist coefficients $w_{ij}$ defined by $\Omega(u_{r+i})=\sum_j w_{ij}v_{r+j}$. It is straightforward to show that $\{w_{ij}\}_{ij}$ are a $(n-r)\times (n-r)$ orthgonal matrix.

We will define the size of an orthogonal matrix $\Omega$ as $|\Omega-I|_F$ This is a reasonable measure of the "size" of an orthogonal matrix, because it is minimized for the identity matrix. It is also considerably easier to work with than the matrix logarithm.

So we'll formalize the problem as finding the orthogonal matrix that minimizes $|\Omega-I|_F$ subject to the constraints $\Omega(u_i)=v_i, i\leq r$. It's easy to see that this is equivalent to maximizing $Tr(\Omega)$.

Since the trace does not depend on the choice of basis, we have: $$Tr(\Omega)=\sum_i <\Omega(u_i),u_i>=\sum_{i\leq r} <u_i,v_i>+\sum_{i\leq n-r} <w_{ij}v_{r+j},u_{r+i}>$$

Now, the first sum is constant (does not depend on $w$) and the second, with some rearrangement, is equal to $Tr(w^tG)$ where $G_{ij}=<u_{r+i},v_{r+j}>$. Thus finding the minimal $\Omega$ amounts to solving another procrustes problem:

$$argmax_{w\in SO(n-r)} Tr(w^tG)=argmin_{w\in SO(n-r)}|V^{\perp}-U^{\perp}w|$$

This problem is also potentially rank-deficient, but at least we have strictly reduced the dimensionality from n to n-r. So we can implement a solution by induction: Given a rank-deficient problem, we reduce it as above from an optimization over $SO(n)$ to an optimization over the strictly smaller space $SO(n-r)$. If the reduced problem is full-rank, then we know how to solve it, and to transform the solution to the solution of the original problem. Otherwise, if the reduced problem is rank-deficient, then we apply the same procedure to the reduced problem. Eventually, we will either arrive at a full-rank problem or to a problem where $n=2$. So the last step is to solve the rank-deficient problem in the case $n=2$, which I leave as an exercise.

$\endgroup$
1
  • $\begingroup$ Awesome. This works great. Very neat approach :) $\endgroup$
    – const-ae
    Dec 14, 2022 at 16:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.