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)