How do I deal with this Diag operator when differentiating with respect to a matrix? This question is an extension of this question. The arbitrary function $B(\cdot)$ is now as specified as follows:
$$
B(\mu, \sigma^{2}) = \exp \left(\mu + \frac{1}{2}\sigma^{2} \right).
$$
So if I write
$$
\begin{align*}
  s &= X_{i}\beta + Z_{i}\mu_{i} + \frac{1}{2}\operatorname{dg}\left(Z_{i}\Lambda_{i}Z_{i}^{T} \right)\\
  e &= \exp \left(s \right)\\
  E &= \operatorname{Diag}\left(e\right)
\end{align*}
$$
where $\exp$ is applied element wise and $\operatorname{Diag}$ is an operator that creates a diagonal matrix with the elements of the vector inside.
The problem is that I need to compute $\frac{\partial^{2} \underline{\ell}}{\partial \operatorname{vech}\left( \Lambda_{i}\right)\partial\mu_{i} }$, I'm stuck with calculating the differential within the $\operatorname{Diag}$ operator.
$$
\frac{\partial \underline{\ell} }{\partial \operatorname{vech}\left(\Lambda_{i} \right)} = \frac{1}{2}\left(\operatorname{vec} \left(\Lambda_{i}^{-1} - \Sigma^{-1} - Z_{i}^{T}EZ_{i} \right)\right)^{T}D_{K}.
$$
Setting $f= \frac{\partial \underline{\ell} }{\partial \operatorname{vech}\left(\Lambda_{i} \right)}$ and trying to differentiate this again with respect to $\mu_{i}$,
$$
\begin{align*}
  df &= -\frac{1}{2} \left(\operatorname{vec} \left(Z_{i}^{T} dE Z_{i} \right)\right)^{T}D_{K}\\
  dE &= \operatorname{Diag}\left(de \right)\\
 &= \operatorname{Diag}\left(e \circ ds \right)\\
&= \operatorname{Diag} \left(e \circ Z_{i}d\mu_{i} \right)
\end{align*}
$$
the differential is inside the $\operatorname{Diag}$ operator. How do I proceed from here?
 A: Ok, I did some research and got the answer for the question so I will post an answer for my own question. So we start with $df = -\frac{1}{2}\left(\operatorname{vec}\left(Z_{i}^{T}dEZ_{i} \right) \right)^{T}D_{K}$ and $dE = \operatorname{Diag}\left(de\right)$. We will use the following theorem:
$$
\operatorname{vec}\left(A^{T}\operatorname{Diag}\left(b\right)A \right) = \mathcal{Q}\left(A\right)^{T}b
$$
where $\mathcal{Q}\left(A\right) = \left(A \otimes \boldsymbol{1}^{T}\right) \circ \left(\boldsymbol{1}^{T}\otimes A\right)$. $\boldsymbol{1}$ is a $d \times 1$ vector of ones when $A$ is $n \times d$ matrix. ($\otimes$ is the Kronecker product and $\circ$ is the Hadamard product.) Then,
$$
\begin{align*}df &= -\frac{1}{2}\left(\operatorname{vec}\left(Z_{i}^{T}\operatorname{Diag}\left(de\right)Z_{i} \right) \right)^{T}D_{K} \\ &= -\frac{1}{2}\left(de\right)^{T}\, \mathcal{Q}\left(Z_{i}\right)D_{K}\\ &= -\frac{1}{2}\left(e\circ ds\right)^{T} \, \mathcal{Q}\left(Z_{i}\right)D_{K} \\ &= -\frac{1}{2}\left(\operatorname{Diag}\left(e\right) Z_{i} \, d\mu_{i} \right)^{T} \, \mathcal{Q}\left(Z_{i}\right)D_{K} \\ &= -\frac{1}{2}\left(d\mu_{i}\right)^{T}Z_{i}^{T}E\, \mathcal{Q}\left(Z_{i}\right)D_{K}  \end{align*} 
$$
Thus,
$$
\frac{\partial^{2} \ell}{\partial \operatorname{vech}\left(\Lambda_{i}\right) \partial \mu_{i}}=-\frac{1}{2}Z_{i}^{T}E\, \mathcal{Q}\left(Z_{i}\right)D_{K}.
$$
For those of you who want the proof of the theorem I used, please refer to M.P. Wand 2013, "Fully Simplified Multivariate Normal Updates in Non-Conjugate Variational Message Passing".
