Second order moment of multivariate Gaussian (bishop p. 83) When computing the second order moment of the Multivariate Gaussian on p. 83 of Bishop's book, the following derivation is given:

It is not clear to me why the integral on the right-hand side of the middle line vanishes due to symmetry. Does someone have a more clear explanation for this? Thanks for your time! 
 A: There are two ways of seeing this. The first is that if you compute the integral in just $y_i$, you are effectively computing the mean of a gaussian with mean 0 and variance $\lambda_k$. So the answer is 0. 
More generally, you're integrating a function $f(y_1,..,y_n)y_iy_j$ $(i\neq j)$, where $f$ itself is even in each $y_k$, but $y_i$ alone is odd, so even times odd is odd, so the integral must vanish when you integrate in $y_i$.  So in the integral the only terms that survive are when $i=j$, giving you $f(y_1,...,y_n)y_i^2$.
A: Another approach involves the moment generating function of $X\sim N(\mu,\Sigma)$:
$$
M_x(t) = exp(\mu't +\frac{1}{2} t'\Sigma t),
$$
which has Hessian with respect to $t$
$$
\dfrac{\partial M_x(t) }{(\partial t)(\partial t)'} 
=
M_x(t)\{\Sigma + \mu\mu' + G(t)\},
$$
where $G(t)$ is directly proportional to $t$. The second moment is obtained from the above after setting $t=0$.
A: After writing the $\textbf{z}$'s as a combination of orthonormal basis vectors, we have the following form as in the book:
$$
\frac{1}{(2{\pi})^{D/2}} \frac{1}{|\Sigma|^{1/2}} \sum_{i} \sum_{j} u_i u_j^{T} \int exp(-\frac{1}{2} \sum\frac{y_{k}^2}{\lambda_{k}})y_{i}y_{j} d\textbf{y}
$$
Define
$$
f_{D}(y_1,...,y_D) := exp(-\frac{1}{2} \sum\frac{y_{k}^2}{\lambda_{k}})y_{i}y_{j}
$$
and for $i\leq j$, define $f_{i}$ analogously.
From the definition of integration of differential forms (See Baby Rudin Chapter 10), we can calculate the integral w.r.t. $\textbf{y}$ by "dropping" one variable at each step.
Define
$$
f_{D-1}(y_1,...,y_{D-1})=\int_{y_D} f_{D}(y_1,...,y_D) dy_{D}
$$
Now after $D$ steps, we will arrive at a number, which is the integral w.r.t. $\textbf{y}$.
Consider two cases:
$i=j$: Then we have,
$$
f_{i-1}(y_1,...,y_{i-1})= [\prod_{k=i+1}^{D} (\lambda_{k})^{1/2}] (\sqrt{2\pi})^{D-i} exp(-\frac{1}{2}C_{i-1}) \int exp(-\frac{y_i^2}{2\lambda_{i}})y_i^2 dy_{i}
$$
where $\frac{1}{2}C_{i-1}$ is the rest of the terms inside the exponential.
Note the last integral is the unnormalized second moment of a normal distribution with mean $0$ and variance $\lambda{i}$, so is equal to the normalizing constant times the variance.
So eventually we will get $u_i u_i^{T} (\sqrt{2\pi})^D \prod_{k=1}^{D} (\lambda_{k})^{1/2} \lambda_{i}$.
Now if wlog $i > j$:
$$
f_{i-1}(y_1,...,y_{i-1})= y_{j} [\prod_{k=i+1}^{D} (\lambda_{k})^{1/2}] (\sqrt{2\pi})^{D-i} exp(-\frac{1}{2}C_{i-1}) \int exp(-\frac{y_i^2}{2\lambda_{i}})y_i dy_{i}
$$
where the last integral is the expectation of a normal distribution with mean $0$, so the value of $f_{i-1}$ is $0$, and eventually the integral vanishes.
