\chapter{Asymmetric Principal Component Analysis}
\label{ch:APCA}


\section{Introduction}

As we saw in the previous chapters
Principal Component Analysis is a statistical method
for finding those subspace directions that contain
the maximum signal information, and as such
it is useful for many applications in 
signal/image processing, data compression,
pattern classification, recognition, etc
\cite{Deprettere88}\cite{Vaccaro91}.
In this chapter we expand this method
by distinguishing two basic categories
of PCA problems: Symmetric and Asymmetric
PCA (APCA). Different Neural Network models can be used to solve
these problems as described below.

\begin{enumerate}
\item {\bf Symmetric (ordinary) PCA and Self-Supervised networks.}

{\em Symmetric} --i.e. ordinary-- PCA is concerned
with the Singular Value Decomposition (SVD) of the autocorrelation
matrix ${\bf R}_x$ of the input $\{{\bf x}_k\}$.
It has been shown that a {\em Self-Supervised} two-layer
network with linear units trained using back propagation, produces
the optimal subspace spanned by the eigen-components \cite{Bourlard88}.
The APEX model (introduced in Chapter \ref{ch:PCA})
and other neural models based on the Hebbian rule
have been proposed for extracting
the components themselves rather than their subspaces
\cite{Oja82}\cite{Sanger89b}\cite{Kung90a}.


\item {\bf Asymmetric PCA and Hetero-Supervised networks.}

The {\em Asymmetric PCA} (APCA) problem
involves two vector stochastic processes:
the input process $\{{\bf x}_k\}$, and the teacher process $\{{\bf y}_k\}$,
with autocorrelation matrices ${\bf R}_x$ and ${\bf R}_y$ respectively
and cross-correlation matrix ${\bf R}_{xy}=E\{{\bf x}_k {\bf y}_k^T\}$.
Two problems are defined to belong to this class of PCA:

$\bullet${\bf Linear Approximation APCA:}

Find the most representative components of $\{{\bf x}_k\}$
so that they may be used to best recover $\{{\bf y}_k\}$.
As shown in Section \ref{sec:linear_approximation_APCA}
if we use the least-squares-error
as the optimality criterion, 
then  the APCA problem can be reduced  to
the problem of finding the principal left (right) singular vectors
of the matrix ${\bf R}_x^{-1/2} {\bf R}_{xy}$, if ${\bf R}_x$
is invertible, or the Generalized SVD (GSVD) of the
matrix pair $({\bf X}^T,{\bf Y}{\bf X}^T)$ where ${\bf X}$
and ${\bf Y}$ are the input and output data matrices
as explained below. The problem is also
equivalent to reduced rank Wiener filtering \cite{Scharf91}
or reduced rank linear regression.

$\bullet${\bf Cross-correlation APCA:}

Find the left and right singular vectors for ${\bf R}_{xy}$, corresponding
to the largest $p$ singular values, namely produce the
optimal $p$-rank approximation of the cross-correlation matrix.

\end{enumerate}
It's worth mentioning that if ${\bf y}_k = {\bf x}_k$
then the ordinary PCA and the two APCA problems
discussed above become the same.



\section{Linear Approximation APCA}
\label{sec:linear_approximation_APCA}

In this chapter (except when explicitly noted) we are dealing
with linear two-layer networks.
Although a two-layer linear network is always equivalent
to some single-layer linear network in terms of the overall mapping, 
the presence of the hidden units guards the rank of the mapping.
Thus, if the network has $n$ inputs, $m$ outputs and $p$
hidden units and if $\ul{\bf W} \in {\cal R}^{n \times p}$,
$\ol{\bf W} \in {\cal R}^{m \times p}$ denote the
lower and upper layer weight matrices respectively,
then the overall mapping ${\bf W} = \ol{\bf W} \ul{\bf W}^T$
is such that $rank({\bf W}) \leq p$. A reduced hidden
layer network is defined as one with $p \leq \min\{m,n\}$.
It was first observed that in the auto-associative case
the optimal least-squares weights are related
to the {\em Singular Value Decomposition (SVD)}
of the input data matrix \cite{Bourlard88}.
Any learning rule that can find the optimal solution
for the sum-of-squares error criterion --like for example,
the linear version of the Back-Propagation
algorithm \cite{Rumelhart86}-- will produce weights
that span the same space as the $p$ principal components
of the input sequence. Later it was found that in the
hetero-associative case the result is related to the
eigenvalue decomposition of the classical least-squares
regression matrix \cite{Baldi89} or equivalently, the SVD
of the cross-correlation between the output-data and the
pre-whitened input data \cite{Kung91,Scharf91}.
In all these publications however, it was assumed that the
input data matrix has full rank and equivalently the
input autocorrelation matrix is non-singular.

In this Section
we'll show that the invertibility condition for the input
autocorrelation matrix can be relaxed and the solution
obtained by the network when minimizing the least-squares
error criterion is related to the {\em Generalized Singular
Value Decomposition (GSVD)} of the input-teacher
cross-correlation matrix and the input data matrix.
Then all previous results can be treated as special cases \cite{Diamantaras92}.

\subsection{Solution to the Reduced-Rank Approximation Problem}

The reduced-rank linear approximation problem is formulated as follows:
Given $N$ input/teacher pairs $({\bf x}_k , {\bf y}_k)$
we want to minimize the cost
\begin{equation}
J = \sum_{k=1}^N \| {\bf y}_k - {\bf W} {\bf x}_k \|^2
= \| {\bf Y} - {\bf W} {\bf X} \|_F^2
\label{eq:lin_approx_error_function}
\end{equation}

where $\|\cdot\|_F^{~}$ stands for the Frobenious norm,
${\bf Y} \in {\cal R}^{m \times N}$,
${\bf X} \in {\cal R}^{n \times N}$,
${\bf W} \in {\cal R}^{m \times n}$,
and $rank({\bf W}) = p \leq q \equiv \min\{m,n\}$.
Assume that\footnote{This condition is necessary for the existence
of the GSV Decomposition to be used later. It is not restricting
however since we can always pad as many 0-columns in ${\bf X}$
and ${\bf Y}$ as needed without changing the error function.}
$N \geq n$ and let
\begin{eqnarray}
{\bf U}^T {\bf X}^T {\bf Q} &=&
{\bf D}_1 = diag\{\alpha_1, \ldots, \alpha_n\}
\label{eq:gsvd_x}
\\
{\bf V}^T {\bf Y} {\bf X}^T {\bf Q} &=&
{\bf D}_2 = diag\{\beta_1, \ldots, \beta_q\}
\label{eq:gsvd_yx}
\end{eqnarray}

be the GSVD of the matrix pair $({\bf X}^T,{\bf Y}{\bf X}^T)$
\cite{VanLoan76}. Thus ${\bf U} \in {\cal R}^{N \times N}$
and ${\bf V} \in {\cal R}^{m \times m}$ are orthonormal matrices,
${\bf Q} \in {\cal R}^{n \times n}$ is invertible,
${\bf D}_1 \in {\cal R}^{N \times n}$,
${\bf D}_2 \in {\cal R}^{m \times n}$ are diagonal
(but not square in general),
$\alpha_1, \ldots, \alpha_r > 0$, $\alpha_{r+1} = \ldots =
\alpha_n = 0$ for some $r$, $1 \leq r \leq n$,
and $\beta_1, \ldots , \beta_s > 0$, $\beta_{s+1} =
\ldots = \beta_q \geq 0$ for some $s$, $1 \leq s \leq q$.
Obviously, $rank({\bf X}^T) = r$, $rank({\bf Y}{\bf X}^T) = s$
and $s \leq r$.
Let's define $\lambda_i \equiv \displaystyle
\frac{\beta_i}{\alpha_i}$, for
$i=1, \ldots , \hat{r} \equiv \min\{q,r\} = \min\{m,r\}$.
The numbers $\lambda_i^2$ are
referred to as the finite {\em generalized eigenvalues}
of the matrix pair $({\bf X}{\bf Y}^T{\bf Y}{\bf X}^T , {\bf X}{\bf X}^T)$.
We will assume, without loss of generality, that the diagonal
entries of ${\bf D}_1$ and ${\bf D}_2$ are ordered so that
$\lambda_1 \geq \lambda_2 \geq \ldots \geq \lambda_s >
\lambda_{s+1} = \ldots = \lambda_{\hat{r}} = 0$.
Given the above observations we can write
\begin{eqnarray}
{\bf D}_1 &=& \left[ \begin{array}{cc} {\bf D}_{1r} & 0 \\ 0 & 0 \end{array}
\right]
\\
{\bf D}_2 &=& [ {\bf D}_{2r} ~ \tilde{\bf D}_2 ]
\\
{\bf U} &=& [ {\bf U}_r ~ \tilde{\bf U} ]
\end{eqnarray}

where ${\bf D}_{1r}$ is a diagonal non-singular $r \times r$ matrix,
and each of the matrices ${\bf D}_{2r}$,
${\bf U}_r$ have $r$ columns (and of course appropriate number of rows).
Then
\begin{equation}
{\bf X}^T {\bf Q} = {\bf U}_r [ {\bf D}_{1r} ~ 0 ].
\end{equation}

\begin{theorem}
\label{th:APCA-optimal-solution}
Let ${\bf v}_1, \ldots, {\bf v}_m$
and ${\bf q}_1, \ldots, {\bf q}_n$ denote the
columns of ${\bf V}$ and ${\bf Q}$ respectively.
Then $J$ in Eq. \ref{eq:lin_approx_error_function} is minimized for
\begin{equation}
{\bf W}^{*} = {\bf V}_p {\bf \Lambda}_p {\bf D}_{1p}^{-1} {\bf Q}_p^T
+ {\bf \Phi} {\bf Q}_1^T
\label{eq:optimal_W}
\end{equation}
where ${\bf V}_p \equiv [{\bf v}_1 \ldots {\bf v}_p]
\in {\cal R}^{m \times p}$,
${\bf Q}_p \equiv [{\bf q}_1 \ldots {\bf q}_p]
\in {\cal R}^{n \times p}$,
${\bf D}_{1p} \equiv diag\{\alpha_1, \ldots , \alpha_p\}
\in {\cal R}^{p \times p}$ and
${\bf \Lambda}_p = diag\{\lambda_1, \ldots , \lambda_p\}
\in {\cal R}^{p \times p}$,
${\bf Q}_1$ is a matrix formed by the last $n-r$ columns
of ${\bf Q}$, and ${\bf \Phi}$ is any $m \times (n - r)$ matrix.
Furthermore, the minimum error is
\begin{equation}
J_{min} = \left\| {\bf V}^T {\bf Y} \tilde{\bf U} \right\|_F^2
+ \sum_{i=p+1}^s \lambda_i^2
\label{eq:optimal_error}
\end{equation}
\end{theorem}

{\bf Proof.}
Using the fact that pre- or post-multiplication with an orthonormal matrix
does not affect the Frobenious norm we have
\begin{eqnarray}
J = \| {\bf Y} - {\bf W} {\bf X} \|_F^2 &=&
\left\| {\bf V}^T {\bf Y} {\bf U} - {\bf V}^T {\bf W} {\bf Q}^{-T} {\bf Q}^T
{\bf X} {\bf U} \right\|_F^2 = \nonumber
\\
&=& \left\| {\bf V}^T {\bf Y} [{\bf U}_r ~ \tilde{\bf U}] -
{\bf V}^T {\bf W} {\bf Q}^{-T}
\left[ \begin{array}{cc} {\bf D}_1^T & 0 \\ 0 & 0 \end{array}
\right] \right\|_F^2 = \nonumber
\\
&=& \left\| {\bf V}^T {\bf Y} {\bf U}_r - {\bf V}^T {\bf W} {\bf Q}^{-T}
\left[ \begin{array}{c} {\bf D}_{1r} \\ 0 \end{array} \right] \right\|_F^2 +
\left\| {\bf V}^T {\bf Y} \tilde{\bf U} \right\|_F^2
\end{eqnarray}
where we used the fact that ${\bf D}_{1r}^T = {\bf D}_{1r}$.
The residue
\[
R = \left\| {\bf V}^T {\bf Y} \tilde{\bf U} \right\|_F^2
\]
does not
depend on ${\bf W}$ thus it can be ignored. This also implies
that $J_{min} \geq R$, where $J_{min}$ is the minimum error.
Thus ${\bf W}$ minimizes
\begin{eqnarray}
\tilde{J} &=&
\left\| {\bf V}^T {\bf Y} {\bf U}_r - {\bf V}^T {\bf W} {\bf Q}^{-T}
\left[ \begin{array}{c} {\bf D}_{1r} \\ 0 \end{array} \right] \right\|_F^2
\nonumber
\\
&=& \left\| {\bf V}^T {\bf Y} {\bf U}_r \left[ {\bf D}_{1r} ~ 0 \right]
\left[ \begin{array}{c} {\bf D}_{1r}^{-1} \\ 0 \end{array} \right] -
{\bf V}^T {\bf W} {\bf Q}^{-T}
\left[ \begin{array}{c} {\bf D}_{1r} \\ 0 \end{array} \right] \right\|_F^2
\nonumber
\\
&=& \left\| {\bf V}^T {\bf Y} {\bf X}^T {\bf Q}
\left[ \begin{array}{c} {\bf D}_{1r}^{-1} \\ 0 \end{array} \right] -
{\bf V}^T {\bf W} {\bf Q}^{-T}
\left[ \begin{array}{c} {\bf D}_{1r} \\ 0 \end{array} \right] \right\|_F^2
\nonumber
\\
&=& \left\| {\bf D}_{2r} {\bf D}_{1r}^{-1} - {\bf V}^T {\bf W} {\bf Q}^{-T}
\left[ \begin{array}{c} {\bf D}_{1r} \\ 0 \end{array} \right] \right\|_F^2
\nonumber
\\
&=& \left\| {\bf \Lambda}_{\hat{r}} - {\bf V}^T {\bf W} {\bf Q}^{-T}
\left[ \begin{array}{c} {\bf D}_{1r} \\ 0 \end{array} \right] \right\|_F^2
\label{eq:APCA-theorem-J}
\end{eqnarray}
where ${\bf \Lambda}_{\hat{r}} = diag\{\lambda_1, \ldots, \lambda_{\hat{r}}\}
= diag\{\lambda_1, \ldots, \lambda_s, 0, \ldots, 0\}
\in {\cal R}^{m \times r}$. Assuming that $p \leq s$
%and $\lambda_1 > \ldots > \lambda_p > \lambda_{p+1} \geq \ldots \geq \lambda_s$
the optimal rank-$p$ matrix ${\bf W}^{*}$, should satisfy the equation
\begin{equation}
{\bf V}^T {\bf W}^{*} {\bf Q}^{-T}
\left[ \begin{array}{c} {\bf D}_{1r} \\ 0 \end{array} \right] =
\left[ \begin{array}{cc} {\bf \Lambda}_p & 0 \\ 0 & 0 \end{array} \right] 
\in {\cal R}^{m \times r}
\label{eq:optimal_W_condition}
\end{equation}

Then Eq. \ref{eq:optimal_W_condition} is solved for
\[
{\bf W}^{*} = {\bf V}_p {\bf \Lambda}_p {\bf D}_{1p}^{-1} {\bf Q}_p^T
+ {\bf \Phi} {\bf Q}_1^T
\]
as claimed by the theorem. Also the minimum error
is easily computed from Eq. \ref{eq:APCA-theorem-J} to be:
\[
J_{min} = R + \sum_{i=p+1}^s \lambda_i^2 ~.
\]
\hfill \qed

In the above proof
the value $R$, is the mismatch of the linear regression between ${\bf X}$
and ${\bf Y}$, and it is the absolute minimum error achievable when there
is no rank restriction on ${\bf W}$.
By introducing rank restrictions on ${\bf W}$
we create the excess error term $\displaystyle \sum_{i=p+1}^s \lambda_i^2$.
Scharf finds a similar result in \cite{Scharf91}
under the special case when ${\bf X} {\bf X}^T$ is non-singular.


\subsection{Comments}

\begin{itemize}
\item The optimal solution ${\bf W}^{*}$ in Eq.~\ref{eq:optimal_W}
is obviously not unique because of the involvement of the arbitrary
matrix ${\bf \Phi}$, but an additional source of non-uniqueness is
created when there is repetition of $\lambda_p$ or any eigenvalue
less than $\lambda_p$, since in this case,
the matrices ${\bf V}$ and ${\bf Q}$ are not unique as well.

\item Although reduced-rank approximation introduces excess error
when minimizing $J$, compared to full-rank approximation,
it might in fact perform better than the full-rank case with respect
to some other, more meaningful error criterion $J'$. For example,
if ${\bf y}_k = {\bf signal}_k + {\bf noise}_k$ it would make
sense, from a filtering application point of view, to minimize
the cost
\[
J' = \sum_k \| {\bf signal}_k - {\bf W} {\bf x}_k \|^2
\]
rather than the standard $J = \sum_k \| {\bf y}_k -
{\bf W} {\bf x}_k \|^2$. This would be the case if the ${\bf signal}$
has low rank $p$, but ${\bf y} = {\bf signal} + {\bf noise}$ has
higher rank $p' > p$. Then the full-rank filter
will be modeling the noise in at least
$p' - p$ dimensions. On the other hand, if we have an a priori
knowledge of $p$ then using a rank-$p$ approximation we get
a good approximation of the signal subspace and get rid of the
noise in all other dimensions. For low signal-to-noise ratios
the power of the noise in these dimensions may be significant
so we can get superior performance (with respect to $J'$)
over full-rank filtering. For high SNR the
performances of these two filters,
with respect to $J'$,
should approach each other asymptotically.
\end{itemize}

\subsection{Special Cases}

\begin{itemize}
\item {\bf Linear Regression: ${\bf X}$ has full rank, ${\bf W}$
has no rank restriction.}

If $r=n$ and in addition there is no
rank constraint on ${\bf W}$ then from Eq.~\ref{eq:optimal_W}
we see that an optimal solution is
\begin{eqnarray}
{\bf W}^{*} &=& {\bf V} {\bf \Lambda}_{\hat{r}} {\bf D}_{1n}^{-1} {\bf Q}^T
= \nonumber
\\
&=& {\bf V} {\bf D}_2 {\bf D}_{1n}^{-1} {\bf D}_{1n}^{-1} {\bf Q}^T
= \nonumber
\\
&=& {\bf V} {\bf D}_2 {\bf Q}^{-1} ({\bf X} {\bf X}^T)^{-1}
= \nonumber
\\
&=& {\bf V} {\bf V}^T {\bf Y} {\bf X}^T ({\bf X} {\bf X}^T)^{-1}
= \nonumber
\\
&=& {\bf Y} {\bf X}^T ({\bf X} {\bf X}^T)^{-1}
\label{eq:wiener}
\end{eqnarray}

which is the classical solution of the linear regression
problem. Assuming that
$\frac{1}{N}{\bf Y} {\bf X}^T \approx {\bf R}_{yx} = {\bf R}_{xy}^T$,
and $\frac{1}{N}{\bf X} {\bf X}^T \approx {\bf R}_x$,
one also recognizes Eq.~\ref{eq:wiener} as the
Wiener filter solution ${\bf R}_{yx}{\bf R}_x^{-1}$.

\item {\bf Reduced Rank Wiener Filtering: ${\bf X}$ has full rank,
${\bf W}$ has reduced rank.}

In this case $r=n$ so the matrix
${\bf X} {\bf X}^T$ is invertible. From Eqs. \ref{eq:gsvd_x}, \ref{eq:gsvd_yx}
we have
\begin{eqnarray}
{\bf Q}^T {\bf X} {\bf X}^T {\bf Q} &=& {\bf D}_1^T {\bf D}_1
\label{eq:QXXQ.eq.D1D1}
\\
{\bf Q}^T {\bf X} {\bf Y}^T {\bf Y} {\bf X}^T {\bf Q} &=&
{\bf D}_2^T {\bf D}_2
\label{eq:QXYYXQ.eq.D2D2}
\end{eqnarray}

and since $\alpha_i > 0$, $i=1, \ldots , n$ and $N \geq n$,
${\bf D}_1^T {\bf D}_1$ is non-singular therefore we can write
\begin{equation}
{\bf X} {\bf Y}^T {\bf Y} {\bf X}^T {\bf Q} =
{\bf X} {\bf X}^T {\bf Q} {\bf \Lambda}_n^2
\label{eq:ge}
\end{equation}

where ${\bf \Lambda}_n^2 = ({\bf D}_1^T {\bf D}_1)^{-1} {\bf D}_2^T {\bf D}_2 =
diag\{\lambda_1^2, \ldots, \lambda_s^2, 0, \ldots, 0\}
\in {\cal R}^{n \times n}$. Therefore, the columns of ${\bf Q}$ are the
generalized eigenvectors of the pair $({\bf X} {\bf Y}^T {\bf Y} {\bf X}^T ,
{\bf X} {\bf X}^T)$. One equivalent way of writing Eq.~\ref{eq:ge}
when ${\bf X}{\bf X}^T$ is invertible is
\begin{equation}
({\bf X}{\bf X}^T)^{-1/2} {\bf X} {\bf Y}^T {\bf Y} {\bf X}^T
({\bf X}{\bf X}^T)^{-1/2} \tilde{\bf Q} =
\tilde{\bf Q} {\bf \Lambda}_n^2
\end{equation}

where $\tilde{\bf Q} = ({\bf X}{\bf X}^T)^{1/2} {\bf Q}$.
Since the matrix $({\bf X}{\bf X}^T)^{-1/2}$ is the pre-whitening
operator for the input sequence,
$\tilde{\bf Q}$ is the eigenvector matrix of the cross-correlation
between the pre-whitened input data and the output data, in accordance
to the results published in \cite{Kung91,Scharf91}.

Yet another way to see the same solution is by using Eq.~\ref{eq:gsvd_yx}
to obtain
\begin{eqnarray}
{\bf Y} {\bf X}^T ({\bf X} {\bf X}^T)^{-1} {\bf X} {\bf Y}^T {\bf V} &=&
{\bf Y} {\bf X}^T ({\bf X} {\bf X}^T)^{-1} {\bf Q}^{-T} {\bf D}_2^T
\nonumber
\\
&=& {\bf Y} {\bf X}^T {\bf Q} ({\bf D}_1^T {\bf D}_1)^{-1} {\bf D}_2^T
\nonumber
\\
&=& {\bf V} {\bf D}_2 ({\bf D}_1^T {\bf D}_1)^{-1} {\bf D}_2^T
\nonumber
\\
&=& {\bf V} {\bf \Lambda}_m^2
\end{eqnarray}

where ${\bf \Lambda}_m = diag\{\lambda_1, \ldots, \lambda_s, 0, \ldots, 0\}
\in {\cal R}^{m \times m}$. Therefore, ${\bf V}$ is the eigenvector matrix
of ${\bf Y} {\bf X}^T ({\bf X} {\bf X}^T)^{-1} {\bf X} {\bf Y}^T$
in accordance to the results published in \cite{Baldi89}.

A possible application example could be an echo-cancellation problem
where the interference $\nu_k$ that corrupts the signal
$s_k$, has inherent low-rankness and it is correlated with some
observed signal $\nu_k'$.
In that case there is reason to
prefer a low-rank approximation filter rather
than the standard Wiener filter for estimating $\nu$ (Figure \ref{fig:echo}).
Another application is system identification when the input data
have low rank. We will further elaborate on this problem in
Chapter \ref{ch:Applications}.

\begin{figure}
\centerline{\psfig{figure=echo.sl.ps,width=4in}}
\caption{When the interference has low rank properties then
a reduced rank filter is preferred.}
\label{fig:echo}
\end{figure}

\item {\bf Auto-association: ${\bf Y} = {\bf X}$.}

In this case $m = q = n$, $s = r$ and Eq.~\ref{eq:gsvd_yx} yields
\[
{\bf V}^T {\bf X} {\bf U} {\bf D}_1 = {\bf D}_2 \Rightarrow
{\bf V}^T {\bf X} [{\bf U}_r ~ \tilde{\bf U}]
\left[ \begin{array}{cc} {\bf D}_{1r} & 0 \\ 0 & 0 \end{array} \right] =
\left[ \begin{array}{cc} {\bf D}_{2r} & 0 \\ 0 & 0 \end{array} \right]
\]

Since ${\bf X} \tilde{\bf U} = 0$ we can write
\[
{\bf V}^T {\bf X} [{\bf U}_r ~ \tilde{\bf U}]
\left[ \begin{array}{cc} {\bf D}_{1r} & 0 \\ 0 & {\bf I} \end{array} \right] =
\left[ \begin{array}{cc} {\bf D}_{2r} & 0 \\ 0 & 0 \end{array} \right] \Rightarrow
{\bf V}^T {\bf X} {\bf U} =
\left[ \begin{array}{cc} {\bf D}_{2r} {\bf D}_{1r}^{-1} & 0 \\ 0 & 0 \end{array}
\right]
\]

Therefore, ${\bf U}$ and ${\bf V}$ are the SVD matrices for ${\bf X}^T$.
Furthermore, we can choose ${\bf Q} = {\bf V}$ so that
Eq.~\ref{eq:gsvd_x} also becomes the SVD of ${\bf X}^T$, and Eq.~\ref{eq:gsvd_yx}
becomes the Eigenvalue Decomposition of ${\bf X} {\bf X}^T$.
This is in agreement with the results published in \cite{Bourlard88}.

\end{itemize}


\subsection{Neural Networks for Linear Approximation APCA}

The error $J$, in Eq.~\ref{eq:lin_approx_error_function},
applied to a linear network with $p$ hidden units
can be written in the following form
\begin{eqnarray}
J_{net}(\ol{\bf W},\ul{\bf W}) &=& \| {\bf Y} - \ol{\bf W} \ul{\bf W}^T {\bf X} \|_F^2 =
\nonumber \\
&=& \left\| {\bf \Lambda}_{\hat{r}} - {\bf V}^T \ol{\bf W} \ul{\bf W}^T
{\bf Q}^{-T} \left[ \begin{array}{c} {\bf D}_{1r} \\ 0 \end{array}
\right] \right\|_F^2 + R
\label{eq:net_J}
\end{eqnarray}

where $\ol{\bf W} \in {\cal R}^{m \times p}$, $\ul{\bf W} \in
{\cal R}^{n \times p}$ are the upper-, lower-layer weights respectively.
Then we can show the following result relevant
to any training algorithm minimizing \ref{eq:net_J}:

\begin{theorem}
\label{th:APCA-net-local-minima}
Under the assumption \label{page:condition_A}
$$
\lambda_1 > \lambda_2 > \ldots > \lambda_{\hat{r}} > 0
\mbox{~~(so $s = \hat{r}$)} \eqno(A)
$$
the error function \ref{eq:net_J} contains no local minima
(all critical points are saddle points).
\end{theorem}

{\bf Proof.}
Any critical point $(\ol{\bf W}, \ul{\bf W})$ of
Eq.~\ref{eq:net_J} is defined as one where
$\partial J_{net} / \partial \ol{\bf W} = 0$ and
$\partial J_{net} / \partial \ul{\bf W}^T = 0$,
with the derivative of a scalar function $f({\bf P})$,
over a matrix ${\bf P}$
defined as the matrix with elements $\partial f / \partial p_{ij}$.
It is easy to show that $\partial \| {\bf G} - {\bf K} {\bf P} {\bf L} \|_F^2
/ \partial {\bf P} = {\bf K}^T {\bf G} {\bf L}^T - {\bf K}^T {\bf K}
{\bf P} {\bf L} {\bf L}^T$, for any conformal set of matrices
${\bf G}, {\bf K}$ and ${\bf L}$, therefore any critical point
of $J_{net}$ satisfies the equations
\begin{eqnarray*}
{\bf V} {\bf \Lambda}_{\hat{r}} [{\bf D}_{1r} ~ 0] {\bf Q}^{-1} \ul{\bf W}
&=& \ol{\bf W} \ul{\bf W}^T {\bf Q}^{-T} \left[ \begin{array}{c}
{\bf D}_{1r} \\ 0 \end{array} \right] [{\bf D}_{1r} ~ 0]
{\bf Q}^{-1} \ul{\bf W}
\\
{\ol{\bf W}}^T {\bf V} {\bf \Lambda}_{\hat{r}} [{\bf D}_{1r} ~ 0]
{\bf Q}^{-1} &=& {\ol{\bf W}}^T \ol{\bf W} \ul{\bf W}^T {\bf Q}^{-T}
\left[ \begin{array}{c} {\bf D}_{1r} \\ 0 \end{array} \right]
[{\bf D}_{1r} ~ 0] {\bf Q}^{-1}
\end{eqnarray*}

Equivalently,
\begin{eqnarray*}
{\bf \Lambda}_{\hat{r}} [{\bf D}_{1r} ~ 0] {\bf Q}^{-1} \ul{\bf W}
&=& {\bf V}^T \ol{\bf W} \ul{\bf W}^T {\bf Q}^{-T} \left[ \begin{array}{c}
{\bf D}_{1r} \\ 0 \end{array} \right] [{\bf D}_{1r} ~ 0]
{\bf Q}^{-1} \ul{\bf W}
\\
{\ol{\bf W}}^T {\bf V} {\bf \Lambda}_{\hat{r}}
&=& {\ol{\bf W}}^T \ol{\bf W} \ul{\bf W}^T {\bf Q}^{-T}
\left[ \begin{array}{c} {\bf D}_{1r} \\ 0 \end{array} \right]
\end{eqnarray*}

or
\begin{eqnarray}
{\bf \Lambda}_{\hat{r}} \ul{\bf F} &=& \ol{\bf F} \ul{\bf F}^T \ul{\bf F}
\label{eq:Lambda_ulF}
\\
\ol{\bf F}^T {\bf \Lambda}_{\hat{r}} &=& \ol{\bf F}^T \ol{\bf F} \ul{\bf F}^T
\label{eq:olF_Lambda}
\end{eqnarray}

where $\ol{\bf F} = [ \ol{\bf f}_1 \ldots \ol{\bf f}_m ]^T
\equiv {\bf V}^T \ol{\bf W} \in {\cal R}^{m \times p}$ and
$\ul{\bf F} = [ \ul{\bf f}_1 \ldots \ul{\bf f}_r ]^T
\equiv [{\bf D}_{1r} ~ 0] {\bf Q}^{-1} \ul{\bf W} \in {\cal R}^{r \times p}$.
Manipulating the above expressions we obtain
\begin{eqnarray*}
Eq.~(\ref{eq:Lambda_ulF}) \Rightarrow {\bf \Lambda}_{\hat{r}}^T
({\bf \Lambda}_{\hat{r}}  \ul{\bf F}) \ol{\bf F}^T &=&
{\bf \Lambda}_{\hat{r}}^T \ol{\bf F} \ul{\bf F}^T \ul{\bf F} \ol{\bf F}^T
\\
&=& \ul{\bf F} \ol{\bf F}^T \ol{\bf F} \ul{\bf F}^T \ul{\bf F} \ol{\bf F}^T
\\
Eq.~(\ref{eq:olF_Lambda}) \Rightarrow \ul{\bf F} (\ol{\bf F}^T
{\bf \Lambda}_{\hat{r}}) {\bf \Lambda}_{\hat{r}}^T &=&
\ul{\bf F} \ol{\bf F}^T \ol{\bf F} \ul{\bf F}^T {\bf \Lambda}_{\hat{r}}^T
\\
&=& \ul{\bf F} \ol{\bf F}^T \ol{\bf F} \ul{\bf F}^T \ul{\bf F} \ol{\bf F}^T
\end{eqnarray*}

therefore
\[
{\bf \Lambda}_{\hat{r}}^T {\bf \Lambda}_{\hat{r}}
\ul{\bf F} \ol{\bf F}^T = \ul{\bf F} \ol{\bf F}^T
{\bf \Lambda}_{\hat{r}} {\bf \Lambda}_{\hat{r}}^T
\]

Using assumption $(A)$ we conclude that the matrix
${\bf \Lambda} \equiv \ol{\bf F} \ul{\bf F}^T \in {\cal R}^{m \times r}$ is diagonal, i.e.
\begin{equation}
{\bf \Lambda} = \ol{\bf F} \ul{\bf F}^T =
{\bf V}^T \ol{\bf W} \ul{\bf W}^T {\bf Q}^{-T}
\left[ \begin{array}{c} {\bf D}_{1r} \\ 0 \end{array} \right] =
diag\{ l_1 , \ldots , l_{\hat{r}} \}
\label{eq:L_eq_diag}
\end{equation}

In addition, from Eq.~\ref{eq:Lambda_ulF} (or Eq.~\ref{eq:olF_Lambda})
we obtain that either $l_i=0$ or $l_i = \lambda_i$ for all $i$.
Furthermore, from these equations and Eq.~\ref{eq:L_eq_diag} we get
$l_i = 0 \Leftrightarrow \ol{\bf f}_i = 0
\Leftrightarrow \ul{\bf f}_i = 0$. It follows that there is an
ordered set of indices ${\cal S} = \{i(1), \ldots , i(\hat{p}) \}
\neq \{ 1, \ldots , p \}$ so that $l_j \neq 0$ iff $j \in {\cal S}$.
Since the rank of ${\bf \Lambda}$ has to be at most $p$ we have
$\hat{p} \leq p$. From the above observations it follows that
there exist two matrices $\ol{\bf C} \in {\cal R}^{\hat{p} \times p}$,
$\ul{\bf C} \in {\cal R}^{\hat{p} \times p}$ so that
\begin{eqnarray*}
\ol{\bf F} &=& {\bf I}_{\cal S}^{m,\hat{p}} \ol{\bf C}
\\
\ul{\bf F} &=& {\bf I}_{\cal S}^{r,\hat{p}} \ul{\bf C}
\end{eqnarray*}

where ${\bf I}_{\cal S}^{m,\hat{p}} = \left[ {\bf e}_{i(1)}^{(m)} \ldots
{\bf e}_{i(\hat{p})}^{(m)} \right] \in {\cal R}^{m \times \hat{p}}$,
${\bf I}_{\cal S}^{r,\hat{p}} = \left[ {\bf e}_{i(1)}^{(r)} \ldots
{\bf e}_{i(\hat{p})}^{(r)} \right] \in {\cal R}^{r \times \hat{p}}$
and ${\bf e}_i^{(m)}$ (${\bf e}_i^{(r)}$)
is the $i$-th column of the $m \times m$ ($r \times r$) identity matrix.
Notice that ${\bf \Lambda} = \ol{\bf F} \ul{\bf F}^T =
{\bf I}_{\cal S}^{m,\hat{p}} {\bf \Lambda}_{\cal S}
\{ {\bf I}_{\cal S}^{r,\hat{p}} \}^T$, where ${\bf \Lambda}_{\cal S} =
diag \{ \lambda_{i(1)} \ldots \lambda_{i(\hat{p})} \} \in
{\cal R}^{\hat{p} \times \hat{p}}$, thus substituting for
$\ol{\bf F}$ and $\ul{\bf F}$ we obtain
$\ol{\bf C} \ul{\bf C}^T = {\bf \Lambda}_{\cal S}$.
Since ${\cal S} \neq \{ 1 , \ldots , p \}$ there exists an index
$i \not\in {\cal S}$, with $1 \leq i \leq p$.
In the sequel we will show that for any critical point
$(\ol{\bf W} , \ul{\bf W})$, which is not a global minimum,
there is an infinitesimal perturbation of the form
\begin{eqnarray*}
\ol{\bf W}' &=& \ol{\bf W} + \epsilon {\bf v}_i \ol{\bf z}^T
\\
\ul{\bf W}' &=& \ul{\bf W} + \epsilon {\bf q}_i \ul{\bf z}^T
\end{eqnarray*}
for some vectors $\ol{\bf z}$, $\ul{\bf z}$,
which achieves less error than $(\ol{\bf W} , \ul{\bf W})$
thus no critical point can be a local minimum.
We will treat the two cases $\hat{p} = p$, and $\hat{p} < p$,
separately. First assume $\hat{p} = p$, and let
$\ol{\bf z} = \ol{\bf C}^T {\bf e}_p^{(p)}$ and
$\ul{\bf z} = \alpha_i^{-1} \ul{\bf C}^T {\bf e}_p^{(p)}$. Then
\begin{eqnarray}
\lefteqn{J_{net}(\ol{\bf W}',\ul{\bf W}') - J_{net}(\ol{\bf W},\ul{\bf W}) =
\left\| {\bf \Lambda}_{\hat{r}} - {\bf V}^T \ol{\bf W}' {\ul{\bf W}'}^T {\bf Q}^{-T}
\left[ \begin{array}{c} {\bf D}_{1r} \\ 0 \end{array} \right]
\right\|_F^2 - \| {\bf \Lambda}_{\hat{r}} - {\bf \Lambda} \|_F^2 =}
\nonumber \\
&=& \left\| {\bf \Lambda}_{\hat{r}} - \left( \ol{\bf F}
+ \epsilon {\bf e}_i^{(m)} \{ {\bf e}_p^{(p)} \}^T \ol{\bf C} \right)
\left( \ul{\bf F}^T
+ \epsilon \ul{\bf C}^T {\bf e}_p^{(p)} \{ {\bf e}_i^{(r)} \}^T
\right) \right\|_F^2 -
\| {\bf \Lambda}_{\hat{r}} - {\bf \Lambda} \|_F^2 =
\nonumber \\
&=& \left\| {\bf \Lambda}_{\hat{r}} - \left( {\bf I}_{\cal S}^{m,p}
+ \epsilon {\bf e}_i^{(m)} \{ {\bf e}_p^{(p)} \}^T \right)
{\bf \Lambda}_{\cal S} \left( \{ {\bf I}_{\cal S}^{r,p} \}^T
+ \epsilon {\bf e}_p^{(p)} \{ {\bf e}_i^{(r)} \}^T
\right) \right\|_F^2
- \| {\bf \Lambda}_{\hat{r}} - {\bf \Lambda} \|_F^2 =
\nonumber \\
&=& \left\| {\bf \Lambda}_{\hat{r}} - {\bf \Lambda}
- \epsilon \lambda_{i(p)} {\bf I}_{\cal S}^{m,p}
{\bf e}_p^{(p)} \{ {\bf e}_i^{(r)} \}^T
- \epsilon \lambda_{i(p)} {\bf e}_i^{(m)} \{ {\bf e}_p^{(p)} \}^T
\{ {\bf I}_{\cal S}^{r,p} \}^T
\right.
\nonumber \\
&& \left. ~~
- \epsilon^2 \lambda_{i(p)} {\bf e}_i^{(m)} \{ {\bf e}_p^{(p)} \}^T
{\bf e}_p^{(p)} \{ {\bf e}_i^{(r)} \}^T
\right\|_F^2
- \| {\bf \Lambda}_{\hat{r}} - {\bf \Lambda} \|_F^2
\label{eq:JminusJ}
\end{eqnarray}

Since for any matrix ${\bf M}$, $\| {\bf M} \|_F^2 =
trace \{ {\bf M} {\bf M}^T \}$ and
$trace \{ {\bf M} \} = trace \{ {\bf M}^T \}$, we obtain
\begin{eqnarray*}
J_{net}(\ol{\bf W}',\ul{\bf W}') - J_{net}(\ol{\bf W},\ul{\bf W}) &=&
\| {\bf \Lambda}_{\hat{r}} - {\bf \Lambda} \|_F^2
\\
&& + \epsilon^2 \lambda_{i(p)}^2 \left\| {\bf I}_{\cal S}^{m,p}
{\bf e}_p^{(p)} \{ {\bf e}_i^{(r)} \}^T \right\|_F^2
\\
&& + \epsilon^2 \lambda_{i(p)}^2 \left\| {\bf e}_i^{(m)} \{ {\bf e}_p^{(p)} \}^T
\{ {\bf I}_{\cal S}^{r,p} \}^T \right\|_F^2
\\
&& - 2 \epsilon \lambda_{i(p)} ~trace \left\{ {\bf I}_{\cal S}^{m,p}
{\bf e}_p^{(p)} \{ {\bf e}_i^{(r)} \}^T
({\bf \Lambda}_{\hat{r}} - {\bf \Lambda})^T \right\}
\\
&& - 2 \epsilon \lambda_{i(p)} ~trace \left\{ {\bf e}_i^{(m)}
\{ {\bf e}_p^{(p)} \}^T \{ {\bf I}_{\cal S}^{r,p} \}^T
({\bf \Lambda}_{\hat{r}} - {\bf \Lambda})^T \right\}
\\
&& - 2 \epsilon^2 \lambda_{i(p)} ~trace \left\{
{\bf e}_i^{(m)} \{ {\bf e}_p^{(p)} \}^T
{\bf e}_p^{(p)} \{ {\bf e}_i^{(r)} \}^T
({\bf \Lambda}_{\hat{r}} - {\bf \Lambda})^T \right\}
\\
&& + {\cal O}(\epsilon^3)
- \| {\bf \Lambda}_{\hat{r}} - {\bf \Lambda} \|_F^2
\end{eqnarray*}

Notice that $\{ {\bf I}_{\cal S}^{r,p} \}^T
({\bf \Lambda}_{\hat{r}} - {\bf \Lambda})^T = 0$ and also
$trace \left\{ {\bf I}_{\cal S}^{m,p}
{\bf e}_p^{(p)} \{ {\bf e}_i^{(r)} \}^T
({\bf \Lambda}_{\hat{r}} - {\bf \Lambda})^T \right\} = $ \newline
$\lambda_i ~trace \left\{ {\bf e}_{i(p)}^{(p)}
\{ {\bf e}_i^{(m)} \}^T \right\} = 0$, since $i(p) \neq i$. So
\begin{equation}
J_{net}(\ol{\bf W}',\ul{\bf W}') - J_{net}(\ol{\bf W},\ul{\bf W}) =
2 \epsilon^2 \lambda_{i(p)} ( \lambda_{i(p)} - \lambda_i )
+ {\cal O}(\epsilon^3)
\label{JminusJ_final}
\end{equation}

Now, the ordered set ${\cal S}$ has $p$ elements
yet ${\cal S} \neq \{ 1 , \ldots , p \}$, so
$i \leq p < i(p)$, therefore $\lambda_{i(p)} < \lambda_i$
and $\lambda_{i(p)} - \lambda_i < 0$.
As $\epsilon \rightarrow 0$ the corresponding negative term
dominates in Eq.~\ref{JminusJ_final}, thus
$J_{net}(\ol{\bf W}',\ul{\bf W}') < J_{net}(\ol{\bf W},\ul{\bf W})$.
On the other hand, if $\hat{p} < p$, the matrices
$\ol{\bf C}$, $\ul{\bf C}$ have more columns than rows. Set
$\ol{\bf z}$ equal to any non-zero vector such that
$\ul{\bf C} \ol{\bf z} = 0$. Then $\ol{\bf z}$ is linearly independent
with respect to the rows of $\ol{\bf C}$;
indeed if this was not true then
there would exist a non-zero vector ${\bf g}$ such that
$\ol{\bf z}^T = {\bf g}^T \ol{\bf C} \Rightarrow
0 = \ol{\bf z}^T \ul{\bf C}^T = {\bf g}^T \ol{\bf C} \ul{\bf C}^T =
{\bf g}^T {\bf \Lambda}_{\cal S} \Rightarrow
{\bf g} = 0$, contradiction. It follows that $\ol{\bf z}
\not\perp \em{Null} \{ \ol{\bf C} \}$, hence there exist a
non-zero vector $\ul{\bf z} \in \em{Null} \{ \ol{\bf C} \}$
so that $\ol{\bf z}^T \ul{\bf z} > 0$.
Then Eq.~\ref{eq:JminusJ} becomes
\begin{eqnarray*}
J_{net}(\ol{\bf W}',\ul{\bf W}') - J_{net}(\ol{\bf W},\ul{\bf W}) &=&
\left\| {\bf \Lambda}_{\hat{r}} - {\bf \Lambda}
- \epsilon^2 \alpha_i {\bf e}_i^{(m)} \ol{\bf z}^T
\ul{\bf z} \{ {\bf e}_i^{(r)} \}^T \right\|_F^2
- \| {\bf \Lambda}_{\hat{r}} - {\bf \Lambda} \|_F^2 =
\\
&=& ( \lambda_i - \epsilon^2 \alpha_i \ol{\bf z}^T \ul{\bf z} )^2 - \lambda_i^2
\end{eqnarray*}

Again as $\epsilon \rightarrow 0$, we have
$( \lambda_i - \epsilon^2 \alpha_i \ol{\bf z}^T \ul{\bf z} )^2 - \lambda_i^2 < 0$,
so $J_{net}(\ol{\bf W}',\ul{\bf W}') < J_{net}(\ol{\bf W},\ul{\bf W})$.
Hence $(\ol{\bf W} , \ul{\bf W})$ is not a local minimum.

\hfill \qed

From the discussion in the proof it follows that
there exist two $p \times p$ matrices
$\ol{\bf C}$, $\ul{\bf C}$ such that $\ol{\bf C} \ul{\bf C}^T =
{\bf \Lambda}_p$ and ${\bf V}^T \ol{\bf W} = {\bf I}_p^m \ol{\bf C}$,~
$[{\bf D}_{1r} 0] {\bf Q}^{-1} \ul{\bf W} = {\bf I}_p^r \ul{\bf C}$~
where ${\bf I}_p^m$ ~(${\bf I}_p^r$) is a matrix comprised of the first
$p$ columns of the $m \times m$ ~($r \times r$) identity matrix.
The matrices $\ol{\bf C}$, $\ul{\bf C}$ are invertible
(since their product has full rank) so
$\ul{\bf C} = {\bf \Lambda}_p \ol{\bf C}^{-T}$,
and we can write the general solution for the upper- and lower-layer
weights as follows
\begin{eqnarray}
\ol{\bf W}^{*} &=& {\bf V}_p \ol{\bf C}
\\
\ul{\bf W}^{*} &=& ( {\bf Q}_p {\bf D}_{1p}^{-1} +
{\bf Q}_1 {\bf \Phi} ) {\bf \Lambda}_p \ol{\bf C}^{-T}
\end{eqnarray}

where ${\bf \Phi}$ is any $(n-r) \times p$ matrix.
If the input has full rank then $r=n$ and the optimal solution is simplified:
\begin{eqnarray}
\ol{\bf W}^{*} &=& {\bf V}_p \ol{\bf C}
\\
\ul{\bf W}^{*} &=& {\bf Q}_p {\bf D}_{1p}^{-1}
{\bf \Lambda}_p \ol{\bf C}^{-T} =\nonumber
\\
&=& {\bf Q}_p ~diag\Bigl\{\frac{\beta_1}{\alpha_1^2},\ldots,
\frac{\beta_p}{\alpha_p^2}\Bigr\}~
\ol{\bf C}^{-T}
\end{eqnarray}

\subsubsection{Linear Back-Propagation for extracting the component subspace}

Theorem \ref{th:APCA-net-local-minima} implies that
any gradient descent algorithm
for minimizing Eq.~\ref{eq:net_J} in a network with $p$ hidden units,
will produce the global minimum solution under assumption $(A)$.
An example of such an algorithm is
the back-propagation learning rule, the linear version of which is as follows:
\begin{eqnarray}
\ol{\bf W}_{k+1} &=& \ol{\bf W}_k +
\beta ~[{\bf y}_k - \ol{\bf W}_k {\bf a}_k] ~{\bf a}_k^T
\label{eq:BPp_ol} \\
\ul{\bf W}_{k+1} &=& \ul{\bf W}_k +
\beta ~{\bf x}_k ~[{\bf b}_k - \ol{\bf W}_k^T \ol{\bf W}_k {\bf a}_k]^T
\label{eq:BPp_ul}
\end{eqnarray}
where ${\bf a}_k = \ul{\bf W}_k^T {\bf x}_k$
and ${\bf b}_k = \ol{\bf W}_k^T {\bf y}_k$.
Therefore, $\ol{\bf W}_k \rightarrow \ol{\bf W}^{*}$ and
$\ul{\bf W}_k \rightarrow \ul{\bf W}^{*}$.
The {\em canonical solution}
for the weights would be the case were $\ol{\bf C} = {\bf I}$,
but the structure of a standard feed-forward network
does not provide any way of distinguishing it
from the other solutions, thus there is no reason why the system would
prefer it. Nevertheless, some extensions of the network structure,
using the concept of deflation or
by introducing lateral connections in the hidden layer, can
enforce preference to the canonical solution as discussed below.
The advantage of such a
structure is the flexibility of gradually increasing or decreasing
the rank of the network without retraining the old neurons.

\subsubsection{Linear BP for extracting the exact components}

\paragraph{Obtaining the first component.}

The linear BP rule for a two-layer network
with one hidden unit (Figure \ref{fig:bpmodel}a)
is capable of extracting the first component.
In this case the standard BP method becomes
\begin{eqnarray}
\ol{\bf w}_{k+1} &=& \ol{\bf w}_k +
\beta ~[{\bf y}_k - \ol{\bf w}_k a_k] ~a_k
\label{eq:bp_ol} \\
\ul{\bf w}_{k+1} &=& \ul{\bf w}_k +
\beta ~{\bf x}_k ~[b_k - \|\ol{\bf w}_k\|^2 a_k]
\label{eq:bp_ul}
\end{eqnarray}
where $a_k = \ul{\bf w}_k^T {\bf x}_k$
and $b_k = \ol{\bf w}_k^T {\bf y}_k$.
\begin{quote}
In the following we shall assume that the input has full rank
for simplicity. A similar discussion can be carried out for the reduced
rank case.
\end{quote}
Then according to Theorem \ref{th:APCA-net-local-minima}
the global minimum is achieved at the points
$(\ul{\bf w},~\ol{\bf w})=(\mu_1 {\bf q}_1,~\rho_1 {\bf v}_1)$
where $\mu_1 \rho_1 = \beta_1 / \alpha_1^2$.

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\psfig{figure=bp1.ps,height=3in} & \hs{0.2in}
\psfig{figure=bp_many.ps,height=3in}
\\
(a) & (b)
\end{tabular}
\end{center}
\caption{(a) The linear BP network with 1 hidden unit.
(b) The proposed model for many hidden units with
Lateral Orthogonalization Network.}
\label{fig:bpmodel}
\end{figure}

\paragraph{Special (symmetric) case: Self-Supervised BP}

A special case of the above result is the autoassociative
BP network, where ${\bf x}_k = {\bf y}_k$.
The fact that there is
only one hidden neuron forces the network to extract
the most efficient representation of the input so as to optimally
reconstruct it at the output, which leads the weights of the network
into a plain PCA solution. Here we'll show how this result
can be derived as a special case of our analysis.
Indeed, since ${\bf X}={\bf Y}$, Eqs.~\ref{eq:gsvd_x} and
\ref{eq:gsvd_yx} become the SVD of ${\bf X}$ and
we have ${\bf Q} = {\bf V}$, $\beta_i = \alpha_i^2$.
Thus ${\bf v}_1$ is
the principal component of ${\bf x}_k$ and
the net converges to the values
$\ul{\bf w}=\mu_1{\bf v}_1$, $\ol{\bf w}=\rho_1{\bf v}_1$,
where $\rho_1 \mu_1=\beta_1/\alpha_1^2=1$.

What is even more interesting is that the BP-learning rule
for the upper layer:
\begin{equation}
\ol{\bf w}_{k+1} = \ol{\bf w}_k+\beta~[{\bf x}_k - a_k \ol{\bf w}_k]~a_k
\label{eq:BP-sym}
\end{equation}

is very similar to the
normalized Hebbian rule,
except that $a_k = \ul{\bf w}^T {\bf x}_k$ rather than
$a_k = \ol{\bf w}^T {\bf x}_k$.
A similar observation has also been mentioned in \cite{Sanger89b}.
Nevertheless, since $\ul{\bf w}$ and $\ol{\bf w}$ are going to converge
to the same final vector (up to an arbitrary scaling factor),
we are motivated to change the BP rule by forcing
$\ul{\bf w}_k=\ol{\bf w}_k={\bf w}_k$ for all times $k$, and by updating
both upper and lower layer weights by the same rule
\begin{equation}
{\bf w}_{k+1} = {\bf w}_k+\beta ~[{\bf x}_k - a_k {\bf w}_k] ~a_k
\label{eq:oja}
\end{equation}

where now $a_k = {\bf w}_k^T {\bf x}_k$.
This rule is now exactly the same as Oja's rule.


\paragraph{Obtaining multiple components.}

Deflation is again the key for extracting the remaining components
after the first. This concept is extended in the APCA and it is
defined below.

\begin{definition}
Let ${\bf U} {\bf D}_1 {\bf Q}^{-1}$, ${\bf V} {\bf D}_2 {\bf Q}^{-1}$,
be the GSVD of a matrix pair $({\bf X}^T,{\bf Y}{\bf X}^T)$
and let the generalized eigenvalues be
$\lambda_1 \geq \lambda_2 \geq \ldots \geq \lambda_{\hat{r}} \geq 0$.
Then the mapping
\begin{equation}
{\bf Y}{\bf X}^T \rightarrow
\tilde{\bf Y}{\bf X}^T \equiv \left( {\bf Y} - \frac{\beta_1}{\alpha_1^2}
{\bf v}_1 {\bf q}_1^T {\bf X} \right) {\bf X}^T
\end{equation}
or
\begin{equation}
{\bf Y}{\bf X}^T \rightarrow
{\bf \bar{Y}}{\bf X}^T \equiv \left( {\bf Y} - {\bf v}_1 {\bf v}_1^T {\bf Y}
\right) {\bf X}^T
\end{equation}
or
\begin{equation}
{\bf Y}{\bf X}^T \rightarrow
{\bf Y}\tilde{\bf X}^T \equiv {\bf Y} \left( {\bf X}^T -
\frac{1}{\alpha_1^2}{\bf X}^T{\bf q}_1{\bf q}_1^T{\bf X}{\bf X}^T \right)
\end{equation}
where ${\bf v}_1$, ${\bf q}_1$ are the first columns of
${\bf V}$ and ${\bf Q}$ respectively, is called {\bf deflation
transformation} of $({\bf X}^T,{\bf Y}{\bf X}^T)$.
\end{definition}

This is an extension of the deflation transformation
introduced in Chapter \ref{ch:PCA}. Let us first
examine the transformation
${\bf Y}{\bf X}^T \rightarrow \tilde{\bf Y}{\bf X}^T$.
From Eq. \ref{eq:QXXQ.eq.D1D1} we get
\begin{eqnarray}
{\bf V}^T \tilde{\bf Y} {\bf X}^T {\bf Q} &=&
{\bf V}^T {\bf Y} {\bf X}^T {\bf Q} -
\frac{\lambda_1}{\alpha_1} {\bf V}^T {\bf v}_1 {\bf q}_1^T {\bf X}
{\bf X}^T {\bf Q} \nonumber
\\
&=& {\bf D}_2 - \frac{\beta_1}{\alpha_1^2}
\left[ \begin{array}{c} 1 \\ 0 \\ \vdots \\ 0 \end{array} \right]
[ 1~ 0~ \ldots ~ 0 ] {\bf D}_1^T {\bf D}_1
= {\bf \hat{D}}_2
\end{eqnarray}
where ${\bf \hat{D}}_2 = diag[0~\beta_2~\ldots~\beta_q]$. Then clearly
$({\bf U} {\bf D}_1 {\bf Q}^{-1},{\bf V} \hat{\bf D}_2 {\bf Q}^{-1})$
is a GSVD of the pair $({\bf X}^T,\tilde{\bf Y}{\bf X}^T)$
from which it follows that the generalized eigenvalues now are:
$\hat{\lambda}_1 = 0$, $\hat{\lambda}_2=\lambda_2$,
$\ldots$, $\hat{\lambda}_{\hat{r}}=\lambda_{\hat{r}}$.
So as in Chapter \ref{ch:PCA}, after the deflation,
the second component becomes dominant and it can be
extracted using exactly the same learning rule used for the first
component (applied now on the deflated data).
Note that the recursive nature of the above
argument lends itself easily to the extension to the
$p$-th component case ($p>1$).
Analogously, for the transformations
${\bf Y}{\bf X}^T \rightarrow {\bf \bar{Y}}{\bf X}^T$ and
${\bf Y}{\bf X}^T \rightarrow {\bf Y}\tilde{\bf X}^T$,
we can show again that
\begin{equation}
{\bf U}^T {\bf X}^T {\bf Q} = {\bf D}_1~\mbox{ and }~
{\bf V}^T {\bf \bar{Y}}{\bf X}^T {\bf Q} = {\bf \hat{D}}_2
\end{equation}
\begin{equation}
{\bf U}^T {\bf X}^T {\bf Q} = {\bf D}_1~\mbox{ and }~
{\bf V}^T {\bf Y}\tilde{\bf X}^T {\bf Q} = {\bf \hat{D}}_2
\end{equation}
which manifests deflation for the pairs
$({\bf X}^T,{\bf \bar{Y}}{\bf X}^T)$ and
$({\bf X}^T,{\bf Y}\tilde{\bf X}^T)$. Since 
all three transformations have the same GSVD
it follows that

\begin{equation}
\left( {\bf Y} - \frac{\beta_1}{\alpha_1^2}
{\bf v}_1 {\bf q}_1^T {\bf X} \right) {\bf X}^T
= \left( {\bf Y} - {\bf v}_1 {\bf v}_1^T {\bf Y} \right) {\bf X}^T
= {\bf Y} \left( {\bf X}^T - \frac{1}{\alpha_1^2}
{\bf X}^T{\bf q}_1{\bf q}_1^T{\bf X}{\bf X}^T \right)
\label{eq:deflation_equivalence}
\end{equation}

We can apply now the above results
straightforwardly into the back-propagation rule.
For the $p$-th neuron we assume that
the previous $p\!-\!1$ principal components
have already been extracted, so that
\[\ul{\bf w}_i = \mu_i {\bf q}_i, \]
\[\ol{\bf w}_i = \rho_i {\bf v}_i ~\mbox{ and} \]
\[\mu_i \rho_i = \beta_i / \alpha_i^2. \]
Then the modified BP rule for this neuron

\begin{eqnarray}
\ol{\bf w}_{p,k+1} &=& \ol{\bf w}_{pk} + \beta
\left[ {\bf y}_k - \sum_{i=1}^{p-1} \ol{\bf w}_i\ul{\bf w}_i^T{\bf x}_k
- \ol{\bf w}_{pk} a_{pk} \right] a_{pk}
\label{eq:bp_deflation_ol1}
\\
\ul{\bf w}_{p,k+1} &=& \ul{\bf w}_{pk} + \beta {\bf x}_k
\left[ {\bf y}_k - \sum_{i=1}^{p-1} \ol{\bf w}_i\ul{\bf w}_i^T{\bf x}_k
- \ol{\bf w}_{pk} a_{pk} \right]^T \ol{\bf w}_{pk}
\label{eq:bp_deflation_ul1}
\end{eqnarray}
where $a_{pk} = \ul{\bf w}_{pk}^T {\bf x}_k$,
is equivalent to a single-hidden-unit
BP rule with a modified teacher:

\[
\tilde{\bf y}_k = {\bf y}_k - \sum_{i=1}^{p-1} \ol{\bf w}_i \ul{\bf w}_i^T
{\bf x}_k.
\]
It is a simple exercise to show that the data matrix
$\tilde{\bf Y} = [\tilde{\bf y}_1, \ldots, \tilde{\bf y}_N]$,
for the modified teacher is the $(p\!-\!1)$-times deflation of
${\bf Y}$ hence the $p$-th component is now dominant and
will be extracted by the network.
So $\ul{\bf w}_p$ and $\ol{\bf w}_p$ will converge to
$\mu_p {\bf q}_p$ and $\rho_p {\bf v}_p$
with $\mu_p \rho_p = \beta_p / \alpha_p^2$.
Another version (Version 2) of the modified BP rule corresponding to
the deflation transformation of the second type is

\begin{eqnarray}
\ol{\bf w}_{p,k+1} &=& \ol{\bf w}_{pk} + \beta
\left[ {\bf y}_k - \sum_{i=1}^{p-1} \ol{\bf w}_i \ol{\bf w}_i^T {\bf y}_k
/ \|\ol{\bf w}_i\|^2
- \ol{\bf w}_{pk} a_{pk} \right] a_{pk}
\label{eq:bp_deflation_ol2}
\\
\ul{\bf w}_{p,k+1} &=& \ul{\bf w}_{pk} + \beta {\bf x}_k
\left[ {\bf y}_k - \sum_{i=1}^{p-1} \ol{\bf w}_i \ol{\bf w}_i^T {\bf y}_k
/ \|\ol{\bf w}_i\|^2
- \ol{\bf w}_{pk} a_{pk} \right]^T \ol{\bf w}_{pk}
\label{eq:bp_deflation_ul2}
\end{eqnarray}

Obviously the modified teacher

\begin{equation}
{\bf \bar{y}}_k = {\bf y}_k - \sum_{i=1}^{p-1}
\frac{\ol{\bf w}_i \ol{\bf w}_i^T}{\|\ol{\bf w}_i\|^2} {\bf y}_k
= {\bf y}_k - \sum_{i=1}^{p-1} {\bf v}_i {\bf v}_i^T {\bf y}_k
\end{equation}

is equal to the $(p\!-1\!)$-times deflated original teacher.
It follows that the weights will again converge to the $p$-th
asymmetric principal components.
Finally, Version 3 of the modified BP

\begin{eqnarray}
\ol{\bf w}_{p,k+1} &=& \ol{\bf w}_{pk} + \beta
\left[ {\bf y}_k{\bf x}_k^T\ul{\bf w}_{pk} - {\bf y}_k{\bf x}_k^T
\sum_{i=1}^{p-1} \frac{\ul{\bf w}_i{\bf w}_i^T{\bf X}{\bf X}^T}
{\ul{\bf w}_i^T{\bf X}{\bf X}^T\ul{\bf w}_i} \ul{\bf w}_{pk}
- \ol{\bf w}_{pk} a_{pk}^2 \right]
\label{eq:bp_deflation_ol3}
\\
\ul{\bf w}_{p,k+1} &=& \ul{\bf w}_{pk} + \beta \left[
\left( {\bf x}_k - \sum_{i=1}^{p-1} \frac{{\bf X}{\bf X}^T
\ul{\bf w}_i\ul{\bf w}_i}{\ul{\bf w}_i^T{\bf X}{\bf X}^T\ul{\bf w}_i}
{\bf x}_k \right) {\bf y}_k^T - {\bf x}_k \|\ul{\bf w}_{pk}\|^2
a_{pk} \right]
\label{eq:bp_deflation_ul3}
\end{eqnarray}
corresponds to the deflation of the third type.

Version 1 (Eqs.~\ref{eq:bp_deflation_ol1},
\ref{eq:bp_deflation_ul1}) is more aesthetically pleasing
since the modified teacher $\tilde{\bf y}$,
is nothing but the original teacher after subtracting the
reconstruction that has been achieved so far by the previous
$p\!-\!1$ neurons.
This is intuitively in agreement with the concept of deflation
which is perceived as a subtraction process of the previous
components.


\section{The Lateral Orthogonalization Net}

It has become apparent so far that orthogonality is a valuable
property in all PCA problems, standard, oriented or asymmetric.
It has also been demonstrated that in all these cases
orthogonality can be achieved via a deflation transformation
assuming that all previous components have been extracted.
Sanger followed this approach in the standard PCA case \cite{Sanger89a}.
However we showed in Chapter \ref{ch:PCA} that deflation can
be essentially implemented in the network using {\em
lateral connections} trained by the normalized Hebbian rule.
This approach not only replaces the artificial explicit
deflation with a more natural in-network deflation,
but also achieves computational advantage over it.
Furthermore, it turns out that the {\em Lateral Orthogonalization
Network (LON)} can be extended to the APCA case as well.
Two kinds of learning rules can be used to train the LON:

\begin{description}
\item[Local Orthogonalization Rule.]
Let $a_{1k}$ and $a_{2k}$ be two stochastic processes
(e.g. the activation sequences of two different neurons)
the {\em inner product} of which is defined as
\[
<a_1,a_2> \equiv E \{a_{1k} a_{2k}\}
\]
The role of the lateral network is to find a weight strength $c$,
such that the modified process
\begin{equation}
a_{2k}' = a_{2k} - c~a_{1k}
\label{eq:olr1}
\end{equation}
is orthogonal to $a_{1k}$, i.e.
\begin{equation}
<a_2',a_1> = 0
\end{equation}
Assuming stationarity the value $c$ is not a function of time.
By Gram-Schmidt orthogonalization \cite{Golub89}
the solution is
\begin{equation}
c^{*}=\frac{<a_2,a_1>}{<a_1,a_1>}=\frac{<a_1,a_2>}{<a_1^2>}
\end{equation}
which means that the component needed to be removed
from $a_2$ in order to achieve orthogonality is exactly
the projection of $a_2$ onto $a_1$ (see Fig.~\ref{fig:LON}).

The following orthogonalization learning rule is then proposed
for tracking the optimal value $c^{*}$:
\begin{equation}
c_{k+1} = c_k + \beta ~\left[ a_{1k} a_{2k}
 - c_k a_{1k}^2 \right]
\label{eq:olr-adaptive}
\end{equation}
Assuming that the system with which this rule is combined is stable,
then from stochastic approximation theory\footnote{
Refer to \cite{Ljung77,Kushner78} and also to the APEX
Proof in Chapter \ref{ch:PCA}.} it follows that
Eq.~\ref{eq:olr-adaptive} approximates the differential equation
\begin{equation}
dc'/dt=E\{a_1 a_2\} - c'~E\{a_1^2\}
\end{equation}
which in turn tracks the value of $c^{*}=\frac{E\{a_1 a_2\}}{E\{a_1^2\}}$.
We call this orthogonalization rule local, because the values
$a_1$ and $a_2$ utilized for updating the weight strength $c$,
are local to the lateral synaptic connection: they are just the
activations of the neurons at both ends of the connection.

\begin{figure}
\centerline{\psfig{figure=orthogonal.ps,width=2.5in}
\hspace{0.5in}
\psfig{figure=projection.ps,width=2.5in}}
\caption{The lateral orthogonalization network is based on the principle of
removing any duplication of old components from a new component.
The part to be removed in order to achieve orthogonality may be depicted as 
the projection from the new component to the old component.
As shown, after the removal, the residue is orthogonal to the old component.}
\label{fig:LON}
\end{figure}

\item[Immediate Orthogonalization Rule.]

In some cases the lateral orthogonalization weights can be trained
by a specialized rule which depends on the
dynamics of the specific system. No general equation can be prescribed
for such a rule since it is usually different
for each different system it is applied to.
However such a rule is usually easy to produce once
the quantities to become orthogonal have been identified
and the system dynamics have been set up.
One example is the APEX lateral orthogonalization
network described in Chapter \ref{ch:PCA}.
The rule proposed there for the lateral weight $c$, was of the form
\begin{equation}
c_{k+1}=c_k+\beta~\left[a_{1k} a_{2k}-c_k a_{2k}^2 \right]
\label{eq:APEX_LON}
\end{equation}
for example, between neurons 1 and 2. In this case $c$
is equal to ${\bf w}_{2k}^T {\bf w}_1$ while
the adaptive rule for ${\bf w}_2$ is
\[
{\bf w}_{2,k+1}={\bf w}_{2k}+\beta~\left[
a_{2k} {\bf x}-{\bf w}_{2k} a_{2k}^2 \right]
\]
therefore a simple multiplication from the left by ${\bf w}_1^T$
yields Eq.~\ref{eq:APEX_LON}. Clearly, had the updating rule
for ${\bf w}_2$ been different we would have obtained a different
adaptive rule for $c$. So this rule is case-dependent.
It is also immediate because $c_k$ equals
the inner product ${\bf w}_{2k}^T {\bf w}_1$ at all times
(assuming $c_0={\bf w}_{2,0}^T {\bf w}_1$)
rather than tracking it through some sort of dynamic equation.
\end{description}

\subsection{Deflation using LON}

\subsubsection{Explicit deflation and stochastic approximation}

Under the stochastic approximation formulation introduced in
Chapter \ref{ch:PCA} the back-propagation learning
Eqs.~\ref{eq:bp_ol}, \ref{eq:bp_ul}
correspond to the following differential equations
\begin{eqnarray}
\frac{d\ol{\bf w}(t)}{dt} &=& {\bf R}_{yx} \ul{\bf w}(t)
- \ol{\bf w}(t) \ul{\bf w}(t)^T {\bf R}_x \ul{\bf w}(t)
\label{eq:bp_ode_ol}
\\
\frac{d\ul{\bf w}(t)}{dt} &=& {\bf R}_{yx}^T \ol{\bf w}(t)
- \| \ol{\bf w}(t) \|^2 {\bf R}_x \ul{\bf w}(t)
\label{eq:bp_ode_ul}
\end{eqnarray}
where ${\bf R}_{yx}=E\{{\bf y}{\bf x}^T\}$ and
${\bf R}_x=E\{{\bf x}{\bf x}^T\}$. If there are finite many training
samples ($N<\infty$) then we can construct the infinite sequences
${\bf x}_k$, ${\bf y}_k$,
by periodically repeating the data. In this case
\[
{\bf R}_x=\frac{1}{N}{\bf X}{\bf X}^T,
~~~ {\bf R}_{yx}=\frac{1}{N}{\bf Y}{\bf X}^T
\]
On the other hand in the asymptotic case $N \rightarrow \infty$
\[
{\bf R}_x=\lim_{N\rightarrow \infty} \frac{1}{N}{\bf X}{\bf X}^T,
~~~ {\bf R}_{yx}=\lim_{N\rightarrow \infty} \frac{1}{N}{\bf Y}{\bf X}^T
\]
In either case from Eqs.~\ref{eq:QXXQ.eq.D1D1}, \ref{eq:gsvd_yx}
we have
\begin{eqnarray}
{\bf Q}^T {\bf R}_x {\bf Q} &=& {\bf \Sigma}_1 =
\frac{1}{N} {\bf D}_1^T {\bf D}_1
\label{eq:QRxQ.eq.S1}
\\
{\bf V}^T {\bf R}_{yx} {\bf Q} &=& {\bf \Sigma}_2 =
\frac{1}{N} {\bf D}_2 ~.
\label{eq:VRyxQ.eq.S2}
\end{eqnarray}
${\bf \Sigma}_1 \in {\cal R}^{n\times n}$ and
${\bf \Sigma}_2 \in {\cal R}^{m\times n}$ are diagonal matrices with
diagonal entries $\bar{\alpha}_i^2=\frac{1}{N}\alpha_i^2$, $i=1,\ldots,n$,
and $\bar{\beta}_i=\frac{1}{N}\beta_i$, $i=1,\ldots,\min\{m,n\}$
(apply limits when appropriate).

By construction of the back-propagation rule,
Eqs.~\ref{eq:bp_ode_ol}, \ref{eq:bp_ode_ul},
form a gradient descent procedure on the surface
of the error function
\[
J=E\{\|{\bf y}-\ol{\bf w}\ul{\bf w}^T {\bf x}\|^2\}
\]
Applying the results of the previous section asymptotically
we can say that $J$ has no local minima (under assumption $(A)$)
and the minimum point is $(\ul{\bf w}^{*},\ol{\bf w}^{*})$=
$(\mu_1{\bf q}_1,\rho_1{\bf v}_1)$ with
$\mu_1 \rho_1=\bar{\beta}_1/\bar{\alpha}_1^2$.
Furthermore, gradient descent is guaranteed to converge to
the global minimum so
$(\ul{\bf w}(t),\ol{\bf w}(t)) \rightarrow (\ul{\bf w}^{*},
\ol{\bf w}^{*})$ as $t \rightarrow \infty$.

For the $p$-th neuron we assume as usual that
\[\ul{\bf w}_i = \mu_i {\bf q}_i, \]
\[\ol{\bf w}_i = \rho_i {\bf v}_i ~\mbox{ and} \]
\[\mu_i \rho_i = \bar{\beta}_i / \bar{\alpha}_i^2. \]
Then Versions 1, 2 and 3 of the modified BP correspond
to the following differential equations

Version 1:
\begin{eqnarray}
\frac{d\ol{\bf w}_p(t)}{dt} &=& \left[ {\bf R}_{yx}
-\sum_{i<p} \ol{\bf w}_i \ul{\bf w}_i^T {\bf R}_x \right] \ul{\bf w}_p(t)
- \ol{\bf w}_p(t) \ul{\bf w}_p(t)^T {\bf R}_x \ul{\bf w}_p(t)
\label{eq:bp_deflation_ode_ol1}
\\
\frac{d\ul{\bf w}_p(t)}{dt} &=& \left[ {\bf R}_{yx}^T
-{\bf R}_x \sum_{i<p} \ul{\bf w}_i \ol{\bf w}_i^T \right]
\ol{\bf w}_p(t)
- \| \ol{\bf w}_p(t) \|^2 {\bf R}_x \ul{\bf w}_p(t)
\label{eq:bp_deflation_ode_ul1}
\end{eqnarray}

Version 2:
\begin{eqnarray}
\frac{d\ol{\bf w}_p(t)}{dt} &=& \left[
{\bf I}-\sum_{i<p} \ol{\bf w}_i \ol{\bf w}_i^T/\|\ol{\bf w}_i\|^2
\right] {\bf R}_{yx} \ul{\bf w}_p(t)
- \ol{\bf w}_p(t) \ul{\bf w}_p(t)^T {\bf R}_x \ul{\bf w}_p(t)
\label{eq:bp_deflation_ode_ol2}
\\
\frac{d\ul{\bf w}_p(t)}{dt} &=& {\bf R}_{yx}^T \left[
{\bf I}-\sum_{i<p} \ol{\bf w}_i \ol{\bf w}_i^T/\|\ol{\bf w}_i\|^2
\right] \ol{\bf w}_p(t)
- \| \ol{\bf w}_p(t) \|^2 {\bf R}_x \ul{\bf w}_p(t)
\label{eq:bp_deflation_ode_ul2}
\end{eqnarray}

Version 3:
\begin{eqnarray}
\frac{d\ol{\bf w}_p(t)}{dt} &=& \left[ {\bf R}_{yx}
- {\bf R}_{yx} \sum_{i<p} \frac{\ul{\bf w}_i \ul{\bf w}_i^T {\bf R}_x}
{\ul{\bf w}_i^T {\bf R}_x \ul{\bf w}_i} \right] \ul{\bf w}_p(t)
- \ol{\bf w}_p(t) \ul{\bf w}_p(t)^T {\bf R}_x \ul{\bf w}_p(t)
\label{eq:bp_deflation_ode_ol3}
\\
\frac{d\ul{\bf w}_p(t)}{dt} &=& \left[ {\bf R}_{yx}^T
- \sum_{i<p} \frac{{\bf R}_x \ul{\bf w}_i \ul{\bf w}_i^T}
{\ul{\bf w}_i^T {\bf R}_x \ul{\bf w}_i} {\bf R}_{yx}^T
\right] \ol{\bf w}_p(t)
- \| \ol{\bf w}_p(t) \|^2 {\bf R}_x \ul{\bf w}_p(t)
\label{eq:bp_deflation_ode_ul3}
\end{eqnarray}

Extending the equivalence shown in
Eq.~\ref{eq:deflation_equivalence} we have
\[
{\bf R}_{yx}
-\sum_{i<p} \ol{\bf w}_i \ul{\bf w}_i^T {\bf R}_x
= \left[
{\bf I}-\sum_{i<p} \ol{\bf w}_i \ol{\bf w}_i^T/\|\ol{\bf w}_i\|^2
\right] {\bf R}_{yx}
= {\bf R}_{yx}
- {\bf R}_{yx} \sum_{i<p} \frac{\ul{\bf w}_i \ul{\bf w}_i^T {\bf R}_x}
{\ul{\bf w}_i^T {\bf R}_x \ul{\bf w}_i}
\]
so
\[
Eq.~\ref{eq:bp_deflation_ode_ol1} \Longleftrightarrow
Eq.~\ref{eq:bp_deflation_ode_ol2} \Longleftrightarrow
Eq.~\ref{eq:bp_deflation_ode_ol3}
\]
\[
Eq.~\ref{eq:bp_deflation_ode_ul1} \Longleftrightarrow
Eq.~\ref{eq:bp_deflation_ode_ul2} \Longleftrightarrow
Eq.~\ref{eq:bp_deflation_ode_ul3}
\]
thus the three sets of differential equations are completely
equivalent. It follows by reverse argument that any one
of Eqs.~\ref{eq:bp_deflation_ol1}, \ref{eq:bp_deflation_ol2}
or \ref{eq:bp_deflation_ol3},
combined with any one of Eqs.~\ref{eq:bp_deflation_ul1},
\ref{eq:bp_deflation_ul2} or \ref{eq:bp_deflation_ul3},
corresponds to the same set
of differential equations, thus any of the 8 combinations
should extract the $p$-th component.

\subsubsection{LON deflation}

Contrary to the standard PCA problem and APEX model,
for the linear approximation APCA we will
need two lateral connections for each pair of hidden neurons:
one for controlling the orthogonality of the upper-layer weights
and one for the lower-layer weights (see Figure \ref{fig:bpmodel}b).
Assuming as usual that the first
$p\!-\!1$ components have already been extracted the proposed model
for the $p$th neuron is defined below

\begin{eqnarray}
\ol{\bf w}_{p,k+1} &=& \ol{\bf w}_k + \beta
\left[ {\bf y}_k a'_{pk} - \ol{\bf w}_{pk} a_{pk}^2 \right]
\label{eq:bp_LON_ol}
\\
\ul{\bf w}_{p,k+1} &=& \ul{\bf w}_k + \beta~
{\bf x}_k \left[ b'_{pk} - \|\ol{\bf w}_{pk}\|^2 a_{pk} \right]
\label{eq:bp_LON_ul}
\end{eqnarray}
where
\begin{eqnarray}
a'_{pk} = a_{pk} - \sum_{i<p} \ul{c}_{pi} a_{ik},
& ~~ &
a_{ik} = \ul{\bf w}_i^T {\bf x}_k, ~i<p,
~a_{pk} = \ul{\bf w}_{pk}^T {\bf x}_k
\label{eq:a_prime}
\\
b'_{pk} = b_{pk} - \sum_{i<p} \ol{c}_{pi} b_{ik}
& ~~ &
b_{ik} = \ol{\bf w}_i^T {\bf y}_k, ~i<p,
~b_{pk} = \ol{\bf w}_{pk}^T {\bf y}_k
\label{eq:b_prime}
\end{eqnarray}
Indeed, expanding we obtain the following associated differential
equations

\begin{eqnarray}
\frac{d\ol{\bf w}_p(t)}{dt} &=& {\bg R}_{yx}\ul{\bf w}_p(t)
- {\bf R}_{yx}\ul{\bf w}_i\ul{\bf c}_{pi}(t)
- \ol{\bf w}_p(t)\ul{\bf w}_p(t)^T{\bf R}_x\ul{\bf w}_p(t)
\\
\frac{d\ul{\bf w}_p(t)}{dt} &=& {\bg R}_{yx}^T\ol{\bf w}_p(t)
- {\bf R}_{yx}^T\ol{\bf w}_i\ol{\bf c}_{pi}(t)
- \|\ol{\bf w}_p(t)\|^2 {\bf R}_x \ul{\bf w}_p(t)
\end{eqnarray}

Ideally, $\ol{c}_{pi}$ and $\ul{c}_{pi}$ should be equal to
\begin{equation}
\ol{c}_{pi}^{*} = \frac{\ol{\bf w}_i^T \ol{\bf w}_p}
{\|\ol{\bf w}_i\|^2} ~~~\mbox{and}~~~
\ul{c}_{pi}^{*} = \frac{\ul{\bf w}_i^T {\bf R}_x \ul{\bf w}_p}
{\ul{\bf w}_i^T {\bf R}_x \ul{\bf w}_i}
\end{equation}
so that the above differential equations
will become the same as
Eqs.~\ref{eq:bp_deflation_ode_ol2} and \ref{eq:bp_deflation_ode_ul3}.
As we already saw this implies that the $p$-th component
is extracted from the network. However, the optimal values
$\ol{c}_{pi}^{*}$, $\ul{c}_{pi}^{*}$ are not available immediately
and therefore must be tracked\footnote{
More precisely,
the values of $\ol{\bf w}_i$ and $\ol{\bf w}_{pk}$ are
immediately available so $\ol{c}_{pi}=\ol{\bf w}_i^T\ol{\bf w}_{pk}/
\|\ol{\bf w}_i\|^2$ can be adapted
by an immediate orthogonalization rule. Indeed,
multiplying Eq.~\ref{eq:bp_LON_ol} from the left by
$\ol{\bf w}_i^T/\|\ol{\bf w}_i\|^2$ one obtains:
\[
\ol{c}_{pi,k+1} = \ol{c}_{pi,k} + \beta~\left[
b_{ik} a_{pk}'/\|\ol{\bf w}_i\|^2 - \ol{c}_{pi,k} a_{pk}^2 \right]
\]
However the value of $\ul{c}_{pi}^{*}$ is not readily available
because of the involvement of ${\bf R}_x$ which is never explicitly
computed by the network. A time average approach is needed
to track $\ul{c}_{pi}^{*}$ thus making the local orthogonalization
rule necessary. For sake of uniformity we'll adopt the local
rule for both $\ol{c}_{pi}$ and $\ul{c}_{pi}$.}.
That is done by the local type orthogonalization rules

\begin{eqnarray}
\ol{c}_{pi,k+1} &=& \ol{c}_{pi,k} + \beta~
(\ol{\bf w}_i^T \ol{\bf w}_{pk} - \ol{c}_{pi,k} \|\ol{\bf w}_i\|^2)
\label{eq:bp_LON_olc}
\\
\ul{c}_{pi,k+1} &=& \ul{c}_{pi,k} + \beta~
(a_{ik} a_{pk} - \ul{c}_{pi,k} a_{ik}^2)
\label{eq:bp_LON_ulc}
\end{eqnarray}
with associated differential equations
\begin{eqnarray*}
\frac{d\ol{c}_{pi}(t)}{dt} &=& \ol{\bf w}_i^T \ol{\bf w}_p(t)
- \ol{c}_{pi}(t) \|\ol{\bf w}_i\|^2
\\
\frac{d\ul{c}_{pi}(t)}{dt} &=& \ul{\bf w}_i^T {\bf R}_x \ul{\bf w}_p(t)
- \ul{c}_{pi}(t) \ul{\bf w}_i^T {\bf R}_x \ul{\bf w}_i
\end{eqnarray*}
which make evident that
the local orthogonalization rules for
$\ol{c}_{pi}$ and $\ul{c}_{pi}$ indeed track the values $\ol{c}_{pi}^{*}$
and $\ul{c}_{pi}^{*}$ as required. We'll show in the next
section via simulation
that the LON network with Eqs.~\ref{eq:bp_LON_ol},
\ref{eq:bp_LON_ul}, \ref{eq:bp_LON_olc} and \ref{eq:bp_LON_ulc},
indeed extracts the $p$-th component (assuming the previous
components have been already extracted) as theoretically predicted.


\subsubsection{Simulation Results}

Figure \ref{fig:bp-sim}a shows a simulation experiment
on multiple components using lateral orthogonalization
network described above. The components are extracted sequentially
The $y$-axis corresponds to the component estimation error
while the $x$-axis corresponds to the number of iterations.
The data used are 100 artificially produced
(${\bf x},{\bf y}$) pairs
from two colored random sequences,
repeated cyclically in sweeps
where the input dimension is $n=4$ and the output dimension is $m=3$.
Unlike APEX the optimal value of $\beta$ is not identified
so we experimentally found $\beta=0.01$ to be good for our
data and this is the value we use in our experiments.
The curves in Fig.~\ref{fig:bp-sim}a are superimposed
and it is understood that the second unit is trained 
after the first one, the third unit after the second, etc.
We found that for larger input and output dimensions
the convergence of
the network can be considerably slower than depicted in Figure
\ref{fig:bp-sim}. However, even if slowly,
the network eventually always converges to the appropriate
solution.

\begin{figure}
\begin{center}
\begin{tabular}{c}
\psfig{figure=bplon_up.ps,width=4in}
\\
(a)
\\
\psfig{figure=bplon_lo.ps,width=4in}
\\
(b)
\end{tabular}
\end{center}
\caption{Convergence of the linear approximation APCA model
with 3 hidden units using LON. (a) Plot of the
component estimation error
$\protect\| {\bf v}_p-\protect\frac{
\protect\ul{\bf w}\protect_p}{\protect\|
\protect\ul{\bf w}\protect_p \protect\|} \protect\|^2$ of the upper layer.
(b) Plot of the component estimation error
$\protect\| \protect\frac{{\bf q}\protect_p}{\protect\|
{\bf q}\protect_p \protect\|}-\protect\frac{
\protect\ul{\bf w}\protect_p}{\protect\|
\protect\ul{\bf w}\protect_p \protect\|} \protect\|^2$
of the lower layer.}
\label{fig:bp-sim}
\end{figure}


\section{Neural Nets for Cross-Correlation  APCA Problem}

As with the linear approximation APCA problem the cross-correlation
APCA problem involves
two random signals ${\bf x}_k \in {\cal R}^n$ and
${\bf y}_k \in {\cal R}^m$,
with at least wide sense stationary statistics and autocorrelation matrices
${\bf R}_x$ and ${\bf R}_y$ respectively,
however, now
we are interested in extracting the principal singular vectors
of the cross-correlation matrix ${\bf R}_{xy}=E\{{\bf x}{\bf y}^T\}$. Let
\begin{equation}
{\bf R}_{yx}={\bf U} {\bf \Sigma} {\bf V}^T
\end{equation}
be the SVD of ${\bf R}_{yx}$
where ${\bf \Sigma}=diag[\sigma_1,\ldots,\sigma_q]$, $q=\min\{m,n\}$,
${\bf U}=[{\bf u}_1,\ldots,{\bf u}_m]$, ${\bf V}=[{\bf v}_1,
\ldots,{\bf v}_n]$.
Our goal is to find the left and right singular vectors ${\bf u}_i$
and ${\bf v}_i$, $i=1,2,\ldots$ that correspond to the largest
singular values $\sigma_1,\sigma_2,\ldots$.


\subsubsection{Single Component Case}

Like any other PCA problem we have discussed so far
our analysis will start from the
single component case before extending it to the general case.
Figure \ref{fig:crosscor}(a) depicts the proposed network.
The objective is to maximize the correlation
\begin{equation}
\mbox{maximize } J=\frac{E\{a_k b_k\}}{\|\ol{\bf w}\|\|\ul{\bf w}\|}
\end{equation}
between the outputs of the two units
\begin{eqnarray}
a_k &\equiv& \ul{\bf w}_k^T {\bf x}_k
\\
b_k &\equiv& \ol{\bf w}_k^T {\bf y}_k
\end{eqnarray}
The problem is equivalent to maximizing the ratio
\begin{equation}
J=\frac{\ol{\bf w}^T {\bf R}_{yx} \ul{\bf w}}{\|\ol{\bf w}\|\|\ul{\bf w}\|}
\end{equation}
The solution is known to be $(\ol{\bf w},\ul{\bf w})=
(\rho {\bf u}_1,\mu {\bf v}_1)$ where $\rho$ and $\mu$ are any scalars
other than zero \cite{Golub89}, namely the optimal weights
are the principal left and right singular vectors of ${\bf R}_{yx}$
probably not normalized. We obviously have $\|\ol{\bf w}\|^2=\rho^2$
and $\|\ul{\bf w}\|^2=\mu^2$ but the maximum value of $J$
does not depend on $\rho$ or $\mu$ and it is equal to the maximum
singular value
\begin{equation}
J_{max}=\sigma_1.
\end{equation}

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\psfig{figure=cross1.ps,width=2in} & \hs{0.2in}
\psfig{figure=cross_many.ps,width=2in}
\\
(a) & (b)
\end{tabular}
\end{center}
\caption{(a) A single hidden unit cross-correlation network.
(b) A multiple hidden units cross-correlation network with
lateral inhibition parameters.}
\label{fig:crosscor}
\end{figure}

We propose the following dual set of learning rules for extracting
the principal cross-correlation components:
\begin{eqnarray}
\ol{\bf w}_{k+1} &=& \ol{\bf w}_k + \beta ~\left[ {\bf y}_k
- \ol{\bf w}_k b_k \right] ~a_k
\label{eq:crosscor_ol}
\\
\ul{\bf w}_{k+1} &=& \ul{\bf w}_k + \beta ~\left[ {\bf x}_k
- \ul{\bf w}_k a_k \right] ~b_k
\label{eq:crosscor_ul} 
\end{eqnarray}

We will show that the algorithm converges with probability 1 to
the desired solution and both final vectors are normalized.
This is a useful feature for the deflation procedure discussed
later for extracting the remaining components.

\begin{theorem}
\label{th:APCA_cross_correlation}
Consider the algorithm described by Eqs.~\ref{eq:crosscor_ol}
and \ref{eq:crosscor_ul} under the assumptions discussed above
and further assume that $\sigma_1 > \sigma_2 \geq \sigma_3 \geq \ldots
\sigma_q \geq 0$. Then with probability 1,
$\ol{\bf w}_k\rightarrow \pm {\bf u}_1$ and
$\ul{\bf w}_k\rightarrow \pm {\bf v}_1$ as $k\rightarrow \infty$.
\end{theorem}
{\bf Proof.}
As in Chapter \ref{ch:PCA} we refer to the stochastic approximation
theory the results of which allow us to associate the stochastic
algorithm with a differential equation in continuous time.
The associated equations in this case are
\begin{eqnarray}
\frac{d\ol{\bf w}(t)}{dt} &=& {\bf R}_{yx} \ul{\bf w}(t)
- \ol{\bf w}(t) \left( \ol{\bf w}(t)^T {\bf R}_{yx} \ul{\bf w}(t)
\right)
\label{eq:crosscor_ode_ol}
\\
\frac{d\ul{\bf w}(t)}{dt} &=& {\bf R}_{yx}^T \ol{\bf w}(t)
- \ul{\bf w}(t) \left( \ol{\bf w}(t)^T {\bf R}_{yx} \ul{\bf w}(t)
\right)
\label{eq:crosscor_ode_ul} 
\end{eqnarray}

We will show the global asymptotic stability of the set
$\{({\bf u}_1,{\bf v}_1),(-{\bf u}_1,-{\bf v}_1)$
by establishing an appropriate Lyapunov function for the
dynamic equations at hand. Let
\begin{equation}
V(\ol{\bf w},\ul{\bf w}) = \sigma_1 - J =
\sigma_1-
\frac{\ol{\bf w}^T{\bf R}_{yx}\ul{\bf w}}{\|\ol{\bf w}\|\|\ul{\bf w}\|}
\label{eq:cross_Lyapunov}
\end{equation}
Clearly $V(\ol{\bf w},\ul{\bf w})\geq 0$ for all
$\ol{\bf w}$, $\ul{\bf w}$. Furthermore,

\begin{eqnarray}
\frac{dV}{dt} &=& -\frac{d\ol{\bf w}^T}{dt}\nabla_{\ol{w}}J 
-\frac{d\ul{\bf w}^T}{dt}\nabla_{\ul{w}}J = \nonumber
\\
&=& -\Bigl( \ul{\bf w}^T{\bf R}_{yx}^T-(\ol{\bf w}^T{\bf R}_{yx}\ul{\bf w})
\ol{\bf w}^T \Bigr)
\left( \frac{{\bf R}_{yx}\ul{\bf w}}{\|\ol{\bf w}\|\|\ul{\bf w}\|}
-\frac{(\ol{\bf w}^T{\bf R}_{yx}\ul{\bf w})\ol{\bf w}}
{\|\ol{\bf w}\|^3\|\ul{\bf w}\|}\right) \nonumber
\\
&& -\Bigl( \ol{\bf w}^T{\bf R}_{yx}-(\ol{\bf w}^T{\bf R}_{yx}\ul{\bf w})
\ul{\bf w}^T \Bigr)
\left( \frac{{\bf R}_{yx}^T\ol{\bf w}}{\|\ol{\bf w}\|\|\ul{\bf w}\|}
-\frac{(\ol{\bf w}^T{\bf R}_{yx}\ul{\bf w})\ul{\bf w}}
{\|\ol{\bf w}\|\|\ul{\bf w}\|^3}\right) = \nonumber
\\
&=& -\left[ \frac{\ul{\bf w}^T{\bf R}_{yx}^T{\bf R}_{yx}\ul{\bf w}}
{\|\ol{\bf w}\|\|\ul{\bf w}\|}-
\frac{\Bigl( \ol{\bf w}^T{\bf R}_{yx}\ul{\bf w}\Bigr)^2}
{\|\ol{\bf w}\|^3\|\ul{\bf w}\|} \right] \nonumber
\\
&& -\left[ \frac{\ol{\bf w}^T{\bf R}_{yx}{\bf R}_{yx}^T\ol{\bf w}}
{\|\ol{\bf w}\|\|\ul{\bf w}\|}-
\frac{\Bigl( \ol{\bf w}^T{\bf R}_{yx}\ul{\bf w}\Bigr)^2}
{\|\ol{\bf w}\|\|\ul{\bf w}\|^3} \right] = \nonumber
\\
&=& -\frac{1}{\|\ol{\bf w}\|\|\ul{\bf w}\|} \left[
\ul{\bf w}^T{\bf R}_{yx}^T \left( {\bf I}-\frac{\ol{\bf w}\ol{\bf w}^T}
{\|\ol{\bf w}\|^2} \right) {\bf R}_{yx} \ul{\bf w}
+ \ol{\bf w}^T{\bf R}_{yx} \left( {\bf I}-\frac{\ul{\bf w}\ul{\bf w}^T}
{\|\ul{\bf w}\|^2} \right) {\bf R}_{yx}^T \ol{\bf w} \right]
\end{eqnarray}

The matrices $\left( {\bf I}-\ol{\bf w}\ol{\bf w}/\|\ol{\bf w}\|^2 \right)$
and $\left( {\bf I}-\ul{\bf w}\ul{\bf w}/\|\ul{\bf w}\|^2 \right)$
are semipositive definite (they have one zero eigenvalue and the rest
are all equal to 1), thus
\begin{equation}
\frac{dV}{dt} \leq 0
\end{equation}
for all $t$. In addition, the only critical points of
Eqs.~\ref{eq:crosscor_ode_ol}, \ref{eq:crosscor_ode_ul} are
$\ol{\bf w}=\pm {\bf u}_i$ and $\ul{\bf w}=\pm {\bf v}_i$ $i=1,\ldots,q$.
However we will show that all cases where $i\neq 1$ are unstable
points.
Indeed the energy for any $i\neq 1$ is
\[
V(\pm{\bf u}_i,\pm{\bf v}_i) = \sigma_1 - \sigma_i > 0
\]
However, an infinitesimal perturbation 
\begin{equation}
\ol{\bf w}'=\pm {\bf u}_i+\epsilon {\bf u}_1,
~~~\ul{\bf w}'=\pm {\bf v}_i+\epsilon {\bf v}_1
\end{equation}
will decrease the energy to
\begin{equation}
V'=\sigma_1-\frac{(\pm {\bf u}_i+\epsilon {\bf u}_1)^T{\bf R}_{yx}
(\pm {\bf v}_i+\epsilon {\bf v}_1)}
{\|\pm {\bf u}_i+\epsilon {\bf u}_1\|\|\pm {\bf v}_i+\epsilon {\bf v}_1\|}
= \sigma_1-\frac{\sigma_i+\epsilon^2\sigma_1}{1+\epsilon^2}
= \frac{\sigma_1-\sigma_i}{1+\epsilon^2}
< V(\pm{\bf u}_i,\pm{\bf v}_i)
\end{equation}
therefore all points corresponding to $i>1$ are saddle points.
By similar perturbation analysis we can show that
the only (globally stable) attractor is $(\ol{\bf w},\ul{\bf w})=
(\pm {\bf u}_1,\pm {\bf v}_1)$.
The probability that the stochastic variables $(\ol{\bf w}_k,\ul{\bf w}_k)$
equal the nonstable critical points is zero.
Therefore, the stochastic
system of Eqs.~\ref{eq:crosscor_ol}, \ref{eq:crosscor_ul}
converges with probability 1 to the principal singular vectors
of ${\bf R}_{yx}$.
\hfill \qed


\subsubsection{Extracting Multiple Components}

In order to extend to the multiple component case
yet another deflation
transformation needs to be introduced.
This transformation is in the same spirit as all the
deflations introduced so far and constitutes
a generalization of the standard PCA deflation,
so that when ${\bf y}_k={\bf x}_k$ the two transformations
become the same.

\begin{definition}
The mapping
\begin{equation}
{\bf R}_{yx} \rightarrow
\tilde{\bf R}_{yx} \equiv {\bf R}_{yx} ({\bf I} - {\bf v}_1 {\bf v}_1^T)
= ({\bf I} - {\bf u}_1 {\bf u}_1^T) {\bf R}_{yx}
\end{equation}
is called a {\bf deflation transformation} of the matrix ${\bf R}_{yx}$.
\end{definition}

As with all kinds of deflations that we have examined
so far one can show that
\[
\tilde{\bf R}_{yx}={\bf U} {\bf \hat{\Sigma}} {\bf V}^T
\]
where ${\bf \hat{\Sigma}}=diag[0~\sigma_2~\sigma_3~\ldots~\sigma_q]$,
i.e. the second component becomes dominant.
The deflation can be achieved by a transformation on the data.
More specifically, one may either transform ${\bf y}_k$:
\begin{equation}
\tilde{\bf y}_k = {\bf y}_k - {\bf u}_1 {\bf u}_1^T {\bf y}_k
\end{equation}
or ${\bf x}_k$:
\begin{equation}
\tilde{\bf x}_k = {\bf x}_k - {\bf v}_1 {\bf v}_1^T {\bf x}_k
\end{equation}
therefore one can extract the second component (and by induction
all components) using Eqs.~\ref{eq:crosscor_ol}, \ref{eq:crosscor_ul}
on either $(\tilde{\bf x}_k,{\bf y}_k)$, $({\bf x}_k,\tilde{\bf y}_k)$
or $(\tilde{\bf x}_k,\tilde{\bf y}_k)$.
For example, the generalization for the
$p$-th component, assuming that the previous $p\!-\!1$ have
been already extracted, is
\begin{eqnarray}
\ol{\bf w}_{p,k+1} &=& \ol{\bf w}_{pk} + \beta ~\left[ {\bf y}_k
- \ol{\bf w}_{pk} b_{pk} \right] ~a_{pk}'
\label{eq:crosscor_lon_ol}
\\
\ul{\bf w}_{p,k+1} &=& \ul{\bf w}_{pk} + \beta ~\left[ {\bf x}_k
- \ul{\bf w}_{pk} a_{pk} \right] ~b_{pk}'
\label{eq:crosscor_lon_ul} 
\end{eqnarray}

where $a_k'$ and $b_k'$ are defined as in Eqs.~\ref{eq:a_prime}
\ref{eq:b_prime}. For convenience we repeat here
\begin{eqnarray}
a'_{pk} = a_{pk} - \sum_{i<p} \ul{c}_{pi} a_{ik},
& ~~ &
a_{ik} = \ul{\bf w}_i^T {\bf x}_k, ~i<p,
~a_{pk} = \ul{\bf w}_{pk}^T {\bf x}_k
\\
b'_{pk} = b_{pk} - \sum_{i<p} \ol{c}_{pi} b_{ik}
& ~~ &
b_{ik} = \ol{\bf w}_i^T {\bf y}_k, ~i<p,
~b_{pk} = \ol{\bf w}_{pk}^T {\bf y}_k
\end{eqnarray}
The lateral connection strengths $\ol{c}_{pi}$ and $\ul{c}_{pi}$,
should be equal to
\begin{equation}
\ol{c}_{pi}^{*} = \ol{\bf w}_i^T \ol{\bf w}_{pk}
~~~\mbox{and}~~~
\ul{c}_{pi}^{*} = \ul{\bf w}_i^T \ul{\bf w}_{pk}
\end{equation}
Under these definitions the corresponding differential equations
for the system become
\begin{eqnarray}
\frac{d\ol{\bf w}(t)}{dt} &=& {\bf R}_{yx}
({\bf I}-{\bf v}_1{\bf v}_1^T) \ul{\bf w}(t)
- \ol{\bf w}(t) \left( \ol{\bf w}(t)^T {\bf R}_{yx}({\bf I}-{\bf v}_1{\bf v}_1^T)
\ul{\bf w}(t) \right)
\\
\frac{d\ul{\bf w}(t)}{dt} &=& {\bf R}_{yx}^T
({\bf I}-{\bf u}_1{\bf u}_1^T) \ol{\bf w}(t)
- \ul{\bf w}(t) \left( \ol{\bf w}(t)^T ({\bf I}-{\bf u}_1{\bf u}_1^T)
{\bf R}_{yx} \ul{\bf w}(t) \right)
\end{eqnarray}
and since both matrices ${\bf R}_{yx}({\bf I}-{\bf v}_1{\bf v}_1^T)$
and $({\bf I}-{\bf u}_1{\bf u}_1^T){\bf R}_{yx}$ are equal to the
deflation matrix $\tilde{\bf  R}_{yx}$ according to Theorem
\ref{th:APCA_cross_correlation} the $p$-th component is going to be extracted.

As before one can also employ the Lateral Orthogonalization Rule for
estimating $\ol{\bf c}$ and $\ul{\bf c}$
(Fig.~\ref{fig:crosscor}(b)). The immediate orthogonalization
rule is applicable in this case: premultiplying Eq.~\ref{eq:crosscor_lon_ol}
by $\ol{\bf w}_i^T$ and Eq.~\ref{eq:crosscor_lon_ul} by
$\ul{\bf w}_i^T$ we obtain
\begin{eqnarray}
\ol{\bf c}_{pi,k+1} &=& \ol{\bf c}_{pi,k} + \beta
[b_{ik}-\ol{\bf c}_{pi,k}b_{pk}] a_{pk}'
\\
\ul{\bf c}_{pi,k+1} &=& \ul{\bf c}_{pi,k} + \beta
[a_{ik}-\ul{\bf c}_{pi,k}a_{pk}] b_{pk}'
\end{eqnarray}

\subsection{Simulations}

Figure \ref{fig:cross-sim} shows the convergence of the
network, using 200 random pairs $({\bf x}_k,{\bf y}_k)$ processed
in sweeps. The input dimension is 10 and the output dimension 8.
while the step-size parameter $\beta$ is 0.01 which is derived
empirically. We see that the weights converge to the
corresponding normal singular vectors as theoretically predicted.
Figure~\ref{fig:cross-sim-lyapunov} depicts the behavior of the Lyapunov
function $V$, in Eq.~\ref{eq:cross_Lyapunov}. Since the system
only approximates the differential equations
\ref{eq:crosscor_ode_ol} and \ref{eq:crosscor_ode_ul}
we don't expect the Lyapunov energy to decrease monotonically,
however we clearly see the monotonic trend {\em on the average}
which makes the energy eventually vanish.


\begin{figure}
\begin{center}
\begin{tabular}{c}
\psfig{figure=cross_up.ps,width=4in}
\\
(a)
\\
\psfig{figure=cross_lo.ps,width=4in}
\\
(b)
\end{tabular}
\end{center}
\caption{Simulation results for the cross-correlation network.
(a) The component estimation error
$\protect\|{\bf u}\protect_p-\protect\ol{\bf w}\protect_{pk}\protect\|^2$
for the upper layer. (b) The error
$\protect\|{\bf v}\protect_p-\protect\ul{\bf w}\protect_{pk}\protect\|^2$
for the lower layer.}
\label{fig:cross-sim}
\end{figure}

\begin{figure}
\centerline{\psfig{figure=cross_lyapunov.ps,width=4in}}
\caption{Cross-correlation network simulation:
the evolution of the Lyapunov energy $V$.}
\label{fig:cross-sim-lyapunov}
\end{figure}
