Covariance of a rectified (relu) Gaussian Given a normal random vector $$X\sim N(\mu,\Sigma)$$ for spd $\Sigma$, I'm interested in the covariance matrix $K=\mathrm{cov}(Y)$ of the variable $$Y = \mathrm{relu}(X)$$ where the relu is performed elementwise $Y_i = \mathrm{Max}(0,x_i)$, so $Y$ is distributed according to the rectified Gaussian distribution.
Given I know everything about $\Sigma$, how can I compute $K$? 
The mean and variance of each $Y_i$ has been covered in other questions on this site, but the off-diagonal elements of $K$ seem pretty challenging to compute, and I haven't found anything on SO or elsewhere online about it.
I'm actually after the eigenvectors of $K$, so if anyone can relate the eigenvectors between $\Sigma$ and $K$ without directly computing $K$, that would be even more interesting.
Thanks!
Edit: Just to note there is a similar question asked here and thoroughly answered, but only in the scalar (or diagonal multivariate) case. For multi-dimensional $X$ with correlations, this seems much more challenging.
 A: In the 2d case,
let
$$
\Sigma=\begin{pmatrix}\sigma_x^2 & r\sigma_x\sigma_y \\ r\sigma_x\sigma_y & \sigma_y^2\end{pmatrix}.
$$
Then by Mathematica (or this paper, table 1),
$\mathrm{E}[\max(X,0)\cdot\max(Y,0)] = \sigma_x\sigma_y f(r),$
where
$$f(r) = \frac{2\sqrt{1-r^2}+r(\pi+2\arctan(\tfrac{r}{\sqrt{1-r^2}}))}{4\pi}.$$
We also have $E[\max(X,0)]=\frac{\sigma_x}{\sqrt{2\pi}}$.
Back to the general case, let $\Sigma \in \mathbb{R}^{n\times n}$ be the covariance matrix, and $\sigma$ be the diagonal of $\Sigma^{1/2}$.
\begin{align}\mathrm{V}(\mathrm{relu}(X))
&=\mathrm{E}(\mathrm{relu}(X)\mathrm{relu}(X)^T)
-\mathrm{E}(\mathrm{relu}(X))\mathrm{E}(\mathrm{relu}(X)^T)
\\&=\sigma\sigma^T\circ f(\Sigma/(\sigma\sigma^T)) - \frac{\sigma\sigma^T}{2\pi}.
\end{align}

We use the bounds $r/4+1/(2\pi)\le f(r)\le (1+r)/4$,
we get
\begin{align}\mathrm{V}(\mathrm{relu}(X))
&\le\sigma\sigma^T\circ (1+\Sigma/(\sigma\sigma^T))/4 - \frac{\sigma\sigma^T}{2\pi}
\\&=
\sigma\sigma^T(1/4-1/(2\pi)) + \Sigma/4.
\end{align}
and
\begin{align}\mathrm{V}(\mathrm{relu}(X))
&\ge\sigma\sigma^T\circ (1/(2\pi)+\Sigma/(4\sigma\sigma^T)) - \frac{\sigma\sigma^T}{2\pi}
\\&=
\Sigma/4.
\end{align}
Or in other words
$$0 \le \mathrm{V}(\mathrm{relu}(X)) - \Sigma/4 \le 0.1 \sigma\sigma^T.$$
In my own experiments, $\Sigma/4$ seems like a pretty good approximation.

In all of the above we have assumed $\mu=0$.
If this is not the case things get even more hairy.
I would assume, that if $(\mu_x, \mu_y)$ is in the positive quadrant, then most of the normal distribution will be within the non-zero region of the relu, so you can basically ignore the relu.
Meanwhile, if $(\mu_x, \mu_y)$ is in another quadrant, then only a very small part of the normal distribution will be in the non-zero region. Because of the sharp tails of the distribution, we can mostly focus on what happens on the border region of the $(\mu_x, \mu_y)$ quadrant and the positive one.
