Packing both missing data and filter into a vector

missing data Now let us examine the theory and coding behind the above examples. Define a roughening filter $A(\omega )$ and a data signal $Y(\omega )$at some stage of interpolation. The fitting goal is $0 \approx A(\omega ) Y(\omega )$where the filter $A(\omega )$ has at least one time-domain coefficient constrained to be nonzero and the data contains both known and missing values. Think of perturbations $\Delta A$ and $\Delta Y$.We neglect the nonlinear term $\Delta A\,\Delta Y$ as follows:

$$\begin{eqnarray} 0 &\approx & (A \ +\ \Delta A)( Y\ +\ \Delta Y) \\ 0 &\approx &... ...lta Y \\ 0 &\approx & A\,\Delta Y \ +\ Y\,\Delta A \ +\ AY\end{eqnarray}$$ (37) (38) (39)

Let us use matrix algebraic notation to rewrite the fitting goals (39). For this we need mask matrices mask matrix (diagonal matrices with ones on the diagonal where variables are free and zeros where they are constrained i.e., where $\Delta a_i=0$ and $\Delta y_i=0$). The free-mask matrix for missing data is denoted $\bold J$and that for the PE filter is $\bold K$.The fitting goal (39) becomes  

$$\bold 0 \quad \approx \quad \bold r \eq \bold A \bold J \Del... ...ta \bold a + ( \bold A \bold y {\rm \ or \ } \bold Y \bold a )$$ (40)

Defining the original residual as $\bar {\bold r} = \bold A\bold y$ this becomes

$$\bold 0 \quad\approx\quad \bold r \eq \left[ \begin{array... ... y \\ \Delta \bold a \end{array} \right] \ +\ \bar {\bold r}$$ (41)

For a 3-term filter and a 7-point data signal, the fitting goal (40) becomes  

$$\left[ \begin{array} {cccccccccc} a_0& . & . & . & . & . &... ...r r_7 \\ \bar r_8 \end{array} \right] \quad \approx \ \bold 0$$ (42)

Recall that $\bar r_t$ is the convolution of a_t with y_t, namely, $\bar r_0=y_0 a_0$ and $\bar r_1=y_0 a_1 + y_1 a_0$, etc. To optimize this fitting goal we first initialize $\bold a= (1,0,0,\cdots )$and then put zeros in for missing data in $\bold y$.Then we iterate over equations (43) to (47).  

$$\bold r \quad\longleftarrow\quad \bold A \bold y$$ (43)

$$\left[ \begin{array} {c} \Delta \bold y \\ \Delta \bold a... ...\bold A' \\ \bold K' \bold Y' \end{array} \right] \ \bold r$$ (44)

$$\Delta \bold r \quad\longleftarrow\quad \left[ \begin{arr... ...ay} {c} \Delta \bold y \\ \Delta \bold a \end{array} \right]$$ (45)

$$\begin{eqnarray} \bold y &\longleftarrow& {\rm cgstep}( \bold y, \Delta \bold y) \\ \bold a &\longleftarrow& {\rm cgstep}( \bold a, \Delta \bold a)\end{eqnarray}$$ (46) (47)

This is the same idea as all the linear fitting goals we have been solving, except that now we recompute the residual $\bold r$ inside the iteration loop so that as convergence is achieved (if it is achieved), the neglected nonlinear term $\Delta A \Delta Y$ tends to zero.

My initial research proceeded by linearization like (39). Although I ultimately succeeded, I had enough difficulties that I came to realize that linearization is dangerous. When you start ``far enough'' from the correct solution the term $\Delta A\,\Delta Y$might not actually be small enough. You don't know how small is small, because these are not scalars but operators. Then the solution may not converge to the minimum you want. Your solution will depend on where you start from. I no longer exhibit the nonlinear solver missif until I find a real data example where it produces noticibly better results than multistage linear-least squares.

The alternative to linearization is two-stage linear least squares. In the first stage you estimate the PEF; in the second you estimate the missing data. If need be, you can re-estimate the PEF using all the data both known and missing (downweighted if you prefer).

If you don't have enough regression equations because your data is irregularly distributed, then you can use binning. Still not enough? Try coarser bins. The point is that nonlinear solvers will not work unless you begin close enough to the solution, and the way to get close is by arranging first to solve a sensible (though approximate) linearized problem. Only as a last resort, after you have gotten as near as you can, should you use the nonlinear least-squares techniques.