Iteratively reweighted least squares

IRLS is used to find the maximum likelihood estimates of a generalized linear model, and in robust regression to find an M-estimator, as a way of mitigating the influence of outliers in an otherwise normally distributed data set, for example, by minimizing the least absolute errors rather than the least square errors.

One of the advantages of IRLS over linear programming and convex programming is that it can be used with Gauss–Newton and Levenberg–Marquardt numerical algorithms.

Examples

L1 minimization for sparse recovery

IRLS can be used for ℓ₁ minimization and smoothed ℓ_p minimization, p < 1, in compressed sensing problems. It has been proved that the algorithm has a linear rate of convergence for ℓ₁ norm and superlinear for ℓ_t with t < 1, under the restricted isometry property, which is generally a sufficient condition for sparse solutions.

Lp norm linear regression

To find the parameters β = (β₁, …,β_k)^T which minimize the L^p norm for the linear regression problem, $${\displaystyle \underset{\mathit{\beta }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\Vert \mathbf{y}-X\mathit{\beta }{\Vert }_p=\underset{\mathit{\beta }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\sum _{i=1}^n{\left|y_i-X_i\mathit{\beta }\right|}^p,}$$ the IRLS algorithm at step t + 1 involves solving the weighted linear least squares problem $${\displaystyle {\mathit{\beta }}^{(t+1)}=\underset{\mathit{\beta }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\sum _{i=1}^n{w}_i^{(t)}{\left|y_i-X_i\mathit{\beta }\right|}^2=(X^{\mathrm{T}}{W}^{(t)}X{)}^{-1}{X}^{\mathrm{T}}{W}^{(t)}\mathbf{y},}$$ where W^(t) is the diagonal matrix of weights, usually with all elements set initially to $${\displaystyle w_i^{(0)}=1}$$ and updated after each iteration to $${\displaystyle w_i^{(t)}=|y_i-X_i{\mathit{\beta }}^{(t)}{|}^{p-2}.}$$

In the case p = 1, this corresponds to least absolute deviation regression (in this case, the problem would be better approached by use of linear programming methods, so the result would be exact) and the formula is $${\displaystyle w_i^{(t)}=\frac{1}{|y_i-X_i{\mathit{\beta }}^{(t)}|}.}$$

To avoid dividing by zero, regularization must be done, so in practice the formula is $${\displaystyle w_i^{(t)}=\frac{1}{max\left\{\delta ,\left|y_i-X_i{\mathit{\beta }}^{(t)}\right|\right\}}.}$$ where ${\displaystyle \delta }$ is some small value, like 0.0001. Note the use of ${\displaystyle \delta }$ in the weighting function is equivalent to the Huber loss function in robust estimation.