Iteratively Reweighted Least Squares for Huber Regression

From Huber Loss to Weighted Least Squares

Least squares gives a large residual a large squared contribution to the fit. Huber regression changes that contribution: residuals near zero are treated quadratically, while the cost of increasingly extreme residuals grows only linearly. Iteratively reweighted least squares (IRLS) turns this loss into a sequence of familiar weighted least-squares calculations.

The weights are not assigned once and then kept. They are recalculated from the current residuals, so each fit supplies the weights for the next fit. Understanding this cycle is important: a case with a large residual gets less weight, but its weight depends on the current fitted line and on the chosen residual scale.

Definition: For observations \((x_i,y_i)\), let \(x_i\) denote the predictor vector for observation \(i\), including a column for the intercept if the model has one. Write the fitted value as \(x_i^\mathsf{T}\beta\), the residual as \(r_i=y_i-x_i^\mathsf{T}\beta\), and let \(s>0\) be a fixed residual scale. The standardized residual is \(u_i=r_i/s\). For a chosen tuning constant \(c>0\), Huber’s loss is \(\rho_c(u)=u^2/2\) when \(|u|\le c\), and \(\rho_c(u)=c|u|-c^2/2\) when \(|u|>c\).

The corresponding Huber regression estimate minimizes \(\sum_{i=1}^n \rho_c(r_i/s)\) over the coefficient vector \(\beta\). For the IRLS derivation, \(s\) and \(c\) are held fixed while the coefficients are updated. This tutorial concerns the coefficient-fitting steps; how a robust scale might be estimated is a separate issue.

The standard setup is a linear predictor \(x_i^\mathsf{T}\beta\) for each observed response. The fitting calculation itself does not require normally distributed errors. For the weighted least-squares update to have a unique solution, the design matrix must have full column rank; in particular, the data must contain enough independent predictor information to estimate all coefficients. IRLS is an algorithm for finding the Huber fit, not by itself a procedure for statistical inference.

Deriving the Residual Weights

The derivative of Huber’s loss is its score function, \(\psi_c(u)\). Differentiating each part of the loss gives a linear score for residuals inside the cutoff and a capped score outside it:

\psi_c(u)= \begin{cases} u, & |u|\le c,\\ c\,\operatorname{sign}(u), & |u|>c. \end{cases}

For \(u\ne 0\), define the weight as the score divided by the standardized residual. At zero, define the weight to be 1, which agrees with the limiting value. Substituting the two score expressions gives:

w(u)= \begin{cases} 1, & |u|\le c,\\ c/|u|, & |u|>c. \end{cases} \qquad\text{or equivalently}\qquad w(u)=\min\left(1,\frac{c}{|u|}\right),

Thus the weight is 1 for residuals within the cutoff and decreases inversely with the residual magnitude beyond it. In original response units, the cutoff is \(cs\): a residual with \(|r_i|\le cs\) receives weight 1, and one with \(|r_i|>cs\) receives weight \(cs/|r_i|\).

Key result: At a current coefficient estimate, calculate \(u_i=r_i/s\) and set \(w_i=w(u_i)\). Since \(\psi_c(u_i)=w_i u_i\), the Huber estimating equation \(\sum_i x_i\psi_c(u_i)=0\) is equivalent, for fixed \(s\), to the weighted least-squares normal equation \(\sum_i w_i x_i r_i=0\).

To see the connection, differentiate the objective \(\sum_i \rho_c(r_i/s)\) with respect to \(\beta\). Because \(r_i=y_i-x_i^\mathsf{T}\beta\), its derivative with respect to \(\beta\) is \(-x_i\). The resulting gradient is \(-\frac{1}{s}\sum_i x_i\psi_c(u_i)\). Setting it equal to zero yields the Huber estimating equation. Replacing \(\psi_c(u_i)\) by \(w_i r_i/s\) gives the weighted normal equation; the common factor \(1/s\) does not change its zero.

\sum_{i=1}^n w_i x_i\left(y_i-x_i^\mathsf{T}\beta_{\mathrm{new}}\right)=0 \quad\Longleftrightarrow\quad X^\mathsf{T}WX\beta_{\mathrm{new}}=X^\mathsf{T}Wy,

Here \(X\) is the design matrix, \(y\) is the response vector, and \(W\) is the diagonal matrix whose diagonal entries are the current weights \(w_i\). Solving this equation is a weighted least-squares step. The weights are computed using the current coefficients, whereas \(\beta_{\mathrm{new}}\) is the result of the weighted fit.

The IRLS Cycle

IRLS repeatedly applies the same sequence. Begin with coefficients, calculate their residuals, convert those residuals into weights, and solve the weighted least-squares problem. Then use the new coefficients to calculate fresh residuals and weights. Stop when successive coefficient estimates (and, in practice, the objective) have changed by less than a chosen tolerance.

Choose fixed settings and a start.
Specify the tuning constant \(c\), the scale \(s\), and starting coefficients \(\beta^{(0)}\).

Calculate residuals and weights.
At iteration \(t\), calculate \(r_i^{(t)}=y_i-x_i^\mathsf{T}\beta^{(t)}\) and \(w_i^{(t)}=\min(1,cs/|r_i^{(t)}|)\), with weight 1 when the residual is zero.

Fit the weighted model.
Solve \(X^\mathsf{T}W^{(t)}X\beta^{(t+1)}=X^\mathsf{T}W^{(t)}y\).

Check and repeat.
Use \(\beta^{(t+1)}\) to recalculate residuals and weights. Repeat until the coefficients and objective are sufficiently stable.

Worked Example: Converting Residuals Into Huber Weights

Suppose the fixed scale is \(s=2\) response units and the tuning constant is \(c=1.5\). At one iteration, four residuals are \(0.4,\ 2,\ 5,\) and \(-8\) response units. Find their standardized residuals and weights.

The cutoff on the original residual scale is \(cs=1.5(2)=3\) response units. Standardizing gives \(u_i=r_i/2\), so the four values are \(0.2,\ 1,\ 2.5,\) and \(-4\). The first three have absolute standardized residuals at most \(1.5\), and therefore each gets weight 1.

For the residual of \(-8\), the standardized residual is \(-4\), so its weight is \(c/|u_i|=1.5/4=0.375\). In residual units the same calculation is \(cs/|r_i|=3/8=0.375\). The weights, in order, are therefore \(1,\ 1,\ 1,\) and \(0.375\). The last observation still contributes to the next fit, but with less weight.

Worked Example: Two Weighted Least-Squares Updates

Consider four hypothetical observations: \((x_i,y_i)=(0,0),(1,1),(2,2),(3,10)\). Use an intercept and slope, fixed \(s=1\), and \(c=1\). Start at \(\beta^{(0)}=(0,1)\), meaning intercept 0 and slope 1. Calculate two IRLS updates.

At the starting fit, the predicted responses are \(0,1,2,3\), so the residuals are \(0,0,0,7\). Since \(cs=1\), the weights are \(1,1,1,\) and \(1/7\). For a weighted straight-line fit, write \(S_0=\sum w_i\), \(S_x=\sum w_ix_i\), \(S_y=\sum w_iy_i\), \(S_{xx}=\sum w_ix_i^2\), and \(S_{xy}=\sum w_ix_iy_i\). With these weights:

S_0=\frac{22}{7},\quad S_x=\frac{24}{7},\quad S_y=\frac{31}{7},\quad S_{xx}=\frac{44}{7},\quad S_{xy}=\frac{65}{7}.

The weighted slope is \((S_0S_{xy}-S_xS_y)/(S_0S_{xx}-S_x^2)\). Substitution gives \((1430-744)/(968-576)=686/392=1.75\). The intercept is \((S_y-1.75S_x)/S_0=(31-42)/22=-0.5\). Thus the first update is \(\beta^{(1)}=(-0.5,1.75)\).

Using that updated line, the fitted values are \(-0.5,\ 1.25,\ 3,\ 4.75\), and the residuals are \(0.5,\ -0.25,\ -1,\ 5.25\). Their weights are \(1,\ 1,\ 1,\) and \(1/5.25=4/21\). Refit using these weights. The weighted sums are \(S_0=67/21,\ S_x=75/21,\ S_y=103/21,\ S_{xx}=141/21,\) and \(S_{xy}=225/21\). Therefore the slope is

\frac{67(225)-75(103)}{67(141)-75^2} =\frac{7350}{3822} =\frac{1225}{637} \approx 1.923.

The intercept is \((103-(1225/637)(75))/67=-8/13\approx-0.615\). The second update is therefore \(\beta^{(2)}\approx(-0.615,\ 1.923)\). This is a second iterate, not evidence that the algorithm has converged; another round would calculate weights from these new residuals and fit again.

What the Weights Mean—and Do Not Mean

A Huber weight is a computational expression of the loss’s changing slope. For a small residual, the score grows in proportion to the residual, as in least squares. For a residual beyond the cutoff, the score’s magnitude is \(c\), rather than growing without bound. Dividing that capped score by the residual produces the decreasing IRLS weight.

Worked Example: A Capped Score Contribution

Take \(c=1\), \(s=1\), and consider an observation whose predictor value is \(x_i=2\) and whose residual is \(r_i=10\). Its standardized residual is \(u_i=10\), so its Huber weight is \(1/10=0.1\). In the slope component of the estimating equation, its contribution is \(x_i\psi_c(u_i)=2(1)=2\).

Under least squares, the corresponding score contribution would be \(x_i u_i=2(10)=20\). Huber’s contribution is capped with respect to residual magnitude for this fixed predictor value; the observation is not discarded. The calculation also shows why the weight should not be read as a probability or as a claim that the observation is erroneous.

These are residual weights, not leverage corrections. They are based on how far a response lies from its current fitted value; they do not directly reduce the effect of an unusual predictor configuration. A case can have high leverage and a modest residual, and hence receive weight 1. Likewise, a low weight is not a rule for automatic deletion. The weights depend on the current fit, scale, and tuning constant, so the final report should make clear what settings were used.

Common Pitfall and Practical Checks

Using the wrong residual scale. The threshold is \(cs\) in response units, not \(c\), unless \(s=1\). Standardize residuals or use the equivalent original-unit formula consistently.
Keeping the first weights fixed. IRLS is iterative: after each weighted fit, recalculate residuals and weights from the new coefficients.
Confusing weights with deletion. A weight below 1 reduces an observation’s contribution; it does not remove the observation from the data.
Assuming a small residual weight handles leverage. The weight responds to the residual, not to how unusual the predictor values are.

A practical implementation should check that the weighted design remains numerically stable and should use a stated convergence tolerance. Since the weights are positive for finite residuals, a full-rank design retains a unique weighted least-squares solution in exact arithmetic. Very small weights can nevertheless make the calculation poorly conditioned in finite-precision arithmetic.

Key takeaway: Huber IRLS sets \(w_i=1\) inside the residual cutoff and \(w_i=cs/|r_i|\) outside it. Each iteration solves a weighted least-squares problem using weights from the current residuals, then updates those weights from the new fit.

Check Your Understanding

Use the definitions and update sequence above to answer the following.

For \(c=2\), \(s=3\), and residual \(r_i=9\), find the standardized residual and the IRLS weight.
Explain why the weight for a zero residual is defined as 1.
At an iteration, why must weights be calculated from the current coefficients rather than from the coefficients being solved for?
With fixed \(s\), show how \(\psi_c(u_i)=w_i u_i\) leads to the weighted least-squares normal equation.
Why does a small Huber residual weight not guarantee that an observation has low leverage?