How to optimize for the best fit nonuniform-scale-rotation to a given 3×3 matrix?

Question

Given a matrix $L\in \mathbb{R}^{3 \times 3}$, I'm looking for a method to find the closest (in a least squares sense) product of a non-uniform scaling matrix and a rotation matrix:

$$ \min_{s\in\mathbb{R}^3,\ R\in SO(3) } \left\| L - \text{diag}(s)\,R \right\|_F^2$$

So far, the best method that I've come up with is to alternate between optimizing $R$ and $s$ in a block coordinate descent fashion. Given initial guess of $s$ (e.g., $s=\mathbb{1}$), repeating:

1. Update $R$: holding $s$ constant, the optimal $R$ is found as the solution to the resulting orthogonal procrustes problem.
2. Update $s$: holding $R$ constant, the optimal $s$ is found as via a linear system solve (and actually just extracting the diagonal: $s = \text{diag}(LR^T)$).

until convergence.

Trying this out on gazillions of random $L$s, I find it works pretty well in most cases, though occasionally convergence is very slow. I've also noticed that the initial guess doesn't appear to affect the resulting minimizers. Could it be proveable that there's only a single (global) minimum (under some mild assumptions of $L$)?

Is there a closed form solution?

If not, is there a superior or best-practice method for optimizing this constrained minimization problem?

A student found "Orthogonal, but not Orthonormal, Procrustes Problems", Everson 1997 which includes the algorithm above ("tandem" in their language) and some other descent approaches (I can't see if they show one's better than the other).

Answer 1

OK, sorry for the delay: yesterday was hectic and then I got distracted by another question. Here is a small piece of theory promised.

First, let us notice that $\|L-SR\|_F=\|LR^*-S\|_F$ so, as you correctly noticed yourself, $S$ must be the diagonal of $LR^*$ or you can improve and also our task can be restated as maximizing the sum of squares of the diagonal entries of $LR^*$ (because the full sum of squares does not depend on $R$). Also, since we can throw the minus sign on any entries of $S$ and the corresponding row of $R$ any time without changing anything, we are free to use the full orthogonal group for $R$ instead of $SO_n(\mathbb R)$ initially requested. In this case the Proctustes step just maximizes $\rm{Tr}[L(SR)^*]=\rm{Tr}[(SL)R*]$ and results in the symmetric positive definite $T=SLR^*$ (by using one SVD). When your algorithm reaches the limit, it should have $\rm{diag T}\approx S$ (with required precision) or you can improve $S$ and the functional further by a noticeable amount.

What we will show now that if $\rm{diag T}=S$, indeed, then the sum of squares of $LR^*$ cannot be raised any further. To this end, notice that the condition that $T=SLR^*$ is symmetric positive definite and $\rm{diag(T)}=S$ can be restated as $LR^*=AS$ where $A$ is a symmetric positive definite matrix with all $1$s on the diagonal. Morally, $A=S^{-1}TS^{-1}$ though $S$ may have zeroes, in which case the corresponding row of $T$ is identically $0$ and so, by symmetry, is the corresponding column, so that block of zeroes can be treated separately as $I\cdot 0$ and in the remaining block the division is legitimate. Our aim is to show that for every orthigonal $U$, the sum $\sum_i (ASU)_{ii}^2\le \sum_i S_{ii}^2$.

Write $A=\sum_k x^{(k)}\otimes x^{(k)}$ with $x^{(k)}=[x^{(k)}_1,\dots,x^{(k)}_n]$. All $1$'s on the diagonal of $A$ mean that $$ \sum_k (x^{(k)}_i)^2=1\tag{$*$} $$ for all $i$. If $u_i$ are the columns of $u$, then $$ (ASU)_{ii}^2=\left[\sum_k x^{(k)}_i (x^{(k)}Su_i)\right]^2 \\ \le \left[\sum_k (x^{(k)}_i)^2\right]\left[\sum_k (x^{(k)}Su_i)^2\right]=\sum_k (x^{(k)}Su_i)^2 $$ by Cauchy-Schwarz. Summing over $i$ and using that $u_i$ form an orthonormal basis, wi get $$ \sum_i (ASU)_{ii}^2\le \sum_k\|x^{(k)}S\|^2=\sum_i S_{ii}^2 $$ as required (here we used $(*)$ once again).

Thus, our task can be reduced to finding a decomposition $$ L=ASR $$ where $A$ is positive definite with ones on the diagonal, $S$ is diagonal (say, with non-negative entries, because we can throw all minuses on $R$), and $R$ is orthogonal.

It is not immediately clear to me how to approach it using SVD or some other standard machinery, but one may notice that if such a decomposition holds (and we just showed that your algorithm yields one in the limit, though we haven't discussed the convergence speed yet; I still have to think here), then $$ AS^2A=LL^*\Longleftrightarrow SAS=\sqrt{SLL^*S} $$ in the sense of positive definite matrices. Since we don't care too much about the off-diagonal entries of $A$ but the correct $S$ leads to the full solution in one Procrustes step, we can just as well rewrite it as $$ S=\sqrt{\rm{diag}\sqrt{SLL^*S}} $$ which is exactly the problem that was posted on MO a few years ago with the question whether we can just use Picard iteration method with the initial guess $S_0=\sqrt{\rm{diag}LL^*}$ to solve it fairly efficiently (the $\sqrt{}$ on positive definite matrices is also one SVD, and you use it for Procrustes anyway, so a single iteration takes the same time here but the convergence speed may be much faster). The answer is still unknown but it is, probably, time to revisit that old question (I don't remember the link now but I'll add it in the comments when I find it). Still, as I said, you may run your theoretically guaranteed method in parallel with this unproved one and have the best of the two worlds (like it is normally done with sure but slow bisection and fast but unreliable Newton for finding a root on an interval).

I hope to be able to say more later but this should give you some food for both thinking and experimenting for a while now. Apologies for the delay again :-)

Answer 2

Here I'll present a few remarks and one reasonably efficient (IMHO) algorithm. Unlike my previous answer the conclusions here are experimental more than theoretical.

1. You are absolutely right that the convergence speed of your iterations is just as fast as that of the square root iterations. However, both are frustratingly slow on the matrix $$ \begin{bmatrix} 1&2&3\\ 1&2.001&3\\ 1&2&3.001 \end{bmatrix} $$ or any other nearly degenerate matrix. This forced me to think of some way to speed it up and for a while I was playing with various convergence acceleration techniques, but could find nothing that was both fast and stable, so I'll skip the description of that story.

2. As I said, the problem is equivalent to maximizing the sum of the squares of the diagonal elements of $LR^*$ where $R$ is orthogonal. $S$ can that be read as the diagonal of that maximizing matrix (up to entry signs). The key idea that finally occurred to me was that one can solve the 2D problem exactly in one simple step. Indeed, the sum of squares of the diagonal elements of $$ \begin{bmatrix} x&y\\ z&t \end{bmatrix} \begin{bmatrix} \cos a&-\sin a\\ \sin a&\cos a \end{bmatrix} $$ can be written as $$ \frac{x^2+y^2+z^2+t^2}2+\frac{x^2+t^2-y^2-z^2}2\cos(2a)+[xy-zt]\sin(2a) $$ so the maximum is achieved when $[\cos(2a), \sin(2a)]$ is just the unit vector in the direction $[\frac{x^2+t^2-y^2-z^2}2, xy-zt]$ (if this latter vector is $[0,0]$, then all $a$ give the same result so you can just choose your favorite one, otherwise the optimal angle can be found using a couple trigonometric functions).

That leads to the idea to make one cycle through three $2D$ coordinate planes optimizing at each step. Let's call the resulting mapping $L\mapsto F(L)$.

Now the natural thing to do would be just to iterate the assignment $L=F(L)$. The iterations indeed, converge reasonably fast (10-15 of them are enough to get to the situation when you have the sum of the squares of the diagonal elements stabilizing within the standard $10^{-15}$ machine precision (note that the matrix itself is not obliged to stabilize that precisely because at the stationary point any deviation by $\sqrt\delta$ from the optimizer results just in the deviation $\delta$ in the value of the functional, so, technically, the fluctuations of the size up to $10^{-8}$ in the matrix elements may continue to occur). Unfortunately, the limit is only guaranteed to be a maximizer with respect to each 2D coordinate rotation but not with respect to the full orthogonal group. It is clear for the matrix $$ \begin{bmatrix} 1& 0.9& -0.9\\ 0.9 & 1 & 0.9\\ -0.9 & 0.9 & 1 \end{bmatrix} $$ where no 2D coordinate rotation can improve anything but the 3D can. However, there are more interesting examples where you arrive at a stable limit (this one was unstable, really) that is far from the global maximum. The computer found, for example, $$ \begin{bmatrix} -0.01 & 0.37 & 0.27\\ -0.18 & 0.28 & -0.24\\ -0.2 & -0.14 & -0.45 \end{bmatrix} $$ If you only iterate $F$ here, you get the limiting value for the sum of squares only of $0.392442207640134$ while the truth is $0.446704013852864$.

To get unstuck, just apply one iteration of your SVD map now and then. Note that while I do not know how many times one can hit the wrong limit, we do know how to recognize that we are at the true $\max$: the matrix $(\rm{diag} L)L$ should be symmetric non-negative definite, so you can just check that property and shake only if needed. In practice, I prefer to shake 4 times with 6 iterations of $F$ in between. Out of a million random examples, none failed to yield a true maximum.

The pseudo-code is below:

real[][] F(real[][] L)
{
real[][] LL=copy(L);
for(int k=0;k<3;++k)
for(int j=k+1;j<3;++j)
{
real A=(LL[k][k]^2+LL[j][j]^2-LL[k][j]^2-LL[j][k]^2)/2, 
B=LL[k][k]*LL[k][j]-LL[j][j]*LL[j][k];
complex P=A+Bi; if(P==0) P=1; else P/=abs(P);
P=sqrt(P);
real[][] U=copy(IdentityMatrix);
U[k][k]=P.x; U[k][j]=-P.y;
U[j][j]=P.x; U[j][k]=P.y;
LL=LL*U;
}
return LL;
}



real[][] MAP(real[][] L)
{
real[][] d=diag(L); 
for(int k=0;k<n;++k)
{
struct T {real[][] U; real[][] D; real[][] V}
T f=SVD(d*L); // so d*L=U*D*V is the singular value decomposition
return L*transpose(f.V)*transpose(f.U);
}

int M=6, K=5;
for(int m=0;m<M;++m)
{
for(int k=0;k<K;++k) L=F(L);
L=MAP1(L);
}

I leave writing the accurate check of the non-negativity of a quadratic form to you. For non-degenerate $L$, it suffices to check that $(\rm{diag}L)L$ has positive determinant but degenerate matrices require more care. Or we can think a bit of why such a check is not really necessary. That gives you $S$. As to $R$, you can either trace it through the process or, better, just compute it once at the end where you know $S$ (just store the original copy of $L$ to use it at that point).

You may want to replace the for loop over m by the while loop until you reach the desired precision. I guess $K=5$ is a fair equivalent of one SVD step just because bringing a 3 by 3 symmetric matrix to the diagonal form by 2D coordinate rotations usually takes 4-5 cycles over coordinate pairs, but even $K=1$ looks better than two pure MAP1 steps.

Try this and let me know if that gives you improved performance. My experiments show that the combination that I used is 2-3 times more efficient than pure SVD iterations in any reasonable metric, but I don't know if it is sufficient for your purposes. :-)

Answer 3

I'm still not sure whether a closed form solution exists. But I recognized that most implementations of SVD will be iterative (the simplest being power iterations), so the block-coordinate descent in the original question can be understood as

while not converged
    update s ← diag(LRᵀ)
    while not converged
       improve svd [U,∑,V] = diag(s) L
    update R ← UVᵀ

The idea here (presented without proof of convergence but with convincing randomized testing) is to merge these while-loops and update $s$ during the iterative singular value decomposition of $\text{diag}(s) L$:

s ← √diag(LLᵀ), initial guess
V ← L, initial guess of eigen vectors of Lᵀ diag(s)² L
while not converged
   V ← Lᵀ diag(s)² L V, power iteration
   V ← gram-schmidt on columns so that Vᵀ V = I
   d ← diag(Vᵀ Lᵀ diag(s)² L V), eigen values
   z ← √d,  adjust so that U diag(z) V = diag(s) L
   U ← diag(s) L V diag(z)⁻¹
   R ← U Vᵀ, update R
   s ← diag(L Rᵀ), update s

There's no need to track $R$ if we can decide to stop based on $s$ alone, so it might be possible to avoid some multiplications in there.

My experiments indicate that this converges just as quickly with respect to (outer) iteration count as the block-descent. It's my understanding Jacobi method is superior to power iterations, so I wonder if the same idea could be applied there.