Regularization is model styling

Operator approach to missing data

For the operator approach to the fitting goal $-\bold F_k \bold m_k \approx \bold F_u \bold m_u$, we rewrite it as $-\bold F_k \bold m_k \approx \bold F \bold J \bold m$ where

$$-\bold F_k \bold m_k \approx \ \left[ \begin{array}{cccc... ...m_5 \\ m_6 \end{array} \right] = \bold F \bold J \bold m$$ (3)

Notice the introduction of the new diagonal matrix $\bold J$, called a masking matrix or a constraint-mask matrix because it multiplies constrained variables by zero leaving freely adjustable variables untouched. Experience shows that a better name than ``mask matrix'' is ``selector matrix'' because what comes out of it, that which is selected, is a less-confusing name for it than which is rejected. With a selector matrix, the whole data space seems freely adjustable, both the missing data values and known values. We see that the CD method does not change the known (constrained) values. In general, we derive the fitting goal (3) by:

$$\bold 0 \approx \bold F \bold m$$ (4)

$$\bold 0 \approx \bold F( \bold J + (\bold I-\bold J) ) \bold m$$ (5)

$$\bold 0 \approx \bold F\bold J\bold m + \bold F(\bold I-\bold J)\bold m$$ (6)

$$\bold 0 \approx \bold F\bold J\bold m + \bold F\bold m_{\rm known}$$ (7)

$$\bold 0 \quad\approx\quad \bold r = \bold F\bold J\bold m + \bold r_0$$ (8)

As usual, we find a direction to go $\Delta \bold m$ by the gradient of the residual energy.

$$\Delta\bold m = \frac{\partial }{\partial \bold m\T} \bold ... ...T + \bold r\T_0) \right) \bold r = \bold J\T \bold F\T \bold r$$ (9)

We begin the calculation with the known data values where missing data values are replaced by zeros, namely $(\bold I-\bold J)\bold m$. Filter this data, getting $\bold F(\bold I-\bold J)\bold m$, and load it into the residual $\bold r_0$. With this initialization completed, we begin an iteration loop. First we compute $\Delta \bold m$ from equation (9).

$$\Delta\bold m \quad\longleftarrow\quad \bold J\T \bold F\T \bold r$$ (10)

$\bold F\T$ applies a crosscorrelation of the filter to the residual and then $\bold J\T$ sets to zero any changes proposed to known data values. Next, compute the change in residual $\Delta\bold r$ from the proposed change in the data $\Delta \bold m$.

$$\Delta\bold r \quad \longleftarrow \quad \bold F \bold J \Delta \bold m$$ (11)

Equation (11) applies the filtering again. Then, use the method of steepest descent (or conjugate direction) to choose the appropriate scaling (or inclusion of previous step) of $\Delta \bold m$ and $\Delta\bold r$, and update $\bold m$ and $\bold r$, accordingly, and iterate.

The subroutine to find missing data is mis1(). It assumes that zero values in the input data correspond to missing data locations. It uses our convolution operator tcai1(). You can also check the book Index for other operators and modules.

user/gee/mis1.c

void mis1(int niter         /* number of iterations */, 
	  float *xx         /* data/model */, 
	  const bool *known /* mask for known data */,
	  const char *step  /* solver */) 
/*< interpolate >*/
{
    switch (step[1]) {
	case 'g': /* conjugate gradients */
	    sf_solver (tcai1_lop, sf_cgstep, nx, ny, xx, zero, 
		       niter, "x0", xx, "known", known, "end");
	    sf_cgstep_close();
	    break;
	case 'd': /* conjugate directions */
	    sf_cdstep_init();
	    sf_solver (tcai1_lop, sf_cdstep, nx, ny, xx, zero, 
		       niter, "x0", xx, "known", known, "end");
	    sf_cdstep_close();
	    break;
	default:
	    sf_error("%s: unknown step %s",__FILE__,step);
	    break;
    }
}

I sought reference material on conjugate gradients with constraints and did not find anything, leaving me to fear that this chapter was in error, and I had lost the magic property of convergence in a finite number of iterations. I tested the code, and it did converge in a finite number of iterations. The explanation is that these constraints are almost trivial. We pretended we had extra variables, and computed a $\Delta \bold m =\bold g$ for each. Using $\bold J$ we then set the gradient $\Delta \bold m =\bold g$ to zero, therefore making no changes to anything, like as if we had never calculated the extra $\Delta \bold m$s.

Regularization is model styling