\documentclass[11pt]{article}
\topmargin=0cm
\oddsidemargin=0truecm
\evensidemargin=0truecm
\headheight=0truecm
\headsep=0truecm
\textheight=23 truecm
\textwidth=15.8 truecm
\renewcommand{\baselinestretch}{1.7}
\newcommand{\mbi}[1]{\mbox{\boldmath $#1$}}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{remark}{Remark}[section]
\newtheorem{corollary}{Corollary}[section]
\newtheorem{example}{Example}[section]
\newtheorem{lemma}{Lemma}[section]
\makeatletter
%\renewcommand{\theequation}{%
%\thesection.\arabic{equation}}
%\@addtoreset{equation}{section}
%\makeatother
\begin{document}
\title{Mathematical Neuroscience --- An Overview}
\author{Shun-ichi Amari\\
RIKEN Brain Science Institute}
\date{}

\maketitle

\begin{abstract}
Mathematical neuroscience studies the fundamental principles of
information processing that has realized in the brain through long
history of evolution.  We use simple and tractable models of neural
networks and establish rigorous mathematical theories to study the
performances of the models.  This reveals the secrets why the brain
works well by using a large number of stochastically fluctuating neurons
and parallel dynamics of mutual interactions.  We give a number of
examples of mathematical neuroscience approaches, including statistical
neurodynamics, associative memory, neural field theory and dynamics of
learning in a hierarchical system including singularity.
\end{abstract}

\section{Introduction}

The brain is a very complex system consisting of a huge number of
neurons.  It forms a system composed of complicated networks, storing
information as memories and processing information to make decision.
One may say that it is a most sophisticated system that the nature has
ever created.  There are many approaches of research to the brain.

Brain science searches for the mechanisms and functions of the brain.
Researchers study, from the microscopic point of view, the role of genes
and proteins to understand the structure and function of neurons by
using the techniques of molecular biology.  They also use molecular
imaging techniques.  They search for the molecular mechanisms of
synapses and their plastic changes.  Systems neuroscientists observe the
activities of neurons of behaving animals to study how information is
represented in the brain and how it is processed.

Recent developments of measurement technology makes it possible to
observe the activities of a human brain by using fMRI, PET, NIRS, EEG,
MEG and even multi-electrodes.  Cognitive neuroscience makes use of
these new technologies.  Computational neuroscience builds models of
neural networks and studies their behaviors by theory or computer
simulations.  Researchers try to build realistic models by taking
experimental details into account, so that they can explain the
behaviors of actual networks by theory.  See, for example, Dayan and
Abott, 2001.

The brain is an organ which a living animal is equipped with for
information processing.  A living object is material, but is special in
the sense that it has an ability of self-reproduction.  It has emerged
through the long history of evolution.  Information, expressed in terms
of DNA molecules, is necessary for this.  Hence, a living object keeps
information in material.  Further, the brain processes information in
the environment and decides actions.  Higher-order information
processing ability as well as structured memory are necessary for this
purpose.

The brain has been developed through long years of evolution, and it
processes information in neuronal networks.  Material has become to be
able to keep information and to process it.  From this point of view,
one may say that information utilizes material to develop itself.  The
brain sits at the cross point of material and information.  It is
important to study the material basis of the brain.  At the same time,
it is necessary to study the principles of information which have
realized in material.

Computational neuroscience builds realistic models as faithful to
experimental details as possible.  This is necessary for explaining
actual processes taking place in the brain.  The results can then be
compared with experimental findings.  This is the orthodox way of
theory.

However, there would be another approach.  The principles of brain
information processing are obviously very different from those of
computers.  Information is distributed over a huge number of neurons and
processed by their dynamic interactions.  We need to understand the
principles which guarantee efficient, reliable and swift information
processing in such a distributed system with dynamic interactions.  The
nature has found such principles through long years of evolution and the
brain is a realization of the principles.

We need abstract models of neural networks in order to reveal the
principles of parallel and dynamic information processing, which are
simple enough to be able to analyze their capabilities deeply.
Mathematical neuroscience is an area that studies the possibilities,
limitations and performances of information processing by using
distributed parallel dynamics (Amari, 1990).  It aims at establishing
mathematical theories, where models are as simple as possible so long as
they do not miss essential features.  The nature has realized these
principles by using realistic neurons, so that the brain is far more
complex than the mathematical models, but they share common principles.
If we understand such principles, we can then proceed further to search
for realistic models and to find evidences of realization of such
principles in the brain.

The present paper shows examples of mathematical neurosciences and its
methods briefly.

\section{Statistical Neurodynamics (1): Simple Random Networks}

Let us consider a very simple model of mathematical neuron, called the
McCulloch-Pitts neuron.  Let us further consider a network composed of
$n$ neurons, and let $x_i(t), \cdots, x_n(t)$ be the outputs of neurons
at discrete time $t$.  We assume that their values are $1$ or $-1$,
representing firing or non-firing of the neurons.  Neuron $i$ receives
inputs from other neurons, which are $x_1(t-1), \cdots, x_n(t-1)$.  Let
$w_{ij}$ be the connection weight from neuron $j$ to neuron $i$.  We
neglect the threshold term for simplicity and also external inputs.
Then, the state of each neuron changes over time as
\begin{equation}
 x_i(t) = \mbox{sgn} \left\{
 \sum_j w_{ij}x_j(t-1) \right\},
\end{equation}
where sgn is the signature function,
\begin{equation}
 \mbox{sgn}(u) = \left\{
 \begin{array}{rr}
  1, & u \ge 0, \\
  -1, & u<0.
 \end{array}
 \right.
\end{equation}

Let us use abbreviated notations:  We use vector
\begin{equation}
 {\mbi{x}}(t) = \left[x_1(t) \cdots x_n(t)\right]
\end{equation}
for the state of the network at time $t$, matrix $W=\left(w_{ij}\right)$
for connections, and $T_W$ for state transition operator which is
nonlinear, so that (1) is rewritten as
\begin{equation}
 \label{eq:am4}
 {\mbi{x}}(t) = T_W{\mbi{x}}(t-1).
\end{equation}
We study the dynamics (\ref{eq:am4}) described by the state transition
operator $T_W$.

The properties of $T_W$ obviously depends on the matrix $W$.
Statistical neurodynamics (Amari, 1971; Amari, 1974; Amari, Yoshida and
Kanatani, 1977) assumes that $W$ is randomly generated from a
probability distribution $p(W)$, and searches for the properties that
hold for almost all networks having $W$ generated subject to $p(W)$ when
$n$ is large.  We give two examples.

When $w_{ij}$ are iid (identical and independently distributed) random
variables subject to a distribution $p(w)$, it is easy to apply the
central limit theorem and the law of large numbers, when $n$ is large.
Here, we assume that $p(w)$ is the Gaussian distribution with mean
$\bar{w}/\sqrt{n}$ and variance $\sigma^2_w$.  We can easily treat a
general distribution by a similar method.  What are common aspects of
dynamics of these networks?  To answer this question, we introduce
macroscopic variables.  The first one is the activity level (Amari,
1971),
\begin{equation}
 X({\mbi{x}})= \frac 1n \sum x_i.
\end{equation}
Then, the macroscopic variable at time $t$ is
\begin{equation}
 X(t) = \frac 1n \sum x_i(t)
\end{equation}
and we have a macroscopic state transition law of $X$, if we could find
a function $F$ which satisfies
\begin{equation}
 X(t)= F \left\{X(t-1)\right\}.
\end{equation}
We search for such a function $F$.  To this end, we use the fact that,
for a given non-random $x_1, \cdots, x_n$,
\begin{equation}
 u_i = \sum_j w_{ij}x_j
\end{equation}
are independently and identically distributed Gaussian random variables
with mean and variance,
\begin{equation}
 E \left[u_i \right] = \sqrt{n}\bar{w}X, \quad
 V \left[u_i \right]= n \sigma^2.
\end{equation}
Hence, by the law of large numbers, we have
\begin{equation}
 F(X) = \mbox{Prob}\left\{u_i>0 \right\}-\mbox{Prob}
 \left\{u_i<0\right\}.
\end{equation}
This is calculated as
\begin{equation}
 F(X) = \Phi \left(\frac{\bar{w}X}{\sigma}\right),
\end{equation}
where
\begin{equation}
 \Phi(u) = \frac 1{\sqrt{2 \pi}} \int^u_{-u} \exp
 \left\{-\frac{t^2}2\right\}dt.
\end{equation}

It is easy to prove that, given ${\mbi{x}}(0)$ at time $0$ of which
macroscopic variable $X(0)$ is given by
\begin{equation}
 X(0) = \frac 1n \sum x_i(0),
\end{equation}
the next macroscopic variable $X(1)$ is
\begin{equation}
 X(1) = \frac 1n \sum x_i(1) = F \left\{X(0)\right\}
\end{equation}
with a probability close to $1$ for almost all randomly generated
networks as $n$ tends to infinity.  However, there is difficulty for
proving that
\begin{equation}
 \label{eq:am15}
 X(2) = F \left\{X(1)\right\},
\end{equation}
or more generally
\begin{equation}
 \label{eq:am16}
 X(t) = F \left\{X(t-1)\right\}.
\end{equation}
This is because, even we show, for a fixed ${\mbi{x}}$,
\begin{equation}
 \label{eq:am17}
 X'=X \left(T_W {\mbi{x}} \right) = F \left(X({\mbi{x}})\right),
\end{equation}
this does not automatically guarantees that
\begin{equation}
 \label{eq:amX}
 X'' = X \left(T_W T_W {\mbi{x}} \right) = F \left\{F(X)\right\}.
\end{equation}
In order to derive (\ref{eq:am15}) or (\ref{eq:am17}), we used the fact
that $w_{ij}$ are independent.  However when we calculate
(\ref{eq:amX}) by using $T_W T_W {\mbi{x}}$, the two operators are no
more independent.  In other words, when we evaluate
\begin{equation}
 v_i(t) = \sum w_{ij}x_j (t),
\end{equation}
$x_j(t)$ are no more constants but depend on $W$.  Therefore, when we
deal with a macroscopic dynamics of a random networks, we need to verify
if the macroscopic equation of type (\ref{eq:amX}) is valid or not.
Computer stimulants show that the marcoscopie dynamics (\ref{eq:am16})
holds.  This is discussed in Rozonoer (1969, A, B, C) and Amari et al.
Amari (1974), Amari et al (1977).  Its theoretical justification is
still an unsolved problem.

We can extend the macroscopic dynamics to a network consisting of a
number of subnetworks, each of which consists of random networks of
different types of neurons.  The macroscopic variables are then
${\mbi{X}}=\left(X_1, \cdots X_m \right)$, each $X_{\alpha}$, $\alpha=1,
\cdots, m$, denoting the macroscopic activity of subnetwork $\alpha$.
We have the following macroscopic equation
\begin{equation}
 X_{\alpha} (t+1) = \Phi \left(
 \frac{\sum \bar{w}_{\alpha \beta}X_{\beta}}{\sigma_{\alpha}}
 \right).
\end{equation}
We can prove that these dynamical equations show monostable,
multistable, oscillatory and chaotic behaviors depending on the
probability distributions of $w^{(\alpha)}_{ij}$ within a subnetwork and
of $w^{(\alpha \beta)}_{ij}$ across subnetworks (Amari, 1971, 1972A).

It is interesting to search for other types of macroscopic variables.
Let us study the stability of the state transition $T_W$ of a simple
random network.  Assume that microscopic state ${\mbi{x}}$ is mapped to
${\mbi{y}}$ by $T_W$
\begin{equation}
 {\mbi{y}}' = T_W{\mbi{x}}.
\end{equation}
When ${\mbi{x}}$ is disturbed and changes to ${\mbi{x}}'$, how much
change is caused in the resultant output ${\mbi{y}}'=T_W{\mbi{x}}'$?
This is the problem of stability of $T_W$.  We introduce the Hamming
distance of two microscopic states
\begin{equation}
 D({\mbi{x}}, {\mbi{x}}') = \frac 1{2n} \sum \left|
  x_i-x'_i \right|
\end{equation}
to show their difference.  The problem is to find the relation
\begin{equation} 
  D' = K(D)
\end{equation}
where
\begin{equation}
 D' = D \left(T_W{\mbi{x}}, T_W{\mbi{x}}' \right).
\end{equation}
We can obtain this relation explicitly when $\bar{w}=0$ and
$\sigma^2=1$.  The answer given by Amari (1974) is
\begin{equation}
 K(D) = \frac 2{\pi} \sin^{-1} \sqrt{D}
\end{equation}
The state transition of a random network has a very specific property:
When $D$ is small,
\begin{equation}
 K(D) \approx \frac 2{\pi} \sqrt{D}.
\end{equation}
Since the derivative of $K(D)$ at $D=0$ is infinity, this means that a
small deviation is enlarged extraordinary.  So this type of network is
convenient for checking small deviations of inputs, since the difference
is expanded in the outputs.

If the macroscopic dynamics
\begin{equation}
 D_{t+1} = K \left(D_t \right)
\end{equation}
holds for two serieses ${\mbi{x}}_t=T^t_W{\mbi{x}}_0$ and ${\mbi{x}}'_t
= T^t_W{\mbi{x}}'_0$ and
\begin{equation}
 D_t=D \left({\mbi{x}}_t, {\mbi{x}}'_t \right),
\end{equation}
we can study the state transition diagram of a random network.  Since
the number of states is $2^n$, we have a state-transition graph having
$2^n$ nodes which correspond to $2^n$ ${\mbi{x}}$'s.  From each node
${\mbi{x}}$, there emerges one branch and its destination is
$T_W{\mbi{x}}$.  Hence, the graph has $2^n$ directed branches.  The
state transition graph is highly different from the random state
transition graph where the state transition branches are randomly
assigned to each node.  First, a random neural network has the
small-world property, where a small number of nodes monopolize branch
outlets.  Let $p_i$ be the probability that a node has $i$ ancestors,
\begin{equation}
 p_i = \mbox{Prob} \left\{ \left|T^{-1}_W{\mbi{x}}\right|=i \right\}. 
\end{equation}
Then, the distribution is subject to the power law
\begin{equation}
 p_i \propto \frac 1{i^{3-\alpha}},
\end{equation}
for some $0< \alpha < 1$.  See Amari (1974).  The small-world network
has a number of interesting properties, and one of them is that the
transient period of dynamics is very short, and the dynamics of
\begin{equation}
 {\mbi{x}} (t+1) = T_W{\mbi{x}}(t)
\end{equation}
converges quickly to its attractor states.

\section{Statistical Neurodynamics (2): Associative Memory}

The conventional model of associative memory networks can also be
analyzed by statistical neurodynamics.  See Amari (1972B), Hopfield
(1982) and many others.  But it has quite different behaviors from the
previous random networks, because the connections $w_{ij}$ are random
but not independent.  Memorized patterns in this network are represented
by the states ${\mbi{x}}$ which satisfy
\begin{equation}
 T_W{\mbi{x}} = {\mbi{x}}.
\end{equation}
Let ${\mbi{\xi}}_1,
\cdots, {\mbi{\xi}}_m$ be $m$ patterns to be memorized in a network such
that they are attractors of the dynamics
\begin{equation}
 T_W {\mbi{\xi}}_{\alpha} = {\mbi{\xi}}_{\alpha}, \quad
  \alpha=1, \cdots, m.
\end{equation} 
Starting at any initial state ${\mbi{x}}_0$, it is expected that the
state reaches one of the memorized pattern ${\mbi{\xi}}_{\alpha}$ by the
state transition, finding a memorized pattern.  When these patterns are
generated randomly, we compose the connection matrix in such a way that
\begin{equation}
 \label{eq:am34}
 w_{ij} = \frac 1m \sum \xi_{\alpha i} \xi_{\alpha j}
\end{equation}
This is a standard way of associative memory.  When these patterns are
orthogonal, it is easy to see that the patterns are attractors,
satisfying
\begin{equation}
 T_W {\mbi{\xi}}_{\alpha} = \mbox{sgn} \left(
 \frac 1m \sum_{\beta} {\mbi{\xi}}_{\alpha} \left(
 {\mbi{\xi}}_{\alpha} \cdot {\mbi{\xi}}_{\beta} \right)
 \right)= {\mbi{\xi}}_{\alpha}.
\end{equation}
When ${\mbi{\xi}}'_{\alpha}$'s are randomly generated subject to
independent and identical probabilities with $\mbox{Prob}
\left\{\xi_{\alpha i}=1 \right\}= \mbox{Prob} \left\{\xi_{\alpha
i}=-1\right\}=1/2$, they are asymptotically orthogonal.

The connection weight matrix $W$ is random, but the component $w_{ij}$
are not independent, because they are composed of (\ref{eq:am34}).  The
law of stability is interesting in this case.  The statistical analysis
shows that
\begin{equation}
 D' = K(D), \quad K(D)= 1-\frac 12 \Phi \left\{
 -\frac{(1-2D)n}m \right\}.
\end{equation}
Therefore, when $D$ is small, $D=K(D)$ is nearly equal to $0$.  Hence
the state transition of recalling a memorized pattern is very stable.
It is interesting to show that the state transition graph of an
associative memory model has a fractal structure (Amari and Maginu,
1988).  Here the basin of attractor has very complicated strange
boundaries.

\section{Pattern Dynamics in Neural Field}

A neural field is a continuous version of neural networks, where neurons
are located in a $2$-dimensional field like cortices.  A neural field
consists of layers including different types of neurons.  We give a
general equation describing the dynamics of excitation in the field,
\begin{equation}
 \label{eq:am37}
 \frac{\partial u_i ({\mbi{r}}, t)}{\partial t} = -u_i
 ({\mbi{r}}, t)+ \sum \int w_{ij} \left({\mbi{r}}, {\mbi{r}}' \right)
  f \left\{u_j \left({\mbi{r}}', t \right)\right\}d{\mbi{r}}'
  + I_i ({\mbi{r}}, t).
\end{equation}
Here, ${\mbi{r}}= \left(r_1, r_2 \right)$ is coordinates of the field,
denoting the locations of neurons, $u_i ({\mbi{r}}, t)$ is the average
membrane potential of neurons at location ${\mbi{r}}$ in the $i$-th
layer at time $t$,
\begin{equation}
 z_i({\mbi{r}}, t) = f \left[u_i({\mbi{r}}, t)\right]
\end{equation}
is the output of neurons whose average membrane potential is $u_i$ and
$f$ is called the activation function.  The functions
$w_{ij}\left({\mbi{r}}, {\mbi{r}}'\right)$ are connection weights from
the neurons at location ${\mbi{r}}'$ in the $j$-th layer to the neurons
at location ${\mbi{r}}$ in the $i$-th layer.  $I_i ({\mbi{r}}, t)$ are
external inputs to the neurons of $i$-th layer at position ${\mbi{r}}$
at time $t$ minus threshold.

This types of dynamics was proposed by Willson and Cown (1973).  A
simplified version of neural field used by Amari (1977) is a
one-dimensional field of one layer, where the activation function is the
step function
\begin{equation}
 f(u)= \left\{
 \begin{array}{rr}
  1, & u \ge 0, \\
  0, & u<0,
 \end{array}
 \right.
\end{equation}
and the connection weight function is homogeneous and isotropic, having
lateral-inhibition type,
\begin{equation}
 w \left(r, r' \right) = w \left(\left|r-r'\right|\right),
\end{equation}
that is, $w(r)$ is positive when $r$ is small and negative otherwise.

Given initial stimuli, we study how an initial excitation pattern $u(r,
0)$ develops in the dynamics.  When $u(r)>0$ holds in a region $R$ of
the field, neurons in the region are excited, $z(r)=1$, $r \in R$.  We
call such a region an excited region, that is a region
\begin{equation}
 R= \left\{r|u(r)>0\right\}.
\end{equation}
When an excited region is an interval $R=\left(a_1, a_2 \right)$, it is
called a localized excitation pattern or a bump.  We first focus on the
dynamics of localized excitation patterns.

Let us consider the case where $I$ is constant over the field.  In order
to analyze how the excited interval develops over time, we remark its
boundaries $a_1$ and $a_2$.  They change over time so that we write
$a_1(t)$ and $a_2(t)$ for the boundaries of the excited interval at time
$t$.  The boundaries satisfy
\begin{equation}
 u \left\{a_1(t), t \right\}= u \left\{a_2(t), t \right\} = 0
\end{equation}
at $t$.  After a short time $dt$, the boundaries change and they satisfy
\begin{equation}
 u \left\{a_i(t+dt), t+dt \right\}=0, \quad i=1, 2.
\end{equation}
We have, by Taylor expansion,
\begin{equation}
 \frac{\partial u \left(a_i, t\right)}{\partial t}+
  \frac{\partial u \left(a_i, t\right)}{\partial r}
 \frac{da_i}{dt}=0.
\end{equation}
The time derivative of $u \left(a_i, t \right)$ satisfies
(\ref{eq:am37}), and the integration part is written as
\begin{equation}
 \int w \left(r-r'\right)f \left\{u \left(r', t \right)\right\}dr'=
 \int^{a_2}_{a_1}w \left(r-r' \right)dr',
\end{equation}
since the activation function is $1$ in the excited interval $\left(a_1,
a_2 \right)$ and $0$ otherwise.  We define the integration of the
connection function by
\begin{equation}
 W(r) = \int^r_0 w(s)ds.
\end{equation}
Then, we have
\begin{equation}
 \int^{a_2}_{a_1}w \left(a_2-r' \right)dr'=W \left(a_2-a_1 \right).
\end{equation}
We put the derivatives of the waveform of pattern $u(r, t)$ at the
boundaries,
\begin{equation}
 c_1= \frac{\partial u \left(r_1, t \right)}{\partial r}, \quad
 c_2= -\frac{\partial u \left(r_2, t\right)}{\partial t}.
\end{equation}
Taking these into account, we have the equations describing the dynamics
of the boundary points,
\begin{equation}
 \frac{d r_i(t)}{dt} = -\frac 1{c_i}W \left(r_2-r_1\right)+I, \quad
  i=1, 2.
\end{equation}

Let
\begin{equation}
 a(t) = a_2(t)-a_1(t),
\end{equation}
which is the length of the excited interval.  It develops as
\begin{equation}
 \frac{da(t)}{dt} = \frac 1c \left\{
 W(a)+I \right\}.
\end{equation}
The equilibrium $a$ is the solution of
\begin{equation}
 W(a) + I = 0.
\end{equation}
Its stability is analyzed (Amari, 1977), and we have the following
results:

When
\begin{equation}
 W(a)+I=0
\end{equation}
has a solution, it is an equilibrium of the length of an excited region.
It is stable when
\begin{equation}
 w(a)>0
\end{equation}
is satisfied.

We show that, under a certain condition, we can find a stable
equilibrium solution.  The field can retain a local excitation pattern
without further inputs in such a case.  Since the field is homogeneous,
such a pattern can exist at any position.  Hence the field has an
infinite number of attractors (stable equilibrium solutions), which is
called a line attractor (Seung, 1988).  Assume that an initial stimulus
arrives at around position $b$.  The stimulus may be very noisy,
distributed around $b$.  Initial excitations are arisen at around $b$
and smoothed by the dynamic interactions in the field.  Even after the
stimulus disappears, a local excitation pattern remains stably at the
position $b$ when the field admits a local excitation pattern.  Hence,
this mechanism can be used as a working memory.  Such a mechanism is
believed to be used in memorizing input patterns and processing them.
One example is a neural field detecting orientations of input bars and
keeping the orientations for processing them further.  The neural field
theory explains the mechanism of the tuning curve in the visual cortex.

More general neural fields give richer behaviors (Ermentrout and Cowan,
1979).  See review paper by Coombes and Owen, 2005.  When a field
consists of an excitatory and inhibitory layers, the field admits a
stable travelling region.  Its velocity can be controlled by stimuli
given outside so that such phenomena will be useful for dynamic
information processing.

A two-dimensional neural field has much richer properties.  See, for
example, Coombes and Owen, 2005.  First, it possesses a localized
excitation pattern or a bump as well.  The stability analysis of the
two-dimensional bump was given in a Japanese book in 1978 (Amari), and
later proved by some others recently.  It also has interesting static
patterns consisting of a number of bumps or of a ring shape.  Their
dynamic behaviors are much more interesting.  It has an expanding ring
pattern, spiral pattern, breathing pattern and others.  It is shown
recently that it has also moving bumps.  The collision of moving bumps
shows very interesting dynamical behaviors.

The neural field model is an important component of computational
neuroscience.  Moreover, there are lots of applications of dynamic
neural field theory to psychological phenomena and robot navigation.

\section{Dynamics of Learning in Multilayer Perceptron and Gaussian Mixture}

We study the dynamical properties of on-line learning in a hierarchical
system such as multilayer perceptrons and Gaussian mixture radial basis
systems.  It is known that the behavior of a student system trained by
using examples given by the teacher system converges to the optimal
value, although one cannot avoid a local optimum.  However, it is also
known that the dynamics of such a learning system suffers from the
so-called plateau phenomenon, that the trajectory of learning is
attracted to a plateau which is not optimal.  It takes long time for the
system to get rid of such a plateau.  The present section analyzes the
dynamical behavior of on-line learning in the neighborhood of a plateau,
based on Wei et al. 2008, Wei and Amari, 2008, Cousseau, Ozeki and
Amari, 2008.

Plateaus are given rise to by singularities of the parameter space on
which learning takes place.  This is caused by the symmetric structure
of a hierarchical system and we cannot avoid them (Amari, Park and
Ozeki, 2006).  

Let us consider a function
\begin{equation}
 y=k({\mbi{x}}),
\end{equation}
where ${\mbi{x}}$ is an $n$-dimensional vector ${\mbi{x}}= \left(x_1
\cdots x_n \right)$ and $y$ is a scalar.  We assume that ${\mbi{x}}$ is
an input to a system and $y$ is its output.  In order to realize this
function, we use a family of functions
\begin{equation}
 F=\left\{f({\mbi{x}}, {\mbi{\theta}})\right\}
\end{equation}
parameterized by
\begin{equation}
 {\mbi{\theta}}= \left(\theta_1 \cdots \theta_n \right).
\end{equation}

Our problem is to estimate ${\mbi{\theta}}$ from a series of
input-output examples $D=\left\{\left({\mbi{x}}_1, y_1 \right), \cdots,
\left({\mbi{x}}_N, y_N \right)\right\}$ which we call the training
examples.  Here, ${\mbi{x}}_t$ is generated subject to a probability
distribution $q(\mbi{x})$ each time independently, and $y_t$ is the
output assumed to be generate by
\begin{equation}
 \label{eq:am58}
 y_t = k({\mbi{x}}_t)+n_t,
\end{equation}
that is, $y_t$ is the output from the true system contaminated by noise
$n_t$ subject to the standard Gaussian distribution.  We assume it is
subject to $N(0, 1)$ without loss of generality, by rescaling $y$
dividing it by the standard deviation of the noise.

When the true function $k(\mbi{x})$ from which the training examples are
generated belongs to $F$,
\begin{equation}
 k({\mbi{x}}) = f \left({\mbi{x}}, {\mbi{\theta}}^{\ast}\right), 
\end{equation}
the problem is to estimate the true value ${\mbi{\theta}}^{\ast}$ from
the training examples.  When $k(\mbi{x})$ does not belong to $F$, we
search for the parameter ${\mbi{\theta}}^{\ast}$ by which $k(\mbi{x})$
is approximated optimally by $f\left({\mbi{x}},
{\mbi{\theta}}^{\ast}\right)$.  Hence, this is a regression problem for
the data $D$.

We use an on-line learning approach such that an estimator
${\mbi{\theta}}_t$ at time $t$ is modified to ${\mbi{\theta}}_{t+1}$ by
using a new example ${\mbi{x}}_t, y_t$ at each discrete time $t, t=1\;
2, \cdots$.  This is a method of sequential estimation or regression,
called on-line learning.  The rule of modifying the parameter is written
as
\begin{equation}
 {\mbi{\theta}}_{t+1}= {\mbi{\theta}}_t-\eta {\mbi{g}}
 \left({\mbi{x}}_t, y_t, {\mbi{\theta}}_t \right)
\end{equation}
by using an adequate function ${\mbi{g}}$.  Here, $\eta$ is called a
learning constant, which may depend on $t$.  We study the properties of
dynamics of the above on-line learning method.

We consider the following two types of the parametric families of
functions as $F$.  One is multilayer perceptrons, written as
\begin{equation}
 f({\mbi{x}}, {\mbi{\theta}}) = \sum^m_{j=1} v_j \varphi
 \left(\sum^n_{i=1} w_{ji}x_i \right).
\end{equation}
Here, the system consists of two layers.  The first is the hidden layer
and the output of the $j$-th hidden neuron is a nonlinear function of
weighted sum of inputs ${\mbi{x}}$,
\begin{equation}
 z_j = \varphi \left(\sum_i w_{ji} x_i \right) =
 \varphi \left({\mbi{w}}_j \cdot {\mbi{x}}\right),
\end{equation}
where $\varphi$ is a sigmoidal function and and we use
\begin{equation}
 \varphi(u) = \frac 1{\sqrt{2 \pi}} \int^u_{-\infty} \exp
 \left\{-\frac 12 s^2 \right\}ds,
\end{equation}
as an example.  The second layer consists of a single neuron, called the
output neuron and the output is a linear summation of the outputs of the
hidden neurons,
\begin{equation}
 z=f({\mbi{x}}, {\mbi{\theta}}) = \sum v_j z_j.
\end{equation}
The parameters to be estimated are
\begin{equation}
 {\mbi{\theta}} = \left(v_1, \cdots, v_m \;;\; w_{11}, \cdots,
		   w_{mn}\right). 
\end{equation}

The second example uses Gaussian radial basis functions given by
\begin{equation}
 f({\mbi{x}}, {\mbi{\theta}}) = \sum v_j \psi
  \left({\mbi{x}}-{\mbi{w}}_j \right)
\end{equation}
where
\begin{equation}
 \psi({\mbi{x}}) = \frac 1{\sqrt{2 \pi}^n} \exp 
 \left\{-\frac{|{\mbi{x}}|^2}2 \right\}
\end{equation}
is the Gaussian function.  This is similar to the first system of
multilayer perceptrons.  The parameters to be estimated are
\begin{equation}
 {\mbi{\theta}} = \left(u_1, \cdots, v_m \;;\; {\mbi{w}}_1, \cdots,
		   {\mbi{w}}_m \right).
\end{equation}

We use the stochastic gradient learning method, minimizing a cost
function each time.  Here, the cost function $l({\mbi{x}}, y,
{\mbi{\theta}})$ is given by a half of the square of the difference of
the actual output $y$ and the hypothetical output $f({\mbi{x}},
{\mbi{\theta}})$ which the system with parameter ${\mbi{\theta}}$ emits
for given input ${\mbi{x}}$.  This is easily calculated for given
parameter ${\mbi{\theta}}$ when an input-output example $({\mbi{x}}, y)$
is given,
\begin{equation}
 \label{eq:am69}
 l({\mbi{x}}, y, {\mbi{\theta}}) = \frac 12
 \left\{
 y-f({\mbi{x}}, {\mbi{\theta}}) \right\}^2.
\end{equation}
The difference is called the error function,
\begin{equation}
 e \left({\mbi{x}}, y, {\mbi{\theta}}\right)=
 y-f({\mbi{x}}, {\mbi{\theta}}).
\end{equation}

The true system of emitting $y$ can be described by a conditional
probability distribution of emitting $y$ from $f({\mbi{x}},
{\mbi{\theta}})$ with additive noise $n$ (\ref{eq:am58}).  Then, the
conditional probability distribution of $y$ under the condition that
input ${\mbi{x}}$ is given is
\begin{equation}
 p(y|{\mbi{x}}, {\mbi{\theta}})= \frac 1{\sqrt{2 \pi}} \exp
 \left[-\frac 12 \left\{y-f({\mbi{x}}, {\mbi{\theta}})\right\}^2
 \right]
\end{equation}
when the parameter is ${\mbi{\theta}}$.  The cost function
(\ref{eq:am69}) is interpreted as the negative of $\log$ likelihood
\begin{equation}
 \label{eq:am72}
 l({\mbi{x}}, y, {\mbi{\theta}}) = -\log p
 \left(y|{\mbi{x}}, {\mbi{\theta}}\right).
\end{equation}
Hence, minimizing the cost function is equivalent to maximizing the
$\log$ likelihood.  This gives the maximum likelihood estimator.

The dynamics of the stochastic gradient on-line learning method is
written as
\begin{equation}
 {\mbi{g}}({\mbi{x}}, y, {\mbi{\theta}})= \nabla l
 \left({\mbi{x}}, y, {\mbi{\theta}}\right) = -e
 \left({\mbi{x}}, y, {\mbi{\theta}}\right) \nabla
 f({\mbi{x}}, y, {\mbi{\theta}})
\end{equation}
where $\nabla$ is the gradient with respect to ${\mbi{\theta}}$,
\begin{equation}
 \nabla = \left(
  \frac{\partial}{\partial \theta_1} \cdots
 \frac{\partial}{\partial \theta_n}
 \right).
\end{equation}
In order to analyze the behavior of dynamics of on-line learning, we use
the average learning equation
\begin{equation}
 {\mbi{\theta}}_{t+1}= {\mbi{\theta}}_t-\eta \left<
 \nabla l ({\mbi{x}}, y, {\mbi{\theta}}) \right>
\end{equation}
where $<\quad>$ denotes the average over all possible input ${\mbi{x}}$
by using the probability distribution (\ref{eq:am72}).  We further use
the continuous time approximation, and analyze the continuous time
average learning equation,
\begin{equation}
 \frac{d{\mbi{\theta}}(t)}{dt} = -\left<
 \nabla l ({\mbi{x}}, y, {\mbi{\theta}}) \right>.
\end{equation}

Let $S=\left\{{\mbi{\theta}}\right\}$ be the parameter space.  When
\begin{equation}
 f \left({\mbi{x}}, {\mbi{\theta}}_1 \right) \ne
 f \left({\mbi{x}}, {\mbi{\theta}}_2 \right) \; \mbox{for}\;
 {\mbi{\theta}}_1 \ne {\mbi{\theta}}_2
\end{equation}
holds, the system is identifiable, because only one ${\mbi{\theta}}$
corresponds to a function $f({\mbi{x}}, {\mbi{\theta}})$.  However, when
there exist at least two points for which
\begin{equation}
 f \left({\mbi{x}}, {\mbi{\theta}}_1 \right) = f
 \left({\mbi{x}}, {\mbi{\theta}}_2 \right), \quad
 {\mbi{\theta}}_1 \ne {\mbi{\theta}}_2
\end{equation}
holds, the system is not identifiable.  Let $R({\mbi{\theta}})$ be a
subset of $S$ in which all the output function is identical,
\begin{equation}
 f({\mbi{x}}, {\mbi{\theta}}) = f \left({\mbi{x}}, {\mbi{\theta}}'
				  \right), \quad
 {\mbi{\theta}}' \in R ({\mbi{\theta}}).
\end{equation}
Then, the all parameter points ${\mbi{\theta}}'$ in $R({\mbi{\theta}})$
are equivalent from the point of view of the input-output behavior.
There are many such subsets $R({\mbi{\theta}})'$s.  We divide the space
$S$ by this equivalence relation.  Then all the points in a
$R({\mbi{\theta}})$ reduces to one behavior point.  The topological
structure of $S$ is destroyed by this division, because every
$R({\mbi{\theta}})$ reduces to a single point.  This gives rise to
singularities.  We call $R({\mbi{\theta}})$ a singular region.  We give
two typical examples of singularity in a hierarchical system.

{\textbf{1) Eliminating singularity:}}

When $v_i=0$, whatever ${\mbi{w}}_i$ is, $v_i \varphi \left({\mbi{w}}_i
	  \cdot {\mbi{x}}\right)=0$, so that
	  $R=\left\{{\mbi{\theta}}|v_i=0, {\mbi{w}}_i\;
	  \mbox{arbitrary} \right\}$ is a singular region.  Since
	  neuron $i$ is eliminated because of $v_i=0$, it is called an
	  eliminating singularity region.

{\textbf{2) Overlapping singularity:}} 

When ${\mbi{w}}_i={\mbi{w}}_j ={\mbi{w}}$, so long as $v_i+v_j=v$ holds,
\begin{equation}
 v_i \varphi \left({\mbi{w}}_i \cdot {\mbi{x}}\right) + v_j
 \varphi \left({\mbi{w}}_j \cdot {\mbi{x}}\right) = 
 v \varphi ({\mbi{w}} \cdot {\mbi{x}}).
\end{equation}
So we cannot identify each of $v_i$ and $v_j$.  The region
\begin{equation}
 R= \left\{{\mbi{\theta}}\left|{\mbi{w}}_i = {\mbi{w}}_j = {\mbi{w}},
 v_i+v_i=v \right. \right\}
\end{equation}
is a singular region, depending on ${\mbi{w}}$ and $v$.  Since the two
neurons become identical in this case, we call it an overlapping
singularity.

There are other singularities where three neurons overlap and so on, but
the above two are typical and other singularities are intersections of
the above singularity regions.  Therefore, we treat a simple system
consisting of two hidden neurons, where ${\mbi{\theta}}= \left(v_1, v_2,
{\mbi{w}}_2, {\mbi{w}}_2 \right)$ and study the dynamical behavior of
learning in such a system.

In the critical region $R(v, {\mbi{w}})$, only one neuron suffices,
where ${\mbi{w}}_1 = {\mbi{w}}_2={\mbi{w}}$ or one of $v_1$ and $v_2$ is
$0$.  In this case, its input-output function is written by using a
single neuron,
\begin{equation}
 \label{eq:am18}
 f ({\mbi{x}}, {\mbi{\theta}}) =  v \varphi \left({\mbi{w}}\cdot
					 {\mbi{x}}\right).
\end{equation}
The situation is the same in the Gaussian mixture case.  

We introduce a new coordinate system ${\mbi{\xi}}= \left({\mbi{v}}, w,
{\mbi{u}}, z \right)$, defined by
%\begin{eqnarray}
\begin{equation}
 \begin{array}{ll}
  {\mbi{u}}= {\mbi{w}}_2-{\mbi{w}}_1, & 
 {\mbi{w}}= {\displaystyle{\frac 12}}
  \left({\mbi{w}}_1+{\mbi{w}}_2 \right), \\
  z = {\displaystyle{\frac{v_1-v_2}{v_1+v_2}}}, &  v=v_1+v_2.
 \end{array}
%\end{eqnarray}
\end{equation}
Then, the singular region is represented by ${\mbi{u}}=0$ or $z= \pm 1$.
 Given $(v, {\mbi{w}})$, singular region $R(v, {\mbi{w}})$ exists.  A
 singular region ${\cal{R}}$ consists of three parts, one given by
 ${\mbi{u}}=0$, that is, ${\mbi{w}}_1={\mbi{w}}_2$, which we call the
 overlapping singular region, and the other is $z=\pm 1$, that is
 $w_1=0$ or $w_2=0$, in which one neuron is eliminated.  We have $R =
 R_1 \cup R_2$,
\begin{eqnarray}
 R_1 \left({\mbi{v}}, w \right) &=& \left\{ {\mbi{\xi}}
      \left|{\mbi{u}}=0, z: \mbox{arbitrary}; {\mbi{v}}\; \mbox{and}\; w \;
      \mbox{fixed} \right.\right\}, \\
 R_2 \left({\mbi{v}}, w \right) &=& \left\{{\mbi{\xi}} \left| {\mbi{u}}:
						       \mbox{arbitrary};
						       z= \pm 1, {\mbi{v}}\;
						       \mbox{and}\; w \;
						       \mbox{fixed} 
						       \right. \right\}.
\end{eqnarray}
In the subspace in which ${\mbi{v}}$ and $w$ are fixed constants, $R_1$
is a line in which $z$ changes, and $R_2$ consists of two
$n$-dimensional surfaces.  The singular region $R$ is a composite of one
line and two surfaces.

We study the dynamical flow near the singularity $R_1 \left({\mbi{v}}, w
\right)$, by using the Taylor expansion with respect to ${\mbi{u}}$,
where $|{\mbi{u}}|$ is small.  The average learning equation is expanded
as
\begin{eqnarray}
 \dot{\mbi{u}} 
& = & \eta \frac w 2(1-z^2) 
       \left\langle e \varphi''(\mbi{v}\cdot\mbi{x})(\mbi{u}\cdot\mbi{x})\mbi{x} \right\rangle
             +O(|\mbi{u}|^2), \label{eq:am22}\\
 \dot{z} & = & 
      -\frac{1+z^2}{w}\eta\left\langle\varphi'(\mbi{v}\cdot\mbi{x})(\mbi{u}\cdot\mbi{x})\right\rangle 
           +O(|\mbi{u}|^2). \label{eq:am23}
\end{eqnarray}

We can integrate (\ref{eq:am22}) and (\ref{eq:am23}) in the subspace in
which $S^{\ast}$ where ${\mbi{v}}={\mbi{v}}^{\ast}$ and $w=v^{\ast}$ are
the optimal values, that is $v^{\ast} \varphi \left({\mbi{w}}^{\ast}
\cdot {\mbi{x}}\right)$ is the optimal approximation of the target
foundation $k({\mbi{x}})$ by using only a single neuron.  This gives the
flows or trajectories of dynamics in the neighborhood of ${\mbi{u}}=0$
of $S^{\ast}$.  They are given by
\begin{eqnarray}
\frac 1 2 \mbi{u}^T\cdot \mbi{u} = \frac {2w^{\ast 2}}{3} \log
   \left[ \frac{(z^2+3)^2}{|z|}\right] + c
\end{eqnarray}
where $c$ is an arbitrary constant.  


\section{Milnor attractor}

The line $R_1 \left({\mbi{v}}^{\ast}, w^{\ast} \right)$ in $S^{\ast}$ is
a critical line of the dynamics of learning.  Under a certain condition,
we show that the line is divided into two parts such that one part is
stable (attractive) and the other is unstable (repulsive).  In such a
case, there are trajectories attracted to the stable part, and the basin
of attraction has a finite measure.  This is the Milnor type attractor
(Milnor, 1985), different from a saddle point.  More strongly, the basin
of attraction includes a neighborhood of this part, so that all the
points near this part are once attracted to it.  After being attracted
to it, the parameters still move randomly in $R_1
\left({\mbi{v}}^{\ast}, w^{\ast}\right)$ by stochastic fluctuations
until they reach the unstable part (Fukumizu and Amari, 2000).  Then,
the parameters leave $R_1$ eventually.

We can prove that $R^{\ast}_1$ is unstable (saddle) when $A=
\left(A_{ij}\right)$, $A_{ij}= \left<e \phi_{ij} \right>$, where
\begin{equation}
 \phi_{ij} = \frac{\partial^2}{\partial w_i \partial w_j} 
 \varphi ({\mbi{w}}\cdot {\mbi{x}})
\end{equation}
includes both
positive and negative eigenvalues.  When $w^{\ast}A$ is
positive-semi-definite, the interval $1<z^2$ is attractive, and the
remaining part $1>z^2$ is repulsive.  Hence, the interval $1<z^2$
behaves like a Milnor attractor.  This is the plateau phenomenon, taking
long time before getting rid of it.  When $w^{\ast}A$ is
negative-semi-definite, the part $1>z^2$ is attractive and the other
part is repulsive.

\section{Conclusions}

We have shown examples of mathematical neuroscience approaches.  They
include 1) the macro-dynamics of randomly connected networks and the
characteristics of state transition are elucidated.  2) The dynamics of
associative memory networks is elucidated from the point of view of
macro-dynamics. 3) Pattern formations and pattern interactions in a
neural field are investigated by using integro-differential equations.
4) It is studied how the singularity existing in a hierarchical system
affects dynamics of learning.

Mathematical approaches are necessary for finding the principles of
information processing which the brain utilizes.  We need to use simple
and tractable models for carrying mathematical analysis.  The models are
different from the real biological neural networks, but we hope the
mathematical theory clarifies the principles from which we can construct
realistic models.  Mathematical neuroscience has not been fully
established and much more efforts should be devoted for constructing
such theories.


\begin{thebibliography}{}
 \bibitem{DA2001} P. Dayan and L.F. Abbott, Theoretical Neuroscience:
	 Computational and Mathematical Modeling of Neural Systems, MIT
	 Press, 2001.
%
 \bibitem{amari1990} S. Amari, Mathematical Foundations of
	 Neurocomputing, {\em Proceedings of the IEEE}, 78, 9,
	 1443--1463, 1990.
%
 \bibitem{amari1971} S. Amari, Characteristics of Randomly Connected
	 Threshold-Element Networks and Network Systems, {\em
	 Proc. IEEE.}, 59, 1, 35--47, 1971.
%
 \bibitem{amari1974} S. Amari, A Method of Statistical Neurodynamics,
	 {\em Kybernetik}, 14, 201--215, (Heft 4) 1974.
%
 \bibitem{AYK1977} S. Amari, K. Yoshida and K. Kanatani, A Mathematical
	 Foundation for Statistical Neurodynamics, {\em SIAM
	 J. Appl. Math.}, 33, 95--126, 1977.
%
 \bibitem{rozonoer69A} L.I. Rozonoer, Random logical nets;
	 I, {\em Avtomat. Telemekh.}, 5, 137--147, 1969A.
%
 \bibitem{rozonoer69b} L.I. Rozonoer, Random logical nets;
	 II, {\em Avtomat. Telemekh.}, 6, 99--109, 1969B.
%
 \bibitem{rozonoer69c} L. I. Rozonoer, Random logical nets;
	 III. {\em Avtomat. Telemekh.}, 7, 127--136, 1969C.
%
 \bibitem{amari1972} S. Amari, Characteristics of Random Nets of Analog
	 Neuron-Like Elements, {\em IEEE Trans. Systems, Man and
	 Cybernetics}, SMC-2, 5, 643--657, November 1972A.
%
 \bibitem{amari1972} S. Amari, Learning Patterns and Pattern Sequences by
	 Self-Organizing Nets of Threshold Elements, {\em IEEE
	 Trans. Computers}, C-21, 11, 1197--1206, 1972B.
%
 \bibitem{hopfield1982} J.J. Hopfield, Neural network and
	 physical systems with emergent collective computational
	 abilities, {\em Proceedings of the National Academy of Sciences of
	 the United States of America}, 79,8, 2554--8, 1982
%
 \bibitem{AM1988} S. Amari and K. Maginu, Statistical neurodynamics of
	 associative memory, {\em Neural Networks}, 1, 1, 63--73, 1988.
%
 \bibitem{WC1973} H.R. Wilson and J.D. Cowan, A mathematical
	 theory of the functional dynamics of cortical and thalamic
	 nervous tissue, {\em Kybernetik}, 13, 55--80, 1973.
%
 \bibitem{amari1977} S. Amari, Dynamics of Pattern Formation in
	 Lateral-Inhibition Type Neural Fields, {\em Biological
	 Cybernetics}, 27, 77--87, 1977.
%
 \bibitem{seung1998} H. S. Seung, Continuous attractors and oculomotor
	 control, {\em Neural Networks}, 11, 1253--58, 1998.
%
 \bibitem{EC1979} G.B. Ermentrout and J.D. Cowan, A mathematical theory
of visual hallucinations patterns, {\em Biological Cybernetics}, 34,
137--150, 1979.
%
 \bibitem{CO2005} S. Coombes and M.R. Owen, Bumps, breathers, and waves
	 in a neural network with spike frequency adaptation, {\em Physical
	 Review Letters}, 94, 148102, 2005.
%
 \bibitem{WZCOA08} H. Wei, J. Zhang, F. Cousseau, T. Ozeki and S. Amari,
	 Dynamics of Learning Near Singularities in Layered Networks,
	 {\em Neural Computation}, 20, 813--843, 2008.
%
 \bibitem{WA2008} H. Wei and S. Amari, Dynamics of learning near
	 singularities in radial basis function networks, {\em Neural
	 Networks}, 21, 989--1005, 2008.
%
 \bibitem{COA08} F. Cousseau, T. Ozeki and S. Amari, Dynamics of
	 Learning in Multilayer Perceptrons Near Singularities, {\em
	 IEEE Transactions on Neural Networks}, 19, 8, 1313--1328, 2008.
%
 \bibitem{APO06} S. Amari, H. Park and T. Ozeki, Singularities Affect
	 Dynamics of Learning in Neuromanifolds, {\em Neural
	 Computation}, 18, 1007--1065, 2006.
%
 \bibitem{milnor1985} J. Milnor, On the concept of attractor, {\em
	 Commum. Math. Phys.}, 99, 177--195, 1985.
%
 \bibitem{FA2000} K. Fukumizu and S. Amari, Local minima and plateaus in
	 hierarchical structures of multilayer perceptrons, {\em Neural
	 Networks}, 13, 317--327, 2000.
\end{thebibliography}



\end{document}