# Minimize ||AX - Y|| for a matrix A that lies in a special orthogonal group

Let $$X$$ and $$Y$$ be two given $$k\times n$$ real matrices. If $$A$$ is a $$k\times k$$ real matrix then $$AX - Y$$ is a $$k\times n$$ real matrix. Applying the Frobenius norm $$\| AX - Y \|$$, we get a non-negative real number. Let $$\mathcal M_k$$ be the space of all $$k\times k$$ real matrices. We have a continuous map $$\mathcal M_k\to\mathbb R$$, given by $$A\mapsto \|AX-Y\|$$. If we restrict $$A$$ to lie in the special orthogonal group $$SO(k)\subset \mathcal M_k$$, we obtain a continuous map $$SO(k) \to \mathbb R$$. Since $$SO(k)$$ is compact, this must attain its non-negative minimum, for some $$A\in SO(k)$$.

What I am asking for is a numerical algorithm with the input of two matrices $$X$$ and $$Y$$ and the output of an arbitrarily close approximation to such an $$A$$ (It is easy to construct examples with an infinite number of arg-minima, but usually, the minimum will be unique).

• en.wikipedia.org/wiki/Orthogonal_Procrustes_problem#Generalized/… Commented Oct 5, 2022 at 13:20
• The objective function you consider, $\|AX-U\|_F^2$ is quadratic in the elements of $A$. The constraint you have is that you want to optimize subject to $AA^T=I$, which is also quadratic. So you have a "quadratically constrained, quadratic problem (QCQP)": en.wikipedia.org/wiki/… Commented Oct 5, 2022 at 15:21
• The additional constraint that $Xe=0$ ($e$ is a vector of all ones) makes this more complicated than a traditional Procrustes problem. The requirement that $\det A=1$ is also a complicating constraint. It might help if you gave a mathematical formulation of the problem that made these constraints clear and if you explained the background of the problem in more detail. Commented Oct 6, 2022 at 2:05
• @FedericoPoloni gives a crucial and very helpful link, to the Wikipedia article "Orthogonal Procrustes". This explains a solution, established in 1924, for O(k), instead of my SO(k). I have been learning more about this, and I believe (but I'm not yet certain) that the proof given there can be adapted to my case. I can report back when I understand the problem better, but I currently believe there will be a reasonably fast solution, based on singular decomposition Commented Oct 18, 2022 at 20:58
• @DavidEpstein Look at the specific section I linked in the Wikipedia page; there is a solution also for SO(k). Commented Oct 18, 2022 at 21:48

Just to give this a (link-only, unfortunately) answer: this is a variant of the orthogonal Procrustes problem. The solution restricted to SO(k) is detailed in a later section of the same page.

$$\def\o{{\tt1}}\def\p{\partial} \def\LR#1{\left(#1\right)} \def\op#1{\operatorname{#1}} \def\trace#1{\op{Tr}\LR{#1}} \def\qiq{\quad\implies\quad} \def\grad#1#2{\frac{\p #1}{\p #2}} \def\c#1{\color{red}{#1}} \def\CLR#1{\c{\LR{#1}}} \def\fracLR#1#2{\LR{\frac{#1}{#2}}}$$Use an unconstrainted matrix $$U$$ to construct $$A$$ \eqalign{ S &= \LR{U-U^T} &\qiq dS = \LR{dU-dU^T} \\ A &= \LR{I+S}^{-1}\LR{I-S} &\qiq dA= -\LR{I+S}^{-1}dS\LR{I+A} \\ } where a Cayley Transform has been used to map the skew matrix to $$SO(k)$$

For later convenience, define the matrices \eqalign{ M &= {\LR{Y-AX}X^T} \\ N &= \LR{I-S}^{-1}M\LR{I+A^T} \qquad\qquad \\ }

The Frobenius product is a concise notation for the trace \eqalign{ X:Y &= \sum_{i=1}^k\sum_{j=1}^n X_{ij}Y_{ij} \;=\; \trace{X^TY} \\ X:X &= \|X\|^2_F \\ } This is also called the double-dot or double contraction product.

The properties of the summation (or trace function) allow the terms in a Frobenius product to be rearranged in many different ways, e.g. \eqalign{ X:Y &= Y:X \\ X:Y &= X^T:Y^T \\ W:\LR{XY} &= \LR{WY^T}:X \\&= \LR{X^TW}:Y \qquad\qquad\qquad \\ } As with the Hadamard product, the matrix on each side of the multiplication symbol $$(:)$$ must have exactly the same dimensions.

Write the objective function using the above notation and calculate its gradient \eqalign{ \phi &= \tfrac 12\LR{Y-AX} : \LR{Y-AX} \\ d\phi &= \LR{Y-AX} : \LR{-dA\,X} \\ &= \LR{\LR{Y-AX}X^T} : \LR{\LR{I+S}^{-1}dS\LR{I+A}} \\ &= \LR{\LR{I-S}^{-1}M\LR{I+A^T}} : {dS} \\ &= N : \LR{dU-dU^T} \\ &= \LR{N-N^T} : dU \\ \grad{\phi}{U} &= \LR{N-N^T} \:=\: G \qquad\quad \big\{{\rm gradient\:matrix}\big\} \\ } Now you can use $$G$$ in your favorite gradient descent iteration, e.g. Barzilai-Borwein \eqalign{ G_k &= G(U_k) \\ U_{k+1} &= U_k - \left[ \frac{\LR{U_k-U_{k-1}}:\LR{G_k-G_{k-1}}}{\LR{G_k-G_{k-1}}:\LR{G_k-G_{k-1}}} \right] G_{k} \\ k &= k+1 \\ } After calculating the optimal $$U,\,$$ calculate the corresponding $$\,A=A(U)$$

• Interesting: I wonder if your process can be related to the proof given in the Wikipedia Orthogonal Procruster article which uses the singular value decomposition of YX^T. Commented Nov 6, 2022 at 22:09
• I don't know how people quantify speed for numerical algorithms, But are you able to say anything about the relative speed of the two approaches? Commented Nov 6, 2022 at 22:14