Operator approach to missing data

For the operator approach to the fitting goal $-\bold A_k \bold m_k \approx \bold A_u \bold m_u$we rewrite it as $-\bold A_k \bold m_k \approx \bold A \bold J \bold m$ where

 

$$-\bold A_k \bold m_k \quad \approx \quad \left[ \begin{arra... ... m_6 \end{array} \right] \quad =\quad \bold A \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. In this new formulation, the whole data space seems freely adjustable, both the missing data values and known values, but we will soon see that the CD method will not be changing the known (constrained) values. In general, we derive the fitting goal (3) by

$$\begin{eqnarray} \bold 0 &\approx& \bold A \bold m \\ \bold 0 &\approx& \bold A... ...0 \quad\approx\quad \bold r &=& \bold A\bold J\bold m + \bold r_0\end{eqnarray}$$ (4) (5) (6) (7) (8)

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

$$\Delta\bold m \eq {\partial \over\partial \bold m'}\ \bold r... ... A' + \bold r'_0) \right) \bold r \eq \bold J' \bold A' \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 A(\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 pass the residual through the adjoint

$$\Delta\bold m \quad\longleftarrow\quad \bold J' \bold A' \bold r$$ (10)

$\bold A'$ applies a crosscorrelation of the filter to the residual and then $\bold J'$ 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 A \bold J \Delta \bold m$$ (11)

This 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.

I could have passed a new operator $\bold A \bold J$into the old solver, but found it worthwhile to write a new, more powerful solver having built-in constraints. To introduce the masking operator $\bold J$ into the solver subroutine , I introduce an optional input parameter known, which contains a logical array of the model size. The values of known are .true. for the known data locations, and .false. for the unknown data. Two lines in the simple_solver module

stat = oper( T, F, g, rr)                          # dm = A' r
stat = oper( F, F, g, gg)                          # dr = A  dm

become three lines in the standard library module solver_mod. (Here I changed the comment lines to correspond to the notation of the missing-data application; $\Delta m$ is g and $\Delta r$ is gg.)

stat = oper( T, F, g, rr)                          # dm  = A' r
if ( present( known))  where(known)   g = 0.       # dm  = JJ' dm
stat = oper( F, F, g, gg)                          # dr =  A   dm

The new solver initializes its residual exactly as did the original solver, namely,

$$\bold r_0 \eq \bold A \bold m_0 \ -\ \bold d$$ (12)

Using the original solver, we could allow any $\bold m_0$,but to get the right behavior with constraints, we would need to set $\bold d =-\bold A\bold m_k$.Using the new solver, we specify $\bold d =\bold 0$, but we must specify a nonzero $\bold m_0$ that matches the known data. Mathematically this is

$$(\bold I-\bold J)(\bold m_0-\bold d_k)=\bold 0$$ (13)

In Fortran77, I simply multiplied by $\bold J$.In Fortran90 subroutines can have optional arguments. The expression if( present( known)) says, ``if the argument known, a logical array, is an argument of the call.'' The expression, where( known) dm = 0 means that each component of the logical array known is examined; if the value of that component is .true. then the corresponding component of $\Delta m$ is set to zero. In other words, we are not going to change the given, known data values.

The subroutine to find missing data is mis1(). It assumes that zero values in the input data correspond to missing data locations.  

module mis_mod {  
  use tcai1
  use cgstep_mod
  use solver_mod
contains
# fill in missing data on 1-axis by minimizing power out of a given filter.
  subroutine mis1 ( niter, mm, aa) {
    integer,             intent (in)      :: niter   # number of iterations
    real, dimension (:), pointer          :: aa      # roughening filter
    real, dimension (:), intent (in out)  :: mm      # in - data with zeroes
                                                     # out - interpolated
    real, dimension (size(mm) + size(aa)) :: zero    # filter output

call tcai1_init( aa, size( mm));         zero = 0.

call solver( tcai1_lop, cgstep, mm, zero, niter,
                              x0 =  mm,    known = (mm /= 0.) )
    call cgstep_close ()
    }
}

Here known is declared in solver_mod to be a logical vector. The call argument known=(xx/=0.) sets the mask vector components to .true. where components of $\bold x$ are nonzero and sets it to .false. where the components are zero.

I sought reference material on conjugate gradients with constraints and didn't find anything, leaving me to fear that this chapter was in error and that 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 m$ for each of them. Then we set each $\Delta m$ to zero, making no changes to anything, like as if we had never calculated the extra $\Delta m$'s.

EXERCISES:

1. Figure 1-4 seem to extrapolate to vanishing signals at the side boundaries. Why is that so, and what could be done to leave the sides unconstrained in that way?
2. Show that the interpolation curve in Figure 2 is not parabolic as it appears, but cubic. (HINT: First show that $(\nabla^2)'\nabla^2 u = \bold 0$.)
3. Verify by a program example that the number of iterations required with simple constraints is the number of free parameters. 1