Scaling Predictors for Stable Optimization

Why Predictor Units Matter for Optimization

A GLMM can be statistically well specified yet numerically awkward to optimize. If one continuous predictor is recorded in units near 100,000 and another in units near 0.01, their coefficient coordinates can have very different scales. That mismatch can make the objective function steep in some directions and comparatively flat in others. An optimizer may then take inefficient steps or have difficulty deciding that it has reached a stable solution.

Rescaling changes the numerical coordinates used by the optimizer. It does not, by itself, add information or alter the scientific model. The key is to treat the scaling as an exact reparameterization and to interpret diagnostics in the coordinates in which they were calculated.

Definition: For a continuous predictor \(x\), choose a positive scale constant \(c\) and define \(z=x/c\). In a linear predictor containing \(\beta x\), write the same contribution as \(\gamma z\), where \(\gamma=c\beta\). Thus \(\beta\) is the change in linear predictor per one unit of \(x\), while \(\gamma\) is the change per \(c\) units of \(x\). Scaling changes the predictor’s unit, not its origin.

Equivalent GLMMs and Coordinate Changes

For example, a random-intercept logistic GLMM may be written as

Y_{ij}\mid u_j \sim \operatorname{Bernoulli}(p_{ij}), \qquad \operatorname{logit}(p_{ij})=\beta_0+\beta_1x_{ij}+u_j,\qquad u_j\sim N(0,\tau^2).

Here \(i\) indexes observations within group \(j\), and the conditional responses are independent given the random effects. Replacing \(x\) by \(z=x/c\) and \(\beta_1\) by \(\gamma_1=c\beta_1\) leaves every conditional linear predictor unchanged because \(\beta_1x=\gamma_1z\). Consequently, the fitted probabilities and the likelihood are unchanged when all other parts of the model are held fixed.

The same principle applies to other GLMM likelihoods. If there is a random slope for \(x\), its coefficient and covariance parameters must also be transformed consistently: the random-slope contribution must remain the same on the linear-predictor scale. Merely rescaling the fixed effect while leaving a random slope’s distribution unchanged generally changes the model. Likewise, if the fitting objective includes penalties or priors, they must be transformed appropriately for the objective to represent the same model.

The optimizer sees an objective as a function of parameter coordinates. Let \(\theta\) denote the original coordinates, including the coefficient for \(x\), and let \(\eta\) denote the coordinates after scaling. For a linear coordinate transformation \(\theta=A\eta\), define \(f_{\mathrm{new}}(\eta)=f_{\mathrm{old}}(A\eta)\), where \(f\) is the objective being minimized, such as a negative log likelihood. At corresponding parameter values, the chain rule gives

g_{\mathrm{new}}=A^{T}g_{\mathrm{old}}, \qquad H_{\mathrm{new}}=A^{T}H_{\mathrm{old}}A,

where \(g\) is the gradient and \(H\) is the Hessian matrix. The Hessian formula follows by differentiating the transformed gradient once more; it holds because \(A\) is constant. If only \(x\) is replaced by \(x/c\), the corresponding original coefficient is \(\beta=\gamma/c\). Thus the entry of \(A\) for that coefficient is \(1/c\), while unchanged parameters have diagonal entries of 1.

Justification: Since \(\theta=A\eta\), a small change \(d\eta\) induces \(d\theta=A\,d\eta\). The first-order change in the objective is \(df=g_{\mathrm{old}}^T d\theta=g_{\mathrm{old}}^T A\,d\eta\), so the new gradient is \(A^Tg_{\mathrm{old}}\). Differentiating this gradient with respect to \(\eta\), with \(A\) constant, gives \(H_{\mathrm{new}}=A^TH_{\mathrm{old}}A\). This also shows why the Hessian’s eigenvalues cannot usually be obtained by simply rescaling the old eigenvalues: the transformed matrix depends on the full Hessian, including its off-diagonal entries.

Worked Calculations

Worked Example: Changing a Coefficient’s Units

Suppose a hypothetical logistic random-intercept model has the fixed-effect contribution \(-2.0+0.0008x\), where \(x\) is a biomarker measured in units that typically reach several hundred. The random intercept is \(u_j\), so the full linear predictor is \(-2.0+0.0008x+u_j\). Rescale using \(z=x/100\).

Because \(x=100z\), the same contribution is

-2.0+0.0008x+u_j =-2.0+(100\times0.0008)z+u_j =-2.0+0.08z+u_j.

The coefficient is now \(0.08\) per 100 units of the biomarker rather than \(0.0008\) per one unit. For an observation with \(x=250\), \(z=2.5\), and the predictor contribution from the biomarker is \(0.0008(250)=0.20\) in the original units and \(0.08(2.5)=0.20\) after rescaling. The fitted linear predictor and probability are therefore identical, provided the rest of the model is unchanged. This is a change in numerical scale, not evidence that the effect became larger.

Worked Example: Transforming the Hessian Correctly

Consider a hypothetical local Hessian for two parameters: \(\beta\), the coefficient of a large-unit predictor, and \(\alpha\), another unchanged parameter. Suppose the Hessian of the objective at the fitted point is

H_{\mathrm{old}}= \begin{pmatrix} 80 & 0.4\\ 0.4 & 0.02 \end{pmatrix}.

This matrix is positive definite because its first diagonal entry is positive and its determinant is \(80(0.02)-(0.4)^2=1.44>0\). Now use \(z=x/100\). The new coefficient is \(\gamma=100\beta\), so \((\beta,\alpha)^T=A(\gamma,\alpha)^T\), with \(A=\operatorname{diag}(0.01,1)\). Apply the transformation to the full matrix:

H_{\mathrm{new}}=A^TH_{\mathrm{old}}A = \begin{pmatrix} 0.01 & 0\\ 0 & 1 \end{pmatrix} \begin{pmatrix} 80 & 0.4\\ 0.4 & 0.02 \end{pmatrix} \begin{pmatrix} 0.01 & 0\\ 0 & 1 \end{pmatrix} = \begin{pmatrix} 0.008 & 0.004\\ 0.004 & 0.02 \end{pmatrix}.

For the original matrix, the trace is \(80.02\), the determinant is \(1.44\), and the eigenvalues are approximately \(80.002\) and \(0.01800\). Its condition number, the ratio of the largest to smallest eigenvalue for this positive-definite matrix, is approximately \(80.002/0.01800=4{,}445\).

For the transformed matrix, the trace is \(0.028\) and the determinant is \(0.008(0.02)-(0.004)^2=0.000144\). Its eigenvalues are

\lambda_{\pm} =\frac{0.028\pm\sqrt{(0.008-0.02)^2+4(0.004)^2}}{2} =\frac{0.028\pm\sqrt{0.000208}}{2},

which gives approximately \(0.02121\) and \(0.006789\). Their product is \(0.000144\), matching the determinant, and their ratio is approximately \(3.124\). In this hypothetical example the coordinate change substantially improves the Hessian condition number. The improvement is not inferable from the original eigenvalues alone: it follows from transforming the actual Hessian, including its off-diagonal term.

Worked Example: Reading Gradient Diagnostics in Both Scales

Suppose, at corresponding trial parameter values, the original-coordinate gradient for the same scaled coefficient and an unchanged parameter is \(g_{\mathrm{old}}=(0.40,\,0.006)^T\). With \(c=100\), the coefficient-coordinate transformation is \(A=\operatorname{diag}(0.01,1)\), so

g_{\mathrm{new}}=A^Tg_{\mathrm{old}} =\begin{pmatrix}0.01&0\\0&1\end{pmatrix} \begin{pmatrix}0.40\\0.006\end{pmatrix} =\begin{pmatrix}0.004\\0.006\end{pmatrix}.

The raw Euclidean gradient norms are \(\sqrt{0.40^2+0.006^2}\approx0.4000\) and \(\sqrt{0.004^2+0.006^2}\approx0.007211\). Their difference does not mean the rescaled fit is closer to a minimum: a raw gradient norm depends on parameter units. Compare gradients using a stated scaling convention and the same objective, and inspect them together with curvature and the optimizer’s convergence information. At a corresponding interior stationary point, a nonsingular coordinate transformation preserves the condition \(g=0\); scaling changes numerical component sizes, not whether the exact gradient is zero.

Using Scaling as a Diagnostic

A useful comparison holds the statistical model fixed and changes only the units of one or more continuous predictors. Fit the model in each parameterization, translate coefficients back to a common scale, and check whether the objective values, fitted linear predictors, and predictions agree to numerical tolerance. Then compare the gradient and Hessian diagnostics in each fit’s own coordinates. This helps distinguish a numerical improvement from a substantive change or an unintended model change.

Scaling is most informative when it addresses a known mismatch in predictor units. A divisor might be 100, 1,000, or a scientifically interpretable unit such as 10 years, chosen so the rescaled values are of a manageable order. A standard-deviation scale is another possible choice, but scaling is not the same operation as centering: here \(z=x/c\), with no subtraction of a mean.

For an interior minimum, a positive-definite Hessian is expected locally. Its eigenvalues describe curvature in the chosen coordinates, and a large condition number signals that some directions are much less curved than others. Since a change of units changes those coordinates, neither the eigenvalues nor the condition number is an invariant property of the statistical model. A smaller condition number after scaling can indicate more favorable numerical geometry, but it does not prove that all parameters are well identified.

A scaling check is less reassuring when the two fits do not yield equivalent objective values or predictions, when a transformed coefficient is inconsistent with its original-scale counterpart, or when a gradient or Hessian diagnostic remains problematic. These discrepancies can arise from a coding error, a changed penalty or prior, incomplete optimization, or an ill-conditioned problem that scaling did not resolve. A flat likelihood direction, weak design information, an invalid Hessian, or a variance estimate at a boundary is not automatically repaired by changing predictor units.

Common Mistakes and Alternatives

Comparing raw gradients across units. A smaller raw gradient norm after rescaling is not, by itself, evidence of a better fit. Gradient components change under the coordinate transformation; interpret them using the stated scaling convention and other convergence diagnostics.
Guessing transformed eigenvalues from the old eigenvalues. The transformation \(H_{\mathrm{new}}=A^TH_{\mathrm{old}}A\) acts on the matrix, not independently on its eigenvalues. Compute the transformed matrix and its eigenvalues, as in the example.
Treating rescaling as a change in the scientific effect. The coefficient changes because its unit changes. Translate it back before comparing effect estimates.
Assuming scaling solves every convergence problem. It can improve numerical conditioning, but it cannot create information or guarantee a unique, well-behaved optimum.

Rescaling predictors is a direct, model-preserving way to improve numerical coordinates. Other responses to a difficult fit include examining the model specification, design information, and optimizer behavior. These are not substitutes for checking equivalence after rescaling: the first question is whether the same model has been represented and fitted. Scaling is most useful as a controlled diagnostic, not as a reason to dismiss a persistent convergence warning.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

If \(z=x/50\) and the original coefficient is \(\beta=0.012\) per unit of \(x\), what is the coefficient per unit of \(z\)?
For a two-parameter rescaling with \(A=\operatorname{diag}(0.1,1)\), what is the transformation rule for the Hessian?
Why can a raw gradient norm decrease after rescaling even if the optimizer has not moved closer to a stationary point?
Does a smaller Hessian condition number after rescaling establish that the model is well identified? Explain.

Solutions. (1) Since \(x=50z\), \(\gamma=50\beta=50(0.012)=0.600\) per unit of \(z\). (2) Use \(H_{\mathrm{new}}=A^TH_{\mathrm{old}}A\); both rows and columns associated with the first coordinate are multiplied by \(0.1\), so its diagonal curvature is multiplied by \(0.01\). (3) Gradient components transform as \(g_{\mathrm{new}}=A^Tg_{\mathrm{old}}\), so their numerical magnitudes, and hence the raw norm, depend on the units of the parameter coordinates. (4) No. The condition number depends on coordinates; improved conditioning can help optimization, but does not establish that the data identify every parameter or rule out flat likelihood directions.