2 added 286 characters in body

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.

1

# 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?