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.