TWO-STAGE LINEAR LEAST SQUARES

To fill empty bins by optimization, we need a filter. That filter should be a prediction-error filter (PEF). Thus, finding missing data is a two-stage process using linear least squares (LS) at each stage.

Our subroutines for finding the PEF would have the pitfall that they would use the missing data to find the PEF, and hence they would get the wrong PEF. We can use a PEF finder that uses weighted least squares. Then we can use a weighting function that vanishes on those fitting equations that involve missing data (and is unity elsewhere). The hconest module is a convolution module that we will use. Instead of multiplying a convolution output by a zero, it simply does not compute the convolution output.  

module hconest {             # masked helix convolution, adjoint is the filter.
  real,    dimension( :), pointer :: x
  logical, dimension( :), pointer :: bad	# output depends on bad inputs
  integer, dimension( :), pointer :: lag
#% _init( bad, x, lag)
#% _lop( a, y)
    integer  ia, ix, iy
    do ia = 1, size( a) {
         do iy = 1  + lag( ia), size( y) {    if( bad( iy)) cycle  
            ix = iy - lag( ia)
             if( adj) 
                    a( ia) = a( ia) + y( iy) * x( ix)
             else
                    y( iy) = y( iy) + a( ia) * x( ix)
         }
    }
}

Suppose that y₄ were a missing or a bad data value in the fitting goal (28). That would spoil the 4th, 5th and 6th fitting equations. Thus we would want to be sure that w₄, w₅, and w₆ were zero.  

$$\bold 0 \ \approx\ \bold W \bold r \ =\ \left[ \begin{arra... ... \begin{array} {c} 1 \\ a_1 \\ a_2 \end{array} \right]$$ (28)

When missing y_i are sprinkled around, which of the w_i should vanish? A simple solution (that we cannot use) is to put zeros in place of missing data and then take the product of every element in a row of the $\bold Y$ matrix. The product will vanish if any element in the row vanishes. We cannot use this solution because it would require us to have the matrix, when actually we have only an operator, a program pair that applies the matrix. A useful solution uses complementary logic: We prepare a template like $\bold y$ with ones where data is missing and zeros where the data is known. We prepare a template for $\bold a$ where all values are ones. We put the templates into (28) and apply the operator. Where the outputs are nonzero we have defective fitting equations and that is where we want the weights to be zero.