expectation of the product of Gaussian kernels and their input

Question

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.

Answer 1

$\newcommand\R{\mathbb R}$Letting $I:=\mathbf I$, $a:=\mathbf xi$, $b:=\mathbf xj$, $N:=(4M+I)^{1/2}$, $c_1:=a^\top Ma+b^\top Mb$, $m:=2N^{-1}M(a+b)$, $c_2:=m^\top m/2=2(a+b)^\top M(4M+I)^{-1}M(a+b)$, completing the squares, using the substitutions $z=Ny$ and $y-m=v$, and doing a bit of algebra, we see that the expectation in question is $$\frac{e^{c_2-c_1}}{\det(4M+I)^{1/2}}(N^{-2}+N^{-1}mm^\top N^{-1}+ab^\top-am^\top N^{-1}-N^{-1}mb^\top). \tag{1}\label{1}$$ Here one can get rid of $N$ and $m$ by noting that $N^{-2}=(4M+I)^{-1}$, so that $N^{-1}m=2(4M+I)^{-1}M(a+b)$ and hence $m^\top N^{-1}=2(a+b)^\top M(4M+I)^{-1}$.

---

Details: The expectation in question is $$J:=(2\pi)^{-n/2}\int_{\R^n}dz\,(z-a)(z-b)^\top e^{g(z)}, \tag{2}\label{2}$$ where \begin{align} g(z)&:=-(z-a)^\top M(z-a)-(z-b)^\top M(z-b)-z^\top z/2 \\ &=-z^\top(2M+I/2)z+2(a+b)^\top Mz-c_1 \\ &=-c_1-y^\top y/2+m^\top y \\ &=-c_1+c_2-v^\top v/2. \tag{3}\label{3} \end{align} Also, \begin{align} (z-a)(z-b)^\top&=(N^{-1}(v+m)-a)(N^{-1}(v+m)-b)^\top \\ &=h(v):=N^{-1}(v+m)(v+m)^\top N^{-1}+ab^\top \\ &\qquad\qquad-a(v+m)^\top N^{-1}-N^{-1}(v+m)b^\top \\ &=N^{-1}vv^\top N^{-1}+N^{-1}mm^\top N^{-1} \\ &+N^{-1}vm^\top N^{-1}+N^{-1}mv^\top N^{-1} \\ &+ab^\top -a(v^\top+m^\top) N^{-1}-N^{-1}(v+m)b^\top \\ \tag{4}\label{4} \end{align} and $dz=dy/\det N=dv/\det N=dv/\det(4M+I)^{1/2}$.

So, using \eqref{2}, \eqref{3}, and \eqref{3}, and letting $V$ denote a standard normal random vector in $\R^n$, we see that \begin{align} J&=\frac{e^{-c_1+c_2}}{\det(4M+I)^{1/2}}Eh(V). \tag{5}\label{5} \end{align} Noting finally that $EV=0$ and $EVV^\top=I$, from \eqref{5} we get the expression \eqref{1} for the expectation $J$ in question.