Let $X$ be a real $n\times n$ positive semidefinite matrix of rank $m\le n$ and let $Y\in\mathbb{R}^{m\times n}$ be the unique matrix satisfying (i) $X=Y^\top Y$, and (ii) $Y\, [I\, |\, 0]^\top = L$ with $L\in\mathbb{R}^{m\times m}$ being upper triangular with positive diagonal entries. (Notice that when $n=m$, $Y$ coincides with the (unique) Cholesky factor of $X$.) Let $A\in\mathbb{R}^{m\times m}$ be a positive definite matrix and $B\in\mathbb{R}^{n\times m}$ be such that $YB$ is nonsingular. Consider the following matrix-valued function
$$
f(X)=(YB)^{-1}A (YB)^{-\top}.
$$



> 
> I'm looking for a closed-form expression of the matrix differential $\mathrm{d}\, f(X)$.
> 
> 
> 




---


**My attempt.** A simple observation is that, in case $A=\alpha I$, $\alpha\in \mathbb{R}$, the sought differential reduces to
$$
\mathrm{d}\, f(X) = -\alpha(B^\top XB)^{-1} {B^{\top} \mathrm{d}\, X\, B} (B^\top XB)^{-1}.
$$


For the case of general positive definite $A$, I describe below my attempt that deals with vectorized ($\mathrm{vec}$ operation) differentials. First, using chain rule,
\begin{equation}\tag{1}\label{eq:1}
\mathrm{vec}(\mathrm{d}\, f(X)) = \frac{\mathrm{vec}(\mathrm{d}\, f(X))}{\mathrm{vec}(\mathrm{d}\, Y)}\frac{\mathrm{vec}(\mathrm{d}\, Y)}{\mathrm{vec}(\mathrm{d}\, X)} \mathrm{vec}(\mathrm{d}\, X).
\end{equation}
Concerning the first term, we have
$$
\frac{\mathrm{vec}(\mathrm{d}\, f(X))}{\mathrm{vec}(\mathrm{d}\, Y)}=-\left((YB)^{-1}A(YB)^{-\top}B^\top \otimes (YB)^{-1}\right)-\left((YB)^{-1}\otimes (YB)^{-1}A(YB)^{-\top}B^\top \right)K^{n,m},
$$
where $K^{n,m}$ is an $nm\times nm$ [commutation matrix](https://en.wikipedia.org/wiki/Commutation_matrix) and $\otimes$ denotes Kronecker product. For the second term, using the same argument of the accepted answer to [this MO question](https://mathoverflow.net/questions/150427/the-derivative-of-the-cholesky-factor?noredirect=1&lq=1), we have
$$
\frac{\mathrm{vec}(\mathrm{d}\, Y)}{\mathrm{vec}(\mathrm{d}\, X)} = \left((Y^\top\otimes I)K^{n,m}+(I\otimes Y^\top)\right)^{-L},
$$
where $(\cdot)^{-L}$ denotes left-inversion. Plugging the above-derived expressions into \eqref{eq:1}, yields a (vectorized) formula for the sought differential. 


At this point, provided that my calculations are correct, I wonder whether Eq. \eqref{eq:1} can be further simplified and written in matrix form (if possible). I conjecture (or, more frankly, I hope) that the final expression is similar to the scalar ($A=\alpha I$) case, i.e., something of the form
$$
\mathrm{d}\, f(X) = -(YB)^{-1}A (YB)^{-\top} {B^{\top} \mathrm{d}\, X\, B} (YB)^{-1}A (YB)^{-\top}.
$$




---


**Edit.** I believe that the hard part of my question regards the computation of the differential of $Y$ w.r.t. $X$. Of course, this can be accomplished via vectorization, as described above, and subsequent "matricization". Nevertheless, I wonder whether a more "genuine" matrix expression for $\mathrm{d}Y$ exists.




---


Any comment/suggestion is very appreciated.



Closed-form expression for differential of matrix function
It suffices to consider the case $n=2$, $m=1$. Namely, write $Y=[Y\_1|Y\_2]$ etc, then $L=Y\_1$ is upper triangular with positives on the diagonal, and 
\begin{align\*}
X &= Y^\top Y = \begin{pmatrix} Y\_1^\top Y\_1 & Y\_1^\top Y\_2 \\ Y\_2^\top Y\_1 & Y\_2^\top Y\_2\end{pmatrix}\,,
\\
dX &= \begin{pmatrix} (dY\_1)^\top. Y\_1 + Y\_1^\top.dY\_1 & (dY\_1)^\top. Y\_2 + Y\_1^\top.dY\_2 \\ 
(dY\_2)^\top. Y\_1 + Y\_2^\top. dY\_1 & (dY\_2)^\top. Y\_2 + Y\_2^\top.dY\_2
\end{pmatrix} 
\\
YB &= Y\_1B\_1 + Y\_2B\_2
\end{align\*}
Moreover,
\begin{align\*}
df(X) =& -(YB)^{-1}.dY.B.(YB)^{-\top} - (YB)^{-1}.A.(YB)^{-\top}.(dY.B)^\top.(YB)^{-\top}
\\=&
-(YB)^{-1}.dY.B.f(X) - f(X)((YB)^{-1}.dY.B)^\top
\\
dY.B =& dY\_1.B\_1 + dY\_2.B\_2 
\end{align\*}


Continued:
==========


Now you can compute $dY\_1$ from $dX\_{1,1}$ as in the invertible situation. Let us drop subindices and go to the invertible situation. 
$C\mapsto C^{-\top}$ is the Cartan involution on the reductive Lie group $GL^+(m)$. Consider the Iwasawa decomposition $GL^+(m) = SO(m).A.N$, $A$ the diagonal matrices with positive entries, and $N$ the upper unipotent matrices (which equals here the Gram-Schmidt orthonormalisation procedure with the coefficients arranged in $A.N$).
First note that $Y:S\_+(m)\to AN$ is a smooth map into a Lie group, so $dY$, better $TY: TS\_+(m)\to T(AN)$, and the right logarithmic derivative
$\delta Y:= TY.Y^{-1}:TS\_+(m)\to \mathfrak{an}$ is a Lie algebra valued 1-form, $\delta Y\in \Omega^1(S\_+(m);\mathfrak{an})$.
You can get back $Y$ from the 1-form $\delta Y$ by Cartan development. Namely, $\delta Y$ describes a flat principal connection on the trivial principal $AN$-bundle $S\_+(m)\times AN \to S\_+(m)$, and any horizontal leaf of it is a right translate of the mapping $Y$. 
Moreover $Z\mapsto Z^\top$ restricts to the mapping $\mathfrak{an}\to \mathfrak{an}^\*$ corresponding to the inner product
$\operatorname{Trace}(U^\top.V)$ on $\mathfrak{gl}(m)$.
We have
\begin{gather\*}
(dY)^\top.Y + Y^{\top}.dY = dX = dX^{\top}
\\
\delta Y + (\delta Y)^\top = dY.Y^{-1} + Y^{-\top}.(dY)^\top = Y^{-\top}.dX.Y^{-1} = (dX.Y^{-1})^{\top}.Y^{-1}
\end{gather\*}
Now, $\delta Y + (\delta Y)^\top$ allows in a simple way to compute $\delta Y$ (take the upper triangular part and 1/2 of the diagonal entries). 


