Normal equations

The iterative reweighted least-squares (IRLS) algorithm provides a mean by which linear systems can be solved by minimizing the L^p norm of the residuals ($1\leq\!p\!\leq 2$). Let us suppose that we try to solve the set of equations:

$$A.x = y \mbox{\hspace{2.0cm}} x\in {\cal R}^m,\; y\in {\cal R}^n,\; A\in {\cal R}^{n\times m}\:.$$

To find the L^p solution of this system we minimize the L^p norm of the residuals:

$${\cal N}_p(x) = (\sum_{i=1}^n \vert y_i - \sum_{j=1}^m A_{ij}x_j\vert^p)^{{1\over p}} = {\cal L}_p(y-A.x)\:.$$

${\cal L}_p$ is not a norm for p<1, so, for the moment, I will only consider $p\geq 1$ . I also omit the exponent 1/p . To establish the normal equations, we set:

$${\partial{\cal N}_p \over\partial x_k} = 0 \mbox{\hspace{2.0cm} for k=1 to m}\:.$$

The partial derivatives can be expressed as follows:

$${\partial{\cal N}_p \over\partial x_k} = -p\sum_{i=1}^n A_{ik}{\rm sign}(r(i))\:\vert r(i)\vert^{p-1}\:.$$

where $r(i)=y_i-\sum_{j=1}^mA_{ij}x_j$ are the residuals. Then:

$${\partial{\cal N}_p \over\partial x_k} = -p\sum_{i=1}^n A_{ik}r(i) \vert r(i)\vert^{p-2}\:.$$ (1)

To form the normal equations we can now replace r(i) by its expression:

$$\sum_{i=1}^n A_{ik}\vert r(i)\vert^{p-2}(\sum_{j=1}^mA_{ij}x_j) = \sum_{i=1}^n A_{ik}\vert r(i)\vert^{p-2}y_i$$

or, after some reorganization:

$$\sum_{j=1}^m(\sum_{i=1}^nA_{ik}\vert r(i)\vert^{p-2}A_{ij})x_j = \sum_{i=1}^n A_{ik}\vert r(i)\vert^{p-2}y_i\:.$$ (2)

There are m such equations (k=1 to m). They can be expressed shortly in a matrix formulation:  

$$W(i) = \vert r(i)\vert^{p-2} \mbox{\hspace{1.0cm} for i=1 to n ,}$$ (3)

$$W = {\rm diag}(W(1),\cdots,W(n))\;,$$

and finally:  

$$A^TWA.x = A^TW.y\;.$$ (4)

This formulation is implicit and non-linear since the residuals r(i), and so the matrix W, depend on the unknown vector x. However, for the classical least-squares inversion (p=2), W is the identity matrix, and equation (4) is the usual least-squares equation:  

$$A^TA.x = A^T.y\;.$$ (5)