Is there a closed-form formula for the derivative of the orthogonal polar factor of a matrix?

Question

$\newcommand{\psym}{\text{Psym}_n}$ $\newcommand{\sym}{\text{sym}}$ $\newcommand{\Sym}{\operatorname{Sym}}$ $\newcommand{\Skew}{\operatorname{Skew}}$ $\newcommand{\SO}{\operatorname{SO}_n}$ $\renewcommand{\skew}{\operatorname{skew}}$ $\newcommand{\GLp}{\operatorname{GL}_n^+}$

Let $\psym$ be the space of real symmetric positive-definite $n \times n$ matrices, and $\GLp$ be the group of real $n \times n$ invertible matrices with positive determinant.

Let $O:\GLp \to \SO$ be the orthogonal polar factor map, defined by requiring $A= O(A)P$ for some symmetric positive-definite $P$. Note that $O(A)=A(\sqrt{A^TA})^{-1}$.

Is there a nice "closed-form algebraic formula" for the differential $dO_A$? If not, perhaps there is a formula for $\langle dO_A(B),C \rangle $? (similarly to what happens with the Levi-Civita connection, where we have an implicit characterization of $\nabla_XY$ in terms of the Koszul formula).

I am fine with using positive matrix square roots and inverses, but not with using integral formulas or vectorization operations like here or here. (I also don't want to use explicitly the singular values of $A$).

Here are some partial results:

1. Since $dO_{QA}(QB)=QdO_A(B)$ (and $dO_{AQ}(BQ)=dO_A(B)Q$), the question can be reduced to the case where $A \in \psym$. (The "dual" orthogonal case is easy: $dO_{Id}(B)=\skew(B)$, and for every $Q \in \SO$, $dO_{Q}(QB)=Q\skew B$).

2. For every $B \in \skew, dO_A(OBP)=OB$, i.e. $dO_A(X)=XP^{-1}$ if $O^TXP^{-1} \in \skew$: Set $\alpha(t)=Oe^{tB}P$. Then $O(\alpha(t))=Oe^{tB}$, so $dO_A(OBP)=OB$.

The problem with calculating $dO_A(OBP)$ when $B \in \sym$, is that $e^{tB}P$ does not need to be symmetric, even though $e^{tB},P$ are both positive-definite. If this was the case, then it would imply $dO_A(OBP)=O\skew B$, which is false in general (see below).

---

One could conjecture that perhaps $dO_A(OBP)=OB$ holds also for $B \in \sym$, or equivalently that $dO_A(X)=XP^{-1}$ for every $X \in M_n$. This is false since $dO_A$ is not injective, due to dimensional incompatibility.

Another possible conjecture is $dO_A(OBP)=O\skew B$. However, this is also false:

Indeed, for $A=P \in \psym$, this reduces to $dO_P(BP)=\skew B$. Suppose that $C:=BP \in \sym$. Then must have $dO_P(BP)=0$. Indeed, by differentiating $A=OP$ we obtain $$ \dot A=\dot O P+O\dot P, \tag{1}$$

which for $A=P,O=Id$ becomes $ \dot A=\dot O P+\dot P$. Note that $C \in \sym \Rightarrow dP_P(C)=C$ (since $P_{\psym}=Id_{\psym}$ and $C \in T_P{\psym}$), i.e. $\dot P=\dot A$, which implies $\dot O=0$.

Thus, we proved that $BP \in \sym$ implies $dO_P(BP)=0$. However, this is incompatible with $dO_P(BP)=\skew B$ in general: $BP \in \sym \iff BP=PB^T$.

For $P=\text{diag}(\sigma_1,\sigma_2), B=\begin{pmatrix} 1 & b \\\ c & 1 \end{pmatrix}$, this happens if and only if $\sigma_2b=\sigma_1c$, so if $\sigma_1 \neq \sigma_2$, then $B$ is not symmetric. Thus, $dO_P(BP)=0 \neq \skew B$.

---

Comment: Equation $(1)$ implies that one can equivalently focus upon the "dual" problem, of computing $dP_A$, instead of $dO_A$. Here is a previous attempt of mine to go in this direction.

Answer 1

Define the matrices $$\eqalign{ F &= A(A^TA)^{-1/2} \cr Q &= I-AA^+ = \,Q^T \cr V^T &= A^+ \cr }$$ where $A^+$ denotes the pseudoinverse.

Then $$\eqalign{ FF^T &= A(A^TA)^{-1}A^T \cr&= AA^+ \cr d(FF^T) &= d(AA^+) \cr dF\,F^T+F\,dF^T &= Q\,dA\,A^+ + V\,dA^T\,Q \cr }$$ At this point, you have stressed that you don't want to use vectorization, so let's try solving it using 4th order tensors.

There are 3 isotropic tensors $({\mathcal H},{\mathcal J},{\mathcal K})$ and their components can be expressed in terms of Kronecker symbols. $$\eqalign{ {\mathcal H}_{ijkl} &= \delta_{ik}\,\delta_{jl} \cr {\mathcal J}_{ijkl} &= \delta_{ij}\,\delta_{kl} \cr {\mathcal K}_{ijkl} &= \delta_{il}\,\delta_{jk} \cr }$$ We can use these to rearrange that last differential $$\eqalign{ ({\mathcal H}F + F{\mathcal K}):dF &= (Q{\mathcal H}V + V{\mathcal H}Q:{\mathcal K}):dA \cr }$$ Now we need to calculate a 4th order tensor ${\mathcal P}$ which is the inverse (under the double-contraction product) of the tensor on the LHS, i.e. $$\eqalign{ {\mathcal H}&= {\mathcal P}:\big({\mathcal H}F + F{\mathcal K}\big)\cr }$$ This allow us to isolate $dF$ and calculate the gradient as $$\eqalign{ dF &= \Big({\mathcal P}:(Q{\mathcal H}V + V{\mathcal H}Q:{\mathcal K})\Big):dA \cr \frac{\partial F}{\partial A} &= {\mathcal P}:(Q{\mathcal H}V + V{\mathcal H}Q:{\mathcal K}) \cr }$$ In the above, juxtaposition indicates a single-contraction product, while a colon represents double-contraction. Here are some examples in component form $$\eqalign{ {\mathcal A} &= {\mathcal B}:{\mathcal C} &\implies {\mathcal A}_{ijmn} &= \sum_{kl}{\mathcal B}_{ijkl}\,{\mathcal C}_{klmn} \cr {\alpha} &= {B}:{C} &\implies {\alpha} &= \sum_{ij}{B}_{ij}\,{C}_{ij} \cr {A} &= BC &\implies {A}_{ik} &= \sum_{j}{B}_{ij}\,{C}_{jk} \cr {\mathcal A} &= B\,{\mathcal C} &\implies {\mathcal A}_{iklm} &= \sum_{j}{B}_{ij}\,{\mathcal C}_{jklm} \cr }$$