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.
It looks like you are using an ANN for function approximation where your target is a continuous variable(?). If so, the uniformity of input data spread over all the dimensions and their range and scale can have a strong impact on how well your ANN works. It has been noticed that, even if there is dense coverage of datapoints over the space of the inputs, an ANN can still have a problem making a good prediction. (I have actually had to smooth the input data $\rightarrow$ i.e., cubic splines, to improve uniformity of datapoint coverage over the range and scale of inputs, before the MSE stabilized and decreased substantially). Therefore, straightforward use of raw input data may not help the situation. If you were performing classification analysis, you would need to use the softmax function on the output side -- but it looks like you are merely performing function approximation. You might also try a linear function (identity) on the output side, and what happens on the output side, because the distribution of your target $y$-variable $[y=f(x_1,x_2,...,x_p)]$ can affect the MSE as well. You could also consider an RBF network, or SOM-ANN, which will reduce dimensionality of your inputs. Lastly, correlation between input features degrades learning speed, since an ANN will waste time learning the correlation between inputs. This is why many use PCA on inputs first, and then input e.g. the 10 PCs associated with the greatest eigenvalues -- effectively decorrelating features so they are orthogonal.
Update(4/27/2016):
A way to more evenly distribute randomly sampled points for inputs to a function approximation problem using an ANN is to employ Latin Hypercube Sampling (LHS) from $\{x,y\}$ in order to predict $\hat{z}$. To begin, split up the range of $x$ and $y$ into $M=10$ uniformly spaced, non-overlapping bins -- the result is a $10\times 10$ square grid with 100 elements (cells) -- call this a ``range grid.'' Next, sample one of the 100 cells, and from this cell draw a random pair of $\{x,y\}$ values from within the range of the bin walls for $x$ and $y$ (of the selected cell), and then block that row and column out from further selection. Next, draw another random element from a row and column that hasn't been sampled from yet, and draw a pair of random $\{x,y\}$ from within that cell. Continue until all rows and column have been selected once. The 10 samples of $\{x,y\}$ will provide pairs of points with no overlap, which is a good way to feed $\{x,y\}$ to an ANN for an $\{x,y,z\}$ problem, or multiple feature problem $\{x_1,z_2,\ldots,x_p\}$.
If you want 100 pairs of $\{x,y\}$, you can start with the combination $\{1,2,3,4,5,6,7,8,9,10\}$. Next, identify 10 permutations for this combination to generate a $10 \times 10$ ``row'' matrix $\mathbf{R}$:
$\{3,2,10,4,1,5,7,9,8,6\}$, $\{2,4,3,1,5,6,10,8,9,7\}$,...,$\{9,1,2,10,5,6,7,8,3,4\}$.
which will give 100 integer values for sampling rows.
Next, generate a $10 \times 10$ ``column'' matrix $\mathbf{C}$ using another set of 10 different permutations:
$\{5,8,10,4,9,3,7,1,2,6\}$, $\{3,7,4,1,5,6,10,8,9,2\}$,...,$\{6,9,2,7,5,1,10,8,3,4\}$
which will provide 100 integers for sampling columns.
The first random draw using the above matrices would be from row 3 and col 5 in the original 10 by 10 grid ``range grid'' of 100 bins for $x$ and $y$. This is another form of LHS.
If you need more than 100 $\{x,y\}$ pairs, then just increase the number of permutations used, and don't be stingy as there are 10! permutations.
Best Answer
There has been quite a lot of theoretical work on Neural Network initialization in the last 5 years, which apparently hasn't still propagated to the wider Deep Learning community. While it's true that there isn't still an initialization which works for all architectures and for all activation functions (and most likely there will never be, from what we have understood so far about the dynamics of deep neural networks) , in practice this isn't a huge limit because most users use two or three activation functions (ReLU for CNNs, tanh and sigmoid for LSTMs, AKA the only RNNs used by most people) and two or three architectures (i.e., ResNets for image classification and LSTMs for sequence prediction or time series forecasting). For these things we do have some powerful results. Granted, they won't work for the Universal Transformer just presented at ICML 2018, but frankly, right now there are much more people trying to apply "standard" architectures such as ResNets to interesting business problems, who need better initializations than the Xavier one, than there are people inventing bleeding-edge architectures, who will always have to use the classic, infallible initialization strategy: elbow grease, also known as "long and boring computational experiments, supported by careful book-keeping".
Deep Linear Networks
One of the first main results is shown in the seminal work "Exact solutions to the nonlinear dynamics of learning in deep linear neural networks" of Andrew Saxe et al., 2014, on optimal activation for deep linear network, which used results from random matrix theory and in particular from free probability theory (AKA the theory of probability for noncommutative random variables, such as random matrices, precisely), to show that using random orthogonal weights provides far superior results than the usual scaled Normal initializations. The theoretical concept behind this result is that of dynamical isometry, i.e., having the product of Jacobians associated with error signal backpropagation act as a near isometry, up to some overall global $O(1)$ scaling, on a subspace (of the weight space) of as high a dimension as possible. This is equivalent to having as many singular values of the product of Jacobians as possible within a small range around an $O(1)$ constant, and is closely related to the notion of restricted isometry in compressed sensing and random projections.
Deep Nonlinear Networks
However, the results from Saxe don't easily translate to "real" NNs, with nonlinear activation functions. One can try and apply the random orthogonal matrix initialization, and there are "regimes", such as "the edge of chaos", where it will indeed work incredibly well for an architecture otherwise horribly hard to train (a fully connected neural network with 100 layers, 1000 units per layer and tanh activations (!!!)). But outside of this regime, there's no guarantee: results can even be worse than Xe initialization, so you're back to "try and see".
However, 3 years after, two different papers appeared which extended Saxe's work to nonlinear neural networks by trying to achieve dynamical isometry for nonlinear networks. One is the work by Pennington et al, "Resurrecting the sigmoid in deep learning through dynamical isometry: theory and practice", 31st Conference on Neural Information Processing Systems (NIPS 2017), Long Beach, CA, USA. This paper fails to find a useful initialization for ReLU networks, while it finds a very good one for tanh networks. Since you're interested in the latter, I'll skip this one, just listing it here for completeness.
The other is Balduzzi et al.'s well-known Shattered Gradients paper (Balduzzi et al., "The Shattered Gradients Problem: if resnets are the answer, then what is the question?", arXiv:1702.08591v2, 2018) which proposes to Saxe's random orthogonal weights with the so-called Looks-Linear (LL) activation. The approach is as follows:
substitute all ReLU with CReLU, i.e., concatenated rectifiers: basically, these are very similar to a ReLU unit, but instead than
$$ \sigma(x) = \max(0,x)$$
we have $$ \boldsymbol{\rho}(x) = (\max(0,x), \max(0,-x)) $$
(note that $\boldsymbol{\rho}(x):\mathbb{R}\to\mathbb{R}^2$). While you're here, you would probably like to reduce the number of units by 2, to keep the same number of parameters as before.
now, each layer has a weight matrix $W_l$, which has twice as many elements as before, if you didn't reduce the number of units in each layer to achieve parameter parity with your initial architecture. In any case, it's a matrix with an odd number of columns, because of the CReLU activation function, thus $W_l=[W_{1l}, W_{2l}]$ where $W_{1l}$ and $W_{2l}$ have the same shape. Now, for each layer sample a matrix $W'_l$ with orthogonal columns, of the same shape as $W_{1l}$, and initialize $W_l$ as $W^0_l=[W'_l, -W'_l]$. Clearly, at initialization you now have a linear network because
$$ \boldsymbol{\rho}(W^0_l \mathbf{x})=W'_l\sigma(\mathbf{x})-W'_l\sigma(-\mathbf{x})=W'_l\mathbf{x} $$
which is why we call this initalization LL (looks-linear).
The LL-init can be "extended" easily to CNNs (see the cited paper for details). It does have the disadvantage of forcing you to change your architecture, though it's admittedly a simple change.
Finally, the most impressive result, based on a mix of theory and experiments is the Delta-Orthogonal initialization of Xiao et al.,"Dynamical Isometry and a Mean Field Theory of CNNs: How to Train 10,000-Layer Vanilla Convolutional Neural Networks", arXiv:1806.05393v2, 2018. This initializer has obtained amazing results, such as allowing successful training of a 10000 layers vanilla CNN with tanh activations, without nearly any regularization techinque (no dropout, no residual connections, no Batch Norm, no weight decay and no learning rate decay: the network relies only on SGD with momentum for regularization). The initialization is even included in Tensorflow as the ConvolutionOrthogonal initializer in the suite of initialization operators of Tensorflow.
Sadly, this initializer only works its magic on tanh CNNs, while it's not guaranteed to deliver its amazing speedups for ReLU CNNs. And tanh CNNs do suck at object classification: the current SOTA for a tanh CNN on CIFAR-10 is a test error of more than 10%, while with modified ReLU ResNets we go below 3%.