# Definition and Convergence of Iteratively Reweighted Least Squares

I've been using iteratively reweighted least squares (IRLS) to minimize functions of the following form,

$J(m) = \sum_{i=1}^{N} \rho \left(\left| x_i - m \right|\right)$

where $N$ is the number of instances of $x_i \in \mathbb{R}$, $m \in \mathbb{R}$ is the robust estimate that I want, and $\rho$ is a suitable robust penalty function. Let's say it's convex (though not necessarily strictly) and differentiable for now. A good example of such a $\rho$ is the Huber loss function.

What I've been doing is differentiating $J(m)$ with respect to $m$ (and manipulating) to obtain,

$\frac{dJ}{dm}= \sum_{i=1}^{N} \frac{\rho'\left( \left|x_i-m\right|\right) }{\left|x_i-m\right|} \left( x_i-m \right)$

and iteratively solving this by setting it equal to 0 and fixing weights at iteration $k$ to $w_i(k) = \frac{\rho'\left( \left|x_i-m{(k)}\right|\right) }{\left|x_i-m{(k)}\right|}$ (note that the perceived singularity at $x_i=m{(k)}$ is really a removable singularity in all $\rho$'s I might care about). Then I obtain,

$\sum_{i=1}^{N} w_i(k) \left( x_i-m{(k+1)} \right)=0$

and I solve to obtain, $m(k+1) = \frac{\sum_{i=1}^{N} w_i(k) x_i}{ \sum_{i=1}^{N} w_i(k)}$.

I repeat this fixed point algorithm until "convergence". I will note that if you get to a fixed point, you are optimal, since your derivative is 0 and it's a convex function.

1. Is this the standard IRLS algorithm? After reading several papers on the topic (and they were very scattered and vague about what IRLS is) this is the most consistent definition of the algorithm I can find. I can post the papers if people want, but I actually didn't want to bias anyone here. Of course, you can generalize this basic technique to many other types of problems involving vector $x_i$'s and arguments other than $\left|x_i-m{(k)}\right|$, providing the argument is a norm of an affine function of your parameters. Any help or insight would be great on this.
2. Convergence seems to work in practice, but I have a few concerns about it. I've yet to see a proof of it. After some simple Matlab simulations I see that one iteration of this is not a contraction mapping (I generated two random instances of $m$ and computing $\frac{\left|m_1(k+1) - m_2(k+1)\right|}{\left|m_1(k)-m_2(k)\right|}$ and saw that this is occasionally greater than 1). Also the mapping defined by several consecutive iterations is not strictly a contraction mapping, but the probability of the Lipschitz constant being above 1 gets very low. So is there a notion of a contraction mapping in probability? What is the machinery I'd use to prove that this converges? Does it even converge?

Any guidance at all is helpful.

Edit: I like the paper on IRLS for sparse recovery/compressive sensing by Daubechies et al. 2008 "Iteratively Re-weighted Least Squares Minimization for Sparse Recovery" on the arXiv. But it seems to focus mostly on weights for nonconvex problems. My case is considerably simpler.

• Looking at the wiki page on IRWLS I struggle to the difference between the procedure you describe and IRWLS (they just use $|y_i-\pmb x_i'\pmb\beta|^2$ as their particular $\rho$ function). Can you explain in what ways you think the algorithm you propose is different from IRWLS? Commented May 24, 2014 at 7:38
• I never stated that it was different, and if I implied it, I didn't mean to. Commented May 27, 2014 at 6:48

As for your first question, one should define "standard", or acknowledge that a "canonical model" has been gradually established. As a comment indicated, it appears at least that the way you use IRWLS is rather standard.

As for your second question, "contraction mapping in probability" could be linked (however informally) to convergence of "recursive stochastic algorithms". From what I read, there is a huge literature on the subject mainly in Engineering. In Economics, we use a tiny bit of it, especially the seminal works of Lennart Ljung -the first paper was Ljung (1977)- which showed that the convergence (or not) of a recursive stochastic algorithm can be determined by the stability (or not) of a related ordinary differential equation.

(what follows has been re-worked after a fruitful discussion with the OP in the comments)

Convergence

I will use as reference Saber Elaydi "An Introduction to Difference Equations", 2005, 3d ed. The analysis is conditional on some given data sample, so the $x's$ are treated as fixed.

The first-order condition for the minimization of the objective function, viewed as a recursive function in $m$, $$m(k+1) = \sum_{i=1}^{N} v_i[m(k)] x_i, \;\; v_i[m(k)] \equiv \frac{w_i[m(k)]}{ \sum_{i=1}^{N} w_i[m(k)]} \qquad [1]$$

has a fixed point (the argmin of the objective function). By Theorem 1.13 pp 27-28 of Elaydi, if the first derivative with respect to $m$ of the RHS of $[1]$, evaluated at the fixed point $m^*$, denote it $A'(m^*)$, is smaller than unity in absolute value, then $m^*$ is asymptotically stable (AS). More over by Theorem 4.3 p.179 we have that this also implies that the fixed point is uniformly AS (UAS).
"Asymptotically stable" means that for some range of values around the fixed point, a neighborhood $(m^* \pm \gamma)$, not necessarily small in size, the fixed point is attractive , and so if the algorithm gives values in this neighborhood, it will converge. The property being "uniform", means that the boundary of this neighborhood, and hence its size, is independent of the initial value of the algorithm. The fixed point becomes globally UAS, if $\gamma = \infty$.
So in our case, if we prove that

$$|A'(m^*)|\equiv \left|\sum_{i=1}^{N} \frac{\partial v_i(m^*)}{\partial m}x_i\right| <1 \qquad [2]$$

we have proven the UAS property, but without global convergence. Then we can either try to establish that the neighborhood of attraction is in fact the whole extended real numbers, or, that the specific starting value the OP uses as mentioned in the comments (and it is standard in IRLS methodology), i.e. the sample mean of the $x$'s, $\bar x$, always belongs to the neighborhood of attraction of the fixed point.

We calculate the derivative $$\frac{\partial v_i(m^*)}{\partial m} = \frac {\frac{\partial w_i(m^*)}{\partial m}\sum_{i=1}^{N} w_i(m^*)-w_i(m^*)\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}}{\left(\sum_{i=1}^{N} w_i(m^*)\right)^2}$$

$$=\frac 1{\sum_{i=1}^{N} w_i(m^*)}\cdot\left[\frac{\partial w_i(m^*)}{\partial m}-v_i(m^*)\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}\right]$$ Then

$$A'(m^*) = \frac 1{\sum_{i=1}^{N} w_i(m^*)}\cdot\left[\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}x_i-\left(\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}\right)\sum_{i=1}^{N}v_i(m^*)x_i\right]$$

$$=\frac 1{\sum_{i=1}^{N} w_i(m^*)}\cdot\left[\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}x_i-\left(\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}\right)m^*\right]$$

and

$$|A'(m^*)| <1 \Rightarrow \left|\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}(x_i-m^*)\right| < \left|\sum_{i=1}^{N} w_i(m^*)\right| \qquad [3]$$

we have

\begin{align}\frac{\partial w_i(m^*)}{\partial m} = &\frac{-\rho''(|x_i-m^*|)\cdot \frac {x_i-m^*}{|x_i-m^*|}|x_i-m^*|+\frac {x_i-m^*}{|x_i-m^*|}\rho'(|x_i-m^*|)}{|x_i-m^*|^2} \\ \\ &=\frac {x_i-m^*}{|x_i-m^*|^3}\rho'(|x_i-m^*|) - \rho''(|x_i-m^*|)\cdot \frac {x_i-m^*}{|x_i-m^*|^2} \\ \\ &=\frac {x_i-m^*}{|x_i-m^*|^2}\cdot \left[\frac {\rho'(|x_i-m^*|)}{|x_i-m^*|}-\rho''(|x_i-m^*|)\right]\\ \\ &=\frac {x_i-m^*}{|x_i-m^*|^2}\cdot \left[w_i(m^*)-\rho''(|x_i-m^*|)\right] \end{align}

Inserting this into $[3]$ we have

$$\left|\sum_{i=1}^{N}\frac {x_i-m^*}{|x_i-m^*|^2}\cdot \left[w_i(m^*)-\rho''(|x_i-m^*|)\right](x_i-m^*)\right| < \left|\sum_{i=1}^{N} w_i(m^*)\right|$$

$$\Rightarrow \left|\sum_{i=1}^{N}w_i(m^*)-\sum_{i=1}^{N}\rho''(|x_i-m^*|)\right| < \left|\sum_{i=1}^{N} w_i(m^*)\right| \qquad [4]$$

This is the condition that must be satisfied for the fixed point to be UAS. Since in our case the penalty function is convex, the sums involved are positive. So condition $[4]$ is equivalent to

$$\sum_{i=1}^{N}\rho''(|x_i-m^*|) < 2\sum_{i=1}^{N}w_i(m^*) \qquad [5]$$

If $\rho(|x_i-m|)$ is Hubert's loss function, then we have a quadratic ($q$) and a linear ($l$) branch,

$$\rho(|x_i-m|)=\cases{ (1/2)|x_i- m|^2 \qquad\;\;\;\; |x_i-m|\leq \delta \\ \\ \delta\big(|x_i-m|-\delta/2\big) \qquad |x_i-m|> \delta}$$

and

$$\rho'(|x_i-m|)=\cases{ |x_i- m| \qquad |x_i-m|\leq \delta \\ \\ \delta \qquad \qquad \;\;\;\; |x_i-m|> \delta}$$

$$\rho''(|x_i-m|)=\cases{ 1\qquad |x_i-m|\leq \delta \\ \\ 0 \qquad |x_i-m|> \delta}$$

$$\cases{ w_{i,q}(m) =1\qquad \qquad \qquad |x_i-m|\leq \delta \\ \\ w_{i,l}(m) =\frac {\delta}{|x_i-m|} <1 \qquad |x_i-m|> \delta}$$

Since we do not know how many of the $|x_i-m^*|$'s place us in the quadratic branch and how many in the linear, we decompose condition $[5]$ as ($N_q + N_l = N$)

$$\sum_{i=1}^{N_q}\rho_q''+\sum_{i=1}^{N_l}\rho_l'' < 2\left[\sum_{i=1}^{N_q}w_{i,q} +\sum_{i=1}^{N_l}w_{i,l}\right]$$

$$\Rightarrow N_q + 0 < 2\left[N_q +\sum_{i=1}^{N_l}w_{i,l}\right] \Rightarrow 0 < N_q+2\sum_{i=1}^{N_l}w_{i,l}$$

which holds. So for the Huber loss function the fixed point of the algorithm is uniformly asymptotically stable, irrespective of the $x$'s. We note that the first derivative is smaller than unity in absolute value for any $m$, not just the fixed point.

What we should do now is either prove that the UAS property is also global, or that, if $m(0) = \bar x$ then $m(0)$ belongs to the neighborhood of attraction of $m^*$.

• Thanks for the response. Give me some time to analyze this answer. Commented May 27, 2014 at 6:49
• Certainly. After all, the question waited 20 months. Commented May 27, 2014 at 8:49
• Yeah, I was reminded of the problem and decided to put up a bounty. :) Commented May 27, 2014 at 17:48
• Lucky me. I wasn't there 20 months ago - I would have taken up this question, bounty or not. Commented May 27, 2014 at 21:49
• Thanks so much for this response. It's looking like, so far, that you've earned the bounty. BTW, your indexing on the derivative of $v_i$ w.r.t $m$ is notationally weird. Couldn't the summations on the second line of this use another variable, such as $j$? Commented May 28, 2014 at 1:16

Is this the standard IRLS algorithm?

IRLS algorithms can be generally said to find/approach a solution to a minimization problem by using an iterative process that generates a sequence of solutions to a weighted least squares problem.

There is no standard, but your approach seems to fit this general description.

What is maybe atypical, different, is that a lot IRLS problems have a different cost function or distance measure. A common case is when the cost function is a $$\ell_1$$, $$\ell_\infty$$, or other $$\ell_p$$ norm instead of a general function $$\rho$$, and the functions in which the approximation is to be found can be multidimensional and non-linear (replacing your $$m$$ with a non-linear function of regressor variables and unknown parameters).

Convergence seems to work in practice, but I have a few concerns about it. I've yet to see a proof of it.

Convergence will depend on the function $$\rho(z)$$ and convexity is not sufficient. A counter example is when $$\rho(z) = z^p$$ with $$p>3$$. That will make the weights as following (where $$z$$ is the difference $$x-m$$)

$$w(z) = \frac{pz^{p-1}}{z} = pz^{p-2}$$

The rate at which the weights increase is larger for the points further away. So a change of $$m$$ will shift the values $$z_i$$ and make the weights of the furthest point larger and make the weight of that point more important in the next itteration. The solution will diverge and end up switching back and forth between two points.

Below is an example when your cost function is a power of 4. You see that the algorithm does not converge.

In addition the example show the convergence for an alternative algorithm (in fact, it is the first printed IRLS algorithm from 1961 by Lawson which is described in The Lawson Algorithm and Extensions Rice and Usow 1968)

    layout(matrix(1:4,2))

set.seed(1)
xi = runif(10)*20-10
xi

new_weights = function(wi) {
m = sum(wi*xi)
ri = xi-m
wi = ri^2/sum(ri^2)
}

new_weights_lawson = function(wi) {
m = sum(wi*xi)
ri = abs(xi-m)
wi = (wi*ri)^(2/3)/sum((wi*ri)^(2/3))
}

#### itterated reweight
wi = rep(0.1,10) ### start
mv = sum(wi*xi)
for (i in 1:30) {
wi = new_weights(wi)
mv = c(mv,sum(wi*xi))
}

plot(mv, type = "l", main = expression(w[i+1,j] == over("|"*r[i,j]*"|"^(p-2), sum("|"*r[i,j]*"|"^(p-2),j=1,n))))

#### itterated reweight lwawson
wi = rep(0.1,10) ### start
mv = sum(wi*xi)
for (i in 1:30) {
wi = new_weights_lawson(wi)
mv = c(mv,sum(wi*xi))
}

plot(mv, type = "l", main = expression(w[i+1,j] == over("|"*w[i,j]*r[i,j]*"|"^((p-2)/(p-1)), sum("|"*w[i,j]*r[i,j]*"|"^((p-2)/(p-1)),j=1,n))))

cost = function(m) {
ri = abs(xi-m)
cost = sum(ri^4)
}
cost = Vectorize(cost)
m = seq(0.6,0.8,0.001)
plot(m,cost(m), type = "l")


### Some doodle

Here is some doodling and graph that I used while figuring out what is going on

set.seed(1)

### cost function
Cost = function(x) {
if (x>0) {
0.1*1/3*x^3
}
else {
1/1.5*abs(x)^1.5
}
}
Cost = Vectorize(Cost)

### derivative of cost function
dCost = function(x) {
if (x>0) {
0.1*x^2
}
else {
-abs(x)^0.5
}
}
dCost = Vectorize(dCost)

### cost as function of m
costm = function(m) {
sum(Cost(xi-m))
}
costm = Vectorize(costm)

### compute the map from old m to new m
solvem = function(m) {
w = dCost(xi-m)/(xi-m)
mn = sum(w*xi)/sum(w)
mn
}
solvem = Vectorize(solvem)

layout(matrix(1:4,2))

x = seq(-10,10,0.01)
dJ = dCost(x)
plot(x,dJ, type = "l", xlab = "xi-m", main = "derivative of cost function")

xi = runif(10)*20-10
dJi = dCost(xi)
points(xi,dJi, pch = 20)

m = seq(-10,10,0.001)
mn = solvem(m)
plot(m,costm(m) , main = "cost as function of m")

plot(m,mn, type = "l", main = "mapping of old to new m", xlab = "old m", ylab = "new m")
lines(c(-10,10),c(-10,10), col = 2)

sel = which.min(costm(m))

lines(c(-10,10),c(10,-10)+2*m[sel], col = 2)

plot(x,dCost(x)/x, type = "l", xlab = "z", main = "weights as function of z = x-m", ylab = "dCost(z)/z")

m[sel]
xi
min(costm(m))

• Amazing response. Thank you for the counter-example and introduction to the Lawson algorithm. I will check correctness soon but assigning the bounty now for the effort and new knowledge, thank you. Commented Nov 9, 2022 at 19:11
• Now that I look at what you're saying, this is interesting. It's the exact opposite of the cases I was worried about not converging, i.e. those that were more like L1 in nature, like Huber losses. But one counterexample is enough. Well-deserved bounty on this one. Commented Nov 9, 2022 at 22:47
• @ChrisA. I have added some off my notes (I called it a doodle) that I made while trying to figure out an answer. Before getting to the power of >3 I used some other functions to see how it looks like. It appears like you get that mapping that is not contractive because of the peaks that occur when a $x_i$ is close to $m$. But you do get that the next $m$ is always at least closer than before (within the region bounded by the red lines). I could not find a way to proof convergence under certain conditions, so I started finding a way to get a situation when it does not work. Commented Nov 9, 2022 at 23:09