I was wondering if anybody knows how to solve: $$\mathbb{E}{\mathbf{z} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})}\left[ (\mathbf{x}{i} - \mathbf{z})(\mathbf{x}{j} - \mathbf{z})^\top \exp\left( - (\mathbf{x}{i} - \mathbf{z})^\top M (\mathbf{x}{i} - \mathbf{z}) -  (\mathbf{x}{j} - \mathbf{z})^\top M (\mathbf{x}{j} - \mathbf{z}) \right) \right].$$  Here $M$ is positive-semidefinite, and $\mathbf{x}i, \mathbf{x}j$ are constant vectors. I tried to solve this by just doing the integral $n$ times. However, you need to be able to solve $$ \int_{-\infty}^\infty (x^n \cdot\exp(-ax - bx^2) \, dx. $$ You can do this by completing the square and making a change of variables, so you get $$ \int_{-\infty}^\infty \left(z-\frac{a}{2b}\right)^n \cdot\exp\left(-bz^2+\frac{a^2}{4b}\right) \,dz $$ which in theory you could expand solve by using binomial expansion and substituting in the $i$th moment of normal distribution. However, doing these repeatedly where a and b are some functions of the other variables don't seem like a promising lead to a closed form solution. 

I'm asking on behalf of a friend, and this is his thought process. I want to do something with the exponential first if I were to work on the problem.