Why is it hard to train deep neural networks? I'm fairly new to the field of ML and I often see in the intro of papers a sentence on how deep networks are hard to train or that it's difficult for the training signal to reach through multiple hidden layers.
Is there some source that actually rigorously explains the problem and analyses the situation and how the chosen activation function affects this? I would assume that it has something do with vanishing gradients, but would actually want to see a rigorous analysis of this and how the number of layers impacts the gradients.
I'm a math Ph.D. so I'm looking for something more precise and less handwavy.
 A: Resources

*

*The chapter Why are deep neural networks hard to train? (in the book "Neural Networks and Deep Learning" by Michael Nielsen) is probably the best answer to your question that I encountered, but hopefully my answer would contain the gist of the chapter.

*The paper On the difficulty of training recurrent neural networks contains a proof that some condition is sufficient to cause the vanishing gradient problem in a simple recurrent neural network (RNN). I would give an explanation which is similar to the proof, but for the case of a simple deep feedforward neural network.

*The chapter How the backpropagation algorithm works (in the same book by Nielsen) explains clearly and rigorously how backpropagation works, and I would use its notations, definitions and conclusions in my explanation.

Unstable Gradient Problem
Nielsen claims that when training a deep feedforward neural network using Stochastic Gradient Descent (SGD) and backpropagation, the main difficulty in the training is the "unstable gradient problem". Here is Nielsen's explanation of this problem:

[...] the gradient in early layers is the product of terms from all the later layers. When there are many layers, that's an intrinsically unstable situation. The only way all layers can learn at close to the same speed is if all those products of terms come close to balancing out. Without some mechanism or underlying reason for that balancing to occur, it's highly unlikely to happen simply by chance. In short, the real problem here is that neural networks suffer from an unstable gradient problem. As a result, if we use standard gradient-based learning techniques, different layers in the network will tend to learn at wildly different speeds.

Next, we would use equations that Nielsen proved to show that "gradient in early layers is the product of terms from all the later layers".
For that, we need some notations and definitions:

*

*Layer $1$ is the input layer.

*Layer $L$ is the output layer.

*$x$ is a vector of inputs in a single training example.

*$y$ is a vector of desired outputs in a single training example.

*$a^l$ is a vector of the activations of the neurons in layer $l$.

*$C\equiv\frac{1}{2}||y-a^{L}||^{2}$ is the cost function with regard to a single training example $(x, y)$. (This is a simplification. In a real implementation, we would use mini-batches instead.)

*$w^l$ is a matrix of weights for the connections from layer $l-1$ to layer $l$.

*$b^l$ is a vector of the biases used while computing the weighted inputs to the neurons in layer $l$.

*$z^{l}\equiv w^{l}a^{l-1}+b^{l}$ is a vector of the weighted inputs to the neurons in layer $l$.

*$\sigma$ is the activation function.

*$a^l\equiv \sigma(z^l)$, while $\sigma$ is applied element-wise.

*$\delta^{l}\equiv\frac{\partial C}{\partial z^{l}}$

*$\Sigma'\left(z^{l}\right)$ is a diagonal matrix whose diagonal is $\sigma'(z^l)$ (while $\sigma'$ is applied element-wise).


Nielsen proved the following equations:

*

*(34): $\delta^{l}=\Sigma'\left(z^{l}\right)\left(w^{l+1}\right)^{T}\delta^{l+1}$

*(30): $\delta^{L}=\left(a^{L}-y\right)\odot\sigma'\left(z^{L}\right)$, which is equivalent to $\delta^{L}=\Sigma'\left(z^{L}\right)\left(a^{L}-y\right)$
Thus: $$\delta^{l}=\Sigma'\left(z^{l}\right)\left(w^{l+1}\right)^{T}\cdots\Sigma'\left(z^{L-1}\right)\left(w^{L}\right)^{T}\delta^{L}\\\downarrow\\\delta^{l}=\Sigma'\left(z^{l}\right)\left(w^{l+1}\right)^{T}\cdots\Sigma'\left(z^{L-1}\right)\left(w^{L}\right)^{T}\Sigma'\left(z^{L}\right)\left(a^{L}-y\right)$$
Nielsen also proved:

*

*(BP3): $\frac{\partial C}{\partial b_{j}^{l}}=\delta_{j}^{l}$

*(BP4): $\frac{\partial C}{\partial w_{jk}^{l}}=\delta_{j}^{l}a_{k}^{l-1}$
Therefore (this is my notation, so don't blame Nielsen in case it is ugly):
$$\frac{\partial C}{\partial b^{l}}\equiv\left(\begin{gathered}\frac{\partial C}{\partial b_{1}^{l}}\\
\frac{\partial C}{\partial b_{2}^{l}}\\
\vdots
\end{gathered}
\right)=\delta^{l}$$
$$\frac{\partial C}{\partial w^{l}}\equiv\left(\begin{matrix}\frac{\partial C}{\partial w_{11}^{l}} & \frac{\partial C}{\partial w_{12}^{l}} & \cdots\\
\frac{\partial C}{\partial w_{21}^{l}} & \frac{\partial C}{\partial w_{22}^{l}} & \cdots\\
\vdots & \vdots & \ddots
\end{matrix}\right)=\delta^{l}\left(a^{l-1}\right)^{T}$$
From these conclusions, we deduce the components of the gradient in layer $l$:
$$\frac{\partial C}{\partial b^{l}}=\Sigma'\left(z^{l}\right)\left(w^{l+1}\right)^{T}\cdots\Sigma'\left(z^{L-1}\right)\left(w^{L}\right)^{T}\Sigma'\left(z^{L}\right)\left(a^{L}-y\right)\\\frac{\partial C}{\partial w^{l}}=\frac{\partial C}{\partial b^{l}}\left(a^{l-1}\right)^{T}$$
Indeed, both components (i.e. partial derivatives with regard to weights and biases) of the gradient in layer $l$ are products that include all of the weight matrices of the next layers, and also the derivatives of the activation function of the next layers.
Vanishing Gradient Problem
If you are still not convinced that the "unstable gradient problem" is real or that it actually matters, we would next show why the "vanishing gradient problem" is probable in a deep feedforward neural network.
As in the proof in the paper, we can use vector norms and induced matrix norms to get a rough upper bound on $\left|\left|\frac{\partial C}{\partial b^{l}}\right|\right|$ and $\left|\left|\frac{\partial C}{\partial w^{l}}\right|\right|$.
In the case of induced matrix norms, both $\left|\left|ABx\right|\right|\le\left|\left|A\right|\right|\cdot\left|\left|B\right|\right|\cdot\left|\left|x\right|\right|$ and $\left|\left|AB\right|\right|\le\left|\left|A\right|\right|\cdot\left|\left|B\right|\right|$ hold for any matrices $A,B$ and vector $x$ such that $ABx$ is defined.
Therefore:
$$\begin{gathered}\left|\left|\frac{\partial C}{\partial b^{l}}\right|\right|=\left|\left|\Sigma'\left(z^{l}\right)\left(w^{l+1}\right)^{T}\cdots\Sigma'\left(z^{L-1}\right)\left(w^{L}\right)^{T}\Sigma'\left(z^{L}\right)\left(a^{L}-y\right)\right|\right|\le\\
\left|\left|\Sigma'\left(z^{l}\right)\right|\right|\left|\left|\left(w^{l+1}\right)^{T}\right|\right|\cdots\left|\left|\Sigma'\left(z^{L-1}\right)\right|\right|\left|\left|\left(w^{L}\right)^{T}\right|\right|\left|\left|\Sigma'\left(z^{L}\right)\right|\right|\left|\left|a^{L}-y\right|\right|\\
\downarrow\\
\left|\left|\frac{\partial C}{\partial b^{l}}\right|\right|\le\overset{L}{\underset{r=l}{\prod}}\left|\left|\Sigma'\left(z^{r}\right)\right|\right|\cdot\overset{L}{\underset{r=l+1}{\prod}}\left|\left|\left(w^{r}\right)^{T}\right|\right|\cdot\left|\left|a^{L}-y\right|\right|
\end{gathered}
$$
and also:
$$\begin{gathered}\left|\left|\frac{\partial C}{\partial w^{l}}\right|\right|\le\left|\left|\frac{\partial C}{\partial b^{l}}\right|\right|\left|\left|\left(a^{l-1}\right)^{T}\right|\right|\\
\downarrow\\
\left(*\right)\\
\left|\left|\frac{\partial C}{\partial w^{l}}\right|\right|\le\left|\left|\frac{\partial C}{\partial b^{l}}\right|\right|\left|\left|a^{l-1}\right|\right|
\end{gathered}
$$
It turns out that $||A||=||A^T||$ for any square matrix $A$, as shown here (which uses what is shown here).
Thus:
$$\left|\left|\frac{\partial C}{\partial b^{l}}\right|\right|\le\overset{L}{\underset{r=l}{\prod}}\left|\left|\Sigma'\left(z^{r}\right)\right|\right|\cdot\overset{L}{\underset{r=l+1}{\prod}}\left|\left|w^{r}\right|\right|\cdot\left|\left|a^{L}-y\right|\right|$$
Let $\gamma\equiv\text{sup}\left\{ \sigma'\left(\alpha\right)\,:\,\alpha\in\mathbb{R}\right\} $.
The norm of a diagonal matrix is the largest absolute value of the elements in the matrix. (This is quite immediate from the claim that the norm of a symmetric matrix is equal to its spectral radius.)
So $\left|\left|\Sigma'\left(z\right)\right|\right|\le\gamma$ for any $z$, and thus:
$$\begin{gathered}\left|\left|\frac{\partial C}{\partial b^{l}}\right|\right|\le\overset{L}{\underset{r=l}{\prod}}\gamma\cdot\overset{L}{\underset{r=l+1}{\prod}}\left|\left|w^{r}\right|\right|\cdot\left|\left|a^{L}-y\right|\right|\\
\downarrow\\
\left(**\right)\\
\left|\left|\frac{\partial C}{\partial b^{l}}\right|\right|\le\gamma^{L-l+1}\cdot\overset{L}{\underset{r=l+1}{\prod}}\left|\left|w^{r}\right|\right|\cdot\left|\left|a^{L}-y\right|\right|
\end{gathered}
$$
Now, consider the derivatives of sigmoid (green) and $\text{tanh}$ (red).

In case $\sigma$ is the sigmoid function, $\gamma=0.25$, and so from $(*)$ and $(**)$ we can deduce that $\left|\left|\frac{\partial C}{\partial b^{l}}\right|\right|$ and $\left|\left|\frac{\partial C}{\partial w^{l}}\right|\right|$ would probably be very small for a high $L-l$. I.e. for an early layer in a deep network with many layers, the gradient would be quite small.
$(*)$ and $(**)$ won't help much in showing that the vanishing gradient problem is also probable for the case that $\sigma$ is $\text{tanh}$, but using the same approach and some approximations would work.
