Proof verification for differentiation through complex SVD I have been working on this particular problem which comes up in research and it would be great if I could have the following proof verified.

Define a zero-mean complex-valued matrix $\mathbf{X}\in\mathbb{C}^{D\times K}$ and another matrix $\mathbf{\Phi}(\epsilon)\in\mathbb{C}^{D\times K}$, the latter depending on a vector $\epsilon \in \mathbb{R}^{D\times1}$  and the covariance matrix
  $$\mathbf{C}=\frac{1}{K}(\mathbf{X}\circ\mathbf{\Phi})(\mathbf{X}\circ\mathbf{\Phi})^H$$
  which is Hermitian. $\circ$ is the Hadamard product operator. $^H$ is the Hermitian (conjugate) transpose operator, $^\ast$ is the complex conjugate, and $^T$ is the usual transpose.
Let
  $$\mathbf{U}\mathbf{\Sigma}\mathbf{V}^H=\mathbf{C}$$
  be the SVD of $\mathbf{C}$. Further, since $\mathbf{C}$ is Hermitian, $\mathbf{U}=\mathbf{V}$, thus
  $$\mathbf{U}\mathbf{\Sigma}\mathbf{U}^H=\mathbf{C}$$
We now seek to find the derivative of 
  $$\mathbf{W}=\mathbf{U}\mathbf{\Sigma}^{-1/2}\mathbf{U}^H$$
  with respect to $\epsilon$.


We begin with
\begin{align}
\renewcommand{\d}{\partial}
\d\mathbf{C}
&= \d\mathbf{U}\mathbf{\Sigma}\mathbf{U}^H
    + \mathbf{U}\d\mathbf{\Sigma}\mathbf{U}^H
    + \mathbf{U}\mathbf{\Sigma}\d\mathbf{U}^H\\
&= \mathbf{U}\d\mathbf{\Sigma}\mathbf{U}^H
    + \d\mathbf{U}\mathbf{\Sigma}\mathbf{U}^H
    + \left(\d\mathbf{U}\mathbf{\Sigma}\mathbf{U}^H\right)^H
\end{align}
Letting $\d\mathbf{P} = U^H\d\mathbf{U}$ thus $\d\mathbf{P}^H = \d\mathbf{U}^HU$,
\begin{align}
\d\mathbf{C}
&= \mathbf{U}\d\mathbf{\Sigma}\mathbf{U}^H
    + \mathbf{U}\d\mathbf{P}\mathbf{\Sigma}\mathbf{U}^H
    + \mathbf{U}\mathbf{\Sigma}\d\mathbf{P}^H\mathbf{U}^H\\
\mathbf{U}^H\d\mathbf{C}\mathbf{U}
&= \d\mathbf{P}\mathbf{\Sigma}
    + \d\mathbf{\Sigma}
    + \mathbf{\Sigma}\d\mathbf{P}^H
\end{align}
Morever, we have
\begin{align}
    \mathbf{U}\mathbf{U}^H =\mathbf{U}^H\mathbf{U} &= \mathbf{I}
\end{align}
thus
\begin{align}
    \mathbf{U}^H\d\mathbf{U} + \d\mathbf{U}^H\mathbf{U} &= \mathbf{0}_D\\
    \d\mathbf{P} + \d\mathbf{P}^H &= \mathbf{0}_D
\end{align}
that is, $\d\mathbf{P}$ is anti-Hermitian and have purely imaginary diagonals. 
From earlier, we define auxiliary variables
\begin{align}
\d\mathbf{G}
    &=\mathbf{U}^H\d\mathbf{C}\mathbf{U}\\
\d\mathbf{G}^H
    &=\mathbf{U}^H\d\mathbf{C}\mathbf{U}\\
%\mathbf{H}\circ\d\mathbf{G} + \Im\left(\mathbf{I}\circ\d\mathbf{G}\right)\\
\d\mathbf{H}
    &= \d\mathbf{G}-\d\mathbf{\Sigma}\\
    &= \d\mathbf{P}\mathbf{\Sigma} + \mathbf{\Sigma}\d\mathbf{P}^H\\
    &= \d\mathbf{P}\mathbf{\Sigma} - \mathbf{\Sigma}\d\mathbf{P}\\
    &= \d\mathbf{H}^H
\end{align}
Hence  $\d\mathbf{G}$ and $\d\mathbf{H}$ are also Hermitian. Hence, it follows that the diagonal entries of $\d\mathbf{G}$ are real and the diagonal entries of $\d\mathbf{H}$ are zero.
Now
\begin{align}
\d\mathbf{H}\mathbf{\Sigma}
    &= \d\mathbf{P}\mathbf{\Sigma}^2 -\mathbf{\Sigma}\d\mathbf{P}\mathbf{\Sigma}\\
\mathbf{\Sigma}\d\mathbf{H}
    &= \mathbf{\Sigma}\d\mathbf{P}\mathbf{\Sigma} -\mathbf{\Sigma}^2\d\mathbf{P}\
\end{align}
Considering element-wise
\begin{align}
\d h_{ij}\sigma_j + \sigma_i\d h_{ij}
    &= \d{p}_{ij}\sigma_j^2 - \sigma_i^2\d{p}_{ij}\\
(\sigma_i + \sigma_j)\d h_{ij}
    &= (\sigma_j^2 - \sigma_i^2)\d{p}_{ij}\\
\d{p}_{ij}
    &= \dfrac{\sigma_i + \sigma_j}{\sigma_j^2 - \sigma_i^2}\d h_{ij}\\
    &= \dfrac{\d h_{ij}}{\sigma_j - \sigma_i}
\end{align}
for $i\ne j$, leaving the diagonal of $\d\mathbf{P}$ unknown.
Then I have
\begin{align}
\mathbf{W} 
    &= \mathbf{U}\mathbf{\Sigma}^{-1/2}\mathbf{U}^H\\
\d\mathbf{W} 
    &= \d\mathbf{U}\mathbf{\Sigma}^{-1/2}\mathbf{U}^H
        + \mathbf{U}\d\mathbf{\Sigma}^{-1/2}\mathbf{U}^H
        + \mathbf{U}\mathbf{\Sigma}^{-1/2}\d\mathbf{U}^H\\
    &= \mathbf{U}\d\mathbf{P}\mathbf{\Sigma}^{-1/2}\mathbf{U}^H
        + \mathbf{U}\d\mathbf{\Sigma}^{-1/2}\mathbf{U}^H
        + \mathbf{U}\mathbf{\Sigma}^{-1/2}\d\mathbf{P}^H\mathbf{U}^H\\
\mathbf{U}^H\d\mathbf{W}\mathbf{U}
    &= \d\mathbf{P}\mathbf{\Sigma}^{-1/2}
        + \d\mathbf{\Sigma}^{-1/2}
        + \mathbf{\Sigma}^{-1/2}\d\mathbf{P}^H\\
\d\mathbf{W}
    &= \mathbf{U}\left[\d\mathbf{P}\mathbf{\Sigma}^{-1/2}
        + (\d\mathbf{P}\mathbf{\Sigma}^{-1/2})^H
        + \d\mathbf{\Sigma}^{-1/2}\right]\mathbf{U}^H\\
    &= \mathbf{U}\left[\d\mathbf{P}\mathbf{\Sigma}^{-1/2}
        - \mathbf{\Sigma}^{-1/2}\d\mathbf{P}
        - \dfrac{1}{2}\mathbf{\Sigma}^{-3/2}(\d\mathbf{\Sigma})\right]\mathbf{U}^H
\end{align}
Note also that
\begin{align}
\left[\d\mathbf{P}\mathbf{\Sigma}^{-1/2}     
    -\mathbf{\Sigma}^{-1/2}\d\mathbf{P}\right]_{ij}
    &= \d{p}_{ij}\sigma_j^{-1/2} - \sigma_i^{-1/2}\d{p}_{ij}
\end{align}
On the diagonal
\begin{align}
    \left[\d\mathbf{P}\mathbf{\Sigma}^{-1/2}     
    -\mathbf{\Sigma}^{-1/2}\d\mathbf{P}\right]_{ii}
    &= \d{p}_{ii}\sigma_j^{-1/2} - \sigma_i^{-1/2}\d{p}_{ii}\\
    &= 0
\end{align}
and elsewhere
\begin{align}
\left[\d\mathbf{P}\mathbf{\Sigma}^{-1/2}     
    -\mathbf{\Sigma}^{-1/2}\d\mathbf{P}\right]_{ij}
    &= \d{p}_{ij}\sigma_j^{-1/2} - \sigma_i^{-1/2}\d{p}_{ij}\\
    &= \dfrac{\sigma_j^{-1/2}- \sigma_i^{-1/2}}{\sigma_ j- \sigma_i}\d h_{ij}
\end{align}
Hence,
\begin{align}
\d\mathbf{P}\mathbf{\Sigma}^{-1/2}     
    -\mathbf{\Sigma}^{-1/2}\d\mathbf{P}
    &= \mathbf{M}\circ\d\mathbf{H}\\
    &= \mathbf{M}\circ(\d\mathbf{G} - \d\mathbf{\Sigma})\\
    &= \mathbf{M}\circ\d\mathbf{G}\\
    &= \mathbf{M}\circ(\mathbf{U}^H\d\mathbf{C}\mathbf{U})
\end{align}
where
\begin{align}
m_{ij}
    &= \begin{cases}
    \dfrac{\sigma_j^{-1/2}- \sigma_i^{-1/2}}{\sigma_j - \sigma_i}  & i\ne j\\
    0   & i = j
    \end{cases}
\end{align}
In all, 
\begin{align}
\d\mathbf{W}
    &= \mathbf{U}\left[\mathbf{M}\circ\d\mathbf{G}
        -\frac{1}{2}\mathbf{\Sigma}^{-3/2}\d\mathbf{\Sigma}\right]\mathbf{U}^H
\end{align}
 A: Define some new variables for later convenience.
$$\eqalign{
&\lambda^2 = K,\quad B = X\circ\Phi \\
&x = {\rm vec}(X),\quad \phi={\rm vec}(\Phi),\quad G={\rm Diag}(x) \\
&b={\rm vec}(B)=x\circ\phi = G\phi,\quad J=\frac{\partial \phi}{\partial \varepsilon} \\
&C = \lambda^{-2}BB^H = U\Sigma U^H \\
}$$
Any analytic function of $C$ can be written in terms of the SVD and vice versa 
$$\eqalign{
f(C) &= U\,f(\Sigma)\,U^H \\
W &= U\bigg(\frac{1}{\sqrt{\Sigma}}\bigg)U^H = C^{-1/2}\\
W^2 &= C^{-1} = \lambda^2(BB^H)^{-1}
}$$
Calculate the differential and vectorize it (note that $C$ and $W$ are both hermitian).
$$\eqalign{
W\,dW + dW\,W &= -\lambda^2(BB^H)^{-1}\Big(B\,dB^H+dB\,B^H\Big)(BB^H)^{-1} \\
-\lambda^{-2}BB^H\big(W\,dW + dW\,W\big)BB^H &= B\,dB^H+dB\,B^H \\
-\lambda^{2}\big(CW\,dW\,C + C\,dW\,WC\big) &= B\,dB^H+dB\,B^H \\
-K\big(C^T\otimes CW  + C^TW^T\otimes C\big)\,dw
  &= (I\otimes B)L\,db^* + (B^*\otimes I)\,db \\
-K\big(C^*\otimes CW  + C^*W^*\otimes C\big)\,dw
  &= (I\otimes B)LG^*\,d\phi^* + (B^*\otimes I)G\,d\phi \\
}$$
where $L$ is the commutator matrix associated with vectorization, i.e.
$$
{\rm vec}(X^T)=L\;{\rm vec}(X) \\
{\rm vec}(X^H)=L\;{\rm vec}(X^*)$$
Compact the differential by introducing names for the Kronecker matrices.
Then find the gradient with respect to the $\varepsilon$-vector. 
$$\eqalign{
P\,dw &= Q\,d\phi + R\,d\phi^* \\
dw &= M\,d\phi + N\,d\phi^* \\
\frac{\partial w}{\partial \varepsilon} &= MJ + NJ^* \\
}$$
