Expected value of error in neural network I wanted to take a look at the properties of the error vector that is propagating during backpropagation. The error vector $\boldsymbol{\delta}$ at layer $i$ is nothing more than the derivative of some loss function $L$ with regard to the inputs $\boldsymbol{z}$ at that layer:
$$\boldsymbol{\delta}^{(i)} = \frac{\partial L(\boldsymbol{y}, \boldsymbol{g})}{\partial \boldsymbol z^{(i)}}$$
where $\boldsymbol{y}$ is the expected output and $\boldsymbol{g}$ is the output of the network.
This formula can be decomposed by means of the chain rule, assuming $L$ hidden layers (layer $0$ is the input layer and layer $L+1$ is the output layer)
$$\boldsymbol{\delta}^{(i)} = \frac{\partial L(\boldsymbol{f}, \boldsymbol{g})}{\partial \boldsymbol{g}} \frac{\partial \boldsymbol{g}}{\partial \boldsymbol{z}^{(L+1)}} \frac{\partial \boldsymbol{z}^{(L+1)}}{\partial \boldsymbol{z}^{(L)}} \ldots \frac{\partial \boldsymbol{z}^{(i+2)}}{\partial \boldsymbol{z}^{(i+1)}} \frac{\partial \boldsymbol{z}^{(i+1)}}{\partial \boldsymbol z^{(i)}}$$
which makes it easy to see that following recursion holds:
$$\boldsymbol{\delta}^{(i)} = \boldsymbol{\delta}^{(i+1)} \frac{\partial \boldsymbol{z}^{(i+1)}}{\partial \boldsymbol{z}^{(i)}}$$
Also note that 
$$\frac{\partial \boldsymbol{z}^{(i+1)}}{\partial \boldsymbol{z}^{(i)}} = \frac{\partial \boldsymbol{z}^{(i+1)}}{\partial \boldsymbol{h}^{(i)}}\frac{\partial \boldsymbol{h}^{(i)}}{\partial \boldsymbol{z}^{(i)}} = {\mathbf{W}^{(i)}}^T f'(\boldsymbol{z}^{(i)})$$
with $f : \mathbb{R}^n \to \mathbb{R}^n$ an (element-wise) activation function and $\boldsymbol{h}^{(i)} = f(\boldsymbol{z}^{(i)})$ the output of layer $i$. $\mathbf{W}^{(i)}$ is a matrix holding the weights between layer $i$ and layer $i+1$.
Now I would be interested in computing 
$$\mathrm{E}\left[\delta_j^{(i)}\right] = \mathrm{E}\left[f'(z_j^{(i)}) \cdot {\boldsymbol{w}^{(i)}}^T_j \cdot \boldsymbol{\delta}^{(i+1)}\right] $$
but due to the dependence of $\boldsymbol{\delta}^{(i+1)}$ with all inputs and weights, there is no way to simplify this and I would need to find an expression of $\delta$ as a function of $z$ and integrate over the whole thing, which I am not sure is possible.
I also tried to compute the expectation as
$\begin{align*}
  \mathrm{E}\left[\delta_j^{(i)}\right] & = \mathrm{E}\left[\frac{\partial L}{\partial z_j^{(i)}}\right] \\
  & = \int_{-\infty}^{\infty} \frac{\partial L}{\partial z} p(z; \mu, \sigma^2) \mathrm{d}z \tag{$z = z_j^{(i)}$} \\
  & = \frac{1}{\sqrt{2 \pi} \sigma} \int_{-\infty}^{\infty} \frac{\partial L}{\partial z} \exp\left(-\frac{(z - \mu)^2}{2 \sigma^2}\right) \mathrm{d}z \tag{Gaussian assumption} \\
  & = \frac{1}{\sqrt{2 \pi} \sigma} \left[ \left. L(\boldsymbol{y}, \boldsymbol{g}) \cdot \exp\left(-\frac{(z - \mu)^2}{2 \sigma^2}\right) \right|^\infty_{-\infty} \right. \\
  & \qquad \left. + \frac{1}{\sigma^2} \int_{-\infty}^{\infty} L(\boldsymbol{y}, \boldsymbol{g}) \cdot (z - \mu) \exp\left(-\frac{(z - \mu)^2}{2 \sigma^2}\right) \mathrm{d}z \right]
\end{align*}$
but this leads to an expression that does not seem to be very integrable, because the network output looks something like $$\boldsymbol{g} = f(\mathbf{W}^{(L)} \cdot f(\mathbf{W}^{(L-1)} \cdot f( \cdots f(\boldsymbol{z}) \cdots )))$$
I am not really a true mathematician, so I hope I did not make too much mistakes in my explanation/notation. This is probably also the reason why I am not completely sure how I could proceed and whether it even makes sense to proceed. Is there any chance I could compute an expectation like this or should I just give up? Are there other tricks I could consider to try to compute this integral?
Any help or hint would be greatly appreciated (also if it concerns my question not being clear enough)
 A: This problem isn't simple to approach. The reason that backpropagation works well for neural networks is because we can pass the inputs through the network once, record the activation value of each neuron, and use them to compute the derivatives of each neuron.
When you integrate over $z$ you are changing the activation value of a neuron which will change every neuron in every subsequent layer. This is very unlike backpropagation where we get the partial derivatives for a single activation value for each neuron.
There is no reason to expect that this would have a nice solution so my honest advice to you would be to give up.
A: Eventually I ended up assuming that everything is independent, i.e.,
$$\mathop{\mathrm{E}}[\delta_j^{(i)}] = \mathop{\mathrm{E}}[f'\big(z_j^{(i)}\big) \cdot \boldsymbol{w}_j^{(i)} \cdot \boldsymbol{\delta}^{(i+1)}] \approx \mathop{\mathrm{E}}[f'\big(z_j^{(i)}\big)] \cdot \mathop{\mathrm{E}}[\boldsymbol{w}_j^{(i)}] \cdot \mathop{\mathrm{E}}[\boldsymbol{\delta}^{(i+1)}].$$
Of course, this is rather inaccurate, but the theoretical results I obtained using this approach turned out to be good enough in empirical experiments.
In the meantime I also went to the struggle of computing the expectation with dependencies for the simplest possible case: the delta for the inputs to the last layer in a network that is to be trained with the MSE loss function at initialisation time. In this case, $\boldsymbol{\delta}^{(i + 1)} = \boldsymbol{z}^{(i)} - \boldsymbol{y} = \boldsymbol{W}^{(i)} \cdot f(\boldsymbol{z}^{(i)}) - \boldsymbol{y}$. As a result, the expectation can be computed as follows:
$$\begin{align}
  \mathop{\mathrm{E}}[\delta_j^{(i)}] &= \mathop{\mathrm{E}}[f'\big(z_j^{(i)}\big) \boldsymbol{w}_j^{(i)} \cdot \boldsymbol{\delta}^{(i+1)}] \\
&= \mathop{\mathrm{E}}[f'\big(z_j^{(i)}\big) \boldsymbol{w}_j^{(i)} \cdot \big(\boldsymbol{W}^{(i)} \cdot f(\boldsymbol{z}^{(i)}) - \boldsymbol{y}\big)] \\
&= \mathop{\mathrm{E}}[f'\big(z_j^{(i)}\big) \boldsymbol{w}_j^{(i)} \cdot \boldsymbol{W}^{(i)} \cdot f(\boldsymbol{z}^{(i)})] - \mathop{\mathrm{E}}[f'\big(z_j^{(i)}\big) \boldsymbol{w}_j^{(i)} \cdot \boldsymbol{y}] \\
&= \mathop{\mathrm{E}}[\boldsymbol{w}_j^{(i)} \cdot \boldsymbol{W}^{(i)}] \cdot \mathop{\mathrm{E}}[f'\big(z_j^{(i)}\big) f(\boldsymbol{z}^{(i)})] - \mathop{\mathrm{E}}[f'\big(z_j^{(i)}\big)] \mathop{\mathrm{E}}[\boldsymbol{w}_j^{(i)}] \cdot \mathop{\mathrm{E}}[\boldsymbol{y}]
\end{align}$$
If we furthermore assume i.i.d. weights, such that $\mathop{\mathrm{E}}[w_{ij}] = \mu_w$ and $\mathop{\mathrm{E}}[w_{ij} w_{kl}] = \delta_{ik}\delta_{jl} \sigma_w^2 + \mu_w^2$, this becomes
$$\begin{align}
  \mathop{\mathrm{E}}[\delta_j^{(i)}] &= M (\sigma_w^2 \boldsymbol{e}_i + \mu_w^2 \boldsymbol{1}) \cdot \mathop{\mathrm{E}}[f'(z_j^{(i)}) f(\boldsymbol{z}^{(i)})] - \mu_w \boldsymbol{1} \cdot \mathop{\mathrm{E}}[\boldsymbol{y}] \mathop{\mathrm{E}}[f'\big(z_j^{(i)}\big)] \\
  &= M \sigma_w^2 \mathop{\mathrm{E}}[f'(z_j^{(i)}) f(z_j^{(i)})] + \mu_w \left(\sum_{k=1}^M \bigg[\mu_w \sum_{l=1}^H \mathop{\mathrm{E}}[f'(z_j^{(i)}) f(z_l^{(i)})] - \mathop{\mathrm{E}}[f'(z_j^{(i)})] \mathop{\mathrm{E}}[y_k]\bigg]\right),
\end{align}$$
with $M$ and $H$ the number of outputs and hidden units, respectively. I.e., $\boldsymbol{y} \in \mathbb{R}^M$ and $\boldsymbol{z} \in \mathbb{R}^H$.
Note that only if $M \sigma_w^2 \mathop{\mathrm{E}}[f'(z_j) f(z_j)] \approx 0$ and
$$\mu_w \left(\frac{\mathop{\mathrm{E}}[f'(z_j) f(z_j)]}{\mathop{\mathrm{E}}[f'(z_j)]} + \sum_{l=1, l \neq j}^H \mathop{\mathrm{E}}[f(z_l^{(i)})]\right) - \mathop{\mathrm{E}}[y_k] \approx \mathop{\mathrm{E}}[\delta_k^{(i+1)}],$$
the independence assumption gives a good approximation.
