# expectation of a quadratic function of a matrix variate normal distribution

I want to compute the following expectation term:

$$E[{\bf{XA}}{{\bf{X}}^T}]$$

where $${\bf X} \in R^{M \times M}$$ and its elements are normal random variables such that

$$vec\left( {\bf{X}} \right)\sim \cal N\left( {\boldsymbol \mu ,\bf \Sigma } \right)$$

$$\bf A$$ is a positive definite matrix with proper dimensions and $$vec(.)$$ is the vectorization operator. Any hint on how I can derive a nice formula?

$$\newcommand{\X}{\mathbf X} \newcommand{\si}{\sigma}$$ Suppose that $$A:=(a_{ij})_{i,j=1}^m$$ is an $$m\times m$$ matrix and $$\X=(X_{ij})_{i,j=1}^m$$ is a random $$m\times m$$ matrix with $$EX_{ij}=\mu_{ij}$$ and $$Cov(X_{ij},X_{kl})=\si_{ij,kl}$$. Then the $$il$$-entry of the matrix $$E\X A\X^T$$ is
$$(E\X A\X^T)_{il}=\sum_{j,k}EX_{ij}a_{jk}X_{lk} =\sum_{j,k}a_{jk}(\mu_{ij}\mu_{kl}+\si_{ij,lk}) =(MAM^T)_{il}+\sum_{j,k}a_{jk}\si_{ij,lk},$$ where $$M:=(\mu_{ij})_{i,j=1}^m$$. So, $$E\X A\X^T=MAM^T+R,$$ where $$R:=(r_{il})_{i,l=1}^m$$ with $$r_{il}:=\sum_{j,k}a_{jk}\si_{ij,lk}$$.