Iteratively reweighted least squares

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.[2][3] However, in most practical situations, the restricted isometry property is not satisfied.

To find the parameters β = (β₁, …,β_k)^T which minimize the L^p norm for the linear regression problem,

${\displaystyle {\underset {\boldsymbol {\beta }}{\operatorname {arg\,min} }}{\big \|}\mathbf {y} -X{\boldsymbol {\beta }}\|_{p}={\underset {\boldsymbol {\beta }}{\operatorname {arg\,min} }}\sum _{i=1}^{n}\left|y_{i}-X_{i}{\boldsymbol {\beta }}\right|^{p},}$

the IRLS algorithm at step t + 1 involves solving the weighted linear least squares problem:[4]

${\displaystyle {\boldsymbol {\beta }}^{(t+1)}={\underset {\boldsymbol {\beta }}{\operatorname {arg\,min} }}\sum _{i=1}^{n}w_{i}^{(t)}\left|y_{i}-X_{i}{\boldsymbol {\beta }}\right|^{2}=(X^{\rm {T}}W^{(t)}X)^{-1}X^{\rm {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)}={\big |}y_{i}-X_{i}{\boldsymbol {\beta }}^{(t)}{\big |}^{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,[5] so the result would be exact) and the formula is:

${\displaystyle w_{i}^{(t)}={\frac {1}{{\big |}y_{i}-X_{i}{\boldsymbol {\beta }}^{(t)}{\big |}}}.}$

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}{\boldsymbol {\beta }}^{(t)}\right|\right\}}}.}$

where ${\displaystyle \delta }$ is some small value, like 0.0001.[5] Note the use of ${\displaystyle \delta }$ in the weighting function is equivalent to the Huber loss function in robust estimation. [6]