The equation in the OP posting is a particular case of the general equation
$$\begin{align}AXB=C\tag{0}\label{zero}\end{align}$$
Where $A, B, C$ are known matrices, $X$ is unknown, and all compatible with respect to matrix multiplication. Here I present first a simple analysis of these types of equations and then focus on the OP's.
The following result summarizes what it is known about \eqref{zero}
Lemma: Given $A\in\operatorname{mat}_\mathbb{C}(m,n)$ and $B\in\operatorname{mat}_\mathbb{C}(p,q)$, and $C\in\operatorname{mat}_\mathbb{C}(m,q)$. The equation
$$\begin{align}AXB=C\end{align}\tag{1}\label{one}$$ has a solution iff
$$\begin{align}AA^+CB^+B=C\tag{2}\label{two}\end{align}$$
(here $A^+$ and $B^+$ are the Moore Penrose pseudoinverses of $A$ and $B$ respectively). In either case, all solutions to \eqref{one} are of the form
$$X=A^+CB^+ +Y-A^+AYBB^+$$
for arbitrary (yet dimension compatible) matrix $Y$.
Proof of Lemma:
Suppose (1) has a solution $X$. Then
$$AA^+CB^+B=AA^+(AXB)B^+B=(AA^+A)X(BB^+B)=AXB=C$$
Conversely, if (2) holds, then $X=A^+CB^+$ solves \eqref{one}.
For the last statement, notice that for any $Y\in\operatorname{mat}_{\mathbb{C}}(n, p)$,
$A(Y-A^+AYBB^+)B=0$.
In the context of the OP, the equation
\begin{align}
B^*XB=A\tag{3}\label{three}
\end{align}
has solution iff
$$B^*(B^*)^+AB^+B=A$$
It is easy to check that for any matrix $B$,
$(B^*)^+=(B^+)^*$. Also, using some properties of the orthogonal projections, it can be shown that
$$B^+=(B^*B)^+B^*$$
Hence, any matrix of the form
\begin{align}
X &= (B^*)^+AB^+ + Y- (B^*)^+B^*YBB^+\\
&=(BB^*)^+BA(B^*B)^+B^*+ Y-(B^*)^+B^*YBB^+,
\end{align}
where $Y$ is an arbitrary matrix compatible with the product, is a solution to \eqref{three}. The choice $Y=0_{n\times p}$ corresponds to the "apparent" solution that the OP posted.
A few comments about the MP-invrse: Recall that $A^+$ is the matrix that satisfies
- $AXA=A$
- $XAX=X$
- $(AX)^*=AX$
- $(XA)^*=XA$
where for any matrix, $M^*$ is the transpose complex conjugate of $M$.
$A^+$ exists and is unique. Furthermore, $A^+$ satisfies
$$A^+=(A^*A)^+A^*$$
This last property is related to the problem of least square errors: Given $A\in\operatorname{mat}_\mathbb{C}(m,n)$ and $b\in\mathbb{C}^m$, find $\beta \in \mathbb{C}^n$ such that
$$\beta=\operatorname{arg}\min\{\|x\|_2: \|Ax-b\|_2=\min\{\|Ay-b\|_2:y\in\mathbb{C}^n\}\}$$
This problem has solution $\beta=A^+b$
Best Answer
Your assumption that $f$ has no derivative anywhere is wrong.
In fact, the derivative exists and is $0$ almost everywhere.
Every matrix with full rank has a neighborhood of matrices that also have full rank, so in this neighborhood $f$ is constant and thus differentiable.
For matrices $X$ that do not have full rank, $f$ is not differentiable at $X$. (The pseudoinverse is not differentiable at such points either; it isn't even continuous).