\documentstyle[11pt]{article}

\def\be{\begin{equation}}
\def\bel#1{\begin{equation}\label{eq:#1}}
\def\ee{\end{equation}}
\def\bea{\begin{eqnarray}}
\def\beal#1{\begin{eqnarray}\label{eq:#1}}
\def\eea{\end{eqnarray}}

\def\bfA{{\bf A}}
\def\bfC{{\bf C}}
\def\bfS{{\bf S}}
\def\bfN{{\bf N}}
\def\bfM{{\bf M}}
\def\bfE{{\bf E}}
\def\bfV{{\bf V}}
\def\bfy{{\bf y}}
\def\sstar{{s_\ast}}
\def\sstarhat{{{\widehat s}_\ast}}
\def\sigmastar{{{\widehat \sigma}_\ast}}
\def\varstar{{{\widehat \sigma}^2_\ast}}
\def\ybar{{\bar y}}


\begin{document}

\title{Some Notes on Wiener Reconstruction}
\author{George B. Rybicki}
\date{10 Oct 1997}

\maketitle

\section{Stationary Processes}

Wiener reconstruction of a curve from one-dimensional stationary data
$y(t_i)$ at a discrete set of points $t_i$, $i=1,2,\ldots,n$ is
treated in Press, Rybicki, \& Hewitt (\cite{PRH}; PRH) and Rybicki \&
Press (\cite{RP}; RP).  The basic results as given in RP are
\bel{1}
        \sstarhat = \bfS^T_\ast [\bfS+\bfN]^{-1} \bfy
\ee
for the reconstructed value [RP eq.\ (8)], and
\bel{2}
  \varstar=\langle \sstar^2 \rangle -\bfS^T_\ast [\bfS+\bfN]^{-1} \bfS_\ast
\ee
for the variance about this value [RP eq.\ (9)], that is,
$\varstar=\langle (\sstar -\sstarhat)^2 \rangle$.
The notation
is described in detail in RP, but briefly: $\bfy$ is the vector of
measued values; $\bfS$ and $\bfN$ are the correlation matrices of
signal and noise; $\bfS_\ast$ is the correlation vector between
the signal at the reconstructed point and the measured points;
and $\langle \sstar^2 \rangle$ is the variance of the signal at the
reconstructed point.

Using the above quantities, one can find the entire Wiener
reconstruction of the curve by evaluating eq.\ (\ref{eq:1}) for each
value of $t$ in the desired range.  Furthermore, one obtains error
estimates for the reconstruction from eq.\ (\ref{eq:2}).  These may be
graphically presented as a ``snake,'' the region between the two
curves $\sstarhat \pm \sigmastar$.
           
\section{Weakly Nonstationary Processes}

It is also of interest to consider weakly nonstationary processes, for
which the mean is not well defined.  As discussed in PRH and RP the
revelant formulas for this case can be derived by assuming that the
correlation function has an essentially infinite constant part.  For
the purpose of many numerical calculations, one can achieve this by
simply adding a sufficiently large constant to the correlation
function and and then using eqs.\ (\ref{eq:1}) and (\ref{eq:2}).
However, it is of some interest to have special formulas for this case
in which the formal limit has been taken analytically.

The detailed derivation is accomplished by first adding a finite
constant $\lambda$ to the correlation function; this is equivalent to
making the replacements
\beal{3}
      \bfS &\rightarrow & \bfS + \lambda \bfE \bfE^T \nonumber\\
      \bfS_\ast &\rightarrow & \bfS_\ast + \lambda \bfE \nonumber\\
  \langle \sstar^2 \rangle &\rightarrow &  \langle \sstar^2 \rangle +\lambda
\eea
where $\bfE$ is a vector of ones.  The desired formulas are now
obtained by taking the limit $\lambda \rightarrow \infty$.

After performing these operations, we find (see appendix) that formula
(\ref{eq:1}) becomes
\bel{4}
  \sstarhat = \bfS^T_\ast [\bfS+\bfN]^{-1} (\bfy-\ybar \bfE) +\ybar
\ee
where
\bel{5}
         \ybar={\bfE^T [\bfS+\bfN]^{-1} \bfy \over
                 \bfE^T [\bfS+\bfN]^{-1} \bfE}
\ee
This is the result given in RP [RP eqs.\ (19) and (20)].

However, the weakly nonstationary result corresponding to formula
(\ref{eq:2}) is not given in RP.  It can be shown (see the appendix)
that the desired formula is
\bel{6}
                \varstar= \langle \sstar^2 \rangle
    -\bfS^T_\ast [\bfS+\bfN]^{-1} \bfS_\ast
  + {\left(\bfS_\ast^T [\bfS+\bfN]^{-1} \bfE -1 \right)^2
     \over \bfE^T [\bfS+\bfN]^{-1} \bfE}
\ee

It is clear that the formulas (\ref{eq:4}--\ref{eq:6}) for the
weakly nonstationary case must be invariant to
the substitutions (\ref{eq:3}) for any finite $\lambda$,
since such a substitution cannot change the limit.
(This can also be proved directly, but the proof is lengthy.)  
In particular, one may choose $\lambda$ to be equal to
$-\langle \sstar^2 \rangle$, so that,
\beal{7}
      \bfS &\rightarrow & -\bfV \nonumber\\
      \bfS_\ast &\rightarrow & -\bfV_\ast \nonumber\\
  \langle \sstar^2 \rangle &\rightarrow & \phantom{-} 0
\eea
where $\bfV$ is the structure function matrix, etc.  Then eqs.\ 
(\ref{eq:4})--(\ref{eq:6}) become
\bel{8}
   \sstarhat= -\bfV^T_\ast [-\bfV+\bfN]^{-1} (\bfy-\ybar \bfE) +\ybar
\ee
where
\bel{9}
         \ybar={\bfE^T [-\bfV+\bfN]^{-1} \bfy \over
                 \bfE^T [-\bfV+\bfN]^{-1} \bfE}
\ee
and
\bel{10}
   \varstar= -\bfV^T_\ast [-\bfV+\bfN]^{-1} \bfV_\ast  
   +{\left(\bfV_\ast^T [-\bfV+\bfN]^{-1} \bfE +1 \right)^2
     \over \bfE^T [-\bfV+\bfN]^{-1} \bfE}
\ee
Thus the calculation of the Wiener reconstruction for weakly
nonstationary processes can be done using only the structure function.


\appendix

\section{Appendix: Proof of Equations (4) and (6)}

It is convenient here to define matrices $\bfA$ and $\bfC$ by
$\bfC=\bfA^{-1}=\bfS+\bfN$.  After making the substitutions
(\ref{eq:3}) the quantities $\sstarhat$ and $\varstar$
become
\bel{11}
 \sstarhat = [\bfS^T_\ast+\lambda \bfE^T] [\bfC +\lambda \bfE \bfE^T]^{-1} \bfy
\ee
and
\bel{12}
  \varstar=[\langle \sstar^2 \rangle +\lambda ]
   -[\bfS^T_\ast+\lambda \bfE^T] [\bfC +\lambda \bfE \bfE^T]^{-1}
    [\bfS_\ast+\lambda \bfE]  
\ee
The quantity $ [\bfC +\lambda \bfE \bfE^T]^{-1}$ can be expressed
by means of the well-known Woodbury formula
\bel{13}
    [\bfC +\lambda \bfE \bfE^T]^{-1}
 = \bfA  - {\lambda \over 1+\lambda g^{-1}} \bfA\bfE\bfE^T\bfA
 = \bfA  - {g \over 1+g\lambda^{-1}} \bfA\bfE\bfE^T\bfA
\ee
where
\bel{13.5}
         g \equiv ( \bfE^T \bfA \bfE)^{-1} 
\ee
is a scalar.  Equation (\ref{eq:13}) can now be
expanded to second order in the small parameter $\lambda^{-1}$ as follows:
\bel{14}
    [\bfC +\lambda \bfE \bfE^T]^{-1} = \bfM_0 + \lambda^{-1}\bfM_1
              + \lambda^{-2} \bfM_2 + \ldots
\ee
where
\bel{15}
         \bfM_0 = \bfA - g \bfA \bfE \bfE^T \bfA,\qquad
         \bfM_i = (-g)^{i+1}\bfA \bfE \bfE^T \bfA, \qquad  i\ge 1
\ee
We now expand the following quantity, which appears in the both
formulas (\ref{eq:11}) and (\ref{eq:12}),
\bel{16}
[\bfS^T_\ast+\lambda \bfE^T] [\bfC +\lambda \bfE \bfE^T]^{-1}  
= (\bfE^T\bfM_1 + \bfS^T_\ast \bfM_0) 
   +\lambda^{-1}(\bfE^T\bfM_2 + \bfS^T_\ast \bfM_1) +O(\lambda^{-2})
\ee
where we have use the results $\bfE^T \bfM_0 = \bfM_0 \bfE=0$.
Therefore as $\lambda \rightarrow \infty$,
\bel{17}
  \sstarhat = (\bfE^T\bfM_1 + \bfS^T_\ast \bfM_0) \bfy
\ee
Substituting for $\bfM_i$, using the definition of $g$, we find
\beal{18}
  \sstarhat &=& g\bfE^T \bfA \bfy + \bfS^T_\ast \bfA \bfy
                 -g\bfS^T_\ast \bfA \bfE \bfE^T \bfy \nonumber\\
       &=& \ybar + \bfS^T_\ast \bfA y - \bfS^T_\ast \bfA \ybar \bfE \nonumber\\
       &=& \bfS^T_\ast \bfA (\bfy - \ybar \bfE) +\ybar
\eea
which proves eq.\ (\ref{eq:4}).  

Returning to eq.\ (\ref{eq:16}), we find that
\beal{19}
&& [\bfS^T_\ast+\lambda \bfE^T][\bfC +\lambda \bfE \bfE^T]^{-1} 
  [\bfS_\ast+\lambda \bfE]  = \lambda (\bfE^T \bfM_1 \bfE) \\
&&{\hskip 1.0in}+(\bfE^T \bfM_2 \bfE + \bfE^T \bfM_1 \bfS_\ast 
          +\bfS^T_\ast \bfM_1 \bfE +\bfS^T_\ast\bfM_0 \bfS_\ast)+O(\lambda^{-1})
       \nonumber
\eea
Substitution into eq.\ (\ref{eq:12}), we note that the term of order $\lambda$
cancels due to $\bfE^T \bfM_1 \bfE = 1$.  Thus, in the limit $\lambda 
\rightarrow \infty$,
\beal{20}
  \varstar &=& \langle \sstar^2 \rangle -\bfS^T_\ast \bfA \bfS_\ast
   +g-g\bfE^T \bfA \bfS_\ast -g\bfS^T_\ast \bfA \bfE 
   +\bfS^T_\ast \bfA \bfE \bfE^T \bfS_\ast \nonumber\\
    &=& \langle \sstar^2 \rangle -\bfS^T_\ast \bfA \bfS_\ast
          +g(\bfS^T_\ast \bfA \bfE -1)^2 
\eea
which proves eq.\ (\ref{eq:6}).  [Note that to obtain this expression
the expansion of Eq.\ (\ref{eq:14}) was required to second order.]

\begin{thebibliography}{}
\bibitem[1992]{PRH}
Press, W.H., Rybicki, G.B., Hewitt, J.N. 1992, ApJ, 385, 404
\bibitem[1992]{RP}
Rybicki, G.B., Press, W.H. 1992, Astrophys. Journ. 398, 169


\end{thebibliography}
\end{document}


From the definitions of $g$ and the $\bfM_i$, we note a number of 
simple relations such as
$\bfE^T\bfM_1 \bfE = 1$, 
$\bfM_1\bfE= g \bfA \bfE$, etc.  Using these we find



\beal{17}
  [\bfS^T_\ast+\lambda \bfE^T] [\bfC +\lambda \bfE \bfE^T]^{-1} 
    &=& \lambda^1
  +  [-g + g\bfE^T \bfA \bfS_\ast 
          +g\bfS^T_\ast \bfA \bfE +\bfS^T_\ast\bfM_0 \bfS_\ast] \nonumber\\
    &=& \lambda + +\bfS^T_\ast\bfC \bfS_\ast 
          +\left( \bfS^T_\ast \bfA \bfE -1 \right)^2
\eea
