Feature Scaling and the Conditioning of Gradient Descent

Two input columns describe the same physical quantity, one in metres and one in millimetres. Nothing about the underlying problem changed, but a model trained on the second pair of columns will take a different path to its solution, or fail to reach one at all. The reason is not numerical precision. It is that rescaling a single input multiplies one eigenvalue of the loss surface’s Hessian by a million and leaves the others alone, and first-order optimisation is governed by the ratio between the largest and smallest of those eigenvalues.

The quadratic case, where the argument is exact

Take a loss L(w) = (1/2) w^T H w - b^T w with H symmetric positive definite. Gradient descent with step size eta iterates w_{k+1} = w_k - eta (H w_k - b). In the eigenbasis of H the update decouples into one scalar recursion per eigenvalue lambda_i: the error component along eigenvector i is multiplied by 1 - eta * lambda_i each step.

Two constraints follow immediately. Stability requires eta < 2 / lambda_max, and speed requires 1 - eta * lambda_min to be far from one in magnitude, which pushes eta toward 2 / (lambda_max + lambda_min). At that optimal fixed step the contraction factor is

contraction per step: rho = (kappa - 1) / (kappa + 1), where the condition number is kappa = lambda_max / lambda_min

Iterations to reduce the error by a factor of e are therefore of order kappa. A condition number of 10^6 means a million steps where a well-conditioned problem would need a handful. This is the whole content of the scaling argument for smooth convex problems; everything else is an application of it.

Now suppose the first feature is measured in units a thousand times smaller than the second. The Hessian entry for that feature is scaled by 10^6 relative to the other diagonal entry, so kappa acquires a factor of roughly 10^6 for two otherwise symmetric features. Dividing each column by its standard deviation removes exactly that factor from the diagonal. It is a diagonal preconditioner applied to the data rather than to the gradient.

What standardisation does and does not fix

Standardisation subtracts the column mean and divides by the column standard deviation, giving each column unit sample variance. On the diagonal of the Hessian this is close to optimal: every feature contributes curvature of comparable magnitude. Off the diagonal, nothing changes in relative terms. Two perfectly correlated features remain perfectly correlated after standardisation, and the Hessian still has a near-zero eigenvalue along the direction that increases both coefficients together while keeping their difference fixed. Standardisation fixes the units, not the redundancy. Removing that requires decorrelation, which is what whitening does: transform the vector by multiplying with Sigma^(-1/2), using the sample covariance matrix Sigma. Whitening makes the Hessian closer to a multiple of the identity, at the cost of estimating d(d+1)/2 covariance parameters from n rows and inverting a d by d matrix, an O(d^3) operation that is unstable when n is not comfortably larger than d. On wide data the covariance estimate is singular and whitening amplifies noise directions instead of suppressing them.

Min-max scaling maps each column affinely onto the interval from 0 to 1 using the observed minimum and maximum. It bounds the input range, which some models require, but its sensitivity to a single value is total: one outlier at ten times the bulk of the data compresses every other value into the first tenth of the interval, and the effective resolution of the column collapses. The breakdown point of both min-max and standardisation is zero, since one observation can move the statistic arbitrarily. Median and interquartile range scaling has a breakdown point of about 25 percent, which is why it is the safer default on columns whose tails are not understood.

The scaler is part of the model, not a preprocessing step

Scaling is usually described as a transformation applied before training. It is better understood as a change to the hypothesis class. Ridge regression penalises the squared Euclidean norm of the coefficient vector, and the meaning of that penalty depends entirely on the units of the inputs: a coefficient on a column measured in millimetres must be a thousand times larger than the same effect measured in metres, so the penalty on it is a million times larger. Two analysts who scale differently and then tune the penalty to the same value are fitting different priors. This is not a defect of regularisation; it is a reason to fix the scaling convention before tuning, and to apply the same convention at serving time.

Distance-based methods inherit the same dependence more directly. In k-nearest neighbours, k-means, and kernel machines with an isotropic kernel, the distance between two points is dominated by whichever column has the largest variance. A column measured in bytes will swamp a column measured in seconds, and the resulting neighbourhoods reflect the units rather than the problem. Principal components analysis has the same property: without standardisation the first component is approximately the column with the largest variance, and its loadings describe measurement scale.

import numpy as np

def standardise(train, other):
    mu = train.mean(axis=0)
    sigma = train.std(axis=0, ddof=0)
    # A constant column has zero spread; substituting 1 keeps the transform total
    sigma = np.where(sigma == 0.0, 1.0, sigma)
    return (train - mu) / sigma, (other - mu) / sigma

X_train, X_test = standardise(X_train, X_test)

The function takes the training matrix and applies the resulting statistics to the held-out matrix. That signature is the point: the mean and standard deviation are estimated parameters of the model, and estimating them on data that includes the test rows is a leak, discussed in how a feature becomes a leak. The effect is usually small for a scaler, because a mean over thousands of rows barely moves when a few hundred are added, but it is unbounded in principle and it is the same mistake, in kind, as fitting a target encoder on the full table.

Where monotone invariance makes scaling irrelevant

Decision trees partition the axis of each feature at thresholds and are invariant under any strictly increasing transformation of a single column. A gradient-boosted model trained on raw counts and on their logarithms will produce the same predictions up to the floating-point representation of the split points, because the ordering of values within a column is all the split search uses. Scaling therefore has no effect on tree ensembles, on rank-based statistics, or on any model whose decision depends only on per-column order.

The invariance is narrower than it looks. It holds per column and under monotone maps only. A linear combination of two columns is not preserved by scaling them differently, so any model that forms such combinations, including a tree ensemble’s own leaf values when features are correlated, is not invariant in the way the axis-aligned split argument suggests. And the invariance does not extend to what happens after the raw model: if feature importances, partial dependence, or coefficient magnitudes are compared across columns, the comparison is in the units of the input unless the input was standardised first.

Neural networks sit between the two cases. Input standardisation is not required for correctness, since the first layer can in principle learn any rescaling, but the initialisation of that layer is scale-dependent: weights drawn from a fixed distribution produce pre-activations whose variance is proportional to the sum of squared input magnitudes, so a column with large variance saturates activations or leaves them near zero. Normalising the input makes the initial forward pass well-conditioned, which is why the practice appears in the standard backpropagation guidance long before batch normalisation existed. Once training proceeds, adaptive methods such as Adam are invariant to a per-parameter rescaling of the gradient in the long run, but not to the interaction between the input scale and weight decay, which is again a penalty on parameters expressed in input units.

Degenerate columns and the cost of centring

A column that is constant across the training set has zero standard deviation. Dividing by it yields a non-finite value for every row; substituting one, as implementations conventionally do, leaves a column of exact zeros after centring, which contributes nothing and can be dropped. The same substitution is what makes a scaler safe to fit on a batch that happens to contain a single distinct value.

Centring a sparse matrix destroys its sparsity. A matrix with one non-zero per row becomes dense after subtracting the mean, because every zero entry becomes a non-zero mean offset. Implementations therefore offer scalers that divide by a per-column maximum absolute value without centring, preserving the zeros, at the price of leaving the column mean wherever it was. For text-like data with hundreds of thousands of columns this is the difference between a few megabytes and tens of gigabytes, and it is one of the few cases where the memory argument dominates the conditioning argument. The related problem of collapsing an unbounded token vocabulary into a fixed-width vector is treated in feature hashing.

None of this changes the underlying statement: the scaler is a fitted parameter, its statistics belong to the training partition, and the choice of scaler determines which directions of the input space the model treats as equivalent. A model trained on unscaled data is not wrong; it is a model whose implicit prior over coefficients, and whose optimiser trajectory, were chosen by the units the data happened to arrive in.