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