Why do we normalize the pre-synaptic values rather than the output of the activation function when using batch normalization? I'm trying to make sense of this batch normalization (1) paper, in Section 3.2, it says

We could have also normalized the layer inputs u, but since u is likely
the output of another nonlinearity, the shape of its distribution is likely to change during training, and constraining its first and second moments would not eliminate the covariate shift.  In contrast, Wu + b is more likely to have a symmetric, non-sparse distribution, that is “more Gaussian” (Hyvarinen & Oja, 2000); normalizing it is likely to produce activations with a stable distribution.

Why the output of a non-linearity $u$ is likely to change during training (as opposed to $Wu+b$)?
Because the input of the non-linearity is in the form of $Wu+b$ as well, if the input is stable or “more Gaussian”, why is the output likely to change?

1: Ioffe S. and Szegedy C. (2015),
"Batch Normalization: Accelerating Deep Network Training by Reducing
Internal Covariate Shift",
Proceedings of the 32nd International Conference on Machine Learning, Lille, France, 2015.
Journal of Machine Learning Research: W&CP volume 37
 A: The shape of the distribution of $Wu+b$ being not likely to change is because of the central limit theorem, which states that,

given certain conditions, the arithmetic mean of a sufficiently large number of iterates of independent random variables, each with a well-defined (finite) expected value and finite variance, will be approximately normally distributed, regardless of the underlying distribution.

There's a nice illustration on the wikipedia page Illustration of the central limit theorem,




*

*Original probability density function           2. Density of a sum of two variables






*Density of a sum of three variables            4. Density of a sum of four variables



As the shape of distribution of $Wu+b$ is "more Gaussian", after applying batch normalization it will become close to $N(0,1)$, regardless what the distribution of $u$ is. 
So by "mapping" some internal distributions to $N(0,1)$ batch normalization helps addressing the internal covariate shift problem.
On the other hand if we apply batch normalization to $u$ or $\sigma(Wu+b)$
instead, as the shape of the distribution is likely to change, normalizing the first and second moments will not result in a stable distribution.
