\chapter{Applications}
\label{ch:Applications}


\section{Signal Detection: Drifting Phase Sinusoids}
\label{sec:applications-detection}

There are various detection applications in underwater acoustics
and optical communications
where the problem can be formulated using the following two alternative
hypotheses:

\begin{eqnarray*}
H_0: ~ x_k &=& n_k
\\
H_1: ~ x_k &=& s_k + n_k
\\
&=& A~ cos(\omega k + \xi \theta_k + \phi) + n_k
\end{eqnarray*}
where $x_k$ is a scalar observation process, $n_k$ is white
Gaussian noise, $\theta_k$ is a unit Brownian motion\footnote{
$\theta_k$ is a zero mean Gaussian process with autocorrelation
function $E\{ \theta_k \theta_m \} = \min \{k,m\}$.}
and $\phi$ is an independent random initial phase uniformly distributed
in $[0, 2\pi]$. Without loss of generality we can assume that
$n_k$ has unit variance. The above model is suitable for the detection
problem of weak acoustic signals which require relatively long observation
intervals.
Although the signal might be nominally harmonic,
during such long intervals the phase
of the signal cannot be assumed constant \cite{Veeravali91}.
Another example where this problem formulation arises
is in optical communications where $\theta_k$ models the laser phase noise
\cite{Foschini88}.
The parameters $A$ and $\omega$
denote the amplitude and the frequency of the sinusoid respectively,
and $\xi$ is called the {\em bandwidth} of the phase drift.
Obviously when $\xi=0$ then the problem is reduced to
detecting sinusoids with random but constant phase immersed
in white Gaussian noise. It is known that in this case the quadrature
(standard noncoherent) detector is optimal \cite{Whalen71}.
When $\xi>0$ however the optimal Neyman-Pearson detector is different
and difficult to calculate analytically. Given an observation
sequence ${\bf x}=[x_1 \ldots x_N]$ the likelihood ratio is defined as

\begin{equation}
L({\bf x}) = \frac{p_1({\bf x})}{p_0({\bf x})}
\label{eq:likelihood_ratio}
\end{equation}
where $p_0(\cdot)$ and $p_1(\cdot)$ denote probability under hypotheses
$H_0$, $H_1$ respectively (shorthand for $Prob(\cdot | H_0)$,
$Prob(\cdot | H_1)$).
Under the Gaussian noise assumption we have

\begin{eqnarray}
p_0({\bf x}) &=& F~ exp\bigl\{ -1/2 \sum_{k=1}^N x_k^2 \bigr\}
\label{eq:p0_of_x}
\\
p_1({\bf x}|\bg{\theta},\phi) &=&
F~ exp\bigl\{ -1/2 \sum_{k=1}^N [x_k - A cos(\omega k + \xi \theta_k
+ \phi_k) ]^2 \bigr\} \nonumber
\\
p_1({\bf x}) &=& \int_{\bgsub{\theta}} \int_{\phi}
p_1({\bf x},\bg{\theta},\phi) = \nonumber
\\
&=& \int_{\bgsub{\theta}} \int_{\phi}
p_1({\bf x}|\bg{\theta},\phi) p(\bg{\theta}) p(\phi)
d\bg{\theta} d\phi
\label{eq:p1_of_x}
\end{eqnarray}
where $F$ is a normalizing constant and $\bg{\theta}=
[\theta_1 \ldots \theta_N]$ is the phase sequence.
Plugging Eqs.~\ref{eq:p0_of_x} and \ref{eq:p1_of_x} into
Eq.~\ref{eq:likelihood_ratio} we get

\begin{eqnarray}
L({\bf x}) &=& \int_{\bgsub{\theta}} \int_{\phi}
exp \left\{A \sum_{k=1}^N cos(\omega k + \xi \theta_k + \phi)
- \frac{1}{2} A^2 \sum_{k=1}^N cos^2(\omega k + \xi \theta_k + \phi) \right\}
\nonumber
\\
&& ~~~ p(\bg{\theta}) p(\phi) d\phi d\bg{\theta}
\end{eqnarray}

We assume $\omega N \gg 2\pi$\footnote{The frequency is assumed normalized
$0<\omega<2\pi$, so the condition in the text implies $N \gg 1$.}
so we approximate
\[
\sum_{k=1}^N cos^2(\omega k + \xi \theta_k + \phi) \approx N/2
\]
to obtain
\begin{eqnarray}
L({\bf x}) &=& exp\{ -A^2 N / 4 \}
\int_{\bgsub{\theta}} \int_{\phi=0}^{2 \pi}
exp\{ A \sum_{k=1}^N cos(\omega k + \xi \theta_k + \phi) \}
p(\bg{\theta}) p(\phi) d\phi d\bg{\theta} = \nonumber
\\
&=& exp\{ -A^2 N / 4 \} E_{\bgsub{\theta},\phi} \left\{
exp\bigl\{ A \sum_{k=1}^N cos(\omega k + \xi \theta_k + \phi) \bigl\}
\right\} =
\label{eq:likelihood_ratio_2}
\\
&\equiv& K \Lambda({\bf x})
\end{eqnarray}
where $K=exp\{ -A^2 N / 4 \}$
and $\Lambda({\bf x})$ is the expectation
in Eq.~\ref{eq:likelihood_ratio_2}.
Expanding $\Lambda$ into Taylor series we get

\begin{eqnarray}
\Lambda({\bf x}) &=& E_{\bgsub{\theta},\phi} \left\{
1 + \sum_{n=0}^{\infty} \frac{A^n}{n!} \left[
\sum_{k=1}^N x_k cos(\omega k + \xi \theta_k + \phi)
\right]^n \right\} = \nonumber
\\
&=& 1 + \sum_{n=0}^{\infty} \frac{A^n}{n!}
E_{\bgsub{\theta},\phi} \left\{
\left[ \sum_{k=1}^N x_k cos(\omega k + \xi \theta_k + \phi)
\right]^n \right\} = \nonumber
\\
&=& 1 + \sum_{n=0}^{\infty} \frac{A^n}{n!}
E_{\bgsub{\theta},\phi} \left\{
\left[ \sum_{k_1,\ldots,k_n=1}^N \prod_{i=1}^n x_{k_i}
cos(\omega k_i + \xi \theta_{k_i} + \phi)
\right] \right\} = \nonumber
\\
&=& 1 + \sum_{n=0}^{\infty} \frac{A^n}{n!}
\left[ \sum_{k_1,\ldots,k_n=1}^N \prod_{i=1}^n x_{k_i}
E_{\bgsub{\theta}} E_{\phi} \left\{ \prod_{i=1}^n
cos(\omega k_i + \xi \theta_{k_i} + \phi) \right\}
\right]
\label{eq:likelihood_ratio_expand}
\end{eqnarray}

Let us denote the quantity $E_{\bgsub{\theta}} E_{\phi} \left\{ \prod_{i=1}^n
cos(\omega k_i + \xi \theta_{k_i} + \phi) \right\}$ by $G_n$.
Notice that for $n=odd$, $G_n=0$ thus Eq.~\ref{eq:likelihood_ratio_expand}
becomes
\begin{equation}
\Lambda({\bf x}) = 1 + \sum_{n=even} \frac{A^n}{n!}
\left[ \sum_{k_1,\ldots,k_n=1}^N \prod_{i=1}^n x_{k_i}
G_n(k_1,\ldots,k_n) \right]
\end{equation}

The optimal Neyman-Pearson detector compares the likelihood ratio
with a threshold $\tau$, determined by the probability of false alarm
\begin{equation}
L({\bf x}) ~\lege~\tau
\end{equation}
or equivalently
\begin{equation}
\Lambda({\bf x}) ~\lege~\tau'
\label{eq:optimal_detector}
\end{equation}

Thus the optimal Neyman-Pearson detector contains only even-power terms in $x$
-- as expected since the signal has zero mean -- while the significance
of these terms varies in proportion to $A^n/n!$ therefore
very large terms have practically no effect.
As an example let us calculate the kernels of
the first two significant terms, $G_2$ and $G_4$:

$n=2$,
\begin{eqnarray}
G_2 &=& E_{\bgsub{\theta}} E_{\phi}\left\{
cos(\omega k_1 + \xi \theta_{k_1} + \phi)
cos(\omega k_2 + \xi \theta_{k_2} + \phi) \right\} = \nonumber
\\
&=& E_{\bgsub{\theta}} \left\{ \frac{1}{2 \pi} \frac{1}{2} \int_{\phi=0}^{2 \pi}
\left[ cos\bigl( \omega (k_1+k_2) + \xi (\theta_{k_1}+\theta_{k_2}) +
2 \phi \bigr) \right. \right. \nonumber
\\
&& \hs{0.5in} \left. \left.
+cos\bigl( \omega (k_1-k_2) + \xi (\theta_{k_1}-\theta_{k_2})
\bigr) \right] d\phi \right\} = \nonumber
\\
&=& \frac{1}{2} E_{\bgsub{\theta}} \left\{
cos\bigl( \omega (k_1-k_2) + \xi (\theta_{k_1}-\theta_{k_2}) \bigr)
\right\} = \nonumber
\\
&=& \frac{1}{2} E_{\bgsub{\theta}} \left\{
cos\bigl( \omega (k_1-k_2) \bigr)
cos\bigl( \xi (\theta_{k_1}-\theta_{k_2}) \bigr) -
sin\bigl( \omega (k_1-k_2) \bigr)
sin\bigl( \xi (\theta_{k_1}-\theta_{k_2}) \bigr)
\right\} = \nonumber
\\
&=& \frac{1}{2} cos\bigl( \omega (k_1-k_2) \bigr)
E_{\bgsub{\theta}} \left\{
cos\bigl( \xi (\theta_{k_1}-\theta_{k_2}) \bigr)
\right\} = \nonumber
\\
&=& \frac{1}{2} cos\bigl( \omega (k_1-k_2) \bigr)
e^{-\xi |\theta_{k_1}-\theta_{k_2}|}
\label{eq:likelihood_second_term}
\end{eqnarray}

$n=4$,
\begin{eqnarray}
G_4 &=& E_{\bgsub{\theta}} E_{\phi}
\left\{ cos(\omega k_1 + \xi \theta_{k_1} + \phi)
cos(\omega k_2 + \xi \theta_{k_2} + \phi) \right. \nonumber
\\
&& \hs{0.3in} \left.
cos(\omega k_3 + \xi \theta_{k_3} + \phi)
cos(\omega k_4 + \xi \theta_{k_4} + \phi) \right\} = \nonumber
\\
&=& \frac{1}{4} E_{\bgsub{\theta}} E_{\phi} \left\{
\left[ cos\bigl( \omega (k_1+k_2) + \xi (\theta_{k_1}+\theta_{k_2}) +
2 \phi \bigr) + cos\bigl( \omega (k_1+k_2) + \xi (\theta_{k_1}-\theta_{k_2})
\bigr) \right] \right. \nonumber
\\
&& \hs{0.3in} \left. 
\left[ cos\bigl( \omega (k_3+k_4) + \xi (\theta_{k_3}+\theta_{k_4}) +
2 \phi \bigr) + cos\bigl( \omega (k_3+k_4) + \xi (\theta_{k_3}-\theta_{k_4})
\bigr) \right] \right\} = \nonumber
\\
&=& \frac{1}{4} E_{\bgsub{\theta}} E_{\phi} \left\{
\left[ cos\bigl( \omega (k_1+k_2) + \xi (\theta_{k_1}+\theta_{k_2}) +
2 \phi \bigr) + cos\bigl( \omega (k_1+k_2) + \xi (\theta_{k_1}-\theta_{k_2})
\bigr) \right] \right. \nonumber
\\
&& \hs{0.3in} \left.
\left[ cos\bigl( \omega (k_3+k_4) + \xi (\theta_{k_3}+\theta_{k_4}) +
2 \phi \bigr) + cos\bigl( \omega (k_3+k_4) + \xi (\theta_{k_3}-\theta_{k_4})
\bigr) \right] \right\} = \nonumber
\\
&=& \frac{1}{4} E_{\bgsub{\theta}} E_{\phi} \left\{
cos\bigl( \omega (k_1+k_2) + \xi (\theta_{k_1}+\theta_{k_2}) + 2 \phi \bigr)
cos\bigl( \omega (k_3+k_4) + \xi (\theta_{k_3}+\theta_{k_4}) + 2 \phi \bigr)
\right. \nonumber
\\
&& \hs{0.3in} \left.
+ cos\bigl( \omega (k_1-k_2) + \xi (\theta_{k_1}-\theta_{k_2}) \bigr)
cos\bigl( \omega (k_3-k_4) + \xi (\theta_{k_3}-\theta_{k_4}) \bigr)
\right\} = \nonumber
\\
&=& \frac{1}{8} E_{\bgsub{\theta}} \left\{
cos\bigl( \omega (k_1+k_2-k_3-k_4) +
\xi (\theta_{k_1}+\theta_{k_2}-\theta_{k_3}-\theta_{k_4}) \bigr)
+ \right. \nonumber
\\
&& \hs{0.3in} \left.
sin\bigl( \omega (k_1-k_2-k_3+k_4) +
\xi (\theta_{k_1}-\theta_{k_2}+\theta_{k_3}-\theta_{k_4}) \bigr)
\right\}
\label{eq:likelihood_fourth_term}
\end{eqnarray}

It is interesting to notice that the term corresponding to $G_2$
in Eq.~\ref{eq:likelihood_ratio_expand} can be written as follows:

\begin{equation}
\sum_{k,m} x_k G_2 x_m = \frac{1}{2}
\sum_{k,m} x_k cos\bigl(\omega (k-m)\bigr) exp\bigl(-\xi |k-m|\bigr) x_m
= \frac{1}{2} {\bf x}^T {\bf R}_s {\bf x}
\label{eq:likelihood_second_term_matrix}
\end{equation}
where ${\bf R}_s$ is the autocorrelation matrix of $s_k$
with entries $r_{ij}=cos\bigl(\omega (i-j)\bigr) exp\bigl(-\xi |i-j|\bigr)$.
In \cite{Picinbono88} it was shown that
Eq.~\ref{eq:likelihood_second_term_matrix} is
the optimal linear-quadratic (second-order) test statistic for this problem.

The difficulty of calculating the likelihood ratio in
Eq.~\ref{eq:likelihood_ratio} makes the optimal detector
in Eq.~\ref{eq:optimal_detector} impractical to use.
Even the approximation \ref{eq:likelihood_second_term_matrix}
is very computationally costly and has the additional disadvantage
that it is difficult to implement in an analog circuit
which could alleviate the computational cost problem.
Various authors have investigated alternatives to the optimal solution.
In \cite{Georghiades85} the approach of estimating $\theta_k$
via an Extended Kalman Filter is considered. Once $\theta_k$ is
approximated, a coherent detector is used for testing the
two hypotheses. In \cite{Veeravali91} the authors
compare the performance of the
$m$-th order noncoherent detector (proposed in \cite{Foschini88})
which uses the test statistic

\begin{equation}
\Lambda_m({\bf x}) = \sum_{i=0}^{m-1} \left[ \left(
\sum_{k=\frac{Ni}{m}+1}^{\frac{N(i+1)}{m}} cos(\omega k) x_k \right)^2 +
\left( \sum_{k=\frac{Ni}{m}+1}^{\frac{N(i+1)}{m}} sin(\omega k) x_k \right)^2
\right]
\end{equation}
with the standard noncoherent detector (i.e. the quadrature detector)
and the optimal quadratic detector
(Eq.~\ref{eq:likelihood_second_term_matrix}).
The authors derived the optimal value of $m$ for
the $m$-th order detector and found that
for this value, it has comparable performance
to the optimal quadratic detector and is quite a lot simpler to implement,
while it is significantly better than the standard noncoherent detector.

In this section we propose to use principal component analysis
for the detection problem discussed above. The idea stems from the
fact that in the coherent case, all eigenvalues of ${\bf R}_s$,
except two, are equal to 0 while the principal component subspace
corresponding to the non-zero ones is spanned by the
quadrature signals $cos(\omega k)$ and $sin(\omega k)$.
The optimal detector then projects the observation signal
onto the optimal subspace following with
a projection-energy threshold test.
In the drifting phase case the signal spectrum has a nonzero bandwidth
proportional to $\xi$ and thus the eigenvalues of ${\bf R}_s$
are all nonzero in general. However, if $\xi$ is not too large
then the first two eigenvalues clearly dominate, so in the
extreme case where $\xi=0$ (coherent case)
only the first two eigenvalues are nonzero as discussed above.
When $\xi$ is non-zero however the space spanned
by $cos$ and $sin$ is only as a rough approximation
to the 2-dimensional principal component subspace.
This motivates us to compare
the performance of the $m$-th order noncoherent detector
with the $m$-order detector produced after substituting
$cos$ and $sin$ with the first two eigenvectors ${\bf e}_1$
and ${\bf e}_2$:

\begin{equation}
\Lambda_m'({\bf x}) = \sum_{i=0}^{m-1} \left[ \left(
\sum_{k=\frac{Ni}{m}+1}^{\frac{N(i+1)}{m}} e_{1k} x_k \right)^2 +
\left( \sum_{k=\frac{Ni}{m}+1}^{\frac{N(i+1)}{m}} e_{2k} x_k \right)^2
\right]
\label{eq:APEX_detector}
\end{equation}
In our comparison experiments we used the same value of $m$
as proposed in \cite{Veeravali91} for both $\Lambda_m$
and $\Lambda_m'$.  The eigenvectors are estimated
using a 2-neuron Parallel APEX Network trained on the incoming data.
We window the observation signal using a window of size
$N/m$ in order to produce vector data
${\bf x}_i=[x_i,\ldots,{\bf x}_{i-N/m+1}]$
which are consequently used for training the network.
Then the sum terms in Eq.~\ref{eq:APEX_detector} can be
denoted by the inner products
${\bf e}_1^T {\bf x}_{Ni/m}$ and
${\bf e}_2^T {\bf x}_{Ni/m}$. Notice that
these values correspond to the outputs of the two neurons of the APEX
network assuming that the weights of the network
store the corresponding eigenvectors. We'll estimate
these values with the actual outputs of the neurons
$y_{1,Ni/m}={\bf w}_1^T {\bf x}_{Ni/m}$ and
$y_{2,Ni/m}={\bf w}_2^T {\bf x}_{Ni/m}$ and our statistic
then becomes
\[
\tilde{\Lambda}_m=\sum_{i=0}^{m-1} (y_{1,Ni/m}^2+y_{2,Ni/m}^2)
\]

In Figure \ref{fig:detection_miss} we show the simulation comparison
between $\tilde{\Lambda}_m$ and $\Lambda_m$ using artificially created
data, with $\xi=10$ and corresponding $m=5$. The probability
of false alarm was set at $0.001$ and the appropriate threshold
was determined by 10000 simulation experiments for each detector.
Using this threshold we then computed, again via simulation,
the miss probability for the two detectors. We plot the result
in Fig.~\ref{fig:detection_miss} versus the signal to noise ratio
defined below
\[
SNR = 10 \log_{10} \left( N \frac{E\{s_k^2\}}{E\{n_k^2\}} \right)
\]

\begin{figure}
\centerline{\psfig{figure=exp1.ps,width=4in}}
\caption{Probability of miss comparison between (a)
the \protect{$m$}-th order noncoherent detector and (b) eigenvalue
based detector.}
\label{fig:detection_miss}
\end{figure}

\begin{figure}
\centerline{\psfig{figure=exper2_defl.ps,width=4in}}
\caption{Deflection ratio comparison between
(a) the \protect{$m$}-th order noncoherent detector and (b) the eigenvalue
based detector.}
\label{fig:detection_deflection_ratio}
\end{figure}
We see that the APEX detector clearly outperforms the
$m$-th order noncoherent detector by $0.5$ to $1db$
in the range of SNR between $-15db$ and $-10db$.

Another useful criterion for comparison is the deflection ratio
\[
H=\frac{(E_1\{\Lambda\}-E_0\{\Lambda\})^2}{var_0\{\Lambda\}}
\]
for a test statistic $\Lambda$. In our particular case
$H$ measures the separation of mean power of the projection
of the observation signal onto the subspaces specified by each
statistic under the two hypotheses
normalized by the variance of the projection power under the
null hypothesis.
Obviously the better the separation of the means the
better the performance of the detector. The comparison
of the two detectors under this criterion is depicted
in Fig.~\ref{fig:detection_deflection_ratio}. We see
that the eigenvector based method significantly outperforms
the $m$-order noncoherent detector which is consistent
with the results of Fig.~\ref{fig:detection_miss}.
This is to be expected since the PCA solution, by definition
maximizes the projection power of the signal
onto the low dimensional subspace
as discussed in Chapter \ref{ch:PCA}.

The price for the superior performance of the APEX detector
comes at the extra cost of training the neural network
to identify the principal subspace. Basically, this
amounts to one inner product ${\bf w}^T{\bf x}$, per
iteration which costs $O(N/m)$ operations per iteration.
The noncoherent detector requires $O(1)$ operations per
iteration so it is computationally simpler
than APEX but at the same time less attractive performance-wise.
So as with many other applications,
there is a trade-off between performance
and computational cost which can be balanced by the specific
needs of the user. Notice however that APEX has the additional advantage
over the noncoherent detector in that the number of components
used doesn't have to be 2. Using more components
in $\Lambda_m'$ will result in improved performance
but the exact amount of improvement is not known. We
propose to investigate this topic in future research.
Another advantage of the APEX-based scheme lies in 
its ability to track slowly changing parameters of the problem.

Finally, notice that the discrete nature of the APEX detector
makes it suitable for low frequency applications, e.g.
acoustics, but not suitable for laser optical communication
applications. If an analog implementation of the
algorithm can be viable then the application domain can be
extended even to such problems. However, the analog implementation
of APEX is not a topic that we shall explore further in this
thesis.




\section{System Identification}
\label{sec:applications-identification}

Consider the identification problem of a
Single-Input Single-Output linear system\footnote{
A similar formulation as the following one arises for a general
MIMO system. We present a SISO example here for simplicity.
}
\begin{eqnarray}
{\bf s}_{k+1} &=& {\bf F} {\bf s}_k + {\bf b} x_k \\
\hat{y}_k &=& {\bf c}^T {\bf s}_k
\end{eqnarray}

where ${\bf F} \in {\cal R}^{p \times p}$ is a non-singular matrix.
If the input sequence has finite length then we can write:
\begin{eqnarray}
\lefteqn{\hspace{-0.5in}\left[ \begin{array}{cccc}
y_{\ell} & y_{\ell+1} & \ldots & y_{\ell+N-1}
\\
y_{\ell+1} & y_{\ell+2} & \ldots & y_{\ell+N}
\\
\vdots & \vdots & & \vdots
\\
y_{\ell+m-1} & y_{\ell+m} & \ldots & y_{\ell+N+m-2}
\end{array} \right] =
\left[ \begin{array}{cccc}
\hat{y}_{\ell} & \hat{y}_{\ell+1} & \ldots & \hat{y}_{\ell+N-1}
\\
\hat{y}_{\ell+1} & \hat{y}_{\ell+2} & \ldots & \hat{y}_{\ell+N}
\\
\vdots & \vdots & & \vdots
\\
\hat{y}_{\ell+m-1} & \hat{y}_{\ell+m} & \ldots & \hat{y}_{\ell+N+m-2}
\end{array} \right] + noise =}
\nonumber \\
&=& \left[ \begin{array}{cccc}
h_{\ell+N-1} & \ldots & h_1 & h_0
\\
h_{\ell+N} & \ldots & h_2 & h_1
\\
\vdots & & \vdots & \vdots
\\
h_{\ell+N+m-2} & \ldots & h_m &  h_{m-1}
\end{array} \right]
\left[ \begin{array}{cccc}
x_0 & 0 & \ldots & 0
\\
\vdots & x_0 & \ldots & 0
\\
x_{\ell-1} & \vdots & \ddots & \vdots
\\
x_{\ell} & x_{\ell-1} & \ddots & x_0
\\
0 & x_{\ell} & & \vdots
\\
\vdots & \vdots & \ddots & x_{\ell-1}
\\
0 & 0 & \ldots & x_{\ell}
\end{array} \right] + noise
\label{eq:linear_system}
\end{eqnarray}

or in matrix form
\begin{equation}
{\bf Y} = \hat{\bf Y} + {\cal N} = {\bf H} {\bf X} + {\cal N}
\end{equation}

where ${\bf X} \in {\cal R}^{n \times N}$ is the input data matrix,
${\bf H} \in {\cal R}^{m \times n}$ is the system Hankel matrix
and ${\bf Y} \in {\cal R}^{m \times N}$ is the output
data matrix corrupted with noise ${\cal N}$. Notice that $n = \ell + N > N$
for $\ell > 0$, and since we assume $x_{\ell} \neq 0$,
${\bf X}$ has full rank $r = N$, because there is no linear combination
of the columns of ${\bf X}$ which is equal to 0 except the trivial one.
The problem then is to find an approximation
${\bf H}^{*}$ of ${\bf H}$, so that the square error
$\| {\bf Y} - {\bf H}^{*} {\bf X} \|_F^2$ is minimum.
From basic linear system theory
we know that ${\bf H}$ can be factored into the product
\begin{equation}
{\bf H}  = \left[ \begin{array}{l} {\bf c}^T \\ {\bf c}^T {\bf F} \\
\vdots \\ {\bf c}^T {\bf F}^{m-1} \end{array} \right]
\left[ {\bf b} ~|~ {\bf F} {\bf b} ~|~ \ldots ~|~ {\bf F}^n {\bf b}
\right] = {\cal O} ~ {\cal C}
\end{equation}

where ${\cal O} \in {\cal R}^{m \times p}$ and ${\cal C} \in {\cal R}^{p \times n}$
are the system observability and controllability
matrices respectively.
From this factorization it follows that
${\bf H}$ has rank $p$, therefore a
rank-$p$ approximation is sufficient to achieve a minimum error $E=0$
in the noise free case.
In some applications (e.g. speech processing)
we might have a priori information about the order of the
underlying system so a reduced rank approximation approach is
justified in this case.
Furthermore, since $n > N$, we cannot assume invertibility of the matrix
${\bf X}{\bf X}^T$ and therefore we cannot use the results
of previous published works on the subject.
Nevertheless, we can pad $\ell$ zero-columns
to both ${\bf X}$ and ${\bf Y}$ without changing the error function:
\[
E = \| {\bf Y} - {\bf W} {\bf X} \|_F^2 =
\left\| ~ [{\bf Y} ~ 0] - {\bf W} [{\bf X} ~ 0] ~ \right\|_F^2
\]

so the new matrix ${\bf X}' = [{\bf X} ~ 0]$ has the same number of rows
as columns and we can apply
the results of Section \ref{sec:linear_approximation_APCA} in
Chapter \ref{ch:APCA} on the problem.
Assuming that $\lambda_1 > \lambda_2 > \ldots > \lambda_s$
(condition true in general), then $\hat{r} = s$ is sufficient
to guarantee absence of local minima in the error function $E$
(cf. condition $(A)$ in page \pageref{page:condition_A})\footnote{
We have not shown that condition $(A)$ is necessary for the
absence of local minima. Nevertheless, simulation experiments
we have conducted seem to support this conjecture.}.
But $\hat{r} = \min \{ m,r \} = \min \{ m,N \}$
and $s = rank \{ {\bf Y} {\bf X}^T \}$. So
in the noise-free case, since $s = p$ and $m > p$ we can guarantee
freedom from local minima if $N = p$.
In the noisy case, $s = \min \{ m,N,n \} = \min \{ m,N \}$
(since $n \geq N$), so we always have $s = \hat{r}$;
thus freedom from local minima
is guaranteed for any value of $m$ and $N$. Of course it is implicitly
assumed that $m, n, N \geq p$ so $p \leq s$ as required
in the theoretical analysis of Section \ref{sec:linear_approximation_APCA}
in Chapter \ref{ch:APCA}.
If the error surface of the cost function contains no
local minima then any learning rule minimizing
\[
\left\| [\hat{\bf Y} ~ 0] - \ol{\bf W} \ul{\bf W}^T [{\bf X} ~ 0] \right\|_F^2
\]

will produce weights such that
\begin{equation}
\hat{\bf Y} = {\cal O} ~ {\cal C} ~ {\bf X} =
\ol{\bf W} ~ \ul{\bf W}^T ~ {\bf X}
\label{eq:OCX_eq_olWulWX}
\end{equation}

in the noise-free case. The matrices
${\cal C} {\bf X}$ and $\ul{\bf W}^T {\bf X}$ are both full rank
so there exists an invertible matrix ${\bf \Psi}$ such that
\begin{equation}
{\cal O} = \ol{\bf W} {\bf \Psi}
\label{eq:O_eq_olW}
\end{equation}

Thus, $\ol{\bf W}$ is the observability matrix of the system
after the similarity transformation $(\tilde{\bf F} , \tilde{\bf b} , \tilde{\bf c}^T )
= ( {\bf \Psi} {\bf F} {\bf \Psi}^{-1} , {\bf \Psi} {\bf b} ,
{\bf c}^T {\bf \Psi}^{-1} )$ and the first row
of $\ol{\bf W}$ is equal to $\tilde{\bf c}^T = {\bf c}^T {\bf \Psi}^{-1}$.
Let $\ol{\bf W}^{\uparrow}$ and $\ol{\bf W}^{\downarrow}$ denote the
matrices resulting from $\ol{\bf W}$ by removing the first and
the last row respectively. Then
\begin{equation}
\ol{\bf W}^{\uparrow} = \ol{\bf W}^{\downarrow}
{\bf \Psi} {\bf F} {\bf \Psi}^{-1} = \ol{\bf W}^{\downarrow} \tilde{\bf F}
\end{equation}

and $\tilde{\bf F}$ can be estimated by \cite{Kung83}
\begin{equation}
\tilde{\bf F} = \left\{ \ol{\bf W}^{\downarrow} \right\}^L
~ \ol{\bf W}^{\uparrow}
\label{eq:F_est}
\end{equation}

where the superscript $L$ denotes the left pseudoinverse.
Furthermore, from Eq.~\ref{eq:OCX_eq_olWulWX} and the fact that ${\cal O}$
has not less rows than columns, we obtain
\begin{equation}
{\cal C} {\bf X}  = {\bf \Psi}^{-1} \ul{\bf W}^T {\bf X}
\end{equation}

Equating the first columns of the two sides of the above equation we get
\begin{equation}
\sum_{i=0}^{\ell} x_i {\bf F}^{\ell-i} {\bf \Psi}^{-1}
{\bf \Psi} {\bf b} = {\bf \Psi}^{-1} \ul{\bf W}^T {\bf x}_1
\end{equation}

so $\tilde{\bf b}$ is computed as
\begin{equation}
\tilde{\bf b} = {\bf \Psi} {\bf b} = \left[ \sum_{i=0}^{\ell} x_i
\tilde{\bf F}^{\ell-i} \right]^{-1}
\ul{\bf W}^T {\bf x}_1
\label{eq:b_est}
\end{equation}

Equations \ref{eq:F_est} and \ref{eq:b_est} along with the fact
that the first row of $\ol{\bf W}$ is equal to $\tilde{\bf c}^T$
complete the solution procedure for the identification problem
in the noise free case. The system matrix ${\bf H}$,
can then be easily reconstructed by computing the Markov parameters
$h_i = \tilde{\bf c}^T {\tilde{\bf F}}^i \tilde{\bf b}$.
In the presence of noise however, this approach becomes an
approximation the accuracy of which deteriorates
with the increase of the noise power.
%%Some preprocessing stage is required
%%for removing the noise, before the actual system identification
%%takes place, in order to increase performance. However,
%%the design of such a preprocessing filter is beyond the scope
%%of this paper.

We simulated a 2-layer linear back-propagation network with $n=19$ inputs,
$m=20$ outputs and $p=10$ hidden units. The input to
the linear system is a finite random sequence of length
$\ell + 1 = 10$ and the output is corrupted with
white gaussian noise of adjustable power. The network is trained
by presenting pairs of columns from the 
matrices ${\bf X}$ and ${\bf Y}$, the
columns of ${\bf X}$ being inputs to the network and the corresponding
columns of ${\bf Y}$ being teachers. There are a total of
$N = 10 ~(= n - \ell)$ such input-teacher column pairs.
In order to produce ${\bf X}$ and ${\bf Y}$
a system of order 10 was simulated by creating random matrices
$({\bf F}, {\bf b}, {\bf c}^T),$ making sure that stability is preserved.
In Fig.~\ref{fig:H_est} we show the estimate of the impulse response
of the system produced by the method described above, after the
2000 sweeps over the training data.
The simulation was done on a SUN SPARC-station
and consumed approximately 7 minutes of execution time.
In Fig.~\ref{fig:spectral_est} we show the estimates of the eigenvalues
of the system, through the eigenvalues of $\tilde{\bf F}$.
For each plot, 100 separate experiments were run, and each
plot corresponds to a different noise power.
It is interesting to notice that identification methods 
performing standard linear regression
between ${\bf X}$ and ${\bf Y}$ of Eq.~\ref{eq:linear_system}
would be able to perfectly fit the data since the
problem is underdetermined. However, these methods
would perform worse than linear-approximation APCA
despite (or because of) the that fact.
The reason is, as already mentioned in Section
\ref{sec:linear_approximation_APCA} of Chapter \ref{ch:APCA},
that not knowing that the signal actually lives in
$p<N$ dimensions the linear regression models will simply
fit the noise in the remaining dimensions so the difference
between actual system output and estimated output would be
maximum.

\begin{figure}
\centerline{\psfig{figure=Hhat.ps,width=6in}}
\caption{The actual impulse response of the system (solid line)
vs. the neural network estimated response (dashed line)
after 2000 sweeps using linear back-propagation algorithm.}
\label{fig:H_est}
\end{figure}

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\psfig{figure=spectral_30.ps,width=2.5in,height=2.5in} &
\psfig{figure=spectral_31.ps,width=2.5in,height=2.5in}
\\
(a) & (b)
\\
\psfig{figure=spectral_32.ps,width=2.5in,height=2.5in} &
\psfig{figure=spectral_33.ps,width=2.5in,height=2.5in}
\\
(c) & (d)
\end{tabular}
\end{center}
\caption{Estimates of the complex poles of a simulated linear
system for various noise powers: (a) 0.05 (b) 0.1
(c) 0.2 (d) 0.4. The noise is white gaussian, and
each subfigure contains 100 experiments with different realizations
of the noise. The actual poles of the system are denoted by
``x'' and the estimated ones by ``\protect{$\cdot$}''.}
\label{fig:spectral_est}
\end{figure}



\section{Image Compression}
\label{sec:applications-compression}

Data compression is an essential task in many image processing
applications related to either picture database storage
or transmission. The necessity of compression emerges from the
fact that images are 2-dimensional entities -- even without considering
moving pictures where
time adds a third dimension -- so the amount of data
involved can often be quite large. Consequently, in moving pictures
the amount of data for either storage or transmission can be overwhelming
(compared to 1-dimensional applications, such as speech for example).
Given these facts it comes as no surprise that 
the success of most commercial products in imaging depend critically
on the amount of bit-rate reduction achievable
through data compression techniques
without significant deterioration of picture quality.
As a result, image compression
has become a focal issue of research in recent years.
Chapter \ref{ch:Motion} will study related image processing problems
in more detail. In this section we wish to demonstrate the
applicability and advantages
of using principal component techniques in bit-rate reduction.

{\em Transform Coding} is a major technique used for image compression
\cite{Netravali88}. It hinges upon the observation that for most pictures,
after a certain transformation of the image,
\begin{enumerate}
\item
the energy is concentrated
in few coefficients in the transform domain and
\item
the coefficients are approximately uncorrelated.
\end{enumerate}
The first part of the
observation is useful because it implies that we may store/transmit
only a few coefficients of the transformed image and still retrieve
a good approximation of the original using the inverse transform.
The second part of the observation tells us that when coding
a coefficient we don't have to worry about repeating information
included in other coefficients since there is no correlation between them.
It implies that coding the transform coefficients independently
is as good as coding them in blocks, a property that simplifies
the coding (and decoding) process.

To be more specific,
let us represent a (sampled) image of size $N_1\times N_2$
by the matrix ${\bf X} \in {\cal R}^{N_1\times N_2}$,
or alternatively, by a vector ${\bf x}$,
of size $N=N_1N_2$ produced by stacking
the columns of ${\bf X}$ on top of each other. For simplicity,
we are looking for linear transformations of the form

\begin{equation}
{\bf y} = {\bf W} {\bf x}
\label{eq:image_vector_transform}
\end{equation}
where ${\bf W} \in {\cal R}^{N\times N}$.
Our target is to approximate ${\bf x}$ as close as possible
using only few, say $m$ ($m<N$)
of the elements of ${\bf y}$. If the approximation
criterion is quadratic and if we consider ${\bf x}$ as an instance
of a random vector with autocorrelation ${\bf R} = E\{ {\bf x}{\bf x}^T \}$
then according to the discussion in Chapter \ref{ch:PCA}
the optimal solution is ${\bf W}_m = [{\bf e}_1~\ldots ~{\bf e}_m]^T$
comprised of the $m$ principal components of ${\bf x}$.
This transformation  satisfies both properties described above

\begin{enumerate}
\item we have
\begin{equation}
var_i = E\{y_i^2\} = E\{{\bf e}_1^T {\bf x}{\bf x}^T {\bf e}_1\}
= \lambda_i
\end{equation}
where $\lambda_1,\lambda_2,\ldots,\lambda_N$ are the eigenvalues
of ${\bf R}_x$ arranged in decreasing order, thus the first elements
of ${\bf y}$ have the largest energy (i.e. variance) while the energy
decreases steadily for smaller components. For a given number of
components the PC coefficients concentrate the maximum energy,
among all other linear transformations with the same
number of components.

\item the different transform coefficients (components) are uncorrelated:
\begin{equation}
E\{y_i y_j\} = E\{{\bf e}_i^T {\bf x}{\bf x}^T {\bf e}_j\} = 0,
~~~ i\neq j
\end{equation}
\end{enumerate}

Unfortunately for typical images a linear transformation of the whole
image can be very computationally intensive. In practice a picture
is divided into equal size blocks, while a 2-dimensional transform
(e.g. DCT, DFT, etc) is applied to each block. If ${\bf X}_n$ is the
$n$-th block of the image ${\bf X}$ then we define a  2-D transform
to be a transformation of the type
\begin{equation}
{\bf Y}_n = {\bf A} {\bf X}_n {\bf B}
\label{eq:2-D_transform}
\end{equation}
where ${\bf A}$ and ${\bf B}$ are two matrices (in the DCT and DFT
transforms these matrices are orthonormal and one is the transpose of
the other). It can be shown that if we stack the columns of ${\bf Y}_n$
and ${\bf X}_n$ on top of each other to produce the vectors
${\bf y}_n$ and ${\bf x}_n$ then Eq.~\ref{eq:2-D_transform} becomes
\[
{\bf y}_n = ({\bf A} \otimes {\bf B}^T) {\bf x}_n
\]
where $\otimes$ denotes the Kronecker product of two matrices.
It can also be shown that if ${\bf A}$ and ${\bf B}$ are
orthonormal then so is their Kronecker product.
Since not all matrices in general can be decomposed
into a Kronecker product
the transform \ref{eq:2-D_transform} is a special case of
\ref{eq:image_vector_transform}. It follows that the best
2-D transformation cannot be better than the best transformation
of the type \ref{eq:image_vector_transform}. However it has been
observed that the DCT transform comes very close to the performance
of the eigenvector transform and it is usually preferred in practice
because of its computational simplicity. The major disadvantage
of PC analysis is that the transform basis (i.e. the eigenvectors)
is image-dependent so not only does one have to compute the eigenvectors
but one also has to make them known to the decoder, by storing or
transmitting them, an act which implies additional bit rate.
However, if the decoder can estimate the eigenvectors
by using a similar picture as reference then the estimation
error can be significantly reduced compared to the DCT transform.

Consider a sequence of pictures where each one
is divided into equal size square blocks, then each block
is coded and transmitted to the decoder.
We will compare the following two methods

\begin{description}
\item[Method 1.] The blocks are coded using 2-D DCT.

\item[Method 2.] The image is divided into two halves, left and right,
while the blocks in each half are independently
transformed by an estimate of the eigenvectors
for the corresponding half. The estimates are obtained from
the eigenvector decomposition run on the right and left
halves of the previous picture in the sequence. Assuming
that the previous picture is not much different than the
current one we expect that its eigenvectors can serve
as a good basis for transformation of the current picture.
The advantage of this scheme is that the decoder can run
its own eigenvector estimator using the previous picture
which has already been decoded thus eliminating the need
for the eigenvector transmission.
In order to avoid computing the components from scratch
for every new picture in the sequence we
propose to use the APEX model for the eigenvector estimation.
\end{description}

Table~\ref{tab:dct_pca_mse_compare} shows the comparison
of methods 1 and 2 above with respect to the mean square
reconstruction error when using 8 components in both methods.
In method 1 the 8 components are the ones corresponding
to the highest variance of the transform coefficients and
are the ones shown in Table \ref{tab:8_high_dct_coefficients}.
For method 2 we simply transmit the 8 highest principal components
as given from the output of the APEX network. We see that we
can achieve approximately 15\% error reduction using the PCA method
over the DCT method. Fig.~\ref{fig:dct_pca_pictures} 
shows the output of the two methods for two sample pictures.
It is worth mentioning that the PCA method can be further improved
by partitioning the image into more than just left and right parts.
Indeed, there is nothing particular into dividing the picture
in two parts: dividing into 4 or 8 parts and computing customized
bases for each part based on the previous picture
will improve the performance.
However we don't expect that by partitioning further we will
keep improving the performance indefinitely because at a certain level the
eigenvectors of some part in the previous picture will become
poor estimators of the eigenvectors of the same part in the current picture. 
In \cite{Martinelli92} it was shown that a similar eigenvector
technique based on APEX, for representing images of similar statistics
results in improved performance versus DCT. The images considered
in that paper are still pictures which are clustered according
to eigenvector similarity. The authors found that their
clustering technique is ``...sufficiently rapid, since the
computation of the KLT in the successive steps requires only
a few adjustments of the step weights as the clustering processes.
Without the adaptive neural network used for computing KLT,
the procedure itself would require an excessive computational cost
since all the KLT should have been performed from the beginning.''

\begin{table}
\begin{center}
\begin{tabular}{||c|c|c||} \hline
Frame \# & DCT MSE & APEX MSE \\ \hline \hline
1 & 205.666 & 175.058 \\
2 & 202.678 & 173.766 \\
3 & 198.581 & 170.852 \\ \hline
\end{tabular}
\caption{Comparison of the MSE for method 1 and 2.}
\label{tab:dct_pca_mse_compare}
\end{center}
\end{table}

\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|} \hline
1 & 2 & 4 & 7 & & & & \\ \hline
3 & 5 & & & & & & \\ \hline
6 & & & & & & & \\ \hline
8 & & & & & & & \\ \hline
 & & & & & & & \\ \hline
 & & & & & & & \\ \hline
 & & & & & & & \\ \hline
 & & & & & & & \\ \hline
\end{tabular}
\caption{The 8 highest DCT coefficients.}
\label{tab:8_high_dct_coefficients}
\end{center}
\end{table}

\begin{figure}
\begin{center}
\begin{tabular}{c}
\psfig{figure=tennisd2.ps,width=4.5in}
\\
(a)
\\
\psfig{figure=tennisp2.ps,width=4.5in}
\\
(b)
\end{tabular}
\end{center}
\caption{(a) A sample reconstructed picture using method 1.
(b) Same using method 2.}
\label{fig:dct_pca_pictures}
\end{figure}

\section{Other Applications}

Constrained principal components have been used for image processing
problems as well. In \cite{Kung91b} it was demonstrated
that combining DCT and CPC can improve the mean square error
compared with using the same total number of DCT components.
The idea is to use the first, say 8, DCT components which 
correspond to 8 orthonormal sinusoidal vectors as a constraint subspace
for the CPC problem. Then using APEX one can extract the constraint
components that will further reduce the error more than
equal number of additional DCT components can.

In addition,
cross-correlation APCA networks
can be potentially used for filtering applications.
Remember that cross-correlation APCA
maximizes
$E\{\ul{\bf w}^T{\bf x}_k{\bf y}_k^T\ol{\bf w}\}$
under $\|\ul{\bf w}\|=\|\ol{\bf w}\|=1$,
which is equivalent to minimizing the cost
\begin{equation}
J = E\{(\ul{\bf w}^T{\bf x}_k-\ol{\bf w}^T{\bf y}_k)^2\}
\label{eq:APCA_mse_criterion}
\end{equation}
assuming that ${\bf x}_k$ and ${\bf y}_k$ are white random sequences,
in which case
\[
J=2-2 E\{(\ul{\bf w}^T{\bf x}_k{\bf y}_k^T\ol{\bf w}\}.
\]
If the vector sequences ${\bf x}_k$, ${\bf y}_k$ are
produced by time windowing two scalar signals $x_k$ and $y_k$,
then we can think of $\ul{\bf w}$, $\ol{\bf w}$ as the
tap coefficients of two linear filters
$\ul{W}(z)$, $\ol{W}(z)$. Then the criterion
\ref{eq:APCA_mse_criterion} corresponds to
minimizing the mismatch between the output of the two
filters.

Fig.~\ref{fig:Cross-Cor_APCA-echo} shows the block
diagram of a filter using the above analysis.
\begin{figure}
\centerline{\psfig{figure=crosscor_echo.ps,height=3in}}
\caption{Noise cancellation using cross-correlation.}
\label{fig:Cross-Cor_APCA-echo}
\end{figure}
Let $x_k=s_k+n_{1k}$ and $y_k=n_{2k}$,
where $s_k$ is some random signal we wish to enhance,
and $n_{1k}$ and $n_{2k}$ are two white noise
signals uncorrelated to $s_k$ but correlated to each other
whose correlation is difficult to compute analytically,
so we prefer to track it using a cross-correlation
APCA neural model discussed in Chapter \ref{ch:APCA}.
If the signal $x_k$ is approximately white
(as for example, in the case of very low SNR)
then the output error $e_k$, between the two filters is minimized.
But since $s_k$ and $n_{2k}$ are uncorrelated this implies
also that the difference between the output noise of the
two filters is minimized
\[
E\{(\ul{\bf w}^T{\bf n}_{1k} - \ol{\bf w}^T{\bf n}_{2k})^2\}
\]
In the ideal case we have
\[
\ul{W}(z)n_1(z)=\ol{W}(z)n_2(z)
\]
so
\[
e(z)=\ul{W}(z)s(z)
\]
i.e. the noise has totally disappeared but
$e(z)$ is a filtered version of $s(z)$.
In detection applications one may be interested only in
the power of $s$ which can be estimated by the
power of $e$ given knowledge of $\ul{W}(z)$.

Another possible application would result
from feeding $e(z)$ into an inverse filter
$1/\ul{W}(z)$ that would restore the original signal $s(z)$.
Notice that the filter structure in Figure~\ref{fig:Cross-Cor_APCA-echo}
followed by an inverse filter $1/\ul{W}(z)$,
is equivalent to IIR filter $\ol{W}(z)/\ul{W}(z)$, in other words
under the white noise assumption, {\em cross-correlation APCA
derives the optimal ARMA model for the echo cancellation problem}.
This is an extension to the Widrow LMS echo-cancellation scheme
\cite{Widrow75} which only considered the MA case alone.
However, the ARMA model has a serious disadvantage compared to
the simple MA model: stability is not guaranteed for the
APCA network since the model does not control the roots
of the polynomial $\ul{W}(z)$ which can often lie outside
the unit circle.
A possible way to come around this problem
is to incorporate some
additional constraint on the network to secure stability
but this issue extends beyond the scope of this discussion
and we propose to look at such problems in future research.
