I am working with two random matrices, $Z$ and $H$:

1. $Z$ is an $n \times K$ matrix with entries sampled i.i.d. from a Bernoulli distribution: $Z_{ij} \sim \mathrm{Bernoulli}(p)$.
2. $H$ is a $K \times K$ matrix with entries sampled i.i.d. from a Gaussian distribution.

I've defined a matrix product as:
$ \kappa = ZHZ^T $

Given that $Z$ and $H$ are independent, I'm interested in computing $E\kappa^2$.

Any suggestions on how to approach this problem would be greatly appreciated. Thank you!

**Update**: I came across [this paper][1] that seems somewhat related to my problem, especially since I'm using an Indian buffet process prior $P(Z)=IBP(\alpha)$. However, the approximate posterior distribution applied for computation of expectation is a Bernoulli with variable $p$. Here's the progress I've made so far, but I'm uncertain of its correctness. I'd love to get feedback on whether this is a structured way to compute the expectation: 
I wrote $\mathrm{vec}(\kappa^2)=((Z \otimes Z) \mathrm{vec}(H) \otimes (Z \otimes Z) \mathrm{vec}(H)) \mathrm{vec}(I)$. Then I have $ \mathrm{vec}(H) \otimes \mathrm{vec}(H) = \mathrm{vec}(H \otimes H)$ 
\begin{equation}
    \mathrm{vec}(\kappa^2) = (Z \otimes Z \otimes Z \otimes Z) \mathrm{vec}(H \otimes H)
\end{equation}
The fourth order moment is
\begin{equation}
    \begin{split}
        \mathbb{E}_{Z}\Big[Z \otimes Z \otimes Z \otimes Z\Big]&=\boldsymbol{p}\otimes\boldsymbol{p}\otimes\boldsymbol{p}\otimes\boldsymbol{p}+\mathfrak{S}_6\big[\mathrm{diag}(\boldsymbol{p}-\boldsymbol{p}^2)\otimes\boldsymbol{p}\otimes\boldsymbol{p}\big]\\
        &+\mathfrak{S}_3\big[\mathrm{diag}(\boldsymbol{p}-\boldsymbol{p}^2)\otimes\mathrm{diag}(\boldsymbol{p}-\boldsymbol{p}^2)\big]\\
&+\mathfrak{S}_4\big[\mathrm{diag}\big(\boldsymbol{p}-3\boldsymbol{p}^2+2\boldsymbol{p}^3\big)\otimes\boldsymbol{p}\big]\\
&+\mathrm{diag}\big(\boldsymbol{p}-7\boldsymbol{p}^2+12\boldsymbol{p}^3-6\boldsymbol{p}^4\big)
    \end{split}
\end{equation}
where $\mathfrak{S}_3\big[A\big]_{ikj}=A_{ijk}+A_{jki}+A_{kij}$ and $\mathfrak{S}_6\big[A\big]_{ikj}=A_{ijk}+A_{kij}+A_{jki}+A_{jik}+A_{kji}+A_{ikj}$ if $A=B\otimes C$. I could not find the definition of $\mathfrak{S}_4$ in the paper. Could anyone clarify the definition of $\mathfrak{S}_4$? Furthermore, is my computation above on track? Thanks. 


  [1]: https://arxiv.org/pdf/1704.00003.pdf