Gaussian expectation of an exponentiated outer product

Question

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?

Answer 1

The following solves the algebra for element-wise calculations by brute-force, and extracts a 'matrix-form solution' of sorts.

For $i = j$

We seek $$ E\left[ \exp\left( X_i^2 \right) \right] = \frac{1}{\Sigma_i\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{X_i^2} e^{-\frac{\left(X_i - \mu_i\right)^2}{2\Sigma_i^2}} dX_i $$ Observe $$ \begin{align*} X_i^2-\frac{\left(X_i - \mu_i\right)^2}{2\Sigma_i^2} &= \frac{2\Sigma_i^2}{2\Sigma_i^2} X_i^2 - \frac{X_i^2 - 2\mu_i X_i + \mu_i^2}{2\Sigma_i^2} \\ &= -\frac{\left( 1 - 2\Sigma_i^2 \right)X_i^2 - 2\mu_i X_i + \mu_i^2}{2\Sigma_i^2} \\ &= -\frac{\frac{1 - 2\Sigma_i^2}{\Sigma_i^2}X_i^2 - 2\frac{\mu_i}{\Sigma_i^2} X_i + \frac{\mu_i^2}{\Sigma_i^2}}{2} \\ &= -\frac{\alpha_i^2 X_i^2 - 2\frac{\mu_i}{\alpha_i\Sigma_i^2} \alpha_i X_i + \frac{\mu_i^2}{\Sigma_i^2}}{2} \ \text{where $\alpha_i^2 = \frac{1 - 2\Sigma_i^2}{\Sigma_i^2}$} \\ &= -\frac{\left(\alpha_i X_i - \frac{\mu_i}{\alpha_i\Sigma_i^2} \right)^2 - \left( \frac{\mu_i}{\alpha_i\Sigma_i^2} \right)^2 + \frac{\mu_i^2}{\Sigma_i^2}}{2} \\ &= -\frac{Y_i^2}{2} + \gamma_i \end{align*} $$ where $Y_i = \alpha_i X_i - \frac{\mu_i}{\alpha_i\Sigma_i^2}$ and $$ \begin{align*} \gamma_i &= \frac{ \frac{\mu_i^2}{\alpha_i^2\left(\Sigma_i^2\right)^2} - \frac{\mu_i^2}{\Sigma_i^2}}{2} \\ &= \frac{\mu_i^2}{2\Sigma_i^2}\left( \frac{1}{\alpha_i^2\Sigma_i^2} - 1 \right) \\ &= \frac{\mu_i^2}{2\Sigma_i^2}\left( \frac{1}{1 - 2\Sigma_i^2} - 1 \right) \\ &= \frac{\mu_i^2}{2\Sigma_i^2}\left( \frac{1 - 1 + 2\Sigma_i^2}{1 - 2\Sigma_i^2} \right) \\ &= \frac{\mu_i^2}{1 - 2\Sigma_i^2} \end{align*} $$ Then by change of variables $$ \begin{align*} E\left[ \exp\left( X_i^2 \right) \right] &= \frac{1}{\Sigma_i\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-\frac{Y_i^2}{2}}e^{\gamma_i} \frac{1}{\alpha_i} dY_i \\ &= \frac{e^{\gamma_i}}{\alpha_i\Sigma_i}\cdot\frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-\frac{Y_i^2}{2}} dY_i \\ &= \frac{e^{\gamma_i}}{\sqrt{1 - 2\Sigma_i^2}} \end{align*} $$

For $i \not= j$

Recall that under iterated expectations $$ E\left[ \exp\left( X_i X_j \right) \right] = E\left[ E\left[ \exp\left( X_i X_j \right) \vert\ X_j \right] \right] $$ So we first seek $$ E\left[ \exp\left( X_i X_j \right) \vert\ X_j \right] = \frac{1}{\Sigma_i\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{X_i X_j} e^{-\frac{\left(X_i - \mu_i\right)^2}{2\Sigma_i^2}} dX_i $$ Observe $$ \begin{align*} X_i X_j - \frac{\left(X_i - \mu_i\right)^2}{2\Sigma_i^2} &= \frac{2\Sigma_i^2}{2\Sigma_i^2} X_i X_j - \frac{X_i^2 - 2\mu_i X_i + \mu_i^2}{2\Sigma_i^2} \\ &= -\frac{X_i^2 - 2\left(\mu_i - \Sigma_i^2 X_j \right)X_i + \mu_i^2}{2\Sigma_i^2} \\ &= -\frac{\frac{1}{\Sigma_i^2}X_i^2 - 2\frac{\mu_i - \Sigma_i^2 X_j}{\Sigma_i^2} X_i + \frac{\mu_i^2}{\Sigma_i^2}}{2} \\ &= -\frac{\beta_i^2 X_i^2 - 2\frac{\mu_i - \Sigma_i^2 X_j}{\beta_i\Sigma_i^2} \beta_i X_i + \frac{\mu_i^2}{\Sigma_i^2}}{2} \ \text{where $\beta_i^2 = \frac{1}{\Sigma_i^2}$} \\ &= -\frac{\left(\beta_i X_i - \frac{\mu_i - \Sigma_i^2 X_j}{\beta_i\Sigma_i^2} \right)^2 - \left( \frac{\mu_i - \Sigma_i^2 X_j}{\beta_i\Sigma_i^2} \right)^2 + \frac{\mu_i^2}{\Sigma_i^2}}{2} \\ &= -\frac{W_i^2}{2} + \delta_i \\ \end{align*} $$ where $W_i = \beta_i X_i - \frac{\mu_i - \Sigma_i^2 X_j}{\beta_i\Sigma_i^2}$ and $$ \begin{align*} \delta_i &= \frac{\frac{\left( \mu_i - \Sigma_i^2 X_j \right)^2}{\beta_i^2\left(\Sigma_i^2\right)^2} - \frac{\mu_i^2}{\Sigma_i^2}}{2} \\ &= \frac{\left( \mu_i - \Sigma_i^2 X_j \right)^2 - \mu_i^2}{2\Sigma_i^2} \\ &= \frac{\mu_i^2 - 2\mu_i\Sigma_i^2 X_j + \left(\Sigma_i^2\right)^2 X_j^2 - \mu_i^2}{2\Sigma_i^2} \\ &= \frac{ \Sigma_i^2 X_j^2 - 2\mu_i X_j}{2} \\ \end{align*} $$ By change of variables $$ \begin{align*} E\left[ \exp\left( X_i X_j \right) \vert\ X_j \right] &= \frac{1}{\Sigma_i\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-\frac{W_i^2}{2} } e^{\delta_i} \frac{1}{\beta_i} dW_i \\ &= \frac{e^{\delta_i}}{\beta_i\Sigma_i} \cdot \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-\frac{W_i^2}{2}} dW_i \\ &= e^{\delta_i} \end{align*} $$ so $$ E\left[ \exp\left( X_i X_j \right) \right] = E\left[ e^{\delta_i} \right] = E\left[ \exp\left( \frac{ \Sigma_i^2 X_j^2 - 2\mu_i X_j}{2} \right) \right] $$ Observe $$ \begin{align*} \frac{ \Sigma_i^2 X_j^2 - 2\mu_i X_j}{2} - \frac{\left(X_j - \mu_j \right)^2}{2\Sigma_j^2} &= \frac{ \Sigma_i^2\Sigma_j^2 X_j^2 - 2\mu_i\Sigma_j^2 X_j}{2\Sigma_j^2} - \frac{X_j^2 - 2\mu_j X_j + \mu_j^2}{2\Sigma_j^2} \\ &= -\frac{\left( 1 - \Sigma_i^2\Sigma_j^2 \right)X_j^2 - 2\left(\mu_j - \mu_i\Sigma_j^2\right)X_j + \mu_j^2}{2\Sigma_j^2} \\ &= -\frac{\frac{1 - \Sigma_i^2\Sigma_j^2}{\Sigma_j^2}X_j^2 - 2\frac{\mu_j - \mu_i\Sigma_j^2}{\Sigma_j^2}X_j + \frac{\mu_j^2}{\Sigma_j^2}}{2} \\ &= -\frac{\eta_j^2X_j^2 - 2\frac{\mu_j - \mu_i\Sigma_j^2}{\eta_j\Sigma_j^2}\eta_jX_j + \frac{\mu_j^2}{\Sigma_j^2}}{2} \ \text{where $\eta_j^2 = \frac{1 - \Sigma_i^2\Sigma_j^2}{\Sigma_j^2}$} \\ &= -\frac{\left(\eta_j X_j - \frac{\mu_j - \mu_i\Sigma_j^2}{\eta_j\Sigma_j^2}\right)^2 - \left( \frac{\mu_j - \mu_i\Sigma_j^2}{\eta_j\Sigma_j^2} \right)^2 + \frac{\mu_j^2}{\Sigma_j^2}}{2} \\ &= -\frac{V_j^2}{2} + \zeta_j \\ \end{align*} $$ where $V_j = \eta_j X_j - \frac{\mu_j - \mu_i\Sigma_j^2}{\eta_j\Sigma_j^2}$ and $$ \begin{align*} \zeta_j &= \frac{\frac{\left( \mu_j - \mu_i\Sigma_j^2\right)^2}{\eta_j^2\left(\Sigma_j^2\right)^2} - \frac{\mu_j^2}{\Sigma_j^2}}{2} \\ &= \frac{\frac{\left( \mu_j - \mu_i\Sigma_j^2\right)^2}{1 - \Sigma_i^2\Sigma_j^2} - \mu_j^2}{2\Sigma_j^2} \\ &= \frac{\left( \mu_j - \mu_i\Sigma_j^2\right)^2 - \mu_j^2\left(1 - \Sigma_i^2\Sigma_j^2\right)}{2\Sigma_j^2\left(1 - \Sigma_i^2\Sigma_j^2\right)} \\ &= \frac{\mu_j^2 - 2\mu_i\mu_j\Sigma_j^2 + \mu_i^2\left(\Sigma_j^2\right)^2 - \mu_j^2 + \mu_j^2\Sigma_i^2\Sigma_j^2}{2\Sigma_j^2\left(1 - \Sigma_i^2\Sigma_j^2\right)} \\ &= \frac{ - 2\mu_i\mu_j + \mu_i^2\Sigma_j^2 + \mu_j^2\Sigma_i^2}{2\left(1 - \Sigma_i^2\Sigma_j^2\right)} \\ &= \frac{ \mu_i^2\Sigma_j^2 - \mu_i\mu_j }{2\left(1 - \Sigma_i^2\Sigma_j^2\right)} + \frac{ \mu_j^2\Sigma_i^2 - \mu_j\mu_i }{2\left(1 - \Sigma_i^2\Sigma_j^2\right)} \end{align*} $$ Thus with a final change of variables $$ \begin{align*} E\left[ \exp\left( X_i X_j\right) \right] &= \frac{1}{\Sigma_j\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-\frac{V_j^2}{2}}e^{\zeta_j} \frac{1}{\eta_j} dV_j \\ &= \frac{e^{\zeta_j}}{\eta_j\Sigma_j}\cdot\frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-\frac{V_j^2}{2}} dV_j \\ &= \frac{e^{\zeta_j}}{\sqrt{1 - \Sigma_i^2\Sigma_j^2}} \end{align*} $$

For matrix-based calculations

In MATLAB I would proceed as follows: Given column vectors muV and sigmaV for $\mathbf{\mu}$ and $\mathbf{\Sigma}$, we calculate ExxM, the matrix of expected values.

Common elements

n = length(muV);
   muSqrV =    muV.^2;
sigmaSqrV = sigmaV.^2;

For the diagonal elements

oneMinusTwoSigmaSqrV = 1 - 2*sigmaSqrV;
gammaV = muSqrV ./ oneMinusTwoSigmaSqrV;
ExxVondiag = exp(gammaV) ./ sqrt(oneMinusTwoSigmaSqrV);

For the non-diagonal elements

oneMinusSigmaSqrSigmaSqrM = 1 - sigmaSqrV*sigmaSqrV.';
muSqrsigmaSqrProdM = muSqrV*sigmaSqrV.';
mumuProdM = muV*muV.';
etaM = (muSqrsigmaSqrProdM + muSqrsigmaSqrProdM.' - 2*mumuProdM) ./ (2*oneMinusSigmaSqrSigmaSqrM);
ExxMoffdiag = exp(etaM) ./ sqrt(oneMinusSigmaSqrSigmaSqrM);

Then put together as

I = eye(n);
ExxM = diag(ExxVondiag) + (1 - I).*ExxMoffdiag;

Answer 2

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).