Neural Networks Lecture Notes Overview
Neural Networks Lecture Notes Overview
Lecture Notes
V 1.26, Sep 10, 2024
Bibliography 223
3
Fasten your seatbelts – we enter the world of neu-
ral networks
Let’s pick up some obvious facts from the surface and take a look at the bottomless
voids that open up underneath them.
What makes you YOU is your brain. Your brain is a neural network. A
neural network is a network made of neurons which connect to each other by
synaptic links. Thus, one of the big question of science and life is how YOU are a
NETWORK of interconnected NEURONS.
The answer seemed clear enough, 80 years ago, for the pioneers of what we
now call computational neuroscience. In the kick-start work on neural networks
(McCulloch and Pitts, 1943), a neuron x was cast as binary switch that could
have two states — call them 0 and 1, or false and true — and this neuron x be-
comes switched depending on the 0-1 states of other neurons x1 , . . . , xk which have
synaptic links to x. A brain thus was seen as a Boolean circuit. The final sentence
in that paper is “Thus in psychology, introspective, behavioristic or physiological,
the fundamental relations are those of two-valued logic.” In other words, brains
(and you) are digital computers.
But.
What followed is 80 years of scientific and philosophical dispute, sometimes
fierce, and with no winners to the present day. The more we have been learning
about neurons and brains and humans and robots and computers, the more con-
fusing the picture has become. As of now (the year 2023), neither of what is a
NEURON, a NETWORK, or a YOU is clear:
4
is almost impossible to observe experimentally. The Human Connectome
Project ([Link]
one of the largest national U.S. research programs in the last years, invested
a gigantic concerted effort to find out more and found that progress is slow.
And, on top of that, let me add that it is also becoming unclear again in these
days what a COMPUTER is (I’ll say more about that in the last lecture of this
course).
Given that the scientific study of neural networks is so closely tied up with
fundamental questions about ourselves, it is no wonder that enormous intellectual
energies have been spent in this field. Over the decades, neural network research
has produced a dazzling zoo of mathematical and computational models. They
range from detailed accounts of the functioning of small neural circuits, compris-
ing a few handfuls of neurons (celebrated: the modeling of a 30-neuron circuit in
crustaceans (Marder and Calabrese, 1996)), to global brain models of grammati-
cal and semantic language processing in humans (Hinaut et al., 2014); from low-
detail models of single neurons (1963 Nobel prize for Alan Hodgkin and Andrew
Huxley for a electrical engineering style formula describing the overall electrical
dynamics of a neuron in good approximation (Hodgkin and Huxley, 1952)), to su-
percomplex geometrical-physiological compartment models on the high-detail side
(Gouwens and Wilson, 2009)); from statistical physics oriented models that can
“only” explain how a piece of brain tissue maintains an average degree of activity
(van Vreeswijk and Hansel, 2001) to AI inspired models that attempt to explain
every millisecond in speech understanding (Shastri, 1999); from models that aim
to capture biological brains but thereby become so complex that they give satis-
factory simulation results only if they are super delicately fine-tuned (Freeman,
1987) to the super general and flexible and robust neural learning architectures
that have made deep learning the most powerful tool of modern machine learning
applications (Goodfellow et al., 2016); or from models springing from a few most
elegant mathematical principles like the Boltzmann machine (Ackley et al., 1985)
5
to complex neural architectures that are intuitively put together by their inventors
and which function well but leave the researchers without a glimmer of hope for
mathematical analysis (for instance the “neural Turing machine” of Graves et al.
(2014)).
In this course I want to unfold for you this world of wonder. My goal is to
make you aware of the richness of neural network research, and of the manifold
perspectives on cognitive processes afforded by neural network models. A student
of AI should, I am convinced, be able to look at “cognition” from many sides —
from application-driven machine learning to neuroscience to philosophical debates
— and at many levels of abstraction. All neural network models spring from the
same core idea, namely that intelligent information processing emerges from the
collective dynamics in networks of simple atomic processing units. This gives the
field a certain coherence. At the same time, it is amazing into how many different
directions one can step forward from this basic idea. My plan is to present a quite
diverse choice of neural network models, all classics of the field and must-know’s
for any serious AI/cognitive science/machine learning disciple. You will see that
the world of neural networks has so much more to offer than just “deep learning”
networks, which in these days outshine all other kinds in the public perception.
This said, I will nonetheless start the course with an introduction to that
currently most visible kind of neural networks, feedforward neural networks, from
the classical, simple Perceptrons via multilayer perceptrons to a number of deep
learning models. I position this material at the beginning of the course because
these models lend themselves to practical programming projects most directly,
such that you can swiftly start working on the practical project that accompanies
the theory lectures.
But, please, be aware that scientific tides come and go, and the day will come
when deep learning methods will recede in favor and some other magic will move
up front. Then the other themes and methods that you have learnt about in this
course will help you to connect with whatever else comes next. In this vein, it may
interest you to learn that the now-dominant paradigm of deep learning directly
emerged from quite another brand of neural networks, the Boltzmann machine
which you will get to know in a few weeks. The celebrated paper which today is
often seen as the starter for deep learning (Hinton and Salakuthdinov, 2006) in fact
was written from the perspective of the statistical physics which rules Boltzmann
machines, and the option to use them as an initialization submechanism for what
today are called deep networks was only mentioned in passing. A few years later,
the Boltzmann machine theme receded into the background and deep deep deep
became the motto of the day. Such shifts in focus will happen again! And I will
dare to give a forecast in the last lecture of this course.
6
use Excel!), which are just matrices if seen with the eye of a machine learner.
Furthermore, neural networks are shaped by connecting neurons with weighted
“synaptic links”, and these weights are again naturally sorted in matrices. And
the main operation that a neural network actually does is formalized by a matrix-
vector multiplication. So it’s matrices and vectors all over the place, no escape
possible. You will need at least a basic, robust understanding of linear algebra to
survive or even enjoy this course. We will arrange a linear algebra crash refresher
early in the course. A good free online resource is the book “Mathematics for
Machine Learning” (Deisenroth et al., 2019).
Furthermore, to a lesser degree, also some familiarity with statistics and prob-
ability is needed. You find a summary of the must-knows in the appendix of these
lecture notes, and again a tutorial exposition in Deisenroth et al. (2019).
Finally, a little (not much) calculus is needed to top it off. If you are familiar
with the notion of partial derivatives, that should do it. In case of doubt - again
it’s all (and more) in Deisenroth et al. (2019).
7
1.1 Training data
A supervised learning tasks starts from labelled training data, that is a sample
S = (ui , yi )i=1,...,N of input-output pairs (ui , yi ). For simplicity we assume that
the inputs and output are real-valued vectors, that is ui ∈ RK , yi ∈ RM , although
other formats are possible.
Brief note on notation: In most ML textbooks the inputs are denoted as x, not
u. However, we will later often deal with the internal processing states of neural
networks. Following traditions in signal processing and dynamical systems maths,
I will use the symbol x for these internal state vectors of a signal processing device
(for us: a neural network), and use u for external inputs.
Two examples:
• Image classification tasks: the input patterns ui are vectors whose entries are
the red, green, blue intensity values of the pixels of photographic images. For
instance, if the images are sized 600 × 800 pixels, the input pattern vectors
are of dimension K = 600 · 800 · 3 = 1, 440, 000. The output vectors might
then be binary “one-hot encodings” of the picture classes to be recognized.
For example, if the task is to recognize images of the ten sorts of handwritten
digits 0, 1, 2, ..., 9, the output vector yi for a “1” input image would be the
10-dimensional vector (0, 1, 0, . . . , 0)′ that has a 1 in the second component
and is zero everywhere else.
Quick note in passing: vectors come as row or column vectors. I always
use column vectors. Thus, when I write ui for the input vector representing
an image, it is a 1, 440, 000 dimensional column vector. To turn a column
vector into a row vector, or vice versa, one applies the transpose operation,
which I denote by a prime. Thus u′i is a row vector, and (0, 1, 0, . . . , 0)′ is a
column vector.
8
1.2 Training objectives
Given a training sample S = (ui , yi )i=1,...,N , the task of supervised learning is to
train a model of the input-output function f that gave rise to the observed training
data. Such an ML model is an algorithm fˆ which accepts inputs from RK and
computes outputs in RM . (A note on terminology: In statistical modeling, fˆ is an
estimate of the true function f . The ˆ· notation is often used to denote estimates
of something. Also, in statistics, functions that map input values to output values
— here we named them f or fˆ — are often called decision functions, a naming
that I adopted in the Machine Learning lecture notes.)
A machine learning algorithm is a computational procedure which gets a train-
ing sample S as input and “learns” (that is, computes — and a stastician would
say, estimates) a model fˆ : RK → RM . In machine learning, the model fˆ will
be an executable algorithm (for instance, a neural network). Thus, a machine
learning algorithm is an algorithm that transforms data into algorithms!
In order to decide whether an estimated model fˆ : RK → RM is “good” or
“bad”, one needs a well-defined measure to quantify the goodness of a model. This
is achieved by a loss function
L : RM × RM → R≥0 . (1)
The idea is that a loss function measures the “cost” of a mismatch between
the correct value y and the model output fˆ(u), given a correct input-output
pair (u, y) from the training data or from testing data. Higher cost means lower
quality of fˆ. By convention, loss functions are non-negative and higher values
means worse fit between true and estimated output. Many loss functions are in
use. An experienced machine learning expert may spend much time and care on
tailoring the loss function to the learning problem that s/he is dealing with. Three
basic examples:
• A loss that counts misclassifications in pattern classification tasks: assume
that one is dealing with a pattern classification set-up where the target
outputs y in the training/testing data and the outputs returned by a model
fˆ are m-dimensional binary “one-hot” class encoding vectors. Let fˆ : RK →
{0, 1}M be any candidate model. Then the loss
ˆ 0, if fˆ(u) = y
L(f (u), y) = (2)
1, if fˆ(u) 6= y
counts misclassifiation errors. It is sometimes called the counting loss.
• A very popular loss penalizes quadratic errors of vector-valued targets:
L(fˆ(u), y) = kfˆ(u) − yk2 . (3)
This loss is often just called “quadratic loss” or “L2 loss”, and it is one of
the most frequently used ones in regression tasks.
9
• In binary classification tasks (that is, the input patterns come from two
classes only, call them class A and class B), one wishes a model h that on
input of some pattern u returns an estimate of the probability fˆ(u) = P̂ (u ∈
A) that this pattern belongs to class A. To train such model behavior, one
usually uses the (binary) cross-entropy loss, which is given by
ˆ − log fˆ(u), if u is indeed in class A
L(f (u), y) = ˆ (4)
− log(1 − f (u)), if u is in class B.
This loss can be straightforwardly generalized to multi-class classification
tasks. It becomes minimal if the model-generated probabilities fˆ(u) =
P̂ (u ∈ A) are the true probabilities P (u ∈ A).
In order to understand the nature of supervised learning tasks, one has to
frame it in the context of probability theory. You find a summary of the necessary
probability concepts and notation in the appendix. We assume that the training
and testing input-output data u, y are obtained from random variables U, Y .
Learning algorithms should minimize the expected loss, that is, a good learning
algorithm should yield a model fˆ whose risk
R(fˆ) = E[L(fˆ(U ), Y )] (5)
is small. The expectation here is taken with respect to the true joint distribution
PU,Y of the data-generating random variables U and Y . For example, in a case
where U and Y are numerical RVs and their joint distribution is described by a
pdf p, the risk of a candidate model h would be given by
Z
R(h) = L(h(u), y) p(u, y) d(u, y).
RK ×RM
However, the true distribution PU,Y and its pdf p are unknown in real-world learn-
ing scenarios. The mission to find a model fˆ which minimizes (5) is, in fact,
hopeless. The only access to PU,Y that the learning algorithm affords is the scat-
tered reflection of PU,Y in the training sample (ui , yi )i=1,...,N .
A natural escape from this impasse is to tune a learning algorithm such that
instead of attempting to minimize the risk (5) it tries to minimize the empirical
risk
X
N
Remp (h) = 1/N L(h(ui ), yi ), (6)
i=1
which is just the mean loss averaged over the training examples. Minimizing this
empirical risk is an achievable goal, and a host of optimization algorithms for
all kinds of supervised learning tasks exist which do exactly this. In machine
learning, these optimization algorithms are called learning algorithms. That is,
such learning algorithms find
X
N
fˆ = hopt = argmin 1/N L(h(ui ), yi ). (7)
h∈H i=1
10
The set H is the hypothesis space – the search space within which a learning
algorithm may look for an optimal model h.
It is important to realize that every learning algorithm comes with a specific
hypothesis space. For instance, when one uses linear regression to solve (7), H
is the set of all linear functions from RK to RM (I assume you have learnt about
linear regression in your first year course on linear algebra or statistics). Or,
if one sets up a neural network learning algorithm, H is typically the set of all
neural networks that have a specific connection structure (number of neuron layers,
number of neuros per layer); the networks in H then differ from each other only
by the weights associated with the synaptic connections.
The empirical risk is often – especially in numerical function approximation
tasks – also called the training error.
Here is an interim take-home summary:
• The ultimate goal for supervised learning algorithms is to estimate a model
fˆ which has a low risk (5), that is, which on average gives low-loss (“good”)
outputs on “testing” data drawn from the distribution PU,Y . That is, the
ultimate goal is to solve the optimization problem
fˆ = hopt = argmin E[L(fˆ(U ), Y )]. (8)
h∈H
• The only source of information that the learning algorithm has is the training
sample S = (ui , yi )i=1,...,N .
• Thus, it appears that the best one can do is to design a learning algorithm
A which minimizes the empirical risk (“training error”), that is, upon input
S the learning algorithm should return the solution of the minimization
problem (7):
XN
A(S) = argmin 1/N L(h(ui ), yi ). (9)
h∈H i=1
In the next subsections we will see that the situation is more involved. First,
in the kind of complex real-world learning tasks that neural networks usually are
used for, algorithms that find exact minimal-training-error solutions do not exist.
One only can design learning algorithms that find approximate solutions. Second,
if one reduces the learning problem to finding minimal-training-error solutions, one
will almost always run into the problem of overfitting. This second complication
is by far more painful and important than the first one, and I will address it in
the following subsection.
11
The overfitting problem is connected to certain properties of learning algo-
rithms which I will collectively refer to as the flexibility of a learning algorithm.
The flexibility of a learning algorithm can be defined in several ways, and there
are several methods to steer the flexibility of a learning algorithm. Flexibility is
not a single, well-defined concept; it is an entire bundle of aspects which have
been addressed in machine learning in many ways. But the general idea is always
the same and can be stated in intuitive terms as “a learning algorithm A is more
flexible than another learning algorithm B if A can fit its computed models fˆ more
closely to the training data than B can do”. More flexible learning algorithms
thus can compute models with lower training error. Maximizing flexibility (by de-
signing learning algorithms that can make hopt (ui ) come very close to the teacher
outputs yi ) can lead to doing too much of a good thing: a super flexible learning
algorithm may even give zero training error, while performing very poorly on new
testing data, rendering the found model fˆ absolutely useless.
Because overfitting is such a fundamental challenge in supervised machine
learning, I illustrate its manifestations with four examples. They are copied with
slight adaptation of notation from the machine learning lecture notes.
p(u) = w0 + w1 u + · · · + wk uk . (10)
12
!
Figure 1: An example of training data (red squares) obtained from a noisy ob-
servation of an underlying “correct” function sin(2 π u) (broken blue
line).
At this moment we don’t bother how this task is solved computationally but simply
rely on the Matlab function polyfit which does exactly this job for us: given data
points (ui , yi ) and polynomial order k, find the coefficients wj which minimize this
MSE. Figure 2 shows the polynomials found in this way for k = 1, 3, 10.
If we compute the MSE’s for the three orders k = 1, 3, 10, we get MSEtrain =
0.4852, 0.0703, 0.0000 respectively. Some observations:
• For k = 3, we get a polynomial that hits our target sine apparently quite
well.
13
• For k = 10, we get a polynomial that perfectly matches the training data,
but apparently misses the target sine function (overfitting).
1.5 6
150
1 4
1
100
1 1
0.5 2
50
0 0.5
0.5 0 0.5
1 0
1
1
0.5 0.5 0.5
0 0 0 0 0 0
Figure 3: Estimating a pdf from 6 data points. Model flexibility grows from left
to right. Note the different scalings of the z-axis: the integral of the pdf
(= the volume under the surface) is 1 in each of the three cases.
Let us consider the task of estimating a 2-dimensional pdf over the unit square
from 6 given training data points {ui }i=1,...,6 , where each ui is in [0, 1]×[0, 1]. This
is an elementary unsupervised learning task, the likes of which frequently occur
14
as a subtask in more involved learning tasks, but which is also of interest in its
own right. Figure 3 shows three pdfs which were obtained from three different
learning runs with models of increasing flexibility (I don’t explain the learning
algorithms here — for the ones who know about it: simple Gaussian Parzen-
window models where the degree of flexibility was tuned by kernel width). Again
we witness the fingerprints of under/overfitting: the low-flexibility model seems too
“unbending” to resolve any structure in the training point cloud (underfitting),
the high-flexibility model is so volatile that it can accomodate each individual
training point (presumably overfitting).
But again, we don’t really know...
15
Figure 5: A nightmare case of overfitting. Picture spotted
by Yasin Cibuk (2018 ML course participant) on
[Link] (link now dead),
designed and crafted by artist Dominic Wilcox, 1999. Quote from the
artist’s description of this object: “I used my own body as a template
for the mattress”. From a ML point of view, this means a size N = 1
training data set.
• On the left end (very inflexible learning algorithms), both training and test-
ing error are high – due to the low flexibility of the algorithm, the obtained
models are just too simple to capture relevant structure in the data distri-
bution.
• On the right end (super flexible learning algorithms), the training error is
small — often, indeed, it can be pushed down to zero — because the high
flexibility of the learning algorithm allowed it to fit itself around each in-
dividual training point. In human psychology terms: the model has learnt
the data by heart. But, the testing error is high: memorizing teaching ma-
terial by heart is stupid; a good learner extracts the underlying regularities
16
and laws from the data and can transfer that valuable, extracted, condensed
knowledge to apply in new (testing) situations.
risk /
empirical
risk
risk,
“testing error”
Empirical risk,
“training error”
Figure 6: The generic, universal, core challenge of machine learning: finding the
right model flexibility which gives the minimal risk.
In order to practically find the sweet spot of optimal flexibility, two technical
conditions must be satisfied:
1. One must have an effective method to change the flexibility of learning al-
gorithms.
2. One must have an effective method to estimate the risk (red curve in the
figure).
If one has these two mechanisms available, one can find the sweet spot (green
line in our figure) by a systematic sweep through flexibilities, learning models for
each flexibility, estimate the risk, and settle for the flexibility that has minimal
estimated risk. In the next two subsections I outline how these two mechanisms
can be instantiated.
17
1.4.1 Tuning learning flexibility through model class size
In the polynomial curve fitting example from Section 1.3.1, the model parameters
were the monomial coefficients w0 , . . . , wk (compare Equation 10). After fixing the
polynomial order k, the polynomial p(u) Pkwith minimal training error was selected
from the set Hk = {p : R → R | p(x) = j=0 wj u }, that is, the learning algorithm
j
It is clear that H1 ⊂ H2 ⊂ . . .. It is also clear that the training error can only
shrink when k grows, because the set of candidate solutions Hk grows with k, thus
there are more candidate solutions to pick from for the optimization algorithm.
Generalizing from this example, we can see one way to obtain a sequence of
learning algorithms of increasing flexibility: A model class inclusion sequence is a
sequence H1 ⊂ H2 ⊂ . . . of sets of candidate models. If one has an algorithm A
that can solve the optimization problem (9) for each k, and if one can restrict the
search space for this algorithm to Hk by a control input k, then the flexibility of
this algorithm can be tuned by this input k. For example, the polyfit algorithm
supplied by Matlab can (in fact, must) be given an argument that specifies the
polynomial order that the algorithm will consider to find the best polynomial fit.
There are many ways how one can set up a sequence of learning algorithms
which pick their respective optimal models from such a model class inclusion se-
quence. In most cases this will just mean to admit larger models with more tune-
able parameters for model classes higher up in the sequence. In polynomial curve
fitting this meant to admit polynomials of increasing order. In neural network
training this means to run neural network training algorithms on neural networks
of increasing size.
18
and size of the considered model class should be rich and large enough to be able
to overfit (!) the available training data. Thus one can be sure that the “right”
model is contained in the search space H.
The models in H are typically characterized by a set of trainable parameters.
In polynomial curve fitting these parameters are the monomial coefficients, and
for a fixed neural network structure, it would be the set of all synaptic weights.
Following the traditional notation in the machine learning literature we denote
this collection of trainable parameters by θ. This is a vector that has as many
components as there are trainable parameters in the chosen kind of model. We
assume that we have D tuneable parameters, that is θ ∈ RD .
Such a high-flexibility model type would inevitably lead to overfitting when
an “optimal” model would be learnt using the basic learning equation (7) which I
repeat here for convenience:
X
N
fˆ = hopt = argmin 1/N L(h(ui ), yi ).
h∈H i=1
X
N
hopt = argmin 1/N L(h(ui ), yi ) + α2 reg(θh ). (12)
h∈H i=1
1. In the polynomial fit task from Section 1.3.1 one might consider for H all
10th order polynomials, but penalize “oscillations” like the ones seen in the
right panel of Figure 2, that is, penalize such 10th order polynomials that
exhibit strong oscillations. The degree of “oscillativity” can be measured,
for instance, by the integral over the (square of the) second derivative of the
polynomial p,
Z 1 2
′ d2 p(u)
reg(θ) = reg((w0 , . . . , w10 ) ) = du.
0 du2
19
Investing a little calculus and linear algebra (good exercise! not too difficult),
it can be seen that this integral resolves to a quadratic form reg(θ) = θ′ C θ
where C is an 11 × 11 sized positive semi-definite matrix. That format is
more convenient to use than the original integral version.
2. A popular regularizer that often works well is just the squared sum of all
model parameters, X
reg(θ) = w2 .
w∈θ
This regularizer favors models with small absolute parameters, which often
amounts to “geometrically soft” models. This regularizer is popular among
other reasons because it supports simple algorithmic solutions for minimiz-
ing risk functions that contain it. It is called the L2 -norm regularizer (or
simply the L2 regularizer) because it penalizes the (squared) L2 -norm of
the parameter vector θ. In the deep learning literature and its fabulous
power-toolboxes, this regularizer is also known as weight decay.
Computing a solution to the minimization task (12) means to find a set of pa-
rameters which simultaneously minimizes the original risk and the penalty term.
The factor α2 in (12) controls how strongly one wishes the regularization to
“soften” the solution. Increasing α2 means downregulating the model flexibility.
For α2 = 0 one returns to the original un-regularized empirical risk (which would
likely mean overfitting). For very large α2 → ∞ the regularization term dominates
the model optimization and one gets a model which does not care anymore about
the training data but instead only is tuned to have minimal regularization penalty.
In case of the L2 norm regularizer this means that all model parameters are zero
– the ultimate wiggle-free model; one should indeed say the model is dead.
When regularization is used to steer the degree of model flexibility, the x-axis
in Figure 6 would be labelled by α2 (highest α2 on the left, lowest at the right end
of the x-axis).
Using regularizers to vary model flexibility is often computationally more con-
venient than using different model sizes, because one does not have to tamper
with differently structured models. One selects a model type with a very large
(unregularized) flexibility, which typically means to select a big model with many
parameters (maybe a neural network with hundreds of thousands of synaptic con-
nections).
20
(ui , yi )i=1,...,N by adding some noise vectors to each input pattern ui : for
each ui , add l variants ui + νi1 , . . . , ui + νil of this pattern to the train-
ing data, where the νij are i.i.d. random vectors (for instance, uniform or
Gaussian noise). This increases the number of training patterns from N to
(l +1)N . The more such noisy variants are added and the stronger the noise,
the more difficult will it be for the learning algorithm to fit all data points,
and the smoother the optimal solution becomes – that is, the more one steers
to the left (underfitting) side of Figure 6. Adding noise to training patterns
is a very common strategy. It is one way of doing data augmentation: create
artificial, “distorted” versions of the training data and add them to the data
used for training.
21
there is no special word for the subset T . It is also called “training set” in the
literature, which is confusing because it is a subset of the originally given complete
training set. I will call T the reduced training set.
A bit more formally, let the flexibility axis in Figure 6 be parametrized by some
appropriate flexibility measure r, where small r means “strong regularization”,
that is, “go left” on the flexibility axis. Let the range of r be r = 1, . . . , l.
For each setting of flexibility r, the data in T is used to train an optimal
model hopt r . The test generalization performance on “new” data is then tested
on the validation set. It is determined which model hopt ropt performs best on
the validation data, that is, has the least average loss over the validation data.
Its flexibility ropt is then taken to be the flexibility that marks the “sweet spot”
indicated by the green bar in Figure 6. After this screening of degrees of flexibility
for the best test performance, a model within the found optimal flexibility ropt is
then finally trained on the original complete training data set S = T ∪ V .
This whole procedure is called cross-validation. Notice that nothing has been
said so far about how to split S into T and V .
A clever way to answer this question is to split S into many subsets Sj of
roughly equal size (j = 1, ..., k). Then, for each flexibility value r carry out k
complete screening runs via cross validation, where in the j-th run the subset Sj
is withheld as a validation set, and the remaining k − 1 sets joined together make
for the reduced training set. After these k runs, average the validation errors
obtained in the k runs. This average is then taken as an estimate for the risk of
learning a model with flexibility r. This is called k-fold cross-validation. Here is
the procedure in detail:
22
Given: A set (ui , yi )i=1,...,N of training data, and a loss function L.
Also given: Some method which allows one to steer the model learning
through different degrees r of flexibility. The weakest regularization
should be weak enough to allow overfitting.
Step 1. Split the training data into k disjoint subsets Sj = (ui , yi )i∈Ij of
roughly equal size N ′ = N/k.
Step 2. Repeat for r = 1, . . . , l:
Step 2.1 Repeat for j = 1, . . . , k:
Step 2.2.1 Designate Sj as validation set Vj and the union of the
other Sj ′ as reduced training set Tj .
Step 2.2.1 Compute the model with minimal training error on
data input Tj
X
hopt r j = argmin 1/|Tj | L(h(ui ), yi ),
h∈Hr
(ui ,yi )∈Tj
Step 2.2 Average the k validation risks Rrvalj obtained from the k “folds”
carried out for this r, obtaining
X
Rrval = 1/k Rrvalj .
j=1,...,k
Step 3. Find the optimal flexibility by looking for that r which minimizes the
averaged validation risk:
Step 4. Compute hropt using the complete original training data set:
X
hropt = argmin 1/N L(h(ui ), yi ).
h∈Hropt
i=1,...,N
This procedure contains two nested loops and looks expensive. For economy,
one starts with the low-end r and increases it stepwise, assessing the generaliza-
23
tion quality through cross-validation for each regularization strength r, until the
validation risk starts to rise. The strength ropt reached at that point is likely to
be about the right one.
The best assessment of the optimal class is achieved when the original training
data set is split into singleton subsets — that is, each Sj contains just a single
training example. This is called leave-one-out cross-validation. It looks like a
horribly expensive procedure, but yet it may be advisable when one has only a
small training data set, which incurs a particularly large danger of ending up with
poorly generalizing models when a wrong model flexibility were used.
K-fold cross validation is widely used. It is a factual standard procedure in
supervised learning tasks when the computational cost of learning a model is
affordable.
I should also mention that cross-validation is not only used for finding the right
degree of regularization, but can similarly be used for tuning hyperparameters of a
learning procedure A. The term “hyperparameter” is generally used for all kinds of
“knobs to play around with” when tuning a complex learning algorithm. Modern
deep learning algorithms have frightfully many tuning options.
24
1.7 How to deal with testing and validation in “real life”
and in your semester projects
There are fine and important distinctions between the concepts of “testing” and
“validation”. A professional machine learning expert is highly aware of these dis-
tinctions and handles and names them always correctly (otherwise he/she isn’t
an expert). But beginners in the game often enter and leave the battlefield of
model optimization and testing in a confused state of mind — no wonder be-
cause the best and correct ways how to handle “testing” and “validation” come in
context-dependent variations. Here is a summary:
25
his/her ultimate model to a final “real-world” test by running it on the offi-
cial testing data Stest . Importantly, these test data Stest may only be visited
once, namely for the final performance assessment of the optimized model
— in no way may data from Stest be used during the model optimization
process.
The majority of free candidate datasets for your semester project will just
be that, plain datasets S, without a predefined train/test split. If you want
to use such a dataset, there are two main ways to proceed:
• If you want to play it by the textbook gold standard method, you first
split S into S = Strain ∪ Stest , creating your own “bemchmark” kind
of dataset. By what percentage you do this split is your decision - a
standard split is an 80:20 percent ratio. You put Stest away in a secret
vault. Then do the obvious: you split Strain further into a reduced
training set and a validation set Strain = T ∪ V , optimize your model
by cross-validation, train a final best model on all of Strain , then unlock
the vault, put your final best model to the final cruel test on the Stest
data, and report the performance metrics that you obtain.
• A lazy man’s shortcut procedure, which is only the second best way
to proceed but which is often found even in the reports of respectable
expert authors, skips the split S = Strain ∪ Stest . No testing data are
defined or used. Instead, the available dataset S is directly split into
a reduced training set T and a validation set V . Then cross-validation
experimentation is carried out as in all the scenarios before in order
to find a training scheme that gives the best validation performance.
When that optimal scheme is found, the validation performance is re-
ported as the final result of the project - work done. In this shortcut
procedure, no final best model needs to be trained on all of S - it would
be pointless because there would be no data to test it on.
26
Rosenblatt (1958)). It is a computational neural architecture inspired by what
was known about the human visual processing system at that time, and its work-
ings were demonstrated on visual character classification tasks. Figure 7 gives a
diagram of the Perceptron’s architecture.
u1
&$ A „neuron“ in today‘s
NN terminology
u2
&)
%
! &" '"
uj "0"
&" s
"#$
"1"
uK-1
&%($
"Response" unit:
uK
&% • Computes binary function s of
%
input
! &" '"
Input "retina" Trainable connection "#$
weigths • Typical choice: threshold unit
27
0-1 valued function on them. In modern machine learning language, these
associator units compute local random features from the retinal image.
• The activation values (denoted uj in Figure 7) of the associator units are
propagated forward into a single “summation unit”, where they are summed
up — but not before they have been scaled by “synaptic” weights wj . Thus,
the
P summation unit together with these weights yields a linear combination
j wj uj of the associator activation values.
• Finally, the output value of the summation unit is passed through a “re-
sponse unit”, which delivers the ultimate output of the Perceptron by ap-
plying a thresholding operation which results in a 0 or 1 value.
• The entire processing pipeline of the Perceptron thus transforms retinal sen-
sor patterns into a binary output — which is the format of an image classi-
fication system. The Perceptron was thus probed on a binary classification
task, for instance letters “A” versus letters “C” (= output 0 versus 1).
The Perceptron was revolutionary in that it was a learning system. In order
to solve the requested image classification task, the synaptic weights must be the
right ones. It is not difficult to find a mathematical formula which computes a set
of well-working synaptic weights if such a set exists at all. Rozenblatt’s perceptron
however “learnt” the right synaptic weights in an incremental, iterative training
process, which gives this model some biological plausibility. In this process, the
network was presented with a series of what we have learnt to call training pat-
terns (uraw raw
i , yi )i=1,...,N , where the ui are vectors containing the response values of
the photocells, and yi ∈ {0, 1} is a binary indicator, for instance yi = 0 if the input
pattern came from an “A” image and yi = 1 if the input came from a “C” image.
When during the training process the Perceptron was shown some input and its
output was correct, the synaptic weights remained unchanged; if the output was
wrong, they were adapted by an incremental weight adaptation rule, the Percep-
tron learning rule. Rosenblatt could prove that this learning rule would converge
to a perfectly performing network if a perfectly working set weights existed at all.
The Perceptron (in its Mark 1 realization) was a remarkable mix of analog
signal processing hardware, electrical engineering, biological inspiration, and ele-
mentary mathematics. It featured many innovations that today are commonplace
for neural networks in machine learning (random nonlinear features, computing
outputs by passing the results of a linear combination through a nonlinearity,
using image classification as a benchmark demonstration).
In plain maths, if one strips off all biologism and does not include the random
Boolean operations of the associator units in the learning story (because they
remain fixed and are not trainable), the Perceptron simply learns a function
fˆ : RK → {0, 1}
P
(u1 , . . . , uK )′ 7→ σ( Kj=1 wj uj ), (13)
28
A B C
Cycle through training data and adapt weights: Present the training pat-
terns to the Perceptron one after the other, starting again from the beginning
when the last pattern has been used. At each presentation of a pattern ui ,
Stop when in one full cycle through all training examples no weight changes were
triggered (that is, all training examples are correctly classified).
29
Today, even a beginner in machine learning would immediately criticize that
this method only minimizes the training error, inviting overfitting; and that this
seems a very weak and restricted learning algorithm, because it is only applicable
to 2-class classification problems, and the search space H for possible solutions fˆ
(see Equation 13) contains only very simple functions, namely affine linear maps
followed by a thresholding (an affine linear map is a linear map plus some optional
shift – the bias).
But, at the time when Rosenblatt presented the Perceptron, machine learning
as a field did not exist at all (the Perceptron can be seen as one of the starters).
The Perceptron bundled an entire collection of 100% new and ingenious ideas, and
it looked like a little brain. Here is a sniplet from the New York Times (“New
Navy Device Learns by Doing”, NYT July 8, 1958), after a press conference held
by Rosenblatt at the US Office of Naval Research on July 7, 1958 (cited after
Olazaran (1996)):
“The Navy revealed the embryo of an electronic computer today that it expects
will be able to walk, talk, see, write, reproduce itself and be conscious of its existence.
Later perceptrons will be able to recognize people and call out their names and
instantly translate speech [...], it was predicted.”
Before you feel immensely superior to the 1958 authors and readership of the
New York Times, think about what you, today, think about the powers of deep
learning. Brain-like, eh? ... let’s wait a few years ...
The hype about the Perceptron did not last long. In a book titled “Percep-
trons”, famous AI pioneers Minsky and Papert (1969) pointed out the obvious:
with an architecture that boils down to a simple, thresholded linear combination
of input values, one can only learn to classify patterns that are linearly separa-
ble. That is, in the input space RK where patterns come from, there must exist
an n − 1 dimensional linear hyperplane, dividing RK into two half-spaces, such
that all patterns of class 0 lie on one side of the hyperplane and all class-1 pat-
terns on the other side. Specifically, Minsky and Papert shone a flashlight on the
embarassing fact that Perceptrons cannot learn the XOR function, because this
function is not linearly separable (see Figure 9). As a consequence of this simple
insight (“perceptrons can’t even learn the simplest Boolean functions”), neural
networks become disreputable as an approach to realize machine intelligence, and
research on neural networks came to a dead halt and remained frozen — until
things happened in the early 1980’s, which will be related in the next section.
I must add that Rosenblatt was very well aware of this limitation of the ele-
mentary Perceptron. In a summary documentation of his work on neural learning
(Rosenblatt, 1961), he described extended versions of the Perceptron which had
several layers of neurons and which today we would call multi-layer perceptrons
(next subsection). He also envisioned learning algorithms for these. What a pio-
neer! Yet, the mutilated message that diffused from the Minsky and Papert book
into the general perception of the young AI community was that “perceptrons
can’t even learn XOR”. History is a beast that does not know fair play.
30
u1 u2 XOR(u1, u2) u2
0 0 0 1
0 1 1
1 0 1
1 u1
1 1 0
Figure 9: The XOR stopper for early neural network research. The training
dataset consists of merely 4 examples (u1 , y1 ) = ((0, 0), 0), (u2 , y2 ) =
((0, 1), 1), (u3 , y3 ) = ((1, 0), 1), (u4 , y4 ) = ((1, 1), 0), but classifying the
(0, 0), (1, 1) patterns in one class and the (0, 1), (1, 0) patterns into an-
other is impossible for a Perceptron, because there is no way to draw a
line (like the two thin blue ones shown) in the pattern space R2 which
separates the two classes, that is, all green points lie on one side of the
line and all red ones on the other (right panel).
31
ton and Terrence Sejnowski (Hinton and Sejnowski, 1983) (subsequently
much more clearly explained in Ackley et al. (1985)), is a universal, unsu-
pervised neural learning system which can, in principle, learn every kind of
distribution from training data. It is mathematically very transparent but
computationally extraordinarily expensive and thus of limited practical use.
Besides producing many other groundbraking innovations to neural network
science, Hinton continued to think about the Boltzmann machine and finally
found a way to cut down the computational cost in the Restricted Boltzmann
Machine (RBM), presented to a wider scientific community in an article in
Science (Hinton and Salakuthdinov, 2006). In passing, in this article it is
mentioned how RBMs can be used to initialize the training of feedforward
neural networks for pattern recognition — in retrospect, these inconspicuous
remarks turned out to be the starting shot for deep learning. I will devote
one session to Boltzmann machines in this course.
Besides these three neural network models, many others were conceived in the
1970 / early 1980 years, often inspired by physics, neuroscience, or cognitive science
(not inspired by AI, since that field was too successful in its own non-neural ways
in those days and too proud to reconsider the disreputed neural networks). These
innovative neural network models helped to understand principles of biological
and cognitive neural information processing.
But none of them was really practically useful.
The tides turned in the year 1986, when a two-volume collection of neural
network articles was published, “Parallel Distributed Processing” (Rumelhart and
McClelland, 1986), of which the first volume soon became known (and still is) as
the PDP book. This volume contained a number of well-written, detailed intro-
ductions to the innovative NN approaches that had been modestly but steadily
sprouting in the shadows outside the glaring light of AI research. Suddenly, neural
networks were back on stage with a fanfare.
The immense and lasting impact of the PDP book is largely due to only one of
its chapters. In this chapter, Rumelhart et al. (1986) give a transparent introduc-
tion how neural networks could be designed and trained for pattern classification
tasks in ways that are much more powerful than what could be achieved with the
Perceptron. The dreaded XOR learning task suddenly was a piece of cake. The
added powers of the multi-layer Perceptrons (MLPs), as the neural networks de-
scribed in this chapter became generally called, arose from two innovations. First,
while the classical Perceptron only has a single “summation-response” unit, in
MLPs many of such units are coupled together in parallel and serial connectivity
patterns. Second, the original Perceptron learning rule (which can’t be applied
to MLPs) was replaced by the error backpropagation (or simply “backprop”) al-
gorithm, which made it possible to iteratively adapt the synaptic weights in an
MLP in order to minimize a given loss function — for instance, in order to classify
images. It must be mentioned that the backpropagation algorithm was not freshly
invented by the authors of the PDP book chapter, but had been described by Paul
32
Werbos much earlier in his 1974 PhD thesis — which did not attract attention and
remained ignored; and in other fields of engineering, similar algorithms had been
discovered even earlier. Schmidhuber (2015), Section 5.5, provides a historical
overview.
At any rate, after 1986 neural networks lifted off again. And again, there was
a hype that soon withered (I will relate this in the lecture). But the decline of
that new hype wasn’t as dramatic as what happened to Perceptrons. In the two
decades from 1986 to 2006, MLPs remained a respected citizen of the world of
AI and machine learning (this field now had come into being). But not more
than that — other methods, like support vector machines or Bayesian networks,
were more prominent in machine learning. MLPs worked well as far as it went,
but that wasn’t too far after all: roughly speaking, they worked fairly well in not
too nonlinear signal processing and pattern recognition tasks. Their economical
impact was not striking.
After this short historical review (I find history fascinating), I will describe
MLPs in sober detail. Much of the material in this section is taken from my
machine learning lecture notes.
MLPs are used for the supervised learning of vectorial input-output tasks,
based on training samples S = (ui , yi )i=1,...,N , where u ∈ RK , y ∈ RM are drawn
from a joint distribution PU,Y .
The MLP is trained to produce outputs y ∈ RM upon inputs u ∈ RK in a way
that this input-output mapping is similar to the relationships ui 7→ yi found in
the training data. Similarity is measured by a suitable loss function. Several loss
functions are today standardly used.
An MLP is a neural network equipped with K input units and M output
units. A K-dimensional input pattern u can be sent to the input units, then
the MLP does some interesting internal processing, at the end of which the M -
dimensional result vector of the computation can be read from the M output units.
An MLP N with K input units and M output units thus instantiates a function
N : RK → RM .
Just like how it worked with the Perceptron, an MLP is defined by its structure
(often called architecture) and by the values of its “synaptic connection” weights.
The architecture is invented and fixed by the human experimenter / data engineer
in a way that should be matched to the given learning task. The architecture is
(mostly) not changed during the subsequent training. MLPs have a more complex
architecture than the Perceptron. A given neuron in an MLP is typically receiving
input from, and sending its output to, many other neurons along synaptic connec-
tions. Like in the Perceptron, each such connection is characterized by a weight.
These weights are iteratively and incrementally adapted in the training process
until the network function N : RK → RM comes close to what one desired.
It is customary in machine learning to lump all trainable parameters of a
machine learning model together in one parameter vector, which is standardly
denoted by θ. For a neural network, θ is thus the vector of all the trainable
33
synaptic connection weights. If there are D such weights, θ ∈ RD . Since the
network function is determined by θ, one also writes Nθ : RK → RM if one wishes
to emphasize the dependance of N ’s functionality on its weights.
The learning task is defined by a loss function L : RM × RM → R≥0 . A
convenient and sometimes adequate choice for L is the quadratic loss L(N (u), y) =
kN (u) − yk2 , but other loss functions are also widely used. Chapter 6.2 in the
deep learning bible Goodfellow et al. (2016) gives an introduction to the theory of
which loss functions should be used in which task settings.
Given the loss function, the goal of training an MLP is to find a weight vector
θopt which minimizes the risk, that is
1 X
N
θopt = argmin L(Nθ (ui ), yi ). (15)
θ∈Θ N i=1
34
Time series prediction: inputs are vector encodings of a past history of a temporal
process, outputs are vector encodings of future observations of the process.
Denoising, restoration and pattern completion: inputs are patterns that are cor-
rupted by noise or other distortions, outputs are cleaned-up or repaired or
completed versions of the same patterns.
Data compression: Inputs are high-dimensional patterns, outputs are low-dimensional
encodings which can be restored to the original patterns using a decoding
MLP. The encoding and decoding MLPs are trained together.
Process control: In control tasks the objective is to send control inputs to a tech-
nological system (called “plant” in control engineering) such that the system
performs in a desired way. The algorithm which computes the control inputs
is called a “controller”. Control tasks range in difficulty from almost trivial
(like controlling a heater valve such that the room temperature is steered to
a desired value) to almost impossible (like operating hundreds of valves and
heaters and coolers and whatnots in a chemical factory such that the chemi-
cal production process is regulated to optimal quality and yield). The MLP
instantiates the controller. Its inputs are settings for the desired plant be-
havior, plus optionally observation data from the current plant performance.
The outputs are the control inputs which are sent to the plant.
35
x1k x2 k ... xLkk
output
neurons
wijk
last hidden
1 x1k-1 x2 k-1 ... xL k-1k-1 layer of
neurons
first hidden
1 x11 x21 ... xL11 layer of
neurons
wij1
2. For κ < k, assume that the activations xκ−1 j of units in layer κ − 1 have
already been computed (or have been externally set to the input values, in
the case of κ − 1 = 0). Then the activation xκi is computed from the formula
!
X
Lκ−1
κ κ−1
xκi = σ wij κ
xj + wi0 . (16)
j=1
That is, xκi is obtained from linearly combining the activations of the lower
layer with combination weights wij κ κ
, then adding the bias wi0 ∈ R, then
wrapping the obtained sum with the activation function σ. The activation
function is a nonlinear, “S-shaped” function which I explain in more detail
κ
below. It is customary to interpret the bias wi0 as the weight of a synaptic
link from a special bias unit in layer κ − 1 which always has a constant
activation of 1 (as shown in Figure 10).
Equation 16 can be more conveniently written in matrix form, and most
conveniently in a version that treats the bias weights as synaptic weights
from bias units. Let xκ = (xκ1 , . . . , xκLκ )′ be the activation vector in layer
κ, let bκ = (w10
κ
, . . . , wLκ κ 0 )′ be the vector of bias weights, and let Wκ =
36
κ
(wij )i=1,...,Lκ ; j=0,...,Lκ−1 be the connection weight matrix for links between
layers κ−1 and κ (including the bias weights, which become the first column
in W). Then (16) becomes
xκ = σ Wκ [1; xκ−1 ] , (17)
where [1; xκ−1 ] is the vector xκ−1 prepended with a leading 1 and the acti-
vation function σ is applied component-wise to the vector Wκ xκ−1 .
y = xk = Wk xk−1 , (18)
that is, in the same way as it was done in the other layers except that no
activation function is applied. The output activation vector y is the result
y = N (u).
In the area of deep learning a drastically simplified “sigmoid” is often used, the
rectifier function defined by r(a) = 0 for a < 0 and r(a) = a for a ≥ 0. The rectifier
has somewhat less pleasing mathematical properties compared to the classical
sigmoids but can be computed much more cheaply. This is of great value in deep
learning scenarios where the neural networks and the training samples both are
often very large and the training process requires very many evaluations of the
sigmoid.
In intuitive terms, the operation of an MLP can be summarized as follows.
After an input vector u is written into the input units, a wave of activation sweeps
forward through the layers of the network. The activation vector xκ in each
37
1
0.5
ï0.5
ï1
ï3 ï2 ï1 0 1 2 3
Figure 11: The tanh (blue), the logistic sigmoid (green), and the rectifier function
(red).
layer κ is directly triggered by the activations xκ−1 according to (17). The data
transformation from xκ−1 to xκ is a relatively “mild” one: just an affine linear map
Wκ xκ−1 followed by a wrapping with the gentle sigmoid σ. But when several such
mild transformations are applied in sequence, very complex “foldings” of the input
vector u can be effected. Figure 12 gives a visual impression of what a sequence
of mild transformations can do.
MLPs are a member of a larger class of neural networks, the feedforward neural
networks (FFNs). FFNs all have an input layer, an output layer, and in between a
directed neural connection network which lets the input activation vector spread
forward through the network, undergoing all sorts of transformations, until it
reaches the output layer. Importantly, there is never a synaptic connection inside
the network that would feed some intermediate activation back into earlier layers.
Mathematically speaking, the connection graph of an FNN is cycle free (also called
acyclical). All FFNs are therefore representing a function from input vectors to
output vectors. Within the class of FNNs, the MLPs are the most simple ones,
with all layers having the same basic structure. For a given task, this generic simple
MLP architecture may be suboptimal (it very likely is!). Much more sophisticated
feedforward architectures have been developed for specific tasks. At the end of
this section I will highlight one of these, the convolutional neural networks whose
architecture is optimized for image recognition.
It is also possible to design neural networks that host cyclical synaptic connec-
tion pathways. In such recurrent neural networks (RNNs), a neuron’s activation
can ultimately feed back on the same neuron. Mathematically speaking, RNNs are
not representing functions but dynamical systems. This is an upper league sector
of mathematics, far more complex, fascinating, and powerful than just functions.
Biological brains are always recurrent, and intelligent reasoning has memory —
38
methods/[Link]
[Link]
Figure 12: Illustrating the power of iterating a simple transformation. The baker
transformation (also known as horseshoe transformation) takes a 2-
dimensional rectangle, stretches it and folds it back onto itself. The
bottom right diagram visualizes a set that is obtained after numer-
ous baker transformations (plus some mild nonlinear distortion). —
Diagrams on the right taken from Savi (2016).
which means it feeds back on itself recurrently. Most of this course will be about
RNNs! Only this Section 2 is about FNNs.
39
However, these constructions have little practical value because the constructed
MLPs N are far too large for practical implementation. You can find more de-
tails concerning such approximation theorems and related results in my legacy
ML lecture notes [Link]
Section 8.1, and a very instructive online mini-tutorial on how these function-
approximation networks can be constructed at [Link]
com/[Link] (pointed out to me by Satchit Chatterji).
Even when the function f that one wants to train into an MLP is very complex
(highly nonlinear and with many “folds”), it can be in principle approximated with
1-hidden-layer MLPs. However, when one employs MLPs that have many hidden
layers, the required overall size of the MLP (quantified by total number of weights)
is dramatically reduced (Bengio and LeCun, 2007). Even for super-complex target
functions f (like photographic image caption generation), MLPs of feasible size
exist when enough layers are used. It is not uncommon for professionally designed
and trained deep networks to have 20 or more layers. This is the basic insight and
motivation to consider deep networks, which is just another word for “many hidden
layers”. Unfortunately it is not at all easy to train deep networks. Traditional
learning algorithms had made non-deep (“shallow”) MLPs popular since the 1980-
ies. But these shallow MLPs could only cope with relatively well-behaved and
simple learning tasks. Attempts to scale up to larger numbers of hidden layers
and more complex data sets largely failed, due to numerical instabilities, very
slow convergence, or poor model quality. Since about 2006 an accumulation of
clever “tricks of the trade” plus the availability of affordable powerful (GPU-based)
computing hardware has overcome these hurdles.
1. Get a clear idea of the formal nature of your learning task. Do you want
a model output that is a probability vector? or a binary decision? or a max-
imally precise transformation of the input? how should “precision” best be
40
measured? and so forth. Only proceed with using MLPs if they are really
looking like a suitable model class for your problem.
2. Decide on a loss function. Go for the simple quadratic loss if you want a
quick baseline solution but be prepared to invest in other loss functions if
you have enough time and knowledge (read Chapter 6.2 in Goodfellow et al.
(2016)).
5. Fix an MLP architecture. Decide how many hidden layers the MLP shall
have, how many units each layer shall have, what kind of sigmoid is used
and what kind of output function and loss function. The structure should be
rich enough that data overfitting becomes possible and your regularization
method can kick in.
8. Do the job. Enjoy the powers, and marvel at the wickedness, of MLPs.
41
You see that “neural network training” is a multi-faceted thing and requires
from you to consider all the issues that always jump at you in supervised ma-
chine learning. It will not miraculously give good results just because it’s “neural
networks inside”. The actual “learning” part, namely solving the optimization
task (14), is only a subtask, albeit a conspicuous one because it is done with an
algorithm that has risen to fame.
Given: Training data (ui , yi )i=1,...,N for a supervised learning task, as usual, where
ui ∈ RK , yi ∈ RM .
42
adaptation procedures which lead to a decreasing risk sequence. Examples are the
family of Expectation-Maximization algorithms (explained in my machine learn-
ing lecture notes if you are interested), or iteration schemes that exploit stability
conditions for fixed points of a map; or various sorts of general-purpose stochastic
search algorithms like genetic algorithms or simulated annealing (we will study
the latter in Section 6.4). Or, finally — gradient descent algorithms. That is the
kind of iterative model optimization algorithm that is typically used with neu-
ral networks (although, as we will see later in this course, by no means the only
one, and also most likely not the one used by Mother Nature to make our brains
learn). Gradient descent algorithms come in many degrees of sophistication. Here
I discuss the plain vanilla kind.
The gradient ∇R(θ) is the vector which points from θ in the direction of
the steepest ascent (“uphill”) of the performance surface. The negative gradient
−∇R(θ) is the direction of the steepest descent (Figure 14).
43
Figure 13: A (2-dimensional cross-section of) a performance surface for a neural
network. The performance landscape shows the variation of the loss
when (merely) 2 weights in the network are varied. Try to imagine the
loss function if it were “plotted” across not only two, but thousands or
millions of parameters! Source: [Link]
16/neural-network-loss-visualization/
The adaptation rate (or learning rate) µ is set to a small positive value.
An obvious weakness of this elegant and natural approach is that the final
model θ(∞) depends on the choice of the initial model θ(0) . In complex risk land-
scapes (as the one shown in Figure 13) there is no hope of guessing an initial
model which guarantees to end in the global minimum. This circumstance is gen-
erally perceived and accepted. There is a substantial mathematical literature that
amounts to “if the initial model is chosen with a good heuristic, the local minimum
that will be reached will be a rather good one with high probability”.
44
R
$2
$ −"R($)
$1
Figure 14: A performance surface for a 2-dimensional model family with parame-
ters θ1 , θ2 , with its contour plot at the bottom. For a model θ (yellow
star in contour plot) the negative gradient is shown as black solid ar-
row. It marks the direction of steepest descent (broken black arrow)
on the performance surface.
45
𝜃(0)
𝜃(1)
𝜃(2)
𝜃(∞)
Figure 15: A gradient descent itinerary, re-using the contour map from Figure 14
and starting from the initial point shown in that figure. Notice the
variable jump length and the sad fact that from this initial model θ(0)
the global minimum is missed. Instead, the itinerary slides toward a
local minimum at θ(∞) . The blue arrows show the negative gradient
at raster points. They are perpendicular to the contour lines and their
size is inversely proportional to the spacing between the contour lines.
46
oping the architecture of convolutional neural networks on his way (LeCun et al.,
1998), and being better able than anybody else to train these beasts by gradient
descent. However, the techniques that he applied were so subtle and not 100%
documented that other researchers were not generally able to reliably reproduce
his achievements.
I will now outline one of the most painful show stoppers for simpleminded
gradient descent.
First I note the obvious: differentiable performance surfaces are ususally not
linear but curved. This simple fact opens the doors for trouble.
In order to understand why curvature makes gradient descent difficult, I con-
sider a special situation where the main argument comes to the surface most
clearly. Assume that your gradient descent optimization already has brought you
into the close neighborhood of a local minimum, like the point θ(∞) indicated in
Figure 15. In the neighborhood of such a local minimum, the curved shape of
the performance surface R(θ) can be approximated by the first and second order
terms of the Taylor expansion of R(θ) around θ(∞) . Allow me to skip some maths
here (you find it detailed in the ML lecture notes, Section 11.3.4) which amounts
to a sequence of coordinate transformations. These coordinate transformations
simplify the picture to the scenario shown in Figure 16.
This graphics shows a contour plot of a performance surface over a two-
dimensional parameter space Θ = R2 after coordinate transformations that move
the local minimum θ(∞) to the origin and rotate the entire coordinate system such
that the main axes of the “trough” of the surface are aligned with the drawing
axes.
Assume that during a gradient descent run you have arrived at some parameter
vector θ(n) . In order to determine the next parameter vector θ(n+1) = θ(n) −
µ ∇R(θ(n) ), you have to decide on a learning rate µ. If you choose µ too large
(red arrow in Figure 16), you will jump forward along the direction of the negative
gradient so far that you end up “on the other side of the valley” higher than you
started — that is, the risk R(θ(n+1) ) at the next point is larger, not smaller, than
at θ(n) . If you would continue to use this too large learning rate even through
the next iterations, you would propel yourself upwards out of the “trough”. The
gradient descent algorithm has become instable.
In order to prevent this instability and divergence, the learning rate must be
set to a sufficiently small value. Then (indicated by the green arrow in the figure)
the gradient descent will continue moving downwards, stably reducing the risk at
every step.
The price to pay for stability is slow convergence. If the safe small learning rate
is kept for future iterations, the speed of approach to θ(∞) becomes very slow as
soon as the iterations reach the bottom of the “valley” (shaded dots in the figure).
But you cannot switch to a larger learning rate: the algorithm would immediately
become instable.
What does it mean, concretely, for a learning rate to be “too large” such that
47
0.3
0.2
0.1
!2
0
−0.1
!1
Figure 16: Gradient descent in the neighborhood of a local minimium θ(n) (cen-
tered at the origin of a shifted coordinate system). Red arrow: gradient
descent overshoots with too large learning rate. Green arrow: learning
rate small enough for stability. Magenta arrow: weight change direc-
tion when a second-order gradient descent method would be used. For
more explanation see text.
gradient descent becomes instable? The maths here are actually not very involved
(detailed in the machine learning lecture notes). Here I give a summary of the
main insights:
• After the coordinate transformations that were used to produce the situation
shown in Figure 16, the geometric shape of the performance landscape (in
its Taylor approximation up to order 2) is given by the formula
X
D
R(θ) = λi θi2 ,
i=1
where θ = (θ1 , . . . , θD )′ and the quadratic terms λi θi2 give the shape of the
parabolic cross-sections through the performance surface along the axes θi .
All λi are non-negative.
• In order to prevent instability, the learning rate must be set to a value smaller
than 1/λmax , where λmax is the largest among the coefficients λi .
48
• Given a stability-ensuring setting of the learning rate, once the iterations
θ(n) have reached the “bottom” of the “valley”, the final approach toward
the minimum value θ(∞) progresses geometrically, that is, the distance from
θ(n+1) to θ(∞) is smaller than the distance from θ(n) to θ(∞) by a constant
factor 0 < β < 1. This approach factor is equal to
λmax /λmin − 1
β= .
λmax /λmin + 1
The ratio λmax /λmin dictates the boundary conditions for successful gradient
descent optimization (close to a local minimum, and up to the order-two
Taylor approximation). If λmax = λmin one gets β = 0 and with an optimally
selected learning rate one achieves convergence in a single step. This is
possible if the performance surface around the minimum has a perfectly
circular shape. If conversely the ratio λmax /λmin is very large (very long and
sharp valley shape of the performance surface), β approaches 1 and gradient
descent becomes so slow that it amounts to a stillstand.
Bad news: the more layers a network has, the more common it is that the
ratio λmax /λmin is very large. In a study of convergence properties of gradient
descent minimization (for hyperparameter optimization of a recurrent neural net-
work (Jaeger et al., 2007)) I once numerically measured values for this ratio in the
order of 1e+14.
Besides this special case of convergence close to a minimum, there are other
geometric scenarios at other places on a performance surface (e.g., saddle points
of the performance landscape) which lead to similar inherent conflicts between
stability and speed of convergence.
One escape from the instability – slowness dilemma is to go for second-order
gradient descent methods. These methods determine the direction vector for the
next weigth adaption θ(n) → θ(n+1) based not only on the gradient, but also on
the curvature of the performance surface at the current model θ(n) . This curva-
ture is given by the D × D sized Hessian matrix which contains the second-order
partial derivatives ∂ 2 R /∂θi ∂θj . In ideal scenarios (when the second-order Taylor
approximation is precise), second-order gradient descent points into exactly the
direction to the target minimum (magenta arrow in Figure 16). However, comput-
ing the Hessian is expensive and may be subject to numerical problems, and the
preconditions for making second-order methods work well (good approximation of
surface by second-order Taylor expansion; being close to local minimum) may be
badly violated. Numerous variants of second-order methods have been developed,
aiming at reduced costs or more robust coping with local geometry. There is un-
fortunately no general rule of when or whether it is beneficial to use what type of
second-order method.
These difficulties are severe. For two decades they had made it impossible
to effectively train MLPs on advanced “cognitive” tasks in image interpretation
49
and language processing. Finding heuristic mechanisms to modify plain gradient
descent in ways that keeps the iterations stable, while maintaining a sufficient
speed of convergence, has been one of the enabling factors for deep learning.
Iterate: Compute a series (θ(n) )n=0,1,... of models of decreasing empirical loss (aka
training error) by gradient descent. Concretely, with
1 X
Remp (Nθ(n) ) = L(Nθ(n) (ui ), yi ) (22)
N i=1,...,N
In these iterations, make sure that the adaptation rate µ is small enough to
guarantee stability.
Stop when a stopping criterion chosen by you is met. This can be reaching a
maximum number of iterations, or the empirical risk decrease falling under
a predermined threshold, or some early stopping scheme.
50
even just BP if you want to have a feeling of belonging to a select community) is a
specific algorithmic scheme to compute this gradient in a computationally efficient
way.
Every student of machine learning must have understood it in detail at least
once in his/her life, even if later it’s just downloaded from a toolbox in some more
sophisticated fashioning. Thus, brace yourself and follow along!
Let us take a closer look at the empirical risk (22). Its gradient can be written
as a sum of gradients
!
1 X 1 X
∇Remp (Nθ ) = ∇ L(Nθ (ui ), yi ) = ∇L(Nθ (ui ), yi ),
N i=1,...,N N i=1,...,N
and this is also how it is actually computed: the gradient ∇L(Nθ (ui ), yi ) is evalu-
ated for each training example (ui , yi ) and the obtained N gradients are averaged.
This means that at every gradient descent iteration θ(n) → θ(n+1) , all training
data points have to be visited individually. In MLP parlance, such a sweep through
all data points is called an epoch. In the neural network literature one finds
statements like “the training was done for 120 epochs”, which means that 120
average gradients were computed, and for each of these computations, N gradients
for individual training example points (ui , yi ) were computed.
When training samples are large — as they should be — one epoch can clearly
be too expensive. Therefore one often takes resort to minibatch training, where
for each gradient descent iteration only a subset of the total training sample S is
used.
The backpropagation algorithm is a subroutine in the gradient descent game.
It is a particular algorithmic scheme for calculating the gradient ∇L(Nθ (ui ), yi )
for a single data point (ui , yi ). Naive highschool calculations of this quantity
incur a cost of O(D2 ) (where D is the number of network weights). When D is
not extremely small (it will almost never be extremely small — a few hundreds
of weights will be needed even for very simple tasks, and easily a billion for deep
networks), this cost O(D2 ) is too high for practical exploits (and it has to be paid
N times in a single gradient descent step!). The backprop algorithm is a clever
scheme for computing and storing certain auxiliary quantities which cuts down
the cost from O(D2 ) to O(D).
Here is how backprop works in order to compute the loss gradient ∇ L(Nθ (u), y)
for a training example (u, y).
1. BP works in two stages. In the first stage, called the forward pass, the
current network Nθ is presented with the input u and the output ŷ = Nθ (u)
is computed using the “forward” formulas (17) and (18). During this forward
pass, for each unit xκi which is not a bias unit and not an input unit the
quantity X
κ κ−1
aκi = wij xj (24)
j=0,...,Lκ−1
51
is computed and stored – this is sometimes referred to as the potential of
unit xκi , that is its internal state before it is passed through the sigmoid.
Define
∂L(Nθ (u), y)
δiκ = . (26)
∂aκi
Using (24) we find
∂aκi
κ
= xκ−1
j . (27)
∂wij
Combining (26) with (27) we get
∂L(Nθ (u), y)
κ
= δiκ xκ−1
j . (28)
∂wij
3. Computing the δ’s for output units. Output units xki are typically set up
differently from hidden units, and their corresponding δ values must be com-
puted in ways that depend on the special architecture. For concreteness here
I stick with the simple linear units introduced in (18). The potentials aki are
thus identical to the output values ŷi and we obtain
∂L(Nθ (u), y)
δik = . (29)
∂ ŷi
This quantity is thus just the partial derivative of the loss with respect to
the i-th output, which is usually simple to compute. For the quadratic loss
L(Nθ (u), y) = kNθ (u) − yk2 , for instance, we get
4. Computing the δ’s for hidden units. In order to compute δiκ for 1 ≤ κ < k
we again make use of the chain rule. We find
52
which is justified by the fact that the only path by which aκi can affect
L(Nθ (u), y) is through the potentials aκ+1
l of the next higher layer. If we
substitute (26) into (31) and observe (24) we get
X ∂aκ+1
δiκ = δlκ+1 l
∂aκi
l=1,...,Lκ+1
P
X ∂ κ+1
j=0,...,Lκ wlj σ(aκj )
= δlκ+1
∂aκi
l=1,...,Lκ+1
X ∂ wliκ+1 σ(aκi )
= δlκ+1
∂aκi
l=1,...,Lκ+1
X
= σ ′ (aκi ) δlκ+1 wliκ+1 . (32)
l=1,...,Lκ+1
This formula describes how the δiκ in a hidden layer can be computed by
“back-propagating” the δlκ+1 from the next higher layer. The formula can
be used to compute all δ’s, starting from the output layer (where (29) is used
— in the special case of a quadratic loss, Equation 30), and then working
backwards through the network in the backward pass of the algorithm.
When the logistic sigmoid σ(a) = 1/(1 + exp(−a)) is used, the computation
of the derivative σ ′ (aκi ) takes a particularly simple form, observing that for
this sigmoid it holds that σ ′ (a) = σ(a) (1 − σ(a)), which leads to
53
• Gradient-descent techniques on performance landscapes can only find a local
minimum of the risk function. This problem can be addressed by various
measures, all of which are computationally expensive. For deep networks
of large size the local minimum problem seems not particularly problem-
atic. These networks afford of such a brute overfitting potential that on the
downhill slide along the negative gradient, the point of overfitting and thus
“early stopping” is reached earlier than the local minimum one is heading
to. Or, in other words: most local minima represent overfitting models, thus
one does not want to reach them.
54
properly acknowledged by LeCun (LeCun et al., 1998) but today is rarely men-
tioned. The Neocognitron was not trained in a supervised way (the backpropo
algorithm was no public knowledge then) but with a biologically inspired unsu-
pervised learning algorithm. This condition, combined with the limited compute
powers of those days and the fact that Fukushima’s scientific goals were in biologi-
cal modeling (not machine learning), effected that the Neocognitron has remained
under-appreciated in the machine learning community. In turn, Fukushima drew
inspiration, and adopted some terminology, from Hubel and Wiesel’s Nobel prize
winning research on the early visual processing in the mammalian brain in the
early 1960’s.
• The first processing layer C1 after the retina is made of 6 feature maps. Each
feature map is a 2-dimensional array of feature detector neurons. Each such
feature detector neuron receives its input from a 5 × 5 sized receptive field in
the retina. Thus, each such neuron has 25 (plus 1 bias) weights that need to
be trained. Within each of the 6 feature maps, every feature detector neuron
however has the same weights as all other neurons in that feature map. This
weight sharing (aka weight tying) is imposed during the backpropagation
training. Mathematically, the transformation from the raw input image to a
55
feature map is a convolution of the image array with the 5 × 5 weight filter,
giving this kind of networks its name.
• The next layer C3 is again a convolutional layer. Each of its feature maps
(of which there are 16 in LeNet5) takes its input from several 5 × 5 sized
receptive fields, located in several of the 6 subsampling layers. In Figure 17,
I indicated one neuron in a C3 feature map receiving input from two 5 × 5
sized receptive fields.
• Then, like what was done from C1 to S2, these 16 feature maps are subsam-
pled into 16 subsampled S4 maps.
• the weight sharing reduces the number of trainable parameters, which leads
to a regularization effect in the sense of machine learning;
• the weight sharing also favors translational invariance for recognition per-
formance (for instance, a local feature which detects a sharp edge like at
the top of the letter A would become activated insensitive to shifts of the
position on the input letter A on the retina);
Convnets that were later developed in the DL era had more convolutional layers
and subsequent processing stages which are more complex than in their ancestor
LeNet5. Das (2017) is an online overview of later CNN architectures.
CNNs can also be designed for 1-dimensional “images”. This makes sense for
input data vectors u = (u1 , . . . , uK )′ where the ordering of vector elements has a
geometrical meaning. This would be expected, for instance, when the inputs u are
energy profiles over a frequency spectrum, or flow velocity measurements along a
56
flow sensor line as in a recent cum laude PhD thesis written in the AI department
(Wolf, 2020). In contrast, it would not make any sense whatsoever, for instance,
to throw a 1-dim CNN at inputs that are bag-of-words count vectors!
57
explaining how important and helpful it is to pre-process “raw” data for dimension
reduction, extracting features which are specifically appropriate for the task at
hand. Using raw, high-dimensional data is the way to failure, for a number of
good reasons that I unfold in those lecture notes.
So, — whom should you believe and follow?
Answer: both views are right — but in different sitations, not at the same
time. The critical issue is overfitting.
Flooding a deep learning system with high-dimensional raw data and training
it “end to end” will work if (and only if) one has large quantities of data. Together
with the clever modern regularization methods for deep networks (plus, possibly,
GAN methods), large data volumes minimize the dangers of overfitting. Such
convenient conditions (very large data sets, plus high-performance computing fa-
cilities and a lot of professional expertise) is what Google, Facebook and their likes
can enjoy.
However, the ordinary mortal machine learning professional most often will
only dispose of scarce training data, while nonetheless wanting to solve a complex
task. A typical non-Google setting is medical image processing. It is a notori-
ous problem in this field that there are never enough training examples. Yet,
the targeted input-output tasks (like cancer diagnostics from liver CT scans) are
eminently complex. Training a deep CNN “end to end” with a few hundreds of
unpreprocessed CT images will not work. Complex data processing pipelines, only
some of them being gradient-descent trained neural networks and others informed
by human insight for appropriate feature definition, will be necessary to come up
with good solutions (for example, Chlebus et al. (2018)).
2.4 A super reading for getting the full flavor of deep learn-
ing
If you want to get the fully brutal-gigantic picture of how DL feels from the inside
at the professional high end, I recommend the project report Le Scao et al. (2023).
This is a very well-written report on the BLOOM project. In this project, a multi-
lingual (40 languages, plus 11 programming languages!) large-scale foundation
model was designed, trained, tested, benchmarked against other large language
models, and documented — in a collaboration of more than 1000 volunteers from
many countries, with the support of French governmental agencies who donated
the compute infrastructure. The distinctive thing about BLOOM is that it is
open-source, that a full documentation of the training dataset is available, and
that it used novel methods to save energy (carbon footprint only 1/20 of what
GPT-3 needed). The report is long (74 pages) but well-structured and as easy
to read as can be, and it gives an overview of the complete project with all the
sub-activities and coordination efforts that were involved. The report has a large
number of co-authors — I did not count them but there are about 400 of them!
58
3 A short visit in the wonderland of dynamical
systems
Note: In previous issues of this course, I gave a very broad introduction to dy-
namical systems, covering both discrete-state and continuous-state systems. This
was too much. This time (Fall 2024) I will only present continuous-state systems,
because most neural networks are of that kind. These lecture notes are still de-
scribing both sorts, but I render in blue font all the parts that I will now skip. Read
these blue parts only if you are interested (it’s good stuff to know) — this material
will not be queried in exams.
1
4 x1 x2 –1
–4
We carry out three simulation runs where we iterate (33) for 1000 timesteps.
In the first run we use W1 , in the second run we use
13.2 −13.2
W2 = 3.3 ∗ W1 = ,
3.3 −3.3
59
in the third run we use W3 = 4.3∗W1 , plugging them into (33). Each run is started
from some arbitrary random initial state x(0). After a short initial transient where
the arbitrary initial state is “forgotten” (or “washed out”), the temporal sequence
x(n) stabilizes into a pattern which characteristically depends on W1 , W2 or W3 ,
respectively.
Figure 18 shows what we get:
• In the run using W1 , all activity dies out and the network gets locked in an
stable fixed-point state x(n) = x(n + 1) = . . ..
• With W2 , the network starts to oscillate, here with a period 4, that is x(n) =
x(n+4) = x(n+8) = . . .. Other periods (whose length are the higher powers
of 2) could also be obtained by using other scalings of W1 , not shown.
• When the weight matrix is scaled above some critical value (which here is
about 3.55), all regularity breaks down and the behavior becomes chaotic.
There are no periodicities; patterns will never repeat exactly.
Biological neural networks exhibit the same kinds of characteristics. Stable
states occur, for instance, in your motor control circuits of your neck muscles
when you keep your head steadily upright — or in an epileptic seizure. Periodic
activation sequences occur extremely often — for instance in the neural ganglia,
embedded on the surface of your heart, which generate heartbeat signals; or in
central pattern generator neural networks in your spine that generate a basic walk-
ing pattern for your leg muscles; or when you memorize a number by mentally
repeating it. Numerous occasions have been hypothetically proposed for chaotic
activation patterns, for instance as being the carrier of complex sensory represen-
tations or as the ’mind-wandering’ processes in your brain when you shut your
eyes and let your thinking flow freely.
The dizzy dancing of just two coupled neurons already is hard to follow even
with a mathematician’s eye. But that’s only 2 (two) neurons. You own about 85
billion neurons (or do they own you?). The door to understanding their dance
is the study of recurrent neural networks. And the key to open that door is the
theory of dynamical systems (DS). In this section I will treat you to a crash course
in DS essentials.
60
1
0.8
W1
0.6
0.4
0 5 10 15 20 25 30 35 40 45 50
3.3 · W1 −0.5
−1
0 5 10 15 20 25 30 35 40 45 50
4.3 · W1 −0.5
−1
0 5 10 15 20 25 30 35 40 45 50
0.4631 0 0
0.4631
0.4631
0.4631 −1 −1
−2 0 2 −1 −0.5 −1 −0.5 0
Figure 18: Running the 2-neuron network with weight matrices W1 (blue, top plot
and first plot in bottom row), W2 = 3.3 · W1 (green, second plot and
second panel in bottom row), and W3 = 4.3 · W1 (red, third plot and
last panel in bottom row). The last 50 time steps from the 1000 step
simulation are shown. First three plots show the development of the
activations of the two neurons plotted against time. The panels in the
bottom row give another visualization of the same traces: here the two
axes of the panel correspond to the two components x1 , x2 of the state
vectors x(n), giving a 2-dimensional representation of the temporal
evolution.
61
Symbolic Numerical
Texts, event and action Physiological models, psychometrical
sequences, DNA, conceptual measurements, motor control
reasoning
Non-deterministic
Language competence models,
Deterministic grammatical sequence generation
Electrodynamics, artificial
neural networks, mean-field Stochastic
spike trains, speech, language
models
performance models
Autonomous
Sleep dynamics (?), central Non-autonomous / input driven
pattern generator models (?), well, ...almost every real-life system
circadian clocks (?)
low-dimensional high-dimensional
Hodgkin-Huxley or FitzHugh-
network-level modeling, modeling of
Nagumo model of neurons,
cognitive processes
oscillator models
linear non-linear
"classical" analysis of neural neural pattern generators, chaotic
dynamics as signals dynamics, coupled oscillators
homogeneous non-homogeneous
fully or sparsely connected modular or hierarchical neural
neural network circuits and architectures
non-stationary
stationary speech (short timescale), adaptation
neural noise, speech (long processes
timescale) evolutionary
language evolution, ontogenesis,
cell differentiation, learning
Figure 19: Some (only some) sorts and aspects of dynamical systems. Each item
comes with its own collection of dedicated mathematical formalisms.
62
Physicists and mathematicians have been developing formalisms for modeling
dynamical systems (DSs) since several hundred years, with Newton setting the
stage, followed in the 20th century by theoretical biologists, control engineers,
economists, neuroscientists and machine learners. Since everything under the sun
is a dynamical system, it is clear that there cannot be a single unified mathe-
matical theory of dynamical systems — it would be a theory of everything and
as mystifying as the world around us. Formal models of dynamical systems thus
always focus on a particular class of real-world (or computer-simulated) systems
and analyse a particular aspect of them. Figure 19 gives a coarse overview of main
aspects that have been studied.
Since brains are made to understand the world, brains have to try to capture
all aspects of the world; and since brains are made from recurrent neural networks,
recurrent neural networks are — or ultimately will be, if evolution carries on —
as complex as the real world. Thus, all the aspects of dynamical systems listed in
Figure 19, plus many more, are relevant for the study of recurrent neural networks.
Close your eyes and look inside...
It is clear that in this course we can only scratch the surface. I start with
introducing some basic concepts and terminology (summarized in Figure 20):
Time. There are two main formalizations of the arrow of time. Discrete time
jumps forward in “time steps” or “increments”, mathematically captured by
the integers Z. A discrete-time DS model can be specified for the entire
timeline Z with infinite past and future, or it can be started (typically at
time n = 0) and run forever which yields a timeline equal to N, or it can
be modeled for a finite interval [0, 1, . . . , N ]. In mathematical formulas one
63
mostly uses the symbol n to denote discrete time points. Figure 20 is based
on a discrete-time view. The main alternative is continuous time where time
points are real numbers, typically denoted by the symbol t ∈ R.
One often says that a DS evolves, and one speaks of the evolution of a DS.
This is a figure of speech that has nothing to do with biological evolution.
It simply means that the system changes its states as time goes on. Even
if one monitors a simple pendulum, physicists or mathematicians will speak
of the temporal evolution of this system, meaning nothing more fancy than
that the pendulum simply swings.
Update operator. This is the mathematical rule which specifies how a system
state changes as time moves on. There are many namings for update opera-
tors, for instance transition function, system law, or (in logic-based symbolic
AI) inference rule. Basic formats for discrete time update operators are
64
which describes a deterministic temporal evolution of states by an iterated
map T , which can be any function T : X → X , or
If you are not familiar with the ẋ notation, be patient, it will be explained
later in this section. I just mention in passing that when states are functions,
one needs partial differential equations (PDEs). They would become neces-
sary if, for instance, one wishes to model the propagation of neural activation
waves on the cortical surface. ODEs and PDEs can be generalized to cap-
ture stochastic dynamics, yielding stochastic [partial] differential equations
— just to drop the words. In theoretical neuroscience, all of these math-
ematical formalisms are used (and theoretical neuroscientists often have a
theoretical physics background).
Inputs. Many DSs of interest are receiving external input signals. This is partic-
ularly relevant for recurrent neural networks: after all, their main purpose
is to process incoming information. Brains are bombarded with sensory
information. Thus in neural network research one almost always consider
scenarios where the DS of interest — an RNN — is continually “driven” by
an input signal. There is no standard notation for input signals. I will follow
the common usage in the field of signals and systems and usually write the
symbol u for input signals. Inputs u(n) can be anything, for instance just a
single number per timestep (example: a voltage fed to an electric motor), or
a vector (example: a million-dimensional vector made from the pixels in the
image frames of a video input stream), or a symbol (like the letters of a text
that is letter-wise read into an online text translation RNN). The evolution
equations (34), (35), (36) become extended by input terms and now appear
as
65
In mathematical parlance, DSs that have no input are called autonomous
dynamical systems. One has to be a bit watchful about who uses this term,
since in the field of robotics and intelligent agents, when one speaks of “au-
tonomous systems” or “autonomous agents”, the intention is to highlight
that these systems develop their action planning all by themselves. But of
course these robots or “agents” receive sensory input all the time, so they are
not autonomous systems in the terminology of dynamical systems theory.
Another point worth mentioning is that physicists are traditionally inter-
ested in exploring their systems (a pendulum or the solar system or the
helium atom) in laboratory conditions where they go to great lengths to iso-
late their experimental system from external influences. They want to under-
stand the “pure” system just as it is for itself and use the rather deprecatory
word perturbation for the undesired external input. Since the mathematical
theory-building of DS has been motivated and often pioneered by physicists,
still most textbook treatments of DSs concentrate on autonomous systems.
This is a pity for modeling neural networks as DSs, because the bulk of
available theory is limited to autonomous systems only, and it is not at all
trivial to extend that huge, traditional body of theory toward input-driven
systems, RNNs in particular. Mathematical research for non-autonomous
DSs is in its infancy.
In contrast, in the field of control engineering one always has considered
input-driven systems. This was by necessity: if one wants to control a system
(the latter being called “plant” in control engineering), one has to give control
input to the system — if you want to steer your car (the plant) you have to
turn the steering wheel. A very rich body of theoretical knowledge about
input-driven dynamical systems has thus grown over the last hundred years
or so. But this theory was almost exclusively built around linear models of
plants, that is system models where the functions T in (37), (39) are linear.
Thus, even today we still have an unsatisfactory theory coverage of input-
driven DS. Physics and mathematics give us beautiful formal theories of non-
linear autonomous systems, and control engineering of linear non-autonomous
ones. But brains and RNNs are nonlinear and non-autonomous.
Outputs and observations. The states of a real-world DS are often not directly
or fully accessible. One cannot experimentally measure all the electric poten-
tials of all neurons in a brain (which are part of the full brain state), nor can
one monitor every molecule of air in meteorological modeling of the earth’s
atmosphere. All that one can do is measure or observe an embarrassingly
small number of measurables (synonym: observables). For instance, one can
plant a few hundred measurement electrodes into a brain (state of the art
today – but the human brain has billions of neurons); or one can monitor
the atmosphere with a few thousands of weather observation stations. Some
technical or computational systems have built-in measurement mechanisms,
66
in which case one often speaks of system output. In particular, artificial neu-
ral networks have designated output neurons. In mathematical formalism,
measurements (which are “pulled out” from the system one might say) and
outputs (which are yielded by the system “voluntarily” one might say) are
treated in the same way. I again follow the standard notation in the signals
and systems tradition and will usually reserve the symbol y to denote out-
puts. In mathematical formalism, in order to capture outputs, a DS model
is extended by an output function O which at any time point extracts the
measurement from the system state, yielding pairs of equations
x(n + 1) = T (x(n), u(n))
y(n) = O(x(n)), (40)
P (Xn+1 = sj | Xn = si , Un = a) = Ta (i, j)
P (Yn = yk | Xn = si ) = O(i, k), (41)
67
3.2 The zoo of standard finite-state discrete-time dynam-
ical systems
When scientists (of all sorts, in all fields) are faced with a very complex system
that they want to understand (brains, societies, molecules, RNNs), a good strategy
is to start analyses with simple and small formal models of the respective target
system. It is clear that highly simplified models cannot capture the full richness
of the system under scrutiny, but one can hope that simplified, highly abstracted
models do capture some of the essential, fundamental properties of the target
system. In this section I present the most commonly used formalisms for DS
modeling which are simple in the sense that
• discrete time is used,
• the state space is assumed to be finite: X = {x1 , . . . , xL },
• possible input and output values likewise come from finite sets which I will
denote by U and Y respectively.
Such models are often called finite-state models. Note that “finite” does not mean
small. A digital computer can be modeled (and is modeled in theoretical computer
science) as a finite-state system. A state x(n) of a digital computer is a binary
vector which assigns the on- or off-switching states to all the Boolean gates and
memory cells that make up the computer’s circuitry. I have no idea how many
gates a modern PC has, but certainly many millions. Say, one million. Then there
are 21,000,000 different possible states: |X | = 21,000,000 , a number that by far exceeds
the number of atoms in the universe. But it’s finite.
I will now take you on a fast guided tour through the zoo of finite-state models.
This is a long tour that will go over many pages, but I would think that you have
already seen much of this material before and you’ll be meeting some old friends.
68
A
q States X = {p, q, r}
p Inputs U = {A, B}
B
A A
B
r
Lookup table for the
B Transition graph transition function
T: X × U → X
example trajectory: A B
p q r
A B B A B
q q p
... ...
r q r
q q p r q p
69
A
q States X = {p, q, r}
p
Inputs U = {A, B}
A,B
B A,B
B
A A
r Lookup table for the
B transition function
T: X × U → Pot X
example trajectories:
A B
A B B A p {q, r} ∅
... ...
q {p, q} {p, q, r}
q p dead end
r {q} {p, r}
p r r p q
the same starting state, many trajectories are possible. It may also happen that
T (x(n), u(n)) = ∅. Then the trajectory comes to a halt in a dead end, since there
is no possible next state (as happens in Figure 22 with T (p, B) = ∅).
NFAs yield simplified models of RNNs as in the following example. The original
state space of an RNN with L neurons that uses the system equation
70
is believed by some (including me) to be the deepest unsolved problem of mathe-
matics, period (Cook, 2000).
I emphasize that non-determinism is not the same as stochasticity. NFAs
and other non-deterministic models of DSs only make statements about which
trajectories are possible at all, not how probable they are.
States X = {p, q, r}
.2 .5
q Initial state distribution p0 Prob(X)
p
.9 .1
0.0 Transition kernel T: X → Prob(X)
.3 .3
.3
r
.4 Markov transition matrix for T
p q r
p 0.2 0.5 0.3 rows sum to 1:
p0 =
q 0.9 0.1 0.0 „Markov matrix“ or
0.3 0.3 „stochastic matrix“
r 0.4
71
2. Later times: if the trajectory point x(n) = xi ∈ X at time n has already
been determined, the next state x(n + 1) = xj is picked from X according
to the probabilities P (Xn+1 = xj | Xn = xi ) which can be found in the i-th
row of the transition matrix.
This generation process can be continued for an arbitrary duration. Among
all initial trajectories of length N , a specific trajectory xi0 , . . . , xiN −1 has the prob-
ability of being generated given by
P (X0 = xi0 , . . . , XN −1 = xiN −1 ) =
= p0 (i0 ) P (X1 = xi1 | X0 = xi0 ) · · · P (XN −1 = xiN −1 | XN −2 = xiN −2 )
= p0 (i0 ) T (i0 , i1 ) · · · T (iN −2 , iN −1 ). (45)
A Markov chain is a particularly simple kind of a stochastic process. Generally,
a (discrete-time) stochastic process is a sequence (Xn )n∈N of random variables
which all take their values in the same state space X . Markov chains are special
in that the probability of what is the next state x(n + 1) only depends on the
previous state x(n). In general stochastic processes, the probability of x(n + 1)
to appear in a trajectory may depend on the entire initial part of the trajectory,
requiring some formalism (of which there are many) that can capture all the
conditional probabilities
P (Xn+1 = x(n + 1) | X0 = x(0), X1 = x(1), . . . , Xn = x(n)).
Stochastic processes where the probabilities for the next state only depend on
the previous state are generally said to have the Markov property; they are called
Markov processes. Markov processes can have infinite or continuous state spaces,
and they can also be defined for continuous time. They play a big role in physics,
because by that dogma of physics that I mentioned earlier, all that can be inferred
or known about the future of a physical system is contained in its present state
— which is just another wording for the Markov property. According to physics,
in particular quantum mechanics, the world is a Markov process, though a more
complex one than what we see in Figure 23.
In the field of stochastic processes, one often calls trajectories realizations or
paths of the process.
72
• in Mealy machines: by emitting the outputs not from the states, but from
the transitions of the underlying DFA, specified by an output function of
type O : X × U → Y.
Figure 24 shows a simple example.
States X = {p, q, r} Inputs U = {A, B} Outputs Y = {0, 1}
Moore: O: X → Y Mealy: O: X × U → Y
A 0 0
A
q q
p p 0
B B
A A B 1 A A
1 B 1
r 0 r
0 B
B
1
a Moore trajectory a Mealy trajectory
A B B A B B B A B B
q q p r q p q p r q p
0 1 0 0 1 1 0 1 1 0
Figure 24: A 3-state Moore and a 3-state Mealy machine.
73
3.2.5 Probabilistic finite-state output-observed DS — aka hidden Markov
models
In the same way as DFAs can be augmented by adding an output mechanism, one
can add such an item to Markov chains. Since here the general perspective is a
stochastic one, the output mechanism also takes on a probabilistic format. Thus
a hidden Markov model (HMM) is defined by a “hidden” (unobservable, invisible)
Markov chain defined as in Section 3.2.3, plus, for every state x ∈ X , conditional
probabilities
P (Yn = y(n) | Xn = x(n))
which declare how probable it is to observe (measure) a “visible” outcome y(n) ∈
Y = {y1 , . . . , yM } at time n when the ongoing random stochastic process is in state
x(n) (see Figure 25). These probabilities are called emission probabilities. They
are usually all collected in an emission matrix E of size |X | × |Y| whose i-th row
collects the emission probabilities of all possible observations y ∈ Y if the process
is in state xi .
example trajectory:
q p r q p unobservable, „hidden“
0 1 0 0 1 observable, „visible“
Here and in Figure 25 I used a finite set of possible observation values. HMMs
can be generalized to continuous observation spaces Y ⊆ Rn by introducing con-
tinuous emission distributions on Y, conditional on the hidden states.
In machine learning and in the natural sciences at large, HMMs are a model
of choice when one wishes to capture a situation where the states of a stochastic
real-world DS cannot be directly or completely measured, and all that one has
available are a few measurement devices that can observe the ongoing process.
74
For example in a chemical factory, what happens inside a reactor vessels will be
monitored only by a temperature and a pressure sensor, which does not at all
fully capture all the chemical reactions, turbulent flows, chemical gradients etc.
inside the reactor. Or, closer to the theme of our lecture, what happens inside a
brain can be observed by the few dozen channels of an EEG scalp electrode array.
The measurement process itself is typically plagued by stochastic uncertainties
(“measurement noise”), hence the need to introduce emission probabilities.
In machine learning, HMMs have been playing a central role in the develope-
ment of speech recognition technologies. The underlying hidden Markov process
models what happens in a brain of a speaker, and the observable signal is the
speech signal coming out of the mouth of the speaker owning the brain. Since
about 2010 HMMs have been superseded by recurrent neural networks in speech
recognition applications. The role of HMMs in the natural sciences however is
irreplaceable by neural networks because HMM models allow a scientist to rigor-
ously analyze the obtained models, for instance making probabilistic predictions
about the long-term future evoluation of a process. HMMs furthermore are still
un-replaceable in computer linguistics and computational linguistics, because they
allow the modeler to connect assumed cognitive mechanisms (the hidden Markov
chain part) with observable speech or text output in a transparent way. RNN
models do not allow this because they are blackbox models which are so complex
and so devoid of mathematical theory that an analysis of system properties is
virtually impossible.
From a machine learning perspective, the learning task is the following: Given
a sequence y(0), y(1), . . . , y(N ) of observations made from an invisible “hidden”
Markov process, learn (estimate, identify) a HMM which can explain the observed
sequence of measurement values. The tutorial text of Rabiner (1990) gives a beau-
tiful introduction to the theory and practical algorithms of HMMs (31K Google
cites, as of April 2022). This classic piece of tutorial scientific literature has been
a major booster for the widespread use of HMMs.
75
States X = {p, q, r}
sound hard, it is hard and needs large volumes of training data and considerable
computing resources.
In the field of reinforcement learning, a subfield of machine learning, controlled
MCs are further expanded by assigning a reward to every transition xi → xj under
an input uk . The underlying interpretation here is that inputs u(n) are actions of
an intelligent agent carried out in a randomly behaving environment (modeled by
the underlying Markov chain, whose states x(n) correspond to external situations
the agent can be in). The objective in reinforcement learning is to let the agent
learn a policy, that is an action-decision-making rule π : X → U which tells the
agent which action to choose next when it finds itself in a situation x(n). A policy
is good to the extent that when it is applied at every time step, the accumulated
reward is maximized. In this context (controlled Markov chain plus reward plus
policy) one speaks of Markov decision processes. Markov decision processes, in
which the policy functions π : X → U are learnt by deep neural networks, play
a leading role in today’s deep reinforcement learning research. You will likely
have heard of machine learning systems beat human champions in chess, Go, and
even advanced computer games — those successes are based on Markov decision
processes combined with deep neural networks.
76
servables, behaves in time when it is steered, controlled, influenced, modulated or
perturbed by external inputs.
0 1 0 0 1 observable
Figure 27: A 3-state Markov chain controlled by two inputs and observed by two
outputs
Stochastic models of this kind play a leading role in designing and training
of artificial intelligent agents like autonomous robots, software agents, or game
characters. These models are the agent’s world model: the states x ∈ X are the
states (situations, conditions) of the external environment, for instance locations
(for a mobile robot), or tool loads and tool tip poses (for a manipulator robot
arm), or the mental states of opponents (for competitive game creatures). The
inputs u ∈ U are the possible actions that the agent can execute. The observations
y ∈ Y are what the agent can see, hear or otherwise sense in its environment. In
the fields of reinforcement learning, rewards and policies are added to the picture,
and the whole thing is called a Partially Observable Markov Decision Processes, or
POMDP for short (pronounce “Pomm-Dee-Pee”). My favorite classical (accessible
but not reflecting the deep learning state of the art) tutorial texts are Basye et al.
(1995) for the case without rewards and policies and Kaelbling et al. (1996) for
the full reinforcement learning picture. It is always a good idea to read seemingly
outdated classical texts because when the field was fresh and young, authors still
explained things that in later texts would be taken for granted. The word /
acronym POMDP is also used for the input-driven, output-observed, finite-state
DS without the reward and policy ingredients; the usage of terminology is not
entirely consistent. Like with Markov decision processes, the discrete input and
77
observation spaces that I displayed here can be extended to continuous-valued
ranges, which today is typically effected by installing (deep) neural networks in
the overall machinery. This leads out of the scope of our lecture.
How it works
Defining components:
1. Initial configuration: x(0): Z → Q • A finite set of local cell states, e.g.
Q = {red, green} (visualized as colors)
x(0) ... ... • a local transition function
Tlocal: Q×Q×Q → Q
2. Apply local transition function to all
local neighborhoods (cell triplets),
get next configuration x(1) Example transition function:
x(0) ... ...
x(1) ... ...
3. Iterate, obtain trajectory of
configurations x(0), x(1), x(2),...
Figure 28: A 1-dimensional CA with two local states (red, green) and a neighbor-
hood width of three.
78
1-dimensional example shown in Figure 28 I use a small neighborhood which com-
prises cell i itself, plus its neighbors i − 1 and i + 1. In formal terms, if the chosen
neighborhood comprises l cells, Tlocal is a function from Ql to Q. Applying Tlocal
in parallel to all cells and their neighborhoods in x(0) gives the next configuration
x(1). To get x(2), repeat this procedure, etc. The same local transition function
is used for all cells. A CA can thus be specified quite economically by a single
simple lookup table for Tlocal .
The evolution of 1-dimensional CAs can be nicely represented in a graphical
format where each row of an evolution graph corresponds to one configuration,
and time progresses from the top to the bottom. Figure 29 shows four (famous)
examples which exhibit four different qualitative kinds of CA rules, today called
the Wolfram classes of CAs after their discoverer, Stephen Wolfram. In all the
four plots, the first configuration is a random binary pattern. In class 1 CAs,
all “activity” quickly dies out. In class 2, a stationary repetitive (cyclic in time)
pattern evolves. Class three CAs generate evolutions that look “just like noise”,
and in fact such CAs were used by Wolfram as random number generators in
early versions of the Mathematica software package that he developed. Class
4 includes CAs whose evolution leads to the “self-organization” of increasingly
globally structured patterns. Universal Turing machines can be embedded in such
class 4 CAs, leading to the arguably simplest know models of universal computing
systems (Cook, 2004), see also the nice Wikipedia page on this subject.
Figure 29: The four classes of 1-dimensional CAs identified by Steven Wolfram.
Image taken from Wolfram (2002).
In two dimensions, one cannot plot entire evolutions but must be satisfied with
snapshots of individual configurations. Figure 30 shows a few.
The evolution of spatial patterns in continuous substrates is usually modeled
with partial differential equations (PDEs). By discretizing the state space, a
PDE can be approximated by a CA. Cellular automata thus can be used as com-
putationally cheap and mathematically analyzable tools to study spatial pattern
formation in general. One of my personal favorites is the simulation of the geomet-
rical patterns that appear on the surface of exotic seashells as the shell is growing
(Figure 31). Note that a seashell is growing through expanding its “lip front” and
winding it around its shell body while it is growing. The geometrical patterns are
79
Cellular automata (CA) – visual demo
Figure 30: Some configurations occurring in some 2-dimensional CAs. Pictures re-
trieved from [Link]/modernca/[Link]
Retrieved from [Link] --- Not online anymore. A beautiful some
Browser-based CA experimentation lab is at [Link]
years ago — not online anymore. But there are other online CA sim-
ulation engines that invite playing with these amazing systems.
80
CAs, comments
81
CNN. In each global update cycle, their locally computed new pixel values are
used to give the new global ’picture’. The fancy thing about this architecture is
that the CNN can be trained from data. For instance, if the training data would
be a video of a burning fire, after training, the neural cellular automaton would
operate as a generator of new burning fire videos. Take a look at the splendid
demos in [Link]
You see, everything in the world of science seems somehow connected to ev-
erything else, and this is one of the reasons that this subsection grew so long.
Another reason is that I have a personal relationship with CAs: they were the
topic of my diploma thesis.
82
saw a used hardcover exemplar offered on Amazon for more than 1000 Euro). I
purchased copies several times and gave them away as presents on very special
occasions. You can “read” (not the right word for bathing in graphical art) it on
one weekend like any graphical novel.
These equations look a bit frightening maybe — we will soon see how they can be
dramatically simplified, but for the moment let us stick to this basic formulation.
In these equations, ẋ and ẏ are the customary shorthands for dxdt
and dy
dt
, the
derivatives of the state variables x, y with respect to time — in other words, the
temporal change rates of these state variables. One could also correctly write these
equations as
bringing to the surface that these equations hold at any time t in the evolution of
the system. However, it is common usage to omit the t’s.
The right-hand sides of (46), (47) each are both functions (I called them f and
g) of x as well as of y: the change rate ẋ(t) at any time t depends on both the
current value x(t) and the current value y(t). One says that the state variables
x, y are coupled.
Yet another way to write these equations is to present them in vector form:
ż = h(z),
where
′ f (x, y)
h(z) = h((x, y) ) = .
g(x, y)
Now let us go graphical. Note that h : R2 → R2 is a function which assigns
to every vector z ∈ R2 another vector h(z) ∈ R2 . Such functions are called vector
fields.
83
1
0.8
0.15
0.6
0.4
0.2
&̇ (0)
0
(x(0), &(0))′
−0.2 0.1
−0.4 !̇ (0)
−0.6
−0.8
−1
0.05
−1 −0.5 0 0.5 1 0.46 0.48 0.5 0.52 0.54 0.56
Figure 32: The vector field yielded by Equations 46, 47. Left: A large-area view.
The blue arrows show the vectors resulting from h(z) for some points
z on a grid. Three trajectories, each started at time t = 0 in the points
z(0) marked by green circles, are evolved for a duration of 2 time units,
that is for t ∈ [0, 2]. Right: a zoom around the starting point of one of
the trajectories.
84
it by its Euclidean coordinates, z = (x, y)′ , which is what I did above. Another
representation is to specify z by its distance r from the origin (its “radius”, hence
the symbol r); and the rotation angle φ of the line leading from the origin to z
(Figure 33).
z = (x, y) y
1
= (r, j)
0.9
0.8
0.7
#̇
0.6
!̇
0.5
r y
0.4
0.3
0.2
j
0.1
0
0 0.2 0.4 0.6 0.8 1
x x
Figure 33: Left: Representing a point z ∈ R2 in Euclidean (green) and polar co-
ordinates (red). Right: According to the decoupled ODEs (48), (49) in
polar coordinates, the motion of a state z is the sum of two indepen-
dent motion components in orthogonal directions, one along the radius
beam and the other along the circle on which z lies.
Now, if in our Whirl system a state z(t) changes its x and y coordinates with
the change rates given by (46), (47), also its polar coordinates r, φ are changing.
Omitting the derivation, the equivalent laws of change for r, φ are given by the
ODEs
ṙ = −r3 + r (48)
φ̇ = 1. (49)
These formulas look certainly simpler than their equivalents in (x, y) coordinates.
The main advantage of the new version over the old is not a simpler look, but the
fact that the new coordinate variables (r, φ) are decoupled. The rate of change ṙ
of the radius only depends on the current radius, and the rate of angular change
doesn’t even depend on anything, since it is constant. In this decoupled system
we can analyze the dynamics in the direction of the radius independently from
the angular dynamics. A trajectory point z(t) moves in the r direction by coming
closer to the origin or moving away from it. And simultaneously, z(t) rotates
anti-clockwise around the origin, with constant angular velocity φ̇ = 1. It’s as if
z(t) sits on a circling radar beam like on a rail on which z(t) can slide inbound or
outbound (Figure 33 right).
85
3.3.2 A cyclic attractor and a repellor
Now that we have made friends with the Whirl system, we will investigate it more
closely (using the polar coordinate representation) and learn a lot.
The vector field shown in Figure 32 has some interesting geometrical properties
which determine the long-term evolution of trajectories.
1 1 1
0 0 0
−1 −1 −1
A B C
Figure 34: A: Two long-term trajectories (black) and two trajectories evolved for
a shorter time (red). Starting points are marked by small green cir-
cles. B: Four short-term trajectories started close to the origin. C:
An elementary phase portrait rendering of the Whirl system, showing
attractors in green and repellors in red, plus two exemplary trajecto-
ries.
Figure 34A shows two trajectories that are started from different initial points
and which were computed for a long timespan (black), and two trajectories that
were computed for a shorter timespan (red). Observation: all these trajectories
converge to the unit circle, whether started from its inside or outside. A state
z(t) that exactly sits on the unit circle will forever rotate along that circle coun-
terclockwise with an angular velocity of φ̇ = 1. One says that the unit circle is an
attractor of this DS, or more specifically, a cyclic attractor. Cyclic attractors indi-
cate the presence of oscillatory (or periodic) behavior of a DS. Our Whirl system
is a case of an oscillator system.
Generally and somewhat loosely speaking, in an n-dimensional dynamical sys-
tem ż = h(z), where z ∈ Rn , an attractor is any subset A ⊆ Rn which has the
properties (i) that trajectories started “near” A will converge to A if evolved for
infinite time t → ∞, and (ii) that any trajectory started in A will remain forever
in A, and (iii) that A is minimal, that is, no subset of A has properties (i) and (ii).
From a purely mathematical perspective, attractors are sets of states, and one
may speak of attracting sets if one wishes to emphasize this fact. Attracting sets
can have many other geometrical shapes besides circles — we will see some other
types later. A general, precise definition of attractors in dynamical systems needs
86
concepts from topology and is beyond our scope, but this intuition of “letting all
nearby trajectories converge” will be good enough.
Now consider Figure 34B. It shows a few trajectories that are started close
to the origin. They all move away from it (and in the long run will converge to
our cyclic attractor). The origin acts as a repellor. Like attractors, repellors R
are sets of states — here, the repellor set contains a single point only, namely the
origin R = {(0, 0)′ }. It is a point repellor. Repellors are the opposite of attractors.
They are (again a little loosely speaking) minimal state sets characterized by the
property that any trajectory that is started near the repellor will move away from
it.
Figure 34C graphically collects the essentials of the Whirl system, by showing
the repellor (red) and the attractor (green) and a few instructive trajectories,
without rendering the vector field. Such graphical “thumbnails” of dynamical
systems are called phase portraits. There is no precise definition of what a phase
portrait is — you may call any graphics that highlights the dynamics of a DS a
“phase portrait”. Attractors and repellors and a few trajectory lines are typically
drawn in phase portraits, but further elements can also be added.
0.5
!̇
0
−0.5 motion of r
−1
−1.5
! ̇ = − r3 + r
−2
0 0.5 1 1.5
r
Figure 35: The radial component of the Whirl dynamics.
Now we take a closer look at the Whirl equations (48) and (49) and investigate
why, exactly, we get the dynamical behavior that we have qualitatively described
above. The rotational component φ̇ = 1 needs no deep thinking — it tells us
that every trajectory is rotating around the origin with constant angular velocity.
Everything that is interesting happens in the radial component ṙ = −r3 +r. Figure
35 plots ṙ against r. We see that ṙ is equal to zero when r = 0 or r = 1. That
is, when a state z has a radial component of 0 or 1 — that is, when it lies in the
origin or on the unit circle — then ṙ is zero, which means that z will not move
in the r direction. Or said in educated terms, r = 0 and r = 1 are fixed points of
the radial dynamics ṙ. And whenever 0 < r < 1, the motion ṙ is positive: when
z has a radial component between 0 and 1, its radial component will grow with
time — the trajectory moves closer to the unit circle. Conversely, when r > 1,
87
then ṙ < 0 which means that states outside the unit circle will also move closer to
it by reducing their radial component.
Furthermore, the fixed point r = 0 is repelling: the slightest perturbation by
pushing the trajectory from 0 to any 0 + ε will initiate a state motion away from
the origin. And the fixed point r = 1 is attracting: all values of r around it move
toward r = 1. For fixed points, instead of calling them repellors or attractors,
one also says that they are instable or stable fixed points. — Once this radial
ṙ is understood, the Whirl phase portrait (Figure 34) becomes “obvious” (start
running when a mathematician says this word! ,) — it’s just the ṙ dynamics on
a rotating radius beam.
0.5
1 1
!̇
0.5 0.5
0 0 0
−0.5 −0.5
−1
−1
−0.5 −1 −0.5 0 0.5 1
0 0.3 0.6 0.9
r −1 −0.5 0 0.5 1
A B C
Figure 36: Another ṙ and the resulting phase portrait. A: The radial rate of
change function. B: Resulting phase portrait. C: The two basins of
attraction for the central point attractor and the cyclic attractor with
radius 0.6
Once one has the knack of it, it is easy to create more complex DSs in the same
spirit. Figure 36 shows a ṙ function that passes through the zero line four times,
at r = 0, 0.3, 0.6, 0.9 and the resulting phase portrait. The constant rotational
motion φ̇ = 1 was kept. This DS has a point attractor in the origin and a cyclic
attractor with diameter 0.6, and two cyclic repellors with diameters 0.3 and 0.9.
No further explanation should be necessary.
Two more broadly used concepts: Given an attractor A, its basin of attraction
(or attractor basin) is the set basin(A) ⊆ Rn of all states z which will eventually
converge to A. Figure 36C shows the basins of the two attractors in this system.
The basin for the point attractor in the origin is the set {z ∈ R2 | kzk < 0.3} and
the basin for the cyclic attractor with diameter 0.6 is {z ∈ R2 | 0.3 > kzk < 0.9}.
— Segments of trajectories that are not close to an attractor are called transients.
This term is used in a loose manner because it is not mathematically clear what
“close” means. It’s more an intuitive term used when one wants to emphasize
that a state motion is not (yet) well described by the controlling influence of an
attractor by which the trajectory might eventually be captured.
88
Attractors are, I think, the most important and decisive phenomenon in math-
ematical models of real-world dynamical systems. The reason: if you look at any
real-world DS — for instance a stone (yes! it’s a dynamical system! its atoms
oscillate like mad!), or your wristwatch, or your brain — then this system most
likely has already existed some time before you look at it. But a system whose
dynamics has evolved for some time will be in, or close to, an attractor. That is
the very definition of attractors. Thus what we observe when we look around are
all attractors, and never repellors. While from a pure-math perspective, repellors
and attractors are symmetric concepts (a repellor in ż = h(z) becomes an attrac-
tor in ż = −h(z) and vice versa, think about it!), their real-world manifestations
are fundamentally different and not interchangeable.
The Whirl system has simple, low-dimensional system equations, and only
a single attractor. This may have given you the impression that attractors are
somehow special and few. But DSs that are not as simple and small as Whirl may
host very large numbers of attractors. This holds in particular for recurrent neural
networks. In Section 5 we will meet with a kind of RNN called Hopfield networks in
which we can embed very many attractors by training the network. Each attractor
corresponds to an input image that the network has learnt to memorize and that
it can remember and retrieve when it is shown a new test image that is similar
enough to the stored image.
1
ṙ = −r3 + r + sin(3φ) sin(r) (50)
2
1
φ̇ = 1 + sin(4r). (51)
2
We have added another vector field to our original one — we have perturbed
our good old Whirl, even quite strongly. The two state variables r, φ are now
coupled, and it seems impossible to predict how the changed system will behave.
We plot a phase portrait of the new system and compare it to the portrait of our
original Whirl system (Figure 37).
You should be surprised! The perturbed system doesn’t look very different
from the original one! It also has a cyclic attractor — which looks a bit wobbled
but it’s still a cyclic attractor with an unstable fixed point in the middle. The main
structural elements have been preserved: the two phase portraits are qualitatively
equivalent. Intuitively, if you imagine that the original phase portrait had been
printed on a rubber sheet, you could smoothly stretch that rubber sheet and
thereby deform the original portrait until it exactly (up to some constant scaling
factor of the vector arrow lengths) matches the perturbed one. More precisely,
89
1.5 1.5
1 1
0.5 0.5
0 0
−0.5 −0.5
−1 −1
−1.5 −1.5
−1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5
Figure 37: The original Whirl system, drawn once again (left) and the perturbed
system (right).
ṙ = 0, (52)
φ̇ = 1. (53)
You can easily see from the system equations (52) & (53) how the phase portrait
will look like: all trajectories are perfect cycles (Figure 38 left). But this needs
extraordinary precision in the vector field: after one rotation, a trajectory must
90
1.5 1.5
1 1
0.5 0.5
0 0
−0.5 −0.5
−1 −1
−1.5 −1.5
−1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5
Figure 38: Left: the pure circling dynamics of the system (52) & (53). Right:
Almost the same system, but slightly perturbed.
hit its own tail exactly. A slight perturbation of the vector field will make the
returning trajectory miss itself, resulting in inward or outward spiraling. Figure
38 shows the latter case.
The two phase portraits in Figure 38 are not rubber-sheet equivalent: one can-
not pull and shear a rubber canvas without rupturing it such that circles become
spirals. Notice that besides the circling vs. spiraling behavior, there are other
qualitative differences between the two systems. The pure circling system has no
attractor and no repellor, while the perturbed system shown in the figure has a
repellor at the origin.
If a DS D has the property that the slightest change of its vector field will
change qualitative properties of its phase portrait, one says that D is structurally
instable.
The mathematical study of structural (in)stability of DSs is involved and since
about half a century an important research subject, with a full picture far from
being known. Some important classes of dynamical systems have been shown to
be generically structurally stable. This means, loosely speaking, that if one writes
down a randomly invented system equation for a DS from that class, it will be
structurally stable with probability 1. If you want to dig deeper, check out the
Wikipedia article on structural stability, or even deeper-digging, the Scholarpedia
article on Morse-Smale systems, an important large class of generically structurally
stable DS which includes many real-world physical systems.
I find structural stability not only mathematically intriguing, but also very
reassuring indeed. We would not exist, and neither would the universe exist, if
structural stability were not a very typical property of DSs. If most real-world DS
were structurally instable, the slightest perturbations would change their behavior
dramatically. Imagine what would happen if the many dynamical systems in our
universe — from the hydrogen atom to planetary systems, with plants and animals
in between — would not be structurally stable — pouuffff! bouummmm!
91
Structural instability is splendidly described in a classical work of German
literature, the Simplicius Simplicissimus by Hans Jakob Christoffel von Grim-
melshausen, written 1668. This novel recounts the life and adventures of the
vagabond Simplicius in the Thirty Year’s War. In one episode (Book 6, Chapter
9) the hero meets a ... a what? That is the question, because that strange being
met by Simplicus is structurally instable:
I was once walking around in the forest and listening to my idle thoughts when
I found a life-sized stone statue lying on the ground. [...] it began to move by
itself and said: “Leave me alone. I’m Soonchanged.” [...] Then he took the book
which I happened to have with me and, after he had changed himself into a scribe,
he wrote the following words in it: “I am the beginning and the end, and I am
valid everywhere. [...]” After he’d written this, he became a large oak tree, then a
sow and then quickly a sausage, and then some peasant’s dung. The he changed
himself into a beautiful meadow of clover and, before I could turn around, into a
cow-pie; then he became a beautiful flower or sprout, a mulberry tree and then a
beautiful silk rug and so on till he finally changed back into human form [...] Then
he changed himself into a bird and flew quickly away. (Grimmelshausen (2012),
translation by Monte Adair)
3.3.4 A bifurcation
We are now ready to call on stage another key actor in the DS show: bifurcations.
We again consider the Whirl equations, but add a scaling parameter a to the radial
component, leaving the constant rotation unchanged:
• When a < 0, phase portraits show a point attractor in the origin, with
trajectories whirling toward it. I mention without proof that these systems
are all structurally stable.
92
1.5 0.5
1
0
0.5
−0.5
a = −1 0
−1
−0.5
−1
−1.5 a < 0:
structurally
stable with
−1.5 −2
−1.5 −1 −0.5 0 0.5 1 1.5 0 0.5 1 1.5
a = − 0.2 0
−1
−0.5
−1.5
−1
−1.5 −2
−1.5 −1 −0.5 0 0.5 1 1.5 0 0.5 1 1.5
1.5 0.5
1
0
0.5
a = 0:
−0.5
structurally
a=0
instable with
0
−1
−0.5
one point
−1
−1.5
attractor of a
−1.5
−1.5 −1 −0.5 0 0.5 1 1.5
−2
0 0.5 1 1.5
special kind
1.5 0.5
1
0
0.5
−0.5
a = 0.2 0
−1
−0.5
a > 0:
−1
−1.5
structurally
−1.5
−1.5 −1 −0.5 0 0.5 1 1.5
−2
0 0.448 1 1.5
stable with
one point
1.5 0.5
repellor and
1
0 one cyclic
0.5 attractor
−0.5
a=1 0
−1
−0.5
−1.5
−1
−1.5 −2
−1.5 −1 −0.5 0 0.5 1 1.5 0 0.5 1 1.5
Figure 39: Whirl phase portraits and ṙ functions for several settings of the control
parameter a.
93
• When a > 0, phase portraits show a point repellor in the origin and a
cyclic attractor that lies on a perfect circle with a diameter equal to the zero
crossing point of the ṙ curve. Again all of these systems are structurally
stable.
• When a = 0, and only for exactly this unique value of the control parameter,
we find again a single point attractor in the origin. However, there is a
qualitative difference between the case a = 0 and the cases a < 0 which is
not immediately visually apparent. Namely, in the cases a < 0, when one
would measure the length of any trajectory, from some arbitrary starting
point until the trajectory line meets the point attractor, then this length is
finite, and when the trajectory spirals toward the origin from some starting
point, it revolves around the center only finitely many times. In contrast,
the length of any trajectory in the a = 0 case, measured from some starting
point toward the center, is infinite; and trajectories revolve around the origin
inifinitely often as they come closer. The trajectories shown in Figure 39 for
a = 0 have been computed for a very long but not infinite time, thus they
don’t reach the center and the plot leaves a little unvisited area around it.
The mathematical cause for this behavior is that in the a = 0 case, the slope
of ṙ at r = 0 is zero. — Finally I note that the a = 0 system is structurally
instable. The slightest perturbation of this system will change its behavior
either to a < 0 or the a > 0 type, i.e. one would get either a finite-length
spiraling-in or a repellor + cyclic attractor type of system.
• The a = 0 system is not rubber-sheet equivalent to the a < 0 systems,
because, intuitively speaking, one would have centrally “unwind” the a = 0
rubber canvas infinitely often to get a 1-1 mapping to a a < 0 phase portrait,
and that cannot be done by a continuous map (it would have a singularity
in the origin).
What we see here is an example of a bifurcation. Bifurcations occur in many
ways, and there are many types of bifurcations. The general scenario for bifurca-
tions looks as follows:
• Consider a system ż = h(z, a) with control parameters.
• Typically, when the control parameters a are changed by small amounts, the
resulting systems will be structurally similar, and each is structurally stable.
• However, for some isolated, critical values a∗ , the system ż = h(z, a∗ ) will
be structurally instable.
• If one considers a sequence of control parameters a0 , a1 , . . . , a∗ , . . . , aN which
passes through a∗ , then systems with a “on the left” of a∗ will be structurally
stable and similar to each other, and so will be the systems with control
parameters “on the right” of a∗ . But the systems on the left vs. on the right
are not structurally similar.
94
This can be summed up a bit loosely as when a control parameter passes
through a critical value, we observe an instantaneous change in qualitative struc-
tural properties of the DS. Such abrupt changes of system behavior are called
bifurcations.
Bifurcations occur in real-world systems everywhere. They are not confined to
systems described by ODEs. Here are two examples:
My old-fashioned wall-clock that I inherited from my grandmother will stop
ticking at the point when the metal spring that I should have wound up
(but forgot) unwinds to a critical slackness. The control parameter here is
the force of the spring, and the bifurcation is one between a system with
a cyclic “ticking” attractor (spring force > a∗ ) and a system with a stable
fixed point (spring force < a∗ , clock stands still).
A dripping faucet. This is my favourite example — a classic in the DS litera-
ture. You can try it out for yourself, but you need a faucet that can regulate
the water very finely, and the outlet nipple should be small — it doesn’t
work satisfactorily with bulky rusty old faucets.
Figure 40: A dripping faucet harnessed for the service of science (picture taken
from Martien et al. (1985))
The control parameter is the degree by which you open the faucet. When you
open it just a tiny bit, you will observe a regular, slow, periodic formation
of water drops coming out in equal time intervals: drip – drip – drip
... : a periodic attractor. When you slowly open the faucet a little further,
the dripping will get faster but stay regular periodic — until a bifurcation
suddenly appears: when the faucet is opened across a critical value, the
dripping sound will suddenly change to a double periodic pattern (your ears
will immediately notice the difference): drip-drop – drip-drop – drip-
drop ... : you witness a period-doubling bifurcation. Now continue (needs
95
a really finely operating faucet) and slooowly open the faucet further. All of
a sudden, the double periodic pattern will double again: drip-drop - dribble-
plop – drip-drop - dribble-plop – drip-drop - dribble-plop ... .
If you continue opening the faucet beyond another certain critical point,
suddenly all periodic regularity in the dropping rhythm will be lost, and
you have entered the chaotic regime of the faucet system (see Section 3.3.6).
The physical reality of this system is very complex, involving water volume
shapes, surface tension, physical interactions between the material of the
faucet nozzle and water, etc., and only approximate mathematical models
are within reach. This system has fascinated DS researchers a lot, and
numerous experimental studies paired with mathematical analyses revealed
an enormously complex picture (Martien et al., 1985; Tufaile et al., 1999).
Not even the most innocent water faucet can escape the bite of science.
we saw that this system changes its qualitative behavior if a is changed. In fact,
this RNN behaves like a perfect water faucet, going through a sequence of period
doublings and entering chaos if a is increased.
Now brace yourself for some disturbing news. RNN models have huge numbers
of control parameters: namely, every single synaptic weight can be regarded a
control parameter. In the simplest version of an RNN with L neurons, x(n + 1) =
tanh(W x(n)) with x ∈ [−1, 1]L , all the elements in the connection matrix W
are control parameters. It is clear that, if already an RNN as simple and small
as (55) can go through a splendid bifurcation cascade under the influence of a
single control parameter, much stranger things can happen in larger RNNs with
hundreds, thousands or (in modern Deep Learning RNNs) even billions of control
parameters.
During the process of training an RNN (we will learn the details of that soon
in our course) these parameters are all gradually adapted by gradient descent
algorithms akin to the ones that we saw in Section 2.2.3. One starts from some
randomly initialized RNN with a weight matrix W init . Usually the initialization
is done by setting all the parameters in W init to random small values. When all
connection weights in W init are very small (and the tanh sigmoid is used), the
96
resulting network dynamics is simple: all state trajectories converge to the all-
zero state x(n) = (0, . . . , 0)′ . The L-dimensional phase portrait of the freshly
created network shows a single point attractor. But at the end of the training
process, depending on the task that one wants to serve with the RNN, one wishes
to see a much richer phase portrait in the final RNN with weights W trained . For
instance, an RNN trained as a motor pattern generator for a legged robot should
host a number of different attractors for motor patterns like walking, running,
hopping, jumping, etc. In the gradual training process, the initial outfit of control
parameters W init is continuously adapted toward W trained . On this way W init →
W trained , the RNN must undergo bifurcations, because the final RNN dynamics is
not structurally similar to the simple point attractor dynamics of the initial RNN.
When the training processes tries to push the adapting weights through the critical
point of a bifurcation, the gradient of the loss function may exhibit discontinuous
jumps. This disrupts learning and can lead to sudden deteriorations (instead
of the expected gradual improvement) of the learning curve. Gradient-descent
training of RNNs may not converge (Doya, 1992). In today’s Deep Learning these
problems are resolved by replacing the RNN with an approximate feedforward
network during training, as we will see in Section 4.5.
3.3.6 Chaos
There are three basic kinds of attractors. The two simpler ones we have already
met: point attractors and cyclic attractors. Now I will cast a glimpse at the third
kind, chaotic attractors, also sometimes called strange attractors. The mathemat-
ical study of chaotic attractors has given rise to an entire subbranch of DS theory,
chaos theory. It is a very rich (and indeed strange) field, and here I can only
present a few basic concepts in a rather informal manner.
I will tell my story using the super-classical demo example, the Lorenz attrac-
tor. It is defined by the 3-dimensional ODE
ẋ = σ(y − x),
ẏ = x(ϱ − z) − y, (56)
ż = xy − βz. (57)
97
20
A C
x1, x2, x3
10
−10
45
40 −20
0 2 4 6 8 10 12
35
2
D
log10(|Dx|)
30 20
z
1
25
z
10
0
20
0 −1
15
−2
10 −10
5 x
x
−3
30 20 10 −20
0 −10 −20 −30 −4
0 2 4 6 8 10 12
yy
t
20
B
10
−10
−20
0 10 20 30 40 50 60 70 80 90 100
then for some reason loses hope and switches to the other bulb... and so on and
on.
This attractor is, like all attractors, defined as a subset of points A ⊂ R3 . The
geometry of this attractor set A is indeed ... strange. It consists of infinitely many
curved lines which fold into the volume cube in an infinitely finely organized way
that defies every plotting resolution. If I would have run the simulation shown in
Figure 41A infinitely long, I would have gotten a line of infinite length contained
in A which would never exactly meet itself — and this infinite line would only
be one among infinitely many other such lines, all belonging to the set A. The
geometry of A cannot be described with the familiar concepts of geometry. It needs
a new sort of geometry theory, the theory of fractal geometry. The study of fractal
geometries was boosted by the graphics powers of modern digital computers and
gives rise to beautiful pictures — and allows computer game engineers to generate
virtual sceneries full of trees, clouds and rippled waters, not to speak of broccoli,
familiar objects around us all of which have fractal geometries.
98
The defining hallmark of chaos is however not its fractal geometry but its ef-
fective unpredictability. Consider the plot of 200 time units length of the x state
component shown in Figure 41B. The spiraling rotations in the two lobes of the
attractor are nicely visible as “upper” and “lower” oscillation episodes. The num-
ber of revolutions that the trajectory takes within one of the lobes seems random.
In the shown plot, picking the “upper” oscillations, we count 4 - 1 - 5 - 1 - 2 -
1 - 8 - 1... In fact, using such apparently random features of chaotic trajecto-
ries, one can design random number generating algorithms, and chaotic dynamics
have been proposed for data encryption schemes. But wait... isn’t the evolution
of an ODE system deterministic? Yes it is. If one fixes the initial starting point
(x(0), y(0), z(0))′ , the future evolution ((x(t), y(t), z(t))′ )t∈[0,∞) is uniquely and per-
fectly determined, without any randomness. And here is the catch: one cannot
effectively determine the starting (x(0), y(0), z(0))′ with the infinite precision that
real numbers need. In all practical computations, one only has finite-precision
numerics; and in all empirical measurements one can measure a physical quantity
only up to a residual measurement error. This necessarily remaining indeterminacy
of the initial condition, regardless how small, renders the long-term continuation
of the trajectory effectively unpredictable. This phenomenon is brought to the
surface in the plots in Figure 41C,D. Panel C shows three x(t) trajectories that
were started from almost, but not quite, the same starting states. The three
starting states were about 0.004 space units distant from each other, that is about
1/10,000 of the diameter of the entire attractor. Up to a simulation time around 5
time units, the three trajectories stay so close together that their separation from
each other is not visible in the plotting resolution. Thereafter, they visibly sepa-
rate from each other and continue in apparently entirely individual futures. Panel
C shows how the pairwise distances between the three shown trajectories grows
as time goes on. Note the log scale of the plotted distances: the mutual distance
grows exponentially on average (red shaded arrow) until complete separation (blue
shaded arrow) is reached. This exponential separation of nearby trajectories is the
essential fingerprint of chaotic dynamics.
Many real-world systems are apparently chaotic. A large class of real-world
systems whose chaotic properties have been intensely studied are turbulent flows in
fluids and gases. Turbulence in flows manifests itself in many ways, the most visible
being the creation of eddies (local “whirls”) of all sizes. The earth’s atmosphere
and oceans are turbulent flows. This means that long-term weather forecasting
is impossible: a complete, precise determination of the global atmosphere’s state,
down to every molecule, is obviously impossible. Here is a thought experiment, the
famous butterfly effect, you will probably have heard of it: consider two complete,
totally precise simulations of our planet’s atmosphere, both started from almost
the same starting state. The only difference between the two starting states is that
in one of the two simulations, a butterfly somewhere in Chine flaps its wings once,
and it doesn’t do that in the other simulation. This tiny difference will blow up,
exponentially invading the future, and lead to entirely different weather conditions
99
after a few weeks, with a tornado striking Groningen in one simulation and calm
sunshine in the other.
Chaos is easy to create in continuous-time RNNs ẋ = tanh(a W x) or in
discrete-time RNNs x(n + 1) = tanh(a W x(n)) when the recurrent weight ma-
trix W is fixed and the single scaling control parameter a is increased. Intuitively,
when a grows beyond some critical value, the network-internal re-circulation of
activation gets so strong that the network over-excites itself and “goes crazy”.
Figure 18 shows that this can happen in RNNs with merely two neurons.
Chaos apparently plays important roles in biological brains. The respective lit-
erature is bottomless. Chaos has been detected or hypothesized on all levels of the
anatomic and functional neural processing hierarchy, from a single spiking neuron
to daydreaming and schizophrenia. In two successive articles, Korn and Faure
(2001) first give a methodological introduction to chaos concepts and the difficul-
ties of its experimental verification (more detailed than this lecture note section
— very recommendable reading), and then survey the wide range of theories and
findings (Korn and Faure, 2003).
Signal processing and control. This large, classical and powerful engineering
100
discipline has developed its own rigorous methods arsenal to deal with time-
series, which are called signals in that field. A main tool there is to de-
compose signals into frequency mixtures by Fourier or Laplace transforms,
and do all analyses in this transformed frequency domain. These methods
are broadly used in theoretical neuroscience and whenever artificial neural
networks get into the hands of engineers.
Stochastic processes. I only briefly pinpointed stochastic DSs by the discrete-
time, discrete-value models of Markov chains, hidden Markov models and
POMDPs. But stochasticity is everywhere in biological neural networks and
is present in almost all real-world signals that one encounters as inputs in
practical applications of RNNs. At some point, an RNN modeler will need
to apply general methods of stochastic process theory.
Information theory. Artificial and biological neural systems use their dynamics
to “process information” in some way — that’s what they are made for. This
will sooner or later force an RNN modeler to read a textbook of information
theory.
Network theory is a rather young, interdisciplinary research field which con-
siders dynamics that evolve on graphs (called networks in that context, not
to be confounded with neural networks — “networks” is a far more general
concept). Many important real-world systems can be best modeled as a
collection of active elements which exchange signals, messages, materials or
forces along connecting channels. Standard examples are power grids (active
elements: power plants, transformation stations, end-users; connecting chan-
nels: power cables); social networks (active elements: people; connecting
channels: communicated messages); metabolic networks (active elements:
proteins and other biomolecules in living cells; connecting channels: chem-
ical reaction pathways); or, yes!, neural networks. The specific modeling
powers of network theory come from the combination of graph theory (the
global connectivity patterns between the active elements) with DS theory
(the local dynamics of the active elements), which opens new ways to study
the emergence of global, network-wide dynamical phenomena.
Non-autonomous dynamical systems. I presented ABC theory only for ODE
systems that are not driven by external input. But neural networks normally
operate by processing input data (recall Equations 37-39). Input-driven DS,
which are a case of what mathematicians call non-autonomous DS (because
their dynamics cannot be understood only from within the system itself),
are quite hard to analyze. In particular, developing adequate generalizations
of ABC phenomena that work in input-driven scenarios is a mountain that
mathematicians have only begun to climb. This is in my view the most
painful gap in today’s math toolboxes when one wishes to really understand
brains and artificial RNNs.
101
4 Recurrent neural networks in deep learning
The deep learning (DL) revolution began with feedforward neural networks, es-
pecially CNNs. Wait... that is not quite true; history isn’t that simple. The
pioneering paper of Hinton and Salakuthdinov (2006), which I (and others) con-
sider a kickoff event for DL, was about a special sort of RNNs, the restricted
Boltzmann machine (RBM). We will treat RBMs later in this course. In that pio-
neering work, the powers of RBMs were demonstrated on data compression tasks.
Only in passing it was also mentioned that RBMs can be used to initialize the
backprop learning for MLPs — the rest is deep learning history. After creating
excitation and much respect (because it isn’t easy at all to work with RBMs) in
the early years of DL, RBMs vanished from the focus of attention, and today are
not playing a leading role in DL anymore.
When it comes to recurrent neural networks, the DL stage is nowadays reserved
for Long Short-Term Memory (LSTM) networks and their derivatives. Like con-
volutional neural networks, LSTM networks have a history that began before the
DL revolution. Essential ideas were introduced already in the diploma thesis of
Hochreiter (1991) and became fully worked out in Hochreiter and Schmidhuber
(1997). However, LSTM networks only rose to the domineering role that they have
today after reliable backpropagation algorithms became available through DL re-
search. In our present time, virtually all advanced machine learning solutions that
are based on recurrent neural networks are LSTM networks or close relatives of
them.
In order to understand LSTM networks, one has to make friends with two sepa-
rate algorithmic and architectural techniques: (i) a generalization of the backprop
algorithm from feedfoward networks to recurrent networks, called backpropagation
through time (BPTT), and (ii) a special kind of complex neuron model, the LSTM
unit. I will treat both topics in turn, but before I do that, I will explain what it
means to carry out a supervised learning task with temporal data. I restrict the
treatment in this section to discrete-time dynamics, where time is represented by
integers . . . , n − 1, n, n + 1, . . .. RNNs in machine learning on digital computers
always use discrete time. In contrast, biological neural systems, and some of the
recent developments in neuromorphic computing which I will hint out in the last
session of this course, are based on continuous time.
where
102
• n = 0, 1, 2, . . . , nmax or n = 0, 1, 2, . . . are the time steps, which may run
until an end time nmax or forever,
• Equation 58 specifies how the network activation state is updated from one
timestep to the next, and Equation 59 specifies how the output vector is
computed at time n,
• x(n) is the vector of activations of the neurons inside the RNN — I will
generally denote the number of neurons in an RNN by L, so x(n) ∈ RL ,
• W is the L × L matrix containing the synaptic connection weights wij ∈ R
giving the strength of the connection from neuron j to neuron i,
• W in is an L × K dimensional real-valued matrix containing the weights from
K input neurons into the RNN — I will generally denote the dimension of
the input signal u(n) by K,
• b ∈ RL is the vector of biases,
• σ is a sigmoid function (like in MLPs, typically the tanh, the logistic sig-
moid, or the rectifier function), which is applied element-wise on the vector
W x(n − 1) + W in u(n) + b,
• y(n) is the output signal generated at time n — I will generally denote the
dimension of the output signal by M ,
• W out is an M × L sized output weight matrix containing the weights of
synaptic connections from RNN neurons to output neurons,
• f is a wrapper function applied to the linear “readout” W out x(n); often
again a sigmoid, but often also just the identity function.
Sometimes one wishes to feed back the generated output signal y(n) into the
RNN. Equation 58 is then extended to the format
where W fb is an L × M sized matrix for feeding back the generated output signal
y(n) into the RNN. Such output feedback is required when the task is not an input-
to-output transformation where the output is determined by the input signal, but
where instead the output is actively generated by the RNN, even in the absence
of input. Such tasks occur specifically in robot motion generation.
In order to get the dynamics (58), (59) started, at time n = 0 the recurrent
network state x(0) is set by the user to some predetermined initial state, often
the all-zero state. For later timepoints, Equation 58 is used to compute network
states x(1), x(2), . . ..
Figure 42 illustrates the wiring schema of such a basic RNN. Many variations
of the system equations (58), (59) are in use. Also, more complex architectures
103
Activation signal
of an internal
neuron xi(n)
Wfb
Figure 42: Schema of a basic RNN. Here an example is shown that has one-
dimensional input and output signals, and hence only one input neu-
ron and one output neuron. This particular RNN has been trained
to generate a sinewave output whose frequency is controlled by the
input: high input value gives high frequency output. Output feed-
back is needed here. Activation signals from a few internal neurons are
shown. Image taken from [Link]
Echo_state_network.
than the one shown in the figure are used. For instance, in hierarchical RNNs,
several recurrent sub-modules are stacked on top of each other. Or the signals
from inside the RNN are propagated through a MLP in order to obtain a more
complex “readout function” than what you can get with (59).
104
where x(n) is the data point at time n. Right-infinite timeseries start at time 0
and continue (potentially) for an endless future duration:
105
phonological frequency spectrum that is characteristically different at differ-
ent positions in the word (speech recognition operates frequency spectrum
signals). The frequency spectrum of the “f” phoneme at the beginning of
“faster” is different from the phonological spectrum of the “r” at the word’s
end, and all intermediate phonemes likewise have their individual spectral
fingerprints. Training data for keyword recognition systems consist in large
numbers of short recordings of individual words. For instance, in the robot
control scenario, the word “faster” would be spoken by many speakers many
times, and each of these utterances yields a new, short timeseries consisting
of spectral measurements. Similarly, the other keywords of interest would
be recorded from a large number of human speaker utterances, each of which
gives another short training timeseries. Different from the stationary mod-
eling case, when you get a datapoint that shows the “f” frequency spectrum,
you can infer that it is located at the beginning of an utterance of “faster”,
or toward the end of an utterance of “lift”. The temporal position of local
observations carries important information, and a given sort of datapoint
(like the “f” spectrum reading) can be found in the data at different tempo-
ral positions with different probabilities. Training datasets for nonstationary
tasks necessarily consist of large numbers of short timeseries, each of which
is ’nonstationary’ in the sense that it looks characteristically different at dif-
ferent times in its brief realization. Formally, such training datasets would
need double indices, one for time and the other indicating the individual
short sequence:
(ui (n), yi (n))i=1,...,N ;n=0,...,nimax ,
where n indicates time, and i the individual training example. Note that
different examples i may have different durations nimax .
The stationary vs. non-stationary distinction is not clear-cut. Whether a sys-
tem is considered stationary or non-stationary depends on how long it is observed.
For instance, when a pair of dancers performing a waltz is monitored for an ex-
tended time, the different waltz figures repeat and the signal can be considered
stationary. However, the entire waltzing performance can be segmented into short
figures (which are individually trained in dancing lessons), where each figure con-
sists in a specific non-stationary sequence of motions. The same holds for speech
(individual words are non-stationary, but a long monologue of a boring person who
repeats his jokes again and again can be considered stationary), weather (during
a day there will be a non-stationary evolution of temperature and wind, but over
the years everything repeats — ignoring climate change), and many other inter-
esting signals. In order to learn RNN models for such mixed signals (short-term
non-stationary, long-term stationary) one may pursue two strategies: either try to
learn one large, complex RNN which is capable of changing its “mode” in shorter
time intervals; or try to learn several smaller, simpler RNNs, one for each typical
non-stationary subpattern, and then combine them into a compound multi-RNN
system where the individual RNNs are activated / deactivated in turn.
106
Machine learning set-ups and architectures differ substantially for stationary
versus nonstationary tasks. Not only are the formats of training data different,
but also the learning mechanism must be adapted to the stationary / nonsta-
tionary character of the task. In the latter case, the learning mechanism must be
sensitive to the temporal location of incoming new datapoints in a sequence. Solv-
ing nonstationary learning tasks is typically more difficult than solving stationary
ones.
Each of these input windows is paired with a single output data point, often
(but not necessarily) from the last time point of the respective window, obtaining
input-output pattern training pairs
107
A consequence of this window-based input feeding is that a trained model can
never incorporate memory effects longer than l steps. If an output value y(n)
has some dependency on inputs u(n′ ) that came earlier than the window (that is,
n′ < n − l + 1), then this modulating impact on the output y(n) cannot be realized
by the trained network. If one needs to account for long memory effects, one must
choose long windows. The required window lengths can become very large in some
technical applications. My favorite example are so-called channel equalization fil-
ters, which essentially are perceptrons that sit close to the receiving antenna in
the high-frequency sector of a mobile phone. Their task is to disentangle the
various ’echos’ in the incoming radio waves, which result from the waves taking
several different travel paths from the transmitter to the phone. These chan-
nel equalization perceptrons have window lengths of tens of thousands. If you
are interested — these filters are described in my machine learning lecture notes
([Link] in Section 11.2.
For timeseries prediction tasks, the output data point comes from the next
timepoint after the window:
This is how the currently popular large language models (like GPT-3) are set
up. They are modeled with feedforward networks. Each input data point u(n) is a
vector (a semantic word embedding vector) which encodes a word (more precisely,
a ’token’ — which is almost the same as a word) from the model’s vocabulary. The
input window is a segment of a text, consisting of l consecutive word embedding
vectors. The target output in the training data is (a vector encoding of) the word
(more precisely the token) in the text that comes next after the input window
text segment. I don’t know for sure, but I have heard from reliable experts that
GPT-3 uses a window length of 200. Thus, a word in a text generated by GPT-3
will typically make much apparent sense in the context of the previous 200 words,
but may fail to blend well into what was expressed in the generated text earlier.
Now we turn to RNNs. Recalling Equation (58), we find that RNNs do not take
windowed inputs but only a single input data point at a time. This is visualized in
Figure 43 (right). At each internal state update x(n − 1) → x(n), the RNN inputs
the single data point u(n) and generates the single output data point y(n). Despite
this single-point input, RNNs can in principle have unbounded long memory, as
we will discuss in the next subsection 4.4. The input sequence must be fed to
the RNN in exactly the temporal order of the input sequence, because otherwise
the memory information from preceding inputs, which are encoded in the current
network state x(n), will get confounded. This is different from the operations of
108
input input point
window of at time n
length l
Figure 43: Input-output set-up for feedforward and recurrent networks. This
schematic shows the case where input and output datapoints u(n) and
y(n) are each single real numbers. Left: A feedforward network with
l input units transforms input windows of length l to a single out-
put point, here placed at the last time of the input window. Right:
An RNN inputs one datapoint u(n) and outputs a datapoint y(n) at
the same time. When the RNN state is updated (Equation (58)), the
next input point u(n + 1) becomes the new input, and the temporal
successor output point y(n + 1) is outputted.
109
tion of a spoken word depends on the previously spoken context, and an intelligent
aircon controller will adapt the gain (response strength) of the power setting to
the swiftness in the room temperature response which the controller has observed
earlier. Figure 44 illustrates how ’memory’ in an RNN lets output y(n) at some
time n be influenced by earlier input u(n′ ), where n′ < n.
input u(n)
RNN
outputs y(n)
Figure 44: Demo of memory effects in RNNs. A ’spike’ in the input signal may
lead to long-term output responses of various sorts, for instance another
single spike that signals the ending of a (trainable) waiting time, or
long-lasting, slowly fading, complex ’repercussion’ responses that set
in directly after the trigger input.
How can an RNN realize such memory effects? How can it preserve input
information from the past to make it affect the present output? The answer is:
past input leaves its traces (’encodings’) in the current network state. This is
very simple maths, as becomes clear when we expand the network state update
equation (58) a little (we drop the bias to make equations shorter):
x(n) =
= σ(W x(n − 1) + W in u(n)) (60)
in in
= σ(W (σ(W x(n − 2) + W u(n − 1))) + W u(n)) (61)
in in in
= σ(W (σ(W (σ(W x(n−3)+W u(n−2)))+W u(n−1)))+W u(n)) (62)
= ...
110
To get these formulas, all that I did was to re-insert the expression for x(n)
again and again, with the right timeshifts applied. For instance, in order to expand
(60) to (61), I replaced the x(n − 1) term in (60) by σ(W x(n − 2) + W in u(n − 1)),
etc. This expansion can be iterated as often as one wants, leading to arbitrarily
deeply nested expressions, which reveal that a network state x(n) carries in it
some traces of all earlier inputs u(n), u(n − 1), u(n − 2), . . .. These traces can be
interpreted as memory traces of the input. In non-mathematical intuitive words,
one could also say that an input pulse u(n) inserted into the RNN at time n will
continue to ’circulate’, ’reverberate’ inside the RNN and become a part of future
states x(n + 1), x(n + 2), . . ..
This sort of input memory effect is not restricted to neural networks - it occurs
in all dynamical systems, whether they are discrete-time iterated maps x(n +
1) = T (x(n), u(n)) or continuous-time ODE systems ẋ(t) = T (x(t), u(t)) (as we
introduced them in (37) and (39)).
This kind of memory ’washes out’ with time — in Equation (62), the contri-
bution of u(n−2) to the state x(n) is smaller than the contribution of u(n−1),
which in turn has a smaller impact on x(n) than the current input u(n). It is good
(and even physically correct) to think of these fading-out input traces as “echos”
or earlier input, and the RNN as a sort of echo chamber. Like real acoustic echos,
these traces fade out with time — in fact they must do so in order to ’make place’
for the freshly incoming information from the next inputs.
There are many sorts of ’memory’ in dynamical systems. It is worth a few
words to take a look that the dynamical systems that are us — our brains. Neuro-
scientists have identified a whole range of neural mechanisms that serve different
memory tasks on different timescales. On the fastest timescale, our brain state dy-
namics gives us the same fading dynamical memory that every dynamical system
has. You can experience it for instance in the visual domain when you cast a quick
look to the side and immediately direct your gaze forward again: for a few tens of
milliseconds you can quite completely recall what your quick sideward glance has
let you see. This state-based, fading dynamical memory is one sort of short-term
memory in brains. Neuroscience has found also other sorts of short-term memory,
thus this is an umbrella concept.
On a longer timescale our brain employs a quite involved set of structures and
mechanisms to achieve working memory. The classical (pre-cellphone era) exam-
ple is to look up a phone number in a phonebook in some room A, then go to
another room B where the phone is installed, and do the call. Doing this requires
the caller to memorize the looked-up number while going from room A to room B.
Humans do this by an active, conscious process — they repeat the number in their
mind (maybe even with moving lips) until it can be used for dialling. This method
of cylcic rehearsal is a classical model of cognitive psychology (Baddeley, 2003).
Characteristic properties of working memory are that only a finite number of dis-
crete memory items can be kept ’alive’ in working memory (you may have heard
of the famous “7 plusminus 2” working memory capacity hypothesis), that stor-
111
ing, preserving, and recalling are active and conscious processes (which require a
complicated neural machinery for attention, evaluation, encoding, decoding), and
that the lifespan of an actively stored memory item is unbounded in principle (if
you concentrate on it, you can keep a phone number in your working memory
for hours). Mathematical models in theoretical neuroscience often use attractors
to explain working memory functionalities: according to these models, when you
actively put an item into your working memory, an attractor is generated in your
brain-dynamical-system, and this attractor encodes the memory item. In the case
of cyclic rehearsal, this would be a periodic attractor. However, empirical neu-
roscience cannot tell us yet how, exactly, working memory is neurally organized,
and theoretical neuroscientists cannot tell us how to best model it in terms of
dynamical systems, and machine learners cannot tell us what is the best way to
endow artificial RNNs with working memory. This is all a super relevant and
fascinating collection of open questions, because for artificial general intelligence
one needs working working memory mechanisms. The most sophisticated models
that I am aware of are the neural Turing machine models of Alex Graves and
colleagues, where RNNs are paired with a new kind of neural processing modules
that roughly correspond to the working tape of a Turing machine (Graves et al.,
2014, 2016). Some of my own research was about working memory mechanisms
in RNNs in a machine learning setting (Pascanu and Jaeger, 2011). The simplest
model that I am aware of are the “long short-term memory” (LSTM) systems,
which are the current standard method in deep learning and which we will treat
in more detail later in this section.
Beyond short-term and working memory, biological and artificial RNNs also
have a sort of memory on the very long timescale of learning. When you ex-
perience a memorable life event, your brain encodes and stores it in a dry and
safe place where you can retrieve it from at any later time in your life (well, no!
humans do not retrieve what they have stored in memory, they recall from their
memory altered versions of past experiences and the alterations can be very strik-
ing — as documented in a classical, famous long-term study of psychologist F.
C. Bartlett who (using personal friends and colleagues as subjects for his study)
tracked the ever-changing versions of recalled memory items over several decades
(Bartlett, 1932). The memory items that are stored in this (lifetime) long-term
memory in RNNs are not stored in the dynamical activation patters like short-
term and working memory items are. Instead, when new items are inserted into
an RNN’s long-term memory by some learning process, the very structure of the
RNN is changed. For more than half a century, it has been firmly believed (both
in neuroscience and machine learning) that these structural changes amount to
permanent changes in synaptic connection weights. This is still the unquestioned
view in machine learning — after all, training an artificial neural network means
to compute its synaptic weights. However, in neuroscience, the view that memory
items (and other knowledge) is encoded in stable synaptic weight patterns has
recently become quite doubtful. With new, highly sophisticated imaging methods
112
it has become possible to “videotape” individual synapses in live, freely behaving
animals (mice) over weeks. Disturbingly, it turns out that synapses grow, shrink,
and change their position on the receiving neuron quite fast, often within hours
(Deitch et al., 2021). It seems that our brains are incessantly rewiring themselves...
at least in some places. If that impression will become confirmed in the next years,
the question how we encode and store long-term memory items in our bio-brains
becomes a big mystery again. Stay tuned.
Some kind of short-term memory is needed for most temporal processing tasks,
and if one wishes to have high-performing neural network solutions for these tasks,
one must make sure that the architecture that one uses can support the right sort
of short-term memory. Designing clever, trainable memory mechanisms is a major
enabler for neural networks in temporal machine learning tasks. This becomes
strikingly clear when we take a look at ChatGPT and other large language models.
They draw their amazing performance from transformers (Phuong and Hutter,
2022), which can be regarded as a working memory mechanism paired with an
attention mechanism, where the latter decide which items in past (text) input
should become stored in working memory. Interestingly, transformer networks
are feedforward networks, not recurrent ones. They emulate temporal processing
of incoming text sequences by windowing. Their maximal memory span is thus
limited by the chosen window length.
113
The training data for pattern generation tasks consists just of the desired
output sequence, S = (y(n))n=1,...,nmax . The objective is to train a RNN such
that its output signal ŷ(n) is “similar” to the training signal.
It is not obvious however how “similarity” is measured when there is no input
in the testing phase. When the training signal S = (y(n))n=1,...,nmax is chaotic
or stochastic, the signal generated by the trained RNN ŷ(n) can not directly
be compared with a “known correct” output. Figure 45 illustrates this dif-
ficulty. The teacher pattern is a chaotic signal generated by solving the
Mackey-Glass equation (an equation describing the temporal change of the
amount of white blood cells in leukemic patients (Mackey and Glass, 1977);
this has become a popular benchmark signal for RNN training). There is
no general recipe for quantifying the accuracy of pattern generation learn-
ing. The only way is to define measurable characteristics of the signal (for
instance mean amplitude, Fourier spectra, or conditional probabilities for
continuations of the signal given an previous few values) and compare the
values of these characteristics measured in the original training signal vs.
the generated ones.
Pattern detection tasks. Sometimes one wants to identify certain patterns when
they occur within an ongoing observed process. For example, in cardiological
monitoring, one wants to detect moments when a certain pathological ab-
normality appears in a ECG signal; or in speech-based car control, the car’s
driver voice recognition system must identify moments when the driver says
“brake”; or the control system of a nuclear fission reactor must trigger an
alarm action when the plasma shows signs of instability (a very tricky task
where machine learning methods are considered a key for making nuclear fu-
sion based energy possible — query “nuclear fusion plasma disruption neural
network” on Google Scholar!).
The training data here often is prepared in the format S = (x(n), y(n))n=1,...,nmax ,
where x(n) is a signal of measurements from the system in which one wants
to identify a certain pattern, and y(n) ∈ {0, 1} is a binary indicator signal
which is 0 when the pattern is not present and which jumps to 1 when the
pattern occurs. This basic scheme comes in a number of more sophisticated
variants. For instance, the indicator flag jumps to 1 only at the end of the
pattern, because sometimes a pattern can be reliably identified only after it
has come to its end — e.g. identifying the word “stop” in a speech signal
can only be done when the “p” has been uttered, in order to not confuse this
word with “stomach”, “stolid”, “stochastic”, etc.). Or the indicator jumps
to 1 a certain time before the target pattern occurs; this is the functionality
needed for giving early warning of an impending system malfunction. A
point in case is the early warning of an approaching epileptic seizure from
EEG recordings (Figure 46) — query Google Scholar with “seizure predic-
tion neural network” to see how much research has been spent on this theme,
114
1
0.5
0.5
0.5
0
0 20 40 60 80 100 120 140 160 180 200
0.5
0
0 0.5 1 0 0.5 1 0 0.5 1
Figure 45: Illustrating the difficulty of assessing the accuracy of a learnt pat-
tern generator. Green: Original (training) pattern. Magenta, blue:
patterns generated by two trained RNNs. Top three panels: pattern
signal against time. Bottom: time delay embeddings. The blue pat-
tern generator matches the teacher better than the magenta one. This
comes out visually clearly in the delay embedding plots, while it would
be hard to judge the matching accuracy from the plots against time.
which cannot be considered solved even today; see Litt and Echauz (2002)
for a very instructive survey.
Timeseries prediction tasks. Sometimes one wants to look into the future.
This task is particularly popular among people who want to get rich by
forecasting financial timeseries, and also among students searching for thesis
topics which will help them to become the kind of people who get rich. Other
applications exist, too, for instance weather forecasting, windspeed predic-
tion for wind energy farming (a theme with an extensive machine learning
literature, check out Google Scholar on “wind speed prediction neural net-
work”), or in epidemology where one wants to understand the spreading
dynamics of a pandemic.
The concrete kind of training data and best suited RNN architectures vary
115
Figure 46: The seizure prediction scenario. The input is a multi-channel EEG
recording. The task is to predict the seizure on the basis of an earlier
observation (gray block). Image taken from Litt and Echauz (2002).
116
The training data consists in a number of example timeseries
S = (u(i) (n), y(i) (n))i=1,...,N ;n=1,...,ni = (d(i) (n − 1), d(i) (n))i=1,...,N ;n=1,...,ni ,
where i is the index of the training example and d(i) (n) is the n-th value in the
i-th training timeseries (there were 111 such timeseries in the example task
illustrated in Figure 47). The output timeseries desired from the network is
thus equal to the input timeseries shifted one step ahead. An RNN is set
up with a single input and a single output node. After training (when the
network has learnt to predict all the training series as well as it can), the
network is used to forecast the future by iterated one-step prediction in the
following way:
The last obtained output ŷ(k) = f (W out x(k)) is saved. This should be
a prediction of the next timeseries value, ŷ(k) ≈ d(k + 1).
Phase 2: iterated prediction: Feed network outputs back as input, that
is, compute for n = k + 1, . . . , k + h
For a serious effort to obtain good predictions, this basic scheme must be
considerably refined (check out Ilies et al. (2007) for a case study). But that
is a common machine learning wisdom: on any nontrivial learning problem,
better and better results are achieved by investing more and more work,
insight, and computing power.
117
Figure 48: Modeling the response of an aircraft (top three panels showing the
aircraft accelerations in the three spatial dimensions) to the pilot’s
action on elevator, aileron and rudder (bottom). The neural network
predictions (red lines in upper three panels) show a good agreement
with actual aircraft flight dynamics (blue lines, measured from actual
flight experiments). Figure taken from Roudbari and Saghafi (2014).
input. When the target system behavior is very nonlinear and has memory
effects, RNNs are a good candidate to yield simulation models. Figure 48
shows an example of an aircraft model realized by an RNN.
Other intensely investigated examples of modeling the responses of com-
plex systems to input are oceanic flows, local weather dynamics (and more
generally, all sorts of turbulent fluid dynamics), all sorts of industrial manu-
facturing pipelines, robot motion, power grids, etc., etc. — basically, every
complex system that one wants to control, understand, or predict through
simulations is a candidate for RNN modeling.
The training data here are of the generic format S = (u(n), y(n))n=1,...,nmax
or S = (u(i) (n), y(i) (n))i=1,...,N ;n=1,...,ni , depending of whether one deals with
stationary or nonstationary system behavior. These data are obtained either
from physical measurements, or from (expensive) first-principle simulations.
In the aircraft modeling example, u(n) would be the three pilot command
signals (for controlling the elevator, aileron and rudder of the aircraft) and
y(n) the measured accelerations of the aircraft.
118
data sources are temporal; no wonder that our brain, the best world modeling
engine that we know, is recurrent.
119
For MLPs, step 1 is done with the backpropagation algorithm. This algo-
rithm crucially depends on the fact that the network topology is unidirectionally
feedforward, that is, there are no connection cycles.
For RNNs, finding weights which minimize the training error is more difficult.
A number of algorithms are known which are based on different mathematical
principles. The “best” of these algorithms does a true gradient descent on the
performance surface and is suitable for online learning (adapting the RNN con-
tinuously while a never-ending stream of training data is arriving). This real-
time recurrent learning (RTRL) algorithm has been introduced 40 years ago by
Williams and Zipser (1989). Sadly, RTRL is too expensive for most practical ex-
ploits, having a cost of O((K + L + M )4 ) per update step. Yet, interest in RTRL
has recently been rekindled in a subfield of deep learning called continual learning,
where the objective is to find learning schemes that enable “life-long” training of
RNNs. If you are interested, a tutorial introduction to RTRL and other RNN
training algorithms is Jaeger (2002).
Today, the training error minimization task (63) is almost always solved by the
backpropagation through time (BPTT) algorithm. The idea is to “unfold in time”
a recurrent neural network into a feedforward neural network, by assigning an
identical copy of the RNN to each timestep and re-wiring the internal connections
(in the weight matrix W ) such that they feed forward into the next copy in time.
Figure 49 illustrates this idea.
…
x(n+1)
u(n+1) y(n+1)
u y
u(n) x(n) y(n)
x
x(n –1)
…
Figure 49: Basic idea of the BPTT algorithm: a recurrent net (left) is identically
replicated for every time step and network-internal connections are
rewired forward to the next copy (right).
120
The BPTT scheme looks straightforward, but it comes with its own new prob-
lems.
The first problem is that the stack of temporal copies (right side in Figure 49)
must have a finite depth to enable the backpropagation algorithm. Call this depth
h (for “horizon”). The unfolded RNN then will be a finite, cycle-free network which
yields an input-to-output mapping from input sequences u(n), . . . , u(n + h − 1) to
output sequences ŷ(n), . . . , ŷ(n + h − 1). If the training data has a length nmax
of more than h steps, the training algorithm involves an averaging over gradient
increments collected from shifting the depth-h network forward through nmax −h+2
length-h segments of the training data ((u(n), y(n))n=j,...,j+h−1 )j=0,...,nmax −h+1 . This
implies that any effects from input u(n) on the output that are delayed by more
than h timesteps cannot be learnt. (Note. This is not entirely and always true and
the story to be told here is involved and certainly not fully understood. If there are
statistical dependencies between the input signals at different times which span
further than the training horizon h, some information about input u(n−h−g) from
times earlier than n−h+1 is encoded in the inputs u(n−h+1), . . . , u(n) available
to the learning algorithm, and effects from these encoded earlier on the current
output can be learnt by BPTT. In the pure memory learning task that I will use
below for demonstration purposes however there are no statistical dependencies
between inputs at different times, and the memory limit of h is strict.) This
is a significant problem in applications where there are input-output effects with
arbitrary long delays — dynamics of that sort cannot be perfectly learnt by BPTT.
Many real-world dynamical systems have this property. Examples are all turbulent
fluid flows (where “eddies” of increasingly large sizes impose effects on increasingly
large timescales), or human language (where what is said now in a text may have
an influence on the text meaning at arbitrary distant points in the future).
For all we know, the human brain solves this multi-timescale problem by ex-
ploiting a host of different physiological and architectural mechanisms which yield
a large compendium of different memory mechanisms, supported by different neu-
ronal circuits and physiological effects, ranging from ultra-short term memories in
the millisecond range, over a spectrum of short-term and working memory loops
from seconds to minutes to hours, to long-term memory mechanisms that oper-
ate in the range of the human lifetime. Neuroscience, cognitive science, AI and
machine learning so far has only given us a very partial understanding of such
memory subsystem cascades. This is an active, interdisciplinary research area in
which I am personally involved (check out the EU project MeM-Scales, “Mem-
ory technologies with multi-scale time constants for neuromorphic architectures”,
[Link]
The second problem is known as the vanishing gradient problem. It also oc-
curs in a mirror version as exploding gradient problem. This problem comes to
the surface when error gradients are back-propagated through many layers. The
deeper a network, the more serious the problem. Unfolding RNNs in time tends
to end up with a particularly large number of time-slice layers (up to several hun-
121
dreds), making these time-unfolded networks much deeper than commonly used
deep feedforward MLPs. Thus vanishing gradients are particularly disruptive in
RNN training.
I will demonstrate the mathematical nature of vanishing gradients with a super-
simple RNN. It has a single input unit with activations u(n), a single output unit
with activations y(n), and the main RNN itself consists of a single linear unit
with activation x(n), and no bias. This leads to the following embryonic RNN
equations:
x(n + 1) = w x(n) + win u(n + 1), (64)
y(n) = wout x(n),
where x, y, u, w, win , wout are just numbers, not vectors. The task which we will
consider is a pure memory task. The training data is given by (u(n), y(n))n=0,...,nmax =
(u(n), u(n − h))n=0,...,nmax , that is, the desired output y(n) is the input from h
timesteps earlier. The input signal u(n) is an i.i.d. signal sampled from the uni-
form distribution on [−1, 1], that is, at every time n the signal value u(n) is freshly
randomly chosen from the interval [−1, 1].
We now want to find values for w, win , wout which minimize the quadratic loss,
by running a BPTT gradient descent with an unfolded version of our embryonic
RNN. We unfold the RNN to depth h, that is we create h + 1 time-slice copies,
which for this task is the minimal needed depth.
For this simple RNN all gradients can be computed easily without invoking the
backpropagation algorithm. During the gradient descent, the unfolded network is
aligned with length-h windows of the input/output teacher signal, from times n to
n + h. The only task-relevant signal propagation pathway goes from the teacher
input u(n) at time n to the teacher output y(n+h) = u(n) at time n+h. All other
input-unit to output-unit pathways in the unfolded network lead an input signals
with an uncorrelated output teacher signal, and their contributions to the gradient
will average to zero as the unfolded network is moved forward through successive
time windows [n, n + h] of the training data. The only gradient component that
may not average to zero is the one on the pathway from the earliest input node
to the latest output node in the unfolded network (see Figure 50).
For the quadratic loss, this gradient component of the loss with respect to the
weigth w is equal to
∂
(y(n + h) − ŷ(n + h))2 =
∂w
∂
= OGCs + (y(n + h) − u(n) win wh wout )2
∂w
= OGCs + 2 (y(n + h) − ŷ(n + h)) u(n) win wout h wh−1
where y(n + h) is the teacher output, ŷ(n + h) the network output, and OGCs are
all the Other Gradient Components stemming from other input-output pathways
and which give zero contributions in temporal averaging.
122
wout
u(n+3) y(n+3) = u(n)
w
u(n+2) y(n+2)
w
u(n+1) y(n+1)
w
u(n) y(n)
win
The critical term in this gradient is wh−1 . You can smell the danger: if w < 1
and h is large, this term will exponentially shrink toward zero — the gradient
vanishes; and when w > 1, the gradient will exponentially explode with the depth
h.
In fully grown-up RNNs the analysis of vanishing/exploding gradients is not
so simple, but the basic mechanism is the same: since all time-slices of the un-
folded RNN are identical copies, gradient components arising from pathways that
span large temporal horizons are repeatedly either quenched or expanded at every
timestep. If the learning task includes the exploitation of long-term delayed input-
to-output effects (long memory), the vanishing / exploding gradient problem will
make it practically impossible to let the network find and encode these long-term
effects during learning.
You can find a more detailed presentation of the vanishing gradient problem
in the (by now, classical) deep learning textbook of Goodfellow et al. (2016),
Section 10.7. The conclusion drawn by the authors at the end of that section is
“... as we increase the span of the dependencies that need to be captured, gradient-
based optimization becomes increasingly difficult, with the probability of successful
training of a traditional RNN via SGD [= stochastic gradient descent] rapidly
reaching 0 for sequences of only length 10 or 20. [...] the problem of learning
long-term dependencies remains one of the main challenges in deep learning.”
123
This statement was made, however, six years ago — ages in the fast times
of deep learning. Long-term temporal dependencies (we called them long-range
memory effects) can now be trained into neural networks. The solution that is
currently preferred in the deep learning community is to abandon RNNs and use
feedforward (’transformer’) networks with a large input window size. However,
very recently, RNNs have returned on stage. New RNN architectures and special-
ized training schemes let them reach the same performance level as transformer
networks while being computationally much cheaper at use-time. This is advanced
stuff; if you are interested and daring, Orvieto et al. (2023) presents the latest re-
sults and tentative insights.
124
Cho et al. (2014)) are the dominating sort of RNNs used in deep learning.
The basic ideas behind LSTM networks are intuitive, but I am not aware of
a transparent mathematical analysis of why and how, exactly, LSTM networks
do function so well as they do. The explanations given in Section 10 in the deep
learning “bible” of Goodfellow et al. (2016) are somewhat handwaving. In what
follows I will try to explain the LSTM mechanism through a worked-out baby
LSTM which extends the embryo RNN example from Figure 50. The extended
version is shown in Figure 51.
r(n+3) wout
1 s(n+3) u(n+3) c(n+3) y(n+3) = u(n)
0
1
r(n+2) s(n+2) y(n+2) = 0
0 0 u(n+2) c(n+2)
1
r(n+1) s(n+1) y(n+1) = 0
0 0 u(n+1) c(n+1)
1
r(n) s(n) u(n) c(n) y(n) = 0
0 1
win
inputs
125
has little in common with biological systems. The memory cell together with
its surrounding control circuitry is called a memory block in the core historical
paper Gers and Schmidhuber (2001). In later literature by other authors, the
entire circuitry is often called an LSTM “cell”. I will stick to the terminology
of Gers and Schmidhuber and use the term “memory cell” only for the central,
memory-preserving neuron, and use “memory block” for the entire circuit.
I will give full equations for LSTM memory blocks later and for the time being
only present a simplified set of equations for our baby LSTM network, which is
made from a single (simplified) memory block.
The memory block has a memory cell with activation state c(n) ∈ R at time
n. The update equation for c(n) and the output equation are
c(n + 1) = 1 · c(n) + s(n) · (win u(n) − c(n)), (65)
y(n) = r(n) · wout c(n). (66)
The “store” signal s(n) acts as a multiplicative gate which only lets the input
u(n) affect the memory cell at the rare times when s(n) = 1. At such times, the
memory cell state is set to win u(n) at the next time step n + 1.
Similarly, the “read” signal r(n) allows the output unit to read out from the
memory cell an output value of wout c(n) only at times when r(n) = 1; at other
times the output reading is zero. We assume that the read signal jumps to 1 only
at the time points when the input h steps earlier was nonzero, that is at times
when the desired output is y(n) = u(n − h) 6= 0.
Furthermore, as long as there is no “store” input signal, that is, as long as
s(n) = 0, the memory cell retains it previous value c(n + 1) = 1 · c(n).
The only two trainable weights in this system are win and wout . If one uses again
the quadratic loss, one will find that the squared error at times when r(n) = 0 is
zero (because the network output and the teacher output are both zero, hence no
error), and at times where r(n) = 1, is
ε2 (n) = (y(n) − r(n) wout c(n))2
= (u(n − h) − wout win u(n − h))2
= u2 (n − h) (1 − wout win )2 .
At these times, the gradient of the error with respect to the two trainable weights
is
∂ 2
ε (n) = −2 u(n − h)2 (1 − wout win ) wout
∂win
∂
out
ε2 (n) = −2 u(n − h)2 (1 − wout win ) win
∂w
Obviously, gradient descent along this gradient will lead to a point where wout win =
1. The important insight here is that this gradient does not vanish or explode when
h gets larger. In fact, the gradient is independent of h.
126
The key to this independence of the error gradient from the memory depth h
is that, as long as s(n) = 0, the activation state c(n) is identically copied to the
next time step by virtue of the update rule c(n + 1) = 1 · c(n) + 0 · . . . in (65). The
weight w that we had in the embryo RNN (64), which made gradients vanish or
explode if not equal to 1, is fixed to be equal to 1 in (65).
LSTMs are, by and large, a sophisticated extension of this trick to carry mem-
ory information forward in time through a linear memory cell update c(n + 1) =
w c(n), where w is equal or close to 1.
Obviously, to make this work in a specific task with specific memory recall
demands, the save and read control signals s(n), r(n) must do a quite sophisticated
job. They acquire this functionality through learning — almost magically the
BPTT algorithm trains the activation dynamics for s(n) and r(n) such that they
best serve the learning task.
An LSTM network is an RNN which consists of several memory blocks and
possibly other, “normal” RNN units. The memory blocks and in particular their
save and read neurons receive trainable input from other blocks/units within the
network through trainable synaptic connections.
Here are the classical equations. LSTM memory blocks are made from 5 spe-
cialized neurons which have different update equations. The central neuron is the
memory cell with state c(n). The remaining four neurons are
an input neuron with state u(n). This corresponds to the input neuron u(n)
shown in our baby LSTM. The input neuron may receive external input (as
in the baby LSTM) and/or input from other sites within the network. The
input neuron is a “normal” unit with an update function
u(n + 1) = f (W u xu (n) + bu ),
an input gate neuron with state g input (n). In our baby LSTM, the role of the
input gate neuron was played by the external “store” input s(n). The input
gating neuron has the update equation
input input input
g input (n + 1) = σ(W g xg (n) + bg ),
input
where σ is always the logistic sigmoid and xg (n) is a vector of signals
composed of external inputs and/or outputs from other blocks or normal
neurons in the network;
an output gate neuron with state g output (n). In our baby LSTM, the role of
the output gate neuron was played by the external “read” input r(n). The
127
output gate neuron is similar to the input gate neuron and has the update
equation
output output output
g output (n + 1) = σ(W g xg (n) + bg ),
output
where again σ is always the logistic sigmoid and xg (n) is a vector of
signals composed of external inputs and/or outputs from other blocks or
normal neurons in the network; and
a forget gate neuron with state g forget (n). In our baby LSTM, the role of the
forget gate neuron was played by the constant 1 factor in (65). The forget
gate neuron is again similar to the input gate neuron and has the update
equation
forget forget forget
g forget (n + 1) = σ(W g xg (n) + bg ),
forget
where again σ is always the logistic sigmoid and xg (n) is is a vector of
signals composed of external inputs and/or outputs from other blocks or
normal neurons in the network.
The central element in a memory block is the memory cell c(n). Its update
equation, which corresponds to (65), is
The main extensions compared to the baby LSTM are that the memory cell
does not necessarily preserve its previous value with a factor of 1, but may “leak”
some of it by the multiplication c(n + 1) = g forget (n + 1) · c(n) . . . with the forget
neuron value.
The output y(n) of a memory block is given (similar to (66)) by
This output of an LSTM block is not necessarily external output, but may be
used internally as input to other memory blocks and “normal” units in the LSTM
network.
Notice that the three gate units must have the logistic sigmoid as their activa-
tion function. This ensures that the gating values g input (n), g output (n), g forget (n) all
range between 0 and 1. With a value of 1, the “gate” is fully opened, with a value
of 0 it is completely shut.
Figure 52 illustrates an LSTM block.
input output forget
LSTM networks are trained with BPTT. All the weights W u , W g , W g ,Wg
and biases in all memory blocks of an LSTM network are trainable. During train-
ing, the various gating cells (almost magically - but it’s just the effect of gradient
descent) learn to let input pass or not into the memory cell in a way that is
maximally helpful to reduce the training error.
Numerous variations of this memory block architecture are in use. For instance,
Gers and Schmidhuber (2001) describe memory blocks that may contain several
128
y
goutput
×
external input
and/or signals gforget
from other × c
locations in
the network
×
ginput
u
Figure 52: A standard LSTM block layout. For explanation see text.
memory cells, and they also admit the signals c(n) as inputs to the gating neurons
of the same block. Significantly simplified versions of LSTM blocks, called “gated
recurrent units”, are also explored and used. A discussion of these variations is
given in Goodfellow et al. (2016), Section 10.10. There appear to exist no clear
universally best winners among the assortment of proposed variations.
Today many high-level programming toolboxes for deep learning are available.
They all incorporate ready-made LSTM modules with automated BPTT adapta-
tions, such that the end-user does not have to care about programming the details.
One just plugs together LSTM units with normal neurons as one deems fit, hits the
“train” button, and hopes for the best. It often works. When it doesn’t — which
also happens to happen — it is good to understand how LSTM blocks function,
which now you do.
5 Hopfield networks
An obvious core functionality of biological brains is long-term memory (LTM).
You can remember your own name as well as details from your parent’s house as
well as impressions from the beach of your last summer vacation — and more. In
fact you carry with you a rich repository of memories from your life’s experiences
— you carry with you all the things that make you. The more you start thinking
seriously about all the things that you “know” (= have in your memory), the more
you will find that you know very very much indeed. How can your little brain store
all of that?
129
The question of human long-term memory has kept psychologists, cognitive
scientists, neuroscientists, psychiatrists, AI researchers and even physicists busy
since the beginnings of these sciences. It is one of the most heavily thought-about
scientific riddles, and it is far from being solved.
One thing is clear: neural memorizing is different from storing items in a shelf.
The car of your parents, which probably you can very well remember in much
detail, is not placed as a little model car somewhere in your brain. A surgeon
would not find it after opening your skull. Here are some of the research riddles
that surround long-term memory:
• How is LTM distinguished from other forms of memory, like short-term
memory, long-short-term memory, working memory?
• Is there a unique dedicated mechanism in the human brain for LTM, or is
it a complex system with many functional modules? Connected with this
question, what are commonalities / differences between LTMs in humans vs.
dogs vs. frogs vs. honeybees vs. little worms? They all have their peculiar
LTM capacities!
• To the extent that long-term memory “items” are “stable”, there must be
something physical in brains that is stably persisting through the tides of
time. For decades, the dogma in neuroscience (and machine learning too, for
that matter) had it that LTM memory traces are physically realized through
synaptic weights. Memorizing something for good means that some synaptic
weights are set and not changed thereafter. Well, this venerable dogma is
dissolving in these days. It can’t be that simple:
– Human memories change over time. A classical, super fascinating ex-
perimental psychology study (Bartlett, 1932) (25K citations, as of May
2022) reveals that over the decades of a human life, what one thinks
one had clearly “memorized” for good — continually changes, even dra-
matically. My wife and I are both aware of this and sometimes marvel
that, when she is super-convinced that the skirt her sister wore at our
wedding was blue and long, and I am dead sure that it was red and
short, we look at photos and find it was yellow and of medium length. If
you sometimes spend time recalling your childhood days, you might be
interested in a report (Nivison et al., 2021) of a long-time study where
parent-children interactions were recorded over decades, and where it
was found that, whether grown-up children remember their childhood
relations with their parents as positive or negative, this depends more
on their current mood and current relationship with their parents than
on the factual emotional experience twenty years before.
– On the microanatomical level, it is becoming clear that neural synapses
are incessantly changing, even becoming deleted or re-growing, at amaz-
ingly high rates (Castello-Waldow et al., 2019).
130
– A long-standing conundrum in artificial NN / machine learning research
is that if an MLP is first successfully trained on some task A, leading to
a specific formation of synaptic weights, and then subsequently trained
again on another task B, it will learn task B allright — but in the pro-
cess modify synaptic weights such that task A is no longer mastered.
This is called the problem of catastrophic forgetting. Attempts to over-
come it have only very recently become halfway successful, forming the
active field of continual learning within deep learning (review: Parisi
et al. (2019)). Again, what is being found out in this field is incompat-
ible with the dogma of that LTM memory traces are realized through
fixing stable synaptic weights.
• If your memory of a car isn’t like putting a little toy car on your brain
shelves, the question is how are memory items neurally encoded? If you
think of your grandmother, does this mean that a specific grandmother neu-
ron is activated? Like always in neuroscience, there is solid evidence in
support of this hypothesis, and equally solid evidence against it (https:
//[Link]/wiki/Grandmother_cell). An alternative view is
that memories are encoded in a distributed way: thinking of your grand-
mother is the effect of a complex neural activation pattern that involves
large portions of your brain.
• Ok., assuming it is understood how memory items are encoded, and how
they are written, how are they retrieved when you “recall” them? Recalling
leads to two subproblems: addressing and decoding:
– How to you mentally “point to” memory items? how do you “know”
what to retrieve? In a digital computer, addressing is done by pointers
to physical memory registers. But equivalents of C pointers are unlikely
to be implemented in biological brains. Instead, it seems more plausible
that brains use content addressing. This means that in order to access
what you have stored about your grandmother, you need to start with
some fragments of your grandmother memories — for instance, her
name, or you think of her house, or of a family gathering. This way
of thinking about mental addressing is known as associative memory
models. Building theories about how associative memories function has
131
been a mainstream activity in neural networks research for decades.
It has led to a mountain of models documented in a wide-spanning
literature.
– The decoding problem is the twin of the encoding problem.
• Finally, it is unclear what a “memory item” is. When you remember your
grandmother, you will be recalling different aspects of her in different situa-
tions. You never think of your grandma in exactly the same way. Apparently
it is quite a simplification to think of well-circumscribed memory “items”.
This immediately leads to the highly disputed problem of what concepts
are and how they are neurally represented. Again, a vast literature on this
topic exists in psychology, AI and philosophy, with no consensus in sight.
A most instructive, influential, and readable book is Lakoff (1987) (30K
Google Scholar cites), and it has a nice title too: “Women, fire and danger-
ous things” The title outlines the semantic contents of a single concept in
an Australian Aboriginal culture and language.
• In deep learning, the problem of continual learning, i.e. the problem of or-
ganizing the incremental growth of the representational repertoire of an ar-
tificial neural network, is not satisfactorily solved. The magnificent achieve-
ments of deep learning come in the form of networks specialized each on a
specific task. A face recognition network cannot be further trained to also
recognize cars, let alone to control a robot arm. Partial solutions are emerg-
ing in these days — if you are interested, read the intro chapter in the PhD
thesis of He (2023).
132
orders of magnitude of energy and time. So far, a generally useful way for
in-memory computing has not been found.
Among the many, many models of neural long-term memories, there is one
which stands out: Hopfield networks (Hopfield, 1982). The Hopfield model of
a long-term memory is simple, mathematically transparent, deeply analyzed and
almost completely understood. In machine learning it spun off an entire family of
energy-based neural network models, among them the Boltzmann machine which
in turn was instrumental in getting deep learning off the ground. In the cognitive
neurosciences it still serves as a foundational reference model for associative neural
memories.
133
Figure 53: Cueing a HN by corrupted patterns leads to a retrieval of the uncor-
rupted, stored pattern. Left: corruption by noise (leading to pattern
restauration upon retrieval). Right: corruption by omission (leading
to pattern completion functionality). Images taken from Hertz et al.
(1991).
134
4. Recall is started by presenting a (corrupted) input pattern u, which is set
to be the initial state of a state trajectory which evolves according to the
recurrent dynamics of the HN. The trajectory leads through a state sequence
which at every step reduces the energy and thereby necessarily ends in a local
minimum — a fundamental memory is retrieved.
E
energy E
u x2
x1 x3
Figure 54: Schematic of energy landscape over the pattern/state space of a HN.
The pattern space is here artificially rendered as a 2-dimensional Eu-
clidean space; this is only for an intuitive visualization as the pat-
tern space {−1, 1}L does not have an Euclidean topology. Funda-
mental memories ξ1 , ξ2 , ξ3 mark local minima of the landscape. Upon
input of a (corrupted) pattern u, the recurrent dynamics of the HN
leads through a state sequence which at every step reduces the en-
ergy, until the nearest local minimum is reached. Image retrieved from
[Link]/∼antonych/[Link] (no longer accessible).
135
u
u
u
u x2
x2
u’
u
x1 x3
Figure 55: Schematic of recall trajectories in a HN. The energy landscape is ren-
dered as a contour plot (blue = low energy). Three fundamental pat-
terns ξ1 , ξ2 , ξ3 here correspond to three pixel images of digits 6, 3, 2.
Since there are three local minima, that is, three point attractors, there
are three basins of attraction. They are marked by orange boundary
lines. Left: cues u lying in the basin of attraction of ξ2 are attracted
toward ξ2 . Right: a concrete 7-step trajectory leading from cue u to
ξ2 . Note that the clean gradient-descent (red) lines would look much
more jittery in the actual HN dynamics because it is stochastic. The
Hopfield dynamics is not a gradient descent dynamics – there is no
gradient in the discrete state space of HNs. Pixel images taken from
Haykin (1999).
136
diagonal. Here is a little example with L = 4:
0 −3 0.5 0.2
−3 0 1 0.1
W = 0.5
(69)
1 0 −2
0.2 0.1 −2 0
The most noteworthy thing here is that the weight matrix is symmetric. This
means that the connection weight of the synaptic link between neurons i and j
is the same in both directions, wij = wji . Connections are undirected in HNs —
which is biologically unrealistic but opens the mathematical doors to energy-based
computing.
A state of a HN is an L-dimensional binary vector with entries from {−1, 1}.
For example, (−1, 1, 1, −1)′ would be a state of our little HN example. There are
2L possible states.
A key idea about HNs and the Boltzmann machine (which we will meet later
in the course) is that each state is assigned to a real-valued quantity that is called
the energy of this state. In HNs, given a state x = (x1 , . . . , xL )′ ∈ {−1, 1}L , the
energy of this state is defined by
X 1
E(x) = − wij xi xj = − x′ W x. (70)
i,j=1,...,L
2
i<j
Note that, while the term x′ W x looks like a quadratic form (remember our
discussion in Section 2.2.3), it isn’t one, because the matrix W is usually not
positive semi-definite. If it would be, then the energy landscape would look like
in Figure 16, and it would have only one local (and global) minimum — useless
for storing many different patterns.
Neural network models in which an energy of states is defined lead to a tight
and mathematically rigorous connection between neural network dynamics and
statistical thermodynamics. Such models and the way of thinking behind all of
them leads to a subfield of neural network models called energy-based models. I
will not dig deeper into the links between HNs and theoretical physics at this point
but leave that for the next section in these lecture notes.
A side fact: Hopfield is a theoretical physicist. While his work certainly has
had a strong impact in theoretical neuroscience and machine learning, it also has
triggered a whole school of neural network research within theoretical physics. This
line of investigation is hardly perceived outside physics, especially not in neural
networks research — one of these strange cases of non-communication between
different communities in the sciences.
The rule for the state update dynamics is stochastic. If at time n the state is
x(n) = (x1 (n), . . . , xL (n))′ , the next state x(n + 1) is obtained as follows:
137
2. Compute its activation value at time n + 1 by
X
xi (n + 1) = sign( wij xj (n)). (71)
j̸=i
where “sign” is the signum function sign :PR → {−1, 1}, sign(z) = 1 if z >
0 else sign(z) = −1. In the rare case that j̸=i wij xj (n) = 0, set xi (n+1) =
xi (n).
3. Set x(n + 1) = (x1 (n), . . . , xi−1 (n), xi (n + 1), xi+1 (n), . . . , xL (n))′ , that is,
update only the activation of the neuron xi .
L=1
L=2
L=3
L=4
138
"Rolling downhill" in a hypercube
Figure
Just 57you
to give attempts to visualize
an impression... a 6-dimensional
Pattern hypercube,
states in a N = with a few state
6 HN (slightly
energies and pattern images included.
misleading though, images have N = 120!)
Indicates
energy at state
(hypercube
corner)
Single-neuron
updates
(flipping one
component of
state vector)
L=6
Figure 57: The hypercube state spaces for L = 6. Some (non-negative) energies of
states are schematically indicated. In any pair of neighboring states,
one of the two states has a lower energy than the other. The HN
update rule always moves to a neighboring state with lower energy, or
repeats the current state. Three “de-noising” sequences starting from
noisy versions of the digit “3” are indicated for intuition, but note that
these digit images have 10 × 12 pixels and thus should be correctly
shown in an L = 120 dimensional hypercube. Furthermore the graphic
is incorrect in that per update step (one move along one edge in the
graph) should change only a single pixel.
5.4 Training a HN
The learning problem for a HN is this:
Given N L-dimensional training patterns ξ1 , . . . , ξN , find a weight matrix W
which creates an energy landscape that has the training patterns located at the local
minima, and every local minimum corresponds to one of the training patterns.
Achieving this goal is not always possible as we will see, but the conditions
when it is possible are well understood.
There are two methods of finding a weight matrix W satisfying the learning
objective (if a solution exists). The first method is very fast and simple: there is an
analytical formula that directly computes W from the training patterns ξ1 , . . . , ξN .
The second method is iterative-incremental and may appear unnecessarily time-
consuming, but it is biologically plausible (well... as far as it goes...) and could
139
be used by real brains, whose neurons cannot compute the analytical formula of
the first method (well...). I will present the formalism of both methods and after
that explain in intuive terms why they work.
Initialization: create a random initial weight matrix W (0) (symmetric with zeros
on diagonal).
140
Loop: • At update step k, present one of the training patterns to the net-
work (picked at random or in cyclic order), say ξ = (ξ (1) , . . . , ξ (L) )′ is
presented. The weight matrix before entering step k is W (k − 1).
• Update all weights wij (k − 1), where i 6= j, by
Stop when you reach a condition of a previously defined stopping criterion. For
instance, you can stop when the largest weight (in absolute value) hits a
predefined ceiling; or when a test set of corrupted input patterns is recalled
correctly; or when the energies of the training patterns are smaller than
the energies of all patterns that are similar to training patterns except at
one flipped vector entry; or (best) when the pairwise ratios of weights in the
sequence W (k) appear to converge up to a predefined, small residual change.
Notice that the incremental update (73) can be written in matrix form as
W (k) = W (k − 1) + λ (ξ ξ ′ − I) . (74)
From this finding it is easy to conclude that this incremental learning rule will
converge to a weight matrix that is the same as the one obtained from the ana-
lytical solution (72), up to a scaling factor which grows larger and a residual error
in matrix entry ratios which decreases the more the longer you run the iterative
computation. The local formulation given in (73) is meant to underline biological
plausibility: in order to effect the change of a weight wij , only information that is
locally available at this synaptic link is needed.
A general note: in any neural network architecture any learning rule that
needs non-local information for computing its update for a weight wij (that is,
information about neurons other than neurons i and j), is considered biologically
impossible. Biological synapses can be adapted only on the basis of information
that is locally available at that very synapse. Specifically, the backpropagation
algorithm is non-local.
141
The learning rule (73) also instantiates a general learning principle that is
believed to be ubiquitously effective in biological brains. This is called Hebbian
learning. Because Hebbian learning is mentioned in the NN literature in many
places, I will expand a little on this topic.
Donald O. Hebb (1904-1985) was a psychologist / neuroscientist (he started
out as a teacher of English) who gave us one of the main guiding principles for
understanding how neural circuits can represent conceptual information. In his
book The Organization of Behavior (Hebb, 1949), he developed a theory of bio-
logical neural learning mechanisms that could explain how the human brain can
memorize, recall and re-generate perceptual and conceptual patterns to which it is
repeatedly exposed at learning time. To explain this capacity, he developed a the-
ory of cell assemblies. A cell assembly can be thought of as a group of neurons that
are mutually exciting each other through positive (“excitatory”) synaptic connec-
tions. These mutually exciting connections develop in a learning process. If some
perception (for instance, a child seeing its mother’s face) is repeatedly made, and
at each presentation a certain subset of neurons in, say, a visual processing area
of the brain is simultaneously activated by this perception, then these repeatedly
co-activated neurons will form mutually excitatory links. In Hebb’s wording, a
memory trace is formed: the perceptual experience becomes encoded in the cell
assembly. This assembly can then function very much like a trained HN: if some of
its neurons are excited by sensory input or input from other brain areas, the entire
assembly tends to self-excite. Hebb stated this learning principle in a paragraph
that has become one of the most often cited sentences in neuroscience:
“When an axon of cell A is near enough to excite a cell B and repeatedly or
persistently takes part in firing it, some growth process or metabolic change takes
place in one or both cells such that A’s efficiency, as one of the cells firing B, is
increased.”
In the folklore of neuroscience, this principle has been shortened to the catch-
phrase
“Cells that fire together, wire together.”
Hebb was, in the first place, a psychologist, and did not use mathematical
formalism to make this learning mechanism more precise. He wanted to model
the biological brain and tried to relate his principle to what was known at his
time about neural microanatomy. He also tried to relate his principle to a large
spectrum of findings in cognitive and developmental psychology that were known
at his time. He certainly would never have thought of simplifying a neuron to a
mathematical variable that can only take the values −1 and +1, and also the idea
to build artificial neural systems was alien to him — that idea would start its rise
but 10 years later, in the form of the Perceptron.
The Hopfield network is nonetheless an obvious instantiation of “Hebbian
learning”. I conclude this excursion by emphasizing that today there is a large
spectrum of rigorously formalized neural learning mechanisms which are Hebbian
in their core. They give a spectrum of different answers to problems which arise
142
when one starts working out the consequences of the cell assembly idea. I mention
just two:
• The mutual excitatory connections of cell assemblies must in some way be
complemented and compensated by inhibitory connections — otherwise the
entire brain would burst into flames by global spread of self-excitation.
• Neurons can only excite one another by sending spikes. A spike lasts only
about one millisecond. But sensory impressions or the activation of con-
ceptual items in one’s mind may last seconds — three orders of magnitude
longer. This mismatch in timescales makes it necessary to develop inge-
nious, nontrivial schemes to account for fine-grained temporal relationships
between the spike events in neurons within an assembly.
The first problem is solved in the HN model by giving equal citizenship rights
to positive and negative neural activations, and positive and negative synaptic
weights. The second problem does not arise in HNs because there are no spikes in
the simple neuron model used in HNs.
5.5 Limitations
Above I inserted the cautionary clause “... if possible ...” at various places. In
fact, it is not always possible to store N training patterns in an L-neuron HN,
such that after learning there are exactly N local minima in the energy landscape
which correspond to the precise training patterns. Several things can go wrong:
• The fundamental memories which are shaped by the learning process may
turn out to be only approximately equal to the original training patterns,
that is, they will differ from the training patterns in a few bits.
• Not all training patterns can be stored — the storage capacity of a HN
is limited. Specifically, consider the ratio N/L of the number of training
patterns over the network size. This fraction is called the load of training a
HN. The following facts are known for HNs:
1. For a load N/L > 0.138 (I give only first three significant digits, pre-
cise number is known), storing training patterns using (72) breaks down
entirely: none of the training patterns will be stored, not even approx-
imately.
2. A HN works really well only for loads 0.03 < N/L < 0.05, in the sense
that the training patterns end up in local minima that correspond to
perfect copies or only slightly altered versions of the original patterns.
• Besides at the locations of the training patterns, other local minima are cre-
ated in the energy landscape which do not correspond to training patterns.
Such “false memories” are called spurious states in the HN literature.
143
There is one thing that always “goes wrong”: if ξ is a fundamental memory,
then also the sign-inverted pattern −ξ is a fundamental pattern (why? find out
for yourself, it is an easy one).
All of this is not so good news, inasmuch as these limitations more or less
render HNs useless for practical applications.
On the other hand, the reasons why these limitations occur, and which of
them strike how badly at what load levels, are almost completely known. The
corresponding mathematical analyses have been carried out, mostly by theoretical
physicists, in the 1980’s and now form a classical body of rigorous insight into the
conditions of storing information in a neural network.
In the remainder of this subsection I document these findings in some more
detail. This is optional reading which will not be queried in exams.
144
a 0
Figure 58: The error function Φ(a) gives the area under the pdf of the standard
normal distribution on the left-hand side of a. It is the cumulative
density function of the standard normal distribution.
training pattern is entirely unstable and is not located near a local minimum of
the energy landscape. In the worst case, all bits in ξ become randomly flipped,
resulting in a maximal Pflip = 0.5.
How bad this avalanching becomes depends on the load. In the limit of large
L, at a load of N/L ≈ 0.138, every training pattern becomes maximally unstable.
Avalanches start to occur for loads n/L > 0.05. Figure 59 shows the growth of
the avalance effect with the load.
Figure 59: Percent of changed pattern bits under iterated network update, vs. load
(here denoted by �). Figure and result from Amit et al. (1985).
145
right panel of Figure 55 were copied from the Haykin book). This makes L = 120.
A small set of training patterns consisting of merely N = 8 images was used. It
contained pixelized versions of the digits 0,1,2,3,4,6,9, plus a pattern showing a 5
× 6 sized black square in the upper left corner (Figure 60).
Figure 60: The 8 training patterns from a HN computer demo presented in Haykin
(1999). The slight jitter visible at some of the black-white boundaries
are artefacts from my postprocessing of photocopies from the Haykin
book.
When testing the trained HN with about 43,000 different cues obtained by
randomly flipping one fourth of the pixels in one of the training patterns, the
network frequently ended in a local energy minimum (attractor) state which was
not one of the 8 training patterns. Such “wrong memories” are called spurious
states in the HN terminology. Figure 61 shows 108 such spurious states. Most
likely there are many more — it is not possible to exhaustively search a pattern
space of 2120 many different patterns! Haykin cites without further explanation
findings of Amit (1989) to the effect that there are three kinds of spurious states
(see Figure 61), one of which is always present: namely, the sign-inverted training
patterns always become local energy minima if the original training patterns do
so.
• For all loads N/L: stable spin glass states exist — these are states that look
“just random” and which are uncorrelated to all training patterns.
• For N/L > 0.138: spin glass states are the only stable ones.
• For 0 < N/L < 0.138: stable states close to desired fundamental patterns
exist.
146
[Link] states (occur
necessarily)
Figure 61: Some spurious states found in the demo in the Haykin book.
• For 0 < N/L < 0.05: pattern-related stable states have lower energy than
spin glass states.
• For 0.05 < N/L < 0.138: spin glass states dominate (some of them have
lower energy than pattern-related states).
• For 0 < N/L < 0.03: additional mixture states exist, with energies not quite
as low as the pattern-related states.
In summary, a HN works really well only for a load 0.03 < N/L < 0.05.
147
A dramatically different update dynamics would result, with no guarantee
to lead to local energy minima.
Such a decisive difference between local stochastic and global deterministic
state update rules is not confined to HNs, but is a common finding in all
kinds of dynamical systems that have states made of discrete bits.
Specifically, digital computer chips use a global deterministic update rule.
This necessitates a global clock which on your notebook computer runs at a
Gigahertz rate. In contrast, in the emerging field of neuromorphic microchips
(I will say more about them in the final lecture), no global clock is available
and the update dynamics is necessarily local and stochastic, driven by local
physical laws. This is one of the reasons why the computational mathemat-
ics in such unconventional microchips is entirely different, and much less
understood, than the mathematical theory for classical digital computing.
Variations. The HN model that I described in this section is but one specific
version among many other formalizations. For instance, one could use state
values {0, 1} instead of {−1, 1}. This is a minor variation which leads to
an entirely equivalent theory. A much stronger modification is to admit
continuous values of the neurons, say from the continuous interval [0, 1].
The quality of findings is similar to what we saw in this section, though
details and formalism differ. The Haykin book discusses both discrete and
continuous HNs.
148
would be a circular sequence ξ (1) , . . . , ξ (N ) = ξ (1) . The learning formula (74)
would be replaced by
′
W (k) = W (k − 1) + λ ξ (i+1) ξ (i) , (75)
and (if convergence is achieved) the weight matrix would become non-symmetrical.
A substantial literature on such heteroassociative networks exists (of course).
HN pro’s. Many things are good and influential and insightful about HNs:
149
6 Moving toward Boltzmann machines
One of the central models of a “cognitive” neural network is the Boltzmann ma-
chine (BM). The Boltzmann machine is a must-know for machine learners, AI
afficionados, cognitive scientist and neuroscientists alike because
• it helped paving the way for what we call deep learning today (and for
some years before 2010, Boltzmann machines were the main driver of deep
learning),
This sounds too good to be true - there must be a drawback. Indeed there
is. Training and using a Boltzmann machine is computationally very expensive.
This has barred a widespread use in machine learning (except for the early years
of deep learning, say between 2006 and 2012, when Boltzmann machines were
competing with MLPs and CNNs for realizing deeply layered neural learning sys-
tems). Although in those exciting years the youthfully emerging DL field was using
a computationally streamlined version, the restricted Boltzmann machine (RBM),
ultimately the computational costs (among other reasons) led to a depreciation of
150
BM/RBMs in machine learning — though there is still ongoing research on them in
the DL world. However, biological brains may be able to do the core computations
which are so expensive in the BM (namely, sampling operations) at no extra cost,
just by exploiting neural noise. Therefore, the BM and related models retain their
fascination for theoretical neuroscience, the theory of cognitive computing, and
maybe for future non-digital neural microchips. Indeed, the PhD project of Fabian
Ijpelaar in my research group [Link] is addressing ex-
actly this challenge: novel, Boltzmann-like, stochastic ’computational mechanics’
for use with new nanoscale materials that are developed in a partner project in
the materials science group of Beatriz Noheda [Link]
zernike/nanostructures-of-functional-oxides/noheda-group/noheda-group.
In order to understand BMs, one must know two formal concepts from statis-
tics and statistical physics, namely the Boltzmann distribution (which gave the
Boltzmann machine its name) and the notion of a sampling algorithm. These con-
cepts are also forming the basis for the general, large class of energy-based models.
It is thus a well-invested effort to make friends with these two concepts. The class
of energy-based models include, for instance,
151
of deep learning. This line of modeling is therefore not dimmed by the deep
learning revolution, a fate that so many other fields in machine learning have
suffered.
Before we can start enjoying the Boltzmann machine, we must equip ourselves
with a fair understanding of the Boltzmann distribution and the idea of sampling
algorithms. This is what this section will give you.
152
strives to achieve: mathematically deduce the macroscopic laws of thermodynam-
ics from assumed mechanisms of how microscopic particles interact in large num-
bers, such that the macroscopic observables and laws can be derived as properties
of statistical distributions.
Brains are also macroscopic lumps of matter in which the interactions of large
numbers of microscopic “particles” (the neurons) give rise to macroscopic observ-
ables (for instance, the words coming out of your mouth). It is a naturally inviting
idea to describe brains with the tools of statistical physics. The Boltzmann ma-
chine is one way of doing exactly this.
A (rather, the) founding father of statistical thermodynamics was Ludwig
Boltzmann (1844-1906, [Link]
an Austrian professor of physics and philosophy (!).
In order to start getting familiar with the Boltzmann distribution, let us con-
sider the textbook example of a vessel filled with water.
The macroscopic description of this system is simple: volume, temperature and
pressure are enough to characterize it. These can be measured with macroscopic
instruments.
The microscopic description is based on the notion of a microstate. A mi-
crostate is (in a first approximation, treating water molecules as elastically re-
bounding balls) a specification of all the 3D position coordinates and 3D velocity
components of all individual H2 O “balls” in that vessel. This is a very high-
dimensional real-valued vector, say of dimension d, which completely character-
izes the molecular-level state of affairs inside the vessel at given moment in time.
Teaser question: what is d if we are modeling one liter of water? Wolfram Alpha
says: one liter of water contains about 3.33 × 1025 molecules, each of which has 3
position and 3 velocity components, thus d ≈ 2 × 1026 .
Boltzmann asked, and answered, the following fundamental question: what is
the probability distribution of these microstates? — and from that distribution he
inferred the laws and values of macroscopic observables.
First let us understand that this distribution is not uniform. Some kinds of
microstates are more probable to occur than others. The water molecules bounce
against each other and the vessel’s walls, exchanging impulses, in a wild stochastic
dance. Boltzmann assumed that the vessel is submerged in a heat bath, which you
can visualize as an infinite ocean of water that has the same temperature as the
water in the vessel and the vessel’s walls. Ocean molecules bounce against the
vessel walls, vessel wall molecules transmit these impulses to water molecules inside
the vessel, and conversely there are stochastic transmissions of molecular impulses
from the inside to the outside. Just by stochastic coincidences, at some moments
large amounts of impulse energy will have found their way from the outside to
the inside or vice versa; and at other moments, less of it. A microstate’s (kinetic)
energy is the sum of all the kinetic energies of the individual water molecules in
the vessel. We don’t have to understand how the kinetic energy of a “ball” with
a given mass and velocity is defined (you remember from high-school physics?).
153
The bottom line is that the energy of microstates is wildly fluctuating all the time
due to the random exchange of energy with the heat bath.
Side remark: this energy fluctuation is due to the embedding of the vessel in
a heath bath. If the vessel would be perfectly isolated, then the energy of all
possible microstates would be constant due to the law of conservation of energy.
The statistical physics of open versus isolated systems are eminently different!
Let us denote microstates by s and the energy of a microstate by E(s). Now
we have put our foot on the doorstep of one of the grandest and most far-reaching
principles of modern physics. Boltzmann reasoned that the probability of a mi-
crostate depends only on the microstate’s energy E(s) and the temperature T of
the vessel/heatbath system. Note that physicists measure absolute temperature
(in Kelvin), where zero is the lowest possible temperature. The probabilities of
microstates are described by a simple formula which today is called the Boltz-
mann distribution. Since the position&velocity vectors describing microstates are
continuous-valued vectors, this distribution is written down as a pdf which assigns
a probability density value p(s) to each microstate:
1 −E(s)
p(s) = exp( ). (76)
Z T
In this equation, Z is the normalizing factor that ensures that the pdf p inte-
grates to unity: Z
−E(s)
Z= exp( ) ds, (77)
s∈S T
where S is the space of all possible microstates. I am not a physicist but I would
think that S = Rd≥0 comprises all non-negative real-valued vectors of dimension
d, where d = 6N and N is the number of particles in the vessel and the factor 6
comes from the three location and velocity coordinates of each particle.
Note that both p and Z depend on the temperature T , so we sometimes write
p(s, T ) and Z(T ). Z is called the partition function in statistical physics (“func-
tion” because it depends on the temperature). The partition function plays a
central role in physics. In most cases it cannot be calculated analytically. One
needs supercomputing power to even only approximately estimate it, and the pos-
sibility to actually do this has changed the face of modern physics (computing
Z is one of the reasons why physicists need supercomputing facilities so badly).
Computing estimates of Z is also of importance in certain applications of machine
learning, deep learning in particular. The entire Chapter 18 of the deep learning
“bible” (Goodfellow et al., 2016) is devoted to estimation algorithms for Z. Luck-
ily, in many applications, among them the Boltzmann machine, Z cancels out and
need not be computed.
The Boltzmann distribution can also be defined on spaces of discrete mi-
crostates (finitely or countably many). Then the pdf from (76) turns into a prob-
154
ability mass function (pmf) and the integral in (77) into a sum:
1 −E(s)
P (s) = exp( ), (78)
Z T
X −E(s)
Z = exp( ). (79)
s∈S
T
The shape of the pdf (or pmf) changes quite dramatically when the temper-
ature parameter is varied. Figure 62 shows an example. There are two notewor-
thy extreme cases. When the temperature is very high, the terms exp( −E(s) T
) all
approach 1, which (after normalization by division with Z) gives a uniform dis-
tribution. At high temperatures, all microstates become almost equally probable.
When the system is cooled down toward zero, the probabilities of the low-energy
states grows relative to the probability of the high-energy states. Ultimately,
the distribution converges to a point distribution where the (single) lowest-energy
state has a probability of 1. Cooling a physical system down toward zero will
localize its distribution at the lowest energy state! It should be noted however
that this cooling down must be done veeery slowly, in a process called adiabatic
cooling in order to see this effect in real physics experiments.
The fact that a (sufficiently slow) cooling of a Boltzmann system leads to the
global minimum of the energy landscape is the basis of the simulated annealing
algorithm. This is a general-purpose optimization algorithm which can find, in
principle, the global (!!!) minimum of any cost surface. In this respect it is vastly
more powerful than gradient-descent optimization algorithms which we discussed
in Section 2.2.2. The simulated annealing algorithm is really good stuff to know
and I will explain it in more detail in Section 6.4.
Two facts about the Boltzmann distribution worth knowing:
155
Energy
Boltzmann pmf
Figure 62: A Boltzmann distribution over a finite space of 100 microstates. The
distribution is given by its pmf (red). The underlying energy function,
which is not changing with temperature, is rendered by a blue broken
line.
While the Boltzmann distribution has been found and explored in physics, its
ideas and maths transfer to any other application domain where there are macro-
scopic systems which can switch between large numbers of microstates, each of
which has an “energy”. In particular, neural networks (artificial ones and real
brains) can be regarded as macroscopic systems, with vectors of activations of all
neurons being the microstates. If one defines some sort of “energy” for the activa-
tion vectors, the Boltzmann machinery can be launched. The “energy” which one
defines can be any function from microstates to the reals; physical considerations
can be ignored.
156
6.2.1 What is it?
A sampling algorithm, or for short a sampler, is a computational procedure which
generates “random” examples from a given distribution. Samplers there are many.
Some are highly specialized and can only generate random examples from a single
distribution, others are generic and can “sample from” any distribution which is
given by a pdf or pmf or energy function.
You all know a sampling algorithm that samples from the uniform distribution
over the interval [0, 1]. It comes as a ready-made function with all programming
languages that I know, including MS Word. It is typically named rand. Every
time you evaluate this function, it generates a new “random” number from the
interval [0, 1].
In fact, the algorithms that sit behind rand are deterministic, and the outputs
from rand are only pseudorandom numbers. Many algorithms are known which
deterministically generate numbers between [0, 1] in a way that can hardly be dis-
tinguished from true randomness. It is a quite sophisticated corner of maths where
such “pseudorandom generators” are cleverly thought out. We just use them as
if they generated truly random outputs and don’t think twice about the math
miracles behind the curtain which turns digital determinism into (almost) ran-
domness. If you want to get true randomness, for instance for unbreakable codes,
you’d need to build a physical, unconventional microchip that gets its randomness
from quantum fluctuations.
Formalizing the concept of a “sampler” in mathematical rigor requires some
stochastic processes theory and is beyond the scope of our course. But I think the
idea is intuitively clear. Consider a pdf as a “landscape” with hills and valleys (as
the two-hill landscape mapped in Figure 63). A sampler is a mechanism which,
each time it is executed, lets fall down a grain of sand on a plane that was flat
and empty at the beginning. In the long run, these grains of sand should pile up
to a landscape whose profile is the same as the pdf landscape.
Note that a sampler need not “jump around” wildly. Modern sampling algo-
rithms indeed have some kind of memory. They linger for a while in the vicinity
of the place where the last grain of sand was dropped: the next grain will typ-
ically fall down somewhat close to the previous one. Such samplers perform a
random walk over the sampling space, dropping a grain at each step, and the av-
erage forward motion that walk can be slow. The entire landscape is covered
only in the long run. The sampler that we will be using for the Boltzmann
machine is of this kind. The general theory behind the design of such sam-
plers is called Markov Chain Monte Carlo (MCMC) sampling. Once powerful
computing hardware became available, MCMC changed the face first of physics,
then of other sciences, because these methods (and only these methods) made it
possible to simulate complex stochastic systems in physics, chemistry, biology,
economics and the social sciences. We cannot dig deeper here. If you are in-
terested: a classical tutorial is Neal (1993), and in my legacy lecture notes on
“Algorithmic and statistical modeling” I give a more in-depth outline (Section 4.5
157
Figure 63: Sampling from a 2-dimensional pdf. The pdf is here rendered by a con-
tour line plot; it consists of two “hills” separated from each other by a
“valley”. Each blue dot is one output of a sampling algorithm. Graphic
taken from the online course notes of IEOR E4703: Monte-Carlo
Simulation, Columbia University [Link]
io/files/MonteCarlo/MCMC_Bayes.pdf where they were taken from
the book Bayesian Reasoning and Machine Learning by D. Barber.
in [Link]
158
to get an intuitive visual idea of the shape of a complex distribution. It is, as per
today, likewise impossible for mathematical formalism to characterize the overall
“geometry” of most real-life pdfs or pmfs.
The only thing that usually can be computed is just the individual pdf or pmf
value of a given point in the state space of the distribution. That is, for any point
s in a high-dimensional state space S, it is possible to compute the corresponding
pdf value p(s) or the pmf value P (s).
Oh... no. Not even that is usually possible. There are two obstacles. The
first obstacle is numerical underflow. For a demonstration, consider a probability
distribution over the set S of possible binary health states of the world popu-
lation. Assume for simplicity that 10 billion humans live on our planet. Each
of them can be ill or healthy — 1 or 0. A global health status vector is thus a
binary vector of dimension 10e10. This makes for a large but finite state space
S = {0, 1}10,000,000,000 . Since S is finite, a probability distribution over this space
has to be described by a pmf. The probability P (s) of a state s is, on average,
2−10,000,000,000 . Written to base 2, which is what digital computers do internally,
this would be written as 0.0 . . . 01, with 10, 000, 000, 000 digits after the dot. But
the machine precision on your 64-bit computer allows only for 53 bits of precision.
Any number smaller than 2−53 is treated as zero by your computer.
A standard escape from such numerical underflow issues is to always work not
with the raw probabilities P (s) but with their log values log(P (s)). In base 2,
the average log probability then is log2 (P (s)) = −10, 000, 000, 000, an order of
magnitude (rather, minitude) that is within convenient precision reach of your
computer. If you read math-oriented papers in machine learning you will see log
probabilities all over the place.
The second obstacle is the partition function. It is unknown and essentially
uncomputable in most cases of interest. What one normally only does have is
knowlege about probability ratios, like “the probability of microstate si is three
times higher than that of microstate sj ”. Knowing the absolute, not only the rel-
ative, probability values (pmf or pdf) amounts to knowing the partition function.
Luckily, the powerful MCMC sampling algorithms only need relative probabilities,
that is, probability ratios — or, if one works with log probabilities, log probability
differences.
Ok., let us return to our main thread. I said above that the only thing that can
be computed for a distribution given by a pdf or pmf are the values p(s) or P (s),
or rather their logs. Let us return to the global health status vector example. An
epidemiologist might want to know, what is the probability Q that more than one
tenth of the world population is ill? This turns out to be a computational show
stopper. Let SQ ⊂ S be the set of all health state vectors that have more than 10
percent 1’s entries. The mathematical formula that defines Q is
X
Q= P (s).
s∈SQ
This is not computable for two reasons: numerical underflow, and the gigantic
159
size of SQ .
But now assume you had a way to get a “fair sample” of manageable size,
drawn from the distribution P . Say, you have “drawn” 1,000 examples s1 , . . . , s1000 .
You count how many of these example vectors have more than 10 percent 1’s in
them. Your count number is q. Then, thinking about it, you see that q/1000 is an
estimate of Q! This estimate will become more precise if you collect larger samples,
converging to the correct probability Q as your sample size goes to infinity.
This example demonstrates one of two major things that sampling is good for:
sampling can be used to get estimates of probabilities for events that a researcher
is interested in. Similarly, sampling can be used to get estimates of other sta-
tistical quantities, like expectations, variances, partition functions, or all sorts of
integrals over functions over the sample space. And there is no other way than
sampling to get access to these quantities of interest. I think you can divine what
a game-changer in the sciences it was to afford of both general-purpose sampling
algorithms and computing hardware with enough bandwidth.
There is an angle to this success story which is worth knowing and think-
ing about. Both the mathematical development of general-purpose sampling al-
gorithms, and the computational exploit of them on the first powerful enough
digital computers, were done in the Los Alamos Labs in the context of devel-
oping the hydrogen bomb. See [Link]
computing-and-manhattan-project (authors, 2014) for a historical outline and
Metropolis et al. (1953) for the landmark scientific publication (48K Google Scholar
cites as of May 2022).
The other major thing that sampling is good for: it’s just the sampling itself.
It can be used to generate examples from a statistical model of some interesting
part of the world. When the Boltzmann machine is used after it has been trained,
this is the way how it is used. In cognitive terms: the random walk of an artificial
brain state sampler creates a “stream of thought” which in the long visits all the
places and themes that the brain knows about. The technical term for this process,
used both in psychology and machine learning, is confabulation.
160
or a computable pmf value P (s) (for discrete spaces). In fact, even less need be
given: it is enough to have the pdf or pmf only up to some unknown normalization
factor. The Metropolis sampler only needs the ratios p(s)/p(s′ ) or P (s)/P (s′ ) to
run, and these ratios remain the same if the pdf or pmf is scaled by some constant
factor. This also makes it unnecessary to compute partition functions, because
they cancel in these ratios.
I will present the Metropolis sampler for the case of a given pmf, because that
is the situation we will meet in the Boltzmann machine. The pdf case is entirely
analogous and you can easily translate the pmf recipes into pdf recipes.
So, here is the scenario. We are given
P a finite state space S and a non-negative
≥0
function F : S → R whose sum s∈S F (s) is finite. I call this function a
proto-pmf
P because it could be turned into a pmf P by scaling it with 1/Z =
1/ s∈S F (s), but this normalization is not needed (and often not feasible), so we
stick with the proto-pmf F .
The task is to generate a potentially endless sequence s1 , s2 , . . . of states from
S such that, in the long run, this sequence of sampling points would re-model the
pmf landscape in the intuitive sense of our “dropping grains of sand” metaphor.
I do not give the mathematical derivation of why the Metropolis sampler does
its job, but just describe the algorithmic procedure. If you are interested in the
mathematical derivation, you can find it in Section 4.7 in my legacy lecture notes
on Algorithmical and Statistical modeling, online at [Link]
minds/uploads/LN_AlgMod.pdf, or the MCMC tutorial of Neal (1993).
The Metropolis algorithm generates the sample point sequence s1 , s2 , . . . by
implementing a Markov process mechanism, that is, every newly generated point
sn depends (only) on the predecessor sn−1 (note: a Markov process generalizes
discrete, finite-state Markov chains to continuous state spaces). In order to get
the whole sequence started, you have to “guess” the initial point s1 . You can pick
it arbitrarily.
Now assume sn has already been computed. In order to compute sn+1 , the
Metropolis algorithm executes two steps, each of which contains a cheap-to-compute
random decision:
Step 1: randomly propose a candidate s∗ for sn+1 . For this step one needs
to design and implement a mechanism to sample from a proposal distribution.
Each instance of a Metropolis algorithm comes with such a mechanism.
Mathematically, a proposal distribution is a conditional distribution over S
which gives the probability to choose s∗ given sn . Let us write Pprop (s∗ | sn )
for the probability to pick s∗ as candidate for sn+1 .
The proposal distribution should be chosen such that one can sample from it
cheaply. For instance, it could be the uniform distribution on a small hyper-
cube centered on sn , or a multidimensional Gaussian distribution centered
on sn . Efficient samplers exist for such elementary distributions.
161
Step 2: randomly accept or reject the candidate s∗ . If s∗ is accepted, it be-
comes the next output of the sampler, that is sn+1 = s∗ . If s∗ is rejected, it
is discarded and the sampler repeats the previous value in its next output,
that is sn+1 = sn .
This requires a subroutine for deciding whether to accept or reject the can-
didate. This decision is again random. It works by first computing an accep-
tance probability, then accept the candidate with this probability (and re-
ject it if it is not accepted). Several such decision-making procedures, called
acceptance functions, are known which result in sampling point sequences
s1 , s2 , . . . which asymptotically re-model the “landscape” of the proto-pmf
F . Here are the two most common and famous ones:
F (s∗ )
Paccept (s∗ | sn ) = . (80)
F (s∗ ) + F (sn )
that is, the Metropolis acceptance function accepts the proposed candi-
date with certainty if its proto-pmf is larger than the one of the previous
value, and it accepts it with probability F (s∗ )/F (sn ) if the F-value of
the candidate is lower than the one of the previous sample point.
Both acceptance functions can be computed solely on the basis of the ratio
F (s∗ )
r = F (sn ) . The Boltzmann acceptance function can be re-written as
r
Paccept (s∗ | sn ) = ,
r+1
and the Metropolis acceptance function as
∗ 1, if r ≥ 1
Paccept (s | s ) =
n
r if r < 1.
162
I conclude this subsection with a few remarks to round off the picture.
• In order to guarantee that the Metropolis algorithm, as outlined above,
indeed yields a sampling point sequence that correctly re-shapes the F -
landscape in the long run, some additional conditions must be satisfied.
– A necessary condition is that the “random walk” process which yields
the sampling point sequences is ergodic. This is an involved concept
from information theory which is beyond the scope of this course. In
intuitive terms this means that every point in S can be reached by some
sample point sequence at some time, regardless of where the sequence
was started. This condition would be violated, for instance, if the F
landscape is zero everywhere except at two “hills” H1 , H2 which are
separated from each other by a zero-F “lowland” of width w, and the
proposal distribution is chosen such that it always proposes a candidate
s∗ which is closer to sn then w. If the sampling sequence is started in
the hill H1 , it can never cross the flatland toward H2 . Checking whether
a Metropolis sampler that one has designed has the ergodicity property
can be nontrivial (sometimes it is easy) and there is no general recipe
to assure this property.
– Not every proposal distribution will yield a valid Metropolis sampler.
Assuming ergodicity, a sufficient condition to obtain a proposal distri-
bution that makes a valid Metropolis sampler for F is to require that
Pprop is symmetric, that is Pprop (s∗ | sn ) = Pprop (sn | s∗ ). The two simple
proposal distributions that I mentioned in Step 1 have this property.
• The proposal distribution is the heart of a Metropolis algorithm. If it is
designed poorly, the sampling process will take a long time to “cover all the
grounds” of the probability landscape. Some thinking is needed to design
a good proposal distribution. If Pprop (s∗ | sn ) is allowed to often propose
candidates s∗ that are far away from sn , the danger is that one lands in a
low-probability zone of the landscape, which means that the candidate is
rarely accepted and the sample point sequence has many repeated points.
The net effect is slow coverage and long required sampling times. If the
proposal distribution mostly suggests candidates from the close vicinity of
the previous point, the resulting random walk will also be too slow. Thus, the
art lies in finding a proposal distribution which jumps far — but mostly hits
candidates that have a high probability. Finding such a proposal distribution
needs insight into the nature of the probability distribution that one wants
to sample from.
• Often one wishes to create a sampler which produces an independently, iden-
tically distributed (i.i.d.) sequence of sample points s1 , s2 , . . .. That is, the
choice of point sn+1 should be statistically independent of the previous point
sn — the sampling point sequence should have no “memory”. But the
163
proposal / acceptance mechanism of Metropolis sampling (and any other
MCMC technique) makes the choice of sn+1 depend on sn . The solution for
this problem is to subsample the sequence generated by the sampler, that is,
instead of using the original sequence s1 , s2 , . . ., retain only each h-th point,
recording only sh , s2h , . . ..
This is (please remember!) the form of the neural network training objective
in supervised learning, where the states s would be parameter vectors θ of a NN
and the cost function R would be the (empirical) risk.
164
But the minimization problem (83) is absolutely general and could, for in-
stance, mean the task to minimize financial loss in stock market transactions; or
to find a protein folding which minimizes the energy of the resulting 3D molecular
structure (which is what nature does, and which gives rise to the proteins you are
made of, and which is a major task in biochemistry research).
Solving (83) analytically is out of the question in many real-world optimization
tasks. If the cost function is differentiable one can try to solve (83) by gradient
descent, which is computationally often quite feasible but which can spot only local
cost minima. If one wants to find the global minimum, one has to “search” through
the entire search space S. A systematic grid search is infeasible if the states s are
high-dimensional vectors because the number of grid points in a search sub-volume
of Rd explodes exponentially with d. In this situation one needs a “clever” random
search method which explores the search space in a way that low-cost candidates
are (much) more often tried than high-cost ones, while making sure that the search
does not become trapped in some subvolume of the search space S but will visit
all parts of it.
Two particularly popular families of such “clever” stochastic search techniques
are evolutionary optimization with the special case of genetic algorithms, and
simulated annealing. Both approaches are inspired by nature:
165
material.
I’ll now give a more detailed description of SA, again for the case of discrete
search spaces where the Boltzmann distribution is characterized by a pmf. Here
is an outline of the complete process:
2. Identify the cost function, of which a global minimum should be found, with
the energy function E(s).
3. Start with some medium or high temperature T0 and consider the Boltzmann
distribution P (s, T0 ), which will be close to uniform (as in the first panel
of Figure 62). Start sampling from this distribution with the Metropolis
algorithm. The sequence of created samples will cover the search space
almost uniformly.
5. Continue until Tn ≈ 0. The hope is that then the cooling process has guided
you toward the global minimum and that the samples that you now get are
closely scattered around that global minimum.
A natural question at this point is, why not start immediately at low temper-
atures (e.g. in a situation like that shown in the 4th panel of Figure 62), wouldn’t
that just save the time of “bouncing aimlessly” around in the search space at high
temperatures, and instead directly lead you to the desired minimum, which is well
pronounced at low temperatures? The answer is, if one starts at low temperatures
— or, for that matter, if one cools to rapidly — one is likely to get “frozen” in
a very suboptimal local minimum far from the best one, from which one cannot
escape. This becomes intuitively clearer if we re-interpret the Metropolis sampling
of P (s, Tn ) in physical terms of “jumping around” in the energy landscape E(s)
directly.
To see this, we consider the two cases when (A) the Metropolis algorithm
accepts with certainty, and (B) when it accepts with probability Paccept (s∗ | sn ) =
P (s∗ ,T )
P (sn ,T )
, and see how we can translate this into energies.
Case A: This case occurs when P (s∗ , T ) ≥ P (sn , T ). This is equivalent to the
condition E(s∗ ) ≤ E(sn ). Thus, whenever the energy of the proposed state
is lower than the energy of the previous state, accept with certainty.
166
∗
Case B: If P (s∗ , T ) < P (sn , T ), then Paccept (s∗ | sn ) = PP (s
(s ,T )
n ,T ) . Rewriting the log
= exp(∆E/T ), (84)
where ∆E < 0 is the energy difference between sn and s∗ .
Summarizing we see that in terms of energy, a new proposed microstate is
accepted with certainty if its corresponding energy jump goes “downhill”, and
if it goes uphill, it is accepted with probability Paccept = exp(∆E/T ). That is,
the greater the energy increase, the accept (exponentially) more unlikely is such
a step taken; however, this may be compensated by a proportional increase in
temperature. In other words, on the average, at higher temperatures we may take
higher jumps uphill.
Equipped with this re-interpretation of the Metropolis algorithm in terms of
an energy-based acceptance function, we can better understand why slow cooling
is important for a final landing in a good local minimum of the energy landscape.
We can now intuitively interpret the SA search process as a random jump sequence
of a “search ball” in the energy landscape E(s), where the temperature determines
the ability of the ball to (randomly) climb uphill and in this way overcome “energy
barriers”. Figure 64 illustrates the different behavior of the SA search process at
different temperatures.
Figure 64: SA seen as an energetic ball game. At high temperatures, the ball
can more easily jump high and overcome energy barrier than at low
temperatures.
We are now aware that the cooling process is important for the success of
running an SA algorithm. A widely used, quick and dirty cooling scheme is ex-
ponential cooling: Put Tn+1 = k Tn for some k < 1 close to 1. Update T after
167
every single sample point. Clearly the size of k is important — a typical way to
optimize it would be just to experiment.
However, such simple cooling schemes, although widely used, may yield un-
satisfactory results. During a SA search run, one may encounter periods where a
particularly slow cooling is required, while at other periods, one may cool faster.
I will first explain this fact intuitively and then give a formal account.
Here are two intuitive examples that illustrate the necessity for slower-than-
others cooling periods. The temperatures where particularly slow cooling is re-
quired are associated with phase transitions, like when water freezes at tempera-
ture T = 273.15 Kelvin.
168
quite some effort (= mental re-heating) to escape. If you cooled to quickly
(= decided for a solution strategy too early, too quickly) your approach is
likely to fail (= lead to a quite suboptimal minimum); if you spent time to
consider different solution options (= hover around the critical temperature)
and then sloooowly decided, your chances of hitting a good solution strategy
are much higher.
F (T ) = ET [E] − T S(T ),
Thus the free energy relates the average energy at temperature T with the
entropy. Intuitively, the free energy of a system is the “useable” part of its energy,
energy that could be exploited at a macroscopic scale (for example, the free energy
of a volume of gas would be the energy that one could exploit by expanding the
gas in a piston, plus the energy that one might gain from cooling the volume).
Phase transitions are defined in physics as discontinuities in the free energy (or
one of its derivatives) as temperature (or another macroscopic variable) passes a
critical value (check out [Link]
For instance, as a volume of water is cooled from some ε value above zero Celsius
to some ε value below, such that it freezes, one has to extract a certain amount of
energy (the melting heat) from that volume — one may in fact exploit this energy;
it is part of the free energy of the volume of water. Therefore, the free energy of
the volume of water just above zero jumps discontinuously to a lower value as the
water is cooled to a temperature just below freezing.
Similarly, when running SA for an optimization problem, one can in principle
compute, at every step n, the free energy Fn of the system, and make the cooling
rate depend on the development of the free energy: cool slowly when Fn shows
signs of changing rapidly, or in other words, cool in a fashion such Fn decreases
smoothly.
The free energy can be computed from the partition function Z by F (T ) =
−T log(Z(T )). So the question is, how can one compute the partition function,
169
which is a gigantic sum (for discrete systems) or an intractable integral (for con-
tinuous systems)? The brutal answer is: use sampling (!) for an approximate
evaluation of this integral. There are a number of specialized sampling procedures
for the partition function, surveyed in Neal (1993) and Chapter 18 of Goodfel-
low et al. (2016). This method of steering the cooling is obviously very expensive:
within an SA run (which may have millions of steps), we repeatedly have to squeeze
in complete auxiliary sampling runs.
It can be mathematically shown that if the cooling is done slower, on average
across different SA runs one ends in lower-energy minima of the energy landscape.
In the limit of infinitesimally slow cooling, SA is guaranteed to find a global
minimum.
I conclude this section with two examples. The first one is in many ways
representative for many SA applications in combinatorial optimization tasks; the
second is just for fun.
170
chips is X aij X aij X aij
(si − sj )2 = − si sj + .
i>j
4 i>j
2 i>j
2
The second sum term is independent of the circuit placement and can be
dropped from the energy function, because it does not affect the location of its
minima.
For (ii), the objective function should grow with the degree of P imbalance of
2
circuits
P assigned to
P the two chips. The squared imbalance score ( s
i i ) is equal
to 2 i>j si sj + i s2i . Again, the second term is independent of the placement
and can be dropped.
Assembling these two cost contributions and replacing the constant 2 by a
weighting factor λ one gets a cost/energy function of the form
X X
E(s) = λ s i sj − aij si sj .
i>j i>j
171
Figure 65: Distributions of total number of pins at various temperatures. The
arrow inserted at the x-axis indicates the average number of pins ob-
tained by a greedy search algorithm (Metropolis sampling with T =
0). Figure taken from the Kirkpatrick et al paper.
172
Figure 66: Unscrambling a photo that was shattered into 1938 fragments. Left:
the mess that I gave out to students. Right: a low-energy re-ordering
found by running SA for 24 hrs on a PC. Right picture taken from the
project report of Ivaylo Enchev and Corneliu Prodescu (unpublished).
nitive systems? How can they be addressed or activated? And for what purposes
are the used, and how do they connect to other concepts? These questions brings
us to one of the big riddles of cognitive science. Numerous and diverse conceptions
of the nature of ’concepts’ have been developed in cognitive science, linguistics,
psychology, AI, philosophy and neuroscience. I highlight the complexity of this
question by pointing out some of the most important views on how concepts are
represented in brains or artificial AI systems. I explain this with with the example
of the concept of the digit pattern of the number four:
In order to recognize new incoming instances of the digit “four”, like this
one:
the stored prototype is matched against the new input, which is classified as
“digit pattern four” if the match is close enough. The prototype view is one
of the classical models of concept representations in cognitive science. The
173
’fundamental memory’ patterns stored in the energy minima of a Hopfield
network are a perfect case of prototype representations.
174
network models were a dominant school of thinking in neural networks. A
good example is the SHRUTI model of Shastri (1999). This linguist spent
his lifetime on developing a detailed connectionist model of human language
understanding.
Figure 67: A sample from a distribution of the pattern “four” (taken from the
widely used MNIST digit benchmark dataset).
175
are forever changing their shape, going through all sorts of variations, and these
variations ultimately give an exact representation of the probability distribution of
the digit-4 concept encoded in the BM. Again, if you think about it, you will find
that indeed your brain has this capability of confabulating ever fresh examples of
digit-4 patterns — and because you can do that, your brain must somehow have
a representation of a distribution of such digit-4 patterns. Agree?
I spent many years during my Phd and postdoc times thinking about the right
way how to scientifically think about the concept of concepts. I still am totally
confused. It’s your task as a new generation of thinkers to sort out these things
and settle the question of what is a ’concept’ for good!
7.1 Architecture
A BM is a recurrent NN whose neurons are all binary, that is, every neuron can
have an activation of 0 or an activation of 1. Neurons can be either visible or they
can be hidden. As we will shortly see, the visible neurons can be very flexibly used
in various ways (in the same BM) for input and output, whereas the hidden units
add internal computing power to the achievable input/output mappings.
We will use the following notation. For a BM with L visible units, v =
(v1 , . . . , vL )′ ∈ {0, 1}L is an activation vector of the visibles; and if the BM has M
hidden units, h = (h1 , . . . , hM )′ ∈ {0, 1}M is the activation vector of the hidden
units. Often we do not want to distinguish between these two sorts and write
s = (v1 , . . . , vL , h1 , . . . , hM )′ =: (s1 , . . . , sL+M )′ ∈ {0, 1}L+M for the entire network
state. Always the first L members of an entire network state will be reserved for
the visible units.
In a BM, there is an undirected synaptic link between any two visible and/or
hidden units si , sj : a BM is totally connected. Every link has a real-valued weight
wij = wji ∈ R. Special case: a zero weight wij = 0 amounts to “no link between
units i and j”. Self-connections wii 6= 0 are not allowed.
You will recognize the similarities with the Hopfield network: binary states,
symmetric weights, no self-connections. I remark without going into details that
this connects BMs to statistical physics and energy interpretations.
Thus, in summary, a BM is fully characterized (i) by its symmetric weight
matrix W of size (L+M )×(L+M ) with a zero diagonal and (ii) by the specification
which of the units are visible, that is by the number L.
A BM can be used for different purposes. Among others, it can be used for
the same tasks as MLPs, namely the supervised training of pattern classification
tasks. I now work out this use-case a little more.
In the supervised learning setting, the training data are (ui , yi )i=1,...,N where
the ui ∈ {0, 1}Lu are binary pattern vectors and the teacher outputs yi ∈ {0, 1}Ly
are binary classification vectors in one-hot encoding. The training objective is the
same as we know it from MLPs: Upon input of a new testing pattern utest from
pattern class j, the j-th output unit should become activated to a value of 1 and
176
the other output units should stay at zero activation. The inputs ui and outputs
yi together make for the visible units of the BM, which means that L = Lu + Ly .
And here is the first amazing new thing about BMs that makes them so differ-
ent from MLPs. A BM, when it was trained on a pattern classification task, can
also be run backwards in a generative confabulation mode. For instance, after it
has been trained to classify the ten classes of handwritten digits 1, 2, . . . , 9, 0, one
can “clamp” the ten classification neurons in one class, for instance in the class of
the digit 4, by fixing the output neurons to the state (0, 0, 0, 1, 0, 0, 0, 0, 0, 0)′ .
In other words, the ten classification units are now used as input units. Then, ex-
amples of the digit 4 patterns will be generated on what was the input layer for
classification, but is now the output layer. This generation of examples will be
driven by Metropolis sampling, that is, you will see a random sequence of various
patterns “4” appearing on the “input” layer. This sequence of patterns “4” is a
(Metropolis generated) sample from the learnt distribution of this pattern class.
The examples shown in Figure 67 might well have been collected from the input
layer while the classification layer was clamped to (0, 0, 0, 1, 0, 0, 0, 0, 0, 0)′ .
The “input” layer can thus be interpreted/used sometimes as an input “retina”
when the BM is run in classification mode, or at other times as an output projection
screen when it is run in confabulation mode. Likewise, the classification layer with
its ten class-coding neurons can be seen as output layer in classification mode, and
as input layer in the confabulation mode. This is why it is common thinking and
terminology in the BM world to drop the distinction between input and output
layers. Instead, the neurons in these layers are just called visible neurons. Figure
68 gives a schematic view of the structure of a BM that is structured as indicated
in this digit classification / generation set-up.
In plain English, this energy function gives low values if pairs of neurons
si , sj which both have an activation of 1 are connected by large positive weights.
Conversely, the energy is high when such si , sj are connected by strong negative
weights.
If during Metropolis sampling a unit si jumps from an activation of 0 to an
activation of 1, while all other units sj retain their activation, the energy (85)
177
“1”
“2”
“3”
“4”
.
.
.
Figure 68: Schematic of a Boltzmann machine. Only some of the all-to-all con-
nections are drawn.
This energy function gives rise to the Boltzmann distribution over the set
S = {0, 1}L+M of all possible states, which has the pmf
1 −E(s)
P (s) = exp( ). (87)
Z T
In most usages of a BM, the temperature is not changed and can be fixed at an
arbitrary value, typically T = 1, which simplifies (87) to P (s) = 1/Z exp(−E(s)).
Thus each particular fixed setting of the weights defines a particular energy
landscape and hence, a particular probability distribution over the possible bi-
nary network states (the microstates in statistical physics terminology). We write
EW , PW for the energy / probability distribution induced by a weight matrix W .
Sampling from the distribution PW is achieved with a special version of the
Metropolis sampler using the Boltzmann acceptance function. This sampler cycles
through the components si of the states s, updating only the selected component
to get a new sample point. This ansatz results in the following update rule:
When unit si is chosen for update from iteration n to n + 1, set it to a value of
178
1 in generation n + 1, regardless of its value at generation n, with the probability
1
P (sn+1
i = 1 | sn ) = , (88)
1 + exp(−∆Ei /T )
where −∆Ei is the energy increment from (86). Proving that this rule is an
instantiation of the Boltzmann acceptance function is a recommended elemen-
tary exercise (hint: start from considering the conditional probabilities P (si =
1 | state of all other units), P (si = 0 | state of all other units); note that the ratio
of these two probabilities is the same as the ratio P (s′ )/P (s) where s′ is the same
as s except for unit i where s′ has a value of 1 and s has a value of 0; exploit that
P (si = 1 | state of all other units) + P (si = 0 | state of all other units) = 1.)
179
(for a change, not for class “4”). The sampling process should then lead
to a sequence of random variations of handwritten-like “2” patterns on the
“retina” part of the visible units.
where pij is the average (over training samples) probability that the units i and
j are both active (that is, si = sj = 1) when the visible units are clamped to the
180
training data point v, and qij is the probability that units i and j are simultane-
ously active in a “free-running” sampling mode with no external clamping. This
yields the following update rule for weights:
where wij (n) is the value of the weight wij at the n-th step of gradient descent
and λ is a learning rate.
A single weight update wij (n) → wij (n + 1) thus involves the following opera-
tions:
1. Estimation of pij : for each training sample point vk , clamp the visible units
to vk . While the visible units remain clamped to vk , run the BM Metropo-
lis sampler (88) until a reasonably representative sample of network states
under this clamping condition has been collected. Use this collection to es-
timate the probability pkij of co-activity of units i and j in this clamping
condition. Do this for all training examples k = 1, . . . , N . Finally, set pij
to the average of all of these pkij . All of this is sometimes called the “wake”
phase of the BM learning algorithm (when the BM has its eyes open and
“sees” the visible input).
2. Estimation of qij : Similar, only without clamping the visible units. This is
sometimes called the “sleep” phase (the BM has its eyes closed).
181
On the other hand, according to the neural sampling hypothesis (Buesing et al.,
2011), biological brains might be carrying out sampling processes all the time, us-
ing random spikes as microstate update operators. This would be massively more
efficient than running the Metropolis sampler on a digital computer, because (i)
physical randomness comes for free, in contrast to producing pseudorandom num-
bers on a digital machine, which needs significant resources, and (ii) the microstate
updates at different neurons i happen in parallel, in contrast to the sequential
cycling through all neurons with the BM sampling algorithm. The Boltzmann
machine model might come amazingly close to what really happens in brains if
the neural sampling hypothesis holds true.
Furthermore, Boltzmann machines may enjoy a practically useful second life in
the emerging new scientific field called neuromorphic computing (it also has other
names, like physical computing or in-materio computing — the field is so new that
its name is still evolving like the field itself). In certain materials, which are being
developed by materials scientists, the atoms of the material can randomly switch
between two ’spin’ values (call them ’up’ and ’down’, or 0 and 1), and these switch
events are governed by ... the Boltzmann distribution! Such a chunk of spin-
flipping atoms looks and feels like a Boltzmann machine, with as many neurons as
there are atoms in the chunk, and with a physically built-in sampling ’algorithm’
that is ultra fast, fully parallel and ultra energy-efficient. The engineering challenge
is to make the energy landscape trainable, which amounts to making the physical
interaction strengths between neighboring atoms adaptable. This problem is not
yet solved for materials with atom-small ’spin-flipping’ units. But there is an
active line of research (google ’Ising machines’) where instead of atoms, larger
multi-atom nanoscale substructures are fabricated that admit a variant of spin-
flipping, and that admit trainable energy landscapes. Research of this kind is done
in the CogniGron research center, and indeed one of my PhD students (Fabian
IJpelaar) is working on the computational-theory side of enterprise.
• DBNs are layered neural networks, where each layer corresponds to one BM.
The hidden units of one such layer BM make the visible units of the next-
higher layer BM.
182
• The connectivity of each participating layer BM is very much reduced: there
are no within-layer connections, only connections between adjacent layers
exist. This led to the name “restricted” Boltzmann machines.
• Each sampling subroutine for the estimation of the probabilities pij , qij is
condensed to only two state updates of one layer BM, using a shortcut
approximate algorithm called contrastive divergence.
These innovations together have brought BMs back on stage with a flourish,
which in turn has triggered the deep learning revolution. I can’t possibly describe
DBNs better (nor more concisely) than Hinton & Salakhutinov did in their cele-
brated Science paper, so if you are interested in digging deeper into the deep roots
of deep learning, that paper is a must-read for you (only 3 pages).
Hinton and Salakhutinov did not foresee or plan what is now called deep learn-
ing. Their paper concentrates on a use of RBMs for data compression. Only in
a 14-line paragraph they mention in passing that they could use RBMs for pre-
training (initializing) the weights of a multilayer perceptron, giving a good starting
point for the subsequent application of the standard backpropagation algorithm.
The rest is history, and you pluck the fruit whenever you use Google translate or
play with ChatGPT.
DBNs remained fashionable for some years after 2006, and there were many
sophisticated attempts to lift them to a generally useful, stand-alone learning
approach for real-world modeling tasks. In fact, — at least, in my personal rec-
ollection — the entire machine learning community was thrilled about DBNs at
that time (check out the 2007 youtube video [Link]
v=AyzOUbkUf3M to see Hinton himself presenting the BM in a Google talk). In
the end, however, deep learning schemes which used cheaper initialization schemes
than via DBNs won over.
183
8 Reservoir computing
Reservoir computing (RC) is an umbrella term for a number of closely related
approaches for designing and training recurrent neural networks. RC methods are
in many ways complementary to the RNN methods which we met in Section 4.
While the computational and learning algorithms differ, the tasks solved by RC
networks are the same as those that are solved by the RNNs and LSTM networks
from Section 4: supervised timeseries-in, timeseries-out learning tasks of all sorts
and for all kinds of applications. I repeat the basic set-up from Section 4 for
convenience. The generic formulas for RNNs and supervised learning tasks are:
Training data: One or several, long or short pairs of discrete-time input and
output timeseries S = (u(n), y(n))n=1,...,nmax , where u(n) ∈ RK is the input
and y(n) ∈ RM is the output vector at time n.
Network equations:
Re-read Section 4.1 if you are unsure about what the symbols in these equa-
tions mean.
Learning task: Find weights W, W in , W out , b, all lumped together in one global
parameter vector θ, such that some loss function L defined to measure the
mismatch between the teacher output signal Ytrain = (y(n))n=1,...,nmax and
the RNN output Ŷθtrain = (ŷ(n))n=1,...,nmax of a network with weights θ is
minimized, that is, solve the minimization problem
The big difference between training RNNs in the “normal” way by gradient
descent via BPTT, and how training is done in RC, is that BPTT optimizes all
parameters in θ = {W, W in , W out , b}, whereas in RC only the readout weights
W out are trained. The other weights contained in {W, W in , b} are randomly set
at network design time and are then kept fixed forever. Figure 69 highlights this
difference.
The recurrent neural network that lies between the input and the output neu-
rons is called the reservoir in this field, and the weights W out from the reservoir
to the output neurons are called the readout weights or just the readouts.
Only training the readouts seems to be a very strong simplification and one
might expect that it will lead to a drastic reduction in achievable performance.
There are a number of reasons why reservoir computing methods are nonetheless
being utilized besides the BPTT training schemes which dominate in deep learning:
184
BPTT: all weights are trained RC: train only output weights
y(n) y(n)
...
+–
error
error
Figure 69: Highlighting the difference between the “normal” BPTT training and
RC training of an RNN. Arrows shown in red are trainable.
2. The training algorithm is numerically robust and there are no local minima
problems.
3. For many tasks RC networks yield solutions that are on a par with what
one gets from BPTT-trained RNNs. In some applications RC even system-
atically outperforms BPTT-based RNN training schemes.
5. RC can work with RNNs whose equations are not differentiable, which pre-
cludes BPTT training. In fact, one does not need a “neural network” be-
tween input and output at all — any kind of nonlinear dynamical system
can be plugged in for the “reservoir”. Specifically, one can use exotic, non-
digital and even non-electronic microchips as reservoirs. This has made RC
a leading paradigm in recent research in optical computing and other phys-
ical substrates which use non-electric nanoscale phenomena. RC is a major
line of approach in the emerging field of neuromorphic computing, which I
also briefly mentioned at the end of the Boltzmann machine Section 7.4.
A note on history and significance: The core idea of RC, namely to use a fixed,
non-trainable RNN and only train readouts, has been independently discovered
a few times. The earliest publication known to me is Kirby (1991), a contribu-
tion to a local low-key AI conference, which was immediately and totally forgotten,
followed after 1994 by a series of papers from Peter F. Dominey, a cognitive neuro-
scientist who identified RC-like circuits in the human brain (for instance, Dominey
(1995)). He used a biologically inspired algorithm for training the readouts which
185
had a low statistical efficiency, with the effect that this work was not taken up
in the machine learning quarters. In the year 2001 the RC principle was again
re-discovered (by myself) within a machine learning context, this time with an
efficient learning algorithm, and branded under the name of echo state networks
(ESNs) (Jaeger, 2001). At the same time it was also independently re-disovered
by Wolfgang Maass in a theoretical neuroscience context, based on biologically
detailed, spiking neuron models, and published under the name of liquid state
machines (LSMs) (Maass et al., 2002a). In those years, BPTT training of RNNs
was not very practical because numerical instabilities and vanishing gradient prob-
lems were not yet under control. ESNs became popular in those years in machine
learning, especially after the publication of Jaeger and Haas (2004) where ESNs
achieved accuracy levels on benchmark tasks of that time which were up to five
orders of magnitude more accurate than the state of the art.
The term “reservoir computing” established itself as an umbrella term for
ESNs, LSMs and some variants. In machine learning contexts, the term “echo
state networks” is still common. When it is used, it is implied that simple neural
networks with equations like (93) are used. When the word “liquid state machine”
is used, this usually means that the author treats a neuroscience modeling topic
and uses more involved, biologically motivated network models with spiking neu-
rons. I will concentrate on the machine learning aspects of RC and therefore use
the word “echo state network” in the next subsections.
In the decade until about 2015, the successful harnessing of BPTT in the deep
learning field diminshed the interest in RC in the machine learning community,
while it continued to be explored in neuroscience.
Since about 2015, both in academia and industry one could witness a quickly
growing interest in developing “brain-inspired” computing microchips. Digital
computing technologies will soon hit ultimate limits in miniaturization, and fur-
thermore the energy demands of digital IT technologies are becoming prohibitive
— it is estimated that today more than 10% of the world’s energy budget is eaten
up by digitial computing hardware. The biological brain is estimated to have an
energy efficiency that is four orders of magnitude better than what can be realized
with classical digital microchip designs. This is a strong economical and ecological
motif to explore non-digital, “brain-like” neuromorphic hardware solutions. The
rise of neuromorphic computing research has pulled RC back into the focus of
attention.
186
while the input pattern is of class C1 and is 0 when the input pattern is of type C2 .
In our demo, class C1 is a rectangular wave while class C2 is a sinewave. Figure
70 shows a portion of the training data.
0.5
u(n) 0
−0.5
−1
y(n) 0.5
In this demo I use a reservoir made of L = 100 neurons. Since the input and
output signals are one-dimensional (K = M = 1), there is one input neuron and
one output neuron. The L × L sized reservoir-internal weight matrix W , the L × 1
input weight matrix W in and the L-dimensional bias vector b from Equation 93
are filled with random values sampled from uniform distributions around 0, that
is, positive and negative values occur roughly equally often.
by the training input and has output weights W out , and let Ytrain = (y(n))n=1,...,nmax
train train
be the teacher output. Note that ŶW out and Y are (finite timeseries) vectors of
length nmax . Using the quadratic loss (which is the most common choice in RC),
the learning objective is to solve
out
Wopt train
= argmin kŶW out − Y
train 2
k. (96)
W out
187
1
0
8.1.2 Step 1: state harvesting
−1
Solving (96) is done in two steps. In the first step, the
1 reservoir network (which
has been randomly created) is driven by the teacher input,
0 that is, for a duration
of nmax steps. While it is being driven, we record the−1 activations xi (n) of each of
the reservoir neurons. This gives L timeseries of length
1
1 nmax . Figure 71 illustrates
this step, which is sometimes called the harvesting of0 reservoir states.
0
−1
−1
1
1
0 x1 (n)
0
−1
−1
1
1
0.5 1
0 0 x2 (n)
0
−0.5
−1
−1 −1
1
1
u(n) 0
0 x3 (n)
−1
−1
1
0
...
−1
1
1
0 xL(n)
0
−1
−1
1
−1
In words: the L vectors xi have to be linearly combined such that their weighted
sum best approximates the teacher vector Ytrain in the least mean square error
188
1
sense. This
0 is just a case of computing a linear regression. Every mathematical
or statistical
−1 programming toolbox offers a choice of ready-made algorithms for
computing a linear regression.
1
If you
0 are not familiar with the concept of linear regression, I can recommend
Section−13.1 in my lecture notes for the Machine Learning course. (In fact manda-
tory reading
1
1
if you need to refresh your understanding of linear regression. Nobody
should leave
0 a neural networks course without knowing what linear regression is!)
0
−1
−1
1
1
x1 0
0 w1out
−1
−1
1
1
x2 0
w2out
0
−1 1
+
−1
1
1
0
w3 out
≈ 0.5
x3 0 0
−1
−1
Ytrain
1
0
...
−1
1
1
xL 0 wLout
0
−1
−1
1
Figure 72:
0 Step 2 in the RC training scheme: compute optimal readout weights
−1 wiout which give the best mean-square error approximation to the
1 teacher Ytrain .
0
−1
Figure
1 72 gives a graphical impression of this situation. After this step 2, the
out
training0 is finished. The found weights Wopt are now inserted into the network
architecture,
−1
the training is finished and the network ready for testing and use.
8.1.4 Testing
For testing, the trained network is driven with fresh input data utest (n) for which
a correct reference output signal ytest (n) is known, and the network-generated
output ŷtest (n) can be compared to the reference signal. Figure 73 shows this for
our little demo example.
189
1
1
0.5
0 0.5
−0.5
−1 0
Figure 73: Testing the trained network (which now has the readout weights in-
stalled) with fresh input data. The reference output ytest (n) is shown
in black and the network output in red.
the case), the random initialization of all the fixed weights comes at a cost of
O(L2 ). The state harvesting needs nmax network updates, each of which has a cost
of O(L2 ). The linear regression costs O(L3 ). All of these together give a total cost
of O(L2 +nmax L2 +L3 ) = O(nmax L2 +L3 ). Normally one should have more training
data time points than reservoir neurons (because if one has more neurons than
training data points, one will get a zero-error solution from the linear regression,
which would mean overfitting in all except the most trivial tasks). Then the cost
is dominated by the harvesting phase and becomes O(nmax L2 ). That is, the cost
of training is essentially the same as running the network once with the training
input data. It cannot be cheaper.
8.2 RC in practice
The demo example illustrated all that needs to be done in training a RNN in
reservoir computing style. Seems simple... But if you want to squeeze good per-
formance out of ESNs, a number of design decisions need to be made appropriately,
and this is not so simple. After all, we are dealing with high-dimensional nonlin-
ear dynamical systems, and these are never easy to handle. In my experience it
takes some months for a RC novice, working full-time, to gain the insight and
routine necessary to handle RC techniques adequately. The “tricks of the trade”
are explained in the detailed practical RC tutorial by Lukosevicius (2012). Here I
mention some things that need to be considered. This summary account cannot
replace reading the tutorial of Lukosevicius (2012) if you seriously want to get
started with RC in practice.
Discard washouts. For state harvesting, the reservoir must be started in an ini-
tial state x(0). This state is arbitrary and unrelated to the learning task,
and traces of it will remain in the next few networks states. Figure 74 shows
this “initial state memory” effect. Our demo ESN was here driven twice
190
with the same input, but started from two different, randomly chosen ini-
tial network states. The top panel in Figure 74 shows the initial 10 step
traces of four reservoir neurons, with the set of traces from the first ini-
tial state in solid lines and the traces arising from the other initial state in
broken lines. One can see that for the first few steps these state sequences
differ from each other, but converge toward each other. The rate of con-
vergence is typically exponential on average, as can be seen in the bottom
plot. This plot was computed as follows. Let x(n), x∗ (n) denote the two
reservoir state sequences. The bottom plot shows the development of the
Euclidean distance between x(n) and x∗ (n) in log10 scaling, that is, it plots
log10 kx(n) − x∗ (n)k. The initial state differences are “forgotten”, in RC
terminology: they are washed out. How fast this washing-out happens de-
pends on many factors; it can be much slower than in this demo where the
differences become invisible in the top plot after three or four steps already.
0.5
−0.5
−1
−2
−4
−6
0 1 2 3 4 5 6 7 8 9 10
Figure 74: The initial state memory effect. For explanation see text.
In order not to pollute the linear regression in step 2 of the RC training pro-
cedure with the un-systematic effects of the arbitrary initial state setting,
one simply starts the harvesting after an initial washout period. An appro-
priate length of this washout time is most easily determined by computing
a diagnostic plot like our Figure 74.
Scaling of initial weights. This is the most important point. The geometrical
and dynamical properties of the harvested state sequences depend crucially
on the average absolute sizes of the weights in W in , W and the bias b. Figure
75 illustrates this.
Concretely, in this scaling exercise I always started from the same reference
in in
matrices Wref , Wref , bref . The entries of Wref and b were sampled from the
uniform distribution in the range [−0.5, 0.5]. In order to get Wref , another
method was used which is typical for practical work with ESNs. A prelimi-
nary version of Wref was created again by sampling weights from [−0.5, 0.5].
191
1
0
−1
1
A 0 m|m|m
−1
1
B 0 l | s |s
−1
1
C 0 l | l|m
−1
1
D 0
s |s | l
−1
1
E 0
s | l |m
−1
0 10 20 30 40 50 60 70 80
Figure 75: The effects of weight scaling. Top panel: the input signal. The other
plots A – E show traces of four neurons in runs with that input, where
the weights were scaled with different scaling combinations. From top
to bottom the three matrices Wref in
| Wref | bref were scaled large (l),
medium (m) or small (s) as indicated on the right hand side.
Then the absolute value |λmax | of the largest eigenvalue of that preliminary
matrix was calculated. This number, which can be computed for any square
matrix, is called the spectral radius of the matrix. It plays an important
role in the theory and practice of RC. The preliminary matrix was scaled by
1/|λmax | to give Wref , which thus had a spectral radius of 1.
in
Given these reference matrices Wref , Wref , bref , I scaled each of them with a
factor that was either small (order of 0.1) or medium (order of 1.0) or large
(order of 10). The exact scalings are not important here — I finetuned them
a little to get visually appealing graphics in Figure 75.
Here are some observations and comments. I refer to the cases A – E in that
figure:
Case A: When all weights are scaled to an intermediate range, the network
states behave in a way that normally works well in RC: their amplitudes
span much of the possible value range [−1, 1]; they are clearly influenced
by the driving input but the reservoir-internal interactions make them
192
markedly different from each other. Note that while such a behavior
is often good for getting accurate results, it is not universally the best
for all sorts of tasks.
Case B: When the input weights are large and all other weights are small,
the activation of a reservoir neurons will be dominated by the input
term in (93) and (almost) all neurons will exhibit activation traces that
look like scaled versions of the input. This is not a desirable behavior
in most applications, except sometimes when the desired output y(n)
at time n depends only on the input u(n) at the same time (no memory
effects needed). But in such cases, a feedforward network would be a
better choice than an RNN.
Case C: If both input weights and reservoir-internal weights are large, the
activations of neurons will typically be pushed toward the −1, +1 limits
of the tanh sigmoid and an almost binary “switching” dynamics inside
the reservoir results. This may be appropriate in extremely nonlinear,
quasi Boolean input-output learning tasks.
Case D: If the bias yields the dominating weight components, the reservoir
dynamics degrades toward constant values in each neuron. I cannot
imagine any interesting task where this would be beneficial.
Case E: If the reservoir-internal weights are large and the input and bias
not, then there is danger that the network dynamics falls prey to a
wild recurrent self-excitation which is no longer modulated in a useful
way by the input: chaotic dynamics (in the strict mathematical sense)
emerge. This is certainly useless because the “echo state property”, to
which I will devote a separate subsection below, is violated.
in
In summary, you see that the absolute and relative scalings of Wref , Wref , bref
exert a dramatic influence on the reservoir dynamics, and getting those
scalings right is crucial for a high-quality RC performance.
Unfortunately there is no general rule of how to set these scalings optimally.
It depends on the task and also on the dimensions K and L of the input
and reservoir. Like so often in practical work with neural networks, the
beginner has to spend a lot of time experimenting, and experienced users
will benefit from their well-honed intuitions. In any case, you should always
create indicative plots of state dynamics as in Figure 75 to get a “feeling” for
what is happening inside your reservoir. It’s like a doctor doing an X-ray.
Output feedback for signal generation tasks. Some tasks demand that the
trained network should generate an output signal. The basic example is to
train a network that has no input and a single output unit which should
yield a generated signal, for instance a sinewave oscillation. Such signal
generation tasks require that the generated output signal is fed back into
193
the reservoir through feedback weigths W fb . Using a linear output unit, the
network update equations are
x(n + 1) = σ(W x(n) + W fb y(n) + b), (98)
y(n + 1) = W out x(n + 1). (99)
The weights W, W fb , b are fixed at design time, the readouts W out are
trained. Figure 76 illustrates the set-up.
1
Wfb
0
−1 training
20 40 60 80 100
Wout
1 exploitation
0
−1
20 40 60 80 100
Figure 76: Training an ESN as a sinewave oscillator. At training time the output
weights, shown in dotted red lines, are not there yet. For explanation
see text.
The teacher signal y(n) is a sample of the desired generated signal, here a
sinewave (black signal in Figure 76). For state harvesting, the teacher signal
is written into the output node in a version that is delayed by one timestep:
in update cycle n the value written into the output unit is y(n − 1). This
time lag accounts for the relative time difference of the y signal that you
witness in (98) versus (99). The feedback weights W fb assume a role as input
weights in the state harvesting phase.
After the output weights are computed as usual with linear regression, the
network is ready for use. If all worked out well, started from a random initial
state the network will settle into an oscillation mode which after an initial
transient, where the random initial reservoir state is washed out, settles into
the desired oscillation. Figure 77 shows this.
One can also add separate input signals which modulate the generated output
signals. Figure 42 shows an example where the task was to generate a
sinewave output whose frequency is set by an input signal. When the input
signal has a high value, the generated sine should have a high frequency.
The demo illustrated in Figure 42 was based on an ESN.
ESNs are particularly well suited for signal generation tasks. Specifically,
they are unrivalled when it comes to generate chaotic signals. The paper
which popularized ESNs (Jaeger and Haas, 2004) included demos of gener-
ating chaotic signals whose precision was essentially machine precision. In
194
1
Some reservoir
0
neuron traces
Wfb
−1
0 20 40 60 80 100
Wout
1
network output and
0 aligned version of
reference sinewave
−1
0 20 40 60 80 100
Figure 77: After training, the network functions as an oscillator. The output unit
will (after an initial washout resulting from the arbitrary initial reser-
voir state) generate the kind of signal that was used as teacher. The
network output (red) is overlaid with a suitably phase-shifted version
of the training sine (thin black line). The reservoir neurons display all
kinds of oscillations whose shapes result from the nonlinear interactions
of the reservoir neurons (four traces shown).
Leaky integrator neurons. (This part in blue font is optional reading) We have
seen in Figure 75 that reservoir dynamics can have quite different properties.
For good RC results, the dynamical and geometrical properties of the reser-
voir should match the properties of the task. Such properties which should
match are, for example, the degree of nonlinearity, the length of memory
timespans, or the sheer complexity of the task.
One of the most important properties is something that one could call the
“speed” of the system from which the training data come. There are fast sys-
tems and there are slow systems. A gigahertz computer clock is faster than
the tides of the ocean. While it is intuitively clear that there are fast and slow
systems, there exists no universal best mathematical way to define or mea-
sure the “speed” of a dynamical system. I will not dig deeper here although
I would love to do so — coming to terms with timescales was the core of my
research in the European Project MeM-Scales (“Memory Technologies with
multi-scale time constants for neuromorphic architectures”, [Link]),
and Jaeger and Catthoor (2023) is a survey report on timescales in neural
networks, microchips, and machine learning.
In practical applications of RC a crucial factor for success is to adjust the
reservoir’s “speed” to the task’s “speed”. To make this possible, one needs
195
a way to design slow or fast reservoirs.
The standard approach is to use a neuron model where each neuron’s dy-
namics is governed by the leaky integration update equations. In fact, RC
experts will almost always use leaky integrator neurons. Such neurons can
be slowed down or sped up as one wishes.
They are best explained by starting from a continuous-time neuron model.
Consider a reservoir that runs in continuous time t ≥ 0. The L-dimensional
activation vector at time t is x(t) (notice that often one uses symbol n ∈ Z
for discrete time points and t ∈ R for continuous time). The activation xi (t)
of neuron i at time t is a continuous-time signal. For a leaky integrator
neuron, it is governed by the ODE
!!
1 XL XK
ẋi (t) = −xi + σ wij xj (t) + in
wik uk (t) + bi , (100)
c j=1 k=1
where the wij are the synaptic connection weights of the incoming connec-
tions of neuron i, uk (t) is the k-th component of the input vector u(t) at
in
time t, the wik are the input connection weights leading into neuron i, bi is
this neuron’s bias, and — the most interesting bit in this equation — c is
the time constant of this ODE.
If this ODE is integrated with different settings of the time constant c, the
obtained trajectory that one could plot in a phase portrait will look the
same. Remember that the trajectories follow tangentially the vectors of the
vector field given by the right hand side of the ODE. The effect of changing
the time constant c only scales the length of these vectors but not their
direction. If the time constant grows, the length of the vectors in the vector
field shrinks in inverse proportion. This means that the “forward speed” of
the point xi (t) along the trajectory slows down when c grows. By setting
c one can create slow or fast continuous-time reservoir dynamics at one’s
discretion.
This dynamical law
P is called the leaky integration model
of a neuron be-
L PK in
cause the term σ j=1 wij xj (t) + k=1 wik uk (t) + bi integrates the input
to this neuron over time, while the term −x always lets the activation x(t)
diminish (“leak”) at a rate that is proportional to the current activation.
This is, by the way, an effect that is also active in biological neurons: due
to elementary electrophysics, their electric potential (measured in millivolt)
likewise would dissipate at a rate proportional to its current level, because
the cell membrane is not a perfect insulator. Therefore, leaky integrator
models are a much better fit to biological neurons than the simple neuron
model that we have been using in this section so far. In computational
neuroscience one almost always uses leaky integrator models of various sorts.
196
However, for practical machine learning applications run on digital com-
puters one needs a discrete-time neuron model. Instead of a continuous
trajectory (xi (t))0≤t one needs a timeseries which advances in discrete steps
with a chosen stepsize ∆, that is one wants to have a discrete-time version
(x̃i (n ∆))n=0,1,2,... of the continuous trajectory. At times t = n ∆ the two
trajectories should be (approximately) equal, xi (n ∆) ≈ x̃i (n ∆).
Finding discretization methods which allow such a discretization of ODEs
with a good balance between approximation accuracy and computational
cost is a main subject of numerical mathematics. Whenever you simulate an
ODE on a digital machine — which is the daily bread and butter in all the
natural sciences — a numerical ODE solver is invoked. Every mathematical
toolbox offers a choice of such solvers.
The simplest of all ODE solvers is called the Euler method (see https:
//[Link]/wiki/Euler_method for an easy introduction). After
some point x̃i (n ∆) has been computed, the next point x̃i ((n + 1) ∆) is
computed by following the direction of the vector given by the right-hand
side of the ODE for a timespan of ∆:
Using ∆ = 1, renaming 1/c to a and joining the two x̃i terms leads to the
simpler looking version
!
X
L X
K
x̃i (n + 1) = (1 − a) x̃i (n) + a σ wij x̃j (n) + in
wik uk (n) + bi , (101)
j=1 k=1
which is the form that you will most commonly find in the RC literature
(without the tilde). Now you can control the “speed” of neurons by setting
a ∈ [0, 1]: the larger this is set, the faster the neuron. In the extreme case
a = 1 one recovers our accustomed simple update equation (93). The slowest
“dynamics” is obtained with a = 0: then nothing happens — the network
state remains frozen in its initial state. The number a is often called the
leaking rate of the neuron.
One can set leaking rates individually for different reservoir neurons. This is
advisable if one wants to have a reservoir that can respond to complex input
signals that have both fast and slow components. Often it is good enough
to tune a single leaking rate uniformly for all reservoir neurons. This gives
the following matrix form for the state update equation:
197
Many real-world physical systems and all truly cognitive systems operate
on several timescales simultaneously. For instance, in atmospheric dynamics
small whirls of air (like gust eddies between houses) have a typical time
constant of a few seconds, while the large whirl of a low over central Europe
evolves over several days. Or, for a cognitive dynamics example, while you
are reading this lecture notes section, your brain at any moment has to
integrate information bits that come just from the preceding syllable (read a
few milliseconds before) with information from the beginning of this section
(maybe an hour ago). As of today, neither deep learning LSTM networks
trained with BPTT, nor ESNs are capable of integrating information across
many timescales. Extending the multi-timescale capabilities of RNNs is a
major topic of current research. In ESNs, one approach is to design reservoirs
with leaky integrator neurons, where different submodules of the reservoir
have different values for the inverse time constant a. Typically, fast modules
or processing layers are closer to the input than slow layers. An example is
Gallicchio et al. (2018) where “deep” ESNs are constructed in this way for
speech recognition and music composition.
Regularization and reservoir size. Like any machine learning method, RC is
susceptible to overfitting. In order to cope with this problem, the flexi-
bility of a RC training scheme must be adapted to the available training
data volume and task complexity by experimentation in a cross-validation
scheme (recall this from Section 1.4). This needs a way to tune the degree
of modeling flexibility. The recommended approach with ESNs is
• use a large reservoir, choosing a size that would allow overfitting,
• then use a regularized version of linear regression for the readout weight
calculation, called ridge regression.
Ridge regression, also known as Tikhonov regularization ([Link]
org/wiki/Tikhonov_regularization) should always be used in machine
learning when a linear regression has to be carried out — which means, all
over the place, not only in a RC context. It is a very valuable thing to know,
and I will explain it in a little detail and present an example.
Let us first rehearse the maths of linear regression, framed in an ESN train-
ing scenario. Assume we want to compute the 1 × L-dimensional readout
weight vector w for a single output neuron from harvested L-dimensional
reservoir states x(1), . . . , x(nmax ) such that the mean squared training error
is minimized, that is, we want to solve the problem
X
nmax
wopt = argmin (w x(n) − y(n))2 , (102)
w
n=1
where y(n) is the teacher for the output unit (note that w is a row vector).
The analytical solution for this problem is the well-known solution formula
198
for linear regression problems,
′
wopt = (X′ X)−1 X′ y, (103)
where X is the nmax × L sized matrix containing the state vectors x(n) in
its rows and y is the nmax × 1 vector y = (y(1), . . . , y(nmax ))′ containing the
teacher output values.
In ridge regression, the objective function in (102) is augmented by a regu-
larization term which penalizes the sum of squared weights,
! !
X
nmax
wopt = argmin (w x(n) − y(n))2 + α w′ w , (104)
w
n=1
199
1
training
0
signal
−1
1
0 ! = 10000
−1
1
0
! = 20
−1
1
!=0
0
−1
0 20 40 60 80 100
Figure 78: The training signal (top, blue lines) and the network output (red) in
three regularization conditions. The thin black sines drawn into the
panels are the clean sinewaves; they are shown for visual intuition only
and the training and testing procedures did not use this information.
200
1 test MSE
0.8
0.6
train MSE
0.4
0.2 weight size
0
10000 1000 100 20 10 1 0
!
Figure 79: THE paradigmatic textbook behavior of training and testing errors as
they vary with the degree of regularization. For explanation see text.
The black curve shows how the mean absolute output weights grow
as the regularization strength is decreased. For α = 0 the average
absolute weight size is about 5.2 (not drawn).
• The reservoir is continually fed with the input data stream u(n), leading to
a reservoir state sequence x(n).
• Assume that at time n, an output weight vector wout (n) is in place which at
the time around n produces outputs ŷ(n) which give a good approximation
(in the mean square error sense) to the teacher signal y(n). In the update
n → n+1, the online algorithm uses the next teacher value y(n+1) to adjust
wout (n) a little, by doing a small step down the error gradient, that is, add
a little correction vector to wout (n) which leads it in that direction which
most strongly reduces the squared error (y(n + 1) − ŷ(n + 1))2 . Skipping the
maths (derivation is straightforward, if interested you find it in Section 11.3
in my Machine Learning lecture notes, or in hundred other textbooks and
201
online tutorials), this leads to the update equation
That’s it. This very cheap and simple rule sits behind most of the signal
processing procedures which make the radio receiver module in your smartphone
work; neuroscientists say that it is biologically plausible; and it enables reservoir
computing to function well in online adaptive task settings, where deep learning
methods are severely challenged. If you have a very good memory of the earlier
parts of these lecture notes you will recognize that the Perceptron learning rule is
in fact just a version of this algorithm with λ = 1.
Like so many powerful ideas in engineering and science, this rule has been
discovered independently several times in different disciplines, where it is known
under different names. In signal processing and control it is called the LMS al-
gorithm, in the neurosciences and sometimes in artificial neural network research
it is called the Widrow-Hoff rule or the Delta-rule, and a plain mathematician
would likely refer to it as stochastic gradient descent on a quadratic error surface.
No introductory course on neural networks would be complete without it, and I
seized the opportunity to introduce it in an ESN context.
202
1
0.5
spectral radius r = 0.95:
0
reservoir has ESP
−0.5
−1
0.5
−1
0 5 10 15 20 25 30 35 40 45 50
Figure 80: The echo state property: to have or have not. Top: with a scaling of the
reservoir weight matrix under a critical value, the ESP is granted and
state sequences converge to each other after some washout time, when
the reservoir is started from different initial states with the same input.
Four reservoir neurons are plotted, with the two runs distinguished by
solid / dashed lines. The input was a sinewave (not shown). Bottom:
the ESP is lost when the reservoir weight matrix is scaled up beyond
a critical value. Here the reservoir engages in a chaotic dynamics that
has lost the connection to the input signal.
There is no known way to predict the critical spectral radius ϱcrit analytically.
A lot of effort has been spent and is being spent on the mathematical study of the
ESP. There are two reasons why the ESP is attracting so much attention:
• For machine learning applications of RC it is obviously important because
RC training only works for input-reservoir combinations that exhibit the
ESP.
• In the general research on cognitive neurodynamics, a currently popular hy-
pothesis states that, roughly speaking, the human brain works at its best
when it increases its self-excitiation level just below the point where it bi-
furcates into chaos. Google ``edge of chaos'' neural network to get a
glimpse of the lively research on this topic, both in machine learning and
neuroscience quarters.
203
Since a few years, interest in reservoir computing has been re-kindled. The rea-
son is that RC is one of the few computational approaches which do not need dig-
ital computers. Non-digital, “neuromorphic”, “brain-inspired”, “unconventional”
computing microchips are gaining relevance due to the reasons that I briefly men-
tioned in the last paragraph before subsection 8.1. The Groningen Cognitive Sys-
tems and Materials Center (CogniGron) ([Link]
cognitive-systems-and-materials/), which was founded at RUG in the year
2018, strives to become a European pioneer in this field. The principle of reservoir
computing, which can be summarized as
“use an input signal to drive a nonlinear excitable medium — the “reservoir”
— and generate a desired output signal by combining many of the local response
signals that can be observed in the medium,”
can be applied to many kinds of “reservoirs” other than neural networks simulated
on digital machines. The recent survey of Tanaka et al. (2019) gives an overview.
If the reservoir is a real, physical piece of material which can be excited into inter-
esting dynamical responses by some physical driver input, one speaks of physical
reservoir computing. The potential benefits are inherent parallelism, a low energy
budget (in physical nanoscale devices), high speed, extreme high-dimensionality
of the reservoir states (even, in principle, infinite-dimensional states in continu-
ous materials), and, hopefully, low cost. As of today, all of this is still academic
research. I will not go into details but just illustrate the flavor of this kind of
research with a few examples. I present them by figures with detailed captions.
204
Figure 81: Reservoir computing in a bucket. Wolfgang Maass, the inventor of
the “liquid state machine” version of RC, used to call the reser-
voir “the liquid” (Maass et al., 2002b). This was taken literally by
students Chrisantha Fernando and Sampsa Sojakka at the School
of Cognitive and Computer Sciences, University of Sussex. They
filled a transparent acryl basin with real water, excited it with eight
Lego-made mechanical pushrods whose oscillations were derived from
speech signals, optically recorded the states of the water surface rip-
ples, and used these states as reservoir states. Their “liquid brain”
could solve the infamous XOR task and classify spoken “Zero” ver-
sus “One”. Their paper at the ECAL 2003 (Fernando and Sojakka,
2003) won the highest impact paper award. Chrisantha, now a Se-
nior Research Scientist at Google DeepMind, explains it on youtube
([Link]
Appendix
A Elementary mathematical structure-forming op-
erations
A.1 Pairs, tuples and indexed families
If two mathematical objects O1 , O2 are given, they can be grouped together in a
single new mathematical structure called the ordered pair (or just pair) of O1 , O2 .
It is written as
(O1 , O2 ).
In many cases, O1 , O2 will be of the same kind, for instance both are integers. But
the two objects need not be of the same kind. For instance, it is perfectly possible
to group integer O1 = 3 together with a random variable (a function!) O2 = X7
in a pair, getting (3, X7 ).
The crucial property of a pair (O1 , O2 ) which distinguishes it from the set
{O1 , O2 } is that the two members of a pair are ordered, that is, it makes sense to
speak of the “first” and the “second” member of a pair. In contrast, it makes not
sense to speak of the “first” or “second” element of the set {O1 , O2 }. Related to
205
Figure 82: Reservoir computing in a plastic worm. In robotics, the control of
body motion needs (among other items) a “forward model” of how
the body limbs will react to motor or muscle action. While in clas-
sical industrial robotics this forward model can be calculated analyt-
ically with high precision, this is not possible with soft robot bodies
or body parts — snakes, worms, trunks, tongues. One way to get
such a model nonetheless is to use the very physical body itself as a
reservoir. Its states are observed by sensors placed on it. The twofold
charm of this approach is, first, that this enables almost delay-less on-
line computing, and second, that this reservoir naturally has exactly
the right dynamical properties to “model” itself. At the University
of Zurich, a team around Rolf Pfeifer realized this idea with a plastic
worm. Their paper Nakajima et al. (2015) comes with a youtube demo
([Link] The idea to use
robot bodies themselves as computational modules within “algorithms”
needed for motion control has become branded as morphological com-
putation — google it if you find this interesting.
206
Figure 83: Reservoir computing in mechanical silicon microchips. Julien Sylvestre
and his group at the Department of Mechanical Engineering, Université
de Sherbrooke, Canada, explores how the mechanical oscillatory dy-
namics of freely suspended microscale silicon beams can be exploited in
microchips which combine mechanical sensing with RC signal process-
ing. The microbeams (marked red in image) can be etched into the sili-
con wafer with standard microchip fabrication technologies (Coulombe
et al., 2017). Several such microbeams on the same chip interact non-
linearly with each other by mechanical couplings.
207
Figure 84: Reservoir sensing. A potentially quite promising future application for
RC is biochemical and environmental sensing. The reservoir is here
a carrier plate coated with some chemical or biochemical material or
mix of materials which change their properties when the surface is ex-
posed to (traces of) chemical or biological substances whose presence
or concentration has to be measured. The property changes induced
in the active coating can be amplified and “dynamified” by additional
electrical impulses given to the plate. The resulting spatiotemporal
dynamics are recorded from the plate in some way, for instance elec-
trically, and used as reservoir states. The desired measurement signal
is trained by the RC principles. This line of sensor engineering was ex-
plored in a European FET-OPEN project (RECORD-IT, 2015-2018)
coordinated by Zoran Konkoli from the Department of Microtechnol-
ogy and Nanoscience, Chalmers University of Technology, Gothenburg,
Sweden. An interim report is Konkoli (2016).
208
Figure 85: Optical reservoir computing. Optical computing — that is, computing
with light instead of with electricity — is a large field of communication
engineering, comparable in importance and promises with quantum
computing. The potential benefits of optical versus electrical comput-
ing are (i) speed: the natural time constants of optical devices are many
orders of magnitude smaller than of electronic devices; (ii) 3D wiring:
unlike electrical signals, which need wires that lead to headaches in
microchip design because they must not cross each other, light beams
can cross each other without interference, (iii) potentially extremely
low energy consumption. These potential benefits have boosted optical
computing research at a large scale, but breakthroughs are still missing
— same as in quantum computing. One approach in this field is op-
tical reservoir computing. This has developed into the currently most
important branch among the many versions of physical RC (google
optical reservoir computing). A wide spectrum of optical effects
and reservoir architectures is being explored. The image shows a mi-
crophotograph of an optical microchip (real size 16 mm2 ) developed
at the University of Gent which implements an all-optical reservoir
with 16 “neurons”. It was demonstrated in Vandoorne et al. (2014)
that with this reservoir one could realize 5-bit header recognition in
internet packages. There are two aspects of this chip that fascinate
me. First, the coils that you see in the image are long spirals of silicon
waveguides whose function is to slow down the reservoir dynamics by
inserting lengthy light travel paths between the “neurons”. Without
these slow-down coils, the native processing speed of this chip would
be orders of magnitude too fast for feeding and analysing I/O signals
with electronic lab equipment. Second, this chip is entirely passive: it
needs no extra energy besides the energy in the incoming light signals.
— Interestingly, very recently, also the quantum computing field has
discovered RC as a potential venue for progress.
209
this is the fact that the two members of a pair can be the same, for instance (2, 2)
is a valid pair. In contrast, {2, 2} makes no sense.
A generalization of pairs is N -tuples. For an integer N > 0, an N -tuple of N
objects O1 , O2 , . . . , ON is written as
(O1 , O2 , . . . , ON ).
1-tuples are just individual objects; 2-tuples are pairs, and for N > 2, N -tuples
are also called lists (by computer scientists that is; mathematicians rather don’t
use that term). Again, the crucial property of N -tuples is that one can identify its
i-th member by its position in the tuple, or in more technical terminology, by its
index. That is, in an N -tuple, every index 1 ≤ i ≤ N “picks” one member from
the tuple.
The infinite generalization of N -tuples is provided by indexed families. For
any nonempty set I, called an index set in this context,
(Oi )i∈I
Oi
to denote the ith “component” of (Oi )i∈I . Writing Oi is a shorthand for applying
the ith projection function on (Oi )i∈I , that is, Oi = πi ((Oi )i∈I ).
S1 × . . . × SN = {(s1 , . . . , sN ) | si ∈ Si }.
210
Q Q
If all the sets Si are the same, say S, then the product i∈I Si = i∈I S is
also written as S I .
An important special case of infinite products is obtained when I = N. This
situation occurs universally in modeling stochastic processes with discrete time.
The elements n ∈ N are the points in time when the amplitude of some signal is
measured. The amplitude is a real number, so at any time n ∈ N, one records an
amplitude value an ∈ Sn = R. The product set
Y [
Sn = {φ : N → Sn | ∀n ∈ I : φ(n) ∈ Sn } = {φ : N → R}
n∈N n∈N
is the set of all right-infinite real-valued timeseries (with discrete time points start-
ing at time n = 0).
f1 ⊗ . . . ⊗ fN : D → S1 × . . . × SN
d 7→ (f1 (d), . . . , fN (d)).
211
P (X = gi , Y = aj ) a1 a2 a3 a4
g1 0.005 0.3 0.2 0.04
(107)
g2 0.005 0.15 0.15 0.04
g3 0.0 0.05 0.05 0.01
The cell (i, j) in this 3 × 4 table contains the probability that a customer with
gender gi falls into the age bracket aj . This is the joint probability of the two
observation values gi and aj . Notice that all the numbers in the table sum to 1.
The mathematical tool to formally describe a category of an observable value is
a random variable (RV). We typically use symbols X, Y, Z, . . . for RVs in abstract
mathematical formulas. When we deal with concrete applications, we may also
use “telling names” for RVs. For instance, in Table (107), instead of P (X =
gi , Y = aj ) we could have written P (Gender = gi , Age = aj ). Here we have two
such observation categories: gender and age bracket, and hence we use two RVs
X and Y for gender and age, respectively. In order to specify, for example, that
female customers in the age bracket 11-20 occur with a probability of 0.3 in the
shop’s customer reservoir (the second entry in the top line of the table), we write
P (X = g1 , Y = a2 ) = 0.3.
Some more info bits of concepts and terminology connected with RVs. You
should consider a RV as the mathematical counterpart of a procedure or apparatus
to make observations or measurements. For instance, the real-world counterpart of
the Gender RV could be an electronic questionnaire posted by the online shop, or
more precisely, the “what is your age?” box on that questionnaire, plus the whole
internet infrastructure needed to send the information entered by the customer
back to the company’s webserver. Or in a very different example (measuring
the speed of a car and showing it to the driver on the speedometer) the real-
world counterpart of a RV Speed would be the total on-board circuitry in a car,
comprising the wheel rotation sensor, the processing DSP microchip, and the
display at the dashboard.
A RV always comes with a set of possible outcomes. This set is called the
sample space of the RV, and I usually denote it with the symbol S. Mathematically,
a sample space is a set. The sample space for the Gender RV would be the set
S = {m, f, o}. The sample space for Age that we used in the table above was S =
{{0, 1, . . . , 10}, {11, . . . , 20}, {21, . . . , 30}, {31, 32, . . .}}. For car speed measuring
we might opt for S = R≥0 , the set of non-negative reals. A sample space can be
larger than the set of measurement values that are realistically possible, but it
must contain at least all the possible values.
Back to our table and the information it contains. If we are interested only in
the age distribution of customers, ignoring the gender aspects, we sum the entries
in each age column and get the marginal probabilities of the RV Y . Formally, we
compute
212
X
P (Y = aj ) = P (X = gi , Y = aj ).
i=1,2,3
a1 a2 a3 a4
g1 0.005 0.3 0.2 0.04 0.545
g2 0.005 0.15 0.15 0.04 0.345 (108)
g3 0.0 0.05 0.05 0.01 0.110
0.01 0.5 0.4 0.09
Notice that the marginal probabilities of age 0.01, 0.5, 0.4, 0.09 sum to 1, as do
the gender marginal probabilities.
Finally, the conditional probability P (X = gi | Y = aj ) that a customer has
gender gi given that the age bracket is aj is computed through dividing the joint
probabilities in column j by the sum of all values in this column:
P (X = gi , Y = aj )
P (X = gi | Y = aj ) = . (109)
P (Y = aj )
There are two equivalent versions of this formula:
P (X = gi , Y = aj ) = P (X = gi | Y = aj )P (Y = aj ) (110)
where the righthand side is called a factorization of the joint distribution on
the lefthand side, and
P (X = gi , Y = aj )
P (Y = aj ) = , (111)
P (X = gi | Y = aj )
demonstrating that each of the three quantities (joint, conditional, marginal prob-
ability) can be expressed by the respective two others. If you memorize one of
these formulas – I recommend the second one – you have memorized the very
key to master “probability arithmetics” and will never get lost when manipulating
probability formulas.
The factorization (110) can be done in two ways: P (Y = aj | X = gi )P (X =
gi ) = P (X = gi | Y = aj )P (Y = aj ), which gives rise to Bayes’ formula
P (X = gi | Y = aj )P (Y = aj )
P (Y = aj | X = gi ) = , (112)
P (X = gi )
which has many uses in statistical modeling because it shows how one can revert
the conditioning direction.
Joint, conditional, and marginal probabilities are also defined when there are
more than two categories of observations. For instance, the online shop marketing
213
people also record how much a customer spends on average, and formalize this by
a third random variable, say Z. The values that Z can take are spending brackets,
say s1 = less than 5 Euros to s20 = more than 5000 Euros. The joint probability
values P (X = gi , Y = aj , Z = sk ) would be arranged in a 3-dimensional array
sized 3 × 4 × 20, and again all values in this array together sum to 1. Now there
are different arrangements for conditional and marginal probabilities, for instance
P (Z = sk | X = gi , Y = aj ) is the probability that among the group of customers
with gender gi and age aj , a person spends an amount in the range sk . Or P (Z =
sk , Y = aj | X = gi ) is the probability that in the gender group gi a person is aged
aj and spends sk . As a last example, the probabilities P (X = gi , Z = sj ) are the
marginal probabilities obtained by summing away the Y variable:
X
P (X = gi , Z = sj ) = P (X = gi , Y = ak , Z = sj ) (113)
k=1,2,3,4
So far I have described cases where all kinds of observations were discrete, that
is, the respective sample spaces S were finite (for example, three gender values)
or countably infinite (for example, the natural numbers 1, 2, 3, . . .). The function
P : S → [0, 1] which assigns to each possible outcome s ∈ S its probability P (s)
is called a probability mass function (pmf) and we denote P it with an upper-case
P . The sum of the pmf over all possible outcomes is one: s∈S P (s) = 1. If one
behaves mathematically very correctly, the symbol P is indexed with the random
variable that gives rise to the distribution, that is one would write PX (s). This is
often not done out of convenience.
Equally often one faces continuous random values which arise from observations
that yield real numbers – for instance, measuring the body height or the weight
of a person. Since each such RV can give uncountably infinite many different
observation outcomes, their probabilities cannot be represented in a table or array,
and they cannot be summed up. Instead, one uses probability density functions
(pdf’s) to write down and compute probability values. We denote pdfs by lower-
case p.
In order to explain pdfs, let’s start with a single RV, say H = Body Height.
Since body heights are non-negative and, say, never larger than 3 m, the distribu-
tion of body heights within some reference population can be represented by a pdf
f : [0, 3] → R≥0 which maps the interval [0, 3] of possible values to the nonnegative
reals (Figure 86). We will be using subscripts to make it clear which RV a pdf
refers to, so the pdf describing the distribution of body height will be written fH .
A pdf for the distribution of a continuous RV X can be used to calculate the
probability that this RV takes values within a particular interval, by integrating
the pdf over that interval. For instance, the probability that a measurement of
body height comes out between 1.5 and 2.0 meters is obtained by
Z 2.0
P (H ∈ [1.5, 2.0]) = fH (x)dx, (114)
1.5
see the shaded area in Figure 86. Some comments:
214
2
1.5
pdf
0.5
0
0 0.5 1 1.5 2 2.5 3
H (body height)
• Be aware that the values f (x) of a pdf are not probabilities! Pdf’s turn into
probabilities only through integration over intervals.
• Values f (x) can be greater than 1 (as in Figure 86), again indicating that
they cannot be taken as probabilities.
215
0.8
0.6
0.4
pdf
0.2
0
3
1
2 3
Y 0 1
0
X
integrates to 1, i.e.
Z b1 Z bk
... f (x1 , . . . , xk ) dxk . . . dx1 = 1,
a1 ak
where also the cases ai = −∞ and bi = ∞ are possible. A more compact notation
for the same integral is Z
f (u) du,
D
where D denotes the k-dimensional box [a1 , b1 ]×. . .×[ak , bk ] and u denotes vectors
in Rk . Mathematicians speak of k-dimensional intervals instead of “boxes”. The
set of points S = {u ∈ Rk | fX1 ,...,Xk > 0} is called the support of the distribution.
Obviously S ⊆ D.
In analogy to the 1-dim case from Figure 86, probabilities are obtained from
a k-dimensional pdf fX1 ,...,Xk by integrating over sub-intervals. For such a k-
dimensional subinterval [r1 , s1 ] × . . . × [rk , sk ] ⊆ [a1 , b1 ] × . . . × [ak , bk ], we get its
probability by
Z s1 Z sk
P (X1 ∈ [r1 , s1 ], . . . , Xk ∈ [rk , sk ]) = ... f (x1 , . . . , xk ) dxk . . . dx1 . (115)
r1 rk
In essentially the same way as we did for discrete distributions, the pdf’s of
marginal distributions are obtained by integrating away the RV’s that one wishes
to expel. In analogy to (113), for instance, one would get
Z b2
fX1 ,X3 (x1 , x3 ) = fX1 ,X2 ,X3 (x1 , x2 , x3 ) dx2 . (116)
a2
216
And finally, pdf’s of conditional distributions are obtained through dividing
joint pdfs by marginal pdfs. Such conditional pdfs are used to calculate that some
RVs fall into a certain multidimensional interval given that some other RVs take
specific values. We only inspect a simple case analog to (109) where we want to
calculate the probability that X falls into a range [a, b] given that Y is known to
be c, that is, we want to evaluate the probability P (X ∈ [a, b] | Y = c), using pdfs.
We can obtain this probability from the joint pdf fX,Y and the marginal pdf fY
by Rb
fX,Y (x, c) dx
P (X ∈ [a, b] | Y = c) = a . (117)
fY (c)
Rb
The r.h.s. expression a fX,Y (x, c) dx / fY (c) is a function of x, parametrized
by c. This function is a pdf, denoted by fX | Y =c , and defined by
fX,Y (x, c)
fX | Y =c (x) = . (118)
fY (c)
Note that (117) and (118) are only defined if fY (c) > 0.
Let me illustrate this with a concrete example. An electronics engineer is test-
ing a device which transforms voltages V into currents I. In order to empirically
measure the behavior of this device (an electronics engineer would say, in order
to “characterize” the device), the engineer carries out a sequence of measurement
trials where he first sets the input voltage V to a specific value, say V = 0.0. Then
he (or she) measures the resulting current many times, in order to get an idea of
the stochastic spread of the current. In mathematical terms, the engineer wants to
get an idea of the pdf fI | V =0.0 . The engineer then carries on, setting the voltage
to other values c1 , c2 , ..., measuring resulting currents in each case, and getting
ideas of the conditional pdfs fI | V =ci . For understanding the characteristics of this
device, the engineer needs to know all of these pdfs.
Conditional distributions arise whenever cause-effect relationships are being
modeled. The conditioning variables are causes, the conditioned variables describe
effects. In experimental and empirical research, the causes are under the control of
an experimenter and can (and have to) be set to specific values in order to assess
the statistics of the effects – which are not under the control of the experimenter.
In ML pattern classification scenarios, the “causes” are the input patterns and
the “effects” are the (stochastically distributed) class label assignments. Since
research in the natural sciences is very much focussed on determining Nature’s
cause-effect workings, and 90% of the applications in machine learning concern
pattern classification (my estimate), it is obvious that conditional distributions lie
at the very heart of scientific (and engineering) modeling and data analysis.
In this appendix (and in the lecture) I consider only two ways of representing
probability distributions: discrete ones by finite probability tables or probability
tables; continuous ones by pdfs. These are the most elementary formats of repre-
senting probability distributions. There are many others which ML experts readily
217
command on. This large and varied universe of concrete representations of prob-
ability distributions is tied together by an abstract mathematical theory of the
probability distributions themselves, independent of particular representations.
This theory is called probability theory. It is not an easy theory and we don’t at-
tempt an introduction to it. If you are mathematically minded, then you can get an
introduction to probability theory in my graduate lecture notes “Principles of Sta-
tistical Modeling” ([Link] At
this point I only highlight two core facts from probability theory:
argmax φ(a)
a
is that d ∈ D for which φ(d) is maximal among all values of φ on D. If there are
several arguments a for which φ gives the same maximal value, – that is, φ does
not have a unique maximum –, or if φ has no maximum at all, then the argmax
is undefined.
219
spaces are considered, as well as distributions that have no pdf — but we will
focus on this basic scenario of real-valued RVs with pdfs).
The expectation of a RV X with sample space RK and pdf fX is defined as
Z
E[X] = x fX (x) dx, (120)
RK
X
N
mean({x1 , . . . , xN }) = 1/N xi ,
i=1
but this sample mean is NOT the expectation of X. If we would have used another
random sample, we would most likely have obtained another sample mean. In
contrast, the expectation E[X] of X is defined not on the basis of a finite, random
sample of X, but it is defined by averaging over the true underlying distribution.
Since in practice we will not have access to the true pdf fX , the expectation
of a RV X cannot usually be determined in full precision. The best one can do is
to estimate it from observed sample data. The sample mean is an estimator for
the expectation of a numerical RV X. Marking estimated quantities by a “hat”
accent, we may write
XN
Ê[X] = 1/N xi .
i=1
220
which in terms of the pdf fX̄ of X̄ can be written as
Z
2
σ (X) = x2 fX̄ (x) dx.
R
but in fact this estimator is not the best possible – on average (across different
samples) it underestimates the true variance. If one wishes to have an estimator
that is unbiased, that is, which on average across different samples gives the correct
variance, one must use
!2
X
N XN
σ̂ 2 ({x1 , . . . , xN }) = 1/(N − 1) xi − 1/N xj
i=1 j=1
Finally, let us inspect the correlation of two scalar RVs X, Y . Here we have to
be careful because this term is used differently in different fields. In statistics, the
correlation is defined as
Cov(X, Y )
Corr(X, Y ) = . (123)
σ(X) σ(Y )
221
It is easy to show that −1 ≤ Corr(X, Y ) ≤ 1. The correlation in the understanding
of statistics can be regarded as a normalized covariance. It has a value of 1 if X
and Y are identical up to some positive scaling factor, it has a value of −1 if X
and Y are identical up to some negative scaling factor. When Corr(X, Y ) = 0, X
and Y are said to be uncorrelated.
The quantity Corr(X, Y ) is also referred to as (population) Pearson’s correla-
tion coefficient, and is often denoted by the greek letter ϱ(X, Y ) = Corr(X, Y ).
In the signal processing literature (for instance in my favorite textbook Farhang-
Boroujeny (1998)), the term “correlation” is sometimes used in quite a different
way, denoting the quantity
E[X Y ],
that is, simply the expectation of the product of the uncentered RVs X and Y .
Just be careful when you read terms like “correlation” or “cross-correlation” or
“cross-correlation matrix” and make sure that your understanding of the term is
the same as the respective author’s.
There are some basic rules for doing calculations with expectations and covari-
ance which one should know:
2. Expectation is idempotent:
E[E[X]] = E[X].
3.
Cov(X, Y ) = E[X Y ] − E[X] E[Y ].
222
References
R. H. Abraham and C. D. Shaw. Dynamics: The Geometry of Behavior.
Addison-Wesley, Redwood City, 1992. e-book at [Link]
com/[Link], 30 USD.
D. H. Ackley, G. E. Hinton, and T. J. Sejnowski. A learning algorithm for Boltz-
mann machines. Cognitive Science, 9:147–169, 1985.
D. J. Amit. Modeling Brain Function: the World of Attractor Neural Networks.
Cambridge Univ. Press, NY, 1989.
D. J. Amit, H. Gutfreund, and H. Sompolinsky. Storing infinite numbers of pat-
terns in a spin-glass model of neural networks. Phys. Rev. Lett., 55(14):1530 –
1533, 1985.
Atomic Heritage Foundation authors. Computing and the Man-
hattan project. [Link]
computing-and-manhattan-project, 2014. Accessed: 2019.
A. Baddeley. Working memory: looking back and looking forward. Nature Reviews:
Neuroscience, 4(10):829–839, 2003.
F. C. Bartlett. Remembering: a study in experimental and social psychology.
Cambridge University Press, 1932.
K. Basye, Th. Dean, and L. P. Kaelbling. Learning dynamics: System identi-
fication for perceptually challenged agents. Artificial Intelligence, 72:139–171,
1995.
Y. Bengio and Y. LeCun. Scaling learning algorithms towards AI. In Bottou L.,
Chapelle O., DeCoste D., and Weston J., editors, Large-Scale Kernel Machines.
MIT Press, 2007.
C. M. Bishop. Pattern Recognition and Machine Learning. Springer Verlag, 2006.
L. Buesing, J. Bill, B. Nessler, and W. Maass. Neural dynamics as sampling:
A model for stochastic computation in recurrent networks of spiking neurons.
PLoS Comp. Biol., 7(11):e1002211, 2011.
T. P. Castello-Waldow, G. Weston, A. Chenani, Y. Loewenstein, A. Chen, and
Al Attardo. Stability of excitatory structural connectivity predicts the probabil-
ity of CA1 pyramidal neurons to become engram neurons. biorxiv manuscript,
[Link] 2019.
G. Chlebus, A. Schenk, J.H. Moltz, B. van Ginneken, H. K. Hahn, and H. Meine.
Automatic liver tumor segmentation in CT with fully convolutional neural net-
works and object-based postprocessing. Scientific Reports, 8:article number
15497, 2018. doi: 10.1038/s41598-018-33860-7.
223
K. Cho, B. van Merrienboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk,
and Y. Bengio. Learning phrase representations using RNN encoder-decoder
for statistical machine translation. arxiv manuscript [Link]
1406.1078, 2014.
L. O. Chua and L. Yang. Cellular neural networks: Theory. IEEE Transactions
on circuits and systems, 35(10):1257–1272, 1988.
A. Clark. Whatever next? Predictive brains, situated agents, and the future of
cognitive science. Behavioral and Brain Sciences, 36(3):1–86, 2013.
M. Cook. Universality in elementary cellular automata. Complex Systems, 15:
1–40, 2004.
S. Cook. The P versus NP problem. Official problem description
of the fourth millenium problem, Clay Mathematics Institute, 2000.
[Link]
J. C. Coulombe, M. C. A. York, and J. Sylvestre. Computing with networks
of nonlinear mechanical oscillators. PLOS ONE, 12(6), 2017. URL Https:
//[Link]/10.1371/[Link].0178663.
S. Das. CNN architectures: LeNet, AlexNet, VGG, GoogLeNet, ResNet
and more. Online article [Link]
cnns-architectures-lenet-alexnet-vgg-googlenet-resnet-and-more-666091488df5,
2017. Accessed: 14 Mar 2020.
M. Dascalu. Cellular automata hardware implementations – an overview. Science
and Technology, 19(4):360–368, 2016.
A. Deisenroth, A. Faisal, and C. S. Ong. Mathematics for Machine Learning.
Cambridge University Press, 2019. Free online copy at [Link]
[Link]/.
D. Deitch, A. Rubin, and Y. Ziv. Representational drift in the mouse visual cortex.
Current Biology, 31(19):4327–4339, 2021.
S. Demyanov. Regularization Methods for Neural Networks and Related Models.
Phd thesis, Dept of Computing and Information Systems, Univ. of Melbourne,
2015.
P. F. Dominey. Complex sensory-motor sequence learning based on recurrent state
representation and reinforcement learning. Biological Cybernetics, 73:265–274,
1995.
K. Doya. Bifurcations in the learning of recurrent neural networks. In Proceedings
of 1992 IEEE Int. Symp. on Circuits and Systems, Vol. 6, pages 2777–2780,
1992.
224
V. Evans, B. K. Bergen, and J. Zinken, editors. The Cognitive
Linguistics Reader. Equinox Publishing Ltd, 2007. Open ac-
cess copy at [Link]
LINGUISTICS-VYVYAN_EVANS20190529-[Link].
225
I. Goodfellow, Y. Bengio, and A. Courville. Deep Learning. MIT Press, 2016.
Open access version at [Link]
D. O. Hebb. The Organization of Behavior. New York: Wiley & Sons, 1949.
226
A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane cur-
rent and its application to conduction and excitation in nerve. The Journal of
Physiology, 117(4):500–544, 1952.
H. Jaeger. The ”echo state” approach to analysing and training recurrent neural
networks. GMD Report 148, GMD - German National Research Institute for
Computer Science, 2001. URL [Link]
227
D. P. Kingma, D. J. Rezende, S. Mohamed, and M. Welling. Semi-
supervised learning with deep generative models. In Proc. NIPS 2014, 2014.
arXiv:1406.5298v2.
H. Korn and P. Faure. Is there chaos in the brain? II. experimental evidence and
related models. Comptes Rendus de l’Académie des Sciences-Series III-Sciences
de la Vie, 326(9), 2003.
G. Lakoff. Women, fire, and dangerous things: What categories reveal about the
mind. University of Chicago Press, 1987.
228
W. Maass, T. Natschläger, and H. Markram. Real-time computing without stable
states: A new framework for neural computation based on perturbations. Neu-
ral Computation, 14(11):2531–2560, 2002a. URL [Link]
papers/[Link].
M. Minsky. The Society of Mind. New York: Simon & Schuster, 1986.
229
M. Olazaran. A sociological study of the official history of the perceptrons con-
troversy. Social Studies of Science, 26(3):611–659, 1996.
M. A. Savi. Nonlinear dynamics and chaos. In V. Lopes Junior and et al, editors,
Dynamics of Smart Systems and Structures, pages 93–117. Springer Interna-
tional Publishing Switzerland, 2016.
230
J. Schmidhuber. Deep learning in neural networks: An overview. Neural Networks,
61:85–117, 2015. Preprint: arXiv:1404.7828.
231