%
\title{Fast Learning in Networks of Locally-Tuned Processing Units}
%
\author{John Moody\\
Christian Darken\\
{\small \sl Yale Computer Science, Post Office Box 2158,
New Haven, CT 06520}}
%
\markboth{John Moody and Christian Darken}
{Fast Learning in Networks of Locally-Tuned Processing Units}
%
\maketitle{moody}
%
\begin{abstract}
We propose a network architecture which uses a single internal layer of
locally tuned processing units to learn to 
approximate real-valued functions (Moody and Darken 1988).
We formulate an efficient learning method which combines 
self-organized and supervised learning.
For a test problem, predicting a chaotic timeseries,
the network learns of order 1000 times faster in digital simulation time than a
three layer perceptron trained with back propagation, but 
requires of order ten times more training
data to achieve comparable prediction accuracy.
Unlike other efficient methods for timeseries prediction, our
network architecture and learning rules are appropriate for real-time
adaptive signal processing and adaptive control.
This is because they can be formulated to learn from examples and do not,
in principle, require storage of past data.
\end{abstract}

\vskip 1pc

Neurons with locally tuned response characteristics
are found in many parts of nervous systems. 
For example, the stereocilia cells in the cochlea have locally-tuned
frequency response.
The cells in somatosensory cortex respond selectively to stimulation from
localized regions of the body surface. 
The orientation-selective cells in visual cortex 
respond selectively to stimulation which is both local in retinal
position and local in object orientation. 
In this letter, we present a simple, idealized network model based upon
processing units which have locally tuned response functions.
Furthermore, we present an efficient method for training such a network.

A network of $M$ locally tuned units (Moody and Darken 1988) 
(see figure 1a) has an
overall response function:
\begin{eqnarray}
f(\vec{x}) &= \sum_{\alpha=1}^{M} f^\alpha R^\alpha (\vec{x}) 
\; \; , 
\label{eq:net} \\ 
R^\alpha (\vec{x}) &\equiv R(\| \vec{x} - \vec{x}^\alpha \| / \sigma^\alpha).
\label{eq:recept}
\end{eqnarray}

\begin{figure}\label{fig:rnet}
\vspace{2in}
\caption{
A network of locally-tuned processing units (A) and
a three layer perceptron (B). Unit response functions are depicted graphically.
Both networks have linear inputs and outputs. The network of locally tuned
units has a single internal layer, while the perceptron has 
two internal layers of sigmoidal units.
}
\end{figure}

Here, $\vec{x}$ is a real-valued vector in the input space,
$R^\alpha$ is the $\alpha^{{\rm th}}$ unit's response function,
$R$ is a radially symmetric function 
with  a single maximum at the origin and
which drops off rapidly to zero at large radii,
$\vec{x}^\alpha$ and $\sigma^\alpha$
are the center and width of the $\alpha^{{\rm th}}$ unit,
and $f^\alpha$ is the weight or amplitude associated with each unit.
For this work, we have chosen gaussian response 
functions with unit normalization:
\begin{equation}
R^\alpha (\vec{x}) = \exp \left[
- \frac{\| \vec{x} - \vec{x}^\alpha\|^2 } 
{(\sigma^\alpha )^2} \right].
\label{eq:gauss}
\end{equation}
One can think of the functional representation (\ref{eq:net}) as a 
decomposition into the functions which are not
orthonormal, not uniformly distributed over the input space,
and do not have uniform width.
With the addition of lateral connections between the processing
units, the network can dynamically
produce the normalized response function (Moody 1988):
\begin{equation}
f(\vec{x}) = \frac{\sum_\alpha f^\alpha R^\alpha (\vec{x}) }
{\sum_\alpha R^\alpha (\vec{x})}. 
\label{eq:norm}
\end{equation}
This is essentially a weighted average or interpolation between weights
or learned function values
$f^\alpha$ of nearby processing units; it is the functional form
we use in our simulations.

The locality of the unit response functions is important for attaining
fast simulation speeds on serial machines. For any given input, only
the small fraction of processing units with centers very close (in the
input space) to the input
will respond with activations which differ significantly from zero.
Thus, only those with centers
close enough to the input need to be evaluated and trained.
We efficiently identify these units by partitioning the input space with an 
adaptive grid (Omohundro 1987) and doing a short search.

In contrast to our network of localized processing units,
which requires only one internal layer,
the perceptron architecture (figure 1b) is believed
to require three layers (two internal layers) of threshold units to
approximate arbitrary functions on ${\cal R}^n$ ($n \geq 2$) with
reasonable accuracy. While this necessity has not been proven, constructive
arguments have been presented which show that three layers of threshold
units are {\it sufficient\/} to approximate arbitrary classification 
(Lippmann 1987) and real-valued mapping (Lapedes and Farber 1987a) problems.
The constructions are based on using two internal layers to isolate localized
convex regions or ``bump" representations of the input space.

We formulate learning in our network as a two stage
hybrid process. The processing unit centers and widths are determined in
a bottom-up or self-organizing manner, while the amplitudes
are found in a top-down manner using the supervised LMS rule.
The bottom-up component serves to allocate
network resources in a meaningful way by placing the unit centers
in only those regions of the input space where data is present.
It also dramatically reduces the amount
of computational time required by supervised learning, since only the
output weights (unit amplitudes)
must be calculated using an error signal.

Although our overall network response function is nonlinear,
each of our learning rules is based on minimizing
a purely quadratic objective function. 
This reduces the learning procedure to a set of simple linear update rules.
Our approach should be contrasted to the back propagation method, 
which formulates learning as a difficult, highly nonlinear 
optimization problem.
Our use of such simple linear learning
procedures offers another significant boost in speed over back propagation. 

The venerable
$k$-means clustering algorithm is used
to find a set of $k$ processing unit centers which
represent a local minimum of the total squared euclidean distances $E$
between the $k$ centers $\vec{x}_\alpha$ and $N$ exemplars 
(training vectors) $\vec{x}_i$:
\begin{equation}
E = \sum_{\alpha =1,k} \sum_{i =1,N} M_{\alpha i} 
\left( \vec{x}_\alpha - \vec{x}_i \right)^2 .
\label{eq:kmeans}
\end{equation}
Here, $M_{\alpha i}$ is the cluster membership function, which is a
$k \times N$ matrix of $0$'s and $1$'s with exactly one $1$ per column
which identifies to which processing unit a given exemplar belongs.
The (local) minimization of $E$ is usually formulated as an iterative
process using fixed exemplars, although 
a real-time, adaptive formulation of $k$-means 
can be easily formulated (Moody and Darken 1988).

The processing unit widths are
determined using various ``$P$ nearest neighbor" heuristics. 
These heuristics vary the widths in order to
achieve a certain amount of response overlap between each unit
and its neighbors so that they
form a smooth and contiguous interpolation over
those regions of the input space which they represent.
They can be formulated in an adaptive fashion
and can be shown to minimize various quadratic objective functions 
with respect to the $\sigma^\alpha$'s or $(\sigma^\alpha)^2$'s
and are thereby stable.
We have found that the particularly simple ``global first nearest
neighbor" heuristic
gives near-optimal results for some problems.
It uses a uniform average
width $\sigma= \langle \Delta x_{\alpha \beta} \rangle $ for all units
where $\Delta x_{\alpha \beta}$ is the euclidean distance in the input
space between each unit $\alpha$ and its nearest neighbor $\beta$ and
$\langle \rangle$ indicates a global average over all such pairs.

The response function amplitudes 
(output weights $f^\alpha$ in  equations \ref{eq:net} or \ref{eq:norm})
are varied to minimize the total error:
\begin{equation}
E = \frac{1}{2} \sum_{i = 1,N} 
\left( f^*(\vec{x}_i) - f(\vec{x}_i) \right)^2 \; \; ,
\end{equation}
where $\vec{x}_i$ is the $i^{{\rm th}}$ training pattern, $f$ is the total
network output, and $f^*$ is the desired result.
This is the supervised Adaline or LMS learning rule.
We have found that convergence occurs 
very quickly using a quasi-Newton method. This is 
possible because the self-organized learning which precedes the supervised
learning has already done most of the work. An adaptive, real-time version
of LMS is well known.

As a test problem, we consider the well-studied prediction of chaotic
timeseries. As it is usually formulated, this problem requires finding
a real-valued mapping $f: {\cal R}^n \mapsto {\cal R}$ which takes
a sequence of $n$ recent samples of a timeseries and predicts
the value of the timeseries at a future moment. 
It is assumed that the underlying process which generates the timeseries
is unknown.
We  shall compare our
network's learning and generalizing capabilities to a three layer
perceptron studied by Lapedes and Farber. The particular timeseries
we use for comparison results from integrating the Mackey-Glass differential
delay equation:
\begin{equation}
\frac{dx[t]}{dt} = -b\, x[t] + a \, \frac{x[t-\tau]}{1+x[t-\tau]^{10}}.
\end{equation}
Figure~\ref{fig:tseries}
shows the resulting timeseries for $\tau =17$, $a=0.2$, and $b=0.1$; 
note that it is quasi-periodic since no two cycles are the same.  
The characteristic time 
of the series, given by the inverse of the mean of the power 
spectrum, is $t_{char} \approx 50$. 
Note that classical techniques like global linear autoregression 
or Gabor-Volterra-Wiener polynomial 
expansions typically do no better than chance at predicting such a 
timeseries beyond $t_{char}$.

\begin{figure}
\vspace{2in}
\caption{
\label{fig:tseries}
Three hundred successive integer timesteps for the Mackey-Glass chaotic
timeseries with delay parameter $\tau = 17$.}
\end{figure}


For our numerical comparison,
both networks (fig. 1) have four real-valued inputs 
$(x[t]$, $x[t-\Delta]$, $x[t-2\Delta]$, $x[t-3\Delta])$ 
and one real-valued output $x[t+T]$
with $\Delta =6$ and $T=85 > t_{char}$. 
The network of locally-tuned units (fig. 1a) has 
between $100$ and $10,000$ internal units arranged in a single layer,
while the back propagation network (Fig.\ 1b) has two internal
layers each containing $20$ sigmoidal units.
The back propagation network thus has $541$ adjustable parameters
(weights and thresholds) total.

Figure~\ref{fig:num1} 
contrasts the prediction accuracy $E$ (Normalized Prediction Error) 
versus number of internal units for
three versions of our algorithm to the back propagation
benchmark (A) of Lapedes and Farber (Lapedes and Farber 1987b).
The three versions of the learning algorithm are:
\begin{enumerate}
\item
Nearest neighbor prediction. Here, the nearest data point in the
training set is used as a predictor. This behavior is actually a special
case for a network of form in equation \ref{eq:norm} where each input/output
training pair $\{\vec{x}_i,f_i\}$
defines a processing unit $\{\vec{x}^\alpha,f^\alpha\}$ of uniform narrow
width $(\sigma^\alpha \rightarrow 0)$.
\item
Adaptive processing units with one unit per data point.
Here, the amplitudes are determined by LMS, the widths by the global
first nearest neighbors heuristic, and the centers are fixed to be
training data vectors. 
\item
Self-organizing, adaptive processing units.
Similar to 2, 
but the training set has ten times more exemplars
than the network has processing units 
and the processing unit centers are found using
$k$-means clustering. 
\end{enumerate}
The back propagation benchmark
used a training set with $500$ exemplars.
For all methods, prediction accuracies were measured on a 500-member test
set. 

\begin{figure}
\vspace{3in}
\caption{
\label{fig:num1}
Comparison of Prediction Accuracy {\it vs.} Number of Internal Units
for four methods:
First Nearest Neighbor (1),
Adaptive Units (one unit per training vector) (2), 
Self-Organizing Units (ten training vectors per processing unit) (3),
and Back Propagation (A). The methods are described in the text. 
For (A), the abscissa indicates the number of training vectors.
The horizontal line associated with (A) is provided for visual reference
and is not intended to suggest a scaling law. In fact, the scaling law is
not known.}
\end{figure}


Note that method 1 requires storage of past timeseries data and method 2
assigns a processing unit to each data point. Hence, neither of these
methods is appropriate for real-time signal processing with fixed memory.
However, method 3 is fully adaptive in that it can vary a fixed set of
network parameters in response to new data in real time and does not,
in principle, require storage of previous data. 
In order to make a fair comparison to the back propagation benchmark, however,
we optimized our networks on a fixed training set, rather than measure
time-averaged, real-time performance.

As is apparent from figure \ref{fig:num1} (note horizontal reference line), 
method 2 achieves a 
prediction accuracy equivalent to back propagation with about $7$ times as much
training data, while method 3 requires about $27$ times as much data
to reach equivalent accuracy. 
These differences in training data requirements stem from the fact
that a back propagation network learns a global fit to the
function and is thus a more powerful generalizer, while the network of
locally tuned units learns only local representations.

Figure~\ref{fig:stime}
correspondingly shows $E$ versus computational time measured in
Sun 3/60 seconds for the 1000 internal units case.
The Lapedes and Farber benchmark
required about 1 hour of Cray X/MP time at an impressive
90 MFlops (Lapedes and Farber 1987b).
Their implementation used the conjugate gradient  minimization technique
rather than the more traditional gradient descent. Based upon our experience,
conjugate gradient provides convergence rates which are 10 to 100 times
faster than gradient descent.
Our simulations
on the Sun 3/60 probably achieved no more than 90 KFlops 
(the LINPACK benchmark), 
a factor of 1000 slower. 
For display purposes, we converted the Cray time 
into peak Sun 3/60 time by multiplying by 1000. 
This is somewhat conservative as 
we made no attempt to fine-tune our implementation on the Sun,
so our simulations were probably 3 to 4
times slower than peak Sun performance would allow.
 From Fig.~\ref{fig:stime}, 
it is clear that the self-organizing network (3)
achieved comparable
precision about {\bf 2000 times faster} than the back propagation network
implemented with conjugate gradient. 
The relative speed would probably be a factor of at least
20,000 if compared against back propagation implemented with
gradient descent.

\begin{figure}
\vspace{2.5in}
\caption{
\label{fig:stime}
Comparison of Prediction Accuracy {\it vs.} Simulation Time
for four methods:
Nearest-Neighbor (1), 
Adaptive Units (one unit per training vector) (2), 
Self-Organizing Units (3),
and Back Propagation (A). The methods are described in the text.}
\end{figure}


Several other methods and results which use techniques similar to
those described here are worthy of note. Various methods for doing
fast classification using local representations 
have been discussed (Huang and Lippmann 1988; Reilly, Scofield, Elbaum, et al.\
1987). Also, several efficient algorithms for doing multidimensional function
interpolation have recently been studied. These include local linear and
quadratic fitting (Farmer and Sidorowich 1988)
and radial basis functions (Casdagli 1988). These algorithms 
utilize local representations in the input space. However, 
they are not appropriate
for real-time use since they require explicit storage of large amounts
of past data. This is in contrast to the neural net approach in which
the only storage requirements are for the fixed set of
tunable network parameters. 

\section*{Acknowledgments}

We gratefully acknowledge helpful comments from and discussions with
Walter Heiligenberg, Alan Lapedes, Richard Lippmann, 
Bartlett Mel, Terry Sejnowski, John Shewchuk, and John Sidorowich. 
This research was supported by ONR grant 
N00014-86-K-0310, AFOSR grant F49620-88-C0025, and a Purdue Army subcontract.

\section*{References}
\refitem
Casdagli, M. 1988. {\it Nonlinear Prediction of Chaotic Time Series}.
Technical Report, London: Queen Mary College.

\refitem
Farmer, J. D.\ and J. J. Sidorowich. 1988. {\it Exploiting Chaos to Predict the
Future and Reduce Noise}. Technical Report, Los Alamos, NM: Los Alamos National
Laboratory.

\refitem
Huang, W.~Y.\ and R.~P. Lippmann. 1988. Neural Net and Traditional Classifiers.
{\it In:} Neural Information Processing
  Systems, ed.\ D.~Z. Anderson. American Institute of Physics.

\refitem
Lapedes, A.\ and R. Farber. 1987. How Neural Nets Work.
{\it In:}  Neural Information Processing
  Systems, ed.\ D.~Z. Anderson. American Institute of Physics.

\refitem
Lapedes, A.~S.\ and R. Farber. 1987. {\it Nonlinear Signal Processing Using
Neural Networks: Prediction
  and System Modeling}.
Technical Report, Los Alamos, NM: Los Alamos National Laboratory.

\refitem
Lippmann, R.~P. 1987. An Introduction to Computing with Neural Nets. {\it
IEEE ASSP Magazine\/} {\bf 4:4}.

\refitem
Moody, J. 1988. In preparation.

\refitem
Moody, J.\ and C. Darken. Learning with Localized Receptive Fields.
{\it In:\/} Proceedings of the
  1988 Connectionist Models Summer School, eds.\ Touretzky, Hinton, and
Sejnowski, Morgan Kaufmann.

\refitem
Omohundro, S. 1987. Efficient Algorithms with Neural Network Behavior. {\it
Complex Systems} {\bf 1}, 273.

\refitem
Reilly, D.~L., C.~Scofield, C.~Elbaum, and L.N.~Cooper. 1987.
Learning System Architectures Composed of Multiple Learning Modules.
{\it In:\/} IEEE First International
  Conference on Neural Networks, eds.\ M. Caudill and C. Butler, II--495.
\endpaper{moody}




