Return to Answer

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^{-1}mm^\top N^{-1}+ab^\top-am^\top N^{-1}-N^{-1}mb^\top).$$ Here one can get rid of $N$ and $m$ by noting 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}$.

$\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.

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$, using the formula $$\int_{\mathbb R^n}e^{-v^\top v/2}\,dv=(2\pi)^{n/2},$$ and doing a bit of algebra, we see that the expectation in question is $$(2\pi)^{n/2}\frac{e^{c_2-c_1}}{\det(4M+I)^{1/2}}(N^{-1}mm^\top N^{-1}+ab^\top-am^\top N^{-1}-N^{-1}mb^\top).$$ Here one can get rid of $N$ and $m$ by noting 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}$.

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^{-1}mm^\top N^{-1}+ab^\top-am^\top N^{-1}-N^{-1}mb^\top).$$ Here one can get rid of $N$ and $m$ by noting 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}$.

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$, using the formula $$\int_{\mathbb R^n}e^{-v^\top v/2}\,dv=(2\pi)^{n/2},$$ and doing a bit of algebra, we see that the expectation in question is $$(2\pi)^{n/2}\frac{e^{c_2-c_1}}{\det(4M+I)^{1/2}}(N^{-1}mm^\top N^{-1}+ab^\top-am^\top N^{-1}-N^{-1}mb^\top).$$ Here one can get rid of $N$ by noting 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}$.

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$, using the formula $$\int_{\mathbb R^n}e^{-v^\top v/2}\,dv=(2\pi)^{n/2},$$ and doing a bit of algebra, we see that the expectation in question is $$(2\pi)^{n/2}\frac{e^{c_2-c_1}}{\det(4M+I)^{1/2}}(N^{-1}mm^\top N^{-1}+ab^\top-am^\top N^{-1}-N^{-1}mb^\top).$$ Here one can get rid of $N$ and $m$ by noting 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}$.