Neural network: Why can NNs have identical structure in the hidden units and still learn? This is the issue that I am still unclear about in Neural Networks, and I was wondering how much of an open issue it still in in Deep Learning research. How can neural networks learn when the nodes in each layer are essentially identical? I found some similar questions, but I was hoping for a more precise and mathematical answer. 
https://stackoverflow.com/questions/20027598/why-should-weights-of-neural-networks-be-initialized-to-random-numbers
QUESTION SETUP
Say I have a very simple fully connected network with two hidden layers, and an input and output layer, such as in the diagram below, taken from this post. 

Now, all of the nodes in the first hidden layer (blue) have the same exact equation, meaning that they will have the same function $z = W_{ij}*X + bias_i$ where $i \in {1...I}$ layers and $j \in {1...J}$ node in that layer. Each of these nodes also has the same non-linear activation function, such as RELU. 
So the nodes in a given layer are essentially identical in structure because the layer is fully connected. The only difference is that the weights are initialized at different values, and hence it takes time for those weights to adapt through backpropagation.
QUESTION
So each node in a given layer has the same exact equation, but varies only in initialization weights and the trajectory that the optimizer takes on that node. So over time, should not all of the weights for each node arrive at the same exact set of weights? If each node is getting the same input values, then why do these nodes not all converge to the same weight matrix over time? Seem like this is the key point of neural networks, that the nodes will specialize to focus on different parts of the task. But I was not sure how this works in practice when all of the nodes seem identical?
Even in CNNs and RNNs they often have fully connected layers at the top. So how is this divergence happening when the node equations are identical? I mean if I ran a bunch of multiple linear regressions on the same $X$ values, I would expect to get the same answer for each run of that regression equation. 
Really I am looking for any research papers or such that anyone can identify on this topic. I would really like to see the mathematical proofs behind this idea. I think that there is a lot going on here, but I would really like to understand this critical topic in training neural nets. 
 A: Yes, it is true that the functional dependence of nodes in a layer on the nodes in a previous layer is the same for all hidden layers (well you can use different activation functions, but let's assume you are using the same one). However, the nodes in different hidden layers are different functions of the input values. 
The general formula for the node values in layer  $l$ given the node valuers in layer $l-1$ is 
\begin{equation} 
  a^{l}_j = \sigma (z^l_j ) = \sigma\left( \sum_k w^{l}_{jk} a^{l-1}_k + b^l_j \right)
\end{equation}
where $w^{l}_{jk}$ are the weights and $b^l_j$ are the biases. 
Consequently, given a set of input values $x_i$ the node values in the first hidden layer are 
\begin{equation} 
  a^{1}_j = \sigma\left( \sum_k w^{1}_{jk} x_k + b^1_j \right)
\end{equation}
and in the second hidden layer the node values are 
\begin{equation} 
  a^{2}_j = \sigma\left( \sum_k w^{2}_{jk} \sigma\left( \sum_k w^{1}_{jk} x_k + b^1_j \right) + b^2_j \right)
\end{equation}
Clearly, the node values are different functions of the input. 
Furthermore, when weights are updated different updates are used for different layers. The general formula for the weight change in a layer $l$ is calculated using the chain rule (assuming the learning rate $\eta=1$ for simplicity):
\begin{eqnarray} 
  \Delta w^{l}_{jk} =  -\delta^l_j  a^{l-1}_k\\
  \delta^l_j = \delta^{l+1}_j \sigma'(z^l_j)
\end{eqnarray}
Starting with the last output layer $L$ we get:
\begin{equation}
  \Delta w^L_{jk} = - \frac{\partial Loss}{\partial  a^{L}_j } \sigma'(z^L_j) a^{L-1}_k
\end{equation}
and for the layer $L-1$ we get:
\begin{equation}
  \Delta w^{L-1}_{jk} = - \frac{\partial Loss}{\partial  a^{L}_j } \sigma'(z^L_j) \sigma'(z^{L-1}_j)a^{L-2}_k
\end{equation}
So the weights in different layers are updated differently.
By the way, that's why a multi-layer neural network is qualitatively different, and better, than a one layer neural network with the same number of nodes.
