What are good initial weights in a neural network? I have just heard, that it's a good idea to choose initial weights of a neural network from the range $(\frac{-1}{\sqrt d} , \frac{1}{\sqrt d})$, where $d$ is the number of inputs to a given neuron. It is assumed, that the sets are normalized - mean 0, variance 1 (don't know if this matters).
Why is this a good idea?
 A: The idea is that you want to initialize the weights in a way that ensures good forward and backward data flow through the network. That is, you don't want the activations to be consistently shrinking or increasing as you progress through the network.
This image shows the activations of a 5 layer ReLU Multi-Layer Perceptron under 3 different initialization strategies after one pass of MNIST through the network.

In all three cases weights are drawn from a zero-centered normal distribution which is determined by its standard deviation. You can see that if the initial weights are too small (the standard deviation is small) the activations get choked, and that if they are too large the activations explode. The middle value, that is approximately right can be found by setting the weights such that the variance of the activations and gradient updates stays approximately the same as you pass through the network.
I wrote a blog post about weight initialization that goes into more detail, but the basic idea is as follows.
If $x^{(i)}$ denotes the activations of the $i$-th layer, $n_i$ the size of the layer, and $w^{(i)}$ the weights connecting them to the $(i+1)$-st layer, then one can show that for activation functions $f$ with $f'(s) \approx 1$ we have
$$
\text{Var}(x^{(i+1)}) = n_i \text{Var}(x^{(i)}) \text{Var}(w^{(i)})
$$
In order to achieve $\text{Var}(x^{(i+1)}) = \text{Var}(x^{(i)})$ we therefore have to impose the condition
$$
\text{Var}( w^{(i)}) = \frac{1}{n_i}\,.
$$
If we denote $\frac{\partial L}{\partial x_j^{(i)}}$  by $\Delta_j^{(i)}$, on the backward pass we similarly want
$$
\text{Var}(\Delta^{(i)} ) = n_{i+1} \text{Var}(\Delta^{(i+1)}) \text{Var}(w^{(i)})\,.
$$
Unless $n_i = n_{i+1}$, we have to compromise between these two conditions, and a reasonable choice is the harmonic mean
$$
\text{Var}(w^{(i)}) = \frac{2}{n_i+n_{i+1}}\,.
$$
If we sample weights from a normal distribution $N(0, \sigma) $ we satisfy this condition with $\sigma = \sqrt{\frac{2}{n_i + n_{i+1}}} $.
For a uniform distribution $U(-a, a) $ we should take $a  = \sqrt{\frac{6}{n_i+n_{i+1}}} $ since  $\text{Var} \left( U(-a,a) \right) = a^2/3 $.
We have thus arrived at Glorot initialization. This is the default initialization strategy for dense and 2D convolution layers in Keras, for instance.
Glorot initialization works pretty well for trivial and $ \tanh $ activations, but doesn't do as well for $ \text{ReLU} $.
Luckily, since $ f(s) = \text{ReLU}(s) $ just zeroes out negative inputs, it roughly removes half the variance and this is easily amended by multiplying one of our conditions above by two:
$$
\text{Var}(w^{(i)}) = \frac{2}{n_i}\,.
$$
A: I assume you are using logistic neurons, and that you are training by gradient descent/back-propagation. 
The logistic function is close to flat for large positive or negative inputs. The derivative at an input of $2$ is about $1/10$, but at $10$ the derivative is about $1/22000$ . This means that if the input of a logistic neuron is $10$ then, for a given training signal, the neuron will learn about $2200$ times slower that if the input was $2$. 
If you want the neuron to learn quickly, you either need to produce a huge training signal (such as with a cross-entropy loss function) or you want the derivative to be large. To make the derivative large, you set the initial weights so that you often get inputs in the range $[-4,4]$. 
The initial weights you give might or might not work. It depends on how the inputs are normalized. If the inputs are normalized to have mean $0$ and standard deviation $1$, then a random sum of $d$ terms with weights uniform on $(\frac{-1}{\sqrt{d}},\frac{1}{\sqrt{d}})$ will have mean $0$ and variance $\frac{1}{3}$, independent of $d$. The probability that you get a sum outside of $[-4,4]$ is small. That means as you increase $d$, you are not causing the neurons to start out saturated so that they don't learn. 
With inputs which are not normalized, those weights may not be effective at avoiding saturation. 
A: One other technique that alleviates the problem of weight initialization is Batch Normalization. It acts to standardize the mean and variance of each unit in order to stabilize learning as described in the original paper. In practice, networks that use Batch Normalization (BN) are significantly more robust to bad initialization. BN works as follows:
$$
\mu_B = \frac{1}{m}\sum_{i=1}^{M}x_i~~~and~~~ \sigma_{B}^{2} = \frac{1}{m}\sum_{i=1}^{m}(x_i - \mu_B)^{2} \\
\hat{x}_i = \frac{x_i - \mu_B}{\sqrt{\sigma_{B}^{2} + \epsilon}}~~~and~~~BN(x_i) = \gamma \hat{x}_i + \beta
$$
We compute empirical mean and variance for each mini-batch, then we standardize the input $x_i$ and form the output $BN(x_i)$ by scaling $\hat{x}_i$ by $\gamma$ and adding $\beta$ both of which are learned during training. 
BN introduces two extra parameters ($\gamma$ and $\beta$) per activation that allow the $\hat{x}_i$ to have any mean and standard deviation. The reason for that is normalizing $x_i$ can reduce its expressive power. This new parameterization has better learning dynamics: in the old parameterization the mean of $x_i$ was determined by a complicated interaction between parameters of all preceding layers - so small changes to the network parameters amplify as the network becomes deeper. In the new parameterization the mean of $\hat{x}_i$ is determined by $\beta$ that we learn along with $\gamma$ during training. Thus, Batch Normalization stabilizes learning. 
As a result, Batch Normalization enables faster training by using much higher learning rates and alleviates the problem of bad initialization. BN also makes it possible to use saturating non-linearities by preventing the network from getting stuck in saturation modes. In summary, Batch Normalization is a differentiable transform that introduces normalized activations into the network. In practice, a BN layer can be inserted immediately after a fully connected layer. 
A: [1] addresses the question:
First, weights shouldn't be set to zeros in order to break the symmetry when backprogragating:

Biases can generally be initialized to zero but weights need to be initialized carefully to break the symmetry between hidden units of the same layer. Because different output units receive different gradient signals, this symmetry breaking issue does not concern the output weights (into the output units), which can therefore also be set to zero.

Some initialization strategies:


*

*[2] and [3] recommend scaling by the inverse of the square root of the fan-in

*Glorot and Bengio (2010) and the Deep Learning Tutorials use a combination of the fan-in and fan-out:


*

*for hyperbolic tangent units: sample a Uniform(-r, r) with $r=\sqrt{\frac{6}{\text{fan-in}+\text{fan-out}}}$ (fan-in is the number of inputs of the unit).

*for sigmoid units : sample a Uniform(-r, r) with $r=4 \sqrt{\frac{6}{\text{fan-in}+\text{fan-out}}}$ (fan-in is the number of inputs of the unit).


*in the case of RBMs, a zero-mean Gaussian with a small standard deviation around 0.1 or 0.01 works well (Hinton, 2010) to initialize the weights.

*Orthogonal random matrix initialization, i.e. W = np.random.randn(ndim, ndim); u, s, v = np.linalg.svd(W) then use u as your initialization matrix.


Also, unsupervised pre-training may help in some situations:

An important choice is whether one should use
  unsupervised pre-training (and which unsupervised
  feature learning algorithm to use) in order
  to initialize parameters. In most settings
  we have found unsupervised pre-training to help
  and very rarely to hurt, but of course that
  implies additional training time and additional
  hyper-parameters.

Some ANN libraries also have some interesting lists, e.g. Lasagne:
Constant([val]) Initialize weights with constant value.
Normal([std, mean]) Sample initial weights from the Gaussian distribution.
Uniform([range, std, mean]) Sample initial weights from the uniform distribution.
Glorot(initializer[, gain, c01b])   Glorot weight initialization.
GlorotNormal([gain, c01b])  Glorot with weights sampled from the Normal distribution.
GlorotUniform([gain, c01b]) Glorot with weights sampled from the Uniform distribution.
He(initializer[, gain, c01b])   He weight initialization.
HeNormal([gain, c01b])  He initializer with weights sampled from the Normal distribution.
HeUniform([gain, c01b]) He initializer with weights sampled from the Uniform distribution.
Orthogonal([gain])  Intialize weights as Orthogonal matrix.
Sparse([sparsity, std]) Initialize weights as sparse matrix.


[1] Bengio, Yoshua. "Practical recommendations for gradient-based training of deep architectures." Neural Networks: Tricks of the Trade. Springer Berlin Heidelberg, 2012. 437-478.
[2] LeCun, Y., Bottou, L., Orr, G. B., and Muller, K. (1998a). Efficient backprop. In Neural Networks, Tricks of the Trade.
[3] Glorot, Xavier, and Yoshua Bengio. "Understanding the difficulty of training deep feedforward neural networks." International conference on artificial intelligence and statistics. 2010.
A: The following explanation is taken from the book: Neural Networks for Pattern Recognition by Christopher Bishop. Great book!
Assume you have previously whitened the inputs to the input units, i.e. $$<x_{i}> = 0$$ and $$<x_{i}^{2}> = 1$$
The question is: how to best choose the weights?. The idea is to pick values of the weights at random following a distribution which helps the optimization process to converge to a meaningful solution.
You have for the activation of the units in the first layer, $$y = g(a) $$ where $$ a = 
\sum_{i=0}^{d}w_{i}x_{i}$$. Now, since you choose the weights independently from the inputs,
$$<a> = \sum_{i=0}^{d}<w_{i}x_{i}> = \sum_{i=0}^{d}<w_{i}><x_{i}> = 0$$ and $$ <a^2> = \left<\left(\sum_{i=0}^{d}w_{i}x_{i}\right) \left(\sum_{i=0}^{d}w_{i}x_{i}\right)\right> = \sum_{i=0}^{d}<w_{i}^{2}><x_{i}^{2}> = \sigma^{2}d $$ where sigma is the variance of the distribution of weights. To derive this result you need to recall that weights are initialized independently from each other, i.e. $$<w_{i}w_{j}> = \delta_{ij}$$
A: Well just as an update, Delving Deep into Rectifiers: Surpassing Human-Level Performance n ImageNet Classification by He et al introduced an initialization specifically with initialization w = U([0,n]) * sqrt(2.0/n) where n is the number of inputs of your NN. I have seen this initialization used in many recent works (also with ReLU). They actually show how this starts to reduce the error rate much faster than the (-1/n, 1/n) that you mentioned. For the thorough explanation, see the paper but here's how fast it converges:

