# compute expected value of a expression with normal distribution

how compute expected value of the below expression $$E[S^TS]$$

in which: $$S = \sum_{i=1}^{\nu} \mathbb{u}_i \mathbb{u}_i^{\top}$$

and $\mathbb{u}_i \sim \mathcal{N}(0,V); V \in R^{p \times p}$ is covariance matrix

[I have used $p$ differently: as a running index]

First off, $S = S^T$ because of the definition of $S$.

$E[S^TS]=E[SS]=E[(\sum_i u_i u_i')(\sum_j u_ju_j')]$

$=E[\sum_{i,j} u_i u_i' u_j u_j']=\sum_{i,j} E[u_i u_i' u_j u_j']=E[(u_i \cdot u_j) u_i u_j']$

$\Rightarrow E[(u_i u_i' u_j u_j')_{pq}]=\sum_k E[u_{ik}u_{jk} u_{ip} u_{jq}]$

Each scalar here is the $(p,q)$th element of the result you are looking for.

So each scalar we need to calculate is similar to a 4th moment.

Now, let $i \neq j$. Then, by marginal probability:

$p(u_{ik},u_{jk}, u_{ip}, u_{jq})=\int p(u_{i1},\ldots,u_{in}, u_{j1},\ldots,u_{jn}) du_\text{others}$

(applying independence of sample $i$ from $j$)

$=\int p(u_{i1},\ldots,u_{in}) d(\text{not }u_{ik} \text{ or } u_{ip}) \int p(u_{j1},\ldots,u_{jn}) d(\text{not }u_{jk} \text{ or } u_{jq})=P_i(u_{ik},u_{ip})P_j(u_{jk},u_{jq})$

We can use the joint rule for the multivariate normal to evaluate each.

$P_i = \text{normpdf}([u_{ik},u_{ip}]|0,V([k,p],[k,p]))$

$P_j = \text{normpdf}([u_{jk}, u_{jq}]|0,V([k,q],[k,q]))$

I am using MATLAB notation here, so that $V([a,b],[a,b])=\begin{bmatrix}V_{aa} &V_{ab} \\V_{ba} &V_{bb} \end{bmatrix}$

So, now that we have the joint over the 4 variables, we can calculate each expectation and perform the sums:

$E[u_{ik}u_{ip}u_{jk}u_{jq}]=\int u_{ik}u_{ip}u_{jk}u_{jq} P_i(u_{ik},u_{ip}) P_j(u_{jk},u_{jq})du_{ik}du_{ip}du_{jk}du_{jq}$

This would involve the tracking of several cases such as $k=i$, $k=p$, $i=j$, etc. and would take some effort to workout. I suspect that things should simplify as you work them out. But in principle at least, you should be able to calculate the final expression.