Why should weights of Neural Networks be initialized to random numbers? Premise

Ok, I know that this question was asked before on ai.SE, on stats.SE and also on SO. So I did my homework in checking before posting my question, but none of them has an answer that fully satisfies me.

Summary of my knowledge up to now

Suppose you have a layer that is fully connected, and that each neuron performs an operation like

a = g(w^T * x + b)


were a is the output of the neuron, x the input, g our generic activation function, and finally w and b our parameters. 

If both w and b are initialized with all elements equal to each other, then a is equal for each unit of that layer. 

This means that we have symmetry, thus at each iteration of whichever algorithm we choose to update our parameters, they will update in the same way, thus there is no need for multiple units since they all behave as a single one. 

In order to break the symmetry we could randomly initialize the matrix w and initialize b to zero (this is the setup that I've seen more often). 
This way a is different for each unit so that all neurons behave differently.

Of course randomly initializing both w and b would be also okay even if not necessary.

The actual question

Is randomly initializing w the only choice? 
Could we randomly initialize b instead of w in order to break the symmetry?
Is the answer dependent on the choice of the activation function and/or the cost function?

My thinking is that we could break the symmetry by randomly initializing b, since in this way a would be different for each unit and, since in the backward propagation the derivatives of both w and b depend on a(at least this should be true for all the activation functions that I have seen so far), each unit would behave differently. Obviously, this is only a thought, and I'm not sure that is absolutely true.  

Remark

Note that this is a theoretical question, not a practical one. 
This means that I'm not really interested in answers that involves performance, but you're welcome to include them as corollary to the real answer.
Randomising just b sort of works, but setting w to all zero causes severe problems with vanishing gradients, especially at the start of learning. 

Using backpropagation, the gradient at the outputs of a layer L involves a sum multiplying the gradient of the inputs to layer L+1 by the weights (and not the biases) between the layers. This will be zero if the weights are all zero.

A gradient of zero at L's output will further cause all earlier layers(L-1, L-2 etc all the way back to layer 1) to receive zero gradients, and thus not update either weights or bias at the update step. So the first time you run an update, it will only affect the last layer. Then the next time, it will affect the two layers closest to the output (but only marginally at the penultimate layer) and so on.

A related issue is that with weights all zero, or all the same, maps all inputs, no matter how they vary, onto the same output. This also can adversely affect the gradient signal that you are using to drive learning - for a balanced data set you have a good chance of starting learning close to a local minimum in the cost function.

For deep networks especially, to fight vanishing (or exploding) gradients, you should initialise weights from a distribution that has an expected magnitude (after multiplying the inputs) and gradient magnitude that neither vanishes nor explodes. Analysis of values that work best in deep networks is how Xavier/Glorot initialisation were discovered. Without careful initialisation along these lines, deep networks take much longer to learn, or in worst cases never recover from a poor start and fail to learn effectively.

Potentially to avoid these problems you could try to find a good non-zero fixed value for weights, as an alternative to Xavier initialisation, along with a good magnitude/distribution for bias initialisation. These would both vary according to size of the layer and possibly by the activation function. However, I would suspect this could suffer from other issues such sampling bias issues - there are more weights, therefore you get a better fit to desired aggregate behaviour when setting all the weight values randomly than you would for setting biases randomly.