Correlation between two nodes of a single layer MLP for joint-Gaussian input

Question

Let's say you have a jointly gaussian vector random variable $\mathbf{x}$, with mean $\mathbf{M}$ and covariance $\mathbf{S}$. I now transform each scalar element of $\mathbf{x}$ , say $x_j$, with a sigmoid:

$$y_j = 1/(1+\exp(-x_j))$$

I am interested in the expectation between two variables $y_j$ and $y_j'$ of the resulting distribution, that is

$$E\{y_j y_j'\}$$

where $E\{\}$ is the expectation operator. Note:

• The whole PDF of $y$ can be computed just by applying a change of variable. Unfortunately the integrals leading to the expectations are intractable (to the best of my knowledge).

• I'm looking for closed-form approximations, no Markov Chain-Monte Carlo, no variational stuff. That is, approximations to the variable change, to the expectation integrals, to the sigmoid, to the resulting PDF, that allow computing $E\{y_j y_j'\}$.

• Some dead ends: Taylor, often used in papers on the topic, is inaccurate by the mile. Gradshteyn and Ryzhik does not seem to contain the integrals.

Answer 1

The question really concerns pairs of normal variates. Let's call them $x_1$ and $x_2$ with means $\mu_i$, standard deviations $\sigma_i$, and correlation $\rho$. Whence their joint pdf is

$$\frac{1}{2 \pi \sqrt{1 - \rho^2} \sigma_1 \sigma_2} e^{-\frac{1}{1-\rho^2} \left(\frac{(x_1 - \mu_1)^2}{2 \sigma_1^2} + \frac{(x_2 - \mu_2)^2}{2 \sigma_2^2} - \frac{\rho (x_1 - \mu_1)(x_2 - \mu_2)}{\sigma_1 \sigma_2}\right)} dx_1 dx_2\text{.}$$

Let $f(x_1,x_2)$ be the product of this with the $y_i$ (as functions of the $x_i$). The first component of the gradient of $\log(f)$ is

$$\frac{\partial \log(f)}{\partial x_1} = \frac{1}{1 + e^{x_1}} + \frac{\rho(\mu_2 - x_2) \sigma_1 + (x_1 - \mu_1)\sigma_2}{(\rho^2-1)\sigma_1^2 \sigma_2},$$

with a similar expression for the second component (via the symmetry achieved by exchanging the subscripts 1 and 2). There will be a unique global maximum, which we can detect by setting the gradient to zero. This pair of nonlinear equations has no closed form solution. It is rapidly found by a few Newton-Raphson iterations. Alternatively, we can linearize these equations. Indeed, through second order, the first component equals

$$\frac{1}{2} + x_1\left(\frac{-1}{4} + \frac{1}{(\rho^2-1)\sigma_1^2}\right) + \frac{-\rho x_2 \sigma_1 + \rho \mu_2 \sigma_1 - \mu_1 \sigma_2}{(\rho^2 -1)\sigma_1^2 \sigma_2}.$$

This gives a pair of linear equations in $(x_1, x_2)$, which therefore do have a closed form solution, say $\hat{x}_i(\mu_1, \mu_2, \sigma_1, \sigma_2, \rho)$, which obviously are rational polynomials.

The Jacobian at this critical point has 1,1 coefficient

$$\frac{e^\hat{x_1}\left(2 - (\rho^2-1)\sigma_1^2 + 2\cosh(\hat{x_1})\right)}{(1+e^\hat{x_1})^2(\rho^2-1)\sigma_1^2},$$

1,2 and 2,1 coefficients

$$\frac{\rho}{\sigma_1 \sigma_2(1 - \rho^2)},$$

and 2,2 coefficient obtained from the 1,1 coefficient by symmetry. Because this is a critical point (at least approximately), we can substitute

$$e^\hat{x_1} = \frac{(\rho^2-1)\sigma_1^2 \sigma_2}{(\mu_2 - \hat{x_2})\rho \sigma_1 + (\hat{x_1} - \mu_1)\sigma_2} - 1$$

and use that also to compute $\cosh(\hat{x_1}) = \frac{e^\hat{x_1} - e^{-\hat{x_1}}}{2}$, with a similar manipulation for $e^\hat{x_2}$ and $\cosh(\hat{x_2})$. This enables evaluation of the Hessian (the determinant of the Jacobian) as a rational function of the parameters.

The rest is routine: the Hessian tells us how to approximate the integral as a binormal integral (a saddlepoint approximation). The answer equals $\frac{1}{2\pi}$ times a rational function of the five parameters: that's your closed form (for what it's worth!).