IRLS algorithm

At the iteration k+1, the algorithm solves:

 

A^TW^kA.x^k+1 = A^TW^k.y (6)

by taking:

• W⁰ = I_n (Identity matrix), at the first iteration,
• W^k formed with the residuals of iteration k (r^k=y-Ax^k), at the iteration k+1 .

Byrd and Payne (1979) showed that this algorithm is convergent under two conditions:

• W(i) must be non-increasing in |r(i)|,
• W(i) must be bounded for all i.

The first condition is satisfied for $p\leq 2$; to assure the second condition, Huber (1981) suggested replacing equation (3) with:  

$$W(i) = \left\{ \begin{array} {ll} \vert r(i)\vert^{p-2} & ... ...mbox{if $\vert r(i)\vert\leq \varepsilon$} \end{array} \right.$$ (7)

for a small positive constant $\varepsilon$.
Equations (6) are the normal equations of the least-squares minimization of:

$${\cal N}_2^k(x) = (y-Ax)^TW^k(y-Ax) \:.$$

W^k has the role of a weighting matrix, the inverse of a covariance matrix. By giving the large residuals small weights (or large variances), it reduces their role; on the contrary, well predicted values will be given large weights (small variances). Aberrant values, receiving low weights, will have a small influence: this makes the algorithm robust, less sensitive to noise bursts in the data. This is especially true for the extreme case p=1, as it had been predicted by Claerbout and Muir (1973).

However, the L¹ norm is not strictly convex; so, the L¹ minimization problem might have non-unique solutions. The IRLS algorithm will lead us to one of these solutions. The fact that the first iteration (W⁰=I_n) is a least-squares inversion makes me think that this algorithm will lead us to a solution close to the L² solution. The choice of W⁰(i)=|y(i)|^p-2 is not reasonable: during the algorithm, x stays close to 0, and the algorithm converges (if it does) to a vector $x_{\infty}$ also close to zero, the residuals staying close to the original vector y.