Are there any theoretically rigorous justification for why scaling or normalizing data should improve statistical performance? I was wondering if there was a rigorous theoretical justification for why normalizing (or scaling) the feature vectors might make sense?
The reason that I am asking this is because, at least in the case when one decides to center the data (i.e. subtract off the mean), it can be viewed as adding an offset that does not have penalization. This justification makes sense in a rather constrained way (I am not aware of any generalizations), but this justification makes sense in the context of squared loss, 2-norm regularization and linear kernel $K(x,x') = x^Tx' = \langle x, x' \rangle$.
Has there been any similar results in terms of feature scaling or normalizations (features have norm 1 or to normalizing with respect to z-scores $Z = \frac{X - \mu}{\sigma}$)?
I do know a rather hand wavy (but good) justifications for why it might make sense. These explanation can be found in this video. Mainly what I took out of it was that scaling helps so that the features are in the same range or order of magnitude. This helps, specially in the context of linear regression and gradient descent because otherwise, gradient descent might jump around the minimum in an awkward way if the features are at a different scale. Centering helps because even if you scale the features, if they are not around "the same" mean, then the scaling might not actually make sense (because of this offset). The reason why this should improve statistical performance should be clear, some data just isn't generated from a zero offset curve, so its hard to generalize if the offset is not present. These explanations have already been discussed at here. 
Apart from that, I am not sure why these methods should work in general or why they should lead to good generalizations.

For the curious minds I will give a rough outline to hows centering the data is equivalent to not penalizing the offset and doing standard empirical risk minimization. Let the empirical risk $\mathcal{E}_X(w,b)$ (for data set $X$) be defined as follows:
$$ \mathcal{E}_X(w,b) = \frac{1}{n} \sum^n_{i=1} ( w^T x_i + b - y_i)^2 + \lambda \| w \|^2$$
what we want to show is that the solution to empirical risk minimization of the following two optimization problems are equivalent
$$ \min_{w,b} \mathcal{E}_X(w,b) = \min_{w} \mathcal{E}_{X^c}(w)$$
where $X^c$ is the centered training set.
Since the above problem is convex, the global minimum is achieved when the gradients are equal to zero. In particular it can be shown:
$$ \frac{\partial \mathcal{E}_X(w,b) }{\partial b} = 0 \iff b^* = \bar{y} - w^T\bar{x}$$
where $\bar{x} = \frac{1}{n} \sum^n_{i=1} x_i$ and $\bar{y} = \frac{1}{n} \sum^n_{i=1} y_i$ and $b^*$ is the optimal (not penalized) offset. Now, its not hard to show that:
$$ \min_{w,b} \mathcal{E}_X(w,b) = \min_{w} \mathcal{E}_X(w,b^*) = \min_{w} \mathcal{E}_{X^c}(w)$$
which completes our proof.
 A: In the context of kernel methods, this is quite easy to see. Any given data matrix $\mathbf{X}$ corresponds to a kernel matrix $\mathbf{K}$, which in turns corresponds to a certain solution of the training problem which, due to convexity, is unique and guaranteed to be the global optimum.
However, if we change $\mathbf{X}$ into $\mathbf{X}'$, via some transformation, then we get a different kernel matrix $\mathbf{K}'$ and evidently a different solution. I explicitly mention this, because it is important to realize that the solution changes when we change $\mathbf{X}$, for instance by scaling. The solution is still found by solving a convex problem, and hence is still unique and the global optimum of its corresponding training problem (which is different for $\mathbf{K}$ and $\mathbf{K}'$).
The reason we do scaling, then, has nothing to do with solving the optimization problem but rather with defining it. For simplicity, lets assume we use the standard linear kernel:
 $$\kappa(\mathbf{u},\mathbf{v}) = \mathbf{u}^T\mathbf{v}$$
If the features are on different scales, say, the first dimension is in $[0, 10^{99}]$ and the second is in $[0, 1]$, then it is easy to see that in all these kernel evaluations the first feature will completely dominate in the resulting distance estimates (that is, entries in the kernel matrix are almost exclusively based on the first dimension).
When you start building a model, you typically want to give each feature a similar contribution in the model, and for kernel methods that implies that they should be on the same scale. Scaling in no way guarantees better performing models, but it is usually the best "prior" you have. For example assume that the first feature in the contrived example above is informative while the second is not, in this case leaving the first feature on a far larger scale is in fact better.
Another example is available here.
