\chapter{Neural Networks}
This chapter is concerned with neurocomputing which is also called brainlike computation. The motivation to study artificial neural networks, commonly referred to as neural networks, or just neural nets, is in attempt to build computers whose construction would be inspired by the brain. The brain is a highly complex information processing system. It is capable of organizing its building blocks, called neurons, to perform complex computations in a very efficient way. A drawing of a more or less typical biological neuron is given in Figure~\ref{fig:bio_neur}.
\begin{figure}
\vspace{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. The signals have the form of electrical pulses. These pulses travel to the neuron's cell body where they are being processed. We note that the surface of the cell body of the neuron is dotted with numerous receptor sites that receive messages form other neurons in a similar fashion as the dendrites do. 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, 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, whose drawing is shown in Figure~\ref{fig:neur_syn}.
\begin{figure}
\vspace{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. 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, which are delivered down to the synaptic vesicles. The neurotransmitters are stored in the synaptic vesicles until the neuron fires in which case 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. 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. In fact 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. These models are referred to as the artificial neurons. Then, we will analyze neural networks which are composed of the artificial neurons.

\section{Threshold Logic Unit}
We now discuss a simple model of a biological neuron called the {\em threshold logic unit\/}~(TLU). We analyze the TLU's computational capability. A TLU is parameterized by a set of weights and a threshold. A diagram of a TLU is shown in Figure~\ref{fig:tlu_ner}, where the summation operation followed by the sign activation function model the cell body of a biological neuron. The dendrites are modeled by the inputs and their strength are modeled by the weights. The axon corresponds to the output $y$.
\begin{figure}
\vspace{3in}
\hskip1in\special{bmp:c:/mybook/tlu_ner.bmp x=4in y=3in}
\caption{A threshold logic unit~(TLU) 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]^T \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 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 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}

In this section, we study properties of the TLU. In particular, we explore issues involving learning of the TLU. 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
\[
\langle \b w^*,\b x\rangle -t^*=0,
\]
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\/}, and the hyperplane is referred to as a {\em separating hyperplane\/}. If no such hyperplane exists, then the training set is referred to as {\em linearly nonseparable}. Finding a separating hyperplane is equivalent to finding a corresponding weight vector and a threshold. If we find such a weight vector and a threshold, then we say that the TLU correctly classifies the 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^*>0&\mb{if}& d_i=1.
\end{array}
\]
If a training set is linearly separable, then the following training algorithm will determine a separating hyperplane.
\begin{center}
{\bf The TLU Learning Algorithm}
\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( \langle \b w^{(k)},\b x^{(i)} \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 perceptron.
\item {\bf Repeat by Going to Step 2}
\end{enumerate}
\markex
\begin{example}\rm
We used the above algorithm to generate a separting 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]^T\quad\mb{and}\quad t_0=0.7621
\]
The 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]^T$ and threshold $t=-2.2379$.
\begin{figure}
\vspace{3.5in}
\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 separting line.}
\label{fig:lsep_ex}
\end{figure}
\end{example}
\markex
Our goal now is to prove that the 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$,
\[
\langle \b w,\b x^{(i)} \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( \langle \b w,\b x^{(i)} \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, $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.
\] 
Furthermore, we can 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,
\begin{eqnarray*}
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),
\end{eqnarray*}
where $\tilde{\b w}^{T}=\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 equivalently express the training problem of a TLU as follows:
\begin{center}
{\bf TLU Learning 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
We further simplify the formulation of the TLU Learning Problem by 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 their 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}
\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, 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\/} 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.
\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 since 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, 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)T} \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)^{T}\left( \b w^{(k-1)} + \b x^{(i)} \right)\\
&=& \b w^{(k-1)T}\b w^{(k-1)}+\b x^{(i)T} \b x^{(i)}+ 2\b w^{(k-1)T} \b x^{(i)}\\
&\le& \|\b w^{(k-1)}\|^2 +L^2,
\end{eqnarray*}
because $\b w^{(k-1)T}\b x^{(i)} \le 0$ and $\|\b x^{(i)}\|^2 \le L^2$. Suppose now 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 an upper bound on $\|\b w^{(k)}\|^2$. We have,
\begin{equation}
\|\b w^{(k)}\|^2 \le k L^2,
\label{eq:upbounonwk}
\end{equation}
because we assumed that $\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 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 obtain
\[
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 is complete.\\
\hspace*{5.5in} $\Box$

We now present a sufficient condition for linear nonseparability.
\begin{thm}
\label{thm:nonseparability}
If there exists a positive linear combination (PLC) 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$. Then, 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^{T}$ to get
\[
\b w^{T}\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 $\sum_{i=1}^p q_i\b x^{(j)} =\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 above training set represents the 2-input Exclusive-OR (XOR) function. 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 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.
\end{example}
\markex
\section{Backpropagation}
In the previous section, we discussed an iterative learning algorithm for a threshold logic unit. This computational unit kept on iterating until it reached a solution, if a solution existed. 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. Specifically, a typical backpropagation algorithm consists of two passes: the forward pass 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, 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.

We now analyze a single computational unit shown in Figure~\ref{fig:cudiffa}. It is called the {\em perceptron}.
\begin{figure}
\vspace{3in}
\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 wish to train the perceptron to learn the training pair
\[
\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^T\b w$. 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 to minimize $E$, we obtain
\begin{eqnarray}
\b w^{(k+1)}&=& \b w^{(k)}-\eta \nabla E\left(\b w^{(k)} \right)\nonumber\\
&=& \b w^{(k)}+\eta \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)T}\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 can represent~(\ref{eq:desgwfss}) as
\begin{equation}
\fbox{$\displaystyle\b w^{(k+1)}= \b w^{(k)}+\eta\delta^k\b x^{(k)}$}
\label{eq:sinpruad}
\end{equation}

We now generalize the above formula for a three-layer feedforward neural network shown in Figure~\ref{fig:fuzzybp}.
\begin{figure}
\vspace{3.5in}
\hskip1.2in\special{bmp:c:/mybook/fuzzybp.bmp x=4.67in y=3.5in}
\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 will derive rules for adjusting the weights $w_{ji}$ and $w_{kj}$ such that the error function
\[
E(\b w)=\frac{1}{2}\sum_{s=1}^m \left(d_s-y_s\right)^2
\]
is minimized. The minimization is performed over the space of weights
\[
\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. Let
\[
net_s=\sum_{j=1}^l w_{sj}^oz_j.
\]
Using the 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)}-\eta_o \frac{\partial E\left(\b w^{(k)} \right)}{\partial w_{sj}^{o(k)}}\nonumber\\
&=& w_{sj}^{o(k)}+\eta_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 in the hidden layer, we use the Chain Rule to compute $\partial E/\partial w_{ji}^h$, and then apply the gradient descent algorithm. We have
\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}
and
\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)}-\eta_h \frac{\partial E\left(\b w^{(k)} \right)}{\partial w_{ji}^{h(k)}}\nonumber\\
&=& w_{ji}^{h(k)}+\eta_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 is the sigmoid,
\begin{equation}
f(v)=\frac{1}{1+e^{-v}}.
\label{eq:sigmoid_act}
\end{equation}
It is easy to check the above activation function 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.
\begin{center}
{\bf The Backpropagation Training Algorithm}
\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)}-\eta\nabla_{\b w} E\left(\b w^{(k)},\b t^{(k)} \right),\\
\b t^{(k+1)} &=& \b t^{(k)}-\eta\nabla_{\b t} E\left(\b w^{(k)}, \b t^{(k)} \right).
\end{eqnarray*}
\item {\bf Repeat by Going to Step 2}
\end{enumerate}
\markex

\begin{example}\label{ex:backproxor}\rm
We consider the 2-input XOR problem. 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. In this example, we use the unipolar representation of the 2-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 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. 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}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}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 in a mixed order. We used a fixed 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.
\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}
\vspace{4in}
\hskip.5in\special{bmp:c:/mybook/xor_bp.bmp x=5.33 y=4in}
\caption{The function implemented by the trained neural network to solve the 2-input XOR problem in Example~\ref{ex:backproxor}.}
\label{fig:xor_bp}
\end{figure}
\end{example}
\markex


\section{Fuzzy Neural Networks}
In this section, we extend the back-propagation algorithm to the case of fuzzy data. Specifically, we present a learning algorithm, for feedforward neural networks, that operates on fuzzy sets rather than on crisp data. Thus the training patterns in the learning phase consist of pairs of fuzzy input vectors and fuzzy target outputs. In our presentation of the algorithm, we use some of the results that are reported in the paper by Ishibuchi, Fujioka, and Tanaka~\cite{ishibuchi93}.

\subsection{Training Pairs}
Suppose we are given $s$ input-output pairs
\[
\left(\mb{\bm$x$}^{(p)},d_p\right),\quad p=1,2,\ldots ,s,
\]
where,
\[
\mb{\bm$x$}^{(p)}=\left[x_1^{(p)},x_2^{(p)},\ldots,x_n^{(p)}\right]^T \in \mb{\Bb R}^n \quad \mb{and}\quad  d_p \in \mb{\Bb R}.
\]
We are also given $r-s$ fuzzy IF-THEN rules of the form:
\begin{equation}
\mb{IF}\; x_1\;\mb{is}\; A_{q1}\; \mb{AND}\; \cdots\; \mb{AND}\; x_n\; \mb{is}\; A_{qn}\; \mb{THEN}\; y\; \mb{is}\; B_q,
\label{eq:mnsfuzrul}
\end{equation}
where $q=1,2,\ldots,r-s$, and $A_{qi}$ and $B_q$ are linguistic values, that is, they are fuzzy sets. We assume that $A_{qi}$ and $B_q$ are fuzzy numbers. The crisp pairs $\left(\b x^{(p)},d_p\right)$ can be viewed as fuzzy singleton pairs, that is, for $p=1,2,\ldots,s$, $\b A_p=\left[A_{p1},A_{p2},\ldots,A_{pn}\right]=\left[x_{1}^{(p)},x_{2}^{(p)},\ldots,x_{n}^{(p)}\right]$, and $B_p=d_p$. Therefore, we can represent the numerical (crisp) data and the fuzzy IF-THEN rules in the same format as given by~(\ref{eq:mnsfuzrul}), where $q=1,2,\ldots,r$. In conclusion, we have $r$ fuzzy input-output pairs
\[
\left(\mb{\bm$A$}_q,B_q\right)\qquad q=1,2,\ldots,r.
\]
Our problem then is to train a feedforward neural network using fuzzy, rather than crisp, training pairs.

In the next subsection, we review some facts from the interval arithmetic that we will be using in the network training.

\subsection{Review of Some Facts From Interval Arithmetic}
The fact that we will be working with fuzzy numbers instead of crisp data requires a generalization of ordinary arithmetic on real numbers to arithmetic on closed intervals. This type of arithmetic is referred to as the interval arithmetic and it is discussed in some detail in Section~\ref{sec:fuzarith}. Here, we will summarize the basic operations of the interval arithmetic that we will employ in the network's training using IF-THEN rules. In particular, we will be using the operations of addition and multiplication on closed intervals. To this end suppose we have two closed intervals
\[
X=\left[x^L,x^U\right],\quad\mb{and}\quad Y=\left[y^L,y^U\right],
\]
where the superscripts $L$ and $U$ denote lower and upper limits, respectively. We define the operation of addition on the intervals $X$ and $Y$ as
\begin{eqnarray*}
X+Y&=&\left[x^L,x^U\right]+\left[y^L,y^U\right]\\
&=&\left[x^L+y^L,x^U+y^U\right].
\end{eqnarray*}
Let $k \in \mb{\Bb R}$. We define the product $k\cdot X$ to be
\begin{eqnarray*}
k\cdot X&=&k \cdot \left[x^L,x^U\right]\\
&=&\left\{\begin{array}{c}
\left[kx^L,kx^U\right]\quad \mb{if}\quad k\ge0\\
\left[kx^U,kx^L\right]\quad \mb{if}\quad k<0.
\end{array}
\right.
\end{eqnarray*}
We now extend the operation of the activation function from real numbers to closed intervals. The activation function that we will is the sigmoid,
\[
f(x)=\frac{1}{1+e^{-x}}.
\]
However, the argument of $f$ will now be a closed interval rather than a real number. Specifically, let $X=\left[x^L,x^U\right]$ be a closed interval. Then,
\begin{eqnarray*}
f(X)&=&f\left(\left[x^L,x^U\right]\right)\\
&=&\left[f\left(x^L\right),f\left(x^U\right)\right].
\end{eqnarray*}
The above operations on closed intervals will be used when manipulating fuzzy numbers. We will use the following operations of addition and multiplication on fuzzy numbers. Given two fuzzy numbers $A$ and $B$, and their corresponding membership functions $\mu_A$ and $\mu_B$. Then, the membership function of the fuzzy set $A+B$ is
\[
\mu_{A+B}(z)=\max\left\{\mu_A(x)\wedge\mu_B(y)\, : \,z=x+y \right\}.
\]
The operation of multiplication of a fuzzy number by a scalar $k \in \mb{\Bb R}$ is defined as:
\[
\mu_{kA}(z)=\max \left\{\mu_A(x)\, : \,z=kx \right\}.
\]
Finally, the operation of the activation function $f$ on a fuzzy number $A$ is defined as:
\[
\mu_{f\left(A\right)}(z)=\max\left\{\mu_A(x)\, : \,z=f(x)\right\}.
\]
The $h$-level set of a fuzzy number $A$, denoted $\left[A\right]_h$, is
\[
\left[A\right]_h=\left\{x\, : \,\mu_A(x) \ge h,\; 0<h\le 1\right\}.
\]
Note that the $h$-level set of a fuzzy number is a closed interval. Thus, we can operate on fuzzy numbers via their $h$-level sets. It is easy to verify the following relations between the operations on closed intervals and fuzzy numbers:
\begin{eqnarray*}
\left[A+B\right]_h&=&\left[A\right]_h+\left[B\right]_h,\\
\left[kA\right]_h&=&k\left[A\right]_h,\\
\left[f\left(A\right)\right]_h&=&f\left(\left[A\right]_h\right).
\end{eqnarray*}
We are now ready to discuss neural networks that learn from fuzzy IF-THEN rules.

\subsection{Neural Network Architecture}
Consider the three-layered feedforward neural network depicted in Figure~\ref{fig:fuzzybp}. The $n$ inputs to the network are now closed intervals that are $h$-level sets of given fuzzy numbers. We use uppercase letters to denote the inputs to the network. Thus, a closed interval corresponding to the $h$-level set of a fuzzy set $X_i$ is the $i$th input of the network. Hence, the input to the $j$-th neuron in the hidden layer is
\[
\left[V_j\right]_h=\sum_{i=1}^{n}w_{ji}^h \cdot \left[X_i\right]_h+\theta _i,
\]
where $w_{ji}^h$ are the weights of the hidden layer and the biases $\theta _j$ are real parameters that are updated in the training phase. We use the superscript $h$ in $w_{ji}^h$ in order to unambiguously designate the hidden layer weights. This should not be confused with the $h$-level set. The output of the $j$th neuron in the hidden layer, denoted $\left[Z_j\right]_h$, is
\begin{eqnarray*}
Z_j&=&f_j\left(\left[V_j\right]_h\right)\\
&=&f_j^h\left(\sum_{i=1}^{n}w_{ji}^h \cdot\left[X_i\right]_h + \theta _i \right).
\end{eqnarray*}
The output of the $s$-th neuron in the output layer, where $s=1,2,\ldots,m$, is
\begin{eqnarray*}
\left[Y_s\right]_h&=&f_s^o\left(\sum_{j=1}^{l}w_{sj}^{o}\cdot\left[Z_j\right]_h+\vartheta_j\right)\\
&=&f_s^o\left(\sum_{j=1}^{l}w_{sj}^o\cdot f_j^h\left(\sum_{i=1}^{n}w_{ji}^h \cdot \left[X_i\right]_h +\theta _i\right)+\vartheta_j \right)\\
&=&F_s\left(\left[X_1\right]_h,\ldots,\left[X_n\right]_h\right),
\end{eqnarray*}
where $w_{sj}^o$ are the weights of the hidden layer and the biases $\vartheta_j$ are real parameters updated in the training phase. Thus, the mapping that the neural network implements is
\[
\left[\begin{array}{c}
\left[Y_1\right]_h\\
\vdots\\
\left[Y_m\right]_h
\end{array}\right]=\left[\begin{array}{c}
F_1\left(\left[X_1\right]_h,\ldots,\left[X_n\right]_h\right)\\
\vdots\\
F_m\left(\left[X_1\right]_h,\ldots,\left[X_n\right]_h\right)
\end{array}\right].
\]
The above expressions are similar to those when the real numbers were used. But now we deal with fuzzy sets. In the training process, we will use interval vectors corresponding to specified level sets of fuzzy sets. Let
\begin{equation}
\mb{\bm$X$}_p=\left[\begin{array}{c}
X_{p1}\\
X_{p2}\\
\vdots\\
X_{pn}
\end{array}\right]=\left[\begin{array}{c}
\left[x_{p1}^L,x_{p1}^U\right]\\
\left[x_{p2}^L,x_{p2}^U\right]\\
\vdots\\
\left[x_{pn}^L,x_{pn}^U\right]
\end{array} \right]
\label{eq:bfcapxpinve}
\end{equation}
and
\begin{equation}
\mb{\bm$Y$}_{dp}=\left[\begin{array}{c}
Y_{dp1}\\
Y_{dp2}\\
\vdots\\
Y_{dpm}
\end{array}\right]=\left[\begin{array}{c}
\left[y_{dp1}^L,y_{dp1}^U\right]\\
\left[y_{dp2}^L,y_{dp2}^U\right]\\
\vdots\\
\left[y_{dpm}^L,y_{dpm}^U\right]
\end{array} \right]
\label{eq:desinterout}
\end{equation}
be interval vectors composed of the closed intervals corresponding to $h$-level sets of given fuzzy sets. Let
\[
\mb{\bm$Y$}_p=\left[\begin{array}{c}
Y_{p1}\\
Y_{p2}\\
\vdots\\
Y_{pn}
\end{array}\right]=\left[\begin{array}{c}
\left[y_{p1}^L,y_{p1}^U\right]\\
\left[y_{p2}^L,y_{p2}^U\right]\\
\vdots\\
\left[y_{pm}^L,y_{pm}^U\right]
\end{array} \right]
\]
be the actual network output, which depends on the weight vector $\b w$ of the network. The components of $\b w$ are
\[
\mb{\bm$w$}=\left\{w_{ji}^h,w_{sj}^o,\theta_j,\vartheta_s\, : \,i=1,2,\ldots,n,\;j=1,2,\ldots,l, \; s=1,2,\ldots,m\right\}.
\]
Suppose then that we are given a training pair $\left(\mb{\bm$X$}_p,\mb{\bm$Y$}_{dp}\right)$ corresponding to specified $h$-level set. Having the actual output $\mb{\bm$Y$}_p$ and the target, desired, output $\mb{\bm$Y$}_{dp}$, we define a cost function to be minimized, over the space of the network weights, corresponding to the $h$-level sets of the components of $\mb{\bm$Y$}_p$ and $\mb{\bm$Y$}_{dp}$,
\begin{eqnarray}
E_h(\mb{\bm$w$})&=&\frac{1}{2}\sum_{s=1}^{m}\left(\left[\left(Y_s^L\right)_{dp}\right]_h-\left[\left(Y_s^L\right)_p\right]_h\right)^2+\left(\left[\left(Y_s^U\right)_{dp}\right]_h-\left[\left(Y_s^U\right)_p\right]_h \right)^2\nonumber\\
&=&\frac{1}{2}\sum_{s=1}^{m}\left(y^L_{dsh}-y^L_{sh}\right)^2+\left( y^U_{dsh}-y^U_{sh}\right)^2.
\label{eq:costforlevh}
\end{eqnarray}
We can also define a cost function corresponding to level sets. An example of such a cost function is
\[
E(\mb{\bm$w$})=\sum_{h}hE_h(\mb{\bm$w$}).
\]
In the next subsection, we present a training algorithm of the three-layer neural network using interval vectors.

\subsection{Training}
The training of a neural network involves adjusting the weights of the network so that the outputs of the network for a given training pair are as close as possible to the target output in the sense of the given cost function. Thus, we faced with an optimization problem,
\[
\mb{minimize  }E(\mb{\bm$w$})\mb{  over the space of weights}.
\]
The inputs to the network are interval vectors of the form given by~(\ref{eq:bfcapxpinve}) corresponding to the desired outputs given by~(\ref{eq:desinterout}). The output of the $j$-th neuron in the hidden layer is
\[
\left[z_{pj}^L,z_{pj}^U\right]=\left[f_j^h\left(v_{pj}^L\right),f_j^h\left(v_{pj}^U\right)\right],
\]
where
\[
v_{pj}^L=\sum_{{i=1}\atop{w_{ji}^h \ge 0}}^{n}w_{ji}^hx_{pi}^L+\sum_{{i=1}\atop{w_{ji}^h < 0}}
^{n}w_{ji}^hx_{pi}^U+\theta_j,
\]
and
\[
v_{pj}^U=\sum_{{i=1}\atop{w_{ji}^h \ge 0}}^{n} w_{ji}^h x_{pi}^U+\sum_{{ i=1}\atop{w_{ji}^h < 0}}^{n}w_{ji}^hx_{pi}^L+\theta_j.
\]
The output from the $s$-th neuron of the output layer is the interval
\[
\left[y_{ps}^L,y_{ps}^U\right]=\left[f_s^o\left(u_{ps}^L \right),f_s^o\left(u_{ps}^U \right)\right],
\]
where
\[
u_{ps}^L=\sum_{{j=1}\atop{ w_{sj}^o \ge 0}}^{l} w_{sj}^oz_{pj}^L+\sum_{{j=1}\atop{w_{sj}^o < 0}}^{l}w_{sj}^oz_{pj}^U+\vartheta_s,
\]
and
\[
u_{ps}^U=\sum_{{j=1}\atop{ w_{sj}^o \ge 0}}^{l} w_{sj}^oz_{pj}^U+\sum_{{j=1}\atop{w_{sj}^o < 0}}^{l}w_{sj}^oz_{pj}^L+\vartheta_s.
\]
The training procedure of the network involves minimization of a given cost function. We will use $E_h(\mb{\bm$w$})$. To minimize $E_h(\mb{\bm$w$})$, we employ a gradient algorithm with a fixed step size $\eta$. The algorithm takes the form
\[
\mb{\bm$w$}^{(k+1)}=\mb{\bm$w$}^{(k)}-\eta \nabla E_h\left(\mb{\bm$w$}^{(k)}\right).
\]
Thus, we need to compute the partial derivatives of $E_h$ with respect to each component of $\mb{\bm$w$}$. For this, we fix the indices $i$, $j$, and $s$. We first compute the partial derivative of $E_h$ with respect to $w_{sj}^o$. We assume that we have only one training pair, that is, $p=1$, and in further manipulations we will drop the index $p$. There are two cases to analyze. First, we consider the case when $w_{sj}^o \ge 0$. We obtain
\begin{eqnarray}
\frac{\partial E_h}{\partial w_{sj}^o}&=&\frac{\partial}{\partial w_{sj}^o}\left\{\frac{1}{2}\sum_{t=1}^{m}\left(\left[\left(Y_t^L\right)_d\right]_h-\left[Y_t^L\right]_h\right)^2\right\}\nonumber\\
&{}&\quad +\frac{\partial}{\partial w_{sj}^o}\left\{\frac{1}{2}\sum_{t=1}^{m}\left(\left[\left(Y_t^U\right)_d\right]_h-\left[Y_t^U\right]_h\right)^2 \right\}\nonumber\\
&=&\frac{1}{2}\frac{\partial}{\partial Y_s^L}\left(\left[\left(Y_s^L\right)_d\right]_h-\left[Y_s^L\right]_h\right)^2\: \frac{\partial Y_s^L}{\partial u_s^L}\:\frac{\partial u_s^L}{\partial w_{sj}^o}\nonumber\\
&{}&\quad +\frac{1}{2}\frac{\partial}{\partial Y_s^U}\left(\left[\left(Y_s^U\right)_d\right]_h-\left[Y_s^U\right]\right)^2\: \frac{\partial Y_s^U}{\partial u_s^U}\:\frac{\partial u_s^U}{\partial w_{sj}^o}\nonumber\\
&=&-\left(\left[\left(Y_s^L\right)_d\right]_h-\left[Y_s^L\right]\right)f_s^{o'}\left(u_s^L\right)z_j^L\nonumber\\
&{}&\quad -\left(\left[\left(Y_s^U\right)_d\right]_h-\left[Y_s^U\right]\right)f_s^{o'}\left(u_s^U\right)z_j^U,
\label{eq:padeehwrosj}
\end{eqnarray} 
where $f_s^{o'}$ is the derivative of $f_s^o$. Let
\begin{equation}
\delta_s^L=\left(\left[\left(Y_s^L\right)_d\right]_h-\left[Y_s^L\right]\right)f_s^{o'}\left(u_s^L\right)'
\label{eq:deslcno}
\end{equation}
and
\begin{equation}
\delta_s^U=\left(\left[\left(Y_s^U\right)_d\right]_h-\left[Y_s^U\right]\right)f_s^{o'}\left(u_s^U\right).
\label{eq:descauno}
\end{equation}
Then, using the above notation, we represent~(\ref{eq:padeehwrosj}) as
\begin{equation}
\frac{\partial E_h}{\partial w_{sj}^o}\left(\mb{\bm$w$}\right)=-\delta_s^Lz_j^L-\delta_s^Uz_j^U.
\label{eq:deekwosjsim}
\end{equation}
Proceeding in a similar manner for the case when $w_{sj}^o<0$, we obtain
\begin{equation}
\frac{\partial E_h}{\partial w_{sj}^o}\left(\mb{\bm$w$}\right)=-\delta_s^Lz_j^U-\delta_s^Uz_j^L.
\label{eq:wosjnegdeeh}
\end{equation}
We next compute the partial derivatives of $E_h$ with respect to $w_{ji}^h$. We start with the cost function written as
\begin{eqnarray*}
E_h\left(\mb{\bm$w$}\right)&=&\frac{1}{2}\sum_{t=1}^{m}\left(\left[\left(Y_t^L\right)_{d}\right]_h-\left[Y_t^L\right]_h\right)^2 +\frac{1}{2}\sum_{t=1}^{m}\left(\left[\left(Y_t^U\right)_{d}\right]_h-\left[Y_t^U\right]_h \right)^2\\
&=&\frac{1}{2}\sum_{t=1}^{m}\left(\left[\left(Y_t^L\right)_{d}\right]_h-f_t^o\left(\sum_{{q=1}\atop{w_{tq}^o \ge 0}}^{l} w_{tq}^of_q^h\left(v_q^L\right)+\sum_{{q=1}\atop{w_{tq}^o < 0}}^{l} w_{tq}^of_q^h\left(v_q^U\right)+\vartheta_t \right)\right)^2\\
&{}&\quad+\frac{1}{2}\sum_{t=1}^{m}\left(\left[\left(Y_t^U\right)_{d}\right]_h-f_t^o\left(\sum_{{q=1}\atop{w_{tq}^o \ge 0}}^{l} w_{tq}^of_q^h\left(v_q^U\right)+\sum_{{q=1}\atop{w_{tq}^o < 0}}^{l}w_{tq}^o f_q^h\left(v_q^L\right)+\vartheta_t \right)\right)^2.
\end{eqnarray*}
This time we have to consider four cases. We first consider the case when $w_{sj}^o \ge 0$ and $w_{ji}^h \ge 0$. Taking into account the expressions for the inputs to the neurons in the hidden layer and applying the Chain Rule again yields
\begin{eqnarray*}
\frac{\partial E_h}{\partial w_{ji}^h}\left(\mb{\bm$w$}\right)&=&-\sum_{t=1}^{m}\left(\left[\left(Y_t^L\right)_d\right]_h-f_t^o\left(u_t^L\right)\right)f_t^{o'}\left(u_t^L\right)w_{tj}^of_j^{h'}\left(v_j^L\right)x_i^L\\
&{}&\quad -\sum_{t=1}^{m}\left(\left[\left(Y_t^U\right)_d\right]_h-f_t^o\left(u_t^U\right)\right)f_t^{o'}\left(u_t^U\right)w_{tj}^of_j^{h'}\left(v_j^U\right)x_i^U.
\end{eqnarray*}
Using the $\delta_t$ notation, similar to~(\ref{eq:deslcno}) and~(\ref{eq:descauno}), we obtain
\[
\frac{\partial E_h}{\partial w_{ji}^h}\left(\mb{\bm$w$}\right)=-\left(\sum_{t=1}^{m}\delta_t^Lw_{tj}^o\right)f_j^{h'}\left(v_j^L\right)x_i^L
 -\left(\sum_{t=1}^{m}\delta_t^Uw_{tj}^o\right)f_j^{h'}\left(v_j^U\right)x_i^U.
\]
Next, we consider the case when $w_{sj}^o \ge 0$ and $w_{ji}^h < 0$. Using the Chain Rule yields
\[
\frac{\partial E_h}{\partial w_{ji}^h}\left(\mb{\bm$w$}\right)=-\left(\sum_{t=1}^{m}\delta_t^Lw_{tj}^o\right)f_j^{h'}\left(v_j^L\right)x_i^U -\left(\sum_{t=1}^{m}\delta_t^Uw_{tj}^o\right)f_j^{h'}\left(v_j^U\right)x_i^L.
\]
The third case is $w_{sj}^o < 0$ and $w_{ji}^h \ge 0$. Using the Chain Rule again, we obtain
\[
\frac{\partial E_h}{\partial w_{ji}^h}\left(\mb{\bm$w$}\right)=-\left(\sum_{t=1}^{m}\delta_t^Lw_{tj}^o\right)f_j^{h'}\left(v_j^U\right)x_i^U -\left(\sum_{t=1}^{m}\delta_t^Uw_{tj}^o\right)f_j^{h'}\left(v_j^L\right)x_i^L.
\]
Finally, we consider the case when $w_{sj}^o < 0$ and $w_{ji}^h < 0$. For this case, we get
\[
\frac{\partial E_h}{\partial w_{ji}^h}\left(\mb{\bm$w$}\right)=-\left(\sum_{t=1}^{m}\delta_t^Lw_{tj}^o\right)f_j^{h'}\left(v_j^L\right)x_i^L -\left(\sum_{t=1}^{m}\delta_t^Uw_{tj}^o\right)f_j^{h'}\left(v_j^L\right)x_i^U.
\]

We now compute the partial derivatives of $E_h$ with respect to the bias parameters $\vartheta_s$ and $\theta_j$. The partial derivative of $E_h$ with respect to $\vartheta_s$ is
\[
\frac{\partial E_h}{\partial \vartheta_s}\left(\mb{\bm$w$}\right)=-\delta_s^L-\delta_s^U.
\]

We consider two cases when computing partial derivatives of $E_h$ with respect to $\theta_j$. In the case when $w_{sj}^o \ge 0$, we have
\[
\frac{\partial E_h}{\partial \theta_j}\left(\mb{\bm$w$}\right)=-\sum_{t=1}^{m}\delta_t^Lw_{tj}^of_j^{h'}\left(v_j^L\right)-\sum_{t=1}^{m}\delta_t^Uw_{tj}^of_j^{h'}\left(v_j^U\right).
\]
In the case when $w_{sj}^o < 0$,
\[
\frac{\partial E_h}{\partial \theta_j}\left(\mb{\bm$w$}\right)=-\sum_{t=1}^{m}\delta_t^Lw_{tj}^of_j^{h'}\left(v_j^U\right)-\sum_{t=1}^{m}\delta_t^Uw_{tj}^of_j^{h'}\left(v_j^L\right).
\]

We now summarize the gradient algorithm for updating the weights of the fuzzy neural network. We will write the update equations for the two sets of weights $w_{sj}^o$ and $w_{ji}^h$ separately. We have
\[
w_{sj}^{o(k+1)}=w_{sj}^{o(k)}+\left\{\begin{array}{c}
\eta \left(\delta_s^{L(k)} z_j^{L(k)}+\delta_s^{U(k)} z_j^{U(k)}\right) \quad\mb{if}\quad w_{sj}^{o(k)} \ge 0\\
\eta \left(\delta_s^{L(k)} z_j^{U(k)}+\delta_s^{U(k)} z_j^{L(k)}\right) \quad\mb{if}\quad w_{sj}^{o(k)} < 0.
\end{array}
\right.
\]
The update equations for the weights of the neurons in the hidden layer are:
\begin{eqnarray*}
w_{ji}^{h(k+1)}&=&w_{ji}^{h(k)}\\
&{}&+\left\{\begin{array}{c}
\eta\left( \left(\sum_{t=1}^{m}\delta_t^{L(k)}w_{tj}^{o(k)}\right)f_j^{h'}\left(v_j^{L(k)}\right)x_i^L+\left(\sum_{t=1}^{m}\delta_t^{U(k)}w_{tj}^{o(k)}\right)f_j^{h'}\left(v_j^{U(k)}\right)x_i^U\right)\\
\mb{for} \quad w_{sj}^{o(k)} \ge 0 \; w_{ji}^{h(k)} \ge 0,\\
\eta \left( \left(\sum_{t=1}^{m}\delta_t^{L(k)}w_{tj}^{o(k)}\right)f_j^{h'}\left(v_j^{L(k)}\right)x_i^U+\left(\sum_{t=1}^{m}\delta_t^{U(k)}w_{tj}^{o(k)}\right)f_j^{h'}\left(v_j^{U(k)}\right)x_i^L\right)\\
\mb{for} \quad w_{sj}^{o(k)} \ge 0 \; w_{ji}^{h(k)} < 0,\\
\eta\left( \left(\sum_{t=1}^{m}\delta_t^{L(k)}w_{tj}^{o(k)}\right)f_j^{h'}\left(v_j^{U(k)}\right)x_i^U+\left(\sum_{t=1}^{m}\delta_t^{U(k)}w_{tj}^{o(k)}\right)f_j^{h'}\left(v_j^{L(k)}\right)x_i^L\right)\\
\mb{for} \quad w_{sj}^{o(k)} < 0 \; w_{ji}^{h(k)} \ge 0,\\
\eta \left( \left(\sum_{t=1}^{m}\delta_t^{L(k)}w_{tj}^{o(k)}\right)f_j^{h'}\left(v_j^{L(k)}\right)x_i^L+\left(\sum_{t=1}^{m}\delta_t^{U(k)}w_{tj}^{o(k)}\right)f_j^{h'}\left(v_j^{L(k)}\right)x_i^U\right)\\
\mb{for} \quad w_{sj}^{o(k)} < 0 \; w_{ji}^{h(k)} < 0.
\end{array}
\right.
\end{eqnarray*}
The update equations for the bias parameters $\vartheta_s$ are:
\[
\vartheta_s^{(k+1)}=\vartheta_s^{(k)}+\eta\left( \delta_s^{L(k)}+\delta_s^{U(k)}\right).
\]
Finally, the update equations for the bias parameters $\theta_j$ are:
\begin{eqnarray*}
\theta_j^{(k+1)}&=&\theta_j^{(k)}\\
&{}&+\left\{\begin{array}{c}
\eta \left(\sum_{t=1}^{m}\delta_t^{L(k)}w_{tj}^{o(k)}f_j^{h'}\left(v_j^{L(k)}\right)+\sum_{t=1}^{m}\delta_t^{U(k)}w_{tj}^{o(k)}f_j^{h'}\left(v_j^{U(k)}\right)\right) \; \mb{for} \; w_{sj}^{o(k)}\ge0\\
\eta\left(\sum_{t=1}^{m}\delta_t^{L(k)}w_{tj}^{o(k)}f_j^{h'}\left(v_j^{U(k)}\right)+\sum_{t=1}^{m}\delta_t^{U(k)}w_{tj}^{o(k)}f_j^{h'}\left(v_j^{L(k)}\right) \right) \; \mb{for} \; w_{sj}^{o(k)}<0.
\end{array}
\right.
\end{eqnarray*}

\section{Neural Fuzzy Modeling}
In this section, we describe a method that can be used to model nonlinear systems as well as to construct neural-network-based fuzzy logic control strategies for nonlinear/uncertain dynamical systems. In the controller synthesis, employing this method, one can use numerical data rather than fuzzy IF-THEN rules. The above problem was considered by Lin and Lee in~\cite{lin91}, \cite{lin94}, and in~\cite{lin96}, by Wang~\cite{lxwang97}, and by Kosko~\cite{kosko97}. Related approaches to fuzzy rule-based neural network controllers can be found in the papers by Higgins and Goodman~\cite{higgins94} and Nie and Linkens~\cite{nie94}. In the following subsection, we show how fuzzy logic and the backpropagation method can be combined to yield a neural fuzzy identifier.

\subsection{Neural Fuzzy Identifier}
A typical fuzzy logic controller is composed of three basic components as depicted in Figure~\ref{fig:typflc}.
\begin{figure}
\vspace{3.5in}
\hskip1in\special{bmp:c:/mybook/nufu_id.bmp x=4.67in y=3.5in}
\caption{An identification scheme of a nonlinear dynamical system using neural fuzzy identifier.}
\label{fig:typflc}
\end{figure}
The fuzzifier maps crisp inputs into fuzzy numbers. The decision making logic unit, the inference engine, uses fuzzy rules to infer fuzzy control actions. The defuzzifier maps the fuzzy outputs from the decision making logic block into crisp values.

In Chapter~\ref{chp:fuzzysys}, we discussed a method for designing a fuzzy logic controller using given to us fuzzy IF-THEN rules that describe in linguistic terms the given dynamic system operation. The availability of such rules is a basic requirement for implementing a fuzzy controller. However, there may be cases where neither a mathematical model of the process nor its linguistic description are available. In such cases we may try to extract fuzzy IF-THEN rules from example data, and then proceed with an implementation of a suitable controller. We will use the available data to train fuzzy feedforward neural-network-based controller shown in Figure~\ref{fig:hybrnn}.
\begin{figure}
\vspace{3.5in}
\hskip1.4in\special{bmp:c:/mybook/hybrnn.bmp x=4.67in y=3.5in}
\caption{Fuzzy neural-network-based controller structure}
\label{fig:hybrnn}
\end{figure}
The network consists of four layers. In the first layer the input signals $x_i$, $i=1,2,\ldots ,n$, are branched-off to the second, hidden, layer. We refer to this layer as the {\em input membership layer}. It implements the input membership functions. The membership functions are bell-shaped and are described as
\[
\mu_{ji}^h=\exp\left\{-\left(\frac{x_i-m_{ji}^h}{\sigma_{ji}^h}\right)^2\right\},
\]
where $j=1,2,\ldots, p$, the parameter $m_{ji}^h$ is the mean of $\mu_{ji}^h$, and the parameter $\sigma_{ji}^h$, the variance, is responsible for the shape of $\mu_{ji}^h$.\\
The outputs of the input membership function layer feed the {\em rule layer}. This layer contains one node for each rule. Each rule node of the layer performs a product of its inputs.
The outputs of the rule nodes are connected via {\em output membership functions} nodes to the {\em defuzzification layer} that generates crisp output signals of the form
\[
y=\frac{\sum_{k=1}^{r}u_k m_k^o/ \sigma_k^o}{\sum_{k=1}^{r}u_k /\sigma_k^o},
\]
where
\[
u_k=\prod_{j=1}^{p}\mu_{kj}^h
\] 
is the output from the $k$th output membership function node, and
\[
\mu_k^o=\exp \left\{-\left(\frac{v-m_k^o}{\sigma_k^o}\right)^2 \right\}.
\]
Observe that the above network implements the fuzzy logic controller discussed at the beginning of the section.

We now present an approach to training the above neural-network-based controller. We discuss the training process next.
\subsection{Training Algorithm}
To train the network, we will use the backpropagation algorithm. The cost function to be minimized is
\[
E=\frac{1}{2}\left(y_d-y\right)^2,
\]
where $y_d$ is the desired output and $y$ is the actual output from the network. We will use a gradient algorithm to minimize $E$ over all
\[
\mb{\bm$w$}=\left\{m_{ji}^h, \: \sigma _{ji}^h,\: m_l^o,\:\sigma_l^o \right\}.
\]
The gradient algorithm we will use has the form
\[
\mb{\bm$w$}^{(k+1)}=\mb{\bm$w$}^{(k)}-\eta \nabla E\left(\mb{\bm$w$}^{(k)}\right),
\]
where $\eta$ is a fixed step size.

In order to be able to use the above algorithm we need to compute the partial derivatives of $E$ with respect to each component of $\mb{\bm$w$}$. For this we first fix indices $i$, $j$, and $l$. We first compute the partial derivatives of $E$ with respect to $m_s^o$ and $\sigma_s^o$. We obtain
\begin{eqnarray*}
\frac{\partial E}{\partial m_s^o}\left(\mb{\bm$w$}\right)&=&\frac{\partial}{\partial m_s^o}\left(\frac{1}{2}\left(y_d-y\right)^2\right)\\
&=&\frac{\partial}{\partial m_s^o}\left(\frac{1}{2}\left(y_d-\frac{\sum_{l=1}^{r}u_l m_l^o/ \sigma_l^o}{\sum_{l=1}^{r} u_l /\sigma_l^o}\right)^2\right)\\
&=&-\left(y_d-y\right)\frac{u_s/ \sigma_s^o}{\sum_{l=1}^{r}u_l / \sigma_l^o},
\end{eqnarray*}
and
\[
\frac{\partial E}{\partial \sigma_s^o}\left(\mb{\bm$w$}\right)=-\left(y_d-y\right)\frac{-\left(u_sm_s^o/\left(\sigma _s ^o \right)^2\right)\left(\sum_{l=1}^{r}u_l/ \sigma_l^o\right)+u_s/\left(\sigma _s^o\right)^2\sum_{l=1}^{r}u_lm_l^o /\sigma_l^o}{\left(\sum_{l=1}^{r}u_l / \sigma_l^o\right)^2}.
\]
For simplicity of notation, let
\[
\delta_s=\left(y_d-y\right)\frac{-\left(u_sm_s^o/\left(\sigma _s ^o \right)^2\right)\left(\sum_{l=1}^{r}u_l/ \sigma_l^o\right)+u_s/\left(\sigma _s^o\right)^2\sum_{l=1}^{r}u_lm_l^o /\sigma_l^o}{\left(\sum_{l=1}^{r}u_l / \sigma_l^o\right)^2}.
\]
Then
\[
\frac{\partial E}{\partial \sigma_s^o}\left(\mb{\bm$w$}\right)=-\delta_s.
\]
We next compute the partial derivatives of $E$ with respect to $m_{ji}^h$, and $\sigma _{ji}^h$. For this we represent $E$ as
\begin{eqnarray*}
E\left(\mb{\bm$w$}\right)&=&\frac{1}{2}\left(y_d-\frac{\sum_{l=1}^{r}u_l m_l^o / \sigma_l^o}{\sum_{l=1}^{r} u_l / \sigma_l^o}\right)^2\\
&=&\frac{1}{2}\left(y_d-\frac{\sum_{l=1}^{r}m_l^o / \sigma_l^o \prod_{j=1}^{p}\mu_{lj}^h}{\sum_{l=1}^{r}1/\sigma_l^o \prod_{j=1}^{p}\mu_{lj}^h}\right)^2\\
&=&\frac{1}{2}\left(y_d-\frac{\sum_{l=1}^{r}m_l^o / \sigma_l^o \prod_{j=1}^{p} \exp\left\{-\left(\frac{x_{j_i}-m_{lj}^h}{\sigma_{lj}^h}\right) ^2\right\}}{\sum_{l=1}^{r} 1/\sigma_l^o \prod_{j=1}^{p}\exp\left\{-\left(\frac{x_{j_i}-m_{lj}^h}{\sigma_{lj}^h}\right) ^2 \right\}}\right)^2,
\end{eqnarray*}
where the notation $x_{j_i}$ is used to indicate the way $x_i$s enter $\mu_{lj}^h$ in order to extract from example data a meaningful set of fuzzy IF-THEN rules.\\
Using the Chain Rule again, we now compute the partial derivative of $E$ with respect to $m_{st}^h$ to get
\begin{eqnarray*}
\frac{\partial E}{\partial m_{st}^h}\left(\mb{\bm$w$}\right)&=&-\left(y_d-y\right)\frac{\partial}{\partial u_s}\left( \frac{\sum_{l=1}^{r}u_l m_l^o /\sigma _l^o}{\sum_{l=1}^{r}u_l / \sigma_l^o} \right) \frac{\partial u_s}{\partial m_{st}^h}\\
&=&-\left(y_d-y\right)\frac{m_s^o/\sigma_s^o\left(\sum_{l=1}^{r}u_l /\sigma_l^o\right)-\left(1/\sigma_s^o \right)\left(\sum_{l=1}^{r}u_lm_l^o/ \sigma_l^o \right)}{\left(\sum_{l=1}^{r}u_l/ \sigma_l^o \right)^2}\frac{\partial u_s}{\partial m_{st}^h}.
\end{eqnarray*}
Let
\[
\delta_s^h=\left(y_d-y\right)\frac{m_s^o/\sigma_s^o\left(\sum_{l=1}^{r}u_l /\sigma_l^o\right)-\left(1/\sigma_s^o \right)\left(\sum_{l=1}^{r}u_lm_l^o/ \sigma_l^o \right)}{\left(\sum_{l=1}^{r}u_l/ \sigma_l^o \right)^2}.
\]
Then
\[
\frac{\partial E}{\partial m_{st}^h}\left(\mb{\bm$w$}\right)=-\delta_s^h\frac{\partial u_s}{\partial m_{st}^h}.
\]
We continue applying the chain rule to get
\begin{eqnarray*}
\frac{\partial E}{\partial m_{st}^h}\left(\mb{\bm$w$}\right)&=&-\delta_s^h\prod_{{j=1}\atop{j \ne t}}^{p}\mu_{sj} \, \frac{\partial}{\partial m_{st}^h}\left(\exp\left\{-\left(\frac{x_{j_i}-m_{st}^h}{\sigma_{st}^h}\right)^2\right\}\right)\\
&=&-2\delta_s^h u_s \frac{x_{j_i}-m_{st}^h}{\left(\sigma_{st}^h\right)^2}.
\end{eqnarray*}
Using the Chain Rule again, we compute the partial derivative of $E$ with respect to $\sigma_{st}^{h}$ to obtain
\[
\frac{\partial E}{\partial \sigma_{st}^h}\left(\mb{\bm$w$}\right)=-2\delta_s^h u_s \frac{\left(x_{ji}-m_{st}^h \right)^2}{\left(\sigma_{st}^h\right)^3}.
\] 
We are now ready to formulate the gradient algorithm for updating the weights of the fuzzy neural-network-based controller. We have
\begin{eqnarray*}
m_s^{o(k+1)}&=&m_s^{o(k)}+\eta \left(y_d-y^{(k)}\right)\frac{u_s^{(k)} / \sigma_s^{o(k)}}{\sum_{l=1}^{r}u_l^{(k)} /\sigma_l^{o(k)}},\\
\sigma_s^{o(k+1)}&=&\sigma_s^{(k)}+\eta \delta_s^{(k)},\\
m_{st}^{h(k+1)}&=&m_{st}^{h(k)}+2\eta\delta_s^{h(k)} u_s^{(k)} \frac{x_{j_i}-m_{st}^{h(k)}}{\left(\sigma_{st}^{h(k)}\right)^2},\\
\sigma_{st}^{h(k+1)}&=&\sigma_{st}^{h(k)}+2 \eta \delta_s^{h(k)} u_s^{(k)} \frac{\left(x_{ji}-m_{st}^{h(k)} \right)^2}{\left(\sigma_{st}^{h(k)}\right)^3}.
\end{eqnarray*}

\section{Radial-Basis Function Networks}
In this section, we present radial-basis function~(RBF) neural networks and their applications to interpolation and then to identification of unknown nonlinear dynamical systems using only the systems' inputs and outputs data. We first show how an RBF network can be used for interpolation purposes. Then, we will 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 is an adaptive system for which the learning algorithm is developed using Lyapunov techniques. The algorithm allocates resources, by adding new nodes, to obtain the desired identification accuracy.

\subsection{Interpolation Using RBF Networks}
A diagram of a radial-basis function~(RBF) network is shown in Figure~\ref{fig:grbf_ip}.
\begin{figure}
\vspace{3in}
\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 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\sigma_i^2} \right),
\label{eq:grbfunct}
\end{equation}
where $\b x^{(i)}$ is the center of the basis function and $\sigma_i$ is sometimes called the radius of the radial basis function. The function given by~(\ref{eq:grbfunct}) is called the Gaussian radial basis function.

The RBF networks are used to perform interpolation. The interpolation problem 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}
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].
\]
Then, 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 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 for the interpolation data given in Table~\ref{tab:interdata}.
\begin{table}[h]
\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 set $\sigma_i=1$, $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]^T.
\]
The interpolation result is depicted in Figure~\ref{fig:grbf_eip}.
\begin{figure}
\vspace{3.5in}
\hskip.9in\special{bmp:c:/mybook/grbf_eip.bmp x=4.67in y=3.5in}
\caption{Results of interpolating using radial-basis functions.}
\label{fig:grbf_eip}
\end{figure}
In Figure~(\ref{fig:grbf_eip})(a), we used solid line to show the result of interpolating using the Gaussian radial-radial basis functions, while 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}

\subsection{Identification of a Single-Input Single-State System}
We first discuss the identification of the single-input single-state dynamical system using neural network. Let the dynamical system be 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$ to be identified by a neural network. 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 the learning rules. We will now proceed with the construction of a neural identifier of the dynamical system~(\ref{eq:sissdynsid}). The identifier will be 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. 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,u,\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. Suppose now that the basis functions are given to us. The problem of choosing the basis functions will be discussed later. Suppose also that we know the ``optimal'' weights $w_k^*$, $k=1,2,\ldots, K$, that result in the 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}) as
\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, 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 is depicted in Figure~\ref{fig:grbf_id}.
\begin{figure}
\vspace{3in}
\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 next define the state estimation error to be
\[
e_x=x-\hat{x},
\]
while the weights' estimation errors,
\[
\xi_k=w^*_k-w_k,\quad k=1,\ldots,K.
\]
We assume that the state estimation error, $e_x$, is bounded by some finite constant $e_K$, that is, $\left|e_x\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}
Our goal now is to use the second method of Lyapunov to devise a learning law for the RBF identifier that would yield stable identification scheme in the sense that the state error $e_x$ decreases as the time increases and the weights do not diverge. For this purpose, we will 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]^T$.
We now evaluate the time derivative of $V$ on the trajectories of~(\ref{eq:errdyneq}),
\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}
Since $\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 the above that if there is no state estimation error, that is, $e_K=0$, then $\dot{V}$ is negative semi-definite in the $\left[\begin{array}{cc}
e_x & \b \xi ^T
\end{array}\right]^T$-space. In such a case the origin of the $\left[\begin{array}{cc}
e_x & \b \xi ^T
\end{array}\right]^T$-space is stable, which guarantees boundedness of the weights $w_k$s. In the presence of the state estimation error, 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 since $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]^T$. 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 always 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$. But the weights, as well as the estimation error, $e_x$,  are upper-bounded. When any of the two bounds is reached, more basis functions maybe needed to improve the approximation of $g$. New basis functions are added to reduce the error until $\dot{V}<0$. We now describe an algorithm 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 error and the bounds are reached, additional basis functions in the RBF network should be added. 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 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)$ 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, since 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}
\mb{round}\left( x(t)/\delta_x\right)\delta_x&\mb{round}\left( u(t)/\delta_u\right)\delta_u
\end{array}\right]^T,
\label{eq:cenofnbfun}
\end{equation}
where the round operator rounds the elements on which it operates to the nearest integer. For example, let $\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 as 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 radius of the new basis function and the adjusted radii are 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\atop \b m_i\ne \b m_k}\left\{ \left\Arrowvert \b m_k-\b m^*_k\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| x(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 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 condition~(\ref{eq:thcofoabf}) is satisfied, a new basis function is added to reduce the state estimation error.

\begin{example}\rm
The above technique for identification of dynamical systems using growing radial-basis function network is illustrated with the following 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)$. 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 use the following parameter values: $e_0=0.001$, $e_{\max}=0.005$, $\alpha=5$, $\kappa=3$. A simple MATLAB code that was used to perform our simulations appears below.

\begin{verbatim}
t=1;dt=0.01;w(1)=0;
i=0;g(1)=0;xhat(1)=0;
emax=0.005;e0=0.001;
dx=0.05;du=0.05;temp=0;
alpha=5;kappa=3;
r(1)=kappa*dx;
% we use the notation: 2*sigma^2_k=r^2_k, that is, in our 
% subsequent expression for the Gaussian radial-basis function,
% we skip the factor 2
x(1)=0;
m1(1)=x(1);
m2(1)=cos(0);
tfinal=1000;
while t<tfinal
    u(t)=cos(t*dt);
    xdot(t)=1.1*(1-x(t)+2*x(t)*u(t)-x(t)^2)*exp(-0.5*x(t)^2-u(t)^2);
    x(t+1)=x(t)+xdot(t)*dt;
      for z=1:100
        g(t)=0;
        if i>0
           for n=1:i
             g(t)=g(t)+w(n)*exp(-((x(t)-m1(n))^2+(u(t)-m2(n))^2)/(r(n)^2));
           end %end for
        end %end if
        xhatdot(t)=-0.5*xhat(t)+g(t);
        xhat(t+1)=xhat(t)+xhatdot(t)*dt;
        e(t)=x(t)-xhat(t);
        xmin=abs(x(t)-m1(1));
        umin=abs(u(t)-m2(1));
        if i>0
           if e0 <= abs(e(t))
               if norm(w) <= 1.5*sqrt(i)
                   for n=1:i
            wdot=alpha*e(t)*exp(-((x(t)-m1(n))^2+(u(t)-m2(n))^2)/(r(n)^2));
                      w(n)=w(n)+wdot*dt;
                      if abs(x(t)-m1(n))<xmin
                           xmin=abs(x(t)-m1(n));
                      end %end if
                      if abs(u(t)-m2(n))<umin
                           umin=abs(u(t)-m2(n));
                      end
                   end
               end
           end
        end
        if abs(e(t))>emax
           if xmin > (0.5*dx)  
               temp=i+1;
           end
           if umin > (0.5*du)
               temp=i+1;
           end
        end
        if i<temp
           i=temp;
           m1(i)=round(x(t)/dx)*dx;
           m2(i)=round(u(t)/du)*du;
           w(i)=0;
           r(i)=3*dx;
           for a=1:i
               m1min=0;
               m2min=0;
                 if a~=i
                      mk=sqrt((m1(a)-m1(a+1))^2+(m2(a)-m2(a+1))^2);
                      m1min=m1(a+1);
                      m2min=m2(a+1);
                 elseif i~=1
                      mk=sqrt((m1(a)-m1(a-1))^2+(m2(a)-m2(a-1))^2);
                      m1min=m1(a-1);
                      m2min=m2(a-1);
                 end
     for b=1:i
       if a~=b
         if mk>sqrt((m1(a)-m1(b))^2+(m2(a)-m2(b))^2)
             mk=sqrt((m1(a)-m1(b))^2+(m2(a)-m2(b))^2);
             m1min=m1(b);
             m2min=m2(b);
         end
       end
      end
       r(a)=kappa*sqrt((m1(a)-m1min)^2+(m2(a)-m2min)^2);
      end
    w(i)=dt*alpha*e(t)*exp(-((x(t)-m1(i))^2+(u(t)-m2(i))^2)/(r(i)^2));
   end
end
t=t+1;
end  
plot(1:tfinal,x,1:tfinal,xhat)
xlabel('Iteration','Fontsize',16)
ylabel('x and its estimate','Fontsize',16)
\end{verbatim}

We mention that the \verb!z! learning loop from 1 to 100 was added to allow the approximation to learn out the error. Using the above program, we generated a number of plots illustrating the performance of the RBF neural identifier. The last three lines in the MATLAB code show how we obtained the plots of Figure~\ref{fig:grbf_xes}. The plots in the remaining figures were obtained using similar commands. In Figure~\ref{fig:grbf_xes}, plots of the true state and its estimate versus time are shown. 
\begin{figure}
\vspace{3.5in}
\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 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}
\vspace{3.5in}
\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}
\vspace{3.5in}
\hskip1.2in \special{bmp:c:/ee675/grbf_uxs.bmp x=3.97in y=3.5in}
\caption{A plot of the true 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}
\vspace{3.5in}
\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 above concepts 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 that we present 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=\mb{trace}\left(\b{\Gamma}^T\b{\Gamma} \right)=\sum_{i=1}^m\sum_{j=1}^n \gamma_{ij}^2.
\]
If the elements, $\gamma_{ij}$, of $\b{\Gamma}$ are functions of time, then
\begin{equation}
\frac{d}{dt}\left( \mb{trace}\left(\b{\Gamma}^T\b{\Gamma} \right)\right)=2\sum_{i=1}^m\sum_{j=1}^n \gamma_{ij}\dot{\gamma}_{ij}=2\:\mb{trace}\left( \b{\Gamma}^T\dot{\b {\Gamma}}\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 nonlinear 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 right-half plane. Since $\b A$, then by the theorem of Lyapunov, which is 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^T\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 deal with 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)).
\]
Let the $k$-th basis function have the form
\[
\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),
\]
and let
\[
\b{\Phi}_K=\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]^T
\]
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}^T & \b m_{\b u k}^T
\end{array}\right]^T$ is the center of the $k$-th radial basis function in the $\left[\begin{array}{cc}
\b x^T & \b u^T
\end{array}\right]^T$-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{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. We have
\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 use 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}^T\b P\b e_{\b x} +\frac{1}{\alpha}\mb{trace}\left(\b{\Gamma}_K\b{\Gamma}_K^T\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}^T\b P\dot{\b e}_{\b x}+\frac{1}{\alpha}\mb{trace}\left( \b{\Gamma}_K\dot{\b{\Gamma}}_K^T\right)\nonumber\\
&=&\b e_{\b x}^T\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^T\right)\nonumber\\
&=& \frac{1}{2}\b e_{\b x}^T\left(\b A^T\b P+\b P\b A\right)\b e_{\b x}+\b e_{\b x}^T\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^T\right)+\b e_{\b x}^T\b P\b e.
\label{eq:lyadermimser}
\end{eqnarray}
Note that by~(\ref{eq:two_mat_trace}) for non-square matrices, we have
\begin{eqnarray}
\b e_{\b x}^T\b P\b{\Gamma}_K\b{\Phi}_K(\b x,\b u)&=&\mb{trace}\left( \b e_{\b x}^T\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}^T\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}^T\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}^T\b P\right)-\mb{trace}\left( \b{\Gamma}_K\dot{\b{W}}_K^T\right)\right)+\b e_{\b x}^T\b P\b e.
\label{eq:inclyanpg}
\end{equation}
We propose to choose the following learning rule,
\begin{equation}
\fbox{$\displaystyle \dot{\b W}_K=\alpha\b P\b e_{\b x}\b{\Phi}_K^T(\b x,\b u) $}
\label{eq:lerufomacw}
\end{equation}
With the above weight matrix update law, the Lyapunov derivative~(\ref{eq:inclyanpg}) evaluates to
\begin{eqnarray}
\dot{V}&=&-\b e_{\b x}^T\b Q\b e_{\b x}+\b e_{\b x}^T\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}
It follows from the above that 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$. This 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^T(\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}


\section{Hopfield Neural Network Model}

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}
\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}
\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}
\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}
\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_1-u_n}{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.

In further analysis, we will 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}
\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]^T\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^T\b T\b x-\b v^T\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^T$. 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]^T.
\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~(\ref{eq:matreofhom}). We have
\begin{equation}
\fbox{$\displaystyle\b S(\b x)\dot{\b x}(t)=-\nabla E(\b x(t)) $}
\label{eq:alagrsys}
\end{equation}

We will analyze the dynamics of the Hopfield network modeled by~(\ref{eq:alagrsys}) using the computational energy function. We will 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^T$, 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)^T\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}^T(t)\b S(\b x(t))\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
\[
\frac{d g^{-1}(x_i)}{d x_i}=\theta\left(1/x_i+1/(1-x_i) \right)>0
\]
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}).
\section{Analog-to-Digital Converter}
In this section, we show how the Hopfield neural network can be used to solve some engineering problems. Specifically, we 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, we 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.

The first step of rounding a number to its nearest integer by a neural network can be organized as solving an optimization problem by the neural network. The second step of converting the integer representation  from integer to binary can be accomplished using a neural network with a suitable nonlinear activation function mapping real numbers into the interval $[0, 1]$. 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$.

Since A/D conversion is viewed as an optimization problem, we need to formulate a suitable cost, or penalty, function that we will minimize. In fact finding a suitable cost function will be the first step in constructing a neural network based 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 network based 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}
\vspace{4.5in}
\hskip.5in \special{eps:c:/brain/ad_conv.eps x=5.1in}
\caption{Constructing the cost function for a 2-bit A/D converter; plots of the cost function and its components.}
\label{fig:ad_conv}
\end{figure}
Observe that minimizing~(\ref{eq:quacost}) corresponds to choosing $V_1$ and $V_0$ in a least-squares sense. However, constructing a circuit minimizing the cost function~(\ref{eq:quacost}) would not yield good results since $V_0$ and $V_1$ would take value 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 of this landscape, 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 will 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}
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, we have
\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 since $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 used the activation function
\[
g(V_i)=\frac{1}{1+e^{-4V_i}}.
\]
Thus,
\[
\frac{d g^{-1}(V_i)}{d V_i}=\frac{1}{4 V_i(1-V_i)}.
\]
\begin{figure}
\vspace{4.5in}
\hskip.75in\special{eps:c:/brain/adphase.eps x=5.1in}
\caption{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'', which are both valid solutions for $x=1.5$. 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.

Using the above technique, we can easily 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{\partial E}{\partial V_i}\left( \frac{d g_i^{-1}(V_i)}{d V_i}\right)^{-1}.
\]

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} 

We apply the second method of Lyapunov to analyze stability in the sense of Lyapunov of the Hopfield neural network. The neural network will be viewed as a large interconnected dynamical system, that is, as an interconnection of many neurons. We first analyze the stability of a single neuron, and then proceed with the stability analysis of the whole network in terms of the properties of the single neurons.

We first review some, not particularly well known, results from matrix theory that we use in our presentation. In particular, we review selected properties of M-matrices that are useful in the analysis of large scale interconnected systems. Next, we present the Comparison Theorem which we use when estimating trajectory bounds of neural networks. Having reviewed the relevant background material, we start with the analysis of a single neuron. After that, we discuss the stability and instability of the whole network in terms of the properties of the individual neurons.

\subsection{Properties of M-Matrices}

We start with the definition of the M-matrix.
\begin{defn}
A real $n \times n$ matrix $\mb{\bm$A$}$ is an M-matrix if all its off-diagonal elements are non-positive, that is, $a_{ij} \le 0$, $i \ne j$, and all its principal minors are positive.
\end{defn}
Note that if $\mb{\bm$A$}$ is an M-matrix, then its diagonal elements are positive since its principal minors are positive.\\
In establishing the needed properties of the M-matrices, we use the following technical lemma due to Fiedler and Pt\'{a}k~\cite{fiedler62}.
\begin{lemma}
\label{lemma:lu}
The following properties of an arbitrary matrix $\mb{\bm$A$} \in \mb{\Bb R}^{n \times n}$ are equivalent:
\begin{itemize}
\item[(i)] The leading principal minors of $\mb{\bm$A$}$ are positive,
\item[(ii)] There exists a lower triangular matrix $\mb{\bm$L$} \in \mb{\Bb R}^{n \times n}$ and an upper triangular matrix $\mb{\bm$U$} \in \mb{\Bb R}^{n \times n}$ both with positive diagonal elements such that $\mb{\bm$A$}=\mb{\bm$L$}\mb{\bm$U$}$
\end{itemize}
\end{lemma}
{\bf Proof}\quad We first prove the implication (i)$ \Rightarrow $(ii). The proof is by induction. Obviously the implication holds for $n=1$. Suppose the statement holds for all real matrices of order $n-1$. Let us represent $\mb{\bm$A$} \in \mb{\Bb R}^{n \times n}$ as follows:
\[
\mb{\bm$A$} = \left[\begin{array}{cc}
\mb{\bm$A$}_{n-1} & \mb{\bm$a$}\\
\b b^T & a_{nn}
\end{array} \right],
\]
where $\mb{\bm$A$}_{n-1} \in \mb{\Bb R}^{(n-1) \times (n-1)}$, $\mb{\bm$a$} \in \mb{\Bb R}^{(n-1) \times 1}$, $\b b^T \in \mb{\Bb R}^{1 \times (n-1)}$, and $a_{nn} \in \mb{\Bb R}$. From the inductive hypothesis, the matrix $\mb{\bm$A$}_{n-1}$ has positive leading principal minors, $\mb{\bm$A$}_{n-1}=\mb{\bm$L$}_{n-1}\mb{\bm$U$}_{n-1}$, where $\mb{\bm$L$}_{n-1} \in \mb{\Bb R}^{(n-1) \times (n-1)}$ is a lower triangular matrix with positive diagonal elements, and $\mb{\bm$U$}_{n-1} \in \mb{\Bb R}^{(n-1) \times (n-1)}$ is an upper triangular matrix with positive diagonal elements. Observe that
\[
\left[\begin{array}{cc}
\mb{\bm$I$}_{n-1} & \mb{\bm$0$}\\
-\b b^T \mb{\bm$A$}_{n-1}^{-1} & 1
\end{array} \right] \mb{\bm$A$} = \left[\begin{array}{cc}
\mb{\bm$A$}_{n-1} & \mb{\bm$a$}\\
\mb{\bm$0$} & a_{nn} - \b b^T\mb{\bm$A$}_{n-1}^{-1}\mb{\bm$a$}
\end{array}\right],
\]
where from the assumptions $\det \mb{\bm$A$}_{n-1}>0$, and $\det \mb{\bm$A$}>0$. Hence,
\[
0< \det \mb{\bm$A$}=\left( a_{nn} - \b b^T\mb{\bm$A$}_{n-1}^{-1}\mb{\bm$a$} \right) \det \mb{\bm$A$}_{n-1}.
\]
Thus, we have
\[
a_{nn} - \b b^T\mb{\bm$A$}_{n-1}^{-1}\mb{\bm$a$}>0.
\]
Let
\[
\mb{\bm$L$}=\left[\begin{array}{cc}
\mb{\bm$L$}_{n-1} & \mb{\bm$0$}\\
\b b^T \mb{\bm$U$}_{n-1}^{-1} & 1
\end{array} \right],
\]
and
\[
\mb{\bm$U$}=\left[\begin{array}{cc}
\mb{\bm$U$}_{n-1} & \mb{\bm$L$}_{n-1}^{-1}\mb{\bm$a$}\\
\mb{\bm$0$} & a_{nn} - \b b^T \mb{\bm$A$}_{n-1}^{-1}\mb{\bm$a$}
\end{array}\right].
\]
One can check that
\[
\mb{\bm$A$}=\mb{\bm$L$}\mb{\bm$U$},
\]
and that the matrices $\mb{\bm$L$}$ and $\mb{\bm$U$}$ satisfy~(ii).

We now prove the implication (i)$ \Leftarrow $(ii). From the generalized Cauchy's Theorem it follows that the leading principal minors of $\mb{\bm$A$}=\mb{\bm$L$}\mb{\bm$U$}$ are just the products of the corresponding leading principal minors of $\mb{\bm$L$}$ and $\mb{\bm$U$}$. This is because the matrix $\mb{\bm$L$}$ is lower triangular with positive diagonal elements and the matrix $\mb{\bm$U$}$ is upper triangular with positive diagonal elements. This completes our proof of the lemma.\\
\hspace*{5.5in} $\Box$

The following theorem, also due to Fiedler and Pt\'{a}k~\cite{fiedler62}, establishes some properties of the M-matrices in terms of their LU-decomposition.
\begin{thm}
\label{thm:m_matlu}
Let $\mb{\bm$A$} \in \mb{\Bb R}^{n \times n}$ be an M-matrix. Then, there is a lower triangular matrix $\mb{\bm$L$} \in \mb{\Bb R}^{n \times n}$ and an upper triangular matrix $\mb{\bm$U$} \in \mb{\Bb R}^{n \times n}$ both with positive diagonal elements. Their off-diagonal elements are non-positive. Furthermore, both $\mb{\bm$L$}^{-1}$ and $\mb{\bm$U$}^{-1}$ are non-negative which implies that $\mb{\bm$A$}^{-1}=\mb{\bm$U$}^{-1}\mb{\bm$L$}^{-1}$ is non-negative.
\end{thm}
{\bf Proof}\quad We prove this theorem by induction. The theorem holds for $n=1$. Suppose the theorem holds for all matrices of order $n-1$. Consider now the factorization from the proof of Lemma~\ref{lemma:lu}. We have
\[
\mb{\bm$A$}=\mb{\bm$L$}\mb{\bm$U$} = \left[\begin{array}{cc}
\mb{\bm$L$}_{n-1} & \mb{\bm$0$}\\
\b b^T\mb{\bm$U$}_{n-1}^{-1} & 1
\end{array} \right]\left[\begin{array}{cc}
\mb{\bm$U$}_{n-1} & \mb{\bm$L$}_{n-1}^{-1}\mb{\bm$a$}\\
\mb{\bm$0$} & a_{nn} - \b b^T\mb{\bm$A$}_{n-1}^{-1}\mb{\bm$a$}
\end{array}\right],
\]
where $\mb{\bm$L$}_{n-1}$, $\mb{\bm$U$}_{n-1}$, $\mb{\bm$L$}_{n-1}^{-1}$, $\mb{\bm$U$}_{n-1}^{-1}$ and $\mb{\bm$A$}_{n-1}^{-1}=\mb{\bm$U$}_{n-1}^{-1}\mb{\bm$L$}_{n-1}^{-1}$ satisfy the inductive hypotheses. We see that $l_{ij} \le 0$ for $i>j$ and $u_{kl} \le 0$ for $k<l$ which means that the off-diagonal elements of $\mb{\bm$L$}$ and $\mb{\bm$U$}$ are non-positive. Their diagonal elements are positive. The inverse of $\mb{\bm$L$}$ is
\[
\mb{\bm$L$}^{-1}=\left[\begin{array}{cc}
\mb{\bm$L$}_{n-1}^{-1} & \mb{\bm$0$}\\
-\b b^T\mb{\bm$L$}_{n-1}^{-1} & 1
\end{array}\right].
\]
From the inductive assumption $\mb{\bm$L$}_{n-1}^{-1}$ and $\mb{\bm$A$}_{n-1}^{-1}$ are non-negative and hence $\mb{\bm$L$}_{n-1}^{-1}$ is non-negative.\\
The inverse of $\mb{\bm$U$}$ is
\[
\mb{\bm$U$}^{-1} = \left[\begin{array}{cc}
\mb{\bm$U$}_{n-1}^{-1} & -\mb{\bm$A$}_{n-1}^{-1} \frac{\mb{\bm$a$}}{ a_{nn} - \b b^T \mb{\bm$A$}_{n-1}^{-1} \mb{\bm$a$}}\\
\mb{\bm$0$} & \left( a_{nn}-\b b^T\mb{\bm$A$}_{n-1}^{-1}\mb{\bm$a$} \right) ^{-1}
\end{array} \right].
\]
By the inductive assumptions and the fact that $a_{nn}-\b b^T\mb{\bm$A$}_{n-1}^{-1}\mb{\bm$a$}>0$, the matrix $\mb{\bm$U$}^{-1}$ is non-negative, and the proof is complete.\\
\hspace*{5.5in} $\Box$

We now consider a more general class of matrices called Z-matrices.
\begin{defn}
A real $n \times n$ matrix $\mb{\bm$A$}$ whose off-diagonal elements are non-positive is called a Z-matrix.
\end{defn}
Note that M-matrices are Z-matrices.\\
To simplify our presentation, we use the notation $\mb{\bm$B$}>0$, to mean that all the elements of $\mb{\bm$B$}$ are positive, where $\mb{\bm$B$}$ can be a matrix or a vector. Respectively, $\mb{\bm$B$} \ge 0$ means that all the elements of $\mb{\bm$B$}$ are non-negative. We are now ready to state and prove the following theorem concerning properties of Z-matrices. The theorem can be found in Fiedler and Pt\'{a}k~\cite{fiedler62}.
\begin{thm}
\label{thm:z_mat}
Let $\mb{\bm$A$} \in \mb{\Bb R}^{n \times n}$ be a Z-matrix. Then, the following conditions are equivalent:
\begin{itemize}
\item[(i)] There exists a vector $\mb{\bm$x$} \ge 0$ such that $\mb{\bm$Ax$}>0$;
\item[(ii)] There exists a vector $\mb{\bm$y$}>0$ such that $\mb{\bm$Ay$}>0$;
\item[(iii)] There exists a diagonal matrix $\mb{\bm$D$}$ with positive diagonal elements such that $\mb{\bm$ADe$}>0$, where $\mb{\bm$e$}=\left[\begin{array}{ccc}
1, & \ldots, & 1
\end{array}\right]^T \in \mb{\Bb R}^n$;
\item[(iv)] There exists a diagonal matrix $\mb{\bm$D$}$ with positive diagonal elements such that the matrix $\mb{\bm$W=AD$}$ has a dominant positive main diagonal, that is for $i=1,2,\ldots,n$, we have $w_{ii}>\sum_{\stackrel{j=1}{j \ne i}}^n w_{ij}$.
\end{itemize}
\end{thm}
{\bf Proof}\quad We prove this theorem using the following scheme:
\[
\mb{(i)} \Rightarrow \mb{(ii)} \Rightarrow \mb{(iii)} \Rightarrow \mb{(iv)} \Rightarrow \mb{(i)}.
\]
We begin with the implication (i)$ \Rightarrow $(ii). Thus suppose that there exists a vector $\mb{\bm$x$} \ge 0$ such that $\mb{\bm$Ax$} > 0$. Let $\mb{\bm$e$}=\left[\begin{array}{ccc}
1, & \ldots, & 1
\end{array}\right]^T \in \mb{\Bb R}^n$, $\epsilon > 0$, and $\mb{\bm$y$} = \mb{\bm$x$} + \epsilon \mb{\bm$e$} > 0$. Then
\[
\mb{\bm$Ay$} = \mb{\bm$A(x$} + \epsilon \mb{\bm$e$}) = \mb{\bm$Ax$} + \epsilon \mb{\bm$Ae$} > 0
\]
for a sufficiently small $\epsilon > 0$.

We next prove that (ii)$ \Rightarrow $(iii). For this let $\mb{\bm$y$}=\left[\begin{array}{ccc}
y_1, & \ldots, & y_n
\end{array}\right]^T > 0$ be such that $\mb{\bm$Ay$} > 0$. Let
\[
\mb{\bm$D$} = \left[\begin{array}{ccccc}
y_1 & 0 & 0 & \ldots & 0\\
0 & y_2 & 0 & \ldots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
0 & 0 & 0 & \ldots & y_n
\end{array} \right].
\]
It then follows that $\mb{\bm$ADe$} = \mb{\bm$Ay$} > 0$.

To prove (iii)$ \Rightarrow $(iv), let $\mb{\bm$D$} > 0$ be a diagonal matrix with positive diagonal elements such that $\mb{\bm$ADe$} > 0$. We denote
\[
\mb{\bm$AD$}=\mb{\bm$W$} = \left[\begin{array}{cccc}
w_{11} & -w_{12} & \ldots & -w_{1n}\\
-w_{21} & w_{22} & \ldots & -w_{2n}\\
\vdots & \vdots & \ddots & \vdots\\
-w_{n1} & -w_{n2} & \ldots & w_{nn}
\end{array} \right],
\]
where $w_{ij} \ge 0$ for $i \ne j$. Note that $\mb{\bm$ADe$}=\mb{\bm$We$}>0$ which is equivalent to
\[
w_{ii} > \sum_{\stackrel{j=1}{j \ne i}}^{n} w_{ij}, \qquad i=1,2,\ldots ,n.
\]
Since $w_{ij} \ge 0$, we have
\[
w_{ii} > \sum_{\stackrel{j=1}{j \ne i}}^{n} |w_{ij}|.
\]
The above means that $\mb{\bm$W$}$ is a row diagonal dominant matrix with positive main diagonal. Note $a_{ii} > 0$ since $0 < w_{ii} = a_{ii} d_i$ and $d_i > 0$.

To complete the proof of the theorem we now show that (iv)$ \Rightarrow $(i). For this let $\mb{\bm$D$}$ be a diagonal matrix with positive diagonal elements such that $\mb{\bm$W$} = \mb{\bm$AD$}$ has a positive dominant main diagonal. Then
\[
0 < \mb{\bm$We$} = \mb{\bm$ADe$} = \mb{\bm$Ax$},
\]
where
\[
\mb{\bm$x$}=\left[\begin{array}{ccc}
d_1 & \ldots & 0\\
\vdots & \ddots & \vdots\\
0 & \ldots & d_n
\end{array} \right] \left[\begin{array}{c}
1\\
\vdots\\
1
\end{array} \right]=\left[\begin{array}{c}
d_1\\
\vdots\\
d_n
\end{array}\right] > 0.
\]
Hence $(i)$ is satisfied.\\
\hspace*{5.5in} $\Box$

To proceed further, we need the following two technical results.
\begin{lemma}
\label{lemma:hadamard}
An $n \times n$ matrix $\mb{\bm$A$}$ with complex coefficients such that
\[
|a_{ii}| > \sum_{\stackrel{j=1}{j \ne i}}^{n} |a_{ij}|\quad\mb{for}\quad i=1,2,\ldots,n,
\]
is nonsingular, that is, $\det \mb{\bm$A$} \ne 0$.
\end{lemma}
{\bf Proof}\quad We prove the above lemma by contraposition. Suppose that $\det \mb{\bm$A$} = 0$. Then, there exists a nontrivial solution $\mb{\bm$x$} \ne \mb{\bm$0$}$ to $\mb{\bm$Ax=0$}$, that is for $k=1,2,\ldots,n$,
\[
\sum_{j=1}^{n} a_{kj} x_j = 0.
\]
Let
\[
x_i = \max_{1 \le k \le n} |x_k|.
\]
Rearranging $\sum_{j=1}^{n} a_{ij} x_j = 0$ yields
\[
-a_{ii}x_i=a_{i1}x_1 + \cdots + a_{i,i-1}x_{i-1} + a_{i,i+1}x_{i+1} + \cdots + a_{in}x_n.
\]
Taking the absolute values of the both sides of the above equation gives
\[
|a_{ii}||x_i|=|a_{i1}x_1 + \cdots + a_{i,i-1}x_{i-1} + a_{i,i+1}x_{i+1} + \cdots + a_{in}x_n| \le \sum_{\stackrel{j=1}{j \ne i}}^{n} |a_{ij}||x_j|.
\]
Thus
\[
|a_{ii}| \le \sum_{\stackrel{j=1}{j \ne i}}^{n} |a_{ij}| \frac{|x_j|}{|x_i|} \le  \sum_{\stackrel{j=1}{j \ne i}}^{n} |a_{ij}|
\]
which contradicts the assumption that $|a_{ii}| > \sum_{\stackrel{j=1}{j \ne i}}^{n} |a_{ij}|$.\\
\hspace*{5.5in} $\Box$

The above lemma says that if a matrix $\mb{\bm$A$}$ is row diagonal dominant then it is nonsingular. Following the arguments of the proof of the lemma, we can show that if a square, real or complex, matrix is column diagonal dominant then it is also nonsingular. We state this result in the following corollary.
\begin{corollary}
A real or complex $n \times n$ matrix $\mb{\bm$A$}$ such that 
for $j=1,2,\ldots,n$,
\[
|a_{jj}| > \sum_{\stackrel{k=1}{k \ne j}}^{n} |a_{kj}|
\]
is nonsingular. 
\end{corollary}

We now use Lemma~\ref{lemma:hadamard} to prove the following theorem due to Gershgorin.
\begin{thm}
\label{thm:gershgorin}
Given a real or complex $n \times n$ matrix $\mb{\bm$A$}$. Let $D_i$, $i=1,2,\ldots,n$, denotes the $i$-th closed disk in the complex plane centered at $a_{ii}$ whose radius is $r_i=\sum_{\stackrel{j=1}{j \ne i}}^{n} |a_{ij}|$. Then the eigenvalues of $\mb{\bm$A$}$ are located in the union $\bigcup_{i=1}^{n}D_i$. 
\end{thm}
{\bf Proof}\quad The eigenvalues of $\mb{\bm$A$}$ are the roots of the characteristic equation
\[
\det [\lambda \mb{\bm$I$} - \mb{\bm$A$}] = \det \left[\begin{array}{cccc}
\lambda - a_{11} & -a_{12} & \ldots & -a_{1n}\\
-a_{21} & \lambda - a_{22} & \ldots & -a_{2n}\\
\vdots & \vdots & \ddots & \vdots \\
-a_{n1} & -a_{n2} & \ldots & \lambda -a_{nn}
\end{array} \right]=0.
\]
By Lemma~\ref{lemma:hadamard}, the characteristic polynomial is nonzero if for all $i=1,2,\ldots,n$,
\[
|\lambda -a_{ii}| > \sum_{\stackrel{j=1}{j \ne i}}^{n} |a_{ij}|,
\]
that is, for all $i=1,2,\ldots,n$
\[
|\lambda -a_{ii}| > r_i.
\]
Thus the eigenvalues of $\mb{\bm$A$}$ must be in the set $\bigcup_{i=1}^{n} D_i$.\\
\hspace*{5.5in} $\Box$

In the following theorem, due to Fiedler and Pt\'{a}k~\cite{fiedler62}, we collect a number of equivalent conditions characterizing the M-matrices.
\begin{thm}
\label{thm:m_matsum}
Let $\mb{\bm$A$}$ be a Z-matrix. Then the following conditions are equivalent
\begin{itemize}
\item[(i)] There exists a vector $\mb{\bm$x$} \ge 0$ such that $\mb{\bm$Ax$} > 0$;
\item[(ii)] There exists a diagonal matrix $\mb{\bm$D$}$ with positive diagonal such that the matrix $\mb{\bm$W=AD$}$ has a positive main diagonal and is row diagonal dominant;
\item[(iii)] Eigenvalues of $\mb{\bm$W$}$ have positive real parts and principal minors of $\mb{\bm$W$}$ are positive;
\item[(iv)] All principal minors of $\mb{\bm$A$}$ are positive, that is, $\mb{\bm$A$}$ is an M-matrix;
\item[(v)] All leading principal minors of $\mb{\bm$A$}$ are positive;
\item[(vi)] There exists a lower triangular matrix $\mb{\bm$L$}$ and an upper triangular matrix $\mb{\bm$U$}$ both with positive diagonals such that $\mb{\bm$A=LU$}$. Furthermore, $l_{ij} \le 0$ for $i>j$, and $u_{kl} \le 0$ for $k<l$;
\item[(vii)] $\mb{\bm$A$}^{-1}$ exists and $\mb{\bm$A$}^{-1} \ge 0$.
\end{itemize}
\end{thm}
{\bf Proof}\quad The implication (i)$ \Rightarrow $(ii) follows from Theorem~\ref{thm:z_mat}. To prove (ii)$ \Rightarrow $(iii), we represent the matrix $\mb{\bm$W=AD$}$ as follows
\[
\mb{\bm$W=AD$}=\left[\begin{array}{cccc}
w_{11} & -w_{12} & \ldots & -w_{1n}\\
-w_{21} & w_{22} & \ldots & -w_{2n}\\
\vdots & \vdots & \ddots & \vdots\\
-w_{n1} & -w_{n2} & \ldots & w_{nn}
\end{array} \right],
\]
where $w_{ij} \ge 0$ for $i \ne j$, and $w_{ii} > 0$, as follows from Theorem~\ref{thm:z_mat}. Furthermore
\[
w_{ii} > \sum_{\stackrel{j=1}{j \ne i}}^{n}w_{ij} \quad \mbox{for} \quad i=1,2,\ldots,n.
\]
The above implies that each principal minor matrix of $\mb{\bm$W$}$ is also row diagonal dominant. By the Theorem of Gershgorin, the eigenvalues of $\mb{\bm$W$}$ as well as its principal minor matrices are located on the disks with centers at $w_{ii}$ and radii $r_i < w_{ii}$. Hence they must have real positive parts.

We shall now show that (iii)$ \Rightarrow $(iv) holds. First note that since the determinant of a square matrix is the product of its eigenvalues and the product of complex eigenvalues of a real matrix is positive, it follows that all principal minors of the matrix $\mb{\bm$W$}$ are positive. Next, the matrix $\mb{\bm$W$}$ has the form
\[
\mb{\bm$W=AD$}=\left[\begin{array}{cccc}
a_{11} & -a_{12} & \ldots & -a_{1n}\\
-a_{21} & a_{22} & \ldots & -a_{2n}\\
\vdots & \vdots & \ddots & \vdots\\
-a_{n1} & -a_{n2} & \ldots & a_{nn}
\end{array} \right]
\left[\begin{array}{cccc}
d_{11} & \ldots & {} & 0\\
{} & d_{22} & {} & {}\\
\vdots & {} & \ddots & \vdots\\
0 & \ldots & {} & d_{nn}
\end{array} \right],
\]
where $a_{ij} \ge 0$ for $i \ne j$, $a_{ii} > 0$, and $d_{ii} > 0$. Thus if the principal minors of $\mb{\bm$W$}$ are positive, the principal minors of $\mb{\bm$A$}$ must be positive, and vice versa. This is because
\[
\b W\left(\begin{array}{cccc}
i_1 & i_2 & \ldots & i_k\\
i_1 & i_2 & \ldots & i_k
\end{array} \right)=
\b A\left[\begin{array}{cccc}
i_1 & i_2 & \ldots & i_k\\
i_1 & i_2 & \ldots & i_k
\end{array} \right]
\left[\begin{array}{cccc}
d_{i_1 i_1} & \ldots & {} & 0\\
{} & d_{i_2 i_2} & {}& {}\\
{} & {} & \ddots & {}\\
0 & \ldots & {} & d_{i_k i_k}
\end{array}\right],
\]
where
\[
\b W\left(\begin{array}{cccc}
i_1 & i_2 & \ldots & i_k\\
i_1 & i_2 & \ldots & i_k
\end{array} \right)
\]
is a matrix formed from $i_1,i_2,\ldots,i_k$ columns and $i_1,i_2,\ldots,i_k$ rows of the matrix $\mb{\bm$W$}$.

The implication (iv)$ \Rightarrow $(v) follows from the observation that the leading principal minors are also principal minors.

The implications (v)$ \Rightarrow \mb{(vi)} \Rightarrow $(vii) follow from Theorem~\ref{thm:m_matlu}.

Finally, to complete the proof, we show that (vii)$ \Rightarrow $(i). For this note that if $\mb{\bm$A$}^{-1}$ exists and $\mb{\bm$A$}^{-1} \ge \mb{\bm$0$}$, then
\[
\mb{\bm$x=A$}^{-1}\mb{\bm$e$} \ge \mb{\bm$0$},
\]
where $\mb{\bm$e$}=\left[1,1,\ldots,1\right]^T \in \mb{\Bb R}^n$. Hence
\[
\mb{\bm$Ax=e$} > \mb{\bm$0$},
\]
which implies that (i) is satisfied.\\
\hspace*{5.5in} $\Box$


\subsection{A Comparison Theorem}
We now state and prove a theorem that allows one to reduce the study of some problems involving differential equations to the study of corresponding problems using differential inequalities instead.

We consider the differential inequality,
\[
\dot{x}(t) \le \omega(t,x(t)),\qquad x(t_0)=x_0,
\]
where $\omega : \mb{\Bb R}^2 \to \mb{\Bb R}$ is a continuous function in a domain $\Delta \subset \mb{\Bb R}^2$. We assume that in $\Delta$, the function $\omega$ satisfies the Lipschitz condition,
\[
|\omega (t,x) - \omega (t,y)| \le L|x-y|,
\]
where $L\ge 0$ is a Lipschitz constant.\\
Consider now an associated differential equation
\[
\dot{\tilde{x}}(t) = \omega(t,\tilde{x}(t)),\qquad \tilde{x}(t_0)=\tilde{x}_0.
\]
We refer to the above equation as the {\em equation of comparison}.\\
Denote the solution to the equation of comparison by
\[
\tilde{x}(t)=\tilde{x}(t;t_0,\tilde{x}_0).
\]
The following theorem, that can be found, for example, in Corduneanu~\cite[page 49]{corduneanu71}, relates the solutions of the differential inequality and the associated differential equation.
\begin{thm}
\label{thm:comparison}
If $x(t)$ is a differentiable solution to the differential inequality and $\tilde{x}(t)$ is the solution to the associated differential equation with the corresponding initial conditions such that
\[
x_0 \le \tilde{x}_0,
\]
then
\[
x(t) \le \tilde{x}(t)
\]
for all $t \ge t_0$.
\end{thm}
{\bf Proof}\quad We prove the theorem by contradiction. We suppose that $x_0 \le \tilde{x}_0$ but $x(t) \le \tilde{x}(t)$ does not hold. Thus, there exists $t_1 > t_0$ such that
\[
x(t_1) > \tilde{x}(t_1).
\]
Let
\[
y(t)=x(t)-\tilde{x}(t).
\]
Note that $y(t_0)=x(t_0)-\tilde{x}(t_0) \le 0$, and $y(t_1)=x(t_1)-\tilde{x}(t_1) > 0$. Denote by $\tau$ the largest $t \in \left[ t_0, \, t_1 \right]$ such that $y(\tau) = 0$. Such a value of $t$ exists by virtue of the continuity of $y$. Thus,
\[
y(t) > 0 \quad \mbox{for all} \quad t \in \left( \tau, t_1 \right].
\]
On the other hand, for $t \in \left[ \tau, t_1 \right]$, we have
\begin{eqnarray*}
\dot{y}(t) &=& \dot{x}(t)-\dot{\tilde{x}}(t)\\
&\le& \omega (t,x(t)) - \omega (t, \tilde{x}(t))\\
&\le& |\omega (t,x(t)) - \omega (t, \tilde{x}(t))|\\
&\le& L|x(t)-\tilde{x}(t)|.
\end{eqnarray*}
In the interval $\left[ \tau, t_1 \right]$, we have
\[
y(t)=x(t)-\tilde{x}(t) \ge 0.
\]
Hence,
\begin{eqnarray*}
\dot{y}(t) &\le& L|x(t)-\tilde{x}(t)|\\
&=& L(x(t)-\tilde{x}(t))\\
&=& Ly(t).
\end{eqnarray*}
Therefore, in the interval $\left[ \tau, t_1 \right]$, the function $y$ satisfies
\[
y(t) \le y( \tau )e^{L(t - \tau )} = 0
\]
since $y(\tau) = 0$. But this is a contradiction to the previous conclusion that $y(t) > 0$ for all $t \in \left( \tau, t_1 \right]$, which completes the proof.\\
\hspace*{5.5in} $\Box$

\subsection{Lyapunov Analysis of the Hopfield Network}
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 equations describing the Hopfield network by $C_i>0$ to get
\[
\dot{u}_i=-b_iu_i+\sum_{j=1}^{n}A_{ij}g_j(u_j) + U_i,\quad i=1,2,\ldots,n,
\]
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. If $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,
\]
for $\mb{\bm$u$}=\left[u_1,u_2,\ldots,u_n \right]^T \in \mb{\Bb R}^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 know that 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
\[
u_i=y_i+\frac{c_i}{b_i},
\]
and
\[
G_i(y_i)=g_i\left(y_i+\frac{c_i}{b_i}\right).
\]
Substituting the above into the algebraic equations in which the $c_i$ terms are present yields
\[
\dot{y}_i=-b_iy_i+\sum_{j=1}^{n}A_{ij}G_j(y_j),\quad i=1,2,\ldots,n.
\]
Observe that the right side of the above differential equations vanishes at the points $\mb{\bm$y$}$ where
\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(y_j + \frac{c_j}{b_j})=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 that the external inputs are zero.

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$. Indeed, suppose $\mb{\bm$u$}^{*} \ne \mb{\bm$0$}$ is an isolated equilibrium point of interest. 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 $\dot{u}_i=-b_iu_i+\sum_{j=1}^{n}A_{ij}g_j(u_j) + U_i$ 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)+U_i\\
&=&-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\} + U_i.
\end{eqnarray*}
Note that, by the definition of the equilibrium point, the term in the curly bracket vanishes. Thus we have
\[
\dot{x}_i =-b_ix_i+ \sum_{j=1}^{n}A_{ij}h_j(x_j+u^{*}_j)+U_i.
\]
Let
\[
G_j(x_j)=h_j(x_j+u^{*}_j)=g_j(x_j+u^{*}_j)-g_j(u^{*}_j).
\]
Then the above differential equations can be written as
\[
\dot{x}_i =-b_ix_i+ \sum_{j=1}^{n}A_{ij}G_j(x_j)+U_i, \quad i=1,2,\ldots,n,
\]
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$. Observe that the above system of differential equations has an equilibrium point at the origin of $\mb{\bm$x$}^{*}=\mb{\bm$0$}$, $\mb{\bm$U$}^{*}=\mb{\bm$0$}$. We assume that there is no other equilibrium point in a neighborhood of the origin, that is, $\mb{\bm$x^{*}=0$}$ is an isolated equilibrium state of the above system corresponding to $\mb{\bm$U$}=\mb{\bm$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_{\stackrel{j=1}{j \ne i}}^{n}A_{ij}G_j(x_j)+U_i, \quad i=1,2,\ldots,n,
\]
We shall view the terms
\[
\sum_{\stackrel{j=1}{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
\[
\dot{x}_i =-b_ix_i+A_{ii}G_i(x_i)+U_i,\quad i=1,2,\ldots,n
\]
as models of free isolated neurons. In what follows we analyze the stability of an isolated, single, neuron.

{\bf 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
\[
\dot{x}_i =-b_ix_i+A_{ii}G_i(x_i).
\]
Assume 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
\[
\sigma_{i1}x_i^2 < x_i G_i (x_i) < \sigma_{i2} x_{i2}^2,
\]
or equivalently
\[
\sigma_{i1} < \frac{G_i (x_i)}{x_i} < \sigma_{i2}.
\]
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 an $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$
\[
\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}
\vspace{3in}
\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 an $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 $\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}
\vspace{3in}
\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}
If in addition to the sufficient conditions for asymptotic stability, shown above, we also have
\[
-b_i+A_{ii} \delta_i < 0,
\]
where
\[
\delta_i = \left\{ \begin{array}{ccc}
\sigma_{i1} & \mbox{if} & A_{ii} < 0\\
\sigma_{i2} & \mbox{if} & A_{ii} > 0,
\end{array} \right.
\]
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\\
&=& \left( \frac{-b_i + A_{ii} \delta_i}{2} \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\{\frac{-b_i + A_{ii} \delta_i}{2}(t-t_0) \right\} V_i(x_i(t_0))\\
&=& \exp \left\{\frac{-b_i + A_{ii} \delta_i}{2}(t-t_0) \right\}\frac{x_i(t_0)}{2}.
\end{eqnarray*}
Hence
\[
|x_i(t)| \le \exp \left\{\frac{-b_i + A_{ii} \delta_i}{4}(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.

{\bf Stability Analysis of the Network}

We now investigate the stability of the null solution of the whole neural network. As we mentioned before, without loss of generality, we can assume that the external inputs $U_i$ are set to zero. Thus the neural network model we analyze is
\[
\dot{x}_i=-b_ix_i+\sum_{j=1}^{n}A_{ij}G_j(x_j), \quad i=1,2,\ldots,n.
\]
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}: There exists a vector $\mb{\bm$\alpha$}=\left[\alpha_1, \alpha_2,\ldots,\alpha_n \right]^T$, where $\alpha_i >0$ for $i=1,2,\ldots,n$ such that the matrix $\mb{\bm$S$}\in \mb{\Bb R}^{n \times n}$ is negative definite, where the elements of the matrix $\mb{\bm$S$}$ are defined as
\[
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 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}
The null solution of the neural network model
\[
\dot{x}_i=-b_ix_i+\sum_{j=1}^{n}A_{ij}G_j(x_j),\quad i=1,2,\ldots,n,
\]
is 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$}^T \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 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$}|^T \mb{\bm$R$} |\mb{\bm$x$}|,
\end{eqnarray*}
where $|\mb{\bm$x$}|=\left[|x_1|,|x_2|,\ldots,|x_n|\right]^T$, 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$}|^T \mb{\bm$R$} |\mb{\bm$x$}| &=& |\mb{\bm$x$}|^T \left( \frac{\mb{\bm$R$}+\mb{\bm$R$}^T}{2} \right)|\mb{\bm$x$}|\\
&=& |\mb{\bm$x$}|^T \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}. The last inequality follows from Rayleigh's Principle. By 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.
\]
Since $V$ is a quadratic Lyapunov function, the null solution is exponentially stable.\\
\hspace*{5.5in} $\Box$

Note that 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,
\[
-b_i+a_{ii} < 0, \quad i=1,2,\ldots,n.
\]
However this condition is not sufficient. 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
\[
x_iA_{ii}G_i(x_i) \le a_{ii} |x_i|^2.
\]
Applying the sector condition $\sigma_{i1}x_i^2 < x_i G_i (x_i) < \sigma_{i2} x_{i2}^2$ to the above inequalities yields
\[
a_{ii}=A_{ii} \delta_i,
\]
where
\[
\delta_i=\left\{\begin{array}{ccc}
\sigma_{i2} & \mb{if} & A_{ii} > 0\\
\sigma_{i1} & \mb{if} & A_{ii} < 0
\end{array} \right.
\]
Therefore for $i=1,2,\ldots,n$, we have
\[
-b_i + A_{ii} \delta_i < 0
\]
since $-b_i+a_{ii} < 0$ and $a_{ii}=A_{ii} \delta_i$. Thus {\bf A2} can only be satisfied if the null solution of each individual neuron is exponentially stable.\\
We now present another sufficiency condition for the exponential stability of the Hopfield neural network model. Before we state the condition we need to perform some algebraic manipulations. We first combine the assumption {\bf A2} with the sector condition $\sigma_{i1}x_i^2 < x_i G_i (x_i) < \sigma_{i2} x_{i2}^2$. It follows from the sector condition that for $x_j \ne 0$
\[
\frac{G_j(x_j)}{x_j} < \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}||x_i \frac{G_j(x_j)}{x_j} x_j|\\
&\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 $\dot{x}_i=-b_ix_i+\sum_{j=1}^{n}A_{ij}G_j(x_j)$, $i=1,2,\ldots,n$ 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[\alpha_1,\alpha_2,\ldots,\alpha_n \right]^T \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_{j2} \right) & \mb{if} & i \ne j.
\end{array} \right.
\]
We are now ready to state and prove an alternative sufficiency condition for the Hopfield neural network to be exponentially stable. This theorem can be found in Michel et al.\ ~\cite{michel89}.
\begin{thm}
The null solution of the neural network model
\[ 
\dot{x}_i=-b_ix_i+\sum_{j=1}^{n}A_{ij}G_j(x_j),\quad i=1,2,\ldots,n,
\]
is exponentially stable if the assumptions {\bf A3} and {\bf A4} are satisfied.
\end{thm}
{\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 A1} and {\bf A2} we evaluate the time derivative of $V$ along the trajectories of the neural network 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_{\stackrel{j=1}{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_{\stackrel{j=1}{j \ne i}}^{n} |A_{ij}||x_i||x_j| \, \sigma_{j2} \right)\\
&=& \frac{1}{2}|\mb{\bm$x$}|^T\mb{\bm$S$}|\mb{\bm$x$}|\\
&<& 0 \quad \mb{for} \quad \mb{\bm$x$} \ne \mb{\bm$0$}.
\end{eqnarray*}
Thus the null solution of the neural network model is exponentially stable since $V$ is a quadratic Lyapunov function.\\
\hspace*{5.5in} $\Box$

\subsection{M-Matrices in Stability Analysis}
An alternative approach to the stability analysis of the Hopfield neural network is by using properties of M-matrices given in Theorem~\ref{thm:m_matsum}. This approach was proposed by Michel at al.\ ~\cite{michel89}. To proceed we need the following assumption:\\
{\bf A5:} The matrix $\mb{\bm$D$} \in \mb{\Bb R}^{n \times n}$ with its elements defined as
\[
d_{ij}=\left\{\begin{array}{ccc}
b_i - \delta_i A_{ii} & \mb{if} & i=j\\
-|A_{ij}| \sigma _{j2} & \mb{if} & i \ne j
\end{array} \right.
\]
is an M-matrix.

We now prove a sufficiency condition, using certain properties of M-matrices, for an equilibrium state of the Hopfield neural network model to be asymptotically stable. Without loss of generality we state and prove the condition for $\b x=\b 0$.
\begin{thm}
\label{thm:stab_m_mat}
The null solution of the neural network model 
\[
\dot{x}_i=-b_ix_i+\sum_{j=1}^{n}A_{ij}G_j(x_j),\quad i=1,2,\ldots,n,
\]
is asymptotically stable if the assumptions {\bf A3} and {\bf A5} are satisfied.
\end{thm}
{\bf Proof}\quad Consider the following positive definite Lyapunov function candidate
\[
V(\mb{\bm$x$}) = \sum_{i=1}^{n} \alpha_i |x_i|,
\]
where $\alpha_i > 0$, for $i=1,2,\ldots,n$, are weighting factors to be specified later in the proof. The time derivative of $V(\mb{\bm$x$}(t))$ along the solutions of the neural network when $x_i \ne 0$ evaluates to
\[
\dot{V}(\mb{\bm$x$}(t))=\sum_{i=1}^{n} \alpha_i (\mb{sign}\,x_i) \dot{x}_i,
\]
where
\[
\mb{sign}\,x_i = \frac{|x_i|}{x_i} = \frac{x_i}{|x_i|} = \left\{ \begin{array}{ccc}
1 & \mb{if} & x_i > 0\\
-1 & \mb{if} & x_i < 0.
\end{array} \right.
\]
Hence
\begin{eqnarray*}
\dot{V} &=& \sum_{i=1}^{n} \alpha_i (\mb{sign}\, x_i) \left\{-b_ix_i+A_{ii}G_i(x_i) + \sum_{\stackrel{j=1}{j \ne i}}^{n}A_{ij}G_j(x_j) \right\}\\
&=& \sum_{i=1}^{n} \alpha_i \left\{-b_i|x_i|+\frac{x_iA_{ii}G_i(x_i)}{|x_i|} + \sum_{\stackrel{j=1}{j \ne i}}^{n}\frac{x_iA_{ij}G_j(x_j)}{|x_i|} \right\}\\
&\le& \sum_{i=1}^{n} \alpha_i \left\{(-b_i +\delta_i A_{ii})|x_i| + \sum_{\stackrel{j=1}{j \ne i}}^{n}|A_{ij}| \sigma_{j2}|x_j|  \right\}\\
&=& -\mb{\bm$\alpha $}^T \mb{\bm$D$}|\mb{\bm$x$}|,
\end{eqnarray*}
where $|\mb{\bm$x$}|=\left[|x_1|,|x_2|,\ldots,|x_n| \right]^T$, and the matrix $\mb{\bm$D$}$ is defined in {\bf A5}. Since $\mb{\bm$D$}$ is an M-matrix, by Theorem~\ref{thm:m_matsum}, $\mb{\bm$D$}^{-1}$ exists and $\mb{\bm$D$}^{-1} \ge \mb{\bm$0$}$. Let
\[
\mb{\bm$y$}^T=\mb{\bm$\alpha$}^T \mb{\bm$D$}.
\]
Then
\[
\mb{\bm$\alpha$}=\left(\mb{\bm$D$}^{-1}\right)^T\mb{\bm$y$}.
\]
If we choose $\mb{\bm$y$}>\mb{\bm$0$}$ then $\mb{\bm$\alpha$}>\mb{\bm$0$}$ and
\[
\dot{V}(\mb{\bm$x$}(t)) \le -\mb{\bm$y$}^T |\mb{\bm$x$}|
\]
is negative definite in some neighborhood of $\mb{\bm$x$}=\mb{\bm$0$}$.\\
The case when one or more $x_i=0$ on an open time interval can be dealt with in a similar fashion as the previous case when all $x_i \ne 0$. In particular, if $x_k=0$ then we consider the following Lyapunov function candidate
\[
V_{(k)}(\mb{\bm$x$})=\sum_{\stackrel{i=1}{i \ne k}}^{n} \alpha_i |x_i|.
\]
The time derivative of $V_{(k)}(\mb{\bm$x$}(t))$ along the solutions of the neural network evaluates to
\[
\dot{V}(\mb{\bm$x$}(t)) \le -\mb{\bm$\alpha$}_{(k)}^T \mb{\bm$D$}_{(k)}| \mb{\bm$x$}_{(k)}|,
\]
where 
\begin{eqnarray*}
\mb{\bm$\alpha$}_{(k)} &=& \left[\alpha_1,\ldots,\alpha_{k-1},\alpha_{k+1},\ldots,\alpha_n \right]^T \in \mb{\Bb R}^{n-1},\\
|\mb{\bm$x$}_{(k)}| &=& \left[ |x_1|,\ldots,|x_{k-1}|,|x_{k+1}|,\ldots,|x_n| \right]^T,
\end{eqnarray*}
and $\mb{\bm$D$}_{(k)} \in \mb{\Bb R}^{(n-1) \times (n-1)}$ is obtained by deleting the $k$th row and the $k$th column from $\mb{\bm$D$}$.\\
Notice that $\mb{\bm$D$}_{(k)}$ is an M-matrix. Hence we can find a vector $\mb{\bm$y$}_{(k)} \in \mb{\Bb R}^{n-1}$, $\mb{\bm$y$}_{(k)} > \mb{\bm$0$}$ that results in
\[
\dot{V}_{(k)}(\mb{\bm$x$}(t)) \le -\mb{\bm$y$}_{(k)}^T |\mb{\bm$x$}_{(k)}| <0
\]
in some neighborhood of the origin of $\mb{\Bb R}^n$.\\
The case when some $x_i=0$ for time instances of measure zero does not pose any problems.\\
In summary, the null solution of the neural network is asymptotically stable.\\
\hspace*{5.5in} $\Box$

We give yet another sufficiency condition, in terms of M-matrices, for asymptotic stability of an equilibrium state of the neural network. This theorem is also due to Michel et al.\ ~\cite{michel89}. This time instead of using the assumption {\bf A5}, we use the following assumption:\\
{\bf A6:} The matrix $\mb{\bm$D$} \in \mb{\Bb R}^{n \times n}$ with its elements defined as
\[
d_{ij}=\left\{\begin{array}{ccc}
\frac{b_i}{\sigma_{i2}} - A_{ii} & \mb{if} & i=j\\
-|A_{ij}| & \mb{if} & i \ne j
\end{array}\right.
\]
is an M-matrix.
\begin{thm}
\label{thm:a3_a6}
The null solution of the neural network model 
\[
\dot{x}_i=-b_ix_i+\sum_{j=1}^{n}A_{ij}G_j(x_j),\quad i=1,2,\ldots,n,
\]
is asymptotically stable if the assumptions {\bf A3} and {\bf A6} are satisfied.
\end{thm}
{\bf Proof}\quad We consider a similar Lyapunov function candidate as in the previous Theorem, that is,
\[
V(\mb{\bm$x$}) = \sum_{i=1}^{n} \alpha_i |x_i|,
\]
where $\alpha_i > 0$, $i=1,2,\ldots,n$, are to be determined. The time derivative of $V$ along the solutions of the neural network, when $x_i \ne 0$, $i=1,2,\ldots,n$, evaluates to
\begin{eqnarray*}
\dot{V}(\mb{\bm$x$}(t)) &=& \sum_{i=1}^{n} \alpha_i (\mb{sign}\, x_i) \left\{-b_ix_i+A_{ii}G_i(x_i) + \sum_{\stackrel{j=1}{j \ne i}}^{n}A_{ij}G_j(x_j) \right\}\\
&=& \sum_{i=1}^{n} \alpha_i \left\{-b_i|x_i|+\frac{A_{ii}G_i(x_i)|x_i|}{x_i} + \sum_{\stackrel{j=1}{j \ne i}}^{n}\frac{A_{ij}G_j(x_j)x_i}{|x_i|}\right\}\\
&\le& \sum_{i=1}^{n} \alpha_i \left\{\left(-b_i \frac{x_i}{G_i(x_i)}+A_{ii} \right) \frac{G_i(x_i)}{x_i}|x_i|+ \sum_{\stackrel{j=1}{j \ne i}}^{n}|A_{ij}||G_j(x_j)|  \right\}.
\end{eqnarray*}
From the sector condition $\sigma_{i1} < \frac{G_i (x_i)}{x_i} < \sigma_{i2}$, we obtain
\[
\frac{1}{\sigma_{i1}} > \frac{x_i}{G_i(x_i)} > \frac{1}{\sigma_{i2}}.
\]
Furthermore, since $\sigma_{i1} > 0$  in the sector condition, and we assumed $x_i \ne 0$, for $i=1,2,\ldots,n$, we can write
\[
|G_j(x_j)|=\frac{G_j(x_j)}{x_j}|x_j|.
\]
Taking the above into account yields
\begin{eqnarray*}
\dot{V}(\mb{\bm$x$}(t)) &\le& \sum_{i=1}^{n} \alpha_i \left\{ \left( -\frac{b_i}{\sigma_{i2}} + A_{ii} \right) \frac{G_i(x_i)}{x_i} |x_i| + \sum_{\stackrel{j=1}{j \ne i}}^{n}|A_{ij}|\frac{G_j(x_j)}{x_j}|x_j|\right\}\\
&=& -\left[\alpha_1,\alpha_2,\ldots,\alpha_n \right]^T \b D \left[\begin{array}{c}
\frac{G_1(x_1)}{x_1}|x_1|\\
\frac{G_2(x_2)}{x_2}|x_2|\\
\vdots\\
\frac{G_n(x_n)}{x_n}|x_n|
\end{array}\right]\\
&=& -\mb{\bm$\alpha$}^T \mb{\bm$D$} \mb{\bm$w$},
\end{eqnarray*}
where
\[
\b D=\left[\begin{array}{cccc}
\frac{b_1}{\sigma_{12}}-A_{11} & -|A_{12}| & \ldots & -|A_{1n}|\\
-|A_{21}| & \frac{b_2}{\sigma_{22}}-A_{22} & \ldots & -|A_{2n}|\\
\vdots & \vdots & \ddots & \vdots \\
-|A_{n1}| & -|A_{n2}| & \ldots & \frac{b_n}{\sigma_{n2}}-A_{nn}\\
\end{array}\right].
\]
Note that
\[
\mb{\bm$w$}=\left[\frac{G_1(x_1)}{x_1}|x_1|,\frac{G_2(x_2)}{x_2}|x_2|,\ldots,\frac{G_n(x_n)}{x_n}|x_n| \right]^T > \mb{\bm$0$},
\]
and the matrix $\mb{\bm$D$}$ is an M-matrix. Hence there exists an $\mb{\bm$\alpha$} > \mb{\bm$0$}$ such that $\mb{\bm$y$} = \mb{\bm$\alpha$}^T \mb{\bm$D$}> \mb{\bm$0$}$ and thus
\[
\dot{V}(\mb{\bm$x$}(t)) \le \mb{\bm$y$}^T\mb{\bm$w$}<0.
\]
If some $x_i=0$ then we proceed similarly as in the second part of the proof of the previous Theorem.

In summary, the null solution of the neural network is asymptotically stable.\\
\hspace*{5.5in} $\Box$

\subsection{Estimating Trajectory Bounds}
Our next task is to estimate trajectories's bounds near an equilibrium point. The bounds obtained can serve as parameters characterizing qualitative behavior of the network. In the subsequent development we do not require the external inputs $U_i$, $i=1,2,\ldots,n$ to be set to zero. Thus the network we analyze is described by
\[
\dot{x}_i =-b_i x_i+ \sum_{j=1}^{n}A_{ij}G_j(x_j)+U_i, \quad i=1,2,\ldots,n.
\]
We use the assumption {\bf A6} that states that the matrix $\mb{\bm$D$}$ is an M-matrix. Since $\mb{\bm$D$}$ is an M-matrix there exist positive constants $\lambda_i$, $i=1,2,\ldots, n$, and $\varepsilon > 0$ such that
\[
\left( \frac{b_i}{\sigma_{i2}} - A_{ii} \right) - \sum_{\stackrel{j=1}{j \ne i}}^{n} \left( \frac{\lambda_j}{\lambda_i} \right) |A_{ij}| \ge \varepsilon .
\]
We also use another assumption:\\
{\bf A7:} For the neural network model $\dot{x}_i =-b_ix_i+ \sum_{j=1}^{n}A_{ij}G_j(x_j)+U_i$, $i=1,2,\ldots,n$, there exists a positive constant $k$ such that for $t \ge 0$
\[
\sum_{i=1}^{n} \lambda_i |U_i(t)| \le k.
\]
We use the following notation. Let
\[
c = \varepsilon \min_{1 \le i \le n} \lambda_i \sigma_{i1},
\]
and
\[
\mb{\bm$\Phi$} (t;t_0,\mb{\bm$x$}_0)=\left[ \Phi_1 (t;t_0,\mb{\bm$x$}_0),\Phi_2 (t;t_0,\mb{\bm$x$}_0),\ldots,\Phi_n (t;t_0,\mb{\bm$x$}_0) \right]^T
\]
denotes the solution to the neural network model subject to the initial condition $\mb{\bm$x$}_0$ at $t_0 \ge 0$.
\begin{thm}
If for $i=1,2,\ldots,n$ we have
\[
\left( \frac{b_i}{\sigma_{i2}} - A_{ii} \right) - \sum_{\stackrel{j=1}{j \ne i}}^{n} \left( \frac{\lambda_j}{\lambda_i}\right) |A_{ij}| \ge \varepsilon
\]
and {\bf A7} are satisfied then
\begin{eqnarray*}
\|\mb{\bm$\Phi$} (t;t_0,\mb{\bm$x$}_0) \| &\stackrel{\rm def}{=}& \sum_{i=1}^{n} \lambda_i |\Phi_i (t;t_0,\mb{\bm$x$}_0)|\\
&\le& \left( \gamma -\frac{k}{c} \right) \exp (-c(t-t_0)) + \frac{k}{c},\quad t \ge t_0 \ge 0,
\end{eqnarray*}
provided
\[
\gamma > \frac{k}{c},
\]
and
\begin{eqnarray*}
\| \mb{\bm$x$}_0 \| &\stackrel{\rm def}{=}& \sum_{i=1}^{n} \lambda_i |x_{i0}|\\
&\le& \gamma.
\end{eqnarray*}
\end{thm}
{\bf Proof}\quad Consider the following positive definite Lyapunov function candidate
\[
V(\mb{\bm$x$})=\sum_{i=1}^{n} \lambda_i |x_i|,
\]
where $x_i \ne 0$, $i=1,2,\ldots, n$. We evaluate the time derivative of $V(\mb{\bm$x$}(t))$ along the solutions of the neural network model in a similar fashion as in the proof of the previous theorem to get
\[
\dot{V}(\mb{\bm$x$}(t)) \le -\mb{\bm$\lambda$}^T \mb{\bm$D$}\mb{\bm$w$} + \sum_{i=1}^{n}\lambda_i |U_i(t)|.
\]
Note that $\mb{\bm$w$} > \mb{\bm$0$}$. It follows from the assumptions that for $i=1,2,\ldots,n$,
\[
\mb{\bm$\lambda$}^T \left[\begin{array}{c}
-|A_{1i}|\\
\vdots\\
\frac{b_i}{\sigma_{i2}} - A_{ii}\\
\vdots\\
-|A_{ni}|
\end{array}\right] \frac{G_i(x_i)}{x_i}|x_i| \ge \lambda_i \, \varepsilon \,\sigma_{i1} |x_i|.
\]
Hence, using the notation introduced before we stated the theorem, we obtain
\[
\dot{V} \le -cV + u,
\]
where $u=\sum_{i=1}^{n}\lambda_i |U_i(t)|$. Applying now Theorem~\ref{thm:comparison} yields
\begin{eqnarray*}
\| \mb{\bm$\Phi$}(t;t_0,\mb{\bm$x$}_0) \| &\stackrel{\rm def}{=}& \sum_{i=1}^{n} \lambda_i |x_i|\\
&=& \sum_{i=1}^{n} \lambda_i |\Phi_i (t;t_0,\mb{\bm$x$}_0)|\\
&\le& \left(\gamma - \frac{k}{c}\right) \exp(-c(t-t_0) +\frac{k}{c}
\end{eqnarray*}
provided
\[
\| \mb{\bm$\Phi$}(t_0;t_0,\mb{\bm$x$}_0) \| = \sum_{i=1}^{n} \lambda_i |x_{i0}| \le \gamma,
\]
and
\[
\gamma > \frac{k}{c}.
\]
The case when some $x_i=0$ can be dealt with in similar fashion as in the proof of Theorem~\ref{thm:stab_m_mat}.\\
\hspace*{5.5in} $\Box$

\subsection{Instability of the Hopfield Network}
Of course not all of the equilibrium states of the Hopfield neural network model have to be stable. We mow present a sufficiency condition for an equilibrium state to be unstable. In the forthcoming discussion we use the following assumption:\\
{\bf A8:} The matrix $\mb{\bm$D$} \in \mb{\Bb R}^{n \times n}$ with its elements defined as
\[
d_{ij}=\left\{\begin{array}{ccc}
\sigma_i & \mb{if} & i=j\\
-|A_{ij}| & \mb{if} & i \ne j
\end{array}\right.
\]
is an M-matrix, where
\begin{eqnarray*}
\sigma_i &=& \frac{b_i}{\sigma_{i2}} - A_{ii} \quad \mb{if} \quad i \in F_s,\\
\sigma_i &=& \frac{b_i}{\sigma_{i1}} - A_{ii} \quad \mb{if} \quad i \in F_u,\\
F &=& F_s \bigcup F_u,\\
F &=& \{1,2,\ldots,n\}.
\end{eqnarray*}
By Theorem~\ref{thm:m_matsum}, the matrix $\mb{\bm$D$}$ is an M-matrix if all its principal minors are positive. Thus for {\bf A8} to hold we have to have
\[
\begin{array}{ccc}
\frac{b_i}{\sigma_{i2}} - A_{ii} > 0 & \mb{if} & i \in F_s,\\
-\frac{b_i}{\sigma_{i1}} + A_{ii} > 0 & \mb{if} & i \in F_u,
\end{array}
\]
that is, all its first order principal minors must be positive. Observe that if $i \in F_u$ then $b_i/ \sigma_{i1} < A_{ii}$. Since $b_i/ \sigma_{i1} > 0$, hence $A_{ii} > 0$ for $i \in F_u$. One can verify that an individual neuron for which $A_{ii} > 0$ and $b_i/ \sigma_{i1} < A_{ii}$ is unstable. Indeed, let us consider the following Lyapunov function candidate
\[
v_i(x_i)=|x_i|
\]
for the single neuron model $\dot{x}_i =-b_ix_i+A_{ii}G_i(x_i)$. The time derivative of $v_i$ evaluated along the trajectories of the single neuron model, when $x_i \ne 0$, is
\begin{eqnarray*}
\dot{v}_i &=& \frac{x_i}{|x_i|}\dot{x}_i\\
&=& \frac{x_i}{|x_i|}\left( -b_i x_i+A_{ii}G_i(x_i) \right)\\
&\ge& -b_i|x_i| + \frac{1}{|x_i|}\sigma_{i1} A_{ii} x_i^2\\
&=& \sigma_{i1} \left( -\frac{b_i}{\sigma_{i1}} + A_{ii} \right)|x_i|\\
&>& 0.
\end{eqnarray*}
Therefore using the Test for Instability we conclude that the single neuron with the index $i \in F_u$ is unstable. On the other hand, using the same Lyapunov function candidate as above, we can show that the single neurons with the indices $i \in F_s$ are asymptotically stable. We now generalize the above results to the whole network. The following theorem is due to Michel et al.\ ~\cite{michel89}.
\begin{thm}
The null solution of the neural network model 
\[
\dot{x}_i=-b_ix_i+\sum_{j=1}^{n}A_{ij}G_j(x_j),i=1,2,\ldots,n,
\]
is unstable if {\bf A3} and {\bf A8} are satisfied, where $F_u \ne \emptyset$, ($\emptyset$ denotes an empty set)
\end{thm}
{\bf Proof}\quad We use the following Lyapunov function candidate
\[
V(\mb{\bm$x$})=\sum_{i \in F_u} \alpha_i \left( -|x_i| \right)+\sum_{i \in F_s} \alpha_i |x_i|,
\]
where $\alpha_i > 0$, $i=1,2,\ldots,n$, are to be determined. Note that in every neighborhood of $\mb{\bm$x$}=\mb{\bm$0$}$ there is at least one $\mb{\bm$x$}$ such that $V(\mb{\bm$x$}) < 0$. We now evaluate the time derivative of $V$ for $x_i \ne 0$, $i=1,2,\ldots,n$, along the solutions of the neural network model to get
\begin{eqnarray*}
\dot{V}(\mb{\bm$x$}(t)) &=& -\sum_{i \in F_u} \alpha_i \mb{sign}(x_i)\dot{x}_i +\sum_{i \in F_s} \alpha_i \mb{sign}(x_i)\dot{x}_i\\
&=& \sum_{i \in F_u} \alpha_i \left( b_i|x_i| - \frac{A_{ii} G_i(x_i)|x_i|}{x_i} - \sum_{\stackrel{j=1}{j \ne i}}^{n} \frac{A_{ij}G_j(x_j)|x_i|}{x_i} \right)\\
&+& \sum_{i \in F_s} \alpha_i \left(- b_i|x_i| + \frac{A_{ii} G_i(x_i)|x_i|}{x_i} + \sum_{\stackrel{j=1}{j \ne i}}^{n} \frac{A_{ij}G_j(x_j)|x_i|}{x_i} \right)\\
&\le& \sum_{i \in F_u} \alpha_i \left( \left(\frac{b_i}{\sigma_{i1}}-A_{ii} \right) \frac{G_i(x_i)}{x_i}|x_i| + \sum_{\stackrel{j=1}{j \ne i}}^{n}|A_{ij}|\frac{G_j(x_j)}{x_i}|x_i| \right)\\
&-& \sum_{i \in F_s} \alpha_i \left( \left(\frac{b_i}{\sigma_{i2}}-A_{ii} \right) \frac{G_i(x_i)}{x_i}|x_i| + \sum_{\stackrel{j=1}{j \ne i}}^{n}|A_{ij}|\frac{G_j(x_j)}{x_i}|x_i| \right)\\
&=& -\mb{\bm$\alpha$}^T \mb{\bm$D$}\mb{\bm$w$}\\
&<&0
\end{eqnarray*}
for some $\mb{\bm$\alpha$}>\mb{\bm$0$}$. We can find an appropriate $\mb{\bm$\alpha$}$ using the method of the proof of Theorem~\ref{thm:a3_a6}. Using the Test for Instability, we conclude that the null solution of the neural network is unstable.\\
\hspace*{5.5in} $\Box$

\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 modelling 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 Reserach Group~\cite{rumelhart86}, and Anderson~\cite{anderson95}. Impatient readers wishing to get a feel of artificial neural networks may consult tutorials by Lippmann~\cite{lippmann87}, and 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. This 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.
\section{Exercises}
\begin{exercise}\rm
Show
\end{exercise}