\documentstyle{article}

\begin{document}
\begin{center}
{\bf\LARGE\bf Contrastive Hebbian Learning in Interactive Networks}
\bigskip
\\Javier R. Movellan


 Department of Psychology
\\Carnegie Mellon  University
\\Pittsburgh, Pa 15213
 \\email: jm2z+@andrew.cmu.edu
\\ April 1990
\\ Submitted to Neural Computation.
\end{center}

Interactive networks  \cite{Grossberg78} \cite{Hopfield84} \cite{IAC} characterized
by the existence of bidirectional connections, may  have an
advantage over feed-forward architectures because of  their completion properties, and
flexibility in the treatment of units as inputs or outputs. Ackley, Hinton and
Sejnowski \cite{Ackley85} derived a learning rule to train Boltzmann machines, which
are discrete, interactive networks.  Unfortunately, because of the discrete
stochasticity of its units,  Boltzmann learning is intolerably slow. 

Peterson and Anderson \cite{Peterson87} showed that
Boltzmann machines with large number of units can be approximated with deterministic
networks whose logistic  activations represent the average  activation of 
discrete Boltzmann units (Mean Field Approximation). Under these conditions a
learning rule that I will call Contrastive Hebbian Learning (CHL) was shown to be a
good approximation to the Boltzmann weight update rule and to achieve learning speeds
comparable to backpropagation. This modification of the Hebbian learning rule was first
applied by Hopfield \cite{Hopfield83} to improve the storage capacity of content addressable
memories \cite{Hopfield82}. Hinton \cite{Hinton89} showed that for Mean Field networks, CHL is at
least a first order approximation to gradient descent on an error function. The purpose in this
paper is to show that CHL  works with any interactive network with bounded, continuous activation
functions and symmetric weights. The approach taken does not  presume the existence of Boltzmann
machines whose behavior is approximated with mean field networks. It is  also shown that CHL
performs gradient descent on a contrastive function of the same form investigated by  Hinton
\cite{Hinton89}. 

The paper is divided in two sections and one appendix. Section 1 studies the
dynamics of the activations in interactive networks. Section 2 shows how to
modify the weights for the stable states of the network to reproduce desired
patterns of activations. The appendix contains
mathematical details, and some comments on how to implement Contrastive Hebbian
Learning in Interactive Networks.

 

\section{ Stability of Activations}
Since interactive networks have recurrent paths, it is in principle possible that
their activations never stabilize. Fortunately, if some
simple conditions are met, we can guarantee that the activations will settle,
and that at equilibrium they will be at a minimum of an $Energy$ function
\cite{Hopfield84}.    \par
Let the activation vector ${\bf a}^T = [ a_1, ... a_n]$, be regulated by  bounded,
monotonically increasing, differentiable activation functions $f_{i}(.)$. Let ${\bf
W} = [ {\bf w}_1, ... {\bf w}_n ] $ be the matrix of connections, where the $
{\bf{w}}_{i}^{T} = [w_{1,i}, ... w_{n,i}]$ are bounded, fan-in weight column vectors.
Let $d ()/dt$ be  derivatives with respect to time, $ net_i = {\bf a}^{T}_{i} \; {\bf
w}_i$, and  $rest = f(0)$. Define a continuous Hopfield Energy function
\cite{Hopfield84}.

\begin{equation}
F = E + S 
\end{equation}
where 

\begin{equation}
E= - \frac{1}{2}\; {\bf a}^T {\bf W} \:{\bf a} =
-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n} a_i w_{ij} a_j \end{equation}
and

\begin{equation}\label{S}
S = \sum_{i=1}^{n} \;\int_{rest}^{a_i} f^{-1}_{i} (a) \;da 
\end{equation}

\bigskip

Notice that E, which reflects soft constraint satisfaction, tends
to drive activations to extreme values, and S is a penalty
function for departing from the  rest activation value. 
In the appendix I show that if the activation functions are the standard (0-1)
logistic, F becomes the Helmholtz free energy function as defined in \cite{Hinton89}.
\bigskip

Hopfield \cite{Hopfield84} showed that if the network is governed by the
differential equation

\begin{equation}\label{snet}
\frac{d\:f^{-1}_{i}(a_i)}{d t} = \lambda\; ( -f^{-1}_{i} (a_i) + net_i) 
\end{equation}
\smallskip
and the weights are symmetric, the activations stabilize in a minimum of
F. For completeness, I present in the appendix a version of Hopfield's proof and
show that stability in a global minimum can also be achieved with the following
equation, typically used in interactive activation networks \cite{Grossberg78}\cite{IAC} \cite{McClelland89}

\begin{equation}\label{sact}
\frac{d\:a_i}{d t} = \lambda\; \left(\left( - a_i + f_i(net_i\right)\right) 
\end{equation}
Notice that if we apply either equation \ref{snet} or \ref{sact}, on equilibrium (when
the derivatives are zero),
\begin{equation}\label{equi}
f^{-1}_{k} (\breve a_i) =\breve {net_i}
\end{equation}
where $(\;\breve{ }\;)$ represents equilibrium.
I will use these properties to derive the learning algorithm in the next
section.

\section{Contrastive Learning}

Learning  can be viewed as the modification of connections between units so that the
stable states of the network reproduce desired patterns of activations.  We will
see that CHL minimizes a {\it contrastive function}  $J$ \footnote{Based on
different arguments, Hinton \cite{Hinton89} showed that Mean Field Learning minimizes
the equivalent of J when the activation is logistic} and then we will discuss the
conditions under which minimization of $J$ produces learning. 

Define a pattern $p$ as the pair $ p = \{ {\bf a}^{I(+)}, {\bf a}^{O(+)} \}$, where
$I$ stands for input set, O for output set, and (+) indicates that the activation of
these units are fixed by the  pattern. We will say that $p$ has been learned if
clamping the input units to ${\bf a}^{I(+)}$ then \begin{equation}
 {\bf \breve {a}}^{O(-)} = {\bf a}^{O(+)}
\end{equation}
where $(\;\breve{ }\;)$ denotes equilibrium,  and
${\bf\breve a}^{O(-)}$  the activation of the output set when inputs are clamped
and outputs are free. 

Define the {\it contrastive function} $J$ as 

\begin{equation} J = \breve F^{(+)} - \breve F^{(-)}  \end{equation} where $\breve  F^{(+)}$
and $\breve F^{(-)} $ respectively are the values of the energy functions at equilibrium
when the inputs and outputs are clamped {\small (+)}, and when the inputs are clamped and
outputs are free {\small (-)}. Notice that $\breve F^{(-)} $ has the same free parameters
that $F^{(+)}$, the activations of the hidden units, plus some additional free parameters,
the activations of the output units. Assume that, over the working  region of activation
states,  $F$ has a unique minimum.\footnote{  This assumption, which is shared with 
Boltzmann learning, Mean Field Learning, and the Almeida-Pineda algorithm is discussed in
section \ref{J}}. Thus, since the minimum is unique   $\breve  F^{(+)} \geq \breve
F^{(-)}$, and if $J=0$, then $ {\bf \breve a}^{O(-)} = {\bf a}^{O(+)}$. This makes $J$
a potential candidate for learning by gradient descent on weight space.  In
particular, following the derivations detailed in the appendix, we have 

\begin{eqnarray}
\frac{\partial \breve E}{\partial w_{ij}} &= &-   \breve a_i \breve a_j  -
\sum\limits_{k=1}^{n} \breve {net_k}(\frac{\partial \breve a_k}{\partial w_{ij}}) \;
;\;\;   i \neq j \label{wij}\\
\frac{\partial \breve E}{\partial w_{ij}} &= &- \frac{1}{2}  \breve a_i^2  -
\sum\limits_{k=1}^{n} \breve {net_k} (\frac{\partial \breve a_k}{\partial w_{ij}})
\;;\;\; i =j \label{wii}
\end{eqnarray}

\begin{equation}\label{S2}
\frac{\partial \breve S}{\partial w_{ij}} = \sum_{k=1}^{n} f^{-1}_{k} (\breve a_k)
\frac{\partial \breve a_k}{\partial w_{ij}}
\end{equation}

And following equation \ref{equi} it is easy to see that
\begin{eqnarray}
\frac{\partial \breve F}{\partial w_{ij}} &= &-   \breve a_i \breve a_j \;
;\;\;   i \neq j\\
\frac{\partial \breve F}{\partial w_{ij}} &= &- \frac{1}{2}  \breve a_i \breve a_j \;
;\;\;   i = j\\
\end{eqnarray}
making
\begin{equation}
\frac{\partial \breve J}{\partial w_{ij}} \propto \breve a_{i}^{(-)}  \breve
a_{j}^{(-)} -  \breve a_{i}^{(+)}
 \breve a_{j}^{(+)}
\end{equation}
which shows that the contrastive Hebbian learning rule
\begin{equation}
\Delta w_{ij} \propto  \breve a_{i}^{(+)}  \breve a_{j}^{(+)} -  \breve a_{i}^{(-)}
 \breve a_{j}^{(-)}
\end{equation}
descends in the J function.
\section{Discussion}\label{J}

We have seen that CHL minimizes the contrastive function 
\begin{equation}
J = \breve F^{(+)} - \breve F^{(-)} 
\end{equation}
so that after each learning step, the difference in energy at equilibrium between the
clamped and free states becomes smaller. A  way to visualize how CHL works
is to imagine the energy surface over activation space as a membrane with several
minima. CHL pushes down the minimum corresponding to the clamped case and
pulls up the stable states for the free phases. If both stable states are in the
same local attractor, after several learning trials, both minima converge.
Here we study how CHL works in three different cases:

\begin{itemize}
\item {\bf Case 1:} {\it Activations are reset to random numbers after each learning phase.} 
In this case, the starting points for the clamped and free phase are different and, very
likely the stable states will also be  apart. Under this condition, CHL learning does not work
well. 

\item {\bf Case 2:} {\it First settle for the clamped phase, and then, without resetting
activations, free the output units and settle again.}  
This procedure guarantees that  $F^{(+)} \geq \breve F^{(-)} $ and that when $F^{(+)} = \breve
F^{(-)} $ the activations on the free and clamped phases are the same. Gradient descent on $J$
assures that when the minima for the clamped phase is achieved, if the output unit activations are
free, they will not change. This form of learning, which can be used for recognition of familiar
patterns is very rapidly achieved with CHL (it just takes about 3 trials to "recognize" the XOR or
the 4-3-4 encoder patterns). Unfortunately, if we just clamp the input units, there is not
assurance that the activations will converge to the desired minima.

\item {\bf  Case 3:} {\it Settle during the free phase and then, without resetting
activations,  clamp the output units and resettle again.}

This scheme minimizes the probability that the clamped and free phases end up in different
attractors. In general this procedure works very well and achieves learning speeds comparable
to backpropagation. There are two phenomena though that sometimes occur. Once in
a while, the network may settle in a different attractor than the one to which it had
converged in previous trials. This may produce a sudden change in activations and
and a temporary "unlearning" of the weights. It is our experience that
as learning progresses these sudden jumps to other minima tend to diminish. Another problem is
that if the the clamped and free phases stabilize in different regions, the energy of the
clamped phase may become lower than the energy in the free phase. If this happens,
learning usually deteriorates making it advisable to start with a different set of weights. In
spite of these problems, CHL works very well and appears to be a promising, simple to implement
learning method that works for a wide variety of interactive architectures. In table 1 we show the
results of some simulations using CHL with the activation update rule of the Interactive
Activation and Competition model \cite{McClelland89}.  
\end{itemize}

\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|} \cline{2-5}
\multicolumn{1}{c|}{ }&\multicolumn{4}{c|}{\sl stepsize} \\ \hline
{\sl hidden units}   &1.0    &0.1    &0.01    &0.001 \\ \hline
1 &19 (0) &30 (0) &246 (5)  &$\geq$300 (10) \\ \hline
2 &13 (0) &15 (0) &72 (1)  &$\geq$300 (7) \\ \hline
4 &55 (3) &20.5 (0) &40.5 (0)   &$\geq$300 (6) \\ \hline
8 &47 (1) &23 (0)  &21 (0)  &245 (4) \\ \hline
16 &$\geq$300 (10) &98.5 (0) &24.5 (0)  &142(1) \\ \hline
\end{tabular}
\end{center}
\caption{ Results for the XOR problem. Networks were fully connected. There were 10 simulations
per cell with different starting random weights from a uniform (-0.5 to 0.5) distribution.  Outside
parenthesis are median number of epochs until maximum absolute error was smaller than 0.2. In
parenthesis are number of simulations in which learning took longer than 300 epochs. Learning was
on batch mode. Update of activation was done according to the IAC equations with the
following parameter values: max=1.0; min=-1.0; rest=0.0; decay=0.1; lambda= 0.2. Update
of activations was stoped after the change in all the activations was smaller than 0.0001 or the
number of cycles bigger than 100.}\end{table}



 \section{Acknowledgement} This paper would not have been possible without the support and
motivating ideas of James McClelland and his research group at Carnegie Mellon University. I also
thank Conrad Galland and Geoffrey Hinton for their advise in the implementational details of the
Deterministic Boltzmann Machine.

\section{Appendix}  
\subsection{S and Mean Field Entropy}
The entropy term of the energy function proposed by Hinton \cite{Hinton89} for the
Mean Field Algorithm is
\begin{equation}\label{entropy}
- \sum^{n}_{i=1} \left( a_{i} \:log a_{i} + \left( 1- a_{1} \right) log \left( 1 - a_i
\right) \right). 
\end{equation}
If we use the standard logistic (0-1) activation function, 
\begin{equation}
y = \frac{1}{1+e^{ -x/T}}
\end{equation}
whose inverse is
\begin{equation}
x = T\;log \left( \frac{y}{1 - y} \right)
\end{equation}
then
\begin{equation}
\sum^{n}_{i=1} \int_{0.5}^{y} T\; log \left( \frac{y}{1 - y} \right) \;dy = T \sum^{n}_{i=1}\left( y
\;log \:y + \left( 1- y \right) log \left( 1 - y \right) \right)
\end{equation}

\subsection{ S and the Interactive Activation and Competition Model}
The activation update rule of the interactive activation and competition model (IAC) as defined in
\cite{IAC} is 
\begin{equation} 
\Delta a_{i} =  \lambda (max - a_i) net_i - (a_i - rest) decay \;;\;net\geq 0 
\end{equation}
\begin{equation}
\Delta a_{i} =  \lambda ( a_i- min) net_i - (a_i - rest) decay\;;\;net \leq 0
\end{equation}
which can be derived from equation \ref{sact} applied to the following activation
function 
\begin{equation}\label{iac1}
a_i =  \frac{ max \: net_i + rest\: decay}{ net + decay}\;;\;net \geq 0 \\
\end{equation}
\begin{equation}\label{iac2}
a_i =  \frac{ min\: net_i - rest\: decay}{ net - decay}\;;\;net \leq 0 
\end{equation}
where max is the maximum value of the activation, rest  the activation when the
net input is zero, min  the minimum value of the activation, and decay  a
positive constant. And applying equation \ref{S}, it is easy to see that 
\begin{equation}
S = \sum\limits ^n_{i=1} S_i
\end{equation}
with 
\begin{equation}
S_i = decay \left( \left(  max - rest \right)\;log \left( \frac{max -
rest}{max - a_i} \right) - \left( a_i -rest\right) \right) ;\;\; a_i \geq rest 
\end{equation}
\begin{equation}
S_i = decay \left( \left( a_i - rest \right) + \left( min - rest \right)\;log 
\left( \frac{a_i - min}{rest - min} \right)\right) ;\;\; a_i \leq rest  
\end{equation}
with the decay parameter assuming the same function that temperature in
the logistic model.

 \subsection{Hopfield's proof of stability of activations}  
Hopfield
showed that if the network is regulated by equation \ref{snet} it will stabilize.
Here,  I use a similar argument to show that using  equation \ref{sact} the
network also stabilizes.   
To facilitate the calculation of the derivatives of
$E$, we collapse into $Q$ the part of $F$ that does not depend on $a_i$,
\begin{equation} 
E = -\frac{1}{2}\sum_{k=1}^{n}\sum_{l=1}^{n} a_k w_{kl} a_l = 
-\frac{1}{2} \left( a_i w_{ii} a_i + a_i \left( \sum_{{k=1}\atop{k \neq
i,j}}^{n} a_k w_{ik} 
 + \sum_{{l=1}\atop{l \neq
i,j}}^{n} a_l w_{li} \right) + Q \right) 
\end{equation} 
and considering that the weights are symmetric, 
\begin{equation} 
\frac{\partial
E}{\partial a_{i}} = -\frac{1}{2} \left( 2 a_i w_{ii} + 2\sum_{{k=1}\atop{k \neq
i,j}}^{n} a_k
 w_{ik} \right) = - net_i  
\end{equation} 
Regarding S, since the derivative of the
integral of a function is the function itself 
\begin{equation} 
\frac{\partial S}{\partial a_{i}} =  f^{-1}_{i} (a_i) 
\end{equation}
and
\begin{equation}
\frac{\partial F}{\partial a_{i}} = \frac{\partial E}{\partial
a_{i}} + \frac{\partial S}{\partial a_{i}} = - net_i + f^{-1}_{i} (a_i)
\end{equation}

If we make  
\begin{equation} 
\frac{d a_i}{d t} = \lambda\; ( - a_i + f(net_k)) 
\end{equation} 
then, applying the chain rule,  
\begin{equation}\label{timed}
 \frac{d F}{d t} =
\sum_{k=1}^{n} \frac{\partial F}{\partial a_{k}} \frac{d a_{k}}{d t} =
\sum_{k=1}^{n} \lambda \left(-net +f^{-1}_k (a_k)\right) \left( -a_k +
f(net_k)\right)  \end{equation}
but since $f_k$ is monotonic, $\left(-net +f^{-1}_k (a_k)\right)$ has the same sign
as $\left( a_k - f(net_k)\right)$, making \begin{equation} \frac{d F}{d t} \leq 0
\end{equation}

Since F is bounded  and on each time step F decreases then
\begin{equation}
\lim _{t \rightarrow \infty} \frac{d F}{d t} = 0
\end{equation}
and since equation \ref{timed} shows that 
\begin{equation}
\frac{d F}{d t} = 0 \Longleftrightarrow \frac{d a_{k}}{d t} = 0
\Longleftrightarrow \frac{\partial F}{\partial a_k} = 0 \;\;; k=1,...,n
\end{equation} then 
\begin{equation}
\lim _{t \rightarrow \infty} \frac{d a_i}{d t} = 0 \;\;; k=1,...,n 
\end{equation} 
which tells us that the activations will tend to equilibrium as time
progresses and that on equilibrium they are in a minima of F.



\subsection{Smoothing net inputs vs. smoothing  activations} 
Equation \ref{snet} can be
discretized as  
\begin{equation} 
\Delta f^{-1}_{i}(a_{i(t)}) = \lambda\; (
-f^{-1}_{i} (a_{i(t)}) + net_{i(t)}) 
\end{equation} 
or 
\begin{equation}\label{dsnet}
f^{-1}_{i}(a_{i(t+1)}) = (1 - \lambda) f^{-1}_{i}(a_{i(t)}) + \lambda net_{i(t)}
\end{equation} 
equivalently, equation \ref{sact} would be discretized as  
\begin{equation}\label{dsact}
a_{i(t+1)} = (1- \lambda) a_{i(t)} + \lambda f(net_{i(t)}) 
\end{equation}

Equation \ref{dsnet} is an exponential smoothing of the net input. The activation
function is then applied to this smoothed net. Another possibility, represented in
equation \ref{dsact}, is to apply the activation function to the instantaneous net
input and then to exponentially smooth these activations. 
\\
Peterson and Anderson \cite{Peterson87}
and Peterson and Hartman\cite{Peterson89} recommend the use of an "{\it annealing schedule}" in
analogy with the Boltzmann machine that the Mean Field network approximates. It is important to
clarify, that since Mean Field networks are deterministic, this procedure, which is better
described as a {\it sharpening schedule} \cite{Gaussian89} is not truly annealing, and no longer
guarantees convergence in a global minima. The effect of 
the sharpening schedules is to gradually increment the net inputs avoiding excessive changes
in the activations. This same
effect can also be achieved by making 
 $\lambda$ in  equation  \ref{dsact} less than unity . \footnote{The
relationship between annealing schedules and exponential smoothing of the activations was
proposed to me by James McClelland}.

\subsection{ The derivatives of the equilibrium points
of F} First consider the weights that connect different units. Extracting the cross
products with a $w_{ij}$ term. We have \begin{equation} \breve E = -\frac{1}{2}
\left( 2 \breve a_i w_{ij} \breve a_j + \sum\limits_{{k=1}\atop{k
\neq i,j}}^{n}\sum\limits_{{l=1}\atop{l \neq i,j}}^{n} \breve a_k w_{kl} \breve a_l
\right) \end {equation} and considering that $w_{ij}$ is the only weight depending on
$w_{ij}$ \begin{equation} \frac{\partial \breve E}{\partial w_{ij}} =  -\frac{1}{2}
\left( 2 \breve a_i \breve a_j +  2 w_{ij}\breve a_i \frac{\partial \breve
a_j}{\partial w_{ij}} + 2 w_{ij}\breve a_j \frac{\partial \breve a_i}{\partial
w_{ij}} + \sum\limits_{{k=1}\atop{k \neq i,j}}^{n} \sum\limits_{{l=1}\atop{l \neq
i,j}}^{n} w_{kl} \left( \breve
a_k \frac{\partial \breve a_l}{\partial w_{ij}} +\breve a_l \frac{\partial \breve
a_k}{\partial w_{ij}} \right) \right)  \end {equation} Reorganizing terms
considering that the weights are symmetric, \begin{equation} \frac{\partial \breve
E}{\partial w_{ij}} = -\frac{1}{2} \left( 2 \breve a_i \breve a_j + 2\sum\limits_{k=1
}^{n} \frac{\partial \breve a_k}{\partial w_{ij}} \sum\limits_{l=1 }^{n} w_{kl} 
\breve a_l \right) \end{equation}
which easily leads to equation \ref{wij}. Similar arguments can be applied to
derive equation \ref{wii}. Equation \ref{S2} is easy to obtain by applying the chain
rule and the fact that the derivative is the inverse of the integral.

\subsection{ Sketch of main routines in a contrastive learning program}

\begin{enumerate}
\item Get a training pattern.
 
\item Reset activations to rest and net inputs to zero.

\item Clamp inputs to desired pattern.

\item Settle activations according to equations \ref{iac1} and \ref{iac2} .

\item Collect cross products of activations multiplied by a negative constant.

\item Clamp also the  output units to the desired pattern.

\item Settle activations according to equations \ref{iac1} and \ref{iac2}.

\item Collect cross products of activations multiplied by a positive constant.
\end{enumerate}


Termination of settling may be done after a sufficient number of iterations, or
after the changes in activations are smaller than a certain criteria. Bellow is an example of a
settling cycle using the update function of the IAC model
\cite{IAC}. \bigskip

{\it for(i=1;$\;$i$\leq$ number\_of\_ units;$\;$i++)\{

$\;\;\;\;$for(j=1;$\;$j$ \leq$ number\_of\_units;$\;$j++)\{	
	 
$\;\;\;\;\;\;\;\;$net[i] = net[i] + w[j][i]*activation[j] ;

$\;\;\;\;$\}



$\;\;\;\;$if(net[i] $\geq$ 0) act[i] = acti[i] + lambda*((actimax -act[i]) -
decay*(acti[i]-rest));\\ $\;\;\;\;$else acti[i] =acti[i] + lambda*((acti[i] - actimin) -
decay*(acti[i]-rest)); \}\\
}\\
Where actimax is the maximum value of the activations, actimin the minimum value, rest the
activation when net is zero, lambda the stepsize of the activation change (smoothing constant),
and decay a positive constant. We have obtained good results with actimax =1.0, actimin = -1.0,
rest =0.0, decay =0.1, lambda = 0.2.

This is an example of a weight change routine  \\ {\it
for(i=1;$\;\leq$ number\_of\_units;$\;$i++)\

$\;\;\;\;$for(j=1;$\;\leq$ number\_of\_units;$\;$ j++)\{	

$\;\;\;\;\;\;\;\;$weight\_change[i][j] = weight\_change[i][j] + phase *epsilon*
a[i]*a[j] ;

$\;\;\;\;$\}

\}\\
}

Where {\it phase} is (+1) for the clamped case and (-1) for the free case. Epsilon
is the stepsize of the weight change. Weights may be updated after each phase,
each pattern or after a batch of patterns

\begin{thebibliography}{99} 

\bibitem{Ackley85} Ackley D, Hinton G and
Sejnowski T (1985) A learning algorithm for Boltzmann machines. {\it Cognitive
Science}, 9, 147-169.

\bibitem{Gaussian89} Akiyama Y, Yamashita A, Kajiura M, Aiso H (1989) Combinatorial Optimization
with Gaussian Machines, {\it Proceedings of the International Joint Conference on Neural Networks}
1, 533-540.

\bibitem{Hinton89} Hinton G E (1989) Deterministic Boltzmann Learning Performs
Steepest Descent in Weight-Space,  {\it Neural Computation}, {\bf1}, 143-150.

\bibitem{Hopfield82} Hopfield J (1982) Neural Networks and Physical Systems with emergent
collective computational abilities. {\it Proceedings of the National Academy of Science} USA, 79,
2254-2558.

\bibitem{Hopfield83} Hopfield J, Feinstein D, Palmer R (1983) Unlearning has a
stabilizing effect in collective memories. {\it Nature} 304, 158-159.

\bibitem{Hopfield84} Hopfield J  (1984) Neurons
with graded response have collective computational properties like those of two-state neurons. {\it
Proceedings of the National Academy of Sciences} U.S.A., {\bf 81}, 3088-3092.

\bibitem{Grossberg78} Grossberg S A (1978) A theory of visual coding, memory, and
development. in E Leeuwenberg and H Buffart (Eds.) {\it Formal Theories of visual
perception} New York, Willey.

\bibitem{IAC} McClelland J, Rumelhart D (1981) An Interactive Activation Model of
Context Effects in Letter Perception: Part 1. An account of Basic Findings. {\it
Psychological Review}, 88, 5.

\bibitem{McClelland89}Mclelland J, Rumelhart D (1989) Interactive activation and Competition.
Chapter II of {\it Explorations in Parallel Distributed Processing: A Handbook of Models,
Programs, and Exercises} Cambridge, MIT Press.

\bibitem{Peterson87} Peterson C, Anderson J R (1987) A
mean field theory learning algorithm for neural networks. {\it Complex Systems}, {\bf
1}, 995-1019.

\bibitem{Peterson89} Peterson C, Hartman E (1989) Explorations of the Mean Field Theory Learning
Algorithm. {\it Neural Networks}, 2, 475-494.



\end{thebibliography}

\end{document}

