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 from a multivariate 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 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$ and how is this moment computed? Furthermore, is my fist computational trick on track? Thanks.