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

• It may be possible, but efforts to deal with a truncated bivariate normal seem to suggest it may not be easy Commented Mar 22, 2020 at 3:01
• @MisterBlobfish I'm looking into a similar question and wondering whether you found a satisfactory answer in the end? What I'm looking for is to compute $\mathbb{E}(\text{relu}(X)\text{relu}(Y))$ when X, Y are dependent. Commented Feb 21, 2023 at 12:51
• – ABIM
Commented Feb 14 at 1:16

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.