\chapter{Neural Networks}
This chapter is concerned with applying the Lyapunov's method to the stability analysis of a class of dynamical systems called {\em artificial neural networks\/}. 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 analayze two different classes of neural networks. We begin with the so called radial-basis function~(RBF) neural networks and their applications to identification of dynamical systems.

\section{Radial-Basis Function Networks}
\label{sec:rbfnets}
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}{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}
In this example, $\phi_{ij}=\exp\left(-\frac{(x^{(i)}-x^{(j)})^2}{2\sigma_i^2}\right)$. 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 shown 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$. 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{4in}
\hskip.75in \special{eps:c:/brain/ad_conv.eps y=4in}
\caption{Constructing the cost function for a 2-bit A/D converter; plots of the cost function and its components.}
\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{4in}
\hskip.75in\special{eps:c:/brain/adphase.eps y=4in}
\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 Lyapunov's second method to analyze stability in the sense of Lyapunov of the Hopfield neural network. The neural network is viewed as a large interconnected dynamical system, that is, as an interconnection of many neurons. We first analyze stability properties of a single neuron, and then proceed with the stability analysis of the whole network in terms of the properties of the single neurons. Interconnection information is used in conjunction with properties of the individual neurons to obtain sufficient conditions for asymptotic stability of the network.


\subsection{Hopfield Network Model Analysis}
Recall the equations describing the Hopfield neural net. They are given by~(\ref{eq:stahoal}), that is,
\[
C_i\frac{du_i}{dt}+\frac{u_i}{r_i}=\sum_{j=1}^{n}T_{ij}x_j + I_i,\quad i=1,2,\ldots,n,
\]
where $x_i=g_i(u_i)$. In our further analysis, we assume that the activation functions $g_i:\mb{\Bb R} \to (-1,1)$, $i=1,2,\ldots,n$, are continuously differentiable and strictly monotonically increasing. We also assume that
\[
u_ig_i(u_i)>0\quad \mbox{for} \quad u_i \ne 0,
\]
and 
\[
g_i(0)=0.
\]
To proceed further, we divide both sides of the equations describing the Hopfield network by $C_i>0$ to get
\begin{equation}
\dot{u}_i=-b_iu_i+\sum_{j=1}^{n}A_{ij}g_j(u_j) + U_i,\quad i=1,2,\ldots,n,
\label{eq:newfoofhof}
\end{equation}
where $b_i=\frac{1}{r_iC_i}$, $A_{ij}=\frac{T_{ij}}{C_i}$, and $U_i=\frac{I_i}{C_i}$.\\
We are interested in the behavior of solutions of the above model near its equilibrium points. 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 can always transform the algebraic equations with non-zero constant external inputs into the algebraic equations where the external inputs are set to zero. Indeed, let for $j=1,2,\ldots,n$,
\[
u_j=y_j+\frac{c_j}{b_j},
\]
and
\[
G_j(y_j)=g_j\left(y_j+\frac{c_j}{b_j}\right).
\]
Substituting the above into~(\ref{eq:newfoofhof}) and assuming that $U_i=c_i$ yields
\begin{equation}
\dot{y}_i=-b_iy_i+\sum_{j=1}^{n}A_{ij}G_j(y_j),\quad i=1,2,\ldots,n.
\label{eq:afttranho}
\end{equation}
The equilibrium points of the above system of differential equations are obtained by solving the following algebraic equations,
\begin{eqnarray*}
-b_iy_i+\sum_{j=1}^{n}A_{ij}G_j(y_j)
&=&-b_iy_i+\sum_{j=1}^{n}A_{ij}g_j(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 and that the network is modeled by the equations,
\begin{equation}
\dot{u}_i=-b_iu_i+\sum_{j=1}^{n}A_{ij}g_j(u_j),\quad i=1,2,\ldots,n,
\label{eq:newfoofhouzer}
\end{equation}


Throughout, we assume that an isolated equilibrium point of interest is located at the origin of the state-space $\mb{\Bb R}^n$. If this is not the case, then by an appropriate change of variables we can arrange for the equilibrium point of interest to be translated to the origin of $\mb{\Bb R}^n$. Indeed, suppose $\mb{\bm$u$}^{*} \ne \mb{\bm$0$}$ is an isolated equilibrium point of interest of~(\ref{eq:newfoofhof}). Let for $i=1,2,\ldots,n$
\begin{eqnarray*}
x_i &=& u_i-u^{*}_i,\\
h_i(u_i) &=& h_i(x_i + u^{*}_i)=g_i(u_i)-g_i(u^{*}_i).
\end{eqnarray*}
Substituting the above into~(\ref{eq:newfoofhouzer}) yields
\begin{eqnarray*}
\dot{x}_i &=&-b_ix_i - b_iu_i^{*} + \sum_{j=1}^{n}A_{ij}\left(h_j(x_j+u^{*}_j)+g_j(u^{*}_j) \right)\\
&=&-b_ix_i+\sum_{j=1}^{n}A_{ij}h_j(x_j+u^{*}_j)+\left\{-b_iu^{*}_i+\sum_{j=1}^{n}A_{ij}g_j(u^{*}_j)\right\}.
\end{eqnarray*}
Note that, by the definition of the equilibrium point, the expression in the curly brackets vanishes. Thus we have
\[
\dot{x}_i =-b_ix_i+ \sum_{j=1}^{n}A_{ij}h_j(x_j+u^{*}_j).
\]
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
\begin{equation}
\dot{x}_i =-b_ix_i+ \sum_{j=1}^{n}A_{ij}G_j(x_j)+U_i, \quad i=1,2,\ldots,n,
\label{eq:xcoorofhosin}
\end{equation}
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), \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
\begin{equation}
\dot{x}_i =-b_ix_i+A_{ii}G_i(x_i),\quad i=1,2,\ldots,n
\label{eq:ensosinneu}
\end{equation}
as models of free isolated neurons. In what follows we analyze the stability of an isolated single neuron.

\subsection{Single Neuron Stability Analysis}
A single neuron in the Hopfield network, when the external input is set to zero, is modeled by a first order differential equation
\begin{equation}
\fbox{$\displaystyle \dot{x}_i =-b_ix_i+A_{ii}G_i(x_i) $}
\label{eq:sinneueq}
\end{equation}
Consider now, as in Michel et al.~\cite{michel89}, the following positive definite Lyapunov function candidate
\[
V_i(x_i)=\frac{1}{2} x_{i}^2.
\]
Its time derivative evaluated along the trajectories of the single neuron model is
\begin{eqnarray*}
\dot{V}_i (x_i(t)) &=& x_i \dot{x}_i\\
&=& x_i \left( -b_ix_i+A_{ii}G_i(x_i) \right).
\end{eqnarray*}
If there exists 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$ and $x_i \in (-r_i,r_i)$,
\[
\dot{V}_i(x_i(t)) < 0
\]
and $\dot{V}_i(0)=0$. Hence the null solution $x_i=0$ is uniformly asymptotically stable. Figure~\ref{fig:st_nen} illustrates this case.
\begin{figure}
\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 $x_i \in (-r_i,r_i)$, and $\dot{V}_i(0)=0$. Thus the null solution $x_i=0$ is unstable. Figure~\ref{fig:ust_nen} illustrates the case when the null solution of a single neuron is unstable.
\begin{figure}
\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}
Assume now that for some $r_i>0$ there exist constants $\sigma _{i1} > 0$ and $\sigma _{i2} > 0$ such that for all $x_i \in (-r_i,r_i)$, $x_i \ne 0$, we have
\[
\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}.
\]
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.

\subsection{Stability Analysis of the Network}
Lyapunov's second method in conjunction with interconnection information will now be used in the stability analysis of the Hopfield neural model viewed as an interconnected dynamical system. Applications of Lyapunov's second method to a general class of interconnected systems were investigated by Bailey~\cite{bailey66}. Here, we use Lyapunov's second method to investigate the stability of the null solution of the neural network described by
\begin{equation}
\fbox{$\displaystyle \dot{x}_i=-b_ix_i+\sum_{j=1}^{n}A_{ij}G_j(x_j), \quad i=1,2,\ldots,n $}
\label{eq:hopasintersy}
\end{equation}
We impose the following assumptions on the interconnections:\\
{\bf A1}: There exist positive constants $r_i$, $r_j$, and $a_{ij}$ such that for $|x_i| < r_i$ and $|x_j| < r_j$, $i,j=1,2,\ldots,n$, we have
\[
x_iA_{ij}G_j(x_j) \le |x_i|a_{ij}|x_j|.
\]
{\bf A2}: Let $\mb{\bm$S$}$ be an $n\times n$ real matrix whose elements are
\[
s_{ij}=\left\{\begin{array}{ccc}
\alpha_i \left(-b_i + a_{ii} \right) & \mbox{if} & i=j\\
\frac{1}{2}\left(\alpha_i a_{ij} + \alpha_j a_{ji}\right) & \mb{if}& i \ne j.
\end{array}\right.
\]
We assume that there exists a vector $\mb{\bm$\alpha$}=\left[\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$}$ is negative definite.

We now state and prove a theorem concerned with the exponential stability of the neural network model. This theorem can be found in Michel et al.~\cite{michel89}.
\begin{thm}
The null solution of the neural network model~(\ref{eq:hopasintersy}) 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}. 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 {\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~(\ref{eq:hopasintersy}) 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$
\section{Notes}
This chapter offers a glimpse at 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.
