Blocky models via the L1/L2 hybrid norm

BLOCKY LOGS: BOTH FITTING AND REGULARIZATION

Here we set out to find blocky functions, such as well logs. We will do data fitting with a somewhat $L2$-like convex penalty function while doing model styling with a more $L1$-like function. We might define the composite norm threshold residual $R_d$ for the data fitting at the 60th percentile, and that for regularization seeking spiky models (with blocky integrals) as $R_m$ at the 5th percentile.

The data fitting goal and the model regularization goal at each $z$ is independent from that at all other $z$ values. The fitting goal says the reflectivity $m(z)$ should be equal to its measurement $d(z)$ (the seismogram). The model styling goal says the reflectivity $m(z)$ should vanish.

$$\approx r_d(z) = m(z) - d(z)$$ 0 (25)

$$\approx r_m(z) = \epsilon m(z)$$ 0 (26)

These two goals are in direct contradiction to each other. With the L2 norm the answer would be simply $m=d/(1+\epsilon^2)$. With the L1 norm, the answer would be either $m=d$ or $m=0$ depending on the numerical choice of $\epsilon$. Let us denote the convex function and its derivatives for data space at the residual as $(B,B',B'')$ and for model space as $(C,C',C'')$. Remember, $m$ and $d$, while normally vectors, are here scalars (independently for each $z$).

loop over all time points {
$m=d/(1+\epsilon^2)$         # These are scalars!
loop over non-linear iterations {
$r_d = m-d$
$r_m = m$
Get derivatives of hybrid norm $B'(r_d)$ and $B''(r_d)$ for data goal.
Get derivatives of hybrid norm $C'(r_m)$ and $C''(r_m)$ for model goal.
# Plan to find $\alpha$ to update $m=m+\alpha$
# Taylor series for data penalty $N(r_d) = B + B'\alpha + B''\alpha^2/2$
# Taylor series for model penalty $N(r_m) = C + C'\alpha + C''\alpha^2/2$
# $0 = \frac{\partial}{\partial \alpha} (N(r_d) + \epsilon N(r_m))$
$\alpha = - (B'+\epsilon C')/(B''+\epsilon C'')$
$m=m+\alpha$
} end of loop over non-linear iterations
} end of loop over all time points

To help us understand the choice of parameters $R_d$, $R_m$, and $\epsilon$, We examine the theoretical relation between $m$ and $d$ implied by the above code as a function of $\epsilon$ and $R_m$ at $R_d\rightarrow\infty$, in other words, when the data has normal behavior and we are mostly interested in the role of the regularization drawing weak signals down towards zero. The data fitting penalty is $B=(m-d)^2/2$ and its derivative $B'=m-d$. The derivative of the model penalty (from equation (15)) is $C'=m/\sqrt{1+m^2/R_m^2}$. Setting the sum of the derivatives to zero we have

$$0 \quad=\quad B'+\epsilon C' \quad=\quad m - d +\frac{\epsilon m}{\sqrt{1+m^2/R_m^2}}$$ (27)

This says $m$ is mostly a little smaller than $d$, but it gets more interesting near $(m,d)\approx 0$. There the slope $m/d = 1/(1+\epsilon)$ which says an $\epsilon=4$ will damp the signal (where small) by a factor of 5. Moving away from $m=0$ we see the damping power of $\epsilon$ diminishes uniformly as $m$ exceeds $R_m$.

Blocky models via the L1/L2 hybrid norm