Correlation between two nodes of a single layer MLP for joint-Gaussian input 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.
 A: 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!).
