Expected value of $xx^{T}$ for multidimensional Gaussian I need a bit of help understanding a step in the derivation of the expected value of $\bf{x x^{T}}$, that is, $E[\bf{x x^{T}}]$ with a Gaussian distribution.
By definition, using the D-dimensional Gaussian distribution, we have:
$$E[{\bf{x x^{T}}}]=\frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{|\Sigma|^{\frac{1}{2}}} \int \exp{\left\{ -\frac{1}{2}(\bf{x-\mu)^T\Sigma^{-1}(x-}\mu) \right\}}{\bf{x x^{T}}}d\bf{x}$$
We substitute $\bf{x-\mu}$ by $\bf{z}$:
$$E[{\bf{x x^{T}}}]=\frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{|\Sigma|^{\frac{1}{2}}} \int \exp{\left\{ -\frac{1}{2}{\bf{z}}^T\Sigma^{-1}{\bf{z}} \right\}}({\bf{z+\mu}})(\bf{z+\mu})^{T}d{\bf{z}}$$
From this equation, the part that is giving me trouble is the $\bf{zz^{T}}$ term:
$$\frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{|\Sigma|^{\frac{1}{2}}} \int \exp{\left\{ -\frac{1}{2}{\bf{z}}^T\Sigma^{-1}{\bf{z}} \right\}}{\bf{z}}\bf{z}^{T}d{\bf{z}}$$
The book says:

Again, we can make use of the eigenvector expansion of the covariance matrix given by (2.45), together with the completeness of the set of eigenvectors, to write ${\bf{z}} = \sum_{j=1}^{D} y_{j}{\bf{u_{j}}}$, where $y_{j} = \bf{u_{j}^{T}z}$

The equation (2.45) is the eigenvector equation:
$${\bf{\Sigma u_{i}}}=\lambda_{i}{\bf{u_{i}}}$$
I know that we can express $\bf{z}$ using the eigenvectors of $\bf{\Sigma}$ but the summation above doesn't make much sense to me. By the way, $y_{j}=\bf{u_{j}^{T}z}$ comes from the exponent in the Gaussian distribution and works as a change of coordinates (from $x$ to $y$), so at least I understand that part.
Here is the relevant section in the book:

And this is the page referenced in the quote above:

So, can anyone explain that part to me? 
Thanks in advance
 A: I think your central question is how to write $z$ as
$$z = \sum_{j=1}^D (u_j^t z) \, u_j = \sum_{j=1}^D \langle u_j, z \rangle u_j \quad (1)$$
This, you can see as follows. As you have a complete, orthonormal set of EVs $\{ u_j \}$ that span $\mathbb{R}^D$, you can write
$$z = \sum_{j=1}^D \alpha_j u_j \quad (2)$$
for some $\alpha_i$. Applying $u_i$ to (2), you get
$$\begin{align}
\langle u_i, z \rangle &= \sum_{j=1}^D\alpha_j \langle u_i, u_j \rangle \quad \text{using (bi)linearity of } \langle .,.\rangle \\
& = \alpha_i \cdot 1 + \sum_{j \neq i} \alpha_j \cdot 0 \quad \text{using the orthonormality of the EVs}\\
&= \alpha_i,
\end{align}$$
and if you plug this back into (2), you get (1).
Edit: At least I understood that you wanted that particular line explained. The exercises at the end would be mostly unnecessary to write down here. Except for the first (what he calls "Eigenvector equation"), they all are steps in showing the so-called "Eigendecomposition of a (symmetric) matrix (Spectral Theorem)" - symmetry coming in to allow for getting an orthonormal basis. You can write such matrices
$$A = Q \Lambda Q^t,$$
with $Q$ being the matrix having as its columns the EVs of $A$, and $\Lambda$ being a diagonal matrix with the eigenvalues on its diagonal. When you look at it closely, this is (2.49) in your book. It's a standard result in linear algebra; proofs are online, eg here. This leaves only (2.45), which I've seen but if you really need it, I would have to think about first (and hope I'll remember/can re-prove).
