Expectation of exponential of product of correlated multivariate Gaussian variables Let $x$ and $y$ be two multivariate Gaussian variables, $x\sim\mathcal{N}(\mu_x,\Sigma_x)$, $y\sim\mathcal{N}(\mu_y,\Sigma_y)$. What is the value of the expectation.
$$
\operatorname{E}\left(e^{x^\top y}\right)
$$
More generally, $x$ and $y$ might be correlated. In which case
$$
\begin{bmatrix} x \\ y \end{bmatrix}
\sim\mathcal{N}\left(
\begin{bmatrix} \mu_x \\ \mu_y \end{bmatrix},
\begin{bmatrix} \Sigma_{xx} & \Sigma_{xy} \\ \Sigma_{xy}^\top & \Sigma_{yy} \end{bmatrix}
\right)
$$ 
At first glance, this seems similar to the expectation $\operatorname{E}\left(e^{x}\right)$, which can be evaluated by completing the square in a Gaussian integral. However, $\operatorname{E}\left(e^{x^\top y}\right)$ has quadratic terms. Does anyone know if this has a closed-form solution, and if so what a good reference might be?
 A: We can use the fact that
$$E\left( e^{x^\top y} \right) 
= E_x\left( E_y\left( e^{x^\top y} \mid x \right) \right)$$
The conditional distribution of $y$ given $x$ is
$$p(y \mid x) = \mathcal{N}\left( \mu_y + \Sigma_{xy}\Sigma_{yy}^{-1} \mu_x, \Sigma_{yy} - \Sigma_{xy}^\top \Sigma_{xx}^{-1} \Sigma_{xy} \right)$$
and so for fixed $x$
$$x^\top y \sim \mathcal{N}(
    \underbrace{x^\top \overbrace{\left( \mu_y + \Sigma_{xy}\Sigma_{yy}^{-1} \mu_x \right)}^u}_{\mu(x)},
    \underbrace{x^\top \overbrace{\left( \Sigma_{yy} - \Sigma_{xy}^\top \Sigma_{xx}^{-1} \Sigma_{xy} \right)}^{A} x}_{\sigma^2(x)} ).$$
For fixed $x$, $e^{x^\top y}$ is therefore log-normally distributed and we have
$$E_y\left(e^{x^\top y} \mid x \right) = e^{\mu(x) + \sigma^2(x) / 2}.$$
By completing the square and using the known normalization constant of the Gaussian distribution, it is possible to derive a closed-form solution. Perhaps somebody else can provide a more elegant or compact solution than my brute-force attempt: 
\begin{align}
&\, E_x\left( \exp \left( \mu(x) + \sigma^2(x) / 2 \right) \right) \\
=& \frac{1}{|2\pi\Sigma_{xx}|^{1/2}} \int \exp\left( -\frac{1}{2} (x - \mu_x)^\top \Sigma_{xx}^{-1} (x - \mu_x) + x^\top u + \frac{1}{2} x^\top A x \right) \, dx \\
=& \frac{1}{|2\pi\Sigma_{xx}|^{1/2}} \int \exp\left( -\frac{1}{2} x^\top \Sigma_{xx}^{-1} x + x^\top \Sigma_{xx}^{-1} \mu_x - \frac{1}{2} \mu_x^\top \Sigma_{xx}^{-1} \mu_x + x^\top u + \frac{1}{2} x^\top A x \right) \, dx \\
=& \frac{1}{|2\pi\Sigma_{xx}|^{1/2}} \int \exp\left( -\frac{1}{2} x^\top (\Sigma_{xx}^{-1} - A) x + x^\top (\Sigma_{xx}^{-1} \mu_x + u) - \frac{1}{2} \mu_x^\top \Sigma_{xx}^{-1} \mu_x \right) \, dx \\
=& \frac{1}{|2\pi\Sigma_{xx}|^{1/2}} \int \exp( -\frac{1}{2} x^\top \underbrace{(\Sigma_{xx}^{-1} - A)}_{B} x + x^\top \underbrace{(\Sigma_{xx}^{-1} \mu_x + u)}_{Bv} - \frac{1}{2} \mu_x^\top \Sigma_{xx}^{-1} \mu_x ) \, dx \\
=& \frac{1}{|2\pi\Sigma_{xx}|^{1/2}} \int \exp\left( -\frac{1}{2} x^\top B x + x^\top B v - \frac{1}{2} v^\top B v + \frac{1}{2} v^\top B v - \frac{1}{2} \mu_x^\top \Sigma_{xx}^{-1} \mu_x \right) \, dx \\
=& \frac{1}{|2\pi\Sigma_{xx}|^{1/2}} \int \exp\left( -\frac{1}{2} (x - v)^\top B (x - v) + \frac{1}{2} v^\top B v - \frac{1}{2} \mu_x^\top \Sigma_{xx}^{-1} \mu_x \right) \, dx \\
=& \frac{|2\pi B^{-1}|^{1/2}}{|2\pi\Sigma_{xx}|^{1/2}} \exp\left( \frac{1}{2} v^\top B v - \frac{1}{2} \mu_x^\top \Sigma_{xx}^{-1} \mu_x \right)
\end{align}
