What are the theoretical/practical reasons to use normal distribution to initialize the weights in Neural Networks? I'm aware that there are many different practices of initializing the weights when training a neural network. It seems traditionally standard normal distribution is the first choice. Most articles I found argue there are better ways to initialize the weights other than using normal distrubtion, but they did not explain why normal distribution would at least work.
(1) I think restricting the weights to have mean at 0 and std at 1 can make the weights as small as possible, which make it convenient for regularization. Am I understanding it correctly?
(2) On the other hand, what are the theoretical/practical reasons to use the normal distribution? Why not sampling random weights from any other arbitrary distributions? Is it because normal distribution has the maximum entropy given the mean and variance? Having the maximum entropy means it's most possible chaotic and thus making least assumptions about the weights. Am I understanding it correctly?
 A: 
(1) I think restricting the weights to have mean at 0 and std at 1 can
  make the weights as small as possible, which make it convenient for
  regularization. Am I understanding it correctly?

No, setting them all to 0 would make them as small as possible. 

(2) On the other hand, what are the theoretical/practical reasons to
  use the normal distribution? Why not sampling random weights from any
  other arbitrary distributions? Is it because normal distribution has
  the maximum entropy given the mean and variance? Having the maximum
  entropy means it's most possible chaotic and thus making least
  assumptions about the weights. Am I understanding it correctly?

I don't think there's too much logic in that decision, perhaps besides the fact that the gaussian distribution is a good "default prior" as many things follow a gaussian distribution. In fact one popular default initialization scheme by Glorot et. al prescribes a uniform distribution, not a normal distribution.
In fact what probably happens is 1. authors provide a theoretical justification for what variance the distribution of the initial weights should be. 2. they choose an arbitrary distribution with that variance. Of course the normal distribution is then a very natural and easy choice!
A: The other answer is good (+1), but just to add to it:
(1) No, one can easily make the weights even smaller by choosing a smaller $\sigma$ than 1. I do think the fact that $L_1$ and $L_2$ weight decay are very common is related to this, in the sense that initializing with $\mu\ne 0$ would be wasteful, as the weight decay penalty would start off higher than necessary (and the weights would then have to slowly "migrate" towards zero mean anyway!). Of course, you could center the weight decay at another value (e.g. $L(w) = ||W - \mu_M||$ for some non-zero $\mu_M$), but this unnecessarily complicated and aesthetically unpleasant. Nevertheless, I don't think it is the cause of it.
(2) Firstly, not all initializations use the normal distribution. Sometimes they use uniform, or in some cases (resnets, some normalizations, etc...) they use some fixed specialized value.
As for the maximum entropy (ME) assumption, I am not sure if this is related (may well be though). ME is true only for that fixed variance. So the question is still why you would want $\sigma=1$. (and also why ME would be preferred either).
Two things come to mind (mostly just gut feelings/things to look into):
(a) There is plenty of theoretical work linking neural nets to Bayesian neural networks (BNNs), and then those to Gaussian processes. For example, see here. A major application is calibrated uncertainty estimation (e.g. here). The idea is that instead of having a single weight value, you have a distribution over weight values. It is common to use Gaussian distributions to (variationally) approximate these distributions, or (more efficiently) use a combination of regularization and noise to approximate it (e.g., here). Usually these approximate having a (variational, factored) Gaussian distribution over the weights.
I suspect that initializing with $\mathcal{N}(0,1)$ is essentially part of the Bayesian  prior in those cases (as you probably know, $L_2$ weight decay regularization is equivalent to a Bayesian Gaussian weight prior). So perhaps Gaussian initialization can be viewed as related to viewing a standard NN as some kind of approximate BNN. It might be worth looking into. 
(b) Normalization in NNs is now standard, in order to prevent problems with gradients in backpropagation. Such layers (e.g., instance, batch, group, layer, etc... norms) usually operate on activations. But there is also the famous weight normalization method, of course. The idea is that you want a "nice" distribution for backprop: too much variance could mean very uneven gradient updates (instability), while too little is also an issue (could be vanishing gradients, or that the weights are not learning different things). You can get this by controlling the weight values or the activation values, or even using methods like SELU.
I suspect $\sigma=1$ just happens to be a nice value (theoretically and practically).
Of course, for initialization, you want to start out close to this distribution.
Hence why we use it (well, sometimes anyway). Of course, choosing $\sigma=0.9$ or $2.1$ would probably be fine (up to a point). Ultimately, it's probably aesthetics (would loved to be proved wrong though).
