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} (z-\frac{a}{2b})^{n}\cdot\exp(-bz^{2}+\frac{a^{2}}{4b}) \,dz $$ which in theory you could expand solve by using binomial expansion and substituting in the ith 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.