Dix inversion constrained by L1-norm optimization

Dix inversion by conjugate direction $L1$ method

Choosing proper parameters for hybrid $l1/l2$ method and IRLS is still quite empirical, even when we understand their physical meanings. Instead, conjugate direction methods do not need setting parameters. Thus, we develop a similar conjugate direction method in $L1$ sense. The pseudo code of this method is given in Table 1.

Table 1: Pseudo Code - Conjugate direction $L1$ solver using Weighted Median

Initialization :

${\bf m}  = {\bf m}_{init}$

${\bf r}^{(d)} = \boldsymbol{F} {\bf m} - {\bf d}$

${\bf s}   = \boldsymbol{0}$

${\bf s}^{(d)} = \boldsymbol{0}$

Iteration i :

     ${\bf r}^{(t)} = sgn({\bf r}^{(t)})$

     ${\bf g}   = \boldsymbol{F'}{\bf r}^{(t)}                    *$

     ${\bf g}^{(d)} = \boldsymbol{F}{\bf g}                      *$

     $(\alpha, \beta) =  \underset{(\alpha, \beta)}{\operatorname{argmin}} \vert\vert{\bf r}^{(t)}+ \alpha {\bf g}^{(d)} + \beta {\bf s}^{(d)}\vert\vert _1$

     ${\bf s}^{(d)} = \alpha {\bf g}^{(d)} + \beta {\bf s}^{(d)}$

     ${\bf r}^{(d)} = {\bf r}^{(d)}+{\bf s}^{(d)}$

     ${\bf s}   = \alpha {\bf g} + \beta {\bf s}$

     ${\bf m}  = {\bf m} + {\bf s}$

The structure of the conjugate-direction $L1$ method is similar to the $L2$ conjugate-direction solver given by Claerbout (2008). The main difference arises in the part of plane search.

Given the gradient $\bold g^{(d)}$, previous step $\bold s^{(d)}$, and the current residual $\bold r^{(d)}$, we construct the $2\times N$ matrix $\bold B = [\bold g^{(d)}  \bold s^{(d)}]$ and the column vector $[\alpha  \beta]'$. We seek to find $[\alpha  \beta]'$ that minimizes $\bold 0 \approx \bold r^{(d)} + \bold B [\alpha  \beta]'$ in the L1-norm sense. This bivariate regression embedded in the plane search is solved in an iterate manner.

At the ultimate solution of the bivariate regression there will be two basis equations that are exactly satisfied. The first one is found by steepest descent. After the first iteration, we do plane searches using the weighted median solver to choose the best equation to be exactly satisfied. ``Best equation'' is the one that decreases the residual the most, while satisfying the equation chosen by the previous iteration exactly as well.

To do this, suppose the previous equation is $g^{(d)}_k \alpha +s^{(d)}_k \beta + r^{(d)}_k = 0$. We seek a $(\Delta\alpha, \Delta\beta)$ that still satisfies the equation $k$. This requirement gives a solution to $(\Delta\alpha, \Delta\beta)$ as $\gamma (s^{(d)}_k, -g^{(d)}_k)$, where $\gamma$ is a scalar. Then the plane search becomes $\bold 0 \approx \gamma \bold B [s^{(d)}_k, -g^{(d)}_k]' + r^{(d)}$, which is a weighted median problem. Thus, using the weighted median solver, we can solve for $\gamma$ and get a new equation, e.g., equation $j$. Then we can update $(\alpha, \beta)$ and $\bold s^{(d)}$ accordingly, drop the old equation $k$, and keep the new equation $j$. We iterate on this process until the inner loop keeps tracking the same equation. The final results of the inner loop are passed out to update the model, residual and the gradient.

The value of expending more effort to find the best step direction will be supported by the real geophysical applications, because the most computationally expensive part of these iterative methods is applying the forward and adjoint operators (steps starred in Table 1). By doing the sophisticate plane search, we hope to decrease the number of outer-loop iterations required for convergence.

However, conjugate direction L1 regression theory is not perfect for a practical problem. The problem of a flat bottom in $L1$ minimization will cause trouble in geophysical practice. Sometimes even where the bottom is not exactly as flat as the median of an even number of points, the slope of the gradient can be so small that we can never reach convergence in a finite number of iterations.

Dix inversion constrained by L1-norm optimization