INVERSE LINEAR INTERPOLATION

 

data Figure 4 The input data are irregularly sampled.

The first example is a simple synthetic test for 1-D inverse interpolation. The input data were randomly subsampled (with decreasing density) from a sinusoid (Figure 4). The forward operator in this case is linear interpolation. We seek a regularly sampled model that could predict the data with a forward linear interpolation. Sparse irregular distribution of the input data makes the regularization enforcement a necessity. I applied convolution with the simple (1,-1) difference filter as the operator that forces model continuity (the first-order spline). An appropriate preconditioner $\bold S$ in this case is recursive causal integration.

 

conv1
Figure 5 Convergence history of inverse linear interpolation. Left: regularization, right: preconditioning. The regularization operator $\bold A$ is the derivative operator (convolution with (1,-1). The preconditioning operator $\bold S$ is causal integration.

As expected, preconditioning provides a much faster rate of convergence. Since iteration to the exact solution is never achieved in large-scale problems, the results of iterative optimization may turn out quite differently. Bill Harlan points out that the two goals in (8) conflict with each other: the first one enforces ``details'' in the model, while the second one tries to smooth them out. Typically, regularized optimization creates a complicated model at early iterations. At first, the data fitting goal (8) plays a more important role. Later, the regularization goal (8) comes into play and simplifies (smooths) the model as much as needed. Preconditioning acts differently. The very first iterations create a simplified (smooth) model. Later, the data fitting goal adds more details into the model. If we stop the iterative process early, we end up with an insufficiently complex model, not in an insufficiently simplified one. Figure 5 provides a clear illustration of Harlan's observation.

Figure 6 measures the rate of convergence by the model residual, which is a distance from the current model to the final solution. It shows that preconditioning saves many iterations. Since the cost of each iteration for each method is roughly equal, the efficiency of preconditioning is evident.

 

schwab1 Figure 6 Convergence of the iterative optimization, measured in terms of the model residual. The ``p'' points stand for preconditioning; the ``r'' points, regularization.

The module invint2 invokes the solvers to make Figures 5 and 6. We use convolution with helicon for the regularization and we use deconvolution with polydiv for the preconditioning. The code looks fairly straightforward except for the oxymoron known=aa%mis.

 

#$=head1 NAME 
#$
#$invint2 - Inverse linear interpolation;
#$
#$=head1 SYNOPSIS
#$
#$C<call invint1(niter,coord,ord,o1,d1,mm,mmov,eps,aa,doprec)>
#$  
#$=head1 PARAMETERS        
#$  
#$=over 4 
#$  
#$=item niter - integer
#$  
#$      Number of itterations
#$
#$=item coord -  C<real(:)> 
#$     
#$      Coordinates
#$  
#$=item o1,d1 -  real
#$     
#$      First model position and sampling
#$
#$=item ord   -  C<real(:)>
#$
#$      Data values
#$
#$=item mm    -  C<real(:)>
#$
#$      Output model
#$
#$=item mmov  -  C<real(:,:)>
#$  
#$      Model movie
#$
#$=item eps   -  real
#$     
#$      Epsilon value if doing preconditioning
#$  
#$=item aa    -  type(filter)
#$
#$      Preconditioning operator
#$
#$=item prec  -  logical
#$  
#$      Whether or not to do preconditioning
#$  
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$Perform inverse linear interpolation
#$
#$
#$=head1 SEE ALSO
#$
#$L<lint1>,L<helicon>,L<polydiv>,L<solver>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
#$

module invint2 {   # Inverse linear interpolation
  use lint1
  use helicon                           # regularized    by helix   filtering
  use polydiv                           # preconditioned by inverse filtering
  use cgstep_mod + solver_mod
contains
  subroutine invint1( niter, coord,ord, o1,d1, mm,mmov, eps, aa, doprec) {
    logical,                  intent( in)  :: doprec
    integer,                  intent( in)  :: niter
    real,                     intent( in)  :: o1, d1, eps
    real,    dimension( :),   intent( in)  :: ord  
    type( filter),            intent( in)  :: aa
    real,    dimension( :),   intent( out) :: mm             
    real,    dimension( :,:), intent( out) :: mmov         # model movie
    real,    dimension( :),   pointer      :: coord        # coordinate
    call lint1_init( o1, d1, coord)
    if( doprec) {                                          # preconditioning
       call polydiv_init( size(mm), aa)
       call solver_prec( lint1_lop, cgstep, niter = niter, x = mm, dat = ord,  
                    prec = polydiv_lop, nprec = size(mm), eps = eps, 
                    xmov = mmov)
       call polydiv_close()
    } else {                                               # regularization
       call helicon_init( aa)
       call solver_reg(  lint1_lop, cgstep, niter = niter, x = mm, dat = ord, 
                    reg = helicon_lop, nreg = size(mm), eps = eps,
                    xmov = mmov)
    }                                               
    call cgstep_close()
  }
}