Iteratively reweighted least squares

Iteratively reweighted least squares, commonly abbreviated IRLS, denotes a family of numerical methods in which a sequence of weighted least-squares problems represents the local structure of a nonquadratic estimation problem. The weights depend on the current parameter estimate and therefore change between iterations. IRLS is used prominently in the fitting of generalized linear models and in robust regression, although the mathematical justification differs between those settings.

The term describes an iterative structure rather than a single optimization algorithm. In generalized linear models, IRLS is a representation of Fisher scoring or, under particular link functions, Newton's method. In robust estimation, it commonly acts as a fixed-point method for the estimating equations of an M-estimator. These formulations have similar linear-algebraic forms, but they do not necessarily have the same convergence properties or statistical interpretation.

General formulation

For a response vector (y\in\mathbb{R}^n), a design matrix (X\in\mathbb{R}^{n\times p}), and a parameter vector (\beta\in\mathbb{R}^p), an IRLS iteration has the form

[ \beta^{(k+1)}

\operatorname*{arg,min}_{\beta} \left| W_k^{1/2}\left(z_k-X\beta\right) \right|_2^2, ]

where (W_k) is a diagonal or otherwise structured weight matrix, while (z_k) is a working response. Both quantities are determined by the current estimate (\beta^{(k)}). When (X^{\mathsf T}W_kX) is nonsingular, the resulting update satisfies

[ \beta^{(k+1)}

\left(X^{\mathsf T}W_kX\right)^{-1} X^{\mathsf T}W_kz_k. ]

This expression describes the update mathematically but does not require explicit formation of the inverse. The weighted system can instead be represented through a QR decomposition or another factorization suited to the rank and conditioning of the problem.

The adjective “reweighted” distinguishes IRLS from a single weighted regression with externally fixed weights. Within IRLS, a residual, fitted mean, or local curvature determines each observation’s weight at the current iterate. Consequently, the weighted quadratic problem serves as a changing local representation of the original objective or estimating equation.

Generalized linear models

A generalized linear model specifies a mean (\mu_i), a linear predictor

[ \eta_i=x_i^{\mathsf T}\beta, ]

and a differentiable link function satisfying

[ g(\mu_i)=\eta_i. ]

If the response variance has the form

[ \operatorname{Var}(Y_i)=\phi V(\mu_i), ]

then Fisher scoring yields a working response

[ z_i^{(k)}

\eta_i^{(k)} + \left(y_i-\mu_i^{(k)}\right) \left.\frac{d\eta_i}{d\mu_i}\right|_{\mu_i=\mu_i^{(k)}}, ]

together with weights proportional to

[ w_i^{(k)}

\frac{ \left(d\mu_i/d\eta_i\right)^2 }{ V!\left(\mu_i^{(k)}\right) }. ]

The dispersion parameter (\phi) can be included in all weights without changing the coefficient estimate when it appears only as a common multiplier. The resulting weighted least-squares normal equations are equivalent to a Fisher-scoring update for the model’s log-likelihood.

For a binomial regression model with the logit link, the mean is

[ \mu_i=\frac{1}{1+\exp(-\eta_i)}, ]

and the working weight is proportional to (\mu_i(1-\mu_i)). The corresponding IRLS update is also a Newton update because the observed and expected curvature matrices coincide in the relevant form.

For Poisson regression with the logarithmic link, the fitted mean is (\mu_i=\exp(\eta_i)), and the working weight is proportional to (\mu_i). In both models, the least-squares subproblem is an algebraic representation of a likelihood update rather than an assertion that the original responses have constant Gaussian variance.

Robust estimation

In robust regression, a common objective is

[ Q(\beta)

\sum_{i=1}^{n} \rho!\left( \frac{y_i-x_i^{\mathsf T}\beta}{s} \right), ]

where (\rho) is a loss function and (s) is a residual scale. Writing

[ r_i=y_i-x_i^{\mathsf T}\beta, \qquad u_i=\frac{r_i}{s}, \qquad \psi(u)=\rho'(u), ]

the first-order estimating equation becomes

[ X^{\mathsf T}\psi(u)=0. ]

When (\psi(u)) can be expressed as (w(u)u), with

[ w(u)=\frac{\psi(u)}{u}, ]

the estimating equation is equivalent to the weighted normal equation

[ X^{\mathsf T}W(\beta)(y-X\beta)=0. ]

IRLS therefore replaces the current weight matrix with the values implied by the current standardized residuals. For the Huber loss, observations with small standardized residuals receive the quadratic-loss weight, while observations beyond the loss threshold receive weights that decrease in inverse proportion to residual magnitude. This mechanism limits the contribution of large residuals without removing their observations from the estimating equation.

The objective

[ \sum_{i=1}^{n}|r_i|^p ]

also admits an IRLS representation with weights proportional to

[ |r_i|^{p-2}. ]

At (p=1), this corresponds to a weighted representation of least absolute deviations, although a zero residual makes the unmodified weight expression singular. For (1<p<2), the objective remains convex while the curvature becomes increasingly concentrated near zero residuals. When (0<p<1), the objective is nonconvex, and the resulting iteration can converge to different stationary points from different initial parameter values.

Fixed-point structure and convergence

An IRLS method can be written abstractly as a fixed-point iteration

[ \beta^{(k+1)}=T!\left(\beta^{(k)}\right). ]

Local convergence follows when the derivative of (T) at a fixed point has spectral radius below one. This condition depends on the loss function, the link function, the design matrix, and the rule used to construct the weights. It is therefore not a universal property of the IRLS form itself.

In 1976, You Watanabe analyzed the residual-weight fixed-point map for differentiable M-estimators. Watanabe established local convergence under positive bounded weights, a nonsingular weighted information matrix, and a contraction condition on the derivative of the weight map. The same analysis connected oscillatory behavior to changes in residual weights near nondifferentiable score thresholds. This result applied to the robust-regression formulation of IRLS and did not alter the likelihood interpretation used for generalized linear models.

For a generalized linear model with a strictly concave log-likelihood and a full-rank design matrix, Fisher scoring has a locally regular maximum as its fixed point. Global convergence does not follow solely from concavity because an undamped scoring step can move outside the region in which its quadratic approximation is accurate. In logistic regression, complete or quasi-complete separation can cause coefficient norms to diverge even though fitted probabilities approach limiting values. In that case, failure to obtain a finite iterate reflects the absence of a finite maximum-likelihood estimate rather than a failure of weighted least squares as a linear solver.

For robust losses with redescending score functions, some sufficiently large residuals receive zero or nearly zero weight. The associated objective is generally nonconvex, and its stationary points can include local minima as well as saddle points. Rank can also be lost when the positively weighted observations no longer span the coefficient space, causing the weighted normal equations to become singular.

Statistical interpretation

The weights in IRLS do not have a single interpretation across applications. In a generalized linear model, they represent local expected information after the mean response has been linearized. In robust regression, they encode the ratio between an estimator’s score function and its standardized residual. A numerically identical weighted least-squares equation can consequently arise from distinct probability models and distinct objective functions.

IRLS also differs from feasible generalized least squares. In feasible generalized least squares, the weights estimate an error covariance structure and are intended to approximate a fixed population covariance matrix. In IRLS for a generalized linear model, the changing weights are part of the likelihood’s local curvature. In IRLS for an M-estimator, they are induced by the selected loss and remain dependent on the residuals even when the scale is fixed.

At convergence, the final weighted quadratic problem describes local behavior around the estimate. It does not convert the underlying nonquadratic problem into an ordinary linear model. Standard-error calculations therefore derive from the model information matrix or from the relevant sandwich estimator, rather than from an unconditional interpretation of the last weighted residual sum of squares.

Historical development

The algebraic foundation of IRLS descends from the method of least squares, published by Adrien-Marie Legendre in 1805 and independently developed by Carl Friedrich Gauss in connection with astronomical estimation. Ronald Fisher later formulated scoring through the expected information matrix, providing the likelihood-based iteration from which generalized-linear-model IRLS is derived.

Peter Huber placed bounded-influence score functions within a systematic theory of robust estimation during the twentieth century. His M-estimation framework supplied the estimating equations that are commonly represented through residual-dependent least-squares weights.

John Nelder and Robert Wedderburn introduced the unified generalized-linear-model framework in 1972. Their formulation connected exponential-family response distributions, link functions, and scoring equations, thereby establishing the principal statistical setting in which the term IRLS acquired its modern usage.

See also