Gaussian expectation of an exponentiated outer product Given a normal random column vector $\mathbf{x} \sim N(\mu, \Sigma)$, I need the expectation,
$$ E\left[ \exp(\mathbf{xx}^\top)\right]$$
where $\exp(\cdot)$ is element-wise exponential function (not a matrix exponential).
Is there a closed form for this expression?
I know that the inner product form has a closed form:
$$ E\left[ \exp(\mathbf{x}^\top A \mathbf{x})\right] = |I - 2A\Sigma|^{-\frac{1}{2}} \exp\left[ -\frac{1}{2} \mu^\top (I - (I - 2A\Sigma)^{-1})\Sigma^{-1}\mu \right]$$
for a real symmetric matrix $A$. Since each element in the resulting expectation is an exponentiated quadratic function, I feel like there should be a closed form solution, but my Matrix-fu is not strong enough.
(Context: this result is needed to derive a statistical estimator for a state-space model. Eventually, I need to numerically evaluate this expression.)
EDIT:
Note that
$$ (\mathbf{xx^\top})_{ij} = \mathbf{x^\top}A\mathbf{x}$$
where $A = \frac{1}{2}(J_{(i,j)} + J_{(j,i)})$, and $J_{(i,j)}$ is a matrix with zeros except a 1 at $(i,j)$. So each entry is computable, but can it be simplified to allow matrix form evaluation?
 A: If a 'closed form' solution is allowed to be an infinite series then ...
Let $X_i$ be the random variable for the i-th row of $\mathbf{x}$ where
$$
  X_i \sim N\left( \mu_i, \Sigma_i \right)
$$
We seek 
$ E \left[ \exp\left( X_i X_j \right) \right] $ for each $i, j$.  
Set $Y_{ij} = X_i X_j$ and observe that
$$
  E \left[ \exp\left( Y_{ij} \right) \right] 
= E \left[ \sum_{n=0}^{\infty} \frac{Y_{ij}^n}{n!} \right]
= \sum_{n=0}^{\infty} \frac{E \left[ Y_{ij}^n \right]}{n!} 
$$
Now to obtain $E \left[ Y_{ij}^n \right]$, the suggested strategy is to use moment-generating functions.
For $i = j$
We have $Y_{ii} = X_i^2 $ and
$$
  \frac{X_i^2}{\Sigma_i^2} \sim \chi_1^2 \left( \lambda \right)
$$
where $\chi_k^2$ denotes the noncentral chi-squared distribution with $k$ degrees of freedom and non-centrality parameter $ \lambda = \frac{\mu_i^2}{\Sigma_i^2} $.  Setting $\frac{Y_{ii}}{\Sigma_i^2} = \frac{X_i^2}{\Sigma_i^2} = Q$ with $Q \sim \chi_k^2$ yields
$$
  E \left[ Y_{ii}^n \right] = \left( \Sigma_i^2 \right)^n E \left[ Q^n \right]
$$
where $E \left[ Q^n \right]$ is obtained from the moment-generating function (Wikipedia)
$$ 
  M(t; k; \lambda) =  \frac{\exp\left(\frac{\lambda t}{1 - 2t}\right)}{\left( 1 - 2t \right)^{k/2}} \ \text{if $2t < 1$} 
$$
For $i \not= j$
We have $Y_{ij} = X_i X_j $. Set $X'_i = \frac{X_i - \mu_i}{\Sigma_i}$, $X'_j = \frac{X_j - \mu_j}{\Sigma_j} $ and $Z = X'_i \cdot X'_j$. Then
$$
  Z \sim \text{Product Normal Distribution}
$$
if $X'_i$ and $X'_j$ are independent.  The Product Normal Distribution has characteristic function (Wikipedia)
$$
  \varphi_Z(t) = \frac{1}{\left( 1 + t^2 \right)^{1/2}}
$$
which should provide $E \left[ Z^n \right]$ (after appropriate adjustment from characteristic function to moment-generating function).
