# High-order moments for a Gaussian random variable

I am currently reviewing the book "State Estimation For Robotics" by Timothy D. Barfoot ( document uploaded by the author http://asrl.utias.utoronto.ca/~tdb/bib/barfoot_ser17.pdf).

The second chapter gives a review of probability concepts. I have studied probability subjects during my undergrad but I get a little lost because of the mixture of algebra and these types of topics for extension to more than one dimension.

At this moment I am in the section that defines Isserlis' theorem (Section 2.2.2, page 16 ). This theorem is used to calculate high order moments of a Gaussian random variable. $$$$\boldsymbol{x}=(x_1,x_2,...,x_{2M})\in \mathbb{R}$$$$ $$$$E[x_1x_2x_3...x_{2M}]=\sum\prod E[x_ix_j]$$$$ Example: $$$$E[x_ix_jx_kx_l]=E[x_ix_j]E[x_kx_l]+E[x_ix_k]E[x_jx_l]+E[x_ix_l]E[x_jx_k]$$$$

And after naming the theorem it gives some examples of calculation of moments of high order. $$$$\boldsymbol{x}\sim \mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma})\in \mathbb{R}^N$$$$

$$$$E[\boldsymbol{x}\boldsymbol{x}^\boldsymbol{T}\boldsymbol{x}\boldsymbol{x}^\boldsymbol{T}]=E\left[\left[x_ix_j\left(\sum_{k=1}^{N}x_k^2\right)\right]_{ij}\right]=\left[\sum_{k=1}^N E[x_ix_kx_j^2]\right]_{ij}=\left[\sum_{k=1}^N\left(E[x_ix_j]E[x_k^2]+2E[x_ix_k]E[x_kx_j]\right)\right]_{ij}=\left[E \left[x_ix_j\right]\right]_{ij}\sum_{k=1}^N E\left[x_k^2\right] +\left[\sum_{k=1}^N 2E[x_ix_k]E[x_kx_j]\right]_{ij}=\boldsymbol{\Sigma} tr(\boldsymbol{\Sigma})+2\boldsymbol{\Sigma}^2=\boldsymbol{\Sigma} (tr(\boldsymbol{\Sigma} )\boldsymbol{1} +2\boldsymbol{\Sigma})$$$$

Where $$\boldsymbol{1}$$ is the identity matrix.

$$$$\boldsymbol{x}=\begin{bmatrix}\boldsymbol{x_1}\\\boldsymbol{x_2}\end{bmatrix} \sim \mathcal{N}\left(\boldsymbol{0},\begin{bmatrix}\boldsymbol{\Sigma}_{11}&\boldsymbol{\Sigma}_{12}\\\boldsymbol{\Sigma}_{12}^T&\boldsymbol{\Sigma}_{22}\end{bmatrix}\right)$$$$

Where dim($$\boldsymbol{x_1})=N_1$$ and dim($$\boldsymbol{x_2})=N_2$$

$$$$E[\boldsymbol{x}\boldsymbol{x_1}^\boldsymbol{T}\boldsymbol{x}_1\boldsymbol{x}^\boldsymbol{T}]=E\left[\left[x_ix_j\left(\sum_{k=1}^{N_1}x_k^2\right)\right]_{ij}\right]=\left[\sum_{k=1}^{N_1} E[x_ix_kx_j^2]\right]_{ij}=\left[\sum_{k=1}^{N_1}\left(E[x_ix_j]E[x_k^2]+2E[x_ix_k]E[x_kx_j]\right)\right]_{ij}=\left[E \left[x_ix_j\right]\right]_{ij}\sum_{k=1}^{N_1} E\left[x_k^2\right] +\left[\sum_{k=1}^{N_1} 2E[x_ix_k]E[x_kx_j]\right]_{ij}=\boldsymbol{\Sigma} tr(\boldsymbol{\Sigma}_{11})+2\begin{bmatrix}\boldsymbol{\Sigma}_{11}^2&\boldsymbol{\Sigma}_{11}\boldsymbol{\Sigma}_{12} \\ \boldsymbol{\Sigma}_{12}^\boldsymbol{T}\boldsymbol{\Sigma}_{11}&\boldsymbol{\Sigma}_{12}^\boldsymbol{T}\boldsymbol{\Sigma}_{12}\end{bmatrix}=\boldsymbol{\Sigma} \left(tr(\boldsymbol{\Sigma}_{11} )\boldsymbol{1} +2\begin{bmatrix}\boldsymbol{\Sigma}_{11}&\boldsymbol{\Sigma}_{12} \\ \boldsymbol{0}&\boldsymbol{0}\end{bmatrix}\right)$$$$

The notation $$[\cdot]_{ij}$$ implies populating the matrix $$\boldsymbol{A}=[a_{ij}]$$

The first example I have more or less arrived at the result although I still don't understand things. I have assumed the vector $$\boldsymbol{x}$$ is: $$$$\boldsymbol{x}=\begin{bmatrix}x_1\\x_2\\\vdots\\x_N\end{bmatrix}$$$$ Then I computed $$\boldsymbol{xx^T}$$ $$$$\boldsymbol{xx^T}=\begin{bmatrix}x_1\\x_2\\\vdots\\x_N\end{bmatrix}\begin{bmatrix}x_1&x_2&...&x_N\end{bmatrix}=\begin{bmatrix}x_1x_1&x_1x_2&\cdots&x_1x_N\\x_2x_1&x_2x_2&\cdots&x_2x_N\\\vdots&\vdots&\ddots&\vdots\\x_Nx_1&x_Nx_2&\cdots&x_Nx_N\end{bmatrix}=[x_ix_j]_{ij}$$$$ As $$\boldsymbol{xx^Txx^T}=(\boldsymbol{xx^T})^2$$ then each component of the new final matrix will be: $$$$\begin{bmatrix}x_ix_1&x_ix_2&\cdots&x_ix_N\end{bmatrix}\cdot\begin{bmatrix}x_1x_j&x_2x_j&\cdots&x_Nx_j\end{bmatrix}=x_ix_j\sum_{k=1}^N x_k^2$$$$ By the properties of linearity of expectation and Isserlis' Theorem

$$$$E\left[\left[x_ix_j\left(\sum_{k=1}^{N}x_k^2\right)\right]_{ij}\right]=\left[\sum_{k=1}^N E[x_ix_kx_j^2]\right]_{ij}=\left[\sum_{k=1}^{N}\left(E[x_ix_j]E[x_k^2]+2E[x_ix_k]E[x_kx_j]\right)\right]_{ij}$$$$ From this moment on, I do not understand specifically how he gets to $$\boldsymbol{\Sigma}^2$$ and how did the identity matrix ($$\boldsymbol{1}$$) appeared.

In the case of the second example I have performed the analogous operations and I have problems to understand how to compute the relations that it is doing. I leave as far as I have arrived. $$$$\boldsymbol{x}=\begin{bmatrix}\boldsymbol{x_1}\\\boldsymbol{x_2}\end{bmatrix}=\begin{bmatrix}x_1^1\\x_2^1\\\vdots\\x_{N_1}^1\\x_1^2\\x_2^2\\\vdots\\x_{N_2}^2\end{bmatrix}$$$$ $$$$\boldsymbol{xx_1^T}=\begin{bmatrix}x_1^1\\x_2^1\\\vdots\\x_{N_1}^1\\x_1^2\\x_2^2\\\vdots\\x_{N_2}^2\end{bmatrix}\begin{bmatrix}x_1^1&x_2^1&\cdots&x_{N_1}^1\end{bmatrix}=\begin{bmatrix}x_1^1x_1^1&x_1^1x_2^1&\cdots&x_1^1x_{N_1}^1\\x_2^1x_1^1&x_2^1x_2^1&\cdots&x_2^1x_{N_1}\\\vdots&\vdots&\ddots&\vdots\\x_{N_1}^1x_1^1&x_{N_1}^1x_2^1&\cdots&x_{N_1}^1x_{N_1}^1\\x_1^2x_1^1&x_1^2x_2^1&\cdots&x_1^2x_{N_1}^2\\x_2^2x_1^1&x_2^2x_2^1&\cdots&x_2^2x_{N_1}\\\vdots&\vdots&\ddots&\vdots\\x_{N_2}^2x_1^1&x_{N_2}^2x_2^1&\cdots&x_{N_2}^2x_{N_1}^1\end{bmatrix}$$$$

I am asking this question to understand where these relationships come from but also to ask for advice on how to perform this type of conversion between matrices and "linear" writing. In some other occasion I have tried to read some paper or other books where I also get lost when they start to do operations with matrices but with summations and this kind of notation. In my head I have the feeling that these types of relationships can be understood without as much matrix computations as I have done. I hope someone can help me and give me some advice or recommend literature or bibliography to consult to improve in this type of deductions.

I have transcribed the formulas from the book if something is not correct please consult the book to answer and let me know the error.

For the $$\Sigma^2$$ term, note that since $$E[x_i] = 0$$ for every $$i$$, then $$\text{Cov}(x_i,x_j) = E[x_ix_j]-E[x_i] E[x_j] = E[x_ix_j]$$ for any $$i,j$$. Then $$2 \sum_{k=1}^N E[x_i x_k] E[x_k x_j] = 2 \sum_{k=1}^N \text{Cov}(x_i,x_k) \text{Cov}(x_k x_j)= 2 \sum_{k=1}^N \Sigma_{ik} \Sigma_{kj} = 2 [\Sigma^2]_{ij}.$$
Now, the expression is $$\Sigma \text{Tr}(\Sigma)+ 2 \Sigma^2 = \Sigma \text{Tr}(\Sigma)I + 2 \Sigma\Sigma = \Sigma (\text{Tr}(\Sigma)I + 2 \Sigma)$$ Note that multiplying things by the identity doesn't change anything, and if we don't have the identity then the term $$(\text{Tr}(\Sigma) + 2 \Sigma)$$ makes no sense because the trace is a scalar and you cannot add a scalar to a matrix.
For the second example, note that $$x x_1^T x_1 x^T = x (x_1^T x_1) x^T = (x_1^T x_1) x x^T$$ in other words, $$(x_1^T x_1)$$ is a scalar. We can write $$x_1^T x_1 = x_{11}^2 + x_{12}^2 + \dots x_{1N_{1}}^2 = \sum_{k=1}^{N_1} x_k^2.$$
What this means is that the $$(i,j)$$-th element of the matrix $$xx_1^T x_1 x^T$$ is given by $$x_i x_j \sum_{k=1}^{N_1} x_k^2,$$ so if you compute the expectation of this element you can easily deduce the expectation of the entire matrix.