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

  • 3
    $\begingroup$ en.wikipedia.org/wiki/Orthogonal_Procrustes_problem#Generalized/… $\endgroup$ Commented Oct 5, 2022 at 13:20
  • 1
    $\begingroup$ 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/… $\endgroup$ Commented Oct 5, 2022 at 15:21
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Oct 6, 2022 at 2:05
  • 2
    $\begingroup$ @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 $\endgroup$ Commented Oct 18, 2022 at 20:58
  • 2
    $\begingroup$ @DavidEpstein Look at the specific section I linked in the Wikipedia page; there is a solution also for SO(k). $\endgroup$ Commented Oct 18, 2022 at 21:48

2 Answers 2


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

  • $\begingroup$ 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. $\endgroup$ Commented Nov 6, 2022 at 22:09
  • $\begingroup$ 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? $\endgroup$ Commented Nov 6, 2022 at 22:14

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.