Generalized-norm conjugate direction solver

Algorithm Pseudo-code

The steps of actual solver have a structure very similar steps as the $L_2$ solver. The chief difference lies in the generalization to allow for different norms in the gradient calculation. The generalized form of the gradient for a given norm $C(\cdot)$ is given by

$${\bf g}^{(m)} =  \Delta{\bf m} =  {\bf F}^T  C'({\bf r}) ,$$ (14)

where $F$ is the forward operator, ${\bf m}$ is the model value and ${\bf g}^{(m)}$ is the gradient. In addition, our code allows for an iterative plane-search to update ${\bf r}$, the intermediate estimated residual. The pseudo-code in Algorithm 1 summarizes the implementation for a simple solver (see Claerbout (2009) for more details). The benefit of our framework is that it can be easily modified to allow regularization and preconditioning without extensively changing the main solver algorithm. Such modifications primarily affect the stepper code. Fortran implementations of plane-search stepper function for both conjugate-direction and our generalized norm stepper are given in the appendices. Note that we approximate the data-space value of the gradient throughout the entire plane-search by a constant, to avoid re-evaluating the forward operator every iteration.

$${\bf F g}^{(m)}_i \approx {\bf F g}^{(m)}_0$$ (15)

This assumption forces our plane-search to remain in a local neighborhood around ${\bf r}_1$, which is the desired outcome for this local approximation anyway.

Generalized-norm conjugate direction solver