\chapter{Oriented Principal Component Analysis}
\label{ch:OPCA}

\section{Introduction}
\label{sec:OPCA-introduction}

In this chapter we extend the standard principal component
analysis discussed in the previous chapter
by introducing {\em Oriented Principal Component
Analysis (OPCA)}.
OPCA corresponds to the generalized
eigenvalue problem of two random signals,
and bears the same relationship to Generalized
Eigenvalue Decomposition (GED), as PCA bears to
ordinary Eigenvalue Decomposition (ED).
More precisely, the goal is to find the direction vector ${\bf w}$,
that maximizes the signal to signal ratio
\begin{equation}
\mbox{maximize } J=\frac{E\Bigl\{({\bf w}^T {\bf x})^2 \Bigr\}}
{E\Bigl\{({\bf w}^T {\bf v})^2 \Bigr\}}=
\frac{{\bf w}^T{\bf R}_x{\bf w}}{{\bf w}^T{\bf R}_v{\bf w}}
\label{eq:OPC_criterion}
\end{equation}
where ${\bf R}_x=E\{{\bf x}_k{\bf x}_k^T\}$ and
${\bf R}_v=E\{{\bf v}_k{\bf v}_k^T\}$.
As usual, stationarity is assumed so none of ${\bf R}_x$,
${\bf R}_v$ or $J$ are functions of $k$.
No additional assumptions are made regarding ${\bf x}_k$ and ${\bf v}_k$
except that we will often find it useful to assume that
${\bf R}_v$ is strictly positive definite, hence non-singular.
The optimal solution to \ref{eq:OPC_criterion} will be called
the {\em principal oriented component} of the pair
of stochastic processes $(\{{\bf x}_k\},\{{\bf v}_k\})$.
Referring to Fig.~\ref{fig:OPC_geometry},
the name oriented principal
component is justified by the fact that the principal
component of ${\bf x}_k$ is now
steered by the distribution of the process ${\bf v}_k$:
it will be oriented towards the directions where ${\bf v}_k$ has
minimum energy while at the same time trying to maximize the
projection energy of ${\bf x}_k$. If ${\bf v}_k$ is evenly distributed
in all directions then it gives no ``guidance'' to the
direction vector. In that case ${\bf w}$ will simply try
to maximize the projected energy of ${\bf x}_k$,
which is equivalent to the ordinary PCA problem.
This happens when ~${\bf R}_v\propto {\bf I}$~ whence the process
${\bf v}_k$ is white noise. Thus OPCA is more general than
ordinary PCA, and is different from it only when ${\bf v}_k$
is a colored random process.

\begin{figure}
\centerline{\psfig{figure=opc_geom.ps,width=4in}}
\caption{A visual interpretation of Oriented Principal
Component Analysis: although the principal component
\protect{${\bf w}_{PC}$} is along the major axis of the signal,
the oriented principal component \protect{${\bf w}_{OPC}$}, is ``steered''
by the distribution of the noise.}
\label{fig:OPC_geometry}
\end{figure}

A limited but very interesting viewpoint for interpreting OPC
is by assuming that the vectors ${\bf x}_k$ and ${\bf v}_k$
are formed by successive time-windowing of two stationary
1-dimensional time series of which the first one we shall call signal
and the second one noise. For example,
using a window of size $n$ on the 1-dimensional
signals $x_k$ and $v_k$, we obtain ${\bf x}_k=[
x_k, x_{k-1}, \ldots, x_{k-n+1}]$,
${\bf v}_k=[v_k, v_{k-1}, \ldots, v_{k-n+1}]$.
In that case ${\bf w}=[w_1, w_2, \ldots, w_n]$ can be interpreted as
the coefficient vector of a linear filter, while the inner
products ${\bf w}^T{\bf x}_k$, ${\bf w}^T{\bf v}_k$,
denote the output of this filter at time $k$ when the input
is $x_k$ or $v_k$ respectively.
Consequently, $J$ is the signal to noise
ratio at the output of this filter which we strive to maximize.
This interpretation is equivalent to {\em Matched Filtering}
\cite{Thomas88} a classical tool used in detection and
communications applications.

It is known that the principal oriented
component is also the principal eigenvector of the symmetric generalized
eigenvalue problem \cite{Moor88}
\begin{equation}
{\bf R}_x {\bf w} = \lambda {\bf R}_v {\bf w}.
\end{equation}
From the properties of the generalized eigenvalue problem,
if ${\bf R}_x$ and ${\bf R}_v$ are symmetric then
all the generalized eigenvalues are real and nonnegative
and the eigenvectors are also real. The eigenvectors
are the stationary points of $J$ which is also
referred to as the {\em Rayleigh quotient} of the two matrices.
The second, third, etc, generalized eigenvectors
(ordered by decreasing associated eigenvalue)
will be also called oriented principal components.
So, if the generalized
eigenvalues are $\lambda_1 \geq \lambda_2 \geq \ldots \geq \lambda_n$,
then the associated eigenvectors ${\bf e}_1,~{\bf e}_2,
\ldots {\bf e}_n$ will be called the first, second, ... ,$n$-th
{\em oriented principal component} respectively.
In addition, all the components
maximize \ref{eq:OPC_criterion} but under different constraints, namely
the second component has to be ${\bf R}_v$-orthogonal
and ${\bf R}_x$-orthogonal to the first one,
the third component has to be the same to the first two, etc.
Formally, for the $i$-th component, $i \geq 2$, we have
\begin{eqnarray}
&& \frac{E\left\{({\bf e}_i^T {\bf x})^2 \right\} }
{E\Bigl\{({\bf e}_i^T {\bf v})^2 \Bigr\} } = \mbox{a maximum,} \nonumber
\\
&& \begin{array}{ll} \mbox{under} & \mbox{${\bf R}_v$-orthogonality: }
{\bf e}_i^T {\bf R}_v {\bf e}_j = 0 \\
 & \mbox{${\bf R}_x$-orthogonality: }
{\bf e}_i^T {\bf R}_x {\bf e}_j = 0 \end{array}
~~~\forall j<i.
\end{eqnarray}

In the following we'll show two extensions
of the Hebbian rule that can tackle the OPC problem
in an adaptive formulation.


\section{Model 1 for OPC Extraction}
\label{sec:OPC_model1}

We initially focus on the extraction of the first
component. Later we'll show how to extend to the multiple
component case.
The proposed model \cite{Diamantaras91} is a simplified linear neuron
with adaptable weights (Fig.~\ref{fig:OPC_net})
similar to the model proposed for ordinary PCA by Oja.

\begin{figure}
\centerline{\psfig{figure=unit.ps,height=6cm}}
\caption{The proposed linear neuron unit.}
\label{fig:OPC_net}
\end{figure}

The learning rule is also similar to the normalized Hebbian rule
except that it is duplicated for both signal
and noise and an additional adaptable parameter
$\xi$, is introduced which plays a crucial role in the correct
functioning of the algorithm. More precisely, the system has $n$
inputs $\{ \mbox{in}_1,~\mbox{in}_2, \ldots ,\mbox{in}_n \}$,
driven interchangeably by the two $n$-dimensional vector
stochastic processes $\{{\bf x}_k\}$ and $\{{\bf v}_k\}$,
and one output which we denote by
\begin{equation}
y_k={\bf w}^T{\bf x}_k
\end{equation}
or
\begin{equation}
z_k={\bf w}^T{\bf v}_k
\end{equation}
depending on whether the corresponding input vector is drawn from
process $\{{\bf x}_k\}$ or process $\{{\bf v}_k\}$.
The weights form the $n$-dimensional vector ${\bf w}$.
As already mentioned, we will assume that the input processes
are at least wide sense stationary and
we will further assume that the autocorrelation matrix
${\bf R}_v$ of ${\bf v}_k$ is strictly positive definite
in order to avoid the complications of infinite eigenvalues.
There is no {\em teacher} process involved in the training phase
thus the algorithm will be classified as unsupervised.
Nevertheless, there is an extra guidance for the system 
(implicit teacher) in the sense that it knows whether
each incoming pattern belongs to the process $\{{\bf x}_k\}$
or $\{{\bf v}_k\}$.

In the proposed algorithm, described below, the input data
are grouped in blocks of size $B$, where $B$ is some large
positive integer:

\begin{algorithm}
\label{alg:OPC_model1}
~
\begin{itemize}
\item[I.] Initialization:

Let $\theta$ be an upper bound of the largest eigenvalue of ${\bf R}_v$.
For example, $\theta$ can be equal to $\alpha T$,
where $T = E\{\|{\bf v}\|^2\}$ and $\alpha>1$ is a safety constant. Also let
\begin{eqnarray}
{\bf w}_0 &=& \mbox{any value} \neq 0
\\
\sigma_0^{(x)} ~=~ \sigma_0^{(v)} &=& 0
\\
\xi_0 &=& 0
\label{eq:OPC_model1_init}
\end{eqnarray}

\item[II.] Within block $\tau$:

For each iteration $k$:

\begin{eqnarray}
{\bf w}_{k+1}&=&{\bf w}_k+\beta (y_k {\bf x}_k - y_k^2 {\bf w}_k)
\label{eq:OPC_model1_dw_x}
\\
\sigma_{k+1}^{(x)} &=& \sigma_k^{(x)} + c \beta_k (y_k^2 - \sigma_k^{(x)})
\label{eq:OPC_model1_dsigmax}
\end{eqnarray}
if the input is drawn from $\{{\bf x}_k\}$, or else
\begin{eqnarray}
{\bf w}_{k+1}&=&{\bf w}_k-\xi_{\tau}\beta
( z_k {\bf v}_k - z_k^2 {\bf w}_k ) \nonumber
\\
&& ~~~ +\theta \beta (1-\|{\bf w}_k\|^2) {\bf w}_k
\label{eq:OPC_model1_dw_v}
\\
\sigma_{k+1}^{(v)} &=& \sigma_k^{(v)} + c \beta_k (z_k^2 - \sigma_k^{(v)})
\label{eq:OPC_model1_dsigmav}
\end{eqnarray}
where $c$ is any positive constant.

\item[III.] At the end of block $\tau$:
\begin{equation}
\xi_{\tau+1} = \xi_{\tau} + \gamma \left( \displaystyle
\frac{\sigma_{\tau B}^{(x)}}{\sigma_{\tau B}^{(v)}} - \xi_{\tau} \right),
~~~ 0 < \gamma \leq 1
\label{eq:OPC_model1_dxi}
\end{equation}
\end{itemize}
\end{algorithm}

Clearly, the second terms in the right-hand side
of Eqs.~\ref{eq:OPC_model1_dw_x}, \ref{eq:OPC_model1_dw_v},
are normalized Hebbian rules
corresponding to ${\bf x}_k$ and ${\bf v}_k$ with the term
in Eq.~\ref{eq:OPC_model1_dw_v}
participating with a negative weight $-\xi$. This is
intuitively correct since we want the network to keep
away from the principal component of ${\bf v}_k$ and
close to the principal component of ${\bf x}_k$ with
the value $\xi$ representing the compromise between the two extremes.
The third term in Eq.~\ref{eq:OPC_model1_dw_v} can be interpreted as a bias
term that is needed to keep the matrix ${\bf R}_x\!-\!\xi {\bf R}_v$
positive so that the weights will not converge to 0.
In the following we will show that the proposed algorithm converges to the
principal generalized eigenvector/eigenvalue pair
of the matrix pair $({\bf R}_x,{\bf R}_v)$, assuming that
\begin{itemize}
\item the samples from the two signals are interleaved, i.e.
the input sequence is ${\bf x}_0,{\bf v}_0,{\bf x}_1,{\bf v}_1,\ldots$.
\item the data block size $B$, is large enough.
\item ${\bf R}_v$ is strictly positive definite.
\item the multiplicity of the first generalized eigenvalue
is 1, i.e. $\lambda_1 > \lambda_2 \geq \ldots \geq \lambda_n > 0$.
\end{itemize}
Let us denote by ${\bf e}_i$ the $i$-th
generalized eigenvector of $({\bf R}_x,{\bf R}_v)$.
Then we can establish the following result

\begin{lemma}
If $\rho_1$ be the largest eigenvalue of the matrix
${\bf R}_x-\xi {\bf R}_v$ then
\begin{itemize}
\item[(a)] $\rho_1 > 0$ if $~\xi < \lambda_1$,
\item[(b)] $\rho_1 < 0$ if $~\xi > \lambda_1$,
\item[(c)] $\rho_1 = 0$ if $~\xi = \lambda_1$, and the corresponding
eigenvector is ${\bf e}_1$.
\hfill\qed
\end{itemize}
\label{lemma:OPC}
\end{lemma}

{\bf Proof.}
(a) It is sufficient to show that at least one eigenvalue of
${\bf R}_x-\xi{\bf R}_v$ is positive. Notice that
\[
{\bf e}_1^T ({\bf R}_x-\xi {\bf R}_v) {\bf e}_1 >
{\bf e}_1^T ({\bf R}_x-\lambda_1 {\bf R}_v) {\bf e}_1 = 0
\]
so indeed the matrix in question can be neither negative definite
nor negative semi-definite.

(b) For all generalized eigenvectors
\[
{\bf e}_i^T ({\bf R}_x-\xi {\bf R}_v) {\bf e}_i <
{\bf e}_i^T ({\bf R}_x-\lambda_1 {\bf R}_v) {\bf e}_i \leq 0
\]
Since the vectors ${\bf e}_i$ are linearly independent it follows
that the matrix in question is negative definite.

(c) We have
\[
({\bf R}_x-\xi {\bf R}_v) {\bf e}_1 =
({\bf R}_x-\lambda_1 {\bf R}_v) {\bf e}_1 = 0
\]
while for all $i$
\[
{\bf e}_i^T ({\bf R}_x-\xi {\bf R}_v) {\bf e}_i \leq
{\bf e}_i^T ({\bf R}_x-\lambda_1 {\bf R}_v) {\bf e}_i = 0
\]
so the matrix in question is semi-negative definite and
${\bf e}_1$ is an eigenvector with an
associated eigenvalue equal to 0.
\hfill \qed

The analysis of Algorithm~\ref{alg:OPC_model1} will be done in two phases:
within-block and block-to-block analysis.

{\bf I. Within block analysis.}
If $B\!\rightarrow\!\infty$ then the parameter $\xi$ does not change and
it can be considered constant at all times.
Then, as we saw in Chapter \ref{ch:PCA}
under some additional mild assumptions the evolution of the adapted parameters
${\bf w}_k$, $\sigma_k^{(x)}$ and $\sigma_k^{(v)}$
can be approximated by the evolution of the continuous-time
parameters ${\bf p}(t)$, $\eta_x(t)$ and $\eta_v(t)$ respectively,
which satisfy the following differential equations
\begin{eqnarray}
\frac{d{\bf p}(t)}{dt} &=&
({\bf R}_x - \xi {\bf R}_v + \xi \theta {\bf I}){\bf p}(t) \nonumber
\\
&& -{\bf p}(t)^T ({\bf R}_x - \xi {\bf R}_v + \xi \theta {\bf I}) {\bf p}(t) {\bf p}(t)
\label{eq:OPC_model1_dp(t)}
\\
\frac{d \eta_x(t)}{dt} &=& c({\bf p}(t)^T {\bf R}_x {\bf p}(t) - \eta_x(t))
\label{eq:OPC_model1_detax(t)}
\\
\frac{d \eta_v(t)}{dt} &=& c({\bf p}(t)^T {\bf R}_v {\bf p}(t) - \eta_v(t))
\label{eq:OPC_model1_detav(t)}
\end{eqnarray}
The right-hand sides of the above equations are the
conditional expectations of the right-hand sides of
Eqs. \ref{eq:OPC_model1_dw_x}, \ref{eq:OPC_model1_dw_v},
\ref{eq:OPC_model1_dsigmax}, \ref{eq:OPC_model1_dsigmav},
assuming ${\bf w}$ and $\xi$ constant.
The stability and convergency properties of these equations
can be now studied through the corresponding properties
of Eqs.~\ref{eq:OPC_model1_dp(t)} -- \ref{eq:OPC_model1_detav(t)}.

\begin{proposition}
Under the above assumptions and if ~$0 \leq \xi \leq \lambda_1$, ${\bf p}(t)$
converges to the normalized principal eigenvector (either one of the two)
of the matrix ${\bf R}_x - \xi {\bf R}_v$,
while $\eta_x(t)$ and $\eta_v(t)$ converge to some positive values $\mu_{\infty}^{x}$
and $\mu_{\infty}^{v}$ respectively,
such that $\xi \leq \displaystyle \frac{\mu_{\infty}^{x}}{\mu_{\infty}^{v}}
\leq \lambda_1$,
as $t \rightarrow \infty$.
\end{proposition}
{\bf Proof}

Define ${\bf A} = {\bf R}_x - \xi {\bf R}_v + \xi \theta {\bf I}$, which is
a symmetric, positive definite matrix
since $\theta$ is larger than the largest eigenvalue of ${\bf R}_v$.
We can then apply the theorem given in \cite{Oja82}
(a subcase of Theorem \ref{th:APEX} for $m=1$)
which showed that ${\bf p}(t)$ converges uniformly and asymptotically
to either ${\bf g}_1$ or $-{\bf g}_1$ where ${\bf g}_1$ is the normalized
principal eigenvector of ${\bf A}$. But ${\bf g}_1$ is also
the principal eigenvector of
${\bf R}_x - \xi {\bf R}_v = {\bf A} - \xi \theta {\bf I}$.
Therefore,
\begin{eqnarray}
&& {\bf p}(t)^T {\bf R}_x {\bf p}(t) \rightarrow {\bf g}_1^T {\bf R}_x {\bf g}_1 > 0,
\\
&& {\bf p}(t)^T {\bf R}_v {\bf p}(t) \rightarrow {\bf g}_1^T {\bf R}_v {\bf g}_1 > 0
\end{eqnarray}
as $t \rightarrow \infty$. From this it follows that
Eqs.~\ref{eq:OPC_model1_detax(t)} and
\ref{eq:OPC_model1_detav(t)}
are stable and asymptotically converging to the values
$\mu_{\infty}^x = {\bf g}_1^T {\bf R}_x {\bf g}_1$ and
$\mu_{\infty}^v = {\bf g}_1^T {\bf R}_v {\bf g}_1$ respectively.
In addition, from Lemma~\ref{lemma:OPC}
it follows that if $0 \leq \xi \leq \lambda_1$ then
${\bf g}_1^T {\bf R}_x {\bf g}_1 - \xi {\bf g}_1^T {\bf R}_v {\bf g}_1 \geq 0$, thus
$\xi \leq \displaystyle \frac{\mu_{\infty}^x}{\mu_{\infty}^v}$.
Furthermore, since $\lambda_1$ maximizes \ref{eq:OPC_criterion}, we have finally
$\xi \leq \displaystyle \frac{\mu_{\infty}^x}{\mu_{\infty}^v} \leq \lambda_1$.
\hfill \qed

The main result of \cite{Ljung77} allows us to imply also that
${\bf w}_k \rightarrow {\bf g}_1$, $\sigma_k^{(x)} \rightarrow \mu_{\infty}^x$
and $\sigma_k^{(v)} \rightarrow \mu_{\infty}^v$, as $k \rightarrow \infty$.
The underlying assumption regarding $\beta$ is the same as with
Theorem \ref{th:APEX} but again the assumption can be relaxed
setting $\beta$ equal to some small constant trading off
strict convergence for tracking capability.

{\bf II. Block-to-block analysis.}
In practice $B$ is large but finite.
Assuming that the above proposition is still approximately true
under this new condition, we can say that
${\bf w}_{\tau B} \approx {\bf g}_{1\tau}$,
where ${\bf g}_{1\tau}$ is the normalized principal eigenvector
of the matrix ${\bf R}_x - \xi_{\tau} {\bf R}_v$.
In addition Eq.~\ref{eq:OPC_model1_dxi} yields
\begin{eqnarray}
\xi_{\tau+1} &\approx& \xi_{\tau} + \gamma \left(
\frac{{\bf g}_{1\tau}^T {\bf R}_x {\bf g}_{1\tau}}
{{\bf g}_{1\tau}^T {\bf R}_v {\bf g}_{1\tau}}
- \xi_{\tau} \right)
\label{eq:OPC_dxitau}
\\
& \geq & \xi_{\tau} \nonumber
\end{eqnarray}
The sequence $\xi_{\tau}$ is nondecreasing and bounded from
above by $\lambda_1$, thus it will converge to some value $\rho \leq \lambda_1$.
But the only limit point of \ref{eq:OPC_dxitau} is $\xi = \lambda_1$.
So, $\xi_{\tau} \rightarrow \lambda_1$
as $\tau \rightarrow \infty$ and ${\bf w}_{\tau B}^{~}$ converges to the
normalized principal eigenvector of ${\bf R}_x-\lambda_1{\bf R}_v$ which is
equal to ${\bf e}_1$ as shown in Lemma \ref{lemma:OPC}.


{\bf Simulation Results}

We tested algorithm \ref{alg:OPC_model1} using various artificially
created random sequences for various input dimensions and block
sizes $B$. Figure \ref{fig:OPC_model1} shows a typical run for an
input dimension $n=10$ and block size $B=200$. In order to demonstrate the
convergence ability of the algorithm, we repeat the same set
of $200$ data points in sweeps. The input data were constructed
by time-windowing two signals using a window of size $n=10$.
The signals were colored by passing two white noise sequences through
different AR filters. The value of $\gamma$ was set experimentally
at $0.95$ while the value of $\beta$ was determined according to our
discussion in Section \ref{sec:APEX} in Chapter \ref{ch:PCA}
since within each block the algorithm is equivalent to the Oja rule
which in turn is a special case of APEX. So our discussion regarding
the value of the step-size constant for APEX applies here too.
The parameter $\theta$ was estimated by averaging $\|{\bf v}_k\|^2$
for 100 iterations before starting the algorithm (notice that
in the plots of Figure \ref{fig:OPC_model1} the adaptation
really starts at $k=200$ or $\tau=1$). In Eqs.~\ref{eq:OPC_model1_dsigmax}
and \ref{eq:OPC_model1_dsigmav} we used a constant learning rate
$\frac{1}{5B}$ instead of $c\beta_k$. Any other small constant
value in the range $1/B$ to $1/(10B)$ would do a reasonable job,
however outside this range the convergence of $\xi$ could be very slow.
We observe that the convergence of the algorithm is relatively slow
compared to APEX but it is smooth and still reasonably fast.
In addition, as expected from the fact that algorithm
\ref{alg:OPC_model1} is really a sequence of short APEX's,
the final convergence vector has unit norm.

\begin{figure}
\begin{center}
\begin{tabular}{c}
\psfig{figure=opc_w.ps,width=4in}
\\
(a)
\\
\psfig{figure=opc_xi.ps,width=4in}
\\
(b)
\end{tabular}
\end{center}
\caption{OPC Model 1: (a) Plot of the estimation error
\protect{$\|{\bf e}_1-{\bf w}_k\|^2$}.
(b) The evolution of \protect{$\xi_{\tau}$}.}
\label{fig:OPC_model1}
\end{figure}


 
\section{Model 2 for OPC Extraction}
\label{sec:OPC-model2}

Motivated by reducing the complexity of the previous model
for OPC extraction, we set off to investigate other alternatives
to the problem. First we'll discuss the use of a Lyapunov energy function
for OPC.
The principal oriented component maximizes $J$ in Eq.~\ref{eq:OPC_criterion},
the maximum value being the generalized eigenvalue $\lambda_1$.
Thus the function
\[
V({\bf w}) = \frac{1}{2} \left( \lambda_1 - J({\bf w}) \right)
\]
is positive definite and it can serve as a Lyapunov energy function
for our system. The proper gradient descent algorithm would
correspond to the differential equation

\begin{equation}
\frac{d{\bf w}}{dt} = - \nabla V = \frac{1}{2} \nabla J =
\frac{1}{{\bf w}^T{\bf R}_v{\bf w}} \left(
{\bf R}_x{\bf w} - \frac{{\bf w}^T{\bf R}_x{\bf w}}
{{\bf w}^T{\bf R}_v{\bf w}} {\bf R}_v{\bf w} \right)
\end{equation}
In fact, even if

\begin{equation}
\frac{d{\bf w}}{dt} =
\left( {\bf R}_x{\bf w} - \frac{{\bf w}^T{\bf R}_x{\bf w}}
{{\bf w}^T{\bf R}_v{\bf w}} {\bf R}_v{\bf w} \right)
\end{equation}
we can show that $dV/dt \leq 0$~ for all $t$, since
\begin{equation}
\frac{dV}{dt}=\frac{d{\bf w}^T}{dt}\nabla V
=-\frac{1}{{\bf w}^T{\bf R}_v{\bf w}}
\left\|{\bf R}_x{\bf w} - \frac{{\bf w}^T{\bf R}_x{\bf w}}
{{\bf w}^T{\bf R}_v{\bf w}} {\bf R}_v{\bf w} \right\|^2 \leq 0
\end{equation}
therefore the system would be asymptotically stable.
The critical problem is the estimation of the ratio
${\bf w}^T{\bf R}_x{\bf w}/{\bf w}^T{\bf R}_v{\bf w}$
by our stochastic algorithm since the matrices
${\bf R}_x$ and ${\bf R}_v$ are not readily available.
Our model will have to track this ratio as the data are processed.
This can be achieved by the following simple adaptive rule
\begin{equation}
{\bf w}_{k+1} = {\bf w}_k+\beta (y_k {\bf x}_k -
\|{\bf w}_k\|^2 z_k {\bf v}_k)
\label{eq:OPC_model2_init}
\end{equation}
which is basically formed
by two counterbalanced Hebbian rules (the terms
$y_k {\bf x}_k$ and $z_k {\bf v}_k$). The differential equation
associated with Eq.~\ref{eq:OPC_model2_init} is

\begin{equation}
\frac{d{\bf w}}{dt}={\bf R}_x{\bf w}-\|{\bf w}\|^2{\bf R}_v{\bf w}
\label{eq:OPC_model2_init_ode}
\end{equation}
The dynamics of $\|{\bf w}\|^2$ are determined by premultiplying
Eq.~\ref{eq:OPC_model2_init_ode} by $2{\bf w}^T$ to obtain
\begin{equation}
\frac{d\|{\bf w}\|^2}{dt}=2({\bf w}^T{\bf R}_x{\bf w} -
\|{\bf w}\|^2{\bf w}^T{\bf R}_v{\bf w})
\end{equation}

So the norm of ${\bf w}$ indeed tracks the ratio
$J={\bf w}^T{\bf R}_x{\bf w}/{\bf w}^T{\bf R}_v{\bf w}$,
i.e. $\|{\bf w}\|^2$ increases if less than $J$ while
it decreases if greater.
Although this is not a rigorous proof that Eq.~\ref{eq:OPC_model2_init}
extracts the principal oriented component
it is as close we can come to a proof.
Simulation results show that this rule indeed extracts the 
principal oriented component for the two random sequences.
However we found that the performance can be significantly
improved in terms of the variance of the estimation
if the norm $\|{\bf w}_k\|^2$ in Eq.~\ref{eq:OPC_model2_init} is replaced
by another parameter $\xi_k$ which tracks it and
has smoother evolution in time.
So in the following we will propose a slight modification
of rule \ref{eq:OPC_model2_init} defined below
\begin{eqnarray}
{\bf w}_{k+1} &=& {\bf w}_k+\beta (y_k {\bf x}_k - \xi_k z_k {\bf v}_k)
\label{eq:OPC_model2_w}
\\
\xi_{k+1} &=& \xi_k+\beta(\|{\bf w}_k\|^2-\xi_k)
\label{eq:OPC_model2_xi}
\end{eqnarray}
Again we find it very difficult to rigorously prove that
model \ref{eq:OPC_model2_w}-\ref{eq:OPC_model2_xi}
indeed extracts the required components
and we'll content ourselves with the demonstration of the
convergence via simulation experiments.

{\bf Simulations.}

We compare the performance of model 2 with model 1 running
the algorithm on the same data used in the simulation experiments
of the previous section.
Figure \ref{fig:OPC_model2} shows the performance
of model 2. The value of the step-size constant was set,
after experimentation, at $\beta=0.01$.
Part (a) of the plot shows the estimation error between ${\bf w}_k$
and the actual component ${\bf e}_1$ while part (b)
shows the evolution of the estimator $\xi_k$,
of the eigenvalue $\lambda_1$.
Although $\xi$ converges slower than for model 1 and
it follows a characteristic damped oscillation,
we observe that the convergence of ${\bf w}_k$ is
both exponential and significantly faster
than the convergence of model 1 (notice
that in Figure \ref{fig:OPC_model1}(a) the $x$-axis extends
from 0 to 20000 while here it only extends up to 10000).
Furthermore, model 2 is simpler,
contains fewer parameters and appears more aesthetically appealing.
The disadvantage of this model is that the final
convergence vector is not normalized and in fact its norm is equal
to $\lambda_1$ which can be significantly large
if ${\bf R}_v$ is near-singular. Luckily, when ${\bf v}_k$
is noise the condition number of ${\bf R}_v$, i.e. the ratio
between its largest and smallest eigenvalue, is
not very high and assuming ${\bf x}$ and ${\bf v}$
are reasonably scaled, $\lambda_1$ is a reasonable
size scalar.

\begin{figure}
\begin{center}
\begin{tabular}{c}
\psfig{figure=opc_simple_w.ps,width=4in}
\\
(a)
\\
\psfig{figure=opc_simple_xi.ps,width=4in}
\\
(b)
\end{tabular}
\end{center}
\caption{(a) Plot of the estimation error
\protect{$\|{\bf e}_1-{\bf w}_k\|^2$}.
(b) Plot of \protect{$\xi_k$}.}
\label{fig:OPC_model2}
\end{figure}


\section{Extracting Multiple OPCs}
Once the first generalized eigenvector ${\bf e}_1$, is extracted
a deflation procedure is followed similar to the classical PC case
in order to extract the rest of the components.

\begin{definition}
\label{def:OPC_deflation}
Any one of the following mappings
\begin{eqnarray}
({\bf R}_x,{\bf R}_v) &\rightarrow& \left( \left\{ {\bf I}-
\frac{{\bf R}_x{\bf e}_1{\bf e}_1^T}{{\bf e}_1^T{\bf R}_x{\bf e}_1}
\right\}
{\bf R}_x \left\{ {\bf I}-
\frac{{\bf e}_1{\bf e}_1^T{\bf R}_x}{{\bf e}_1^T{\bf R}_x{\bf e}_1}
\right\},{\bf R}_v \right)
\\
({\bf R}_x,{\bf R}_v) &\rightarrow& \left( \left\{ {\bf I}-
\frac{{\bf R}_v{\bf e}_1{\bf e}_1^T}{{\bf e}_1^T{\bf R}_v{\bf e}_1}
\right\}
{\bf R}_x \left\{ {\bf I}-
\frac{{\bf e}_1{\bf e}_1^T{\bf R}_v}{{\bf e}_1^T{\bf R}_v{\bf e}_1}
\right\},{\bf R}_v \right)
\\
({\bf R}_x,{\bf R}_v) &\rightarrow& \left( {\bf R}_x \left\{ {\bf I}-
\frac{{\bf e}_1{\bf e}_1^T{\bf R}_x}{{\bf e}_1^T{\bf R}_x{\bf e}_1}
\right\}, {\bf R}_v \right)
\\
({\bf R}_x,{\bf R}_v) &\rightarrow& \left( {\bf R}_x \left\{ {\bf I}-
\frac{{\bf e}_1{\bf e}_1^T{\bf R}_v}{{\bf e}_1^T{\bf R}_v{\bf e}_1}
\right\}, {\bf R}_v \right)
\end{eqnarray}
is called a {\bf deflation transformation} of the matrix pair
$({\bf R}_x,{\bf R}_v)$ with respect to the generalized eigenvalue
problem.
\end{definition}

As with the standard PCA case, the matrix resulting after deflation
has the same eigenvalue/vector pairs except that the eigenvalue
associated which ${\bf e}_1$ becomes 0. So, after deflation,
the originally second component becomes dominant. One can easily
verify the above claims by showing that if $(\tilde{\bf R}_x,{\bf R}_v)$
is anyone of the four pairs resulting after the deflation then
\[
\tilde{\bf R}_x {\bf e}_1 = 0
\]
while
\[
\tilde{\bf R}_x {\bf e}_i = {\bf R}_x {\bf e}_i, ~~~i>1
\]
based on the fact that the generalized eigenvectors are
${\bf R}_x$- and ${\bf R}_v$-orthogonal. Therefore, the relationship
\[
\tilde{\bf R}_x {\bf e}_i = \tilde{\lambda}_i {\bf R}_v {\bf e}_i
\]
holds for all $i$ where now $\tilde{\lambda}_1=0$
and $\tilde{\lambda}_i=\lambda_i$ for all $i>1$.

We can achieve deflation directly on the data for example, through
the transformation
\begin{eqnarray}
\tilde{\bf v}_k &=& {\bf v}_k
\\
\tilde{\bf x}_k &=& {\bf x}_k - y_k \frac{E\{y_k{\bf x}_k\}}{E\{y_k^2\}}\nonumber
\\
&=& {\bf x}_k - \frac{{\bf R}_x{\bf w}}{{\bf w}^T{\bf R}_x{\bf w}}{\bf x}_k
\label{eq:OPC_deflation_x}
\end{eqnarray}
assuming that the weight vector ${\bf w}$ is proportional to the
principal eigenvector ${\bf e}_1$. Then the autocorrelation matrices
$\tilde{\bf R}_x=E\{\tilde{\bf x}_k\tilde{\bf x}_k^T\}$,
$\tilde{\bf R}_v=E\{\tilde{\bf v}_k\tilde{\bf v}_k^T\}={\bf R}_v$,
form a deflated pair of the first type in Definition~\ref{def:OPC_deflation}.
The expectations in Eq.~\ref{eq:OPC_deflation_x} can be approximated
by averages computed off-line after the first neuron has converged.
This limits the usefulness of the model since it
forbids the parallel estimation of more than one components.
An alternative would be the use of a lateral connection network
as we did in APEX to enforce ${\bf R}_v$-orthogonality
between the units. We propose to investigate
this interesting alternative in a future work. 

 
\section{OPC and classification, control and tree-perceptrons.}
\label{sec:OPC_relations}

We have already observed the relationship between OPC
and matched filtering which implies that the neural models
described above can also be thought of as adaptive matched-filters.
Below we shall discuss briefly the relations between
OPC and various other techniques used for classification,
control or neural applications.

\subsection{Pattern Classification (Fisher's Linear Discriminant)}

In the pattern classification context the vectors
$\{{\bf x}_k\}$ and $\{{\bf v}_k\}$ are
random samples drawn from two different classes of data
which have different distributions. The
task is to classify observed patterns into their
corresponding classes as correctly as possible given
some knowledge about the distributions or given a set of
training samples from which the classifier can
learn these distributions.
A classical second order classification technique (i.e. based
only on first and second order statistics)
is Fisher's linear discriminant analysis. It is based
on the principal generalized eigenvector of a pair of symmetric
matrices ${\bf S}_w$ and ${\bf S}_b$ called the {\em within-cluster} and
{\em between-cluster} scatter matrices. In particular,
let ${\bf m}_1$, ${\bf m}_2$, be the means of the  two classes
and ${\bf u}_{1k}$, ${\bf u}_{2k}$ be the samples drawn from each class.
Then the Fisher criterion is

\begin{eqnarray}
\mbox{maximize } J &=& \frac{({\bf w}^T{\bf m}_1-{\bf w}^T{\bf m}_2)^2}
{E\left\{\Bigl[{\bf w}^T({\bf u}_1-{\bf m}_1)\Bigr]^2\right\}+
E\left\{\Bigl[{\bf w}^T({\bf u}_2-{\bf m}_2)\Bigr]^2\right\}}
\\
&=& \frac{{\bf w}^T{\bf S}_b{\bf w}}{{\bf w}^T{\bf S}_w{\bf w}}
\end{eqnarray}
where ${\bf S}_b=({\bf m}_1-{\bf m}_2)({\bf m}_1-{\bf m}_2)^T$
and ${\bf S}_w=E\{({\bf u}_1-{\bf m}_1)({\bf u}_1-{\bf m}_1)^T\} +
E\{({\bf u}_2-{\bf m}_2)({\bf u}_2-{\bf m}_2)^T\}$.
Fisher analysis opts to
separate the two classes as much as possible using
1-dimensional projection by maximizing the
distance between the projection means
normalized by the spread of the projection
(for more details see \cite[chapter5]{Duda73}).
The problem is easily transformed into an OPC formulation by setting
${\bf x}_k={\bf m}_1-{\bf m}_2$, and
${\bf v}_k={\bf u}_{2k}-{\bf m}_2$ or
${\bf v}_k={\bf u}_{1k}-{\bf m}_1$
in which case ${\bf R}_x={\bf S}_b$ and ${\bf R}_v={\bf S}_w$.

Another popular group of methods for pattern classification
based on PCA is {\em subspace classifiers}
\cite{Fukunaga70,Kittler78,Oja83}.
In order to battle the curse of dimensionality, subspace classifiers
describe each class by a small dimensional principal component subspace.
The main criticism against these methods
is that the representation of each class is not taking into account
the relative difference with the other classes. For example, in the cases where
the first few principal components of two classes are very similar to each other
subspace classifiers have very poor classification performance.
Many attempts have been carried out for improving the PC-based
methods by introducing coupling terms from the other classes \cite{Oja87}
but these methods are rather heuristic.
OPC analysis offers a PC-like alternative for taking the other-class
information into consideration and also has the advantage
that it is based on a specific optimality criterion.


\subsection{Control and Signal Processing}

Recently GED and GSVD have been found useful in numerous
applications like control system identification \cite{Moonen90},
system modeling via the total least-squares criterion \cite{Willems87},
high resolution sensor-array processing \cite{Roy86a,Roy86b},
mother-fetal ECG separation \cite{Vanderschoot87}, etc.
For a comprehensive summary see \cite{Moor88}. Due to the relation
between OPC and GED/GSVD the OPC methods proposed in this chapter
are applicable to these problems as well. However in order for our
models to be practical we need to identify
formulations of the above problems that need recursive
solutions, so that solving the problem from scratch every
time new data arrive becomes undesirable. Such formulations
remain an open question for future research.


\subsection{Training Decision Tree Networks}

In \cite{Koutsogeras88} the authors proposed a tree-structured neural model
for pattern classification based on the Fisher's linear
discriminant analysis that clusters the input space
using hyperplanes. Thus, the proposed OPC models
can thus be employed for training the
individual units of the network which correspond to
these hyperplanes.
Similar tree-networks have been proposed
also in \cite{Utgoff88,Utgoff90}.

