\chapter{Neural Networks}
\begin{flushright}
\parbox{5.5in}{We tend to regard forgetting as a failure of memory, the result of an inadequacy in the machinery for remembering, but a moment's reflection shows that that cannot always be right. If we were to remember every telephone number we have ever used, every scene we have ever looked at, every conversation we have ever heard, every sentence we have ever read, we should be overwhelmed with memories. In fact we transfer to our long-term memory stores only a tiny fraction of the contents of our short-term stores. We retain only what is likely to be useful.~\cite[page~331]{glynn99}}
\end{flushright}

\bigskip

%This chapter is concerned with neurocomputing which is also called the brainlike computation.
The motivation behind studying artificial neural networks,\index{artificial neural networks|see{neural networks}} commonly referred to as neural networks,\index{neural networks} or just neural nets,\index{neural nets|see{neural networks}} is to build computers whose construction would mimic the organization of the brain. The brain is the most complex natural information processing system.\index{brain!natural information processing system} It is capable of organizing its building blocks, called neurons,\index{neurons!brain's building blocks} to perform computations in a very efficient way. A drawing of a typical biological neuron is given in Figure~\ref{fig:bio_neur}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{bio_neur.eps}
%\vspace{2.5in}
%\centerline{\psfig{figure=bio_neur.eps,height=2.5in}}
%\hskip.5in\special{eps:c:/mybook/bio_neur.eps y=2.5in}
\caption{A typical biological neuron.}
\label{fig:bio_neur}
\end{figure}
The neuron collects signals from other neurons through its dendrites.\index{dendrites} The signals have the form of electrical pulses. These pulses travel to the neuron's cell body\index{neuron's cell body} where they are being processed. The surface of the cell body of the neuron is dotted with numerous receptor sites that can receive messages from other neurons in a fashion similar to the dendrites. The neuron may produce its own pulses in response to excitation from other neurons. The neuron's own pulses travel through a long, thin strand called an axon,\index{axon} which splits into many end branches. The neuron's message to the other neurons is delivered via the end branches. The key to the transmission of the neuron's message is the synapse,\index{synapse} shown in Figure~\ref{fig:neur_syn}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{neur_syn.ps}
%\vspace{2.5in}
%\centerline{\psfig{figure=neur_syn.ps,height=2.5in}}
%\hskip1in\special{ps:c:/mybook/neur_syn.ps y=2.5in}
\caption{The synapse is a junction point where neuron meets neuron.}
\label{fig:neur_syn}
\end{figure}
The end branches of the neuron's axon terminate in synaptic knobs.\index{synaptic knobs} The synaptic knob of the first neuron is separated from the receptor site of another neuron by a microscopic distance. The cell body of the first neuron produces chemical substances called neurotransmitters,\index{neurotransmitters} which are delivered down to the synaptic vesicles. The neurotransmitters are stored in the synaptic vesicles\index{synaptic vesicles} until the neuron fires and a burst of neurotransmitters is released by the vesicles. The neurotransmitters then flow across the synaptic cleft and act on the second neuron. The neurotransmitters of the first neuron may stimulate or inhibit the second neuron activity.\index{neurotransmitters!stimulate or inhibit} The first neuron has many other synaptic knobs elsewhere, while the second neuron has many other receptor sites that receive messages from many other neurons. Each neuron is receiving messages from hundreds or thousands of other neurons. The neuron sets off the nervous pulse, that is, it fires off its own nervous pulse only when it is properly stimulated. It was observed that, in general, the neuron will not fire in response to a single message. To fire, the neuron needs many messages arriving simultaneously or in quick succession from other neurons. In addition, the stimulating messages to fire must outweigh the messages that inhibit the neuron from firing.

In the following sections, we discuss a few mathematical models of the biological neuron.\index{mathematical models!of biological neuron} These models are referred to as the artificial neurons.\index{artificial neurons} Then, we analyze neural networks that are composed of the artificial neurons, and finally we discuss applications of the artificial neural networks in the controller design of dynamical systems.

\section{Threshold Logic Unit}
\label{sec:tlu}
A simple model of a biological neuron is the so called the {\em threshold logic unit\/}~(TLU),\index{threshold logic unit (TLU)}\index{TLU|see{threshold logic unit}}\index{TLU|see{perceptron}} shown in Figure~\ref{fig:tlu_ner}. The TLU is also called the {\em perceptron\/}.\index{perceptron} A TLU is parameterized by a set of weights and a threshold (bias). The summation operation\index{summation operation} is followed by the sign activation function,\index{sign activation function} which models the cell body of a biological neuron. The dendrites are modeled by the inputs and their strengths are modeled by the weights. The axon\index{axon} corresponds to the output $y$.
\begin{figure}
\centering
\includegraphics[scale=0.5]{tlu_ner.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/tlu_ner.bmp}
%\hskip1in\special{bmp:c:/mybook/tlu_ner.bmp x=4in y=3in}
\caption{A threshold logic unit~(TLU), also called the perceptron, as a model of a biological neuron.}
\label{fig:tlu_ner}
\end{figure}
Suppose that we are given a TLU with $n$ weights, $w_1,w_2,\ldots,w_n$, and a threshold $t$. Given an input vector $\b x^{(i)} =\left[\begin{array}{cccc}
x_{1}^{(i)} & x_{2}^{(i)} & \cdots & x_{n}^{(i)} \end{array}\right]^{\top} \in \mb{\Bb R}^n$, then the corresponding output of the TLU is
\begin{equation}
y_i=\mb{sign}\left(\sum_{j=1}^n w_jx_j^{(i)} -t \right),
\label{eq:tluoutyi}
\end{equation}
where
\[
\mb{sign}(z)=\left\{\begin{array}{ccc}
-1&\mb{if}& z < 0\\
1&\mb{if}& z\ge0.
\end{array}\right.
\]
We note that sometimes in the literature the sign function\index{sign function} is undefined at 0, or it is defined as a point to a set mapping, that is, $\mb{sign}(0)\in [-1,1]$. Here, we define, $\mb{sign}(0)=1$. Using the scalar product\index{scalar product} notation, we represent~(\ref{eq:tluoutyi}) as
\[
y_i=\mb{sign}\left( \langle \b w,\b x^{(i)}\rangle -t\right).
\]
Note that $y_i \in \{-1,1\}$. Specifically,
\begin{equation}
y_i=\left\{\begin{array}{ccc}
-1&\mb{if}& \langle \b w,\b x^{(i)}\rangle -t<0\\
1&\mb{if}& \langle \b w,\b x^{(i)}\rangle -t\ge0
\end{array}\right. .
\label{eq:outoftlulel}
\end{equation}
The goal of the TLU (perceptron)\index{perceptron!to classify} is to correctly classify the set of given input vectors $\b x^{(1)}, \ldots,\b x^{(p)}$ into one of the two classes, $\mb{\Sc C}_1$ or $\mb{\Sc C}_2$. In other words, if a given vector, say $\b x^{(i)}$, belongs to class $\mb{\Sc C}_1$, then the corresponding output $y_i$ should be $+1$. On the other hand, if $\b x^{(i)}$ belongs to class $\mb{\Sc C}_2$ then the TLU's output should be $-1$. Thus, we are faced with a two-class pattern classification problem.\index{two-class pattern!classification problem} The TLU's role is to find a hyperplane that separates two decision regions.\index{TLU's role!find a separating hyperplane} The separating hyperplane\index{separating hyperplane} (decision boundary)\index{separating hyperplane|see{decision boundary}}\index{decision boundary|see{separating hyperplane}} has the form
\begin{equation}
\langle \b w^*,\b x\rangle -t^*=\sum_{i=1}^n w_i^*x_i -t^*=0.
\label{eq:decsi_boun}
\end{equation}
In the case of two variables $x_1$ and $x_2$ the decision boundary is a straight line. An illustration of such a case is given in Figure~\ref{fig:lsep_ex}. In this case a point that lies above the boundary line is assigned to class $\mb{\Sc C}_1$, while a point below the boundary line is assigned to class $\mb{\Sc C}_2$. We present an iterative algorithm  for the decision boundary weights and threshold provided that a decision boundary exists. We assume that we are given $p$ training pairs
\[
\left(\b x^{(1)},d_1\right),\left(\b x^{(2)},d_2\right),\ldots,\left(\b x^{(p)},d_p\right),
\]
where for $i=1,2,\ldots,p$, the input vectors $\b x^{(i)} \in \mb{\Bb R}^n$ and $d_i \in \{-1,1\}$ represent the desired responses. We would like to determine a hyperplane~(\ref{eq:decsi_boun}) such that each input vector lies on a pre-assigned side of the hyperplane. If such a hyperplane exists for the given training set, then the set of input vectors is referred to as {\em linearly separable\/},\index{linearly separable} and the hyperplane is referred to as a {\em separating hyperplane\/}.\index{separating hyperplane} If no such hyperplane exists, then the training set is referred to as {\em linearly nonseparable}.\index{linearly nonseparable} Finding a separating hyperplane is equivalent to finding a corresponding weight vector and a threshold. If we find such a weight vector\index{weight vector!TLU} and a threshold,\index{threshold!TLU} then we say that the TLU correctly classifies the training set.\index{training set} Thus, designing a TLU that correctly classifies the given training set amounts to finding a weight vector $\b w^*$ and a threshold $t^*$ such that for $i=1,2,\ldots,p$, we have
\[
\begin{array}{ccc}
\langle \b w^*,\b x^{(i)}\rangle -t^*<0&\mb{if}& d_i=-1\\
\langle \b w^*,\b x^{(i)}\rangle -t^*\ge0&\mb{if}& d_i=1.
\end{array}
\]
In our discussion, we assume the threshold value, $t$, to be equal zero. We can achieve this by increasing the dimension of every input vector, $\b x^{(i)}$, by augmenting it with an entry equal to $-1$, and by increasing the dimension of the weight vector to $n+1$ by augmenting it with $t$. Indeed,
\[
y_i=\mb{sign}\left( \langle \b w,\b x^{(i)}\rangle -t\right)= \mb{sign}\left( \left\langle \tilde{\b w},\tilde{\b x}^{(i)}\right\rangle\right),
\]
where $\tilde{\b w}^{\top}=\left[\begin{array}{ccccc}
w_1&w_2&\cdots&w_n&t
\end{array}\right]$ and $\tilde{\b x}^{(i){T}}=\left[\begin{array}{ccccc}
x_{1}^{(i)}&x_{2}^{(i)}&\cdots&x_{n}^{(i)}&-1
\end{array}\right]$. In our further discussion, we assume that the given input vectors have been already augmented. In view of the above, we can formulate the training problem of a TLU\index{TLU!training problem} as follows:
\begin{center}
{\bf TLU Training Problem}
\end{center}
Given $p$ training pairs $\left(\b x^{(1)},d_1\right),\ldots,\left(\b x^{(p)},d_p\right)$, where $\b x^{(i)} \in \mb{\Bb R}^{n+1}$ and $d_i \in \{-1,1\}$. Determine a weight vector $\b w \in \mb{\Bb R}^{n+1}$, if there exists one, such that for $i=1,\ldots,p$,
\[
\mb{sign}\left( \left\langle \b w,\b x^{(i)}\right\rangle\right)=d_i.
\]
\markex
Assume that the given vectors $\b x^{(1)},\ldots, \b x^{(p)}$ are linearly separable, that is, there exists a weight vector $\b w^*$ such that
\[
\left\langle\b w^*,\b x^{(i)}\right\rangle\ge 0\quad\mb{for every input vector $\b x^{(i)}$ from class}\quad \mb{\Sc C}_1
\]
and
\[
\left\langle\b w^*,\b x^{(i)}\right\rangle< 0\quad\mb{for every input vector $\b x^{(i)}$ from class}\quad \mb{\Sc C}_2.
\]
We can now formulate the algorithm for adapting the weight vector\index{algorithm!for adapting weight vector} resulting in a separating hyperplane\index{separating hyperplane} of a given set of linearly separable points.\index{linearly separable points} The idea is very simple and we proceed as follows. If the given vector $\b x^{(i)}$ is classified correctly, do not update the weights, that is,
\[
\b w(k+1)=\b w(k).
\]
Otherwise,
\begin{equation}
\b w(k+1)=\b w(k)-\b x^{(i)}\quad\mb{if}\quad \b w^{\top}(k)\b x^{(i)}\ge 0\quad\mb{and}\quad\b x^{(i)} \in \mb{\Sc C}_2
\label{eq:firupofwe}
\end{equation}
and
\begin{equation}
\b w(k+1)=\b w(k)+\b x^{(i)}\quad\mb{if}\quad \b w^{\top}(k)\b x^{(i)}< 0\quad\mb{and}\quad\b x^{(i)} \in \mb{\Sc C}_1.
\label{eq:secweup}
\end{equation}
Note that if $\b w^{\top}(k)\b x^{(i)}\ge 0$, that is, $y(k)=1$ and $\b x^{(i)} \in \mb{\Sc C}_2$, that is, $d_i=-1$, then the error between the actual TLU output and the desired TLU output\index{error!between actual and desired TLU outputs} is $d_i-y(k)=-2$. Therefore, we can represent the weight update in~(\ref{eq:firupofwe}) as
\begin{equation}
\b w(k+1)=\b w(k)+\frac{1}{2}(d_i-y(k))\b x^{(i)}.
\label{eq:modweiupfor}
\end{equation}
Analogously, if $\b w^{\top}(k)\b x^{(i)}< 0$, that is, $y(k)=-1$, and $\b x^{(i)} \in \mb{\Sc C}_1$, then $d_i-y(k)=2$. Here, as in the previous case, we can represent the weight update in~(\ref{eq:secweup}) as~(\ref{eq:modweiupfor}). We now summarize the above considerations in the form of the following algorithm.
\begin{algorithm}\rm
\begin{center}
{\bf TLU (Perceptron) Learning}
\end{center}
Given a set of training pairs $\left\{(\b x^{(1)},d_1),\ldots,(\b x^{(i)},d_i),\ldots,(\b x^{(p)},d_p)  \right\}$.
\begin{enumerate}
\item {\bf Initialize Weight Vector and Threshold}\\
Set the values of the components of the initial weight vector $\b w(0)$ and threshold $t(0)$ to zero or to small random values.
\item {\bf Present New Input Vector and Desired Output}\\
Present the new input vector, $\b x^{(i)}$ and the corresponding desired output $d_i$.
\item {\bf Calculate Actual Output}
\[
y(k)=\mb{sign}\left( \left\langle \b w(k),\b x^{(i)} \right\rangle-t(k)\right).
\]
\item {\bf Adapt Weight Vector and Threshold}
\begin{eqnarray*}
\b w(k+1)&=&\b w(k)+\frac{1}{2}(d_i-y(k))\b x^{(i)},\\
t(k+1) &=& t(k)-\frac{1}{2}(d_i-y(k)).
\end{eqnarray*}
Note that the weight vector and threshold are unchanged if the correct decision is made by the TLU.
\item {\bf Repeat by Going to Step 2}
\end{enumerate}
\end{algorithm}\markex

\begin{example}\rm
We used the above algorithm to generate a separating line\index{separating line} for eight training pairs. For $i=1,\ldots,4$ and $d_i=1$,
\[
\left[\begin{array}{cccc}
\b x^{(1)} & \b x^{(2)} & \b x^{(3)} & \b x^{(4)}
\end{array}\right]= \left[\begin{array}{cccc}
-5 & -4 & 3 & 6\\
-4  & 6 & 4 & 7
\end{array}\right].
\]
The above points are labeled ``x'' in Figure~\ref{fig:lsep_ex}. The other four training pairs corresponding to $d_i=-1$, $i=5,\ldots,8$, labeled ``o'' are
\[
\left[\begin{array}{cccc}
\b x^{(5)} & \b x^{(6)} & \b x^{(7)} & \b x^{(8)}
\end{array}\right]= \left[\begin{array}{cccc}
-6 & 0 & 3 & 8\\
-8  & -1 & -4 & 3
\end{array}\right]
\]
The weight vector and threshold were initialized randomly, where
\[
\b w(0)=\left[\begin{array}{cc}
-0.4860 & 0.8913
\end{array}\right]^{\top}\quad\mb{and}\quad t_0=0.7621
\]
The TLU learning algorithm\index{TLU!learning algorithm} generated a separating line, shown in Figure~\ref{fig:lsep_ex}, corresponding to the weight vector $\b w^*=\left[\begin{array}{cc}
-5.5140 & 6.8913
\end{array}\right]^{\top}$ and threshold $t^*=-2.2379$.
\begin{figure}
\centering
\includegraphics[scale=0.75]{lsep_ex.eps}
%\vspace{3.5in}
%\centerbmp{4.67in}{3.5in}{c:/mybook/lsep_ex.bmp}
%\hskip1in\special{bmp:c:/mybook/lsep_ex.bmp x=4.67in y=3.5in}
\caption{An example of separable points in two-dimensional space and a separating line.}
\label{fig:lsep_ex}
\end{figure}
\end{example}
\markex

Our goal is to prove that the TLU learning algorithm\index{TLU!learning algorithm} can always find a separating hyperplane for linearly separable input vectors. In our proof, we use some technical results that we present next.
\begin{lemma}
\label{lemma:newt}
Let $\b x^{(1)},\b x^{(2)},\ldots,\b x^{(p)}$, $\b x^{(i)} \in \mb{\Bb R}^n$, be a set of $p$ input vectors, $\b w \in \mb{\Bb R}^n$ be a weight vector of an $n$-input TLU, and $t\in \mb{\Bb R}$ its threshold. Then, there exists a threshold $t^{\prime}$ such that, for $i=1,2,\ldots,p$,
\[
\left\langle \b w,\b x^{(i)} \right\rangle -t^{\prime} \ne 0,
\]
and
\[
\mb{\rm sign}\left(\left\langle \b w,\b x^{(i)} \right\rangle -t \right)=\mb{\rm sign}\left( \left\langle \b w,\b x^{(i)} \right\rangle -t^{\prime} \right).
\]
\end{lemma}
{\bf Proof}\quad Assume, without loss of generality, that the input vectors are sorted into two sets $\left\{\b x^{(1)},\ldots,\b x^{(l)}\right\}$ and $\left\{\b x^{(l+1)},\ldots,\b x^{(p)}\right\}$ such that the desired response of the vectors from the first set is $d_i=-1$, while the desired response of the vectors from the second set is $d_i=1$. Let
\[
\delta=\min\left\{\left| \langle \b w,\b x^{(i)} \rangle -t\right|:1\le i \le l  \right\}.
\]
It follows from~(\ref{eq:outoftlulel}) that $\delta>0$. Let $t^{\prime}=t-\delta/2$. Then,
\[
\langle \b w,\b x^{(i)} \rangle -t^{\prime} =\langle \b w,\b x^{(i)} \rangle -t + \delta/2 \ne 0 \quad \mb{for}\quad i=1,\ldots,p,
\]
and the output of the TLU with the modified threshold,\index{modified threshold}\index{threshold!modified} $t^{\prime}$, is the same as the output of the TLU with the old threshold, $t$.\\
\hspace*{5.5in}$\Box$

By the above lemma, we can always assume that $\langle \b w,\b x^{(i)} \rangle -t^{\prime}\ne 0$, for $i=1,\ldots,p$ and that the output of a TLU can be defined as
\[
y_i=\left\{\begin{array}{ccc}
-1&\mb{if}& \langle \b w,\b x^{(i)}\rangle -t<0\\
1&\mb{if}& \langle \b w,\b x^{(i)}\rangle -t>0.
\end{array}\right.
\]
We further simplify the formulation of the TLU learning problem by preprocessing training pairs.\index{preprocessing!training pairs} Specifically, we replace every training pair, $\left(\b x^{(i)},d_i\right)$ with $d_i=-1$ with their negative, $\left(-\b x^{(i)},-d_i\right)$. After performing the above preprocessing, the TLU learning problem can be reformulated in a simple form that we use while analyzing the TLU learning algorithm.
\begin{center}
{\bf Simplified TLU Learning Problem}
\end{center}
Given a set of input vectors $\left\{ \b x^{(1)},\ldots,\b x^{(p)} \right\}$. Determine a weight vector $\b w$, if there exists one, such that for $i=1,\ldots,p$,
\[
\left\langle \b w,\b x^{(i)}\right\rangle>0.
\]
\markex
In view of the above, for a linearly separable set of vectors, $\left\{\b x^{(i)}\right\}$, after preprocessing there exists a hyperplane passing through the origin in $\mb{\Bb R}^{n+1}$ such that all the vectors lie on one side of the hyperplane. The normal to the hyperplane is the weight vector $\b w$. This is illustrated in Figure~\ref{fig:tlu_sep}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{tlu_sep.eps}
%\vspace{3.5in}
%\hskip1.4in\special{eps:c:/mybook/tlu_sep.eps y=3.5in}
\caption{Linearly separable vectors after their preprocessing.}
\label{fig:tlu_sep}
\end{figure}
Taking into account the simplified TLU learning problem,\index{simplified TLU!learning problem}\index{TLU!simplified learning problem} we can state the TLU learning algorithm equivalently as follows:
\begin{center}
{\bf Simplified TLU Learning Algorithm}
\end{center}
\begin{enumerate}
\item {\bf Initialize Weight Vector}\\
Set the values of the components of the initial weight vector $\b w(0)$ to zero.
\item {\bf Test}\\
If $\left\langle \b w(k),\b x^{(i)}\right\rangle\le 0$, go to {\bf Add}, else {\bf Test} next input vector.
\item {\bf Add}
\[
\b w(k+1)=\b w(k)+\b x^{(i)},\quad k=k+1.
\]
Go to {\bf Test}.
\end{enumerate}
\markex
The algorithm is said to {\em converge\/}\index{algorithm!is said to converge}\index{converge!algorithm} if for some $k<\infty$,
\[
\left\langle \b w(k),\b x^{(i)} \right\rangle >0,\quad \mb{for}\quad i=1,\ldots,p.
\]

The following theorem states that if the set of input vectors is linearly separable, then the TLU learning algorithm converges.\index{TLU convergence theorem}
\begin{thm}[The TLU Convergence Theorem]
If the set of vectors $\left\{ \b x^{(1)},\ldots,\b x^{(p)}\right\}$ is linearly separable, then the simplified TLU learning algorithm determines a weight vector $\b w$ in a finite number of steps, such that $\langle \b w,\b x^{(i)}\rangle>0$ for $i=1,2,\ldots,p$
\end{thm}
{\bf Proof}\quad The proof is adapted from Siu, Roychowdhury, and Kailath~\cite{siu95}. We begin by observing that because the set of input vectors is linearly separable, there exists a weight vector $\b w^*$ and a constant $\delta>0$ such that for $i=1,2,\ldots,p$, we have
\[
\langle \b w^*,\b x^{(i)} \rangle >\delta.
\]
To simplify further manipulations, we assume that $\left\Arrowvert \b w^*\right\Arrowvert=1$. We also assume that $\left\Arrowvert \b x^{(i)}\right\Arrowvert\le L$, for $i=1,2,\ldots,p$, that is, the lengths of all the input vectors are bounded above by $L$. Let $k$ denote the number of times the weight vector $\b w$, in the simplified TLU learning algorithm,\index{simplified TLU!learning algorithm} has been updated, and let $\b w(k)$ be the value of the weight vector, $\b w$, after the $k$-th update. Our goal is to show that $k$ is bounded above. To proceed, we note that if an input vector $\b x^{(i)}$ is used for the $k$-th update in the algorithm, then
\[
\b w(k)=\b w(k-1)+\b x^{(i)},
\]
where $\langle \b w(k-1),\b x^{(i)} \rangle =\b w(k-1)^{\top} \b x^{(i)} \le 0$. We first derive an upper bound on $\|\b w(k)\|$. We have,
\begin{eqnarray*}
\|\b w(k)\|^2 &=& \|\b w(k-1) + \b x^{(i)} \|^2\\
&=& \left( \b w(k-1) + \b x^{(i)} \right)^{\top}\left( \b w(k-1) + \b x^{(i)} \right)\\
&=& \b w(k-1)^{\top}\b w(k-1)+\b x^{(i)\top} \b x^{(i)}+ 2\b w(k-1)^{\top} \b x^{(i)}\\
&\le& \|\b w(k-1)\|^2 +L^2,
\end{eqnarray*}
because $\b w(k-1)^{\top}\b x^{(i)} \le 0$ and $\|\b x^{(i)}\|^2 \le L^2$. Suppose that for the $(k-1)$-st update, the input vector $\b x^{(j)}$ was used, that is,
\[
\b w(k-1)=\b w(k-2)+\b x^{(j)}.
\]
Hence,
\begin{eqnarray*}
\|\b w(k)\|^2 &\le& \|\b w(k-1)\|^2 +L^2\\
&=& \|\b w(k-2)+\b x^{(j)}\|^2+L^2\\
&\le& \|\b w(k-2)\|^2 + 2L^2.
\end{eqnarray*}
Continuing in this manner, we obtain,
\begin{equation}
\|\b w(k)\|^2 \le k L^2,
\label{eq:upbounonwk}
\end{equation}
because by assumption, $\b w(0)=\b 0$.

We next derive a lower bound on $\|\b w(k)\|$. We have,
\begin{eqnarray*}
\left\langle \b w^*,\b w(k) \right\rangle &=& \left\langle \b w^*,\b w(k-1) \right\rangle+ \left\langle \b w^*,\b x^{(i)} \right\rangle\\
&\ge& \left\langle \b w^*,\b w(k-1) \right\rangle + \delta,
\end{eqnarray*}
because we assumed that $\langle \b w^*,\b x^{(i)} \rangle >\delta$. Observing that $\b w(k-1)=\b w(k-2)+\b x^{(j)}$, we get
\begin{eqnarray*}
\left\langle \b w^*,\b w(k) \right\rangle &\ge&\left\langle \b w^*,\b w(k-1) \right\rangle + \delta\\
&\ge& \left\langle \b w^*,\b w(k-2) \right\rangle + 2\delta.
\end{eqnarray*}
Continuing in this manner, we obtain
\begin{equation}
\left\langle \b w^*,\b w(k) \right\rangle \ge k\delta.
\label{eq:lobouonwks}
\end{equation}
Applying to the above the Cauchy-Schwarz inequality\index{Cauchy-Schwarz inequality} yields
\[
\|\b w^*\|\|\b w(k)\| \ge \left\langle \b w^*,\b w(k) \right\rangle \ge k\delta.
\]
By assumption $\|\b w^*\|=1$. Hence,
\[
\|\b w(k) \| \ge k \delta.
\]
Combining the obtained upper and lower bounds on $\|\b w(k)\|^2$, given by~(\ref{eq:upbounonwk}) and~(\ref{eq:lobouonwks}), we have
\[
k^2\delta^2 \le \|\b w(k)\|^2 \le k L^2.
\]
Therefore,
\begin{equation}
k \le \frac{L^2}{\delta^2},
\label{eq:upboonstenu}
\end{equation}
which implies that the TLU learning algorithm converges. The proof of the TLU convergence theorem\index{TLU convergence theorem} is complete.\\
\hspace*{5.5in} $\Box$

We now present a sufficient condition for linear nonseparability.\index{linear nonseparability!sufficient condition for}
\begin{thm}
\label{thm:nonseparability}
If there exists a positive linear combination~(PLC)\index{positive linear combination (PLC)}\index{PLC|see{positive linear combination}} that equals $\b 0$ of the given set of input vectors, then the given set of input vectors is linearly nonseparable.
\end{thm}
{\bf Proof}\quad We prove the theorem by contradiction. Assume that there exists a PLC of the given input vectors that is equal to $\b 0$, that is, there exist nonnegative constants, $q_i$, such that for some $j$, $1\le j\le p$, we have $q_j>0$, and
\[
\sum_{i=1}^p q_i\b x^{(i)} =\b 0.
\]
Suppose now that the given input vectors are linearly separable, that is, there is a weight vector, $\b w$, such that $\left\langle \b w,\b x^{(i)}\right\rangle >0$ for $i=1,\ldots,p$. We premultiply $\sum_{i=1}^p q_i\b x^{(i)}=\b 0$ by $\b w^{\top}$ to get
\[
\b w^{\top}\sum_{i=1}^p q_i\b x^{(i)} = \sum_{i=1}^p q_i\langle\b w,\b x^{(i)} \rangle
> 0
\]
since $\left\langle \b w,\b x^{(i)}\right\rangle >0$ and $q_j>0$. However, we arrived at a contradiction because by assumption $\sum_{i=1}^p q_i\b x^{(i)} =\b 0$, which completes the proof.\\
\hspace*{5.5in} $\Box$
\begin{example}\label{ex:tlu_xor}\rm
Suppose we are given the following training set:
\begin{eqnarray*}
\b x^{(1)}&=&\left[\begin{array}{c}
-1\\
-1
\end{array}\right],\; d_1=-1;\quad \b x^{(2)}=\left[\begin{array}{c}
1\\
1
\end{array}\right],\; d_2=-1;\\
\b x^{(3)}&=&\left[\begin{array}{c}
-1\\
1
\end{array}\right],\; d_3=1;\quad \b x^{(4)}=\left[\begin{array}{c}
1\\
-1
\end{array}\right],\; d_4=1.
\end{eqnarray*}
The preceding training set represents the 2-input Exclusive-OR (XOR) function.\index{Exclusive-OR (XOR) function}\index{XOR|see{Exclusive-OR}} Let us now preprocess the training set so we can apply Theorem~\ref{thm:nonseparability}. By preprocessing we mean eliminating the threshold by augmenting the input\index{eliminating threshold!by augmenting the input vector} vectors with $-1$ and then negating the input vectors for which the desired output is $-1$. The preprocessed input vectors are:
\[
\tilde{\b x}^{(1)}=\left[\begin{array}{c}
1\\
1\\
1
\end{array}\right],\quad \tilde{\b x}^{(2)}=\left[\begin{array}{c}
-1\\
-1\\
1
\end{array}\right],\quad \tilde{\b x}^{(3)}=\left[\begin{array}{c}
-1\\
1\\
-1
\end{array}\right],\quad \tilde{\b x}^{(4)}=\left[\begin{array}{c}
1\\
-1\\
-1
\end{array}\right].
\]
Note that $\sum_{i=1}^4\tilde{\b x}^{(i)}=\b 0$. It follows from Theorem~\ref{thm:nonseparability} that the given input vectors are not linearly separable.\index{linearly separable!vectors}
\end{example}
\markex

%\section{Identification Using Adaptive Linear Element}
%\label{sec:idadaline}
%The adaptive linear element, also called the adaline, is shown in Figure~\ref{fig:adaline}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{adaline.eps}
%%\vspace{3in}
%%\hskip1.5in\special{eps:c:/mybook/adaline.eps y=3in}
%\caption{Adaline.}
%\label{fig:adaline}
%\end{figure}
%We first explain the symbols used in our discussion of its operation. Let $\b x(t)=\left[\begin{array}{ccc}
%x_1(t) & \cdots & x_n(t)
%\end{array}\right]^{\top}$ denote the adaline input. We assume that for all $t$
%\[
%\|\b x(t)\|=\sqrt{x_1^2(t)+\cdots + x_n^2(t)} \le B_x,
%\]
%where $B_x>0$ is known positive constant and $\|\cdot \|=\|\cdot\|_2$, that is, whenever we write $\|\cdot\|$ we mean the Euclidean norm $\|\cdot \|_2$. We also assume that for all $t$,
%\[
%\|\dot{\b x}(t)\|_2=\sqrt{\dot{x}_1^2(t)+\cdots + \dot{x}_n^2(t)} \le B_{\dot{x}},
%\]
%The weight vector $\b w(t)=\left[\begin{array}{ccc}
%w_1(t) & \cdots & w_n(t)
%\end{array}\right]^{\top}$. The bias weight $w_{n+1}(t)$ is connected to a constant input equal to $-1$. We define $\tilde{\b x}(t)$ to be an augmented input vector that includes the constant input to the bias weight, that is,
%\[
%\tilde{\b x}(t)=\left[\begin{array}{cccc}
%x_1(t) & \cdots & x_n(t) & x_{n+1}(t)
%\end{array}\right]^{\top}.
%\]
%Note that
%\[
%\left\Arrowvert \tilde{\b x}(t) \right\Arrowvert_2^2=\tilde{\b x}^{\top}(t)\tilde{\b x}(t)=\left\Arrowvert \b x(t) \right\Arrowvert_2^2 +1,
%\]
%and hence $\tilde{\b x}(t)$ is bounded away from zero at all times. The augmented vector $\tilde{\b x}(t)$ is bounded, for all $t$, if and only if $\b x(t)$ is bounded. We denote a bound of $\tilde{\b x}(t)$ by $B_{\tilde{x}}$, that is, $\left\Arrowvert \tilde{\b x}(t) \right\Arrowvert \le B_{\tilde{x}}$ for all $t$. Observe that $\left\Arrowvert \dot{\tilde{\b x}}(t) \right\Arrowvert = \left\Arrowvert \dot{\b x}(t) \right\Arrowvert$. Thus, we can take $B_{\dot{\tilde{x}}}=B_{\dot{x}}$. Similarly, we define an augmented weight vector,
%\[
%\tilde{\b w}(t)=\left[\begin{array}{cccc}
%w_1(t) & \cdots & w_n(t) & w_{n+1}(t)
%\end{array}\right]^{\top}.
%\]
%We assume that for all $t$,
%\[
%\left\Arrowvert \tilde{\b w}(t) \right\Arrowvert \le B_{\tilde{w}}
%\]
%for some constant $B_{\tilde{w}}>0$. The output, $y(t)$, of the adaline is a linear combination of its inputs. That is,
%\[
%y(t)=\tilde{\b x}^{\top}(t)\tilde{\b w}(t).
%\]
%We denote by $y_d(t)$ the desired output of the adaline. We assume that
%\[
%\left| y_d(t)\right|\le B_y\quad\mb{and}\quad \left| \dot{y}_d(t)\right|\le B_{\dot{y}}
%\]
%for some positive constants $B_y$ and $B_{\dot{y}}$. Let
%\[
%e(t)=y(t)-y_d(t)
%\]
%denote the error signal. The error dynamics are described by
%\begin{eqnarray}
%\dot{e} &=& \dot{y}-\dot{y}_d\nonumber\\
%&=& \dot{\tilde{\b w}}^{\top}\tilde{\b x}+\tilde{\b w}^{\top}\dot{\tilde{\b x}}-\dot{y}_d.
%\label{eq:err_dyna}
%\end{eqnarray}
%Our objective is to devise an adaptation algorithm for the weights $w_i(t)$, $i=1,2,\ldots,n+1$, so that the error $e(t)$ converges to zero in finite time, say $t_h$, and $e(t)=0$ for $t>t_h$. Sira-Ram{\'{\i}}rez and Colina-Morles~\cite{sira95} showed that we can achieve our objective if we adapt the weights according to the adaptation law,
%\begin{equation}
%\fbox{$\displaystyle \dot{\tilde{\b w}}=-\eta\frac{\tilde{\b x}}{\tilde{\b x}^{\top}\tilde{\b x}}\,\mb{sign}(e(t)) $}
%\label{eq:adaptlawforwtil}
%\end{equation}
%where the gain $\eta>0$ is sufficiently large. We now state and prove the theorem of Sira-Ram{\'{\i}}rez and Colina-Morles~\cite{sira95}.
%\begin{thm}
%If the weight vector $\tilde{\b w}$ is adapted according to~(\ref{eq:adaptlawforwtil}), where
%\[
%\eta>B_{\tilde{w}}B_{\dot{x}}+B_{\dot{y}},
%\]
%then for arbitrary initial error $e(0)$ the error signal $e(t)$ converges to zero in finite time, say $t_h$, such that
%\[
%t_h\le \frac{|e(0)|}{\eta-B_{\tilde{w}}B_{\dot{x}}-B_{\dot{y}}}
%\]
%and $e(t)=0$ for $t>t_h$.
%\label{thm:adaforfth}
%\end{thm}
%{\bf Proof}\quad Recall that the error dynamics are governed by~(\ref{eq:err_dyna}), where the vector $\tilde{\b w}(t)$ is adapted according to~(\ref{eq:adaptlawforwtil}). Consider the following Lyapunov function candidate
%\[
%V(e)=\frac{1}{2}e^2(t).
%\]
%The time derivative of $V$ evaluated on the solutions of~(\ref{eq:err_dyna}) is
%\begin{eqnarray}
%\dot{V}(e) &=& e\dot{e}\nonumber\\
%&=& e\left( \dot{\tilde{\b w}}^{\top}\tilde{\b x}+\tilde{\b w}^{\top}\dot{\tilde{\b x}}-\dot{y}_d\right).
%\label{eq:vdotforerrada}
%\end{eqnarray}
%Substituting~(\ref{eq:adaptlawforwtil}) into the above and taking into account that
%\[
%\mb{sign}(e)=\frac{e}{|e|}=\frac{|e|}{e},
%\]
%we obtain
%\begin{eqnarray}
%\dot{V} &=& -\eta|e| + e\left(\tilde{\b w}^{\top}\dot{\tilde{\b x}}-\dot{y}_d \right)\nonumber\\
%&\le& -\eta|e| +|e|\left( B_{\tilde{w}}B_{\dot{\tilde{x}}}+B_{\dot{y}}\right)\nonumber\\
%&=& \left(-\eta + B_{\tilde{w}}B_{\dot{\tilde{x}}}+B_{\dot{y}}\right)|e|\nonumber\\
%&<& 0\quad\mb{for all}\quad e\ne 0\quad\mb{and}\quad \eta > B_{\tilde{w}}B_{\dot{\tilde{x}}}+B_{\dot{y}},
%\label{vdotsmalforallt}
%\end{eqnarray}
%which means that $e(t)\to 0$ as $t$ increases. It follows from our discussion above that we can take $B_{\dot{\tilde{x}}}=B_{\dot{x}}$. We now show that for some finite time $t_h$, we have $e(t_h)=0$ and $e(t)=0$ for $t>t_h$. For this substitute~(\ref{eq:adaptlawforwtil}) into~(\ref{eq:err_dyna}) to get
%\begin{equation}
%\dot{e}=-\eta\,\mb{sign}(e)+\tilde{\b w}^{\top}\dot{\tilde{\b x}}-\dot{y}_d.
%\label{eq:erreqafsub}
%\end{equation}
%Integrating the above from $t=0$ to $t\le t_h$ gives
%\begin{equation}
%e(t)-e(0)=-\eta t\,\mb{sign}(e(0))+\int_0^{t} \left( \tilde{\b w}^{\top}(\tau)\dot{\tilde{\b x}}(\tau)-\dot{y}_d(\tau)\right)d\tau .
%\label{eq:aftintoth}
%\end{equation}
%At $t=t_h$, we have $e(t_h)=0$ and therefore for $t=t_h$ the above equation takes the form
%\begin{equation}
%-e(0)=-\eta t_h \,\mb{sign}(e(0))+\int_0^{t_h} \left( \tilde{\b w}^{\top}(t)\dot{\tilde{\b x}}(t)-\dot{y}_d(t)\right)dt .
%\label{eq:erratth}
%\end{equation}
%Multiplying both sides of the above by $-\mb{sign}(e(0))$ yields
%\begin{eqnarray}
%|e(0)|&=&\eta t_h-\mb{sign}(e(0))\int_0^{t_h} \left( \tilde{\b w}^{\top}(t)\dot{\tilde{\b x}}(t)-\dot{y}_d(t)\right)dt\nonumber\\
%&\ge& \left(\eta-\left( B_{\tilde{w}}B_{\dot{x}}+B_{\dot{y}}\right) \right) t_h.
%\label{eq:afmulbysign}
%\end{eqnarray}
%Thus, the time $t_h$ at which $e(t)$ reaches zero value for the first time satisfies the inequality
%\[
%t_h \le \frac{|e(0)|}{\eta-\left( B_{\tilde{w}}B_{\dot{x}}+B_{\dot{y}}\right)}.
%\]
%Because for $e(t)\ne 0$ we have $\frac{d}{dt}e^2(t)<0$, we are guaranteed that we also have $e(t)=0$ for $t>t_h$. The proof of the theorem is now complete.\\
%\hspace*{5.5in} $\Box$
%
%If the time derivative $\dot{\tilde{\b x}}(t)$ is available, then we can use it in an adaptation law of $\tilde{\b w}$. Using $\dot{\tilde{\b x}}(t)$ in an adaptation law allows us to reduce the gain $\eta$ and the bound on time $t_h$ when the error $e(t)$ becomes zero. One such adaptation strategy that includes $\dot{\tilde{\b x}}(t)$ was proposed by Sira-Ram{\'{\i}}rez and Colina-Morles~\cite{sira95}. We present this strategy in the following theorem.
%\begin{thm}
%If the adaptation law for the augmented weight vector is chosen as
%\begin{equation}
%\fbox{$\displaystyle \dot{\tilde{\b w}}(t)=-\frac{\tilde{\b x}(t)\dot{\tilde{\b x}}^{\top}(t)}{\tilde{\b x}^{\top}(t)\tilde{\b x}(t)}\tilde{\b w}(t)-\eta \frac{\tilde{\b x}(t)}{\tilde{\b x}^{\top}(t)\tilde{\b x}(t)}\,\mb{sign}(e(t)) $}
%\label{eq:adalfwwxd}
%\end{equation}
%where $\eta > B_{\dot{y}}$, then for any given initial error $e(0)$, the error $e(t)$ converges to zero in finite time $t_h$ such that
%\begin{equation}
%t_h \le \frac{|e(0)|}{\eta-B_{\dot{y}}}
%\label{eq:finthwxd}
%\end{equation}
%end $e(t)=0$ for $t>t_h$.
%\label{thm:adawxdo}
%\end{thm}
%{\bf Proof}\quad The proof of this theorem is analogous to the proof of Theorem~\ref{thm:adaforfth}. The only difference is that (\ref{eq:erreqafsub}) becomes
%\begin{equation}
%\dot{e}=-\eta\,\mb{sign}(e)-\dot{y}_d.
%\label{eq:rederreq}
%\end{equation}
%\hspace*{5.5in} $\Box$
%
%When using the adaline for control purposes, the signal to be processed by the adaline is first passed through a prefilter as shown in Figure~\ref{fig:pre_ada}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{pre_ada.eps}
%%\vspace{1.5in}
%%\hskip1in\special{eps:c:/mybook/pre_ada.eps y=1.5in}
%\caption{Prefilter connection to adaline.}
%\label{fig:pre_ada}
%\end{figure}
%The prefilter is described by the equation, $\dot{\b x}=\b A\b x+\b b r$, where $r$ is the signal to be processed. The matrix $\b A$ of the prefilter is chosen to be asymptotically stable. Recall that the inputs to the adaline are the components of the augmented vector $\tilde{\b x}$ where the $(n+1)$-st component is equal to $-1$. Thus, $\dot{x}_{n+1}=0$. Taking this into account, the prefiltered inputs to the adaline can be modeled as
%\begin{equation}
%\dot{\tilde{\b x}}=\tilde{\b A}\tilde{\b x}+\tilde{\b b}r,
%\label{eq:prefinta}
%\end{equation}
%where
%\[
%\tilde{\b A}=\left[\begin{array}{cc}
%\b A & \b 0\\
%\b 0^{\top} & 0
%\end{array}\right]\quad\mb{and}\quad \tilde{\b b}=\left[\begin{array}{c}
%\b b\\
%0
%\end{array}\right].
%\]
%Note that $\dot{\tilde{\b x}}$ is available to us in the architecture shown in Figure~\ref{fig:pre_ada}. Thus, we can employ the adaptation strategy given by~(\ref{eq:adalfwwxd}), which results in the following description of the structure depicted in Figure~\ref{fig:pre_ada}:
%\begin{equation}
%\left.\begin{array}{lll}
%\dot{\tilde{\b x}}&=& \tilde{\b A}\tilde{\b x}+\tilde{\b b}r\\
%\dot{\tilde{\b w}}&=&-\frac{\tilde{\b x}\dot{\tilde{\b x}}^{\top}}{\tilde{\b x}^{\top}\tilde{\b x}}\tilde{\b w} - \eta \frac{\tilde{\b x}}{\tilde{\b x}^{\top}\tilde{\b x}}\,\mb{sign}(e)\\
%y&=& \tilde{\b x}^{\top}\tilde{\b w}
%\end{array}\right\}
%\label{eq:combdescaf}
%\end{equation}
%We now demonstrate, with a numerical example, how the adaline system modeled by~(\ref{eq:combdescaf}) can be applied to forward and inverse dynamics identification of dynamical systems. The nonlinear plant model is taken from Sira-Ram{\'{\i}}rez and Colina-Morles~\cite{sira95}. This plant model has the form
%\begin{equation}
%\left.\begin{array}{lll}
%\dot{\alpha}&=&p+\frac{\sin(\alpha)}{l}u\\
%\dot{p} &=& \left( \frac{g}{l}-\frac{u^2}{l^2}\cos(\alpha)\right)-\frac{u}{l}p\cos(\alpha)\\
%\dot{z}&=&u\\
%y_p &=& \alpha
%\end{array}\right\}
%\label{eq:kapitsamo}
%\end{equation}
%The parameters $g$ and $l$ have the values: $g=9.81$\,\,\,m/sec$^2$ and $l=0.7$\,\,\,m. The plant input is labeled as $u$ while the plant output is $y_p$. We first concern ourselves with forward dynamics identification of the above plant. A block diagram of forward dynamics identification scheme is depicted in Figure~\ref{fig:forward}.
%\begin{figure}
%%\vspace{2in}
%\centerline{\psfig{figure=forward.eps,height=2in}}
%%\hskip2in\special{eps:c:/mybook/forward.eps y=2in}
%\caption{Forward dynamics identifier scheme.}
%\label{fig:forward}
%\end{figure}
%The input to the plant is
%\[
%u(t)=A_1+A_2\cos\left(\frac{t}{2\epsilon}\right)+A_3\sin\left(\frac{t}{\epsilon}\right),
%\]
%where $A_1=0.4$, $A_2=2$, $A_3=3$, and $\epsilon=0.05$. The prefilter equation~(\ref{eq:prefinta}) in our example has the form
%\[
%\dot{\tilde{\b x}}=\left[\begin{array}{cccc}
%0 & 1 & 0 & 0\\
%0 & 0 & 1 & 0\\
%-1 & -3 & -3 & 0\\
%0 & 0 & 0 & 0
%\end{array}\right]\tilde{\b x}+\left[\begin{array}{c}
%0\\
%0\\
%1\\
%0
%\end{array}\right]u.
%\]
%Thus, the adaline has total four weights, including the bias weight. In our simulations we use the following approximation of the sign function,
%\begin{equation}
%\mb{sign}\,(e(t))\approx \frac{e(t)}{|e(t)|+\delta},
%\label{eq:sign_approx}
%\end{equation}
%where $\delta=0.05$. We used the above approximation to eliminate chattering and speed up the simulations. We used $\eta=10$. The performance of the adaline as a forward dynamics identifier of the dynamical system modeled by~(\ref{eq:kapitsamo}) is shown in Figure~\ref{fig:sim_for}.
%\begin{figure}
%%\vspace{3in}
%\centerline{\psfig{figure=sim_for.eps,height=3in}}
%%\hskip1.25in\special{eps:c:/mybook/sim_for.eps y=3in}
%\caption{Simulation results of forward dynamics identification of~(\ref{eq:kapitsamo}).}
%\label{fig:sim_for}
%\end{figure}
%We used the following initial conditions:
%\[
%\tilde{\b x}(0)=\left[\begin{array}{cccc}
%0.1 & 0 & 0 & -1
%\end{array}\right]^{\top},\quad \tilde{\b w}(0)=\left[\begin{array}{cccc}
%0.1 & 0 & 0 & 0
%\end{array}\right]^{\top},
%\]
%and $\left[\begin{array}{ccc}
%\alpha(0) & p(0) & z(0)
%\end{array}\right]=\left[\begin{array}{ccc}
%0.1 & 0 & 0
%\end{array}\right]$. Increasing the gain $\eta$ results in more accurate identification, that is, the error $e(t)$ goes to zero faster. Furthermore, the effects of using the approximation~(\ref{eq:sign_approx}) are reduced for higher gain $\eta$.
%
%In our next experiment, we employed the adaline as an inverse dynamics identifier using the scheme shown in Figure~\ref{fig:inve_id}.
%\begin{figure}
%%\vspace{2in}
%\centerline{\psfig{figure=inve_id.eps,height=2in}}
%%\hskip2in\special{eps:c:/mybook/inve_id.eps y=2in}
%\caption{Inverse dynamics identifier scheme.}
%\label{fig:inve_id}
%\end{figure}
%The results of our simulation, using the same initial conditions as before, are depicted in Figure~\ref{fig:inv_sim}. Here we used $\eta=90$. Increasing the gain $k$ decreases the identification error as well as reduces the effects of the approximation of the sign function.
%\begin{figure}
%%\vspace{3in}
%\centerline{\psfig{figure=inv_sim.eps,height=3in}}
%%\hskip1.25in\special{eps:c:/mybook/inv_sim.eps y=3in}
%\caption{Simulation results of inverse dynamics identification of~(\ref{eq:kapitsamo}).}
%\label{fig:inv_sim}
%\end{figure}

\section{Backpropagation}
In Section~\ref{sec:tlu}, we discussed an iterative learning algorithm for a threshold logic unit.\index{iterative learning algorithm!for TLU} This computational unit keeps on iterating until it reaches a solution, if a solution exists. In this section, we present iterative learning algorithms for computational units with differentiable activation functions. This family of algorithms are called the backpropagation algorithms because they back propagate the error between the desired and actual output of the neural network for the purpose of adjusting the network's weights in order to reduce the learning error.\index{learning error!reduce by adjusting weights} Specifically, a typical backpropagation algorithm\index{backpropagation!algorithm}\index{backpropagation} consists of two passes:\index{backpropagation algorithm!consists of two passes} the forward pass\index{forward pass!of backpropagation algorithm} during which the input vectors are presented to the network and the actual outputs of the network are calculated. Then, the error between the actual and desired outputs are calculated. In the reverse pass,\index{reverse pass!of backpropagation algorithm} the error gradient with respect to the network's weights is calculated by propagating the error backwards through the network. Once the error gradient is calculated, the weights are adjusted using, for example, a descent gradient method given by \eqref{eq:grad_desc_metho} in Subsection~\ref{subsec:gard_desc_meth}.

We first analyze a single computational unit shown in Figure~\ref{fig:cudiffa}. It is called the {\em perceptron}.\index{perceptron}
\begin{figure}
\centering
\includegraphics[scale=0.5]{cudiffa.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/cudiffa.bmp}
%\hskip1.5in\special{bmp:c:/mybook/cudiffa.bmp x=4in y=3in}
\caption{A single perceptron with a differentiable sigmoidal activation function.}
\label{fig:cudiffa}
\end{figure}
We consider a simple case, where we train the perceptron to learn the training pair\index{perceptron!training of}
\[
\left(\b x,d\right),
\]
where $\b x \in \mb{\Bb R}^{n+1}$, $x_{n+1}=\left(\b x \right)_{n+1}=-1$, and $d \in [-1,1]$. Define the error function\
\begin{eqnarray}
E(\b w)&=&\frac{1}{2}\left(d-y \right)^2\nonumber\\
&=& \frac{1}{2}\left(d-f(net) \right)^2,
\label{eq:percrfsper}
\end{eqnarray}
where $net=\b x^{\top}\b w$. Note that $w_{n+1}=t$, where $t$ is a threshold, bias.\index{threshold!perceptron}\index{bias|see{threshold}} The gradient of $E$, with respect to $\b w$, is
\begin{equation}
\nabla E(\b w)=-(d-f(net))f^{'}(net)\b x.
\label{eq:gradjww}
\end{equation}
Using the descent gradient algorithm with a fixed step size\index{descent gradient algorithm!with a fixed step size} to minimize $E$, we obtain
\begin{eqnarray}
\b w(k+1)&=& \b w(k)-\alpha \nabla E\left(\b w(k) \right)\nonumber\\
&=& \b w(k)+\alpha \left(d_k-f\left( net_k\right) \right)f^{'}\left( net_k\right)\b x(k),
\label{eq:desgwfss}
\end{eqnarray}
where $net_k=\b x(k)^{\top}\b w(k)$. Let
\begin{equation}
\delta ^k=\left(d_k-f\left( net_k\right) \right)f^{'}\left( net_k\right).
\label{eq:defofdsp}
\end{equation}
Then, we represent~(\ref{eq:desgwfss}) as
\begin{equation}
\fbox{$\displaystyle\b w(k+1)= \b w(k)+\alpha\delta^k\b x(k)$}
\label{eq:sinpruad}
\end{equation}

We generalize the above formula for a three-layer feedforward neural network\index{three-layer!feedforward neural network}\index{feedforward neural network} shown in Figure~\ref{fig:fuzzybp}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{fuzzybp.eps}
\caption{A three-layer feedforward neural network.}
\label{fig:fuzzybp}
\end{figure}
We wish to train the network to learn the training pairs
\[
\left(\b x^{(1)},\b d^{(1)}\right),\left(\b x^{(2)},\b d^{(2)}\right),\ldots,\left(\b x^{(p)},\b d^{(p)}\right),
\]
where $\b x^{(i)} \in \mb{\Bb R}^{n}$ and $\b d^{(i)} \in \mb{\Bb R}^{m}$. The network has $n$ inputs, which is the dimension of the training vectors, and $m$ outputs. The network's inputs are connected to a layer of hidden units, which in turn is connected to a layer of output units. To construct a network to perform some task, we must set the weights of the connections appropriately. The weights specify the strength of the influence of one unit upon another. We derive rules for adjusting the weights $w_{ji}$ and $w_{kj}$ such that the error function\index{error function!of feedforward neural network}
\begin{eqnarray*}
E(\b w) &=& \frac{1}{2}\sum_{s=1}^m \left(d_s-y_s\right)^2\\
&=& \frac{1}{2}\sum_{s=1}^m \left(d_s-f_{o}(net_s)\right)^2
\end{eqnarray*}
is minimized, where
\begin{equation}
net_s=\sum_{j=1}^l w_{sj}^oz_j.
\label{eq:net_s_defn}
\end{equation}
The minimization is performed over the space of weights\index{space of weights!minimization over}
\[
\left\{w_{ji}^h,w_{sj}^o : i=1,2,\ldots,n,\: j=1,2,\ldots,l,\: s=1,2,\ldots,m\right\}.
\]
We first give formulas for updating the weights connecting the hidden and output layers.

Using the descent gradient algorithm\index{descent gradient algorithm} in a similar way as in~(\ref{eq:sinpruad}), we obtain
\begin{eqnarray}
w_{sj}^{o}(k+1)&=&w_{sj}^{o}(k)-\alpha_o \frac{\partial E\left(\b w(k) \right)}{\partial w_{sj}^{o}(k)}\nonumber\\
&=& w_{sj}^{o}(k)+\alpha_o\left(d_s-y_s \right)f_o^{'}(net_s)z_j(k),
\label{eq:rulefhidw}
\end{eqnarray}
where $f_o^{'}$ is the derivative of $f_{o}(net_s)$ with respect to $net_s$, and the values of $z_j$ are computed by feeding into the network the input vector $\b x$ and then propagating it forward through the network. We have
\[
z_j=f_h\left(\sum_{i=1}^n w_{ji}^h x_i\right)=f_h(net_j),\quad j=1,2,\ldots,l.
\]
To derive the update formulas for the weights\index{update formulas!for the weights} in the hidden layer,\index{hidden layer} we use the chain rule to compute $\partial E/\partial w_{ji}^h$, and then apply the gradient descent algorithm to obtain
\begin{equation}
\frac{\partial E}{\partial w_{ji}^h}=\frac{\partial E}{\partial z_j}\frac{\partial z_j}{\partial net_j}\frac{\partial net_j}{\partial w_{ji}^h},
\label{eq:chrutoewh}
\end{equation}
where
\begin{equation}
\frac{\partial net_j}{\partial w_{ji}^h}=x_i,
\label{eq:lneww}
\end{equation}
\begin{equation}
\frac{\partial z_j}{\partial net_j}=f^{'}_h(net_j).
\label{eq:midtwn}
\end{equation}
Finally, using~\eqref{eq:net_s_defn} gives
\begin{eqnarray}
\frac{\partial E}{\partial z_j}&=&\frac{\partial}{\partial z_j}\left(\frac{1}{2}\sum_{s=1}^m \left( d_s-f_o\left( net_s\right)\right)^2 \right)\nonumber\\
&=& -\sum_{s=1}^m \left( d_s-f_o\left( net_s\right)\right)\frac{\partial f_o\left( net_s\right)}{\partial z_j}\nonumber\\
&=& -\sum_{s=1}^m \left( d_s-y_s\right)f_o^{'}\left(net_s \right)w_{sj}^o.
\label{eq:pewzjbi}
\end{eqnarray}
Combining~(\ref{eq:chrutoewh})--(\ref{eq:pewzjbi}), we obtain rules for updating the hidden weights,
\begin{eqnarray}
w_{ji}^{h}(k+1) &=& w_{ji}^{h}(k)-\alpha_h \frac{\partial E\left(\b w(k) \right)}{\partial w_{ji}^{h}(k)}\nonumber\\
&=& w_{ji}^{h}(k)+\alpha_h\left(\sum_{s=1}^m \left( d_s-y_s\right)f_o^{'}\left(net_s \right)w_{sj}^{o}(k)\right)f_h^{'}(net_j)x_i.
\label{eq:rfhiwfo}
\end{eqnarray}

A commonly used activation function\index{activation function!sigmoid} is the sigmoid,\index{sigmoid!as an activation function}
\begin{equation}
f(v)=\frac{1}{1+e^{-v}}.
\label{eq:sigmoid_act}
\end{equation}
It is easy to check that the activation function~(\ref{eq:sigmoid_act}) has the useful property,
\begin{equation}
\frac{d f(v)}{d v}=f(v)(1-f(v)).
\label{eq:useproofsig}
\end{equation}
We summarize the above derivations in the following algorithm.\index{backpropagation!training algorithm}
\begin{algorithm}\rm
\begin{center}
{\bf Backpropagation Training}
\end{center}
Given a set of training pairs $\left\{\left(\b x^{(1)},\b d^{(1)}\right),\ldots,\left(\b x^{(i)},\b d^{(i)}\right),\ldots,\left(\b x^{(p)},\b d^{(p)}\right)  \right\}$.
\begin{enumerate}
\item {\bf Initialize Weights and Thresholds}\\
Initialize the weights and thresholds to small random values.
\item {\bf Present New Input Vector and Desired Output}\\
Present the new input vector, $\b x^{(i)}$ and the corresponding desired output $\b d^{(i)}$. Calculate actual outputs.
\item {\bf Calculate Error Gradients}
\item {\bf Adapt Weight Vectors and Thresholds}
\begin{eqnarray*}
\b w(k+1)&=&\b w(k)-\alpha\nabla_{\b w} E\left(\b w(k),\b t(k) \right),\\
\b t(k+1) &=& \b t(k)-\alpha\nabla_{\b t} E\left(\b w(k), \b t(k) \right).
\end{eqnarray*}
\item {\bf Repeat by Returning to Step 2}
\end{enumerate}
\end{algorithm}\markex

\begin{example}\label{ex:backproxor}\rm
We consider the two-input XOR problem.\index{XOR problem!two-input} We showed, in Example~\ref{ex:tlu_xor}, that this logic function cannot be implemented using a single TLU. Here, we use a two layer neural network with two nodes in the hidden layer and one node in the output layer as shown in Figure~\ref{fig:xor_example}.
\begin{figure}
\centering
\includegraphics[scale=0.45]{xor_example.eps}
\caption{A two-layer feedforward neural network for solving the XOR problem.}
\label{fig:xor_example}
\end{figure}

In this example, we use the unipolar representation\index{unipolar representation!of XOR function}\index{XOR function!unipolar representation of} of the two-input XOR function, that is, the elements of the training pairs are from the set $\{0, 1\}$. Therefore, the input vectors and the corresponding desired outputs have the form
\begin{eqnarray*}
\b x^{(1)}&=&\left[\begin{array}{c}
0\\
0
\end{array}\right],\; d_1=0;\quad \b x^{(2)}=\left[\begin{array}{c}
1\\
1
\end{array}\right],\; d_2=0;\\
\b x^{(3)}&=&\left[\begin{array}{c}
0\\
1
\end{array}\right],\; d_3=1;\quad \b x^{(4)}=\left[\begin{array}{c}
1\\
0
\end{array}\right],\; d_4=1.
\end{eqnarray*}
The error function for this example takes the form
\[
E=\frac{1}{2}(d_i-y_i)^2.
\]
The activation function\index{activation function!sigmoid} for all three neurons is the sigmoid given by~(\ref{eq:sigmoid_act}). The outputs of the hidden neurons are
\begin{eqnarray}
z_1&=&f(w_{11}^h x^{(i)}_1+w_{12}^hx^{(i)}_2-t_1)\label{eq:outfhidd}\\
z_2&=&f(w_{21}^h x^{(i)}_1+w_{22}^hx^{(i)}_2-t_2)\label{eq:outsechidd}
\end{eqnarray}
The network output is
\begin{equation}
y=f(w_{11}^o z_1+w_{12}^o z_2-t_o)
\label{eq:net_outxor}
\end{equation}
The threshold parameters are treated as learning weights.\index{threshold parameters!as learning weights} Thus, the above simple network requires updating of nine weights, where, in particular,
\begin{eqnarray*}
\frac{\partial E}{\partial t_1}&=&(d_i-y)y'w_{11}^o z_1'=(d_i-y)y(1-y)w_{11}^oz_1(1-z_1)\\
\frac{\partial E}{\partial t_2}&=&(d_i-y)y'w_{12}^o z_2'=(d_i-y)y(1-y)w_{12}^oz_2(1-z_2)\\
\frac{\partial E}{\partial t_o}&=&(d_i-y)y'=(d_i-y)y(1-y).
\end{eqnarray*}
In calculating the above, we used~(\ref{eq:useproofsig}). In our simulation the weights of the network were initialized to small random values in the interval $[-0.1, 0.1]$. The four patterns were presented in a mixed order. We used a fixed learning rate\index{learning rate} of 0.75 for all nine weights. After 50,000 iterations we obtained the following weight values:
\[
\b W^h=\left[\begin{array}{cc}
-4.9288 & -4.9334\\
-6.7240 & -6.7526
\end{array}\right],\quad \b W^o=\left[\begin{array}{cc}
10.2447 & -10.7258
\end{array}\right],
\]
and $t_1=-7.3398$, $t_2=-2.7853$, and $t_o=4.8701$. In Table~\ref{tab:xor_res} we show the desired and actual responses of the trained network.\index{trained network}
\begin{table}
\caption{The desired and actual responses of the trained network in Example~\ref{ex:backproxor}.}
\label{tab:xor_res}
\begin{center}
\begin{tabular}{cccc}\hline
$x_1$ & $x_2$ & $d_i$ & $y_i$\\\hline
0 & 0& $0$ & 0.0087\\
1 & 1 & $0$ & 0.0162\\
0 & 1 & $1$ & 0.9870\\
1 & 0 & $1$ & 0.9870\\ \hline
\end{tabular}
\end{center}
\end{table}
In Figure~\ref{fig:xor_bp} we show a plot of the function that the trained network implements.
\begin{figure}[h]
\centerline{\psfig{figure=xor_bp.eps,height=3.5in}}
\caption{The function implemented by the trained neural network to solve the two-input XOR problem in Example~\ref{ex:backproxor}.}
\label{fig:xor_bp}
\end{figure}
\end{example}
\markex

\section{Neural Fuzzy Identifier}
\label{subsec:nfuzid}
In this section, we describe a method that can be used to construct identifiers of dynamical systems\index{identifiers of dynamical systems!neural-network-based} that in turn could be employed to construct neural-network-based fuzzy logic control strategies. The idea behind the method is to apply the backpropagation algorithm to a fuzzy logic system.\index{backpropagation!applied to a fuzzy logic system}

It follows from our discussion in Section~\ref{sec:sam} that a typical fuzzy logic system\index{fuzzy logic system!as nonlinear system} can be thought of as a nonlinear system described by
\begin{equation}
f(x_1,\ldots,x_n)=f(\b x)=\frac{\sum_{l=1}^M v_l\left( \prod_{i=1}^n \mu_{il}(x_i)\right)}{\sum_{l=1}^M  \prod_{i=1}^n \mu_{il}(x_i)},
\label{eq:flcasnlsys}
\end{equation}
where, for example,
\[
\mu_{il}=\exp\left\{ - \frac{\left(x_i-m_{il}\right)^2}{2\sigma_{il}^2}\right\},
\]
$M$ is the number of IF-THEN fuzzy rules,\index{IF-THEN fuzzy rules} $\b x=\left[\begin{array}{cccc}
x_1&x_2&\cdots&x_n
\end{array}\right]^{\top}$, and $m_{il}$ and $\sigma_{il}$ are the parameters of the Gaussian membership functions,\index{Gaussian membership functions}\index{membership functions!Gaussian } while $v_l$ is a weight parameter. These parameters are determined using the backpropagation algorithm\index{backpropagation!algorithm} discussed in the previous section. To facilitate our discussion, we represent~(\ref{eq:flcasnlsys}), following Wang~\cite[page~169]{lxwang97}, as a feedforward network shown in Figure~\ref{fig:flcfeed}.
\begin{figure}
%\vspace{3.5in}
\centerline{\psfig{figure=zakflow.eps,height=3.5in}}
%\hskip1in\special{eps:c:/mybook/zakflow.eps y=3.5in}
\caption{Fuzzy logic system~(\ref{eq:flcasnlsys}) represented as a feedforward network.}
\label{fig:flcfeed}
\end{figure}
The variables that characterize the feedforward network\index{feedforward network} depicted in Figure~\ref{fig:flcfeed} satisfy the following relations:
\begin{equation}
\left.\begin{array}{ccl}
z_l&=&\prod_{i=1}^n\exp\left\{ - \frac{\left(x_i-m_{il}\right)^2}{2\sigma_{il}^2}\right\},\quad l=1,2\ldots,M\\
b&=& \sum_{l=1}^M z_l\\
a&=& \sum_{l=1}^M v_l z_l\\
f&=& a/b
\end{array}\right\}
\label{eq:relfeednet}
\end{equation}
We then use this feedforward network as a component of an identifier of dynamical systems. The identification scheme\index{identification!scheme} involving~(\ref{eq:flcasnlsys}) is depicted in Figure~\ref{fig:nufu}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{nufu_id.eps}
%\vspace{3.5in}
%\centerbmp{4in}{3in}{c:/mybook/nufu_id.bmp}
%\hskip1in\special{bmp:c:/mybook/nufu_id.bmp x=4.67in y=3.5in}
\caption{Neural fuzzy identification scheme of a nonlinear dynamical system.}
\label{fig:nufu}
\end{figure}
The adaptation algorithm for the identifier parameters\index{adaptation algorithm!for the identifier parameters} is derived using the backpropagation algorithm. The cost function\index{cost function} to be minimized is
\[
E=\frac{1}{2}\left(y-\hat{y}\right)^2=\frac{1}{2}\left(y-f\right)^2,
\]
where $y$ is the output of the identified plant and $f$ is the output of the identifier. We use a gradient algorithm to minimize $E$ over the space of parameters
\[
\mb{\bm$w$}=\left\{v_l, \: \sigma _{il},\: m_{il} \right\}.
\]
The gradient algorithm that we use has the form
\begin{equation}
\mb{\bm$w$}(k+1)=\mb{\bm$w$}(k)-\alpha \nabla E\left(\mb{\bm$w$}(k)\right),
\label{eq:fixstegradal}
\end{equation}
where $\alpha$ is a fixed step size. Note that~(\ref{eq:fixstegradal}) is just a simple descent gradient algorithm\index{descent gradient!algorithm} that we used in our discussion of the backpropagation algorithm.\index{backpropagation!algorithm}\index{backpropagation algorithm!descent gradient algorithm}

In order to be able to use the above algorithm, we need to compute partial derivatives of $E$ with respect to each component of $\mb{\bm$w$}$. For this we first fix indices $i$, and $l$. We next compute the partial derivatives of $E$ with respect to $v_j$. Then, the parameters $v_j$ will be adapted accordingly as
\[
v_j(k+1)=v_j(k)-\alpha \frac{\partial E\left(\b w(k)\right)}{\partial v_j}.
\]
We note that $E$ depends upon $v_j$ only through $a$, where $f=a/b$. Hence, using the chain rule, we obtain
\begin{eqnarray*}
\frac{\partial E}{\partial v_j}\left(\mb{\bm$w$}\right)&=&\frac{\partial}{\partial v_j}\left(\frac{1}{2}\left(y-f\right)^2\right)\\
&=&\frac{\partial}{\partial v_j}\left(\frac{1}{2}\left(y-a/b\right)^2\right)\\
&=&-\frac{\left(y-f\right)}{b}\frac{\partial a}{\partial v_j}\\
&=& \frac{f-y}{b}z_j.
\end{eqnarray*}
We next compute the partial derivative of $E$ with respect to $m_{st}$. Observe that $E$ depends upon $m_{st}$ only through $z_t$. Using the chain rule again gives
\begin{eqnarray*}
\frac{\partial E}{\partial m_{st}}\left(\mb{\bm$w$}\right)&=&-\left(y-f\right)\frac{\partial f}{\partial z_t} \frac{\partial z_t}{\partial m_{st}}\\
&=&-\left(y-f\right)\left( \frac{v_t b-a}{b^2}\right)z_t\frac{\left( x_s-m_{st}\right)}{\sigma_{st}^2}\\
&=&\left(f-y\right)\left( \frac{v_t-f}{b}\right)z_t\frac{\left( x_s-m_{st}\right)}{\sigma_{st}^2}.
\end{eqnarray*}
Hence, the adaptation formula for $m_{st}$ takes the form
\[
m_{st}(k+1)=m_{st}(k)-\alpha\left(f-y\right)\left( \frac{v_t-f}{b}\right)z_t\frac{\left( x_s-m_{st}(k)\right)}{\sigma_{st}^2(k)}.
\]
To find the partial derivative of $E$ with respect to $\sigma_{st}$, we note that $E$ depends upon $\sigma_{st}$ only through $z_t$. Applying the chain rule gives
\begin{eqnarray*}
\frac{\partial E}{\partial \sigma_{st}}\left(\mb{\bm$w$}\right)&=&-\left(y-f\right)\frac{\partial f}{\partial z_t} \frac{\partial z_t}{\partial \sigma_{st}}\\
&=&\left(f-y\right)\left( \frac{v_t-f}{b}\right)z_t\frac{\left( x_i-m_{st}\right)^2}{\sigma_{st}^3}.
\end{eqnarray*}
We are now ready to formulate an algorithm for updating on-line the weights of the neural fuzzy identifier.
\begin{algorithm}\rm
\begin{center}
{\bf On-line Identification}
\end{center}
\begin{enumerate}
\item {\bf Initialize the identifier parameters}\\
Set the values of $v_l$, $m_{il}$, and $\sigma_{il}$ to small random values.
\item {\bf Calculate the output of the identifier at the time instant $k$}\\
Compute
\begin{eqnarray*}
z_l&=&\prod_{i=1}^n\exp\left( - \frac{\left(x_i-m_{il}\right)^2}{2\sigma_{il}^2} \right),\\
b&=& \sum_{l=1}^M z_l,\\
a&=& \sum_{l=1}^M v_l z_l,\\
f&=& a/b;
\end{eqnarray*}
\item {\bf Calculate the identification error}
\[
e(k)=y(k)-f(k);
\]
\item {\bf Adapt the identifier parameters}
\begin{eqnarray*}
v_l(k+1)&=&v_l(k)+\alpha\frac{e(k)}{b(k)}z_l(k),\\
m_{il}(k+1) &=& m_{il}(k)+\alpha\frac{e(k)}{b(k)}\left( v_l(k)-f(k)\right)z_l(k)\frac{\left(x_i(k)-m_{il}(k)\right)}{\sigma_{il}^2(k)},\\
\sigma_{il}(k+1) &=& \sigma_{il}(k)+\alpha\frac{e(k)}{b(k)}\left( v_l(k)-f(k)\right)z_l(k)\frac{\left(x_i(k)-m_{il}(k)\right)^2}{\sigma_{il}^3(k)},\\
k&=&k+1;
\end{eqnarray*}
\item {\bf Repeat by Returning to Step 2}
\end{enumerate}
\end{algorithm}\markex

\begin{example}\label{ex:nufuzid}\rm
The plant to be identified\index{plant!to be identified} is the same as that considered by Wang~\cite[page~175]{lxwang97}, and is described by
\[
y(k+1)=0.3y(k)+0.6y(k-1)+g(u(k)),
\]
where the ``unknown'' function $g(u)$ has the form
\[
g(u)=0.6\sin(\pi u)+0.3\sin(3\pi u)+0.1\sin(5 \pi u),
\]
and $u(k)=\sin(2\pi k/200)$. The identifier has the structure shown in Figure~\ref{fig:nufu} and is described by
\[
\hat{y}(k+1)=0.3 y(k)+0.6y(k-1)+f(u(k)),
\]
where $f$ has the form of~(\ref{eq:flcasnlsys}) with $M=10$. We chose $\alpha=0.5$, $y(0)=y(1)=0$, $f(1)=f(2)=-3$, and we started with $u(1)=0$. The identifier's parameters were initialized in a random fashion using MATLAB's command \verb!rand!.\index{MATLAB's command!rand}
\begin{figure}
%\vspace{3.5in}
\centerline{\psfig{figure=nufuzid.eps,height=3.5in}}
%\hskip1in\special{eps:c:/fuzzy/nufuzid.eps y=3.5in}
\caption{Plots of the outputs of the plant and its identifier during the learning phase in Example~\ref{ex:nufuzid}.}
\label{fig:nufuzid}
\end{figure}
After the identifier was trained, we tested its prediction power\index{prediction power!of neural identifier} using the same input signal as before and the initial conditions $y(0)=y(1)=0$, $f(1)=f(2)=-3$, and $u(1)=0$. The results of this experiment are shown in Figure~\ref{fig:nfu_tst}.
\begin{figure}
%\vspace{3.5in}
\centerline{\psfig{figure=nfu_tst.eps,height=3.5in}}
%\hskip1in\special{eps:c:/fuzzy/nfu_tst.eps y=3.5in}
\caption{Plots of the outputs of the plant and its identifier during the testing phase in Example~\ref{ex:nufuzid}.}
\label{fig:nfu_tst}
\end{figure}
\end{example}\markex

\section{Radial-Basis Function (RBF) Networks}
\label{sec:rbfnets}
In this section, we present radial-basis function~(RBF) neural networks\index{radial-basis function (RBF)}\index{RBF|see{radial-basis function}}\index{radial-basis function!neural networks} and their applications to interpolation and then to identification of unknown nonlinear dynamical systems using only the systems' inputs and outputs data.\index{identification!of unknown systems} We first show how an RBF network\index{RBF network!for interpolation} can be used for interpolation purposes. Then, we concern ourselves with the problem of identification of dynamical systems. Our exposition of using RBF neural networks to identify dynamical systems is based on the paper by Liu, Kadirkamanathan, and Billings~\cite{liu96}. The RBF neural identifier\index{RBF neural identifier!an adaptive system} is an adaptive system for which the learning algorithm\index{learning algorithm!using Lyapunov techniques} is developed using Lyapunov techniques. The algorithm can add new nodes to the neural identifier to obtain the desired identification accuracy.

\subsection{Interpolation Using RBF Networks}
A diagram of a radial-basis function~(RBF) network\index{RBF network} is shown in Figure~\ref{fig:grbf_ip}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{grbf_ip.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/grbf_ip.bmp}
%\hskip1in \special{bmp:c:/mybook/grbf_ip.bmp x=4in y=3in}
\caption{A two-input radial-basis function~(RBF) network.}
\label{fig:grbf_ip}
\end{figure}
The functions $\phi\left(\left\|\b x-\b x^{(i)} \right\| \right)$, $i=1,2,\ldots,N$, are called the radial-basis functions, where $\|\cdot \|$ denotes a $p$-norm\index{p-norm!for example, Euclidean 2-norm} that is usually the Euclidean 2-norm. The real parameters $w_i$, $i=1,2, \ldots, N$, are the weights. A common choice for $\phi$ is
\begin{equation}
\phi\left(\left\|\b x-\b x^{(i)} \right\| \right)=\exp\left(-\frac{\left\| \b x-\b x^{(i)}\right\|^2}{2\sigma_i^2} \right),
\label{eq:grbfunct}
\end{equation}
where $\b x^{(i)}$ is the center of the basis function\index{basis function!center of} and $\sigma_i$ is sometimes called the radius of the radial basis function.\index{radial basis function!radius of}\index{basis function!radial} The function given by~(\ref{eq:grbfunct}) is called the Gaussian radial basis function.\index{Gaussian radial basis function}

The RBF networks are used to perform interpolation.\index{RBF networks!perform interpolation} The interpolation problem\index{interpolation problem!formulation} can be formulated as: For a given set of $N$ different points,
\[
\left\{ \b x^{(i)}:\b x^{(i)}\in \mb{\Bb R}^n,\;i=1,2,\ldots, N\right\}
\]
and a corresponding set of $N$ real numbers
\[
\left\{ d_i:d_i\in \mb{\Bb R},\;i=1,2,\ldots, N\right\}
\]
find a function
\[
F:\mb{\Bb R}^n \to \mb{\Bb R}
\]
such that
\begin{equation}
F\left( \b x^{(i)}\right)=d_i,\;i=1,2,\ldots, N.
\label{eq:interp_con}
\end{equation}
We attempt to solve the interpolation problem using the RBF network.\index{interpolation!using RBF network} The function that the RBF network implements is
\begin{equation}
F(\b x)=\sum_{j=1}^{N}w_j \phi\left(\left\|\b x-\b x^{(j)} \right\| \right).
\label{eq:rbffunimp}
\end{equation}
Substituting the interpolation conditions given by~(\ref{eq:interp_con}) into~(\ref{eq:rbffunimp}) gives
\begin{equation}
\left[\begin{array}{cccc}
\phi_{11} & \phi_{12} & \cdots & \phi_{1N}\\
\phi_{21} & \phi_{22} & \cdots & \phi_{2N}\\
\vdots & \vdots & \ddots & \vdots\\
\phi_{N1} & \phi_{N2} & \cdots & \phi_{NN}
\end{array}\right]\left[\begin{array}{c}
w_1\\
w_2\\
\vdots\\
w_N
\end{array}\right]=\left[\begin{array}{c}
d_1\\
d_2\\
\vdots\\
d_N
\end{array}\right],
\label{eq:matvecinter}
\end{equation}
where $\phi_{ij}=\phi\left(\left\|\b x^{(i)}-\b x^{(j)} \right\| \right)$. Let $\b w=\left[\begin{array}{cccc}
w_1& w_2 &\cdots &w_N
\end{array}\right]$, $\b d=\left[\begin{array}{cccc}
d_1& d_2 &\cdots &d_N
\end{array}\right]$, and
\[
\b{\Phi}=\left[\begin{array}{cccc}
\phi_{11} & \phi_{12} & \cdots & \phi_{1N}\\
\phi_{21} & \phi_{22} & \cdots & \phi_{2N}\\
\vdots & \vdots & \ddots & \vdots\\
\phi_{N1} & \phi_{N2} & \cdots & \phi_{NN}
\end{array}\right].
\]
Using the above notation, we represent~(\ref{eq:matvecinter}) as
\begin{equation}
\b{\Phi}\b w=\b d.
\label{eq:matvecip}
\end{equation}
The matrix $\b{\Phi}$ is called the {\em interpolation matrix\/}.\index{interpolation matrix} If the interpolation matrix is nonsingular, then we can solve~(\ref{eq:matvecip}) for $\b w$ to obtain
\begin{equation}
\b w=\b{\Phi}^{-1}\b d.
\label{eq:weqqinvd}
\end{equation}
\begin{example}\label{ex:iterusrgf}\rm
We construct an RBF network\index{RBF network} that interpolates data given in Table~\ref{tab:interdata}.
\begin{table}
\caption{Interpolation data for Example~\ref{ex:iterusrgf}.}
\begin{center}
\begin{tabular}{cc}
$x^{(i)}$ & $d_i$\\\hline
0 & 0.8\\
6 & 0.2\\
7 & 0.4\\
8 & 0.6\\
10 & 0.3\\\hline
\end{tabular}
\end{center}
\label{tab:interdata}
\end{table}
We have, $\phi_{ij}=\exp\left(-\frac{(x^{(i)}-x^{(j)})^2}{2\sigma_i^2}\right)$. We set $\sigma_i=1$, for $i=1,\ldots,5$, and evaluate~(\ref{eq:weqqinvd}) to obtain
\[
\b w=\left[\begin{array}{ccccc}
0.8000& 0.1459& -0.0390& 0.5737& 0.2227
\end{array}\right]^{\top}.
\]
The interpolation result is shown in Figure~\ref{fig:grbf_eip}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{grbf_eip.eps}
%\vspace{3.5in}
%\centerbmp{4.67in}{3.5in}{c:/mybook/grbf_eip.bmp}
%\hskip.9in\special{bmp:c:/mybook/grbf_eip.bmp x=4.67in y=3.5in}
\caption{Results of interpolating data from Table~\ref{tab:interdata} using radial-basis functions.}
\label{fig:grbf_eip}
\end{figure}
In Figure~(\ref{fig:grbf_eip})(a), we used a solid line to show the result of interpolating using the Gaussian radial-radial basis functions,\index{Gaussian radial-radial basis functions} while a broken line indicates the result of the straight line interpolation. The plots of the basis functions that contribute to the interpolating curve are shown in Figure~(\ref{fig:grbf_eip})(b).
\end{example}
\markex
If the number of radial basis functions was fixed and the rows of the matrix $\b{\Phi}$ and the elements of the vector $\b d$ were obtained sequentially, we could use the recursive least squares algorithm, described in Section~\ref{sec:rls}, to compute the least squares estimate of $\b w$.

\subsection{Identification of a Single-Input Single-State System}
\label{subsec:sistsiin}
Here, we discuss the identification of a single-input single-state dynamical system using neural network\index{identification!using neural network}\index{identification!of a single-input single-state system}. We consider a dynamical system modeled by
\begin{equation}
\dot{x}(t)=f(x(t),u(t)),\quad x(0)=x_0,
\label{eq:sissdynsid}
\end{equation}
where $f:\mb{\Bb R}\times \mb{\Bb R} \to \mb{\Bb R}$ is an unknown, possibly nonlinear, function of $x$ and $u$. We will construct a neural network to identify the function $f$. Adding and subtracting $ax(t)$, where $a>0$ is a design parameter, to the right-hand side of~(\ref{eq:sissdynsid}) gives
\begin{eqnarray}
\dot{x}(t) &=& -ax(t)+ax(t)+f(x(t),u(t))\nonumber\\
&=& -ax(t) +g(x(t),u(t)),
\label{eq:altrepfodyid}
\end{eqnarray}
where $g(x(t),u(t))=ax(t)+f(x(t),u(t))$. The reason for the above manipulations will become apparent when we discuss learning rules\index{neural identifier!learning rules}\index{learning rules!of neural identifier} of the neural network identifier of the dynamical system~(\ref{eq:sissdynsid}). This identifier is described by
\begin{equation}
\dot{\hat{x}}(t)=-a\hat{x}(t)+\hat{g}(x(t),u(t),\b m_k,\sigma_k),\quad \hat{x}(0)=\hat{x}_0,
\label{eq:nesysid}
\end{equation}
where $\hat{x}$ is the state of the identification network, $\b m_k$ and $ \sigma_k$ are adjustable parameters of the identifier.\index{neural identifier!adjustable parameters}\index{adjustable parameters!of neural identifier} The function $\hat{g}$ has the form,
\begin{equation}
\hat{g}(x(t),u(t),\b m_k,\sigma_k)=\sum_{k=1}^K w_k \phi(x,u,\b m_k,\sigma_k),
\label{eq:ghatnnout}
\end{equation}
where for $k=1,2,\dots,K$, we have $\b m_k=\left[\begin{array}{cc}
m_{xk}& m_{uk}
\end{array}\right]$, and
\[
\phi(x(t),u(t),\b m_k,\sigma_k)=\exp\left(-\frac{1}{2\sigma_k^2}\left( \left( x(t)-m_{xk}\right)^2+\left(u(t)-m_{uk} \right)^2\right) \right).
\]
Thus, we postulate that the function $\hat{g}$ be generated using an RBF network.\index{RBF network!neural identifier} 


Suppose that the basis functions are given to us. The problem of choosing basis functions is discussed later. Suppose also that we know the ``optimal'' weights $w_k^*$, $k=1,2,\ldots, K$, that result in the least identification error\index{least identification error} between $g$ in~(\ref{eq:altrepfodyid}) and $\hat{g}$ in~(\ref{eq:nesysid}), where the identification error is defined as
\begin{equation}
e(t)= g(x(t),u(t))- \sum_{k=1}^K w_k^* \phi(x(t),u(t),\b m_k,\sigma_k).
\label{eq:errbgandgha}
\end{equation}
Using the above notation, we represent the dynamical system model~(\ref{eq:altrepfodyid}) in the form
\begin{equation}
\dot{x}(t)=-ax(t)+\sum_{k=1}^K w_k^* \phi(x(t),u(t),\b m_k,\sigma_k)+e(t).
\label{eq:morepuserr}
\end{equation}
The neural identifier,\index{neural identifier!representation of} on the other hand, can be represented as
\begin{equation}
\dot{\hat{x}}(t)=-a\hat{x}(t)+\sum_{k=1}^K w_k \phi(x(t),u(t),\b m_k,\sigma_k).
\label{eq:idneuserr}
\end{equation}
A schematic of the above identification architecture\index{neural identifier!architecture of} is depicted in Figure~\ref{fig:grbf_id}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{grbf_id.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/grbf_id.bmp}
%\hskip1.2in \special{bmp:c:/mybook/grbf_id.bmp x=4in y=3in}
\caption{A configuration for the radial-basis function neural identifier.}
\label{fig:grbf_id}
\end{figure}

\subsection{Learning Algorithm for the RBF Identifier}
We assume  that both the state $x$ and the input $u$ of the dynamical system model~(\ref{eq:sissdynsid}), and hence that of~(\ref{eq:altrepfodyid}), are bounded. We define the state estimation error\index{state estimation error}\index{error!of state estimation} to be
\[
e_x=x-\hat{x},
\]
while the weights' estimation errors are
\[
\xi_k=w^*_k-w_k,\quad k=1,\ldots,K.
\]
We assume that the modeling error,\index{modeling error} $e$, defined in~(\ref{eq:errbgandgha}), is bounded by a finite constant $e_K$, that is, $\left|e\right|\le e_K$. However, $e_K$ is unknown to us. Next, we use~(\ref{eq:morepuserr}) and~(\ref{eq:idneuserr}) to form the differential equation describing the state error dynamics. For simplicity of notation, we omit the argument $t$ of the time-varying functions. We have
\begin{eqnarray}
\dot{e}_x &=& \dot{x}-\dot{\hat{x}}\nonumber\\
&=& -ax+\sum_{k=1}^K w_k^* \phi(x,u,\b m_k,\sigma_k)+e+a\hat{x}-\sum_{k=1}^K w_k \phi(x,u,\b m_k,\sigma_k)\nonumber\\
&=& -ae_x+\sum_{k=1}^K \xi_k \phi(x,u,\b m_k,\sigma_k)+  e.
\label{eq:errdyneq}
\end{eqnarray}
We now use the second method of Lyapunov to devise a learning law for the RBF identifier that would yield a stable identification scheme in the sense that the state estimation error $e_x$ decreases as the time increases and the weights do not diverge. For this purpose, we use the following Lyapunov function candidate,
\begin{equation}
V(e_x,\b{\xi})=\frac{1}{2}\left( e_x^2+\frac{1}{\alpha}\sum_{k=1}^K\xi_k^2\right),
\label{eq:lyafuncarbf}
\end{equation}
where $\alpha>0$ is a design parameter and $\b{\xi}=\left[\begin{array}{cccc}
\xi_1 & \xi_2&\cdots & \xi_K
\end{array}\right]^{\top}$.
The time derivative of $V$ evaluated on the trajectories of~(\ref{eq:errdyneq}) is
\begin{eqnarray}
\dot{V}&=& e_x\dot{e}_x+\frac{1}{\alpha}\sum_{k=1}^K \xi_k\dot{\xi}_k\nonumber\\
&=&e_x\left( -ae_x+\sum_{k=1}^K \xi_k \phi(x,u,\b m_k,\sigma_k)+  e\right)+\frac{1}{\alpha}\sum_{k=1}^K \xi_k\dot{\xi}_k\nonumber\\
&=& -ae_x^2+\frac{1}{\alpha}\sum_{k=1}^K \xi_k\left(\alpha e_x\phi(x,u,\b m_k,\sigma_k)+\dot{\xi}_k\right)+e_x e.
\label{eq:lyaderofsterid}
\end{eqnarray}
Let
\begin{equation}
\dot{\xi}_k=-\alpha e_x\phi(x,u,\b m_k,\sigma_k).
\label{eq:updfoxi}
\end{equation}
Then~(\ref{eq:lyaderofsterid}) simplifies to
\begin{equation}
\dot{V}=-ae_x^2+e_x e.
\label{eq:vdotsimpnid}
\end{equation}
Because $\dot{w}_k=-\dot{\xi}_k$, we can represent~(\ref{eq:updfoxi}) as
\begin{equation}
\fbox{$\displaystyle \dot{w}_k=\alpha e_x\phi(x,u,\b m_k,\sigma_k) $}
\label{eq:weupdnid}
\end{equation}
We further evaluate~(\ref{eq:vdotsimpnid}) to obtain
\begin{eqnarray}
\dot{V} &\le & -ae_x^2+\left\arrowvert e_x\right\arrowvert e_K\nonumber\\
&=& -a\left\arrowvert e_x\right\arrowvert\left( \left\arrowvert e_x\right\arrowvert-\frac{e_K}{a}\right)
\label{eq:vdotevalwek}
\end{eqnarray}
It follows from~(\ref{eq:vdotsimpnid}) that if the state estimation error is zero, that is, $e_K=0$, then
\[
\dot{V}=-ae_x^2.
\]
Thus in this case, $\dot{V}$ is negative semi-definite in the $\left[\begin{array}{cc}
e_x & \b \xi ^{\top}
\end{array}\right]^{\top}$-space, and hence the origin of the $\left[\begin{array}{cc}
e_x & \b \xi ^{\top}
\end{array}\right]^{\top}$-space is stable, which guarantees boundedness of $e_x$ and the weights $w_k$s. We will  show, using the Lyapunov-like lemma stated on page~\pageref{lemma:lyapunov_like}, that when $e=0$, then $e_x(t)\to 0$ as $t\to \infty$. We have to check if the assumptions of the Lyapunov-like lemma are satisfied. First, observe that $V$ is bounded below and positive semi-definite. Next, because $e_x$, $\xi_k$ as well as $\phi(x,u,\b m_k,\sigma_k)$ are bounded, the derivative of the state estimation error, given by~(\ref{eq:errdyneq}) for $e=0$, is bounded as well, which implies that
\[
\ddot{V}=-2 a e_x\dot{e}_x
\]
is bounded. This in turn implies that $\dot{V}$ is uniformly continuous. Hence, by the Lyapunov-like lemma, $e_x(t)\to 0$ as $t\to \infty$.


In the presence of the identification error, $e$, if $\left\arrowvert e_x\right\arrowvert <e_K/a$, then it is possible that $\dot{V}>0$, which implies that the weights $w_k$s may drift to infinity with time. To avoid this problem, Liu, Kadirkamanathan, and Billings~\cite{liu96} propose to set $\dot{w}_k=0$ if $\left\arrowvert e_x\right\arrowvert <e_K/a$. Also, because $e_K$ is unknown, a bound is imposed on the weights as well as on the error. The learning rule becomes
\begin{equation}
\fbox{$\displaystyle\dot{w}_k=\left\{\begin{array}{lcl}
\alpha e_x\phi(x,u,\b m_k,\sigma_k) & \mb{if}& \left\arrowvert e_x\right\arrowvert\ge e_0\quad\mb{and}\quad \left\Arrowvert \b w\right\Arrowvert \le M\\
0 & {}& \mb{otherwise}
\end{array}\right. $}
\label{eq:molearru}
\end{equation}
where the bounds $e_0$ and $M$ are design parameters, and $\b w=\left[\begin{array}{ccc}
w_1 & \cdots & w_K
\end{array}\right]^{\top}$. With the above modified learning law, whenever $\left\arrowvert e_x\right\arrowvert \ge e_0\ge e_K/a$, the Lyapunov derivative $\dot{V}$ is negative semi-definite. However, when $e_0\le\left\arrowvert e_x\right\arrowvert \le e_K/a$, the weights may increase with time because it is possible that $\dot{V}>0$. We impose an upper-bound on the weights. We also impose an upper bound on the estimation error, $e_x$. When any of the two bounds is reached, more basis functions may be needed to improve the approximation of $g$. New basis functions are added to reduce the error until $\dot{V}<0$. In what follows, we describe a procedure for adding radial-basis functions to the identification network.

\subsection{Growing RBF Network}
As we mentioned above, when $\dot{V}>0$, the identification scheme is no longer stable, but the learning rule~(\ref{eq:molearru}) guarantees that the adjustable parameters are bounded. Once the parameter adjustments can no longer reduce the estimation error and the bounds on the weights or on the estimation error are reached, additional basis functions in the RBF network should be added to reduce the modeling error. This continues until the approximation leads to the error reduction so that $\dot{V}<0$.

In the proposed identification scheme, we use the Gaussian radial-basis functions, whose areas of ``influence'' are determined by the functions' centers and radii (widths). As the function to be approximated moves further away from a particular basis function, the contribution of that basis function decreases. This is due to the exponential form of the Gaussian basis functions. As $|z|$ increases $\exp(-z^2)$ decreases. The maximum contribution of a particular basis function to approximating a given function occurs when the center of the basis function and the approximated function are at the same point because in this case we have $\exp(0)=1$. This heuristic argument leads to the following rule for determining the centers of new basis functions to be added at a time instance $t$,
\begin{equation}
\b m_{K+1}=\left[\begin{array}{cc}
m_{x,K+1} & m_{u,K+1}
\end{array}\right]^{\top}=
\left[\begin{array}{cc}
\delta_x\mb{round}\left( x(t)/\delta_x\right)&\delta_u\mb{round}\left( u(t)/\delta_u\right)
\end{array}\right]^{\top},
\label{eq:cenofnbfun}
\end{equation}
where the round operator rounds the elements on which it operates to the nearest integer. For example, if $\b x=\left[\begin{array}{ccc}
-1.8 & -0.3 & 3.7
\end{array}\right]$, then round$(\b x)=\left[\begin{array}{ccc}
-2 & 0 & 4
\end{array}\right]$. The parameters, $\delta_x$ and $\delta_u$ represent the spacing of the centers of the basis functions. We can think of having a grid over the input-state space where the spacing between the vertices of the grid are determined by the parameters $\delta_x$ and $\delta_u$. Each vertex of the grid represents the potential center of a basis function. The parameter $\sigma_k$ of the new basis function is determined from the relation,
\begin{equation}
\sigma_k=\kappa \left\Arrowvert \b m_k-\b m^*_k\right\Arrowvert,
\label{eq:radofbfun}
\end{equation}
where $\kappa>0$ is a design parameter and
\begin{equation}
\b m_k^*=\arg\min_{i=1,\ldots,K+1\atop i\ne k}\left\{ \left\Arrowvert \b m_k-\b m_i\right\Arrowvert\right\}.
\label{eq:defmkstar}
\end{equation}
Thus, $\b m_k^*$ is the nearest center of all basis functions to the center of the $k$-th basis function. With the above background, we can now present the rules for adding a new basis function to the RBF network. Initially, the network contains no basis functions. Then, once the following two conditions are satisfied, a new basis function is added. The first condition is:
\begin{equation}
\min_{k=1,\ldots,K}\left| x(t)-m_{xk}\right|>\delta_x/2\quad\mb{or}\quad \min_{k=1,\ldots,K}\left| u(t)-m_{uk}\right|>\delta_u/2.
\label{eq:cosfxau}
\end{equation}
The above will be satisfied if the state or input, at time $t$, is closer to an available ``unoccupied'' grid vertex than to any of the existing centers of the basis functions. The basis function at the new vertex would be the closest to the existing state or input and would have the greatest influence in the approximation of $g$. The second condition is
\begin{equation}
\left| x(t)-\hat{x}(t)\right|>e_{\max},
\label{eq:thcofoabf}
\end{equation}
where $e_{\max}$ is the maximum allowable state estimation error. If~(\ref{eq:thcofoabf}) is satisfied, a new basis function is added to reduce the state estimation error by reducing the modeling error.

\begin{example}\rm
The above technique for identification of dynamical systems using radial-basis function network is illustrated with the system,
\begin{equation}
\dot{x}=1.1\left(1-x+2xu-x^2 \right)\exp\left( -0.5x^2-u^2\right),\quad x(0)=0,
\label{eq:simsssiex}
\end{equation}
where $u=\cos(t/100)$. This model was adapted from~\cite{liu96}. There are many ways to simulate the above identification scheme. In our simulations, we discretized the above differential equation using the Euler method with the step size $dt=0.01$~sec. The initial weights and the initial state of the identifier were set to zero. We used the following parameter values: $e_0=0.001$, $e_{\max}=0.005$, $\alpha=5$, $\kappa=3$. In Figure~\ref{fig:grbf_xes}, plots of the true state and its estimate versus time are shown.
\begin{figure}
\centering
\includegraphics[scale=0.5]{grbf_xes.eps}
%\vspace{3.5in}
%\centerbmp{3.97in}{3.5in}{c:/ee675/grbf_xes.bmp}
%\hskip1.2in \special{bmp:c:/ee675/grbf_xes.bmp x=3.97in y=3.5in}
\caption{Plots of the true state and its estimate versus time. The estimated state ``hugs'' closely the actual state.}
\label{fig:grbf_xes}
\end{figure}
A plot of the state estimation error, $e_x$, between the true state and its estimate is shown in Figure~\ref{fig:grbf_ese}. Initially, the error exceeds the bound $e_{\max}$. Adding additional basis functions drives the error back down.
\begin{figure}
\centering
\includegraphics[scale=0.5]{grbf_ese.eps}
%\vspace{3.5in}
%\centerbmp{3.97in}{3.5in}{c:/ee675/grbf_ese.bmp}
%\hskip1.2in \special{bmp:c:/ee675/grbf_ese.bmp x=3.97in y=3.5in}
\caption{A plot of the error, $e_x$, between the true state and its estimate.}
\label{fig:grbf_ese}
\end{figure}
The plots shown in Figure~\ref{fig:grbf_uxs} depict the true state and its estimate versus the input $u$.
\begin{figure}
\centering
\includegraphics[scale=0.5]{grbf_uxs.eps}
%\vspace{3.5in}
%\centerbmp{3.97in}{3.5in}{c:/ee675/grbf_uxs.bmp}
%\hskip1.2in \special{bmp:c:/ee675/grbf_uxs.bmp x=3.97in y=3.5in}
\caption{A plot of the true state and its estimate versus the input $u$.}
\label{fig:grbf_uxs}
\end{figure}
Finally, in Figure~\ref{fig:grbf_xxd}, plots of the derivative of the true state versus the true state and the derivative of the estimated state versus the state estimate are shown.
\begin{figure}
\centering
\includegraphics[scale=0.5]{grbf_xxd.eps}
%\vspace{3.5in}
%\centerbmp{3.97in}{3.5in}{c:/ee675/grbf_xxd.bmp}
%\hskip1.2in \special{bmp:c:/ee675/grbf_xxd.bmp x=3.97in y=3.5in}
\caption{Plots of the derivative of the true state versus the true state and the derivative of the estimated state versus the state estimate.}
\label{fig:grbf_xxd}
\end{figure}
\end{example}
\subsection{Identification of Multivariable Systems}
We now extend the results of Subsection~\ref{subsec:sistsiin} to multi-input multi-state dynamical systems. In our derivation of the update law for the weights of the RBF network, we use a result from matrix algebra presented next. Recall the definition of the trace of a square matrix as given by~(\ref{eq:mat_trace}). If $\b{\Gamma}\in \mb{\Bb R}^{m\times n}$, then by~(\ref{eq:frob_trace}) its Frobenius norm satisfies
\[
\|\b{\Gamma}\|_F^2 = \sum_{i=1}^m\sum_{j=1}^n \gamma_{ij}^2 = \mb{trace}\left(\b{\Gamma}\b{\Gamma}^{\top} \right).
\]
If the elements, $\gamma_{ij}$, of $\b{\Gamma}$ are functions of time, then
\begin{equation}
\frac{d}{dt}\left( \mb{trace}\left(\b{\Gamma}\b{\Gamma}^{\top} \right)\right)=2\sum_{i=1}^m\sum_{j=1}^n \gamma_{ij}\dot{\gamma}_{ij}=2\:\mb{trace}\left( \b{\Gamma}\dot{\b {\Gamma}}^{\top}\right),
\label{eq:deroftrac}
\end{equation}
where $\dot{\b{\Gamma}}=\frac{d \b{\Gamma}}{dt}$. With the above background, we can proceed with the analysis of the RBF neural identifier for the multivariable systems.

The multivariable dynamical systems that we consider are modeled by
\begin{equation}
\dot{\b x}(t)=\b f(\b x(t),\b u(t)), \quad \b x(0)=\b x_0,
\label{eq:mimsdysymo}
\end{equation}
where $\b u(t)\in \mb{\Bb R}^m$, $\b x(t)\in \mb{\Bb R}^n$, and $\b f:\mb{\Bb R}^n\times \mb{\Bb R}^m \to \mb{\Bb R}^n$ is the function to be identified. Let $\b A\in \mb{\Bb R}^{n \times n}$ be a Hurwitz matrix, that is, the matrix whose eigenvalues are all in the open left-half plane. Because $\b A$ is a Hurwitz matrix, by the theorem of Lyapunov (Theorem~\ref{Lyapunov}), for any given real symmetric positive definite (r.s.p.d.) matrix ${\b Q}$ the solution ${\b P}$ to the continuous Lyapunov matrix equation,
\begin{equation}
\b A^{\top}\b P+\b P\b A=-\b Q
\label{eq:lyaeqfoagiq}
\end{equation}
is also (real symmetric) positive definite. We use this fact later when we develop a learning law for the RBF neural identifier.

We represent the model~(\ref{eq:mimsdysymo}) in an equivalent form,
\begin{equation}
\dot{\b x}(t)=\b A\b x(t)+\b g(\b x(t),\b u(t)), \quad \b x(0)=\b x_0,
\label{eq:eqreofmimsmo}
\end{equation}
where
\[
\b g(\b x(t),\b u(t))=-\b A\b x(t)+\b f(\b x(t),\b u(t)).
\]
The $k$-th basis function can be represented as
\[
\phi_k(\b x,\b u)=\exp\left( -\frac{1}{2\sigma_k^2}\left\Arrowvert \left[\begin{array}{c}
\b x\\
\b u
\end{array}\right]-\b m_k \right\Arrowvert^2\right).
\]
Let
\[
\b{\Phi}_K(\b x,\b u)=\left[\begin{array}{cccc}
\phi_1(\b x,\b u) & \phi_2(\b x,\b u)&\cdots&\phi_K(\b x,\b u)
\end{array}\right]^{\top}
\]
and
\[
\hat{\b g}(\b x,\b u,\b m_k,\sigma_k)=\b W_K \b{\Phi}_K,
\]
where $\b W_K$ is the $n\times K$ weight matrix and $\b m_k=\left[\begin{array}{cc}
\b m_{\b x k}^{\top} & \b m_{\b u k}^{\top}
\end{array}\right]^{\top}$ is the center of the $k$-th radial basis function in the $\left[\begin{array}{cc}
\b x^{\top} & \b u^{\top}
\end{array}\right]^{\top}$-space. The RBF neural identifier has the form,
\begin{equation}
\dot{\hat{\b x}}(t)=\b A\hat{\b x}(t)+\hat{\b g}(\b x,\b u,\b m_k,\sigma_k),\quad \hat{\b x}(0)=\b x_0.
\label{eq:mimsneuid}
\end{equation}
Let $\b W^*_K$ be the optimal weight matrix that minimizes the modeling error,
\[
\b e(t)={\b g}(\b x(t),\b u(t))-\b W_K^* \b{\Phi}_K(\b x(t),\b u(t)).
\]
Using the above definition, we represent the system model~(\ref{eq:eqreofmimsmo}) as
\begin{equation}
\dot{\b x}(t)=\b A\b x(t)+ \b W_K^* \b{\Phi}_K(\b x(t),\b u(t))+\b e(t).
\label{eq:sysmouswstbph}
\end{equation}
We define the state estimation error and the weight estimation error as
\[
\b e_{\b x}=\b x-\hat{\b x}\quad\mb{and}\quad \b{\Gamma}_K=\b W_K^*-\b W_K.
\]
Next, we form the differential equation describing the state estimation error dynamics,
\begin{eqnarray}
\dot{\b e}_{\b x} &=& \dot{\b x}-\dot{\hat{\b x}}\nonumber\\
&=& \b A\b x+ \b W_K^* \b{\Phi}_K(\b x,\b u)+\b e-\b A\hat{\b x}-\hat{\b g}(\b x,\b u,\b m_k,\sigma_k)\nonumber\\
&=&\b A\left(\b x-\hat{\b x}\right)+\b W_K^* \b{\Phi}_K(\b x,\b u)-\b W_K \b{\Phi}_K(\b x,\b u)+\b e\nonumber\\
&=& \b A\b e_{\b x}+\b{\Gamma}_K\b{\Phi}_K(\b x,\b u)+\b e.
\label{eq:mimserdys}
\end{eqnarray}
We employ the following Lyapunov function candidate,
\begin{equation}
V=V\left(\b e_{\b x},\b{\Gamma}_K\right)=\frac{1}{2}\left( \b e_{\b x}^{\top}\b P\b e_{\b x} +\frac{1}{\alpha}\mb{trace}\left(\b{\Gamma}_K\b{\Gamma}_K^{\top}\right)\right).
\label{eq:lyafcmims}
\end{equation}
Using~(\ref{eq:deroftrac}), we evaluate the time derivative of $V$ on the solutions of the state estimation error equation to get
\begin{eqnarray}
\dot{V} &=& \b e_{\b x}^{\top}\b P\dot{\b e}_{\b x}+\frac{1}{\alpha}\mb{trace}\left( \b{\Gamma}_K\dot{\b{\Gamma}}_K^{\top}\right)\nonumber\\
&=&\b e_{\b x}^{\top}\left( \b P\b A\b e_{\b x}+\b P\b{\Gamma}_K\b{\Phi}_K(\b x,\b u)+\b P\b e\right)+\frac{1}{\alpha}\mb{trace}\left( \b{\Gamma}_K\dot{\b{\Gamma}}_K^{\top}\right)\nonumber\\
&=& \frac{1}{2}\b e_{\b x}^{\top}\left(\b A^{\top}\b P+\b P\b A\right)\b e_{\b x}+\b e_{\b x}^{\top}\b P\b{\Gamma}_K\b{\Phi}_K(\b x,\b u)+\frac{1}{\alpha}\mb{trace}\left( \b{\Gamma}_K\dot{\b{\Gamma}}_K^{\top}\right)\nonumber\\
&{}&\quad+\b e_{\b x}^{\top}\b P\b e.
\label{eq:lyadermimser}
\end{eqnarray}
Using the definition of the trace operator, we obtain
\begin{eqnarray}
\b e_{\b x}^{\top}\b P\b{\Gamma}_K\b{\Phi}_K(\b x,\b u)&=&\mb{trace}\left( \b e_{\b x}^{\top}\b P\b{\Gamma}_K\b{\Phi}_K(\b x,\b u)\right)\nonumber\\
&=&\mb{trace}\left(\b{\Gamma}_K\b{\Phi}_K(\b x,\b u)\b e_{\b x}^{\top}\b P\right).
\label{eq:auxtoevesx}
\end{eqnarray}
Observe that
\[
\dot{\b W}_K=-\dot{\b{\Gamma}}_K.
\]
Using the above and~(\ref{eq:auxtoevesx}) along with the Lyapunov equation~(\ref{eq:lyaeqfoagiq}), we write
\begin{equation}
\dot{V}=-\b e_{\b x}^{\top}\b Q\b e_{\b x}+\frac{1}{\alpha}\left(\mb{trace}\left( \alpha\b{\Gamma}_K\b{\Phi}_K(\b x,\b u)\b e_{\b x}^{\top}\b P\right)-\mb{trace}\left( \b{\Gamma}_K\dot{\b{W}}_K^{\top}\right)\right)+\b e_{\b x}^{\top}\b P\b e.
\label{eq:inclyanpg}
\end{equation}
Consider the following weight matrix update law,
\begin{equation}
\fbox{$\displaystyle \dot{\b W}_K = \alpha\b P\b e_{\b x}\b{\Phi}_K^{\top}(\b x,\b u) $}
\label{eq:lerufomacw}
\end{equation}
With the above update law, the Lyapunov derivative~(\ref{eq:inclyanpg}) evaluates to
\begin{eqnarray}
\dot{V}&=&-\b e_{\b x}^{\top}\b Q\b e_{\b x}+\b e_{\b x}^{\top}\b P\b e\nonumber\\
&\le& -\lambda_{\min}(\b Q)\|\b e_{\b x}\|^2+\lambda_{\max}(\b P)\|\b e\|\|\b e_{\b x}\|\nonumber\\
&=&-\lambda_{\min}(\b Q)\left(\|\b e_{\b x}\|-\frac{\lambda_{\max}(\b P)}{\lambda_{\min}(\b Q)}\|\b e\| \right)\|\b e_{\b x}\|.
\label{eq:lyaderwitupl}
\end{eqnarray}
Utilizing the same type of arguments as in the scalar case (analyzed in Subsection~\ref{subsec:sistsiin}) and the Lyapunov-like lemma stated on page~\pageref{lemma:lyapunov_like}, we can show that when $\b e=\b 0$, then $\b e_x(t)\to \b 0$ as $t\to \infty$.

If the magnitude of the state estimation error is such that $\|\b e_{\b x}\|<\frac{\lambda_{\max}(\b P)}{\lambda_{\min}(\b Q)}\|\b e\|$, then it may happen that $\dot{V}>0$, which in turn may lead to the weights' drift to infinity over time. To prevent this, we modify~(\ref{eq:lerufomacw}) following the analysis of the single-input single-state case. The modified weight update law is
\begin{equation}
\fbox{$\displaystyle \dot{\b W}_K=\left\{\begin{array}{lcl}
\alpha\b P\b e_{\b x}\b{\Phi}_K^{\top}(\b x,\b u) &\mb{if}& \|\b e_{\b x}\| \ge e_0\quad\mb{and}\quad \|\b W_K\|\le M\\
\b O &{}& \mb{otherwise}
\end{array}\right.$}
\label{eq:mowulmims}
\end{equation}

\begin{example}\label{ex:mimo_neural_adaptive}\rm
Consider an uncertain dynamical system model,
\begin{eqnarray*}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2\\
\dot{x}_3
\end{array}\right]&=&\left[\begin{array}{ccc}
0 & 1 & 0\\
0 & 0 & 1\\
-116 & -69 & -14
\end{array}\right] \left[\begin{array}{c}
x_1\\
x_2\\
x_3
\end{array}\right]+\left[\begin{array}{cc}
0 & 0\\
1 & 0\\
0 & 1
\end{array}\right]\left(\left[\begin{array}{c}
u_1\\
u_2
\end{array}\right]+\left[\begin{array}{c}
\zeta _1\\
\zeta _2
\end{array}\right]\right)\\
&=& \b A\b x+\b B(\b u+\b{\zeta}),
\end{eqnarray*}
where the vector modeling the system uncertainties satisfies the norm bound condition,
\[
\left\|\left[\begin{array}{c}
\zeta _1\\
\zeta _2
\end{array}\right]\right\| \le 10.
\]
The uncontrolled nominal system $\dot{\b x}={\b A\b x}$ is asymptotically stable. The solution to the Lyapunov matrix equation ${\b A}^{\top}{\b P+\b P\b A}=-{\b I}_3$ is
\[
{\b P}=\left[\begin{array}{ccc}
9.2833 & 4.8451 & 0.0043\\
4.8451 & 3.9263 & 0.0775\\
0.0043 & 0.0775 & 0.0412
\end{array}\right].
\]
Our objective is to design a stabilizing controller, that contains an {\bf  adaptive neural term} employing 10 radial basis functions, so that the state of the closed-loop system driven by your controller asymptotically converges to the origin.

We propose the following stabilizing controller,
\begin{equation}
\b u=-\b u_{neural}-\eta\frac{\b B^{\top}\b P\b x}{\|\b B^{\top}\b P\b x\|},
\label{neur_contr}
\end{equation}
where
\[
\b u_{neural}=\b W\b{\Phi}
\]
and the gain $\eta$ will be determined later. The vector $\b{\Phi}$ is composed of radial basis functions and is bounded. We assume that the ``optimal'' weights approximate the uncertainties modeled by the function $\b{\zeta}$ so that
\begin{equation}
\b{\zeta}=\b W^{*}\b{\Phi}+\b e,
\label{eq:zetaappropti}
\end{equation}
where $\b e$ is the approximation error. We assume that the approximation error is bounded by a known constant $\rho$, that is,
\[
\|\b e\|\le \rho.
\]
In further discussion, we use the following notation,
\[
\b{\sigma}=\b B^{\top}\b P\b x\;\;\mb{and let}\;\; \b{\Gamma}=\b W^*-\b W.
\]
Employing the above notation, we represent the proposed control strategy~(\ref{neur_contr}) as
\begin{equation}
\b u=-\b W\b{\Phi}-\eta\frac{\b{\sigma}}{\|\b{\sigma}\|}.
\label{eq:cont_sig}
\end{equation}
Consider the following Lyapunov function candidate,
\[
V = \frac{1}{2}\b x^{\top}\b P\b x + \frac{1}{2\alpha} \mb{trace} \left( \b{\Gamma}\b{\Gamma}^{\top}  \right).
\]
The Lyapunov derivative is
\begin{eqnarray*}
\dot{V} &=& \b x^{\top}\b P\dot{\b x} + \frac{1}{\alpha}\mb{trace} \left( \b{\Gamma}\dot{\b{\Gamma}}^{\top}  \right)\\
&=& \b x^{\top}\b P\left(  \b A\b x+\b B(\b u+\b{\zeta})\right)+ \frac{1}{\alpha}\mb{trace} \left( \b{\Gamma}\dot{\b{\Gamma}}^{\top}  \right)\\
&=& -\frac{1}{2}\|\b x\|^2+\b x^{\top}\b P\b B\b u + \b x^{\top}\b P\b B\b{\zeta} + \frac{1}{\alpha}\mb{trace} \left( \b{\Gamma}\dot{\b{\Gamma}}^{\top}  \right).
\end{eqnarray*}
Substituting into the above~(\ref{eq:zetaappropti}) and the control strategy given by~(\ref{eq:cont_sig}) and manipulating gives
\begin{eqnarray*}
\dot{V} &=& -\frac{1}{2}\|\b x\|^2 -\eta \|\b{\sigma}\| -\b{\sigma}^{\top}\b W\b{\Phi}+ \b{\sigma}^{\top}\left(\b W^{*}\b{\Phi}+\b e\right) + \frac{1}{\alpha}\mb{trace} \left( \b{\Gamma}\dot{\b{\Gamma}}^{\top}  \right)\\
&=& -\frac{1}{2}\|\b x\|^2 -\eta \|\b{\sigma}\|+\b{\sigma}^{\top}\left(\b W^*-\b W   \right)\b{\Phi} + \b{\sigma}^{\top}\b e + \frac{1}{\alpha}\mb{trace} \left( \b{\Gamma}\dot{\b{\Gamma}}^{\top}  \right)\\
&\le& -\frac{1}{2}\|\b x\|^2 -\|\b{\sigma}\|\left(\eta  - \|\b e\|\right) +\b{\sigma}^{\top}\b{\Gamma}\b{\Phi}+ \frac{1}{\alpha}\mb{trace} \left( \b{\Gamma}\dot{\b{\Gamma}}^{\top}  \right)\\
&=& -\frac{1}{2}\|\b x\|^2 -\|\b{\sigma}\|\left(\eta  - \|\b e\|\right) +\mb{trace}\left(\b{\sigma}^{\top}\b{\Gamma}\b{\Phi}\right)+ \frac{1}{\alpha}\mb{trace} \left( \b{\Gamma}\dot{\b{\Gamma}}^{\top}  \right)\\
&=& -\frac{1}{2}\|\b x\|^2 -\|\b{\sigma}\|\left(\eta  - \|\b e\|\right) +\mb{trace}\left(\b{\Gamma}\b{\Phi}\b{\sigma}^{\top}+ \frac{1}{\alpha} \b{\Gamma}\dot{\b{\Gamma}}^{\top}  \right).
\end{eqnarray*}
We select
\[
\eta  \ge \|\b e\|\;\;\mb{and}\;\; \dot{\b{\Gamma}}=-\alpha \b{\sigma}\b{\Phi}^{\top}
\]
to obtain
\begin{equation}
\dot{V}\le -\frac{1}{2}\|\b x\|^2,
\label{eq:neg_def_ld}
\end{equation}
which is negative semi-definite in the augmented $(\b x,\b{\Gamma})$-space. Therefore $\b x$ as well as $\b{W}$ are bounded. Hence the proposed controller is also bounded, which implies that $\dot{\b x}$ is bounded. This in turn implies that
\[
\frac{d}{dt}\|\b x\|^2=2\b x^{\top}\dot{\b x}
\]
is bounded and hence $\|\b x\|^2$ is uniformly continuous. Integrating both sides of~(\ref{eq:neg_def_ld}) gives
\[
\int_0^{t}\dot{V}(\tau)d\tau = V(t)-V(0)\le \int_0^{t}\left( -\frac{1}{2}\|\b x\|^2 \right)d\tau
\]
and hence
\[
\int_0^{t} \|\b x\|^2 d\tau\le 2 (V(0)-V(t))\le 2V(0),
\]
which means that
\[
\lim_{t\to \infty}\int_0^{t} \|\b x\|^2 d\tau
\]
exists and is finite. By Barbalat's lemma,
\[
\|\b x\|^2 \to 0 \;\;\mb{as}\;\; t\to \infty,
\]
which implies that
\[
\lim_{t\to \infty} \b x = \b 0.
\]
\end{example}
\markex

\section{A Self-Organizing Network}
Up to now we discussed neural networks that must have a teacher in order to learn presented patterns. This mode of learning is referred to as the supervised mode. In this section, we analyze a neural network that can learn in the unsupervised mode and the network undergoes the self-organization process. Specifically, the network is capable of clustering training patterns. It is assumed that the number of classes, $p$, is known. The network's number of weight vectors is equal to the number of classes in the training set. The weights are initialized and then normalized. After initialization, the training patterns, $\b x^{(1)}, \b x^{(2)}, \ldots, \b x^{(N)}$, are presented. In response to each presented pattern, we compute the winning neuron of the network, that is, we compute the network weight that is closest to the presented pattern. We call this computation the {\em similarity matching\/}. Let $\b x$ be the presented pattern, where, for simplicity, we omitted its index. Then the similarity matching can be described as
\[
\left\| \b x-\b w^{(i)}\right\|=\min_{1\le j \le p}\left\{ \left\| \b x-\b w^{(j)}\right\|\right\},
\]
where $i$ is the index of the winning neuron. Thus, searching for the winning neuron corresponds to finding the maximum among the $p$ scalar products,
\[
\max_{j=1,2,\ldots,p}\left\langle \b w^{(j)}, \b x\right\rangle.
\]
To see that this is indeed the case, note that because the weights and the patterns are normalized, that is, $\left\|\b w^{(j)} \right\|=\left\|\b x \right\|=1$, we have
\begin{eqnarray*}
\left\| \b x-\b w^{(j)}\right\|&=&\left( \left\|\b x \right\|^2-2\left\langle \b w^{(j)},\b x\right\rangle + \left\|\b w^{(j)} \right\|^2\right)^{1/2}\\
&=&\left(2- 2\left\langle \b w^{(j)},\b x\right\rangle\right)^{1/2}.
\end{eqnarray*}
Hence, minimizing $\left\| \b x-\b w^{(j)}\right\|$ is equivalent to maximizing the scalar product $\langle \b w^{(j)},\b x\rangle$. On the other hand,
\[
\langle \b w^{(j)},\b x\rangle=\left\| \b w^{(j)}\right\|\left\| \b x\right\|\cos\left(\b w^{(j)},\b x\right),
\]
which means that the winning neuron is ``more parallel'' to the input pattern $\b x$ than the other neurons.

After the winning neuron is declared, it is updated, while the remaining weights are not changed. This is the reason for the name of the network---the winner-take-all~(WTA) network. This update algorithm is also known as the {\em Kohonen learning rule\/} after Kohonen who first proposed this learning rule. The role of the update is to reduce the distance between the winning neuron, $\b w^{(i)}$, and the presented pattern, $\b x$. One way to accomplish this task is to apply the descent gradient method to the objective function $\frac{1}{2}\left\| \b x-\b w^{(i)}\right\|^2$. The gradient of $\frac{1}{2}\left\| \b x-\b w^{(i)}\right\|^2$ with respect to $\b w^{(i)}$ is
\[
\frac{1}{2}\nabla_{\b w^{(j)}}\left\| \b x-\b w^{(i)}\right\|^2=-\left( \b x-\b w^{(i)}\right).
\]
Applying the descent gradient method to minimize $\frac{1}{2}\left\| \b x-\b w^{(i)}\right\|^2$ with a step size $\alpha_k$ yields
\[
\b w^{(i)}(k+1)=\b w^{(i)}(k)+\alpha_k\left( \b x-\b w^{(i)}(k)\right).
\]
The weights of other neurons are unchanged, that is,
\[
\b w^{(j)}(k+1)=\b w^{(j)}(k),\quad\mb{for}\quad j\ne i.
\]
We summarize our discussion in the form of the following algorithm.
\begin{center}
{\bf Winner-Take-All~(WTA) Learning Algorithm}
\end{center}
\begin{enumerate}
\item {\bf Initialize Weights}\\
Generate the initial weights. Normalize each weight to have unit length: $\frac{\b w}{\|\b w\|}$.
\item {\bf Serially Present Training Patterns}\\
Normalize each training pattern $\b x^{(j)}$ to have unit length. The training patterns can be presented in a random order a fixed number of times.
\item {\bf Select Winning Neuron}\\
For each WTA weight, compute the scalar product with the presented pattern. Select a weight that gives the largest scalar product. Call this the winning weight.
\item {\bf Update Weight of Winning Neuron}\\
Let $i$ be the index of the winning weight. Set
\begin{eqnarray*}
\b w^{(i)}(k+1)&=&\b w^{(i)}(k)+\alpha_k\left(\b x^{(j)}-\b w^{(i)}(k)\right)\\
\b w^{(j)}(k+1)&=&\b w^{(j)}(k),\quad\mb{for}\quad j\ne i
\end{eqnarray*}
Normalize $\b w^{(i)}(k+1)$ to have unit length.
\item {\bf Repeat by Returning to Step 2}\\
Repeat until stopping criterion is met.
\end{enumerate}
\markex
\begin{example}\label{ex:wta_zur}\rm
We take the same set of training patterns as in Zurada~\cite[page~407]{zurada92}:
\[
\left[\begin{array}{ccccc}
\b x^{(1)}&\b x^{(2)}&\b x^{(3)}&\b x^{(4)}&\b x^{(5)}
\end{array}\right]=\left[\begin{array}{ccccc}
0.8  &  0.1736 & 0.707 & 0.342 & 0.6\\
0.6  & -0.9848 & 0.707 & -0.9397 & 0.8
\end{array}\right].
\]
The above patterns are already normalized and are marked in Figure~\ref{fig:wta_ex}. As can be seen from Figure~\ref{fig:wta_ex}, we can classify the training patterns into two classes, that is, $p=2$. The initial normalized weights are
\[
\left[\begin{array}{cc}
\b w^{(1)}_0 & \b w^{(2)}_0
\end{array}\right]=\left[\begin{array}{cc}
1 & -1\\
0 & 0
\end{array}\right].
\]
The initial weights are also marked in Figure~\ref{fig:wta_ex}. The training patterns, $\b x^{(1)}, \b x^{(2)}, \ldots, \b x^{(5)}$ are presented in ascending order and then recycled. After 200 presentations of the training set, we obtained the following final weights of the network:
\[
\left[\begin{array}{cc}
\b w^{(1)} & \b w^{(2)}
\end{array}\right]
=\left[\begin{array}{cc}
0.6628  &  0.2868\\
0.7488  & -0.9580
\end{array}\right].
\]
The final weights are marked in Figure~\ref{fig:wta_ex} with $\ast$.
\begin{figure}
%\vspace{3in}
\centering \epsfig{file=wta_ex.eps,height=3in}
%\hskip1.75in \special{eps:c:/bsb/wta_ex.eps y=3in}
\caption{Training patterns, initial weights and the final weights of the WTA network of Example~\ref{ex:wta_zur}.}
\label{fig:wta_ex}
\end{figure}
\end{example}\markex

\section{Hopfield Neural Network}
\subsection{Hopfield Neural Network Modeling and Analysis}
The analog Hopfield neural network is a dynamical system whose schematic is shown in Figure~\ref{fig:hopfi_1}. It is a nonlinear interconnected system proposed by Hopfield~\cite{hopfield84} as a model of biological neural networks.
\begin{figure}
%\vspace{4.88in}
\centering \epsfig{file=hopfi_1.eps,height=5in}
%\centerbmp{6.5in}{4.88in}{c:/mybook/hopfi_1.bmp}
%\special{bmp:c:/mybook/hopfi_1.bmp x=6.5in y=4.88in}
\caption{Hopfield neural network model.}
\label{fig:hopfi_1}
\end{figure}
Each node, or neuron, in Figure~\ref{fig:hopfi_1} is represented by a circle. The structure of a neuron is depicted in Figure~\ref{fig:hopfi_2}.
\begin{figure}
%\vspace{3.5in}
\centering \epsfig{file=hopfi_2.eps,height=2in}
%\centerbmp{4.67in}{3.5in}{c:/mybook/hopfi_2.bmp}
%\hskip1in\special{bmp:c:/mybook/hopfi_2.bmp x=4.67in y=3.5in}
\caption{Model of a neuron in the Hopfield neural network.}
\label{fig:hopfi_2}
\end{figure}
The nonlinear amplifier in Figure~\ref{fig:hopfi_2} has an input-output characteristic described by a function $g:\mb{\Bb R} \to [0, 1]$. We assume that $g$ is continuously differentiable and strictly monotonically increasing, that is, $g(u)>g(u')$ if $u>u'$. The function $g$ is sometimes called, the {\em activation function\/}. A typical choice of $g$ is the sigmoid function given by
\begin{equation}
x=g(u)=\frac{1}{1+e^{-u/\theta}},
\label{eq:sigmofoho}
\end{equation}
where the positive parameter $\theta$ controls the slope of $g$. A circuit implementation of the Hopfield network is shown in Figure~\ref{fig:hopfi_3}.
\begin{figure}
%\vspace{4.88in}
\centering \epsfig{file=hopfi_3.eps,height=4.5in}
%\special{bmp:c:/mybook/hopfi_3.bmp x=6.5in y=4.88in}
\caption{A circuit realization of the Hopfield neural network.}
\label{fig:hopfi_3}
\end{figure}
To derive equations that model the Hopfield neural network, we use an equivalent representation of the single node of the Hopfield network. This equivalent representation is shown in Figure~\ref{fig:hopfi_4}.
\begin{figure}
%\vspace{3.5in}
\centering \epsfig{file=hopfi_4.eps,height=2in}
%\hskip.75in\special{bmp:c:/mybook/hopfi_4.bmp x=4.67in y=3.5in}
\caption{Equivalent representation of a single node in the circuit realization of the Hopfield neural network.}
\label{fig:hopfi_4}
\end{figure}
Applying Kirchhoff's current law at the input of each amplifier of the equivalent node representation and taking into account that the current into the amplifier is negligibly small due to the amplifier's high input resistance, we get
\begin{equation}
C_i\frac{du_i}{dt}+\frac{u_i}{R_i}=\frac{x_1-u_i}{R_{i1}}+\frac{x_2-u_i}{R_{i2}}+\cdots+\frac{x_n-u_i}{R_{in}}+I_i,\quad i=1,2,\ldots,n,
\label{eq:hoeqkcl}
\end{equation}
where $x_i=g(u_i)$. Let
\begin{equation}
T_{ij}=\frac{1}{R_{ij}},\quad \frac{1}{r_i}=\sum_{j=1}^{n}T_{ij}+\frac{1}{R_i}.
\label{eq:imtoresnt}
\end{equation}
Substituting~(\ref{eq:imtoresnt}) into~(\ref{eq:hoeqkcl}) yields
\begin{equation}
\fbox{$ \displaystyle C_i\frac{d u_i}{dt}+\frac{u_i}{r_i}=\sum_{j=1}^n T_{ij}x_j +I_i,\quad i=1,2,\ldots,n $}
\label{eq:stahoal}
\end{equation}
The system $n$ first-order nonlinear differential equations given by~(\ref{eq:stahoal}) constitutes a model of the Hopfield neural network. The Hopfield neural network can be viewed as an associative memory. This is because it can have asymptotically stable states that attract the neighboring states to develop in time. This behavior can be interpreted as an evolution of an imperfect, corrupted pattern towards the correct, stored pattern---hence the term associative memory. More on neural associative memory can be found in Subsection~\ref{subsec:assoscmem}.

In further analysis, we use an equivalent representation of~(\ref{eq:stahoal}). To obtain this equivalent representation of the Hopfield model, we first note that the activation function $g$ is, by assumption, strictly increasing and differentiable. Therefore $g^{-1}$ exists and is also differentiable---see the inverse function theorem on page~\pageref{invfuntgen}. In Figure~\ref{fig:sigainv}, we show plots of $x=g(u)$ and $u=g^{-1}(x)$ for $\theta=0.2$.
\begin{figure}
\centering
\includegraphics[scale=0.5]{sigainv.eps}
%\vspace{3.5in}
%\hskip1.2in\special{bmp:c:/mybook/sigainv.bmp x=3.97in y=3.5in}
\caption{Plots of $x=g(u)$ and $u=g^{-1}(x)$.}
\label{fig:sigainv}
\end{figure}
We re-write equations~(\ref{eq:stahoal}) as
\begin{equation}
C_i\frac{d g^{-1}(x_i)}{dx_i}\frac{d x_i}{dt}+\frac{g^{-1}(x_i)}{r_i}=\sum_{j=1}^n T_{ij}x_j +I_i,\quad i=1,2,\ldots,n
\label{eq:eq:stahoalalt}
\end{equation}
Let $\b T \in \mb{\Bb R}^{n\times n}$ be the matrix whose elements are $T_{ij}$ and $\b v=\left[\begin{array}{cccc}
I_1&I_2&\cdots I_n
\end{array}\right]^{\top}\in \mb{\Bb R}^n$. Using this notation, we write~(\ref{eq:stahoal}) in matrix form,
\begin{equation}
\left[\begin{array}{ccc}
C_1\frac{dg^{-1}(x_1)}{d x_1}& \cdots & 0\\
\vdots & \ddots & \vdots\\
0 & \cdots & C_n\frac{dg^{-1}(x_n)}{d x_n}
\end{array}\right]\dot{\b x}(t)=\b T\b x+\b v-\left[\begin{array}{c}
\frac{g^{-1}(x_1)}{r_1}\\
\vdots\\
\frac{g^{-1}(x_n)}{r_n}
\end{array}\right].
\label{eq:matreofhom}
\end{equation}
Associated with the Hopfield neural network is its computational energy function, $E:\mb{\Bb R}^n \to \mb{\Bb R}$, given by
\begin{eqnarray}
E(\b x)&=&-\frac{1}{2}\sum_{i=1}^n\sum_{j=1}^n T_{ij}x_i x_j -\sum_{i=1}^n I_i x_i +\sum_{i=1}^n \int_{0}^{x_i}\frac{g^{-1}(s)}{r_i}ds\nonumber\\
&=&-\frac{1}{2}\b x^{\top}\b T\b x-\b v^{\top}\b x +\sum_{i=1}^n \int_0^{x_i}\frac{g^{-1}(s)}{r_i}ds.
\label{eq:enefuhop}
\end{eqnarray}
Suppose now that $\b T=\b T^{\top}$. Then, the gradient of $E$ is
\begin{equation}
\nabla E(\b x)=-\b T\b x-\b v+\left[\begin{array}{ccc}
\frac{g^{-1}(x_1)}{r_1}&\cdots&
\frac{g^{-1}(x_n)}{r_n}
\end{array}\right]^{\top}.
\label{eq:grofcoenfu}
\end{equation}
Let
\begin{equation}
\b S(\b x)=\left[\begin{array}{ccc}
C_1\frac{dg^{-1}(x_1)}{d x_1}& \cdots & 0\\
\vdots & \ddots & \vdots\\
0 & \cdots & C_n\frac{dg^{-1}(x_n)}{d x_n}
\end{array}\right].
\label{eq:defofbs}
\end{equation}
We now combine~(\ref{eq:grofcoenfu}) and~(\ref{eq:defofbs}) to obtain an equivalent form of the Hopfield network modeling equations,
\begin{equation}
\fbox{$\displaystyle\b S(\b x)\dot{\b x}(t)=-\nabla E(\b x(t)) $}
\label{eq:alagrsys}
\end{equation}

We analyze the dynamics of the Hopfield network modeled by~(\ref{eq:alagrsys}) using the computational energy function. We first show that $E$ is bounded. Taking norms of~(\ref{eq:enefuhop}) yields
\begin{equation}
\left\arrowvert E(\b x)\right\arrowvert \le \frac{1}{2}\left\|\b T\right\|\|\b x\|^2 +\|\b v\|\|\b x\|+\sum_{i=1}^n \frac{1}{r_i}\left\arrowvert \int_0^{x_i}g^{-1}(s) ds\right\arrowvert.
\label{eq:abvaofe}
\end{equation}
The components of the vector $\b x$ are the outputs of the nonlinear amplifiers and they take values from the interval $[0, 1]$. Hence, $\|\b x\|=\|\b x\|_2 \le \sqrt{n}$. Thus, to show that $E(\b x)$ is bounded it is enough to show that the integral $\int_0^{x_i}g^{-1}(s) ds$ is bounded. We first find the expression for $g^{-1}$. Recall that $g(u)$ is given by~(\ref{eq:sigmofoho}). Therefore,
\[
\frac{1}{x}=1+e^{-u/\theta},
\]
or
\begin{equation}
\frac{1}{x}-1=e^{-u/\theta},
\label{eq:copinog}
\end{equation}
Taking logarithms of both sides of the above yields
\[
\ln\left(\frac{1-x}{x}\right)=-\frac{u}{\theta}.
\]
Hence,
\begin{equation}
u=g^{-1}(x)=\theta\left( \ln(x)-\ln(1-x)\right).
\label{eq:invgacf}
\end{equation}
Using the symmetry property of $g^{-1}$, we obtain
\begin{eqnarray}
\left\arrowvert \int_0^{x_i}g^{-1}(s) ds\right\arrowvert &\le & \int_0^{x_i}\left|g^{-1}(s)\right| ds\nonumber\\
&\le & -2\int_0^{1/2}g^{-1}(s)ds.
\label{eq:boinogin}
\end{eqnarray}
Recall that
\[
\int \ln (s)ds =s \ln(s)-s.
\]
Hence,
\[
\int\ln(1-s)ds=-(1-s)\ln(1-s)-s.
\]
Using the above, we obtain
\begin{eqnarray}
\int_0^{1/2}g^{-1}(s)ds&=&\int_0^{1/2}\theta\left( \ln(s)-\ln(1-s)\right)ds\nonumber\\
&=& \theta\left.\left( s\ln(s)-s +(1-s)\ln(1-s) +s \right)\right\arrowvert^{1/2}_0 \nonumber\\
&=& \theta\left.\left( s\ln(s)+(1-s)\ln(1-s) \right)\right\arrowvert^{1/2}_0.
\label{eq:inginoh}
\end{eqnarray}
Note that
\[
s\ln(s)=\frac{\ln(s)}{1/s}.
\]
Applying L'Hospital's rule to the above yields
\[
\lim_{s\to 0}\frac{\ln(s)}{1/s}=\lim_{s\to 0}\frac{1/s}{-1/s^2}=0.
\]
Applying the above to~(\ref{eq:inginoh}), we get
\begin{equation}
\int_0^{1/2}g^{-1}(s)ds=\theta\ln(1/2).
\label{eq:bouonigi}
\end{equation}
Combining~(\ref{eq:boinogin}) and~(\ref{eq:bouonigi}) yields
\begin{equation}
\left\arrowvert \int_0^{x_i}g^{-1}(s) ds\right\arrowvert \le -2\theta\ln(1/2),
\label{eq:bouoningi}
\end{equation}
which completes the proof of the statement that $E$ is bounded.

Let $\b x(t)$ denote the trajectory of the Hopfield network model~(\ref{eq:alagrsys}). We will show that if $\b T=\b T^{\top}$, then
\[
\frac{d}{dt}E(\b x(t))\le 0
\]
on the Hopfield network trajectories. Using the chain rule gives
\begin{equation}
\frac{d}{dt}E(\b x(t))=\left( \nabla E(\b x(t))\right)^{\top}\dot{\b x}(t).
\label{eq:uschrutef}
\end{equation}
It follows from~(\ref{eq:alagrsys}) that
\begin{equation}
\nabla E(\b x(t))=-\b S(\b x(t))\dot{\b x}(t).
\label{eq:usprmoxd}
\end{equation}
Substituting~(\ref{eq:usprmoxd}) into~(\ref{eq:uschrutef}), we obtain
\begin{eqnarray}
\frac{d}{dt}E(\b x(t))&=&-\dot{\b x}^{\top}(t)\b S(\b x(t))\dot{\b x}(t)\nonumber\\
&=&-\sum_{i=1}^n C_i\frac{dg^{-1}(x_i(t))}{d x_i}\dot{x}_i^2,
\label{eq:deocoefu}
\end{eqnarray}
where $C_i>0$ and
\begin{equation}
\frac{d g^{-1}(x_i)}{d x_i}=\theta\left(1/x_i+1/(1-x_i) \right)>0
\label{eq:posderofginv}
\end{equation}
because $x_i\in(0,1)$. Hence,
\[
\frac{d}{dt}E(\b x(t))\le 0
\]
on the trajectories of the Hopfield model~(\ref{eq:alagrsys}).

\subsection{Analog-to-Digital Converter}
We will now show how the Hopfield neural network can be used to solve some engineering problems. Specifically, we will discuss the construction of an analog-to-digital~(A/D) converter using the Hopfield neural network. The idea was put forward by Tank and Hopfield~\cite{tank86}. In A/D conversion, an analog value is transformed into its corresponding binary representation. The conversion that we will discuss consists of two steps. The first step is to round the given analog value to the nearest integer. The second step is to convert the integer into its binary equivalent. For example, 3.4 is first rounded to 3, then 3 is converted into its binary equivalent ``11'' if we have a 2-bit converter, ``011'' if we have a 3-bit converter, or ``0011'' if we have a 4-bit converter. Therefore, an $N$-bit A/D converter can be viewed as a device that assigns to a given analog input a corresponding vertex of the hypercube $[0, 1]^N$.

The first step in A/D conversion, rounding of a given number to its nearest integer, can be organized as solving an optimization problem by a neural network. This entails formulation of a suitable cost, or penalty, function to be minimized. In fact, finding a suitable cost function will be the first step in constructing a neural A/D converter.

To better understand the idea behind solving A/D conversion problems using neural networks, we first consider constructing a 2-bit A/D converter. Let $V_1$ be the most significant bit~(MSB), $V_0$ be the least significant bit~(LSB) and $x$ an analog value to be converted. In the actual neural network, $V_1$ and $V_0$ are the outputs of the two neurons in the neural A/D converter. We first consider a very simple cost function,
\begin{equation}
E=\frac{1}{2}\left( V_1 2^1+V_0 2^0-x\right)^2.
\label{eq:quacost}
\end{equation}
A plot of the above function, for an analog input of $x=1.5$, is shown in Figure~\ref{fig:ad_conv}(a). Note that $E$ is defined over the two-dimensional hypercube $[0,1]^2$.
\begin{figure}
\centering
\includegraphics[scale=0.5]{ad_conv.eps}
%\vspace{4in}
%\hskip.75in \special{eps:c:/brain/ad_conv.eps y=4in}
\caption{Constructing the cost function for a 2-bit A/D converter---plots of the cost function and its components; (a) plot of~(\ref{eq:quacost}) for an analog input of $x=1.5$; (b) plot of the term used to cancel the quadratic terms in~(\ref{eq:quacost}); (c) plot of the cost function~(\ref{eq:final_cost}).}
\label{fig:ad_conv}
\end{figure}
However, constructing a circuit minimizing the cost function~(\ref{eq:quacost}) would not yield good results because $V_0$ and $V_1$ would take values from the interval $[0,1]$ rather than the set $\{0, 1\}$. To overcome the problem, we add to~(\ref{eq:quacost}) two terms whose purpose is to penalize for non-integer values of $V_0$ and $V_1$. The modified cost function is
\begin{equation}
E=\frac{1}{2}\left( V_1 2^1+V_0 2^0-x\right)^2+aV_1(V_1-1)+bV_0(V_0-1),
\label{eq:mod_cost}
\end{equation}
where $a$ and $b$ are weighting coefficients. Tank and Hopfield~\cite{tank86} suggest selecting $a$ and $b$ so that the quadratic terms $V_0^2$ and $V_1^2$ are cancelled out. In particular, if we set $a=-2$ and $b=-0.5$, then~(\ref{eq:mod_cost}) becomes
\begin{equation}
E=\frac{1}{2}x^2+2V_0V_1-2V_1x-V_0x+2V_1+\frac{1}{2}V_0.
\label{eq:final_cost}
\end{equation}
The term used to cancel the quadratic terms is plotted in Figure~\ref{fig:ad_conv}(b). The final cost function given by~(\ref{eq:final_cost}) is plotted in Figure~\ref{fig:ad_conv}(c). When an analog input $x$ is converted to its binary equivalent, the cost function~(\ref{eq:final_cost}) is minimized. One can interpret this process as follows. A given analog input $x$ defines a corresponding cost function $E$. We can view $E$ as a landscape. Now, a ball placed somewhere on this landscape will move towards the lowest point, which should correspond to one of the vertices of the hypercube $[0,1]^2$.

Our next task is to construct differential equations that model a neural network implementing an A/D converter. The parameters of the state equations of the network will be chosen so that the time derivative of $E$ evaluated on the solutions of the equations describing the network is negative. We consider a general case of an $N$-bit A/D converter. After this, we derive state equations for a 2-bit A/D converter. Note that
\[
\frac{dE}{dt}=\frac{\partial E}{\partial V_0}\dot{V}_0+\frac{\partial E}{\partial V_1}\dot{V}_1+\cdots + \frac{\partial E}{\partial V_{N-1}}\dot{V}_{N-1},
\]
where $V_i=g_i(u_i)$. Hence, similarly as in~(\ref{eq:deocoefu}), it is enough to set
\begin{equation}
\frac{\partial E}{\partial V_i}=-C_i \dot{u}_i=-C_i \frac{d g_i^{-1}(V_i)}{d V_i}\dot{V}_i,\quad i=0,1,\ldots,N-1
\label{eq:parewrtvi}
\end{equation}
to obtain
\begin{equation}
\frac{d E}{dt}=-\sum_{i=0}^{N-1} C_i \frac{d g_i^{-1}(V_i)}{d V_i}\dot{V}_i^2\le 0,
\label{eq:dedtinad}
\end{equation}
where as in~(\ref{eq:posderofginv}), $\frac{d g_i^{-1}(V_i)}{d V_i}>0$. Thus, if the state equations of an $N$-bit A/D converter satisfy~(\ref{eq:parewrtvi}), then $dE/dt\le 0$ on the trajectories of the converter. Using~(\ref{eq:parewrtvi}), we can write the state equations that model the 2-bit A/D converter consisting of two neurons. Let $V_0=g_0(u_0)$ and $V_1=g_1(u_1)$. Then, substituting~(\ref{eq:final_cost}) into~(\ref{eq:parewrtvi}) and performing required differentiation, we obtain
\begin{equation}
\left.\begin{array}{ccl}
C_0\dot{u}_0 &=& -\frac{\partial E}{\partial V_0}=-\left( 2V_1-x+0.5\right)=-(2g_1(u_1)-x+0.5)\\
C_1\dot{u}_1 &=& -\frac{\partial E}{\partial V_1}=-2\left( V_0-x+1\right)=-2(g_0(u_0)-x+1)
\end{array}\right\}
\label{eq:twobitadcoe}
\end{equation}
For simplicity, we assume that $C_0=C_1=1$ and that $g_0=g_1=g$. Note that because $C$ is present in both state equations, assigning $C_0=C_1=1$ does not distort the dynamic behavior of the network; it only changes the time scale of both equations by the same factor. Observing that for $i=1,2$,
\[
\dot{u}_i=\frac{d g_i^{-1}(V_i)}{d V_i}\dot{V}_i,
\]
we represent~(\ref{eq:twobitadcoe}) as
\begin{equation}
\left.\begin{array}{ccl}
\dot{V}_0 &=& -\left( 2V_1-x+0.5\right)\left( \frac{d g^{-1}(V_0)}{d V_0}\right)^{-1}\\
\dot{V}_1 &=& -2\left( V_0-x+1\right)\left(\frac{d g^{-1}(V_1)}{d V_1}\right)^{-1}
\end{array}\right\}
\label{eq:stsptwobitad}
\end{equation}
We display a plot of the cost function $E$ for the analog value $x=1.5$ in Figure~(\ref{fig:adphase})(a) and a phase portrait of the 2-bit A/D converter in Figure~(\ref{fig:adphase})(b). We chose the activation function
\[
V_i=g(u_i)=\frac{1}{1+e^{-4u_i}}.
\]
Using~(\ref{eq:invgacf}), we find
\[
u_i=g^{-1}(V_i)=\frac{1}{4}\left( \ln(V_i)-\ln(1-V_i)\right).
\]
Thus,
\[
\frac{d g^{-1}(V_i)}{d V_i}=\frac{1}{4 V_i(1-V_i)}.
\]
\begin{figure}
\centering
\includegraphics[scale=0.5]{adphase.eps}
%\vspace{4in}
%\hskip.75in\special{eps:c:/brain/adphase.eps y=4in}
\caption{A plot of the cost function and a phase portrait of the 2-bit A/D converter for the analog value $x=1.5$.}
\label{fig:adphase}
\end{figure}
Note that for $x=1.5$, the cost function has two minimizers in the hypercube $[0,1]^2$ that are separated by a ``ridge''. Initial conditions on each side of this ridge yield a different final value of the state vector $[V_0, V_1]$. Depending on the initial conditions of the state vector of the A/D converter, the analog value $x=1.5$ will either be converted to the binary ``01'' or ``10''. Thus, initial conditions influence the dynamic behavior of the A/D converter. Also, the  activation nonlinearity influences the dynamical behavior of the A/D converter via the multiplicative term $\left( \frac{d g^{-1}(V_i)}{d V_i}\right)^{-1}$ present in the state equations~(\ref{eq:stsptwobitad}). One can see that the multiplicative factor influences the direction in which the converter state trajectory evolves.

Employing the above technique, we can construct higher order A/D converters. For example, to construct a 3-bit A/D converter, we begin with the cost function
\[
E=\frac{1}{2}\left( V_2 2^2+V_1 2^1 +V_0 2^0-x\right)^2+aV_2(V_2-1)+bV_1(V_1-1)+cV_0(V_0-1).
\]
To cancel the quadratic terms, select $a=-8$, $b=-2$, $c=-1/2$. The state equations of the converter are obtained from the equations
\[
\dot{V}_i=-\frac{1}{C_i}\frac{\partial E}{\partial V_i}\left( \frac{d g_i^{-1}(V_i)}{d V_i}\right)^{-1}\quad i=0,1,2.
\]

In our analysis of A/D neural converters, we observed that while the 2-bit converter yielded correct results, the 3-bit and 4-bit converters showed some deficiency in their performance. The errors observed for these two converters were quite significant and increased with the number of bits and consequently with network size.

%\section{Hopfield Network Stability Analysis}
%
%In this section, we apply Lyapunov's second method to analyze stability in the sense of Lyapunov of the Hopfield neural network. The neural network is viewed as a large interconnected dynamical system, that is, as an interconnection of many neurons. We first analyze stability properties of a single neuron, and then proceed with the stability analysis of the whole network in terms of the properties of the single neurons. Interconnection information is used in conjunction with the properties of the individual neurons to obtain sufficient conditions for asymptotic stability of the whole network.
%
%
%\subsection{Hopfield Network Model Analysis}
%Recall the equations describing the Hopfield neural net. They are given by~(\ref{eq:stahoal}), that is,
%\[
%C_i\frac{du_i}{dt}+\frac{u_i}{r_i}=\sum_{j=1}^{n}T_{ij}x_j + I_i,\quad i=1,2,\ldots,n,
%\]
%where $x_i=g_i(u_i)$. In our further analysis, we assume that the activation functions $g_i:\mb{\Bb R} \to (-1,1)$, $i=1,2,\ldots,n$, are continuously differentiable and strictly monotonically increasing. We also assume that
%\[
%u_ig_i(u_i)>0\quad \mbox{for} \quad u_i \ne 0,
%\]
%and
%\[
%g_i(0)=0.
%\]
%To proceed further, we divide both sides of the above equations describing the Hopfield network by $C_i>0$ to get
%\begin{equation}
%\dot{u}_i=-b_iu_i+\sum_{j=1}^{n}A_{ij}g_j(u_j) + U_i,\quad i=1,2,\ldots,n,
%\label{eq:newfoofhof}
%\end{equation}
%where $b_i=\frac{1}{r_iC_i}$, $A_{ij}=\frac{T_{ij}}{C_i}$, and $U_i=\frac{I_i}{C_i}$. We are interested in the behavior of solutions of the above model near its equilibrium points. For $U_i=0$, the equilibrium points are obtained by solving the set of algebraic equations
%\[
%-b_iu_i+\sum_{j=1}^{n}A_{ij}g_j(u_j)=0,\quad i=1,2,\dots,n .
%\]
%In some applications the external inputs $U_i$ are held constant, that is, $U_i=c_i=\mb{constant}$. In such a case the equilibrium points are obtained by solving for $u_i$ the set of the following algebraic equations
%\[
%-b_iu_i+\sum_{j=1}^{n}A_{ij}g_j(u_j) + c_i=0,\quad i=1,2,\dots,n.
%\]
%However, we can always transform the algebraic equations with non-zero constant external inputs into the algebraic equations where the external inputs are set to zero. Indeed, let for $j=1,2,\ldots,n$,
%\[
%u_j=y_j+\frac{c_j}{b_j},
%\]
%and
%\[
%G_j(y_j)=g_j\left(y_j+\frac{c_j}{b_j}\right).
%\]
%Substituting the above into~(\ref{eq:newfoofhof}) and assuming that $U_i=c_i$ yields
%\begin{equation}
%\dot{y}_i=-b_iy_i+\sum_{j=1}^{n}A_{ij}G_j(y_j),\quad i=1,2,\ldots,n.
%\label{eq:afttranho}
%\end{equation}
%The equilibrium points of the above system of differential equations are obtained by solving the following algebraic equations,
%\begin{eqnarray*}
%-b_iy_i+\sum_{j=1}^{n}A_{ij}G_j(y_j)
%&=&-b_iy_i+\sum_{j=1}^{n}A_{ij}g_j\left(y_j + \frac{c_j}{b_j}\right)\\
%&=&0,\quad i=1,2,\ldots,n.
%\end{eqnarray*}
%Thus, after the transformation, the algebraic equations with external constant inputs used for computing the equilibrium points assume the same form as the algebraic equations corresponding to zero external inputs. Hence, without loss of generality, we can assume in our further analysis that the external inputs are zero and that the network is modeled by the equations,
%\begin{equation}
%\dot{u}_i=-b_iu_i+\sum_{j=1}^{n}A_{ij}g_j(u_j),\quad i=1,2,\ldots,n,
%\label{eq:newfoofhouzer}
%\end{equation}
%
%
%Throughout, we assume that an isolated equilibrium point of interest is located at the origin of the state-space $\mb{\Bb R}^n$. If this is not the case, then by an appropriate change of variables we can arrange for the equilibrium point of interest to be translated to the origin of $\mb{\Bb R}^n$. Let $\mb{\bm$u$}=\left[\begin{array}{ccc}
%u_1&\cdots &u_n
%\end{array}\right]^{\top} \in \mb{\Bb R}^n$ and suppose $\mb{\bm$u$}^{*} \ne \mb{\bm$0$}$ is an isolated equilibrium point of interest of~(\ref{eq:newfoofhouzer}). Let for $i=1,2,\ldots,n$,
%\begin{eqnarray*}
%x_i &=& u_i-u^{*}_i,\\
%h_i(u_i) &=& h_i(x_i + u^{*}_i)=g_i(u_i)-g_i(u^{*}_i).
%\end{eqnarray*}
%Substituting the above into~(\ref{eq:newfoofhouzer}) yields
%\begin{eqnarray*}
%\dot{x}_i &=&-b_ix_i - b_iu_i^{*} + \sum_{j=1}^{n}A_{ij}\left(h_j(x_j+u^{*}_j)+g_j(u^{*}_j) \right)\\
%&=&-b_ix_i+\sum_{j=1}^{n}A_{ij}h_j(x_j+u^{*}_j)+\left\{-b_iu^{*}_i+\sum_{j=1}^{n}A_{ij}g_j(u^{*}_j)\right\}.
%\end{eqnarray*}
%Note that, by the definition of the equilibrium point, the expression in the curly brackets vanishes. Thus we have
%\begin{equation}
%\dot{x}_i =-b_ix_i+ \sum_{j=1}^{n}A_{ij}h_j(x_j+u^{*}_j), \quad i=1,2,\ldots,n.
%\label{eq:dendeschop}
%\end{equation}
%Let
%\[
%G_j(x_j)=h_j(x_j+u^{*}_j).
%\]
%Then, differential equations~(\ref{eq:dendeschop}) can be written as
%\[
%\dot{x}_i =-b_ix_i+ \sum_{j=1}^{n}A_{ij}G_j(x_j), \quad i=1,2,\ldots,n,
%%\label{eq:xcoorofhosin}
%\]
%where $G_j:\mb{\Bb R} \to (c_1,c_2) \subseteq (-2,2)$, $c_1 < 0 < c_2$, $G_j(0)=0$, $G_j$ is continuously differentiable and strictly monotonically increasing in $x_j$. Furthermore $x_jG_j(x_j) >0$ for all $x_j \ne 0$. We now represent the above model of the Hopfield neural network in the form
%\[
%\dot{x}_i =-b_ix_i+A_{ii}G_i(x_i)+ \sum_{{j=1}\atop{j \ne i}}^{n}A_{ij}G_j(x_j), \quad i=1,2,\ldots,n,
%\]
%We shall view the terms
%\[
%\sum_{{j=1}\atop{j \ne i}}^{n}A_{ij}G_j(x_j)
%\]
%as the interconnections between the $i$-th neuron and the remaining neurons in the network. We shall refer to the equations
%\begin{equation}
%\dot{x}_i =-b_ix_i+A_{ii}G_i(x_i),\quad i=1,2,\ldots,n
%\label{eq:ensosinneu}
%\end{equation}
%as models of free isolated neurons. In what follows we analyze stability properties of an isolated single neuron.
%
%\subsection{Single Neuron Stability Analysis}
%A single neuron in the Hopfield network, when the external input is set to zero, is modeled by a first order differential equation
%\begin{equation}
%\fbox{$\displaystyle \dot{x}_i =-b_ix_i+A_{ii}G_i(x_i) $}
%\label{eq:sinneueq}
%\end{equation}
%Consider now, as in Michel et al.~\cite{michel89}, the following positive definite Lyapunov function candidate
%\[
%V_i(x_i)=\frac{1}{2} x_{i}^2.
%\]
%Its time derivative evaluated along the trajectories of the single neuron model is
%\begin{eqnarray*}
%\dot{V}_i (x_i(t)) &=& x_i \dot{x}_i\\
%&=& x_i \left( -b_ix_i+A_{ii}G_i(x_i) \right).
%\end{eqnarray*}
%If there exists $r_i > 0$ such that for all $x_i \in (-r_i,r_i)$
%\[
%\left\{ \begin{array}{ccc}
%-b_ix_i+A_{ii}G_i(x_i)>0 & \mbox{if} & x_i < 0\\
%-b_ix_i+A_{ii}G_i(x_i)=0 & \mbox{if} & x_i=0\\
%-b_ix_i+A_{ii}G_i(x_i)<0 & \mbox{if} & x_i > 0,
%\end{array} \right.
%\]
%then for $x_i \ne 0$ and $x_i \in (-r_i,r_i)$,
%\[
%\dot{V}_i(x_i(t)) < 0
%\]
%and $\dot{V}_i(0)=0$. Hence the null solution $x_i=0$ is uniformly asymptotically stable. Figure~\ref{fig:st_nen} illustrates this case.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{st_nen.eps}
%%\vspace{3in}
%%\centerbmp{3.33in}{2.5in}{c:/mybook/st_nen.bmp}
%%\hskip1in\special{bmp:c:/mybook/st_nen.bmp x=4in y=3in}
%\caption{Characteristics of a single neuron model whose null solution is uniformly asymptotically stable.}
%\label{fig:st_nen}
%\end{figure}
%Now, if there exists $r_i > 0$ such that for $x_i \in (-r_i,r_i)$ we have
%\[
%\left\{ \begin{array}{ccc}
%-b_ix_i+A_{ii}G_i(x_i)<0 & \mbox{if} & x_i < 0\\
%-b_ix_i+A_{ii}G_i(x_i)=0 & \mbox{if} & x_i=0\\
%-b_ix_i+A_{ii}G_i(x_i)>0 & \mbox{if} & x_i > 0,
%\end{array} \right.
%\]
%then in addition to $V_i$ being positive definite,
%\[
%\dot{V}_i(x_i(t)) > 0
%\]
%for $x_i \ne 0$ and $x_i \in (-r_i,r_i)$, and $\dot{V}_i(0)=0$. Thus the null solution $x_i=0$ is unstable. Figure~\ref{fig:ust_nen} illustrates the case when the null solution of a single neuron is unstable.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{ust_nen.eps}
%%\vspace{3in}
%%\centerbmp{3.33in}{2.5in}{c:/mybook/ust_nen.bmp}
%%\hskip1in\special{bmp:c:/mybook/ust_nen.bmp x=4in y=3in}
%\caption{Characteristics of single neuron model whose null solution is unstable.}
%\label{fig:ust_nen}
%\end{figure}
%Assume now that for some $r_i>0$ there exist constants $\sigma _{i1} > 0$ and $\sigma _{i2} > 0$ such that for all $x_i \in (-r_i,r_i)$, $x_i \ne 0$, we have
%\begin{equation}
%\sigma_{i1}x_i^2 < x_i G_i (x_i) < \sigma_{i2} x_{i2}^2,
%\label{eq:sigoandt}
%\end{equation}
%or equivalently
%\[
%\sigma_{i1} < \frac{G_i (x_i)}{x_i} < \sigma_{i2}.
%\]
%If in addition to the sufficient conditions for asymptotic stability we also have
%\begin{equation}
%-b_i+A_{ii} \delta_i < 0,
%\label{eq:sing_neur}
%\end{equation}
%where
%\begin{equation}
%\delta_i = \left\{ \begin{array}{ccc}
%\sigma_{i1} & \mbox{if} & A_{ii} < 0\\
%\sigma_{i2} & \mbox{if} & A_{ii} > 0,
%\end{array} \right.
%\label{eq:defofdeli}
%\end{equation}
%then the null solution of the single neuron model is exponentially stable. Indeed, for $x_i \ne 0$ we have
%\begin{eqnarray*}
%\dot{V}_i(x_i(t)) &=& x_i\dot{x}_i\\
%&=&  x_i \left( -b_ix_i+A_{ii}G_i(x_i) \right)\\
%&=& \left( -b_i+A_{ii} \frac{G_i(x_i)}{x_i} \right)x_i^2\\
%&\le& \left(-b_i + A_{ii} \delta_i \right) x_i^2\\
%&=& 2\left( -b_i + A_{ii} \delta_i \right)V_i(x_i(t)).
%\end{eqnarray*}
%Applying Theorem~\ref{thm:comparison} yields
%\begin{eqnarray*}
%V_i(x_i(t)) &=& \frac{x_i^2(t)}{2}\\
%&\le& \exp \left\{2\left(-b_i + A_{ii} \delta_i\right)(t-t_0) \right\} V_i(x_i(t_0))\\
%&=& \exp \left\{2\left(-b_i + A_{ii} \delta_i\right)(t-t_0) \right\}\frac{x_i^2(t_0)}{2}.
%\end{eqnarray*}
%Hence
%\[
%|x_i(t)| \le \exp \left\{\left(-b_i + A_{ii} \delta_i\right)(t-t_0) \right\}|x_i(t_0)|.
%\]
%From the above we conclude that the null solution is not only asymptotically stable but even more, it is exponentially stable.
%
%\subsection{Stability Analysis of the Network}
%Lyapunov's second method in conjunction with interconnection information will now be used in the stability analysis of the Hopfield neural model viewed as an interconnected dynamical system. Applications of Lyapunov's second method to a general class of interconnected systems were investigated by Bailey~\cite{bailey66}. Here, we use Lyapunov's second method to investigate stability properties of the null solution of the neural network described by
%\begin{equation}
%\fbox{$\displaystyle \dot{x}_i=-b_ix_i+\sum_{j=1}^{n}A_{ij}G_j(x_j), \quad i=1,2,\ldots,n $}
%\label{eq:hopasintersy}
%\end{equation}
%We impose the following assumptions on the interconnections:\\
%{\bf A1}: There exist positive constants $r_i$, $r_j$, and $a_{ij}$ such that for $|x_i| < r_i$ and $|x_j| < r_j$, $i,j=1,2,\ldots,n$, we have
%\[
%x_iA_{ij}G_j(x_j) \le |x_i|a_{ij}|x_j|.
%\]
%{\bf A2}: Let $\mb{\bm$S$}$ be an $n\times n$ real matrix whose elements are
%\[
%s_{ij}=\left\{\begin{array}{ccc}
%\alpha_i \left(-b_i + a_{ii} \right) & \mbox{if} & i=j\\
%\frac{1}{2}\left(\alpha_i a_{ij} + \alpha_j a_{ji}\right) & \mb{if}& i \ne j.
%\end{array}\right.
%\]
%We assume that there exists a vector $\mb{\bm$\alpha$}=\left[\begin{array}{ccc}
%\alpha_1 & \cdots &\alpha_n
%\end{array}\right]^{\top}$, where $\alpha_i >0$ for $i=1,2,\ldots,n$ such that the matrix $\mb{\bm$S$}$ is negative definite.
%
%We now state and prove a theorem concerned with the exponential stability of the neural network model. This theorem can be found in Michel et al.~\cite{michel89}.
%\begin{thm}\label{thm:expstann}
%The null solution of the neural network model~(\ref{eq:hopasintersy}) is uniformly exponentially stable if the assumptions {\bf A1} and {\bf A2} are satisfied.
%\end{thm}
%{\bf Proof}\quad Consider the following Lyapunov function candidate
%\begin{eqnarray*}
%V(\mb{\bm$x$}) &=& \frac{1}{2}\sum_{i=1}^{n}\alpha_i x_i^2\\
%&=&\frac{1}{2}\mb{\bm$x$}^{\top} \left[\begin{array}{cccc}
%\alpha_1 & {} & \ldots& 0\\
%{} & \alpha_2 & {} & {}\\
%\vdots & {} & \ddots & \vdots\\
%0 & {} & \ldots & \alpha_n
%\end{array}\right] \mb{\bm$x$},
%\end{eqnarray*}
%where $\alpha_i > 0$, $i=1,2,\ldots,n$ are defined in assumption {\bf A2}. Note that $V$ is real symmetric positive definite. The time derivative of $V$ evaluated along the trajectories of the neural model~(\ref{eq:hopasintersy}) is
%\begin{eqnarray*}
%\dot{V}(\mb{\bm$x$}(t)) &=& \sum_{i=1}^{n} \alpha_i x_i \dot{x}_i\\
%&=& \sum_{i=1}^{n} \alpha_i x_i \left( -b_ix_i + \sum_{j=1}^{n} A_{ij} G_j(x_j)\right)\\
%&\le& \sum_{i=1}^{n} \alpha_i \left( -b_i x_i^2 + |x_i|\sum_{j=1}^{n}a_{ij}|x_j| \right)\\
%&=& |\mb{\bm$x$}|^{\top} \mb{\bm$R$} |\mb{\bm$x$}|,
%\end{eqnarray*}
%where $|\mb{\bm$x$}|=\left[\begin{array}{cccc}
%|x_1|&|x_2|&\cdots &|x_n|
%\end{array}\right]^{\top}$, and the matrix $\mb{\bm$R$} \in \mb{\Bb R}^{n \times n}$ has its elements defined as
%\[
%r_{ij}= \left\{\begin{array}{ccc}
%\alpha_i \left( -b_i + a_{ii} \right) & \mb{if} & i=j\\
%\alpha_i a_{ij} & \mb{if} & i \ne j.
%\end{array}\right.
%\]
%Note that
%\begin{eqnarray*}
%|\mb{\bm$x$}|^{\top} \mb{\bm$R$} |\mb{\bm$x$}| &=& |\mb{\bm$x$}|^{\top} \left( \frac{\mb{\bm$R$}+\mb{\bm$R$}^{\top}}{2} \right)|\mb{\bm$x$}|\\
%&=& |\mb{\bm$x$}|^{\top} \mb{\bm$S$}|\mb{\bm$x$}|\\
%&\le& \lambda_{\max}(\mb{\bm$S$})\|\mb{\bm$x$}\|^2,
%\end{eqnarray*}
%where the matrix $\mb{\bm$S$}$ is defined in {\bf A2}. By the assumption {\bf A2} the matrix $\mb{\bm$S$}$ is real symmetric negative definite. Hence $\lambda_{\max}(\mb{\bm$S$})$ is real and negative. Thus
%\[
%\dot{V}(\mb{\bm$x$}(t)) \le  \lambda_{\max}(\mb{\bm$S$})\|\mb{\bm$x$}(t)\|^2 < 0
%\]
%in a neighborhood of $\mb{\bm$x$}=\mb{\bm$0$}$ of radius
%\[
%r=\min_{1 \le i \le n}r_i.
%\]
%By Proposition~\ref{pro:exponsaofz}, the null solution is uniformly exponentially stable.\\
%\hspace*{5.5in} $\Box$
%
%In our stability analysis of the Hopfield neural network we do not require that the parameters $A_{ij}$, and hence the interconnections $T_{ij}$, be symmetric. The symmetry requirement on the interconnections was usually imposed in the early stability analyses of this network.
%
%Note also that a necessary condition for the matrix $\mb{\bm$S$}$ to be negative definite is that its diagonal elements be negative, that is,
%\begin{equation}
%-b_i+a_{ii} < 0, \quad i=1,2,\ldots,n.
%\label{eq:diagelofsne}
%\end{equation}
%However this condition is not sufficient for the matrix $\b S$ to be negative definite. The above necessary condition, though, can shed some light on the influence of the stability properties of the individual neurons on the stability of the whole network. To show how the stability of the individual neurons may affect the stability of the whole network, we set $i=j$ in {\bf A1} to obtain
%\begin{equation}
%x_iA_{ii}G_i(x_i) \le a_{ii} |x_i|^2.
%\label{eq:afteraone}
%\end{equation}
%Applying now the sector condition, $\sigma_{i1}x_i^2 < x_i G_i (x_i) < \sigma_{i2} x_{i}^2$, to $x_iA_{ii}G_i(x_i)$ gives
%\begin{equation}
%x_iA_{ii}G_i(x_i) \le A_{ii}\delta_i |x_i|^2,
%\label{eq:afterseccon}
%\end{equation}
%where $\delta_i$ is defined in~(\ref{eq:defofdeli}). Thus the sector condition~(\ref{eq:sigoandt}) is sufficient for~(\ref{eq:afteraone}) to hold and we can write
%\[
%a_{ii}=A_{ii} \delta_i,
%\]
%Substituting the above into~(\ref{eq:diagelofsne}) yields
%\[
%-b_i + A_{ii} \delta_i < 0,\quad i=1,2,\ldots,n,
%\]
%which are sufficient conditions for the exponential stability of the individual neurons. Thus, if we first find constants $a_{ii}$, $i=1,2,\ldots,n$, and then use them to construct the matrix $\b S$, then for so formed matrix $\b S$ to be negative definite it is necessary that the above sufficient conditions for the exponential stability of the individual neurons hold.
%
%We now present a corollary to Theorem~\ref{thm:expstann}. Before stating the corollary, we need to perform some algebraic manipulations. We first combine {\bf A2} with the sector condition $\sigma_{i1}x_i^2 < x_i G_i (x_i) < \sigma_{i2} x_{i}^2$. It follows from the sector condition that for $x_j \ne 0$
%\[
%\left\arrowvert \frac{G_j(x_j)}{x_j} \right\arrowvert < \sigma_{j2}.
%\]
%Applying the above when $i \ne j$ yields
%\begin{eqnarray*}
%x_i A_{ij} G_j(x_j) &\le& |A_{ij}||x_iG_j(x_j)|\\
%&=&|A_{ij}|\left\arrowvert x_i \frac{G_j(x_j)}{x_j} x_j\right\arrowvert\\
%&\le& |A_{ij}||x_i||x_j|\, \sigma_{j2}.
%\end{eqnarray*}
%When $i=j$, we obtain
%\[
%x_iA_{ii}G_i(x_i) \le \delta_i A_{ii} x_i^2,
%\]
%where $\delta_i$ is defined as above. We use the results of the above manipulations in the following assumption:\\
%{\bf A3:} There exist positive constants $r_i$ and $r_j$ such that for $|x_i| < r_i$ and $|x_j| < r_j$, the interconnections of the neural network model~(\ref{eq:hopasintersy}) satisfy the conditions
%\begin{eqnarray*}
%x_i A_{ii} G_i(x_i) &\le& \delta_i A_{ii} x_i^2 \quad \mb{for} \quad i=j\\
%x_i A_{ij} G_j (x_j) &\le& |x_i||A_{ij}||x_j| \sigma_{j2} \quad \mb{for} \quad i \ne j.
%\end{eqnarray*}
%We shall use another assumption.\\
%{\bf A4:} There exists a vector $\mb{\bm$\alpha$}=\left[\begin{array}{ccc}
%\alpha_1 & \cdots &\alpha_n
%\end{array}\right]^{\top} \in \mb{\Bb R}^n$, $\alpha_i > 0$, such that the matrix $\mb{\bm$S$} \in \mb{\Bb R}^{n \times n}$ is negative definite, where the elements of $\mb{\bm$S$}$ are defined as
%\[
%s_{ij} = \left\{ \begin{array}{ccc}
%\alpha _i (-b_i + \delta _i A_{ii} )& \mb{if} & i=j\\
%\frac{1}{2} \left(\alpha _i |A_{ij}| \sigma_{j2} + \alpha _j |A_{ji}| \sigma_{i2} \right) & \mb{if} & i \ne j.
%\end{array} \right.
%\]
%We are now ready to state and prove a corollary to Theorem~\ref{thm:expstann}. This result also comes from Michel et al.~\cite{michel89}.
%\begin{corollary}\label{cor:exphop}
%The null solution of the neural network model~(\ref{eq:hopasintersy}) is uniformly exponentially stable if the assumptions {\bf A3} and {\bf A4} are satisfied.
%\end{corollary}
%{\bf Proof}\quad Consider the following quadratic positive definite Lyapunov function candidate
%\[
%V(\mb{\bm$x$})=\frac{1}{2} \sum_{i=1}^{n} \alpha_i x_i^2,
%\]
%where $\alpha_i > 0$, $i=1,2,\ldots,n$, are defined in {\bf A4}. Taking into account the assumptions {\bf A3} and {\bf A4}, we evaluate the time derivative of $V$ along the trajectories of the neural network~(\ref{eq:hopasintersy}) to get
%\begin{eqnarray*}
%\dot{V}(\mb{\bm$x$}(t)) &=& \sum_{i=1}^{n} \alpha_i x_i \left(-b_ix_i+\sum_{j=1}^{n}A_{ij}G_j(x_j) \right)\\
%&=& \sum_{i=1}^{n} \alpha_i x_i \left(-b_i x_i+ A_{ii} G_i(x_i) + \sum_{{j=1}\atop{j \ne i}}^{n}A_{ij} G_j(x_j) \right)\\
%&\le& \sum_{i=1}^{n} \alpha_i \left( \left(-b_i + \delta_i A_{ii} \right) x_i^2 +\sum_{{j=1}\atop{j \ne i}}^{n} |A_{ij}||x_i||x_j| \, \sigma_{j2} \right)\\
%&=& \frac{1}{2}|\mb{\bm$x$}|^{\top}\mb{\bm$S$}|\mb{\bm$x$}|\\
%&<& 0 \quad \mb{for} \quad \mb{\bm$x$} \ne \mb{\bm$0$}.
%\end{eqnarray*}
%Using the above, the fact that $V$ is a quadratic function, and Proposition~\ref{pro:exponsaofz}, we conclude that the null solution is uniformly exponentially stable.\\
%\hspace*{5.5in} $\Box$
%
%
%We illustrate the above corollary with a numerical example involving a neural network consisting of two neurons.
%\begin{example}\label{ex:twodhopnet}\rm
%The neural net that we analyze in this example is described by the equations,
%\begin{equation}
%\left[\begin{array}{c}
%\dot{u}_1\\
%\dot{u}_2
%\end{array}\right]=\left[\begin{array}{c}
%-1.4 u_1+1.2g_2(u_2)\\
%-1.1 u_2+g_1(u_1)
%\end{array}\right],
%\label{eq:netinucoord}
%\end{equation}
%where
%\begin{equation}
%g_i(u_i)=\frac{2}{\pi}\arctan\left(\lambda\frac{\pi}{2} u_i\right),\quad i=1,2.
%\label{eq:actnonho}
%\end{equation}
%The value of the parameter $\lambda$ was set to be equal to $1.4$. In Figure~\ref{fig:atan_ac}, we show a plot of the above activation nonlinearity.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{atan_ac.eps}
%%\vspace{2.5in}
%%\hskip1.75in\special{eps:c:/mybook/atan_ac.eps y=2.5in}
%\caption{A plot of the activation nonlinearity~(\ref{eq:actnonho}).}
%\label{fig:atan_ac}
%\end{figure}
%The network has three equilibrium states. They are:
%\[
%\b u^{(1)}=\left[\begin{array}{c}
%-0.4098\\
%-0.4245
%\end{array}\right], \quad \b u^{(2)}=\left[\begin{array}{c}
%0.4098\\
%0.4245
%\end{array}\right],\quad\mb{and}\quad \b u^{(3)}=\left[\begin{array}{c}
%0\\
%0
%\end{array}\right].
%\]
%When computing the equilibrium points, we began by finding the first coordinates of them. This can be done, for example, by combining the two algebraic equations,
%\[
%-1.4 u_1+1.2g_2(u_2)=0\quad\mb{and}\quad -1.1 u_2+g_1(u_1)=0
%\]
%and then solving them for $u_1$. A plot of the function obtained by substituting the second of the above equations into the first one is shown in Figure~\ref{fig:equi_comp}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{equi_comp.eps}
%%\vspace{2.5in}
%%\hskip1.75in\special{eps:c:/mybook/equi_comp.eps y=2.5in}
%\caption{A plot of the function used to find the first coordinates of the equilibrium points of the neural net of Example~\ref{ex:twodhopnet}.}
%\label{fig:equi_comp}
%\end{figure}
%We can use Lyapunov's first method, as we did in Example~\ref{ex:hopf_lyap}, to show that the equilibrium states $\b u^{(1)}$ and $\b u^{(2)}$ are asymptotically stable, while the equilibrium $\b u^{(3)}$ is unstable.
%
%We now apply Corollary~\ref{cor:exphop} to show that the equilibrium states $\b u^{(1)}$ and $\b u^{(2)}$ are in fact uniformly exponentially stable. We analyze the equilibrium state $\b u^{(1)}$; the analysis of $\b u^{(2)}$ is identical. We first translate the equilibrium state of interest to the origin of $\mb{\Bb R}^2$ using the translation,
%\[
%\b x=\b u-\b u^{(1)}
%\]
%Let
%\[
%G_i(x_i)=g_i\left(x_i+u^{(1)}_i\right)-g_i\left(u^{(1)}_i \right),\quad i=1,2.
%\]
%Taking into account the above relation, the neural network model in the new coordinates takes the form,
%\begin{equation}
%\left[\begin{array}{c}
%\dot{x}_1\\
%\dot{x}_2
%\end{array}\right]=\left[\begin{array}{c}
%-1.4x_1+1.2 G_2(x_2)\\
%-1.1x_2+G_1(x_1)
%\end{array}\right].
%\label{eqn:hoftdinxco}
%\end{equation}
%The equilibrium state $\b x_e=\b 0$ of the above net corresponds to the equilibrium state $\b u_e=\b u^{(1)}$ of the net in the $u$-coordinates. To be able to apply Corollary~\ref{cor:exphop}, we need to find $\sigma_{i2}$ for $i=1,2$ such that
%\[
%0<\frac{G_i(x_i)}{x_i}<\sigma_{i2}\quad\mb{for}\quad |x_i|<r,
%\]
%where $r>0$. In this example, we took $r=\min\left\{\left|u^{(1)}_1\right|,\left|u^{(1)}_2\right| \right\}=0.4098$. Graphs of $\frac{G_i(x_i)}{x_i}$ versus $x_i$, for i=1,2, are shown in Figure~\ref{fig:two_gix}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{two_gix.eps}
%%\vspace{3.5in}
%%\hskip1.2in\special{eps:c:/mybook/two_gix.eps y=3.5in}
%\caption{Plots of $G_i(x_i)/x_i$ in Example~\ref{ex:twodhopnet}.}
%\label{fig:two_gix}
%\end{figure}
%We calculated $\sigma_{12}=1.1394$ and $\sigma_{22}=1.1165$. We next form the matrix $\b S\in \mb{\Bb R}^{2\times 2}$ defined in {\bf A4}. Selecting $\alpha_1=\alpha_2=1$ and taking into account that $A_{11}=0$, $A_{12}=1.4$, and  $A_{21}=1$, we obtain
%\begin{eqnarray*}
%\b S&=&\left[\begin{array}{cc}
%- 1.4\alpha_1 & \left(\alpha_1\left|A_{12}\right|\sigma_{22} + \alpha_2\left|A_{21}\right|\sigma_{12}\right)/2\\
%\left(\alpha_2\left|A_{21}\right|\sigma_{12} + \alpha_1\left|A_{12}\right|\sigma_{22}\right)/2 & -1.1 \alpha_2
%\end{array}\right]\\
%&=&\left[\begin{array}{cc}
%-1.4000  &  1.2396\\
%1.2396  & -1.1000
%\end{array}\right].
%\end{eqnarray*}
%The eigenvalues of S are located at $-2.4987$ and $-0.0013$. Hence, the matrix $\b S$ is negative definite, and therefore the equilibrium state $\b x_e=\b 0$ of the net~(\ref{eqn:hoftdinxco}), and hence $\b u_e=\b u^{(1)}$ of the net~(\ref{eq:netinucoord}), is uniformly exponentially stable. In Figure~\ref{fig:mod_hop}, we show a state plane portrait of the neural network modeled by~(\ref{eq:netinucoord}).
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{mod_hop.eps}
%%\vspace{3.5in}
%%\hskip1.35in\special{eps:c:/mybook/mod_hop.eps y=3.5in}
%\caption{Phase plane portrait of the neural net of Example~\ref{ex:twodhopnet}.}
%\label{fig:mod_hop}
%\end{figure}
%\end{example}\markex
%
%\section{Brain-State-in-a-Box (BSB) Models}
%\subsection{Associative Memories}
%\label{subsec:assoscmem}
%In this section, we analyze a class of neural network models that can be used to construct the so called associative memories. We distinguish between two types of association: autoassociation and heteroassociation. In autoassociation, after storing a set of patterns by a neural network, where patterns are usually represented as vectors, a distorted (noisy) pattern of a stored pattern is subsequently presented to the network. The task of the neural network is to retrieve (recall) the original stored pattern from the noisy pattern. In heteroassociation an arbitrary set of input patterns is paired with another arbitrary set of output patterns. Operation of an associative memory is characterized by two stages: storage phase, where patterns are being stored by the neural network, and recall phase, where memorized patterns are being retrieved in response to a noisy pattern presented to the network. Figure~\ref{fig:patt_assoc} illustrates the recall phase. In this figure, the input pattern represents a noisy version of stored pattern $\b x_j$. The pattern associator produces the output pattern $\b y$ in response to the input pattern $\b x$. For a perfect recall in an autoassociative memory, we should have $\b y=\b x_j$. When $\b y \ne \b x_j$, the associative memory is said to have made a recall error (or error in recall).
%\begin{figure}
%\setlength{\unitlength}{.1in}
%\begin{center}
%\begin{picture}(60,20)
%\thicklines
%\put(11,9){\vector(1,0){11}}
%\put(11,10){Input pattern}
%\put(16,7){$\b x$}
%\put(22,5){\framebox(10,8){\shortstack{Pattern\\\\associator}}}
%\put(32,9){\vector(1,0){12}}
%\put(33,10){Output pattern}
%\put(37,7){$\b y$}
%\end{picture}
%\end{center}
%\caption{Operation of pattern associator.}
%\label{fig:patt_assoc}
%\end{figure}
%We now analyze a static linear autoassociator. For the sake of simplicity, we assume that the given patterns, $\b v^{(j)}$, $j=1,\ldots, r$, are orthogonal, that is, $\left\langle \b v^{(j)},\b v^{(k)}\right\rangle =0$ for $k\ne j$. Suppose that we form the weight matrix as the sum of the outer products of the given patterns, that is,
%\[
%\b W = \sum_{j=1}^r \b v^{(j)} \b v^{(j)T} =  \b V\b V^{\top},
%\]
%where $\b V=\left[\begin{array}{cccc}
%\b v^{(1)} & \b v^{(2)} & \cdots & \b v^{(r)}
%\end{array}\right]$. Let $\b v$ be one of the given patterns. We present $\b v$, to the linear associator. Then, because the patterns $\b v^{(j)}$, $j=1,\ldots, r$, are mutually orthogonal, the output of the associator is
%\[
%\b W\b v = \sum_{j=1}^r \b v^{(j)} \b v^{(j)T} \b v = \b v \b v^{\top} \b v = c \b v,
%\]
%where $c=\b v^{\top} \b v$. Thus, we obtained back the presented pattern.
%
%We will now show that our linear associator can reconstruct the missing part of a stored pattern. Assume that
%\[
%\b v=\b v{'}+\b v{''},
%\]
%that is, the given pattern $\b v$ is composed of two parts. Suppose that $\b v{'}$ and $\b v{''}$ are orthogonal to each other and that $\b v{'}$ is also orthogonal  to the remaining patterns in the set. Taking this into account while presenting $\b v{'}$ to our linear associator gives
%\begin{eqnarray*}
%\b W\b v{'} &=& \sum_{j=1}^r \b v^{(j)} \b v^{(j)T} \b v{'}\\
%&=& \b v\b v^{\top}\b v{'}\\
%&=& \left( \b v{'}+\b v{''}\right)\left(\b v{'}+\b v{''} \right)^{\top}\b v{'}\\
%&=& \left( \b v{'}+\b v{''}\right)\b v{'}^{\top}\b v{'}\\
%&=& d \b v,
%\end{eqnarray*}
%where $d=\b v{'}^{\top}\b v{'}$. Thus the linear associator was able to reconstruct the missing part of the presented pattern. In the above discussion, we assumed that the patterns are mutually orthogonal. If this is not the case, the output of the associator will contain an extra component called ``cross-talk''. This linear static associator has very limited capability. However, we shall use the idea of the outer product to construct neural associative memories.
%
%
%A neural network that can be used as an associative memory is the Brain-State-in-a-Box~(BSB) neural model, which was proposed by Anderson, Silverstein, Ritz, and Jones in 1977~\cite{anderson77}. It can recognize a pattern from its given noisy version. For this reason the BSB neural model is often referred to as an associative memory. The BSB neural model is a discrete dynamical system whose dynamics are described by the difference equation,
%\begin{equation}
%\b x(k+1)=\b g\left(\b x(k)+\alpha \b W\b x(k)\right)
%\label{eq:classic_bsb}
%\end{equation}
%subject to the initial condition
%\[
%\b x(0)=\b x_0,
%\]
%where $\b x(k) \in \mb{\Bb R}^n$ is the state of the BSB network at time $k$, $\alpha>0$ is the step size, and $\b W \in \mb{\Bb R}^{n \times n}$ is a symmetric weight matrix determined during the training procedure. The function $\b g:\mb{\Bb R}^n \to \mb{\Bb R}^n$ is a vector valued function whose i-th component $(\b g)_i$ is defined as:
%\[
%(\b g(\b x))_i=\left\{ \begin{array}{lll}
%1 & \mbox{if}& \quad x_i \ge 1 \\
%x_i & \mbox{if} & \quad -1<x_i<1 \\
%-1 & \mbox{if}& \quad x_i \le -1.
%\end{array}
%\right.
%\]
%The above nonlinear activation function is sometimes referred to as the linear saturating nonlinearity. A plot of such a  function is shown in Figure~\ref{fig:lin_sat}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{sat_non.ps}
%%\vspace{3in}
%%\hskip1.5in\special{ps:c:/mybook/sat_non.ps y=3in}
%\caption{A plot of linear saturating nonlinearity used in the BSB model as an activation function.}
%\label{fig:lin_sat}
%\end{figure}
%The BSB model is named from the fact that the state of the network trajectories are constrained to be in the hypercube $H_n=\left[\begin{array}{cc}
%-1, & 1
%\end{array}\right]^n$. In Figure~\ref{fig:bsb_tr2}, we show a state plane portrait of the BSB model~(\ref{eq:classic_bsb}) with the weight matrix
%\begin{equation}
%\b W=\left[\begin{array}{cc}
%7 & -3\\
%-3 & 10
%\end{array}\right]
%\label{eq:bsbweighton}
%\end{equation}
%and the step size is $\alpha=0.1$. The state space of this two-dimensional neural net is the square $H_2=\left[\begin{array}{cc}
%-1, & 1
%\end{array}\right]^2$.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{bsb_tr2.eps}
%%\vspace{3in}
%%\hskip1.75in\special{eps:c:/mybook/bsb_tr2.eps y=3in}
%\caption{Trajectories of the BSB model with the weight matrix $\b W$ given by~(\ref{eq:bsbweighton}).}
%\label{fig:bsb_tr2}
%\end{figure}
%Each vertex of the square has its basin of attraction. In this example all four vertices are asymptotically stable equilibrium states of the BSB neural net. In Figure~\ref{fig:bsb_tr1}, we show a state plane portrait of the BSB model~(\ref{eq:classic_bsb}) with the weight matrix
%\begin{equation}
%\b W=\left[\begin{array}{cc}
%5 & -3\\
%-3 & 2
%\end{array}\right]
%\label{eq:bsbweighttwo}
%\end{equation}
%and the step size $\alpha=0.1$. In this example only two vertices are asymptotically stable equilibrium states of the BSB neural net.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{bsb_tr1.eps}
%%\vspace{3in}
%%\hskip1.75in\special{eps:c:/mybook/bsb_tr1.eps y=3in}
%\caption{Trajectories of the BSB model with the weight matrix $\b W$ given by~(\ref{eq:bsbweighttwo}).}
%\label{fig:bsb_tr1}
%\end{figure}
%When the BSB model is used as an associative memory, asymptotically stable equilibrium points of the model represent stored patterns. Noisy patterns are represented by points in the basins of attraction of the stable equilibrium points. Thus, suppose we are given a noisy, or imperfect, pattern represented by an initial condition $\b x_0$ that is located in a basin of attraction of a stable equilibrium point corresponding to the correct pattern. Then, the corresponding trajectory of the BSB model evolves in time to this stable equilibrium point representing the correct stored pattern. An associative memory is able to recall all the information stored in the memory when it is excited with a sufficient large portion of that memory's information. This is in contrast to the address addressable memory (AAM) used in digital computers. In this type of memory, a desired set of information is recalled only when the entire correct address of the information is given. For this reason, we can view associative memories as content addressable memories (CAM).
%
%In summary, realization of associative memories with neural networks requires storage of a set of given patterns as asymptotically stable equilibrium points of the network. The points in basins of attraction about each stable equilibrium point correspond to noisy versions of the stored patterns. When a noisy version of a stored pattern is sufficiently close to it, the network trajectory should converge to the equilibrium point corresponding to the correct pattern. Neural network associative memories should possess the following useful properties:
%\begin{enumerate}
%\item \label{aemem:one} Each given pattern to be stored should be stored as an asymptotically stable equilibrium point of the network.
%\item \label{aemem:two} The number of asymptotically stable equilibrium points of the network that do not correspond to the given stored patterns, that is, the number of spurious states, should be minimal.
%\item \label{aemem:three} A nonsymmetric interconnection structure, which eases difficulties in implementing associative memory using neural networks.
%\item \label{aemem:four} The ability to control the extent of the basin of attraction about the equilibrium point corresponding to each stored pattern.
%\item \label{aemem:five} Learning capability, that is, the ability to add patterns to be stored as asymptotically stable equilibrium points to an existing set of stored patterns without affecting the existing equilibria.
%\item \label{aemem:six} Forgetting capability, that is, the ability to delete specified equilibria from a given set of stored equilibria without affecting the rest of the equilibria.
%\item \label{aemem:seven} A high storage and retrieval efficiency, that is, the ability to efficiently store and retrieve a large number, compared to the order $n$ of the network, of patterns.
%\end{enumerate}
%\subsection{Analysis of BSB Models}
%The BSB neural net has a symmetric weight matrix, which is not a desirable property of a neural associative memory. Hui and \.{Z}ak~\cite{hui92bsb} modified the BSB model to allow for nonsymmetric weight matrix as well as to better control the extent of basins of attraction. A modified form of the BSB neural model, which we refer to as the generalized Brain-State-in-a-Box~(gBSB) neural model, is well suited for implementing associative memories. Its dynamics are described by the equation
%\begin{equation}
%\b x(k+1)=\b g((\b I_n+\alpha \b W)\b x(k)+\alpha \b b),
%\label{eq:gbsbmo}
%\end{equation}
%where $\b I_n$ is the $n \times n$ identity matrix, $\b b \in \mb{\Bb R}^n$, and the weight matrix $\b W\in \mb{\Bb R}^{n \times n}$ need not be symmetric as in the BSB model. Thus, the original BSB model can be viewed as a special case of the gBSB model when the interconnection matrix $\b W$ is symmetric and $\b b=\b 0$. We add that the analysis of the BSB model behavior on the boundary regions of the hypercube $H_n$ can be reduced to the study of the reduced order gBSB type models rather than that of the BSB type.
%
%We now describe an algorithm for computing the weight matrix $\b W$ and the vector $\b b$ that result in the given patterns being stored as asymptotically stable vertices of the hypercube $H_n$. In further discussion, we use the following notation and definitions. Let
%\[
%\b L(\b x)=\left(\b I_n+\alpha \b W\right)\b x+\alpha \b b,
%\]
%and
%\[
%\b T(\b x)=\b g(\b L(\b x)).
%\]
%\begin{defn}
%A point $\b x \in H_n$ is an equilibrium point of $\b x(k+1)=\b T(\b x(k))$ if $\b x=\b T(\b x)$.
%\end{defn}\markdefn
%
%Note that an equilibrium point is a constant solution to the system $\b x(k+1)=\b T(\b x(k))$.
%
%Let $(\b L(\b x))_i$ be the $i$-th component of $\b L(\b x)$ and $v_i^{(j)}$ be the $i$-the component of the vector $\b v^{(j)}$. Let $\b v=\left[\begin{array}{ccc}
%v_1 & \cdots & v_n
%\end{array}\right]^{\top} \in \{-1, 1\}^n$, that is, $\b v$ is a vertex of $H_n$ because for $i=1,2,\ldots, n$, we have $v_i=\pm 1$. For a vertex $\b v$ of $H_n$ to be an equilibrium point of the gBSB net, we must have
%\[
%v_i=\left( \b g\left( \b L(\b v)\right)\right)_i\quad i=1,2,\ldots, n,
%\]
%that is, if $v_i=1$, then we should have
%\begin{equation}
%\left( \b L(\b v)\right)_i \ge 1,
%\label{eq:eqpforvone}
%\end{equation}
%and for $v_1=-1$, we should have
%\begin{equation}
%\left( \b L(\b v)\right)_i \le -1.
%\label{eq:eqpforvmin}
%\end{equation}
%Combining~(\ref{eq:eqpforvone}) and (\ref{eq:eqpforvmin}), we obtain the following necessary and sufficient conditions for the vertex $\b v$ to be an equilibrium point.
%\begin{thm}
%\label{thm:nascforeq}
%A vertex $\b v$ of the hypercube $H_n$ is an equilibrium point of the gBSB neural net if and only if
%\[
%(\b L(\b v))_iv_i \ge 1, \quad i=1,2, \ldots ,n.
%\]
%\end{thm}
%Equivalently, a vertex $\b v$ of the hypercube $H_n$ is not an equilibrium point of the gBSB neural net if and only if
%\[
%(\b L(\b v))_iv_i < 1, \quad \mbox{ for some } \quad i=1,2, \ldots ,n.
%\]
%\begin{example}\rm
%We verify that the vertex $\b v^{(1)}=\left[\begin{array}{cc}
%1 & 1
%\end{array}\right]^{\top}$ of $H_2$ is not an equilibrium point of the BSB net with the weight matrix given by (\ref{eq:bsbweighttwo}) and the step size $\alpha=0.1$. Note that for the BSB net,
%\[
%\b L(\b x)=(\b I_n+\alpha \b W)\b x.
%\]
%We first compute $\b L\left(\b v^{(1)} \right)$ to get
%\[
%\b L\left( \b v^{(1)}\right) = \left(\b I_2+\alpha\b W\right)\b v^{(1)}
%= \left[\begin{array}{c}
%1.2\\
%0.9
%\end{array}\right].
%\]
%Note that because
%\[
%\left(\b L\left(\b v^{(1)}\right)\right)_2v_2^{(1)}=0.9 < 1,
%\]
%the vertex $\b v^{(1)}$ is not an equilibrium point of the above BSB net.
%
%On the other hand, the vertex $\b v^{(2)}=\left[\begin{array}{cc}
%1 & -1
%\end{array}\right]^{\top}$ is an equilibrium point of this BSB model because
%\[
%\left(\b L\left(\b v^{(2)}\right)\right)_1v_1^{(2)}=1.8\ge 1\quad\mb{and}\quad \left(\b L\left(\b v^{(2)}\right)\right)_2v_2^{(2)}=1.5\ge 1,
%\]
%that is, the necessary and sufficient conditions of Theorem~\ref{thm:nascforeq} are satisfied.
%\end{example}\markex
%Using the fact that $v_i^2=1$ and $\alpha >0$, we reformulate Theorem~\ref{thm:nascforeq} as follows.
%\begin{thm}
%\label{thm:refonascfeq}
%A vertex $\b v$ of the hypercube $H_n$ is an equilibrium point of the gBSB neural model if and only if
%\[
%\left(\b W\b v+\b b\right)_iv_i \ge 0, \quad i=1,2, \ldots ,n.
%\]
%\end{thm}
%
%
%In our analysis we use stability definitions in the sense of Lyapunov, which we now recall.
%\begin{defn}
%An equilibrium point $\b x^*$ of $\b x(k+1)=\b T(\b x(k))$ is (uniformly) stable if for every $\varepsilon > 0$ there is a $\delta=\delta(\varepsilon) > 0$ such that if $\|\b x(0)-\b x^*\|< \delta $ then  $\|\b x(k)-\b x^*\|< \varepsilon $ for all $k \ge 0$,
%where $\| \cdot\|$ may be any $p$-norm of a vector, for example the Euclidean norm $\|\cdot\|_2$.
%\end{defn}\markdefn
%
%The above definition says that if an equilibrium is stable then a trajectory starting sufficiently close to the equilibrium, will stay close to the equilibrium for all future times. If the equilibrium is not stable, then it is called unstable.
%\begin{defn}
%The equilibrium solution $\b x^*$ is asymptotically stable if it is stable and if there exists an $\eta > 0$ such that if $\left\Arrowvert\b x(0)-\b x^* \right\Arrowvert < \eta$ then $\lim_{k \rightarrow \infty}\b x(k) = \b x^*$.
%\end{defn}\markdefn
%
%The following theorem gives us a sufficient condition for a vertex of $H_n$ to be an asymptotically stable equilibrium point.
%\begin{thm}
%\label{thm_asymstab}
%Let $\b v$ be a vertex of the hypercube $H_n$. If
%\[
%(\b L(\b v))_i v_i > 1, \quad i=1,2, \ldots ,n,
%\]
%then $\b v$ is an asymptotically stable equilibrium point of the gBSB neural model.
%\end{thm}
%{\bf Proof}\quad Because $(\b L(\b x))_i x_i$ is a continuous function, there exists a neighborhood $N(\b v)$ about $\b v$ such that
%\[
%(\b L(\b x))_i x_i > 1\quad\mb{for all}\quad i=1,2,\ldots, n\quad\mb{and}\quad \b x\in N(\b v).
%\]
%Therefore, for any $\b x\in N(\b v)$, we have
%\[
%\b g(\b L(\b x))=\b v,
%\]
%which means that $\b v$ is indeed asymptotically stable equilibrium point.\\
%\hspace*{5.5in} $\Box$
%
%Taking into account the fact that $v_i^2=1$ and $\alpha>0$, we reformulate the above theorem as follows.
%\begin{thm}
%\label{thm:refofasst}
%Let $\b v$ be a vertex of the hypercube $H_n$. If
%\[
%\left( \b W\b v+\b b\right)_i v_i > 0, \quad i=1,2, \ldots ,n,
%\]
%then $\b v$ is an asymptotically stable equilibrium point of the gBSB neural model.
%\end{thm}
%
%\subsection{Synthesis of Neural Associative Memory}
%
%
%In our further discussion, we use a special class of matrices that are characterized in the following definition.
%\begin{defn}
%A matrix $\mbox{\boldmath$ D$} \in \mb{\Bb R}^{n \times n}$ is said to be row diagonal dominant if
%\[
%d_{ii} \ge \sum_{{j=1}\atop{j \ne i}}^{n}|d_{ij}|, \quad i=1,2,\ldots ,n.
%\]
%It is strongly row diagonal dominant if
%\[
%d_{ii} > \sum_{{j=1}\atop{j \ne i}}^{n}|d_{ij}|, \quad i=1,2,\ldots ,n.
%\]
%\end{defn}\markdefn
%
%Now, suppose we are given a matrix $\mbox{\boldmath$ V$} \in \mb{\Bb R}^{n \times r}$. Its pseudo-inverse, or the Moore-Penrose inverse, is defined as follows.
%\begin{defn}
%Given $\b V \in \mb{\Bb R}^{n \times r}$, a matrix $\b V^{\dagger} \in \mb{\Bb R}^{r \times n}$ is called a pseudo-inverse of $\mbox{\boldmath$ V$}$ if
%\begin{enumerate}
%\item $\b V\b V^{\dagger}\b V=\b V$;
%\item $\b V^{\dagger}\b V \b V^{\dagger}=\b V^{\dagger} $;
%\item $(\b V\b V^{\dagger})^{\top}=\b V\b V^{\dagger}$;
%\item $(\b V^{\dagger}\b V)^{\top}=\b V^{\dagger}\b V$.
%\end{enumerate}
%\end{defn}\markdefn
%
%The pseudo-inverse has the following useful properties:
%\begin{itemize}
%\item $(\b V^{\top})^{\dagger} =(\b V^{\dagger})^{\top}$;
%\item $(\b V^{\dagger})^{\dagger}=\b V$.
%\end{itemize}
%%For further properties of pseudo-inverses, we refer the reader to Campbell %and \newline Meyer~(1991) and Boullion and Odell~(1971).
%
%Let $\b v^{(1)},\ldots ,\b v^{(r)}$ be given vertices of $H_n$ that we wish to store as asymptotically stable equilibrium points of the gBSB neural net, that is, we are interested in building a neural autoassociative memory. We now present a number of technical results that will lead to formulation of a design algorithm for the weight matrix of the gBSB neural associative memory. Let $\b D\in \mb{\Bb R}^{n\times n}$ satisfy the condition
%\begin{equation}
%d_{ii} > \sum_{{k=1}\atop{k \ne i}}^{n}|d_{ik}|+|b_i|\qquad i=1,2,\ldots ,n.
%\label{eq:condonmd}
%\end{equation}
%We now present an expression for constructing the weight matrix $\b W$ ensuring that the given patterns are stored as asymptotically stable equilibrium points of the gBSB neural net.
%\begin{thm}
%\label{thm:simexfow}
%Let $\b V=\left[\begin{array}{ccc}
%\b v^{(1)} & \cdots & \b v^{(r)}
%\end{array}\right] \in \{ -1,1 \}^{n \times r}$. If
%\[
%\b W=\b D\b V\b V^{\dagger},
%\]
%where the matrix $\b D$ satisfies~(\ref{eq:condonmd}), then $\b v^{(1)},\ldots,\b v^{(r)}$ are asymptotically stable points of the gBSB net.
%\end{thm}
%{\bf Proof}\quad Using the first property of the Moore-Penrose inverse of $\b V$, we obtain
%\[
%\b W\b V=\b D\b V\b V^{\dagger}\b V=\b D\b V.
%\]
%Hence, for $j=1,\ldots,r$,
%\begin{equation}
%\b W\b v^{(j)}=\b D\b v^{(j)}.
%\label{eq:prowvj}
%\end{equation}
%Substituting the above into the condition of Theorem~\ref{thm:refofasst} and using~(\ref{eq:condonmd}) gives
%\begin{eqnarray*}
%\left(\b W\b v^{(j)}+\b b \right)_iv_i^{(j)} &=& \left(\b D\b v^{(j)}+\b b \right)_iv_i^{(j)}\\
%&=& \left( d_{ii}v_i^{(j)} + \sum_{{k=1}\atop{k \ne i}}^{n}d_{ik}v_k^{(j)}
%+b_i\right)v_i^{(j)}\\
%&\ge & d_{ii}-\sum_{{k=1}\atop{k \ne i}}^{n}|d_{ik}|-|b_i|\\
%&>&0.
%\end{eqnarray*}
%By Theorem~\ref{thm:refofasst}, $\b v^{(j)}$ is asymptotically stable. The above argument holds for all $j=1,\ldots,r$. The proof is complete.\\
%\hspace*{5.5in} $\Box$
%
%Note that the weight matrix $\b W=\b D\b V\b V^{\dagger}$ does not have to be symmetric thanks to presence of the matrix $\b D$ in the formula for $\b W$. The matrix $\b D$ is a design parameter and therefore it can be chosen to be nonsymmetric so that $\b W$ becomes nonsymmetric.
%\begin{corollary}
%\label{cor:verinimaofv}
%Every vertex that is contained in $\mb{Range}\,(\b V)$ is also stored as an asymptotically stable equilibrium point.
%\end{corollary}
%{\bf Proof}\quad Let $\b v\in\{-1,1\}^n$ be such that
%\[
%\b v\in\mb{Range}\,(\b V).
%\]
%Hence, there is a vector $\b{\gamma}\in \mb{\Bb R}^r$ such that
%\[
%\b v=\b V\b{\gamma}.
%\]
%Using the above and the first property of the pseudo-inverse matrix of $\b V$ yields
%\begin{eqnarray}
%\b W\b v &=& \b D\b V\b V^{\dagger}\b v\nonumber\\
%&=& \b D\b V\b V^{\dagger}\b V\b{\gamma}\nonumber\\
%&=& \b D\b V\b{\gamma}\nonumber\\
%&=& \b D\b v.
%\label{eq:wvisdv}
%\end{eqnarray}
%Substituting the above into the condition of Theorem~\ref{thm:refofasst} and performing manipulation similar to the ones in the proof of Theorem~\ref{thm:simexfow}, we arrive at the result.\\
%\hspace*{5.5in} $\Box$
%\begin{thm}
%Under the conditions of Theorem~\ref{thm:simexfow}, every vertex $\b v$ of $H_n$ that is contained in $\mb{Range}\,(\b V)^{\perp}$ and such that for some $i=1,\ldots,n$
%\[
%b_iv_i <0
%\]
%is not an equilibrium point of the gBSB net.
%\label{thm:notstore}
%\end{thm}
%{\bf Proof}\quad It follows from the definition of the pseudo-inverse that if $\b V^{\dagger}$ is a pseudo-inverse matrix of $\b V\in\mb{\Bb R}^{n\times r}$, then there exists a square matrix $\b M$ such that
%\begin{equation}
%\b V^{\dagger}=\b M\b V^{\top}.
%\label{eq:psismattr}
%\end{equation}
%Suppose that $\b v\in \mb{Range}\,(\b V)^{\perp}$. Then,
%\[
%\b V^{\top} \b v=\b 0,
%\]
%and therefore, by~(\ref{eq:psismattr}),
%\begin{equation}
%\b V^{\dagger}\b v=\b 0.
%\label{eq:vdagvlze}
%\end{equation}
%Hence,
%\[
%\b W\b v=\b D\b V\b V^{\dagger}\b v=\b 0,
%\]
%and
%\begin{equation}
%\left(\b W\b v+\b b \right)_iv_i=b_iv_i,\quad i=1,\ldots,n.
%\label{eq:wvbawvz}
%\end{equation}
%If for some $i$, we have $b_iv_i<0$, then
%\[
%\left(\b W\b v+\b b \right)_iv_i<0,
%\]
%which implies that $\b v$ is not an equilibrium point of the gBSB neural net.\\
%\hspace*{5.5in} $\Box$
%
%The above results are clear and easy to apply to construct the weight matrix $\b W$. However, if the diagonal elements of the matrix $\b D$ are too large, we may end up storing all of the vertices of the hypercube $H_n$. This follows from the following result concerning the weight matrix.
%\begin{thm}
%If the weight matrix $\b W$ of the gBSB neural net is such that for all $i=1,\ldots, n$,
%\[
%w_{ii} > \sum_{{k=1}\atop{k \ne i}}^{n}|w_{ik}|+|b_i|,
%\]
%then all the vertices of $H_n$ are asymptotically stable equilibrium points of the gBSB net.
%\label{thm:wstrdiado}
%\end{thm}
%{\bf Proof}\quad Let $\b v$ be an arbitrary vertex of $H_n$, that is, $\b v\in \{-1, 1\}^n$. Then, using the assumption on the weight matrix, we obtain
%\begin{eqnarray*}
%\left(\b W\b v+\b b\right)_iv_i &=& w_{ii} + \sum_{{k=1}\atop{k \ne i}}^{n}w_{ik}v_k v_i+b_i v_i\\
%&\ge & w_{ii} - \sum_{{k=1}\atop{k \ne i}}^{n}|w_{ik}|-|b_i|\\
%&>& 0.
%\end{eqnarray*}
%It follows from Theorem~\ref{thm:refofasst} that all vertices of $H_n$ are asymptotically stable points of the gBSB net.\\
%\hspace*{5.5in} $\Box$
%
%We now provide an ``improved'' method for constructing the weight matrix. The method is based on the following theorem.
%\begin{thm}
%\label{thm_lindep}
%Let $\b V=\left[\begin{array}{ccc}
%\b v^{(1)}&\cdots &\b v^{(r)}
%\end{array}\right] \in \{ -1,1 \} ^{n \times r}$. Let $\b B =\left[\begin{array}{ccc}
%\b b  & \cdots &\b b
%\end{array}\right] \in \mb{\Bb R}^{n \times r}$ be the matrix  consisting of the column vector $\b b$ repeated $r$ times, $\b \Lambda \in \mb{\Bb R}^{n \times n}$, ${\beta}_j=\sum_{i=1}^{r}\left(\b V^{\dagger}\b v^{(j)}\right)_{i}$, and $\beta=-\max _j|\beta _j|$, where $j=1,2,\ldots,r$. Suppose
%\[
%d_{ii} > \sum_{{k=1}\atop{k \ne i}}^{n}|d_{ik}|+|b_i(1- \beta)| \qquad i=1,2, \ldots ,n.
%\]
%If
%\begin{equation}
%\b W=(\b D\b V-\b B)\b V^{\dagger}+\b \Lambda\left(\b I_n-\b V\b V^{\dagger}\right),
%\label{eq:matwingbsb}
%\end{equation}
%then all the patterns $\b v^{(j)}$, $j=1,2, \ldots ,r$, are stored as asymptotically stable equilibrium points of the gBSB neural net.
%\end{thm}
%{\bf Proof}\quad Using the first property of the Moore-Penrose inverse of $\b V$, we obtain
%\begin{eqnarray*}
%\b W\b V&=& (\b D\b V-\b B)\b V^{\dagger}\b V+\b \Lambda\left(\b I_n-\b V\b V^{\dagger}\right)\b V\\
%&=&\b D\b V-\b B\b V^{\dagger}\b V\\
%&=& \b D\b V-\left[\begin{array}{cccc}
%\b B\b V^{\dagger}\b v^{(1)} & \b B\b V^{\dagger}\b v^{(2)}& \cdots & \b B\b V^{\dagger}\b v^{(r)}
%\end{array}\right]\\
%&=& \b D\b V-\left[\begin{array}{cccc}
%\b b\beta_1 & \b b\beta_2 & \cdots & \b b\beta_r
%\end{array}\right].
%\end{eqnarray*}
%The $j$-th column of the above is
%\[
%\b W\b v^{(j)}=\b D\b v^{(j)}-\b b\beta _j \qquad j=1,2, \ldots ,r.
%\]
%We now show that a sufficient condition for a pattern $\b v^{(j)}$ to be an asymptotically stable equilibrium of $\b x(k+1)=\b T(\b x(k))$ is satisfied. Performing simple manipulations and taking into account the assumptions we obtain for $j=1,2, \ldots ,r$:
%\begin{eqnarray*}
%\left(\b W\b v^{(j)}+\b b\right)_iv_i^{(j)}
%&=& \left(\b D\b v^{(j)}-\b b \beta_j+ \b b\right)_i v_i^{(j)}\\
%&=&d_{ii}+\sum_{{k=1}\atop{k \ne i}}^{n} d_{ik}v_k^{(j)}v_i^{(j)}+b_i(1-\beta_j)v_i^{(j)} \\
%&\ge& d_{ii}-\sum_{{k=1}\atop{k \ne i}}^{n}|d_{ik}|-\left\arrowvert b_i\left(1-\left(-\max_j|\beta_j|\right)\right)\right\arrowvert \\
%&= &d_{ii}-\sum_{{k=1}\atop{k \ne i}}^{n}|d_{ik}|-|b_i(1-\beta)|\\
%&>&0,
%\end{eqnarray*}
%where the last inequality follows from the assumption on the matrix $\b D$. Hence, by Theorem~\ref{thm:refofasst}, the given patterns are asymptotically stable equilibrium points of the gBSB model.\\
%\hspace*{5.5in} $\Box$
%
%
%If the patterns are linearly independent and $r<n$, then $\mbox{rank}\b V =r$ and therefore $\b V^{\dagger}\b V=\b I_r$. In this case $\beta _j=1$ for all $j$ and
%\[
%\b W\b v^{(j)}=\b D\b v^{(j)}-\b b\qquad j=1,2, \ldots ,r,
%\]
%Thus in this case, a sufficient condition for the patterns $\b v^{(j)}$, $j=1,2, \ldots ,r$, to be stored as asymptotically stable equilibrium points of the gBSB model is that for $i=1,2,\ldots ,n$,
%\[
%d_{ii} > \sum_{{k=1}\atop{k \ne i}}^{n}|d_{ik}|.
%\]
%The above condition indeed guarantees that all the desired patterns are stored as asymptotically stable equilibria. As we mentioned before, if $d_{ii}$'s are too large, then we may end up storing many other vertices of $H_n$ in addition to the given patterns that we want to be stored. To minimize the number of spurious patterns stored, we impose the constraints of the form
%\[
%d_{ii}<\sum_{{k=1}\atop{k \ne i}}^{n}|d_{ik}|+|b_i|,\qquad i=1,2,\ldots ,n,
%\]
%on the matrix $\b D$ while synthesizing the weight matrix $\b W$.
%
%We now discuss the role of different terms in the expression~(\ref{eq:matwingbsb}) on the weight matrix in storing the desired patterns and minimizing the number of spurious states. First, the term $(\b D\b V-\b B)\b V^{\dagger}$ is responsible for the given patterns to be stored as asymptotically stable equilibrium points of the gBSB model. This is because for a given pattern, say $\b v^{(j)}$, we have $\b \Lambda (\b I_n-\b V\b V^{\dagger})\b v^{(j)}=\b 0$, where $\b 0$ is a vector of zeros. The role of the term  $\b \Lambda (\b I_n-\b V\b V^{\dagger})$ is to destabilize patterns that are not in $\mb{Span} \left[\begin{array}{ccc}
%\b v^{(1)}& \cdots &\b v^{(r)}
%\end{array}\right]$. Specifically, if a pattern $\b y$ is orthogonal to $\mb{Span}\left[\begin{array}{ccc}
%\b v^{(1)}& \cdots & \b v^{(r)}
%\end{array}\right]$, then the term $\b \Lambda (\b I_n-\b V\b V^{\dagger})$ will not only destabilize the vertex of $H_n$ corresponding to the pattern $\b y$, if it is an equilibrium point of the gBSB model, it will make this point a non-equilibrium point of the gBSB model. To see this, consider the expression
%\[
%\b W\b y=(\b D\b V-\b B)\b V^{\dagger}\b y+\b \Lambda (\b I_n-\b V\b V^{\dagger})\b y.
%\]
%Because
%\[
%\b y \perp \mb{Span} \left[\begin{array}{ccc}
%\b v^{(1)}& \cdots &\b v^{(r)}
%\end{array}\right],
%\]
%we have
%\[
%\b W\b y = \b \Lambda (\b I_n-\b V\b V^{\dagger })\b y =\b \Lambda \b y.
%\]
%Therefore,
%\[
%\left(\b W\b y+\b b\right)_iy_i
%=\left(\b \Lambda \b y+\b b\right)_iy_i =\lambda_{ii}+\sum_{{j=1}\atop{j \ne i}}^{n}{\lambda }_{ij}y_jy_i+b_iy_i.
%\]
%Hence, if we choose $\b \Lambda$ such that for some $i=1,2,\ldots ,n$,
%\[
%{\lambda }_{ii}<-\sum_{{j=1}\atop{j \ne i}}^{n}|{\lambda }_{ij}|-|b_i|,
%\]
%then
%\[
%\left(\b W\b y+\b b\right)_iy_i <0 \quad\mb{for some}\quad i=1,2,\ldots ,n.
%\]
%Therefore, the vertex $\b y$ cannot be an equilibrium point. The inclusion of the matrix $\b B$ in the synthesis formula for the weight matrix $\b W$ was to prevent an automatic storage of the negatives of the given patterns as well as to somehow help to direct the gBSB model trajectories towards the given patterns. For this reason, we select $\b b$ to be in $\mb{Span}\{ \b V\}$. We recommend making $\b b$ a weighted sum of $\b v^{(i)}$, $i=1,2, \ldots ,r$, with positive weight coefficients. The above results can be used as guides for selecting the components of the matrices $\b D$, $\b \Lambda$, the vector $\b b$, and consequently the matrix $\b B$. We summarize the above considerations in the form of an algorithm for synthesizing the weight matrix $\b W$.
%\begin{center}
%{\bf Weight Matrix Construction Algorithm}
%\end{center}
%Given: $\b v^{(j)}$, $j=1,2, \ldots ,r$.
%\begin{itemize}
%\item Form $\b B=\left[\begin{array}{cccc}
%\b b & \b b & \cdots &\b b
%\end{array}\right]$, where
%\[
%\b b=\sum_{j=1}^{r}\epsilon _j\b v^{(j)},\quad \epsilon_j >0,\quad j=1,2,\ldots , r.
%\]
%\item Select $\b D$ such that
%\begin{eqnarray*}
%d_{ii} &>& \sum_{{k=1}\atop{k \ne i}}^{n}|d_{ik}|,\quad\mbox{and}\\
%d_{ii} &<& \sum_{{k=1}\atop{k \ne i}}^{n}|d_{ik}|+|b_i|, \qquad i=1,2,\ldots ,n.
%\end{eqnarray*}
%\item Choose $\b \Lambda$ such that for $i=1,2, \ldots, n$,
%\[
%{\lambda }_{ii}<-\sum_{{j=1}\atop{j \ne i}}^{n}|{\lambda }_{ij}|-|b_i|.
%\]
%\item Compute
%\[
%\b W=(\b D\b V-\b B)\b V^{\dagger}+\b \Lambda(\b I_n-\b V\b V^{\dagger}).
%\]
%\item Identify spurious equilibrium patterns. You can use the condition of Theorem~\ref{thm:refonascfeq} for this purpose.
%\end{itemize}
%\markex
%When the set of patterns to be stored is not linearly independent, we can still use the preceding algorithm to compute the weight matrix $\b W$. However, the second constraint on the matrix $\b D$ may lead to a violation of the condition on $\b D$ in Theorem~\ref{thm_lindep}. In such a case we cannot guarantee that all the given patterns will be stored, nor we will have any control on the number of spurious states. We may then relax this condition. If attempts in synthesizing associative memory for a set of linearly dependent patterns result in unacceptable designs, we may try a different approach. We can start by picking a maximal set of patterns that are linearly independent from the given set of patterns. Let us assume that there are $s$ patterns to be stored and $r$ of them form a maximal linearly independent subset. Suppose, for notational simplicity, that the patterns $\b v^{(q)}$, $q=1, \ldots ,r$, form a linearly independent set. We then apply the above algorithm to generate the weight matrix $\b W$ based on the linearly independent patterns $\b v^{(q)}$, $q=1, \ldots ,r$. After that, we modify the obtained weight matrix and the vector $\b b$ to ensure that the remaining patterns $\b v^{(j)}$, $j=r+1,r+2, \ldots ,s$, will also be stored. We will focus our attention on ensuring that the pattern $\b v^{(r+1)}$ is stored as an asymptotically stable point. The remaining patterns $\b v^{(j)}$, $j=r+2,r+3, \ldots ,s$, can be stored using the same procedure. Because
%\[
%\b v^{(r+1)}\in \mb{Span}\left\{\b v^{(1)},\b v^{(2)},\ldots ,\b v^{(r)}\right\},
%\]
%that is, $\b v^{(r+1)}$ is a linear combination of $\b v^{(i)}$, $i=1,2,\ldots, r$, it can be represented as
%\[
%\b v^{(r+1)}=\sum_{q=1}^{r}\gamma _q\b v^{(q)}=\b V \b \gamma,
%\]
%where $\b \gamma =\left[\begin{array}{ccc}
%\gamma_1  & \cdots & \gamma_r
%\end{array}\right]^{\top} \in \mb{\Bb R}^n$. Then, using~(\ref{eq:matwingbsb}), we obtain
%\begin{eqnarray*}
%\b W\b v^{(r+1)}&=&\b D\b V\b \gamma -\b B\b \gamma\\
%&=&\b D\b v^{(r+1)}-\left( \sum_{q=1}^{r}\gamma _q \right)\b b.
%\end{eqnarray*}
%Hence, if
%\begin{eqnarray*}
%\left(\b W\b v^{(r+1)}+\b b\right)v^{(r+1)}_i &=& \left(\b D\b v^{(r+1)}-\left[\sum_{q=1}^{r}\gamma_q-1\right]\b b\right)_iv^{(r+1)}_i\\
%&>& 0,\qquad i=1,2,\ldots ,n,
%\end{eqnarray*}
%then the vertex $\b v^{(r+1)}$ is asymptotically stable. Observe that because of the assumptions on the matrix $\b D$, we have
%\[
%\mb{sign}\left(\b D\b v^{(r+1)}\right)_i=\mb{sign}\left( v^{(r+1)}_i\right).
%\]
%With the above in mind, we obtain stronger sufficient conditions for $\b v^{(r+1)}$ to be asymptotically stable. They are:
%\[
%\mb{sign}\left[\left(1-\sum_{q=1}^{r}\gamma _q\right)\b b\right]_i=\mb{sign}\left(v^{(r+1)}_i\right) \qquad i=1,2,\ldots ,n.
%\]
%We can represent the above conditions as
%\[
%\mb{sign}\,b_i=\mb{sign}\left(v^{(r+1)}_i\right)\,\mb{sign}\left(1-\sum_{q=1}^{r}\gamma_q\right)\qquad i=1,2,\ldots ,n.
%\]
%Thus, if we adjust the components of $\b b$ so that the above conditions are satisfied, then we are guaranteed that $\b v^{(r+1)}$ is stored as an asymptotically stable equilibrium point. We may try to further adjust the components of $\b b$, if possible, to store the remaining patterns. Then, using the same matrices $\b D$, $\b \Lambda$, and $\b V$ used to compute the original weight matrix $\b W$, we can choose an appropriate $\mb{\bm$b$}$ to generate a new weight matrix $\b W$.
%
%If $\b b$ cannot be adjusted in the above way to accommodate all the patterns, we may then have to use the following procedure. First, recall the stability conditions that the patterns $\b v^{(r+1)},\ldots ,\b v^{(s)}$ should satisfy,
%\[
%\left(\b D\b v^{(k)}-\left[\sum_{q=1}^{r}\gamma_q-1\right]\b b\right)_i v^{(k)}_i>0,\qquad i=1,2,\ldots ,n, \quad k=r+1,r+2,\ldots ,s.
%\]
%Using the requirement for the matrix $\mb{\bm$D$}$ to be diagonally dominant, we arrive at stronger conditions of the form
%\[
%d_{ii}-\sum_{{j=1}\atop{j \ne i}}^{n}|d_{ij}|+\left(1-\sum_{q=1}^{r}\gamma _{q}\right)b_i v^{(k)}_i > 0,\quad i=1,\ldots ,n,\quad
%k=r+1,\ldots ,s.
%\]
%Thus, we may use the above conditions and try to decrease the magnitude of the $b_i$'s until the conditions are met for all $k=r+1,\ldots ,s$. However, we are not guaranteed to always satisfy the above conditions. If the conditions are satisfied, then employing the same matrices $\b D$, $\mb{\bm$\Lambda$}$, and $\b V$ used to compute the original weight matrix $\b W$, we choose an appropriate $\b b$ to generate a new weight matrix $\b W$.
%\begin{example}
%\label{ex_synthes}
%\rm
%We synthesize a gBSB based associative memory to store the following three patterns:
%\[
%\mb{\bm$v$}^{(1)}=\left[\begin{array}{c}
%-1\\
%-1\\
%1\\
%1\\
%1
%\end{array}\right],\qquad \mb{\bm$v$}^{(2)}=\left[\begin{array}{c}
%1\\
%1\\
%-1\\
%-1\\
%-1
%\end{array}\right],\quad \mb{ and }\quad \mb{\bm$v$}^{(3)}=\left[\begin{array}{c}
%1\\
%-1\\
%1\\
%-1\\
%1
%\end{array}\right].
%\]
%We first choose matrices $\mb{\bm$D$}$ and $\mb{\bm$\Lambda$}$. Let
%\[
%\mb{\bm$D$}=\left[\begin{array}{ccccc}
%1 & 0 & 0 & 0 & -0.5\\
%0 & 1 & 0 & 0 & 0\\
%0 & -0.2 & 1 & 0 & 0\\
%0 & 0 & 0 & 1 & 0\\
%0 & 0 & 0 & 0 & 1
%\end{array}\right],\qquad \mb{\bm$\Lambda$}=-2\mb{\bm$I$}_5.
%\]
%Set $\mb{\bm$V$}=\left[\begin{array}{ccc}\mb{\bm$v$}^{(1)} & \mb{\bm$v$}^{(2)} & \mb{\bm$v$}^{(3)} \end{array}\right]$ and $\mb{\bm$B$}=\left[\begin{array}{ccc}\mb{\bm$b$} & \mb{\bm$b$} & \mb{\bm$b$} \end{array}\right]$, where
%\[
%\mb{\bm$b$}=\left[\begin{array}{ccccc} 0.6 & -0.5 & 0.7 & -0.4 & 0.8 \end{array}\right]^{\top}.
%\]
%Note that $\mb{rank}\mb{\bm$V$}=2$. Let the step size $\alpha =0.9$. Using the above, we compute the weight matrix $\mb{\bm$W$}$ to get
%\[
%\mb{\bm$W$}=\left[\begin{array}{ccccc}
%-0.6500 & 0.2667 & -0.2667 & -1.3500 & -0.2667\\
%0.1250 & -1.0833 & -0.9167 & -0.1250 & -0.9167\\
%-0.1750 & -0.9500 & -1.0500 & 0.1750 & 0.9500\\
%-1.4000 & -0.0667 & 0.0667 & -0.6000 & 0.0667\\
%-0.2000 & -0.8667 & 0.8667 & 0.2000 & -1.1333
%\end{array}\right].
%\]
%We then check which vertices of $H_n$ satisfy the sufficiency condition of Theorem~\ref{thm_asymstab} for asymptotic stability. In other words, we wish to identify the spurious states. In this example only the given patterns satisfy the condition of Theorem~\ref{thm_asymstab}.
%\end{example}\markex
%\begin{example}\label{ex:gbsbsnaps}\rm In this example, we use a gBSB neural net with the state in the hypercube $\{-1, 1\}^{25}$. We construct the weight matrix to store the four letters of alphabet whose bit maps are shown in Figure~\ref{fig:bitmaps}.
%\begin{figure}
%\centerline{\epsfig{figure=sn_stor.eps,height=1in}}
%\caption{Bit maps of the characters $A$, $B$, $C$, and $D$ stored in the gBSB neural net of Example~\ref{ex:gbsbsnaps}.}
%\label{fig:bitmaps}
%\end{figure}
%The bit maps are $5\times 5$ with elements from the set $\{-1,1\}$ with $-1$ corresponding to black and 1 to white. The pattern vectors are obtained from the bit maps by stacking their columns. For example, if $A=\left[\begin{array}{ccccc}
%\b a_1,&\b a_2,&\b a_3,&\b a_4,&\b a_5
%\end{array}\right]$, where $\b a_i$, $i=1,2,\ldots,5$ is the $i$-th column of $A$, then the corresponding pattern vector is obtained from $A$ by applying the stacking operator $\b s:\mb{\Bb R}^{5\times 5} \to \mb{\Bb R}^{25}$  to get
%\[
%\b s(A)=A(:)=\left[\begin{array}{c}
%\b a_1\\
%\vdots\\
%\b a_5
%\end{array}\right] \in \{-1,1\}^{25}.
%\]
%The weight matrix, $\b W$, of the gBSB net was constructed using the weight matrix construction algorithm, where
%\[
%\b V=\left[\begin{array}{cccc}
%A(:) & B(:) & C(:) &D(:)
%\end{array}\right].
%\]
%The values of other parameters chosen in this example were: $\b D=1.1\b I_{25}$, $\b{\Lambda}=-2 \b I_{25}$, $\b b=0.3 A(:)+0.1 B(:)$, and $\alpha=0.5$. Snapshots of the network iterates for different initial patterns are shown in Figures~\ref{fig:snap_at}--\ref{fig:snap_dt}.
%\begin{figure}
%\centerline{\epsfig{figure=snap_at.eps,height=2.5in}}
%\caption{Snapshots of the iterations of the gBSB net of Example~\ref{ex:gbsbsnaps} converging towards the letter $A$.}
%\label{fig:snap_at}
%\end{figure}
%\begin{figure}
%\centerline{\epsfig{figure=snap_bt.eps,height=2.5in}}
%\caption{Snapshots of the iterations of the gBSB net of Example~\ref{ex:gbsbsnaps} converging towards the letter $B$.}
%\label{fig:snap_bt}
%\end{figure}
%\begin{figure}
%\centerline{\epsfig{figure=snap_ct.eps,height=2.5in}}
%\caption{Snapshots of the iterations of the gBSB net of Example~\ref{ex:gbsbsnaps} converging towards the letter $C$.}
%\label{fig:snap_ct}
%\end{figure}
%\begin{figure}
%\centerline{\epsfig{figure=snap_dt.eps,height=2.5in}}
%\caption{Snapshots of the iterations of the gBSB net of Example~\ref{ex:gbsbsnaps} converging towards the letter $D$.}
%\label{fig:snap_dt}
%\end{figure}
%
%\end{example}\markex
%\subsection{Learning}
%A pattern is learned by a memory if its noisy or incomplete version presented to the memory is mapped back to this pattern. In our discussion of the proposed iterative learning, we use the following theorem which allows one to iteratively compute the pseudo-inverse of a given matrix. The theorem can be found in Campbell and Meyer~\cite[page~54]{campbell91}.
%\begin{thm}
%\label{thm_pseuditer}
%Let $\mbox{\boldmath$ U$}\in \mb{\Bb R}^{n \times (r+1)}$ be partitioned by slicing off its last column, so that
%\[
%\b U=\left[\begin{array}{ccc}
%\b V & \vdots &\b y
%\end{array}\right].
%\]
%Let
%\begin{eqnarray*}
%\mbox{\boldmath$ k$}&=&\mbox{\boldmath$ V$}^{\dagger}\mbox{\boldmath$y$},\\
%\mbox{\boldmath$ u$}&=&\left(\mbox{\boldmath$ I$}_n-\mbox{\boldmath$ VV$}^{\dagger}\right)\mbox{\boldmath$ y$},\\
%\mbox{\boldmath$ z$}&=&\left\{ \begin{array}{ccc}
%\mbox{\boldmath$ u$}^{\dagger} & \mbox{if}& \mbox{\boldmath$ u \ne 0$}\\
%{\textstyle \frac{\mbox{\boldmath$ k$}^{\top}\mbox{\boldmath$ V$}^{\dagger}}{1+\mbox{\boldmath$ k$}^{\top}\mbox{\boldmath$ k$}}} & \mbox{if}& \mbox{\boldmath$ u=0$}.
%\end{array}
%\right.
%\end{eqnarray*}
%Then,
%\[
%\mbox{\boldmath$ U$}^{\dagger}=\left[ \begin{array}{c}
%\mbox{\boldmath$ V$}^{\dagger}-\mbox{\boldmath$ kz$}\\
%\mbox{\boldmath$ z$}
%\end{array}
%\right].
%\]
%\end{thm}
%
%Learning a new pattern requires updating the weight matrix $\mb{\bm$W$}$. This in turn involves augmenting the matrices $\mb{\bm$V$}$ and $\mb{\bm$B$}$, and recomputing $\mb{\bm$V$}^{\dagger}$. Augmenting  $\mb{\bm$V$}$ and $\mb{\bm$B$}$ is straightforward, and updating $\mb{\bm$V$}^{\dagger}$ can be performed recursively using the above theorem. The above procedure guarantees that the new pattern will be stored as an asymptotically stable equilibrium point if it is linearly independent of the learned patterns. If, however, the pattern to be learned is a linear combination of the stored patterns, then we do not update the matrix $\mb{\bm$V$}$ but rather modify the vector $\mb{\bm$b$}$ as described at the end of the previous subsection. After that we update the weight matrix $\mb{\bm$W$}$. We now illustrate the iterative learning process on a numerical example.
%\begin{example}
%\label{ex_learn}
%\rm
%Consider Example~\ref{ex_synthes}. Suppose that we wish to add
%\[
%\mb{\bm$v$}^{(4)}=\left[\begin{array}{ccccc}
%-1 &-1 & -1 &-1 & 1
%\end{array}\right]^{\top}
%\]
%to the previously designed memory. We augment the pattern matrix to include $\mb{\bm$v$}^{(4)}$ in it. Thus, we get $\mb{\bm$V$}=\left[\begin{array}{cccc}
%\mb{\bm$v$}^{(1)} & \mb{\bm$v$}^{(2)} & \mb{\bm$v$}^{(3)} & \mb{\bm$v$}^{(4)}
%\end{array}\right]$. We also augment the matrix $\mb{\bm$B$}$ by adding extra $\mb{\bm$b$}$ column to it. We then update $\mb{\bm$V$}^{\dagger}$ using Theorem~\ref{thm_pseuditer} and recompute the weight matrix $\mb{\bm$W$}$ to obtain
%\[
%\mb{\bm$W$}=\left[\begin{array}{ccccc}
%0.1714 & 0.8143 & 0.8286 & -0.5286 & -0.8143\\
%0.4643 & -0.8571 & -0.4643 & 0.2143 & -1.1429\\
%0.7786 & -0.3143 & 0.2214 & 1.1286 & 0.3143\\
%-0.8286 & 0.3143 & 0.8286 & -0.0286 & -0.3143\\
%-0.4857 & -1.0571 & 0.4857 & -0.0857 & -0.9429
%\end{array}\right].
%\]
%Next, we check which vertices of $H_n$ satisfy the condition of Theorem~\ref{thm_asymstab} for asymptotic stability. We find that only the given four patterns satisfy the condition of Theorem~\ref{thm_asymstab}.
%\end{example}\markex
%\subsection{Forgetting}
%A previously stored pattern is forgotten or deleted from the memory if a stimulus that is a perturbed version of that pattern is presented to the memory and is not mapped back to the stored pattern. We now discuss the case when we wish to update the weight matrix in such a way that the specified patterns are removed from memory. We refer to the process of removing patterns from the associative memory as forgetting. Our discussion uses the following theorem which can be considered as the inverse of Theorem~\ref{thm_pseuditer}. A general version of the theorem appears in Boullion and Odell~\cite[page~14]{boullion71}.
%\begin{thm}
%\label{thm_pseudinv}
%Let $\mbox{\boldmath$ U$}\in \mb{\Bb R}^{n \times (r+1)}$ be partitioned by slicing off its last column, so that
%\[
%\b U=\left[\begin{array}{ccc}
%\b V & \vdots & \b y
%\end{array}\right],
%\]
%and let its pseudoinverse $\mbox{\boldmath$U$}^{\dagger} \in \mb{\Bb R}^{(r+1) \times n}$ be partitioned by slicing off its last row, so that
%\[
%\mbox{\boldmath$U$}^{\dagger}=\left[\begin{array}{c}
%\mbox{\boldmath$S$}\\
%\mbox{\boldmath$z$}
%\end{array}\right].
%\]
%Then,
%\[
%\mbox{\boldmath$V$}^{\dagger}=\left\{\begin{array}{cc}
%\mbox{\boldmath$S$}+\mbox{\boldmath$S$}\frac{\mbox{\boldmath$yz$}}{1-\mbox{\boldmath$zy$}} & \mbox{if}\quad \mbox{\boldmath$zy$}\ne 1\\
%\mbox{\boldmath$S$}-\mbox{\boldmath$S$}\mbox{\boldmath$z$}^{\dagger}\mbox{\boldmath$z$} & \mbox{if}\quad \mbox{\boldmath$zy$}=1.
%\end{array}\right.
%\]
%\end{thm}
%{\bf Proof}\quad Referring to Theorem~\ref{thm_pseuditer}, we first show that $\mbox{\boldmath$u$}=\mbox{\boldmath$0$}$ if and only if $\mbox{\boldmath$zy$}\ne 1$, or equivalently, $\mbox{\boldmath$u$}\ne \mbox{\boldmath$0$}$ if and only if $\mbox{\boldmath$zy$}= 1$. Indeed, if $\mbox{\boldmath$u$}=\mbox{\boldmath$0$}$, then
%\[
%\mbox{\boldmath$z$}=\frac{\mbox{\boldmath$k$}^{\top}\mbox{\boldmath$V$}^{\dagger}}{1+\mbox{\boldmath$ k$}^{\top}\mbox{\boldmath$ k$}}.
%\]
%Hence,
%\[
%\mbox{\boldmath$zy$}=\frac{\mbox{\boldmath$k$}^{\top}\mbox{\boldmath$V$}^{\dagger}\mbox{\boldmath$y$}}{1+\mbox{\boldmath$ k$}^{\top}\mbox{\boldmath$ k$}}=\frac{\mbox{\boldmath$ k$}^{\top}\mbox{\boldmath$ k$}}{1+\mbox{\boldmath$ k$}^{\top}\mbox{\boldmath$ k$}}<1.
%\]
%On the other hand, if $\mbox{\boldmath$u$}\ne\mbox{\boldmath$0$}$, then
%\[
%\mbox{\boldmath$z$}=\mbox{\boldmath$u$}^{\dagger}=\left\{(\mbox{\boldmath$I$}_n-\mbox{\boldmath$VV$}^{\dagger})\mbox{\boldmath$y$}\right\}^{\dagger}.
%\]
%Taking into account the fact that $\mbox{\boldmath$u$}\ne \mbox{\boldmath$0$}$ and that $\mbox{\boldmath$u$}\in \mb{\Bb R}^n$, its pseudoinverse is
%\[
%\mbox{\boldmath$u$}^{\dagger}=(\mbox{\boldmath$u$}^{\top}\mbox{\boldmath$u$})^{-1}\mbox{\boldmath$u$}^{\top}=\left(\mbox{\boldmath$y$}^{\top}\left(\mbox{\boldmath$I$}_n-\mbox{\boldmath$VV$}^{\dagger}\right)\mbox{\boldmath$y$}\right)^{-1}\mbox{\boldmath$y$}^{\top}\left(\mbox{\boldmath$I$}_n-\mbox{\boldmath$VV$}^{\dagger}\right),
%\]
%because $\left(\mbox{\boldmath$I$}_n-\mbox{\boldmath$VV$}^{\dagger}\right)$ is an orthogonal projector and hence,
%\[
%\left(\mbox{\boldmath$I$}_n-\mbox{\boldmath$VV$}^{\dagger}\right)^2=\left(\mbox{\boldmath$I$}_n-\mbox{\boldmath$VV$}^{\dagger}\right).
%\]
%Therefore,
%\[
%\mbox{\boldmath$zy$}=\mbox{\boldmath$u$}^{\dagger}\mbox{\boldmath$y$}=1.
%\]
%We now prove the expressions for $\mbox{\boldmath$V$}^{\dagger}$. We start with the case when $\mbox{\boldmath$zy$}\ne 1$. It follows from Theorem~\ref{thm_pseuditer} that
%\begin{equation}
%\mbox{\boldmath$S$}=\mbox{\boldmath$V$}^{\dagger}-\mbox{\boldmath$V$}^{\dagger}\mbox{\boldmath$yz$}=\mbox{\boldmath$V$}^{\dagger}\left(\mbox{\boldmath$I$}_n-\mbox{\boldmath$yz$}\right).
%\label{eq:afapptfo}
%\end{equation}
%Because $\mbox{\boldmath$zy$}\ne 1$, the matrix $\left(\mbox{\boldmath$I$}_n-\mbox{\boldmath$yz$}\right)$ is invertible. Using the matrix inversion formula~(\ref{eq:specialsmw}), from the Appendix, gives
%\begin{equation}
%\left(\mbox{\boldmath$I$}_n-\mbox{\boldmath$yz$}\right)^{-1}=\mbox{\boldmath$I$}_n+\frac{\mbox{\boldmath$yz$}}{1-\b z \b y}.
%\label{eq:simpappofminvf}
%\end{equation}
%Applying formula~(\ref{eq:simpappofminvf}) to (\ref{eq:afapptfo}) yields
%\[
%\mbox{\boldmath$V$}^{\dagger}=\mbox{\boldmath$S$}+\mbox{\boldmath$S$}\frac{\mbox{\boldmath$yz$}}{1-\mbox{\boldmath$zy$}}.
%\]
%
%We now prove the expression for $\mbox{\boldmath$V$}^{\dagger}$ for the case when $\mbox{\boldmath$zy$}=1$. It follows from Theorem~\ref{thm_pseuditer} that
%\begin{equation}
%\mbox{\boldmath$S$}=\mbox{\boldmath$V$}^{\dagger}-\mbox{\boldmath$kz$},
%\label{eq:exfos}
%\end{equation}
%and hence, because $\b z\ne \b 0$,
%\begin{equation}
%\mbox{\boldmath$S$}\mbox{\boldmath$z$}^{\dagger}=\mbox{\boldmath$V$}^{\dagger}\mbox{\boldmath$z$}^{\dagger}-\mbox{\boldmath$k$}.
%\label{eq:szpsinv}
%\end{equation}
%But in this case, $\mbox{\boldmath$u$}\ne \mbox{\boldmath$0$}$, and therefore
%\[
%\mbox{\boldmath$z$}^{\dagger}=\left(\mbox{\boldmath$u$}^{\dagger}\right)^{\dagger}=\mbox{\boldmath$u$}=\left(\mbox{\boldmath$I$}_n-\mbox{\boldmath$VV$}^{\dagger}\right)\mbox{\boldmath$y$}.
%\]
%Substituting the above expression for $\mbox{\boldmath$z$}^{\dagger}$ into~(\ref{eq:szpsinv}), we get
%\begin{eqnarray*}
%\mbox{\boldmath$S$}\mbox{\boldmath$z$}^{\dagger}&=&\mbox{\boldmath$V$}^{\dagger}\left(\mbox{\boldmath$I$}_n-\mbox{\boldmath$VV$}^{\dagger}\right)\mbox{\boldmath$y$}-\mbox{\boldmath$k$}\\
%&=&\left(\mbox{\boldmath$V$}^{\dagger}-\mbox{\boldmath$V$}^{\dagger}\right)\mbox{\boldmath$y$}-\mbox{\boldmath$k$}\\
%&=&-\mbox{\boldmath$k$}.
%\end{eqnarray*}
%We now substitute the above into~(\ref{eq:exfos}) to obtain
%\[
%\mbox{\boldmath$V$}^{\dagger}=\mbox{\boldmath$S$}-\mbox{\boldmath$S$}\mbox{\boldmath$z$}^{\dagger}\mbox{\boldmath$z$},
%\]
%which completes the proof of the theorem.\\
%\hspace*{5.5in} $\Box$
%
%We can delete a pattern from the memory by first removing it from $\mb{\bm$V$}$, then updating $\mb{\bm$V$}^{\dagger}$ using Theorem~\ref{thm_pseudinv}, removing one column from $\mb{\bm$B$}$, and recomputing the weight matrix $\mb{\bm$W$}$. We illustrate the iterative forgetting process on the following numerical example.
%\begin{example}
%\label{ex_forget}
%\rm
%We consider  Example~\ref{ex_learn}. Suppose that we wish to delete
%\[
%\mb{\bm$v$}^{(3)}=\left[\begin{array}{ccccc}
%1 &-1 &1 &-1 & 1
%\end{array}\right]^{\top}
%\]
%from the previously designed memory. We can use Theorem~\ref{thm_pseudinv} to update $\mb{\bm$V$}^{\dagger}$, that is, to find the pseudoinverse of $\mb{\bm$V$}=\left[\begin{array}{ccc}
%\mb{\bm$v$}^{(1)} & \mb{\bm$v$}^{(2)} & \mb{\bm$v$}^{(4)}
%\end{array}\right]$. We also update the matrix $\mb{\bm$B$}$ by deleting from it one column. We leave the remaining data unchanged. We recompute the weight matrix $\mb{\bm$W$}$ to obtain
%\[
%\mb{\bm$W$}=\left[\begin{array}{ccccc}
%-0.7333 & 1.2667 & 0.1500 & 0.1500 & -1.2667\\
%0.9167 & -1.0833 & -0.1250 & -0.1250 & -0.9167\\
%0.0500 & 0.0500 & -0.3250 & 1.6750 & -0.0500\\
%-0.0667 & -0.0667 & 1.4000 & -0.6000 & 0.0667\\
%-0.8667 & -0.8667 & 0.2000 & 0.2000 & -1.1333
%\end{array}\right].
%\]
%We then check which vertices of $H_n$ satisfy the condition of Theorem~\ref{thm_asymstab} for asymptotic stability. We find that only $\mb{\bm$v$}^{(1)}$, $\mb{\bm$v$}^{(2)}$, and $\mb{\bm$v$}^{(4)}$  satisfy the condition of Theorem~\ref{thm_asymstab}, which is what we desire.
%\end{example}\markex

%\section{Principal Components Analysis}
\section{Notes}
This chapter gives an introduction to computing using artificial neural
 networks from a dynamical system perspective. The disclaimer that we found
  in the book by Hertz, Krogh, and Palmer~\cite[page~xix]{hertz91} also applies to this chapter as well as to the whole field of artificial neural networks, ``The adjective ``neural'' is used because much of the inspiration for such networks comes from neuroscience, {\em not\/} because we are concerned with networks of real neurons. Brain modeling is a different field and, though we sometimes describe biological analogies, our prime concern is with what {\em artificial\/} networks can do, and why.'' The readers interested in principles underlying the organization of biological neurons into circuits and networks within different regions of the brain may consult the book edited by Shepherd~\cite{shepherd90}. For textbooks with an engineering approach to neural networks, we recommend Hassoun~\cite{hassoun95}, Hagan, Demuth, and Beale~\cite{hagan96}, Zurada~\cite{zurada92}, Haykin~\cite{haykin99}, Fausett~\cite{fausett94}, Khanna~\cite{khanna90}, or Harvey~\cite{harvey94}. Applications of artificial neural networks to signal processing and solving optimization problems are presented in the book by Cichocki and Unbehauen~\cite{cichocki93}. A psychological slant to the subject is presented in the book by Rumelhart, McClelland, and the PDP Research Group~\cite{rumelhart86}, and Anderson~\cite{anderson95}. Impatient readers wishing to get a feel of artificial neural networks may consult a tutorial by Lippmann~\cite{lippmann87}, and also by Hush and Horne~\cite{hush93}.

According to Haykin~\cite[page~38]{haykin99}, ``The modern area of neural
 networks began with the pioneering work of McCulloch and Pitts (1943). McCulloch was a psychiatrist and neuroanatomist by training; he spent some 20 years thinking about the representation of an event in the nervous system. Pitts was a mathematical prodigy, who joined McCulloch in 1942.'' In their paper~\cite{mcculloch43}, McCulloch and Pitts, assuming that the activity of the neuron is an ``all-or-none'' process, developed a logical calculus of neural events. They argued that, with addition of memory, a neural network consisting of a sufficiently large number of their neurons can compute anything that can be computed on any digital computer. McCulloch and Pitts' paper influenced John von~Neumann who used idealized switch-delay elements derived from the McCulloch-Pitts neuron model in the construction of the EDVAC~(Electronic Discrete Variable Automatic Computer) that was developed out of the ENIAC~(Electronic Numerical Integrator and Computer). The ENIAC was the first general purpose computer and it was invented by J. Presper Eckert and John Mauchly.


More recent investigators of the computational processes in the nervous
 system compare the brain to a hybrid, analog-digital computer~\cite{koch97}. According to Koch~\cite[page~210]{koch97}, individual nerve cells communicate between each other via streams of digital pulses. Information is, on the other hand, processed in the analog domain. The resulting signal is then converted back into digital pulses and conveyed to the following neurons.

%For other analog learning algorithms for adaptive linear elements
%(adalines) in the spirit of Section~\ref{sec:idadaline}, see~\cite{sira00}.
%
%As we mentioned before, the Brain-State-in-a-Box~(BSB) neural model
% was proposed by Anderson et al.\ in 1977~\cite{anderson77}.
%  The dynamics of the BSB model and its variants were further analyzed
%   by Greenberg~\cite{greenberg88}, Golden~\cite{golden86,golden93},
%    Grossberg~\cite{grossberg88}, Hui and \.{Z}ak~\cite{hui92bsb},
%     Hui, Lillo, and \.{Z}ak~\cite{hui93b},
%      Schultz~\cite{schultz93}, Lillo et al.~\cite{lillo94},
%       Anderson~\cite{anderson95}, Perfetti~\cite{perfetti95},
%       Hassoun~\cite{hassoun95}, \.{Z}ak, Lillo, and Hui~\cite{zak96},
%        Chan and \.{Z}ak~\cite{chan97,chan98}, Park, Cho, and Park~\cite{parkcho99}.
%         A continuous counterpart of the BSB model referred to as the linear
%          system operating in a saturated mode was analyzed in~Li, Michel,
%          and Porod~\cite{jhli89}.
%
%
%One of the first techniques used to design neural network associative
% memories is the {\em outer product method}. This method was first
%  used by Hopfield in~\cite{hopfield82} and~\cite{hopfield84}. The
%  outer product method gives a network with a symmetric interconnection
%   structure that does not necessarily store the desired patterns as
%    equilibrium points of the network. The network designed using this
%     method can efficiently store up to $0.15n$ arbitrary patterns and
%     has learning capabilities. The {\em projection learning rule\/}
%     proposed by Personnaz, Guyon, and Dreyfus in~\cite{personnaz85}
%     and in~\cite{personnaz86} gives a network with a symmetric
%     interconnection structure that always stores the desired
%     patterns as equilibrium points although not necessarily asymptotically
%     stable. The network designed with this method can effectively store up to $0.5n$
%     desired patterns and has learning capabilities. The {\em eigenstructure method\/}
%     used by Li et al.~\cite{jhli89} yields a neural network associative memory capable
%      of storing any given pattern as an asymptotically stable equilibrium point.
%      The network designed with this method is capable of effectively storing a number
%       of patterns that may exceed the order of the network. However, the resulting
%       network will have a symmetric interconnection structure. Furthermore, there are
%        no provisions for learning for networks designed by the eigenstructure method.
%         In a series of papers, Michel, Farrell, and Porod~\cite{michel89}, Michel and
%          Farrell~\cite{michel90a}, Michel, Farrell, and Sun~\cite{michel90b}, and Yen
%and Michel~\cite{yen91} modified the eigenstructure method. The modified method can
%be used to design networks capable of storing on order of $0.5n$ patterns as
%asymptotically stable equilibrium points. Networks designed by this method need not
%have symmetric interconnecting structures. Furthermore, learning is incorporated
%into the method. Networks designed by the previous methods use energy function
%techniques. A benefit of using energy functions is that the resulting networks are
%globally stable, that is, the networks' trajectories always tend to some equilibrium
%points. On the other hand, if the modified method yields a nonsymmetric
%interconnecting structure, the resulting network is not guaranteed to be globally
%stable.

\section{Exercises}
\begin{exercise}\rm
Write a MATLAB program that implements the TLU~(Perceptron) Learning Algorithm for generating a separating line for points on a plane. Test your program on the following data. For $i=1,\ldots,4$ and $d_i=1$,
\[
\left[\begin{array}{cccc}
\b x^{(1)} & \b x^{(2)} & \b x^{(3)} & \b x^{(4)}
\end{array}\right]= \left[\begin{array}{cccc}
5 & -7 & 8 & -1\\
5  & 6 & -3 & 2
\end{array}\right].
\]
The above points constitute class $\mb{\Sc C}_1$. You can label the above points using, for example, ``x''. Class $\mb{\Sc C}_2$ consists of the points corresponding to $d_i=-1$. They are,
\[
\left[\begin{array}{ccc}
\b x^{(5)} & \b x^{(6)} & \b x^{(7)}
\end{array}\right]= \left[\begin{array}{ccc}
-6 & -4 & 7\\
1  & -4 & -5
\end{array}\right]
\]
You can label these points using ``o''. The initial weight vector and threshold are
\[
\b w(0)=\left[\begin{array}{cc}
-5 & -1
\end{array}\right]^{\top}\quad\mb{and}\quad t_0=5.
\]
Plot the initial separating line as well as the final separating line.
\end{exercise}

\begin{exercise}\label{ex:wtawthrcl}\rm
Write a MATLAB program that implements the WTA algorithm. Test your program on the following data. The initial weights of the WTA network are:
\[
\left[\begin{array}{ccc}
\b w_0^{(1)} & \b w_0^{(2)} & \b w_0^{(3)}
\end{array}\right]=\left[\begin{array}{ccc}
0.4403  &  0.7411 &  -0.8866\\
0.8979 &  -0.6714 &   0.4625
\end{array}\right],\quad \alpha=0.5.
\]
The training patterns are:
\[
\b X=\left[\begin{array}{ccccccccc}
-0.2&0.2&-.1&0.98&0.98&.8&-0.6&-0.8&-0.7\\
0.98&0.98&.9&0.2&-0.2& 0.1&-0.8&-0.6&-0.7
\end{array}\right]
\]
Present the training set 200 times. Find the weights of the network after training.
\end{exercise}
\begin{exercise}\label{exs:nfuzid}\rm
Use the identification scheme depicted in Figure~\ref{fig:nufu} and the neural fuzzy on-line identification algorithm from Section~\ref{subsec:nfuzid} to identify the plant described by
\[
y(k+1)=\frac{y(k)(y(k-1)+2)y(k)+2.5)}{8.5+y^2(k)+y^2(k-1)}+u(k),
\]
where
\[
u(k)=2\cos(2\pi k/100)\quad\mb{for}\quad k\le 200
\]
and
\[
u(k)=1.2\sin(2\pi/20)\quad\mb{for}\quad k>200.
\]
The identifier is to be described by
\[
\hat{y}(k+1)=f(k)+u(k),
\]
where $f$ has the form of~(\ref{eq:flcasnlsys}) with $M=100$. Choose $\alpha=0.5$, $y(0)=y(1)=0$, $f(1)=f(2)=0$, and start with $u(1)=0$. (The above plant was also used by Lin and Lee~\cite[page~731]{lin96} to test their identification scheme.)
\end{exercise}
\begin{exercise}\rm
Write a MATLAB program that implements the learning algorithm for the RBF identifier discussed in Section~\ref{sec:rbfnets}.
\end{exercise}
%\begin{exercise}\label{ex:parkcholem}\rm
%Prove the following lemma that can be found in~\cite{parkcho99}.
%\begin{lemma}
%\label{lem:parkchopark}
%Let $\b v\in \{-1, 1\}^n$ be an asymptotically stable equilibrium point of the gBSB neural net. Let $\b x$ be a vertex of $H_n$ that differs from $\b v$ in the component $x_i$ for some $i=1,\ldots, n$. If
%\[
%\left(\b W\b x+\b b\right)_iv_i >0,
%\]
%then $\b x$ is not an equilibrium point of the gBSB neural net.
%\end{lemma}
%\end{exercise}
%\begin{exercise}\label{ex:parkchoprop}\rm
%Use Lemma~\ref{lem:parkchopark} from Exercise~\ref{ex:parkcholem} to prove the following proposition that can be found in~\cite{parkcho99}.
%\begin{proposition}
%\label{prop:parkchopark}
%Let $\b v\in \{-1, 1\}^n$ be an asymptotically stable equilibrium point of the gBSB neural net. If $w_{ii}=0$ for $i=1,\ldots,n$, then none of the vertices of $H_n$ that differ from $\b v$ in just one component is an equilibrium point of the gBSB neural net.
%\end{proposition}
%\end{exercise}
%\begin{exercise}\label{ex:parkchoproptwo}\rm
%Use Lemma~\ref{lem:parkchopark} from Exercise~\ref{ex:parkcholem} to prove the following proposition that can be found in~\cite{parkcho99}. Before stating the proposition, we define the {\em Hamming distance} between two vectors whose components are $-1$ or 1 as the number of components in which the two vectors differ. Given two vectors $\b x$ and $\b y$, the Hamming distance between them is denoted $HD(\b x,\b y)$.
%\begin{proposition}\label{prop:parkchoparktwo}
%Let $\b v\in \{-1, 1\}^n$ be an asymptotically stable equilibrium point of the gBSB neural net and $k$ an integer from the set $\{1,2,\ldots,n\}$. If
%\[
%\left(\sum_{j=1}^n w_{ij}x_j+b_i \right)v_i > 2 k \max_j |w_{ij}|
%\]
%then a vertex $\b x$ of $H_n$ such that, $0<HD(\b x,\b v)\le k$, is not an equilibrium point of the gBSB neural net.
%\end{proposition}
%\end{exercise}
