TL;DR
This paper provides comprehensive lecture notes on neural networks, covering theoretical foundations, supervised and unsupervised learning, and recent advances in deep learning architectures for various applications.
Contribution
It offers an organized overview of neural network principles, including Hopfield networks, supervised learning, and unsupervised techniques, with insights into recent developments like convolutional and recurrent networks.
Findings
Explains fundamental learning rules for neural networks.
Highlights recent successes of convolutional and recurrent networks.
Provides a unified view of supervised and unsupervised learning methods.
Abstract
These are lecture notes for a course on machine learning with neural networks for scientists and engineers that I have given at Gothenburg University and Chalmers Technical University in Gothenburg, Sweden. The material is organised into three parts: Hopfield networks, supervised learning of labeled data, and learning algorithms for unlabeled data sets. Part I introduces stochastic recurrent networks: Hopfield networks and Boltzmann machines. The analysis of their learning rules sets the scene for the later parts. Part II describes supervised learning with multilayer perceptrons and convolutional neural networks. This part starts with a simple geometrical interpretation of the learning rule and leads to the recent successes of convolutional networks in object recognition, recurrent networks in language processing, and reservoir computers in time-series analysis. Part III explains what…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 18
Figure 19
Figure 20
Figure 21
Figure 22
Figure 23
Figure 24
Figure 25
Figure 26
Figure 27
Figure 28
Figure 29
Figure 30
Figure 1
Figure 32
Figure 33
Figure 34
Figure 35
Figure 36
Figure 37
Figure 38
Figure 39
Figure 40| [5976, 1543, 1319, 1120, 42] |
| [4513, 2823, 2057, 607] |
| [4513, 1543, 1319, 1120, 607, 514, 342, 42] |
| [8479, 4868, 3696, 2646, 169, 142] |
| [11968, 5026, 1081, 1050, 691, 184] |
| [8479, 4167, 2646, 1081, 881, 859, 701, 691, 184, 169, 142] |
| [9979, 9348, 8022, 4020, 2693, 1892, 1714, 1371, 510, 451] |
| [9492, 8453, 7749, 7365, 2292, 2180, 1023, 959, 278, 124, 85] |
| [7042, 5608, 5464, 4371, 3884, 3121, 1901, 1768, 1590, 959, 899, 707, 702, 510, 451, 412, |
| [ 278, 124, 124, 85] |
| output | target | classification | correct? | |
|---|---|---|---|---|
| [0.4, 0.5, 0.4] | [0,1,0] | versicolor | yes | |
| [0.4, 0.3, 0.5] | [0,0,1] | setosa | yes | |
| [0.6, 0.5, 0.4] | [1,0,0] | virginica | yes |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Code & Models
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Machine learning with neural networks
Bernhard Mehlig
(2021)
\documenttype\department
Department of Physics
Acknowledgements.
This textbook is based on lecture notes for the course Artificial Neural Networks that I have given at Gothenburg University and at Chalmers Technical University in Gothenburg, Sweden. When I prepared my lectures, my main source was Introduction to the theory of neural computation by Hertz, Krogh, and Palmer [1]. Other sources were Neural Networks: a comprehensive foundation by Haykin [2], Horner’s lecture notes [3] from Heidelberg, Deep learning by Goodfellow, Bengio & Courville [4], and the online book Neural Networks and Deep Learning by Nielsen [5]. I thank Martin Čejka for typesetting the first version of my hand-written lecture notes, Erik Werner and Hampus Linander for their help in preparing Chapter 8, Kristian Gustafsson for his detailed feedback on Chapter 11, Nihat Ay for his comments on Section 4.5, and Mats Granath for discussions about autoencoders. I would also like to thank Juan Diego Arango, Oleksandr Balabanov, Anshuman Dubey, Johan Fries, Phillip Gräfensteiner, Navid Mousavi, Marina Rafajlovic, Jan Schiffeler, Ludvig Storm, and Arvid Wenzel Wartenberg for implementing algorithms described in this book. Many Figures are based on their results. Oleksandr Balabanov, Anshuman Dubey, Jan Meibohm, and in particular Johan Fries and Marina Rafajlovic contributed exam questions that became exercises in this book. Finally, I would like to express my gratitude to Stellan Östlund, for his encouragement and criticism. Last but not least, a large number of colleagues and students – past and present – pointed out misprints and errors, and suggested improvements. I thank them all. The present version does not contain exercises (copyright owned by Cambridge University Press). The complete book is available from Cambridge University Press. The cover image shows an input pattern designed to maximise the output of neurons corresponding to one feature map in a given convolution layer of a deep convolutional neural network [129,130]. See also page 8.7. Image by Hampus Linander. Reproduced with permission.
\coverfigure
\begin{overpic}[scale={1.28}]{figures/cover_image/out_7_980.png} \end{overpic}
\covercaptionThe cover image shows an input pattern designed to maximise the output of neurons corresponding to one feature map in a given convolution layer of a deep convolutional neural network [129,130]. See also page 8.7. Image by Hampus Linander. Reproduced with permission.
\makecoverpage\maketitlepage\makeacknowledgementspage
Contents
Chapter 1 Introduction
The term neural networks historically refers to networks of neurons in the mammalian brain. Neurons are its fundamental units of computation, and they are connected together in networks to process data. This can be a very complex task. The dynamics of such neural networks in response to external stimuli is therefore often quite intricate. Inputs and outputs of each neuron vary as functions of time in the form of spike trains, but also the network itself changes over time: we learn and improve our data-processing capabilities by establishing new connections between neurons.
Neural-network algorithms for machine learning are inspired by the architecture and the dynamics of networks of neurons in the brain. The algorithms use highly idealised neuron models. Nevertheless, the fundamental principle is the same: artificial neural networks learn by changing the connections between their neurons. Such networks can perform a multitude of information-processing tasks.
Neural networks can for instance learn to recognise structures in a data set and, to some extent, generalise what they have learnt. A training set contains a list of input patterns together with a list of corresponding labels, or target values, that encode the properties of the input patterns the network is supposed to learn. Artificial neural networks can be trained to classify such data very accurately by adjusting the connection strengths between their neurons, and can learn to generalise the result to other data sets – provided that the new data is not too different from the training data. A prime example for a problem of this type is object recognition in images, for instance in the sequence of camera images taken by a self-driving car. Recent interest in machine learning with neural networks is driven in part by the success of neural networks in visual object recognition.
Another task at which neural networks excel is machine translation with dynamical, or recurrent networks. Such networks take sentences as inputs. As one feeds word after word, the network outputs the words in the translated sentence. Recurrent networks can be efficiently trained on large training sets of input sentences and their translations. Google translate works in this way. Recurrent networks have also been used with considerable success to predict chaotic dynamics. These are all examples of supervised learning, where the networks are trained to associate a certain target, or label, with each input.
Artificial neural networks are also good at analysing large sets of unlabeled, often high-dimensional data – where it may be difficult to determine a priori which questions are most relevant and rewarding to ask. Unsupervised-learning algorithms organise the unlabeled input data in different ways. They can for instance detect familiarity and similarity (clusters) of input patterns and other structures in the input data. Unsupervised-learning algorithms work well when there is redundancy in the input data, and they are particularly useful for large, high-dimensional data sets, where it may be a challenge to detect clusters or other data structures by inspection.
Many problems lie between these two extremes of supervised and unsupervised learning. Consider how an agent may learn to navigate a complex environment, in order to get from one location to another as quickly as possible, or expending as little energy as possible. The method of reinforcement learning allows the agent to do just that, by optimising its behaviour in response to environmental cues in the shape of penalties and rewards. In short, the agent learns to act in such a way that it receives positive feedback (reward) more often than a penalty.
The tools for machine learning with neural networks were developed long ago, most of them during the second half of the last century. In 1943, McCulloch and Pitts [6] analysed how networks of neurons can process information. Using an abstract model for a neuron, they demonstrated how such units can be coupled together to represent logical functions (Figure 1.1). Their analysis and conclusions are formulated using Carnap’s logical syntax [7], not in terms of algebraic equations as we are used to today. Nevertheless, their neuron model is essentially the binary threshold unit, closely related to the fundamental building block of most neural-network algorithms for machine learning to this date. In this book we therefore refer to this model as the McCulloch-Pitts neuron. The purpose of this early research on neural networks was to explain neuro-physiological mechanisms [8]. Perhaps the most significant advance was Hebb’s learning principle, describing how neural networks learn by strengthening connections between neurons that are active simultaneously. The principle is described in Hebb’s book The organization of behavior: A neuropsychological theory [9], published in 1949.
About ten years later, research in artificial neural networks had intensified, sparked by the influential work of Rosenblatt. In 1958 he formulated a learning rule for the McCulloch-Pitts neuron, related to Hebb’s rule, and demonstrated that the rule converges to the correct solution for all problems this model can solve [10]. He coined the term perceptron for layered networks of McCulloch-Pitts neurons, and showed that such neural networks could in principle solve tasks that a single McCulloch-Pitts neuron could not. However, there was no general learning rule for perceptrons at the time. The work of Minsky and Papert [11] emphasised the geometric aspects of learning. This allowed them to prove which kinds of problems perceptrons could solve, and which not. In 1969 they summarised these results in their elegant book Perceptrons. An Introduction to Computational Geometry.
Perceptrons are classifiers that output a label for each input pattern. A perceptron represents a mathematical function, an input-output mapping. A breakthrough in perceptron learning was the paper by Rumelhart et al. [12]. The authors demonstrated in 1986 that perceptrons can be trained by gradient descent. This means that the connection strengths between the neurons are iteratively changed in small steps, to eventually minimise the output error. For a single McCulloch-Pitts neuron, this gives essentially Hebb’s rule. The important point is that gradient descent allows to efficiently train perceptrons with many layers (backpropagation for multi-layer perceptrons). A second contribution of Rumelhart et al. is the idea to use local feature maps for object recognition with neural networks. The corresponding mathematical operation is a convolution. Therefore such architectures are now known as convolutional networks.
The work of Hopfield popularised an entirely different approach, also based on Hebb’s rule. In 1982, Hopfield analysed the properties of a dynamical, or recurrent, network that functions as a memory [13]: the dynamics is designed to find stored patterns by converging to a corresponding steady state. Such Hopfield networks were especially popular amongst physicists because there are close connections to the statistical physics of spin glasses that made it possible to derive a precise mathematical understanding of such artificial neural networks. Hopfield networks are the basis for important developments in computer science. More general recurrent networks, for example, are trained like perceptrons for language processing. Hopfield networks with hidden neurons, so-called Boltzmann machines [14] are generative models that allow to sample from a distribution the neural network learned. The training algorithm for Boltzmann machines with many hidden layers [15],
published in 1986, is one of the first algorithms for training networks with many layers, so-called deep networks.
An important problem in behavioural psychology is to understand how we learn from experience. One hypothesis is that desirable behaviours are reinforced by positive feedback. Around the same time as researchers began to analyse perceptrons, a different line of neural-network research emerged: to find learning rules that allow neural networks to learn by reinforcement. In many problems of this kind, the positive feedback or reward is not immediate but is received at some time in the future, as in a board game for example. Therefore it is necessary to understand how to estimate the future reward for a certain behaviour, and how to find strategies that optimise the future reward. Reinforcement learning [16] is designed for this purpose. In 1995, an early application of this method demonstrated how a neural network could learn to play the board game backgammon [17].
A related research field originated from the neuro-physiological question: how do we learn to map visual or sensory stimuli to spatio-temporal patterns of neural activity in the brain? In 1992 Kohonen described a self-organising map [18] that suggests how neurons might create meaningful geometric representations of inputs. At the same time, Kohonen’s algorithm is one of the first methods for non-linear dimensionality reduction for large data sets.
There are many connections between neural-network algorithms for machine learning and methods used in mathematical statistics, such as for instance Markov-chain Monte-Carlo algorithms and simulated-annealing methods. Certain unsupervised learning algorithms are related to principal-component analysis, others to clustering algorithms such as -means clustering. Supervised learning with deep networks is essentially regression analysis, trying to fit an input-output function to the training data. In other words this is just function fitting – and usually with a very large number of fitting parameters. Recent convolutional neural networks have millions of parameters. To determine so many parameters requires very large and accurate data sets. This makes it clear that neural networks are not a solution of everything. One of the difficult problems is to understand when machine learning with neural networks is called for, and when not. Therefore we need a detailed understanding of how the algorithms work, and in particular when and how they fail.
There were some early successes of machine learning with neural networks, but these methods were not widely used in the last century. During the past decade, by contrast, machine learning with neural networks has become increasingly successful and popular. For many applications, neural-network based algorithms are now regarded as the method of choice, for example for predicting how proteins fold [19]. What caused this paradigm shift? After all, the methods are essentially those developed forty or more years ago. A reason for the new success is perhaps that industry, in acute need of machine-learning algorithms for large data sets, has invested time and effort into generating larger and more accurate training sets than previously available. Computer hardware has improved too, so that networks with many layers containing many neurons can now be efficiently trained, making the recent progress possible.
This book is based on notes for lectures on artificial neural networks I have given for more than 15 years at Gothenburg University and Chalmers University of Technology in Gothenburg, Sweden. When I prepared these lectures, my primary source was Introduction to the theory of neural computation by Hertz, Krogh, and Palmer [1]. The material is organised into three parts: Hopfield networks, supervised learning of labeled data, and learning for unlabeled data sets (unsupervised and reinforcement learning). One reason for starting with Hopfield networks is that there is an elegant mathematical theory that describes how these neural networks learn, making it possible to understand their strengths and weaknesses from first principles. This is not the case for most of the other algorithms discussed in this book. The analysis of Hopfield networks sets the scene for the later parts of the book. Part II describes supervised learning with multilayer perceptrons and convolutional neural networks, starting from the simple geometrical picture emphasised by Minsky and Papert, and leading to the recent successes of convolutional networks in object recognition, and recurrent networks in language processing. Part III explains what neural networks can learn about data that is not labeled, with particular emphasis on reinforcement learning. The overall goal is to explain the fundamental principles that allow neural networks to learn, emphasising ideas and concepts that are common to all three parts.
1.1 Neural networks
Different regions in the mammalian brain perform different tasks. The cerebral cortex is the outer part of the mammalian brain, one of its largest and most developed segments. We can think of the cerebral cortex as a thin sheet (about 2 to 5 mm thick) that folds upon itself to form a layered structure with a large surface area that can accommodate large numbers of nerve cells, neurons. The human cerebral cortex contains about neurons. They are linked together by nerve strands (axons) that branch and end in synapses. These synapses are the connections to other neurons. The synapses connect to dendrites, branches extending from the neural cell body that are designed to collect input from other neurons in the form of electrical signals. A neuron in the human brain may have thousands of synaptic connections with other neurons. The resulting network of connected neurons in the cerebral cortex is responsible for processing of visual, audio, and sensory data.
Figure 1.2 shows neurons in the cerebral cortex. This drawing was made by Santiago Ramón y Cajal more than 100 years ago. By microscope he studied the structure of neural networks in the brain and documented his observations by ink-on-paper drawings like the one reproduced in Figure 1.2. One can distinguish the cell bodies of the neural cells, their axons (f), and their dendrites. The axons of some neurons connect to the dendrites of other neurons, forming a neural network (see Ref. [21] for a slightly more detailed description of this drawing).
A schematic image of a neuron is drawn in Figure 1.3. Information is processed from left to right. On the left are the dendrites that receive signals and connect to the cell body of the neuron where the signal is processed. The right part of the Figure shows the axon, through which the output is sent to other neurons. The axon connects to their dendrites via synapses.
Information is transmitted as an electrical signal. Figure 1.4 shows an example of the time series of the electric potential for a pyramidal neuron in fish [22]. The time series consists of an intermittent series of electrical-potential spikes. Quiescent periods without spikes occur when the neuron is inactive, during spike-rich periods we say that the neuron is active.
1.2 McCulloch-Pitts neurons
In artificial neural networks, the ways in which information is processed and signals are transferred are highly idealised. McCulloch and Pitts [6] modelled the neuron, the computational unit of the neural network, as a binary threshold unit.
It has only two possible outputs, or states: active or inactive. To compute the output, the unit sums the weighted inputs. If the sum exceeds a given threshold, the state of the neuron is said to be active, otherwise inactive. A slightly more general model than the original one is illustrated in Figure 1.5. The model performs repeated computations in discrete time steps . The state of neuron number at time step is denoted by
[TABLE]
Given the states , neuron number computes
[TABLE]
Here is the signum function (Figure 1.5):
[TABLE]
The argument of the signum function,
[TABLE]
is called the local field.
We see that the neuron performs a weighted average of the inputs . The parameters are called weights. Here the first index, , refers to the neuron that does the computation, and labels the neurons that connect to neuron . In general weights between different pairs of neurons assume different numerical values, reflecting different strengths of the synaptic couplings. Weights can be positive or negative, and we say that there is no connection when .
In this book we refer to the model described in Figure 1.5 as the McCulloch-Pitts neuron, although their original model had additional constraints on the weights. The threshold111In the deep-learning literature [4], the thresholds are called biases, defined as the negative of , with a plus sign in Equation (1.4). In this book we use the convention (1.4), with the minus sign. for neuron is denoted by .
Finally, note that the computation (1.2) is performed for all neurons in parallel, given the states at time step . The outputs are the inputs to all neurons at the next time step. Therefore the outputs have the time argument . These steps are repeated many times, resulting in time series of the activity levels of all neurons in the network.
1.3 Activation functions
The McCulloch-Pitts model approximates the patterns of spiking activity in Figure 1.4 in terms of two states, and , representing the inactive and active periods shown in the Figure. For many computation tasks this is sufficient, and for our purposes it does not matter that the dynamics of electrical signals in the cortex is quite different in detail. The aim after all is not to model the neural dynamics in the brain, but to construct computation models inspired by this dynamics.
It will become apparent later that the simplest model described above must be generalised somewhat for certain tasks and questions. For example, the jump in the signum function at may cause large fluctuations in the activity levels of a network of neurons, caused by infinitesimal changes of the local fields across . To dampen this effect, one allows the neuron to respond continuously to its inputs, replacing Eq. (1.2) by
[TABLE]
Here is a continuous activation function. It could just be a linear function, . But we shall see that many tasks require non-linear activation functions, such as (Figure 1.7). When the activation function is continuous, the neuron states assume continuous values too, not just the discrete values and given in Equation (1.1).
Alternatively one may use a piecewise linear activation function (Figure 1.8). This is motivated in part by the response curve of the leaky integrate-and-fire neuron, a model for the relation between the electrical current through the cell membrane into the neuron cell and the membrane potential
. The simplest model for the dynamics of the membrane potential represents the neuron as a capacitor. In the leaky integrate-and-fire neuron, leakage is modelled by adding a resistor in parallel with the capacitor , so that
[TABLE]
For a constant current, the membrane potential grows from zero as a function of time, , where is the time constant of the model. One says that the neuron produces a spike when the membrane potential exceeds a critical value, . Immediately after, the membrane potential is set to zero (and begins to grow again). In this model, the firing rate is thus given by , where is the solution of . It follows that the firing rate exhibits a threshold behaviour. In other words, the system works like a rectifier:
[TABLE]
This response curve is illustrated in Figure 1.8 (a). The main point is that there is a threshold below which the response is strictly zero. The response function looks qualitatively like the piecewise linear function
[TABLE]
shown in panel (b). Neurons with this activation function are called rectified linear units, and the activation function (1.8) is called the ReLU function.
1.4 Asynchronous updates
Equations (1.2) and (1.5) are called synchronous update rules,
because all neurons are updated in parallel: at time step all inputs are stored. Then all neurons are simultaneously updated using the stored inputs. An alternative is to update only a single neuron per iteration, for example the one with index :
[TABLE]
This is called an asynchronous update rule. Different schemes for choosing neurons are used in asynchronous updating. One possibility is to arrange the neurons into a two-dimensional array and to update them one by one, in a certain order. In the typewriter scheme, for example, one updates the neurons in the top row of the array first, from left to right, then the second row from left to right, and so forth. A second possibility is to choose randomly which neuron to update.
If there are neurons, then one synchronous step corresponds to asynchronous steps, on average. This difference in time scales is not the only difference between synchronous and asynchronous updating. The asynchronous dynamics can be shown to converge to a definite state in certain cases, while the synchronous dynamics may fail to do so, resulting in periodic cycles that persist forever.
1.5 Summary
Artificial neural networks use a highly idealised model for the fundamental computation unit: the McCulloch-Pitts neuron (Figure 1.5) is a binary threshold unit, very similar to the model introduced originally by McCulloch and Pitts [6]. The units are linked together by weights , and each unit computes a weighted average of its inputs. The network performs these computations in sequence. Most neural-network algorithms are built using the model described in this Chapter.
1.6 Further reading
Two accounts of the history of artificial neural networks are especially recommended. First, the early history of the field is summarised in the Introduction to the second edition of Perceptrons. An introduction to computational geometry by Minsky and Papert [11], which came out in 1988. This book also contains a concise bibliography with comments by Minsky and Papert. Second, in a short note, Kanal [23] reviews the work of Rosenblatt and puts it into context.
Part I Hopfield networks
Hopfield networks [24, 13] are artificial neural networks that can recognise or reconstruct images, for instance the binary images of digits shown in Figure 2.1. The images are stored in the network by choosing its weights according to Hebb’s rule [9]. One feeds a distorted digit (Figure 2.2) to the network by assigning the initial states of its neurons to the bits of the distorted digit. The idea is that the neural-network dynamics converges to the closest stored digit. In this way the network can recognise the input as a distorted version of the correct pattern, it can retrieve the correct digit. Hopfield networks recognise patterns with many bits quite reliably, and in the past such networks were used to perform pattern-recognition tasks. Today there are more efficient algorithms for this purpose (Chapter 8).
Nevertheless, Hopfield networks exemplify fundamental principles of machine learning with neural networks. For a start, most neural-network algorithms discussed in this book are built from similar building blocks and use learning rules related to Hebb’s rule. Moreover, Hopfield networks are examples of recurrent networks, their neurons are updated following a dynamical rule. Widely used algorithms for machine translation and time-series prediction are based on this principle.
Furthermore, restricted Boltzmann machines are closely related to Hopfield networks. These machines use hidden neurons to learn distributions of input patterns.
This makes it possible to generate image textures and to complete partially obscured images [25]. Generalisations of these machines, deep-belief networks, are examples of the first deep network architectures for machine learning. Restricted Boltzmann machines have been developed into more efficient generative models, Helmholtz machines, to sample new patterns similar to those in a given data distribution. The training algorithm for recent generative models, variational autoencoders, is based on the same principles as the learning rule for Helmholtz machines.
The dynamics of Hopfield networks is closely related to stochastic Markov-chain Monte-Carlo algorithms used in a wide range of problems in the natural sciences. Finally,Hopfield networks highlight the role of stochasticity in neural-network dynamics. A certain degree of noise, not too much, can substantially improve the performance of the Hopfield network. In engineering problems it is usually better to avoid stochasticity, when it is due to multiplicative or additive noise that diminishes the performance of the system. In neural-network dynamics, by contrast, stochasticity is often helpful, because it allows to explore a wider range of configurations or actions and thus helps the dynamics to converge to a better solution. In general it is challenging to analyse the stochastic dynamics of neural networks. But for the Hopfield network much is known. The reason is that Hopfield networks are closely related to stochastic systems studied in statistical physics, so-called spin glasses [26, 27].
Chapter 2 Deterministic Hopfield networks
2.1 Pattern recognition
As an example for a pattern-recognition task, consider images (patterns), for instance the digits shown in Figure 2.1. The different patterns are indexed by . Here is the number of patterns ( in Figure 2.1). The bits of pattern are denoted by . The index labels the different bits of a given pattern, and is the number of bits per pattern ( in Figure 2.1). The bits are binary: they can take only the values and . To determine the generic properties of the algorithm, one often turns to random patterns where each bit is chosen randomly, taking either value with probability equal to . It is convenient to gather the bits of a pattern in a column vector like this
[TABLE]
In this book, vectors are written in bold math font, as in Equation (2.1).
The task of the neural network is to recognise distorted patterns, to determine for instance that the pattern on the right in Figure 2.2 is a perturbed version of the digit on the left in this Figure. To this end, one stores patterns in the network and presents it with a distorted version of one of these patterns, and asks the network to find the stored pattern that is most similar to the distorted one.
The formulation of the problem makes it necessary to define how similar a distorted pattern is to any of the stored patterns, \mbox{\boldmathx}^{(\nu)} say. One possibility is to quantify the distance between the patterns and \mbox{\boldmathx}^{(\nu)} in terms of the mean squared error:
[TABLE]
The prefactor is chosen so that the distance equals the fraction of bits by which two -patterns differ. Note that the distance (2.2) does not refer to distortions by translations, rotations, or shearing. An improved measure of distance might take the minimum distance between the patterns subject to all possible translations, rotations, and so forth.
2.2 Hopfield networks and Hebb’s rule
Hopfield networks [24, 13] are networks of McCulloch-Pitts neurons designed to solve the pattern-recognition task described in the previous Section. The bits of the patterns to be recognised are encoded in a particular choice of weights called Hebb’s rule, as explained in the following.
All possible states of the McCulloch-Pitts neurons in the network,
[TABLE]
form the configuration space of the network. The components of the states are updated either with the synchronous rule (1.2):
[TABLE]
or with the asynchronous rule
[TABLE]
To recognise a distorted pattern, one feeds its bits into the network, by assigning the initial states of the neurons to the pattern bits,
[TABLE]
The idea is to choose a set of weights so that the network dynamics (2.4) or (2.5) converges to the correct stored pattern. The choice of weights must depend on all patterns, \mbox{\boldmathx}^{(1)},\ldots,\mbox{\boldmathx}^{(p)}. We say that these patterns are stored in the network by assigning appropriate weights. For example, if is a distorted version of \mbox{\boldmathx}^{(\nu)} (Figure 2.2), then we want the network to converge to this pattern:
[TABLE]
Equation (2.7) means that the network corrects the errors in the distorted pattern . If this works, the stored pattern \mbox{\boldmathx}^{(\nu)} is said to be an attractor of the network dynamics.
But convergence is not guaranteed. If the initial distortion is too large, the network may converge to another pattern. It may converge to some other state that bears no or little relation to the stored patterns, or it may not converge at all. The region around pattern \mbox{\boldmathx}^{(\nu)} in which all patterns converge to \mbox{\boldmathx}^{(\nu)} is called the region of attraction of \mbox{\boldmathx}^{(\nu)}. The size of this region depends in an intricate way upon the ensemble of stored patterns, and there is no general convergence proof.
Therefore we ask a simpler question first: if one feeds one of the undistorted patterns, for instance \mbox{\boldmathx}^{(\nu)}, does the network recognise it as one of the stored, undistorted ones? The network should not make any changes to \mbox{\boldmathx}^{(\nu)} because all bits are correct:
[TABLE]
How can we choose weights and thresholds to ensure that Equation (2.8) holds? Let us consider a simple case first, where there is only one pattern, \mbox{\boldmathx}^{(1)}, to recognise. In this case a suitable choice of weights and thresholds is
[TABLE]
This learning rule reminds of a relation between activation patterns of neurons and their coupling, postulated by Hebb [9] more than 70 years ago:
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.
This is a mechanism for learning through establishing connections: the coupling between neurons tends to increase if they are active at the same time. Equation (2.25) expresses an analogous principle. Together with Equation (2.7) it says that the coupling between two neurons is positive if they are both active (); if their states differ, the coupling is negative. Therefore the rule (2.25) is called Hebb’s rule. Hopfield networks are networks of of McCulloch-Pitts neurons with weights determined by Hebb’s rule.
Does a Hopfield network recognise the pattern \mbox{\boldmathx}^{(1)} stored in this way? To check that the rule (2.9) does the trick, we feed the pattern to the network by assigning , and evaluate Equation (2.4):
[TABLE]
The second equality follows because can only take the values , so that . Using , we obtain
[TABLE]
Comparing Equation (2.11) with the update rule (2.4) shows that the bits of the pattern \mbox{\boldmathx}^{(1)} remain unchanged under the update, as required by Equation (2.8). In other words, the network recognises the pattern as a stored one, so Equation (2.9) does what we asked for. Note that we obtain the same result if we leave out the factor of in Equation (2.9).
But does the network correct errors? Is the pattern \mbox{\boldmathx}^{(1)} an attractor [Eq. (2.7)]? This question cannot be answered in general. Yet in practice Hopfield networks work often quite well. It is a fundamental insight that neural networks may perform well although no proof exists that their dynamics converges to the correct solution.
To illustrate the difficulties, consider an example, a Hopfield network with neurons (Figure 2.3). Store the pattern \mbox{\boldmathx}^{(1)} shown in Figure 2.3 by assigning the weights using Equation (2.9). Now consider a distorted pattern that has a non-zero distance to \mbox{\boldmathx}^{(1)} [Figure 2.4(a)]:
[TABLE]
To feed the pattern to the network, we set . Then we iterate the dynamics using the synchronous rule (2.4). Results for different distorted patterns are shown in Figure 2.4. We see that the first two distorted patterns converge to the stored pattern, cases (a) and (b). But the third distorted pattern does not [case (c)].
To understand this behaviour, we analyse the synchronous dynamics (2.4) using the weight matrix
[TABLE]
Here \mbox{\boldmathx}^{(1)\textsf{T}} denotes the transpose of the column vector \mbox{\boldmathx}^{(1)},
so that \mbox{\boldmathx}^{(1)\textsf{T}} is a row vector. The standard rules for matrix multiplication apply also to column and row vectors, they are just and matrices. So the product on the r.h.s. of Equation (2.13) is an matrix. In the following, matrices with elements or are written as , , and so forth. The product in Equation (2.13) is also referred to as an outer product. If we change the order of \mbox{\boldmathx}^{(1)} and \mbox{\boldmathx}^{(1)\textsf{T}} in the product, we obtain a number instead:
[TABLE]
The product (2.14) is called scalar product. It is also denoted by \mbox{\boldmatha}\cdot\mbox{\boldmathb}=\mbox{\boldmatha}^{\textsf{T}}\mbox{\boldmathb} and equals |\mbox{\boldmatha}||\mbox{\boldmathb}|\cos\varphi, where is the angle between the vectors and , and |\mbox{\boldmatha}| is the magnitude of . We use the same notation for multiplying a transposed vector with a matrix from the left: \mbox{\boldmathx}\cdot\mathbb{A}=\mbox{\boldmathx}^{\textsf{T}}\mathbb{A}.
An excellent source for those not familiar with these terms from Linear Algebra (Figure 2.5) is Chapter 6 of Mathematical methods of physics by Mathews and Walker [29].
Using Equation (2.14) we see that projects onto the vector \mbox{\boldmathx}^{(1)},
[TABLE]
In addition, it follows from Equation (2.14) that the matrix (2.13) is idempotent:
[TABLE]
Equations (2.15) and (2.16) together with imply that the network recognises the pattern \mbox{\boldmathx}^{(1)} as the stored one. This example illustrates the general proof, Equations (2.10) and (2.11).
Now consider the distorted pattern (a) in Figure 2.4. We feed this pattern to the network by assigning
[TABLE]
To compute one step in the synchronous dynamics (2.4), we apply to \mbox{\boldmaths}(t=0). This is done in two steps, using the outer-product form (2.13) of the weight matrix. We first multiply \mbox{\boldmaths}(t=0) with \mbox{\boldmathx}^{(1)\textsf{T}} from the left
[TABLE]
and then we multiply this result with \mbox{\boldmathx}^{(1)}. This gives
[TABLE]
Taking the signum of the -th component of the vector \mathbb{W}\mbox{\boldmaths}(t=0) yields :
[TABLE]
We conclude that the state of the network converges to the stored pattern, in one synchronous update. Since is idempotent, the network stays there: the pattern \mbox{\boldmathx}^{(1)} is an attractor. Case (b) in Figure 2.4 works in a similar way.
Now look at case (c), where the network fails to converge to the stored pattern. We feed this pattern to the network by assigning \mbox{\boldmaths}(t=0)=[-1,1,-1,-1]^{\sf T}. For one iteration of the synchronous dynamics we evaluate
[TABLE]
It follows that
[TABLE]
Using the update rule (2.4) we find
[TABLE]
Equation (2.16) implies that
[TABLE]
Thus the network shown in Figure 2.3 has two attractors,
the pattern \mbox{\boldmathx}^{(1)} as well as the inverted pattern -\mbox{\boldmathx}^{(1)}. As we shall see in Section 2.5, this is a general property of Hopfield networks: if \mbox{\boldmathx}^{(1)} is an attractor, then the pattern -\mbox{\boldmathx}^{(1)} is an attractor too. In the next Section we discuss what happens when more than one patterns are stored in the Hopfield network.
2.3 The cross-talk term
When there are more patterns than just one, we need to generalise Equation (2.9). One possibility is to simply sum Equation (2.9) over the stored patterns [13]:
[TABLE]
Equation (2.25) generalises Hebb’s rule to patterns. Because of the sum over , the relation to Hebb’s learning hypothesis is less clear, but we nevertheless refer to Equation (2.25) as Hebb’s rule. At any rate, we see that the weights are proportional to the second moments of the pattern bits. It is plausible that a neural network based the rule (2.25) can recognise properties of the patterns \mbox{\boldmathx}^{(\mu)} that are encoded in the two-point correlations of their bits.
Note that the weights are symmetric, . Also, note that the prefactor in Equation (2.25) is not important. It is chosen to simplify the large- analysis of the model (Chapter 3). An alternative version of Hebb’s rule [13, 2] sets the diagonal weights to zero:
[TABLE]
In this Section we use the form (2.26) of Hebb’s rule. If we assign the weights in this way, does the network recognise the stored patterns? The question is whether
[TABLE]
holds or not. To check this, we repeat the calculation described in the previous Section. As a first step we evaluate the local field
[TABLE]
Here we split the sum over the patterns into two contributions. The first term corresponds to , where refers to the pattern that was fed to the network, the one we want the network to recognise. The second term in Equation (2.28) contains the sum over the remaining patterns. Condition (2.27) is satisfied if the second term in (2.28) does not affect the sign of the r.h.s. of this Equation. This second term is called cross-talk term.
Whether adding the cross-talk term to \mbox{\boldmathx}^{(\nu)} affects or not, depends on the stored patterns. Since the cross-talk term contains a sum over , we expect that this term does not matter if is small enough. The fewer patterns we store in the network, the more likely it is that they are all recognised. Furthermore, by analogy with the example described in the previous Section, it is plausible that the stored patterns are then also attractors, so that slightly distorted patterns converge to the correct stored pattern.
For a more quantitative analysis of the effect of the cross-talk term we store patterns with random bits (random patterns). Different bits are assigned independently with equal probability:
[TABLE]
This means in particular that , and that the covariance of different bits vanishes:
[TABLE]
Here denotes an average over many realisations of random patterns, and is the Kronecker delta, equal to unity if but zero otherwise. In summary, different bits are uncorrelated.
Given an ensemble of random patterns, what is the probability that the cross-talk term changes ? In other words, what is the probability that the network produces a wrong bit in one asynchronous update, if all bits were initially correct? The magnitude of the cross-talk term does not matter when it has the same sign as . If it has a different sign, then the cross-talk term may matter. To determine when this is the case, one defines
[TABLE]
If then the cross-talk term has same sign as , so that the cross-talk term does not make a difference: adding this term does not change the sign of . If it does not matter either, only when . The network produces an error in bit of pattern if (we approximated in Equation (2.28), assuming that is large).
2.4 One-step error probability
The one-step error probability is defined as the probability that an error occurs in one attempt to update a bit, given that initially all bits were correct:
[TABLE]
Since patterns and bits are identically distributed, does not depend on or . Therefore does not carry any indices.
How does depend on the parameters of the problem, and ? When both and are large, we can use the central-limit theorem [29, 30] to answer this question. Since different bits/patterns are independent, we can think of as a sum of independent random numbers that take the values and with equal probabilities,
[TABLE]
There are terms in the sum on the r.h.s. because terms with are excluded, and also those with [Equation (2.26)]. If we use the rule (2.25) instead, then there is a correction to Equation (2.33) from the diagonal weights. For this
correction is small.
When and are large, the sum contains a large number of independently identically distributed random numbers with mean zero and variance unity. It follows from the central-limit theorem that is Gaussian distributed with mean zero and variance .
Since the central-limit theorem plays an important role in the analysis of neural-network algorithms, it is worth discussing this theorem in a little more detail. To begin with, note that the sum equals , where is the number of occurrences of in the sum. Choosing randomly to equal either or is called a Bernoulli trial [30], and the probability of drawing times and times is given by the binomial distribution [30]. In our case the probability of equals , so that
[TABLE]
Here denotes the number of ways in which occurrences of can be distributed over places.
We want to show that approaches a Gaussian distribution for large , with mean zero and with variance . Since the variance diverges as , it is convenient to use the variable . The central-limit theorem implies that is Gaussian with mean zero and unit variance in the limit of large . To prove that this is the case, we substitute into Equation (2.34) and take the limit of large using Stirling’s approximation
[TABLE]
Expanding to leading order in assuming that remains of order unity gives . Now one changes variables from to . This stretches local neighbourhoods to . Conservation of probability implies that . It follows that , so that . In other words, the distribution of is Gaussian with zero mean and unit variance, as we intended to show.
Returning to Equation (2.33), we conclude that is Gaussian distributed
[TABLE]
with zero mean, as illustrated in Figure 2.6, and with variance
[TABLE]
Here we used for large and .
Another way to compute this variance is to square Equation (2.33) and to average over random patterns:
[TABLE]
Here denotes the average over random realisations of . Since the random numbers are independent for different indices and because , we have that . So only the diagonal terms in the double sum contribute, summing to . This yields Equation (2.37).
To determine [Equation (2.32)] we must integrate the distribution of from to :
[TABLE]
Here is the error function defined as [31]
[TABLE]
Since erf increases monotonically as increases, we conclude that increases as increases, or as decreases. This is expected: it is more difficult for the network to distinguish stored patterns when there are more of them. On the other hand, it is easier to differentiate stored patterns if they have more bits. We also see that the one-step error probability depends on and only through the combination
[TABLE]
The parameter is called the storage capacity of the network. Figure 2.7 shows how depends on the storage capacity. For for example, the one-step error probability is slightly larger than 1%.
In the derivation of Equation (2.39) we assumed that the stored patterns are random with independent bits. Realistic patterns are not random. We nevertheless expect that describes the typical one-step error probability of the Hopfield network when and are large. However, it is straightforward to construct counter examples. Consider for example orthogonal patterns:
[TABLE]
For such patterns, the cross-talk term vanishes in the limit of large (Exercise 2.2), so that .
More importantly, the error probability defined in this Section refers only to the initial update, the first iteration. What happens in the next iteration, and after many iterations? Numerical experiments show that the error probability can be much higher in later iterations, because one error tends to increase the probability of making another error later on. So the estimate is only a lower bound for the probability of observing errors in the long run.
2.5 Energy function
Consider the long-time limit . Does the Hopfield dynamics converge, as required by Equation (2.7)? This is an important question in the analysis of neural-network algorithms, because an algorithm that does not converge to a meaningful solution is useless.
The standard way of analysing convergence of neural-network algorithms is to define an energy function H(\mbox{\boldmaths}) that has a minimum at the desired solution, \mbox{\boldmaths}=\mbox{\boldmathx}^{(1)} say. We monitor how the energy function changes as we iterate, and keep track of the smallest values of encountered, to find the minimum. If we store only one pattern, , then a suitable energy function is
[TABLE]
This function is minimal when \mbox{\boldmaths}=\mbox{\boldmathx}^{(1)}, i.e., when for all . It is customary to insert the factor , this does not change the fact that is minimal at \mbox{\boldmaths}=\mbox{\boldmathx}^{(1)}.
A crucial point is that the asynchronous McCulloch-Pitts dynamics (2.5) converges to a minimum of [13]. This follows from the fact that cannot increase under the update rule (2.5). To prove this important property, we begin by evaluating the expression on the r.h.s. of Equation (2.43):
[TABLE]
Using Hebb’s rule (2.9) we find that the energy function (2.43) becomes
[TABLE]
This function has the same form as the energy function (or Hamiltonian) for certain physical models of magnetic systems consisting of interacting spins [32], where the interaction energy between spins and is . Note that Hebb’s rule (2.9) yields symmetric weights, , and . Note also that setting the diagonal weights
to zero does not change the fact that is minimal at \mbox{\boldmaths}=\mbox{\boldmathx}^{(1)} because . The diagonal weights just give a constant contribution to , independent of .
The second step is to show that cannot increase under the asynchronous McCulloch-Pitts dynamics (2.5). In this case we say that the energy function is a Lyapunov function, or loss function. To demonstrate that the energy function is a Lyapunov function, choose a neuron and update it according to Equation (2.5). We denote the updated state of neuron by :
[TABLE]
All other neurons remain unchanged. There are two possibilities, either or . In the first case remains the same, . Here refers to the value of the energy function after the update (2.46). When , by contrast, the energy function changes by the amount
[TABLE]
The sum goes over all neurons that are connected to the neuron , the one to be updated in Equation (2.46). Now if the weights are symmetric, equals
[TABLE]
Since the sign of is that of and if , it follows that
[TABLE]
In other words,the value of must decrease when the state of neuron changes, . In summary,111The derivation outlined here did not use the specific form of Hebb’s rule (2.9), only that the weights are symmetric, and that . However, the derivation fails when . In this case it is still true that assumes a minimum at \mbox{\boldmaths}=\mbox{\boldmathx}^{(1)}, but can increase under the update rule, so that convergence is not guaranteed. We therefore require that the diagonal weights are not negative. either remains constant under the asynchronous McCulloch-Pitts dynamics (), or its value decreases (). Note that this does not hold for the synchronous dynamics (2.4), see Exercise 2.9. Since the energy function cannot increase under the asynchronous McCulloch-Pitts dynamics, it must converge to a minimum of the energy function. For the energy function (2.43) this implies that the dynamics must either converge to the stored pattern or to its inverse. Both are attractors.
We assumed the thresholds to vanish, but the proof also works when the thresholds are not zero, in this case for the energy function
[TABLE]
in conjunction with the update rule .
Up to now we considered only one stored pattern, . If we store more than one pattern [Hebb’s rule (2.25)], the proof that (2.45) cannot increase under the McCulloch-Pitts dynamics works in the same way because no particular form of the weights was assumed, only that they must be symmetric, and that the diagonal weights must not be negative.
Therefore it follows in this case too that the minima of the energy function must correspond to attractors, as illustrated schematically in Figure 2.8.
The configuration space of the network, corresponding to all possible choices of \mbox{\boldmaths}=[s_{1},\ldots s_{N}]^{\sf T}, is drawn as a single axis, the -axis. But when is large, the configuration space is really very high dimensional.
However, some stored patterns may not be attractors when . This follows from our analysis of the cross-talk term in Section 2.2 . If the cross-talk term causes errors for a certain stored pattern, then this pattern is not located at a minimum of the energy function. Another way to see this is to combine Equations (2.25) and (2.45) to give:
[TABLE]
While the energy function defined in Equation (2.43) has a minimum at \mbox{\boldmathx}^{(1)}, Equation (2.51) need not have a minimum at \mbox{\boldmathx}^{(1)} (or at any other stored pattern), because a maximal value of \,\big{(}\sum_{i=1}^{N}s_{i}x_{i}^{(1)}\big{)}^{2} may be compensated by terms stemming from other patterns. This happens rarely when is small (Section 2.2).
Conversely there may be minima that do not correspond to stored patterns. Such states are referred to as spurious states. The network may converge to spurious states. This is undesirable but it occurs even when there is only one stored pattern, as we saw in Section 2.2: the McCulloch-Pitts dynamics may converge to the inverted pattern. This follows also from Equation (2.51): if \mbox{\boldmaths}=\mbox{\boldmathx}^{(1)} is a local minimum of , then so is \mbox{\boldmaths}=-\mbox{\boldmathx}^{(1)}. This is a consequence of the invariance of under \mbox{\boldmaths}\to-\mbox{\boldmaths}. There are other types of spurious states besides inverted patterns. An example are mixed states, superpositions of an odd number of patterns [1]. For , for example, the bits of a mixed state read:
[TABLE]
The number of mixed states increases as increases. There are mixed states that are superpositions of out of patterns, for (Exercise 2.4). Mixed states such as (2.52) are sometimes recognised by the network (Exercise 2.5), therefore it may happen that the network converges to these states. Finally, there are spurious states that are not related in any way to the stored patterns . Such spin-glass states are discussed in detail in Refs. [33, 34, 27], and also by Hertz, Krogh and Palmer [1].
2.6 Summary
Hopfield networks are networks of McCulloch-Pitts neurons that recognise, or retrieve, patterns (Algorithm 1). Their layout is defined by connection strengths, or weights, chosen according to Hebb’s rule. The weights are symmetric, and the network is in general fully connected. Hebb’s rule ensures that stored patterns are recognised, at least most of the time if the number of patterns is not too large. Convergence of the McCulloch-Pitts dynamics is analysed in terms of an energy function, which cannot increase under this dynamics.
A single-step estimate for the error probability of the network dynamics was derived in Section 2.2. If one iterates several steps, the error probability is usually much larger, but it is difficult to evaluate in general. For stochastic Hopfield networks the steady-state error probability can be estimated more easily, because the dynamics converges to a steady state.
Chapter 3 Stochastic Hopfield networks
Two related problems became apparent in the previous Chapter. First, the Hopfield dynamics may get stuck in spurious minima. In fact, if there is a local minimum downhill from a given initial state, between this state and the correct attractor, then the dynamics arrests in the local minimum, so that the algorithm fails to converge to the correct attractor. Second, the energy function usually is a strongly varying function over a high-dimensional configuration space. Therefore it is difficult to predict the first local minimum encountered by the down-hill dynamics of the network.
Both problems are solved by introducing an element of stochasticity into the dynamics. This is a trick that works for many neural-network algorithms. In general, however, it is quite challenging to analyse the stochastic dynamics. For the Hopfield network, by contrast, much is known. The reason is that the stochastic Hopfield network is closely related to systems studied in statistical mechanics, so-called spin glasses. Like these systems – and many other physical systems – the stochastic Hopfield network exhibits an order-disorder transition. This transition becomes sharp in the limit of a large number of neurons. This has important consequences. Suppose that the network produces satisfactory results for a given number of patterns with a certain number of bits. If one tries to store just one more pattern, the network may fail to recognise anything. The goal of this Chapter is to explain why this occurs, and how it can be avoided.
3.1 Stochastic dynamics
The asynchronous update rule (2.5) is called deterministic, because a given set of states determines the outcome of the update of neuron . To introduce noise, one replaces the rule (2.5) by an asynchronous stochastic rule [36]:
[TABLE]
The function is plotted in Figure 3.1.
The parameter is the noise parameter. When is large, the noise level is small. As tends to infinity, the function approaches zero if is negative, and it tends to unity if is positive. So for , the stochastic update rule
(3.1) converges to the deterministic rule (2.5). In the opposite limit, when , the update probability simply equals . In this case is updated to or randomly, with equal probability. The dynamics does not depend upon the stored patterns contained in the local field .
The idea is to keep a small but finite noise level. Then the network dynamics is very similar to the deterministic Hopfield dynamics analysed in the previous Chapter. But the noise allows the system to escape spurious minima. However, since the dynamics is stochastic, we must rephrase the convergence criterion (2.7). This is discussed next.
3.2 Order parameters
If we feed one of the stored patterns, \mbox{\boldmathx}^{(1)} for example, then we want the stochastic dynamics to stay in the vicinity of \mbox{\boldmathx}^{(1)}. This can only work if the noise is weak enough, and even then it is not guaranteed. At time step , bit is correct if . All bits are correct when , otherwise the sum takes a value smaller than . One measures success by averaging over the asynchronous stochastic dynamics of the network from to , for given bits :
[TABLE]
If we feed pattern \mbox{\boldmathx}^{(1)} to the network, we have initially. We want that remains close to unity, so that the network recognises the pattern \mbox{\boldmathx}^{(1)}. In practice, the quantity settles into a steady state, where it fluctuates around a mean value with a definite distribution that becomes independent of the iteration number . If the network works well, the finite-time average converges to a value of order unity after an initial transient (Figure 3.2). The limiting value,
[TABLE]
is called the order parameter. Since there is noise, the order parameter is usually smaller than unity. The last equality in Equation (3.2b) defines the time average over the stochastic network dynamics.
Figure 3.2 also illustrates a subtlety. For finite values of , the order parameter depends upon the stored patterns. Different realisations \mbox{\boldmathx}^{(1)},\ldots,\mbox{\boldmathx}^{(p)} of random patterns yield different values of . In the limit of this problem does not occur, the order parameter becomes independent of the stored patterns. We say that the system is self averaging in this limit. To obtain a definite value for the order parameter for finite values of , one usually averages over different realisations of random patterns stored in the network (thick solid line in Figure 3.2). The dashed line in Figure 3.2 shows this average, .
The other components, for , are expected to be small. This is certainly true for random patterns with many independent bits. If , the individual terms in the sumer over in Equation (3.2b) cancel approximately upon summation, because the bits of the patterns \mbox{\boldmathx}^{(2)} to \mbox{\boldmathx}^{(p)} are independent from those of \mbox{\boldmathx}^{(1)}. In summary, if we feed pattern \mbox{\boldmathx}^{(1)} and if the network works well, we expect in the limit of large :
[TABLE]
Whether this is the case or not depends on the values of , , and . In the next Sections we determine how depends on these parameters.
3.3 Mean-field theory
The order parameter is defined as an average over the stochastic dynamics of the network in its steady state (Figure 3.2). It is a challenging task to compute this average because all neurons interact with each other in a nonlinear fashion. Consider neuron number . The fate of is determined by its local field , through Equation (3.1). The difficulty is that the local field in turn depends on the states of all other neurons in the network:111We set the thresholds to zero, as assumed in Hebb’s rule (2.26).
[TABLE]
When is large, we may assume that remains essentially constant in the steady state, independent of , because fluctuations of average out when summing over :
[TABLE]
Since is given by a sum over many random numbers, we appeal to the central-limit theorem and argue that the fluctuations of are of order . Since , we ignore the fluctuations in the limit of large and write
[TABLE]
using Hebb’s rule (2.26) for given patterns \mbox{\boldmathx}^{(\mu)}. The time-averaged local field is called the mean field. Theories that neglect the fluctuations in Equation (3.5) are called mean-field theories. They require a self-consistent solution, because the average on the r.h.s. of Equation (3.6) depends on the mean field. Using the stochastic update rule (3.1) we find:
[TABLE]
Equations (3.6) and (3.7) yield a set of non-linear self-consistent equations for ,
[TABLE]
Recall that the averages are time averages, evaluated for given patterns \mbox{\boldmathx}^{(\mu)}.
An equivalent yet slightly different derivation of the mean-field equations (3.8) is this: suppose we average over the dynamics (3.1) at fixed , and then we average all over the dynamics. This gives . Comparing with Equation (3.8), we see that the mean-field approximation corresponds to approximating .
Now, in order to calculate the order parameters (3.9),
[TABLE]
we must solve the mean-field equations (3.8) to obtain the time averages in Equation (3.9). To this end we express the mean field in terms of the order parameters :
[TABLE]
The last equality is only approximate because the -sum in the definition of contains the term . Whether or not to include this term makes only a small difference in the limit of large .
Let us first calculate assuming Equation (3.3), neglecting terms with in Equation (3.10). To make sure that these small -terms do not add up to a substantial correction to the first term, the storage capacity must be small enough. For large values of we assume
[TABLE]
In this case it is sufficient to keep only the first term on the r.h.s. of Equation (3.10) [37]. This approximation yields together with Equation (3.8):
[TABLE]
Applying the definition (3.9) of the order parameter, one finds
[TABLE]
Using that as well as the fact that the bits can only assume the values , one obtains:
[TABLE]
This is a self-consistent equation for . For , it has the solution . This is not the desired solution because means that \mbox{\boldmathx}^{(1)} is not recognised. For , by contrast, there are three solutions, . Figure 3.3 shows results of the numerical evaluation of Equation (3.14) for intermediate values of . For larger than the critical value,
[TABLE]
the three solutions persist. The solution is unstable . This can be shown by computing the derivatives of the free energy of the Hopfield network [1]. In other words, if we start with an initial condition that corresponds to , the network dynamics does not stay there. The other two solutions are stable: when the network is initialised close to \mbox{\boldmathx}^{(1)}, then it converges to , as long as .
The symmetry of the problem dictates that there must also be a solution with . This solution corresponds to the inverted pattern -\mbox{\boldmathx}^{(1)} (Section 2.5). If we start in the vicinity of \mbox{\boldmathx}^{(1)}, then the network is unlikely to converge to -\mbox{\boldmathx}^{(1)}, provided that is large enough. The probability of the dynamical transition \mbox{\boldmathx}^{(1)}\to-\mbox{\boldmathx}^{(1)} vanishes very rapidly as increases and as the noise level decreases. If this transition happens in a simulation in this limit, the network then stays near -\mbox{\boldmathx}^{(1)} for a very long time. Consider the limit where tends to at a finite but large value of . Then the network jumps back and forth between \mbox{\boldmathx}^{(1)} and -\mbox{\boldmathx}^{(1)} at a very small rate. As a result, the order parameter averages to zero. This shows that the limits of large and large do not commute:
[TABLE]
In practice the interesting limit is the left one, that of a large network run for a time much longer than the initial transient, but not infinite. This is precisely where the mean-field theory applies. It corresponds to taking the limit first, at finite but large . This describes simulations where the transition \mbox{\boldmathx}^{(1)}\to-\mbox{\boldmathx}^{(1)} does not occur.
In summary, Equation (3.14) predicts that the order parameter converges to a definite value, , independent of the stored patterns in the limit . When is finite, the limiting value of the order parameter does depend on the stored patterns (Figure 3.2). In this case one averages also over different realisations of the stored patterns, as mentioned above. The value of this average, , determines the average
error probability in the steady state, the average fraction of wrong bits. The steady-state average number of correct bits is given by
[TABLE]
because if is correct, and equal to zero otherwise. The outer average is over different realisations of random patterns (the inner average is over the network dynamics). The second equality follows from Equation (3.2b). Since the l.h.s. of Equation (3.17) equals times , we deduce that
[TABLE]
Since as , the steady-state error probability tends to zero in this limit. This is expected since the stored patterns \mbox{\boldmathx}^{(\mu)} are recognised for small enough values of in the deterministic limit, when the cross-talk term is negligible. But note that the stochastic dynamics slows down as the noise level tends to zero. The lower the noise level, the longer the network remains stuck in local minima, so that it takes longer time to reach the steady state, and to sample the steady-state statistics of . In the opposite limit, , the steady-state error probability tends to , because .
In this noise-dominated limit the stochastic network ceases to function. If one were to assign bits entirely randomly, then half of them would be correct, on average, .
It is important to note that noise can help in another way: it may prevent the network dynamics from converging to mixed states (Section 2.5). This can be seen as follows [33, 1]. To derive the above mean-field result we assumed and for . Mixed states correspond to solutions where an odd number of components of is non-zero, for example:
[TABLE]
Neglecting the cross-talk term, the mean-field equation reads
[TABLE]
In the limit of , the averages converge to the mixed states (2.52) when \mbox{\boldmathm}^{({\rm mix})} is given by Equation (3.19). Averaging over the bits of the random patterns one finds:
[TABLE]
The numerical solution of Equation (3.21) shows that there is a non-zero solution for . Yet this solution is unstable for [33]. In other words, the mixed states have a lower critical noise level than the stored patterns, equal to . For noise levels larger than that, but still smaller than unity, the network can recognise the stored patterns, and it does not converge to mixed states.
However, these results were obtained assuming that only one (or a few) order parameters are not zero. This corresponds to the limit of , where the cross-talk term (Section 2.3) is negligible. The next Section describes a mean-field theory that remains accurate for larger values of .
3.4 Critical storage capacity
The analysis in the preceding Section replaced the sum (3.10) by its first term, . This can only work when is small enough. Now we discuss how to proceed when is not small.
Note that the analysis in Section 2.2 did not assume that is small, but it yielded only the one-step error probability , and we discussed the storage capacity in relation to the one-step error probability. As the network dynamics is iterated, however, the number of errors tends to increase, at least when is large enough so that the cross-talk term matters. Now we describe how to compute for general values of the storage capacity , in order to demonstrate how the errors multiply as one iterates, causing the network to fail.
As before, we store patterns in the network using Hebb’s rule (2.26) and feed pattern \mbox{\boldmathx}^{(1)} to the network. The aim is to determine the order parameter and the corresponding error probability in the steady state for , so that remains finite as . In this case we can no longer approximate the sum in Equation (3.10) just by its first term, because the other terms for may sum up to a contribution that is of the same order as . Instead we must evaluate all to compute the mean field .
The relevant calculation is summarised in Chapter 4 of Ref. [38]. It is also outlined in Section 2.5 of Hertz, Krogh and Palmer [1]. The remainder of this Section follows this outline quite closely. One starts by rewriting the mean-field equations (3.8)
in terms of the order parameters . Using
[TABLE]
we find
[TABLE]
This coupled set of non-linear equations is equivalent to the mean-field equations (3.8).
Now feed pattern \mbox{\boldmathx}^{(1)} to the network. We assume that the network stays close to the pattern \mbox{\boldmathx}^{(1)} in the steady state, so that remains of order unity. The other remain small. When is large, however, we cannot simply approximate the sum over on the r.h.s. of Equation (3.23) by its first term only, because the sum of the remaining (small) terms might not be negligible. Therefore we need to estimate these terms, the other order parameters for .
The trick is to assume that the pattern bits are random, uncorrelated with mean zero [Equations (2.29) and (2.30)]. In this case the order parameters , , become random numbers that fluctuate around zero with variance (this average is over random patterns). We use Equation (3.23) to compute the variance approximately.
In the -sum on the r.h.s of Equation (3.23) we must treat the term separately, because the index appears also on the l.h.s. of this equation. Also the term must be treated separately, as before, because is the index of the pattern that is fed to the network. As a consequence, the calculations of and for proceed slightly differently. We begin with the first case. Using that , and that is an odd function, Equation (3.23) simplifies to:
[TABLE]
The next steps are similar to the analysis of the cross-talk term in Section 2.2. One assumes that the patterns are random, that their bits are independently and identically distributed. In the limit of large and , the sums in Equation (3.24) can then be estimated using the central-limit theorem . For random patterns, the variable
[TABLE]
is a sum of many independent, identically distributed random numbers with mean zero and finite variance. The variable is therefore approximately Gaussian distributed, with mean zero. As a consequence, the distribution of is entirely determined by its variance , and it is indepedent of .
Returning to Equation (3.24), one approximates the sum as an average over the Gaussian distributed variable . This yields:
[TABLE]
Equation (3.26) is the desired result, a self-consistent equation for replacing the mean-field equation (3.14).
In order to determine , we need to estimate the variance featuring in Equation (3.26). To this end, one squares Equation (3.25) and averages the resulting double sum. Since the bits and are independent when , only the diagonal terms in this double sum contribute to the average:
[TABLE]
Here we assumed that is large, and approximated . To evaluate the variance further, it is necessary to estimate the remaining order parameters. One starts again from Equation (3.23) and writes for
[TABLE]
Consider the three terms in the argument of . The term ① is of order unity, it is independent of . The term ③ may be of the same order, because the sum over contains terms. The term ②, by contrast, is small for large values of . Therefore it is a good approximation to Taylor-expand tanh(…) as follows:
[TABLE]
Using one obtains
[TABLE]
Using the fact that and thus , this expression simplifies:
[TABLE]
The goal is now to solve for . Approximating the sum in the second line as an average over the Gaussian distributed variable [Equation (3.25)] gives:
[TABLE]
Defining the parameter
[TABLE]
one can write Equation (3.32) as
[TABLE]
Returning to Equation (3.31), we see that it takes the form
[TABLE]
Solving for one finds for :
[TABLE]
This expression allows us to compute the variance , defined by Equation (3.27). Equation (3.36) shows that the average contains a double sum over the bit index, . Since the bits are independent, only the diagonal terms contribute, so that
[TABLE]
for , but otherwise independent of . The numerator is just , from Equation (3.33). So the variance evaluates to
[TABLE]
In summary there are three coupled equations, for , , and , Equations (3.26), (3.34), and (3.38). They must be solved together to determine how depends on and .
In order to compare with the results described in Section 2.2, we must take the deterministic limit, . In this limit, approaches unity, yet remains finite [1]. Setting in Equation (3.38) but retaining one finds:
[TABLE]
Recall expression (3.18) for the steady-state error probability. Inserting Equation (3.39c) for into this expression we find in the same limit:
[TABLE]
Compare this with Equation (2.39) for the one-step error probability in the deterministic limit. That equation was derived for only one step of the network dynamics, while Equation (3.40) describes the limit of many steps, the long-time or steady-state limit.
Yet it turns out that Equation (3.40) reduces to (2.39) as . To see this, one solves the set of Equations (3.39) by introducing the variable . One obtains the following one-dimensional equation for [1, 34]:
[TABLE]
The relevant solutions are those satisfying , because the order parameter is restricted to this range (transitions to do not occur in the limit ). Figure 3.4 shows the steady-state error probability obtained from Equations (3.40) and (3.41). Also shown is the one-step error probability
[TABLE]
derived in Section 2.2. As stated above, approaches for small . We conclude: in this limit, for small , the error probability does not increase significantly as one iterates the network dynamics. Errors in earlier iterations have little effect on the probability that later errors occur.
The situation is different at larger values of . In that case, significantly underestimates the steady-state error probability. In the hashed region, errors in the dynamics increase the probability of errors in subsequent steps, giving rise to error avalanches. Figure 3.4 illustrates that the steady-state error probability tends to as the parameter increases beyond a critical value, . Equation (3.41) yields
[TABLE]
for the critical storage capacity . When , the network produces just noise. When , by contrast, the network works well. The smaller the storage capacity , the better the network performs.
Figure 3.4 shows that the steady-state error probability changes very abruptly near . Suppose we store 137 patterns with 1000 bits in a Hopfield network. In this case the network can retrieve the patterns with a comparatively small error probability. However, if we try to store one or two more patterns, the network fails to produce output meaningfully related to the stored patterns. This rapid change is an example of a phase transition. In many physical systems one observes similar transitions between ordered and disordered phases [32].
What happens at higher noise levels? The numerical solution of Equations (3.34), (3.26), and (3.38) shows that the critical storage capacity decreases as the noise level increases (smaller values of ). This is shown schematically in Figure 3.5. Below the solid line the error probability is smaller than , so that the network operates reliably (although less so as one approaches the phase-transition boundary). Outside this region the the error probability equals . In this region the network fails. In the limit of small the critical noise level is . In this regime the network is described by the theory explained in Section 3.3, Equation (3.14).
Alternatively these two different phases of the Hopfield network are characterised in terms of the order parameter . We see that below the solid line, while above it, in the limit of large .
3.5 Beyond mean-field theory
The theory summarised in this Chapter rests on a mean-field approximation for the local field, Equation (3.6). The main result is the phase diagram shown in Figure 3.5, derived in the limit . For smaller values of one expects the transition to be less sharp, so that is non-zero also for values of larger than the critical storage capacity .
But even for large values of , the question remains how reliable the mean-field theory really is.
To answer this question, one uses a more accurate theory, based on the so-called replica trick. One starts from the steady-state distribution of for fixed patterns \mbox{\boldmathx}^{(\mu)}. In Chapter 4 we will see that the steady-state distribution for the McCulloch-Pitts dynamics is the Boltzmann distribution
[TABLE]
(the proof in Chapter 4 assumes that the diagonal weights are set to zero). The normalisation factor is called the partition function
[TABLE]
In order to compute the order parameter, one adds a threshold term to the energy function (2.45)
[TABLE]
Then the order parameter is obtained by taking a derivative w.r.t :
[TABLE]
The outer average is over different realisations of random patterns. The logarithm of is averaged using the replica trick. The idea is to represent the average of the logarithm as
[TABLE]
The function looks like the partition function of copies of the system, hence the name replica trick. If one assumes that all copies yield the same order parameter, one obtains the mean-field solution described in Section 3.4. If one allows different copies to have different order parameters (replica-symmetry breaking),
one obtains a more accurate solution for the critical storage capacity [39],
[TABLE]
The mean-field result (3.42) differs only slightly from Equation (3.48). The most precise Monte-Carlo simulations (Section 4.2) for finite values of [40] yield upon extrapolation to
[TABLE]
This is close to, yet significantly different from the best theoretical estimate, Equation (3.48), and also different from the mean-field result (3.42).
To put these results into context, note that for other systems mean-field theories tend to give results much worse than here. Usually, mean-field theories yield at best a qualitative description of a phase transition. For the Hopfield network, by contrast, the mean-field theory works very well because every neuron is connected with every other neuron. This helps to average out the fluctuations in Equation (3.6). In physical systems with local interactions, mean-field theories tend to work better in higher dimensions, because there are more neighbours to average over (Exercise 3.5).
3.6 Correlated and non-random patterns
In the two previous Sections we assumed that the stored patterns are random with independently identically distributed bits. This allowed us to calculate the storage capacity of the Hopfield network using the central-limit theorem. The hope is that the result describes what happens for typical, non-random patterns, or for random patterns with correlated bits. Correlations affect the distribution of the cross-talk term, and thus the storage capacity of the Hopfield network. It has been argued that the storage capacity increases when the patterns are more strongly correlated, while others have claimed that the capacity decreases in this limit (see Ref. [41] for a discussion).
For a set of definite patterns (no randomness to average over), the situation seems to be even more challenging. Yet there is a way of modifying Hebb’s rule to deal with this problem, at least when the patterns are linearly independent. The recipe is explained by Hertz, Krogh, and Palmer [1]. One simply incorporates the overlaps
[TABLE]
into Hebb’s rule. To this end one defines the overlap matrix with elements and writes:
[TABLE]
For orthogonal patterns (), this modified Hebb’s rule is identical to Equation (2.25). For non-orthogonal patterns, the rule (3.51) ensures that all patterns are recognised. Equation (3.51) requires that the matrix is invertible: its columns must be linearly independent (and this implies that the rows are linearly independent too). This limits the number of patterns one can store with the rule (3.51), because implies linear dependence.
For linearly independent patterns one can find the weights iteratively, by successive improvement from an arbitrary starting point. We can say that the network learns the task through a sequence of weight changes. This is the idea used to solve classification tasks with perceptrons (Part II).
3.7 Summary
In this Chapter we analysed the stochastic dynamics of Hopfield networks. We asked under which circumstances the network dynamics can reliably recognise stored patterns. If the stored patterns are random, the performance of the Hopfield network depends on their number, on the number of bits per pattern, and upon the noise level. The storage capacity equals the ratio of the number of stored patterns to the number of bits per pattern. The network operates reliably when this ratio is small, and provided the noise level is not too large. A mean-field analysis of the -limit shows that there is a phase transition in the parameter plane of the Hopfield network (Figure 3.5): when exceeds the critical storage capacity , the network ceases to function.
Hopfield networks share many properties with the networks discussed later on in this book. The most important point is perhaps that introducing noise in the dynamics allows to study the convergence and performance of the network: in the presence of noise there is a well-defined steady state that can be analysed. Without noise, in the deterministic limit, the network dynamics arrests in local minima of the energy function, and may not reach the stored patterns. Naturally the noise must be small enough for the network to function accurately. Finally, the building blocks of Hopfield networks are McCulloch-Pitts neurons and Hebb’s rule for the weights. Many of the algorithms discussed in the coming Chapters use these elements in some form.
3.8 Further reading
The statistical mechanics of Hopfield networks is explained in Introduction to the theory of neural computation by Hertz, Krogh, and Palmer [1]. Starting from the Boltzmann distribution, Chapter 10 in this book summarises how to compute the order parameters, and how to evaluate the stability of the corresponding solutions. For more details on the replica trick, see the books by Müller, Reinhard and Strickland [37] and by Engel and van den Broeck [42], as well as the review article [43].
Chapter 4 The Boltzmann distribution
In Chapter 2 we saw that the deterministic dynamics (2.5) of Hopfield networks admits the Lyapunov function
[TABLE]
if the weights are symmetric, and . In this Chapter111In this Chapter we set the diagonal weights to zero. we show that the asynchronous stochastic McCulloch-Pitts dynamics (3.1) converges to a steady state where the state vector follows the Boltzmann distribution
[TABLE]
The stochastic dynamics (3.1) is closely related to that of Markov-chain Monte-Carlo algorithms, designed to efficiently sample from the Boltzmann distribution.
We also discuss how to solve optimisation tasks by Monte-Carlo simulation: one assigns a suitable energy to each configuration , so that the function H(\mbox{\boldmaths}) has global minimum for the optimal configuration \mbox{\boldmaths}_{\rm min}. The stochastic dynamics finds low-energy configurations (but not necessarily \mbox{\boldmaths}_{\rm min}), in particular if one iteratively decreases the noise level by increasing (simulated annealing [44]).
Last but not least we look at Boltzmann machines [15, 45, 14, 46, 47], stochastic Hopfield networks with hidden neurons that are neither used for input nor for output. Boltzmann machines can be trained to learn the properties of a distribution P_{\rm data}(\mbox{\boldmathx}) of binary input patterns . The idea is to iteratively change the weights in Equation (4.1) until the Boltzmann distribution represents the input distribution. This idea, to iterate the weights until the network learns the input distribution , is used in a slightly different form in supervised learning (Part II). Boltzmann machines are closely related to Hopfield networks. Without hidden neurons, both models learn
to represent two-point correlations of pattern bits.
When important information about the inputs is encoded in higher-order correlations, one can use hidden neurons to represent
these correlations. Generally Boltzmann machines are hard to train, in particular if they have many hidden neurons. Restricted Boltzmann machines are neural networks with hidden neurons, but with fewer connections: only those between visible and hidden neurons are allowed. These neural networks can be fairly efficiently trained and can solve a number of different tasks. Apart from learning a distribution of input patterns, they can for instance be trained to recognise incomplete input patterns, and to classify inputs [25].
4.1 Convergence of the stochastic dynamics
We begin by showing that the stochastic dynamics (3.1) has a steady state where is distributed according to the Boltzmann distribution (4.2). To this end, we consider an alternative yet equivalent formulation of the network dynamics. It consists of two parts. First, choose a neuron randomly, number say. Second, change to with probability
[TABLE]
To explore the relation between the stochastic rules (4.3) and (3.1), we use that
[TABLE]
with local field . To derive Equation (4.4), we assumed that the weights are symmetric, and that the diagonal weights vanish. The result is obtained with a calculation similar to the one leading to Equation (2.48), except that we have non-zero thresholds here. To proceed, we break the rule (4.3) up into different cases. The state of neuron changes with probability
[TABLE]
Comparing with Equation (3.1) we conclude that the two schemes (3.1) and (4.3) are equivalent under the assumptions made ( and ). Note that Equation (4.3) is more general than the stochastic Hopfield dynamics, because it does not require the energy function to be of the form (4.1). In particular it is neither needed that the weights are symmetric, nor that the diagonal weights vanish.
Equations (3.1) and (4.3) are not equivalent if these conditions are not satisfied (Exercise 4.1).
The rule (4.3) defines a Markov chain of states
[TABLE]
As before, the index counts the iteration steps. A Markov chain is a memoryless random sequence of states defined by transition probabilities p(\mbox{\boldmaths}^{\prime}|\mbox{\boldmaths}) from state to \mbox{\boldmaths}^{\prime} [48]. The transition probability p(\mbox{\boldmaths}^{\prime}|\mbox{\boldmaths}) connects arbitrary states. One distinguishes between local moves where only one neuron may change, as above, and global moves where many neurons may change their states in a single step.
In both cases, an update consists of two parts. First, a new state \mbox{\boldmaths}^{\prime} is suggested with probability q(\mbox{\boldmaths}^{\prime}|\mbox{\boldmaths}). Second, the new state \mbox{\boldmaths}^{\prime} is accepted with acceptance probability
[TABLE]
As result, the transition probability is given by a product of two factors
[TABLE]
These steps are repeated many times, creating the chain of states (4.6).
The Markov chain defined by the transition probability (4.8) has the Boltzmann distribution (4.2) as a steady-state distribution if the detailed-balance condition is satisfied:
[TABLE]
Note that this is a sufficient condition, not a necessary one [49]. There are Markov chains that do not satisfy detailed balance but still have a steady state (Exercise 4.4). Usually detailed balance implies not only that the Markov chain has P_{\rm B}(\mbox{\boldmaths}) as a steady-state distribution, but also that the distribution of states generated by the sequence (4.6) converges to P_{\rm B}(\mbox{\boldmaths}), see Ref. [48] for details.
To prove that the detailed-balance condition (4.9) holds for the transition probability (4.8), assume that a single neuron is picked randomly with uniform probability
[TABLE]
where is the number of neurons in the network. Since does not depend on either or \mbox{\boldmaths}^{\prime}, the probability of suggesting a new state is clearly symmetric. Equations (4.2), (4.7) then imply:
[TABLE]
This demonstrates that the Boltzmann distribution is a steady state of the Markov chain defined by (4.7), (4.8), and (4.10). As a consequence, the Boltzmann distribution is a steady state of the Markov chain. If the simulation converges to the steady state (as it usually does), then states visited by the Markov chain are distributed according to the Boltzmann distribution. This also means that the steady-state distribution for the Hopfield model is the Boltzmann distribution, as stated in Section 3.5.
It is important to stress that Equation (4.9) is a condition for the transition probability p(\mbox{\boldmaths}^{\prime}|\mbox{\boldmaths})=q(\mbox{\boldmaths}^{\prime}|\mbox{\boldmaths})p_{\rm a}(\mbox{\boldmaths}^{\prime}|\mbox{\boldmaths}), not just for the acceptance probability p_{\rm a}(\mbox{\boldmaths}^{\prime}|\mbox{\boldmaths}). For the local moves discussed above, is a constant, so that p(\mbox{\boldmaths}^{\prime}|\mbox{\boldmaths})\propto p_{\rm a}(\mbox{\boldmaths}^{\prime}|\mbox{\boldmaths}). In this case it is sufficient to check the detailed-balance condition for the acceptance probability. In general, and in particular for global moves, it is necessary to include q(\mbox{\boldmaths}^{\prime}|\mbox{\boldmaths}) in the detailed-balance check [50].
4.2 Monte-Carlo simulation
The Markov chain described in the previous Section is the basis for the Markov-chain Monte-Carlo algorithm. This method is widely used in statistical physics and in mathematical statistics. It is therefore important to understand the connections between the different formulations.
The Boltzmann distribution describes the probabilities of observing configurations of a large class of physical systems in their steady states [32]. The statistical mechanics of systems with energy function (also called Hamiltonian) shows that their configurations are distributed according to the Boltzmann distribution in thermodynamic equilibrium at a given temperature (in this context where is the Boltzmann constant), and free from any other constraints. If we denote the configuration of a system by the vector , then the Boltzmann distribution takes the form (4.2). The normalisation factor Z=\sum_{\mbox{\boldmaths}}{\rm e}^{-\beta H(\mbox{\boldmaths})} is also called partition function. For systems with a large number of interacting degrees of freedom, the partition function can be very expensive to compute, because the sum over contains many terms. Therefore, instead of computing the distribution directly one generates a Markov chain of states with a suitable transition probability, for instance (4.3).
In practice one often uses a slightly different form of the transition probability (Metropolis algorithm [51]). Assuming that is constant, one takes:
[TABLE]
with \Delta H=H(\mbox{\boldmaths}^{\prime})-H(\mbox{\boldmaths}) as before. That the Metropolis rates obey the detailed-balance condition (4.9) can be seen as follows:
[TABLE]
The Metropolis algorithm is summarised in Algorithm 2. It provides an elegant way of computing the average of an observable A(\mbox{\boldmaths}) over the Boltzmann distribution of :
[TABLE]
This particular way of evaluating the average is a special case of the more general method of importance sampling [52]. The central-limit theorem implies that the error of this estimate for decreases as increases. The prefactor is determined by the correlations between subsequent terms in the sum (4.20): the states in the sequence (4.6) are correlated, in particular when the moves are local, because then subsequent configurations are similar. Generating many quite strongly correlated samples from a distribution is not a very efficient way of sampling this distribution. Sometimes it may be more efficient to suggest global moves instead, in order to avoid that subsequent
states in the Markov chain are similar. But it is not guaranteed that global moves lead to weaker correlations. For global moves, may be more likely to assume large positive values, so that fewer suggested moves are accepted. As a consequence the Markov chain may stay in certain states for a long time, increasing correlations in the sequence. Usually a compromise is most efficient, moves that are neither local nor global. In summary, the convergence of Monte-Carlo sampling is quite slow. This motivated Sokal to begin his lecture notes on Monte-Carlo simulation with the warning [49]
Monte Carlo is an extremely bad method; it should be used only when all alternative methods are worse.
Monte-Carlo algorithms are very widely used, and the original reference for the Metropolis algorithm [51] is generally considered one of the most significant scientific papers in computational physics. Sokal’s point is of course that many problems cannot be solved in any other way, so that Monte-Carlo simulation is the only option. But we should be aware of the shortcomings of the method. The same caution applies more generally to the topic of this book, machine-learning algorithms with neural networks.
4.3 Simulated annealing
Combinatorial optimisation problems admit or configurations -
too many to find the optimal one by complete enumeration when is large.
An alternative strategy is to assign an energy H(\mbox{\boldmaths}) to each configuration so that is minimal at the optimal configuration \mbox{\boldmaths}_{\rm min}. One minimises H(\mbox{\boldmaths}) by Monte-Carlo simulation, using that the Monte-Carlo dynamics tends to decrease when the temperature is low, Figure 4.1. A common strategy is to lower the temperature on the fly. In the beginning of the simulation, the temperature is high, so that the dynamics first explores the rough features of the energy landscape. When the temperature is lowered, the dynamics perceives finer and finer features of H(\mbox{\boldmaths}). The hope is that it ends up in the global minimum H_{\rm min}=H(\mbox{\boldmaths}_{\rm min}) at zero temperature, where P_{\rm B}(\mbox{\boldmaths})=0 when H(\mbox{\boldmaths})>H_{\rm min} and P_{\rm B}(\mbox{\boldmaths})>0 only for H(\mbox{\boldmaths})=H_{\rm min}. This method is called simulated annealing [44], see also Section 10.9 in Numerical Recipes [53]. Slowly lowering the temperature during the simulation mimics the slow cooling of a physical system. It passes through a sequence of quasi-equilibrium Boltzmann distributions with lower and lower temperatures, until the system hopefully finds the global minimum .
For a number of combinatorial optimisation problems one can write
down energy functions that have the same form as Equation (2.50) with symmetric weights [54]. Since , one can always assume that the diagonal weights
vanish, because they make only a constant contribution to . In short, one can use the Hopfield dynamics (3.1) to minimise . The travelling-salesman problem has been solved in this way [54, 1], gradually reducing the noise level as one iterates
the stochastic dynamics.
It is by no means necessary to use a Hopfield model for this purpose. Instead we can just use the stochastic dynamics (4.3) or the Metropolis algorithm (4.12) to solve combinatorial optimisation problems by simulated annealing.
Nevertheless, a crucial step is to find a suitable energy function.
As an example, consider the double-digest problem. It arose when sequencing the human genome [55, 56]. The human genome sequence was first assembled by piecing together overlapping DNA segments in the right order by making sure that overlapping segments share the same DNA sequence. To this end it is necessary to uniquely identify the DNA segments. The actual DNA sequence of a segment is a unique identifier. But it is sufficient and more efficient to identify a DNA segment by a fingerprint, for example the sequence of restriction sites. These are short subsequences (four or six base pairs long) that are recognised by enzymes that cut (digest) the DNA strand precisely at these sites. A DNA segment is identified by the types and locations of restriction sites that it contains, the so-called restriction map.
When a DNA segment is cut by two different enzymes one can experimentally determine the lengths of the resulting fragments. Is it possible to determine how the cuts were ordered in the DNA sequence of the segment from the fragment lengths, to find the restriction map? This is the double-digest problem [55]. In a double-digest experiment, a given DNA sequence is first digested by one enzyme ( say). Assume that this results in fragments with lengths (). Second, the DNA sequence is digested by another enzyme, . In this case fragments are found, with lengths . Third, the DNA sequence is digested with both enzymes and , yielding fragments with lengths , see Table 4.1 for examples. The task is now to determine all possible orderings of the - and -cuts that result in fragments with lengths . Since the solutions of the double-digest problem are degenerate, an important question is to determine how many distinct solutions there are (Exercise 4.5).
To write down an energy function, denote the ordered set of fragment lengths produced by digesting with enzyme by , where . Similarly () for fragment lengths produced by enzyme , and () for fragment lengths produced by digesting first with and then with . Permutations and of the sets and result in a set of -fragments that we call . Solutions of the double-digest problem correspond to permutations that yield . A suitable energy function is therefore
[TABLE]
and configuration space is the space of all permutation pairs \mbox{\boldmaths}=[\sigma,\mu]. Local moves in configuration space correspond to inversions of short subsequence of and/or . One can show that the corresponding q(\mbox{\boldmaths}^{\prime}|\mbox{\boldmaths}) is symmetric (Exercise 4.5). As mentioned above, this is necessary for the stochastic dynamics to converge in its simplest form, Equation (4.12) and Algorithm 2.
For the simulation one chooses a larger temperature to begin with, so that the stochastic dynamics explores the rough features of the energy landscape at first. As the simulation proceeds, the temperature is gradually reduced. This allows the dynamics to learn finer features of the landscape, as described above.
4.4 Boltzmann machines
Boltzmann machines are generalised Hopfield networks that can learn to approximate data distributions of binary input patterns.
Boltzmann machines differ from Hopfield networks
in two essential ways. First, instead of using Hebb’s rule, the weights are adjusted until the Boltzmann machine approximates the data distribution precisely. The weights are iteratively refined to minimise the difference between the data distribution and the model (the Boltzmann distribution). Nevertheless, this procedure is closely related to Hebb’s rule, as we shall see. Second, to represent higher-order correlations between bits of input patterns, Boltzmann machines employ hidden neurons.
We begin with Boltzmann machines without hidden neurons (Figure 4.2), because they are simpler to analyse. Then we discuss why hidden neurons are necessary to learn the properties of general input distributions P_{\rm data}(\mbox{\boldmathx}) of binary inputs . The training algorithm for Boltzmann machines with hidden neurons is described in Section 4.5.
The goal of the training algorithm is to find weights so that the Boltzmann distribution
[TABLE]
approximates the distribution P_{\rm data}(\mbox{\boldmathx}) as precisely as possible. Here and in the remainder of this Chapter we set . The input patterns have binary bits [Equation (2.1)] with values . The weight matrix is symmetric, , and its diagonal elements are set to zero, . In this Section we also set the thresholds to zero.
The Boltzmann machine is trained by iteratively adjusting the weights , using a sequence of input patterns \mbox{\boldmathx}^{(\mu)} (), independently sampled from the data distribution P_{\rm data}(\mbox{\boldmathx}). This is achieved by maximising the likelihood \mathscr{L}=\prod_{\mu=1}^{p}P_{\rm B}(\mbox{\boldmaths}=\mbox{\boldmathx}^{(\mu)}) that the Boltzmann machine produces the sequence \mbox{\boldmathx}^{(1)},\ldots,\mbox{\boldmathx}^{(p)} of input patterns. Any pattern may appear more than once in the sequence, with frequency proportional to P_{\rm data}(\mbox{\boldmathx}). Maximising therefore corresponds to approximating the data distribution as accurately as possible. Usually one maximises the logarithm of the likelihood, the log-likelihood function
[TABLE]
The logarithm is a monotonic function, so the log-likelihood has its maximum at the same weight values as the likelihood. Taking the logarithm simplifies the analysis of the learning algorithm, because \log P_{\rm B}(\mbox{\boldmaths}=\mbox{\boldmaths}^{(\mu)}) is simply a quadratic function of . Also, a learning algorithm based on the log-likelihood is usually more stable numerically.
A different reasoning behind maximising the log-likelihood starts from the Kullback-Leibler divergence, defined as
[TABLE]
Terms in the sum with P_{\rm data}(\mbox{\boldmathx}^{(\mu)}){=0} are set to zero, and is defined to equal infinity when there are patterns for which but . The Kullback-Leibler divergence is a measure of the difference between the two distributions: is non-negative, and it assumes its global minimum for P_{\rm data}(\mbox{\boldmathx}^{(\mu)})=P_{\rm B}(\mbox{\boldmaths}=\mbox{\boldmathx}^{(\mu)}), see Exercise 4.6. We infer from Equation (4.24) that minimising corresponds to maximising .
To find the global maximum of the log-likelihood, we use gradient ascent: we repeatedly change the weights by adding increments
[TABLE]
The small parameter is the learning rate. The gradient points in the steepest uphill direction of . The idea is to take many small uphill steps until one hopefully (but not necessarily) reaches the global maximum. Since the likelihood is a product of many possibly quite small factors, can become very small. This can lead to numerical instabilities. Maximising instead of can be more stable because it yields an additional factor in the gradient: .
To evaluate the gradient of we start from Eq. (4.23)
[TABLE]
This expression assumes that the diagonal weights vanish, just like Equation (4.22). The first step is to evaluate the derivative of
[TABLE]
To compute we use the chain rule together with
[TABLE]
This relation is valid for symmetric weights and provided that and . In Equation (4.28), is the Kronecker delta, if and zero otherwise (Chapter 2). In particular, only if and . Otherwise the product of Kronecker deltas equals zero. Equation (4.28) is illustrated by the following story (a modification of a well-known maths joke):
The linear function, , and the constant function are going for a walk. When they suddenly see the derivative approaching, the constant function gets worried. "I’m not worried" says the function confidently, "I’m not put to zero by the derivative." When the derivative comes closer, it says "Hi! I’m . How are you?"
The moral is: since and are independent variables, . Equation (4.28) reflects the same principle: the weights and are independent variables unless their indices agree. Equation (4.28) is valid for off-diagonal weights, and there are two terms on the r.h.s. because the weights are symmetric.
Returning to the derivative of with respect to , one finds using Equation (4.28):
[TABLE]
The last equality defines the two-point correlations of the model, , computed using the steady-state distribution (4.22) of the Boltzmann machine. Evaluating the derivative of the second term in Equation (4.26) gives:
[TABLE]
In summary,
[TABLE]
Here is the two-point correlation of the input data. Using (4.25), the learning rule becomes:
[TABLE]
where we dropped a factor of that only affects the numerical value of the learning rate . The weight increments are determined by the two-point pattern correlations, just like Hebb’s rule (2.25). The first term on the r.h.s. of Eq. (4.32) has precisely the same form as Equation (2.25), a sum over two-point correlations of the input patterns. The second average is over the steady-state distribution (4.22) of the Boltzmann machine. The learning rule takes the form of the difference between two-point correlations because the task is to minimise the difference between two distributions. It is plausible that the learning rule may converge
because the weight increments vanish when the model correlations equal the data correlations.
The average can be approximated by numerical simulation of the McCulloch-Pitts dynamics
[TABLE]
with and . One must iterate Equation (4.33) until the system has reached its steady state, long enough so that any initial transient becomes negligible.
The training algorithm can be summarised as follows. One initialises all weights and computes from the given sequence of input patterns. One estimates by numerical simulation of the dynamics of the Boltzmann machine, and changes the weights using (4.32). This step is iterated, either with a sequence of new inputs, or with the same inputs but in permuted sequence. In each iteration one must compute again, because the weights changed. This procedure is quite slow, because it usually takes long simulations to estimate accurately, in each iteration of the learning algorithm.
There is a more fundamental problem [47]. Like Hebb’s rule, the learning rule (4.32) relies entirely upon two-point correlations of the input bits. This means that the Boltzmann machine cannot learn higher-order correlations between inputs. However, two-point correlations may not be sufficient to represent the information encoded in the input data. To illustrate this point, consider the Boolean XOR function (Exercise 2.13 and Chapter 5). It can be encoded in the four patterns , , , and . The first two components represent the input to the XOR function. The third component represents the output, which depends on both input variables as prescribed by the XOR function. Let us define an input distribution that reflects these three-point correlations by assigning to the four patterns, and setting otherwise. A Boltzmann machine with three neurons cannot represent this input distribution, because there is no energy function of the form (4.1) that has four minima at these patterns. So the three-point correlations encoded in the four patterns cannot be represented in terms of a Boltzmann machine in its simplest form.
Also Hopfield networks fail for the XOR function: the four states are not attractors of a Hopfield network with three neurons (Exercise 2.13). One could consider neural networks with third- or higher-order interactions [47],
[TABLE]
(Exercise 2.7). But the number of weights proliferates as the order increases, rendering the training very slow.
An alternative is to use Boltzmann machines with hidden neurons, that are neither input nor output units. The idea is that the hidden neurons can learn to represent such correlations [47]. The learning rule for the Boltzmann machines with hidden neurons is very similar to Equation (4.32), but when the number of hidden neurons is large, the Boltzmann machine is very slow to train. It is more efficient to remove all weights between visible neurons, and between hidden neurons. This is described in the next Section.
4.5 Restricted Boltzmann machines
Restricted Boltzmann machines [57] consist of visible and hidden neurons arranged in an undirected bipartite graph (Figure 4.3): the only connections are between neurons of different kinds, there are no connections between visible neurons, no connections between hidden neurons either. So the energy function for a restricted Boltzmann machine for visible neurons and hidden neurons takes the form
[TABLE]
with weights and thresholds and . The McCulloch-Pitts dynamics reads
[TABLE]
The diagonal weights are assumed to vanish, but the weight matrix is not required to be symmetric (Exercise 4.9). Since most often , it is usually not even a square matrix.
The learning rule for the weights of the restricted Boltzmann machine is derived using gradient ascent on the log-likelihood for a single pattern \mbox{\boldmathx}^{(\mu)}:
[TABLE]
Proceeding as in the previous Section one finds:
[TABLE]
The first average,
[TABLE]
can be evaluated further, using the fact that there are no connections between the hidden units. Making use of the update rule (4.36a) we find
[TABLE]
just like Equation (3.7). For the other sums in Equation (4.5) we use the normalisation condition 1=\sum_{h_{k}=\pm 1}P(h_{k}|\mbox{\boldmathv}=\mbox{\boldmathx}^{(\mu)}) to obtain:
[TABLE]
The second average on the r.h.s. of Equation (4.38) simplifies to
[TABLE]
The average is computed by Monte-Carlo sampling, using the McCulloch-Pitts dynamics (4.36) to generate the Markov chain
[TABLE]
In the limit , the steady state of this sequence is distributed according to the model distribution, the Boltzmann distribution
with energy function (4.35). In general only the asynchronous McCulloch-Pitts dynamics can be proven to converge (Sections 2.5 and 4.1). Here, however, the Markov chain can be generated more efficiently by updating all hidden neurons \mbox{\boldmathh}_{t} at the same time, given \mbox{\boldmathv}_{t}, because the components of \mbox{\boldmathh}_{t} are independent from each other since there are no connections between them. In the same way the visible neurons \mbox{\boldmathv}_{t} are updated in parallel. To speed up the computation further, one usually only iterates for a finite number of steps, up to say, and initialises the chain with \mbox{\boldmathv}_{t=0}=\mbox{\boldmathx}^{(\mu)}. After steps one approximates
[TABLE]
This algorithm is called contrastive-divergence or CD- algorithm (Algorithm 3). Since the average over the model distribution is approximated [Equation (4.43)], this algorithm does not precisely correspond to gradient ascent.
In summary,
[TABLE]
The analogous learning rules for the thresholds read:
[TABLE]
The derivation of Equations (4.45) is left as an exercise (Exercise 4.10). Restricted Boltzmann machines may have 0/1 neurons with state values [math] and instead of and . For 0/1 neurons, the CD- algorithm is slightly different (Exercise 4.11).
Figure 4.4 illustrates how a restricted Boltzmann machine can learn to complete patterns, using the bars-and-stripes data set [47, 25] as an example. To begin with, the restricted Boltzmann machine is trained using the CD- algorithm. Then consider a partially obscured pattern. Assume for instance that only the upper row of its bits is known: (), (), and (). The remaining bits are obscured, their states are set to zero as shown in Figure 4.4(b). To complete the pattern, one samples from the Boltzmann distribution keeping fixed (clamping these neurons), and iterates the McCulloch-Pitts dynamics for the remaining ones. Panel (b) shows how the machine outputs the correct completed pattern.
This requires hidden neurons, because the data distribution has non-zero three-point correlations [47] (Exercise 2.13). In general a restricted Boltzmann machine can approximate a distribution of binary input data better with more hidden neurons. How many are needed [58, 59]? The answer is not known in general, but it is plausible hidden neurons are sufficient, because each hidden neuron can encode one of the binary input patterns (winning neuron, Section 7.1). More precisely, it can be shown that hidden neurons are sufficient to reach arbitrarily small Kullback-Leibler divergence [60]. An upper bound for the Kullback-Leibler divergence was derived in Refs. [61, 62]:
[TABLE]
Here denotes the integer part. Figure 4.5 illustrates this result. It demonstrates how well a restricted Boltzmann machine approximates the
XOR distribution introduced in Section 4.4. The Figure shows how the Kullback-Leibler divergence depends on the number of hidden neurons (Exercise 4.7).
In this example there are inputs. We see that three hidden neurons are sufficient to allow the restricted Boltzmann machine to approximate the data distribution very precisely, consistent with Equation (4.46). In general, however, the CD- algorithm is not guaranteed to converge to the optimal solution corresponding to the estimate (4.46).
Restricted Boltzmann machines are generative models, they can be used to sample from a distribution the machine has learned [25]. In this way, the machine can complete missing information, as illustrated in Figure 4.4. Restricted Boltzmann machines can also learn to classify patterns, by learning a distribution of binary inputs together with their labels. To this end one splits the visible neurons into input neurons and output neurons with labels or targets. This is a supervised-learning task, the subject of Part II. Recently, restricted Boltzmann machines were used to represent and analyse ground-state wave functions of quantum many-body systems [63].
4.6 Summary
This Chapter dealt with the Boltzmann distribution. Two main points are, first, that the stochastic McCulloch-Pitts dynamics (3.1) has the Boltzmann distribution as a steady state. Second, the update rule (3.1) is a special case of the Markov-chain Monte-Carlo algorithm, for Hopfield models with the energy function (2.45). Since this algorithm tends to decrease the energy function, it can be used to solve complex optimisation problems. In simulated annealing one gradually reduces the noise level as the simulation proceeds. This mimics the slow cooling of a physical system, usually an efficient way of bringing the system into its global optimum.
Boltzmann machines are generalisations of Hopfield networks that can learn distributions of binary data by iteratively changing the weights and thresholds until the corresponding Boltzmann distribution approximates the data distribution. The learning rule is derived using gradient ascent on a target function, in this case the log-likelihood. A related idea is used for training deep neural networks with stochastic gradient descent (Part II). To learn general input distributions of binary patterns requires hidden neurons, also this is a central topic of Part II. Since Boltzmann machines with many hidden neurons are hard to train, one removes connections that are not needed. Restricted Boltzmann machines have connections only between visible and hidden neurons.
4.7 Further reading
Older but still good references for Monte-Carlo methods in statistical physics are the book Monte Carlo methods in Statistical Physics edited by Binder [52], and Sokal’s lecture notes [49]. Some historical notes are found in Ref. [64].
For a concise introduction to Boltzmann machines, refer to Information theory, inference and learning algorithms by MacKay [47], or to Machine learning: a probabilistic perspective by Murphy [65]. Ref. [66] is a more mathematical review of restricted Boltzmann machines.
How many hidden neurons should one allow for in a restricted Boltzmann machine? Little is known apart from the upper bound (4.46) for the Kullback-Leibler divergence, and simulations [60] show that one can get very precise approximations of with less hidden neurons than stipulated by Equation (4.46).
Deep-belief networks consist of layers of restricted Boltzmann machines [2]. Contrastive-divergence training for such deep architectures (networks with many layers) is one of the first examples of deep-learning algorithms [67] (Chapter 7).
Helmholtz machines [68, 69] are generalisations of Boltzmann machines designed as more efficient generative models. They consist of two networks, encoder and decoder, just like variational autoencoders (Section 10.6). The encoder (called recognition model in Ref. [69]) generates a compressed representation of the data distribution, and the decoder (generative model) generates patterns from the compressed representation.
Part II Supervised learning
The Hopfield network described in Part I recognises patterns stored using Hebb’s rule. Its neurons act as inputs and outputs. After feeding a distorted pattern into the network, the network dynamics runs until it reaches a steady state which hopefully corresponds to the stored pattern closest to the distorted one. In this case, the network classifies the distorted pattern by associating it with the closest one amongst the stored patterns.
Part II describes supervised learning, a different way of solving classification tasks with neural networks using labeled data sets. The machine-learning repository [70] at the University of California Irvine contains a number of such data sets. An example is the iris data set which lists certain properties of 150 iris plants. For each plant, four attributes are given (Figure 5.1): its sepal length, sepal width, petal length, and petal width. Each entry in the data set contains a label (or target) that says which class the plant belongs to: iris setosa, iris versicolor, or iris virginica. This data set was described by the geneticist R. A. Fisher [71].
The machine-learning task is to adjust weights and thresholds of a neural network so that it correctly determines the class of each plant from its attributes. To this end one uses a training set of labeled data. Each set of attributes is an input pattern to the network. The neural network is supposed to output the correct label (or target), in this case whether the plant is an iris setosa, iris versicolor, or iris virginica. One compares the network output with the corresponding target, for all input patterns in the training set, and changes the weights and thresholds until the network computes the correct output for each input pattern. The crucial question is whether the trained network can generalise: does it find the correct labels for an input pattern not contained in the training set?
The networks used for supervised learning are called perceptrons [10]. They consist of layers of McCulloch-Pitts neurons: usually some layers
of hidden neurons, and an output layer. We briefly discussed the idea of hidden neurons in connection with restricted Boltzmann machines (Section 4.5), but perceptrons have different layouts, and they are trained in a different way. The layers are usually arranged from the left (input) to the right (output). All connections are one-way, from neurons in one layer to neurons in the layer immediately to the right. There are no connections between neurons in a given layer, or back to layers on the left. This arrangement ensures convergence of the training algorithm (stochastic gradient descent). During training with this algorithm, the network parameters are changed iteratively. In each step, an input is applied, and weights and thresholds of the network are updated to reduce the output error. Loosely speaking, each step corresponds to adding a little bit of Hebb’s rule to the weights. This is repeated until the network classifies the training set correctly.
Stochastic gradient descent for multilayer perceptrons has received much attention recently, after it was realised that networks with many hidden layers can be trained to reliably recognise and classify image data (deep learning).
Chapter 5 Perceptrons
In 1958 Rosenblatt [10] suggested to connect McCulloch-Pitts neurons into layered feed-forward networks to process information. He referred to these networks as perceptrons. The layout is illustrated in in Figure 5.2. The leftmost layer consists of input terminals, drawn in black in Figure 5.2. To the right follow two layers of McCulloch-Pitts neurons. The rightmost layer consists of output neurons. The intermediate layer is a hidden layer. The states of its neurons are not read out.
All connections are one-way: every neuron feeds forward, only to neurons in the layer immediately to the right. There are no connections within layers, no back connections, no connections that skip a layer. There are input terminals. As in Part I, we denote the input patterns
by
[TABLE]
The index labels the different input patterns. The hidden neurons compute
[TABLE]
with weights and thresholds . The function is an activation function, and its argument is called local field (Section 1.2).
The output neurons of the network shown in Figure 5.2 perform the computation
[TABLE]
The index labels the output neurons with weights , and with thresholds .
A classification problem is given by a training set of input patterns \mbox{\boldmathx}^{(\mu)} and the corresponding target vectors
[TABLE]
The idea is to choose all weights and thresholds so that the network produces the desired output:
[TABLE]
In the Hopfield networks described in Part I, the weights were assigned using Hebb’s rule (2.26). Perceptrons, by contrast, are trained by iteratively updating their weights and thresholds until Equation (5.5) is satisfied. This is achieved by repeatedly adding small multiples of Hebb’s rule to the weights (Section 5.2). An alternative approach is to define an energy function, a function of the weights of the network, that has a global minimum when Equation (5.5) is satisfied. The network is trained by taking small steps in weight space that reduce the energy function (gradient descent, Section 5.3).
5.1 A classification problem
To illustrate how perceptrons can solve classification problems, consider the simple example shown in Figure 5.3. There are ten patterns, each has two real-valued components:
[TABLE]
In Figure 5.3 the patterns are drawn as points in the - plane, the input plane. There are two classes of patterns, with targets :
[TABLE]
A single neuron suffices to classify these patterns, a binary threshold unit with activation function , consistent with the possible target values. Since there is only one neuron, we can arrange the weights into a weight vector
[TABLE]
The network performs the computation
[TABLE]
Here \mbox{\boldmathw}\cdot\mbox{\boldmathx}=w_{1}x_{1}+w_{2}x_{2} is the scalar product between the vectors and (Chapter 2).
This example allows us to find a geometrical interpretation of the classification problem. We see in Figure 5.3 that the patterns fall into two clusters: to the left and to the right. We can classify the patterns by drawing a line that separates the two clusters, so that everything to the right of the line has , while the patterns to the left of the line have . This line is called the decision boundary. To find the geometrical significance of Equation (5.9), let us put the threshold to zero for a moment, so that
[TABLE]
Then the classification problem takes the form
[TABLE]
To evaluate the scalar product, we write the vectors as
[TABLE]
Here |\mbox{\boldmathw}|=\sqrt{w_{1}^{2}+w_{2}^{2}} denotes the norm of the vector , and and are the angles of the vectors with the -axis. Then \mbox{\boldmathw}\cdot\mbox{\boldmathx}=|\mbox{\boldmathw}||\mbox{\boldmathx}|\cos(\alpha-\beta)=|\mbox{\boldmathw}||\mbox{\boldmathx}|\cos\varphi, where is the angle between the two vectors. When is between and , the scalar product is positive, otherwise negative. As a consequence, the network classifies the patterns in Figure 5.3 correctly if the weight vector is orthogonal to the decision boundary, as shown in Figure 5.4.
What is the role of the threshold ? Equation (5.9) implies that the decision boundary is parameterised by \mbox{\boldmathw}\cdot\mbox{\boldmathx}=\theta, or
[TABLE]
Therefore the threshold determines the intersection of the decision boundary with the -axis (equal to ). This is illustrated in Figure 5.5.
The decision boundary – the straight line orthogonal to – should divide inputs with positive and negative targets. If such a line can be found, then the problem can be solved with a single neuron. We say that the problem is linearly separable. Conversely, if no such line exists, the problem not linearly separable. This can occur only when . Figure 5.6 shows two problems. The left one is linearly separable, the right one is not.
Other examples are Boolean functions. A Boolean function takes binary inputs and has one binary output. The Boolean AND function (two inputs) is illustrated in Figure 5.7.
The value table of the function is shown on the left. The graphical representation is shown on the right of the Figure ( corresponds to and to ). Also shown is the decision boundary of a binary threshold unit and its weight vector .
It is important to note that the decision boundary is not unique, neither are the weight vector and threshold value that solve the problem. The norm of the weight vector, in particular, is arbitrary. Figure 5.8 illustrates that the Boolean XOR function is
not linearly separable [11].
There are 16 different Boolean functions of two variables. Only two are not linearly separable (Exercise 5.2), XOR (Figure 5.8) and XNOR.
Up to now we discussed only one single neuron. If the classification problem requires several output neurons, each has its own weight vector \mbox{\boldmathw}_{i} and threshold . We can group the weight vectors into a weight matrix as in Part I, so that the row vectors \mbox{\boldmathw}_{i}^{\sf T} are the rows of the weight matrix .
5.2 Iterative learning algorithm
In the previous Section we determined the weights and threshold for the Boolean AND function by inspection (Figure 5.7). Now we discuss an algorithm that finds the weights iteratively. It is illustrated in Figure 5.9. In panel (a), the pattern \mbox{\boldmathx}^{(8)} () is on the wrong side of the decision boundary.
In order to correct this error, one turns the decision boundary anti-clockwise. To this end, one adds a small multiple of the pattern vector \mbox{\boldmathx}^{(8)} to the weight vector
[TABLE]
The parameter is called the learning rate. It must be small, so that the decision boundary is not rotated too far. The result is shown in panel (b). Panel (c) shows another case, where pattern \mbox{\boldmathx}^{(4)} () is on the wrong side of the decision boundary. In order to turn the decision boundary in the right way, anti-clockwise, one subtracts a small multiple of \mbox{\boldmathx}^{(4)}:
[TABLE]
These two learning rules combine to the learning rule of Rosenblatt [10]:
[TABLE]
For more than one neuron, the rule reads
[TABLE]
This rule is reminiscent of Hebb’s rule (2.9), except that here inputs and outputs are associated with distinct units. Therefore we have instead of . One applies (5.17) iteratively for a sequence of randomly chosen patterns , until the problem is solved. This corresponds to adding a little bit of Hebb’s rule in each iteration. To ensure that the algorithm stops when the problem is solved, one can use the learning rule [1]
[TABLE]
5.3 Gradient descent for linear units
In this Section, the learning algorithm (5.18) is derived in a different way, by minimising an energy function using gradient descent. This requires differentiation, therefore we must choose a differentiable activation function. The simplest choice is a linear activation function, . We set , so that the network computes:
[TABLE]
A neuron with a linear activation function is called a linear unit. The outputs assume continuous values, but not necessarily the targets . For linear units, the classification problem
[TABLE]
has the formal solution
[TABLE]
This can be verified by inserting Equation (5.21) into (5.19). Here is the overlap matrix with elements
[TABLE]
(Section 3.6). For the solution (5.21) to exist, the matrix must be invertible. As mentioned in Section 3.6, this requires that , because otherwise the input-pattern vectors are linearly dependent, and thus also the columns (and rows) of . If the matrix has linearly dependent columns or rows, it cannot be inverted.
Let us assume that the input patterns are linearly independent, so that the solution (5.21) exists. In this case we can find the solution iteratively. To this end one defines the energy function
[TABLE]
This function is non-negative, and it vanishes when all ouputs equal the corresponding targets, for all patterns.
The energy function (5.23) is regarded as a function of the weights , unlike the energy function in Part I which is a function of the state-variables of the neurons. The goal is now to find weights that minimise . If the input patterns are linearly independent, vanishes at the global miminum, corresponding to the desired solution of the problem (Exercise 5.1). Let us use gradient descent to minimise ,
[TABLE]
with learning rate . This is analogous to Equation (4.25), apart from the minus sign. In Section 4.4 the goal was to maximise the target function, here we want to minimise by taking many downhill steps in search of the global minimum. The derivatives in Equation (5.24) are evaluated with the chain rule, together with Equation (4.28) which takes the form
[TABLE]
for asymmetric weights. This yields the weight increments
[TABLE]
This learning rule is very similar to Equation (5.18). One difference is that Equation (5.26) contains a sum over all patterns. It is important to keep in mind also that the activations functions are different, while Equation (5.18) was derived for , the learning rule (5.26) was derived for . An advantage of the rule (5.26) is that it is derived from an energy function.
This helps to analyse the convergence of the algorithm, as we have seen in Chapter 2.
Linear units [Equation (5.19)] are special. The Boolean AND problem (Figure 5.7) does
not admit the solution (5.21), even though the problem is linearly separable. Since the pattern vectors \mbox{\boldmathx}^{(\mu)} are linearly dependent, the solution (5.21) does not exist. Shifting the patterns or introducing a threshold does not change this fact.
In Section 5.5 we discuss how to solve problems that are not linearly separable using a hidden layer of neurons with non-linear activation functions. Note that introducing hidden layers with linear units does not help,
because the resulting input-output mapping is still linear if all neurons have linear activation functions, so that only problems with can be solved. This is the main reason for using hidden layers with non-linear activation functions.
There are four points to keep in mind. First, if the the patterns are linearly independent, then we can use gradient descent to determine suitable weights (and thresholds) of linear units. Second, in general hidden layers with non-linear units are required, because a single neuron with a continuous non-linear and monotonous activation function can only solve problems with linearly independent patterns (Exercise 5.11). Third, for gradient descent for non-linear units we must require that the activation function is differentiable, or at least piecewise differentiable. Fourth, in this case we calculate the gradients using the chain rule, resulting in factors of derivatives . This is the origin of the vanishing-gradient problem (Chapter 7).
5.4 Classification capacity
In Chapter 3 we analysed the storage capacity of Hopfield networks. The analogous question for the classification problem described in Section 5.1 is: how many patterns can a single neuron with activation function classify? As in the case of Hopfield networks, one can find a general answer for random binary classification problems.
Consider points with coordinate vectors \mbox{\boldmathx}^{(\mu)} in -dimensional input space, and assign random targets:
[TABLE]
This random classification problem is homogeneously linearly separable if we can find an -dimensional weight vector , so that \mbox{\boldmathw}\cdot\mbox{\boldmathx}=0 is a valid decision boundary that goes
through the origin:
[TABLE]
So homogeneously linearly separable problems are binary classification problems that are linearly separable by a hyperplane that contains the origin. Problems with this property can be solved by a binary threshold unit with threshold .
Now assume that the points (including the origin) are in general position (Figure 5.10). In this case Cover’s theorem [72] gives an expression for the probability that the random binary classification problem of patterns in dimension is homogeneously linearly separable:
[TABLE]
Here are the binomial coefficients, for . Equation (5.29) is proven by recursion, starting from a set of points in general position. Assume that the number of homogeneously linearly separable classification problems given these points is known. After adding one more point, one can compute the in terms of , and recursion yields Equation (5.29). Figure 5.11 shows this result as a function of for different values of . For , any random classification problem is homogeneously linearly separable. In this case the pattern vectors are linearly independent, so that the problem can also be solved by a linear unit (Section 5.3). But a neuron with activation function can classify problems with more than patterns. In the limit of , the function approaches a step function (Exercise 5.12). In this limit the maximal classification capacity is therefore .
What is the expected classification capacity for finite values of ? To answer this question, consider a random sequence of patterns \mbox{\boldmathx}^{(1)},\mbox{\boldmathx}^{(2)},\ldots and targets and ask [72]: what is the distribution of the largest integer so that the problem \mbox{\boldmathx}^{(1)},\mbox{\boldmathx}^{(2)},\ldots,\mbox{\boldmathx}^{(n)} is separable in dimension , but \mbox{\boldmathx}^{(1)},\mbox{\boldmathx}^{(2)},\ldots,\mbox{\boldmathx}^{(n)},\mbox{\boldmathx}^{(n+1)} is not? is the probability that patterns are linearly separable in -dimensional input space. We can write where is the conditional probability that patterns are linearly separable if the patterns are. Then the probability that patterns are not separable (but patterns are) reads . We can interpret the right-hand side of this Equation as a distribution of the random variable , the maximal number of separable patterns in dimension :
[TABLE]
It follows that the expected maximal number of separable patterns is
[TABLE]
So the expected classification capacity is twice the input dimension:
[TABLE]
This quantifies the notion that it is easier to separate patterns in higher-dimensional input space. As an illustration, consider the XOR problem which is not linearly separable in two-dimensional input space. The problem becomes separable when we embed the points in three-dimensional space, for instance by assigning to the patterns and to the patterns (Figure 5.12).
5.5 Multi-layer perceptrons
In Sections 5.1 and 5.2 we discussed how to solve linearly separable problems [Figure 5.13(a)]. The aim of this Section is to show that non-separable problems like the one in Figure 5.13(b) can be solved by a perceptron with one hidden layer. A network that does the trick for the classification problem in Figure 5.13(b) is depicted in Figure 5.14.
As in the previous Section, all neurons have the signum function as activation function, with possible outputs :
[TABLE]
Each of the three neurons in the hidden layer has its own decision boundary.
The idea is to choose the weights and the thresholds in such a way that the three decision boundaries partition the input plane into distinct regions, so that each region contains either only patterns or patterns [3].
How this construction works is shown in Figure 5.15.
The left part of the Figure shows the three decision boundaries with their weight vectors, and how they divide the input plane into different regions which contain either only or only . Each region bears a three-digit code made out of the symbols and . The codes are determined by the states of the hidden neurons. A sign in the -th entry of the code means that . So the region in question is on the weight-vector side of the decision boundary . A sign, by contrast, corresponds to . In this case the region is on the other side of the decision boundary, the one opposite the weight vector. The value table shows the targets associated with each region, together with the code of the region.
The weights and the threshold of the output neuron are chosen so that it associates the correct target value with each region. A graphical representation of the output problem is shown in Figure 5.16. This problem is linearly separable (Exercise 5.3).
The following function computes the correct output for each region:
[TABLE]
This solves the binary classification problem described in Figure 5.15, but note that the solution is not unique. There is a range of different weights and thresholds that solve the problem, and there are other solutions based on different network layouts. Nevertheless, the solution illustrates how non-linearly separable classification problems can be solved by adding a hidden layer to the network layout. The neurons in the hidden layer define segments of a piecewise linear decision boundary. More hidden neurons are needed if the decision boundary is very wiggly.
Figure 5.17
shows another example, how to solve the Boolean XOR problem with a perceptron that has two neurons in a hidden layer, with activation functions , thresholds and , and all weights equal to unity. The output neuron has weights and and unit threshold:
[TABLE]
Minsky and Papert [11] proved in 1969 that all Boolean functions can be represented by multilayer perceptrons, but that at least one hidden neuron must be connected
to all input terminals. This means that not all neurons in the network are locally connected (the neurons have only a few incoming weights). Since fully connected networks are much harder to train than locally connected ones, this was considered a shortcoming at the time. Now, almost 50 years later, the perspective has changed. Convolutional networks (Chapter 8) have only local connections to the inputs and can be trained to recognise objects in images with high accuracy.
In summary, perceptrons are trained on a training set [\mbox{\boldmathx}^{(\mu)},\mbox{\boldmatht}^{(\mu)}], , by moving the decision boundaries into the correct positions. This is achieved by repeatedly applying Hebb’s rule to adjust all weights. A related learning rule is obtained by gradient-descent on the energy function (5.23). Also, we have not discussed how to update the thresholds yet, but it is clear that they too can be updated with gradient-descent learning.
Once all decision boundaries are in the right places we must ask:
what happens if we apply the trained network to a new data set? Does it classify the new inputs correctly? In other words, can the network generalise? An example is shown in Figure 5.18. Panel (a) shows the result of training the network on a training set. The decision boundary separates patterns from patterns, so that the network classifies all patterns in the training set correctly. In panel (b) the trained network is applied to patterns in a validation set. We see that most patterns are correctly classified, save for one error. This means that the energy function (5.23) is not exactly zero for the validation set. Nevertheless, the network does quite a good job. Usually it is not a good idea to try to precisely classify all patterns near the decision boundary, because real-world data sets are subject to noise (Section 6.4).
5.6 Summary
Perceptrons are layered feed-forward networks that can learn to classify data in a training set [\mbox{\boldmathx}^{(\mu)},\mbox{\boldmatht}^{(\mu)}]. For each input pattern \mbox{\boldmathx}^{(\mu)}, the network finds the correct target vector \mbox{\boldmatht}^{(\mu)}. We discussed the learning algorithm for a simple example: real-valued patterns with just two components, and one binary target. This allowed us to represent the classification problem graphically, and to see how linearly separable classification problems can be solved by a simple perceptron. There are three different ways of understanding how the perceptron learns. First, geometrically, the perceptron learn by moving decision boundaries into the correct locations. Second, this can be achieved by repeatedly adding a little bit of Hebb’s rule. Third, these rules are similar to the learning rule derived from gradient descent on the energy function (5.23). Cover’s theorem quantifies the capacity of a simple perceptron to separate patterns with binary targets. Finally we discussed how to solve non-linearly separable classification problems with perceptrons with a hidden layer.
5.7 Further reading
As mentioned in the Introduction, a short account of the history of perceptron research is the review by Kanal [23]. The remarkable book by Minsky and Papert explains the geometry of perceptron learning in great depth, and in a very elegant fashion. For a proof of Cover’s theorem see Ref. [73].
Chapter 6 Stochastic gradient descent
In Chapter 5 we discussed how a hidden layer helps to classify problems that are not linearly separable. We explained how the decision boundary
in Figure 5.15 is represented in terms of the weights and thresholds of the hidden neurons, and introduced a training algorithm based on gradient descent. In this Section, the training algorithm is discussed in more detail.
Figure 5.2 shows the layout of the network to be trained. There are input patterns \mbox{\boldmathx}^{(\mu)} with components each, as before. The output of the network has components:
[TABLE]
to be matched to the target vector \mbox{\boldmatht}^{(\mu)}. The network shown in Figure 5.2 computes
[TABLE]
Equation (6.2) shows that the outputs are obtained in terms of nested activation functions. They must be differentiable (or at least piecewise differentiable). Apart from that there is no need to specify them further at this point.
6.1 Chain rule and error backpropagation
The network in Figure 5.2 is trained by gradient-descent learning in the same way as in Section 5.3. The weight increments are given by:
[TABLE]
with energy function
[TABLE]
The small parameter in Equation (6.3) is the learning rate, as in Section 5.3. The derivatives of the energy function are evaluated with the chain rule. For the weights connecting to the output layer we apply the chain rule once
[TABLE]
Here is the derivative of the activation function with respect to the local field , and is the Kronecker delta: =1 if and zero otherwise.
An important point is that the states of the neurons in the hidden layer do not depend on , because these neurons do not have incoming connections with these weights, a consequence of the feed-forward layout of the network. In summary we obtain for the increments of the weights connecting to the output layer:
[TABLE]
is a weighted output error : it vanishes when . The weights connecting to the hidden layer are adjusted in a similar fashion, by applying the chain rule four times:
[TABLE]
Here we used Equation (5.25). With the definition of the output error, , Equation (6.3) yields:
[TABLE]
The last equality defines weighted errors,
[TABLE]
associated with the hidden layer. Note that the vanish when the output errors are zero. Equation (6.9) shows that the errors are determined recursively. The neuron states are also updated recursively, Equation (6.2), but there is an important difference between Equations (6.9) and (6.2). The feed-forward structure of the layered network implies that the neurons are updated from left to right. Equation (6.9), by contrast, says that the errors are updated from right to left, from the output layer to the hidden layer. The term backpropagation refers to this difference: the neurons are updated forward, the errors are updated backwards.
In terms of the errors and ,
the weight increments have the same form for both layers:
[TABLE]
The rule (6.10) is also called -rule [1]. The thresholds are adjusted in a similar way:
[TABLE]
So, the general form for the threshold increments is analogous to Equation (6.10)
[TABLE]
but without the state variables of the neurons (or the inputs). A way to remember the difference between Equations (6.10) and (6.12) is to note that the formula for the threshold increments looks like the one for the weight increments if one sets the state values of the neurons to . This follows from Equation (6.2).
The backpropagation rules (6.10) and (6.12) contain sums over patterns. This corresponds to feeding all patterns at the same time to compute the increments of weights and thresholds (batch training). Alternatively one may choose a single pattern, update the weights by backpropagation, and then continue to iterate these training steps many times (sequential training). One iteration corresponds to feeding a single pattern, iterations are called one epoch (in batch training, one iteration corresponds to one epoch). If one chooses the patterns randomly, then sequential training results in stochastic gradient descent:
[TABLE]
Since the sum over pattern is absent, the steps do not necessarily decrease the energy function. Their directions fluctuate, but the average weight increment (averaged over all patterns) points downhill. The result is a stochastic path through parameter space, less prone to getting stuck in local minima (but see Section 7.8).
6.2 Stochastic gradient-descent algorithm
The stochastic-gradient descent formulae derived in the previous Section were derived for a network with one hidden layer. This Section describes the details of the stochastic-gradient algorithm for deep networks with many hidden layers. To this end we need to adapt our notation, as described in Figure 6.1. We label the layers by the index . The layer of input terminals has label , while layer denotes the layer of output neurons. The state variables for the neurons in layer are , the weights connecting into these neurons from the left are , the errors associated with layer are denoted by . In this notation, Equations (6.2) read:
[TABLE]
Repeating the steps outlined in the previous Section, we arrive at the update formulae
[TABLE]
with errors
[TABLE]
where is the local field of . It involves the matrix-vector product between the weight matrix and the vector \mbox{\boldmathV}^{(\ell-1)}. Evaluating the gradients with the chain rule, one obtains the recursion
[TABLE]
with initial condition . For one hidden layer, Equation (6.17) is equivalent to (6.9). The result of the recursion (6.17) is a vector \mbox{\boldmath\delta}^{(\ell-1)} with components , obtained by component-wise multiplication of [{\mathbb{W}^{(\ell)}}^{\sf T}\mbox{\boldmath\delta}^{(\ell)}]_{j} with . Component-wise multiplication of vectors is sometimes called Schur or Hadamard product [75], denoted by \mbox{\boldmatha}\odot\mbox{\boldmathb}=[a_{1}b_{1},\ldots,a_{N}b_{N}]^{\small\sf T}. It does not have a geometric meaning like the scalar product or the cross product of vectors, and therefore there is little point in using it. Also, note that the vector \mbox{\boldmath\delta}^{(\ell)} is multiplied by the transpose of the weight matrix, , rather than by the weight matrix itself. We return to this point in Section 9.1.
The stochastic-gradient algorithm is summarised in Algorithm 4. One feeds an input \mbox{\boldmathx}^{(\nu)}, updates the weights using (6.15), and iterates these steps until the energy function (5.23) is deemed sufficiently small. Note that the resulting weights and thresholds are not unique. In Figure 5.17 all weights for the Boolean XOR function are equal to . But the training algorithm
(6.10) corresponds to repeatedly adding weight increments. This may cause the weights to grow.
In practice, the stochastic gradient-descent dynamics may be too noisy. In this case it is better to average over a small number of randomly chosen patterns. Such a set is called mini batch, of size say. In stochastic gradient descent with mini batches one replaces Equations (6.10) and (6.12) by
[TABLE]
Sometimes the mini-batch rule is quoted with prefactors of before the sums. The factors can just be absorbed in the learning rate, but when comparing learning rates for different implementations one needs to check whether or not there are factors of in front of the sums in Equation (6.18).
How does one select which inputs to include in a mini batch? This is discussed below, in Section 6.3: at the beginning of each epoch, one randomly shuffles the sequence of the input patterns in the training set. Then the first mini batch contains patterns , and so forth.
Common choices for the activation functions are the sigmoid function or tanh:
[TABLE]
In both cases, the derivatives can be expressed in terms of the function itself:
[TABLE]
The second equality was used in Section 3.4. The following short-hand notation for the derivative of the activation function is common: .
As illustrated in Figure 6.2, the activation functions (6.19) saturate at large values of : their derivatives tend to zero. Since the backpropagation rule (6.16) contains factors of , this implies that the algorithm slows down if becomes too large. For this reason, the initial weights and thresholds should be chosen so that the local fields are not too large in magnitude, to avoid that becomes too small. A standard procedure is to take all weights to be initially randomly distributed, for example Gaussian with zero mean, and with a suitable variance. The performance of networks with many hidden layers (deep networks) can be sensitive to the initialisation of the weights (Section 7.2).
It is sometimes argued that the initial values of the thresholds are not so critical. The idea is that they are learned more rapidly than the weights, at least initially, and a common choice is to initialise the thresholds to zero. Section 7.2 summarises a mean-field argument that comes to a different conclusion.
6.3 Preprocessing the input data
It can be useful to preprocess the input data, although any preprocessing may remove information from the data. Nevertheless, it is usually advisable to shift the data so the mean of each component over all patterns vanishes:
[TABLE]
There are several reasons for this. First, large mean values can cause large gradients in the energy function (Exercise 6.9) that are difficult to navigate with gradient descent. Different input-data variances in different directions have a similar effect. Therefore one scales the inputs so that the input-data distribution has the same variance in all directions (Figure 6.3), equal to unity for instance:
[TABLE]
Second, to avoid that the neurons connected to the inputs saturate, their local fields must not be too large (Section 6.2). If one initialises the weights to Gaussian random numbers with mean zero and unit variance, large activations are quite likely if the distribution of input patterns has a large mean or a large variance. Third, enforcing zero input mean by shifting the input data avoids that the weights of the neurons in the first hidden layer must decrease or increase together [76]. Equation (6.18) shows that the components of \delta\!\mbox{\boldmathw}_{m}\propto\delta_{m}\mbox{\boldmathx} into hidden neuron are likely to have the same signs if the input data has a large mean. This makes it difficult for the network to learn to differentiate. In summary, it is advisable to shift and scale the input-data distribution so that it has mean zero and unit variance, as illustrated in Figure 6.3.
The same transformation (using the mean values and scaling factors determined for the training set) should be applied to any new data set that the network is supposed to classify after it has been trained on the training set.
Figure 6.4 shows a distribution of inputs that falls into two distinct clusters. The difference between the clusters is sometimes called covariate shift, here covariate is just another term for input. Imagine feeding just inputs from one of the clusters to the network. It will learn local properties of the decision boundary, instead of its global features. Such global properties are efficiently learned if the network is more frequently confronted with unfamiliar data. For sequential training (stochastic gradient descent) this is not a problem, because the sequence of input patterns presented to the network is random. However, if one trains with mini batches, the mini batches should contain randomly chosen patterns in order to avoid covariate shifts. To this end one randomly shuffles the sequence of the input patterns in the training set, at the beginning of each epoch.
It is also recommended [76] to observe the output errors during training. If the errors are similar for a number of subsequent learning steps, the corresponding inputs appear familiar to the network. Larger errors correspond to unfamiliar inputs, and Ref. [76] suggests to feed such inputs more often.
When the input data is very high dimensional, many input terminals are needed. This usually means that one should use many neurons in the hidden layers. This can be problematic because it increases the risk of overfitting the input data. To avoid this as far as possible, one can reduce the dimensionality of the input data by principal-component analysis. This method allows to project high-dimensional data to a lower dimensional subspace (Figure 6.5).
The data shown on the left of Figure 6.5 falls approximately onto a straight line, the principal direction \mbox{\boldmathu}_{1}. We see that the coordinate orthogonal to the principal direction is not useful in classifying the data. Consequently this coordinate can be disregarded, reducing the dimensionality of the data set. The idea of principal component analysis is to rotate the basis in input space so that the variance of the data along the first axis of the new coordinate system, \mbox{\boldmathu}_{1}, is maximal. One keeps the input components corresponding to \mbox{\boldmathu}_{1}, discarding those corresponding to \mbox{\boldmathu}_{2} (Figure 6.5).
To determine the maximal-variance direction, consider the data variance along a unit direction vector (|\mbox{\boldmathv}|=1):
[TABLE]
Here
[TABLE]
is the data covariance matrix. The variance is maximal when points in the direction of the leading eigenvector of the covariance matrix . This can be seen as follows. The covariance matrix is symmetric, therefore its eigenvectors \mbox{\boldmathu}_{1},\ldots,\mbox{\boldmathu}_{N} form an orthonormal basis of input space. This allows us to express the matrix as
[TABLE]
The eigenvalues are non-negative. This follows from Equation (6.24) and the eigenvalue equation \mathbb{C}\mbox{\boldmathu}_{\alpha}=\lambda_{\alpha}\mbox{\boldmathu}_{\alpha}. We arrange the eigenvalues by magnitude, . Using Equation (6.25) we can write for the variance
[TABLE]
with v_{\alpha}=\mbox{\boldmathv}\cdot\mbox{\boldmathu}_{\alpha}. We want to show that is maximal for \mbox{\boldmathv}=\pm\mbox{\boldmathu}_{1} subject to the constraint that is normalised to unity,
[TABLE]
To ensure that this constraint is satisfied as the are varied, one introduces a Lagrange multiplier (Exercises 6.10 and 6.11). The constraint (6.27) is multiplied with and added to the target function (6.26). The function to maximise reads
[TABLE]
To find the maximum of , we determine its singular points, defined by . This yields . The maximum of is obtained for , where is the maximal eigenvalue of with eigenvector \mbox{\boldmathu}_{1}. We conclude that all components must vanish, except one which must equal unity. This shows that the variance is maximised by the principal direction.
In more than two dimensions there is commonly more than one direction along which the data varies significantly. These principal directions correspond to the eigenvectors of with the largest eigenvalues. This can be shown recursively. One projects the data to the subspace orthogonal to \mbox{\boldmathu}_{1} by applying the projection matrix \mathbb{P}_{1}=\mathbb{1}-\mbox{\boldmathu}_{1}\mbox{\boldmathu}_{1}^{\sf T}. Then one repeats the procedure outlined above, and finds that the data varies maximally along \mbox{\boldmathu}_{2}. Upon iteration, one obtains the principal directions \mbox{\boldmathu}_{1},\ldots,\mbox{\boldmathu}_{k}. Often there is a gap between the largest eigenvalues and the small ones (all close to zero). Then one can safely project the data onto the subspace spanned by the principal directions. If there is no gap then it is less clear what to do.
The data set shown on the right of Figure 6.5 illustrates another problem. This data set is much harder to classify if we use only the principal component alone. In this case we lose important information by projecting the data on its principal component.
6.4 Overfitting and cross validation
The goal of supervised learning is to generalise from a training set to new data. Only general properties of the training set are of interest, not specific ones that are particular to the training set in question. A neural network with more neurons may classify the input data better, because it more accurately represents all specific features of the given data set. But a different set of patterns from the same input distribution can look quite different in detail, in which case the decision boundary may not classify the new data very well (Figure 6.6). In other words, the network may fit too fine details (for instance noise in the training set) that have no general meaning. This problem, illustrated in Figure 6.6, is called overfitting. The tendency to overfit is larger for networks with more neurons. In general, we should look for a compromise, reducing the tendency of the network to overfit at the expense of training accuracy.
One way of avoiding overfitting is to use cross validation and early stopping. One splits the data into two sets: a training set and a validation set. The idea is that these sets share the general features to be learnt. But although training and validation sets are drawn from the same distribution, they may differ in details that are not of interest.
While the network is trained on the training set, one monitors not only the energy function for the training set, but also the energy function evaluated using the validation data. As long as the network learns general features of the input distribution, both training and validation energies decrease. But when the network starts to learn specific features of the training set, then the validation energy saturates, or may start to increase. At this point the training is stopped. The scheme is illustrated in Figure 6.7.
Often the possible state values of the output neurons are continuous while the targets assume only discrete values. In this case one may also monitor the classification error of the validation set. The definition of the classification error depends on the type of the classification problem. For one single output neuron with targets , the classification error is defined as
[TABLE]
As a third example, consider a classification problem where inputs must be classified into mutually exclusive classes, such as the MNIST data set of hand-written digits (Section 8.3) where . Another example is given in Table 6.1, with classes. In both cases, one of the targets equals unity while all others equal zero. As a consequence, the targets sum to unity: . Now assume that the network has sigmoid outputs, . To classify input \mbox{\boldmathx}^{(\mu)} from the network outputs we define
[TABLE]
In all cases, the classification accuracy is defined as , it is usually quoted in percent.
The classification error determines the fraction of inputs that are classified wrongly. However, it contains less information than the energy function, which is in fact a mean-squared error of the outputs. This is illustrated in Table 6.1. All three inputs are classified correctly, but there is a substantial mean-squared error. This indicates that the classification is not very reliable.
6.5 Adaptation of the learning rate
It is tempting to choose larger learning rates, because they enable the network to escape more efficiently from shallow minima. But this can lead to problems when the energy function varies rapidly, causing the training to fail. To avoid this, one uses an adaptive learning rule, such as:
[TABLE]
Here labels the iteration number. We see that the increment at step depends not only on the instantaneous gradient, but also on the weight change of the previous iteration. We say that the dynamics becomes inertial, the weights gain momentum. The parameter is called momentum constant. It determines how strong the inertial effect is. We see that corresponds to the usual backpropagation rule. When is positive, then how does inertia change the learning process? Iterating Equation (6.31) yields
[TABLE]
Here and in the following we use the short-hand notation
[TABLE]
Equation (6.32) shows that is a weighted average of the gradients encountered during training. Now assume that the training is stuck in a shallow minimum. Then the gradient remains roughly constant through many time steps. To illustrate what happens, let us assume that for . In this case we can write
[TABLE]
In this situation, convergence is accelerated when is close to unity. We also see that it is necessary that for the sum in Equation (6.33) to converge.
The other limit to consider is that the gradient changes rapidly from iteration to iteration. How is the learning rule modified in this case? As an example, let us assume that the gradient remains of the same magnitude, but that its sign oscillates, for . Inserting this into Equation (6.32), we obtain:
[TABLE]
Here the increments are much smaller compared with those in Equation (6.33). This shows that introducing inertia can substantially accelerate convergence without sacrificing accuracy. The disadvantage is, of course, that there is yet another parameter to choose, namely the momentum constant .
Nesterov’s accelerated gradient method [77] is another way of implementing momentum. The algorithm was developed for smooth optimisation problems, but it has been suggested to use the method when training deep neural networks with gradient descent [78]:
[TABLE]
A suitable sequence of coefficients is defined by recursion [78]. The coefficients approach unity from below as increases.
Nesterov’s accelerated-gradient method is more accurate than the simple momentum method, because the accelerated-gradient method evaluates the gradient at an extrapolated point, not at the initial point. Figure 6.8 illustrates a situation where Nesterov’s method converges more rapidly. Nesterov’s method is not much more difficult to implement than Equation (6.31), and it is not much more expensive in terms of computational cost.
There are other ways of adapting the learning rate during training, described in Section 4.10 in Haykin’s book [2]. Finally, the learning rate need not be the same for all neurons. If the weights of neurons in different layers change at very different speeds (Section 7.2), one could define a layer-dependent learning rate that is larger for neurons with smaller gradients.
6.6 Summary
Backpropagation is an efficient algorithm for stochastic gradient-descent on the energy function (6.4) in weight space, because it refers only to quantities that are local to the weight to be updated. Networks with many hidden neurons have many free parameters (their weights and thresholds). This increases the risk of overfitting, which reduces the power of the network to generalise. Deep networks with many hidden layers are particularly prone to overfitting (Chapter 7). The tendency of networks to overfit can be reduced by cross validation.
6.7 Further reading
The backpropagation algorithm is explained in Section 6.1 of Hertz, Krogh and Palmer [1], and in Chapter 4 of Haykin’s book [2]. The paper [76] by LeCun et al. predates deep learning, but it is still a very nice collection of recipes for making backpropagation more efficient.
One of the first papers on error backpropagation is the one by Rumelhart et al. [12] from 1986. The authors provide an elegant explanation and summary of the backpropagation algorithm. They also describe results of different numerical experiments, and one of them introduces convolutional networks (Chapter 8) to learn to tell the difference between the letters T and C (Figure 6.9).
Chapter 7 Deep learning
7.1 How many hidden layers?
In Chapter 5 we saw why it is sometimes necessary to have a hidden layer: in order to solve problems that are not linearly separable. Under which circumstances is one hidden layer sufficient? Are there problems that require more than one hidden layer? Even if not necessary, may additional hidden layers improve the performance of the network?
To understand how many hidden layers suffice, it is useful to view the
classification problem as an approximation problem [79]. Consider the classification problem [\mbox{\boldmathx}^{(\mu)},t^{(\mu)}] for . This problem defines a target function t(\mbox{\boldmathx}). Training a network to solve this task corresponds to approximating the target function t(\mbox{\boldmathx}) by the output function O(\mbox{\boldmathx}) of the network, from -dimensional input space to one-dimensional outputs.
How many hidden layers are necessary or sufficient to approximate a given set of functions to a certain accuracy, by choosing weights and thresholds? The answer depends on the nature of the target function. Is it real-valued, perhaps continuous, or does it assume only discrete values?
We start with real-valued inputs and a single output [80, 1]. Consider the network drawn in Figure 7.2. The neurons in the hidden layers have sigmoid activation functions . The output is a linear unit, with linear activation function .
With two hidden layers one tries to approximate the function t(\mbox{\boldmathx}) by
[TABLE]
In the simplest case the inputs are one-dimensional (Figure 7.3). The training set consists of pairs that encode the target function . The task is then to approximate by the network output :
[TABLE]
To this end one uses linear combinations of the basis functions shown in Figure 7.4. Any reasonable real-valued function can be approximated by sums of such basis functions, each suitably shifted and scaled. Furthermore, these basis functions can be expressed as scaled differences of sigmoid activation functions
[TABLE]
Comparing with Equation (7.1) shows that one hidden layer is sufficient to construct the function in this way.
Now consider two-dimensional inputs. In this case, a suitable basis function is (Figure 7.5):
[TABLE]
So two hidden layers are sufficient for two input dimensions. For each basis function we require four neurons in the first hidden layer and one neuron in the second hidden layer. The construction is analogous for more than two input dimensions. For each basis function we need neurons in the first and one neuron in second layer. In conclusion, two hidden layers are sufficient to approximate a real-valued target function.
Yet it is not always necessary to have two layers for real-valued functions. For continuous functions, one hidden layer is sufficient. This is ensured by the universal approximation theorem [2]. It says any continuous function can be approximated to arbitrary accuracy by a network with a single hidden layer, for sufficiently many neurons in the hidden layer.
In Chapter 5 we considered discrete Boolean functions.
Any Boolean function with -dimensional inputs can be represented by a network with one hidden layer, using neurons in the hidden layer. An example for such a network is discussed in Ref. [1]:
[TABLE]
A difference compared with the Boolean-function representations in Section 5.5 is that here the inputs take the values . The reason is that this simplifies the proof, which is by construction [1]. For each hidden neuron one assigns the weights as follows
[TABLE]
with (see below). The thresholds of all hidden neurons are the same, equal to . The idea is that each input pattern turns on exactly one neuron in the hidden layer, the winning neuron). The weights feeding into the output neuron are determined as follows. If the output for the pattern represented by neuron is , let , otherwise . The threshold is set to .
To show how this construction works, consider the Boolean XOR function as an example (Figure 7.6).
To confirm that only the corresponding winning neuron gives a positive signal, consider pattern \mbox{\boldmathx}^{(1)}=[-1,-1]^{\sf\small T}. It activates the first neuron in the hidden layer (). To see this, compute the local fields of the hidden neurons:
[TABLE]
If we choose then the output of the first hidden neuron gives a positive output (), the other neurons produce negative outputs, for . In conclusion, output neuron is the winning neuron for this pattern. Now consider \mbox{\boldmathx}^{(3)}=[-1,+1]^{\sf\small T}. In this case
[TABLE]
Now the third hidden neuron gives a positive output, while the others yield negative values. It works in the same way for the other two patterns, \mbox{\boldmathx}^{(2)} and \mbox{\boldmathx}^{(4)}. In summary, there is a unique winning neuron for each pattern.111That pattern gives the winning neuron is of no importance, it is just a consequence of how the patterns are ordered in the value table in Figure 7.6 Figure 7.7 shows how the four decision boundaries corresponding to partition the input plane.
According to the scheme outlined above, the output neuron computes
[TABLE]
with . For \mbox{\boldmathx}^{(1)} and \mbox{\boldmathx}^{(4)} we find the correct result . The same is true for \mbox{\boldmathx}^{(2)} and \mbox{\boldmathx}^{(3)}, in this case we obtain . In summary, this example illustrates how an -dimensional Boolean function is represented by a network with one hidden layer, with neurons. The problem is of course that this network is expensive to train for large values of because the number of hidden neurons is very large.
There are more efficient layouts if one uses more than one hidden layer. As an example, consider the parity function for -dimensional binary inputs with bits equal to 0 or 1. It measures the parity of sequences of input bits. The function evaluates to unity if there is an odd number of ones in the input, otherwise to zero. A construction similar to the above yields a network layout with neurons in the hidden layer. If one instead wires together the XOR networks, one can solve the parity problem with neurons [81] (Figure 7.8). When is a power of two, this network has neurons. To see this, set the number of inputs to . Figure 7.8 shows that the number of neurons satisfies the recursion with . The solution of this recursion is .
This example also illustrates a second reason why it may be useful to have more than one hidden layer. To design a neural network for a certain task it is often convenient to build the network from well-studied building blocks. One wires them together, often in a hierarchical fashion. In Figure 7.8 there is only one building block, the XOR network from Figure 5.17. Other examples are convolutional networks for image analysis. Here the fundamental building blocks are so-called feature maps, they recognise different geometrical features in the image, such as edges or corners (Chapter8).
7.2 Vanishing and exploding gradients
This Section describes an inherent instability in the training of deep networks with stochastic gradient descent, the vanishing- or exploding-gradient problem.
In Chapter 6 we saw that learning slows down when the factors in the recursion (6.17) become small. When the network has several hidden layers, like the one shown in Figure 7.9, potentially small factors of are multiplied, aggravating the problem. As a consequence, the weights of hidden neurons close to the input layer change only by small amounts, the smaller the more hidden layers the network has. This is the vanishing-gradient problem.
Figure 7.10 quantifies the problem. The Figure shows that the r.m.s. errors averaged over different realisations of random initial weights, \delta_{\rm rms}^{(\ell)}\equiv\big{(}\langle N^{-1}\sum_{j=1}^{N}[\delta^{(\ell)}_{j}]^{2}\rangle\big{)}^{1/2}, tend to be very small during initial training. To explain this phenomenon, consider the simple example discussed in Ref. [5]: a long chain of neurons with only one neuron per layer (Figure 7.11). The output is given by nested activation functions
[TABLE]
Let us compute the errors using Equation (6.16). The partial derivative in (6.16) is evaluated using the chain rule:
[TABLE]
where is the local field for neuron . This yields the following expression for :
[TABLE]
Inserting this expression into Equation (6.16), we find:
[TABLE]
One can also obtain this result by applying the recursion from Algorithm 4, .
Now consider the early stages of training [5]. For the activation functions (6.19), the maximum of is and , respectively, and becomes exponentially small if is large. If one initialises the weights as described in Section 6.2, to Gaussian random variables with mean zero and variance say, then the factors tend to be smaller than unity. In this case, Equation (7.12) implies that the error or gradient vanishes quickly as decreases. The reason is simply that the number of small factors in the product (7.12) increases when becomes smaller, and multiplying many small numbers gives a very small product. As a result, the training slows down. As mentioned above, this is the vanishing-gradient problem (phase I in Figure 7.10).
What happens at later times? Figure 7.10 indicates that the network continues to learn slowly. For the particular example shown in Figure 7.10, the effect persists for about 20 epochs. Then the first layers begin to learn faster (phase II). There is to date no mathematical theory describing how this transition occurs. Much later in training, the errors decay as the learning converges (phase III in Figure 7.10).
There is a second, equivalent, point of view [5]: the learning is slow in a layer far from the output because the output is not very sensitive to the state of these neurons. The effect of a given neuron on the output is measured by Equation (7.12), which describes how the output of the network changes when changing the state of a neuron in a particular layer. At any rate, Equation (7.13) demonstrates that hidden layers far from the output learn slowly, at least initially when the weights are still random.
Suppose we try to combat the vanishing-gradient problem by increasing the weight variance . The problem is that this may cause the factors to become larger than unity. As a consequence, the gradients increase exponentially instead (exploding gradients). In conclusion, the training dynamics is fundamentally unstable. This is due to the multiplicative nature of the recursion for the errors. Taking the logarithm of the product in Equation (7.12) and assuming that the weights are independently distributed random numbers, the central-limit theorem (Chapter 2) implies that the distribution of is Gaussian. In other words, the distribution of the errors is lognormal, implying that very small and very large values of occur with high probability.
In networks like the one shown in Figure 7.9 the principle is the same. Assume that all layers have the same number of neurons. When , one multiplies matrices, instead of numbers. The product (7.12) of random numbers becomes a product of random matrices. Using the chain rule we find:
[TABLE]
With the update rule
[TABLE]
we can evaluate each factor:
[TABLE]
Substituting this result into Equation (7.14), we see that the partial derivatives can be computed in the form of a matrix product. The matrix with elements is given by:
[TABLE]
Here is the matrix of weights feeding into layer , and
[TABLE]
is the diagonal matrix with entries . The matrix product (7.17) determines the error dynamics, just
like Equation (7.13):
[TABLE]
Does the magnitude of \big{|}\mbox{\boldmath\delta}^{(\ell)}\big{|}^{2}={\mbox{\boldmath\delta}^{(\ell)}}^{\sf T}\mbox{\boldmath\delta}^{(\ell)} of the errors shrink or grow as they propagate through the layers? This is determined by the eigenvalues of the left Cauchy-Green matrix
, with . This matrix is symmetric, and its eigenvalues are non-negative. Their square roots are the singular values of :
[TABLE]
It is customary to sort the singular values by their magnitudes, as in Equation (7.20). When there are many layers, the number of factors in Equation (7.17) is large. In this case, the maximal singular value either decreases or increases exponentially as a function of [82]. The corresponding rate
[TABLE]
is called the maximal Lyapunov exponent. A negative maximal Lyapunov exponent indicates that the errors vanish exponentially. The eigenvectors of are called backward Lyapunov vectors. They describe how the errors change as they propagate through the network.
How small differences between the inputs change, is determined by the forward Lyapunov vectors, the eigenvectors of . Since and have the same eigenvalues, the rate of decay or increase of the magnitude of input differences is the same as that of the error magnitudes.
The concept of a maximal Lyapunov exponent is borrowed from chaos theory [83, 84, 85], where implies that small perturbations of the initial conditions grow exponentially as a function of time. The iterated map (7.10) is a dynamical system.
The transition in Figure 7.10 is triggered by a change of the Lyapunov exponent from negative values to [86]. In summary, the unstable-gradient problem in deep networks is due to the fact that the maximal singular value either increases or decreases exponentially as one moves away from the output layer, depending on whether the maximal Lyapunov exponent is negative or positive.
Pennington et al. [87] suggested to combat the unstable-gradient problem by initialising weights and thresholds in such a way that the maximal Lyapunov exponent is close to zero, in order to make sure that the errors neither grow nor shrink exponentially. Consider the network shown in Figure 7.9, with neurons per hidden layer, and initialise the weights to independent Gaussian random numbers with mean zero and variance . The thresholds are initialised in the same way, with variance . In the limit of one can use a mean-field theory [87], just as in Chapter 3, to estimate the maximal Lyapunov exponent.
Following Ref. [87], the first step is to compute how the errors propagate through the network. We assume uncorrelated random input patterns, Equation (2.29), and random weights with mean zero and variance . When , the errors are sums of many random numbers [Equation (6.17)]. Invoking the central-limit theorem (Chapter 2), one concludes that the errors are approximately Gaussian distributed, with mean zero and with variance
[TABLE]
The last approximation neglects possible correlations with the local fields. The variance does not depend on , so that the sum just gives a factor of .
In the limit of large , the central-limit theorem ensures that the local fields are Gaussian distributed too, with mean zero and variance
[TABLE]
This allows us to estimate
[TABLE]
Equation (7.23) describes how the distribution of local fields narrows or broadens as one iterates. For , it was shown in Ref. [87] that approaches a fixed point, , under certain conditions on the variances of the weights and thresholds, and . If is well approximated by , Equation (7.22) simplifies to: . This results in a mean-field estimate of the the maximal Lyapunov exponent,
[TABLE]
The network parameters should be adjusted so that this exponent is as close to zero as possible. This means, in particular, that one should take
[TABLE]
see also Refs. [88, 89]. But we must keep in mind that Equation (7.25) relies on taking the limit . It is expected that the assumptions underlying Equation (7.25) break down when is finite, causing the mean-field theory to fail. The tails of the error distribution, for example, are expected to become heavier as decreases, as indicated by the results for described above.
Note also that assumes rank zero with a small but non-zero probability when is finite. In this case .
There are a number of other tricks that help to cope with unstable gradients in practice, to some extent at least. First, it is sometimes argued that activation functions which do not saturate at large , such as the ReLU function , help against the vanishing-gradient problem. Second, batch normalisation (Section 7.6.5) may reduce the unstable-gradient problem. Third, introducing connections that skip layers (residual networks) can alleviate the unstable-gradient problem. This is discussed in Section 7.4.
Finally, there is an important aspect of the problem that we did not discuss: unstable gradients limit the extent to which information can propagate through the network in a meaningful way. This is explained in Ref. [90].
In this Section we assumed all along that the weights are random numbers. When the network starts to learn, this is no longer the case. The question is how the singular values of change when correlations between different factors in the product (7.17) develop.
7.3 Rectified linear units
Glorot et al. [91] suggested to use a piecewise activation function, the ReLU function222Since the derivative of the ReLU function is discontinuous at , a common convention is to set the derivative to zero at . (Chapter 1). What is the point of using ReLU neurons? When training
a deep network with ReLU activation functions, many of the hidden neurons produce output zero. This means that the network of active neurons (non-zero output) is sparsely connected. It is sometimes argued that sparse networks have desirable properties; at least sparse representations of a classification problem tend to be easier to learn because they are more likely to be linearly separable (Section 5.4). Figure 7.12 illustrates that for a given input pattern, only a certain fraction of hidden neurons is active. For these neurons the computation is linear, yet different input patterns give different sets of active neurons. The product in Equation (7.17) acquires a particularly simple structure: the matrices are diagonal with 0/1 entries. But while the weight matrices are independent initially, they become correlated as the training proceeds. Also the -matrices develop correlations: which elements vanish depends on which pattern is clamped to the input terminals.
A hidden layer with only one or very few active neurons might act as a bottleneck preventing efficient backpropagation of output errors which could in principle slow down training. For the examples given in Ref. [91], this does not occur. To describe information propagation through a network with ReLU neurons, one should compute the probability that a given number of singular values of the matrix vanish (Section 7.2).
The ReLU function is unbounded for large positive local fields. Therefore, the vanishing-gradient problem (Section 7.2) is thought to be less severe in networks made of rectified linear units. However, since the ReLU function does not saturate, the weights tend to increase. Glorot et al. [91] suggested to use -regularisation (Section 7.6.1) to make sure that the weights do not grow.
Finally, using ReLU functions instead of sigmoid functions speeds up the training,
because the ReLU function has piecewise constant derivatives. Such functions are faster to evaluate than non-linear activation functions and their derivatives.
7.4 Residual networks
One way of reducing the vanishing-gradient problem is to introduce short cuts, connections that skip layers [92]. Empirical evidence shows that networks with such short cuts are easier to train than standard multilayer perceptrons. The likely reason is that the vanishing-gradient problem is less severe in networks with short cuts, because error propagation in such networks is determined by the matrix product with the smallest number of factors.
This Section explains how to train networks with short cuts [93].
The layout is illustrated schematically in Figure 7.13. Black arrows stand for usual feed-forward connections, and the gray arrow indicates a connection that skips a layer. The notation in Figure 7.13 differs somewhat from that of Algorithm 4. The weights from layer to are denoted by , and those from layer to by (gray arrow in Figure 7.13). Note that the superscripts are ordered in the same way as the subscripts: the right index refers to the layer on the left. According to Figure 7.13, neuron in layer computes
[TABLE]
As usual, the argument of the activation function is the local field . The weights of all connections are trained in the usual fashion, by stochastic gradient descent.
To illustrate the structure of the resulting formulae, consider a chain of neurons, just one neuron per layer, with short cuts that skip one neuron (Figure 7.14). We calculate the weight increments using Equations (6.15) and (6.16). The recursion (6.17) applies only to standard feed-forward networks without skipping layers. In order to determine how to update the weights for the network shown in Figure 7.14, we need to evaluate the gradients . To begin with, consider the learning rule for . Using Equations (6.15) and (6.16) one finds
[TABLE]
as in Algorithm 4. In the same way one obtains
[TABLE]
Now consider the learning rule for . Using gives
[TABLE]
as before. But the update for is different, because now the short cuts come into play. The connection from layer to gives rise to an extra term:
[TABLE]
Evaluating the partial derivatives yields:
[TABLE]
Iterating further in this way, one finds the following error-backpropagation rule:
[TABLE]
with for . for . The first term is the same as in the error recursion in Algorithm 4. The second term is due to the skipping connections. These connections reduce the vanishing-gradient problem. To see this, note that we can write the error as
[TABLE]
where the sum is over all paths back through the network. The structure of the general formula, for networks with more than only one neuron per layer, is analogous to Equation (7.34). According to this equation, the smallest errors, or gradients, in networks with many layers are dominated by the product corresponding to the path with the smallest number of steps (factors). Therefore short cuts tend to increase small gradients.
Finally, the network described in Ref. [92] used unit weights for the skipping connections. In this case, the local field of takes the form
[TABLE]
assuming that the hidden layers have the same number of neurons. Here is a residual contribution to the local field (when , the inputs are passed right through to ). Therefore such networks are called residual networks [92]. But note that the networks described in Ref. [92] use convolution layers (Section 8.1).
7.5 Outputs and energy functions
Up to now we discussed networks that have the same activation functions for all neurons in all layers, either sigmoid or activation functions
[Equation (6.19)], or ReLU functions (Sections 1.3 and 7.3). In the output layer, one often uses neurons with a different activation function, so-called softmax outputs:
[TABLE]
Here are the local fields in the output layer. In the limit , we see that where is the index of the winning output neuron, the one with the largest value (Chapter 10.1). For , Equation (7.36) is a soft version of this maximum criterion, thus the name softmax. We set to unity from now on.
Two important properties of softmax outputs are, first, that . Second, the values of the outputs sum to unity,
[TABLE]
Therefore the outputs of softmax units can be interpreted as probabilities. Consider classification problems where the inputs must be assigned to one of classes. In this case, the output of softmax unit is assumed to represent the probability that the input \mbox{\boldmathx}^{(\mu)} is in class (in terms of the targets: while for ). If , we assume that the network is quite certain that input \mbox{\boldmathx}^{(\mu)} is in class . On the other hand, if all , we interpret the network output as uncertain. But note that neural networks may fail like humans sometimes do: their output can be very certain yet wrong (Section 8.6).
Softmax units are used in conjunction with a different energy function,
[TABLE]
Here and in the following stands for the natural logarithm. The function (7.38) is minimal when (Exercise 7.5). To find the correct backpropagation formula for the energy function (7.38), we need to evaluate
[TABLE]
Here the labels denoting the output layer were omitted, and in the following equations the index that refers to the input pattern is dropped as well. Using the identities
[TABLE]
one obtains
[TABLE]
So
[TABLE]
since for the type of classification problem where each input belongs to precisely one class. The corresponding learning rule for the thresholds reads
[TABLE]
Equations (7.42) and (7.43) highlight a further advantage of softmax output neurons (apart from the fact that they allow the output to be interpreted in terms of probabilities). The weight and threshold increments for the output layer derived in Section 6 [Equations (6.6a) and (6.11a)] contain factors of derivatives . As noted earlier, these derivatives tend to zero when the activation function saturates, slowing down the learning. But here the rate at which the neuron learns is simply proportional to the output error, , without any possibly small factor . Softmax units are normally only used in the output layer, because the learning speedup is coupled to the use of the energy function (7.38), and because it is customary to avoid dependencies between the neurons within a hidden layer.
There is an alternative form of the energy function that is very similar to the above, but works with sigmoid activation functions and 0/1 targets.
Instead of Equation (7.38) one chooses:
[TABLE]
with , , and where denotes the sigmoid function (6.19a). To compute the weight increments, we apply the chain rule:
[TABLE]
Using Equation (6.20) we obtain
[TABLE]
identical to Equation (7.42). The thresholds are adjusted in an analogous fashion, Equation (7.43). But now the interpretation of the outputs is slightly different, since the values of the softmax units in the output layers sum to unity, while those with sigmoid activation functions do not. In either case one can use the definition (6.30) for the classification error.
To conclude this Section, we briefly discuss the meaning of the energy functions (7.38) and (7.44). In Chapter 6 we saw how deep neural networks are trained to fit input-output functions (Section 7.1) by minimising the quadratic energy function (6.4). This reminds of regression analysis in mathematical statistics, where the predictive accuracy of a model is improved by minimising the sum over the squared errors. Now consider the energy function (7.44) for a single sigmoid output with targets and . In this case, the network output is interpreted as the probability O^{(\mu)}=\mbox{Prob}(t^{(\mu)}=1|\mbox{\boldmathx}^{(\mu)}) of observing .
The corresponding likelihood is the joint probability of observing the outcomes for independent inputs \mbox{\boldmathx}^{(\mu)}:
[TABLE]
under the model determined by the weights and thresholds of the network.
Minimising the negative log-likelihood (Section 4.4) corresponds to minimising (7.44). This is just binary logistic
regression [94] to predict a binary outcome or . The case corresponds to a multivariate regression problem [94] with possibly correlated outcome variables .
When the targets describe mutually exclusive categorical outcomes, with , the softmax output is interpreted as the probability of observing . An example is the problem of classifying hand-written digits (Section 8.3). Training the network then corresponds to multinomial regression [94] with the log-likelihood (7.38). Note that Equation (7.44), for , is equivalent to (7.38) for , because and . At any rate, these remarks motivate why the energy functions (7.38) and (7.44) are sometimes called log-likelihoods. Equation (7.44) is also referred to as cross entropy, because it has the same for as the cross entropy [65] characterising the difference between two Bernoulli distributions: the network output , and the target .
7.6 Regularisation
Deeper networks have more neurons, so the problem of overfitting (Figure 6.6) tends to be more severe for deeper networks. Regularisation schemes limit the tendency to overfit. Apart from cross validation (Section 6.4), a number of other regularisation schemes have proved useful for deep networks: weight decay, pruning, drop out, expansion of the training set, and batch normalisation . This Section summarises the most important aspects of these methods.
7.6.1 Weight decay
Recall Figure 5.17 which shows a solution of the classification problem defined by the Boolean XOR function.
In the solution illustrated in this Figure, all weights equal , and also the thresholds are of order unity. If one uses the backpropagation algorithm to find a solution to this problem, one may find that the weights continue to grow during training. As mentioned above, this can be problematic because it may imply that the local fields become so large that the activation functions saturate. Then training slows down, as explained in Section 7.2.
To prevent the weights from growing, one can reduce them by some factor during training, either at each iteration or in regular intervals, for , or
[TABLE]
This is achieved by adding a term to the energy function:
[TABLE]
Gradient descent on gives:
[TABLE]
with . One can include a corresponding term for the thresholds. The scheme summarised here is sometimes called -regularisation. An alternative scheme is -regularisation. It amounts to
[TABLE]
This gives the learning rule
[TABLE]
The discontinuity of the learning rule at is cured by defining . Comparing Equations (7.50) and (7.52), we see that -regularisation puts more weights to zero, compared with the -scheme [5].
An alternative to these two methods is max-norm regularisation [95], where the weights are constrained to remain smaller than a given constant: . If a exceeds the positive constant , then is rescaled so that .
These weight-decay schemes are referred to as regularisation schemes because they tend to help against overfitting. How does this work? Weight decay adds a constraint to the problem of minimising the energy function. The result is a compromise [5] between a small value of and small weight values. The idea is that a network with smaller weights is more robust to the effect of noise. When the weights are small, then small changes in some of the patterns do not give a substantially different training result. When the network has large weights, by contrast, it may happen that small changes in the input yield significant differences in the training result that are difficult to generalise.
7.6.2 Pruning
The term pruning refers to removing unnecessary weights or neurons from the network, to improve its efficiency. The simplest approach is weight elimination by weight decay [96]. Weights that tend to remain very close to zero during training are removed by setting them to zero and not updating them anymore. Neurons that have zero weights for all incoming connections are effectively removed (pruned). Pruning is a regularisation method: by removing unnecessary weights, one reduces the risk of overfitting. As opposed to drop out (Section 7.6.3), where hidden neurons are only temporarily ignored, pruning refers to permanently removing hidden neurons. The idea is to train a large network, and then to prune a large fraction of neurons to obtain a much smaller network. It is usually found that such pruned networks generalise better than small networks that were trained without pruning. Up to of the hidden neurons can be removed in some cases. In general, pruning is an excellent way to create efficient classifiers for real-time applications.
An efficient pruning algorithm is based on the idea to remove weights in such a way that the effect upon the energy function is as small as possible [97]. The idea is to find the optimal weight, to remove it, and to change the other weights in such a way that the energy function increases as little as possible. The algorithm works as follows. Assume that the network was trained, so that it reached a (local) minimum of the energy function . One expands the energy function around this minimum. To second order, the expansion of reads:
[TABLE]
The term linear in \delta\!\mbox{\boldmathw} vanishes because we expand around a local minimum. The matrix is the Hessian, the matrix of second derivatives of the energy function.
For the next step it is convenient to adopt the following notation [97]. One groups all weights in the network into a long weight vector (as opposed to grouping them into a weight matrix as we did in Chapter 2). A particular component is extracted from the vector as follows:
[TABLE]
Here is the Cartesian unit vector in the direction , with components . In this notation, the elements of are . Now, eliminating the weight amounts to setting
[TABLE]
To minimise the damage to the network we should eliminate the weight that has least effect upon , changing the other weights at the same time so that increases as little as possible (Figure 7.15). This is achieved by minimising
[TABLE]
The constant term was dropped because it does not matter. Now we first minimise w.r.t. \delta\!\mbox{\boldmathw}, for a given value of . The linear constraint is incorporated using a Lagrange multiplier as in Section 6.3, to form the Lagrangian
[TABLE]
A necessary condition for a minimum [\delta\!\mbox{\boldmathw},\lambda] satisfying the constraint is
[TABLE]
We denote the solution of these Equations by \delta\!\mbox{\boldmathw}^{\ast} and . It is obtained by solving the linear system
[TABLE]
If is invertible, then the top rows of Eq. (7.59) give
[TABLE]
Inserting this result into \hat{\textbf{e}}_{q}^{\sf T}\delta\!\mbox{\boldmathw}^{\ast}+w_{q}=0 we find
[TABLE]
We see that \hat{\textbf{e}}_{q}\cdot\delta\!\mbox{\boldmathw}^{\ast}=-w_{q}, so that the weight is eliminated. The other weights are also changed (black arrow in Figure 7.15). The final step is to find the optimal by minimising
[TABLE]
The Hessian of the energy function is expensive to evaluate, and so is the inverse of this matrix. Usually one resorts to an approximate expression for [97]. One possibility is to set the off-diagonal elements of to zero [98]. But in this case the other weights are not adjusted, because \hat{\textbf{e}}_{q^{\prime}}\cdot\delta\!\mbox{\boldmathw}_{q}^{\ast}=0 for if is diagonal. In this case it is necessary to retrain the network after weight elimination.
The algorithm is summarised in Algorithm 5. It succeeds better than elimination by weight decay in removing the unnecessary weights in the network [97]. Weight decay eliminates the smallest weights. One obtains weight elimination of the smallest weights by substituting in the algorithm described above [Equation (7.62)]. Since small weights are often needed to achieve a small training error, this is usually not a good approximation.
To illustrate the effect of pruning for neural networks with hidden layers, consider the XOR function. Recall that it can be represented by a hidden layer with two neurons (Figure 5.17). For random initial weights, backpropagation takes a long time to find a valid solution, and networks with many more hidden neurons tend to perform much better [99]. The numerical experiments in Ref. [99] indicate that with two hidden neurons, only about 49% of the networks learned the task in 10 000 training steps of stochastic gradient descent, and networks with more neurons in the hidden layer learn more easily (98.5 % for hidden neurons). The data in Ref. [99] also shows that pruned networks, initially trained with hidden neurons, still show excellent training success (83.3 % if only hidden neurons remain). The networks were pruned iteratively during training, removing the neurons with the largest average magnitude. After training, the weights and threshold were reset to their initial values, the values before training began.
One can draw three conclusions from the numerical experiments described in Ref. [99]. First, iterative pruning during training singles out neurons in the hidden layer that had initial weights and thresholds resulting in the correct decision boundaries (lottery-ticket effect [99]). Second, the pruned network with two hidden neurons has much better training success than the network that was trained with only two hidden neurons. Third, despite pruning more than 50% of the hidden neurons, the network with hidden neurons performs almost as well as the one with hidden neurons (97.9 % training success). When training deep networks it is common to start with many neurons in the hidden layers, and to prune up to 90% of them. This results in small trained networks that can classify efficiently and reliably.
7.6.3 Drop out
In this regularisation scheme, some hidden neurons are ignored during training. In each step of the training algorithm (for each mini batch, or for each individual pattern), one disregards at random a fraction of neurons from each hidden layer by setting their outputs to zero, and by updating only the weights and thresholds of the remaining neurons. For the next step in the training algorithm, the ignored neurons at put back, and another set of hidden neurons is disregarded. Once the training is completed, all hidden neurons are activated, but their outputs are multiplied by , to ensure that the local fields are independent of , on average.
This method is motivated by noting that the performance of machine-learning algorithms is usually improved by combining the results of several learning attempts [5, 95], for instance by separately training networks with different layouts, and averaging over their outputs. However, for deep networks this is computationally very expensive. Drop out is an attempt to achieve the same goal more efficiently. The idea is that dropout corresponds to effectively training a large number of different networks. If there are hidden neurons, then there are different combinations of neurons that are turned on or off. The hope is that the network learns more robust features of the input data in this way, and that this reduces overfitting. In practice the method is applied together with max-norm regularisation (Section 7.6.1).
7.6.4 Expanding the training set
If one trains a network with a fixed number of hidden neurons on larger training sets, one observes that the network generalises with higher accuracy (better classification success). The reason is that overfitting is reduced when the training set is larger. Thus, a way of avoiding overfitting is to expand or augment the training set.
It is sometimes argued that the recent success of deep neural networks in image recognition and object recognition is in large part due to larger training sets. One example is ImageNet, a database of more than 107 hand-classified images, into more than 20 000 categories [100].
It is expensive to augment training sets because it requires manual annotation. An alternative is to expand a training set artificially. For digit recognition (Figure 2.1), for example, one could create more input patterns by by shifting, rotating, and shearing the digits, or by adding noise.
7.6.5 Batch normalisation
Batch normalisation [101] can significantly speed up the training of deep networks with backpropagation. The idea is to shift and normalise the input data for each hidden layer, not only for the input patterns (Section 6.3). This is done separately for each mini batch, and for each component of the inputs into the given layer (Algorithm 6). Denoting the states of the neurons feeding into the layer in question by , , one calculates the average and variance over each mini batch
[TABLE]
subtracts the mean from the , and divides by . The parameter is added to the denominator to avoid division by zero when evaluates to zero. There are two additional parameters in Algorithm 6, and . They are are learnt by backpropagation, just like the weights and thresholds. In general the new parameters are allowed to differ from layer to layer, and .
Batch normalisation was originally motivated by arguing that it reduces possible covariate shifts faced by hidden neurons in layer : as the parameters of the neurons in the preceding layer change, their outputs shift thus forcing the neurons in layer to adapt. However in Ref. [102] it was argued that batch normalisation does not reduce the internal covariate shift, but that it speeds up the training by effectively smoothing the energy landscape.
Batch normalisation helps to combat the vanishing-gradient problem because it prevents local fields of hidden neurons to grow. This makes it possible to use sigmoid functions in deep networks, because the distribution of inputs remains normalised. It is sometimes argued that batch normalisation has a regularising effect, and it has been suggested [101] that batch normalisation can replace drop out (Section 7.6.3). It is also argued that batch normalisation may help the network to generalise better, in particular if each mini batch contains randomly picked inputs. Then batch normalisation corresponds to randomly transforming the inputs to each hidden neuron (by the randomly changing means and variances). This may help to make the learning more robust. There is no theory that proves either of these claims, but it is an empirical fact that batch normalisation often speeds up the training.
7.7 Summary
Neural networks with many layers of hidden neurons are called deep networks. Error backpropagation in deep networks suffers from the vanishing-gradient problem. It can be reduced by using ReLU units, by initialising the weights in certain ways, and with networks containing connections that skip layers. Yet vanishing or exploding gradients remain a fundamental difficulty, slowing learning down in the initial phase of training. Nevertheless, deep neural networks have become immensely successful in object recognition, outperforming other algorithms significantly.
Since deep networks contain many free parameters, deep networks tend to overfit the training data, so that the networks must be regularised. Apart from cross validation, there are other ways of regularising the problem: weight decay, drop out, pruning, and data-set augmentation.
7.8 Further reading
Deep networks suffer from catastrophic forgetting: when we train a network on a new input distribution that is quite different from the one the network was originally trained on, then the network tends to forget what it learned initially. A good starting point for further reading is Ref. [103].
The stochastic-gradient descent algorithm (with or without minibatches) samples the input-data distribution uniformly randomly. As mentioned in Section 6.3, it may be advantageous to sample those inputs more frequently that initially cause larger output errors. More generally, the algorithm may use other criteria to choose certain input data more often, with the goal to speed up learning. It may even suggest how to augment a given training set most efficiently, by asking to specifically label certain types of input data (active learning) [104].
Another question concerns the structure of the energy landscape for multilayer perceptrons. It seems that local minima are perhaps less important for deep networks than for Hopfield networksindexHopfield network, because the energy functions of deep networks tend to have more saddle points than minima [105], just like Gaussian random functions [106]. A recent study explores the relation between the multilayer layout of the perceptron network and the properties of the energy landscape [107].
Finally, training tends to work best when all input patterns appear with roughly the same frequency in the training set. Unlike humans, neural networks may struggle with rare input patterns. Special techniques, however, allow networks to recognise rare patterns, by comparison with features of the input distribution that are well represented (few-shot learning [108]). Standard algorithms for few-shot learning use elements of unsupervised learning (Chapter 10.1).
Chapter 8 Convolutional networks
Convolutional networks have been around since the 1980’s. They became widely used after Krizhevsky et al. [109] won the ImageNet challenge (Section 8.5) with a convolutional net. One reason for the recent success of convolutional networks is that they have fewer connections than fully connected networks with the same number of neurons. This has two advantages. Firstly, such networks are obviously cheaper to train. Secondly, as pointed out above, reducing the number of connections regularises the network, it reduces the risk of overfitting.
Convolutional neural networks are designed for object recognition and image classification. They take images as inputs (Figure 7.1), not just a list of attributes (Figure 5.1). Convolutional networks have important properties in common with networks of neurons in the visual cortex of the human brain [4]. First, there is a spatial array of input terminals. For image analysis this is the two-dimensional array of bits shown in Figure 8.1(a). Second, neurons are designed to detect local features of the image (such as edges or corners for instance). The maps learned by such neurons, from inputs to output, are referred to as feature maps. Since these features occur in different parts of the image, one uses the same kernel (or filter) for different parts of the image, always with the same weights and thresholds. Since these kernels are local, and since they act in a translational-invariant way, the number of neurons from the two-dimensional input array is greatly reduced, compared with fully connected networks. Feature maps are obtained by convolution of the kernel with the input image. Therefore, layers consisting of a number of feature maps corresponding to different kernels are also referred to as convolution layers, Figure 8.1(b).
Convolutional networks can have many convolution layers. The idea is that the additional layers can learn more abstract features (Section 8.7). Apart from feature maps, convolutional networks contain other types of layers. Pooling layers perform local averages over the output of the convolution layers, to speed up learning by reducing the number of variables. Convolutional networks may also contain fully connected layers.
8.1 Convolution layers
Figure 8.1(a) illustrates how a feature map is obtained by convolution of the input image with a kernel which reads a part of the input image [5]. In analogy with the terminology used in neuroscience, this array is called the local receptive field of the kernel. The outputs of the kernel from different parts of the input image make up the feature map, here an array of hidden neurons: neuron connects to the area in the upper left-hand corner of the input image. Neuron connects to a shifted area, as illustrated in Figure 8.1(a), and so forth. Since the input has pixels, the dimension of the feature map is in this example. The important point is that the neurons and , and all other neurons in this convolution layer, share their weights and the threshold. In the example shown in Figure 8.1(a) there are thus only nine independent weights, and one threshold. Since the different hidden neurons share
weights and thresholds, their computation rule is a discrete convolution [4]:
[TABLE]
In Figure 8.1(a) the local receptive field is shifted by one pixel at a time. Sometimes it is useful to use a larger stride , to shift the receptive field by pixels horizontally and by pixels vertically. Also, the local receptive regions need not have size . If we assume that their size is , and that , the rule (8.1) takes the form
[TABLE]
Figure 8.1(a) depicts a two-dimensional input array. For colour images there are three colour channels, in this case the input array is three-dimensional, and the input bits are labeled by three indices: two for position and the last one for colour, . Usually one connects several feature maps with different kernels to the input layer, as shown in Figure 8.1(b). The different kernels detect different features of the input image, one detects edges for example, and another one detects corners, and so forth. To account for these extra dimensions, one groups weights (and thresholds) into higher-dimensional arrays (tensors). The convolution takes the form:
[TABLE]
(see Figure 8.2). All neurons in a given convolution layer have the same threshold. The software package TensorFlow [110] is designed to efficiently perform tensor operations as in Equation (8.3).
If one couples several convolution layers together, the number of neurons in these layers decreases as one moves to the right. To avoid this, one can pad the image (and the convolution layers) by adding rows and columns of bits set to zero [4]. In Figure 8.1(a), for example, one obtains a convolution layer of the same dimension as the original image by adding one column each on the left-hand and right-hand sides of the image, as well as two rows, one at the bottom and one at the top. In general, the numbers of rows and columns need not be equal, so the amount of padding is specified by four numbers, .
Convolution layers are trained with backpropagation. Consider the simplest case, Equation (8.1). As usual, we use the chain rule to evaluate the gradients:
[TABLE]
with local field . The derivative of is evaluated by applying rule (5.25):
[TABLE]
In this way one can train networks with several stacked convolution layers too. It is important to keep track of the summation boundaries. To that end it helps to pad out the image and the convolution layers, so that the upper bounds remain the same in different layers.
Details aside, the fundamental principle of feature maps is that the map is applied in the same form to different parts of the image (translational invariance). In this way, each weight in a given feature map is trained on different parts of the image. This effectively increases the training set for the feature map and combats overfitting.
8.2 Pooling layers
Pooling layers process the output of convolution layers. A neuron in a pooling layer takes the outputs of several neighbouring feature maps and compresses their outputs into a single number [5]. There are no weights or thresholds associated with the pooling layers. Max-pooling units,
for example, take the maximum over several nearby feature-map outputs. Instead, one may compute the root-mean square of the map values (-pooling). Just as for convolution layers, we need to specify stride and padding for pooling layers. Other ways of pooling are discussed in Ref. [4].
Usually several feature maps are connected to the input. Pooling is performed independently for each feature map [5]. The network layout looks like the one shown schematically in Figure 8.3. In this Figure, the pooling layers connect to a number of fully connected hidden layers that feed into the output neurons. There are as many output neurons as there are classes to be recognised. This layout is similar to the layout used by Krizhesvky et al. [109] in the ImageNet challenge (Section 8.5).
8.3 Learning to read handwritten digits
Figure 8.4 shows patterns from the MNIST data set of handwritten digits [111]. The data set derives from a data set compiled by the National Institute of Standards and Technology (NIST), of digits handwritten by high-school students and employees of the United States Census Bureau. The data contains 60 000 images of digits, each with pixels, and a test set of 10 000 digits. The images are grayscale with 8-bit resolution, so each pixel contains a value ranging from 0 to 255. The images in the database were preprocessed. The procedure is described on the MNIST home page. Each original binary image from the National Institute of Standards and Technology was represented as a 2020 gray-scale image, preserving the aspect ratio of the digit. The resulting image was placed in a 2828 image so that the centre-of-mass of the image coincided with its geometrical centre. These preprocessing steps improve the performance of the learning algorithm.
The goal of this Section is to show how the principles described up to now allow neural networks to learn the MNIST data with low classification error, following Ref. [5]. As described in Chapter 6, one divides the data set into a training set and a validation set, here with 50 000 digits and 10 000 digits, respectively [5]. The validation set is used for cross validation. The test data set allows to measure the classification error after training. For this purpose one must use a data set that was not involved in the training. As described in Section 6.3, the inputs are preprocessed further by subtracting the mean image averaged over the whole training set from each input image [Equation (6.21)].
To find good parameter values and network layouts is one of the main difficulties when training a neural network, and it usually requires a fair deal of experimenting. There are recipes for finding certain parameters [112], but the general approach is still trial and error [5]. Consider first a network with one hidden layer with ReLU activation functions (Section 7.3), and a softmax output layer (Section 7.5) with ten outputs and energy function (7.38). Output is interpreted as the probability that the pattern fed to the network falls into category . The networks are trained with stochastic gradient descent with momentum, Equation (6.31). The learning rate is set to , and the momentum constant to . The mini-batch size [Equation (6.18)] equals 8192. Cross validation and early stopping is implemented as follows: during training, the algorithm keeps track of the smallest validation error observed so far. Training stops when the validation error becomes larger than the minimum for a specified number of times, equal to in this case.
Figure 8.5 shows how the training and the validation energies decrease during training, for networks with and hidden neurons [5]. One epoch corresponds to applying patterns or iterations (Section 6.1). The energies are a little lower for the network with hidden neurons. But one observes overfitting in both cases: after many training steps the validation energy is much higher than the training energy. Early stopping caused the training of the larger network to abort after epochs, this corresponds to iterations. The resulting classification accuracy is about 97.2% for the network with 100 hidden neurons.
It is difficult to increase the classification accuracy by adding more hidden layers, most likely because the network overfits the data (Section 6.4). This problem becomes more acute as one adds more hidden neurons. The tendency of the network to overfit is reduced by regularisation (Section 7.6). For the network with one hidden layer with 100 ReLU neurons, -regularisation improves the classification accuracy to almost 98%.
Convolutional networks can be optimised to yield higher classification accuracies than those quoted above. A convolutional network with one convolution layer with 20 feature maps, a max-pooling layer, and a fully connected hidden layer with 100 ReLU neurons, similar to the network shown schematically in Figure 8.6, gives classification accuracy only slightly above 98% after training for 60 epochs. Adding a second convolution layer and batch normalisation (Section 7.6.5) gives a classification accuracy is 98.99% after 30 epochs (this layout is similar to a layout described in MathWorks [113]). The accuracy can be improved further by tuning parameters and network layout, and by using ensembles of convolutional neural networks [111]. The best classification accuracy found in this way is 99.77% [114]. Several of the MNIST digits are difficult to classify for humans too (Figure 8.7), so we conclude that convolutional networks really work very well.
The above examples show also that it takes much experimenting to find the right parameters and network layout, as well as long training times to reach the best classification accuracies. It could be argued that one reaches a stage of diminishing returns as the classification error falls below a fraction of a percent.
8.4 Coping with deformations of the input distribution
How well does a MNIST-trained convolutional network classify your own hand-written digits? Figure 8.8(a) shows examples of digits drawn by colleagues at the University of Gothenburg, preprocessed in the same way as the MNIST data.
Using a MNIST-trained convolutional network on these digits yields a classification accuracy of about 90%, substantially lower than the classification errors quoted in the previous Section.
A possible cause is that the digits in Figure 8.8(a) have a more slender stroke than those in Figure 8.4. It was suggested in Ref. [115] that differences in line thickness can confuse algorithms designed to read hand-written text [116]. There are different methods for normalising the line thickness of hand-written text. Applying the method proposed in Ref. [116] to our digits results in Figure 8.8(b). The algorithm has a free parameter, , that specifies the line thickness. In Figure 8.8(b) it was taken to be , close to the average line thickness of the MNIST digits, which is approximately . If we run a MNIST-trained convolutional network on a data set of 60 digits with normalised line thickness, it fails on only two digits. This corresponds to a classification accuracy of roughly , not so bad – yet not as good as the best results in Section 8.3. But note that this estimate of the classification accuracy is not very precise, because the test set had only 60 digits. To obtain a better estimate, more test digits are needed.
A question is of course whether there are other significant differences between our non-MNIST hand-written digits and those in the MNIST data. At any rate, the results of this Section raise a point of fundamental importance. We have seen that convolutional networks can be trained to represent a distribution of input patterns with very high accuracy. But the network may not work as well on a data set with a slightly different input distribution, perhaps because the patterns were preprocessed differently, or because they were slightly deformed in other ways.
8.5 Deep learning for object recognition
Deep learning has become so popular in the last few years because deep convolutional networks are good at recognising objects in images. Figure 8.9 shows a frame from a movie taken by a data-collection vehicle. A convolutional network was trained to recognise objects, and to localise them in the image by means of bounding boxes around the objects.
Convolutional networks excel at this task, as demonstrated by the ImageNet large-scale visual recognition challenge (ILSVRC) [117], a competition for object recognition and localisation in images, based upon the ImageNet database [100]. The challenge is based on a subset of ImageNet. The training set contains more than images manually classified into one of 1000 classes. There are approximately 1000 images for each class. The validation set contains 50 000 images.
The ILSVRC challenge consists of several tasks. One task is image classification, to list the object classes found in the image. A common measure for accuracy is the so-called top-5 error for this classification task. The algorithm lists the five object classes it with the highest softmax outputs. The result is defined to be correct if the annotated class is among these five. The error equals the fraction of incorrectly classified images. Why does one not simply judge whether the most probable class is the correct one? The reason is that the images in the ImageNet database are annotated by a single-class identifier. Often this is not unique. The image in Figure 8.9, for example, shows not only a car but also trees, yet the image is annotated with the class label car. The resulting classification ambiguity is reduced by considering the top five softmax outputs, and checking whether the annotated class is among them.
The tasks in the ILSVRC challenge are significantly more difficult than the digit recognition described in Section 8.3. One reason is that the ImageNet classes are organised into a deep hierarchy of subclasses. This results in highly specific sub classes that can be difficult to tell apart. The algorithm must be very sensitive to small differences between similar sub classes. We say that the algorithm must have high inter-class variability [120]. Different images in the same sub class, on the other hand, may look quite different. The algorithm should nevertheless recognise them as similar, belonging to the same class. In other words, the algorithm should have small intra-class variability [120].
Since 2012, algorithms based on deep convolutional networks won the ILSVRC challenge. Figure 8.10 shows that the error has significantly decreased until 2017, the last year of the challenge in the form described above. We saw in previous Sections that deep networks are difficult to train. So how can these algorithms work so well? It is generally argued that the recent success of deep convolutional networks is mainly due to three factors.
First, there are now much larger and better annotated training sets available. ImageNet is an example. Excellent training data is now recognised as one of the most important factors. Companies developing software for self-driving cars and systems that help to avoid accidents understand that good training sets are indispensable. At the same time, it is a challenge to create high-quality training data, because one must manually collect and annotate the data (Figure 8.11). This is costly, also because it is important to have as large data sets as possible, in order to reduce overfitting. In addition one must aim for a large variability in the collected data. Second, the hardware is much better today. Deep networks are nowadays implemented on single or multiple GPUs. There are even dedicated chips for this purpose. Third, improved regularisation techniques (Section 7.6) help to fight overfitting, and skipping connections (Section 7.4) render the networks less susceptible to the vanishing-gradient problem (Section 7.2).
The winning algorithm for 2012 was based on a network with five convolution layers and three fully connected layers, using drop out, ReLU activation functions, and data-set augmentation [109]. The algorithm was implemented on GPU processors. The 2013 ILSVRC challenge was also won by a convolutional network [121], with 22 layers. Nevertheless, the network has substantially fewer free parameters (weights and thresholds) than the 2012 network: 4106 instead of 60 106. In 2015, the winning algorithm [92] had 152 layers. One significant new element in the layout was the idea to allow connections that skip layers (Section 7.4). The best algorithms in 2016 [122] and 2017 [119] used ensembles of convolutional networks, where the classification is based on the ensemble average of the outputs.
8.6 Summary
Convolutional networks can be trained to recognise objects in images with high accuracy. An advantage of convolutional networks is that they have fewer weights than fully connected networks with the same number of neurons, and that the weights of a given feature map are trained on different parts of the input images, effectively increasing the size of the training set. This helps against overfitting. Another view is that the hidden neurons are forced to agree on a particular choice of weights, they must compromise. This yields a more robust training result.
It is sometimes stated that convolutional networks are now better than humans, in that they recognise objects with lower classification errors than humans [123]. This and similar statements refer to an experiment showing that the human classification error in recognising objects in the ImageNet database is about 5.1% [124], worse than the most recent convolutional neural-network algorithms (Figure 8.10).
This notion is not unproblematic, for several reasons. To begin with, the article [123] refers to the 2015 ILSVRC competition, where the top scores were quite similar, and it has been debated whether interpreting the rules of the competition in different ways allowed competitors to gain an advantage. Second, and more importantly, it is clear that these algorithms learn in quite a different way from humans. The algorithms can detect local features, but since these convolutional networks rely on translational invariance, they do not easily understand global features, and can mistake a leopard-patterned sofa for a leopard [125]. It may help to include more leopard-patterned sofas in the training set, but the essential difficulty remains: translational invariance imposes constraints on what convolutional networks can learn [125]. More fundamentally one may argue that humans learn differently, by abstraction instead of going through very large training sets.
We have also seen that convolutional networks are sensitive to small changes in the input data. Convolutional networks excel at learning the properties of a given input distribution, but they may have difficulties in recognising patterns sampled from a slightly different distribution, even if the two distributions appear to be very similar to the human eye. Note also that this problem cannot be solved by cross validation, because training and validation sets are drawn from the same input distribution, but here we are concerned with what happens when the network is applied to a input distribution different from the one it was trained on.
Here is another example illustrating this point: the authors of Ref. [126] trained a convolutional network on perturbed grayscale images from the ImageNet data base, adding a little bit of noise independently to each pixel (white noise) before training. This network failed to recognise images that were weakly perturbed in a different way, by setting a small number of pixels to white or black. But when we look at the images we have no difficulties seeing through the noise.
Refs. [127, 128] illustrate intriguing failures of convolutional networks [5]. Szegedy et al. [127] demonstrate that the way convolutional networks partition input space can lead to unexpected results. The authors took an image that the network classifies correctly with high confidence, and perturbed it slightly. The perturbation was not random, but specifically designed to push the input pattern over a decision boundary. The difference between the original and perturbed images (adversarial images) is undetectable to the human eye, yet the network misclassifies the perturbed image with high confidence [127]. This reflects the fact that decision boundaries are always close in high-dimensional input space.
Figure 1 in Ref. [128] shows images that are completely unrecognisable to the human eye. Yet a convolutional network classifies these images with high confidence. This illustrates that there is no telling what a network may do if the input is far away from the training distribution. Unfortunately the network can sometimes be highly confident yet wrong. Nevertheless, despite these problems, deep convolutional networks have enjoyed tremendous success in image classification during the past years, and they have found widespread use in industry and science.
Finally, the fundamental mechanisms of deep learning are quite well understood, but many open questions remain. It is fair to say that the theory of deep learning has somewhat lagged behind the associated practical successes, although some progress has been made in recent years.
8.7 Further reading
The online book of Nielsen [5] is an excellent introduction to convolutional neural networks, and guides the reader through all the steps required to program a convolutional network to recognise hand-written digits. Nielsen’s chapter Deep Learning [5] is the main source for Section 8.3.
What do the hidden layers in a convolutional layer actually compute? Feature maps that are directly coupled to the inputs detect local features, such as edges or corners. Yet it is unclear precisely how hidden convolutional layers help the network to learn. To which input features do the neurons of a certain hidden layer react most strongly? Input patterns chosen to maximise the outputs of neurons in a given layer [129, 130] reveal intricate geometric structures that defy straightforward interpretation. An example is shown on the cover of this book, see also Exercise 8.7.
It has been suggested that more general models, normally used for natural-language processing, may outperform convolutional nets in image-processing tasks when there is enough data [131]. An advantage is that these models do not rely on translational invariance, unlike convolutional networks.
Chapter 9 Supervised recurrent networks
The layout of the perceptrons analysed in the previous Chapters is special. All connections are one way, and only to the layer to the right, so that the update rule for the -th neuron in layer becomes, for example,
[TABLE]
The backpropagation algorithm relies on this feed-forward layout. It means that the derivatives vanish. This ensures that the outputs are nested functions of the inputs, which in turn implies the simple iterative structure of the backpropagation algorithm (Chapter 6).
In some cases it is necessary or convenient to use networks that do not have this simple layout. The Hopfield networks discussed in part I are examples where all connections are symmetric. More general networks may have a feed-forward layout with feedbacks, as shown in Figure 9.1. Such networks are called recurrent networks. There are many different ways in which the feedbacks can act: from the output layer to hidden neurons for example, or there could be connections between the neurons in a given layer. Neurons and in Figure 9.1 are output neurons, they are associated with targets just as in Chapters 5 to 7. The layout of recurrent networks is very general, but because of the feedbacks we must consider how such networks can be trained.
Unlike multilayer perceptrons that represent an input-to-output mapping in terms of nested activation functions, recurrent networks are used as dynamical networks, where the iteration index replaces the layer index :
[TABLE]
See Figure 9.1 for the definition of the different weights, and the parameters are thresholds. Equation (9.2) is analogous to the deterministic McCulloch-Pitts dynamics of Hopfield networks and Boltzmann machines [c.f. Equation (1.5)]. As in the case of Hopfield networks (Exercise 2.10), one may also consider a continuous network dynamics:
[TABLE]
with time constant . We shall see in a moment why it is convenient to assume that the dynamics is continuous in , as in Equation (9.3).
Recurrent networks can learn in different ways. One possibility is to use a training set of pairs [\mbox{\boldmathx}^{(\mu)},\mbox{\boldmathy}^{(\mu)}] with . To avoid confusion with the iteration index , the targets are denoted by in this Chapter. One feeds a pattern from this set and runs the dynamics (9.2) or (9.3) for the given \mbox{\boldmathx}^{(\mu)} until it reaches a steady state \mbox{\boldmathV}^{\ast} (if this does not happen, the training fails). Then one adjusts the weights by one gradient-descent step using the energy function
[TABLE]
The asterisk in this Equation indicates that all variables are evaluated in the steady state, at \mbox{\boldmathV}=\mbox{\boldmathV}^{\ast}. Iterating these steps, one feeds another pattern \mbox{\boldmathx}^{(\mu)}, finds the steady state \mbox{\boldmathV}^{\ast}, adjusts the weights and so forth. Instead of defining the energy function in terms of the mean-squared output errors, one could also use the negative log-likelihood function (7.44). One continues to iterate these steps until the steady-state outputs yield the correct targets for all input patterns. This is reminiscent of the algorithms discussed in Chapters 5 to 7.
Another possibility is that inputs and targets change as functions of time while the network dynamics runs. This allows to solve temporal classification tasks. The network is trained on a set of input sequences \mbox{\boldmathx}(t) and corresponding target sequences \mbox{\boldmathy}(t). In this way, recurrent networks can translate written text or recognise speech. The network can be trained by unfolding its dynamics in time as explained in Section 9.2, although this algorithm suffers from the vanishing-gradient problem discussed in Chapter 7.
9.1 Recurrent backpropagation
This Section summarises how to generalise Algorithm 4 to recurrent networks with feedback connections. Recall the recurrent network shown in Figure 9.1. The neurons have smooth activation functions, and they are connected by weights . Several neurons may be linked to inputs , with weights . Other neurons are output units with associated target values .
One takes the dynamics to be continuous in time, Equation (9.3), and assumes that \mbox{\boldmathV}(t) runs into a steady state
[TABLE]
Equation (9.3) implies
[TABLE]
and it is assumed that \mbox{\boldmathV}^{\ast} is a linearly stable steady state of the dynamics (9.3), so that small perturbations \delta\mbox{\boldmathV} away from \mbox{\boldmathV}^{\ast} decay with time.
The synchronous discrete dynamics (9.2) can exhibit undesirable stable periodic solutions [133], as mentioned in Section 1.3. This is a reason for using the continuous dynamics (9.3), yet convergence to the steady state is not guaranteed in this case either.
Equation (9.6) is a non-linear self-consistent condition for the components of \mbox{\boldmathV}^{\ast}, in general difficult to solve. However, if the steady state \mbox{\boldmathV}^{\ast} is stable, we can use the dynamics (9.3) to automatically pick out the steady-state solution \mbox{\boldmathV}^{\ast}. This solution depends on the pattern \mbox{\boldmathx}^{(\mu)}. Note that the superscript is left out in Equations (9.5) and (9.6) and also in the remainder of this Section.
The goal is to find weights so that the outputs give the correct target values in the steady state. To this end one uses gradient descent on the energy function (9.4). Following Ref. [1], consider first how to adjust the weights :
[TABLE]
One calculates the gradients of \mbox{\boldmathV}^{\ast} by differentiating Equation (9.6):
[TABLE]
Here is the local field in the steady state. Equation (9.8) is a self-consistent equation for the gradient, as opposed to the explicit expressions we found in Chapter 6. The reason for the difference is that the recurrent network has feedbacks.
Since Equation (9.8) is linear in the gradients, it can be solved by matrix inversion, at least formally. In terms of the matrix with elements
[TABLE]
Equation (9.8) can be written as
[TABLE]
If is invertible, one applies to both sides. Using the fact that one finds:
[TABLE]
Inserting this result into (9.7) one obtains:
[TABLE]
This learning rule can be written in the form of the backpropagation rule (6.10) by introducing the error
[TABLE]
Then the learning rule (9.12) takes the form
[TABLE]
If there are no recurrent connections, then . In this case Equation (9.13) reduces to the standard expression (6.6b), Exercise 9.1.
The learning rule for the weights is derived in an analogous fashion. The result is:
[TABLE]
The learning rules (9.14) and (9.15) are well-defined only if the matrix is invertible. Otherwise the solution (9.11) does not exist. Also, matrix inversion is an expensive operation. As described in Chapter 5, one can try to avoid the problem by finding the inverse iteratively. The trick [1] is to write down a dynamical equation for that has a steady state at the solution of Equation (9.13):
[TABLE]
It is left as an exercise (Exercise 9.2) to verify that the dynamics (9.16) has a steady state satisfying Equation (9.13). Equation (9.16) is written in a form to stress that (9.16) and (9.3) exhibit the same duality as Algorithm 4, between forward propagation of states of neurons and backpropagation of errors. The sum in Equation (9.16) has the same form as the recursion for the errors in Algorithm 4, except that there are no layer indices here.
Equation (9.16) admits the steady state (9.13). But does converge to ? For convergence it is necessary that the steady state is linearly stable. Whether or not this is the case is determined by linear stability analysis [85]. One asks: does a small deviation from the steady state increase or decrease under Equation (9.16)? In other words, if one writes
[TABLE]
do \delta\!\mbox{\boldmathV}(t) and \delta\!\mbox{\boldmath\Delta}(t) grow in magnitude? To answer this question, one inserts this ansatz into (9.3) and (9.16), and linearises:
[TABLE]
Equation (9.18a) shows: whether or not the norm of \delta\mbox{\boldmathV}(t) grows is determined by the eigenvalues of the matrix . We say that \mbox{\boldmathV}^{\ast} is a linearly stable steady state of Equation (9.3) if all eigenvalues of have negative real parts. In this case |\delta\mbox{\boldmathV}(t)|\to 0. If at least one eigenvalue has a positive real part then |\delta\!\mbox{\boldmathV}| grows. In this case we say that \mbox{\boldmathV}^{\ast} is linearly unstable. Since the matrix with elements has the same eigenvalues as , \mbox{\boldmath\Delta}^{\ast} is a stable steady state of (9.16) if \mbox{\boldmathV}^{\ast} is a stable steady state of (9.3). If the steady states are unstable, the algorithm does not converge.
In summary, recurrent backpropagation is analogous to backpropagation (Algorithm 4) for layered feed-forward networks, save for two differences. First, the non-linear network dynamics is no longer a simple input-to-output mapping with nested activation functions, but a non-linear dynamics that may (or may not) converge to a steady state. Second, the feedbacks give rise to linear self-consistent equations for the
steady-state gradients , which can be viewed as steady-state conditions for a dual dynamics of the errors.
The main conclusion of this Section is that convergence of the training is not guaranteed if the network has feedback connections (for a layered feed-forward network without feedbacks, recurrent backpropagation simplifies to stochastic gradient descent, Algorithm 4, see Exercise 9.1). This explains why stochastic gradient descent is used mostly for multi-layer networks with feed-forward layouts. The algorithm tends to fail for networks with feedbacks. However, it is possible to get rid of the feedbacks in recurrent networks by unfolding the dynamics in time.
This is described in the next Section.
9.2 Backpropagation through time
Recurrent networks can be used to learn sequential inputs, as in speech recognition and machine translation. The training set consists of time sequences [\mbox{\boldmathx}(t),\mbox{\boldmathy}(t)] of inputs and targets. The network is trained on the sequences and learns to predict the targets. In this context the layout differs from the one described in the previous Section. There are two main differences. Firstly, the inputs and targets depend on , and one uses a discrete-time update rule. Secondly, separate output neurons are added to the layout. The update rule takes the form
[TABLE]
The activation function of the output neurons can be different from that of the hidden neurons . One possibility is to use the softmax function for the outputs [134, 135]. For the hidden neurons one often uses tanh activations.
To train recurrent networks with time-dependent inputs and with the dynamics (9.19), one uses backpropagation through time. The idea is to unfold the network in time to get rid of the feedbacks. The price paid is that one obtains large networks in this way, with as many copies of the original neurons as there are time steps.
The procedure is illustrated in Figure 9.2 for a recurrent network with one hidden neuron, one input terminal, and one output neuron. The unfolded network has inputs and outputs. It can be trained in the usual way with stochastic gradient descent. The errors are calculated using backpropagation as in Algorithm 4, but here the error is propagated back in time, not from layer to layer. The energy function is the squared error summed over all time steps
[TABLE]
One could use the negative log-likelihood function (7.38), but here we use the squared output-error function (9.20). There is only one hidden neuron in our example, and the inputs and outputs are also one-dimensional. Here and in the following we write the time argument as a subscript, instead of and so forth, because there is no risk of confusing the time index with other subscripts.
Consider first how to adjust the weight . Gradient descent (5.24) yields
[TABLE]
Here
[TABLE]
is an output error, and is the local field of the output neuron at time [Equation (9.19)]. Equation (9.21a) is similar to the learning rule for recurrent backpropagation, Equations (9.7) and (9.8), but the derivative is evaluated differently. Equation (9.19a) yields the recursion
[TABLE]
for . Since , Equation (9.23) implies:
[TABLE]
Equation (9.21a) says that we must sum over . Regrouping the terms in this sum yields:
[TABLE]
To write the learning rule in the usual form, we define errors recursively:
[TABLE]
Then the learning rule (9.21a) takes the form
[TABLE]
just like Equation (6.9), or like the error recursion in Algorithm 4. The factor in the recursion (9.24) gives rise to a product of many such factors in when is large, exactly as described in Section 7.2 for multilayer perceptrons. This means that the training of recurrent networks suffers from unstable gradients, as backpropagation of multilayer perceptrons does: if the factors are smaller than unity, then the errors become very small when becomes small (vanishing-gradient problem). This means that the early states of the hidden neuron no longer contribute to the learning, causing the network to forget what it has learned about early inputs. When , on the other hand, exploding gradients make learning impossible. In summary, unstable gradients in recurrent neural networks occur much in the same way as in multilayer perceptrons. The resulting difficulties for training recurrent neural networks are discussed in more detail in Section 9.3, see also Ref. [136].
A slight variation of the above algorithm (truncated backpropagation through time) suffers less from the exploding-gradient problem. The idea is that the exploding gradients are tamed by truncating the memory. This is achieved by limiting the error propagation backwards in time, errors are computed back to and not further, where is the truncation time [2]. Naturally this implies that long-time correlations cannot be learnt.
The learning rules for the weights are obtained in a similar fashion. Equation (9.19a) yields the recursion
[TABLE]
This looks just like Equation (9.23), except that is replaced by . As a consequence we have
[TABLE]
The learning rule for is simpler to derive. From Equation (9.19b) we find by differentiation w.r.t. :
[TABLE]
How are the thresholds adjusted? Going through the above derivation we see that we must replace in Equation (9.25) by . It works in the same way for the output threshold.
In order to keep the formulae simple, we derived the algorithm for a single hidden neuron, a single output neuron, and one-component inputs, so that we could leave out the indices referring to different hidden neurons, and different input and output components. If we consider several hidden and output neurons and multi-dimensional inputs, the structure of the Equations remains exactly the same, except for a number of extra sums over those indices:
[TABLE]
The second term in the recursion for is analogous to the error recursion in of Algorithm 4. The time index here plays the role of the layer index in Algorithm 4. A difference is that the weights in Equation (9.29) are the same for all time steps. The algorithm is summarised in Algorithm 7.
In conclusion we see that backpropagation through time for recurrent networks is similar to backpropagation for multilayer perceptrons. After the recurrent network is unfolded to get rid of the feedback connections, it can be trained by backpropagation. The time index takes the role of the layer index . Backpropagation through time is the standard approach for training recurrent networks, despite the fact that it suffers from the vanishing-gradient problem. The next Section describes how improvements to the layout make it possible to more efficiently train recurrent networks.
9.3 Vanishing gradients
Hochreiter and Schmidhuber [137] suggested to replace the hidden neurons of the recurrent network with computation units that are specially designed to reduce the vanishing-gradient problem. The method is referred to as long-short-term memory (LSTM). The basic ingredient is the same as in residual networks (Section 7.4): short cuts reduce the vanishing-gradient problem. For our purposes we can think of LSTMs as units that replace the hidden neurons. For a detailed description of LSTMs see Ref. [138].
Gated recurrent units [139] serve the same purpose as LSTMs, and they function in a similar way. It has been argued that LSTMs outperform gated recurrent units for certain tasks, but since they are simpler than LSTMs, the remainder of this Section focuses on gated recurrent units. As illustrated in Figure 9.3, these units replace the hidden neurons of a recurrent neural network [with update rule (9.19a)] by a new rule:
[TABLE]
The first two Equations are referred to as gates because they regulate how the values of the hidden state variables are passed through the unit. Here is the sigmoid function (6.19a). If for all , and for all , Equation (9.30) coincides with the standard update rule (9.19a), save for the thresholds which were left out in Equation (9.30). As explained above, the resulting recurrent network suffers from the vanishing-gradient problem. This means that states in the past history, \mbox{\boldmathV}(0),\mbox{\boldmathV}(1),\ldots, have little effect upon the present state \mbox{\boldmathV}(t) for . The recurrent network forgets early inputs, so that it cannot learn from them.
If by contrast for all , then the input is passed right through the unit. Since in this case, the gradients do not decrease as the dynamics explores the history. However, since \mbox{\boldmathV}(t-1)=\mbox{\boldmathV}(t), the recurrent network reproduces previous states. This is analogous to skipping layers in a residual network, although comparison with Equation (7.35) reveals some differences in detail.
For the recurrent network to learn in a meaningful way from past inputs, the weights in Equation (9.30) (and the thresholds) are adjusted so that the gated recurrent unit operates between these two extreme limits. This is achieved by including the weights and thresholds of the gated recurrent unit in the gradient-descent minimisation of the energy function (9.20). The learning rules for the weights (and thresholds) are calculated in the same as before. Using Equation (9.19b) one has:
[TABLE]
where stands for . The derivatives are evaluated using the chain rule and Equations (9.30).
It is instructive to inspect the values of and when the recurrent network operates after training. Suppose that a unit assumes small values of and . This means that the update of the state variable is determined entirely by the instantaneous inputs . Since the unit does not refer to the past history of the hidden-state variables, it truncates the dynamical memory. In the opposite limit, when , the unit can contribute to building up long-term dynamical memory. These arguments suggest that a unit with just one gate may achieve the same goal [140]:
[TABLE]
This unit is easier to train because it has fewer parameters than the standard gated recurrent unit (9.30). Yet, the additional parameters in Equation (9.30) may help to represent and exploit correlations on different time scales. LSTMs have even more parameters. How this tradeoff between ease of training and accurate representation of time correlations works out may well depend on the problem at hand. In the following Section we describe recurrent networks with LSTM units, following Refs. [134, 135].
9.4 Recurrent networks for machine translation
Recurrent networks are used for machine translation [134, 135]. How does this work? The networks are trained using backpropagation through time. The vanishing-gradient problem is dealt with by using LSTMs (Section 9.3).
How are the network inputs and outputs coded? For machine translation one respresents all words in a given dictionary in terms of a code. The conceptually simplest code is one where 100…represents the first word in the dictionary, 010…the second word, and so forth. The drawback of this scheme is that it does not account for the fact that two given words might be more or less closely related to each other. Other encoding schemes are described in Ref. [135].
Each input to the recurrent network is a vector with as many components as there are words in the dictionary. A sentence corresponds to a sequence \mbox{\boldmathx}_{1},\mbox{\boldmathx}_{2},\ldots,\mbox{\boldmathx}_{T}. Each sentence ends with an end-of-sentence tag, <EOS>. This tag tells the network when the input sentence ends. This is necessary because the number of words per sentence is not fixed. Now suppose that a possible translation reads \mbox{\boldmathx}_{1}^{\prime},\mbox{\boldmathx}_{2}^{\prime},\ldots,\mbox{\boldmathx}_{T^{\prime}}^{\prime}. The task of the network is to determine the probability p(\mbox{\boldmathx}_{1}^{\prime},\ldots,\mbox{\boldmathx}_{T^{\prime}}^{\prime}|\mbox{\boldmathx}_{1},\ldots\mbox{\boldmathx}_{T}) that the translation is correct. The idea is to estimate this probability recursively as
[TABLE]
Sutskever et al. [134] describe how to achieve this with a recurrent network with two hidden LSTMs. The network uses softmax outputs \mbox{\boldmathO}_{\!t}, where -th component of \mbox{\boldmathO}_{\!t} is interpreted as the probability that the -th component of \mbox{\boldmathx}_{t}^{\prime} is the correct word at position in the translated sentence. As shown in Figure 9.4, the first LSTM processes the input sentence \mbox{\boldmathx}_{1},\ldots,\mbox{\boldmathx}_{T}, encoding its contents in the hidden states. When the <EOS> tag appears, the second LSTM takes over, using the information encoded by the first LSTM as an input. The second LSTM recursively outputs the translated sentence word by word, using Equation (9.33).
There is a large number of recent papers on machine translation with recurrent neural networks. Most studies are based on the training algorithm described in Section 9.2, backpropagation through time. The different approaches mainly differ in their network layouts. Google’s machine translation system uses a deep network with several layers of hidden units arranged in a bidrectional layout [141]. In such bidirectional networks, different hidden units are unfolded forward as well as backwards in time, as shown schematically in Figure 9.5.
For several hidden and output neurons and multidimensional inputs, the bidirectional network has the dynamics
[TABLE]
It is natural to use bidirectional networks for machine translation because correlations go either way in a sentence, forward and backwards. In German, for example, the finite verb form is usually at the end of the sentence.
Different schemes for scoring the accuracy of a translation are described by Lipton et al. [135]. One difficulty is that there are often several different valid translations of a given sentence, and the score must compare the machine translation with all of them. More recent papers on machine translation usually use the so-called BLEU score to evaluate the translation accuracy. The acronym stands for bilingual evaluation understudy. The scheme was proposed by Papieni et al. [142]. It is argued to evaluate the accuracy of a translation not too differently from humans.
9.5 Reservoir computing
An alternative to backpropagation through time for recurrent networks is reservoir computing [143]. This method has been used with success to predict chaotic dynamics [144, 145] and rare transitions in stochastic bi-stable systems [146].
Consider input data in the form of a time series \mbox{\boldmathx}(0),\ldots\mbox{\boldmathx}(T-1) of -dimensional vectors \mbox{\boldmathx}(t), and a corresponding series of -dimensional targets \mbox{\boldmathy}(t). The goal is to train the recurrent network so that its outputs \mbox{\boldmathO}(t) approximate the targets as precisely as possible, by minimising the energy function , where is the output error, and represents an initial transient that is disregarded.
Figure 9.6 shows the layout for this task. There are input terminals. They are connected with weights
to a reservoir of hidden neurons with state variables . The reservoir is linked to linear output units with weights . The reservoir itself is a large recurrent network with weights . The update rule is similar to Equation (9.19). There are many different versions that differ in detail [147]. One possibility is [146]
[TABLE]
for with initial conditions .
The main difference to the training algorithms described in the previous Sections of this Chapter is that the input weights and the reservoir weights are randomly initialised and then kept constant. Only the output weights are trained. The idea is that the dynamics of a sufficiently large reservoir finds nonlinear, high-dimensional representations of the input data [143], not unlike sparse representations of binary classification problems embedded in a high-dimensional space that become linearly separable in this way (Section 5.4).
In addition, and this is a difference to the problem described in Section 5.4, the reservoir is a dynamical memory.
This requires that the reservoir states faithfully represent the input sequence: similar input sequences should yield similar reservoir activations, provided one iterates it long enough. However, for random weights the recurrent reservoir dynamics can be chaotic [85]. In this case, the state of the reservoir after many iterations bears no relation to the input sequence. To avoid this, one requires that the reservoir dynamics is linearly stable.
Linearising the reservoir dynamics (9.35a) gives
[TABLE]
where is a diagonal matrix with entries , and where . Whether or not \delta\mbox{\boldmathr} grows is then determined by the singular values of , as in Section 7.2. The singular values of are denoted by . At large times, when driven with a stationary input series, the maximal Lyapunov exponent must be negative to ensure that the reservoir dynamics is stable:
[TABLE]
Sometimes the stability criterion is quoted in terms of the maximal eigenvalue of . If one uses tanh activation-functions and if the local fields remain small, then the diagonal elements of remain close to unity. In this case the stability condition for the reservoir dynamics is given by the weight matrix alone. In general the singular values of are different from its eigenvalues, but what matters here is that the maximal singular value of approaches , where is the eigenvalue of with largest modulus (Exercise 9.8).
For inputs with long time correlations, the reservoir must not decay too quickly, so that it can represent the dynamical correlations in the input sequence. There is no precise mathematical theory that says how to optimise the reservoir. In practice one adjusts the maximal Lyapunov exponent by trial and error. Its optimal value depends on the properties of the input series, for instance on its time correlations.
There are many different recipes for how to set up a reservoir. Usually the reservoir is sparse, with only a small fraction of weights non-zero. The elements of the resulting weight matrix are rescaled to adjust [145]. The weight matrix is commonly taken to be a full matrix, and its elements are drawn from the same distribution as those of the reservoir. Lukosevicius [147] gives a practical overview over different schemes for setting up reservoir computers.
For time-series prediction, one trains the network on an input series \mbox{\boldmathx}(0),\ldots, \mbox{\boldmathx}(T-1) with targets \mbox{\boldmathy}(t)=\mbox{\boldmathx}(t). After training, one continues to iterate the network dynamics with inputs \mbox{\boldmathx}(T+k)=\mbox{\boldmathO}(T+k) to predict \mbox{\boldmathx}(T+k+1), for . Using Equation (9.35), we see that the reservoir dynamics takes the form r_{i}(t+1)=g\big{(}\sum_{j}\big{[}w_{ij}+\sum_{k=1}^{N}w^{({\rm in})}_{ik}w^{({\rm out})}_{kj}]r_{j}(t)\big{)} during the prediction phase. In order to represent complex spatio-temporal patterns, Pathak et al. [144] found it necessary to employ several parallel reservoirs. Lim et al. [146] used a chain of reservoirs, replacing Equation (9.35a) by a set of nested update rules.
Tanaka et al. [148] describe different physical implementations of reservoir computers, based on electronic RC-circuits, optical cavities or resonators, spin-torque oscillators, or mechanical devices.
9.6 Summary
It is sometimes said that recurrent networks learn dynamical systems, while multilayer perceptrons learn input-output maps. This emphasises a difference in how these networks are usually used, but we should bear in mind that they are trained in similar ways, by backpropagation. Neither is it given that the tasks must differ: recurrent networks are also used to learn time-independent data. It is true, however, that tools from dynamical-systems theory help to analyse the dynamics of recurrent networks [149, 136].
Recurrent neural networks can be trained by stochastic gradient descent after unfolding the network in time, to get rid of feedback connections. This algorithm suffers from the vanishing-gradient problem. To overcome this difficulty, the hidden neurons in the recurrent network are replaced by composite units that are trained to sometimes act as residual connections, passing the signal right through, and sometimes as non-linear units that can learn correlations in a meaningful way. There are different versions, long-short-term memory units and gated recurrent units. They all work in similar ways. Succesful layouts for machine translation use deep bidirectional networks with layers of LSTMs.
An alternative scheme is reservoir computing, where a large reservoir of hidden neurons is used to represent correlations in the input data, and a set of linear output units is trained to learn the original sequence from such representations. The idea is that it is easier to learn intricate features of an input sequence from a high-dimensional, sparse representation of the data.
9.7 Further reading
The training of recurrent networks is discussed in Chapter 15 of Ref. [2], see also Refs. [150, 151]. Recurrent backpropagation is described by Hertz, Krogh and Palmer [1], for a very similar network layout. How LSTMs combat the vanishing-gradient problem is explained in Ref. [138]. For a recent review of recurrent neural networks, see Ref. [135]. This webpage [152] gives a very enthusiastic overview about what recurrent networks can do. A more pessimistic view is expressed in this blog. For a review of reservoir computing, see Ref. [143].
Part III Learning without labels
Chapters 5 to 9 describe supervised learning of labeled data with neural networks. The network is trained to reproduce the correct labels (targets) for each input pattern.
The analysis of unlabeled data requires different methods. Machine learning can be applied with success to large data sets of high-dimensional unlabeled data. The machine can for instance mark patterns that are typical for the given distribution, or detect outliers. Other tasks are to detect similarity, to find clusters in the data (Figure 10.1), and to determine non-linear, low-dimensional representations of high-dimensional data. More recently, such unsupervised learning algorithms have been used to generate synthetic data, patterns that resemble those in a certain data set. One possible application is data-set augmentation for supervised learning.
Learning without labels is called unsupervised learning , because there are no targets that tell the network whether it has learnt correctly or not. There is no obvious function to fit, or dynamics to learn. Instead the network organises the input data in relevant ways. This requires redundancy in the input data. It is sometimes said that unsupervised learning corresponds to learning without a teacher, implying that the network itself discovers suitable ways of organising the input data. This is inaccurate, because unsupervised networks usually operate with a pre-determined learning rule, like Hopfield networks.
Part III of this book is organised as follows. Chapter 10.1 describes unsupervised-learning algorithms, starting with unsupervised Hebbian learning to detect familiarity and similarity of input patterns (Sections 10.1 and 10.2). Related algorithms can be used to find low-dimensional non-linear projections of high-dimensional input data (self-organising maps, Section 10.3). In Section 10.4, these algorithms are compared and contrasted with a standard unsupervised clustering algorithm, -means clustering. Section 10.5 introduces radial basis-function networks, they learn using a hybrid algorithm with supervised and unsupervised elements. Section 10.6 explains how to use layered feed-forward networks for unsupervised learning.
Chapter 11 deals with learning tasks that lie in between supervised and unsupervised learning, problems where the machine receives partial feedback on its performance in the shape of a penalty or a reward. Such tasks can be solved by reinforcement-learning algorithms that allow a neural network or more generally an agent to learn to reproduce outputs that tend to give positive rewards. Several algorithms for reinforcement learning are described, the associative reward-penalty algorithm (Section 11.1), temporal difference learning (Section 11.2), and -learning (Section 11.3). The -learning algorithm is illustrated by demonstrating how it allows two players to learn to compete in the board game tic-tac-toe.
Chapter 10 Unsupervised learning
10.1 Oja’s rule
A simple example for an unsupervised-learning algorithm uses a single McCulloch-Pitts neuron with linear activation function (Figure 10.2). The neuron computes111In this Chapter we follow a common convention [1] and denote the output of unsupervised-learning algorithms by . y=\mbox{\boldmathw}\cdot\mbox{\boldmathx} with weight vector \mbox{\boldmathw}=[w_{1},\ldots,w_{N}]^{\sf T}. Now consider a distribution P_{\rm data}(\mbox{\boldmathx}) of input patterns \mbox{\boldmathx}=[x_{1},\ldots,x_{N}]^{\sf T} with continuous-valued components . Patterns are drawn from this distribution at random and fed one after another to the net. For each pattern , the weights are adjusted as follows:
[TABLE]
This rule is also called Hebbian unsupervised learning rule [1], because it is reminiscent of Hebb’s rule (Chapter 2). As usual, is the learning rate.
What can this rule learn about the input distribution P_{\rm data}(\mbox{\boldmathx})? Since we keep adding multiples of the pattern vectors to the weights (just as described in Section 5.2), the magnitude of the output becomes the larger the more often the input pattern occurs in the distribution P_{\rm data}(\mbox{\boldmathx}). So the most familiar pattern produces the largest output. In this way the network can detect how familiar certain input patterns are.
A problem is that the components of the weight vector continue to grow as we keep adding. This means that the simple Hebbian learning rule (10.1) does not converge to a steady state. To analyse learning outcomes we want the learning to converge. This is achieved by adding a weight-decay term (Section 7.6.1) with coefficient proportional to to Equation (10.1):
[TABLE]
Making use of y=\mbox{\boldmathw}\cdot\mbox{\boldmathx}=\mbox{\boldmathw}^{\sf\small T}\mbox{\boldmathx}=\mbox{\boldmathx}^{\sf\small T}\mbox{\boldmathw}, Equation (10.2) can be rewritten in the following form:
[TABLE]
This learning rule is called Oja’s rule [155]. Equation (10.3) ensures that remains normalised. To see why, consider an analogy: a vector that obeys the differential equation
[TABLE]
For a general matrix , the norm |\mbox{\boldmathq}| may increase or decrease, depending on the singular values of . We can ensure that remains normalised by adding a term to Equation (10.4):
[TABLE]
The vector turns in the same way as , and if we set |\mbox{\boldmathw}|=1 initially, then remains normalised, \mbox{\boldmathw}=\mbox{\boldmathq}/|\mbox{\boldmathq}| (Exercise 10.1). Equation (10.5) describes the dynamics of the normalised orientation vector of a small rod in turbulence [156], where is the matrix of fluid-velocity gradients.
Returning to Equation (10.2), we note that the dynamics of (10.2) and (10.5) is the same in the limit of small learning rates . Therefore we conclude that remains normalised under (10.3) when the learning rate is small enough. Oja’s algorithm is summarised in Algorithm 8. One draws a pattern from the distribution P_{\rm data}(\mbox{\boldmathx}) of input patterns, applies it to the network, and updates the weights as prescribed in Equation (10.2). This is repeated many times. In the following we denote the average over input patterns as .
While the rule (10.1) does not have a steady state, Oja’s rule (10.3) does. For zero-mean input data, its steady state \mbox{\boldmathw}^{\ast} corresponds to the principal component of the input data, as illustrated in Figure 10.3. This can be seen by analysing the steady-state condition
[TABLE]
Here \langle\cdots\rangle_{\small\mbox{\boldmathw}^{*}} is an average over iterations of the learning rule (10.3) at fixed \mbox{\boldmathw}^{\ast}, the steady state. Equation (10.6) says that the weight increments \delta\!\mbox{\boldmathw} must average to zero in the steady state, to ensure that the weights neither grow nor decrease in the long run. Equation (10.6) is a condition upon \mbox{\boldmathw}^{\ast}. Using the learning rule (10.3), it can be written as:
[TABLE]
Equation (10.7) shows that \mbox{\boldmathw}^{*} must be an eigenvector of the matrix222For zero-mean input data, equals the data-covariance matrix, Equation (6.24). , normalised to unity, |\mbox{\boldmathw}^{\ast}|=1. But which one?
We denote the eigenvectors and eigenvalues of by \mbox{\boldmathu}_{\alpha} and , and investigate the stability of \mbox{\boldmathw}^{\ast}=\mbox{\boldmathu}_{\alpha} for different values of by linear stability analysis, just as in Section 9.1. To this end, consider a small perturbation \mbox{\boldmath\varepsilon}_{\!t} away from \mbox{\boldmathw}^{\ast}=\mbox{\boldmathu}_{\alpha}:
[TABLE]
A difference to the analysis in Section 9.1 is that the dynamics is discrete in time. The perturbation at the next time step, \mbox{\boldmath\varepsilon}_{t+1}, is defined by \mbox{\boldmathw}_{t+1}=\mbox{\boldmathu}_{\alpha}+\mbox{\boldmath\varepsilon}_{t+1}. A second difference is that the sequence of weight increments depends on the randomly chosen input patterns. In order to determine the linear stability one should iterate and then linearise the dynamics (10.3), to see whether \mbox{\boldmath\varepsilon}_{\!t} grows or not. However, in the limit of small learning rate it is sufficient to average over before iterating (Exercise 10.4). To linear order in \mbox{\boldmath\varepsilon}_{\!t} one finds:
[TABLE]
where the last equality sign defines the matrix . The steady state \mbox{\boldmathw}^{\ast}=\mbox{\boldmathu}_{\alpha} is linearly stable if all eigenvalues of have real parts with magnitudes smaller than unity.333For time-continuous dynamics (Section 9.1), linear stability is ensured when all eigenvalues have negative real parts, for discrete dynamics their magnitudes must all be smaller than unity [85]. To determine the eigenvalues of , we use the fact that has the same eigenvectors as . Since is symmetric, these eigenvectors form an orthonormal basis, \mbox{\boldmathu}_{\alpha}\cdot\mbox{\boldmathu}_{\beta}=\delta_{\alpha\beta}. As a consequence, the eigenvalues of are simply given by
[TABLE]
Since is a positive-semidefinite matrix (its eigenvalues cannot be negative), Equation (10.10) shows that there are eigenvalues with unless \mbox{\boldmathw}^{\ast} is the leading eigenvector of , the one corresponding to its largest eigenvalue. This means that Algorithm 8 finds the principal component of zero-mean data, and it also implies that the algorithm maximises over all with |\mbox{\boldmathw}|=1, see Section 6.3. Note that for zero-mean input data.
Now consider inputs with non-zero mean. In this case Algorithm 8 still finds the maximal-eigenvalue direction of . But for inputs with non-zero mean, this direction is different from the principal direction. Figure 10.4 illustrates this difference.
The Figure shows three data points in a two-dimensional input plane. The elements of \mathbb{C}^{\prime}=\langle\mbox{\boldmathx}\mbox{\boldmathx}^{\sf\small T}\rangle are
[TABLE]
with eigenvalues and eigenvectors
[TABLE]
So the the maximal-eigenvalue direction of is \mbox{\boldmathu}_{1}. To compute the principal direction of the data we must determine the data-covariance matrix , Equation (6.24). Its maximal-eigenvalue direction is \mbox{\boldmathu}_{2}, and this is the principal component of the data shown in Figure 10.4.
Oja’s rule can be generalised to determine principal components
of zero-mean input data using output neurons that compute y_{i}=\mbox{\boldmathw}_{i}\cdot\mbox{\boldmathx} for :
[TABLE]
This is called Oja’s M-rule [1]. For , Equation (10.13) simplifies to Oja’s rule.
10.2 Competitive learning
Oja’s -rule (10.13) results in neurons that are activated simultaneously. Any input usually causes several outputs to assume non-zero values at the same time. In Sections 4.5 and 7.1 we encountered the notion of a winning neuron where the weights are trained in such a way that each pattern activates only a single neuron, and different patterns activate different winning neurons. This allows to represent a distribution of input patterns with a neural network.
Unsupervised learning algorithms can categorise or cluster input data in this way: similar inputs are classified to belong to the same category, and activate the same winning neuron. This is called competitive learning [1]. Figure 10.5(a) shows an example, input patterns on the unit circle that cluster into two distinct clusters.
The idea is to find weight vectors \mbox{\boldmathw}_{i} that point into the direction of the clusters. To this end we take linear output units with weight vectors \mbox{\boldmathw}_{i}, . We feed a pattern from the distribution P_{\rm data}(\mbox{\boldmathx}) and define the winning neuron as the one that has minimal angle between its weight and the pattern vector . This is illustrated in Figure 10.5(b), where . Then only this weight vector is updated by adding a little bit of the difference \mbox{\boldmathx}-\mbox{\boldmathw}_{i_{0}} between the pattern vector and the weight of the winning neuron. The other weights remain unchanged:
[TABLE]
In other words, only the winning neuron is updated, \mbox{\boldmathw}_{\!i_{0}}^{\prime}=\mbox{\boldmathw}_{i_{0}}+\delta\!\mbox{\boldmathw}_{i_{0}}. Equation (10.14) is called competitive-learning rule.
The learning rule (10.14) has the following geometrical interpretation: the weight of the winning neuron is drawn towards the pattern . Upon iterating (10.14), the weight vectors are drawn to clusters of inputs. If the input patterns are normalised as in Figure 10.5,
the weights end up normalised on average, even though |\mbox{\boldmathw}_{i_{0}}|=1 does not imply that |\mbox{\boldmathw}_{i_{0}}+\delta\mbox{\boldmathw}_{i_{0}}|=1, in general. The algorithm for competitive learning is summarised in Algorithm 9. When weight and input vectors are normalised, then the winning neuron is the one with the largest scalar product \mbox{\boldmathw}_{i}\cdot\mbox{\boldmathx}. For linear output units y_{i}=\mbox{\boldmathw}_{i}\cdot\mbox{\boldmathx} (Figure 10.2) this is simply the unit with the largest output. Equivalently, the winning neuron is the one with the smallest distance |\mbox{\boldmathw}_{i}-\mbox{\boldmathx}|. Output units with \mbox{\boldmathw}_{i} that are very far away from any pattern may never be updated (dead units). There are several strategies to avoid this [1]. One possibility is to initialise the weights to directions found in the inputs. Also, how to choose the number of weight vectors is a matter of trial and error. Clearly it is better to start with too many rather than too few.
Finally, consider the relation between the competitive learning rule (10.14) and Oja’s rule (10.13). If we define
[TABLE]
then the rule (10.14) can be written in the form of Oja’s -rule:
[TABLE]
Equation (10.16) is reminiscent of Hebb’s rule (Chapter 2) with weight decay.
10.3 Self-organising maps
In order to analyse high-dimensional data it is often useful to map the high-dimensional input patterns to a low-dimensional output space, to obtain a low-dimensional representation of the input distribution. Principal-component analysis (Section 6.3) does just that. However, it does not necessarily preserve distance. To visualise clusters or other arrangements of the input patterns, similar patterns or patterns that are close in input space should be mapped to nearby points in output space, and patterns that are far apart should be mapped to outputs that are far from each other. Maps that achieve this are called semantic or topgraphic maps.
Moreover, principal-component analysis is a linear method. As explained in Section 6.3, it projects the data to the space spanned by the leading eigenvectors of the correlation matrix. In many cases, however, the data may not lie in a linear subspace, as illustrated in Figure 10.6. In order to project the data onto the non-linear principal manifold (solid line), a non-linear map is needed.
In neuroscience, the term topographic map refers to the relation between the spatial arrangement of stimuli and the activation patterns in certain parts of the mammalian brain. Similar patterns of visual stimuli on the retina, for instance, activate close-by regions in the visual cortex [157]. Other cognitive stimuli, auditory and sensory, are mapped in analogous ways. The complex neural networks in the mammalian cortex contain large numbers of such maps, arranged in a hierarchical fashion. They represent local stimuli in terms of spatially localised neural activation. How did this complex structure arise? One possibility is that the mappings are coded in the genetic sequence, that the connections are hard wired, so to speak. However, it is observed that such maps can change over time [158], leading to the hypothesis that they are learned, and that the our DNA merely encodes a set of fairly simple learning rules.
This motivated Kohonen [158, 159] and others to propose and analyse learning rules for topographic maps. The term self-organising map [18, 160] emphasises that the mapping develops in response to the stimuli it maps, that it learns in an unsupervised fashion.
Kohonen’s model for a non-linear self-organising map relies on an ordered array of output neurons, as illustrated in Figure 10.7. The map learns to activate nearby output neurons for similar inputs. This is achieved using a competitive learning rule, similar to the learning rule (10.14) described in the previous Section. In order to represent the proximity or similarity of inputs, the rule is endowed with the notion of distance in the output array, by updating not only the winning neuron, but also those that are neighbours in the output array. To this end one replaces the competitive-learning rule (10.14) by
[TABLE]
where is the index of the winning neuron, the one with weight vector closest to the input . The neighbourhood function function depends on the distance of the neurons and in the output array. The neighbourhood function has a maximum at and decreases as the distance between and increases. One possibility is to assign decreasing values to for nearest neighbours, next-nearest neighbours, and so forth (Figure 10.8).
Another possibility is to use a Gaussian function of the Euclidean distance |\mbox{\boldmathr}_{i}-\mbox{\boldmathr}_{i_{0}}| in the output array [1]:
[TABLE]
Here \mbox{\boldmathr}_{i} is the position of neuron in the output array (Figure 10.7). Different normalisations of the Gaussian [2] can be subsumed in different learning rates.
Kohonen’s rule has two parameters: the learning rate , and the width of the neighbourhood function.
Usually one adjusts these parameters as the learning proceeds. Typically one begins with large values for and (ordering phase), and then reduces these parameters as the elastic network evolves (convergence phase): quickly at first and then in smaller steps, until the algorithm converges [158, 1, 2].
According to Equations (10.17) and (10.18), similar patterns activate nearby neurons in output space, and their weight vectors change in similar ways. Kohonen’s rule drags the winning weight vector \mbox{\boldmathw}_{\!i_{0}} towards , just as the competitive learning rule (10.14), but it also drags the neighbouring weight vectors along. Figure 10.9 illustrates a geometrical interpretation of Kohonen’s rule [18]. We can think of the weight vectors as pointing to the nodes of an elastic net that has the same layout as the output array. As one feeds patterns from the input distribution, the weights are updated, causing the nodes of the network to move. This changes the shape of the elastic network until it resembles the shape defined by the distribution of input patterns. Figure 10.6 shows another example where the dimensionality of the output array (one-dimensional) is lower than that of the input space (two-dimensional). The algorithm finds a non-linear approximation to the data, the principal manifold. As opposed to the principal direction in principal-component analysis, the principal manifold need not be linear. Therefore it can approximate the data more precisely, leading to a smaller residual variance (Exercise 10.7).
In summary, Kohonen’s algorithm learns by distributing the weight vectors of the output neurons to reflect the distribution of input patterns. In general this works well, but problems occur at the boundaries. Why this happens is quite clear (Figure 10.9): since the density of patterns outside the parallelogram is low, the elastic network cannot be drawn very close to the boundary. To analyse how the boundaries affect learning for Kohonen’s rule, consider the steady-state condition
[TABLE]
This condition is more complicated than it looks at first sight, because depends on the weights and on the patterns, as mentioned above. The steady-state condition (10.19) is very difficult to analyse in general. One of the reasons is that global geometric information is difficult to learn. It is usually much easier to learn local structures. This is particularly true in the continuum limit where we can analyse local learning progress using Taylor expansions.
The analysis of condition (10.19) in the continuum limit
is due to Ritter and Schulten [161], and it is described in detail by Hertz, Krogh, and Palmer [1]. One assumes that there is a very dense network of output neurons, so that one can approximate i\to\mbox{\boldmathr}, i_{0}\to\mbox{\boldmathr}_{0}, \mbox{\boldmathw}_{i}\to\mbox{\boldmathw}(\mbox{\boldmathr}), h(i,i_{0})\to h\big{(}\mbox{\boldmathr}-\mbox{\boldmathr}_{0}(\mbox{\boldmathx})\big{)}, and \tfrac{1}{T}\sum_{t}\to\int{\rm d}\mbox{\boldmathx}P_{\rm data}(\mbox{\boldmathx}). In this continuum limit, Equation (10.19) reads
[TABLE]
This is a condition for the steady-state learning outcome, the function \mbox{\boldmathw}^{\ast}(\mbox{\boldmathr}).
In the continuum limit the position \mbox{\boldmathr}_{0}(\mbox{\boldmathx}) of the winning neuron in the output array for pattern is given by
[TABLE]
We use this relation to write Equation (10.20) as:
[TABLE]
Equation (10.21) defines a mapping \mbox{\boldmathr}_{0}(\mbox{\boldmathx}) from input space to output space, the self-organising map (Figure 10.7). Assuming that this mapping is one-to-one, we change integration variable from to \mbox{\boldmathr}_{0}:
[TABLE]
where Q(\mbox{\boldmathr}_{0})\equiv P_{\rm data}\big{(}\mbox{\boldmathx}(\mbox{\boldmathr}_{0})\big{)}, and where the determinant represents the volume element of the variable transformation. Using Equation (10.21), the Jacobian of the transformation has elements
[TABLE]
The neighbourhood function is sharply peaked at \mbox{\boldmathr}=\mbox{\boldmathr}_{0},
and this makes it possible to evaluate the steady-state condition (10.23) approximately,
expanding the integrand in \delta\mbox{\boldmathr}=\mbox{\boldmathr}_{0}-\mbox{\boldmathr}, assuming that \mbox{\boldmathw}^{\ast}(\mbox{\boldmathr}) is a smooth function. This is illustrated in Figure 10.10, for one-dimensional inputs and outputs. We consider this special case not only to simplify the notation, but also because it is one of the few cases that admit mathematical analysis (Exercise 10.9). Expanding as shown in Figure 10.10 yields
[TABLE]
The other factors in Equation (10.23) are expanded in a similar way:
[TABLE]
Inserting these expressions into Equation (10.20), discarding terms of order higher than , and changing the integration variable to , one finds
[TABLE]
where we introduced the short-hand notation , dropped the asterisk, and used that the neighbourhood function (10.18) is symmetric, .
Since the integral in Equation (10.27) is non-zero, we must either have
[TABLE]
The first solution can be excluded because it corresponds to a singular weight distribution that does not contain any geometrical information about the input distribution . The second solution gives
[TABLE]
In other words, , and this means that . The density of output weights can be computed as
[TABLE]
where is the Dirac -function [162]. Changing variables in the -function
[TABLE]
and assuming that the function is one-to-one, one finds
[TABLE]
This tells us that the self-organising map learns the input distribution in the following way: the distribution of output weights in the steady state reflects the distribution of input patterns. Equation (10.32) shows that the two distributions are not equal (equality would have been a perfect outcome). The distribution of weights is instead proportional to \big{[}P_{\rm data}(w)\big{]}^{\frac{2}{3}}. Little is known in higher dimensions, but the general idea is that the elastic network has difficulties reaching the corners and edges of the domain where the input distribution is non-zero.
The output of a self-organising map can be interpreted in different ways. For a low-dimensional inputs and outputs, one can simply plot the map w^{\ast}(\mbox{\boldmathr}), as in Figure 10.7. Dense regions of weights point to regions in input space with a high density of inputs.
Often the output dimension is taken to be much lower than the dimension of input space. In this case the self-organising map performs non-linear dimensionality reduction, and it can be used to find clusters in high-dimensional input data [163]. The analysis proceeds in two steps. First, one runs Kohonen’s algorithm until the map has converged to a steady state. Second, one feeds all inputs into the net, and for each input one determines the location of the winning neuron in the output array. The spatial activation patterns in the output array represent clusters of similar inputs. This is illustrated in Figure 10.11, which shows how a self-organising map represents handwritten digits from the MNIST data set. To reveal the semantic map, the Figure labels clusters of outputs that correspond to the same digits (as determined by the labels in the training set). We see that the self-organising map groups the same digits together, but it has some difficulty distinguishing the digits 3 and 8, and also 4 and 9.
10.4 -means clustering
Sections 10.2 and 10.3 described different ways of finding clusters in input data. In particular, it was shown how self-organising maps can find clusters in high-dimensional input data, and represent them in a low-dimensional, non-linear projection. -means clustering [2] is an alternative unsupervised-learning algorithm for finding clusters in the input data. Let us compare and contrast this algorithm with Kohonen’s self-organising map. The goal is to cluster input patterns \mbox{\boldmathx}^{(\mu)}, into clusters. Usually is much smaller than the number of inputs, , and than the input dimension .
A solution of the clustering task is a mapping that associates each input \mbox{\boldmathx}^{(\mu)} with one of the clusters . The function is determined by minimising the energy function
[TABLE]
The second sum is over all values of that satisfy . The vector \mbox{\boldmathw}_{k} becomes the average of all pattern vectors in cluster , and the expression in the parentheses is the variance associated with this cluster:
[TABLE]
In other words, measures the sum of the cluster variances . A solution to the clustering problem corresponds to a local minimum of . To determine the cluster vectors \mbox{\boldmathw}_{k} and the corresponding variances , one starts from an initial guess for . For each cluster, one begins by adjusting \mbox{\boldmathw}_{k} to minimise the cluster variance:
[TABLE]
In a second step, one optimises the encoding function
[TABLE]
given the vectors \mbox{\boldmathw}_{k}. These steps are repeated until a satisfactory solution is found. The solution is not unique, usually the algorithm converges to a local minimum of . In practice one should try different random initialisations to find the best local minimum.
All three algorithms, competitive learning, the self-organising map, and -means clustering move weight vectors towards clusters in input space. A difference between the self-organising map and the other two algorithms is that
the self-organising map uses a neighbourhood function (so that similar inputs activate close-by neurons in the output array), and updates their weight vectors in similar fashion. In this way, a self-organising map with a large output array can find a smooth parameterisation of the principal manifold. If we shrink the neighbourhood function to the centre point in Figure 10.8, all geometric information is lost, and the self-organising map becomes equivalent to competitive learning (Algorithm 9). Essentially, competitive learning and -means clustering are sequential and batch versions of the same algorithm [1]. So the self-organising map becomes equivalent to -means clustering when the neighbourhood range tends to zero.
10.5 Radial basis functions
Problems that are not linearly separable can be solved by perceptrons with hidden layers,
as we saw in Chapter 5. Figure 5.13(b), for example, shows a piecewise linear decision boundary parameterised by hidden neurons.
Separability can also be achieved by a non-linear transformation of input space. Figure 10.13 shows how the XOR problem can be transformed into a linearly separable problem by the transformation
[TABLE]
The Figure shows the non-separable problem in the - plane, and in the new coordinates and .
The problem is homogeneously linearly separable in the - plane. We can solve it by a single McCulloch-Pitts neuron with weights and zero threshold, parameterising the decision boundary as \mbox{\boldmathW}\cdot\mbox{\boldmathu}(\mbox{\boldmathx})=0.
It is even better to map the patterns (non-linearly) to a space of higher dimension, because Cover’s theorem (Section 5.4) says that it is easier to separate the patterns there: consider a set \mbox{\boldmathu}(\mbox{\boldmathx})=[u_{1}(\mbox{\boldmathx}),\ldots,u_{m}(\mbox{\boldmathx})]^{\sf T} of polynomial functions of finite order that embed -dimensional input space in an -dimensional space. Then the probability that a problem with points \mbox{\boldmathx}^{(\mu)} in -dimensional input space is separable by a polynomial decision boundary is given by [Equation (5.29)] [72, 2]. Note that this probability is independent of the dimension of input space.
The question is of course how to find the non-linear mapping \mbox{\boldmathu}(\mbox{\boldmathx}). One possibility is to use radial basis functions. The idea is to parameterise the functions u_{j}(\mbox{\boldmathx}) in terms of weight vectors \mbox{\boldmathw}_{j}, and to use an unsupervised-learning algorithm to find weights that separate the input data. A common choice [2] is to use radial basis functions of the form:
[TABLE]
Note that these functions are not of the finite-order polynomial form that was assumed above. So strictly speaking we cannot invoke Cover’s theorem. In practice the mapping u_{j}(\mbox{\boldmathx}) works nevertheless quite well. The parameters parameterise the widths of the radial basis functions. In the simplest version of the algorithm they are set to unity. Hertz, Krogh, and Palmer [1] discuss radial basis-function networks with normalised radial basis functions
[TABLE]
Other choices for radial basis functions are given by Haykin [2].
Figure 10.14 shows a radial basis-function network for and . The four neurons in the hidden layer stand for the four radial basis functions (10.38) that map the inputs to four-dimensional -space. The network looks like a perceptron (Chapter 5). But here the hidden layers work in a different way. Perceptrons have hidden McCulloch-Pitts neurons that compute non-local outputs \sigma(\mbox{\boldmathw}_{j}\cdot\mbox{\boldmathx}-\theta). The output of radial basis functions u_{j}(\mbox{\boldmathx}), by contrast, is localised in input space [Figure 10.15 (left)]. We saw in Section 7.1 how to make localised basis functions out of McCulloch-Pitts neurons with sigmoid activation functions , but
one needs two hidden layers to do that [Figure 10.15 (right)].
Radial basis functions produce localised outputs with a single hidden layer, they divide up input space into localised regions, each corresponding to one radial basis function. Imagine for a moment that we have as many radial basis functions as input patterns. In this case we can simply take \mbox{\boldmathw}_{\nu}=\mbox{\boldmathx}^{(\nu)} for . The linear output in Figure 10.14 computes O^{(\mu)}=\mbox{\boldmathW}\cdot\mbox{\boldmathu}(\mbox{\boldmathx}^{(\mu)}), and so the classification problem in -space takes the form
[TABLE]
with U_{\mu\nu}=u_{\nu}(\mbox{\boldmathx}^{(\mu)}). If all patterns are pairwise different, \mbox{\boldmathx}^{(\mu)}\neq\mbox{\boldmathx}^{(\nu)} for , then the matrix is invertible [2]. In this case the solution of the classification problem reads
[TABLE]
where is the symmetric matrix with elements .
In practice one can get away with fewer radial basis functions by choosing their weights to point in the directions of clusters of input data. To this end one uses unsupervised competitive learning (Algorithm 10), where the index of the winning neuron is defined to be the one with largest . How are the widths determined? The width of radial basis function u_{j}(\mbox{\boldmathx}) is taken to be equal to the minimum distance between \mbox{\boldmathw}_{j} and the centers of the surrounding radial basis functions. Once weights and widths of the radial basis functions are found, the weights of the output neuron are determined by minimising
[TABLE]
with respect to . This works even if is not invertible. An approximate solution can be obtained by stochastic gradient descent on , keeping the parameters of the radial basis functions fixed. Cover’s theorem indicates that the problem is more likely to be separable if the embedding dimension is higher.
Radial basis-function networks are similar to the perceptrons described in Chapters 5 to 7, in that they are feed-forward networks designed to solve classification problems. A fundamental difference is that the parameters of the radial basis functions are determined by unsupervised learning, whereas perceptrons are trained using supervised learning for all units. While McCulloch-Pitts neurons compute weights to minimise their output from given targets, the radial basis functions compute weights by maximising the output as a function of . The algorithm for finding the weights of the radial basis functions is summarised in Algorithm 10. Further, as opposed to the deep networks from Chapter 7, radial basis-function networks have only one hidden layer, and a linear output neuron. In summary, radial basis-function networks learn using a hybrid scheme: unsupervised learning for the parameters of the radial basis functions, and supervised learning for the weights of the output neuron.
10.6 Autoencoders
Multi-layer perceptrons, layered feed-forward networks, were developed for supervised learning, as described in Part II. Such layouts can also be used for unsupervised learning. Examples are autoencoders and generative adversarial networks.
Autoencoders employ layered feed-forward networks for unsupervised learning of an unlabeled data set of input patterns, using the inputs as targets, \mbox{\boldmatht}^{(\mu)}=\mbox{\boldmathx}^{(\mu)}. The layout is illustrated in Figure 10.16. The network consists of two main parts, an encoder (on the left), and a decoder (on the right). The encoder consists for instance of several fully connected or convolutional layers and maps the inputs to a bottleneck layer with a small number of neurons, significantly smaller than the input dimension, . We denote the states of the bottleneck neurons by . The encoder corresponds to a non-linear mapping \mbox{\boldmathz}=\mbox{\boldmathf}_{\!\rm e}(\mbox{\boldmathx}). The decoder maps the bottleneck (or latent) variables back to the inputs, \mbox{\boldmathx}=\mbox{\boldmathf}_{\!\rm d}(\mbox{\boldmathz}). One adjusts weights and thresholds by backpropagation until the network learns to approximate the inputs as
[TABLE]
The energy function reads:
[TABLE]
where |\mbox{\boldmathx}|^{2}=\mbox{\boldmathx}^{\sf T}\mbox{\boldmathx}. In other words, the autoencoder learns the identity function. The point is that the identity is represented in terms of two non-linear functions, the encoder \mbox{\boldmathf}_{\!\rm e} and the decoder \mbox{\boldmathf}_{\!\rm d}. While the identity function is trivial, the encoding and decoding functions need not be. The bottleneck ensures that the network does not simply learn \mbox{\boldmathf}_{\!\rm e}(\mbox{\boldmathx})=\mbox{\boldmathf}_{\!\rm d}(\mbox{\boldmathx})=\mbox{\boldmathx}.
The latent variables may encode interesting properties of the input patterns. If the number of neurons is much smaller than the number of pattern bits, as indicated by the term bottleneck, the encoder is a low-dimensional (compressed) representation of the input data. In this way autoencoders can perform non-linear dimensionality reduction, like self-organising maps (Section 10.3). If both encoder and decoder are linear functions with zero thresholds, then H=\frac{1}{2}\sum_{\mu}|\mbox{\boldmathx}^{(\mu)}-\mathbb{W}_{\rm d}\mathbb{W}_{\rm e}\mbox{\boldmathx}^{(\mu)})|^{2}. In this case, z_{1}(\mbox{\boldmathx}),\ldots,z_{M}(\mbox{\boldmathx}) are simply the first principal components
of zero-mean input data [164] (Exercise 10.14).
Sparse autoencoders [165] have a large number of neurons in the bottleneck, possibly more than the number of pattern bits. But only a small number of bottleneck neurons are allowed to be active at the same time. The idea is that sparse representations of input data are more robust than dense ones, and generalise more reliably. At least high-dimensional but sparse representations of binary classification problems are more likely to be linearly separable (Section 5.4). There are different ways of enforcing sparsity, for instance using - or -regularisation (Section 7.6.1). An alternative [165] is to ensure that the average activation of each bottleneck neuron with sigmoid activation function,
[TABLE]
remains small. This is achieved by adding the term
[TABLE]
to the energy function, with Lagrange multiplier (Section 6.3). This term penalises any deviation of from , where is a sparsity parameter. Each term in the sum is non-negative and vanishes when for all , because each term can be interpreted as the Kullback-Leibler divergence (Section 4.4) between two Bernoulli distributions with parameters and .
Variational autoencoders [166, 167, 168] have layouts similar to the one shown schematically in Figure 10.16, but their purpose is quite different. Variational autoencoders are generative models (Section 4.5): just like restricted Boltzmann machines they approximate a data distribution of inputs P_{\rm data}(\mbox{\boldmathx}), and allow to sample from it. As an example consider the MNIST data set of handwritten digits. The patterns define a data distribution that encodes the properties of the digits in terms of covariances and higher-order correlations. The question is how to generate new digits from this distribution, different from those in the data set, yet with their defining properties. In other words, how can a machine learn to generate images that look like handwritten digits?
The idea of variational autoencoders is to represent the data distribution in terms of a Gaussian distribution P_{\rm L}(\mbox{\boldmathz}) of latent variables , using the fact that one can approximate any given data distribution P_{\rm data}(\mbox{\boldmathx}) in terms of P_{\rm L}(\mbox{\boldmathz}) by a suitable non-linear transformation \mbox{\boldmathf}(\mbox{\boldmathz}) [167]. Variational autoencoders are trained not unlike neural networks, but an essential difference to the algorithms described in Part II is that variational autoencoders learn probabilities rather than deterministic input-output mappings.
Given the Gaussian distribution P_{\rm L}(\mbox{\boldmathz}) of the latent variables, the goal is to maximise the log-likelihood (Section 4.5):
[TABLE]
Here P(\mbox{\boldmathx}|\mbox{\boldmathz}) is the probability to generate given . In the simplest case, this distribution is assumed to be Gaussian with mean \mbox{\boldmath\mu}_{P}(\mbox{\boldmathz})=\mbox{\boldmathf}(\mbox{\boldmathz}), and with correlation matrix \mathbb{C}_{P}{(\mbox{\boldmathz})}. The decoder represents these functions in terms of a multilayer perceptron or a convolutional neural network. Weights and thresholds are determined to maximise by gradient ascent. To this end we must find an efficient way of computing and its gradients. One possibility is Monte-Carlo sampling (Section 4.2), but this is not very efficient because most values of drawn from P_{\rm L}(\mbox{\boldmathz}) result in unlikely patterns , with only negligible contributions to . To get around this problem, one needs to know which values of are likely to produce a given pattern . The idea is to learn a second approximate distribution Q(\mbox{\boldmathz}|\mbox{\boldmathx}) of given . We can think of Q(\mbox{\boldmathz}|\mbox{\boldmathx}) as an encoder. So Q(\mbox{\boldmathx}|\mbox{\boldmathz}) corresponds to the encoder \mbox{\boldmathf}_{\!\rm e} discussed above, while P(\mbox{\boldmathz}|\mbox{\boldmathx}) corresponds to the decoder \mbox{\boldmathf}_{\!\rm d}. An important difference is that and are probabilities, not deterministic functions.
To determine a good approximation Q(\mbox{\boldmathz}|\mbox{\boldmathx}), we minimise the difference between Q(\mbox{\boldmathz}|\mbox{\boldmathx}) and the unknown exact distribution P(\mbox{\boldmathz}|\mbox{\boldmathx}),
[TABLE]
A first trick is to rewrite this expression using Bayes’ theorem [29],
[TABLE]
This gives:
[TABLE]
The second trick is to note that the l.h.s of Equation (10.50) is a suitable target function to maximise. We want to maximise subject to the constraint that the unknown function Q(\mbox{\boldmathz}|\mbox{\boldmathx}) approximates the probability of encoding the pattern . Usually one takes Q(\mbox{\boldmathz}|\mbox{\boldmathx}) to be a Gaussian with mean \mbox{\boldmath\mu}_{Q} and correlation matrix . The task is then to determine the functions \mbox{\boldmath\mu}_{P}{(\mbox{\boldmathz})},\mathbb{C}_{P}{(\mbox{\boldmathz})},\mbox{\boldmath\mu}_{Q}{(\mbox{\boldmathz})}, and \mathbb{C}_{Q}{(\mbox{\boldmathz})} by adjusting the weights of two neural networks, encoder and decoder.
This task, maximising the r.h.s. of Equation (10.50), is not as straightforward as it may seem, because the target function (10.50) involves an average over the distribution which in turn depends on the weights. The question is how to move the derivative inside the average , to obtain an unbiased expression for the weight updates . In other words, the goal is to ensure that the average weight increments are proportional to the gradients of the target function (10.50). A related problem occurs when training binary stochastic neurons (Section 11.1). Similarities and differences are described in Ref. [169].
One solution is to use stochastic backpropagation [168]. In its simplest form, this algorithm makes use of a relation for the gradient of the average of a test function F(\mbox{\boldmathz}) with a Gaussian probability Q(\mbox{\boldmathz};w_{ij}) that depends on the weights :
[TABLE]
Here and are the gradient and the Hessian of the function F(\mbox{\boldmathz}). A challenge is that these derivatives tend to be difficult to compute reliably. Suitable approximations are described in Ref. [168]. The expression inside the average in Equation (10.51) is the unbiased weight increment. Iterating the learning rule allows one to determine the parameters of P(\mbox{\boldmathx}|\mbox{\boldmathz}) and Q(\mbox{\boldmathz}|\mbox{\boldmathx}). This allows tio efficiently sample by sampling the latent variables and then applying the decoder.
Generative adversarial networks [170] are generative models based on learning rules similar to that described above for variational autoencoders, but there are some differences in detail. Generative adversarial networks consist of two multilayer perceptrons, a generator and a discriminator. The generative network produces new outputs from a given data distribution (fakes), and the task of the discriminator is to classify these outputs into two classes: real or fake data. Generator and discriminator are trained together. The weights of the generator are adjusted to maximise the classification error of the discriminator, while those of the discriminator are trained to minimise this error [171].
10.7 Summary
The unsupervised-learning algorithms described in Sections 10.1 and 10.2 are based on Hebb’s rule. These algorithms can learn different features of unlabeled input data: they can detect the familiarity of inputs, perform principal-component analysis, and identify clusters in the input data.
Self-organising maps also rely on Hebb’s rule. An important difference is that the outputs are arranged in an array, and that output neurons that are close-by in the output array are updated in similar ways. Self-organising maps can therefore represent topographic and semantic maps, where close-by or similar inputs are mapped to nearby outputs. When the dimension of the output array is much lower than the input dimension, self-organising maps perform non-linear dimensional reduction.
Radial basis-function networks are classifiers, just like multilayer perceptrons. Their output neurons are trained in the same way, using labeled input data. However, the decision boundaries of radial basis-function networks are polynomial functions (not just hyperplanes), and their parameters are determined by unsupervised learning.
Autoencoders are multilayer perceptrons. They can learn to encode non-linear features of unlabeled input data by using the input patterns as targets. Finally, generative adversarial networks do not require labeled inputs, so they can be considered unsupervised-learning machines. They are used to generate synthetic data in order to expand training sets for supervised learning, and pose an ethical dilemma because they can be used to generate deep fakes [172], manipulated videos where someone’s facial expression and speech are replaced by another person’s.
In summary, the simple algorithms described in this Chapter provide a proof of concept: how machines can learn without labels.
10.8 Further reading
The primary source for Sections 10.1 and 10.2 is the book by Hertz, Krogh, and Palmer [1]. A good reference for self-organising maps is Kohonen’s book [159]. Radial basis-function networks are discussed by Haykin in Chapter 5 of his book [2]. It has been argued that radial-basis function networks do not generalise as well as perceptrons do [173]. To solve this problem, Poggio and Girosi [174] suggested to determine the parameters \mbox{\boldmathw}_{j} of the radial basis function by supervised learning, using stochastic gradient descent.
Autoencoders can generate non-linear, low-dimensional representations of an input distribution. The relation to principal component analysis is discussed in Refs. [164, 175].
The recommended introduction to variational autoencoders is the tutorial by Doersch [167]. He also mentions that the underlying mathematics for variational autoencoders is similar to that of Helmholtz machines (Section 4.7), although the two machines learn in quite different ways.
Variational autoencoders are used for a number of different purposes. Ref. [176] suggests to employ a variational autoencoder for active learning (Section 7.8). The idea is to represent the input distribution in terms of lower-dimensional latent variables, and to use -means clustering (Section 10.4) to identify groups of patterns that should be labeled. Variational autoencoders have also been used for outlier detection [177] and language generation [178].
Chapter 11 Reinforcement learning
Supervised learning requires labeled data, where each input comes with a target the network is supposed to learn. Unsupervised learning, by contrast, does not require labeled data. Reinforcement learning
lies between these extremes. The term reinforcement describes the principle of learning by means of a reward function. This function assigns penalties or rewards to the network output, depending on how the output relates to the learning goal. For a neural network with a vector of outputs, the reward function could be
[TABLE]
The goal is to learn to produce outputs that receive a reward more frequently than those that trigger a penalty. We say that rewarded outputs are reinforced. The feedback may be random, given by a distribution initially unknown to the network.
The reward function reflects the learning goal.The training process as well as the learning outcome depend crucially on this function. Suppose one replaces the reward function (11.1) by the more lenient alternative: if at least one output is correct, and otherwise. Naturally this leads to more errors, possibly not a good idea if the goal is to teach a robot to fly.
One distinguishes two different types of reinforcement problems, associative and non-associative tasks [179]. An example for a non-associative task is the -armed bandit problem [16]. Imagine slot machines with different reward distributions, initially unknown to the player. Given a finite amount of money, the question is in which order to play the machines so as to maximise the overall profit. The dilemma is whether to stick with a machine that yields a good reward, or whether to try out other machines that may yield a low reward initially, but could give much higher rewards eventually (exploit-versus-explore dilemma). In this type of problem, the player receives only the reinforcement signal, no other inputs. In associative tasks, by contrast, the agent receives inputs, or stimuli, and it should learn to associate with each stimulus the output that yields the highest reward. Such tasks occur for instance in behavioural psychology, where the problem is to discriminate between different stimuli, and to associate the right behaviour with each stimulus.
In general, such associative tasks can be described as sequential decision processes (Figure 11.1), where an agent explores a sequence of states \mbox{\boldmaths}_{0},\mbox{\boldmaths}_{1},\mbox{\boldmaths}_{2},\ldots through a sequence of actions \mbox{\boldmatha}_{0},\mbox{\boldmatha}_{1},\mbox{\boldmatha}_{2},\ldots. Consider for instance a motile microorganism in the turbulent ocean that should swim to the water surface as quickly as possible [180]. It determines its state by observing the local environment. The microorganism might measure local strain and vorticity of the flow. The environment provides a reinforcement signal (the distance to the surface for example), and the organism determines which action to take, given its state and the reinforcement signal. Should it turn, stop to swim, or accelerate? The organism learns to associate actions with certain states that maximise the reward. This sounds quite similar to associating optimal outputs with stimuli. A conceptual difference is that the action of the agent modifies the environment: its actions take it to a different place in the turbulent flow, with different vorticity and different strain. A second point is that the reward to an action may not be immediate. In this case the challenge is to credit actions that optimise the expected future reward, given the information collected so far. This is the credit-assignment problem [181].
There are two different kinds of associative tasks: continuous and episodic ones. In continuous tasks, the intertwined sequences of states and actions have no natural end, so that one must either terminate the sequence in an ad-hoc fashion or introduce a weighting factor to ensure that the expected future reward remains finite. In episodic tasks, by contrast, the learning is divided into episodes that terminate after a finite number of steps. An example is to learn a strategy for winning a board game. In this case, each episode corresponds to a round of the game, and the reward is incurred at the end of each episode. The number of steps per round, the episode length , may vary from round to round. In order to estimate the expected reward one usually needs many episodes. A second, very simple example is the stimulus problem described above, where the states (stimuli) are independent from the actions. Each episode consists of only one step, so that . In response to a randomly chosen state \mbox{\boldmaths}_{0}, the agent learns to perform the action \mbox{\boldmatha}_{0} that maximises the immediate reward. This can be achieved by the associative reward-penalty algorithm (Section 11.1). It uses stochastic neurons with weights that are trained by gradient ascent to maximise the expected immediate reward.
To estimate the expected future reward when , one must use a different method, usually temporal difference learning. It allows to estimate the expected future reward, after steps, by breaking up the learning into time steps . The idea is that it is better to adjust the prediction of a future reward as one iterates, rather than waiting for iterations before updating the prediction. In temporal difference learning one expresses the reward at time in terms of differences at time steps and [182].
Temporal difference learning builds up a lookup table that summarises the best actions for each state, the -table Q(\mbox{\boldmaths},\mbox{\boldmatha}). Given states and possible actions , Q(\mbox{\boldmaths},\mbox{\boldmatha}) is a table. Its elements contain the expected future reward for each state-action pair. To implement temporal difference learning, one must adopt a policy or strategy. It specifies for each state which action is taken. In total there are possibilities of assigning actions to states. When there are many states and actions it quickly becomes impossible to determine the best strategy by simple sampling, because there are too many to consider.
The advantage of temporal difference learning and related algorithms is that they do not rely on the evaluation of all possible policies. Instead, the policy is updated using iterated estimates of the table. We write for the estimate at time step . There are different ways of deriving a policy from a -table. The greedy policy is a deterministic policy, it corresponds to choosing the action that corresponds to the largest -element in a given row of the -table: \mbox{\boldmatha}=\mbox{argmax}_{\footnotesize\mbox{\boldmatha}^{\prime}}Q_{t}(\mbox{\boldmaths},\mbox{\boldmatha}^{\prime}). This policy maximises the current estimate of the future reward.
Stochastic policies are often better, in particular for non-stationary or stochastic environments, because they allow the agent to explore potentially better alternatives. Also, for a deterministic environment, a deterministic policy may lead to cycles. This can be avoided with a stochastic policy. One example for a stochastic policy is the -greedy policy. With probability , it chooses the greedy action \mbox{\boldmatha}=\mbox{argmax}_{\footnotesize\mbox{\boldmatha}^{\prime}}Q_{t}(\mbox{\boldmaths},\mbox{\boldmatha}^{\prime}), but with a small probability it takes a suboptimal action. Another example is the softmax policy, where argmax is replaced by the softmax function (Section 7.5). The softmax policy can handle actions described in terms of continuous variables.
In general the policy can change as the algorithm is iterated. A common choice is to reduce the parameter in the -greedy policy as one iterates, so that the algorithm converges to the optimal deterministic policy.
-learning is an approximation to the temporal difference algorithm. In -learning, the -table is updated assuming that the agent always follows the greedy policy, even though it might actually follow a different policy. -learning allows agents to learn to play strategic games [183]. A simple example is the game of tic-tac-toe (Section 11.3). Games such as chess or go require to keep track of a very large number of states, so large that -learning in its simplest form becomes impractical. An alternative is to represent the state-action mapping by a deep neural network [184].
11.1 Associative reward-penalty algorithm
The associative reward-penalty algorithm uses stochastic neurons that are trained to maximise the average immediate reward.
In Chapters 5 to 9 the output neurons were deterministic functions of their inputs. For reinforcement learning, by contrast, it is better to use stochastic neurons. The idea is the same as in Chapters 3 and 4: stochastic neurons can explore a wider range of possible states which may in the end lead to a better solution. The state of neuron is given by the stochastic update rule (3.1):
[TABLE]
where b_{i}=\mbox{\boldmathw}_{i}\cdot\mbox{\boldmathx} is the local field (no thresholds), and . Recall that the parameter is the noise level. Since the outputs can assume only two values, , Equation (11.2) describes a binary stochastic neuron.
To illustrate the associative reward-penalty algorithm for a single binary stochastic neuron, consider an agent experiencing different stimuli drawn with equal probability from a distribution of inputs. Upon receiving stimulus , the stochastic neuron outputs either or . Given and , the environment provides a stochastic reward r(\mbox{\boldmathx},y)=\pm 1 drawn from a reward distribution p_{\rm reward}(\mbox{\boldmathx},y):
[TABLE]
The goal is to adjust the weights so that the neuron produces outputs that are rewarded with high probability. Figure 11.2(a) shows an example with just two stimuli, {\mbox{\boldmathx}^{(1)}}=[1,0]^{\sf T} and {\mbox{\boldmathx}^{(2)}}=[1,1]^{\sf T}. The numerical values of p_{\rm reward}(\mbox{\boldmathx},y) indicate that the expected reward is maximised when the neuron outputs in response to {\mbox{\boldmathx}^{(1)}}, and in response to {\mbox{\boldmathx}^{(2)}}. Since {\mbox{\boldmathx}^{(1)}} and {\mbox{\boldmathx}^{(2)}} occur with equal probability, the maximal expected reward is
[TABLE]
Here we used \langle r(\mbox{\boldmathx},y)\rangle_{\rm reward}=p_{\rm reward}(\mbox{\boldmathx},y)-[1-p_{\rm reward}(\mbox{\boldmathx},y)], as well as the numerical values for p_{\rm reward}(\mbox{\boldmathx},y) given in Figure 11.2(a).
Figure 11.2 (b) shows the contingency space [185] of the problem, representing the inputs in a plane with coordinates p_{\rm reward}(\mbox{\boldmathx},+1) and p_{\rm reward}(\mbox{\boldmathx},-1). It is easier to learn to associate the correct output with inputs that lie in the shaded regions where p_{\rm reward}(\mbox{\boldmathx},+1)>\tfrac{1}{2} and p_{\rm reward}(\mbox{\boldmathx},-1)<\tfrac{1}{2}, or vice versa. In this case one can solve the problem by fixing and sampling p_{\rm reward}(\mbox{\boldmathx},+1) for all . If p_{\rm reward}(\mbox{\boldmathx},+1)>\tfrac{1}{2} then is the optimal output for , otherwise it is . This strategy cannot be used outside the shaded region. For example, if both reward probabilities are larger than one half, it is necessary to sample both p_{\rm reward}(\mbox{\boldmathx},-1) and p_{\rm reward}(\mbox{\boldmathx},+1) sufficiently often in order to determine which one is larger: one must find the greater of two goods according to Barto [185]. This illustrates the fundamental dilemma of reinforcement learning: an output that appears at first to yield a high reward may not be the optimal one in the long run. To find the optimal output it is necessary to estimate both reward probabilities precisely. This means that one must try all possible outputs frequently, not only the one that appears to be optimal at the moment.
To derive a learning rule we need a cost function. One possibility is to use the average of the immediate reward for a given stimulus ,
[TABLE]
Here is an average over the response of the environment, determined by the reward distribution p_{\rm reward}(\mbox{\boldmathx},\mbox{\boldmathy}). It is assumed that the reward distribution is stationary. Furthermore,
[TABLE]
is the probability that the network produces the output \mbox{\boldmathy}=[y_{1},\ldots,y_{M}]^{\sf T} given the local field .
To find the maximum of one uses gradient ascent on , analogous to maximising the log-likelihood for Boltzmann machines (Section 4.4), and to gradient descent on the energy function for perceptrons in supervised learning (Chapter 6).
The gradient is computed by applying the chain rule, as usual. The calculation is similar to the one for Boltzmann machines (Chapter 4). After some algebra (Exercise 11.2) one finds for given and that the derivative of P(\mbox{\boldmathy}|\mbox{\boldmathx}) with respect to equals P(\mbox{\boldmathy}|\mbox{\boldmathx})\,\beta[y_{m}-\tanh(\beta b_{m})]x_{n}. We conclude that
[TABLE]
with , as before. The average is taken over the output of the network and over the reward distribution, just as in Equation (11.5).
Now we seek a learning rule that increases the expected immediate reward . In other words, we require that the weight increment is an unbiased of the gradient of the expected immediate reward [179]:
[TABLE]
Comparison with Equation (11.7) leads to:
[TABLE]
with . This learning rule belongs to a set of more general rules derived by Williams [179]. It is plausible that the rule (11.9) converges to a steady state, because the weight increments approach zero as the network learns to produce the output \max_{\small\mbox{\boldmathy}}\{p_{\rm reward}(\mbox{\boldmathx},\mbox{\boldmathy})\}, independent of , so that \mbox{\boldmathy}-\langle\mbox{\boldmathy}\rangle averages to zero. But there is no proof of convergence.
An alternative is the associative reward-penalty rule [185]:
[TABLE]
with . For , the learning rules (11.9) and (11.10) give the same weight increment, but for the increments are different. With rule (11.10) the agent learns primarily from positive feedback. One advantage of this asymmetric rule is that it can be proven to converge in the limit of [185]. In general, however, the convergence becomes quite slow when is small. Figure 11.2(c) shows simulation results for the immediate reward, averaged over independent realisations of the learning process, versus the iteration number of the rule (11.10). We see that the average immediate reward approaches a steady state. The steady-state average of the immediate reward is smaller than , but as expected it approaches as decreases.
The averaged learning curves still exhibit substantial fluctuations. They reflect significant variations within and between individual realisations. Furthermore, the convergence proof assumes that the input patterns are linearly independent. This means that the number of patterns cannot exceed the input dimension . Associative reinforcement problems with linearly dependent inputs can be solved by embedding the input patterns in a higher-dimensional input space (Section 5.4).
The associative reward-penalty rule illustrates how an agent can use a reinforcement signal to maximise the expected immediate reward. The algorithm could for instance be a model for how an animal learns to respond in different ways to different stimuli.
Yet there are many problems where the reward is not immediate. When we play chess, the reward comes at the end of the game, for example if we won, if we lost, and if the game ended in a draw. More generally, an agent navigating a complex environment should not only consider immediate rewards, but also how a certain action may affect possible future rewards. One way of estimating future rewards for such tasks is temporal difference learning, discussed next.
11.2 Temporal difference learning
How does temporal difference learning allow an agent to optimise its expected future reward? For an episodic task, given an episode with steps, the agent visits the finite sequence of states \mbox{\boldmaths}_{0},\ldots,\mbox{\boldmaths}_{T-1}, and collects the rewards . The future reward is defined as
[TABLE]
Continuous tasks, by contrast, do not have defined end points. Since the sum in (11.11) might diverge as , it is customary to introduce a weighting factor in the sum over rewards:
[TABLE]
The weighting factor reduces the contribution of the far future to the estimate. Smaller values of give more weight to the immediate future, and the limit corresponds to . The sum in Equation (11.12) is called future discounted reward.
We use a neural network with input \mbox{\boldmaths}_{t} to estimate . In general, the network output is a non-linear function of the inputs, parameterised by weights that could be arranged into several layers of hidden neurons (Part II). The simplest choice is to use a single linear unit, just as in Equation (5.19):
[TABLE]
The components of the weight vector are determined so that the network output O(\mbox{\boldmaths}_{t}) approximates . This can be achieved by minimising the energy function
[TABLE]
using gradient descent. The corresponding learning rule reads:
[TABLE]
The idea of temporal difference learning [182] is to express the error R_{t}-O(\mbox{\boldmaths}_{t}) as a sum of temporal differences:
[TABLE]
where O(\mbox{\boldmaths}_{T}) is defined to be zero, O(\mbox{\boldmaths}_{T})\equiv 0. Using the gradient-descent rule (11.15) one obtains
[TABLE]
The terms in this double sum can be summed in a different way, as illustrated in Figure 11.3:
[TABLE]
Exchanging the summation variables and introducing a weighting factor gives:
[TABLE]
The purpose of the weighting factor is to reduce the weight of past states in the sum [186]. Alternatively one may update (and hence ) at each time step, with increment [182]
[TABLE]
This is the temporal-difference learning rule, also called TD()
[186]. Temporal difference learning allows a machine to learn the board game backgammon [17, 187], using a deep layered network and backpropagation (Section 6.1) to determine the weights.
The rule TD([math]) is similar to the learning rule (6.6a) with target r_{t+1}+O({\mbox{\boldmathw}_{t},}\mbox{\boldmaths}_{t+1}). It allows to learn one-step prediction of the time series. Using Equation (11.13), we see that the TD([math])-learning rule corresponds to the following learning rule for the output :
[TABLE]
The subscript in emphasises that the output function is updated iteratively . The learning rule (11.21) applies to estimating the future reward (11.11) for episodic tasks. If the environment is stationary, one may average over many consecutive episodes, using the final weights from episode as initial weight values for episode . For continuous tasks, the corresponding rule for estimating the future discounted reward (11.12) reads:
[TABLE]
Returning to the problem outlined in the beginning of this Chapter, consider an agent exploring a complex environment. The task might be to get from location to location as quickly as possible, or expending as little energy as possible. At time the agent is at position \mbox{\boldmathx}_{t} with velocity \mbox{\boldmathv}_{t}. These variables as well as the local state of the environment are summarised in the state vector \mbox{\boldmaths}_{t}. Given \mbox{\boldmaths}_{t}, the agent can act in certain ways: it might for example slow down, speed up, or turn. These possible actions are summarised in a vector \mbox{\boldmatha}_{t}. At each time step, the agent takes the action \mbox{\boldmatha}_{t} that optimises the expected future discounted reward (11.12), given its present state \mbox{\boldmaths}_{t}. The estimated expected future reward for any state-action pair is summarised in a table: the -table with elements Q_{t}(\mbox{\boldmaths}_{t},\mbox{\boldmatha}_{t}) is the analogue of O_{t}(\mbox{\boldmaths}_{t}). Different rows of the -table correspond to different states, and different columns to different actions. The TD([math]) rule for the table reads:
[TABLE]
This algorithm is called SARSA, because one needs \mbox{\boldmaths}_{t},\mbox{\boldmatha}_{t},r_{t+1},\mbox{\boldmaths}_{t+1}, and \mbox{\boldmatha}_{t+1} to update the -table (Figure 11.4). A difficulty with the rule (11.23) is that it depends not only on the present state-action pair [\mbox{\boldmaths}_{t},\mbox{\boldmatha}_{t}], but also on the next action \mbox{\boldmatha}_{t+1}, and thus indirectly upon the policy. Sometimes this is indicated by writing for the -table given policy .
11.3 -learning
The -learning rule [188] is an approximation to Eq. (11.23) that does not depend on \mbox{\boldmatha}_{t+1}. Instead one assumes that the next action, \mbox{\boldmatha}_{t+1}, is the optimal one:
[TABLE]
regardless of the policy that is currently followed. Although Equation (11.24) does not refer to any policy, the learning outcome nevertheless depends on it, because the policy determines the sequence of states and actions [\mbox{\boldmaths}_{t},\mbox{\boldmatha}_{t}]. For the greedy policy, Eq. (11.24) is equivalent to (11.23), but in general the two algorithms differ, and converge to different solutions. While the -table converges to the expected future reward of the greedy policy in -learning, for SARSA it converges to the expected future reward corresponding to the policy used in training. This can be an advantage when performance during training is important, for example when training an expensive robot that should not crash too often, or for a small bird that learns flying by doing. -learning is simpler, and can be used for problems where the final strategy based on \mbox{argmax}_{\mbox{\boldmatha}^{\prime}}Q(\mbox{\boldmaths},\mbox{\boldmatha}^{\prime}) matters, but where the reward during training is less important. Examples are board games where only the quality of the final strategy counts, not how often one loses during training. In summary, -learning is simpler, and if one takes during training, it yields the optimal strategy as well as SARSA, but it may give lower rewards during training.
The -learning algorithm is summarised in Algorithm 11. Usually one sets the initial entries in the -table to large positive values (optimistic initialisation), because this prompts the agent to explore many different actions, at least in the beginning. If the agent is in state \mbox{\boldmaths}_{t}, it chooses the action \mbox{\boldmatha}_{t} from Q_{t}(\mbox{\boldmaths}_{t},\mbox{\boldmatha}_{t}) according to the given policy. For the -greedy policy, for example, the agent picks a random action from the corresponding row of the -table with probability . With probability it chooses the action \mbox{\boldmatha}_{t} that yields the largest111If several elements in the relevant row have the same maximal value then any one of them is chosen with equal probability. Q_{t}(\mbox{\boldmaths}_{t},\mbox{\boldmatha}_{t}) given \mbox{\boldmaths}_{t}. The choice of action \mbox{\boldmatha}_{t} determines the next state \mbox{\boldmaths}_{t+1}, and this in turn allows to update the -table: given the new state \mbox{\boldmaths}_{t+1} resulting from the action \mbox{\boldmatha}_{t}, one updates Q_{t}(\mbox{\boldmaths}_{t},\mbox{\boldmatha}_{t}) using Equation (11.24). For episodic tasks one puts in (11.24), and one averages over many episodes using the outcome from episode as initial condition for the -table for episode . Each new episode can start with a new initial state \mbox{\boldmaths}_{0}. It helps the exploration process if \mbox{\boldmaths}_{0} is one of the states that are rarely visited by the learning algorithm.
When the sequence \mbox{\boldmaths}_{0},\mbox{\boldmaths}_{1},\mbox{\boldmaths}_{2},\ldots is a Markov chain (Section 4.2), then the -learning algorithm can be shown [189, 186] to converge if one uses a time-dependent learning rate that satisfies
[TABLE]
Often -learning is implemented in combination with the -greedy policy. This policy shares an important property with the associative reward-penalty algorithm with stochastic neurons: stochasticity allows for a wider range of responses, some of which may turn out beneficial in the long run. When is very small, the agent picks the action that appears optimal. As a consequence, suboptimal -elements are sampled less frequently and are therefore subject to larger errors. Therefore it is advantageous to begin with a relatively large value of . It is customary to decrease as the algorithm is iterated, because this accelerates convergence to the greedy policy.
It is important to bear in mind that the learning outcome depends on the reward function, as mentioned above. In general it is a good idea to analyse how the optimal strategy changes as one varies the reward function. Sometimes we are faced with the inverse problem: consider how a microorganism swimming in the turbulent ocean responds to different stimuli. How was this behaviour shaped by genetic evolution? Which quantity was optimised? Is it most important to reduce the energy cost for propulsion? Or is it more important to avoid predation?
Another challenge is to determine suitable states and actions. An agent navigating a complex environment may have a continuous range of positions and velocities, and may experience continuous-valued signals from the environment. To represent the corresponding states in a -table it is necessary to discretise. To this end one must determine suitable ranges and resolutions of these variables, and for the actions. If there are too many states and actions, -learning becomes inefficient. This is referred to as the curse of dimensionality [190, 16].
Let us see how -learning works for a very simple example, for the associative task described in Fig. 11.2. One episode corresponds to computing the output of the neuron given its initial state, so . There is no sequence of states, and the task is to estimate the immediate reward. In this case the learning rule (11.24) simplifies to
[TABLE]
Since each episode consists only of a single time step, we dropped the subscript . Also, the term \max_{\footnotesize\mbox{\boldmatha}}Q(\mbox{\boldmaths}_{t+1},\mbox{\boldmatha}) from Equation (11.24) does not appear in (11.26) since estimates the immediate reward. There are only two states in this problem, \mbox{\boldmaths}={\mbox{\boldmathx}^{(1)}} and \mbox{\boldmaths}={\mbox{\boldmathx}^{(2)}}, and the possible actions are . In other words, in this case. In each round, one of the states is chosen randomly, with equal probability. The action is determined from the current estimate of the immediate reward as {\rm argmax}_{a}Q(\mbox{\boldmaths},a) with probability , and uniformly randomly otherwise. These steps are iterated over many iterations (episodes), using the outcome of episode as initial condition for episode . The rule (11.26) describes exponential relaxation to the target for small learning rate . In this limit, Equation (11.26) is approximated by the stochastic differential equation
[TABLE]
where f_{\varepsilon}(\mbox{\boldmaths},a) is the stationary frequency with which the state-action pair [\mbox{\boldmaths},a] is visited using the -greedy policy:
[TABLE]
The frequency is normalised to unity, 1=\sum_{s,a}f_{\varepsilon}(\mbox{\boldmaths},a). Averaging the solution of Equation (11.27) with initial condition Q_{0}(\mbox{\boldmaths},a)=1 over the reward distribution gives:
[TABLE]
For , Q_{t}(\mbox{\boldmaths},\mbox{\boldmatha}) converges on average to
[TABLE]
Here we used that \langle r(\mbox{\boldmathx},y)\rangle=2p_{\rm reward}(\mbox{\boldmathx},y)-1, as well the reward probabilities in Figure 11.2. Figure 11.5 illustrates how the rate of convergence depends on the value of the parameter . For , all state-action pairs are visited and evaluated equally often, independently of the present elements of the -table. Equations (11.28) and (11.29) show that all elements of the -table converge to their steady-state values at the same rate, equal to . For small values of , by contrast, the algorithm tends to take optimal actions, \mbox{argmax}_{a^{\prime}}Q_{k}(\mbox{\boldmaths},a^{\prime}). Therefore it finds the optimal elements of the -table more quickly, at the rate . However, the other elements converge much more slowly. Initially, the theory (11.29) does not apply because optimal and suboptimal -elements in each row of the -table are not well separated. As a result the decay rates of all elements are similar at first. But once optimal and suboptimal elements are significantly different, the suboptimal ones decay at the rate , as predicted by the theory.
This example illustrates the strength of -learning with the -greedy policy. For small values of , the algorithm tends to converge to the optimal strategy more quickly than a brute-force algorithm that visits every state-action pair equally often. Equation (11.28) shows that this advantage becomes larger for larger values of . The price one pays is that the suboptimal entries of the -table converge more slowly.
The gain is even more significant for episodic tasks with steps. In this case the learning rule for the table depends on \mbox{argmax}_{\footnotesize\mbox{\boldmatha}^{\prime}}Q_{t}(\mbox{\boldmaths},\mbox{\boldmatha}^{\prime}). As mentioned above, there are ways in which the largest elements can be distributed over the rows of the -table. To find the optimal deterministic strategy by complete enumeration, one must run a large number of episodes for each of the possibilities, to find out which one gives the largest reward. This is impractical when either or or both become too large. -learning with a small value of , by contrast, tends to explore actions that appear to yield the largest expected future reward given the current estimates of the -values. Using experience in this way allows the algorithm to simultaneously improve its policy and the -values towards optimality. As a consequence, -learning can find optimal, or at least good strategies when complete enumeration of all possibilities fails.
A second example for -learning is illustrated in Figure 11.6, the board game tic-tac-toe. It is a very simple game where two players take turns in placing their pieces on a board. The player who manages to first obtain three pieces in a row, column, or diagonal wins and receives the reward . A draw gives , and the player receives when the round is lost. The goal is to win as often as possible, to maximise the expected future reward. However, there is a strategy for both players to ensure that they do not lose. If both players try to maximise their expected future reward, then they end up following this strategy. As a consequence, every game must end in a draw [191]. As a result, the game is quite boring.
Nevertheless it is instructive to ask how the players can learn to find this strategy using -learning with the -greedy policy.
To this end we let two agents play many rounds against each other. The state space is the collection of all board configurations. Player starts, and thus always sees a board with an even number of pieces, while the number of pieces is odd for player . Since the players encounter different sets of states, each must keep track of their own -table. The task is episodic, and the number of steps may vary from round to round. Feedback is only obtained at the end of each round.
We use Equation (11.24) with a constant learning rate . We can set since the number of steps in each round is finite.
The -table is a table where each entry is a array. Here is the number of states the player encountered so far. The first row lists the states, each a array with entries (), (), or [math] (empty). The second row contains the -values. Given a certain state of the board, a player can play a piece on any empty field. The corresponding estimate of the expected future reward is stored in the corresponding array in the second row. Since one cannot place a piece onto an occupied field, the corresponding entries in the -table are assigned NaN.
During a round, the -tables of the players are updated in turns, always the one of the player who places a piece. After a round, the -tables of both players are updated. During the first round, the elements of -tables encountered are initialised to zero and remain zero. The first change to occurs in the last step of this round. If player wins, for example, the element Q_{T-1}(\mbox{\boldmaths}_{T-1},a_{T-1}) corresponding to the state-action pair that led to the final state \mbox{\boldmaths}_{T-1}^{\prime} is updated for the winning player, and Q_{T-2}(\mbox{\boldmaths}_{T-2}^{\prime},\mbox{\boldmatha}_{T-2}^{\prime}) is set to for player .
Both players follow the -greedy policy . With probability they take the optimal move (if the maximal -element in the relevant row is degenerate, then one of the maximal elements is chosen randomly). With probability , a random action is chosen. As the players continue to play rounds against each other, the rewards spread to other elements of . Suppose that the state \mbox{\boldmaths}_{T-1} is encountered once more, the one that allowed player to win the first round with \mbox{\boldmatha}_{T-1}. Then the term {\rm max}_{\scriptstyle a}Q_{T-1}(\mbox{\boldmaths}_{T-1},\mbox{\boldmatha}) causes a -element for the previous state to change, the one from which \mbox{\boldmaths}_{T-1} was reached the second time. However as times goes on, this process slows down, because later updates are multiplied with higher powers of the learning rate . Also, if the opponent lost in the previous round, it will try different actions that may block winning moves for the other player.
Figure 11.7 illustrates how the players learn, after playing many rounds against each other. Since both players try to maximise their expected future reward in the steady state of the -learning algorithm, all games end in a draw in this case. The corresponding -tables contain the strategies each player should adopt to maximise their reward. Suppose player places the first piece as shown below:
[TABLE]
How does this game continue? There are several different ways in which player may try to win. The left -table in Equation (11.31) shows that one possibility is to place the piece in the top or bottom row, because this creates the opportunity of creating a bridge in the next move, a configuration that cannot be blocked by the opponent, allowing to win. The right -table shows that player could still lose or end up with a draw if he makes the wrong move. The corresponding -entries have not quite converged to and [math], respectively. -entries corresponding to suboptimal states are not estimated as precisely because they are visited less frequently. Here was chosen quite large. Smaller values of give even less accurate estimates for the suboptimal -elements compared with those in Equation (11.31), after training for the same number of rounds.
As pointed out above, the learning outcome depends on the reward function. If one increases reward for winning, to for instance, the optimal strategy appears to be to take turns in winning. The same learning outcome is expected if one imposes a penalty for a draw, (win), (draw, lose), Exercise 11.8. More examples of reinforcement-learning problems in robotics and in the natural sciences are described in Ref. [192].
The -learning algorithm described above is quite efficient when the number of states and actions is not too large. For very large -tables, the algorithm becomes quite slow. In this case it may be more efficient to replace the -table by an approximate -function that maps states to actions. As explained in Section 11.2, one can use a neural network to represent the -function [184] (deep reinforcement learning). An application of this method is AlphaGo, a machine-learning algorithm that learnt to play the game of go [183]. The -function is represented in terms of a convolutional neural network. This makes it possible to use a variant of -learning despite the fact that the number of states is enormous.
In recent years many proof-of-principle studies have demonstrated the possibilities of reinforcement learning in a wide range of scientific problems. Recent advances in deep reinforcement learning hold promise for the future, for real-world control problems in the engineering sciences.
11.4 Summary
Reinforcement learning lies between unsupervised learning (Chapter 10.1) and supervised learning (Chapters 5 to 9). In reinforcement learning, there are no labeled data sets. Instead, the neural network or agent learns through feedback from the environment in the shape of a reward or a penalty. The goal is to find a strategy that maximises the expected reward. Reinforcement learning is applied in a wide range of fields, from psychology to mechanical engineering, using a large variety of algorithms. The associative reward-penalty algorithm and many versions of temporal difference learning were originally formulated using neural networks. -learning is an approximation to temporal difference learning for sequential decision processes. In its simplest form it does not rely on neural networks. However, when the number of states and actions is large this algorithm becomes slow. In this case it may be more efficient to approximate the -function by a neural network.
11.5 Further reading
The standard reference for reinforcement learning is Reinforcement learning: an introduction by Sutton and Barto [16]. The original reference for the convergence of the -learning algorithm is Ref. [189]. A more mathematical introduction to reinforcement learning is given in Ref. [186]. Examples for reinforcement learning in statistical and non-linear physics are summarised in Ref. [192].
An open question is when and how symmetries can be exploited to simplify a reinforcement problem. For a small microorganism learning to navigate a turbulent flow, some aspects are discussed in Ref. [193], but little is known in general. Another open question concerns the convergence of the -learning algorithm. Convergence to the optimal policy is assured if the sequence of states is a Markov chain. However, most real-world problems are not Markovian, so that convergence is not guaranteed. The algorithm appears to perform well nevertheless (a recurring theme in this book), but it is an open question under which circumstances it may fail.
\nobibliography
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Hertz et al. [1991] Hertz, J, Krogh, A & Palmer, R 1991 Introduction to the Theory of Neural Computation . Addison-Wesley.
- 2Haykin [1999] Haykin, S 1999 Neural Networks: a comprehensive foundation , 2nd edn. New Jersey: Prentice Hall.
- 3Horner [2001] Horner, H , Neuronale Netze, www.tphys.uni-heidelberg.de/~horner , [Last accessed 8-November-2018].
- 4Goodfellow et al. [2016] Goodfellow, I. J, Bengio, Y & Courville, A , Deep learning, www.deeplearningbook.org , [Last accessed 5-September-2018].
- 5Nielsen [2017] Nielsen, M , Neural networks and deep learning, http://neuralnetworksanddeeplearning.com , [Last accessed 13-August-2018].
- 6Mc Culloch & Pitts [1943] Mc Culloch, W & Pitts, W 1943 A logical calculus of the ideas immanent in nervous activity. Bull. Math. Biophys. 5 , 115.
- 7Carnap [1937] Carnap, R 1937 The logical syntax of language . London: K. Paul, Trench, Trubner & Co., Limited.
- 8Mc Culloch & Pitts [1947] Mc Culloch, W & Pitts, W 1947 How we know universals the perception of auditory and visual forms. Bull. Math. Biophys. 9 , 127–147.
