Procrustes Problem for rank-deficient input 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).

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)
    }
  }
}


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.
