Closed-form expression for differential of matrix function 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 and $\otimes$ denotes Kronecker product. For the second term, using the same argument of the accepted answer to this MO question, 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.
 A: 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). 
