Model fitting by least squares

Conjugate-direction theory for programmers

Fourier-transformed variables are often capitalized. This convention is helpful here, so in this subsection only, we capitalize vectors transformed by the $\bold F$ matrix. As everywhere, a matrix, such as $\bold F$, is printed in boldface type but in this subsection, vectors are not printed in boldface. Thus, we define the solution, the solution step (from one iteration to the next), and the gradient by:

$$X = \bold F \ x\ \quad\quad\quad {\rm modeled\ data\ } = \bold F\ {\rm model}$$ (66)

$$S_j = \bold F \ s_j \quad\quad\quad {\rm solution\ change}$$ (67)

$$G_j = \bold F \ g_j \quad\quad\quad \Delta\bold r = \bold F\Delta\bold m$$ (68)

A linear combination in solution space, say $s+g$, corresponds to $S+G$ in the conjugate space, the data space, because $S+G = \bold F s + \bold F g = \bold F(s+g)$. According to equation (51), the residual is the modeled data minus the observed data.

$$R \eq \bold F x \ -\ D \eq X \ -\ D$$ (69)

The solution $x$ is obtained by a succession of steps $s_j$, say:

$$x \eq s_1 \ +\ s_2 \ +\ s_3 \ +\ \cdots$$ (70)

The last stage of each iteration is to update the solution and the residual:

$${\rm solution\ update:} \quad\quad\quad x \ \leftarrow x +\ s$$ (71)

$${\rm residual\ update:} \quad\quad\quad R \ \leftarrow R +\ S$$ (72)

The gradient vector $g$ is a vector with the same number of components as the solution vector $x$. A vector with this number of components is:

$$g = \bold F\T \ R \eq \hbox{gradient}$$ (73)

$$G = \bold F \ g \eq \hbox{conjugate gradient} \ =\ \Delta r$$ (74)

The gradient $g$ in the transformed space is $G$, also known as the conjugate gradient.

What is our solution update $\Delta \bold x=\bold s$? It is some unknown amount $\alpha$ of the gradient $\bold g$ plus another unknown amount $\beta$ of the previous step $\bold s$. Likewise in residual space.

$$\Delta \bold x = \alpha \bold g + \beta \bold s \quad\quad \hbox{model space}$$ (75)

$$\Delta \bold r = \alpha \bold G + \beta \bold S \quad\quad \hbox{data space}$$ (76)

The minimization (56) is now generalized to scan not only in a line with $\alpha$, but simultaneously another line with $\beta$. The combination of the two lines is a plane. We now set out to find the location in this plane that minimizes the quadratic $Q$.

$$Q(\alpha ,\beta ) \eq ( R + \alpha G + \beta S) \ \cdot\ (R + \alpha G + \beta S )$$ (77)

The minimum is found at $\partial Q / \partial \alpha \,=\,0$ and $\partial Q / \partial \beta \,=\,0$, namely,

$$0 = G \ \cdot\ ( R + \alpha G + \beta S )$$ (78)

$$0 = S \ \cdot\ ( R + \alpha G + \beta S )$$ (79)

$$-\ \ \left[ \begin{array}{c} (G\cdot R) \\ (S\cdot R) \... ...eft[ \begin{array}{c} \alpha \\ \beta \end{array} \right]$$ (80)

Equation (81) is a set of two equations for $\alpha$ and $\beta$. Recall the inverse of a $2\times 2$ matrix, equation (111) and get:

$$\left[ \begin{array}{c} \alpha \\ \beta \end{array} \rig... ...egin{array}{c} (G\cdot R) \\ (S\cdot R) \end{array} \right]$$ (81)

The many applications in this book all need to find $\alpha$ and $\beta$ with (81), and then update the solution with (71) and update the residual with (72). Thus, we package these activities in a subroutine named cgstep(). To use that subroutine, we have a computation template with repetitive work done by subroutine cgstep(). This template (or pseudocode) for minimizing the residual $\bold 0\approx \bold r= \bold F \bold x - \bold d$ by the conjugate-direction method is:

		 $\bold r \quad\longleftarrow\quad \bold F \bold x - \bold d$ 


iterate {                                                     

$\Delta \bold x \quad\longleftarrow\quad \bold F\T\ \bold r$ 

$\Delta \bold r\ \quad\longleftarrow\quad \bold F \ \Delta \bold x$ 

$(\bold x,\bold r) \quad\longleftarrow\quad {\rm cgstep} (\bold x,\bold r, \Delta\bold x,\Delta\bold r )$        

} 

where the subroutine cgstep() remembers the previous iteration and works out the step size and adds in the proper proportion of the $\Delta \bold x$ of the previous step.

Model fitting by least squares