\chapter{Principal Component Neural Networks}
\label{ch:PCA}


\section{Introduction}

The problems of data compression and pattern classification have
attracted a lot of attention from many researchers of various fields like
pattern recognition, artificial intelligence, signal processing,
etc.  Both problems rely on finding an efficient representation
of the input data. This representation should extract the essential
information out of the original sequence of patterns.
For data compression this process is
a mapping from a higher dimensional (input) space to a lower
dimensional (representation) space.
The data compression problem is also important in modeling
multilayer perceptrons where the hidden neurons may be regarded as
a layer of nodes corresponding to the most effective representation
of the input patterns. Their activation patterns can be viewed as the
target patterns of the representation transformation that may facilitate the
``classification'' problem tackled by the following layer(s).

\section{Principal Component Analysis.}
\label{sec:PCA-formulation}

A powerful statistical tool for extracting such representations is
Principal Component Analysis (PCA) \cite{Hotelling33} which derives the
optimal (in the sense of the mean square projection error)
subspace ${\cal L}$, on which to project the input to.
In other words,
for a given input vector stochastic process $\{ {\bf x}_k \} \in {\cal R}^n$,
assumed at least wide sense stationary,
and a given subspace dimension $m<n$, we are looking for the
optimal linear transformation:
\[
{\bf y}_k = {\bf W}{\bf x}_k
\]
where ${\bf W}$ is an {\em orthonormal basis} of ${\cal L}$,
so that a certain error is minimal.
The error criterion is that the reconstructed input:
\begin{equation}
\hat{\bf x}_k = {\bf W}^T {\bf y}_k = {\bf W}^T {\bf W} {\bf x}_k
\end{equation}
which is the projection of ${\bf x}_k$ onto ${\cal L}$,
should minimize
\begin{eqnarray}
J &=& E \{ \| {\bf x}_k - \hat{\bf x}_k \|^2 \} =
trace\Bigl( E\{({\bf x}_k-\hat{\bf x}_k)
({\bf x}_k-\hat{\bf x}_k)^T\} \Bigr) \nonumber
\\
&=& E\{\|{\bf x}_k\|^2\}+ trace\Bigl( E\{{\bf W}^T{\bf W}{\bf x}_k{\bf x}_k^T
{\bf W}^T{\bf W}\} \Bigr) -2 trace\Bigl(
E\{{\bf x}_k{\bf x}_k^T{\bf W}^T{\bf W}\} \Bigr) \nonumber
\\
&=& var({\bf x}) + trace\Bigl( {\bf W}^T{\bf W}{\bf R}_x
{\bf W}^T{\bf W} \Bigr) - 2 trace\Bigl( {\bf W}{\bf R}_x{\bf W}^T \Bigr)
\label{eq:PCA_cost}
\end{eqnarray}
where
\[
{\bf R}_x = E\{ {\bf x}_k {\bf x}_k^T \}
\]
is called the {\em auto-correlation matrix} for the input process ${\bf x}_k$,
and we made use of the fact that $trace({\bf A}{\bf B}^T)=trace({\bf B}^T{\bf A})$
for all matrices ${\bf A}$, ${\bf B}$ with same dimensions.
Notice that since we assume stationarity neither ${\bf R}_x$
nor $J$ are functions
of time. Let the eigenvalues of ${\bf R}_x$ be
$\lambda_1, \lambda_2, \ldots, \lambda_n$, arranged in decreasing order,
with corresponding normalized eigenvectors ${\bf e}_1, {\bf e}_2,
\ldots, {\bf e}_n$. Then the minimizer of \ref{eq:PCA_cost} is
\[
{\bf W}^{*} = [{\bf e}_1\,\ldots\,{\bf e}_m]^T
\]
while the corresponding minimal error is
\begin{equation}
J_{min} = \sum_{i=m+1}^{n} \lambda_i.
\label{eq:PCA_Jmin}
\end{equation}

Since ${\bf R}_x$ is an autocorrelation matrix,
it follows that it is symmetric and semi-positive
definite, i.e. all the eigenvalues
are real and nonnegative. Thus the right-hand side
of Eq.~\ref{eq:PCA_Jmin} is a nonnegative quantity
in consistency with the left-hand side which must
also be nonnegative since it is a quadratic form.
The (normalized) eigenvectors of ${\bf R}_x$
corresponding to its highest eigenvalues (i.e. the rows of ${\bf W}^{*}$)
are called {\em principal eigenvectors}.
Since ${\bf R}_x$ is symmetric
the rows of ${\bf W}$ are orthonormal forming a basis
of the {\em principal component subspace}.
The mapping ${\bf y}_k={\bf W}^{*}{\bf x}_k$, is also
called the {\em Karhunen-Lo\`{e}ve Transform} of ${\bf x}_k$
\cite{Watanabe65} although originally both Karhunen \cite{Karhunen46}
and Lo\`{e}ve \cite{Loeve48} investigated statistical
orthogonal expansions of continuous time stochastic processes.
The reconstruction $\hat{\bf x}_k = {\bf W}^T{\bf y}_k$
for $p<n$, is consequently called a truncated K-L expansion of ${\bf x}_k$.

An alternative criterion maximized by the principal eigenvectors
is the constrained output variance
\begin{equation}
\begin{array}{c}
J' = var\{{\bf w}^T{\bf x}\} = E\{({\bf w}^T{\bf x})^2\} =
{\bf w}^T {\bf R}_x {\bf w},
\\
\mbox{under } \|{\bf w}\|=1.
\end{array}
\end{equation}
For ${\bf w} = {\bf e}_1$, the constrained output variance
is maximum, while for ${\bf w} = {\bf e}_p$, $p>1$
the variance is maximized under the additional constraint that
${\bf w}$ is orthogonal to all higher components than $p$,
namely to the vectors ${\bf e}_1,\ldots,{\bf e}_{p-1}$.
Figure \ref{fig:pca_axes} shows a geometric interpretation
of the principal component subspaces: based on the variance
criterion the principal component direction should be the one
where the signal has more energy (i.e. projection variance);
the least principal direction
is the one with the least energy. If the signal is zero-mean then
the maximum energy direction is also the direction of maximum spread
or in information-theoretical parlance, the direction that
contains the most information regarding the signal. It should be
also intuitively clear from the same figure that the principal
direction is the one corresponding to the minimum reconstruction
error. For more thorough discussion of the properties of the
eigenvectors and eigenvalues of symmetric matrices the reader
is referred to \cite{Parlett80}.

\begin{figure}
\centerline{\psfig{figure=linsker.ps,width=3.5in}}
\caption{The principal component: (1) contains most of
the energy of the signal, (2) has the minimum signal projection
error and (3) keeps the maximum information of the signal,
if the latter is zero mean. The least principal axis has just the
opposite properties.}
\label{fig:pca_axes}
\end{figure}

In many situations the autocorrelation matrix is not known in advance
but it may be estimated by
a time-average over a sample sequence of input vectors
\begin{equation}
\tilde{\bf R}_x =
\frac{1}{N} \sum_{k=1}^{N} {\bf x}_k {\bf x}_k^T =
\frac{1}{N} {\bf X}_N {\bf X}_N^T
\end{equation}
where ${\bf X}_N=\left[{\bf x}_1\,{\bf x}_2\,\ldots\,{\bf x}_N\right]$
is the input data matrix.
In this context, there is a straightforward connection between
PCA and {\em Singular Value Decomposition}\footnote{
The SVD Theorem states that for any matrix ${\bf X}$,
there exist two orthonormal matrices ${\bf U}$ and ${\bf V}$,
and a (nonnegative) diagonal matrix ${\bf \Sigma}$, such that
${\bf X} = {\bf U} {\bf \Sigma} {\bf V}^T$ \cite{Eckart39,Golub89}} :
If
\[
{\bf \bar{X}}_N \equiv \frac{1}{\sqrt{N}} {\bf X}_N
= {\bf U}_N {\bf \Sigma}_N {\bf V}_N^T,
\]
is the SVD of ${\bf \bar{X}}_N$, then the Principal Components
of the stochastic process $\{{\bf x}_k\}$
are the columns of ${\bf U}_N$ as $N \rightarrow \infty$.
The diagonal elements $\sigma_i$, of
${\bf \Sigma}_N$ are called the {\em singular values}
of ${\bf \bar{X}}_N$ and are equal (in the limit) to
the square roots of the eigenvalues $\lambda_i$.

The applications of Principal Component Analysis and SVD
arise in modeling of control systems \cite{Kung78},
in high-resolution spectral estimation \cite{Schmidt79,Tufts82,Kung83,Roy86b},
in antenna array processing \cite{Roy86a},
as well as in training multilayer perceptrons \cite{Linsker88}.
The PCA approach can result in improved performance, e.g. accuracy
or resolution of the target signals. More on these topics
is discussed in Chapter \ref{ch:Applications}.

Principal Component Analysis is traditionally done in
a batch mode: first the data are collected and then
an estimate of ${\bf R}_x$ is computed
(very rarely ${\bf R}_x$ is known in advance)
followed by an eigenvalue decomposition.
In contrast, neural network
techniques are adaptive (recursive) by nature, in other words
the components are estimated every time a new input sample
comes along rather than waiting for all the data to come in before
attempting analysis.
Adaptive techniques in principle work for infinitely long
sequences and require less storage than batch methods
since the data are used only as they arrive and need
not be remembered for the future.
This is useful when one wants to track
the principal component subspaces in a slowly changing environment
without having to do the eigenvalue decomposition
every time from scratch.
In the following we'll discuss different neural-net approaches
for the PC extraction problem and we shall propose  
a new neural model which is applicable to the pure
PC problem as well as the extension to this problem which we call
Constrained PC (see Section~\ref{sec:CPC}).  
Other extensions will be discussed in later chapters.

\subsection{The Hebbian Rule and PC Extracting NN's}

In his seminal work, ``{\em The Organization of Behavior}'' \cite{Hebb49}
Donald Hebb proposed the following rule for adjusting the synaptic
weights between neurons:
\begin{quote}
``When unit $A$ and unit $B$ are simultaneously excited,
increase the strength of the connection between them.''
\end{quote}
Although this rule is quite vague to be useful for practical
applications, a lot of artificial neural training methods are based on the
above philosophy or an extension of it. A natural extension,
proposed by McClelland, Rumelhart and Hinton \cite{McClelland86} is:
\begin{quote}
``Adjust the strength of the connection between units $A$ and $B$
in proportion to the product of their simultaneous activation.''
\end{quote}
This extension covers the general case of the neurons having positive
and negative activation values. It came to be known as the {\em Hebbian
rule} although it is not exactly what Hebb originally proposed.
In Hebb's philosophy the synaptic weights are a measure of the
correlation between the activations of the units: if $A$ and $B$
are often simultaneously positive or simultaneously negative then their connection
weight should be excitory (positive). However if their activation levels are
often of opposite sign then their connection should be inhibitory (negative).

Interestingly, extensions of the Hebbian rule
are in the heart of most neural net approaches
proposed recently for solving the standard PC problem
\cite{Oja82,Sanger89a,Sanger89b,Foldiak89,Kung90a}.
In 1982 \cite{Oja82}, Oja proposed a linear single unit neural
network (see Fig.~\ref{fig:Oja_net}) that extracts just the first
principal component adaptively from a given vector stochastic process
$\{ {\bf x}_k \}$.  The output value $y_k$, is linearly related
with the synaptic weights which form the vector ${\bf w}$:
\begin{equation}
y_k = {\bf w}_k^T {\bf x}_k.
\label{eq:y.eq.wx}
\end{equation}
The proposed adaptive learning rule (so called {\em the normalized
Hebbian}, or {\em Oja rule}) is
\begin{equation}
\begin{array}{c}
{\bf w}_{k+1} = {\bf w}_k + \beta (y_k {\bf x}_k - y_k^2 {\bf w}_k),
\\
\mbox{any initial } {\bf w}_0 \neq 0.
\end{array}
\label{eq:Oja_rule}
\end{equation}

It was proved that under stationarity assumption,
${\bf w}_k$ converges to the normalized principal
eigenvector ${\bf e}_1$,
of ${\bf R}_x$ as ~$k\!\rightarrow\!\infty$.
The method however cannot be trivially extended to extract
more than one components: by just replicating this network,
each of the neurons will simply
extract independently the first component. It is clear that some
kind of interaction between the units is necessary if orthogonal
components are to be extracted.

Oja's rule is also a linearized version of the normalized Hebbian
rule, the Hebbian part being the input-output product
$y_k{\bf x}_k$, while the normalization part corresponds to
the term $-y_k^2{\bf w}_k$. Indeed observe that
in its simplest form, the realization of the Hebbian rule is
\begin{equation}
{\bf w}_{k+1} = {\bf w}_k + \beta (y_k {\bf x}_k)
\label{eq:Hebbian_rule}
\end{equation}
which is equal to Oja's rule without the $-y_k^2{\bf w}_k$ part.
However Eq.~\ref{eq:Hebbian_rule} is numerically unstable
when ${\bf x}_k$ and $y_k$ are related by Eq~\ref{eq:y.eq.wx}.
Actually, in this case the Hebbian rule becomes
\[
{\bf w}_{k+1}=({\bf I}+{\bf x}_k{\bf x}_k^T){\bf w}_k
\]
where one eigenvalue of the matrix $({\bf I}+{\bf x}_k{\bf x}_k^T)$
is always larger than 1 (for any $k$) while all the other eigenvalues
are equal to 1. Thus the norm of ${\bf w}$ keeps increasing
monotonically making the system unstable.
The problem can be combatted though by
the {\em normalized Hebbian rule} defined below:
\begin{equation}
\tilde{\bf w}_{k+1} = {\bf w}_k + \beta (y_k {\bf x}_k),
\label{eq:normhebb.1}
\end{equation}
\begin{equation}
{\bf w}_{k+1} = \frac{\tilde{\bf w}_{k+1}}{\|\tilde{\bf w}_{k+1}\|}~.
\label{eq:normhebb.2}
\end{equation}
The normalization at each step is vital for preventing instability,
while it is a common feature in many standard
eigenvalue-related algorithms for example,
the power iteration method
\cite[chapter 7]{Golub89}\cite[chapter 4]{Gourlay73}.
The equivalence of Oja's updating rule (Eq.~\ref{eq:Oja_rule})
and the normalized Hebbian rule (Eqs.~\ref{eq:normhebb.1},
\ref{eq:normhebb.2})
can be established through a {\em linearization} procedure:
from Eqs.~\ref{eq:normhebb.1} and \ref{eq:normhebb.2} we obtain
\[
\| \tilde{\bf w}_{k+1} \|^2 = 1 + 2 \beta y_k^2 + {\cal O}(\beta^2),
\]
\[
{\bf w}_{k+1} = \tilde{\bf w}_{k+1} (1 + 2 \beta y_k^2 +
{\cal O}(\beta^2))^{-1/2}~.
\]
Approximating the inverse square root by a
first order Taylor expansion and ignoring $\beta^2$ terms we get
\[
{\bf w}_{k+1} = \tilde{\bf w}_{k+1} (1 - \beta y_k^2).
\]
Substituting Eq.~\ref{eq:normhebb.1} for $\tilde{\bf w}_{k+1}$ we have
\[
{\bf w}_{k+1} = {\bf w}_k + \beta (y_k {\bf x}_k - y_k^2 {\bf w}_k) -
\beta^2 y_k^3 {\bf x}.
\]
Finally, ignoring the $\beta^2$
term we get the original Oja rule, hence we established
the aforementioned equivalence.

\begin{figure}
\centerline{\psfig{figure=oja.ps,height=3in}}
\caption{Oja's simplified neuron model.}
\label{fig:Oja_net}
\end{figure}

An attempt to extend to the multi-component extraction problem has
been carried out by Foldiak \cite{Foldiak89} in 1989. The proposed
model has $m$ output neurons with feedforward connections from input to
output and full connectivity of the output neurons.
The system learns the same space
spanned by the $m$ principal components, but the final vectors are not
the components themselves.  The disadvantages of the model are the
fact that it has only been proved by simulation,
and that it is less efficient than other algorithms
discussed below.

In 1989 \cite{Sanger89a} Sanger proposed another method which is
more attractive than Foldiak's and is equal to Oja's method
for a single component. The model has $m$ output neurons without any
lateral connections between them. There are only feedforward
connections between input and output
and the updating equations for neuron $i$
($i=1,\ldots,m$) are as follows
\begin{eqnarray}
{\bf w}_{i,k+1} &=& {\bf w}_{ik}+\beta (y_{ik} {\bf x}_k -
y_{ik} \sum\limits_{j \leq i} y_{jk} {\bf w}_{jk})
\label{eq:Sanger}
\\
y_{ik} &=& {\bf w}_{ik}^T {\bf x}_k
\end{eqnarray}
The model has been proved to extract the $m$ principal components of
the input stochastic process. The method has been criticized for using
non-local information for the adaptation of the weights without
specifically putting connections between the neurons \cite{Foldiak89}. If an analogy
has to be drawn with biological neural nets this is a disadvantage.
Sanger tried to overcome this criticism by introducing a local
version of his algorithm \cite{Sanger89b},
which is actually Oja's model followed
by a well known {\em deflation} procedure (see Section~\ref{sec:APEX}).

None of the previous methods can effectively support
a recursive approach for the calculation of the $m$-th principal
component given the first $m\!-\!1$ components.
The motivation behind such an approach is the need
to extract the principal components of a random vector sequence
when the number of required PC's is not known a priori. It is also
useful in environments where ${\bf R}_x$ might be slowly changing with
time (e.g. speech applications). Then the new PC may be added
trying to compensate for that change without affecting
the previously computed PCs. This is similar to the idea
of lattice filtering used in signal processing
-- where for every increased filter order,
one new  lattice filtering section is added to the original structure
but all the old sections remain completely intact.
This idea leads to a new neural model called {\em APEX}, an acronym
standing for {\em Adaptive Principal-component EXtractor}
\cite{Kung90a}.
APEX amends the problems associated with other NN
models described above and at the same time gains a
large advantage in convergence speed along with a considerable
advantage in terms of number of operations per sweep. Furthermore,
the method is more powerful in that it can be applied to the
Constrained PC problem discussed later in this Chapter.
This method is described below.


\section{A New Neural Network called APEX}
\label{sec:APEX}

APEX is a neural PC extractor structure which is able to extract
the $m$-th component assuming that all the previous components
have been stored in the synaptic weights of the network.
The APEX model is depicted in Figure \ref{fig:APEX_model}.
There are $n$ inputs $\{x_1 \ldots x_n\}$,
connected to $m$ outputs $\{y_1 \ldots y_m\}$,
through the feedforward weights $\{w_{ij}\}$.
Additionally, there exist lateral
weights $c_j$, forming vector $\bf c$, that
connect the first $m-1$ units with the $m$-th one.
These connections play a very important role in the model:
they work toward the orthogonalization of the synaptic weights of the
$m$-th neuron versus the extracted principal components
stored in the weights of the previous $m\!-\!1$ neurons. Eventually,
all the weights in the network will span an orthonormal base of the
$m$-dimensional principal component subspace. Interestingly enough,
the lateral weights tend to zero as orthogonalization is achieved.
Thus, we associate the term ``anti-Hebbian connections''
or ``orthogonalization connections'' with $c_j$, in contrast
to the term ``Hebbian connections'' associated with $w_{ij}$.
The mathematical details of the anti-Hebbian rule
(also called {\em Lateral Orthogonalization Rule}) used to train
$c_j$ are to follow.

We assume that the input ${\bf x}_k$, $k=0,1,2,\ldots$,
is a stationary stochastic process
with positive-definite autocorrelation matrix ${\bf R}_x$.
The most important feature of APEX hinges upon the fact that
the $m$-th neuron is able to extract the largest component
which is orthogonal to the first $m\!-\!1$ components
represented by the already trained $m\!-\!1$ neurons.
The activation of each neuron is a linear function
of its inputs
\begin{equation}
{\bf y} = {\bf W} {\bf x}
\label{eq:APEX_y}
\end{equation}
\begin{equation}
y_m = {\bf w}^T {\bf x} - {\bf c}^T {\bf y}
\label{eq:APEX_y.eq.wx-cy}
\end{equation}
where ${\bf x} = [ x_1 \ldots x_n ]^T$,
${\bf y} = [ y_1 \ldots y_{m-1} ]^T$, ${\bf W}=[{\bf e}_1 \ldots
{\bf e}_{m-1}]^T$ is the weight matrix
for the first $m\!-\!1$ neurons, and ${\bf w}$ is the
vector comprised of the synaptic weights $w_i$,
connecting the inputs with the $m$-th neuron.

\begin{figure}
\centerline{\psfig{figure=apex.ps,height=3in}}
\caption{The linear multi-output model. The solid lines denote the
weights \protect{$w_i$}, \protect{$c_j$}
which are trained at the \protect{$m$}-th stage.
The dashed lines correspond to the weights of the already
trained neurons. Note that the lateral weights asymptotically
converge to zero so they don't appear between the already
trained units.}
\label{fig:APEX_model}
\end{figure}

Only ${\bf w}$ and ${\bf c}$ are trained (${\bf W}$ remains fixed).
The $k$-th iteration of the algorithm is
\begin{eqnarray}
{\bf w}_{k+1} &=& {\bf w}_k + \beta_k (y_{mk} {\bf x}_k - y_{mk}^2 {\bf w}_k)
\label{eq:APEX_dw}
\\
{\bf c}_{k+1} &=& {\bf c}_k + \beta_k (y_{mk} {\bf y}_k - y_{mk}^2 {\bf c}_k)
\label{eq:APEX_dc}
\end{eqnarray}
where $\beta_k$ is a positive sequence of step-size parameters.
If we expand the above equations for each
individual weight we get the following formulas
$$
\Delta w_{ik} = \beta_k (y_{mk} x_{ik} - y_{mk}^2 w_{ik}),
~~~i = 1 \ldots n
$$
$$
\Delta c_{ik} = \beta_k (y_{mk} y_{ik} - y_{mk}^2 c_{ik}),
~~~i = 1 \ldots m-1
$$

Notice that both equations are of the Oja type (see Eq.~\ref{eq:Oja_rule}).
The difference
lies on the sign with which ${\bf w}$ and ${\bf c}$
participate in Eq.~\ref{eq:APEX_y.eq.wx-cy}.
The ${\bf w}$ weights intuitively work towards finding most principal
directions of the signal, that's why we call Eq.~\ref{eq:APEX_dw}
the {\em Hebbian part} of the algorithm.
The $\bf c$ weights intuitively
subtract the first $m\!-\!1$ components from the
$m$-th neuron, thus the $m$-th output neuron tends to become
orthogonal (rather than be correlated) to all the previous components.
We call Eq.~\ref{eq:APEX_dc} the {\em orthogonalization rule}
which constitutes the {\em anti-Hebbian part} of the algorithm.
It will be shown below that the combination
of the Hebbian and orthogonalization rules works as described.
This will be proved by the
theoretical analysis in Section \ref{sec:APEX-theory},
and demonstrated by simulations
in Section \ref{sec:APEX-simulations}.
The orthogonal learning rule also has a very important
application to the problem of extracting ``constrained principal
components'' as discussed in Section \ref{sec:CPC} and
further elaborated in \cite{Kung90b}.

\section{APEX Proof and Convergence Rates}
\label{sec:APEX-theory}

In this section we will prove the main convergence theorem
related to APEX. First we state some definitions and results of
matrix theory which will be useful for our analysis.
\begin{definition} \cite{Golub89}
The mapping
\[
{\bf R}_x \rightarrow \tilde{\bf R}_x =
({\bf I} - {\bf e}_1 {\bf e}_1^T) {\bf R}_x ({\bf I} - {\bf e}_1 {\bf e}_1^T)
\]
is called a {\bf deflation transformation} of ${\bf R}_x$.
\end{definition}

{\bf Fact:} The eigenvectors of \label{page:fact}
$\tilde{\bf R}_x$ are ${\bf e}_1, {\bf e}_2,
\ldots, {\bf e}_n$, as for ${\bf R}_x$, but the corresponding eigenvalues
are $0,\lambda_2,\ldots,\lambda_n$.
Also $({\bf I}-{\bf e}_1{\bf e}_1^T){\bf R}_x = \tilde{\bf R}_x$.

Indeed, since the eigenvectors are orthonormal we have
\[
\tilde{\bf R}_x {\bf e}_i =
\left\{ \begin{array}{ll} 0 ~~~\mbox{if } i=1 & \\
\lambda_i({\bf I}-{\bf e}_1{\bf e}_1^T){\bf e}_i =
\lambda_i {\bf e}_i & \mbox{if } i>1 \end{array} \right.
\]
In other words the principal eigenvalue
is ``killed'' while the second component
becomes now the dominant one.
All components smaller than the second one (inclusive)
keep their eigenvalues intact.
Also notice that
\[
({\bf I}-{\bf e}_1{\bf e}_1^T){\bf R}_x {\bf e}_i =
\lambda_i({\bf I}-{\bf e}_1{\bf e}_1^T){\bf e}_i =
\left\{ \begin{array}{ll} 0 & \mbox{if } i<m \\
\lambda_i {\bf e}_i & \mbox{if } i\geq m \end{array} \right.
\]
i.e. the eigenvalue/vector pairs of $\tilde{\bf R}_x$ are also
eigenvalue/vector pairs of $({\bf I}-{\bf e}_1{\bf e}_1^T) {\bf R}_x$
thus the two matrices must be equal:
\begin{equation}
\tilde{\bf R}_x
=({\bf I}-{\bf e}_1{\bf e}_1^T) {\bf R}_x ({\bf I}-{\bf e}{\bf e}_1^T)
=({\bf I}-{\bf e}_1{\bf e}_1^T) {\bf R}_x.
\end{equation}
This also implies that
$({\bf I}-{\bf e}_1{\bf e}_1^T){\bf R}_x$ is symmetric.

Inductively, we can show that the $m\!-1\!$ times deflation
of ${\bf R}_x$
\begin{eqnarray}
\tilde{\bf R}_x &=& \prod_{i=1}^{m-1} ({\bf I}-{\bf e}_i{\bf e}_i^T)
{\bf R}_x \prod_{i=1}^{m-1} ({\bf I}-{\bf e}_i{\bf e}_i^T) \nonumber
\\
&=& \prod_{i=1}^{m-1} ({\bf I}-{\bf e}_i{\bf e}_i^T) {\bf R}_x
\label{eq:PCA_deflation_m}
\end{eqnarray}
``kills'' the $m\!-\!1$ most principal eigenvalues so the
eigenvalue/vector pairs now are $(0,{\bf e}_1)$,$\ldots$,
$(0,{\bf e}_{m-1})$,$(\lambda_m,{\bf e}_m)$,$\ldots$,$(\lambda_n,{\bf e}_n)$.
Clearly the $m$-th component is now dominant. Because of the
orthogonality property ${\bf e}_i^T{\bf e}_j=0$ for $i\neq j$,
Eq.~\ref{eq:PCA_deflation_m} can be rewritten as follows
\begin{eqnarray}
\tilde{\bf R}_x &=& ({\bf I}-\sum_{i<m} {\bf e}_i{\bf e}_i^T)
{\bf R}_x ({\bf I}-\sum_{i<m} {\bf e}_i{\bf e}_i^T) \nonumber
\\
&=& ({\bf I}-{\bf W}^T{\bf W})
{\bf R}_x ({\bf I}-{\bf W}^T{\bf W})
\\
&=& ({\bf I}-{\bf W}^T{\bf W}){\bf R}_x
\end{eqnarray}
where ${\bf W}=[{\bf e}_1,\ldots,{\bf e}_{m-1}]^T$.
These facts will be useful in our analysis below.
Next we formally state the assumptions used by the main theorem:

\begin{enumerate}
\item[A.1] The input sequence $\{{\bf x}_k\}$, is at least
wide sense stationary with autocorrelation matrix ${\bf R}_x$,
whose eigenvalues are distinct, positive and arranged
in descending order: $\lambda_1 > \lambda_2 > \ldots > \lambda_n > 0$.
Let us denote the corresponding normalized eigenvectors (choose
arbitrarily the sign) by ${\bf e}_1, {\bf e}_2, \ldots , {\bf e}_n$.
\item[A.2] The first $m\!-\!1$ principal components,
i.e. the eigenvectors ${\bf e}_1$ through ${\bf e}_{m-1}$, are
stored in the synaptic weights connected to neurons $1$ through $m\!-\!1$.
\item[A.3] The step-size parameter sequence $\beta_k$ is such that
\[
\sum_{k=0}^\infty \beta_k = \infty ~~
\mbox{ and } ~~
\beta_k \rightarrow 0, \mbox{ as } k \rightarrow \infty.
\]
This assumption will be useful for showing the asymptotic convergence
of the algorithm. Later we'll discuss the case where it doesn't hold.
\end{enumerate}


\begin{lemma}
Consider the differential equation
\begin{equation}
d {\bf p}/dt = ({\bf I}-{\bf V}^T{\bf V}) {\bf R} {\bf p} -
({\bf p}^T {\bf R} {\bf p}) {\bf p}
\label{eq:lemma_ode}
\end{equation}
where ${\bf R}$ is symmetric and positive definite
and the rows ${\bf V}$
are orthonormal, i.e. ${\bf V}{\bf V}^T={\bf I}$.
If the initial condition is such that $\|{\bf p}(0)\|>0$ then
\begin{enumerate}
\item[(a)] $\|{\bf p}(t)\|$ is bounded
\item[(b)]
${\bf p}(t)^T {\bf R} {\bf p}(t)>0$, for all $t$.
\end{enumerate}
\label{lemma:APEX}
\end{lemma}

{\bf Proof.}
\begin{enumerate}
\item[(a)] Premultiply Eq.~\ref{eq:lemma_ode} by $2{\bf p}^T$ to obtain
\begin{equation}
\frac{d\|{\bf p}\|^2}{dt} =
2\left[ {\bf p}^T({\bf I}-{\bf V}^T{\bf V}){\bf R}{\bf p}
- ({\bf p}^T{\bf R}{\bf p})\|{\bf p}\|^2 \right]
\label{eq:lemma_ode_pp}
\end{equation}
It is easy to show that the quantity
\[
\mu = \max_{\bf p}
\{{\bf p}^T({\bf I}-{\bf V}^T{\bf V}){\bf R}{\bf p}/\|{\bf p}\|^2\}
\]
is positive. Indeed let ${\bf p}=({\bf I}-{\bf V}^T{\bf V}){\bf a}$
for some nonzero vector ${\bf a}$. Then
\[
{\bf p}^T({\bf I}-{\bf V}^T{\bf V}){\bf R}{\bf p}/\|{\bf p}\|^2=
{\bf a}^T({\bf I}-{\bf V}^T{\bf V}){\bf R}({\bf I}-{\bf V}^T{\bf V})
{\bf a}/\|{\bf p}\|^2>0
\]
since ${\bf R}$ is positive definite. Thus $\mu>0$.
Also the quantity
\[
\rho_n=\min_{\bf p}\{{\bf p}^T{\bf R}{\bf p}/\|{\bf p}\|^2\}
\]
equals the minimum eigenvalue of ${\bf R}$ and so $\rho_n>0$.
From Eq.~\ref{eq:lemma_ode_pp} we get
\[
\frac{d\|{\bf p}\|^2}{dt} \leq 2\mu \|{\bf p}\|^2 - 2\rho_n \|{\bf p}\|^4
\]
From the basic theorem of differential inequalities \cite[chapter 2]{Hale69}
it follows that if $u(t)$ is the solution to the equation
\begin{equation}
\frac{du}{dt} = 2\mu u - 2\rho_n u^2
\label{eq:lemma_ode_riccati}
\end{equation}
with initial condition $u(0)=\|{\bf p}(0)\|^2$ then
\[
\|{\bf p}\|^2<u(t), \mbox{ for all } t
\]
However, Eq.~\ref{eq:lemma_ode_riccati} is a stable differential
equation (Riccati) so $\|{\bf p}\|^2$ is bounded.

\item[(b)] Premultiply Eq.~\ref{eq:lemma_ode} by $2{\bf p}^T{\bf R}$
to obtain
\begin{equation}
\frac{d{\bf p}^T{\bf R}{\bf p}}{dt} =
2\left[ {\bf p}^T{\bf R}({\bf I}-{\bf V}^T{\bf V}){\bf R}{\bf p}
- ({\bf p}^T{\bf R}{\bf p})^2 \right]
\label{eq:lemma_ode_pRp}
\end{equation}
The matrix ${\bf R}({\bf I}-{\bf V}^T{\bf V}){\bf R}$
is semi-positive definite so
\[
\frac{d{\bf p}^T{\bf R}{\bf p}}{dt} \leq
-2 ({\bf p}^T{\bf R}{\bf p})^2
\]
If $s$ is the solution to the the equation
\[
\frac{ds}{dt}=-2s^2
\]
with initial condition $s(0)={\bf p}^T{\bf R}{\bf p}>0$
then $0<s(t)\leq {\bf p}(t)^T{\bf R}{\bf p}(t)$ for all $t$.

\end{enumerate}
\hfill \qed

Now we are ready to prove the main theorem

\begin{theorem}
\label{th:APEX}
Consider the algorithm defined in Eqs. \ref{eq:APEX_dw}, \ref{eq:APEX_dc}
subject to the conditions A.1 -- A.3. Then with probability 1,
${\bf w}_k \rightarrow {\bf e}_m$ (or $-{\bf e}_m$)
and ${\bf c}_k \rightarrow 0$, as $k \rightarrow \infty$.
\end{theorem}

{\bf Proof.}
We'll use the now standard approach for analyzing stochastic approximation
type algorithms which was first proposed in \cite{Ljung77} and later in
\cite{Kushner78}. First we rewrite Eqs. \ref{eq:APEX_dw}, \ref{eq:APEX_dc}
\begin{eqnarray}
{\bf w}_{k+1} &=& {\bf w}_k + \beta_k f_1({\bf x}_k,{\bf w}_k,{\bf c}_k)
\\
{\bf c}_{k+1} &=& {\bf c}_k + \beta_k f_2({\bf x}_k,{\bf w}_k,{\bf c}_k)
\end{eqnarray}
where
\begin{eqnarray}
f_1({\bf x},{\bf w},{\bf c}) &=&
{\bf x} {\bf x}^T {\bf w} -
{\bf x} {\bf x}^T {\bf W}^T {\bf c}
- \bigl[
({\bf w} -{\bf W}^T {\bf c})^T {\bf x} {\bf x}^T
({\bf w} -{\bf W}^T {\bf c}) \bigr]
{\bf w}
\\
f_2({\bf x},{\bf w},{\bf c}) &=&
{\bf W} {\bf x} {\bf x}^T {\bf w} -
{\bf W} {\bf x} {\bf x}^T {\bf W}^T {\bf c}
- \bigl[
({\bf w} -{\bf W}^T {\bf c})^T {\bf x} {\bf x}^T
({\bf w} -{\bf W}^T {\bf c}) \bigr] {\bf c}
\end{eqnarray}

Then for both discrete-time dynamic equations
there is an associated differential equation defined as follows
\begin{eqnarray}
\frac{d {\bf w}'(t)}{dt} &=&
E_x \{ f_1({\bf x},{\bf w}'(t),{\bf c}'(t)) \} = \nonumber
\\
&=& {\bf R}_x {\bf w}'(t) - {\bf R}_x {\bf W}^T {\bf c}'(t)
- \sigma'(t) {\bf w}'(t)
\label{eq:APEX_ode_w}
\\
\frac{d {\bf c}'(t)}{dt} &=&
E_x \{ f_2({\bf x},{\bf w}'(t),{\bf c}'(t)) \} = \nonumber
\\
&=& {\bf W} {\bf R}_x {\bf w}'(t) - {\bf W} {\bf R}_x {\bf W}^T {\bf c}'(t)
- \sigma'(t) {\bf c}'(t)
\label{eq:APEX_ode_c}
\end{eqnarray}
where
\[
\sigma'(t) \equiv {\bf q}'(t)^T {\bf R}_x {\bf q}'(t)
\]
and
\[
{\bf q}'(t) = {\bf w}'(t) - {\bf W}^T {\bf c}'(t).
\]
The time index $t$, of the differential equations and the index $k$,
of the discrete-time equations are related by the formula
\begin{equation}
t(k) = \sum_{\ell=0}^k \beta_\ell
\label{eq:t_k}
\end{equation}
The values $[{\bf w}_k,{\bf c}_k]$ approach the corresponding values
$\Bigl[ {\bf w}'\bigl(t(k)\bigr),{\bf c}'\bigl(t(k)\bigr)\Bigr]$,
as $k$ (and hence $t$) get very large \cite{Ljung77,Ljung83}
under the additional assumption that the invariant set
$D_c = \{({\bf e}_m,0), (-{\bf e}_m,0)\}$ for
the system of differential equations \ref{eq:APEX_ode_w}, \ref{eq:APEX_ode_c},
is asymptotically stable and that $({\bf w}_k,{\bf c}_k)$
belongs to its domain of attraction $D_A$, infinitely often
with probability 1.
To that end, we shall show that $D_A$ equals the whole space ${\cal R}^n$,
except for some lower dimensional invariant subspaces
thus their probability measure is 0. Let us expand
\[
{\bf w}'(t) = \sum_{i=1}^n \theta_i(t) {\bf e}_i
\]
where ${\bf e}_i$ is the normalized $i$-th eigenvector of ${\bf R}_x$
(any one of the two).
Given that the rows of ${\bf W}$ 
are ${\bf e}_1, \ldots, {\bf e}_{m-1},$
Eqs. \ref{eq:APEX_ode_w} and \ref{eq:APEX_ode_c}
can be written as
\begin{eqnarray}
\left[ \begin{array}{l} \frac{d \theta_i(t)}{dt} \\
\frac{d c'_i(t)}{dt}
\end{array} \right] &=&
\left[ \begin{array}{cc} \lambda_i-\sigma'(t) & -\lambda_i \\
\lambda_i  & -\lambda_i-\sigma'(t) \end{array} \right]
\left[ \begin{array}{l} \theta_i(t) \\ c'_i(t) \end{array} \right],
~~~ i=1, \ldots ,m-1
\label{eq:APEX_ode_1tom-1}
\\
\frac{d \theta_i(t)}{dt} &=& \bigl[ \lambda_i-\sigma'(t) \bigr] \theta_i(t),
~~~ i=m, \ldots ,n
\label{eq:APEX_ode_mton}
\end{eqnarray}

Eq. \ref{eq:APEX_ode_1tom-1} has a double eigenvalue at ~$-\sigma'(t)
={\bf q}'(t)^T {\bf R}_x {\bf q}'(t)$.
Since ${\bf R}_x$ is positive definite $\sigma'(t) > 0$ at all times,
unless ${\bf q}'(t)=0$ for some $t$. However, we'll show that
${\bf q}'(t) \neq 0$, $\forall t \geq 0$. Indeed,
premultiplying Eq. \ref{eq:APEX_ode_c} by ${\bf W}^T$, subtracting from
Eq. \ref{eq:APEX_ode_w} we obtain
\begin{equation}
\frac{d {\bf q}'(t)}{dt} = ({\bf I} - {\bf W}^T{\bf W}) {\bf R}_x {\bf q}'(t)
- {\bf q}'(t)^T {\bf R}_x {\bf q}'(t) {\bf q}'(t)
\label{eq:APEX_ode_q}
\end{equation}
From Lemma \ref{lemma:APEX}, $\sigma'(t)>0$ for all $t$
so $\theta_i\rightarrow 0$, $i<m$ and ${\bf c}'(t)\rightarrow 0$
as $t\rightarrow \infty$.

For the components $i>m$, we have
\begin{eqnarray*}
\frac{d}{dt} \left( \frac{\theta_i}{\theta_m} \right) &=&
\frac{1}{\theta_m} \frac{d \theta_i}{dt} -
\frac{\theta_i}{\theta_m^2} \frac{d \theta_m}{dt}
\\
&=& [\lambda_i - \sigma'(t)] \frac{\theta_i(t)}{\theta_m(t)}
- [\lambda_m - \sigma'(t)] \frac{\theta_i(t)}{\theta_m(t)}
\\
&=& [\lambda_i - \lambda_m] \frac{\theta_i(t)}{\theta_m(t)}
\end{eqnarray*}
Since ~$\lambda_i < \lambda_m$ it follows that
\begin{equation}
\frac{\theta_i(t)}{\theta_m(t)} \rightarrow 0,
\mbox{ as } t \rightarrow \infty
\label{eq:ratio_i_goto_0}
\end{equation}
From the boundedness of ${\bf q}'(t)$ and the fact that ${\bf c}'(t)
\rightarrow 0$ it follows that $\theta_m(t)$ is bounded,
therefore Eq.~\ref{eq:ratio_i_goto_0} implies that $\theta_i\rightarrow 0$,
$i>m$. It remains to show that
$\theta_m(t)$ does not vanish as $t \rightarrow \infty$.
Indeed, Eq.~\ref{eq:APEX_ode_mton} for $i=m$ yields
\[
\frac{d\theta_m}{dt}=[\lambda_m-\sigma'(t)]\theta_m
\]
So unless $\theta_m(0)=0$, the magnitude of $\theta_m(t)$ increases
when $\sigma'(t)$ becomes less than $\lambda_m$ thus $\theta_m(t)$
cannot converge to 0. Let us denote the final
nonzero convergence value of $\theta_m$ by $\alpha$.
Plugging the steady state values ${\bf w}'(t)=\alpha {\bf e}_m$
and ${\bf c}'(t)=0$ into Eqs.~\ref{eq:APEX_ode_w}
and \ref{eq:APEX_ode_c} we obtain that $\alpha=\pm 1$.
Therefore under the initial condition
\begin{equation}
{\bf w}'(0)^T {\bf e}_m = \theta_m(0) \neq 0
\label{eq:APEX_ode_initial_condition}
\end{equation}
we have
\begin{equation}
{\bf w}'(t) \rightarrow \pm {\bf e}_m, \mbox{ as } t \rightarrow \infty.
\end{equation}

The domain of attraction of this solution contains the whole space
except for the subspace ${\bf w}^T {\bf e}_m = 0$.
The probability that ${\bf w}_k$ does not belong to the
domain of attraction infinitely often is equal to the probability
that it belongs to this subspace almost always (i.e. except for
finite many times). However the probability of this event is 0
since almost certainly the disturbance $\Delta {\bf w}_k$
has a component orthogonal to this subspace.
Therefore we conclude that
\[
{\bf w}_k \rightarrow \pm {\bf e}_m,
\]
\[
{\bf c}_k \rightarrow 0, ~\mbox{ as } k \rightarrow \infty.
\]
\hfill \qed

The APEX rule for the first unit is equivalent to Oja's
rule, the lateral weights being absent in this case.
We can form the basis of an induction from the fact
that Oja's rule extracts the first principal component while
Theorem \ref{th:APEX} provides the induction step by
stating that given the previous $m\!-1\!$ components
the $m$-th one is extracted by the APEX network.
From the complete induction argument it follows that
the APEX algorithm can sequentially extract all the principal
components $\{{\bf x}_k\}$ provided that infinite time
can be used by the network to converge to the exact solutions.
However, in practice finite amount of time is available
for training the model and the algorithm must stop
without extracting the exact components. That affects the
accuracy of the later units which rely on the accuracy
of the previous components for the corresponding deflation.
So errors can accumulate as units are added. However through
our simulation experience we observed that the first 8--16 components
are usually estimated with quite reasonable accuracy
and without excessive training time.

Although Theorem \ref{th:APEX} assumes that $\beta_k \rightarrow 0$~
as $k \rightarrow \infty$, for practical purposes
this is usually not the case. In fact the well known trade-off between
tracking capability and accuracy \cite{Ljung83} is present here as well
as any other adaptive algorithm. One wants to keep the step-size small,
but non-zero, in order to be able to follow (slow) changes in the statistics
of the input signal. It is known \cite{Kushner81,Benveniste82} that
even if $\beta_k$ does not go to zero
but remains equal to a small constant
then the mean values $\bar{\bf w}_k$ of ${\bf w}_k$
and $\bar{\bf c}_k$ of ${\bf c}_k$ still approximate
the ordinary differential equations \ref{eq:APEX_ode_w}
and \ref{eq:APEX_ode_c} therefore the algorithm still
converges (in the mean) to the solution.
However now in the steady state,
the adapted parameters oscillate around the fixed
point without ever actually converging to it.
In fact the final error approximates a Gaussian
process which is the solution of a {\em diffusion
stochastic differential equation} driven by a standard Brownian motion.
This phenomenon is due to
the variance of the stochastic input that drives the
system. Naturally the amplitude of the oscillation
is proportional to this variance as well as to the size of $\beta$.
However in many practical applications, a small uncertainty
in the final value is tolerable, traded-off for a
valuable property of the algorithm, namely the capability
to track the changing parameters of the observed (input) process
and adapt its internal parameters to these changes.

The above observations hold true for any adaptive algorithm
in general, so they hold true for APEX as well.
Nevertheless, when discussing APEX specifically,
a lot more can be said about the step-size parameter
as it will become clear next.

\subsection{Relationship with RLS}

We saw before that
principal component analysis can be formulated as
mean square minimization problem.
Naturally then, we pursue the investigation
of the relationship between APEX
and the Recursive Least Squares (RLS) algorithm
which is well studied in the literature. By establishing this relationship
we will extract some valuable conclusions that relate to the value of the
step-size constant as well as
the numerical performance of the algorithm.
The first principal component optimizes the reconstruction mean-square-error
by picking the best 1-dimensional subspace to project the signal onto.
The criterion for the first component is
\begin{equation}
\mbox{minimize }
J_1({\bf w}) = E \{ \| {\bf x} - {\bf w} {\bf w}^T {\bf x} \|^2 \}
\end{equation}
Expanding $J_1$ and doing some mathematical manipulations we obtain
\begin{eqnarray}
J_1 &=& E \{ \| {\bf x} \|^2 \} - 2 E \{ 
{\bf w}^T {\bf x} {\bf x}^T {\bf w} \} + trace \left(
E \{ {\bf w} {\bf w}^T {\bf x} {\bf x}^T {\bf w} {\bf w}^T \} \right) = \nonumber
\\
&=& E \{ \| {\bf x} \|^2 \} - 2 {\bf w}^T {\bf R}_x {\bf w}
+ trace \left( {\bf w} {\bf w}^T {\bf R}_x {\bf w} {\bf w}^T \right)
\label{eq:1-component-cost}
\end{eqnarray}
The solution to the above minimization problem is ${\bf w} = \alpha {\bf e}_1$,
where $\alpha\neq 0$ is any scalar.

For the $m$-th component the criterion becomes
\begin{equation}
\mbox{minimize }
J_m({\bf w}) = E \{ \| {\bf x} - {\bf w} {\bf w}^T
({\bf I} - {\bf W}^T {\bf W}) {\bf x} \|^2 \}
\end{equation}
where ${\bf W} = [ {\bf e}_1 \ldots {\bf e}_{m-1}]^T$.
Indeed, expanding $J_m$ as we did for $J_1$ we get
\begin{eqnarray}
J_m &=& E \{ \| {\bf x} \|^2 \}
- 2 {\bf w}^T ({\bf I} - {\bf W}^T {\bf W}) {\bf R}_x {\bf w} \nonumber
\\
&& + trace \left( {\bf w} {\bf w}^T ({\bf I} - {\bf W}^T {\bf W})
{\bf R}_x ({\bf I} - {\bf W}^T {\bf W}) {\bf w} {\bf w}^T \right)
\label{eq:m-component-cost}
\end{eqnarray}

From the facts stated in page \pageref{page:fact},
Eq. \ref{eq:m-component-cost} becomes
\begin{equation}
J_m = E \{ \| {\bf x} \|^2 \}
- 2 {\bf w}^T \tilde{\bf R}_x {\bf w}
+ trace \left( {\bf w} {\bf w}^T \tilde{\bf R}_x {\bf w} {\bf w}^T \right)
\end{equation}
therefore $J_m$ is the same as $J_1$ for the deflated matrix
$\tilde{\bf R}_x$. So the minimizer of $J_m({\bf w})$
is the principal component of $\tilde{\bf R}_x$ which is  ${\bf e}_m$.

In order to relate the above analysis with RLS
we first transform the cost function $J_m({\bf w})$
into a non-stochastic formulation, replacing the expectation by (weighted)
average
\begin{equation}
\tilde{J}_m({\bf w},N) = \frac{1}{N+1} \sum_{k=0}^{N}
\gamma^{N-k}\,\| {\bf x}_k - {\bf w} y_k \|^2
\label{eq:m-component-RLS-cost}
\end{equation}
where $0 < \gamma \leq 1$ is a forgetting factor chosen by the user,
and $y_k = {\bf w}^T {\bf x}_k - {\bf w}^T {\bf W}^T {\bf W} {\bf x}_k$.
We shall define
\begin{equation}
{\bf c} \equiv {\bf W} {\bf w}
\label{eq:c_Ww}
\end{equation}
to get
\begin{equation}
y_k = {\bf w}^T {\bf x}_k - {\bf c}^T {\bf W} {\bf x}_k.
\end{equation}

$\tilde{J}_m({\bf w},N)$ can be minimized iteratively using the RLS algorithm
\cite{Ljung83} which yields the following updating equations:
\begin{eqnarray}
{\bf w}_{k+1} &=& {\bf w}_k + L_k ({\bf x}_k - {\bf w}_k y_k)
\\
L_k &=& \left[ \sum_{\ell=0}^{k} \gamma^{k-\ell} y_\ell^2 \right]^{-1}
y_k~.
\end{eqnarray}

Letting
\begin{equation}
\beta_k = \frac{1}{\sum\limits_{\ell=0}^{k} \gamma^{k-\ell} y_\ell^2}
\label{eq:APEX_beta1}
\end{equation}
the RLS algorithm becomes
\begin{equation}
{\bf w}_{k+1} = {\bf w}_k + \beta_k ({\bf x}_k - {\bf w}_k y_k) y_k
\label{eq:RLS_dw}
\end{equation}
which is the same as APEX equation \ref{eq:APEX_dw}. The only difference
is that now we have a specific choice for the value of $\beta_k$ which
is optimal in the sense of the criterion $\tilde{J}_m({\bf w},N)$.
Moreover, the optimal choice of the
step-size parameter has a profound impact in
the convergence speed of the algorithm as will be discussed
in Section \ref{sec:APEX-simulations}.
Clearly $\beta$ can also be calculated iteratively by
\[
\beta_k=\frac{\beta_{k-1}}{\gamma-y_k^2\beta_{k-1}}
\]
The value of the forgetting factor $\gamma$,
induces an effective time-window of size
$M = 1/(1-\gamma)$, so for example, with $\gamma=0.99$
the number of previous values $y_k$ which are important in
$\beta$ is effectively 100.
Therefore $\gamma$ should be chosen
to be closer to 1 if larger time
window is to be used for the averaging, and further from 1
when the past values are less important.
In practice we
can use a hard-limiting window of size $M$, instead of an
effective one induced by $\gamma$, setting
\begin{equation}
\beta_k = \frac{1}{\sum\limits_{\ell=k-M+1}^{k} y_\ell^2}~.
\label{eq:APEX_beta2}
\end{equation}
This formula is more natural in the case where the number of samples
is $M$, whence the algorithm uses the same data again and again
in a loop (each iteration of the loop is called a {\em sweep}).
This situation is may occur for example, in classification
applications where the number of samples is limited
and stored in a computer database, the size of which
is known a priori. The formula \ref{eq:APEX_beta2} for $\beta$
is very similar to what was originally proposed for APEX
\cite{Kung90a} although in that paper we proposed to keep
the value of $\beta$ constant within each sweep but still
inversely proportional to the output power.
In \cite{Abbas92} the authors propose
to track the value $1/\sigma_k=1/\sum_{k-M+1}^k y_\ell^2$ by following
a gradient descent on the error surface $J=1/2[1-\beta_k\sigma_k]^2$.
However neither of these methods are as general as our current treatment
since they are not applicable when $M$ is not defined,
while they fail to establish the relationship between $\beta$
and the least squares optimality criterion \ref{eq:m-component-RLS-cost}.

We proceed now
to establish the full connection between APEX and RLS,
by showing that the updating equation for ${\bf c}$ in APEX
can also be derived from recursive least squares optimization.
Indeed, premultiplying Eq.~\ref{eq:RLS_dw} by ${\bf W}$
and using Eq.~\ref{eq:c_Ww} we obtain
\begin{eqnarray}
{\bf c}_{k+1} &=& {\bf c}_k + \beta_k
({\bf W} {\bf x}_k - {\bf W} {\bf w}_k y_k) y_k = \nonumber
\\
&=& {\bf c}_k + \beta_k ({\bf y}_k - {\bf c}_k y_k) y_k
\label{eq:RLS_dc}
\end{eqnarray}
which is the same as Eq.~\ref{eq:APEX_dc}.
Thus both updating equations defining APEX can be derived
via the RLS theory.
%The RLS theory will
%be very useful when we'll  discuss the constrained principal
%component problem later as well.

\subsection{Rates of convergence}

As we discussed in the proof of Theorem \ref{th:APEX},
the decay rate for ${\bf c}'(t)$ as well as 
$\theta_i(t)$, $i=1,\ldots,m\!-\!1$, is $\sigma'(t)$.
Given that the discrete time parameters approach the continuous
time ones with the relationship between the two time indices
defined in Eq.~\ref{eq:t_k} we can approximate the decay rate
for the discrete-time parameters
\begin{eqnarray}
\frac{{\bf c}_{k+1} - {\bf c}_k}{t(k+1) - t(k)} &=&
- \sigma'(t(k)) {\bf c}_k \Rightarrow \nonumber
\\
{\bf c}_{k+1} &=& [1 - \beta_k \sigma_k] {\bf c}_k
\label{eq:discr_de_c}
\end{eqnarray}
where $\sigma_k = \sigma'(t(k)) = E\{y_k^2\}$. Similarly,
\begin{eqnarray}
\theta_{i,k+1} &=& [1 - \beta_k \sigma_k] \theta_{ik},
~~~~~~ i=1, \ldots, m-1
\label{eq:discr_de_theta_1tom-1}
\\
\theta_{i,k+1} &=& [1 - \beta_k (\sigma_k - \lambda_i)] \theta_{ik},
~~~ i=m, \ldots, n
\label{eq:discr_de_theta_mton}
\end{eqnarray}
From Eq.~\ref{eq:discr_de_theta_mton} we can derive the rate
of decay of the component ratio
\[
r_i=\theta_i/\theta_m, ~~i=m+1,\ldots, n
\]
\begin{equation}
r_{i,k+1} = \frac{1 - \beta_k (\sigma_k - \lambda_i)}
{1 - \beta_k (\sigma_k - \lambda_m)} r_{ik}
\end{equation}

The above theoretical predictions are justified as shown in
Fig.~\ref{fig:APEX_decay_rates} where a comparison is carried out with
simulation results.

\begin{figure}
\centerline{\psfig{figure=n1.ps,width=5in}}
\caption{Comparison of the theoretical vs. the actual decay rates
via simulation.}
\label{fig:APEX_decay_rates}
\end{figure}

%\begin{figure}
%\vspace{1.8 in}
%\caption{ The convergence speed in terms of the number of sweeps 
%vs.~the computed eigenvalues.}
%\label{sweep}
%\end{figure}


\section{Computational Efficiency and Simulation Results}
\label{sec:APEX-simulations}

Based on the above theoretical analysis, the APEX Algorithm is
formally presented below: 

\begin{algorithm}[APEX, Sequential Version]

For every neuron $m=1$ to $N$, ($N \leq n$)
\begin{enumerate}
\item Initialize ${\bf w}$ and ${\bf c}$ to some random values.

\item Choose $\beta$ either as in Eq.~\ref{eq:APEX_beta1}
or Eq.~\ref{eq:APEX_beta2} selecting
$\gamma$ or $M$ appropriately (see Section \ref{sec:APEX-theory}).

\item Update ${\bf w}_k$ and ${\bf c}_k$ according to Eqs.~\ref{eq:APEX_dw},
\ref{eq:APEX_dc}, until a stopping criterion is met.
\end{enumerate}
\end{algorithm}

One such criterion could be that
the average square output $\sigma_k$, approximately converges to a
fixed value. This value is the $m$-th eigenvalue of ${\bf R}_x$
and the corresponding weights ${\bf w}$ form the $m$-th
normal eigenvector. Another viable criterion is detecting
that $\| {\bf c}_k \|$ is below a certain threshold. This criterion
cannot be used though for the first neuron because of the absence
of the lateral weights in this case. On the other hand if someone
wishes to track slow changes in the principal component subspaces
then probably no terminating criterion is specified but rather
the system is left to follow the changes.

Since the neurons in the APEX model are hierarchically structured
each node does not get affected by any nodes following it.
If the neurons prior to a node have converged reasonably close
to the appropriate components then the node would converge to
its corresponding component. So we heuristically contend that if
we let all the nodes work in parallel following the APEX rule
the network will extract the principal components in parallel
rather than one after the other. Indeed, the first node
is unaffected since it has no prior nodes so it will extract
the first component. In turn, the second neuron will start
converging to the second component no later than the first
one converges to the first component assuming that the final
variance of the steady state of the first neuron
does not affect the convergence of the second neuron.
Similarly the $m$-th neuron will start converging to the
$m$-th component no later than the $m\!-\!1$ neuron has converged.
We expect however the final variance of the $m$-th neuron
to increase as $m$ increases accumulating the variances
of all the previous units.
The network for the parallel APEX version
is shown in Fig.~\ref{fig:APEX_parallel}.
In \cite{Chen91} it was shown that the heuristic
analysis above is indeed true. Formally the parallel version of APEX
is described below.

\begin{algorithm}[APEX, Parallel Version]
For all neurons $m=1$ to $N$, ($N\leq n$) in parallel do:
\begin{enumerate}
\item Initialize the weights ${\bf w}_{m}$, ${\bf c}_m$
for each neuron to some random values.

\item Set $\beta_m$ from Eqs.~\ref{eq:APEX_beta1}
or \ref{eq:APEX_beta2}. Notice that each node has
its own customized step-size parameter.

\item Update ${\bf w}_{mk}$ and ${\bf c}_{mk}$ according
to Eqs.~\ref{eq:APEX_dw}, \ref{eq:APEX_dc}, until a stopping criterion is met.
\end{enumerate}
\end{algorithm}

\begin{figure}
\centerline{\psfig{figure=parAPEX.ps,height=3in}}
\caption{The parallel APEX network model.}
\label{fig:APEX_parallel}
\end{figure}


\subsection{Neural Net models vs. batch methods}

Comparison of the APEX model with batch methods
for eigenvalue decomposition is
not straightforward because they approach the problem
using different assumptions and they are practical for
different applications. Because of storage considerations
batch methods are probably more preferable when
relatively few data are to be processed relatively few times.
Adaptive methods on the other hand,
are preferred when arbitrarily long data
are to be processed, while they 
can provide at any time a reasonable approximation of the
optimal solution.
In addition to requiring less memory for data
storage, adaptive methods can track gradual changes in the optimal
solution rather inexpensively compared to batch methods.
On the negative side, adaptive methods tend to consume a
lot of iterations in order to produce as accurate results
as produced by batch methods. As with APEX, all adaptive schemes
converge to the final solution arbitrarily close only if
$\beta_k \rightarrow 0$ and unfortunately, the optimal
schedule for decreasing $\beta$ is not known.

As an example the Housholder-QR, probably the fastest batch
method for eigenvalue decomposition of non-sparse matrices,
requires $O(n^3)$ floating point operations ({\em flops})
\cite[chapter 8]{Golub89}
where $n$ is the size of the matrix which in this case is
the dimension of ${\bf x}_k$. In addition to that complexity,
the estimation of ${\bf R}_x$ as an average of the outer products
${\bf x}_k{\bf x}_k^T$, would typically require $O(Mn^2)$ flops
where $M$ is the number of patterns ${\bf x}_k$ to be processed.
APEX on the other hand has complexity $O(n)$ per iteration
for the computation of the $m$-th component.
In this we assumed that the activations of the $m\!-\!1$
previous neurons corresponding to each pattern
are stored in memory so they do not need to
be recomputed every time.
In many applications we don't need all the components there are.
For example, in harmonic retrieval problems the signal lies
in a low dimensional space the dimension of which is usually quite
smaller than $n$ therefore the components required are much
less than $n$. Therefore for some applications, $m$ can also
be considered constant and can be dropped from the
argument of $O(\cdot)$. However, 
the number of sweeps $S$, required
for each neuron to extract the corresponding component within
reasonable accuracy is not easily identifiable. Therefore
even for fixed $m$ the
comparison between the total APEX complexity $O(SMn)$,
and the $O(n^3+Mn^2)$ complexity of QR is
difficult. However, simulation evidence supports the proposition
that QR is faster for reasonable size data and
for the same accuracy requirement. On the other hand
using the QR method for computing the eigenvalues every time
a new pattern arrives would incur $O(n^3)$ flops per pattern
or a total of $O(Mn^3)$ flops. This puts the QR method
in disadvantage compared to APEX. The disadvantage of
batch methods over APEX
in adaptive applications has also been observed in
\cite{Martinelli92}.


\subsection{APEX vs. other NN models}

The APEX algorithm also compares
favorably with all the previous NN methods for PCA.
The claim can be supported in the following aspects:

\begin{enumerate}
\item Efficiency in Recursively Computing New PCs.

Every time an order update of the model is required
one can {\em grow} the APEX network by one more
neuron or {\em prune} it by one or more units
without need for retraining the old nodes.
Referring to Eq.~\ref{eq:Sanger},
growing the network by one unit would require
more operations per iteration for the Sanger model
compared with APEX.
This is because Sanger computes the outer product
$\sum_{j\leq m} y_{jk} {\bf w}_{mk}$
that costs $O(mn)$ multiplications per iteration
even if the values $y_{jk}$ are stored in memory,
while APEX utilizes the lateral connection network
to effectively compute something similar to this product
with only $O(n)$ operations. The computational advantage
of APEX is clear over Sanger's model.

\item Convergence Speed-up by Adopting
Optimal Step-size Parameter.

Our study regarding the best step-size parameter
$\beta$, results in impressive convergence speeds.
Such numerical analysis is lacking from all other models.

\item Extensibility to the CPC problem.

The APEX model is even more general in that it can be used
for extracting the constrained principal components of a random process
as will be discussed in Section \ref{sec:CPC}.
\end{enumerate}

\subsection{Simulation Results}

In this section we present two typical sets of simulations
that demonstrate the numerical performance of the APEX
algorithm. The two sets are one for the sequential and
one for the parallel model. The data used are artificially
created colored random vectors of dimension $n=64$. The same
data are used for both the sequential and the parallel model
for purposes of comparison. The autocorrelation matrix
${\bf R}_x$ has eigenvalue distribution shown in Fig.~\ref{fig:APEX_sim_eigs}.
The eigenvalues are not equal so ${\bf R}_x \neq {\bf I}$,
and the sequence ${\bf x}_k$ is not white.
In order to ensure stationarity we have repeated the sequence
periodically, so that the convergence performance of the algorithm
can be studied. Each period of the data is called a sweep
and contains 200 data samples.

\begin{figure}
\centerline{\psfig{figure=apex_sim_eigs.ps,width=4in}}
\caption{The eigenvalue distribution of \protect{${\bf R}_x$}.}
\label{fig:APEX_sim_eigs}
\end{figure}

In Fig.~\ref{fig:APEX_seq_sim} we show the convergence of
the first 4 components using the sequential APEX model.
We used Eq.~\ref{eq:APEX_beta1}
with $\gamma=0.995$ for determining the value of $\beta$\footnote{
Obviously we could have also used Eq.~\ref{eq:APEX_beta2} with $M=200$
achieving similar results.}.
The convergence is exponential
as predicted by our theoretical analysis in Section \ref{sec:APEX-theory}.
Notice that we have superimposed the plots for each neuron
in order to save space:
in fact the second unit starts {\em after} the first one has
converged so the actual time reference for the start of
curve $b$ is $k=5000$ instead of $k=0$. Similarly
the third unit starts {\em after} the second one has converged, etc.
Since the convergence speed for each neuron depends predominantly
upon the difference $\lambda_m-\lambda_{m-1}$ (see the proof 
of Theorem \ref{th:APEX}) and not upon the
order of the unit, it comes as no surprise that the first
neuron has slower convergence compared to the other units.
Only in the parallel APEX model,
the speeds of convergence for the different
neurons are necessarily sorted according to their order.
This is indeed shown in the next set of simulations
which concern the parallel APEX model.

\begin{figure}
\begin{center}
\begin{tabular}{c}
\psfig{figure=apex_seq_sim_a.ps,height=3in}
\\
(a)
\\
\psfig{figure=apex_seq_sim_b.ps,height=3in}
\\
(b)
\end{tabular}
\end{center}
\caption{(a) Plot of the square distance
\protect{$\|{\bf e}_m-{\bf w}_{mk}\|^2$},
between the actual
components and the ones estimated using sequential APEX.
The data are repeated in sweeps (1 sweep =
200 iterations). (b) Average $y_{mk}^2$ over each sweep
for each neuron $m$.}
\label{fig:APEX_seq_sim}
\end{figure}

Fig.~\ref{fig:APEX_par_sim} shows the results of the simulation
on the same data used in Fig.~\ref{fig:APEX_seq_sim}
with the same choice of $\beta$ and $\gamma$.
Notice that in both figures related to the parallel model the
final variance of the components is larger than the variance of
the corresponding figures in the sequential case. As mentioned
before, this is due to the fact that each neuron in the
parallel model receives a noisy estimate of the components
prior to it so this noise adds up to inherent variance that the
algorithm produces. Notice also that in the parallel model,
the first unit converges faster than the second which in turn
converges faster than the third, etc, as expected from the
fact that each unit needs the convergence of all previous
ones prior to its own convergence.

\begin{figure}
\begin{center}
\begin{tabular}{c}
\psfig{figure=apex_par_sim_a.ps,height=3in}
\\
(a)
\\
\psfig{figure=apex_par_sim_b.ps,height=3in}
\\
(b)
\end{tabular}
\end{center}
\caption{(a) Plot of the quantity \protect{$\|{\bf e}_m-{\bf w}_{mk}\|^2$}.
Same data used in Fig.~\protect{\ref{fig:APEX_seq_sim}}.
(b) Average \protect{$y_{mk}^2$} over each sweep.}
\label{fig:APEX_par_sim}
\end{figure}


\section{Constrained Principal Component Analysis (CPC)}
\label{sec:CPC}

In the previous section we proposed the APEX model
for the recursive calculation of the principal components
of a stochastic process.
The new approach offers considerable computational
advantages over the previous approaches. It turns out also
that APEX introduces a new application domain (which can
not be handled by any of the previous known algorithms)
in its ability to extract {\em Constrained Principal Components (CPCs)}.
The problem arises in
applications, where certain subspaces are less
preferred than others, thus affecting the selection of ``best''
components. In CPC, we assume that there is
an undesirable subspace ${\cal L}$, spanned by the
columns of an orthonormal {\em constraint matrix}  ${\bf V}$.
${\cal L}$ may stand for the interference subspace in
interference-cancellation applications or
for a redundant subspace in applications of extracting novel components.

More formally, the CPC problem can be defined as follows:
given an $n$-dimensional stationary stochastic
input vector process $\{{\bf x}_k\}$,
and an $\ell$-dimensional ($\ell < n$)
constraint process, $\{{\bf v}_k\}$, such that
\[
\begin{array}{c} {\bf v}_k = {\bf V}{\bf x}_k, \\
{\bf V}: \mbox{ orthonormal, spanning } {\cal L} \end{array}
\]
find the most representative $m$-dimensional ($\ell +m\leq n$)
subspace ${\cal L}_m$ in the principal component sense
constrained to be orthogonal to ${\cal L}$.
In other words,
as in Section \ref{sec:PCA-formulation}, we are looking for the
optimal linear transformation
\[
{\bf y}_k = {\bf W} {\bf x}_k
\]
where ${\bf W}$ is orthonormal, such that the error
\[
J=E\|{\bf x}-\hat{\bf x}\|^2=E\|{\bf x}-{\bf W}^T {\bf y}\|^2
\]
is minimal under the constraint
\begin{equation}
{\bf W} {\bf V}^T = 0
\label{eq:CPC_orthogonality}
\end{equation}

${\bf y}_k$ is basically the projection of ${\bf x}_k$
on the subspace spanned by the columns of ${\bf W}$
which is orthogonal to ${\cal L}$.
In \cite{Kung90b} it was shown that the optimal solution to the CPC
problem is
\begin{equation}
{\bf W}^{*}=[{\bf e}_1 \ldots {\bf e}_m]^T
\end{equation}
where ${\bf e}_1,\ldots,{\bf e}_m$ are the principal
right eigenvectors of
the input {\em skewed autocorrelation matrix}:
\[
{\bf R}_s=({\bf I}-{\bf V}{\bf V}^T) {\bf R}_x.
\]
while the minimum error is
\[
J_{min} = \sum_{i=m+1}^n {\tilde \lambda_i}
\]
where ${\tilde \lambda_1}, {\tilde \lambda_2}, \ldots, {\tilde \lambda_n}$
are the eigenvalues of ${\bf R}_s$ in decreasing order. 
Much like PC analysis, the components now maximize the output variance,
but under the additional constraint that they are orthogonal to ${\bf V}$.

\subsection{APEX solves CPC}

The following facts can be easily shown
regarding the eigenvalues/vectors of ${\bg R}_s$:
\begin{enumerate}
\item All eigenvalues are real and nonnegative.
To see that notice that the matrix ${\bf I}-{\bf V}^T{\bf V}$
is idempotent, i.e. $({\bf I}-{\bf V}^T{\bf V})^k = 
({\bf I} -{\bf V}^T{\bf V})$, for all $k\geq 1$, and let ${\bf e}$
be a right eigenvector associated with eigenvalue $\lambda$.
Assume for the moment that ${\bf e}$ is complex; then
\[
{\bf e}^{*}{\bf R}_x({\bf I}-{\bf V}^T{\bf V}){\bf R}_x{\bf e} = 
\lambda {\bf e}^{*}{\bf R}_x{\bf e} \Rightarrow
\]
\[
\|({\bf I}-{\bf V}^T{\bf V}){\bf R}_x{\bf e}\|^2 = \lambda \|{\bf R}_x^{1/2}
{\bf e}\|^2
\]
where $^{*}$ denotes conjugate transpose.
Therefore $\lambda$ is not only real but also $\lambda\geq 0$,
and ${\bf e}$ is also real.

\item There are exactly $\ell$ zero eigenvalues of ${\bf R}_s$ and the
rest $n-\ell$ eigenvalues are positive. This is true because the
matrix $({\bf I}-{\bf V}^T{\bf V})$ has a null space of dimension
${\ell}$ and ${\bf R}_x$ was assumed invertible.

\item Eigenvectors corresponding to different nonzero eigenvalues are
orthogonal. Indeed, let
\begin{eqnarray*}
{\bf R}_s {\bf e}_i&=&\lambda_i {\bf e}_i,
\\
{\bf R}_s {\bf e}_j&=&\lambda_j {\bf e}_j,
~~~\lambda_i,\lambda_j >0, \lambda_i\neq \lambda_j.
\end{eqnarray*}
Then
\[
{\bf e}_i^T {\bf R}_x({\bf I}-{\bf V}^T{\bf V}) {\bf R}_x {\bf e}_j
= \lambda_i {\bf e}_i^T {\bf R}_x {\bf e}_j
= \lambda_j {\bf e}_i^T {\bf R}_x {\bf e}_j \Rightarrow
\]
\[
{\bf e}_i^T {\bf R}_x({\bf I}-{\bf V}^T{\bf V}) {\bf R}_x {\bf e}_j=0
\]
But
\[
{\bf e}_i^T {\bf R}_x({\bf I}-{\bf V}^T{\bf V}) {\bf R}_x {\bf e}_j
= \lambda_i \lambda_j {\bf e}_i^T {\bf e}_j = 0 \Rightarrow
\]
\[
{\bf e}_i^T {\bf e}_j=0
\]
\end{enumerate}

We'll show now that the APEX model tackles the CPC problem
defined above. Consider the APEX network in Fig.~\ref{fig:APEX_CPC}(a)
where the first $\ell$ units correspond to the undesirable
components ${\bf v}_1,\ldots,{\bf v}_{\ell}$. We shall call them
{\em constraint units}. We train the first
nonconstraint neuron with the APEX rule
\begin{eqnarray}
{\bf w}_{k+1}&=&{\bf w}_k+\beta (y_k{\bf x}_k-y_k^2{\bf w}_k)
\\
{\bf c}_{k+1}&=&{\bf c}_k+\beta (y_k{\bf v}_k-y_k^2{\bf c}_k)
\end{eqnarray}
where
\[
y_k={\bf w}_k^T{\bf x}_k-{\bf c}_k^T{\bf v}_k=
({\bf w}_k^T-{\bf c}_k^T{\bf V}){\bf x}_k
\]
We can prove the following

\begin{theorem}
Consider the algorithm defined above. Then with probability 1,
the quantity ${\bf q}_k={\bf w}_k-{\bf V}^T{\bf c}_k$
converges to $\pm {\bf e}_1$ as $k\rightarrow \infty$.
\label{th:CPC}
\end{theorem}

{\bf Proof.}
Following the proof of Theorem \ref{th:APEX} we see that
Eq.~\ref{eq:APEX_ode_q} still holds (with ${\bf V}$ substituted
for ${\bf W}$). We repeat this equation here for convenience
\begin{equation}
\frac{d {\bf q}'(t)}{dt} = ({\bf I} - {\bf V}^T{\bf V}) {\bf R}_x {\bf q}'(t)
- \sigma'(t) {\bf q}'(t)
\label{eq:CPC_ode_q}
\end{equation}
where now
\[
{\bf q}'(t) = {\bf w}'(t) - {\bf V}^T {\bf c}'(t)
\]
\[
\sigma'(t)={\bf q}'(t)^T{\bf R}_x{\bf q}'(t).
\]
Expanding ${\bf q}'(t)$ into the basis ${\bf e}_1,\ldots,{\bf e}_n$
\[
{\bf q}'(t)=\sum_{i=1}^n \theta_i(t) {\bf e}_i
\]
we transform Eq.~\ref{eq:CPC_ode_q} into
\begin{equation}
\frac{d\theta_i(t)}{dt}=[\tilde{\lambda}_i-\sigma'(t)] \theta_i(t)
\end{equation}
Lemma \ref{lemma:APEX} is applicable to Eq.~\ref{eq:CPC_ode_q}
so following exactly the reasoning in the proof of
Theorem \ref{th:APEX} we can show that
$\theta_i(t)\rightarrow 0$ as $t\rightarrow \infty$
for all $i>1$,
and $\theta_1$ converges to some finite nonzero value $\alpha$.
Plugging ${\bf q}'(t)=\alpha {\bf e}_1$ into Eq.~\ref{eq:CPC_ode_q}
we obtain that $\alpha=\pm\sqrt{\lambda_1/({\bf e}_1^T{\bf R}_x{\bf e}_1)}$.
However, notice that ${\bf e}_1^T{\bf V}^T=0$ (by definition)
so ${\bf e}_1^T{\bf R}_x{\bf e}_1={\bf e}_1^T({\bf I}-{\bf V}^T{\bf V})
{\bf R}_x{\bf e}_1=\lambda_1$ implying $\alpha=\pm 1$.
Thus finally,
\[
{\bf q}_k \rightarrow
\pm {\bf e}_1, ~~\mbox{as } k \rightarrow \infty
\]
\hfill \qed

In the general case for CPC,
it is not true that ${\bf c}_k\rightarrow 0$.
In the special case where the rows of ${\bf V}$ are principal
eigenvectors of ${\bf R}_x$ then the problem is reduced to the
standard Principal Component Analysis and Theorem \ref{th:APEX}
applies.

\begin{figure}
\centerline{\psfig{figure=oln.ps,height=2.5in}}
\caption{(a) The single-output model. The connections denote the
weights \protect{$w_i$}, \protect{$c_j$} which are trained.
(b) The multiple-output
model.}
\label{fig:APEX_CPC}
\end{figure}

The CPC problem can be solved recursively be the APEX model.
Since the components are orthogonal, once the
first component has been extracted it can be appended to the rows
of ${\bf V}$ to obtain a new orthogonal constraint matrix
\[
{\bf V}' = \left[ \begin{array}{c} {\bf V}\\{\bf e}_1^T \end{array} \right]
\]
which will be used to obtain the next component.

\subsection{Simulation Results}

Figure \ref{fig:cpc_sim} depicts the convergence of the APEX
network for the CPC problem. The data are artificially
created and have dimension $n=64$. There $\ell=8$ constraint 
neurons while the orthonormal constraint matrix ${\bf V}$
is randomly picked. The convergence is exponential
as in the PCA case.

\begin{figure}
\begin{center}
\begin{tabular}{c}
\psfig{figure=cpc_seq_sim_a.ps,height=3in}
\\
(a)
\\
\psfig{figure=cpc_seq_sim_b.ps,height=3in}
\\
(b)
\end{tabular}
\end{center}
\caption{APEX model performance for CPC. We use
200 64-dimensional data vectors repeated in sweeps.
(a) The convergence of the error \protect{$\|{\bf q}_k-{\bf e}_m\|^2$}
for each neuron \protect{$m$}. (b) The output variance:
\protect{$average\{y^2\}$}
over each sweep for each neuron. The final variance equals approximately
to the corresponding eigenvalue \protect{$\tilde{\lambda}_m$}.}
\label{fig:cpc_sim}
\end{figure}
