1803 08823 PDF
1803 08823 PDF
Marin Bukov
Department of Physics,
University of California,
arXiv:1803.08823v2 [[Link]-ph] 1 Feb 2019
Berkeley, CA 94720,
USA†
Charles K. Fisher
[Link], San Francisco,
CA 94108
David J. Schwab
Initiative for the Theoretical Sciences,
The Graduate Center,
City University of New York,
365 Fifth Ave., New York,
NY 10016
Machine Learning (ML) is one of the most exciting and dynamic areas of modern re-
search and application. The purpose of this review is to provide an introduction to the
core concepts and tools of machine learning in a manner easily understood and intuitive
to physicists. The review begins by covering fundamental concepts in ML and modern
statistics such as the bias-variance tradeoff, overfitting, regularization, generalization,
and gradient descent before moving on to more advanced topics in both supervised
and unsupervised learning. Topics covered in the review include ensemble models, deep
learning and neural networks, clustering and data visualization, energy-based models (in-
cluding MaxEnt models and Restricted Boltzmann Machines), and variational methods.
Throughout, we emphasize the many natural connections between ML and statistical
physics. A notable aspect of the review is the use of Python Jupyter notebooks to in-
troduce modern ML/statistical packages to readers using physics-inspired datasets (the
Ising Model and Monte-Carlo simulations of supersymmetric decays of proton-proton
collisions). We conclude with an extended outlook discussing possible uses of machine
learning for furthering our understanding of the physical world as well as open problems
in ML where physicists may be able to contribute.
XVIII. Outlook 109 sity in Fall of 2016. As such, it assumes a level of familiar-
A. Research at the intersection of physics and ML 109 ity with several topics found in graduate physics curricula
B. Topics not covered in review 110
(partition functions, statistical mechanics) and a fluency
C. Rebranding Machine Learning as “Artificial
Intelligence” 111 in mathematical techniques such as linear algebra, multi-
D. Social Implications of Machine Learning 111 variate calculus, variational methods, probability theory,
and Monte-Carlo methods. It also assumes a familiar-
XIX. Acknowledgments 112
ity with basic computer programming and algorithmic
A. Overview of the Datasets used in the Review 112 design.
1. Ising dataset 112
2. SUSY dataset 112
3. MNIST Dataset 113 A. What is Machine Learning?
References 113
Most physicists learn the basics of classical statistics
early on in undergraduate laboratory courses. Classical
I. INTRODUCTION statistics is primarily concerned with how to use data
to estimate the value of an unknown quantity. For in-
Machine Learning (ML), data science, and statistics stance, estimating the speed of light using measurements
are fields that describe how to learn from, and make pre- obtained with an interferometer is one such example that
dictions about, data. The availability of big datasets is relies heavily on techniques from statistics.
a hallmark of modern science, including physics, where Machine Learning is a subfield of artificial intelligence
data analysis has become an important component of di- with the goal of developing algorithms capable of learning
verse areas, such as experimental particle physics, ob- from data automatically. In particular, an artificially in-
servational astronomy and cosmology, condensed matter telligent agent needs to be able to recognize objects in its
physics, biophysics, and quantum computing. Moreover, surroundings and predict the behavior of its environment
ML and data science are playing increasingly important in order to make informed choices. Therefore, techniques
roles in many aspects of modern technology, ranging from in ML tend to be more focused on prediction rather than
biotechnology to the engineering of self-driving cars and estimation. For example, how do we use data from the
smart devices. Therefore, having a thorough grasp of the interferometry experiment to predict what interference
concepts and tools used in ML is an important skill that pattern would be observed under a different experimental
is increasingly relevant in the physical sciences. setup? In addition, methods from ML tend to be applied
The purpose of this review is to serve as an introduc- to more complex high-dimensional problems than those
tion to foundational and state-of-the-art techniques in typically encountered in a classical statistics course.
ML and data science for physicists. The review seeks to Despite these differences, estimation and prediction
find a middle ground between a short overview and a full- problems can be cast into a common conceptual frame-
length textbook. While there exist many wonderful ML work. In both cases, we choose some observable quantity
textbooks (Abu-Mostafa et al., 2012; Bishop, 2006; Fried- x of the system we are studying (e.g., an interference pat-
man et al., 2001; Murphy, 2012), they are lengthy and use tern) that is related to some parameters θ (e.g., the speed
specialized language that is often unfamiliar to physicists. of light) of a model p(x|θ) that describes the probability
This review builds upon the considerable knowledge most of observing x given θ. Now, we perform an experiment
physicists already possess in statistical physics in order to obtain a dataset X and use these data to fit the model.
to introduce many of the major ideas and techniques Typically, “fitting” the model involves finding θ̂ that pro-
used in modern ML. We take a physics-inspired peda- vides the best explanation for the data. In the case when
gogical approach, emphasizing simple examples (e.g., re- “fitting” refers to the method of least squares, the esti-
gression and clustering), before delving into more ad- mated parameters maximize the probability of observ-
vanced topics. The intention of this review and the ing the data (i.e., θ̂ = argmaxθ {p(X|θ)}). Estimation
accompanying Jupyter notebooks (available at https: problems are concerned with the accuracy of θ̂, whereas
//[Link]/~pankajm/[Link]) is to prediction problems are concerned with the ability of the
give the reader the requisite background knowledge to model to predict new observations (i.e., the accuracy of
follow and apply these techniques to their own areas of p(x|θ̂)). Although the goals of estimation and prediction
interest. are related, they often lead to different approaches. As
While this review is written with a physics background this review is aimed as an introduction to the concepts of
in mind, we aim for it to be useful to anyone with some ML, we will focus on prediction problems and refer the
background in statistical physics, and it is suitable for reader to one of many excellent textbooks on classical
both graduate students and researchers as well as ad- statistics for more information on estimation (Lehmann
vanced undergraduates. The review is based on an ad- and Casella, 2006; Lehmann and Romano, 2006; Wasser-
vanced topics graduate course taught at Boston Univer- man, 2013; Witte and Witte, 2013).
4
The last three decades have seen an unprecedented in- Any review on ML must simultaneously accomplish
crease in our ability to generate and analyze large data two related but distinct goals. First, it must convey the
sets. This “big data” revolution has been spurred by an rich theoretical foundations underlying modern ML. This
exponential increase in computing power and memory task is made especially difficult because ML is very broad
commonly known as Moore’s law. Computations that and interdisciplinary, drawing on ideas and intuitions
were unthinkable a few decades ago can now be routinely from many fields including statistics, computational neu-
performed on laptops. Specialized computing machines roscience, and physics. Unfortunately, this means mak-
(such as GPU-based machines) are continuing this trend ing choices about what theoretical ideas to include in the
towards cheap, large-scale computation, suggesting that review. This review emphasizes connections with sta-
the “big data” revolution is here to stay. tistical physics, physics-inspired Bayesian inference, and
This increase in our computational ability has been ac- computational neuroscience models. Thus, certain ideas
companied by new techniques for analyzing and learning (e.g., gradient descent, expectation maximization, varia-
from large datasets. These techniques draw heavily from tional methods, and deep learning and neural networks)
ideas in statistics, computational neuroscience, computer are covered extensively, while other important ideas are
science, and physics. Similar to physics, modern ML given less attention or even omitted entirely (e.g., statis-
places a premium on empirical results and intuition over tical learning, support vector machines, kernel methods,
the more formal treatments common in statistics, com- Gaussian processes). Second, any ML review must give
puter science, and mathematics. This is not to say that the reader the practical know-how to start using the tools
proofs are not important or undesirable. Rather, many and concepts of ML for practical problems. To accom-
of the advances of the last two decades – especially in plish this, we have written a series of Jupyter notebooks
fields like deep learning – do not have formal justifica- to accompany this review. These python notebooks in-
tions (much like there still exists no mathematically well- troduce the nuts-and-bolts of how to use, code, and im-
defined concept of the Feynman path-integral in d > 1). plement the methods introduced in the main text. Luck-
ily, there are numerous great ML software packages avail-
Physicists are uniquely situated to benefit from and able in Python (scikit-learn, tensorflow, Pytorch, Keras)
contribute to ML. Many of the core concepts and tech- and we have made extensive use of them. We have also
niques used in ML – such as Monte-Carlo methods, simu- made use of a new package, Paysage, for energy-based
lated annealing, variational methods – have their origins generative models which has been co-developed by one
in physics. Moreover, “energy-based models” inspired by of the authors (CKF) and maintained by [Link] (a
statistical physics are the backbone of many deep learn- company affiliated with two of the authors: CKF and
ing methods. For these reasons, there is much in modern PM). The purpose of the notebooks is to both familiarize
ML that will be familiar to physicists. physicists with these resources and to serve as a starting
Physicists and astronomers have also been at the fore- point for experimenting and playing with ideas.
front of using “big data”. For example, experiments such ML can be divided into three broad categories: super-
as CMS and ATLAS at the LHC generate petabytes of vised learning, unsupervised learning, and reinforcement
data per year. In astronomy, projects such as the Sloan learning. Supervised learning concerns learning from la-
Digital Sky Survey (SDSS) routinely analyze and release beled data (for example, a collection of pictures labeled
hundreds of terabytes of data measuring the properties of as containing a cat or not containing a cat). Common
nearly a billion stars and galaxies. Researchers in these supervised learning tasks include classification and re-
fields are increasingly incorporating recent advances in gression. Unsupervised learning is concerned with find-
ML and data science, and this trend is likely to acceler- ing patterns and structure in unlabeled data. Examples
ate in the future. of unsupervised learning include clustering, dimensional-
Besides applications to physics, part of the goal of this ity reduction, and generative modeling. Finally, in rein-
review is to serve as an introductory resource for those forcement learning an agent learns by interacting with an
looking to transition to more industry-oriented projects. environment and changing its behavior to maximize its
Physicists have already made many important contribu- reward. For example, a robot can be trained to navigate
tions to modern big data applications in an industrial in a complex environment by assigning a high reward to
setting (Metz, 2017). Data scientists and ML engineers actions that help the robot reach a desired destination.
in industry use concepts and tools developed for ML to We refer the interested reader to the classic book by Sut-
gain insight from large datasets. A familiarity with ML ton and Barto Reinforcement Learning: an Introduction
is a prerequisite for many of the most exciting employ- (Sutton and Barto, 1998). While useful, the distinction
ment opportunities in the field, and we hope this review between the three types of ML is sometimes fuzzy and
will serve as a useful introduction to ML for physicists fluid, and many applications often combine them in novel
beyond an academic setting. and interesting ways. For example, the recent success
5
of Google DeepMind in developing ML algorithms that MNIST dataset, on the other hand, introduces the fla-
excel at tasks such as playing Go and video games em- vor of present-day ML problems. By re-analyzing the
ploy deep reinforcement learning, combining reinforce- same datasets with multiple techniques, we hope readers
ment learning with supervised learning methods based will be able to get a sense of the various, inevitable trade-
on deep neural networks. offs involved in choosing how to analyze data. Certain
Here, we limit our focus to supervised and unsuper- techniques work better when data is limited while others
vised learning. The literature on reinforcement learning may be better suited to large data sets with many fea-
is extensive and uses ideas and concepts that, to a large tures. A short description of these datasets are given in
degree, are distinct from supervised and unsupervised the Appendix.
learning tasks. For this reason, to ensure cohesiveness This review draws generously on many wonderful text-
and limit the length of this review, we have chosen not books on ML and we encourage the reader to con-
to discuss reinforcement learning. However, this omis- sult them for further information. They include Abu
sion should not be mistaken for a value judgement on Mostafa’s masterful Learning from Data, which intro-
the utility of reinforcement learning for solving physical duces the basic concepts of statistical learning theory
problems. For example, some of the authors have used (Abu-Mostafa et al., 2012), the more advanced but
inspiration from reinforcement learning to tackle difficult equally good The Elements of Statistical Learning by
problems in quantum control (Bukov, 2018; Bukov et al., Hastie, Tibshirani, and Friedman (Friedman et al., 2001),
2018). Michael Nielsen’s indispensable Neural Networks and
In writing this review, we have tried to adopt a style Deep Learning which serves as a wonderful introduction
that reflects what we consider to be the best of the to the neural networks and deep learning (Nielsen, 2015)
physics tradition. Physicists understand the importance and David MacKay’s outstanding Information Theory,
of well-chosen examples for furthering our understand- Inference, and Learning Algorithms which introduced
ing. It is hard to imagine a graduate course in statistical Bayesian inference and information theory to a whole
physics without the Ising model. Each new concept that generation of physicists (MacKay, 2003). More com-
is introduced in statistical physics (mean-field theory, prehensive (and much longer) books on modern ML
transfer matrix techniques, high- and low-temperature techniques include Christopher Bishop’s classic Pattern
expansions, the renormalization group, etc.) is applied to Recognition and Machine Learning (Bishop, 2006) and
the Ising model. This allows for the progressive building the more recently published Machine Learning: A Prob-
of intuition and ultimately a coherent picture of statisti- abilistic Perspective by Kevin Murphy (Murphy, 2012).
cal physics. We have tried to replicate this pedagogical Finally, one of the great successes of modern ML is deep
approach in this review by focusing on a few well-chosen learning, and some of the pioneers of this field have writ-
techniques – linear and logistic regression in the case of ten a textbook for students and researchers entitled Deep
supervised learning and clustering in the case of unsu- Learning (Goodfellow et al., 2016). In addition to these
pervised learning – to introduce the major theoretical textbooks, we have consulted numerous research papers,
concepts. reviews, and web resources. Whenever possible, we have
In this same spirit, we have chosen three interest- tried to point the reader to key papers and other refer-
ing datasets with which to illustrate the various algo- ences that we have found useful in preparing this review.
rithms discussed here. (i) The SUSY data set consists However, we are neither capable of nor have we made any
of 5, 000, 000 Monte-Carlo samples of proton-proton col- effort to make a comprehensive review of the literature.
lisions decaying to either signal or background processes, The review is organized as follows. We begin by
which are both parametrized with 18 features. The sig- introducing polynomial regression as a simple example
nal process is the production of electrically-charged su- that highlights many of the core ideas of ML. The next
persymmetric particles, which decay to W bosons and an few chapters introduce the language and major concepts
electrically-neutral supersymmetric particle, invisible to needed to make these ideas more precise including tools
the detector, while the background processes are various from statistical learning theory such as overfitting, the
decays involving only Standard Model particles (Baldi bias-variance tradeoff, regularization, and the basics of
et al., 2014). (ii) The Ising data set consists of 104 Bayesian inference. The next chapter builds on these
states of the 2D Ising model on a 40 × 40 square lat- examples to discuss stochastic gradient descent and its
tice, obtained using Monte-Carlo (MC) sampling at a generalizations. We then apply these concepts to linear
few fixed temperatures T . (iii) The MNIST dataset com- and logistic regression, followed by a detour to discuss
prises 70000 handwritten digits, each of which comes in how we can combine multiple statistical techniques to
a square image, divided into a 28 × 28 pixel grid. The improve supervised learning, introducing bagging, boost-
first two datasets were chosen to reflect the various sub- ing, random forests, and XG Boost. These ideas, though
disciplines of physics (high-energy experiment, condensed fairly technical, lie at the root of many of the advances
matter) where we foresee techniques from ML becom- in ML over the last decade. The review continues with
ing an increasingly important tool for research. The a thorough discussion of supervised deep learning and
6
neural networks, as well as convolutional nets. We then on the training set is called the in-sample error Ein =
turn our focus to unsupervised learning. We start with C(ytrain , f (Xtrain ; θ)) and the value of the cost func-
data visualization and dimensionality reduction before tion on the test set is called the out-of-sample error
proceeding to a detailed treatment of clustering. Our Eout = C(ytest , f (Xtest ; θ)).
discussion of clustering naturally leads to an examina- One of the most important observations we can make
tion of variational methods and their close relationship is that the out-of-sample error is almost always greater
with mean-field theory. The review continues with a than the in-sample error, i.e. Eout ≥ Ein . We explore
discussion of deep unsupervised learning, focusing on this point further in Sec. VI and its accompanying note-
energy-based models, such as Restricted Boltzmann Ma- book. Splitting the data into mutually exclusive train-
chines (RBMs) and Deep Boltzmann Machines (DBMs). ing and test sets provides an unbiased estimate for the
Then we discuss two new and extremely popular model- predictive performance of the model – this is known as
ing frameworks for unsupervised learning, generative ad- cross-validation in the ML and statistics literature. In
versarial networks (GANs) and variational autoencoders many applications of classical statistics, we start with a
(VAEs). We conclude the review with an outlook and mathematical model that we assume to be true (e.g., we
discussion of promising research directions at the inter- may assume that Hooke’s law is true if we are observing
section physics and ML. a mass-spring system) and our goal is to estimate the
value of some unknown model parameters (e.g., we do
not know the value of the spring stiffness). Problems in
II. WHY IS MACHINE LEARNING DIFFICULT? ML, by contrast, typically involve inference about com-
plex systems where we do not know the exact form of the
A. Setting up a problem in ML and data science mathematical model that describes the system. There-
fore, it is not uncommon for ML researchers to have mul-
Many problems in ML and data science starts with tiple candidate models that need to be compared. This
the same ingredients. The first ingredient is the dataset comparison is usually done using Eout ; the model that
D = (X, y) where X is a matrix of independent variables minimizes this out-of-sample error is chosen as the best
and y is a vector of dependent variables. The second is model (i.e. model selection). Note that once we select
the model f (x; θ), which is a function f : x → y of the the best model on the basis of its performance on Eout ,
parameters θ. That is, f is a function used to predict an the real-world performance of the winning model should
output from a vector of input variables. The final ingre- be expected to be slightly worse because the test data
dient is the cost function C(y, f (X; θ)) that allows us to was now used in the fitting procedure.
judge how well the model performs on the observations y.
The model is fit by finding the value of θ that minimizes
the cost function. For example, one commonly used cost B. Polynomial Regression
function is the squared error. Minimizing the squared er-
ror cost function is known as the method of least squares, In the previous section, we mentioned that multiple
and is typically appropriate for experiments with Gaus- candidate models are typically compared using the out-
sian measurement errors. of-sample error Eout . It may be at first surprising that
ML researchers and data scientists follow a standard the model that has the lowest out-of-sample error Eout
recipe to obtain models that are useful for prediction usually does not have the lowest in-sample error Ein .
problems. We will see why this is necessary in the fol- Therefore, if our goal is to obtain a model that is use-
lowing sections, but it is useful to present the recipe up ful for prediction we may not want to choose the model
front to provide context. The first step in the analysis that provides the best explanation for the current obser-
is to randomly divide the dataset D into two mutually vations. At first glance, the observation that the model
exclusive groups Dtrain and Dtest called the training and providing the best explanation for the current dataset
test sets. The fact that this must be the first step should probably will not provide the best explanation for future
be heavily emphasized – performing some analysis (such datasets is very counter-intuitive.
as using the data to select important variables) before Moreover, the discrepancy between Ein and Eout be-
partitioning the data is a common pitfall that can lead to comes more and more important, as the complexity of our
incorrect conclusions. Typically, the majority of the data data, and the models we use to make predictions, grows.
are partitioned into the training set (e.g., 90%) with the As the number of parameters in the model increases,
remainder going into the test set. The model is fit by min- we are forced to work in high-dimensional spaces. The
imizing the cost function using only the data in the train- “curse of dimensionality” ensures that many phenomena
ing set θ̂ = arg minθ {C(ytrain , f (Xtrain ; θ))}. Finally, that are absent or rare in low-dimensional spaces become
the performance of the model is evaluated by computing generic. For example, the nature of distance changes in
the cost function using the test set C(ytest , f (Xtest ; θ̂)). high dimensions, as evidenced in the derivation of the
The value of the cost function for the best fit model Maxwell distribution in statistical physics where the fact
7
y
Training 1.0
−2 Linear
Poly 3 0.5
−4 Poly 10
0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.00 0.25 0.50 0.75 1.00 1.25
x x
4 Test
60 linear
2 3rd order
40 10th order
0
y
Training
−2 Linear 20
Poly 3
−4 Poly 10
0
0.0 0.2 0.4 0.6 0.8 1.0 0.00 0.25 0.50 0.75 1.00 1.25
x x
FIG. 1 Fitting versus predicting for noiseless data. Ntrain = 10 points in the range x ∈ [0, 1] were generated from a
linear model (top) or tenth-order polynomial (bottom). This data was fit using three model classes: linear models (red), all
polynomials of order 3 (yellow), all polynomials of order 10 (green) and used to make prediction on Ntest = 20 new data points
with xtest ∈ [0, 1.2] (shown on right). Notice that in the absence of noise (σ = 0), given enough data points that fitting and
predicting are identical.
that all the volume of a d-dimensional sphere of radius we have, the “noise” in the data generation process, and
r is contained in a small spherical shell around r is ex- our prior knowledge about the system. The goal is to
ploited. Almost all critical points of a function (i.e., the build intuition about why prediction is difficult in prepa-
points where all derivatives vanish) are saddles rather ration for introducing general strategies that overcome
than maxima or minima (an observation first made in these difficulties.
physics in the context of the p-spin spherical spin glass). Before reading the rest of the section, we strongly en-
For all these reasons, it turns out that for complicated courage the reader to read Notebook 1 and complete the
models studied in ML, predicting and fitting are very accompanying exercises.
different things (Bickel et al., 2006). Consider a probabilistic process that assigns a label yi
To develop some intuition about why we need to pay to an observation xi . The data are generated by drawing
close attention to out-of-sample performance, we will samples from the equation
consider a simple one-dimensional problem – polynomial yi = f (xi ) + ηi , (1)
regression. Our task is a simple one, fitting data with
polynomials of different order. We will explore how our where f (xi ) is some fixed (but possibly unknown) func-
ability to predict depends on the number of data points tion, and ηi is a Gaussian, uncorrelated noise variable,
8
y
Training
−2 Linear 0
Poly 3
−5
−4 Poly 10
−10
0.0 0.2 0.4 0.6 0.8 1.0 0.00 0.25 0.50 0.75 1.00 1.25
x x
Training
−2 Linear 0
Poly 3
−5
−4 Poly 10
−10
0.0 0.2 0.4 0.6 0.8 1.0 0.00 0.25 0.50 0.75 1.00 1.25
x x
FIG. 2 Fitting versus predicting for noisy data. Ntrain = 100 noisy data points (σ = 1) in the range x ∈ [0, 1] were
generated from a linear model (top) or tenth-order polynomial (bottom). This data was fit using three model classes: linear
models (red), all polynomials of order 3 (yellow), all polynomials of order 10 (green) and used to make prediction on Ntest = 20
new data points with xtest ∈ [0, 1.2](shown on right). Notice that even when the data was generated using a tenth order
polynomial, the linear and third order polynomials give better out-of-sample predictions, especially beyond the x range over
which the model was trained.
such that we will consider three different model classes: (i) all poly-
nomials of order 1 which we denote by f1 (x; θ1 ), (ii) all
hηi i = 0, polynomials up to order 3 which we denote by f3 (x; θ3 ),
hηi ηj i = δij σ 2 . and (iii) all polynomials of order 10, f10 (x; θ10 ). Notice
that these three model classes contain different number
We will refer to the f (xi ) as the function used to generate of parameters. Whereas f1 (x; θ1 ) has only two parame-
the data, and σ as the noise strength. The larger σ is the ters (the coefficients of the zeroth and first order terms
noisier the data; σ = 0 corresponds to the noiseless case. in the polynomial), f3 (x; θ3 ) and f10 (x; θ10 ) have four
To make predictions, we will consider a family of func- and eleven parameters, respectively. This reflects the
tions fα (x; θα ) that depend on some parameters θα . fact that these three models have different model com-
These functions represent the model class that we are us- plexities. If we think of each term in the polynomial as a
ing to model the data and make predictions. Note that “feature” in our model, then increasing the order of the
we choose the model class without knowing the function polynomial we fit increases the number of features. Using
f (x). The fα (x; θα ) encode the features we choose to a more complex model class may give us better predic-
represent the data. In the case of polynomial regression
9
tive power, but only if we have a large enough sample noise encodes real information. This problem is called
size to accurately learn the model parameters associated “overfitting” and leads to a steep drop-off in predictive
with these extra features from the training dataset. performance.
To learn the parameters θα , we will train our models We can guard against overfitting in two ways: we can
on a training dataset and then test the effectiveness of use less expressive models with fewer parameters, or we
the model on a different dataset, the test dataset. Since can collect more data so that the likelihood that the noise
we are interested only in gaining intuition, we will simply appears patterned decreases. Indeed, when we increase
plot the fitted polynomials and compare the predictions the size of the training data set by two orders of mag-
of our fits for the test data with the true values. As we nitude to Ntrain = 104 (see Figure 3) the tenth order
will see below, the models that give the best fit to existing polynomial clearly gives both the best fits and the most
data do not necessarily make the best predictions even predictive power over the entire training range x ∈ [0, 1],
for a simple task like polynomial regression. and even slightly beyond to approximately x ≈ 1.05.
To illustrate these ideas, we encourage the reader to This is our first experience with what is known as the
experiment with the accompanying notebook to gener- bias-variance tradeoff, c.f. Sec. III.B. When the amount
ate data using a linear function f (x) = 2x and a tenth of training data is limited as it is when Ntrain = 100,
order polynomial f (x) = 2x − 10x5 + 15x10 and ask one can often get better predictive performance by using
how the size of the training dataset Ntrain and the noise a less expressive model (e.g., a lower order polynomial)
strength σ affect the ability to make predictions. Obvi- rather than the more complex model (e.g., the tenth-
ously, more data and less noise leads to better predic- order polynomial). The simpler model has more “bias”
tions. To train the models (linear, third-order, tenth- but is less dependent on the particular realization of the
order), we uniformly sampled the interval x ∈ [0, 1] and training dataset, i.e. less “variance”. Finally we note that
constructed Ntrain training examples using (1). We then even with ten thousand data points, the model’s perfor-
fit the models on these training samples using standard mance quickly degrades beyond the original training data
least-squares regression. To visualize the performance of range. This demonstrates the difficulty of predicting be-
the three models, we plot the predictions using the best yond the training data we mentioned earlier.
fit parameters for a test set where x are drawn uniformly This simple example highlights why ML is so difficult
from the interval x ∈ [0, 1.2]. Notice that the test interval and holds some universal lessons that we will encounter
is slightly larger than the training interval. repeatedly in this review:
Figure 1 shows the results of this procedure for the
noiseless case, σ = 0. Even using a small training set • Fitting is not predicting. Fitting existing data well
with Ntrain = 10 examples, we find that the model class is fundamentally different from making predictions
that generated the data also provides the best fit and the about new data.
most accurate out-of-sample predictions. That is, the
linear model performs the best for data generated from a
• Using a complex model can result in overfitting. In-
linear polynomial (the third and tenth order polynomials
creasing a model’s complexity (i.e number of fitting
perform similarly), and the tenth order model performs
parameters) will usually yield better results on the
the best for data generated from a tenth order polyno-
training data. However when the training data size
mial. While this may be expected, the results are quite
is small and the data are noisy, this results in over-
different for larger noise strengths.
fitting and can substantially degrade the predictive
Figure 2 shows the results of the same procedure for
performance of the model.
noisy data, σ = 1, and a larger training set, Ntrain = 100.
As in the noiseless case, the tenth order model provides
the best fit to the data (i.e., the lowest Ein ). In contrast, • For complex datasets and small training sets, sim-
the tenth order model now makes the worst out-of-sample ple models can be better at prediction than com-
predictions (i.e., the highest Eout ). Remarkably, this is plex models due to the bias-variance tradeoff. It
true even if the data were generated using a tenth order takes less data to train a simple model than a com-
polynomial. plex one. Therefore, even though the correct model
At small sample sizes, noise can create fluctuations in is guaranteed to have better predictive performance
the data that look like genuine patterns. Simple mod- for an infinite amount of training data (less bias),
els (like a linear function) cannot represent complicated the training errors stemming from finite-size sam-
patterns in the data, so they are forced to ignore the pling (variance) can cause simpler models to out-
fluctuations and to focus on the larger trends. Complex perform the more complex model when sampling is
models with many parameters, such as the tenth order limited.
polynomial in our example, can capture both the global
trends and noise-generated patterns at the same time. In • It is difficult to generalize beyond the situations
this case, the model can be tricked into thinking that the encountered in the training data set.
10
y
Training
−2 Linear 0
Poly 3
−5
−4 Poly 10
−10
0.0 0.2 0.4 0.6 0.8 1.0 0.00 0.25 0.50 0.75 1.00 1.25
x x
FIG. 3 Fitting versus predicting for noisy data. Ntrain = 104 noisy data points (σ = 1) in the range x ∈ [0, 1] were
generated from a tenth-order polynomial. This data was fit using three model classes: linear models (red), all polynomials
of order 3 (yellow), all polynomials of order 10 (green) and used to make prediction on Ntest = 100 new data points with
xtest ∈ [0, 1.2](shown on right). The tenth order polynomial gives good predictions but the model’s predictive power quickly
degrades beyond the training data range.
III. BASICS OF STATISTICAL LEARNING THEORY general about the relationship between Ein and Eout ?
Surprisingly, the answer is ‘Yes’. We can in fact say
In this section, we briefly summarize and discuss the quite a bit. This is the domain of statistical learning
sense in which learning is possible, with a focus on su- theory, and we give a brief overview of the main results
pervised learning. We begin with an unknown function in this section. Our goal is to briefly introduce some of
y = f (x) and fix a hypothesis set H consisting of all func- the major ideas from statistical learning theory because
tions we are willing to consider, defined also on the do- of the important role they have played in shaping how we
main of f . This set may be uncountably infinite (e.g. if think about machine learning. However, this is a highly
there are real-valued parameters to fit). The choice of technical and theoretical field, so we will just skim over
which functions to include in H usually depends on our some introductory topics. A more thorough introduction
intuition about the problem of interest. The function to statistical learning theory can be found in the intro-
f (x) produces a set of pairs (xi , yi ), i = 1 . . . N , which ductory textbook by Abu Mostafa (Abu-Mostafa et al.,
serve as the observable data. Our goal is to select a func- 2012).
tion from the hypothesis set h ∈ H that approximates
f (x) as best as possible, namely, we would like to find
h ∈ H such that h ≈ f in some strict mathematical A. Three simple schematics that summarize the basic
sense which we specify below. If this is possible, we say intuitions from Statistical Learning Theory
that we learned f (x). But if the function f (x) can, in
principle, take any value on unobserved inputs, how is it The basic intuitions of statistical learning can be sum-
possible to learn in any meaningful sense? marized in three simple schematics. The first schematic,
The answer is that learning is possible in the restricted shown in Figure 4, shows the typical out-of-sample er-
sense that the fitted model will probably perform approx- ror, Eout , and in-sample error, Ein , as a function of the
imately as well on new data as it did on the training data. amount of training data. In making this graph, we have
Once an appropriate error function E is chosen for the assumed that the true data is drawn from a sufficiently
problem under consideration (e.g. sum of squared errors complicated distribution, so that we cannot exactly learn
in linear regression), we can define two distinct perfor- the function f (x). Hence, after a quick initial drop (not
mance measures of interest. The in-sample error, Ein , shown in figure), the in-sample error will increase with
and the out-of-sample or generalization error, Eout . Re- the number of data points, because our models are not
call from Sec II that both metrics are required due to the powerful enough to learn the true function we are seeking
distinction between fitting and predicting. to approximate. In contrast, the out-of-sample error will
This raises a natural question: Can we say something decrease with the number of data points. As the number
11
Optimum
E out
E out
Error
{
Variance
Variance
}
}
Error
Bias
E in
FIG. 5 Bias-Variance tradeoff and model complexity.
This schematic shows the typical out-of-sample error Eout as
function of the model complexity for a training dataset of fixed
Number of data points size. Notice how the bias always decreases with model com-
plexity, but the variance, i.e. fluctuation in performance due
FIG. 4 Schematic of typical in-sample and out-of- to finite size sampling effects, increases with model complex-
sample error as a function of training set size. The ity. Thus, optimal performance is achieved at intermediate
typical in-sample or training error, Ein , out-of-sample or gen- levels of model complexity.
eralization error, Eout , bias, variance, and difference of errors
as a function of the number of training data points. The
schematic assumes that the number of data points is large (in
fitting and predicting. Models with a large difference be-
particular, the schematic does not show the initial drop in
Ein for small amounts of data), and that our model cannot tween the in-sample and out-of-sample errors are said to
exactly fit the true function f (x). “overfit” the data. One of the lessons of statistical learn-
ing theory is that it is not enough to simply minimize
the training error, because the out-of-sample error can
of data points gets large, the sampling noise decreases still be large. As we will see in our discussion of regres-
and the training data set becomes more representative sion in Sec. VI, this insight naturally leads to the idea of
of the true distribution from which the data is drawn. “regularization”.
For this reason, in the infinite data limit, the in-sample The second schematic, shown in Figure 5, shows the
and out-of-sample errors must approach the same value, out-of-sample, or test, error Eout as a function of “model
which is called the “bias” of our model. complexity”. Model complexity is a very subtle idea
The bias represents the best our model could do if we and defining it precisely is one of the great achieve-
had an infinite amount of training data to beat down ments of statistical learning theory. In many cases, model
sampling noise. The bias is a property of the kind of complexity is related to the number of parameters we
functions, or model class, we are using to approximate are using to approximate the true function f (x)1 . In
f (x). In general, the more complex the model class we the example of polynomial regression discussed above,
use, the smaller the bias. However, we do not generally higher-order polynomials are more complex than the lin-
have an infinite amount of data. For this reason, to get ear model. If we consider a training dataset of a fixed
best predictive power it is better to minimize the out-of- size, Eout will be a non-monotonic function of the model
sample error, Eout , rather than the bias. As shown in complexity, and is generally minimized for models with
Figure 4, Eout can be naturally decomposed into a bias, intermediate complexity. The underlying reason for this
which measures how well we can hypothetically do in the is that, even though using a more complicated model
infinite data limit, and a variance, which measures the always reduces the bias, at some point the model be-
typical errors introduced in training our model due to comes too complex for the amount of training data and
sampling noise from having a finite training set. the generalization error becomes large due to high vari-
The final quantity shown in Figure 4 is the difference ance. Thus, to minimize Eout and maximize our predic-
between the generalization and training error. It mea- tive power, it may be more suitable to use a more bi-
sures how well our in-sample error reflects the out-of-
sample error, and measures how much worse we would
do on a new data set compared to our training data. For 1 There are, of course, exceptions. One neat example in the context
this reason, the difference between these errors is pre- of one-dimensional regression in given in (Friedman et al., 2001),
cisely the quantity that measures the difference between Figure 7.5.
12
" #
X
2
ED, [C(y, f (X; θ̂D ))] = ED, (yi − f (xi ; θ̂D ))
i
" #
X
= ED, (yi − f (xi ) + f (xi ) − f (xi ; θ̂D ))2
i
X
= E [(yi − f (xi ))2 ] + ED, [(f (xi ) − f (xi ; θ̂D ))2 ] + 2E [yi − f (xi )]ED [f (xi ) − f (xi ; θ̂D )]
i
X
= σ2 + ED [(f (xi ) − f (xi ; θ̂D ))2 ], (5)
i
where in the last line we used the fact that our noise has zero mean and variance σ2 and the sum over i applies to all
13
The first term is called the bias However, in practice gradient descent is full of surprises
X and a series of ingenious tricks have been developed by
Bias2 = (f (xi ) − ED [f (xi ; θ̂D )])2 (7) the optimization and machine learning communities to
i
improve the performance of these algorithms.
and measures the deviation of the expectation value of The underlying reason why training a machine learn-
our estimator (i.e. the asymptotic value of our estimator ing algorithm is difficult is that the cost functions we
in the infinite data limit) from the true value. The second wish to optimize are usually complicated, rugged, non-
term is called the variance convex functions in a high-dimensional space with many
X local minima. To make things even more difficult, we
V ar = ED [(f (xi ; θ̂D ) − ED [f (xi ; θ̂D )])2 ], (8) almost never have access to the true function we wish
i
to minimize: instead, we must estimate this function di-
and measures how much our estimator fluctuates due rectly from data. In modern applications, both the size
to finite-sample effects. Combining these expressions, of the dataset and the number of parameters we wish to
we see that the expected out-of-sample error, Eout := fit is often enormous (millions of parameters and exam-
ED, [C(y, f (X; θ̂D ))], can be decomposed as ples). The goal of this chapter is to explain how gradient
descent methods can be used to train machine learning
Eout = Bias2 + V ar + N oise, (9) algorithms even in these difficult settings.
with N oise = i σ .
P 2 This chapter seeks to both introduce commonly used
The bias-variance tradeoff summarizes the fundamen- methods and give intuition for why they work. We
tal tension in machine learning, particularly supervised also include some practical tips for improving the per-
learning, between the complexity of a model and the formance of stochastic gradient descent (Bottou, 2012;
amount of training data needed to train it. Since data LeCun et al., 1998b). To help the reader gain more in-
is often limited, in practice it is often useful to use a tuition about gradient descent and its variants, we have
less-complex model with higher bias – a model whose developed a Jupyter notebook that allows the reader to
asymptotic performance is worse than another model – visualize how these algorithms perform on two dimen-
because it is easier to train and less sensitive to sampling sional surfaces. The reader is encouraged to experi-
noise arising from having a finite-sized training dataset ment with the accompanying notebook whenever a new
(smaller variance). This is the basic intuition behind the method is introduced (especially to explore how changing
schematics in Figs. 4, 5, and 6. hyper-parameters can affect performance). The reader
may also wish to consult useful reviews that cover these
topics (Ruder, 2016) and this blog [Link]
IV. GRADIENT DESCENT AND ITS GENERALIZATIONS optimizing-gradient-descent/.
η =0.1
η =0.5 1
−2.5
η =1 E(θ + v) ≈ E(θ) + ∇θ E(θ)v + vT H(θ)v,
2
η =1.01
−5.0 where H(θ) is the Hessian matrix of second derivatives.
−4 −2 0 2 4 Differentiating this equation respect to v and noting that
x for the optimal value vopt we expect ∇θ E(θ + vopt ) = 0,
yields the following equation
FIG. 7 Gradient descent exhibits three qualitatively
different regimes as a function of the learning rate. 0 = ∇θ E(θ) + H(θ)vopt . (12)
Result of gradient descent on surface z = x2 + y 2 − 1 for
learning rate of η = 0.1, 0.5, 1.01. Notice that the trajectory Rearranging this expression results in the desired update
converges to the global minima in multiple steps for small
rules for Newton’s method
learning rates (η = 0.1). Increasing the learning rate fur-
ther (η = 0.5) causes the trajectory to oscillate around the vt = H −1 (θt )∇θ E(θt ) (13)
global minima before converging. For even larger learning
rates (η = 1.01) the trajectory diverges from the minima. See θt+1 = θt − vt . (14)
corresponding notebook for details.
Since we have no guarantee that the Hessian is well con-
ditioned, in almost all applications of Netwon’s method,
data points, one replaces the inverse of the Hessian H −1 (θt ) by some
suitably regularized pseudo-inverse such as [H(θt )+I]−1
n
X with a small parameter (Battiti, 1992).
E(θ) = ei (xi , θ). (10) For the purposes of machine learning, Newton’s
i=1 method is not practical for two interrelated reasons.
For example, for linear regression ei is just the mean First, calculating a Hessian is an extremely expensive
square-error for data point i; for logistic regression, it is numerical computation. Second, even if we employ first-
the cross-entropy. To make analogy with physical sys- order approximation methods to approximate the Hes-
tems, we will often refer to this function as the “energy”. sian (commonly called quasi-Newton methods), we must
In the simplest gradient descent (GD) algorithm, we store and invert a matrix with n2 entries, where n is the
update the parameters as follows. Initialize the parame- number of parameters. For models with millions of pa-
ters to some value θ0 and iteratively update the param- rameters such as those commonly employed in the neu-
eters according to the equation ral network literature, this is close to impossible with
present-day computational power. Despite these practi-
vt = ηt ∇θ E(θt ), cal shortcomings, Newton’s method gives many impor-
θt+1 = θt − vt (11) tant intuitions about how to modify GD algorithms to
improve their performance. Notice that, unlike in GD
where ∇θ E(θ) is the gradient of E(θ) w.r.t. θ and we where the learning rate is the same for all parameters,
have introduced a learning rate, ηt , that controls how big Newton’s method automatically “adapts” the learning
a step we should take in the direction of the gradient at rate of different parameters depending on the Hessian
time step t. It is clear that for sufficiently small choice of matrix. Since the Hessian encodes the curvature of the
the learning rate ηt this methods will converge to a local surface we are trying to find the minimum of – more
minimum (in all directions) of the cost function. How- specifically, the singular values of the Hessian are in-
ever, choosing a small ηt comes at a huge computational versely proportional to the squares of the local curvatures
cost. The smaller ηt , the more steps we have to take to of the surface – Newton’s method automatically adjusts
reach the local minimum. In contrast, if ηt is too large, the step size so that one takes larger steps in flat di-
we can overshoot the minimum and the algorithm be- rections with small curvature and smaller steps in steep
comes unstable (it either oscillates or even moves away directions with large curvature.
from the minimum). This is shown in Figure 7. In prac- Our derivation of Newton’s method also allows us to
tice, one usually specifies a “schedule” that decreases ηt develop intuition about the role of the learning rate in
at long times. Common schedules include power law and GD. Let us first consider the special case of using GD
exponential decay in time. to find the minimum of a quadratic energy function of
To better understand this behavior and highlight some a single parameter θ (LeCun et al., 1998b). Given the
of the shortcomings of GD, it is useful to contrast GD current value of our parameter θ, we can ask what is
15
• GD is very sensitive to choices of the learning rates. really experiment with different methods in landscapes
As discussed above, GD is extremely sensitive to of varying complexity using the accompanying notebook.
the choice of learning rates. If the learning rate is
very small, the training process takes an extremely
long time. For larger learning rates, GD can di- C. Stochastic Gradient Descent (SGD) with mini-batches
verge and give poor results. Furthermore, depend-
ing on what the local landscape looks like, we have One of the most widely-applied variants of the gra-
to modify the learning rates to ensure convergence. dient descent algorithm is stochastic gradient descent
Ideally, we would “adaptively” choose the learning (SGD)(Bottou, 2012; Williams and Hinton, 1986). As
rates to match the landscape. the name suggests, unlike ordinary GD, the algorithm
is stochastic. Stochasticity is incorporated by approx-
• GD treats all directions in parameter space uni- imating the gradient on a subset of the data called a
formly. Another major drawback of GD is that minibatch 2 . The size of the minibatches is almost al-
unlike Newton’s method, the learning rate for GD ways much smaller than the total number of data points
is the same in all directions in parameter space. For n, with typical minibatch sizes ranging from ten to a
this reason, the maximum learning rate is set by the few hundred data points. If there are n points in total,
behavior of the steepest direction and this can sig- and the mini-batch size is M , there will be n/M mini-
nificantly slow down training. Ideally, we would like batches. Let us denote these minibatches by Bk where
to take large steps in flat directions and small steps k = 1, . . . , n/M . Thus, in SGD, at each gradient descent
in steep directions. Since we are exploring rugged step we approximate the gradient using a single mini-
landscapes where curvatures change, this requires batch Bk ,
us to keep track of not only the gradient but second n
X X
derivatives of the energy function (note as discussed ∇θ E(θ) = ∇θ ei (xi , θ) −→ ∇θ ei (xi , θ). (19)
above, the ideal scenario would be to calculate the i=1 i∈Bk
Hessian but this proves to be too computationally We then cycle over all k = 1, . . . , n/M minibatches one
expensive). at a time, and use the mini-batch approximation to the
gradient to update the parameters θ at every step k. A
• GD is sensitive to initial conditions. One conse- full iteration over all n data points – in other words using
quence of the local nature of GD is that initial con- all n/M minibatches – is called an epoch. For notational
ditions matter. Depending on where one starts, one convenience, we will denote the mini-batch approxima-
will end up at a different local minimum. There- tion to the gradient by
fore, it is very important to think about how one X
initializes the training process. This is true for GD ∇θ E M B (θ) = ∇θ ei (xi , θ). (20)
as well as more complicated variants of GD intro- i∈Bk
duced below. With this notation, we can rewrite the SGD algorithm as
• GD can take exponential time to escape saddle vt = ηt ∇θ E M B (θ),
points, even with random initialization. As we men- θt+1 = θt − vt . (21)
tioned, GD is extremely sensitive to the initial con- Thus, in SGD, we replace the actual gradient over the
dition since it determines the particular local min- full data at each gradient descent step by an approxima-
imum GD would eventually reach. However, even tion to the gradient computed using a minibatch. This
with a good initialization scheme, through random- has two important benefits. First, it introduces stochas-
ness (to be introduced later), GD can still take ex- ticity and decreases the chance that our fitting algorithm
ponential time to escape saddle points, which are gets stuck in isolated local minima. Second, it signifi-
prevalent in high-dimensional spaces, even for non- cantly speeds up the calculation as one does not have
pathological objective functions (Du et al., 2017). to use all n data points to approximate the gradient.
Indeed, there are modified GD methods developed Empirical and theoretical work suggests that SGD has
recently to accelerate the escape. The details of additional benefits. Chief among these is that introduc-
these boosted method are beyond the scope of this ing stochasticity is thought to act as a natural regular-
review, and we refer avid readers to (Jin et al., izer that prevents overfitting in deep, isolated minima
2017) for details. (Bishop, 1995a; Keskar et al., 2016).
In the next few subsections, we will introduce variants
of GD that address many of these shortcomings. These 2 Traditionally, SGD was reserved for the case where you train on
generalized gradient descent methods form the backbone a single example – in other words minibatches of size 1. However,
of much of modern deep learning and neural networks, we will use SGD to mean any approximation to the gradient on
see Sec IX. For this reason, the reader is encouraged to a subset of the data.
17
wt+∆t − 2wt + wt−∆t wt+∆t − wt E. Methods that use the second moment of the gradient
m +µ = −∇w E(w).
(∆t)2 ∆t
(25) In stochastic gradient descent, with and without mo-
Rearranging this equation, we can rewrite this as mentum, we still have to specify a “schedule” for tuning
the learning rate ηt as a function of time. As discussed in
(∆t)2 m the context of Newton’s method, this presents a number
∆wt+∆t = − ∇w E(w) + ∆wt . (26) of dilemmas. The learning rate is limited by the steepest
m + µ∆t m + µ∆t
direction which can change depending on the current po-
Notice that this equation is identical to Eq. (23) if we sition in the landscape. To circumvent this problem, ide-
identify the position of the particle, w, with the parame- ally our algorithm would keep track of curvature and take
ters θ. This allows us to identify the momentum param- large steps in shallow, flat directions and small steps in
eter and learning rate with the mass of the particle and steep, narrow directions. Second-order methods accom-
the viscous damping as: plish this by calculating or approximating the Hessian
and normalizing the learning rate by the curvature. How-
m (∆t)2 ever, this is very computationally expensive for models
γ= , η= . (27) with extremely large number of parameters. Ideally, we
m + µ∆t m + µ∆t
18
would like to be able to adaptively change the step size respectively, and (βj )t denotes βj to the power t. The
to match the landscape without paying the steep compu- parameters η and have the same role as in RMSprop.
tational price of calculating or approximating Hessians. Like in RMSprop, the effective step size of a parameter
Recently, a number of methods have been introduced depends on the magnitude of its gradient squared. To
that accomplish this by tracking not only the gradient, understand this better, let us rewrite this expression in
but also the second moment of the gradient. These terms of the variance σt2 = ŝt − (m̂t )2 . Consider a single
methods include AdaGrad (Duchi et al., 2011), AdaDelta parameter θt . The update rule for this parameter is given
(Zeiler, 2012), RMSprop (Tieleman and Hinton, 2012), by
and ADAM (Kingma and Ba, 2014). Here, we discuss
the last two as representatives of this class of algorithms. m̂t
∆θt+1 = −ηt p . (32)
In RMSprop, in addition to keeping a running average σt2 + m̂2t +
of the first moment of the gradient, we also keep track of
the second moment denoted by st = E[gt2 ]. The update We now examine different limiting cases of this expres-
rule for RMSprop is given by sion. Assume that our gradient estimates are consistent
so that the variance is small. In this case our update
gt = ∇θ E(θ) (29) rule tends to ∆θt+1 → −ηt (here we have assumed that
st = βst−1 + (1 − β)gt2 m̂t ). This is equivalent to cutting off large persis-
gt tent gradients at 1 and limiting the maximum step size
θt+1 = θt − ηt √ ,
st + in steep directions. On the other hand, imagine that the
gradient is widely fluctuating between gradient descent
where β controls the averaging time of the second mo-
steps. In this case σ 2 m̂2t so that our update becomes
ment and is typically taken to be about β = 0.9, ηt is a
∆θt+1 → −ηt m̂t /σt . In other words, we adapt our learn-
learning rate typically chosen to be 10−3 , and ∼ 10−8
ing rate so that it is proportional to the signal-to-noise
is a small regularization constant to prevent divergences.
ratio (i.e. the mean in units of the standard deviation).
Multiplication and division by vectors is understood as
From a physics standpoint, this is extremely desirable:
an element-wise operation. It is clear from this formula
the standard deviation serves as a natural adaptive scale
that the learning rate is reduced in directions where the
for deciding whether a gradient is large or small. Thus,
gradient is consistently large. This greatly speeds up the
ADAM has the beneficial effects of (i) adapting our step
convergence by allowing us to use a larger learning rate
size so that we cut off large gradient directions (and hence
for flat directions.
prevent oscillations and divergences), and (ii) measuring
A related algorithm is the ADAM optimizer. In
gradients in terms of a natural length scale, the stan-
ADAM, we keep a running average of both the first and
dard deviation σt . The discussion above also explains
second moment of the gradient and use this information
empirical observations showing that the performance of
to adaptively change the learning rate for different pa-
both ADAM and RMSprop is drastically reduced if the
rameters. In addition to keeping a running average of the
square root is omitted in the update rule. It is also
first and second moments of the gradient (i.e. mt = E[gt ]
worth noting that recent studies have shown adaptive
and st = E[gt2 ], respectively), ADAM performs an addi-
methods like RMSProp, ADAM, and AdaGrad to gener-
tional bias correction to account for the fact that we are
alize worse than SGD in classification tasks, though they
estimating the first two moments of the gradient using a
achieve smaller training error. Such discussion is beyond
running average (denoted by the hats in the update rule
the scope of this review so we refer readers to (Wilson
below). The update rule for ADAM is given by (where
et al., 2017) for more details.
multiplication and division are once again understood to
be element-wise operations)
gt = ∇θ E(θ) (30) F. Comparison of various methods
mt = β1 mt−1 + (1 − β1 )gt
To better understand these methods, it is helpful to
st = β2 st−1 + (1 − β2 )gt2 visualize the performance of the five methods discussed
mt
m̂t = above – gradient descent (GD), gradient descent with
1 − (β1 )t momentum (GDM), NAG, ADAM, and RMSprop. To
st do so, we will use Beale’s function:
ŝt =
1 − (β2 )t
m̂t f (x, y) = (1.5 − x + xy)2 (33)
θt+1 = θt − ηt √ , 2 2 3 2
ŝt + +(2.25 − x + xy ) + (2.625 − x + xy ) .
(31)
This function has a global minimum at (x, y) = (3, 0.5)
where β1 and β2 set the memory lifetime of the first and and an interesting structure that can be seen in Fig. 9.
second moment and are typically taken to be 0.9 and 0.99 The figure shows the results of using all five methods
19
main idea of statistical mechanics in physics, where we look at the data then we would like to select an unin-
use probability to describe the behavior of large systems formative prior that reflects our ignorance, otherwise we
where we cannot know the positions and momenta of all should select an informative prior that accurately reflects
the particles even if the system itself is fully determinis- the knowledge we have about θ. This review will focus
tic (at least classically). In practice, Bayesian inference on informative priors that are commonly used for ML
provides a set of principles and procedures for learning applications. However, there is a large literature on un-
from data and for describing uncertainty. In this section, informative priors, including reparameterization invari-
we give a gentle introduction to Bayesian inference, with ant priors, that would be of interest to physicists and
special emphasis on its logic (i.e. Bayesian reasoning) and we refer the interested reader to (Berger and Bernardo,
provide a connection to ML discussed in Sec. II and III. 1992; Gelman et al., 2014; Jaynes, 1996; Jeffreys, 1946;
For a technical account of Bayesian inference in general, Mattingly et al., 2018).
we refer readers to (Barber, 2012; Gelman et al., 2014). Using an informative prior tends to decrease the vari-
ance of the posterior distribution while, potentially, in-
creasing its bias. This is beneficial if the decrease in
A. Bayes Rule variance is larger than the increase in bias. In high-
dimensional problems, it is reasonable to assume that
To solve a problem using Bayesian methods, we have many of the parameters will not be strongly relevant.
to specify two functions: the likelihood function p(X|θ), Therefore, many of the parameters of the model will
which describes the probability of observing a dataset be zero or close to zero. We can express this belief
X for a given value of the unknown parameters θ, and using two commonly used priors: the Gaussian prior
the prior distribution p(θ), which describes any knowl- Q q λ −λθ2
p(θ|λ) = j 2π e j is used to express the assump-
edge we have about the parameters before we collect the
tion that many of the parameters will be small, and the
data. Note that the likelihood should be considered as
Laplace prior p(θ|λ) = j λ2 e−λ|θj | is used to express the
Q
a function of the parameters θ with the data X held
assumption that many of the parameters will be zero.
fixed. The prior distribution and the likelihood function
We’ll come back to this point later in Sec. VI.F.
are used to compute the posterior distribution p(θ|X)
via Bayes’ rule:
B. Bayesian Decisions
p(X|θ)p(θ)
p(θ|X) = R . (34)
dθ 0 p(X|θ 0 )p(θ 0 ) The above section presents the tools for computing the
posterior distribution p(θ|X), which uses probability as
The posterior distribution describes our knowledge about
a framework for expressing our knowledge about the pa-
the unknown parameter θ after observing the data X.
rameters θ. In most cases, however, we need to summa-
In many cases, it will not be possible to analytically
rize our knowledge and pick a single “best” value for the
compute the normalizing constant in the R denominator of parameters. In principle, the specific value of the param-
the posterior distribution, i.e. p(X) = dθ p(X|θ)p(θ),
eters should be chosen to maximize a utility function.
and Markov Chain Monte Carlo (MCMC) methods are
In practice, however, we usually use one of two choices:
needed to draw random samples from p(θ|X).
the posterior mean hθi = dθ θp(θ|X), or the poste-
R
The likelihood function p(X|θ) is a common feature
rior mode θ̂MAP = arg maxθ p(θ|X). Often, hθi is called
of both classical statistics and Bayesian inference, and
is determined by the model and the measurement noise. the Bayes estimate and θ̂MAP is called the maximum-a-
Many common statistical procedures such as least-square posteriori or MAP estimate. While the Bayes estimate
fitting can be cast as Maximum Likelihood Estimation minimizes the mean-squared error, the MAP estimate is
(MLE). In MLE, one chooses the parameters θ̂ that max- often used instead because it is easier to compute.
imize the likelihood (or equivalently the log-likelihood
since log is a monotonic function) of the observed data:
C. Hyperparameters
θ̂ = arg max log p(X|θ). (35) The Gaussian and Laplace prior distributions, used to
θ
express the assumption that many of the model parame-
In other words, in MLE we choose the parameters that ters will be small or zero, both have an extra parameter
maximize the probability of seeing the observed data λ. This hyperparameter or nuisance variable has to be
given our generative model. MLE is an important con- chosen somehow. One standard Bayesian approach is to
cept in both frequentist and Bayesian statistics. define another prior distribution for λ – usually using an
The prior distribution, by contrast, is uniquely uninformative prior – and to average the posterior distri-
Bayesian. There are two general classes of priors: if we bution over all choices of λ. This is called a hierarchical
do not have any specialized knowledge about θ before we prior. Computing averages, however, often requires long
21
Markov Chain Monte Carlo simulations that are compu- the columns X:,j ∈ Rn , j = 1, · · · p being measured fea-
tationally intensive. Therefore, it is simpler if we can tures. Bear in mind that this function f is never known
find a good value of λ using an optimization procedure to us explicitly, though in practice we usually presume
instead. We will discuss how this is done in practice when its functional form. For example, in linear regression, we
discussing linear regression in Sec. VI. assume yi = f (x(i) ; wtrue ) + i = wtrue
T
x(i) + i for some
unknown but fixed wtrue ∈ R . p
in response to the applied electric current I. Normally In other words, we are looking to find the w which mini-
we would first make a bunch of measurements labeled mizes the L2 error. Geometrically speaking, the predictor
by i and plot them on a two-dimensional scatterplot, function g(x(i) ; w) = wT x(i) defines a hyperplane in Rp .
(Vi , Ii ). The next step is to assume, either from an oracle Minimizing the least squares error is therefore equivalent
or from theoretical reasoning, some models that might to minimizing the sum of all projections (i.e. residuals)
explain the measurements and measuring their perfor- for all points x(i) to this hyperplane (see Fig. 10). For-
mance. Mathematically, this amounts to finding some mally, we denote the solution to this problem as ŵLS :
function f such that Vi = f (Ii ; w), where w is some pa-
rameter (e.g. the electrical resistance R of the metal slab
in the case of Ohm’s law). We then try to minimize the ŵLS = arg min ||Xw − y||22 , (38)
errors made in explaining the given set of measurements w∈Rp
based on our model f by tuning the parameter w. To
which, after straightforward differentiation, leads to
do so, we need to first define the error function (formally
called the loss function) that characterizes the deviation ŵLS = (X T X)−1 X T y. (39)
of our prediction from the actual response.
Before formulating the problem, let us set up the no- Note that we have assumed that X T X is invertible,
tation. Suppose we are given a dataset with n sam- which is often the case when n p. Formally speak-
ples D = {(yi , x(i) )}ni=1 , where x(i) is the i-th obser- ing, if rank(X) = p, namely, the predictors X:,1 , . . . , X:,p
vation vector while yi is its corresponding (scalar) re- (i.e. columns of X) are linearly independent, then ŵLS is
sponse. We assume that every sample has p features, unique. In the case of rank(X) < p, which happens when
namely, x(i) ∈ Rp . Let f be the true function/model p > n, X T X is singular, implying there are infinitely
that generated these samples via yi = f (x(i) ; wtrue ) + i , many solutions to the least squares problem, Eq. (38).
where wtrue ∈ Rp is a parameter vector and i is some In this case, one can easily show that if w0 is a solution,
i.i.d. white noise with zero mean and finite variance. w0 +η is also a solution for any η which satisfies Xη = 0
Conventionally, we cast all samples into an n × p ma- (i.e. η ∈ null(X)). Having determined the least squares
trix, X ∈ Rn×p , called the design matrix, with the rows solution, we can calculate ŷ, the best fit of our data X,
Xi,: = x(i) ∈ Rp , , i = 1, · · · , n being observations and as ŷ = X ŵLS = PX y, where PX = X(X T X)−1 X T ,
22
c.f. Eq. (37). Geometrically, PX is the projection matrix regularization: the first one employs an L2 penalty and
which acts on y and projects it onto the column space is called Ridge regression, while the second uses an L1
of X, which is spanned by the predictors X:,1 , · · · , X:,p penalty and is called LASSO.
(see FIG. 11). Notice that we found the optimal solu-
tion ŵLS in one shot, without doing any sort of iterative
optimization like that discussed in Section IV. y
span( { X 1 , · · · , X p } )
y − ŷ
ŷ
tion of Ēin and Ēout in Section III, the average in-sample This problem is equivalent to the following constrained
and out-of-sample error can be shown to be optimization problem
p
Ēin = σ 2 1 − (40) ŵRidge (t) = arg min ||Xw − y||22 . (44)
n
p w∈Rp : ||w||22 ≤t
Ēout = σ 2 1 + , (41)
n This means that for any t ≥ 0 and solution ŵRidge in
provided we obtain the least squares solution ŵLS from Eq. (44), there exists a value λ ≥ 0 such that ŵRidge
i.i.d. samples X and y generated through y = Xwtrue + solves Eq. (43), and vice versa4 . With this equivalence, it
3 . Therefore, we can calculate the average generaliza- is obvious that by adding a regularization term, λ||w||22 ,
tion error explicitly: to our least squares loss function, we are effectively con-
p straining the magnitude of the parameter vector learned
|Ēin − Ēout | = 2σ 2 . (42) from the data.
n
To see this, let us solve Eq. (43) explicitly. Differenti-
This imparts an important message: if we have p n ating w.r.t. w, we obtain,
(i.e. high-dimensional data), the generalization error is
extremely large, meaning the model is not learning. Even ŵRidge (λ) = (X T X + λIp×p )−1 X T y. (45)
when we have p ≈ n, we might still not learn well due
to the intrinsic noise σ 2 . One way to ameliorate this
is, as we shall see in the following few sections, to use
4
regularization. We will mainly focus on two forms of Note that the equivalence between the penalized and the con-
strained (regularized) form of least square optimization does not
always hold. It holds for Ridge and LASSO (introduced later),
but not for best subset selection which is defined by choosing a
L0 norm: λ||w||0 . In this case, for every λ > 0 and any ŵBS
3 This requires that is a noise vector whose elements are i.i.d. of that solves the penalized form of best subset selection, there is a
zero mean and variance σ 2 , and is independent of the samples value t ≥ 0 such that ŵBS also solves that constrained form of
X. best subset selection, but the converse is not true.
23
In fact, when X is orthogonal, one can simplify this ex- As in Ridge regression, there is another formulation for
pression further: LASSO based on constrained optimization, namely,
ŵLS ŵLASSO (t) = arg min ||Xw − y||22 . (53)
ŵRidge (λ) = , for orthogonal X, (46)
1+λ w∈Rp : ||w||1 ≤t
where ŵLS is the least squares solution given by Eq. (39). The equivalence interpretation is the same as in Ridge
This implies that the ridge estimate is merely the least regression, namely, for any t ≥ 0 and solution ŵLASSO in
squares estimate scaled by a factor (1 + λ)−1 . Eq. (53), there is a value λ ≥ 0 such that ŵLASSO solves
Can we derive a similar relation between the fitted Eq. (52), and vice versa. However, to get the analytic
vector ŷ = X ŵRidge and the prediction made by least solution of LASSO, we cannot simply take the gradient
squares linear regression? To answer this, let us do a sin- of Eq. (52) with respect to w, since the L1 -regularizer is
gular value decomposition (SVD) on X. Recall that the not everywhere differentiable, in particular at any point
SVD of an n × p matrix X has the form where wj = 0 (see Fig. 13). Nonetheless, LASSO is a
convex problem. Therefore, we can invoke the so-called
X = U DV T , (47)
“subgradient optimality condition" (Boyd and Vanden-
where U ∈ Rn×p and V ∈ Rp×p are orthogonal matrices berghe, 2004; Rockafellar, 2015) in optimization theory
such that the columns of U span the column space of to obtain the solution. To keep the notation simple, we
X while the columns of V span the row space of X. only show the solution assuming X is orthogonal:
D ∈ Rp×p =diag(d1 , d2 , · · · , dp ) is a diagonal matrix with
ŵjLASSO (λ) = sign(ŵjLS )(|ŵjLS | − λ)+ , for orthogonal X,
entries d1 ≥ d2 ≥ · · · dp ≥ 0 called the singular values of
(54)
X. Note that X is singular if there is at least one dj = 0.
where (x)+ denotes the positive part of x and ŵjLS is
By writing X in terms of its SVD, one can recast the
the j-th component of least squares solution. In Fig. 12,
Ridge estimator Eq. (45) as
we compare the Ridge solution Eq. (46) with LASSO
ŵRidge = V (D 2 + λI)−1 DU T y, (48) solution Eq. (54). As we mentioned above, the Ridge
solution is the least squares solution scaled by a factor
which implies that the Ridge predictor satisfies of (1 + λ). Here LASSO does something conventionally
ŷRidge = X ŵRidge called “soft-thresholding" (see Fig. 12). We encourage
interested readers to work out the exercises in Notebook
= U D(D 2 + λI)−1 DU T y
p
3 to explore what this function does.
X d2j
= U:,j 2 U:jT y (49) LASSO Ridge
j=1
dj + λ
≤ U U Ty (50)
scaled by
= X ŷ ≡ ŷLS , (51) λ
λ 1+ λ
where U:,j are the columns of U . Note that in the in-
λ
equality step we assumed λ ≥ 0 and used SVD to sim-
plify Eq. (39). By comparing Eq. (49) with Eq. (51), it is
clear that in order to compute the fitted vector ŷ, both
Ridge and least squares linear regression have to project FIG. 12 [Adapted from (Friedman et al., 2001)] Comparing
y to the column space of X. The only difference is that LASSO and Ridge regression. The black 45 degree line is
Ridge regression further shrinks each basis component j the unconstrained estimate for reference. The estimators are
by a factor d2j /(d2j + λ). We encourage the reader to do shown by red dashed lines. For LASSO, this corresponds to
the exercises in Notebook 3 to develop further intuition the soft-thresholding function Eq. (54) while for Ridge regres-
about how Ridge regression works. sion the solution is given by Eq. (46)
Performance
0.6
cause the regressor contours are in the shape of ovals (it
is quadratic in w), their intersection tends to occur at 0.4
the vertex of the feasibility region, implying the solution
vector will be sparse. 0.2
0.0
2 1 0 1 2
ŵ 10 10 10 10 10
ŵ
w2 w2
FIG. 14 Performance of LASSO and ridge regression on the
t
diabetes dataset measured by the R2 coefficient of determina-
t tion. The best possible performance is R2 = 1. See Notebook
w1 w1 3.
D. Using Linear Regression to Learn the Ising Hamiltonian numbers Ej mean, and the configuration {Sj }L j=1 can be
interpreted in many ways: the outcome of coin tosses,
To gain deeper intuition about what kind of physics black-and-white pixels of an image, the binary represen-
problems linear regression allows us to tackle, consider tation of integers, etc. Your goal is to learn a model that
the following problem of learning the Hamiltonian for predicts Ej from the spin configurations.
the Ising model. Imagine you are given an ensemble of Without any prior knowledge about the origin of the
random spin configurations, and assigned to each state data set, physics intuition may suggest to look for a spin
its energy, generated from the 1D Ising model: model with pairwise interactions between every pair of
L
X variables. That is, we choose the following model class:
H = −J Sj Sj+1 (55)
L X
L
j=1 X
Hmodel [S i ] = − Jj,k Sji Ski , (56)
where J is the nearest-neighbor spin interaction, and j=1 k=1
Sj ∈ {±1} is a spin variable. Let’s assume the data
was generated with J = 1. You are handed the data The goal is to determine the interaction matrix Jj,k by
set D = ({Sj }L
j=1 , Ej ) without knowledge of what the applying linear regression on the data set D. This is a
25
well-defined problem, since the unknown Jj,k enters lin- data, was set to unity, J = 1. Moreover, Jj,k was NOT
early into the definition of the Hamiltonian. To this end, defined to be symmetric (we only used the Jj,j+1 but
we cast the above ansatz into the more familiar linear- never the Jj,j−1 elements). Figure. 17 shows the matrix
regression form: representation of the learned weights Jj,k . Interestingly,
OLS and Ridge regression learn nearly symmetric weights
Hmodel [S i ] = Xi · J. (57) J ≈ −0.5. This is not surprising, since it amounts to tak-
ing into account both the Jj,j+1 and the Jj,j−1 terms, and
The vectors Xi represent all two-body interactions
the weights are distributed symmetrically between them.
j,k=1 , and the index i runs over the samples in
{Sji Ski }L
LASSO, on the other hand, tends to break this symme-
the dataset. To make the analogy complete, we can also
try (see matrix elements plots for λ = 0.01) 5 . Thus,
represent the dot product by a single index p = {j, k},
we see how different regularization schemes can lead to
i.e. Xi · J = Xpi Jp . Note that the regression model does
learning equivalent models but in different ‘gauges’. Any
not include the minus sign. In the following, we apply
information we have about the symmetry of the unknown
ordinary least squares, Ridge, and LASSO regression to
model that generated the data should be reflected in the
the problem, and compare their performance.
definition of the model and the choice of regularization.
In addition to the diabetes dataset in Notebook 3, we
1.0
encourage the reader to work out Notebook 4 in which
linear regression is applied to the one-dimensional Ising
0.8
Train (OLS) model.
0.6 Test (OLS)
R2
FIG. 17 Learned interaction matrix Jij for the Ising model ansatz in Eq. (56) for ordinary least squares (OLS) regression
(left), Ridge regression (middle) and LASSO (right) at different regularization strengths λ. OLS is λ-independent but is shown
for comparison [Link] Notebook 4.
27
readers to (Boyd and Vandenberghe, 2004; Rockafellar, Using the assumption that samples are i.i.d., we can write
2015). the log-likelihood as
Now let us examine the two regularizers we introduced n
earlier. A close inspection reveals that LASSO and Ridge l(θ) ≡ log p(D|θ) =
X
log p(yi |x(i) , θ). (63)
regressions are both convex problems but only Ridge re- i=1
gression is a strictly convex problem (assuming λ > 0).
From convexity theory, this means that we always have a Note that the conditional dependence of the response
unique solution for Ridge but not necessary for LASSO. variable yi on the independent variable x(i) in the likeli-
In fact, it was recently shown that under mild conditions, hood function is made explicit since in regression the ob-
such as demanding general position for columns of X, served value of data, yi , is predicted based on x(i) using
the LASSO solution is indeed unique (Tibshirani et al., a model that is assumed to be a probability distribution
2013). Apart from this theoretical characterization, (Zou that depends on unknown parameter θ. This distribu-
and Hastie, 2005) introduced the notion of Elastic Net to tion, when endowed with θ, can, as we hope, potentially
retain the desirable properties of both LASSO and Ridge explain our prediction on yi . By definition, such distri-
regression, which is now one of the standard tools for bution is the likelihood function we discussed in Sec. V.
regression analysis and machine learning. We refer to Note that this is consistent with the formal statistical
reader to explore this in Notebook 2. treatment of regression where the goal is to estimate the
conditional expectation of the dependent variable given
the value of the independent variable (sometimes called
F. Bayesian formulation of linear regression the covariate) (Wasserman, 2013). We stress that this
notation does not imply x(i) is unknown– it is still part
In Section V, we gave an overview of Bayesian inference of the observed data!
and phrased it in the context of learning and uncertainty Using Eq. (61), we get
quantification. In this section we formulate least squares n 2 n
regression from a Bayesian point of view. We shall see 1 X
yi − wT x(i) − log 2πσ 2
l(θ) = − 2
that regularization in learning will emerge naturally as 2σ i=1 2
part of the Bayesian inference procedure. 1
From the setup of linear regression, the data D used to =− ||Xw − y||22 + const. (64)
2σ 2
fit the regression model is generated through y = xT w +
. We often assume that is a Gaussian noise with mean By comparing Eq. (38) and Eq. (64), it is clear that per-
zero and variance σ 2 . To connect linear regression to the forming least squares is the same as maximizing the log-
Bayesian framework, we often write the model as likelihood of this model.
What about adding regularization? In Section V, we
p(y|x, θ) = N (y|µ(x), σ 2 (x)). (61) introduced the maximum a posteriori probability (MAP)
estimate. Here we show that it actually corresponds to
In other words, our regression model is defined by a con- regularized linear regression, where the choice of prior
ditional probability that depends not only on data x but determines the type of regularization. Recall Bayes’ rule
on some model parameters θ. For example, if the mean
is a linear function of x given by µ = xT w, and the p(θ|D) ∝ p(D|θ)p(θ). (65)
variance is fixed σ 2 (x) = σ 2 , then θ = (w, σ 2 ). Now instead of maximizing the log-likelihood, l(θ) =
In statistics, many problems rely on estimation of some log p(D|θ), let us maximize the log posterior, log p(θ|D).
parameters of interest. For example, suppose we are Invoking Eq. (65), the MAP estimator becomes
given the height data of 20 junior students from a regional
high school, but what we are interested in is the average θ̂MAP ≡ arg max log p(D|θ) + log p(θ). (66)
height of all high school juniors in the whole county. It θ
is conceivable that the data we are given are not rep- In Sec. V.C, we discussed that a common choice for the
resentative of the student population as a whole. It is prior is a Gaussian distribution. Consider the Gaussian
therefore necessary to devise a systematic way to pre- prior6 with zero mean and variance τ 2 , namely, p(w) =
form reliable estimation. Here we present the maximum
likelihood estimation (MLE), and show that MLE for θ
is the one that minimizes the mean squared error (MSE)
used in OLS, see [Link].A. 6 Indeed, a Gaussian prior is the conjugate prior that gives a Gaus-
MLE is defined by maximizing the log-likelihood sian posterior. For a given likelihood, conjugacy guarantees the
preservation of prior distribution at the posterior level. For ex-
w.r.t. the parameters θ: ample, for a Gaussian (Geometric) likelihood with a Gaussian
(Beta) prior, the posterior distribution is still Gaussian (Beta)
θ̂ ≡ arg max log p(D|θ). (62) distribution.
θ
28
N (wj |0, τ 2 ). Then, we can recast the MAP estimator probably aware of the time of day, that it is a Wednes-
Q
j
into day, your friend Alice being in attendance, your friend
Bob being absent with a cold, the country in which you
n n
1 X 1 X are doing the experiment, and the planet you are on, but
θ̂MAP ≡ arg max − 2 (yi − wT x(i) )2 − 2 w2
θ 2σ i=1 2τ j=1 j you almost assuredly haven’t written these down in your
notebook. Why not? The reason is because you entered
1 1 the classroom with strongly held prior beliefs that none of
= arg max − 2 ||Xw − y||2 − 2 ||w||2 . (67)
2 2
θ 2σ 2τ those things affect the physics which takes place in that
room. Even of the things you did write down in an effort
Note that we dropped constant terms that do not de-
to be a careful scientist you probably hold some doubt
pend on the maximization parameters θ. The equiva-
as to their importance to your result and what is serving
lence between MAP estimation with a Gaussian prior
you here is the intuition that probably only a few things
and Ridge regression is established by comparing Eq. (67)
matter in the physics of pendula. Hence again you are
and Eq. (44) with λ ≡ σ 2 /τ 2 . We relegate the analogous
approaching the experiment with prior beliefs about how
derivation for LASSO to an exercise in Notebook 3.
many features you will need to pay attention to in order
to predict what will happen when you swing an unknown
G. Recap and a general perspective on regularizers pendulum. This example might seem a bit contrived, but
the point is that we live in a high-dimensional world of
In this section, we explored least squares linear regres- information and while we have good intuition about what
sion with and without regularization. We motivated the to write down in our notebook for well-known problems,
need for regularization due to poor generalization, in par- often in the field of ML we cannot say with any confi-
ticular in the “high-dimensional limit" (p n). Instead dence a priori what the small list of things to write down
of showing the average in-sample and out-of-sample er- will be, but we can at least use regularization to help us
rors for the regularized problem explicitly, we conducted enforce that the list not be too long so that we don’t end
numerical experiments in Notebook 3 on the diabetes up predicting that the period of a pendulum depends on
dataset and showed that regularization typically leads Bob having a cold on Wednesdays.
to better generalization. Due to the equivalence between Of course, in both LASSO and Ridge regression there
the constrained and penalized form of regularized regres- is a parameter λ involved. In principle, this hyper-
sion (in LASSO and Ridge, but not generally true in cases parameter is usually predetermined, which means that
such as L0 penalization), we can regard the regularized it is not part of the regression process. As we saw in
regression problem as an un-regularized problem but on Fig. 15, our learning performance and solution depends
a constrained set of parameters. Since the size of the al- strongly on λ, thus it is vital to choose it properly. As
lowed parameter space (e.g. w ∈ Rp when un-regularized we discussed in Sec. V.C, one approach is to assume an
vs. w ∈ C ⊂ Rp when regularized) is roughly a proxy for uninformative prior on the hyper-parameters, p(λ), and
model complexity, solving the regularized problem is in average the posterior over all choices of λ following this
effect solving the un-regularized problem with a smaller distribution. However, this comes with a large computa-
model complexity class. This implies that we’re less likely tional cost. Therefore, it is simpler to choose the regular-
to overfit. ization parameter through some optimization procedure.
We also showed the connection between using a reg- We’d like to emphasize that linear regression can be
ularization function and the use of priors in Bayesian applied to model non-linear relationship between input
inference. This connection can be used to develop more and response. This can be done by replacing the input
intuition about why regularization implies we are less x with some nonlinear function φ(x). Note that doing
likely to overfit the data: Let’s say you are a young so preserves the linearity as a function of the parame-
Physics student taking a laboratory class where the goal ters w, since model is defined by the their inner product
of the experiment is to measure the behavior of several φT (x)w. This method is known as basis function expan-
different pendula and use that to predict the formula sion (Bishop, 2006; Murphy, 2012).
(i.e. model) that determines the period of oscillation. Recent years have also seen a surge of interest in un-
In your investigation you would probably record many derstanding generalized linear regression models from a
things (hopefully including the length and mass!) in an statistical physics perspective. Much of this research has
effort to give yourself the best possible chance of deter- focused on understanding high-dimensional linear regres-
mining the unknown relationship, perhaps writing down sion and compressed sensing (Donoho, 2006) (see (Ad-
the temperature of the room, any air currents, if the ta- vani et al., 2013; Zdeborová and Krzakala, 2016) for ac-
ble were vibrating, etc. What you have done is create a cessible reviews for physicists). On a technical level,
high-dimensional dataset for yourself. However you actu- this research imports and extends the machinery of spin
ally possess an even higher-dimensional dataset than you glass physics (replica method, cavity method, and mes-
probably would admit to yourself. For example you are sage passing) to analyze high-dimensional linear models
29
0 1 2 3
data point xi belongs to a category yi = {0, 1} is given Since the cost (error) function is just the negative log-
by likelihood, for logistic regression we find
where θ = w are the weights we wish to learn from the The right-hand side in Eq. (76) is known in statistics as
data. To gain some intuition for these equations, consider the cross entropy.
a collection of non-interacting two-state systems coupled Having specified the cost function for logistic regres-
to a thermal bath (e.g. a collection of atoms that can be sion, we note that, just as in linear regression, in practice
in two states). Furthermore, denote the state of system we usually supplement the cross-entropy with additional
i by a binary variable: yi ∈ {0, 1}. From elementary regularization terms, usually L1 and L2 regularization
statistical mechanics, we know that if the two states have (see Sec. VI for discussion of these regularizers).
energies 0 and 1 the probability for finding the system
in a state yi is:
B. Minimizing the cross entropy
e−β0 1
P (yi = 1) = −β0 = ,
e + e−β1 1 + e−β∆ The cross entropy is a convex function of the weights w
P (yi = 1) = 1 − P (yi = 0). (71) and, therefore, any local minimizer is a global minimizer.
Minimizing this cost function leads to the following equa-
Notice that in these expressions, as is often the case in tion
physics, only energy differences are observable. If the n
difference in energies between two states is given by ∆ = 0 = ∇C(w) =
X
σ(xTi w) − yi xi ,
(77)
xTi w, we recover the expressions for logistic regression. i=1
We shall use this mapping between partition functions
and classification to generalize the logistic regressor to where we made use of the logistic function identity
SoftMax regression in Sec. VII.D. Notice that in terms ∂z σ(s) = σ(s)[1−σ(s)]. Equation (77) defines a transcen-
of the logistic function, we can write dental equation for w, the solution of which, unlike linear
regression, cannot be written in a closed form. For this
P (yi = 1) = σ(xTi w) = 1 − P (yi = 0). (72) reason, one must use numerical methods such as those
introduced in Sec. IV to solve this optimization problem.
We now define the cost function for logistic regression
using Maximum Likelihood Estimation (MLE). Recall,
that in MLE we choose parameters to maximize the prob- C. Examples of binary classification
ability of seeing the observed data. Consider a dataset
D = {(yi , xi )} with binary labels yi ∈ {0, 1} from which Let us now show how to use logistic regression in
the data points are drawn independently. The likelihood practice. In this section, we showcase two pedagogi-
of observing the data under our model is just: cal examples to train a logistic regressor to classify bi-
nary data. Each example comes with a corresponding
n
Y yi 1−yi Jupyter notebook, see [Link] panka-
σ(xTi w) 1 − σ(xTi w)
P (D|w) = jm/[Link].
i=1
(73)
1. Identifying the phases of the 2D Ising model
from which we can readily compute the log-likelihood:
The goal of this example is to show how one can employ
n
X logistic regression to classify the states of the 2D Ising
yi log σ(xTi w) + (1 − yi ) log 1 − σ(xTi w) .
l(w) = model according to their phase of matter.
i=1
(74) The Hamiltonian for the classical Ising model is given
The maximum likelihood estimator is defined as the set by
of parameters that maximize the log-likelihood: X
H = −J Si Sj , Sj ∈ {±1}, (78)
n
X hiji
yi log σ(xTi w)+(1−yi ) log 1 − σ(xTi w) .
ŵ = arg max
θ i=1 where the lattice site indices i, j run over all nearest
(75) neighbors of a 2D square lattice, and J is an interaction
31
10 10 10
20 20 20
30 30 30
0 20 0 20 0 20
FIG. 20 Examples of typical states of the 2D Ising model for three different temperatures in the ordered phase (T /J = 0.75,
left), the critical region (T /J = 2.25, middle) and the disordered phase (T /J = 4.0, right). The linear system dimension is
L = 40 sites.
energy scale. We adopt periodic boundary conditions. to, among other things, critical slowing down of the MC
Onsager proved that this model undergoes a phase tran- algorithm. Perhaps identifying the phases is also harder
sition in the thermodynamic limit from an ordered fer- in the critical region. With this in mind, consider the
romagnet with all spins aligned to a disordered √ phase at following three types of states: ordered (T /J < 2.0),
the critical temperature Tc /J = 2/ log(1 + 2) ≈ 2.26. near-critical (2.0 ≤ T /J ≤ 2.5) and disordered (T /J >
For any finite system size, this critical point is smeared 2.5). We use both ordered and disordered states to train
out to a critical region around Tc . the logistic regressor and, once the supervised training
An interesting question to ask is whether one can procedure is complete, we will evaluate the performance
train a statistical classifier to distinguish between the two of our classification model on unseen ordered, disordered,
phases of the Ising model. If successful, this can be used and near-critical states.
to locate the position of the critical point in more compli- Here, we deploy the liblinear routine (the default for
cated models where an exact analytical solution has so far Scikit’s logistic regression) and stochastic gradient de-
remained elusive (Morningstar and Melko, 2017; Zhang scent (SGD, see Sec. IV for details) to optimize the logis-
et al., 2017b). In other words, given an Ising state, we tic regression cost function with L2 regularization. We
would like to classify whether it belongs to the ordered or define the accuracy of the classifier as the percentage of
the disordered phase, without any additional information correctly classified data points. Comparing the accuracy
other than the spin configuration itself. This categorical on the training and test data, we can study the degree
machine learning problem is well suited for logistic re- of overfitting. The first thing to notice in Fig. 21 is the
gression, and will thus consist of recognizing whether a small degree of overfitting, as suggested by the training
given state is ordered by looking at its bit configurations. (blue) and test (red) accuracy curves being very close to
Notice that, for the purposes of applying logistic regres- each other. Interestingly, the liblinear minimizer outper-
sion, the 2D spin state of the Ising model will be flattened forms SGD on the training and test data, but not on the
out to a 1D array, so it will not be possible to learn in- near-critical data for certain values of the regularization
formation about the structure of the contiguous ordered strength λ. Moreover, similar to the linear regression
2D domains [see Fig. 20]. Such information can be incor- examples, we find that there exists a sweet spot for the
porated using deep convolutional neural networks, see SGD regularization strength λ that results in optimal
Section IX. performance of the logistic regressor, at about λ ∼ 10−1 .
To this end, we consider the 2D Ising model on a 40×40 We might expect that the difficulty of the phase recogni-
square lattice, and use Monte-Carlo (MC) sampling to tion problem depends on the temperature of the queried
prepare 104 states at every fixed temperature T out of sample. Looking at the states in the near-critical region,
a pre-defined set. We furthermore assign a label to each c.f. Fig. 20, it is no longer easy for a trained human eye
state according to its phase: 0 if the state is disordered, to distinguish between the ferromagnetic and the disor-
and 1 if it is ordered. dered phases close to Tc . Therefore, it is interesting to
It is well-known that near the critical temperature Tc , also compare the training and test accuracies to the ac-
the ferromagnetic correlation length diverges, which leads curacy of the near-critical state predictions. (Recall that
32
FIG. 22 The probability of an event being a classified as a signal event for true signal events (left, blue) and background events
(right, red).
there is minimal overfitting, in part because we trained we’ve already come up with a set of discriminating vari-
on such a large dataset (4.5 million events). More impor- ables, including higher order ones derived from theories
tantly, however, is the underlying data we are working about SUSY particles, it’s worth reflecting on whether
with: each input variable is an important feature. there is utility to the increased sophistication of ML. To
show why we would want to use such a technique, recall
that, even to the learning algorithm, some signal events
and background events look similar. We can illustrate
this directly by looking at a plot comparing the pT spec-
trum of the two highest pT leptons (often referred to as
the leading and sub-leading leptons) for both signal and
background events. Figure 24 shows these two distribu-
tions, and one can see that while some signal events are
easily distinguished, many live in the same part of phase
space as the background. This effect can also be seen by
looking at figure 22 where you can see that some signal
events look like background events and vice-versa.
based in the field of HEP) to our regression results, differ- the probability of xi being in class m0 is given by
ent ROC curves have been made for each of the following T
e−xi wm0
cases: logistic regression with just the simple kinematic P (yim0 = 1|xi , {wk }M −1
k=0 ) = PM −1 −xT wm , (79)
variables, logistic regression with the full set of variables, m=0 e
i
and simply putting a requirement on the leading lepton where yim0 ≡ [yi ]m0 refers to the m0 -th component of vec-
pT . Figure 25 shows that there is a clear performance tor yi . This is known as the SoftMax function. There-
benefit from using logistic regression. Note also that in fore, the likelihood of this M -class classifier is simply
the cut-based approach we have only used one variable (cf. Sec. VII.A):
where we could have put requirements on all of them. n M −1
While putting more requirements would indeed increase yim
Y Y
P (D|{wk }M −1
k=0 ) = [P (yim = 1|xi , wm )]
background rejection, it would also decrease signal effi- i=1 m=0
ciency. Hence, the cut-based approach will never yield as 1−yim
× [1 − P (yim = 1|xi , wm )] (80)
strong discrimination as the logistic regression we have
performed. One other interesting point about these re- from which we can define the cost function in a similar
sults is that the higher-order variables noticeably help the fashion:
ML techniques. In later sections, we will return to this n M
X X −1
point to see if more sophisticated techniques can provide C(w) = − yim log P (yim = 1|xi , wm )
further improvement. i=1 m=0
+ (1 − yim ) log (1 − P (yim = 1|xi , wm )) . (81)
As expected, for M = 1, we recover the cross entropy for
logistic regression, cf. Eq. (76).
So far we have focused only on binary classification, One of the most powerful and widely-applied ideas in
in which the labels are dichotomous variables. Here we modern machine learning is the use of ensemble methods
generalize logistic regression to multi-class classification. that combine predictions from multiple, often weak, sta-
One approach is to treat the label as a vector yi ∈ ZM 2 ,
tistical models to improve predictive performance (Diet-
namely a binary string of length M with only one com- terich et al., 2000). Ensemble methods, such as random
ponent of yi being 1 and the rest zero. For example, forests (Breiman, 2001; Geurts et al., 2006; Ho, 1998),
yi = (1, 0, · · · , 0) means data the sample xi belongs to and boosted gradient trees, such as XGBoost (Chen and
class 17 , cf. Fig. 18. Following the notation in Sec. VII.A, Guestrin, 2016; Friedman, 2001), undergird many of the
winning entries in data science competitions such as Kag-
gle, especially on structured datasets 8 . Even in the con-
FIG. 26 Visualization of the weights wj after training a SoftMax Regression model on the MNIST dataset (see Notebook
7). We emphasize that SoftMax Regression does not have explicit 2D spatial knowledge; the model learns from data points
flattened out in a one-dimensional array.
In this section, we give an overview of ensemble meth- The bias-variance tradeoff summarizes the fundamen-
ods and provide rules of thumb for when and why they tal tension in machine learning between the complexity
work. On one hand, the idea of training multiple models of a model and the amount of training data needed to fit
and then using a weighted sum of the predictions of the it (see Sec. III). Since data is often limited, in practice
all these models is very natural. After all, the idea of the it is frequently useful to use a less complex model with
“wisdom of the crowds” can be traced back, at least, to higher bias – a model whose asymptotic performance is
the writings of Aristotle in Politics. On the other hand, worse than another model – because it is easier to train
one can also imagine that the ensemble predictions can and less sensitive to sampling noise arising from having a
be much worse than the predictions from each of the in- finite-sized training dataset (i.e. smaller variance). Here,
dividual models that constitute the ensemble, especially we will revisit the bias-variance tradeoff in the context
when pooling reinforces weak but correlated deficiencies of ensembles, drawing upon the beautiful discussion in
in each of the individual predictors. Thus, it is impor- Ref. (Louppe, 2014).
tant to understand when we expect ensemble methods to A key property that will emerge from this analysis is
work. the correlation between models that constitute the en-
In order to do this, we will revisit the bias-variance semble. The degree of correlation between models9 is
trade-off, discussed in Sec. III, and generalize it to con- important for two distinct reasons. First, holding the en-
sider an ensemble of classifiers. We will show that the semble size fixed, averaging the predictions of correlated
key to determining when ensemble methods work is the models reduces the variance less than averaging uncor-
degree of correlation between the models in the ensemble related models. Second, in some cases, correlations be-
(Louppe, 2014). Armed with this intuition, we will intro-
duce some of the most widely-used and powerful ensem-
ble methods including bagging (Breiman, 1996), boosting 9 For example, the correlation coefficient between the predictions
(Freund et al., 1999; Freund and Schapire, 1995; Schapire made by two randomized models based on the same training set
and Freund, 2012), random forests (Breiman, 2001), but with different random seeds, see Sec. VIII.A.1 for precise
and gradient boosted trees such as XGBoost (Chen and definition.
36
tween models within an ensemble can result in an in- measures the deviation of the expectation value of our es-
crease in bias, offsetting any potential reduction in vari- timator (i.e. the asymptotic value of our estimator in the
ance gained from ensemble averaging. We will discuss limit of infinite data) from the true value. The variance
this in the context of bagging below. One of the most X
dramatic examples of increased bias from correlations is V ar = EL [(ĝL (xi ) − EL [ĝL (xi )])2 ], (86)
the catastrophic predictive failure of almost all derivative i
models used by Wall Street during the 2008 financial cri- measures how much our estimator fluctuates due to
sis. finite-sample effects. The noise term
X
N oise = σ2i (87)
1. Bias-Variance Decomposition for Ensembles
i
We will discuss the bias-variance tradeoff in the con- is the part of the error due to intrinsic noise in the data
text of continuous predictions such as regression. How- generation process that no statistical estimator can over-
ever, many of the intuitions and ideas discussed here also come.
carry over to classification tasks. Before discussing en- Let us now generalize this to ensembles of estimators.
sembles, let us briefly review the bias-variance tradeoff in Given a dataset XL and hyper-parameters θ that param-
the context of a single model. Consider a data set consist- eterize members of our ensemble, we will consider a pro-
ing of data XL = {(yj , xj ), j = 1 . . . N }. Let us assume cedure that deterministically generates a model ĝL (xi , θ)
that the true data is generated from a noisy model given XL and θ. We assume that the θ includes some
random parameters that introduce stochasticity into our
y = f (x) + , (82) ensemble (e.g. an initial condition for stochastic gradient
where is a normally distributed with mean zero and descent or a random subset of features or data points
standard deviation σ . used for training.) Concretely, with a giving dataset L,
Assume that we have a statistical procedure (e.g. least- one has a learning algorithm A that generates a model
squares regression) for forming a predictor ĝL (x) that A(θ, L) based on a deterministic procedure which intro-
gives the prediction of our model for a new data point x duced stochasticity through θ in its execution on dataset
given that we trained the model using a dataset L. This L. We will be concerned with the expected prediction
estimator is chosen by minimizing a cost function which, error of the aggregate ensemble predictor
for the sake of concreteness, we take to be the squared M
error 1 X
A
ĝL (xi , {θ}) = ĝL (xi , θm ). (88)
X M m=1
C(X, g(x)) = (yi − ĝL (xi ))2 . (83)
i For future reference, let us define the mean, variance,
and covariance (i.e. the connected correlation function in
The dataset L is drawn from some underlying distribu-
the language of physics), and the normalized correlation
tion that describes the data. If we imagine drawing many
coefficient of a single randomized model ĝL (x, θm ) as:
datasets {Lj } of the same size as L from this distribu-
tion, we know that the corresponding estimators ĝLj (x) EL,θm [ĝL (x, θm )] = µL,θm (x)
will differ from each other due to stochastic effects aris- EL,θm [ĝL (x, θm ) ] − EL,θm [ĝL (x, θ)]2 = σL,θ
2 2
(x)
ing from sampling noise. For this reason, we can view m
our estimator ĝL (x) as a random variable and define an EL,θm [ĝL (x, θm )ĝL (x, θm0 )] − Eθ [ĝL (x, θm )]2 = CL,θm ,θm0 (x)
expectation value EL in the usual way. Note that the sub- CL,θm ,θm0 (x)
script denotes that the expectation is taken over L. In ρ(x) = 2 . (89)
σL,θ
practice, EL is computed by by drawing infinitely many
different datasets {Lj } of the same size, fitting the corre- Note that the expectation EL,θm [·] is computed over the
sponding estimator, and then averaging the results. We joint distribution of L and θm . Also, by definition, we
will also average over different instances of the “noise” . assume m 6= m0 in CL,θm ,θm0 .
The expectation value over the noise will be denoted by We can now ask about the expected generalization
E . (out-of-sample) error for the ensemble
As discussed in Sec. III, we can decompose the ex- " #
pected generalization error as
X
A A
(xi , {θ}))2 .
EL,,θ C(X, ĝL (x)) = EL,,θ (yi − ĝL
EL, [C(X, g(x))] = Bias2 + V ar + N oise. (84) i
(90)
where the bias, As in the single estimator case, we decompose the error
X into a noise term, a bias-term, and a variance term. To
Bias2 = (f (xi ) − EL [ĝL (xi )])2 , (85) see this, note that
i
37
" #
X
A A
EL,,θ [C(X, ĝL (x))] = EL,,θ (yi − f (xi ) + f (xi ) − ĝL (xi , {θ}))2
i
X
= EL,,θ [(yi − f (xi ))2 + (f (xi ) − ĝL
A
(xi , {θ}))2 + 2(yi − f (xi ))(f (xi ) − ĝL
A
(xi , {θ})]
i
X X
= σ2i + A
EL,θ [(f (xi ) − ĝL (xi , {θ}))2 ], (91)
i i
where in the last line we have used the fact that E [yi ] = f (xi ) to eliminate the last term. We can further decompose
the second term as
A
EL,θ [(f (xi ) − ĝL (xi , {θ}))2 ] = EL,θ [(f (xi ) − EL,θ [ĝL
A A
(xi , {θ})] + EL,θ [ĝL A
(xi , {θ})] − ĝL (xi , {θ}))2 ]
A
= EL,θ [(f (xi ) − EL,θ [ĝL (xi , {θ})])2 ] + EL,θ [(EL,θ [ĝL
A A
(xi , {θ})] − ĝL (xi , {θ}))2 ]
A A A
+ 2EL,θ [(EL,θ [ĝL (xi , {θ})] − ĝL (xi , {θ}))(f (xi ) − EL,θ [ĝL (xi , {θ})])
A
= (f (xi ) − EL,θ [ĝL (xi , {θ})])2 + EL,θ [(ĝL
A A
(xi , {θ}) − EL,θ [ĝL (xi , {θ})])2 ]
≡ Bias2 (xi ) + V ar(xi ), (92)
where we have defined the bias of an aggregate predictor So far the calculation for ensembles is almost iden-
as tical to that of a single estimator. However, since the
aggregate estimator is a sum of estimators, its variance
Bias2 (x) ≡ (f (x) − EL,θ [ĝL
A
(x, {θ})])2 (93) implicitly depends on the correlations between the indi-
vidual estimators in the ensemble. Using the definition
and the variance as of the aggregate estimator Eq. (88) and the definitions in
Eq. (89), we see that
A
V ar(x) ≡ EL,θ [(ĝL A
(x, {θ}) − EL,θ [ĝL (x, {θ})])2 ]. (94)
A A
V ar(x) = EL,θ [(ĝL (x, {θ}) − EL,θ [ĝL (x, {θ})])2 ]
1 X X
= 2 EL,θ [ĝL (x, θm )ĝL (x, θm0 )] − M 2 [µL,θ (x)]2
M 0 i
m,m
1 − ρ(x) 2
2
= ρ(x)σL,θ + σL,θ . (95)
M
This last formula is the key to understanding the power single model
of random ensembles. Notice that by using large ensem-
bles (M → ∞), we can significantly reduce the variance,
and for completely random ensembles where the mod-
Bias2 (x) = (f (x) − EL,θ [ĝL
A
(x, {θ}])2
els are uncorrelated (ρ(x) = 0), maximally suppresses
M
the variance! Thus, using the aggregate predictor beats 1 X
= (f (x) − EL,θ [ĝL (x, θm ])2 (96)
down fluctuations due to finite-sample effects. The key, M m=1
as the formula indicates, is to decorrelate the models as
much as possible while still using a very large ensemble. = (f (x) − µL,θ )2 . (97)
One can be worried that this comes at the expense of a
very large bias. This turns out not to be the case. When
models in the ensemble are completely random, the bias Thus, for a random ensemble one can always add more
of the aggregate predictor is just the expected bias of a models without increasing the bias. This observation lies
behind the immense power of random forest methods dis-
cussed below. For other methods, such as bagging, we
will see that the bootstrapping procedure actually does
increase the bias. But in many cases, this increase in bias
is negligible compared to the reduction in variance.
38
oretical discussion above, we know that this can signifi- a stable procedure the accuracy is limited by the bias
cantly reduce the variance without increasing the bias. introduced by using bootstrapped datasets. This means
While simple and intuitive, this form of aggregation that there is an instability-to-stability transition point
clearly works only when we have enough data in each par- beyond which bagging stops improving our prediction.
titioned set Li . To see this, one can consider the extreme
limit where Li contains exactly one point. In this case,
the base hypothesis gLi (x) (e.g. linear regressor) becomes
extremely poor and the procedure above fails. One way
to circumvent this shortcoming is to resort to empir-
ical bootstrapping, a resampling technique in statis-
tics introduced by Efron (Efron, 1979) (see accompany-
ing box and Fig. 28). The idea of empirical bootstrap-
ping is to use sampling with replacement to create new
“bootstrapped” datasets {LBS 1 , . . . , LM } from our origi-
BS
and
M
X
BS
ĝL (x) = arg max I[gLBS
i
(x) = j]. (101)
j M̂ n ( D )
i=1
small changes in the training set result in large changes for k = 1, · · · , B. These are then used to evaluate the accu-
in predictions (Breiman, 1996). When the procedure is racy of M̂n (D) (see also box on Bootstrapping in main text).
unstable, the prediction error is dominated by the vari- It can be shown that in the n → ∞ limit the distribution
?(k)
ance and one can exploit the aggregation component of of Mn would be a Gaussian centered around M̂n (D) with
bagging to reduce the prediction error. In contrast, for variance σ 2 defined by Eq. (102) scales as 1/n.
40
x2
placement from D to get a new set D?(1) =
?(1) ?(1) 0
{X1 , · · · , Xn }, called a bootstrap sample,
which possibly contains repetitive elements. Then −1
we repeat the same procedure to get in total B
−2
such sets: D?(1) , · · · , D?(B) . The next step is to
use these B bootstrap sets to get the bootstrap −3
−1.0 −0.5 0.0 0.5 1.0
estimate of the quantity of interest. For example, x1
?(k)
let Mn = M edian(D?(k) ) be the sample median
of bootstrap data D?(k) . Then we can construct FIG. 29 Bagging applied to the perceptron learning
the variance of the distribution of bootstrap medi- algorithm (PLA). Training data size n = 500, number of
ans as : bootstrap datasets B = 25, each contains 50 points. Col-
ors corresponds to different classes while the marker indicates
B 2 how these points are labelled: cross for true label and circle
1 X ?(k) for that obtained by bagging. Each gray dashed line indicates
Vd
arB (Mn ) = Mn − M̄n? , (102)
B−1 the prediction made, based on every bootstrap set while the
k=1
dark dashed black line is the average of these.
where
B
1 X ?(k) more ‘democratic’ ways as is done in bagging. In all cases,
M̄n? = Mn (103) the idea is to build a strong predictor by combining many
B
k=1
weaker classifiers.
is the mean of the median of all bootstrap sam- In boosting, an ensemble of weak classifiers {gk (x)}
ples. Specifically, Bickel and Freedman (Bickel and is combined into an aggregate, boosted classifier. How-
Freedman, 1981) and Singh (Singh, 1981) showed ever, unlike bagging, each classifier is associated with a
that in the n → ∞ limit, the distribution of the weight αk that indicates how much it contributes to the
bootstrap estimate will be a Gaussian centered aggregate classifier
around M̂n (D) = M edian(X1 , · · · √
, Xn ) with stan- M
dard deviation proportional to 1/ n. This means gA (x) =
X
αk gk (x), (104)
that the bootstrap distribution M̂n? − M̂n approxi- K=1
mates fairly well the sampling distribution M̂n −M
from which we obtain the training data D. Note where k αk = 1. For the reasons outlined above, boost-
P
that M is the median based on which the true dis- ing, like all ensemble methods, works best when we com-
tribution D is generated. In other words, if we plot bine simple, high-variance classifiers into a more complex
?(k)
the histogram of {Mn }B k=1 , we will see that in
whole.
the large n limit it can be well fitted by a Gaus- Here, we focus on “adaptive boosting” or AdaBoost,
sian which sharp peaks at M̂n (D) and vanishing first proposed by Freund and Schapire in the mid 1990s
variance whose definition is given by Eq. (102) (see (Freund et al., 1999; Freund and Schapire, 1995; Schapire
Fig. 28). and Freund, 2012). The basic idea behind AdaBoost,
is to form the aggregate classifier in an iterative pro-
cess. Importantly, at each iteration we reweight the error
C. Boosting function to "highlight" data points where the aggregate
classifier performs poorly (so that in the next round the
Another powerful and widely used ensemble method is procedure put more emphasis on making those right.) In
Boosting. In bagging, the contribution of all predictors this way, we can successively ensure that our classifier
is weighted equally in the bagged (aggregate) predictor. has good performance over the whole dataset.
However, in principle, there are myriad ways to combine We now discuss the AdaBoost procedure in greater
different predictors. In some problems one might prefer detail. Suppose that we are given a data set L =
to use an autocratic approach that emphasizes the best {(xi , yi ), i = 1, · · · , N } where xi ∈ X and yi ∈ Y =
predictors, while in others it might be better to opt for {+1, −1}. Our objective is to find an optimal hypoth-
41
esis/classifier g : X → Y to classify the data. Let method widely deployed for complex classification tasks.
H = {g : X → Y} be the family of classifiers available A random forest is composed of a family of (randomized)
in our ensemble. In the AdaBoost setting, we are con- tree-based classifier decision trees (discussed below). De-
cerned with the classifiers that perform somehow better cision trees are high-variance, weak classifiers that can
than “tossing a fair coin”. This means that for each clas- be easily randomized, and as such, are ideally suited for
sifier, the family H can predict yi correctly at least half ensemble-based methods. Below, we give a brief high-
of the time. level introduction to these ideas.
We construct the boosted classifier as follows: A decision tree uses a series of questions to hierarchi-
cally partition the data. Each branch of the decision tree
• Initialize wt=1 (xn ) = 1/N, n = 1, · · · , N .
consists of a question that splits the data into smaller
• For t = 1 · · · , T (desired termination step), do: subsets (e.g. is some feature larger than a given num-
ber? See Fig. 30), with the leaves (end points) of the
1. Select a hypothesis gt ∈ H that minimizes the tree corresponding to the ultimate partitions of the data.
weighted error When using decision trees for classification, the goal is
N
to construct trees such that the partitions are informa-
X tive about the class label (see Fig. 30). It is clear that
t = wt (xi )1(gt (xi ) 6= yi ) (105)
more complex decision trees lead to finer partitions that
i=1
give improved performance on the training set. How-
2. Let αt = 21 ln 1− ever, this generally leads to over-fitting10 , limiting the
t , update the weight for each
t
where the index i runs over data points and the index j
Before we turn to applications of these techniques, we runs over decision trees in our ensemble. In XGBoost,
briefly discuss one final class of ensemble methods that the regularization function is chosen to be
has become increasingly popular in the last few years:
Gradient-Boosted Trees (Chen and Guestrin, 2016; Fried- λ
Ω(g) = γT + ||w||22 , (108)
man, 2001). The basic idea of gradient-boosted trees is 2
to use intuition from boosting and gradient descent (in with γ and λ regularization parameters that must be
particular Newton’s method, see Sec. IV) to construct chosen appropriately. Notice that this regularization pe-
43
nalizes both large weights on the leaves (similar to L2 - It is clear that ∆Ctopt measures the in-sample performance
regularization) and having large partitions with many of gt and we should find the decision tree that minimizes
leaves. this value. In principle, one could enumerate all possi-
As in boosting, we form the ensemble iteratively. For ble trees over the data and find the tree that minimizes
(t)
this reason, we define a family of predictors ŷi as ∆Ctopt . However, in practice this is impossible. Instead,
an approximate greedy algorithm is run that optimizes
t
(t)
X (t−1) one level of the tree at a time by trying to find optimal
ŷi = gj (xi ) = ŷi + gt (xi ). (109)
splits of the data. This leads to a tree that is a good
j=1
local minimum of ∆Ctopt which is then added to the en-
Note that by definition yi
(M )
= gA (xi ). The central idea semble. We emphasize that this is only a very high level
is that for large t, each decision tree is a small pertur- sketch of how the algorithm works. In practice, addi-
bation to the predictor (of order 1/T ) and hence we can tional regularization such as shrinkage(Friedman, 2002)
perform a Taylor expansion on our loss function to second and feature subsampling(Breiman, 2001; Friedman et al.,
order: 2003) is also used. In addition, there are many numerical
N
and technical tricks used for the approximate algorithm
X (t−1) and how to find splits of the data that give good decision
Ct = l(yi , ŷi + gt (xi )) + Ω(gt )
trees (Chen and Guestrin, 2016).
i=1
≈ Ct−1 + ∆Ct , (110)
with F. Applications to the Ising model and Supersymmetry
(t−1) 1 Datasets
∆Ct = ai l(yi , ŷi )gt (xi ) + bi gt (xi )2 + Ω(gt ), (111)
2
We now illustrate some of these ideas using two exam-
where ples drawn from physics: (i) classifying the phases of the
ai = ∂ŷ(t−1) l(yi , ŷi
(t−1)
), (112) spin configurations of the 2D-Ising model above and be-
i
low the critical temperature using random forests and (ii)
classifying Monte-Carlo simulations of collision events in
(t−1)
bi = ∂ŷ2(t−1) l(yi , ŷi ). (113) the SUSY dataset as supersymmetric or standard using
i
an XGBoost implementation of gradient-boosted trees.
We then choose the t-th decision tree gt to minimize ∆Ct . Both examples were analyzed in Sec. VII.C using logistic
This is almost identical to how we derived the Newton regression. Here we show that on the Ising dataset, the
method update in the section on gradient descent, see RFs perform significantly better than logistic regression
Sec. IV. models whereas gradient boosted trees seem to yield an
We can actually derive an expression for the param- accuracy of about 80%, comparable to published results.
eters of gt that minimize ∆Ct analytically. To simplify The two accompanying Jupyter notebooks discuss practi-
notation, it is useful to define the set of points xi that get cal details of implementing these examples and the read-
mappedP to leaf j: Ij = {i P : qt (xi ) = j} and the functions ers are encouraged to experiment with the notebooks.
Bj = i∈Ij bi and Aj = i∈Ij ai . Notice that in terms The Ising dataset used for classification by RFs here
of these quantities, we can write is identical to that used to study logistic regression in
T Sec. VII.C. We assign a label to each state according to
X 1 its phase: 0 if the state is disordered, and 1 if it is ordered.
∆Ct = [Bj wj + (Aj + λ)wj2 ] + γT, (114)
j=1
2 We divide the dataset into three categories according to
the temperature at which samples are drawn: ordered
where we made the t-dependence of all parameters im- (T /J < 2.0), near-critical (2.0 ≤ T /J ≤ 2.5) and disor-
plicit. Note that λ comes from the regularization term, dered (T /J > 2.5) (see Figure 20). We use the ordered
Ω(gt ), through Eq.(108). To find the optimal wj , just as and disordered states to train a random forest and eval-
in Newton’s method we take the gradient of the above uate our learned model on a test set of unseen ordered
expression with respect to wj and set this equal to zero, and disordered states (test sets). We also ask how well
to get our RF can predict the phase of samples drawn in the
Bj critical region (i.e. predict whether the temperature of
wjopt = − . (115) a critical sample is above or below the critical tempera-
Aj + λ
ture). Since our model is never trained on samples in the
Plugging this expression into ∆Ct gives critical region, prediction in this region is a test of the
T algorithm’s ability to generalize to new regions in phase
1 X Bj2 space.
∆Ctopt = − + γT. (116)
2 j=1 Aj + λ The results of fits using RFs to predict phases are
44
1.00
Train (coarse) Feature importance
0.95 Test (coarse)
Critical (coarse) axial MET 517
0.90 Train (fine) missing energy magnitude 505
MET_rel 461
Test (fine) dPhi_r_b 443
0.85 Critical (coarse) lepton 2 eta 416
Accuracy
Features
0.75 lepton 1 pT 351
lepton 2 pT 348
0.70 S_R 346
M_R 331
0.65 R 319
10 20 30 40 50 60 70 80 90 100 lepton 1 phi 314
Nestimators MT2 314
lepton 2 phi 290
M_TR_2 283
90 M_Delta_R 273
80 Fine 0 100 200 300 400 500
Coarse F score
70
60 FIG. 33 Feature Importance Scores in SUSY dataset from
applying XGBoost to 100, 000 samples. See Notebook 10 for
Run time (s)
50 more details.
40
30
20 physics-related datasets, we used the XGBoost imple-
mentation of gradient boosted trees to classify Monte-
10
Carlo collisions from the SUSY dataset. With default
0 parameters using a small subset of the data (100, 000 out
10 20 30 40 50 60 70 80 90 100
Nestimators of the full 5, 000, 000 samples), we were able to achieve
a classification accuracy of about 79%, which could be
FIG. 32 Using Random Forests (RFs) to classify Ising Phases. improved to nearly 80% after some fine-tuning (see ac-
(Top) Accuracy of RFs for classifying the phase of samples
companying notebook). This is comparable to published
from the Ising mode for the training set (blue), test set (red),
and critical region (green) using coarse trees with a few leaves results (Baldi et al., 2014) and those obtained using lo-
(triangles) and fine decision trees with many leaves (filled cir- gistic regression in earlier chapters. One nice feature of
cles). RFs were trained on samples from ordered and disor- ensemble methods such as XGBoost is that they auto-
dered phases but were not trained on samples from the critical matically allow us to calculate feature scores (Fscores)
region. (Bottom) The time it takes to train RFs scales lin- that rank the importance of various features for clas-
early with the number of estimators in the ensemble. For the sification. The higher the Fscore, the more important
upper panel, note that the train case (blue) overlaps with the
the feature for classification. Figure 33 shows the feature
test case (red). Here ‘fine’ and ‘coarse’ refer to trees with 2
and 10,000 leaves, respectively. For implementation details, scores from our XGBoost algorithm for the production of
see Jupyter notebooks 9 electrically-charged supersymmetric particles (χ±) which
decay to W bosons and an electrically neutral supersym-
metric particle χ0 , which is invisible to the detector. The
shown in Figure 32. We used two types of RF classifiers, features are a mix of eight directly measurable quanti-
one where the ensemble consists of coarse decision trees ties from the detector, as well as ten hand crafted fea-
with a few leaves and another with finer decision trees tures chosen using physics knowledge. Consistent with
with many leaves (see corresponding notebook). RFs the physics of these supersymmetric decays in the lepton
have extremely high accuracy on the training and test channel, we find that the most informative features for
sets (over 99%) for both coarse and fine trees. How- classification are the missing transverse energy along the
ever, notice that the RF consisting of coarse trees per- vector defined by the charged leptons (Axial MET) and
form extremely poorly on samples from the critical region the missing energy magnitude due to χ0 .
whereas the RF with fine trees classifies critical samples
with an accuracy of nearly 85%. Interestingly, and unlike
with logistic regression, this performance in the critical IX. AN INTRODUCTION TO FEED-FORWARD DEEP
region requires almost no parameter tuning. This is be- NEURAL NETWORKS (DNNS)
cause, as discussed above, RFs are largely immune to
overfitting problems even as the number of estimators in Over the last decade, neural networks have emerged
the ensemble becomes large. Increasing the number of es- as the one of most powerful and widely-used supervised
timators in the ensemble does increase performance but learning techniques. Deep Neural Networks (DNNs) have
at a large cost in computational time (Fig. 32 bottom). a long history (Bishop, 1995b; Schmidhuber, 2015), but
In the second application of ensemble methods to re-emerged to prominence after a rebranding as “Deep
45
Learning” in the mid 2000s (Hinton et al., 2006; Hinton and widely held beliefs not as universal as once imagined
and Salakhutdinov, 2006). DNNs truly caught the atten- (Lee et al., 2017; Zhang et al., 2016). This is especially
tion of the wider machine learning community and indus- true in modern neural networks where results are largely
try in 2012 when Alex Krizhevsky, Ilya Sutskever, and empirical and heuristic and lack the firm footing of many
Geoff Hinton used a GPU-based DNN model (AlexNet) earlier machine learning methods. For this reason, in
to lower the error rate on the ImageNet Large Scale this review we have chosen to emphasize tried and true
Visual Recognition Challenge (ILSVRC) by an incred- fundamentals, while pointing out what, from our current
ible twelve percent from 28% to 16% (Krizhevsky et al., vantage point, seem like promising new techniques. The
2012). Just three years later, a machine learning group field is rapidly evolving and readers are urged to read
from Microsoft achieved an error of 3.57% using an ultra- papers and to implement these algorithms themselves in
deep residual neural network (ResNet) with 152 layers order to gain a deep appreciation for the incredible power
(He et al., 2016)! Since then, DNNs have become the of modern neural networks, especially in the context of
workhorse technique for many image and speech recogni- image, speech, and natural language processing, as well
tion based machine learning tasks. The large-scale indus- as limitations of the current methods.
trial deployment of DNNs has given rise to a number of In physics, DNNs and CNNs have already found
high-level libraries and packages (Caffe, Keras, Pytorch, numerous applications. In statistical physics, they
and TensorFlow) that make it easy to quickly code and have been applied to detect phase transitions in 2D
deploy DNNs. Ising (Tanaka and Tomiya, 2017a) and Potts (Li
Conceptually, it is helpful to divide neural networks et al., 2017) models, lattice gauge theories (Wetzel and
into four categories: (i) general purpose neural net- Scherzer, 2017), and different phases of polymers (Wei
works for supervised learning, (ii) neural networks de- et al., 2017). It has also been shown that deep neural net-
signed specifically for image processing, the most promi- works can be used to learn free-energy landscapes (Sidky
nent example of this class being Convolutional Neural and Whitmer, 2017). At the same time, methods from
Networks (CNNs), (iii) neural networks for sequential statistical physics have been applied to the field of deep
data such as Recurrent Neural Networks (RNNs), and learning to study the thermodynamic efficiency of learn-
(iv) neural networks for unsupervised learning such as ing rules (Goldt and Seifert, 2017), to explore the hy-
Deep Boltzmann Machines. Here, we will limit our dis- pothesis space that DNNs span, make analogies between
cussions to the first two categories (unsupervised learn- training DNNs and spin glasses (Baity-Jesi et al., 2018;
ing is discussed later in the review). Though increas- Baldassi et al., 2017), and to characterize phase transi-
ingly important for many applications such as audio and tions with respect to network topology in terms of er-
speech recognition, for the sake of brevity, we omit a dis- rors (Li and Saad, 2017). In relativistic hydrodynam-
cussion of sequential data and RNNs from this review. ics, deep learning has been shown to capture features of
For an introduction to RNNs and LSTM networks see non-linear evolution and has the potential to accelerate
Chris Olah’s blog, [Link] numerical simulations (Huang et al., 2018), while in me-
08-Understanding-LSTMs/, and Chapter 13 of (Bishop, chanics CNNs have been used to predict eigenvalues of
2006) as well as the introduction to RNNs in Chapter 10 photonic crystals (Finol et al., 2018). Deep CNNs were
of (Goodfellow et al., 2016) for sequential data. used in lensing reconstruction of the cosmic microwave
Due to the number of recent books on deep learning background (Caldeira et al., 2018). Recently, DNNs have
(see for example Michael Nielsen’s introductory online been used to improve the efficiency of Monte-Carlo algo-
book (Nielsen, 2015) and the more advanced (Goodfellow rithms (Shen et al., 2018).
et al., 2016)), the goal of this section is to give a high- Deep learning has also found interesting applications
level introduction to the basic ideas behind of DNNs, in quantum physics. Various quantum phase transi-
and provide some practical knowledge for coding simple tions (Arai et al., 2017; Broecker et al., 2017; Iakovlev
neural nets for supervised learning tasks. This section et al., 2018; van Nieuwenburg et al., 2017b; Suchs-
assumes the reader is familiar with the basic concepts land and Wessel, 2018) can be detected and studied
introduced in earlier sections on logistic and linear re- using DNNs and CNNs, including the transverse-field
gression. Throughout, we strive to provide intuition be- Ising model (Ohtsuki and Ohtsuki, 2017), topological
hind the inner workings of DNNs, as well as highlight phases (Yoshioka et al., 2017; Zhang et al., 2017a,b)
limitations of present-day algorithms. and non-invasive topological quality control (Caio et al.,
The influx of corporate and industrial interests has 2019), and even non-equilibrium many-body localiza-
rapidly transformed the field in the last few years. This tion (van Nieuwenburg et al., 2017a,b; Schindler et al.,
massive influx of money and researchers has given rise to 2017; Venderley et al., 2017) and the characterization of
new dogmas and best practices that change rapidly. As photoexcited quantum states (Shinjo et al., 2019). DNNs
with most intellectual fields experiencing rapid expan- were recently applied in cold atoms to identify critical
sion, many commonly accepted heuristics many turn out points (Rem et al., 2018). Representing quantum states
not to be as powerful as thought (Wilson et al., 2017), as DNNs (Gao et al., 2017; Gao and Duan, 2017; Levine
46
{
timate rates of coherent and incoherent quantum pro- output
layer
cesses (Greplova et al., 2017), to obtain spectra of 1/f -
{
noise in spin-qubit devices (Zhang and Wang, 2018), and input
the recognition of state and charge configurations and layer
auto-tuning in quantum dots (Kalantre et al., 2017). In
quantum information theory, it has been shown that one
can perform gate decompositions with the help of neural
nets (Swaddle et al., 2017). In lattice quantum chromo- FIG. 34 Basic architecture of neural networks. (A)
dynamics, DNNs have been used to learn action param- The basic components of a neural network are stylized neu-
eters in regions of parameter space where principal com- rons consisting of a linear transformation that weights the
ponent analysis fails (Shanahan et al., 2018). Last but importance of various inputs, followed by a non-linear activa-
tion function. (b) Neurons are arranged into layers with the
not least, DNNs also found place in the study of quan- output of one layer serving as the input to the next layer.
tum control (Yang et al., 2017), and in scattering theory
to learn s-wave scattering length (Wu et al., 2018) of po-
tentials. a neuron-specific bias b(i) :
output with respect to the weights tends to zero since hidden layers (hence the “deep” in deep learning). There
∂σ/∂z → 0 for z 1. Such “vanishing gradients” are are many ideas of why such deep architectures are fa-
a feature of any saturating activation function (top row vorable for learning. Increasing the number of layers in-
of Fig. 35), making it harder to train deep networks. In creases the number of parameters and hence the represen-
contrast, for a non-saturating activation function such as tational power of neural networks. Indeed, recent numer-
ReLUs or ELUs, the gradients stay finite even for large ical experiments suggests that as long as the number of
inputs. parameters is larger than the number of data points, cer-
tain classes of neural networks can fit arbitrarily labeled
random noise samples (Zhang et al., 2016). This suggests
2. Layering neurons to build deep networks: network that large neural networks of the kind used in practice can
architecture. express highly complex functions. Adding hidden layers
is also thought to allow neural nets to learn more complex
The basic idea of all neural networks is to layer neurons features from the data. Work on convolutional networks
in a hierarchical fashion, the general structure of which is suggests that the first few layers of a neural network learn
known as the network architecture (see Fig. 34). In the simple, “low-level” features that are then combined into
simplest feed-forward networks, each neuron in the in- higher-level, more abstract features in the deeper layers.
put layer of the neurons takes the inputs x and produces Other works suggest that it is computationally and al-
an output ai (x) that depends on its current weights, see gorithmically easier to train deep networks rather than
Eq. (118). The outputs of the input layer are then treated shallow, wider nets, though this is still an area of major
as the inputs to the next hidden layer. This is usually controversy and active research (Mhaskar et al., 2016).
repeated several times until one reaches the top or output Choosing the exact network architecture for a neural
layer. The output layer is almost always a simple clas- network remains an art that requires extensive numer-
sifier of the form discussed in earlier sections: a logistic ical experimentation and intuition, and is often times
regression or soft-max function in the case of categorical problem-specific. Both the number of hidden layers and
data (i.e. discrete labels) or a linear regression layer in the number of neurons in each layer can affect the per-
the case of continuous outputs. Thus, the whole neural formance of a neural network. There seems to be no
network can be thought of as a complicated nonlinear single recipe for the right architecture for a neural net
transformation of the inputs x into an output ŷ that de- that works best. However, a general rule of thumb that
pends on the weights and biases of all the neurons in the seems to be emerging is that the number of parameters in
input, hidden, and output layers. the neural net should be large enough to prevent under-
The use of hidden layers greatly expands the represen- fitting (see theoretical discussion in (Advani and Saxe,
tational power of a neural net when compared with a sim- 2017)).
ple soft-max or linear regression network. Perhaps, the Empirically, the best architecture for a problem de-
most formal expression of the increased representational pends on the task, the amount and type of data that is
power of neural networks (also called the expressivity) is available, and the computational resources at one’s dis-
the universal approximation theorem which states that a posal. Certain architectures are easier to train, while
neural network with a single hidden layer can approxi- others might be better at capturing complicated depen-
mate any continuous, multi-input/multi-output function dencies in the data and learning relevant input features.
with arbitrary accuracy. The reader is strongly urged Finally, there have been numerous works that move be-
to read the beautiful graphical proof of the theorem in yond the simple deep, feed-forward neural network ar-
Chapter 4 of Nielsen’s free online book (Nielsen, 2015). chitectures discussed here. For example, modern neural
The basic idea behind the proof is that hidden neurons networks for image segmentation often incorporate “skip
allow neural networks to generate step functions with ar- connections” that skip layers of the neural network (He
bitrary offsets and heights. These can then be added et al., 2016). This allows information to directly propa-
together to approximate arbitrary functions. The proof gate to a hidden or output layer, bypassing intermediate
also makes clear that the more complicated a function, layers and often improving performance.
the more hidden units (and free parameters) are needed
to approximate it. Hence, the applicability of the ap-
proximation theorem to practical situations should not B. Training deep networks
be overemphasized. In condensed matter physics, a good
analogy are matrix product states, which can approxi- In the previous section, we introduced the basic ar-
mate any quantum many-body state to an arbitrary ac- chitecture for neural networks. Here we discuss how to
curacy, provided the bond dimension can be increased efficiently train large neural networks. Luckily, the basic
arbitrarily – a severe requirement not met in any useful procedure for training neural nets is the same as we used
practical implementation of the theory. for training simpler supervised learning algorithms, such
Modern neural networks generally contain multiple as logistic and linear regression: construct a cost/loss
48
0 0 0
-1 -1 -1
-5 0 5 -5 0 5 -5 0 5
ReLU Leaky ReLU ELU
6 6 6
4 4 4
2 2 2
0 0 0
-5 0 5 -5 0 5 -5 0 5
FIG. 35 Possible non-linear activation functions for neurons. In modern DNNs, it has become common to use non-linear
functions that do not saturate for large inputs (bottom row) rather than saturating functions (top row).
function and then use gradient descent to minimize the implement regularization (e.g. L1 or L2 regularizers), see
cost function and find the optimal weights and biases. Sec. VI.
Neural networks differ from these simpler supervised pro- For categorical data, the most commonly used loss
cedures in that generally they contain multiple hidden function is the cross-entropy (Eq. (76) and Eq. (81)),
layers that make taking the gradient computationally since the output layer is often taken to be a logistic clas-
more difficult. We will return to this in Sec. IX.D which sifier for binary data with two types of labels, or a soft-
discusses the “backpropagation” algorithm for computing max classifier if there are more than two types of labels.
gradients. The cross-entropy was already discussed extensively in
Like all supervised learning procedures, the first thing earlier sections on logistic regression and soft-max classi-
one must do to train a neural network is to specify a fiers, see Sec. VII. Recall that for classification of binary
loss function. Given a data point (xi , yi ), xi ∈ Rd+1 , data, the output of the top layer of the neural network is
the neural network makes a prediction ŷi (w), where w the probability ŷi (w) = p(yi = 1|xi ; w) that data point
are the parameters of the neural network. Recall that i is predicted to be in category 1. The cross-entropy be-
in most cases, the top output layer of our neural net tween the true labels yi ∈ {0, 1} and the predictions is
is either a continuous predictor or a classifier that makes given by
discrete (categorical) predictions. Depending on whether n
one wants to make continuous or categorical predictions,
X
E(w) = − yi log ŷi (w) + (1 − yi ) log [1 − ŷi (w)] .
one must utilize a different kind of loss function. i=1
For continuous data, the loss functions that are com- More generally, for categorical data, y can take on M
monly used to train neural networks are identical to those values so that y ∈ {0, 1, . . . , M − 1}. For each datapoint
used in linear regression, and include the mean squared i, define a vector yim called a ‘one-hot’ vector, such that
error (
1, if yi = m
1X
n yim = (121)
E(w) = (yi − ŷi (w))2 , (119) 0, otherwise.
n i=1
We can also define the probability that the neural net-
where n is the number of data points, and the mean- work assigns a datapoint to category m: ŷim (w) = p(yi =
absolute error (i.e. L1 norm) m|xi ; w). Then, the categorical cross-entropy is defined
as
1X n M −1
E(w) = |yi − ŷi (w)|. (120) X X
n i E(w) = − yim log ŷim (w)
i=1 m=0
The full cost function often includes additional terms that +(1 − yim ) log [1 − ŷim (w)] . (122)
49
51 [Link](
volutional layers. Every Dense() layer accepts as its first 52 optimizer=[Link](lr=0.01,
required argument an integer which specifies the number momentum=0.9),
of neurons. The type of activation function for the layer 53 loss=[Link].
is defined using the activation optional argument, the categorical_crossentropy,
input of which is the name of the activation function 54 metrics=[’accuracy’]
in string format. Examples include ’relu’, ’tanh’, 55 )
’elu’, ’sigmoid’, ’softmax’, see Fig. 35. In order Training the DNN is a one-liner using the fit()
for our DNN to work properly, we must ensure that the method of the Sequential class. The first two re-
numbers of output and input neurons for consecutive lay- quired arguments are the training input and output
ers match. Therefore, we specify the shape of the input data. As optional arguments, we specify the mini-
in the first layer of the model explicitly using the op- batch_size, the number of training epochs, and the test
tional argument input_shape=(N_features,), see line or validation_data. To monitor the training procedure
39 below. The sequential construction of the model then for every epoch, we set verbose=True.
allows Keras to infer the correct input/output dimensions 57 ##### train model using minibatches
of all hidden layers automatically. Hence, we only need 58 # train DNN
to specify the size of the softmax output layer to match 59 history=[Link](X_train, Y_train,
the number of categories, see line 47. 60 batch_size=64,
61 epochs=10,
35 ##### create deep neural network 62 validation_data=(X_test, Y_test),
36 # instantiate model 63 verbose=True
37 model = Sequential() 64 )
38 # add a dense all-to-all sigmoid layer
51
0.3
the activities of the output layer aL
j . It of course also indi-
rectly depends on all the activities of neurons in lower lay-
ers in the neural network through iteration of Eq. (123).
0.2 Let us define the error ∆L j of the j-th neuron in the L-th
layer as the change in cost function with respect to the
0.1 weighted input zjL
∂E
0 2 4 6 8 ∆L
j = . (125)
epoch ∂zjL
We now derive the final two backpropagation equations of the backpropagation algorithm. We can calculate all
using the chain rule. Since the error depends on neurons the derivatives using a single “forward” and “backward”
in layer l only through the activation of neurons in the pass of the neural network. This computational efficiency
subsequent layer l + 1, we can use the chain rule to write is crucial since we must calculate the gradient with re-
spect to all parameters of the neural net at each step
∂E X ∂E ∂z l+1 of gradient descent. These basic ideas also underly al-
∆lj = = k
∂zjl ∂zk
l+1 ∂z l
j most all modern automatic differentiation packages such
k
as Autograd (Pytorch).
X ∂zkl+1
= ∆l+1
k
k
∂zjl
! 2. Computing gradients in deep networks: what can go wrong
X with backprop?
= ∆l+1 l+1
k wkj σ 0
(zjl ). (III)
k
Armed with backpropagation and gradient descent, it
This is the third backpropagation equation. The final seems like it should be straightforward to train any neural
equation can be derived by differentiating of the cost network. However, until fairly recently it was widely be-
function with respect to the weight wjk
l
as lieved that training deep networks was an extremely dif-
ficult task. One reason for this was that even with back-
∂E ∂E ∂zjl propagation, gradient descent on large networks is ex-
l
= = ∆lj al−1
k (IV) tremely computationally expensive. However, the great
∂wjk ∂zjl ∂wjk
l
advances in computational hardware (and the widespread
Together, Eqs. (I), (II), (III), and (IV) define the four use of GPUs) has made this a much less vexing prob-
backpropagation equations relating the gradients of the lem than even a decade ago. It is hard to understate
activations of various neurons alj , the weighted inputs the impact these advances in computing have had on the
P l l−1 l practical utility of neural networks.
zjl = k wjk ak +bj , and the errors ∆lj . These equations
On a more technical and mathematical note, another
can be combined into a simple, computationally efficient
problem that occurs in deep networks, which transmit
algorithm to calculate the gradient with respect to all
information through many layers, is that gradients can
parameters (Nielsen, 2015).
vanish or explode. This is, appropriately, known as the
The Backpropagation Algorithm problem of vanishing or exploding gradients. This prob-
lem is especially pronounced in neural networks that try
1. Activation at input layer: calculate the activa- to capture long-range dependencies, such as Recurrent
tions a1j of all the neurons in the input layer. Neural Networks for sequential data. We can illustrate
this problem by considering a simple network with one
2. Feedforward: starting with the first layer, exploit
neuron in each layer. We further assume that all weights
the feed-forward architecture through Eq. (123) to
are equal, and denote them by w. The behavior of the
compute z l and al for each subsequent layer.
backpropagation equations for such a network can be in-
3. Error at top layer: calculate the error of the top ferred from repeatedly using Eq. (III):
layer using Eq. (I). This requires to know the ex- L−1 L−1
pression for the derivative of both the cost function Y Y
∆1j = ∆L
j wσ 0 (zj ) = ∆L
j (w)
L
σ 0 (zj ), (126)
E(w) = E(aL ) and the activation function σ(z).
j=0 j=0
We can now see where the name backpropagation and the gradient blow up. On the other hand, if wσ00 < 1
comes from. The algorithm consists of a forward pass the errors and gradients vanish. Only when the weights
from the bottom layer to the top layer where one calcu- satisfy wσ00 ≈ 1 and the neurons are not saturated will
lates the weighted inputs and activations of all the neu- the gradient stay well behaved for deep networks.
rons. One then backpropagates the error starting with This basic behavior holds true even in more compli-
the top layer down to the input layer and uses these errors cated networks. Rather than considering a single weight,
to calculate the desired gradients. This description makes we can ask about the eigenvalues (or singular values) of
clear the incredible utility and computational efficiency the weight matrices wjk l
. In order for the gradients to
53
be finite for deep networks, we need these eigenvalues to important to experiment with different variances.
stay near unity even after many gradient descent steps. The second important thing is to appropriately choose
In modern feedforward and ReLU neural networks, this the learning rate or step-size by searching over five log-
is achieved by initializing the weights for the gradient de- arithmic grid points (Wilson et al., 2017). If the best
scent in clever ways and using non-linearities that do not performance occurs at the edge of the grid, repeat this
saturate, such as ReLUs (recall that for saturating func- procedure until the optimal learning rate is in the middle
tions, σ 0 → 0, which will cause the gradient to vanish). of the grid parameters. Finally, it is common to center
Proper initialization and regularization schemes such as or whiten the input data (just as we did for linear and
gradient clipping (cutting-off gradients with very large logistic regression).
values), and batch normalization also help mitigate the Another important form of regularization that is of-
vanishing and exploding gradient problem. ten employed in practice is Early Stopping. The idea of
Early Stopping is to divide the training data into two por-
tions, the dataset we train on, and a smaller validation
E. Regularizing neural networks and other practical set that serves as a proxy for out-of-sample performance
considerations
on the test set. As we train the model, we plot both the
training error and the validation error. We expect the
DNNs, like all supervised learning algorithms, must
training error to continuously decrease during training.
navigate the bias-variance tradeoff. Regularization tech-
However, the validation error will eventually increase due
niques play an important role in ensuring that DNNs
to overfitting. The basic idea of early stopping is to halt
generalize well to new data. The last five years have seen
the training procedure when the validation error starts to
a wealth of new specialized regularization techniques for
rise. This Early Stopping procedure ensures that we stop
DNNs beyond the simple L1 and L2 penalties discussed in
the training and avoid fitting sample specific features in
the context of linear and logistic regression, see Secs. VI
the data. Early Stopping is a widely used essential tool
and VII. These new techniques include Dropout and
in the deep learning regularization toolbox.
Batch Normalization. In addition to these specialized
regularization techniques, large DNNs seem especially
well-suited to implicit regularization that already takes
2. Dropout
place in the Stochastic Gradient Descent (SGD) (Wilson
et al., 2017), cf. Sec. IV. The implicit stochasticity and
Another important regularization schemed that has
local nature of SGD often prevent overfitting of spurious
been widely adopted in the neural networks literature
correlations in the training data, especially when com-
is Dropout (Srivastava et al., 2014). The basic idea of
bined with techniques such as Early Stopping. In this
Dropout is to prevent overfitting by reducing spurious
section, we give a brief overview of these regularization
correlations between neurons within the network by in-
techniques.
troducing a randomization procedure similar to that un-
derlying ensemble models such as Bagging. Recall that
1. Implicit regularization using SGD: initialization, the basic idea behind ensemble methods is to train an en-
hyper-parameter tuning, and Early Stopping semble of models that are created using a randomization
procedure to ensure that the members of the ensemble
The most commonly employed and effective optimizer are uncorrelated, see Sec. VIII. This reduces the vari-
for training neural networks is SGD (see Sec. IV for other ance of statistical predictions without creating too much
alternatives). SGD acts as an implicit regularizer by in- additional bias.
troducing stochasticity (from the use of mini-batches) In the context of neural networks, it is extremely costly
that prevents overfitting. In order to achieve good per- to train an ensemble of networks, both from the point of
formance, it is important that the weight initialization is view of the amount of data needed, as well as computa-
chosen randomly, in order to break any leftover symme- tional resources and parameter tuning required. Dropout
tries. One common choice is drawing the weights from a circumnavigates these problems by randomly dropping
Gaussian centered around zero with some variance that out neurons (along with their connections) from the neu-
scales inversely with number of inputs to the neuron (He ral network during each step of the training (see Figure
et al., 2015; Sutskever et al., 2013). Since SGD is a lo- 39). Typically, for each mini-batch in the gradient de-
cal procedure, as networks get deeper, choosing a good scent step, a neuron is dropped from the neural network
weight initialization becomes increasingly important to with a probability p. The gradient descent step is then
ensure that the gradients are well behaved. Choosing performed only on the weights of the “thinned” network
an initialization with a variance that is too large or too of individual predictors.
small will cause gradients to vanish and the network to Since during training, on average weights are only
train poorly – even a factor of 2 can make a huge differ- present a fraction p of the time, predictions are made
ence in practice (He et al., 2015). For this reason, it is by reweighing the weights by p: wtest = pwtrain .The
54
Batch Normalization is a regularization scheme that Now that we have gained sufficient high-level back-
has been quickly adopted by the neural network com- ground knowledge about deep neural nets, let us discuss
munity since its introduction in 2015 (Ioffe and Szegedy, how to use them in practice.
2015). The basic inspiration behind Batch Normalization
is the long-known observation that training in neural net-
1. Deep learning packages
works works best when the inputs are centered around
zero with respect to the bias. The reason for this is that
In Sec. IX.C, we demonstrated that the numerical
it prevents neurons from saturating and gradients from
implementation of DNNs is greatly facilitated by open
vanishing in deep nets. In the absence of such center-
source python packages, such as Keras, TensorFlow, and
ing, changes in parameters in lower layers can give rise
Pytorch (and many more). The complexity and learning
to saturation effects in higher layers, and vanishing gra-
curves for these packages differ, depending on the user’s
dients. The idea of Batch Normalization is to introduce
level of familiarity with python. The reader should keep
additional new “BatchNorm” layers that standardize the
mind mind that there are DNN packages written in other
inputs by the mean and variance of the mini-batch.
languages, such as Caffe which uses C++, but we do not
Consider a layer l with d neurons whose inputs are
discuss them in this review for brevity.
(z1l , . . . , zdl ). We standardize each dimension so that
Keras is a high-level framework which does not require
any knowledge about the inner workings of the under-
zkl − E[zkl ]
zkl −→ ẑkl = q , (127) lying deep learning algorithms. Coding DNNs in Keras
Var[zkl ] is particularly simple, see Sec. IX.C, and allows one to
quickly grasp the big picture behind the theoretical con-
where the mean and variance are taken over all sam- cepts which we introduced above. However, for advanced
ples in the mini-batch. One problem with this procedure applications, which may require more direct control over
55
the operations in between the layers, Keras’ high-level 5. Evaluate and study the model performance
design may prove insufficient. on the test data.
If one opens up the Keras black box, one will find that
it wraps the functionality of another package – Tensor- 6. Use the validation data to adjust the hyper-
Flow12 . Over the last years, TensorFlow, which is sup- parameters (and, if necessary, network ar-
ported by Google, has been gaining popularity and has chitecture) to optimize performance for the
become the preferred library for deep learning. It is fre- specific dataset.
quently used in Kaggle competitions, university classes, At this point, a few remarks are in order. While we
and industry. In TensorFlow one constructs data flow treat Step 1 above as consisting mainly of loading and re-
graphs, the nodes of which represent mathematical oper- shaping a dataset prepared ahead of time, we emphasize
ations, while the edges encode multidimensional tensors that obtaining a sufficient amount of data is a typical
(data arrays). A deep neural net can then be thought challenge in many applications. Oftentimes insufficient
of as a graph with a particular architecture. One needs data serves as a major bottleneck on the ultimate per-
to understand this concept well before one can truly un- formance of DNNs. In such cases one can consider data
leash TensorFlow’s full potential. The learning curve can augmentation, i.e. distorting data samples from the ex-
sometimes be rather steep for TensorFlow beginners, and isting dataset in some way to enhance size the dataset.
requires a certain degree of perseverance and devoted Obviously, if one knows how to do this, one already has
time to internalize the underlying ideas. partial information about the important features in the
There are, however, many other open source packages data.
which allow for control over the inter- and intra-layer op- One of the first questions we are usually faced with is
erations, without the need to introduce computational how to determine the sizes of the training and test data
graphs. Such an example is Pytorch, which offers li- sets. The MNIST dataset, which has 10 classification
braries for automatic differentiation of tensors at GPU categories, uses 80% of the available data for training
speed. As we discussed above, manipulating neural nets and 20% for testing. On the other hand, the ImageNet
boils down to fast array multiplication and contraction data which has 1000 categories is split 50% − 50%. As
operations and, therefore, the [Link] library often a rule of thumb, the more classification categories there
does the job of providing enough access and controlla- are in the task, the closer the sizes of the training and
bility to manipulate the linear algebra operations under- test datasets should be in order to prevent overfitting.
lying deep neural nets. Once the size of the training set is fixed, it is common to
For the benefit of the reader, we have prepared Jupyter reserve 20% of it for validation, which is used to fine-tune
notebooks for DNNs using all three packages for the deep the hyperparameters of the model.
learning problems we discuss below. We invite the reader Also related to data preprocessing is the standardiza-
to carefully examine the differences in the code which tion of the dataset. It has been found empirically that if
should help them decide on which package they prefer to the original values of the data differ by orders of magni-
use. tude, training can be slowed down or impeded. This can
be traced back to the vanishing and exploding gradient
problem in backprop discussed in Sec. IX.D. To avoid
2. Approaching the learning problem
such unwanted effects, one often resorts to two tricks:
(i) all data should be mean-centered, i.e. from every data
Let us now analyze a typical procedure for using neu-
point we subtract the mean of the entire dataset, and (ii)
ral networks to solve supervised learning problems. As
rescale the data, for which there are two ways: if the data
can be seen already from the code snippets in Sec. IX.C,
is approximately normally distributed, one can rescale by
constructing a deep neural network to solve ML problems
the standard deviation. Otherwise, it is typically rescaled
is a multiple-stage process. Generally, one can identify a
by the maximum absolute value so the rescaled data lies
set of key steps:
within the interval [−1, 1]. Rescaling ensures that the
1. Collect and pre-process the data. weights of the DNN are of a similar order of magnitude.
The next issue is how to choose the right hyperpa-
2. Define the model and its architecture. rameters to begin training the model with. According
to Bengio, the optimal learning rate is often an order
3. Choose the cost function and the optimizer.
of magnitude lower than the smallest learning rate that
4. Train the model. blows up the loss (Bengio, 2012). One should also keep
in mind that, depending on how ambitious of a problem
one is dealing with, training the model can take a con-
siderable amount of time. This can severely slow down
12 While Keras can also be used with a Theano backend, we do not any progress on improving the model in Step 6. There-
discuss this here since Theano support has been discontinued. fore, it is usually a good idea to play with a small enough
56
accuracy (%)
data set size
10000 44.6% 65.0% 74.6% 77.0% 78.2% 60%
hyperparameter landscape is to use grid searches.
Whereas it is always possible to view Steps 1-5 as
generic and independent of the particular problem we 100000 42.7% 48.8% 74.1% 78.6% 79.6% 40%
are trying to solve, it is only when these steps are put
together in Step 6 that the real benefit of deep learning
is revealed, compared to less sophisticated methods such 200000 56.7% 59.9% 76.0% 78.9% 79.9%
20%
as regression or bagging, see Secs. VI, VII, and VIII. The
optimal choice of network architecture, cost function, and learning rate
optimizer is determined by the properties of the training 0%
and test datasets, which are only revealed when we try
to improve the model. FIG. 40 Grid search results for the test set accuracy of the
While there is no single recipe to approach all ML prob- DNN for the SUSY problem as a function of the learning rate
lems, we believe that the above list gives a good overview and the size of the dataset. The data used includes all high-
and low-level features.
and can be a useful guideline to the layman. Further-
more, as it becomes clear, this ‘recipe’ can be applied
to generic supervised learning tasks, not just DNNs. We cost function and minimize it using SGD with batches of
refer the reader to Sec. XI for more useful hints and tips size 10% of the training data size. We train the DNN for
on how to use the validation data during the training 10 epochs.
process. Figure 40 shows the accuracy of our DNN on the test
data as a function of the learning rate and the size of the
dataset. It is considered good practice to start with a
3. SUSY dataset logarithmic scale to search through the hyperparameters,
to get an overall idea for the order of magnitude of the
As a first example from physics, we discuss a DNN optimal values. In this example, the performance peaks
approach to the SUSY dataset already introduced in the at the largest size of the dataset and a learning rate of 0.1,
context of logistic regression in Sec. VII.C.2, and Bagging and is of the order of 80%. Since the optimal performance
in Sec. VIII.F. For a detailed description of the SUSY is obtained at the edge of the grid, we encourage the
dataset and the corresponding classification problem, we reader to extend the grid size to beat our result. For
refer the reader to Sec. VII.C.2. There is an interest in comparison, in the original study (Baldi et al., 2014),
using deep learning methods to automate the discovery the authors achieved ≈ 89% by using the entire dataset
of collision features from data. Benchmark results using with 5, 000, 000 points and a more sophisticated network
Bayesian Decision Trees from a standard physics pack- architecture, trained using GPUs.
age, and five-layer neural networks using Dropout were
presented in the original paper (Baldi et al., 2014); they
demonstrate the ability of deep learning to bypass the 4. Phases of the 2D Ising model
need of using hand-crafted high-level features. Our goal
here is to study systematically the accuracy of a DNN As a second example from physics, we discuss a DNN
classifier as a function of the learning rate and the dataset approach to the Ising dataset introduced in Sec. VII.C.1.
size. We study the problem of classifying the states of the 2D
Unlike the MNIST example where we used Keras, here Ising model with a DNN (Tanaka and Tomiya, 2017a),
we use the opportunity to introduce the Pytorch package, focussing on the model performance as a function of both
see the corresponding notebook. We leave the discussion the number of hidden neurons and the learning rate.
of the code-specific details for the accompanying note- The discussion is accompanied by a notebook written in
book. TensorFlow. As in the previous example, the interested
To classify the SUSY collision events, we construct a reader can find the discussion of the code-specific details
DNN with two dense hidden layers of 200 and 100 neu- in the notebook.
rons, respectively. We use ReLU activation between the To classify whether a given spin configuration is in the
input and the hidden layers, and a sofmax output layer. ordered or disordered phase, we construct a minimalistic
We apply dropout regularization on the weights of the model for a DNN with a single hidden layer containing
DNN. Similar to MNIST, we use the cross-entropy as a a number of hidden neurons. The network architecture
57
accuracy (%)
10 60.7% 43.1% 28.0% 75.6% 95.7% 99.2% 60%
nians, or work directly in momentum space if the system
is translationally invariant. This basic idea, tailoring our
100 68.9% 70.8% 85.1% 93.3% 99.7% 99.9%
40%
analysis to exploit additional structure, is a key feature of
modern physical theories from general relativity, through
1000 79.2% 53.9% 84.2% 97.9% 99.8% 99.9% gauge theories, to critical phenomena.
20% Like physical systems, many datasets and supervised
learning rate learning tasks also possess additional symmetries and
structure. For instance, consider a supervised learning
0% task where we want to label images from some dataset as
100%
being pictures of cats or not. Our statistical procedure
must first learn features associated with cats. Because
1e-06 1e-05 0.0001 0.001 0.01 0.1
a cat is a physical object, we know that these features
80% are likely to be local (groups of neighboring pixels in the
1 50.2% 50.1% 49.8% 49.2% 53.7% 60.1% two-dimensional image corresponding to whiskers, tails,
eyes, etc). We also know that the cat can be anywhere
hidden neurons
accuracy (%)
10 68.6% 46.4% 39.3% 74.5% 83.7% 84.2% 60% in the image. Thus, it does not really matter where in
the picture these features occur (though relative positions
of features likely do matter). This is a manifestation of
100 67.6% 63.9% 55.5% 67.1% 73.1% 93.7%
40% translational invariance that is built into our supervised
learning task. This example makes clear that, like many
1000 67.6% 51.9% 56.2% 59.9% 74.0% 96.0%
physical systems, many ML tasks (especially in the con-
20%
learning rate text of computer vision and image processing) also pos-
sess additional structure, such as locality and translation
0%
invariance.
The all-to-all coupled neural networks in the previous
FIG. 41 Grid search results for the test set accuracy (top) section fail to exploit this additional structure. For exam-
and the critical set accuracy (bottom) of the DNN for the ple, consider the image of the digit ‘four’ from the MNIST
Ising classification problem as a function of the learning rate dataset shown in Fig. 36. In the all-to-all coupled neural
and the number of hidden neurons. networks used there, the 28 × 28 image was considered
a one-dimensional vector of size 282 = 796. This clearly
throws away lots of the spatial information contained in
thus includes a ReLU-activated input layer, the hidden the image. Not surprisingly, the neural networks com-
layer, and the softmax output layer. We pick the cate- munity realized these problems and designed a class of
gorical cross-entropy as a cost function and minimize it neural network architectures, convolutional neural net-
using SGD with mini-batches of size 100. We train the works or CNNs, that take advantage of this additional
DNN for 100 epochs. structure (locality and translational invariance) (LeCun
Figure 41 shows the outcome of a grid search over a et al., 1995). Furthermore, what is especially interest-
log-spaced learning rate and the number of neurons in the ing from a physics perspective is the recent finding that
hidden layer. We see that about 10 neurons are enough these CNN architectures are intimately related to mod-
at a learning rate of 0.1 to get to a very high accuracy on els such as tensor networks (Stoudenmire, 2018; Stouden-
the test set. However, if we aim at capturing the physics mire and Schwab, 2016) and, in particular, MERA-like
close to criticality, clearly more neurons are required to architectures that are commonly used in physical mod-
reliably learn the more complex correlations in the Ising els for quantum condensed matter systems (Levine et al.,
configurations. 2017).
One of the core lessons of physics is that we should ex- A convolutional neural network is a translationally in-
ploit symmetries and invariances when analyzing physi- variant neural network that respects locality of the in-
58
H Fully
Connected
Layer
W
Convolution Coarse-graining Convolution Coarse-graining
(pooling) (pooling)
FIG. 42 Architecture of a Convolutional Neural Network (CNN). The neurons in a CNN are arranged in three
dimensions: height (H), width (W ), and depth (D). For the input layer, the depth corresponds to the number of channels (in
this case 3 for RGB images). Neurons in the convolutional layers calculate the convolution of the image with a local spatial
filter (e.g. 3 × 3 pixel grid, times 3 channels for first layer) at a given location in the spatial (W, H)-plane. The depth D of the
convolutional layer corresponds to the number of filters used in the convolutional layer. Neurons at the same depth correspond
to the same filter. Neurons in the convolutional layer mix inputs at different depths but preserve the spatial location. Pooling
layers perform a spatial coarse graining (pooling step) at each depth to give a smaller height and width while preserving the
depth. The convolutional and pooling layers are followed by a fully connected layer and classifier (not shown).
put data. CNNs are the backbone of many modern deep square filter of size F is a three-dimensional array of
learning applications and here we just give a high-level size F × F × Dl−1 . The convolution consists of running
overview of CNNs that should allow the reader to delve this filter over all locations in the spatial plane. To
directly into the specialized literature. The reader is also demonstrate how this works in practice, let us a consider
strongly encouraged to consult the excellent, succinct the simple example consisting of a one-dimensional
notes for the Stanford CS231n Convolutional Neural Net- input of depth 1, shown in Fig. 43. In this case, a filter
works class developed by Andrej Karpathy and Fei-Fei Li of size F × 1 × 1 can be specified by a vector of weights
([Link] We have drawn heavily w of length F . The stride, S, encodes by how many
on the pedagogical style of these notes in crafting this neurons we translate the filter by when performing
section. the convolution. In addition, it is common to pad
There are two kinds of basic layers that make up a the input with P zeros (see Fig. 43). For an input
CNN: a convolution layer that computes the convolu- of width W , the number of neurons (outputs) in the
tion of the input with a bank of filters (as a math- layer is given by (W − F + 2P )/S + 1. We invite the
ematical operation, see this practical guide to image reader to check out this visualization of the convolution
kernels: [Link] and procedure, [Link]
pooling layers that coarse-grain the input while main- arithmetic/blob/master/[Link], for a square
taining locality and spatial structure, see Fig. 42. For input of unit depth. After computing the filter, the
two-dimensional data, a layer l is characterized by three output is passed through a non-linearity, a ReLU in Fig.
numbers: height Hl , width Wl , and depth Dl 13 . The 43. In practice, one often inserts a BatchNorm layer
height and width correspond to the sizes of the two- before the non-linearity, cf. Sec. IX.E.3.
dimensional spatial (Wl , Hl )-plane (in neurons), and the These convolutional layers are interspersed with pool-
depth Dl (marked by the different colors in Fig. 42) – ing layers that coarse-grain spatial information by per-
to the number of filters in that layer. All neurons corre- forming a subsampling at each depth. One common pool-
sponding to a particular filter have the same parameters ing operation is the max pool. In a max pool, the spatial
(i.e. shared weights and bias). dimensions are coarse-grained by replacing a small region
In general, we will be concerned with local spatial (say 2×2 neurons) by a single neuron whose output is the
filters (often called a receptive field in analogy with maximum value of the output in the region. In physics,
neuroscience) that take as inputs a small spatial patch this pooling step is very similar to the decimation step
of the previous layer at all depths. For instance, a of RG (Iso et al., 2018; Koch-Janusz and Ringel, 2017;
Lin et al., 2017; Mehta and Schwab, 2014). This gener-
ally reduces the dimension of outputs. For example, if
13 The depth Dl is often called “number of channels”, to distin-
the region we pool over is 2 × 2, then both the height
guish it from the depth of the neural network itself, i.e. the total and the width of the output layer will be halved. Gen-
number of layers (which can be convolutional, pooling or fully- erally, pooling operations do not reduce the depth of the
connected), cf. Fig. 42. convolutional layers because pooling is performed sepa-
59
3 0 1 0 110
0 1 1 1 3 1
100
Accuracy
4 2
90
2 3 2 1
80
4 1 0 1 test
70 critical
. .
.. .. 60 1 5 10 20 50
Depth of hidden layer
2 5 2 1 FIG. 45 Single-layer convolutional network for classi-
fying phases in the Ising mode. Accuracy on test set and
critical samples for a convolutional neural network with sin-
1 1 0 1 5 2 gle layer of varying depth with filters of size 2, max-pool layer
with receptive field of size 2, followed by soft-max classifier.
Notice that the test accuracy is 100% even for a CNN of depth
1 0 1 0 3 1 one with a single set of weights. Accuracy on the near-critical
dataset is significantly below that for the test set.
3 1 0 1
age recognition based tasks, not just the ones they were
originally trained for. In other words, we expect that,
FIG. 44 Illustration of Max Pooling. Illustration of max- since images reflect the natural world, the filters learned
pooling over a 2 × 2 region. Notice that pooling is done at by these CNNs should transfer over to new tasks with
each depth (vertical axis) separately. The number of outputs only slight modifications and fine-tuning. In practice,
is halved along each dimension due to this coarse-graining. this turns out to be true for many tasks one might be
interested in.
There are three distinct ways one can take a pre-
Regarding the SUSY dataset, we stress that the absence trained CNN and repurpose it for a new task. The follow-
of spatial locality in the collision features renders apply- ing discussion draws heavily on the notes from the course
ing CNNs to that problem inadequate. CS231n mentioned in the introduction to this section.
• Use CNN as fixed feature detector at top
layer. If the new dataset we want to train on is
C. Pre-trained CNNs and transfer learning
small and similar to the original dataset, we can
simply use the CNN as a fixed feature detector
The immense success of CNNs for image recognition
and retrain our classifier. In other words, we re-
has resulted in the training of huge networks on enormous
move the classifier (soft-max) layer at the top of the
datasets, often by large industrial research teams from
CNN and replace it with a new classifier (linear sup-
Google, Microsoft, Amazon, etc. Many of these mod-
port vector machine (SVM) or soft-max) relevant
els are known by name: AlexNet, GoogLeNet, ResNet,
to our supervised learning problem. In this proce-
InceptionNet, VGGNet, etc. Most researchers and prac-
dure, the CNN serves as a fixed map from images
titioners do not have the resources, data, or time to train
to relevant features (the outputs of the top fully-
networks on this scale. Luckily, the trained models have
connected layer right before the original classifier).
been released and are now available in standard packages
This procedure prevents overfitting on small, simi-
such as the Torch Vision library in Pytorch or the Caffe
lar datasets and is often a useful starting point for
framework. These models can be used directly as a basis
transfer learning.
for fine-tuning in different supervised image recognition
tasks through a process called transfer learning. • Use CNN as fixed feature detector at inter-
The basic idea behind transfer learning is that the fil- mediate layer. If the dataset is small and quite
ters (receptive fields) learned by the convolution layers different from the dataset used to train the origi-
of these networks should be informative for most im- nal image, the features at the top level might not
61
Yes
Valida;on error high? More data
Overfi9ng Regulariza;on
1. Neural networks as representation learning
New model architecture
No
One important and powerful aspect of the deep learn-
ing paradigm is the ability to learn relevant features
DONE!
of the data with relatively little domain knowledge, i.e.
FIG. 46 Organizing a workflow for Deep Learning.
with minimal hand-crafting. Often, the power of deep
Schematic illustrating a deep learning workflow inspired by learning stems from its ability to act like a black box
navigating the bias-variance tradeoff (Figure based on An- that can take in a large stream of data and find good
drew Ng’s talk at the 2016 Deep Learning School available at features that capture properties of the data we are in-
[Link] In this terested in. This ability to learn good representations
diagram, we have assumed that there in no mismatch be- with very little hand-tuning is one of the most attractive
tween the distributions the training and test sets are drawn properties of DNNs. Many of the other supervised learn-
from.
ing algorithms discussed here (regression-based models,
ensemble methods such as random forests or gradient-
architecture. boosted trees) perform comparably or even better than
neural networks but when using hand-crafted features
If the validation and test sets are drawn from the same with small-to-intermediate sized datasets.
distributions, then good performance on the validation The hierarchical structure of deep learning models is
set should lead to similarly good performance on the thought to be crucial to their ability to represent com-
test set. (Of course performance will typically be slightly plex, abstract features. For example, consider the use
worse on the test set because the hyperparameters were of CNNs for image classification tasks. The analysis of
fit to the validation set.) However, sometimes the train- CNNs suggests that the lower-levels of the neural net-
ing data and test data differ in subtle ways because, for works learn elementary features, such as edge detectors,
example, they are collected using slightly different meth- which are then combined into higher levels of the net-
ods, or because it is cheaper to collect data in one way works into more abstract, higher-level features (e.g. the
versus another. In this case, there can be a mismatch famous example of a neuron that “learned to respond to
between the training and test data. This can lead to cats”) (Le, 2013). More recently, it has been shown that
the neural network overfitting these small differences be- CNNs can be thought of as performing tensor decompo-
tween the test and training sets, and a poor performance sitions on the data similar to those commonly used in
on the test set despite having a good performance on numerical methods in modern quantum condensed mat-
the validation set. To rectify this, Andrew Ng suggests ter (Cohen et al., 2016).
making two validation or dev sets, one constructed from
the training data and one constructed from the test data. One of the interesting consequences of this line of
The difference between the performance of the algorithm thinking is the idea that one can train a CNN on one
on these two validation sets quantifies the train-test mis- large dataset and the features it learns should also be use-
match. This can serve as another important diagnostic ful for other supervised tasks. This results in the ability
when using DNNs for supervised learning. to learn important and salient features directly from the
data and then transfer this knowledge to a new task. In-
deed, this ability to learn important, higher-level, coarse-
B. Why neural networks are so successful: three high-level grained features is reminiscent of ideas like the renormal-
perspectives on neural networks ization group (RG) in physics where the RG flows sep-
arate out relevant and irrelevant directions, and certain
Having discussed the basics of neural networks, we con- unsupervised deep learning architectures have a natural
clude by giving three complementary perspectives on the interpretation in terms of variational RG schemes (Mehta
success of DNNs and Deep Learning. This high-level dis- and Schwab, 2014).
63
Performance
world. One important reason for the success of DNNs Medium NN
is that they are able to exploit the additional signal in
large datasets for difficult supervised learning tasks. Fun- Small NN
damentally, modern DNNs are unique in that they con-
tain millions of parameters, yet can still be trained on
existing hardwares. The complexity of DNNs (in terms Traditional
of parameters) combined with their simple architecture (e.g. logistic reg)
(layer-wise connections) hit a sweet spot between expres-
sivity (ability to represent very complicated functions)
and trainability (ability to learn millions of parameters).
Indeed, the ability of large DNNs to exploit huge
datasets is thought to differ from many other commonly
Amount of Data
employed supervised learning methods such as Support
FIG. 47 Large neural networks can exploit the vast
Vector Machines (SVMs). Figure 47 shows a schematic amount of data now available. Schematic of how neural
depicting the expected performance of DNNs of differ- network performance depends on amount of available data
ent sizes with the number of data samples and compares (Figure based on Andrew Ng’s talk at the 2016 Deep Learn-
them to supervised learning algorithms such as SVMs or ing School available at [Link]
ensemble methods. When the amount of data is small, F1ka6a13S9I.)
DNNs offer no substantial benefit over these other meth-
ods and often perform worse. However, large DNNs seem
to be able to exploit additional data in a way other meth- train much larger DNNs than was thought possible fifteen
ods cannot. The fact that one does not have to hand years ago. Furthermore, many of these computational
engineer features makes the DNN even more well suited gains are quickly incorporated into modern packages with
for handling large datasets. Recent theoretical results industrial resources. This makes it easy to perform nu-
suggest that as long as a DNN is large enough, it should merical experiments on large datasets, leading to further
generalize well and not overfit (Advani and Saxe, 2017). engineering gains.
In the data-rich world we live in (at least in the context
of images, videos, and natural language), this is a recipe
for success. In other areas where data is more limited, C. Limitations of supervised learning with deep networks
deep learning architectures have (at least so far) been less
successful. Like all statistical methods, supervised learning using
neural networks has important limitations. This is es-
pecially important when one seeks to apply these meth-
ods, especially to physics problems. Like all tools, DNNs
3. Neural networks scale up well computationally
are not a universal solution. Often, the same or better
performance on a task can be achieved by using a few
A final feature that is thought to underlie the success
hand-engineered features (or even a collection of random
of modern neural networks is that they can harness the
features). This is especially important for hard physics
immense increase in computational capability that has
problems where data (or Monte-Carlo samples) may be
occurred over the last few decades. The architecture of
hard to come by.
neural networks naturally lends itself to parallelization
Here we list some of the important limitations of su-
and the exploitation of fast but specialized processors
pervised neural network based models.
such as graphical processing units (GPUs). Google and
NVIDIA set on a course to develop TPUs (tensor pro- • Need labeled data.—Like all supervised learn-
cessing units) which will be specifically designed for the ing methods, DNNs for supervised learning require
mathematical operations underlying deep learning archi- labeled data. Often, labeled data is harder to ac-
tectures. The layered architecture of neural networks also quire than unlabeled data (e.g. one must pay for
makes it easy to use modern techniques such as automatic human experts to label images).
differentiation that make it easy to quickly deploy them.
Algorithms such as stochastic gradient descent and the • Supervised neural networks are extremely
use of mini-batches make it easy to parallelize code and data intensive.—DNNs are data hungry. They
64
perform best when data is plentiful. This is doubly sional reduction lies in the fact that low-dimensional rep-
so for supervised methods where the data must also resentations of high-dimensional data necessarily incurs
be labeled. The utility of DNNs is extremely lim- information lost. Below, we introduce some common lin-
ited if data is hard to acquire or the datasets are ear and nonlinear methods for performing dimensional
small (hundreds to a few thousand samples). In reduction with applications in data visualization of high-
this case, the performance of other methods that dimensional data.
utilize hand-engineered features can exceed that of
DNNs.
A. Some of the challenges of high-dimensional data
• Homogeneous data.—Almost all DNNs deal
with homogeneous data of one type. It is very hard Before we begin exploring some specific dimensional
to design architectures that mix and match data reduction techniques, it is useful to highlight some of the
types (i.e. some continuous variables, some discrete generic difficulties encountered when dealing with high-
variables, some time series). In applications beyond dimensional data.
images, video, and language, this is often what is
required. In contrast, ensemble models like random
forests or gradient-boosted trees have no difficulty a. High-dimensional data lives near the edge of sample space.
handling mixed data types. Geometry in high-dimensional space can be counterintu-
itive. One example that is pertinent to machine learning
• Many physics problems are not about
is the following. Consider data distributed uniformly at
prediction.—In physics, we are often not inter-
random in a D-dimensional hypercube C = [−e/2, e/2]D ,
ested in solving prediction tasks such as classifi-
where e is the edge length. Consider also a D-dimensional
cation. Instead, we want to learn something about
hypersphere S of radius e/2 centered at the origin and
the underlying distribution that generates the data.
contained within C. The probability that a data point x
In this case, it is often difficult to cast these ideas in
drawn uniformly at random in C is contained within S
a supervised learning setting. While the problems
is well approximated by the ratio of the volume of S to
are related, it’s possible to make good predictions
that of C : p(kxk2 < e/2) ∼ (1/2)D . Thus, as the di-
with a “wrong” model. The model might or might
mension of the feature space D increases, p goes to zero
not be useful for understanding the physics.
exponentially fast. In other words, most of the data will
Some of these remarks are particular to DNNs, oth- concentrate outside the hypersphere, in the corners of
ers are shared by all supervised learning methods. This the hypercube. In physics, this basic observation under-
motivates the use of unsupervised methods which in part lies many properties of ideal gases such as the Maxwell
circumnavigate these problems. distribution and the equipartition theorem (see Chapter
3 of (Sethna, 2006) for instance).
a) b)
2D 1D
FIG. 48 The “Swiss roll”. Data distributed in a three- FIG. 49 Illustration of the crowding problem. (Left) A two-
dimensional space (a) that can effectively be described on a dimensional dataset X consisting of 3 equidistant points.
two-dimensional surface (b). A common goal of dimensional (Right) Mapping X to a one-dimensional space while try-
reduction techniques is to preserve ordination in the data: ing to preserve relative distances leads to a collapse of the
points that are close-by in the original space are also near-by mapped data points.
in the mapped (latent) space. This is true of the mapping (a)
to (b) as can be seen by inspecting the color gradient.
signal
noise
We are interested in finding a new basis for the data and Cox, 2000). Moving forward, we use the term “dis-
that emphasizes highly variable directions while reducing tance” and “dissimilarity” interchangeably. There are two
redundancy between basis vectors. In particular, we will types of MDS: metric and non-metric. In metric MDS,
look for a linear transformation that reduces the covari- the distance is computed under a pre-defined metric and
ance between different features. To do so, we first per- the latent coordinates Ỹ are obtained by minimizing the
form singular value decomposition (SVD) on the design difference between the distance measured in the original
matrix X, namely, X = U SV T , where S is a diagonal space (dij (X)) and that in the latent space (dij (Y )):
matrix of singular value si , the orthogonal matrix U con- X
tains (as its columns) the left singular vectors of X, and Ỹ = arg min wij |dij (X) − dij (Y )|, (131)
Y
similarly V contains (as its columns) the right singular i<j
vectors of X. With this, one can rewrite the covariance
where wij ≥ 0 are weight values. The weight matrix wij
matrix as
is a set of free parameters that specify the level of confi-
1 dence (or precision) in the value of dij (X). If Euclidean
Σ(X) = V SU T U SV T
N −1 metric is used, MDS gives the same result as PCA and is
S2
usually referred to as classical scaling (Torgerson, 1958).
=V VT Thus MDS is often considered as a generalization of PCA.
N −1
In non-metric MDS, dij can be any distance matrix. The
≡ V ΛV T . (130) objective function is then to preserve the ordination in
the data, i.e. if d12 (X) < d13 (X) in the original space,
where Λ is a diagonal matrix with eigenvalues λi in the
then in the latent space we should have d12 (Y ) < d13 (Y ).
decreasing order along the diagonal (i.e. eigendecompo-
Both MDS and PCA can be implemented using stan-
sition). It is clear that the right singular vectors of X
dard Python packages such as Scikit. MDS algorithms
(i.e. the columns of V ) are principal directions of Σ(X),
typically have a scaling of O(N 3 ) where N corresponds
and the singular values of X are related to the eigenval-
to the number of data points, and are thus very limited
ues of the covariance matrix Σ(X) via λi = s2i /(N − 1).
in their application to large datasets. However, sample-
To reduce the dimensionality of data from p to p̃ < p, we
based methods have been introduce to reduce this scaling
first construct the p× p̃ projection matrix Ṽp0 by selecting
to O(N log N ) (Yang et al., 2006). In the case of PCA,
the singular components with the p̃ largest singular val-
a complete decomposition has a scaling of O(N p2 + p3 ),
ues. The projection of the data from p to a p̃ dimensional
where p is the number of features. Note that the first
space is simply Ỹ = X Ṽp0 . The same idea is central
term N p2 is due to the computation of covariance ma-
to matrix-product-state-like techniques used to compress
trix Eq.(129) while the second, p3 , stems from eigenvalue
the number of components in quantum wavefunctions in
decomposition. Nothe that PCA can be improved to bear
studies of low-dimensional many-body lattice systems.
complexity O(N p2 +p) if only the first few principal com-
The singular vector with the largest singular value (i.e
ponents are desired (using iterative approaches). PCA
the largest variance) is referred to as the first principal
and MDS are often among the first data visualization
component; the singular vector with the second largest
techniques one resorts to.
singular value as the second principal component,
Pp and so
on. An important quantity is the ratio λi / i=1 λi which
is referred as the percentage of the explained variance D. t-SNE
contained in a principal component (see FIG. 51.b).
It is common in data visualization to present the data It is often desirable to preserve local structures in
projected on the first few principal components. This is high-dimensional datasets. However, when dealing with
valid as long as a large part of the variance is explained datasets having clusters delimitated by complicated sur-
in those components. Low values of explained variance faces or datasets with a large number of clusters, preserv-
may imply that the intrinsic dimensionality of the data is ing local structures becomes difficult using linear tech-
high or simply that it cannot be captured by a linear rep- niques such as PCA. Many non-linear techniques such as
resentation. For a detailed introduction to PCA, see the non-classical MDS (Cox and Cox, 2000), self-organizing
tutorials by Shlens (Shlens, 2014) and Bishop (Bishop, map (Kohonen, 1998), Isomap (Tenenbaum et al., 2000)
2006). and Locally Linear Embedding (Roweis and Saul, 2000)
have been proposed and to address this class of problems.
These techniques are generally good at preserving local
C. Multidimensional scaling structures in the data but typically fail to capture struc-
tures at the larger scale such as the clusters in which the
Multidimensional scaling (MDS) is a non-linear dimen- data is organized (Maaten and Hinton, 2008).
sional reduction technique which preserves the pairwise Recently, t-stochastic neighbor embedding (t-SNE) has
distance or dissimilarity dij between data points (Cox emerged as one of the go-to methods for visualizing high-
68
dimensional data. It has been shown to offer insight- t-SNE constructs a similar probability distribution qij
ful visualization for many benchmark high-dimensional in a low dimensional latent space (with coordinates Y =
0
datasets (Maaten and Hinton, 2008). t-SNE is a non- {yi }, yi ∈ Rp , where p0 < p is the dimension of the latent
parametric16 method that constructs non-linear embed- space):
dings. Each high-dimensional training point is mapped
to low-dimensional embedding coordinates, which are op- (1 + ||yi − yj ||2 )−1
qij = P 2 −1
. (134)
timized in a way to preserve the local structure in the k6=i (1 + ||yi − yk || )
data.
When used appropriately, t-SNE is a powerful tech- The crucial point to note is that qij is chosen to be a long
nique for unraveling the hidden structure of high- tail distribution. This preserves short distance informa-
dimensional datasets while at the same time preserv- tion (relative neighborhoods) while strongly repelling two
ing locality. In physics, t-SNE has recently been used points that are far apart in the original space (see FIG.
to reduce the dimensionality and classify spin configu- 52). In order to find the latent space coordinates yi , t-
rations, generated with the help of Monte Carlo sim- SNE minimizes the Kullback-Leibler divergence between
ulations, for the Ising (Carrasquilla and Melko, 2017) qij and pij :
and Fermi-Hubbard models at finite temperatures (Ch’ng X
pij
et al., 2017). It was also applied to study clustering tran- C(Y ) = DKL (p||q) ≡ pij log . (135)
qij
sitions in glass-like problems in the context of quantum ij
short-tail
yi = arg min |p(xi ) q(y)|
y
long-tail
p
q
x0,y0
y1 x2 y2
x1
FIG. 52 Illustration of the t-SNE embedding. xi points cor-
respond to the original high-dimensional points while the yi
points are the corresponding low-dimensional map points pro-
duced by t-SNE. Here we consider two points, x1 , x2 , that are
respectively “close” and “far” from x0 . The high-dimensional
Gaussian (short-tail) distribution p(x) of x0 ’s neighbors is
shown in blue. The low-dimensional Cauchy (fat-tail) distri-
bution q(y) of x0 ’s neighbors is shown in red. The map point FIG. 53 Different visualizations of a Gaussian mixture formed
yi , are obtained by minimizing the difference |q(y) − p(xi )| of K = 30 mixtures in a D = 40 dimensional space. The
(similar to minimizing the KL divergence). We see that point Gaussians have the same covariance but have means drawn
x1 is mapped to short distances |y1 −y0 |. In contrast, far-away uniformly at random in the space [−10, 10]40 . (a) Plot of the
points such as x2 are mapped to relatively large distances first two coordinates. The labels of the different Gaussian is
|y2 − y0 |. indicated by the different colors. Note that in a realistic set-
ting, label information is of course not available, thus making
it very hard to distinguish the different clusters. (b) Random
• Scales are deformed in t-SNE. Since a scale-free dis- projection of the data onto a 2 dimensional space. (c) pro-
tribution is used in the latent space, one should not jection onto the first 2 principal components. Only a small
put too much emphasis on the meaning of the size fraction of the variance is explained by those components (the
of any clusters observed in the latent space. ratio is indicated along the axis). (d) t-SNE embedding (per-
plexity = 60, # iteration = 1000) in a 2 dimensional latent
• t-SNE is computationally intensive. Finally, a space. t-SNE captures correctly the local structure of the
data.
direct implementation of t-SNE has an algorith-
mic complexity of O(N 2 ) which is only appli-
cable to small to medium data sets. Improved t-SNE keeps nearby points close together while repelling
scaling of the form O(N log N ) can be achieved points that are far apart.
at the cost of approximating Eq. (135) by us-
Figure 54 shows t-SNE and PCA plots for the MNIST
ing the Barnes-Hut method (Van Der Maaten,
dataset of ten handwritten numerical digits (0-9). It is
2014) for N -body simulations (Barnes and Hut,
clear that the non-linear nature of t-SNE makes it much
1986). More recently extremely efficient t-
better at capturing and visualizing the complicated cor-
SNE implementation making use of fast Fourier
relations between digits, compared to PCA.
transforms for kernel summations in (136) have
been made available on [Link]
KlugerLab/FIt-SNE (Linderman et al., 2017).
XIII. CLUSTERING
As an illustration, in Figure 53 we applied t-SNE to
a Gaussian mixture model consisting of thirty Gaus- In this section, we continue our discussion of unsuper-
sians, whose means are uniformly distributed in forty- vised learning methods. Unsupervised learning is con-
dimensional space. We compared the results to a random cerned with discovering structure in unlabeled data (for
two-dimensional projection and PCA. It is clear that un- instance learning local structures for data visualization,
like more naïve dimensional reduction techniques, both see section XII). The lack of labels make unsupervised
PCA and t-SNE can identify the presence of well-formed learning much more difficult and subtle than its super-
clusters. The t-SNE visualization cleanly separates all vised counterpart. What is somewhat surprising is that
the clusters while certain clusters blend together in the even without labels it is still possible to uncover and ex-
PCA plot. This is a direct consequence of the fact that ploit the hidden structure in the data. Perhaps, the sim-
70
1 X 10.0
µk = rnk xn . (138)
Nk n 9.5
9.0
2. Maximization: Given a set of cluster means {µk },
find the assignments {rnk } which minimizes C. 8.5
Clearly, this is achieved by assigning each data C
8.0
point to their nearest cluster-mean: 0 10 20
( t
1 if k = arg mink0 (xn − µk0 )2
rnk = (139) FIG. 55 K-means with K = 3 applied to an artificial two-
0 otherwise
dimensional dataset. The cluster means at each iteration are
indicated by cyan star markers. t indicates the iteration num-
K-means clustering consists in alternating between these ber and C the value of the objective function. (a) The algo-
two steps until some convergence criterion is met. Practi- rithm is initialized by randomly partitioning the space into 3
cally, the algorithm should terminate when the change in sectors to generate an initial assignment. (b)-(c) For well sep-
the objective function from one iteration to another be- arated clusters, the algorithm converges rapidly to the true
clusters. (d) The objective function as a function of the it-
comes smaller than a pre-specified threshold. A simple
eration. C converges after t = 18 iterations for this choice of
example of K-means is presented in FIG. 55. random seed (for center initialization).
A nice property of the K-means algorithm is that it
is guaranteed to converge. To see this, one can verify
explicitly (by taking second-order derivatives) that the sively merged to form larger clusters. The merging pro-
expectation step always decreases C. This is also true for cess generates a hierarchy of clusters that can be visu-
the assignment step. Thus, since C is bounded from be- alized in the form of a dendrogram (see FIG. 56). This
low, the two-step iteration of K-means always converges hierarchy can be useful to analyze the relation between
to a local minimum of C. Since C is generally a non- clusters and the subcomponents of individual clusters.
convex function, in practice one usually needs to run the Agglomerative methods are usually specified by defin-
algorithm with different initial random cluster center ini- ing a distance measure between clusters18 . We denote
tializations and post-select the best local minimum. A the distance between clusters X and Y by d(X, Y ) ∈ R.
simple implementation of K-means has an average com- Different choices of distance result in different clustering
putational complexity which scales linearly in the size of algorithms. At each step, the two clusters that are the
the data set (more specifically the complexity is O(KN ) closest with respect to the distance measure are merged
per iteration) and is thus scalable to very large datasets. until a single cluster is left.
As we will see in section XIII.B, K-means is a hard-
assignment limit of the Gaussian mixture model where
all cluster variances are assumed to be the same. This Agglomerative clustering algorithm Agglomerative
highlights a common drawback of K-means: if the true clustering algorithms can thus be summarized as follows:
clusters have very different variances (spreads), K-means
can lead to spurious results since the underlying assump- 1. Initialize each point to its own cluster.
tion is that the latent model has uniform variances.
2. Given a set of K clusters X1 , X2 , · · · , XK , merge
clusters until one cluster is left (K = 1):
2. Hierarchical clustering: Agglomerative methods
d(X, Y )
xi ∈Xi ,xj ∈Xj
0.15
3. Average linkage: average distance between points
of different clusters 0.10
1 X 0.05
d(Xi , Xj ) = ||xi − xj ||2 (142)
|Xi | · |Xj |
xi ∈Xi ,xj ∈Xj
0.00
0 1 4 5 2 3
4. Ward’s linkage: This distance measure is analogous Leaf label
to the K-means method as it seeks to minimize the FIG. 56 Hierarchical clustering example with single linkage.
total inertia. The distance measure is the “error (a) The data points are successively grouped as denoted by the
squared” before and after merging which simplifies colored dotted lines. (b) Dendrogram representation of the
to: hierarchical decomposition. Each node of the tree represents
a cluster. One has to specify a scale cut-off for the distance
|Xi ||Xj | measure d(X, Y ) (corresponding to a horizontal cut in the
d(Xi , Xj ) = (µi − µj )2 , (143)
|Xi ∪ Xj | dendrogram) in order to obtain a set of clusters.
Relative density
expect. Finally, a point xi is said to be density-reachable
if it is in the ε-neighborhood of a core-point. From these
definitions, the algorithm can be simply formulated (see
also Figure 57):
(cluster identity) are a way of representing and abstract- variables”). For instance if we were considering a Gaus-
ing the correlations between datapoints. sian mixture with K = 3, we would have three possible
In this language, we can think of clustering as an algo- values for z ≡ (z1 , z2 , z3 ) : (1, 0, 0), (0, 1, 0) and (0, 0, 1).
rithm to learn the most probable value of a latent variable We cannot directly observe the variable z. It is a latent
(cluster identity) associated with each datapoint. Calcu- variable that encodes the cluster identity of point x. Let
lating this latent variable requires additional assumptions us also denote all the N latent variables corresponding
about the structure of our dataset. Like all unsuper- to a dataset X by Z.
vised learning algorithms, in clustering we must make Viewing the GMM as a generative model, we can write
an assumption about the underlying probability distri- the probability p(x|z) of observing a data point x given
bution from which the data was generated. Our model z as
for how the data is generated is called the generative
K
model. In clustering, we assume that data points are as- Y
p(x|z; {µk , Σk }) = N (x|µk , Σk )zk (148)
signed a cluster, with each cluster characterized by some
k=1
cluster-specific probability distribution (e.g. a Gaussian
with some mean and variance that characterizes the clus- as well as the probability of observing a given value of
ter). We then specify a procedure for finding the value latent variable
of the latent variable. This is often done by choosing
K
the values of the latent variable that minimize some cost Y
p(z|{πk }) = πkzk . (149)
function.
k=1
One common choice for a class of cost functions for
many unsupervised learning problems is Maximum Like- Using Bayes’ rule, we can write the joint probability of
lihood Estimation (MLE), see Secs. V and VI. In MLE, a clustering assignment z and a data point x given the
we choose the values of the latent variables that maximize GMM parameters as
the likelihood of the observed data under our generative
model (i.e. maximize the probability of getting the ob- p(x, z; θ) = p(x|z; {µk , Σk })p(z|{πk }). (150)
served dataset under our generative model). Such MLE
equations often give rise to the kind of Expectation Maxi- We can also use Bayes rule to rearrange this expression
mization (EM) equations that we first encountered in the to give the conditional probability of the data point x
last section in the context of K-means clustering. being in the k-th cluster, γ(zk ), given model parameters
Gaussian Mixtures models (GMM) are a generative θ as
model often used in the context of clustering. In GMM, πk N (x|µk , Σk )
points are drawn from one of K Gaussians, each with its γ(zk ) ≡ p(zk = 1|x; θ) = PK . (151)
own mean µk and covariance matrix Σk , j=1 πj N (x|µj , Σj )
For future reference, let us denote the set of parameters where θ = {µk , Σk , πk }. Once we know the MLEs θ̂, we
(of K Gaussians in the model) {µk , Σk , πk } by θ. could use Eq. (151) to calculate the optimal hard cluster
To see how we can use GMM and MLE to perform assignment arg maxk γ̂(zk ) where γ̂(zk ) = p(zk = 1|x; θ̂).
clustering, we introduce discrete binary K-dimensional In practice, due to the complexity of Eq. (147), it is
latent variables z for each data point x whose k-th com- almost impossible to find the global maximum of the like-
ponent is 1 if point x was generated from the k-th Gaus- lihood function. Instead, we must settle for a local max-
sian and zero otherwise (these are often called “one-hot imum. One approach to finding a local maximum of the
75
likelihood is to use a method like stochastic gradient de- converges to a local maximum of the likelihood (Neal and
scent on the negative log-likelihood, cf. Sec IV. Here, we Hinton, 1998).
introduce an alternative, powerful approach for finding The central observation underlying EM is that it is
local minima in latent variable models using an iterative often much easier to calculate the conditional likelihoods
procedure called Expectation Maximization (EM). Given of the latent variables p̃(t) (Z) = p(Z|X; θ(t) ) given some
an initial guess for the parameters θ(0) , the EM algorithm choice of parameters, and the maximum of the expected
iteratively generates new estimates for the parameters log likelihood given an assignment of the latent variables:
θ(1) , θ(2) , . . .. Importantly, the likelihood is guaranteed to θ(t+1) = arg maxθ Ep(Z|X;θ(t) ) [log p(X, Z; θ)]. To get an
be non-decreasing under these iterations and hence EM intuition for this later quantity notice that we can write
N X
K
(t)
X
Ep̃(t) [log p(X, Z; θ)] = γik [log N (xi |µk , Σk ) + log πk ] , (153)
i=1 k=1
(t)
where we have used the shorthand γik = p(zik |X; θ(t) ) random measurement error for each feature. This in turn
with zik the k-th component of zi . Taking the derivative leads to increased errors for pairwise similarity and dis-
of this equation with
P respect to µk , Σk , and πk (subject tance measures and thus tends to “blur” distances be-
to the constraint k πk = 1) and setting this to zero tween data points (Domingos, 2012; Kriegel et al., 2009;
yields the intuitive equations Zimek et al., 2012). Many clustering algorithms rely on
PN (t) the explicit use of a similarity measure or distance met-
(t+1) i γik xi rics that weigh all features equally. For this reason, one
µk = P (t)
i γik must be careful when using an off-the-shelf method in
PN (t) high dimensions.
(t+1) γ (xi − µk )(xi − µk )T
Σk = i ik P (t) In order to perform clustering on high-dimensional
i γik data, it is often useful to denoise the data before pro-
(t+1) 1 X (t) ceeding with using a standard clustering method such as
πk = γik (154)
N
k
K-means (Kriegel et al., 2009). Figure 54 illustrates an
application of denoising to high-dimensional data. PCA
These are just the usual estimates for the mean and vari- (section XII.B) was used to denoise the MNIST dataset
ance, with each data point weighed according to our cur- by projecting the 784 original dimensions onto the 40 di-
rent best guess for the probability that it belongs to clus- mensions with the largest principal components. The re-
ter k. We can then use our new estimate θ(t+1) to calcu- sulting features were then used to construct a Euclidean
(t+1)
late responsibility γik and repeat the process. This is distance matrix which was used by t-SNE to compute
essentially the K-Means algorithm discussed in the first the two-dimensional embedding that is presented. Using
section. t-SNE directly on original data leads to a “blurring” of
This discussion of the Gaussian mixture model intro- the clusters (the reader is encouraged to test this them-
duces several concepts that we will return to repeatedly selves).
in the context of unsupervised learning. First, it is often
However, simple feature selection or feature denoising
useful to think of the visible correlations between features
(using PCA for instance) can sometimes be insufficient
in the data as resulting from hidden or latent variables.
for learning clusters due to the presence of large vari-
Second, we will often posit a generative model that en-
ations in the signal and noise of the features that are
codes the structure we think exists in the data and then
relevant for identifying the underlying clusters (Kriegel
find parameters that maximize the likelihood of the ob-
et al., 2009). Recent promising work suggests that one
served data. Third, often we will not be able to directly
way to overcome these limitations is to learn the latent
estimate the MLE, and will have to instead look for a
space and the cluster labels at the same time (Xie et al.,
computationally efficient way to find a local minimum of
2016).
the likelihood.
Finally we end the clustering section with a short dis-
cussion on clustering validation, which can be particu-
C. Clustering in high dimensions larly difficult for high-dimensional data. Often clustering
validation, i.e. verifying whether the obtained labels are
Clustering data in high-dimension can be very chal- “valid” is done by direct visual inspection. That is, the
lenging. One major problem that is aggravated in high- data is represented in a low-dimensional space and the
dimensions is the generic accumulation of noise due to cluster labels obtained are visually inspected to make
76
from the underlying distribution (this can be done know- spite the extreme simplicity relative to real-world sys-
ing only the relative probabilities) and then use these tems, Ising models exhibit a high level of intrinsic com-
samples to numerically estimate the partition function. plexity, and the degrees of freedom can become correlated
This is the philosophy behind powerful methods such as in sophisticated ways. Often, spins interact spatially lo-
Markov Chain Monte Carlo (MCMC) (Andrieu et al., cally, and respond to externally applied magnetic fields.
2003) and annealed importance sampling (Neal and Hin- A spin configuration s specifies the values si of the
ton, 1998) which are widely used in both the statistical spins at every lattice site. We can assign an “energy” to
physics and machine learning communities. An alterna- every such configuration
tive approach – which we focus on here – is to approxi- 1X X
mate the the probability distribution p(x) and partition E(s, J ) = − Jij si sj − hi si , (158)
2 i,j
function using a “variational distribution” q(x; θq ) whose i
partition function we can calculate exactly. The varia- where hi is a local magnetic field applied to the spin si ,
tional parameters θq are chosen to make the variational and Jij is the interaction strength between the spins si
distribution as close to the true distribution as possible and sj . In textbook examples, the coupling parameters
(how this is done is the focus of much of this section). J = (J, h) are typically uniform or, in studies of disor-
One of the most-widely applied examples of a varia- dered systems, (Ji , hi ) are drawn from some probability
tional method in statistical physics is Mean-Field Theory distribution (i.e. quenched disorder).
(MFT). MFT can be naturally understood as a procedure The probability of finding the system in a given spin
for approximating the true distribution of the system by configuration at temperature β −1 is given by
a factorized distribution. The deep connection between
1
MFT and variational methods is discussed below. These p(s|β, J ) = e−βE(s,J) ,
variational MFT methods have been extended to under- Zp (J )
stand more complicated spin models (also called graph-
X
Zp (β, J ) = e−βE(s,J) , (159)
ical models in the ML literature) and form the basis of {si =±1}
powerful set of techniques that go under the name of Be-
with {si =±1} denoting the sum over all possible con-
P
lief Propagation and Survey Propagation (MacKay, 2003;
figurations of the spin variables. We write Zp to em-
Wainwright et al., 2008; Yedidia et al., 2003).
phasize that this is the partition function corresponding
Variational methods are also widely used in ML to ap-
to the probability distribution p(s|β, J ), which will be-
proximate complex probabilistic models. For example,
come important later. For a fixed number of lattice sites
below we show how the Expectation Maximization (EM)
N , there are 2N possible configurations, a number that
procedure, which we discussed in the context of Gaus-
grows exponentially with the system size. Therefore, it
sian Mixture Models for clustering, is actually a general
is not in general feasible to evaluate the partition func-
method that can be derived for any latent (hidden) vari-
tion Zp (β, J ) in closed form. This represents a major
able model using a variational procedure (Neal and Hin-
practical obstacle for extracting predictions from physi-
ton, 1998). This section serves as an introduction to this
cal theories since the partition function is directly related
powerful class of variational techniques. For readers in-
to the free-energy through the expression
terested in an in-depth discussion on variational infer-
ence for probabilistic graphical models, we recommend βFp (J ) = − log Zp (β, J ) = βhE(s, J )ip − Hp ,(160)
the great treatise written by Michael I. Jordan and oth- with
ers(Jordan et al., 1999), the more physics oriented dis- X
cussion in (Yedidia, 2001; Yedidia et al., 2003), as well as Hp = − p(s|β, J ) log p(s|β, J ) (161)
David MacKay’s outstanding book (MacKay, 2003). {si =±1}
where hE(s, J )iq is the expectation value of the energy entropy separates into a sum of one-body terms
E(s, J ) with respect to the distribution q(s, θ), and Hq X
is the entropy of q(s, θ). Hq (θ) = − q(s, θ) log q(s, θ)
{si =±1}
Before proceeding further, it is helpful to introduce
a new quantity: the Kullback-Leibler divergence (KL-
X
=− qi log qi + (1 − qi ) log(1 − qi ), (166)
divergence or relative entropy) between two distributions i
p(x) and q(x). The KL-divergence measures the dissim- θi
where qi = 2 cosh
e
θi is the probability that spin si is in the
ilarity between the two distributions and is given by
+1 state. Next, we need to evaluate the average of the
q(x) Ising energy E(s, J ) with respect to the variational dis-
DKL (qkp) = Trx q(x) log , (163)
p(x) tribution q. Although the energy contains bilinear terms,
we can still evaluate this average easily, because the spins
which is the expectation w.r.t. q of the logarithmic dif- are independent (uncorrelated) in the q distribution. The
ference between the two distributions p and q. The trace mean value of spin si in the q distribution, also known
Trx denotes a sum over all possible configurations x. as the on-site magnetization, is given by
Two important properties of the KL-divergence are (i)
positivity: DKL (pkq) ≥ 0 with equality if and only if X eθi si
mi = hsi iq = si = tanh(θi ). (167)
p = q (in the sense of probability distributions), and (ii) 2 cosh θi
si =±1
DKL (pkq) 6= DKL (qkp), that is the KL-divergence is not
symmetric in its arguments. Since the spins are independent, we have
Variational mean-field theory is a systematic way for 1X X
constructing such an approximate distribution q(s, θ). hE(s, J )iq = − Jij mi mj − hi mi . (168)
2 i,j
The main idea is to choose parameters that minimize the i
1. Expectation: Given a set of assignments at iteration be the set of latent or hidden variables that we cannot di-
t, θ (t) , calculate the corresponding magnetizations rectly observe. Denote the underlying probability distri-
m(t) using Eq. (167) bution from which x and z are drawn by p(z, x|θ), with
θ representing all relevant parameters. Given a dataset
2. Maximization: Given a set of magnetizations mt , x, we wish to find the maximum likelihood estimate of
find new assignments θ(t+1) which minimize the the parameters θ that maximizes the probability of the
variational free energy Fq . From, Eq. (170) this observed data.
is just As in variational MFT, we view θ as variational pa-
rameters chosen to maximize the log-likelihood L(θ) =
(t+1) (t)
X
θi =β Jij mj + hi . (171) hlog p(x|θ)iPx , where the expectation is taken with re-
j spect to the marginal distributions of x. Algorithmically,
this can be done by iterating the variational parameters
From these equations, it is clear that we can think of the θ (t) in a series of steps (t = 1, 2, . . . ) starting from some
MFT of the Ising model as an EM-like procedure similar arbitrary initial value θ (0) :
to the one we used for K-means clustering and Gaussian
1. Expectation step (E step): Given the known
Mixture Models in Sec. XIII.
values of observed variable x and the current esti-
As is well known in statistical physics, even though
mate of parameter θt−1 , find the probability distri-
MFT is not exact, it can often yield qualitatively and
bution of the latent variable z:
even quantitatively precise predictions (especially in high
dimensions). The discrepancy between the true physics qt−1 (z) = p(z|θ (t−1) , x) (172)
and MFT predictions stems from the fact that the varia-
tional distribution q we chose cannot capture correlations 2. Maximization step (M step): Re-estimate the
between the spins. For instance, it predicts the wrong parameter θ (t) to be those with maximum likeli-
value for the critical temperature for the two-dimensional hood, assuming qt−1 (z) found in the previous step
Ising model. It even erroneously predicts the existence of is the true distribution of hidden variable z:
a phase transition in one dimension at a non-zero tem-
perature. We refer the interested reader to standard θt = arg maxhlog p(z, x|θ)iqt−1 (173)
textbooks on statistical physics for a detailed analysis θ
FIG. 59 Convergence of EM algorithm. Starting from θ (t) , 1. Expectation step: Construct the approximating
E-step (blue) establishes −Fq (θ (t) ) which is always a lower probability distribution of unobserved z given the
bound of −Fp := hlog p(x|θ)iPx (green). M-step (red) is then values of observed variable x and parameter esti-
applied to update the parameter, yielding θ (t+1) . The up- mate θ (t−1) :
dated parameter θ (t+1) is then used to construct −Fq (θ (t+1) )
in the subsequent E-step. M-step is performed again to up- qt−1 (z) = arg min Fq (θ (t−1) ) (177)
date the parameter, etc. q
the flood of new, large scale datasets has resulted in a Let us illustrate how this works in more detail. Sup-
resurgence of interest in MaxEnt models in many fields pose that we have chosen a set of functions {fi (x)} whose
including physics (especially biological physics), compu- average value we want to fix to some observed values
tational neuroscience, and ecology (Elith et al., 2011; hfi iobs . The Principle of Maximum Entropy states that
Schneidman et al., 2006; Weigt et al., 2009). MaxEnt we should choose the distribution p(x) with the largest
models are often presented as the class of generative mod- uncertainty (i.e. largest Shannon entropy Sp ), subject to
els that make the least assumptions about the underlying the constraints that the model averages match the ob-
data. However, as we have tried to emphasize throughout served averages:
the review, all ML and statistical models require assump- Z
tions, and MaxEnt models are no different. Overlooking hfi imodel = dxfi (x)p(x) = hfi iobs . (180)
this can sometimes lead to misleading conclusions, and
it is important to be cognizant of these implicit assump- We can formulate the Principle of Maximum Entropy
tions (Aitchison et al., 2016; Schwab et al., 2014). as an optimization problem using the method of Lagrange
multipliers by minimizing:
by Shannon in his landmark treatise on information the- The maximum entropy distribution is clearly just
ory (Shannon, 1949). The Shannon entropy quantifies thePusual Boltzmann distribution with energy E(x) =
the statistical uncertainty one has about the value of a − i λi fi (x). The values of the Lagrange multipliers are
random variable x drawn from a probability distribution chosen to match the observed averages for the set of func-
p(x). The Shannon entropy of the distribution is defined tions {fi (x)} whose average value is being fixed:
as
∂ log Z
Z
hfi imodel = dxp(x)fi (x) = = hfi iobs . (182)
Sp = −Trx p(x) log p(x) (179) ∂λi
In other words, the parameters of the distribution can be
where the trace is a sum/integral over all possible val-
chosen such that
ues a variable can take. Jaynes showed that the Boltz-
mann distribution follows from the Principle of Maxi- ∂λi log Z = hfi idata . (183)
mum Entropy. A physical system should be described
by the probability distribution with the largest entropy To gain more intuition for the MaxEnt distribution, it
subject to certain constraints (often provided by measur- is helpful to relate the Lagrange multipliers to the famil-
ing the average value of conserved, extensive quantities iar thermodynamic quantities we use to describe physical
such as the energy, particle number, etc.) The princi- systems (Jaynes, 1957a). Our x denotes the microscopic
ple uniquely specifies a procedure for parametrizing the state of the system, i.e. the MaxEnt distribution is a
functional form of the probability distribution. Once we probability distribution over microscopic states. How-
have specified and learned this form we can, of course, ever, in thermodynamics we only have access to average
generate new examples by sampling this distribution. quantities. If we know only the average energy hE(x)iobs ,
84
the MaxEnt procedure tells us to maximize the entropy 3. Generalized Ising Models from MaxEnt
subject to the average energy constraint. This yields
The form of a MaxEnt model is completely specified
1 −βE(x)
p(x) = e , (184) once we choose the averages {fi } we wish to constrain.
Z One common choice often used in MaxEnt modeling is to
where we have identified the Lagrange multiplier conju- constrain the first two moments of a distribution. When
gate to the energy λ1 = −β = 1/kB T with the (negative) our random variables x are continuous, the corresponding
inverse temperature. Now, suppose we also constrain the MaxEnt distribution is a multi-dimensional Gaussian. If
particle number hN (x)iobs . Then, an almost identical the x are binary (discrete), then the corresponding Max-
calculation yields a MaxEnt distribution of the functional Ent distribution is a generalized Ising (Potts) model with
form all-to-all couplings.
1 −β(E(x)−µN (x)) To see this, consider a random variable x with first
p(x) = e , (185) and second moments hxi idata and hxi xj idata , respectively.
Z
According to the Principle of Maximum Entropy, we
where we have rewritten our Lagrange multipliers in
should choose to model this variable using a Boltzmann
the familiar thermodynamic notation λ1 = −β and
distribution with constraints on the first and second mo-
λ2 = µ/β. Since this is just the Boltzmann distribu-
ments. Let ai be the Lagrange multiplier associated with
tion, we can also relate the partition function in our
hxi idata and Jij /2 be the Lagrange multiplier associated
MaxEnt model to the thermodynamic free-energy via
with hxi xj idata . Using Eq. (182), it is easy to verify that
F = −β −1 log Z. The choice of which quantities to
the energy function
constrain is equivalent to working in different thermo-
dynamic ensembles. X 1X
E(x) = − ai xi − Jij xi xj (186)
i
2 ij
2. From statistical mechanics to machine learning satisfies the above constraints.
Partition functions for maximum entropy models are
The MaxEnt idea also provides a general procedure
often intractable to compute. Therefore, it is helpful to
for learning a generative model from data. The key dif-
consider two special cases where x has different support
ference between MaxEnt models in (theoretical) physics
(different kinds of data). First, consider the case that the
and ML is that in ML we have no direct access to ob-
random variables x ∈ Rn are real numbers. In this case
served values hfi iobs . Instead, these averages must be
we can compute the partition function directly:
directly estimated from data (samples). To denote this
difference, we will call empirical averages calculated from
Z
T 1 T 1 T −1
Z = dx ea x+ 2 x Jx = (2π)n detJ −1 e− 2 a J a .
p
data as hfi idata . We can think of MaxEnt as a statisti-
cal inference procedure simply by replacing hfi iobs by (187)
hfi idata above.
This subtle change has important implications for The resulting probability density function is,
training MaxEnt models. First, since we do not know
p(x) = Z −1 e−E(x)
these averages exactly, but must estimate them from the
data, our training procedures must be careful not to over- 1 1 T −1 T 1 T
=p e 2 a J a+a x+ 2 x Jx
fit to the observations (our samples might not be reflec- (2π) detJ
n −1
with Λ controlling the regularization strength. For Λ = 0, In terms of the operators {Oi (x)}, the gradient of
there is no regularization and we are simply performing Eq. (189) takes the form (Ackley et al., 1987)
MLE. In contrast, a choice of large Λ will force many
parameters to be close to or exactly zero. Just as in ∂L({θi }) D ∂E(x; θi ) E ∂ log Z({θi })
− = +
regression, an L1 penalty enforces sparsity, with many of ∂θi ∂θi data ∂θi
the θi set to zero, and L2 regularization shrinks the size = hOi (x)idata − hOi (x)imodel . (195)
of the parameters towards zero.
One challenge of generative models is that it is often These equations have a simple and beautiful interpre-
difficult to choose the regularization strength Λ. Recall tation. The gradient of the log-likelihood with respect to
that, for linear and logistic regression, Λ is chosen to a model parameter is a difference of moments – one calcu-
maximize the out-of-sample performance on a validation lated directly from the data and one calculated from our
dataset (i.e. cross-validation). However, for generative model using the current model parameters. The data-
models our data are usually unlabeled. Therefore, choos- dependent term is known as the positive phase of the
ing a regularization strength is more subtle and there ex- gradient and the model-dependent term is known as the
ists no universal procedure for choosing Λ. One common negative phase of the gradient. This derivation also gives
strategy is to divide the data into a training set and a val- an intuitive explanation for likelihood-based training pro-
idation set and monitor a summary statistic such as the cedures. The gradient acts on the model to lower the en-
log-likelihood, energy distance (Székely, 2003), or varia- ergy of configurations that are near observed data points
tional free-energy of the generative model on the training while raising the energy of configurations that are far
and validation sets (the variational free-energy was dis- from observed data points. Finally, we note that all infor-
cussed extensively in Sec. XIV ) (Hinton, 2012). If the mation about the data only enters the training procedure
gap between the training and validation datasets starts through the expectations hOi (x)idata and our generative
growing, one is probably overfitting the model even if model is blind to information beyond what is contained
the log-likelihood of the training dataset is still increas- in these expectations.
ing. This also gives a procedure for “early stopping” – To use SGD, we must still calculate the expectation
a regularization procedure we introduced in the context values that appear in Eq. (195). The positive phase of the
of discriminative models. In practice, when using such gradient – the expectation values with respect to the data
regularizers it is important to try many different values – can be easily calculated using samples from the training
of Λ and then try to use a proxy statistic for overfitting dataset. However, the negative phase – the expectation
to evaluate the optimal choice of Λ. values with respect to the model – is generally much more
difficult to compute. We will see that in almost all cases,
we will have to resort to either numerical or approximate
methods. The fundamental reason for this is that it is
D. Computing gradients impossible to calculate the partition function exactly for
most interesting models in both physics and ML.
We still need to specify a procedure for minimizing the There are exceptional cases in which we can calcu-
cost function. One powerful and common choice that late expectation values analytically. When this hap-
is widely employed when training energy-based models pens, the generative model is said to have a Tractable
is stochastic gradient descent (SGD) (see Sec. IV). Per- Likelihood. One example of a generative model with a
forming MLE using SGD requires calculating the gradi- Tractable Likelihood is the Gaussian MaxEnt model for
ent of the log-likelihood Eq. (189) with respect to the real valued data discussed in Eq. (188). The param-
parameters θi . To simplify notation and gain intuition, eters/Lagrange multipliers for this model are the local
it is helpful to define “operators” Oi (x), conjugate to the fields a and the pairwise coupling matrix J. In this case,
parameters θi the usual manipulations involving Gaussian integrals al-
low us to exactly find the parameters µ = −J −1 a and
∂E(x; θi ) Σ = −J −1 , yielding the familiar expressions µ = hxidata
Oi (x) = . (193) and Σ = h(x − hxidata )(x − hxidata )T idata . These are
∂θi
the standard estimates for the sample mean and covari-
Since the partition function is just the cumulant gener- ance matrix. Converting back to the Lagrange multipli-
ating function for the Boltzmann distribution, we know ers yields
that the usual statistical mechanics relationships between
expectation values and derivatives of the log-partition J = −h(x − hxidata )(x − hxidata )T i−1
data . (196)
function hold:
Returning to the generic case where most energy-based
∂ log Z({θi }) models have intractable likelihoods, we must estimate ex-
hOi (x)imodel = Trx pθ (x)Oi (x) = − . (194) pectation values numerically. One way to do this is draw
∂θi
87
samples Smodel = {x0i } from the model pθ (x) and evalu- In practice, it is helpful to supplement this basic proce-
ate arbitrary expectation values using these samples: dure with some tricks that help training. As with dis-
Z criminative neural networks, it is important to initial-
ize the parameters properly and print summary statistics
X
hf (x)imodel = dxpθ (x)f (x) ≈ f (x0i ). (197)
x0i ∈Smodel
during the training procedure on the training and vali-
dation sets to prevent overfitting. These and many other
The samples from the model x0i ∈ Smodel are often re- “cheap tricks” have been nicely summarized in a short
ferred to as fantasy particles in the ML literature and note from the Hinton group (Hinton, 2012).
can be generated using simple MCMC algorithms such A major computational and practical limitation of
as Metropolis-Hasting which are covered in most modern these methods is that it is often hard to draw samples
statistical physics classes. However, if the reader is unfa- from generative models. MCMC methods often have long
miliar with MCMC methods or wants a quick refresher, mixing-times (the time one has to run the Markov chain
we recommend the concise and beautiful discussion of to get uncorrelated samples) and this can result in bi-
MCMC methods from both the physics and ML point- ased sampling. Luckily, we often do not need to know
of-view in Chapters 29-32 of David MacKay’s masterful the gradients exactly for training ML models (recall that
book (MacKay, 2003). noisy gradient estimates often help the convergence of
Finally, we note that once we have the fantasy particles gradient descent algorithms), and we can significantly re-
from the model, we can also easily calculate the gradient duce the computational expense by running MCMC for
of any expectation value hf (x)imodel using what is com- a reasonable time window. We will exploit this observa-
monly called the “log-derivative trick” in ML (Fu, 2006; tion extensively in the next section when we discuss how
Kleijnen and Rubinstein, 1996): to train more complex energy-based models with hidden
variables.
∂ ∂pθ (x)
Z
hf (x)imodel = dx f (x)
∂θi ∂θi
D ∂ log p (x) E
θ
= f (x)
∂θi model
= hOi (x)f (x)imodel XVI. DEEP GENERATIVE MODELS: HIDDEN VARIABLES
AND RESTRICTED BOLTZMANN MACHINES (RBMS)
X
≈ Oi (xj )f (x0j ). (198)
x0j ∈Smodel
The last section introduced many of the core ideas be-
This expression allows us to take gradients of more com- hind energy-based generative models. Here, we extend
plex cost functions beyond the MLE procedure discussed this discussion to energy-based models that include la-
here. tent or hidden variables.
Including latent variables in generative models greatly
enhances their expressive power – allowing the model to
E. Summary of the training procedure
represent sophisticated correlations between visible fea-
tures without sacrificing trainability. By having multiple
We now summarize the discussion above and present a layers of latent variables, we can even construct powerful
general procedure for training an energy based model us- deep generative models that possess many of the same
ing SGD on the cost function (see Sec. IV). Our goal is to desirable properties as deep, discriminative neural net-
fit the parameters of a model pλ ({θi }) = Z −1 e−E(x,{θi }) . works.
Training the model involves the following steps:
We begin with a discussion that tries to provide a sim-
1. Read a minibatch of data, {x}. ple intuition for why latent variables are such a pow-
erful tool for generative models. Next, we introduce a
2. Generate fantasy particles {x0 } ∼ pλ using an powerful class of latent variable models called Restricted
MCMC algorithm (e.g., Metropolis-Hastings). Boltzmann Machines (RBMs) and discuss techniques for
training these models. After that, we introduce Deep
3. Compute the gradient of log-likelihood using these Boltzmann Machines (DBMs), which have multiple layers
samples and Eq. (195), where the averages are of latent variables. We then introduce the new Paysage
taken over the minibatch of data and the fantasy package for training energy-based models and demon-
particles from the model, respectively. strate how to use it on the MNIST dataset and sam-
ples from the Ising model. We conclude by discussing
4. Use the gradient as input to one of the gradient recent physics literature related to energy-based genera-
based optimizers discussed in section Sec. IV. tive models.
88
A. Why hidden (latent) variables? model with this form of the energy function is known as
the Hopfield model (Amit et al., 1985; Hopfield, 1982).
Latent or hidden variables are a powerful yet elegant The Hopfield model has played an extremely important
way to encode sophisticated correlations between observ- role in statistical physics, computational neuroscience,
able features. The underlying reason for this is that and machine learning, and a full discussion of its prop-
marginalizing over a subset of variables – “integrating erties is well beyond the scope of this review [see (Amit,
out” degrees of freedom in the language of physics – in- 1992) for a beautiful discussion that combines all these
duces complex interactions between the remaining vari- perspectives]. Therefore, here we refer to all energy func-
ables. The idea that integrating out variables can lead tions of the form Eq. (200) as (generalized) Hopfield mod-
to complex correlations is a familiar component of many els, even for the case when the Wiµ are continuous vari-
physical theories. For example, when considering free ables.
electrons living on a lattice, integrating out phonons gives We now “decouple” the visible variables vi by intro-
rise to higher-order electron-electron interactions (e.g. su- ducing a set of normally, distributed continuous latent
perconducting or magnetic correlations). More generally, variables hµ (in condensed matter language we perform a
in the Wilsonian renormalization group paradigm, all ef- Hubbard-Stratonovich transformation). Using the usual
fective field theories can be thought of as arising from identity for Gaussian integrals, we can rewrite the Boltz-
integrating out high-energy degrees of freedom (Wilson mann distribution for the generalized Hopfield model as
and Kogut, 1974).
ai vi + 21
P P
vi Wiµ Wjµ vj
Generative models with latent variables run this logic e i ijµ
p(v) =
in reverse – encode complex interactions between visible Z
variables by introducing additional, hidden variables that h2µ −
1
P P
dhµ e− 2 vi Wiµ hµ
P
i a i vi
Q R
e µ
µ i
We note that in the special case when both vi ∈ A Restricted Boltzmann Machine (RBM) is an energy-
{−1, +1} and Wiµ ∈ {−1, +1} are binary variables, a based model with both visible and hidden units where the
89
above). When the the visible units are continuous and (n)
the hidden units are discrete, the RBM is often called a Kµ (t) is defined such that the nth cumulant is κµ =
Gaussian Bernoulli Restricted Boltzmann Machine (Dahl ∂tn Kµ |t=0 .
et al., 2010; Hinton and Salakhutdinov, 2006). It is even The cumulant generating function appears in the
possible to perform multi-modal learning with a mixture marginal free-energy of the visible units, which can be
90
C. Training RBMs where the positive expectation with respect to the data
is understood to mean sampling from the model while
clamping the visible units to their observed values in the
RBMs are a special class of energy-based generative
data. As before, calculating the negative phase of the
models, which can be trained using the Maximum Like-
gradient (i.e. the expectation value with respect to the
lihood Estimation (MLE) procedure described in detail
model) requires that we draw samples from the model.
in Sec. XV. To briefly recap, first, we must choose a cost
Luckily, the bipartite form of the interactions in RBMs
function – for MLE this is just the negative log-likelihood
were specifically chosen with this in mind.
with or without an additional regularization term to pre-
vent overfitting. We then minimize this cost function us-
ing one of the Stochastic Gradient Descent (SGD) meth- 1. Gibbs sampling and contrastive divergence (CD)
ods described in Sec. IV.
The gradient itself can be calculated using Eq. (195). The bipartite interaction structure of an RBM makes it
For example, for the Bernoulli-Bernoulli RBM in possible to calculate expectation values using a Markov
91
Chain Monte Carlo (MCMC) method known as Gibbs as is often the case in ML, CD-n sacrifices the ability
sampling. The key reason for this is that since there are to generalize to some extent in order to make the model
no interactions of visible units with themselves or hidden easier to train.
units with themselves, the visible and hidden units of an Some of these undesirable features can be tempered
RBM are conditionally independent: by using a slightly different variant of CD called Persis-
Y tent Contrastive Divergence (PCD) (Tieleman and Hin-
p(v|h) = p(vi |h) ton, 2009). In PCD, rather than restarting the Gibbs
i sampler from the data at each gradient descent step, we
start the Gibbs sampling at the fantasy particles in the
Y
p(h|v) = p(hµ |v), (212)
µ last gradient descent step (see Fig. 62). Since parameters
change slowly compared to the Gibbs sampling, samples
with that are high probability at one step of the SGD are also
likely to be high probability at the next step. This en-
X
p(vi = 1|h) = σ(ai + Wiµ hµ ) (213)
µ
sures that PCD does not introduce large errors in the
X estimation of the gradients. The advantage of using fan-
p(hµ = 1|v) = σ(bµ + Wiµ vi ) tasy particles to initialize the Gibbs sampler is to allow
i PCD to explore parts of the feature space that are much
and where σ(z) = 1/(1 + e−z ) is the sigmoid function. further from the training dataset than one could reach
Using these expressions it is easy to compute expec- with ordinary CD.
tation values with respect to the data. The input to We note that, in applications using RBMs as a vari-
gradient descent is a minibatch of observed data. For ational ansatz for quantum states, Gibbs sampling is
each sample in the minibatch, we simply clamp the visi- not necessarily the best option for training, and in prac-
ble units to the observed values and apply Eq. (213) using tice parallel tempering or other Metropolis schemes can
the probability for the hidden variables. We then average outperform Gibbs sampling. In fact, Gibbs sampling is
over all samples in the minibatch to calculate expectation not even feasible with complex-valued weights required
values with respect to the data. To calculate expectation for quantum wavefucntions, whereas Metropolis schemes
values with respect to the model, we use (block) Gibbs might be feasible (Carleo, 2018).
sampling. The idea behind (block) Gibbs sampling is
to iteratively sample from the conditional distributions
ht+1 ∼ p(h|vt ) and vt+1 ∼ p(v|ht+1 ) (see Figure 62, 2. Practical Considerations
top). Since the units are conditionally independent, each
step of this iteration can be performed by simply draw- The previous section gave an overview of how to train
ing random numbers. The samples are guaranteed to RBMs. However, there are many “tricks of the trade”
converge to the equilibrium distribution of the model in that are missing from this discussion. Luckily, a succinct
the limit that t → ∞. At the end of the Gibbs sampling summary of these has been compiled by Geoff Hinton and
procedure, one ends up with a minibatch of samples (fan- published as a note that readers interested in training
tasy particles). RBMs are urged to consult (Hinton, 2012).
One drawback of Gibbs sampling is that it may take For completeness, we briefly list some of the important
many back and forth iterations to draw an independent points here:
sample. For this reason, the Hinton group introduced • Initialization.—The model must be initialized.
an approximate Gibbs sampling technique called Con- Hinton suggests taking the weights Wiµ from a
trastive Divergence (CD) (Hinton, 2002; Hinton et al., Gaussian with mean zero and standard deviation
2006). In CD-n, we just perform n iterations of (block) σ = 0.01 (Hinton, 2012). An alternative initial-
Gibbs sampling, with n often taken to be as small as 1 ization scheme proposed by Glorot and Bengio in-
(see Figure 62)! The price for this truncation is, of course, stead chooses the standard deviation to scale with
that we are not drawing samples from the true model dis- √
the size of the layers: σ = 2/ Nv + Nh where Nv
tribution. But for our purpose – using the expectations and Nh are number of visible and hidden units re-
to estimate the gradient for SGD – CD-n has proven to spectively (Glorot and Bengio, 2010). The bias of
work reasonably well. As long as the approximate gra- the hidden units is initialized to zero while the bias
dients are reasonably correlated with the true gradient, of the visible units is typically taken to be inversely
SGD will move in a reasonable direction. CD-n of course proportional to the mean activation, ai = hvi i−1 data .
does come at a price. Truncating the Gibbs sampler pre-
vents sampling far away from the starting point, which • Regularization.—One can of course use an L1 or
for CD-n are the data points in the minibatch. Therefore, L2 penalty, typically only on the weight parame-
our generative model will be much more accurate around ters, not the biases. Alternatively, Dropout has
regions of feature space close to our training data. Thus, been shown to decrease overfitting when training
92
Below, we show how to build and train four different kinds of models: (i) a “Hopfield” type RBM with Gaussian
hidden units and Bernoulli (binary) visible units, (ii) a conventional RBM where both the visible and hidden units are
Bernoulli, (iii) a conventional RBM with an additional L1-penalty that enforces sparsity, and (iv) a Deep Boltzmann
Machine (DBM) with three Bernoulli layers with L1 penalty each. In the following, we demonstrate the simplicity of
using Paysage walking the reader step-by-step through short snippets of code.
We kick off by loading the required packages. Note that Paysage requires Python 3.6 or higher (see additional guides
to install the package in Notebook 17). We also fix the seed of the random number generator to ensure reproducibility
of our numerical experiment.
1 from __future__ import print_function, division
2 import os
3 import paysage
4 from [Link] import Model # model constructor
5 from paysage import optimizers # optimizer
6 from [Link] import BernoulliLayer, GaussianLayer # layers
7 from [Link] import DataShuffler, HDFBatch # data handler
8 from [Link] import ProgressMonitor, SequentialMC, SGD, LayerwisePretrain, pcd
9 from [Link] import PowerLawDecay # hyperparameter schedule
10 from [Link].model_utils import State
11 from [Link] import l1_penalty # regularization
12 # fix random seed to ensure deterministic behavior
13 [Link].set_seed(137)
We want to study the MNIST digit dataset. A preprocessed version of the data is conveniently built into Paysage,
and our first task is to download it. To this end, let us fetch the directory Paysage was installed in and print it:
15 ##### download and preprocess MNIST data in Paysage
16 # fetch paysage directory
17 paysage_path=[Link]([Link](paysage.__file__))
18 mnist_path=[Link](paysage_path, "mnist", "mnist.h5")
19 shuffled_mnist_path=[Link](paysage_path, "mnist", "shuffled_mnist.h5")
20 # print path to Paysage directory
21 print(paysage_path)
To download the data, open up a terminal and navigate to the Paysage directory to run python
mnist/download_mnist.py. We can also check if the data have been successfully downloaded:
94
If this is the first time using the data set, we need to shuffle it. This step is necessary, since we shall shortly
employ SGD-based algorithms in the training process (cf. Sec. IV) which requires using small minibatches of data to
compute the gradient at each step. If the data have an order, then the estimates for the gradients computed from the
minibatches will be biased. Shuffling the data ensures that the gradient estimates are unbiased (though still noisy).
The data can be compressed by setting complevel to a value greater than 0, but we do not use that here.
26 ##### set up minibatch data generator
27 # shuffle data if running for the first time
28 if not [Link](shuffled_mnist_path):
29 DataShuffler(mnist_path,shuffled_mnist_path,complevel=0).shuffle()
Next, we create a python generator, which splits the data into a training and validation sets, and separates them into
minibatches of size batch_size. Before we begin training, we set data to training mode.
30 # batch size
31 batch_size=100
32 # create data generator object with minibathces
33 data=HDFBatch(shuffled_mnist_path,’train/images’, batch_size,
34 transform=[Link].binarize_color,train_fraction=0.95)
35 # reset the data generator in training mode
36 data.reset_generator(mode=’train’)
To monitor the progress of performance metrics during training, we define the variable performance which tells
Paysage to measure the reconstruction error from the validation set. Possible metrics include the reconstruction error
(used in this example) and metrics related to difference in the energy of random samples and samples from the model
(see [Link] in Paysage documentation for a complete list).
37 # the reconstruction error will be computed from the validation set
38 performance=ProgressMonitor(data,metrics=[’ReconstructionError’])
Having loaded and preprocessed the data, we now move on to construct a hopfield model. To do this, we use the
Model class with a visible BernoulliLayer and a hidden GaussianLayer. Note that the visible layer has the same
size as the input data points, which is read off [Link]. The number of hidden units is num_hidden_units.
We also standardize the mean and variance of the Gaussian layer setting them to zero and unity, respectively (the
nomenclature of Paysage here is inspired by the terminology in Variational Autoencoders, cf. Sec. XVII).
40 ##### create hopfield model
41 # hidden units
42 num_hidden_units=200
43 # set up the model
44 hopfield=Model([BernoulliLayer([Link]), # visible layer
45 GaussianLayer(num_hidden_units) # hidden layer
46 ])
47 # set mean and standard deviation of hidden layer to to 0 and 1, respectively
48 [Link][1].set_fixed_params([’loc’, ’log_var’])
We choose to train the model with the Adam optimizer. To ensure convergence, we attenuate the learning_rate
hyperparameter according to a PowerLawDecay schedule: learning_rate(t) =initial/(1+coefficient×t). It will
prove convenient to define the function Adam_optimizer for this purpose.
49 # set up an optimizer method (ADAM in this case)
50 def ADAM_optimizer(initial,coefficient):
51 # define learning rate attenuation schedule
52 learning_rate=PowerLawDecay(initial=initial,coefficient=coefficient)
53 # return optimizer object
54 return [Link](stepsize=learning_rate)
Next, we have to create the model. First, we initialize the model using the initialize function method which accepts
the data as a required argument. We choose the initialization routine glorot, cf. discussion in Sec XVI.C.2. Second,
we define an optimizer calling the function Adam_optimizer, and store the object under the name opt. To create
95
an MCMC sampler, we use the method from_batch of the SequentialMC class, passing the model and the data.
Next, we create an SGD object called trainer to train the model using Persistent Contrastive Divergence (pcd) with
a fixed number of monte_carlo_steps. We can also monitor the reconstruction error during training. Last, we train
the model in epochs (cf. variable num_epochs), calling the train() method of trainer. These steps are universal
for shallow generative models, and it is convenient to combine them in the function train_model, which we shall use
repeatedly.
55 # define function to compile and train model
56 num_epochs=20 # training epochs
57 monte_carlo_steps=1 # number of MC sampling steps
58 def train_model(model,num_epochs,monte_carlo_steps,performance):
59 # make a simple guess for the initial parameters of the model
60 [Link](data,method=’glorot_normal’)
61 # set optimizer
62 opt=ADAM_optimizer(1E-2,1.0)
63 # set up a Monte Carlo sampler
64 sampler=SequentialMC.from_batch(model,data)
65 # use persistent contrastive divergence to fit the model
66 trainer=SGD(model,data,opt,num_epochs,sampler,
67 method=pcd,mcsteps=monte_carlo_steps,monitor=performance)
68 # train model
69 [Link]()
70 # train hopfield model
71 train_model(hopfield,num_epochs,monte_carlo_steps,performance)
Let us now show how to build a few more generative models with Paysage. We can easily create a Bernoulli RBM
and train it using the functions defined above as follows:
73 ##### Bernoulli RBM
74 rbm = Model([BernoulliLayer([Link]), # visible layer
75 BernoulliLayer(num_hidden_units) # hidden layer
76 ])
77 # train Bernoulli RBM
78 train_model(rbm,num_epochs,monte_carlo_steps,performance)
Constructing a Bernoulli RBM with L1 regularization is also straightforward in Paysage, using the add_penalty
method which accepts a dictionary as an input. Some layers may have multiple properties (such as the location and
scale parameters of a Gaussian layer) so the dictionary key specifies which property the penalty should be applied to
80 ##### Bernoulli RBM with L1 regularizer
81 rbm_L1 = Model([BernoulliLayer([Link]), # visible layer
82 BernoulliLayer(num_hidden_units) # hidden layer
83 ])
84 # add an L1 penalty to the weights
85 rbm_L1.weights[0].add_penalty({’matrix’:l1_penalty(1e-3)})
86 # train Bernoulli RBM with L1 regularizer
87 train_model(rbm_L1,num_epochs,monte_carlo_steps,performance)
To define a deep Boltzmann machine (DBM), we just add more layers, and an L1 penalty for every layer.
89 ##### Deep Boltzmann Machine
90 # set up the model
91 dbm_L1 = Model([BernoulliLayer([Link]), # visible layer
92 BernoulliLayer(num_hidden_units), # hidden layer 1
93 BernoulliLayer(num_hidden_units) # hidden layer 2
94 ])
95 # add an L1 penalty to the weights
96 for weight in dbm_L1.weights:
97 weight.add_penalty({’matrix’:l1_penalty(1e-3)})
Recalling the essential trick with layer-wise pre-training to prepare the weights of the DBM, we define a pretrainer
as an object of the LayerwisePretrain class (see code snippet below). This results in a slight modification of the
function train_model.
98 # add pre-training
99 def train_model(model,num_epochs,monte_carlo_steps, performance):
96
100 # make a simple guess for the initial parameters of the model
101 [Link](data,method=’glorot_normal’)
102 # set SGD retrain optimizer
103 opt=ADAM_optimizer(1E-2,1.0)
104 # set up a Monte Carlo sampler
105 sampler = SequentialMC.from_batch(model, data)
106 # check if model is deep
107 is_deep = model.num_layers > 2
108 if is_deep:
109 print("layerwise pretraining")
110 pretrainer = LayerwisePretrain(model, data, opt, num_epochs,
111 method=pcd, mcsteps=monte_carlo_steps,
112 metrics=[’ReconstructionError’])
113 [Link]()
114 # reset the optimizer using a lower learning rate
115 opt = ADAM_optimizer(initial/10.0, coefficient)
116 print("use persistent contrastive divergence to fit the model")
117 # use persistent contrastive divergence to fit the model
118 trainer=SGD(model,data,opt,num_epochs,sampler,
119 method=pcd,mcsteps=monte_carlo_steps,monitor=performance)
120 # train model
121 [Link]()
122 # train DBM
123 train_model(dbm_L1,num_epochs,monte_carlo_steps,performance)
Having trained our models, let us see how they perform by computing some reconstructions and fantasy particles
from the validation data. Recall that a reconstruction v0 of a given data point x is computed in two steps: (i) we
fix the visible layer v = x to be the data, and use MCMC sampling to find the state of the hidden layer h which
maximizes the probability distribution p(h|v). (ii) fixing the same obtained state h, we find the reconstruction v0 of
the original data point which maximizes the probability p(v0 |h). In the case of a DBM, the forward pass continues
until we reach the last of the hidden layers, and the backward pass goes in reverse. A configuration sampled from an
RBM needs to specify the values of both the visible and hidden units. Since the data only specify the visible units,
we need to initialize some hidden unit values. The visible and hidden units are stored in a State object. To compute
reconstructions, we define an MCMC sampler based on the trained model. The stating point for the MCMC sampler is
set using the set_state() method. To compute reconstructions, we need to keep the probability distribution learned
by the generative model fixed which is done with the help of the deterministic_iteration function method, that
takes in its first argument the number of passes (1 for a single v → h → v0 pass), and the state of the sampler
[Link] as required arguments. We can combine these steps in the function compute_reconstructions.
Figure 60 shows the result.
125 ##### compute reconstructions
126 def compute_reconstructions(model, data):
127 """
128 Computes reconstructions of the input data.
129 Input v -> h -> v’ (one pass up one pass down)
130
131 Args:
132 model: a model
133 data: a tensor of shape (num_samples, num_visible_units)
134
135 Returns:
136 tensor of shape (num_samples, num_visible_units)
137
138 """
139 # a configuration sampled from an RBM needs to specify the values
140 # of both the visible and hidden units
141 # since the data only specify the visible units, we need to initialize
142 # some hidden unit values
143 # the visible and hidden units are stored in a State object
144 data_state=State.from_visible(data,model)
145 # define MC sampler
146 sampler=SequentialMC(model)
147 # define a starting point for MC sampler
97
148 sampler.set_state(data_state)
149 # compute reconstructions
150 recons=model.deterministic_iteration(1,[Link]).units[0]
151 #
152 return [Link].to_numpy_array(recons)
Once we have the trained models ready, we can use MCMC to draw samples from the corresponding probability
distributions, the so-called fantasy particles. To this end, let us draw a random_sample from the validation data and
compute the model_state. Next, we define an MCMC sampler based on the model, and set its state to model_state.
To compute the fantasy particles, we do layer-wise Gibbs sampling for a total of n_steps equilibration steps. The
last step (controlled by the boolean mean_field) is a final mean-field iteration [see the tricks discussed in (Hinton,
2012)]. Figure 64 shows the result.
154 ##### compute fantasy particles
155 def compute_fantasy_particles(model,data,num_steps,mean_field=True):
156 """
157 Draws samples from the model using Gibbs sampling Markov Chain Monte Carlo .
158 Starts from randomly initialized points.
159
160 Args:
161 model: a model
162 data: a tensor of shape (num_samples, num_visible_units)
163 num_steps (int): the number of update steps
164 mean_field (bool; optional): run a final mean field step to compute probabilities
165
166 Returns:
167 tensor of shape (num_samples, num_visible_units)
168
169 """
170 # compute random data sample
171 random_sample=[Link](data)
172 # get model state from visible layer
173 model_state=State.from_visible(random_sample,model)
174 # define MC sampler
175 sampler=SequentialMC(model)
176 # change sampler state
177 sampler.set_state(model_state)
178 # does n_steps forward and backward passes
179 sampler.update_state(num_steps)
180 if mean_field: # see Hinton’s 2012 paper: trick (practical guide for training)
181 fantasy_particles=model.mean_field_iteration(1,[Link]).units[0]
182 else:
183 fantasy_particles=[Link][0]
184 #
185 return [Link].to_numpy_array(fantasy_particles)
One can use generative models to reduce the noise in images (de-noising). Let us randomly flip a fraction,
fraction_to_flip, of the black&white bits in the validation data, and use the models defined above to reconstruct
(de-noise) the digit images. Figure 65 shows the result.
187 ##### denoise MNIST images
188 # get validation data
189 examples = [Link](mode=’validate’) # shape (batch_size, 784)
190 # reset data generator to beginning of the validation set
191 data.reset_generator(mode=’validate’)
192 # add some noise to the examples by randomly flipping some pixels 0 -> 1 and 1 -> 0
193 fraction_to_flip=0.15
194 # create flipping mask
195 flip_mask=[Link].rand_like(examples) < fraction_to_flip
196 # compute noisy data
197 noisy_data=(1-flip_mask) * examples + flip_mask * (1 - examples)
198 # define number of digits to display
199 num_to_display=8
200 # compute de-noised images
98
201 hopfield_denoised=compute_reconstructions(hopfield,noisy_data[:num_to_display])
202 rbm_denoised=compute_reconstructions(rbm,noisy_data[:num_to_display])
203 rbmL1_denoised=compute_reconstructions(rbmL1,noisy_data[:num_to_display])
204 dbm_L1_denoised=compute_reconstructions(dbm_L1,noisy_data[:num_to_display])
The full code used to generate Figs. 60, 64 and 65 is available in Notebook 17.
FIG. 66 MC samples, their reconstructions and fantasy particles generated by a Deep Boltzmann Machine in the ordered
phase of the the 2D Ising data set at T /J = 1.75. We used two hidden layers of 1000 and 100 layers, respectively.
FIG. 67 MC samples, their reconstructions and fantasy particles generated by a Deep Boltzmann Machine in the critical
regime of the the 2D Ising data set at T /J = 2.25. We used two hidden layers of 1000 and 100 layers, respectively.
tably, RBM-inspired variational ansatzes were used to XVII. VARIATIONAL AUTOENCODERS (VAES) AND
learn both complex-valued wavefunctions and the real- GENERATIVE ADVERSARIAL NETWORKS (GANS)
valued probability distribution associated with the abso-
lute square of a quantum state (Carleo et al., 2018; Car-
leo and Troyer, 2017; Freitas et al., 2018; Nomura et al., In the previous two sections, we considered energy-
2017; Torlai et al., 2018) and, in this context, RBMs are based generative models. Here, we extend our discus-
sometimes called Born machines (Cheng et al., 2017). sion to two new generative model frameworks that have
Further applications include the detection of order in gained wide appeal in the the last few years: generative
low-energy product states (Rao et al., 2017), and learning adversarial networks (GANs) (Goodfellow, 2016; Good-
Einstein-Podolsky-Rosen correlations on an RBM (Wein- fellow et al., 2014; Radford et al., 2015) and variational
stein, 2017). Inspired by the success of tensor networks in autoencoders (VAEs) (Kingma and Welling, 2013). Un-
physics, the latter have been used as a basis for both gen- like energy-based models, both these generative modeling
erative and discriminative learning (Huggins et al., 2019): frameworks are based on differentiable neural networks
RBMs (Chen et al., 2018) were used to extract the spa- and consequently can be trained using backpropagation-
tial geometry from entanglement (You et al., 2017), and based methods. VAEs, in particular, can be easily im-
generative models based on matrix product states have plemented and trained using high-level packages such as
been developed (Han et al., 2017). Last but not least, Keras making them an easy-to-deploy generative frame-
Quantum entanglement was studied using RBM-encoded work. These models also differ from the energy-based
states (Deng et al., 2017) and tensor product based gen- models in that they do not directly seek to maximize like-
erative models have been used to understand MNIST and lihood. GANs, for example, employ a novel cost function
other ML datasets (Stoudenmire and Schwab, 2016). based on adversarial learning (a concept we motivate and
explain below). Finally we note that VAEs and GANs are
100
FIG. 68 MC samples, their reconstructions and fantasy particles generated by a Deep Boltzmann Machine in the disordered
phase of the the 2D Ising data set at T /J = 2.75. We used two hidden layers of 1000 and 100 layers, respectively.
already starting to make their way into physics (Heimel satisfy the triangle inequality.
et al., 2018; Liu et al., 2017; Rocchetto et al., 2018; Wet- Given two distributions, there are two distinct KL-
zel, 2017) and astronomy (Ravanbakhsh et al., 2017), and divergences we can construct:
methods from physics may prove useful for furthering
our understanding of these methods (Alemi and Abbara, p(x)
Z
DKL (p||q) = dxp(x) log (214)
2017). More generally, GANs have found important ap- q(x)
plications in many artistic and image manipulation tasks Z
q(x)
(see references in (Goodfellow, 2016)). DKL (q||p) = dxq(x) log . (215)
p(x)
The section is organized as follows. We start by moti-
vating adversarial learning by discussing the limitations A related quantity called the Jensen-Shannon divergence,
of maximum likelihood based approaches. We then give
a high-level introduction to the main idea behind gen- 1
p + q
p + q
DJS (p, q) = DKL p
+ DKL q
erative adversarial networks and discuss how they over- 2 2
2
come some of these limitations, simultaneously highlight-
ing both the power of GANs and some of the difficulties. does satisfy all of the properties of a squared metric (i.e.,
We then show how VAEs integrate the variational meth- the square root of the Jensen-Shannon divergence is a
ods introduced in Sec. XIV with deep, differentiable neu- metric). An important property of the KL-divergence
ral networks to build more powerful generative models that we will make use of repeatedly is its positivity:
that move beyond the Expectation Maximization (EM). DKL (p||q) ≥ 0 with equality if and only if p(x) = q(x)
We then briefly discuss VAEs from an information theo- almost everywhere.
retic perspective, before discussing practical tips for im- In generative models in ML, the two distributions
plementing and training VAEs. We conclude by using we are usually concerned with are the model distribu-
VAEs on examples using the Ising and MNIST datasets tion pθ (x) and the data distribution pdata (x). We of
(see also Notebooks 19 and 20). course would like these models to be as similar as possi-
ble. However, as we discuss below, there are many sub-
tleties about how we measure similarities that can have
large consequences for the behavior of training proce-
A. The limitations of maximizing Likelihood
dures. Maximizing the log-likelihood of the data under
the model is the same as minimizing the KL divergence
The Kullback-Leibler (KL)-divergence plays a central
between the data distribution and the model distribu-
role in many generative models. Developing an intuition
tion DKL (pdata ||pθ ). To see this, we can rewrite the KL
about KL-divergences is one of the keys to understanding
divergence as:
why adversarial learning has proved to be such a powerful
method for generative modeling. Here, we revisit the KL- Z
divergence with an eye towards understanding GANs and DKL (pdata ||pθ ) = dxpdata (x) log pdata (x)
motivate adversarial learning. The KL-divergence mea- Z
sures the similarity between two probability distributions − dxpdata (x) log pθ (x)
p(x) and q(x). Strictly speaking, the KL divergence is
not a metric because it is not symmetric and does not = −S[pdata ] − hlog pθ (x)idata (216)
101
∆ = 2.0 ∆ = 5.0
Model Model
Data Data
15
10
0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 −4 −2 0 2 4
∆ x
FIG. 69 KL-divergences between the data distribution pdata
and the model pθ . Data is drawn from a bimodal Gaus- Data-Model
sian distribution with unit variances peaked at ±∆ with 5
Model-Data
∆ = 2.0 and the model pθ (x) is a Gaussian with mean zero 4
and same variance as pθ (x). (Top) pdata and pθ for ∆ = 2.
(Bottom) DKL (pdata ||pθ ) (Data-Model) and DKL (pθ ||pdata ) 3
(Model-Data) as a function of ∆. Notice that DKL (pdata ||pθ )
is insensitive to placing weight in the model distribution in 2
regions where pdata ≈ 0 whereas DKL (pθ ||pdata ) punishes this
harshly. 1
0
0 2 4 6 8 10
Rearranging this equation, we have ∆
hlog pθ (v)idata = −S[pdata ] − DKL (pdata ||pθ ) (217) FIG. 70 KL-divergences between the data distribution pdata
and the model pθ . Data is drawn from a Gaussian mixture
The equivalence follows from the positivity of KL- of the form pdata = 0.25N (−∆) + 0.25 ∗ N (∆) + 0.5N (0)
divergence and the fact that the entropy of the data where N (a) is a normal distribution with unit variance cen-
distribution is constant. In contrast, the original for- tered at x = a. pθ (x) is a Gaussian with σ 2 = 2. (Top)
mulation of GANs minimizes an upper bound on the pdata and pθ for ∆ = 5. (Middle) pdata and pθ for ∆ = 1.
Jensen-Shannon divergence between the model distribu- (Bottom) DKL (pdata ||pθ ) [Data-Model] and DKL (pθ ||pdata )
tion pθ (x) and the data distribution pdata (x) (Goodfellow [Model-Data] as a function of ∆. Notice that DKL (pθ ||pdata )
is insensitive to placing weight in the model distribution in
et al., 2014). regions where pθ ≈ 0 whereas DKL (pdata ||pθ ) punishes this
This difference in objectives underlies the difference in harshly .
behavior between GANs and likelihood based generative
models. To see this, we can compare the behavior of the
two KL-divergences DKL (pdata ||pθ ) and DKL (pθ ||pdata ).
As is evident from Fig. 69 and Fig. 70, though both of to placing weight in the model distribution in regions
these KL-divergences measure similarities between the where pdata ≈ 0 whereas DKL (pθ ||pdata ) punishes this
two distributions, they are sensitive to very different harshly. In other words, DKL (pdata ||pθ ) prefers models
things. DKL (pθ ||pdata ) is insensitive to setting pθ ≈ 0 that have a high probability in regions with lots of train-
even when pdata 6= 0 whereas DKL (pdata ||pθ ) punishes ing data points whereas DKL (pθ ||pdata ) punishes models
this harshly. In contrast, DKL (pdata ||pθ ) is insensitive for putting high probability where there is no data.
102
An alternative approach that avoids these sampling On the other hand, calculating the gradient with re-
issues is to use the variational approach discussed in spect to the parameters φ is more complicated since φ
Sec. XIV. We know from Eq. (162) that we can write also appears in the expectation value Eqφ (z|x) . Ideally, we
the log-likelihood as would like to also use backpropagation to calculate this
as well. It turns out that this can be done by a simple
log p(x) = DKL (qφ (z|x)kpθ (z|x, θ)) − Fqφ (x), (222) change of variables that often goes under the name the
where the variational free energy is defined as “reparameterization trick” (Kingma and Welling, 2013;
Rezende et al., 2014). The basic idea is to change vari-
− Fqφ (x) ≡ Eqφ (z|x) [log pθ (x, z)] − DKL (qφ (z|x)|p(z)). ables so that φ no longer appears in the distribution we
(223) are taking an expectation value with respect to. To do
In writing this term, we have used Bayes rule and this, we express the random variable z ∼ qφ (z|x) as some
Eq. (174). Since the KL-divergence is strictly positive, differentiable and invertible transformation of another
the (negative) variational free energy is a lower-bound on random variable :
the log-likelihood. For this reason, in the VAE literature,
it is often called the Evidence Lower BOund or ELBO. z = g(, φ, x), (226)
Equation (223) has a beautiful interpretation. The where the distribution of is independent of x and φ.
first term in this equation can be viewed as a “recon- Then, we can replace expectation values over qφ (z|x) by
struction error”, where we start with data x, encode it expectation values over p
into the latent representation using our approximate pos-
terior qφ (z|x), and then evaluate the log probability of Eqφ (z|x) [f (z)] = Ep [f (z)]. (227)
the original data given the inferred latents. For binary
variables, this is just the cross-entropy which we first en- Evaluating the derivative then becomes quite straight for-
countered when studying logistic regression, cf. Sec. VII. ward since
The second term acts as a regularizer and encourages the ∇φ Eqφ (z|x) [f (z)] ∼ Ep [∇φ f (z)]. (228)
posterior distributions to be close to p(z). By maximizing
the ELBO, we minimize the KL-divergence between the Of course, when we do this we still need to be able to
approximate and true posterior. By choosing a tractable calculate the Jacobian of this change of variables
qφ (z|x), we make this feasible (see Fig. 72).
∂z
dφ (x, φ) = Det (229)
∂
2. Training via the reparametrization trick
since
VAEs train models by minimizing the variational free log qφ (z|x) = log p() − log dφ (x, φ). (230)
energy (maximizing the ELBO). Training a VAE is some-
what complicated because we must simultaneously learn Since we can calculate gradients, we can now use back-
two sets of parameters: the parameters θ that define our propagation on the full the ELBO objective function (we
generative model pθ (x, z) as well as the variational pa- return to this below when we discuss concrete architec-
rameters φ in qφ (z|x). The basic approach will be the tures and implementations of VAE).
same as for all DNN models: we will use gradient de- One of the problems that commonly occurs when train-
scent with the variational free energy as the objective ing VAEs by performing a stochastic optimization of the
(cost) function. For a dataset L, we can write our cost ELBO (variational free energy) is that it often gets stuck
function as in undesirable local minima, especially at the beginning
X of the training procedure (Bowman et al., 2015; Kingma
Cθ,φ (L) = −Fqφ (x). (224) et al., 2017; Sønderby et al., 2016). The underlying rea-
x∈L son for this is that the ELBO objective function can be
Taking the gradient with respect to θ is easy since only improved in two qualitatively different ways correspond-
the first term in Eq. (223) depends on θ, ing to each of the two terms in Eq. (223): by minimizing
the reconstruction error or by making the posterior dis-
Cθ,φ (x) = Eqφ (z|x) [∇θ log pθ (x, z)] tribution qφ (z|x) to be close to p(z) (Of course, the goal
∼ ∇θ log pθ (x, z) (225) is to do both!). For complex datasets, at the beginning of
training when the reconstruction error is extremely poor,
where in the second line we have replaced the expec- the model often quickly learns to make q(z|x) ≈ p(z) and
tation value with a single Monte-Carlo sample z drawn gets stuck in this local minimum. For this reason, in prac-
from qφ (z|x) (see Fig. XVII.C.2). When pθ (x|z) is ap- tice it is found that it makes sense to modify the ELBO
proximated by a neural network, this can be calculated objective to use an optimization schedule of the form
using backpropagation with the reconstruction error as
the objective function. Eqφ (z|x) [log pθ (x, z)] − βDKL (qφ (z|x)|p(z)) (231)
105
Datapoint
Negative
Variational Free Energy:
(ELBO)
E [log pθ(x|z) - KL(qφ(z|x)||p(z))]
q (z|x)
φ
FIG. 73 Schematic explaining the computational flow of VAEs. Figure based on Kingma’s Ph.D. dissertation Chapter 2.
(Kingma et al., 2017).
where β is slowly annealed from 0 to 1 (Bowman et al., implementing IB is that it requires knowledge of the joint
2015; Sønderby et al., 2016). An alternative regulariza- distribution p(x, y) and that we must be able to compute
tion is the “method of free bits”: modifying the objective the mutual information, a notoriously difficult quantity
function of ELBO to ensure that on average qφ (z|x) has to estimate from samples. Hence, IB has in recent years
at least λ natural units of information about p(z) (see been utilized less than it might otherwise.
Kingma Ph.D thesis (Kingma et al., 2017) for details) . To address these problems, variational approximations
These observations hints at the more general connec- to the IB objective function have been developed (Alemi
tion between VAEs and information theory that we turn et al., 2016; Chalk et al., 2016). These approximations,
to in the next section. when applied to a particular choice of p(x, y) give the
same objective as the variational autoencoder. Here
we follow the exposition from Alemi et al.(Alemi et al.,
3. Connection to the information bottleneck 2016). To see this, consider a dataset of N points, xi .
We set x = i and y = xi in the IB objective, similar
There is a fundamental connection between the vari- to (Slonim et al., 2005; Strouse and Schwab, 2017). We
ational autoencoder objective and the information bot- choose p(i) = 1/N and p(x|i) = δ(x − xi ). That is, we
tleneck (IB) for lossy compression (Tishby et al., 2000). would like to find a compression of the data that preserves
The information bottleneck imagines we have input data information about data point location while reducing in-
x that is correlated with another variable of interest, y, formation about data point identity.
and we are given access to the joint distribution, p(x, y). Imagine that we are unable to directly work with the
Our task is to take x as input and compress it in such a decoder p(x|z). The first approximation replaces the ex-
way as to retain as much information as possible about act decoder inside the logarithm with an approximation,
the relevance variable, y. To do this, Tishby et al. pro- q(x|z). Due to the positivity of KL-divergence, namely,
pose to maximize the objective function
DKL (p(x|z)||q(x|z)) ≥ 0
Z Z
⇒ dx p(x|z) log p(x|z) ≥ dx p(x|z) log q(x|z),
LIB = I(y; z) − βI(x; z) (232)
(233)
over a stochastic encoding distribution q(z|x), where z is
our compression of the input, and β is a tradeoff param- we have
eter that sets the relative preference of compression and
p(x|z)
Z
accuracy, and I(y; z) is the mutual information between y I(x; z) = dxdz p(x)p(z|x) log
p(x)
and z. Note that we choose a slightly different but equiv- Z
alent form of the objective relative to Tishby et al.. This ≥ dxdz p(x)p(z|x) log q(x|z) + Hp (x)
objective is only known to have a closed-form solution Z
when x and y are jointly Gaussian (Chechik et al., 2005). ≥ dxdz p(x)p(z|x) log q(x|z), (234)
Otherwise, the optimization can be performed through a
Blahut-Arimoto type iterative update scheme (Arimoto, where Hp (x) ≥ 0 is the Shannon entropy of x. This
1972; Blahut, 1972). However, this is only guaranteed to quantity can be estimated from data samples (i, xi ) af-
converge to a local optimum. A significant difficulty in ter drawing from p(z|i) = p(z|xi ). Similarly, we can
106
a Monte-Carlo method and HMC was developed physics interpretations in terms of transfer matri-
by physicists and exploits Hamiltonian dynamics ces (see (Mehta et al., 2011) for explicit exam-
to improve proposal distributions. ple of this). Recently, Recurrent Neural Networks
(RNNs) have become an important and powerful
• Statistical physics style theory of Deep tool for dealing with sequence data (Goodfellow
Learning. Many techniques in ML have origins et al., 2016). RNNs generalize many of the ideas
in statistical physics. Yet, a physics-style theory discussed in the DNN section to deal with temporal
of Deep Learning remains elusive. A key question data.
is to ask when and why these models manage to
generalize well. Physicists are only beginning to
ask these questions (Advani and Saxe, 2017; Mehta
• Reinforcement Learning. Many of the most ex-
and Schwab, 2014; Saxe et al., 2013; Shwartz-Ziv
citing developments in the last five years have come
and Tishby, 2017). But right now, it is fair to say
from combining ideas from reinforcement learning
that the insights remain scattered and a cohesive
with deep neural networks (Mnih et al., 2015; Sut-
theoretical understanding is lacking.
ton and Barto, 1998). RL traces its origins to be-
• Biological physics and ML. Biological physics haviourist psychology, when it was conceived as
is generating ever more datasets in fields ranging a way to explain and study reward-based deci-
from neuroscience to evolution and immunology. It sion making. RL was put on solid mathematical
is likely that ML will be an important part of the grounds in the 50’s by Richard Bellman and col-
biophysics toolkit in the future. Many of the au- laborators, and has by now become an insepara-
thors of this review were inspired to engage with ble part of robotics and artificial intelligence. RL
ML for this reason. is a field of Machine Learning, in which an agent
learns how to master performing a specific task
• Using ideas from physics to develop new through an interaction with its environment. De-
ML algorithms. Many of the core ideas of ML pending on the reward it receives, the agent chooses
from Monte-Carlo techniques to variational meth- to take an action affecting the environment, which
ods have their origin in physics. There has been a in turn determines the value of the next received
tremendous amount of recent work developing tools reward, and so on. The long-term goal of the
to understand physical systems that may be of po- agent is to maximise the cumulative expected re-
tential use to ML. For example, in quantum con- turn, thus improving its performance in the longer
densed matter techniques such as DMRG, MERA, run. Shadowed by more traditional optimal con-
etc. have enriched both our practical and con- trol algorithms, Reinforcement Learning has only
ceptual understandings (Stoudenmire and White, recently taken off in physics (Albarran-Arriagada
2012; Vidal, 2007; White, 1992). It will be inter- et al., 2018; August and Hernández-Lobato, 2018;
esting to figure how and if these numerical methods Bukov, 2018; Bukov et al., 2018; Cárdenas-López
can be translated from a physics to a ML setting. et al., 2017; Chen et al., 2014; Chen and Xue, 2019;
There are tantalizing hints that this is likely to be Dunjko et al., 2017; Fösel et al., 2018; Lamata,
a fruitful direction (Han et al., 2017; Stoudenmire, 2017; Melnikov et al., 2017; Neukart et al., 2017;
2018; Stoudenmire and Schwab, 2016). Niu et al., 2018; Ramezanpour, 2017; Reddy et al.,
2016b; Sriarunothai et al., 2017; Zhang et al., 2018).
Of particular interest are biophysics inspired works
B. Topics not covered in review
that seek to use RL to understand navigation and
sensing in turbulent environments (Colabrese et al.,
Despite the considerable length of the review, we have
2017; Masson et al., 2009; Reddy et al., 2016a; Ver-
had to make many omissions for the sake of brevity. It
gassola et al., 2007).
is our hope and belief that after reading this review the
reader will have the conceptual and practical knowledge
to quickly learn about these other topics. Among the
most prominent topics missing from this review are: • Support Vector Machines (SVMs) and Ker-
nel Methods. SVMs and kernel methods are a
• Temporal/Sequential Data. We have not cov- powerful set of techniques that work well when the
ered techniques for dealing with temporal or se- amount of training data is limited (Burges, 1998).
quential data. Here, too there are many connec- The mathematics and theory of SVM are very dif-
tions with statistical physics. A powerful class of ferent from statistical physics and for this reason
models for sequential data called Hidden Markov we chose not to include them here. However, SVMs
Models (Rabiner, 1989) that utilize dynamical and kernel methods have played an extremely im-
programming techniques have natural statistical portant role in ML and are worth understanding.
111
C. Rebranding Machine Learning as “Artificial Intelligence” ity distributions, and (c) finding policies for tasks with
well-defined goals and clear rules. We hope that this
The immense scientific progress in ML has also been review has given the reader enough conceptual tools to
accompanied by a massive public relations effort cen- start forming their own opinions about reality and hype
tered around Silicon Valley. Starting with the success when it comes to modern ML research. As Michael I. Jo-
of ImageNet (the most prominent early use of GPUs ran puts it, “...if the acronym âĂIJAIâĂİ continues to be
for training large models) and the widespread adoption used as placeholder nomenclature going forward, letâĂŹs
of Deep Learning based techniques by the Silicon Val- be aware of the very real limitations of this placeholder.
ley companies, there has been a deliberate move to re- LetâĂŹs broaden our scope, tone down the hype and rec-
brand modern ML as “artificial intelligence” or AI (see ognize the serious challenges ahead "(Jordan, 2018).
graphs in (Katz, 2017)). Recently, computer scientist
Michael I. Jordan (who is famously known for his formal-
ization of variational inference, Bayesian network, and
expectation-maximization algorithm in machine learning D. Social Implications of Machine Learning
research) cautioned that “This confluence of ideas and
technology trends has been rebranded as âĂIJAIâĂİ over The last decade has also seen a systematic increase in
the past few years. This rebranding is worthy of some the use and deployment of Machine Learning techniques
scrutiny”(Jordan, 2018). into new areas of life and society. Some of the readers of
AI, by design, is an ambiguous term that mixes aspi- this review may currently be (or eventually be) employed
rations with reality. It also conflates the statistical ideas in industrial settings that seek to harness ML for practi-
that form the basis of modern ML with the more com- cal purposes. However, caution is in order when applying
monplace notions about what humans and behavioral sci- ML. Without foresight and accountability, the scale and
entists mean by intelligence (see (Lake et al., 2017) for scope of modern ML algorithms can lead to large scale
an enlightening and important modern discussion of this unaccountable and undemocratic outcomes that can re-
distinction from a quantitative cognitive science point of inforce or even worsen existing inequality and inequities.
view as well as (Dreyfus, 1965) for a surprisingly relevant Mathematician and data scientist turned social commen-
philosophy-based critique from 1965). tator Cathy O’Neil has dubbed the indiscriminate use of
Almost all the techniques discussed here rely on op- these Big Data techniques “Weapons of Math Destruc-
timizing a pre-specified objective function on a given tion” (O’Neil, 2017).
dataset. Yet, we know that for large, complex models When ML is used in a social context, abstract statis-
changing the data distribution or the goal can lead to an tical relationships have real social consequences. False
immediate degradation of performance. Deep networks positives can mean the difference between life and death
have poor generalizations to even a slightly different con- (for example in the context of “signature drone strikes”)
text (the infamous Validation-Test set mismatch). This (Mehta, 2015). ML algorithms, like all techniques, have
inability to abstract and generalize is a common criticism important limitations and should be employed with great
lobbied against branding modern ML techniques as AI caution. It is our hope that ML practitioners keep this
(Lake et al., 2017). For all these reasons, we have chosen in mind when working in social settings.
to use the term Machine Learning rather than artificial All algorithms involve inherent tradeoffs in fairness, a
intelligence through out the review. point formalized by computer scientist Jon Kleinberg and
This is far from the first time we have seen the use of collaborators in a very interesting recent paper (Klein-
the term artificial intelligence and the grandiose promises berg et al., 2016). It is far from clear how to make al-
that it implies. In fact, the early 1950’s and 1960’s as gorithms fair for all people involved. This is even more
well as the early 1980’s saw similar AI bubbles (see this true with methods like Deep Learning that are hard to
interesting summary by Luke Muehlhauser for Open Phi- interpret. All ML algorithms have implicit assumptions
lanthropy (Muehlhauser, 2016)). These AI bubbles have and choices reflected in the datasets we use to the kind
been followed by what have been dubbed “AI Winters” of functions we choose to optimize. It is important to
(McDermott et al., 1985). remember that there is no “ view from nowhere” (Adam,
The “Singularity” may not be coming but the advances 2006; Katz, 2017) – all ML algorithms reflect a point of
in computing and the availability of large data sets likely view and a set of assumptions about the world we live
ensure that the kind of statistical learning frameworks in. For this reason, we hope that ML practitioners and
discussed are here to stay. Rather than a general arti- data scientists will take the time to consider the social
ficial intelligence, the kind of techniques presented here consequences of their actions. For example, developing a
seem to be best suited for three important tasks: (a) Hippocratic Oath for data scientists is now being consid-
automating prediction from lots of labeled examples in ered (Simonite, 2018). Doing no harm seems like a good
a narrowly-defined setting (b) learning how to parame- start for making sure that we harness ML for the benefit
terize and capture the correlations of complex probabil- of all members of society.
112
3. MNIST Dataset Amit, Daniel J (1992), Modeling brain function: The world of
attractor neural networks (Cambridge university press).
The MNIST dataset is one of the simplest and most Amit, Daniel J, Hanoch Gutfreund, and Haim Sompolinsky
widely used Machine Learning Datasets. The MNIST (1985), “Spin-glass models of neural networks,” Physical
Review A 32 (2), 1007.
dataset consists of hand-written images of numerical
Andrieu, Christophe, Nando De Freitas, Arnaud Doucet, and
characters 0−9 and consists of a training set of 60,000 ex- Michael I Jordan (2003), “An introduction to mcmc for
amples, and a test set of 10,000 examples (LeCun et al., machine learning,” Machine learning 50 (1-2), 5–43.
1998a). Information about the MNIST database and its Arai, Shunta, Masayuki Ohzeki, and Kazuyuki Tanaka
historical importance can be found at Yann Lecun’s wed- (2017), “Deep neural network detects quantum phase tran-
site: [Link] A brief sition,” arXiv preprint arXiv:1712.00371 .
description from the website: Arimoto, Suguru (1972), “An algorithm for computing the
capacity of arbitrary discrete memoryless channels,” IEEE
The original black and white (bilevel) images Transactions on Information Theory 18 (1), 14–20.
from NIST were size normalized to fit in a Arsenault, Louis-François, Alejandro Lopez-Bezanilla,
20x20 pixel box while preserving their aspect O Anatole von Lilienfeld, and Andrew J Millis (2014),
ratio. The resulting images contain grey lev- “Machine learning for many-body physics: the case of the
anderson impurity model,” Physical Review B 90 (15),
els as a result of the anti-aliasing technique
155136.
used by the normalization algorithm. the im- Arunachalam, Srinivasan, and Ronald de Wolf (2017),
ages were centered in a 28x28 image by com- “A survey of quantum learning theory,” arXiv preprint
puting the center of mass of the pixels, and arXiv:1701.06806 .
translating the image so as to position this August, Moritz, and José Miguel Hernández-Lobato (2018),
point at the center of the 28x28 field. “Taking gradients through experiments: Lstms and mem-
ory proximal policy optimization for black-box quantum
The MNIST is often included by default in many modern control,” arXiv:1802.04063.
ML packages. Baireuther, P, TE O’Brien, B Tarasinski, and CWJ
Beenakker (2017), “Machine-learning-assisted correction
of correlated qubit errors in a topological code,” arXiv
REFERENCES preprint arXiv:1705.07855 .
Baity-Jesi, M, L. Sagun, M. Geiger, S. Spigler, G. Ben Arous,
Abu-Mostafa, Yaser S, Malik Magdon-Ismail, and Hsuan- C. Cammarota, Y. LeCun, M. Wyart, and G. Biroli (2018),
Tien Lin (2012), Learning from data, Vol. 4 (AMLBook “Comparing dynamics: Deep neural networks versus glassy
New York, NY, USA:). systems,” .
Ackley, David H, Geoffrey E Hinton, and Terrence J Se- Baldassi, Carlo, Federica Gerace, Hilbert J Kappen, Carlo
jnowski (1987), “A learning algorithm for boltzmann ma- Lucibello, Luca Saglietti, Enzo Tartaglione, and Riccardo
chines,” in Readings in Computer Vision (Elsevier) pp. Zecchina (2017), “On the role of synaptic stochasticity
522–533. in training low-precision neural networks,” arXiv preprint
Adam, Alison (2006), Artificial knowing: Gender and the arXiv:1710.09825 .
thinking machine (Routledge). Baldassi, Carlo, Federica Gerace, Luca Saglietti, and Ric-
Advani, Madhu, and Surya Ganguli (2016), “Statistical me- cardo Zecchina (2018), “From inverse problems to learning:
chanics of optimal convex inference in high dimensions,” a statistical mechanics approach,” in Journal of Physics:
Physical Review X 6 (3), 031034. Conference Series, Vol. 955 (IOP Publishing) p. 012001.
Advani, Madhu, Subhaneil Lahiri, and Surya Ganguli (2013), Baldi, Pierre, Peter Sadowski, and Daniel Whiteson (2014),
“Statistical mechanics of complex neural systems and high “Searching for exotic particles in high-energy physics with
dimensional data,” Journal of Statistical Mechanics: The- deep learning,” Nature communications 5, 4308.
ory and Experiment 2013 (03), P03014. Barber, David (2012), Bayesian reasoning and machine learn-
Advani, Madhu S, and Andrew M Saxe (2017), “High- ing (Cambridge University Press).
dimensional dynamics of generalization error in neural net- Barnes, Josh, and Piet Hut (1986), “A hierarchical o (n log
works,” arXiv preprint arXiv:1710.03667. n) force-calculation algorithm,” nature 324 (6096), 446.
Aitchison, Laurence, Nicola Corradi, and Peter E Latham Barra, Adriano, Alberto Bernacchia, Enrica Santucci, and
(2016), “Zipfs law arises naturally when there are under- Pierluigi Contucci (2012), “On the equivalence of hopfield
lying, unobserved variables,” PLoS computational biology networks and boltzmann machines,” Neural Networks 34,
12 (12), e1005110. 1–9.
Albarran-Arriagada, F, J. C. Retamal, E. Solano, and Barra, Adriano, Giuseppe Genovese, Peter Sollich, and
L. Lamata (2018), “Measurement-based adapta- Daniele Tantari (2016), “Phase transitions in restricted
tion protocol with quantum reinforcement learning,” boltzmann machines with generic priors,” arXiv preprint
arXiv:1803.05340 . arXiv:1612.03132 .
Alemi, Alexander A, Ian Fischer, Joshua V Dillon, and Kevin Barra, Adriano, Giuseppe Genovese, Daniele Tantari, and Pe-
Murphy (2016), “Deep variational information bottleneck,” ter Sollich (2017), “Phase diagram of restricted boltzmann
arXiv preprint arXiv:1612.00410. machines and generalised hopfield networks with arbitrary
Alemi, Alireza, and Alia Abbara (2017), “Exponential capac- priors,” arXiv preprint arXiv:1702.05882 .
ity in an autoencoder neural network with a hidden layer,” Battiti, Roberto (1992), “First-and second-order methods for
arXiv preprint arXiv:1705.07441 . learning: between steepest descent and newton’s method,”
114
Neural computation 4 (2), 141–166. ing the quantum kapitza oscillator,” Phys. Rev. B 98,
Benedetti, Marcello, John Realpe-Gómez, Rupak Biswas, 224305.
and Alejandro Perdomo-Ortiz (2016), “Quantum-assisted Bukov, Marin, Alexandre G. R. Day, Dries Sels, Phillip Wein-
learning of graphical models with arbitrary pairwise con- berg, Anatoli Polkovnikov, and Pankaj Mehta (2018), “Re-
nectivity,” arXiv preprint arXiv:1609.02542 . inforcement learning in different phases of quantum con-
Benedetti, Marcello, John Realpe-Gómez, and Alejandro trol,” Phys. Rev. X 8, 031086.
Perdomo-Ortiz (2017), “Quantum-assisted helmholtz ma- Burges, Christopher JC (1998), “A tutorial on support vector
chines: A quantum-classical deep learning framework for machines for pattern recognition,” Data mining and knowl-
industrial datasets in near-term devices,” arXiv preprint edge discovery 2 (2), 121–167.
arXiv:1708.09784 . Caio, Marcello D, Marco Caccin, Paul Baireuther, Timo
Bengio, Yoshua (2012), “Practical recommendations for Hyart, and Michel Fruchart (2019), “Machine learning as-
gradient-based training of deep architectures,” in Neural sisted measurement of local topological invariants,” arXiv
networks: Tricks of the trade (Springer) pp. 437–478. preprint arXiv:1901.03346 .
Bennett, Robert (1969), “The intrinsic dimensionality of sig- Caldeira, J, WLK Wu, B Nord, C Avestruz, S Trivedi, and
nal collections,” IEEE Transactions on Information Theory KT Story (2018), “Deepcmb: Lensing reconstruction of the
15 (5), 517–525. cosmic microwave background with deep neural networks,”
Bény, Cédric (2018), “Inferring relevant features: from qft to arXiv preprint arXiv:1810.01483 .
pca,” arXiv preprint arXiv:1802.05756 . Canabarro, Askery, Samuraí Brito, and Rafael Chaves (2018),
Berger, James O, and José M Bernardo (1992), “On the devel- “Machine learning non-local correlations,” arXiv preprint
opment of the reference prior method,” Bayesian statistics arXiv:1808.07069 .
4, 35–60. Cárdenas-López, FA, L Lamata, JC Retamal, and E Solano
Bickel, Peter J, and David A Freedman (1981), “Some asymp- (2017), “Generalized quantum reinforcement learning with
totic theory for the bootstrap,” The Annals of Statistics , quantum technologies,” arXiv preprint arXiv:1709.07848 .
1196–1217. Carleo, Giuseppe (2018), , Private Communication.
Bickel, Peter J, Bo Li, Alexandre B Tsybakov, Sara A van de Carleo, Giuseppe, Yusuke Nomura, and Masatoshi Imada
Geer, Bin Yu, Teófilo Valdés, Carlos Rivero, Jianqing Fan, (2018), “Constructing exact representations of quantum
and Aad van der Vaart (2006), “Regularization in statis- many-body systems with deep neural networks,” arXiv
tics,” Test 15 (2), 271–344. preprint arXiv:1802.09558 .
Bishop, C M (2006), Pattern recognition and machine learning Carleo, Giuseppe, and Matthias Troyer (2017), “Solving the
(springer). quantum many-body problem with artificial neural net-
Bishop, Chris M (1995a), “Training with noise is equivalent to works,” Science 355 (6325), 602–606.
tikhonov regularization,” Neural computation 7 (1), 108– Carrasquilla, Juan, and Roger G Melko (2017), “Machine
116. learning phases of matter,” Nature Physics 13 (5), 431.
Bishop, Christopher M (1995b), Neural networks for pattern Chalk, Matthew, Olivier Marre, and Gasper Tkacik (2016),
recognition (Oxford university press). “Relevant sparse codes with variational information bottle-
Blahut, Richard (1972), “Computation of channel capacity neck,” in Advances in Neural Information Processing Sys-
and rate-distortion functions,” IEEE transactions on Infor- tems, pp. 1957–1965.
mation Theory 18 (4), 460–473. Chamberland, Christopher, and Pooya Ronagh (2018), “Deep
Bottou, Léon (2012), “Stochastic gradient descent tricks,” in neural decoders for near term fault-tolerant experiments,”
Neural networks: Tricks of the trade (Springer) pp. 421– arXiv preprint arXiv:1802.06441 .
436. Chechik, Gal, Amir Globerson, Naftali Tishby, and Yair
Bowman, Samuel R, Luke Vilnis, Oriol Vinyals, Andrew M Weiss (2005), “Information bottleneck for gaussian vari-
Dai, Rafal Jozefowicz, and Samy Bengio (2015), “Gener- ables,” Journal of machine learning research 6 (Jan), 165–
ating sentences from a continuous space,” arXiv preprint 188.
arXiv:1511.06349. Chen, Chunlin, Daoyi Dong, Han-Xiong Li, Jian Chu, and
Boyd, Stephen, and Lieven Vandenberghe (2004), Convex Tzyh-Jong Tarn (2014), “Fidelity-based probabilistic q-
optimization (Cambridge university press). learning for control of quantum systems,” IEEE transac-
Bradde, Serena, and William Bialek (2017), “Pca meets rg,” tions on neural networks and learning systems 25 (5), 920–
Journal of Statistical Physics 167 (3-4), 462–475. 933.
Breiman, Leo (1996), “Bagging predictors,” Machine learning Chen, Jing, Song Cheng, Haidong Xie, Lei Wang, and Tao Xi-
24 (2), 123–140. ang (2018), “Equivalence of restricted boltzmann machines
Breiman, Leo (2001), “Random forests,” Machine learning and tensor network states,” Phys. Rev. B 97, 085104.
45 (1), 5–32. Chen, Jun-Jie, and Ming Xue (2019), “Manipulation of spin
Breuckmann, Nikolas P, and Xiaotong Ni (2017), “Scalable dynamics by deep reinforcement learning agent,” arXiv
neural network decoders for higher dimensional quantum preprint arXiv:1901.08748 .
codes,” arXiv preprint arXiv:1710.09489 . Chen, Tianqi, and Carlos Guestrin (2016), “Xgboost: A scal-
Broecker, Peter, Fakher F Assaad, and Simon Trebst (2017), able tree boosting system,” in Proceedings of the 22nd acm
“Quantum phase recognition via unsupervised machine sigkdd international conference on knowledge discovery and
learning,” arXiv preprint arXiv:1707.00663 . data mining (ACM) pp. 785–794.
Bromley, Thomas R, and Patrick Rebentrost (2018), Cheng, Song, Jing Chen, and Lei Wang (2017), “Information
“Batched quantum state exponentiation and quantum heb- perspective to probabilistic modeling: Boltzmann machines
bian learning,” arXiv:1803.07039 . versus born machines,” arXiv preprint arXiv:1712.04144 .
Bukov, Marin (2018), “Reinforcement learning for au- Ch’ng, Kelvin, Nick Vazquez, and Ehsan Khatami
tonomous preparation of floquet-engineered states: Invert- (2017), “Unsupervised machine learning account of mag-
115
netic transitions in the hubbard model,” arXiv preprint Dreyfus, Hubert L (1965), “Alchemy and artificial intelli-
arXiv:1708.03350 . gence,”.
Ciliberto, Carlo, Mark Herbster, Alessandro Davide Ialongo, Du, Simon S, Chi Jin, Jason D Lee, Michael I Jordan, Aarti
Massimiliano Pontil, Andrea Rocchetto, Simone Severini, Singh, and Barnabas Poczos (2017), “Gradient descent can
and Leonard Wossnig (2017), “Quantum machine learning: take exponential time to escape saddle points,” in Advances
a classical perspective,” . in Neural Information Processing Systems, pp. 1067–1077.
Cohen, Nadav, Or Sharir, and Amnon Shashua (2016), “On Duchi, John, Elad Hazan, and Yoram Singer (2011), “Adap-
the expressive power of deep learning: A tensor analysis,” tive subgradient methods for online learning and stochas-
in Conference on Learning Theory, pp. 698–728. tic optimization,” Journal of Machine Learning Research
Colabrese, Simona, Kristian Gustavsson, Antonio Celani, 12 (Jul), 2121–2159.
and Luca Biferale (2017), “Flow navigation by smart mi- Dunjko, Vedran, and Hans J Briegel (2017), “Machine learn-
croswimmers via reinforcement learning,” Physical review ing and artificial intelligence in the quantum domain,”
letters 118 (15), 158004. arXiv preprint arXiv:1709.02779 .
Cossu, Guido, Luigi Del Debbio, Tommaso Giani, Ava Kham- Dunjko, Vedran, Yi-Kai Liu, Xingyao Wu, and Jacob M
seh, and Michael Wilson (2018), “Machine learning deter- Taylor (2017), “Super-polynomial and exponential im-
mination of dynamical parameters: The ising model case,” provements for quantum-enhanced reinforcement learning,”
arXiv preprint arXiv:1810.11503 . arXiv preprint arXiv:1710.11160 .
Cox, Trevor F, and Michael AA Cox (2000), Multidimen- Efron, B (1979), “Bootstrap methods: another look at the
sional scaling (CRC press). jackknife annals of statistics 7: 1–26,” View Article Pub-
Cristoforetti, Marco, Giuseppe Jurman, Andrea I Nardelli, Med/NCBI Google Scholar.
and Cesare Furlanello (2017), “Towards meaningful physics Efron, Bradley, Trevor Hastie, Iain Johnstone, Robert Tib-
from generative models,” arXiv preprint arXiv:1705.09524 shirani, et al. (2004), “Least angle regression,” The Annals
. of statistics 32 (2), 407–499.
Dahl, George, Abdel-rahman Mohamed, Geoffrey E Hinton, Eisen, Michael B, Paul T Spellman, Patrick O Brown, and
et al. (2010), “Phone recognition with the mean-covariance David Botstein (1998), “Cluster analysis and display of
restricted boltzmann machine,” in Advances in neural in- genome-wide expression patterns,” Proceedings of the Na-
formation processing systems, pp. 469–477. tional Academy of Sciences 95 (25), 14863–14868.
Daskin, Ammar (2018), “A quantum implementation model Elith, Jane, Steven J Phillips, Trevor Hastie, Miroslav Dudík,
for artificial neural networks,” Quanta , 7–18. Yung En Chee, and Colin J Yates (2011), “A statistical
Davaasuren, Amarsanaa, Yasunari Suzuki, Keisuke Fujii, explanation of maxent for ecologists,” Diversity and distri-
and Masato Koashi (2018), “General framework for con- butions 17 (1), 43–57.
structing fast and near-optimal machine-learning-based de- Ernst, Oliver K, Thomas Bartol, Terrence Sejnowski, and
coder of the topological stabilizer codes,” arXiv preprint Eric Mjolsness (2018), “Learning dynamic boltzmann dis-
arXiv:1801.04377 . tributions as reduced models of spatial chemical kinetics,”
Day, Alexandre GR, Marin Bukov, Phillip Weinberg, Pankaj arXiv preprint arXiv:1803.01063 .
Mehta, and Dries Sels (2019), “Glassy phase of opti- Ester, Martin, Hans-Peter Kriegel, Jörg Sander, Xiaowei Xu,
mal quantum control,” Physical Review Letters 122 (2), et al. (1996), “A density-based algorithm for discovering
020601. clusters in large spatial databases with noise.” in Kdd,
Day, Alexandre GR, and Pankaj Mehta (2018), “Validated Vol. 96, pp. 226–231.
agglomerative clustering,” in preparation. Faddeev, Ludvig D, and Victor N Popov (1967), “Feynman
Decelle, Aurélien, Giancarlo Fissore, and Cyril Furtlehner diagrams for the yang-mills field,” Physics Letters B 25 (1),
(2017), “Spectral learning of restricted boltzmann ma- 29–30.
chines,” arXiv preprint arXiv:1708.02917 . Finol, David, Yan Lu, Vijay Mahadevan, and Ankit
Decelle, Aurélien, Giancarlo Fissore, and Cyril Furtlehner Srivastava (2018), “Deep convolutional neural networks
(2018), “Thermodynamics of restricted boltzmann ma- for eigenvalue problems in mechanics,” arXiv preprint
chines and related learning dynamics,” arXiv preprint arXiv:1801.05733 .
arXiv:1803.01960 . Fisher, Charles K, and Pankaj Mehta (2015a), “Bayesian fea-
Dempster, Arthur P, Nan M Laird, and Donald B Rubin ture selection for high-dimensional linear regression via the
(1977), “Maximum likelihood from incomplete data via the ising approximation with applications to genomics,” Bioin-
em algorithm,” Journal of the royal statistical society. Series formatics 31 (11), 1754–1761.
B (methodological) , 1–38. Fisher, Charles K, and Pankaj Mehta (2015b), “Bayesian
Deng, Dong-Ling, Xiaopeng Li, and S Das Sarma (2017), feature selection with strongly regularizing priors maps to
“Quantum entanglement in neural network states,” Physi- the ising model,” Neural computation 27 (11), 2411–2422.
cal Review X 7 (2), 021021. Foreman, Sam, Joel Giedt, Yannick Meurice, and Judah
Dietterich, Thomas G, et al. (2000), “Ensemble methods in Unmuth-Yockey (2017), “Rg inspired machine learning for
machine learning,” Multiple classifier systems 1857, 1–15. lattice field theory,” arXiv preprint arXiv:1710.02079 .
Do, Chuong B, and Serafim Batzoglou (2008), “What is the Fösel, Thomas, Petru Tighineanu, Talitha Weiss, and Flo-
expectation maximization algorithm?” Nature biotechnol- rian Marquardt (2018), “Reinforcement learning with neu-
ogy 26 (8), 897–899. ral networks for quantum feedback,” arXiv:1802.05267 .
Domingos, Pedro (2012), “A few useful things to know about Freitas, Nahuel, Giovanna Morigi, and Vedran Dun-
machine learning,” Communications of the ACM 55 (10), jko (2018), “Neural network operations and susuki-
78–87. trotter evolution of neural network states,” arXiv preprint
Donoho, David L (2006), “Compressed sensing,” IEEE Trans- arXiv:1803.02118 .
actions on information theory 52 (4), 1289–1306. Freund, Yoav, Robert Schapire, and Naoki Abe (1999), “A
116
short introduction to boosting,” Journal-Japanese Society Zhang (2017), “Unsupervised generative modeling using
For Artificial Intelligence 14 (771-780), 1612. matrix product states,” arXiv preprint arXiv:1709.01662 .
Freund, Yoav, and Robert E Schapire (1995), “A desicion- He, Kaiming, Xiangyu Zhang, Shaoqing Ren, and Jian Sun
theoretic generalization of on-line learning and an applica- (2015), “Delving deep into rectifiers: Surpassing human-
tion to boosting,” in European conference on computational level performance on imagenet classification,” in Proceed-
learning theory (Springer) pp. 23–37. ings of the IEEE international conference on computer vi-
Friedman, Jerome, Trevor Hastie, and Robert Tibshirani sion, pp. 1026–1034.
(2001), The elements of statistical learning, Vol. 1 (Springer He, Kaiming, Xiangyu Zhang, Shaoqing Ren, and Jian Sun
series in statistics New York). (2016), “Deep residual learning for image recognition,” in
Friedman, Jerome H (2001), “Greedy function approximation: Proceedings of the IEEE conference on computer vision and
a gradient boosting machine,” Annals of statistics , 1189– pattern recognition, pp. 770–778.
1232. Heimel, Theo, Gregor Kasieczka, Tilman Plehn, and Jen-
Friedman, Jerome H (2002), “Stochastic gradient boosting,” nifer M Thompson (2018), “Qcd or what?” arXiv preprint
Computational Statistics & Data Analysis 38 (4), 367–378. arXiv:1808.08979 .
Friedman, Jerome H, Bogdan E Popescu, et al. (2003), “Im- Higgins, Irina, Loic Matthey, Arka Pal, Christopher Burgess,
portance sampled learning ensembles,” Journal of Machine Xavier Glorot, Matthew Botvinick, Shakir Mohamed, and
Learning Research 94305. Alexander Lerchner (2016), “beta-vae: Learning basic vi-
Fu, Michael C (2006), “Gradient estimation,” Handbooks in sual concepts with a constrained variational framework,”
operations research and management science 13, 575–616. .
Funai, Shotaro Shiba, and Dimitrios Giataganas (2018), Hinton, Geoffrey E (2002), “Training products of experts by
“Thermodynamics and feature extraction by machine learn- minimizing contrastive divergence,” Neural computation
ing,” arXiv preprint arXiv:1810.08179 . 14 (8), 1771–1800.
Gao, Jun, Lu-Feng Qiao, Zhi-Qiang Jiao, Yue-Chi Ma, Cheng- Hinton, Geoffrey E (2012), “A practical guide to training re-
Qiu Hu, Ruo-Jing Ren, Ai-Lin Yang, Hao Tang, Man-Hong stricted boltzmann machines,” in Neural networks: Tricks
Yung, and Xian-Min Jin (2017), “Experimental machine of the trade (Springer) pp. 599–619.
learning of quantum states with partial information,” arXiv Hinton, Geoffrey E, Simon Osindero, and Yee-Whye Teh
preprint arXiv:1712.00456 . (2006), “A fast learning algorithm for deep belief nets,”
Gao, Xun, and Lu-Ming Duan (2017), “Efficient representa- Neural computation 18 (7), 1527–1554.
tion of quantum many-body states with deep neural net- Hinton, Geoffrey E, and Ruslan R Salakhutdinov (2006), “Re-
works,” arXiv preprint arXiv:1701.05039 . ducing the dimensionality of data with neural networks,”
Gelman, Andrew, John B Carlin, Hal S Stern, David B Dun- science 313 (5786), 504–507.
son, Aki Vehtari, and Donald B Rubin (2014), Bayesian Ho, Tin Kam (1998), “The random subspace method for con-
data analysis, Vol. 2 (CRC press Boca Raton, FL). structing decision forests,” IEEE transactions on pattern
Gersho, Allen, and Robert M Gray (2012), Vector quantiza- analysis and machine intelligence 20 (8), 832–844.
tion and signal compression, Vol. 159 (Springer Science & Hopfield, John J (1982), “Neural networks and physical sys-
Business Media). tems with emergent collective computational abilities,”
Geurts, Pierre, Damien Ernst, and Louis Wehenkel (2006), Proceedings of the national academy of sciences 79 (8),
“Extremely randomized trees,” Machine learning 63 (1), 2554–2558.
3–42. Huang, Haiping (2017a), “Mean-field theory of input dimen-
Glorot, Xavier, and Yoshua Bengio (2010), “Understanding sionality reduction in unsupervised deep neural networks,”
the difficulty of training deep feedforward neural networks,” arXiv preprint arXiv:1710.01467 .
in Proceedings of the Thirteenth International Conference Huang, Haiping (2017b), “Statistical mechanics of unsuper-
on Artificial Intelligence and Statistics, pp. 249–256. vised feature learning in a restricted boltzmann machine
Goldt, Sebastian, and Udo Seifert (2017), “Thermodynamic with binary synapses,” Journal of Statistical Mechanics:
efficiency of learning a rule in neural networks,” arXiv Theory and Experiment 2017 (5), 053302.
preprint arXiv:1706.09713 . Huang, Hengfeng, Bowen Xiao, Huixin Xiong, Zeming Wu,
Goodfellow, Ian (2016), “Nips 2016 tutorial: Generative ad- Yadong Mu, and Huichao Song (2018), “Applications
versarial networks,” arXiv preprint arXiv:1701.00160. of deep learning to relativistic hydrodynamics,” arXiv
Goodfellow, Ian, Yoshua Bengio, and Aaron preprint arXiv:1801.03334 .
Courville (2016), Deep Learning (MIT Press) Hubbard, J (1959), “Calculation of partition functions,” Phys-
[Link] ical Review Letters 3 (2), 77.
Goodfellow, Ian, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, Huggins, William, Piyush Patil, Bradley Mitchell, K Birgitta
David Warde-Farley, Sherjil Ozair, Aaron Courville, and Whaley, and E Miles Stoudenmire (2019), “Towards quan-
Yoshua Bengio (2014), “Generative adversarial nets,” in Ad- tum machine learning with tensor networks,” Quantum Sci-
vances in neural information processing systems, pp. 2672– ence and Technology 4 (2), 024001.
2680. Iakovlev, I A, O. M. Sotnikov, and V. V. Mazurenko
Greplova, Eliska, Christian Kraglund Andersen, and Klaus (2018), “Supervised learning magnetic skyrmion phases,”
Mølmer (2017), “Quantum parameter estimation with a arXiv:1803.06682 .
neural network,” arXiv preprint arXiv:1711.05238 . Innocenti, Luca, Leonardo Banchi, Alessandro Ferraro,
Grisafi, Andrea, David M Wilkins, Gábor Csányi, and Sougato Bose, and Mauro Paternostro (2018), “Supervised
Michele Ceriotti (2017), “Symmetry-adapted machine- learning of time-independent hamiltonians for gate design,”
learning for tensorial properties of atomistic systems,” arXiv:1803.07119 .
arXiv preprint arXiv:1709.06757 . Ioffe, Sergey, and Christian Szegedy (2015), “Batch normal-
Han, Zhao-Yu, Jun Wang, Heng Fan, Lei Wang, and Pan ization: Accelerating deep network training by reducing in-
117
ternal covariate shift,” in International Conference on Ma- Kohonen, Teuvo (1998), “The self-organizing map,” Neuro-
chine Learning, pp. 448–456. computing 21 (1-3), 1–6.
Iso, Satoshi, Shotaro Shiba, and Sumito Yokoo (2018), Kong, Qingkai, Daniel T. Trugman, Zachary E. Ross,
“Scale-invariant feature extraction of neural net- Michael J. Bianco, Brendan J. Meade, and Peter Ger-
work and renormalization group flow,” arXiv preprint stoft (2018), “Machine learning in seismology: Turn-
arXiv:1801.07172 . ing data into insights,” Seismological Research Letters
Jarzynski, Christopher (1997), “Nonequilibrium equality for 10.1785/0220180259.
free energy differences,” Physical Review Letters 78 (14), Krastanov, Stefan, and Liang Jiang (2017), “Deep neural
2690. network probabilistic decoder for stabilizer codes,” arXiv
Jaynes, Edwin T (1957a), “Information theory and statistical preprint arXiv:1705.09334 .
mechanics,” Physical review 106 (4), 620. Kriegel, Hans-Peter, Peer Kröger, and Arthur Zimek (2009),
Jaynes, Edwin T (1957b), “Information theory and statistical “Clustering high-dimensional data: A survey on subspace
mechanics. ii,” Physical review 108 (2), 171. clustering, pattern-based clustering, and correlation clus-
Jaynes, Edwin T (1968), “Prior probabilities,” IEEE Trans- tering,” ACM Transactions on Knowledge Discovery from
actions on systems science and cybernetics 4 (3), 227–241. Data (TKDD) 3 (1), 1.
Jaynes, Edwin T (1996), Probability theory: the logic of sci- Krizhevsky, Alex, Ilya Sutskever, and Geoffrey E Hinton
ence (Washington University St. Louis, MO). (2012), “Imagenet classification with deep convolutional
Jaynes, Edwin T (2003), Probability theory: the logic of sci- neural networks,” in Advances in neural information pro-
ence (Cambridge university press). cessing systems, pp. 1097–1105.
Jeffreys, Harold (1946), “An invariant form for the prior prob- Krzakala, Florent, Andre Manoel, Eric W Tramel, and Lenka
ability in estimation problems,” Proceedings of the Royal Zdeborová (2014), “Variational free energies for compressed
Society of London. Series A, Mathematical and Physical sensing,” in Information Theory (ISIT), 2014 IEEE Inter-
Sciences , 453–461. national Symposium on (IEEE) pp. 1499–1503.
Jin, Chi, Praneeth Netrapalli, and Michael I Jordan (2017), Krzakala, Florent, Marc Mézard, François Sausset, YF Sun,
“Accelerated gradient descent escapes saddle points faster and Lenka Zdeborová (2012a), “Statistical-physics-based
than gradient descent,” arXiv preprint arXiv:1711.10456. reconstruction in compressed sensing,” Physical Review X
Jordan, Michael (2018), “Artificial intelligence: The revolu- 2 (2), 021005.
tion hasnâĂŹt happened yet. medium,”. Krzakala, Florent, Marc Mézard, Francois Sausset, Yifan Sun,
Jordan, Michael I, Zoubin Ghahramani, Tommi S Jaakkola, and Lenka Zdeborová (2012b), “Probabilistic reconstruc-
and Lawrence K Saul (1999), “An introduction to vari- tion in compressed sensing: algorithms, phase diagrams,
ational methods for graphical models,” Machine learning and threshold achieving matrices,” Journal of Statistical
37 (2), 183–233. Mechanics: Theory and Experiment 2012 (08), P08009.
Kalantre, Sandesh S, Justyna P Zwolak, Stephen Ragole, Lake, Brenden M, Tomer D Ullman, Joshua B Tenenbaum,
Xingyao Wu, Neil M Zimmerman, MD Stewart, and Ja- and Samuel J Gershman (2017), “Building machines that
cob M Taylor (2017), “Machine learning techniques for learn and think like people,” Behavioral and Brain Sciences
state recognition and auto-tuning in quantum dots,” arXiv 40.
preprint arXiv:1712.04914 . Lamata, Lucas (2017), “Basic protocols in quantum reinforce-
Katz, Yarden (2017), “Manufacturing an artificial intelligence ment learning with superconducting circuits,” Scientific Re-
revolution,” SSRN Preprint . ports 7.
Kerenidis, Iordanis, and Anupam Prakash (2017), “Quan- Larsen, Bjornar, and Chinatsu Aone (1999), “Fast and effec-
tum gradient descent for linear systems and least squares,” tive text mining using linear-time document clustering,” in
arXiv preprint arXiv:1704.04992 . Proceedings of the fifth ACM SIGKDD international con-
Keskar, Nitish Shirish, Dheevatsa Mudigere, Jorge Nocedal, ference on Knowledge discovery and data mining (ACM)
Mikhail Smelyanskiy, and Ping Tak Peter Tang (2016), pp. 16–22.
“On large-batch training for deep learning: Generalization Le, Quoc V (2013), “Building high-level features using large
gap and sharp minima,” arXiv preprint arXiv:1609.04836. scale unsupervised learning,” in Acoustics, Speech and Sig-
Kingma, Diederik P, and Jimmy Ba (2014), “Adam: nal Processing (ICASSP), 2013 IEEE International Con-
A method for stochastic optimization,” arXiv preprint ference on (IEEE) pp. 8595–8598.
arXiv:1412.6980. LeCun, Yann, Yoshua Bengio, et al. (1995), “Convolutional
Kingma, Diederik P, and Max Welling (2013), networks for images, speech, and time series,” The hand-
“Auto-encoding variational bayes,” arXiv preprint book of brain theory and neural networks 3361 (10), 1995.
arXiv:1312.6114. LeCun, Yann, Léon Bottou, Yoshua Bengio, and Patrick
Kingma, DP, et al. (2017), “Variational inference & deep Haffner (1998a), “Gradient-based learning applied to docu-
learning,” PhD thesis 978-94-6299-745-5. ment recognition,” Proceedings of the IEEE 86 (11), 2278–
Kleijnen, Jack PC, and Reuven Y Rubinstein (1996), “Op- 2324.
timization and sensitivity analysis of computer simulation LeCun, Yann, Léon Bottou, Genevieve B Orr, and Klaus-
models by the score function method,” European Journal Robert Müller (1998b), “Efficient backprop,” in Neural net-
of Operational Research 88 (3), 413–427. works: Tricks of the trade (Springer) pp. 9–50.
Kleinberg, Jon, Sendhil Mullainathan, and Manish Raghavan Lee, Jason D, Ioannis Panageas, Georgios Piliouras, Max Sim-
(2016), “Inherent trade-offs in the fair determination of risk chowitz, Michael I Jordan, and Benjamin Recht (2017),
scores,” arXiv preprint arXiv:1609.05807. “First-order methods almost always avoid saddle points,”
Koch-Janusz, Maciej, and Zohar Ringel (2017), “Mutual in- arXiv preprint arXiv:1710.07406.
formation, neural networks and the renormalization group,” Lehmann, Erich L, and George Casella (2006), Theory of
arXiv preprint arXiv:1704.06279 . point estimation (Springer Science & Business Media).
118
Neukart, Florian, David Von Dollen, Christian Seidel, and Ramezanpour, A (2017), “Optimization by a quantum rein-
Gabriele Compostella (2017), “Quantum-enhanced rein- forcement algorithm,” Phys. Rev. A 96, 052307.
forcement learning for finite-episode games with discrete Rao, Wen-Jia, Zhenyu Li, Qiong Zhu, Mingxing Luo, and
state spaces,” arXiv preprint arXiv:1708.09354 . Xin Wan (2017), “Identifying product order with restricted
Nguyen, H Chau, Riccardo Zecchina, and Johannes Berg boltzmann machines,” arXiv preprint arXiv:1709.02597 .
(2017), “Inverse statistical problems: from the inverse ising Ravanbakhsh, Siamak, Francois Lanusse, Rachel Mandel-
problem to data science,” Advances in Physics 66 (3), 197– baum, Jeff G Schneider, and Barnabas Poczos (2017), “En-
261. abling dark energy science with deep generative models of
Nielsen, Michael A (2015), Neural networks and deep learning galaxy images.” in AAAI, pp. 1488–1494.
(Determination Press). Rebentrost, Patrick, Thomas R. Bromley, Christian Weed-
van Nieuwenburg, Evert, Eyal Bairey, and Gil Refael (2017a), brook, and Seth Lloyd (2017), “A quantum hopfield neural
“Learning phase transitions from dynamics,” arXiv preprint network,” arXiv:1710.03599 .
arXiv:1712.00450 . Rebentrost, Patrick, Maria Schuld, Francesco Petruccione,
van Nieuwenburg, Evert PL, Ye-Hua Liu, and Sebastian D and Seth Lloyd (2016), “Quantum gradient descent and
Huber (2017b), “Learning phase transitions by confusion,” newton’s method for constrained polynomial optimization,”
Nature Physics 13 (5), 435. arXiv preprint arXiv:1612.01789 .
Niu, Murphy Yuezhen, Sergio Boixo, Vadim Smelyanskiy, Reddy, Gautam, Antonio Celani, Terrence J Sejnowski, and
and Hartmut Neven (2018), “Universal quantum con- Massimo Vergassola (2016a), “Learning to soar in turbu-
trol through deep reinforcement learning,” arXiv preprint lent environments,” Proceedings of the National Academy
arXiv:1803.01857 . of Sciences 113 (33), E4877–E4884.
Nomura, Yusuke, Andrew Darmawan, Youhei Yamaji, and Reddy, Gautam, Antonio Celani, and Massimo Vergassola
Masatoshi Imada (2017), “Restricted-boltzmann-machine (2016b), “Infomax strategies for an optimal balance be-
learning for solving strongly correlated quantum systems,” tween exploration and exploitation,” Journal of Statistical
arXiv preprint arXiv:1709.06475 . Physics 163 (6), 1454–1476.
Ohtsuki, Tomi, and Tomoki Ohtsuki (2017), “Deep learning Rem, Benno S, Niklas Käming, Matthias Tarnowski, Luca
the quantum phase transitions in random electron systems: Asteria, Nick Fläschner, Christoph Becker, Klaus Seng-
Applications to three dimensions,” Journal of the Physical stock, and Christof Weitenberg (2018), “Identifying quan-
Society of Japan 86 (4), 044708. tum phase transitions using artificial neural networks on
O’Neil, Cathy (2017), Weapons of math destruction: How big experimental data,” arXiv preprint arXiv:1809.05519 .
data increases inequality and threatens democracy (Broad- Rezende, Danilo Jimenez, Shakir Mohamed, and Daan
way Books). Wierstra (2014), “Stochastic backpropagation and approx-
Papanikolaou, Stefanos, Michail Tzimas, Hengxu Song, An- imate inference in deep generative models,” arXiv preprint
drew CE Reid, and Stephen A Langer (2017), “Learning arXiv:1401.4082.
crystal plasticity using digital image correlation: Exam- Rocchetto, Andrea, Scott Aaronson, Simone Severini, Gon-
ples from discrete dislocation dynamics,” arXiv preprint zalo Carvacho, Davide Poderini, Iris Agresti, Marco Ben-
arXiv:1709.08225 . tivegna, and Fabio Sciarrino (2017), “Experimental learn-
Pedregosa, F, G. Varoquaux, A. Gramfort, V. Michel, ing of quantum states,” arXiv preprint arXiv:1712.00127
B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, .
R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cour- Rocchetto, Andrea, Edward Grant, Sergii Strelchuk,
napeau, M. Brucher, M. Perrot, and E. Duchesnay (2011), Giuseppe Carleo, and Simone Severini (2018), “Learn-
“Scikit-learn: Machine learning in Python,” Journal of Ma- ing hard quantum distributions with variational autoen-
chine Learning Research 12, 2825–2830. coders,” npj Quantum Information 4 (1), 28.
Perdomo-Ortiz, Alejandro, Marcello Benedetti, John Realpe- Rockafellar, Ralph Tyrell (2015), Convex analysis (Princeton
Gómez, and Rupak Biswas (2017), “Opportunities university press).
and challenges for quantum-assisted machine learn- Rodriguez, Alex, and Alessandro Laio (2014), “Clustering by
ing in near-term quantum computers,” arXiv preprint fast search and find of density peaks,” Science 344 (6191),
arXiv:1708.09757 . 1492–1496.
Polyak, Boris T (1964), “Some methods of speeding up the Rokach, Lior, and Oded Maimon (2005), “Clustering meth-
convergence of iteration methods,” USSR Computational ods,” in Data mining and knowledge discovery handbook
Mathematics and Mathematical Physics 4 (5), 1–17. (Springer) pp. 321–352.
Qian, Ning (1999), “On the momentum term in gradient de- Roweis, Sam T, and Lawrence K Saul (2000), “Nonlinear
scent learning algorithms,” Neural networks 12 (1), 145– dimensionality reduction by locally linear embedding,” sci-
151. ence 290 (5500), 2323–2326.
Rabiner, Lawrence R (1989), “A tutorial on hidden markov Rudelius, Tom (2018), “Learning to inflate,” arXiv preprint
models and selected applications in speech recognition,” arXiv:1810.05159 .
Proceedings of the IEEE 77 (2), 257–286. Ruder, Sebastian (2016), “An overview of gradient descent
Radford, Alec, Luke Metz, and Soumith Chintala (2015), optimization algorithms,” arXiv preprint arXiv:1609.04747
“Unsupervised representation learning with deep convo- .
lutional generative adversarial networks,” arXiv preprint Rumelhart, David E, and David Zipser (1985), “Feature dis-
arXiv:1511.06434. covery by competitive learning,” Cognitive science 9 (1),
Ramezanali, Mohammad, Partha P Mitra, and Anirvan M 75–112.
Sengupta (2015), “Critical behavior and universality classes Ruscher, Céline, and Jörg Rottler (2018), “Correlations in
for an algorithmic phase transition in sparse reconstruc- the shear flow of athermal amorphous solids: A principal
tion,” arXiv preprint arXiv:1509.08995. component analysis,” arXiv preprint arXiv:1809.06487 .
120
Saito, Hiroki, and Masaya Kato (2017), “Machine learning preprint arXiv:1703.00810.
technique to find quantum many-body ground states of Sidky, Hythem, and Jonathan K Whitmer (2017), “Learn-
bosons on a lattice,” arXiv preprint arXiv:1709.05468 . ing free energy landscapes using artificial neural networks,”
Salazar, Domingos SP (2017), “Nonequilibrium thermody- arXiv preprint arXiv:1712.02840 .
namics of restricted boltzmann machines,” arXiv preprint Simonite, Tom (2018), “Should data scientist adhere to a hip-
arXiv:1704.08724 . pocratic oath?” Wired .
Sander, Jörg, Martin Ester, Hans-Peter Kriegel, and Xiaowei Singh, Kesar (1981), “On the asymptotic accuracy of efron’s
Xu (1998), “Density-based clustering in spatial databases: bootstrap,” The Annals of Statistics , 1187–1195.
The algorithm gdbscan and its applications,” Data mining Slonim, Noam, Gurinder S Atwal, Gasper Tkacik, and
and knowledge discovery 2 (2), 169–194. William Bialek (2005), “Estimating mutual information
Saxe, Andrew M, James L McClelland, and Surya Gan- and multi–information in large networks,” arXiv preprint
guli (2013), “Exact solutions to the nonlinear dynamics of cs/0502017.
learning in deep linear neural networks,” arXiv preprint Sønderby, Casper Kaae, Tapani Raiko, Lars Maaløe,
arXiv:1312.6120. Søren Kaae Sønderby, and Ole Winther (2016), “Ladder
Schapire, Robert E, and Yoav Freund (2012), Boosting: variational autoencoders,” in Advances in neural informa-
Foundations and algorithms (MIT press). tion processing systems, pp. 3738–3746.
Schindler, Frank, Nicolas Regnault, and Titus Neupert Springenberg, Jost Tobias, Alexey Dosovitskiy, Thomas Brox,
(2017), “Probing many-body localization with neural net- and Martin Riedmiller (2014), “Striving for simplicity: The
works,” Phys. Rev. B 95, 245134. all convolutional net,” arXiv preprint arXiv:1412.6806.
Schmidhuber, Jürgen (2015), “Deep learning in neural net- Sriarunothai, Theeraphot, Sabine Wölk, Gouri Shankar Giri,
works: An overview,” Neural networks 61, 85–117. Nicolai Fries, Vedran Dunjko, Hans J Briegel, and Christof
Schneidman, Elad, Michael J Berry II, Ronen Segev, and Wunderlich (2017), “Speeding-up the decision making of a
William Bialek (2006), “Weak pairwise correlations imply learning agent using an ion trap quantum processor,” arXiv
strongly correlated network states in a neural population,” preprint arXiv:1709.01366 .
Nature 440 (7087), 1007. Srivastava, Nitish, Geoffrey E Hinton, Alex Krizhevsky, Ilya
Schoenholz, Samuel S (2017), “Combining machine learning Sutskever, and Ruslan Salakhutdinov (2014), “Dropout: a
and physics to understand glassy systems,” arXiv preprint simple way to prevent neural networks from overfitting.”
arXiv:1709.08015 . Journal of machine learning research 15 (1), 1929–1958.
Schuld, Maria, Mark Fingerhuth, and Francesco Petruc- Stoudenmire, E Miles (2018), “Learning relevant features
cione (2017), “Implementing a distance-based classifier of data with multi-scale tensor networks,” arXiv preprint
with a quantum interference circuit,” arXiv preprint arXiv:1801.00315 .
arXiv:1703.10793 . Stoudenmire, E Miles, and David J Schwab (2016), “Super-
Schuld, Maria, and Nathan Killoran (2018), “Quantum ma- vised learning with tensor networks,” in Advances in Neural
chine learning in feature hilbert spaces,” arXiv:1803.07128 Information Processing Systems, pp. 4799–4807.
. Stoudenmire, EM, and Steven R White (2012), “Studying
Schuld, Maria, and Francesco Petruccione (2017), “Quan- two-dimensional systems with the density matrix renormal-
tum ensembles of quantum classifiers,” arXiv preprint ization group,” Annu. Rev. Condens. Matter Phys. 3 (1),
arXiv:1704.02146 . 111–128.
Schuld, Maria, Ilya Sinayskiy, and Francesco Petruccione Stratonovich, RL (1957), “On a method of calculating quan-
(2015), “An introduction to quantum machine learning,” tum distribution functions,” in Soviet Physics Doklady,
Contemporary Physics 56 (2), 172–185. Vol. 2, p. 416.
Schwab, David J, Ilya Nemenman, and Pankaj Mehta (2014), Strouse, DJ, and David J Schwab (2017), “The deterministic
Physical review letters 113 (6), 068102. information bottleneck,” Neural computation 29 (6), 1611–
Sethna, James (2006), Statistical mechanics: entropy, or- 1630.
der parameters, and complexity, Vol. 14 (Oxford University Suchsland, Philippe, and Stefan Wessel (2018), “Parameter
Press). diagnostics of phases and phase transition learning by neu-
Shanahan, Phiala E, Daniel Trewartha, and William ral networks,” arXiv preprint arXiv:1802.09876 .
Detmold (2018), “Machine learning action parameters Sutskever, Ilya, James Martens, George Dahl, and Geoffrey
in lattice quantum chromodynamics,” arXiv preprint Hinton (2013), “On the importance of initialization and
arXiv:1801.05784 . momentum in deep learning,” in International conference
Shannon, Claude E (1949), “Communication theory of secrecy on machine learning, pp. 1139–1147.
systems,” Bell Labs Technical Journal 28 (4), 656–715. Sutton, Richard S, and Andrew G Barto (1998), Reinforce-
Shen, Huitao, Junwei Liu, and Liang Fu (2018), “Self-learning ment learning: An introduction, Vol. 1 (MIT press Cam-
monte carlo with deep neural networks,” arXiv preprint bridge).
arXiv:1801.01127 . Swaddle, Michael, Lyle Noakes, Liam Salter, Harry Small-
Shinjo, Kazuya, Shigetoshi Sota, Seiji Yunoki, and bone, and Jingbo Wang (2017), “Generating 3 qubit
Takami Tohyama (2019), “Characterization of photoex- quantum circuits with neural networks,” arXiv preprint
cited states in the half-filled one-dimensional extended hub- arXiv:1703.10743 .
bard model assisted by machine learning,” arXiv preprint Székely, GJ (2003), “E-statistics: The energy of statistical
arXiv:1901.07900 . samples,” Bowling Green State University, Department of
Shlens, Jonathon (2014), “A tutorial on principal component Mathematics and Statistics Technical Report 3 (05), 1–18.
analysis,” arXiv preprint arXiv:1404.1100. Tanaka, Akinori, and Akio Tomiya (2017a), “Detection of
Shwartz-Ziv, Ravid, and Naftali Tishby (2017), “Opening the phase transition via convolutional neural networks,” Jour-
black box of deep neural networks via information,” arXiv nal of the Physical Society of Japan 86 (6), 063001.
121
Tanaka, Akinori, and Akio Tomiya (2017b), “Towards reduc- higgsable gauge groups in 4d f-theory,” arXiv preprint
tion of autocorrelation in hmc by machine learning,” arXiv arXiv:1804.07296 .
preprint arXiv:1712.03893 . Wasserman, Larry (2013), All of statistics: a concise course in
Tenenbaum, Joshua B, Vin De Silva, and John C Langford statistical inference (Springer Science & Business Media).
(2000), “A global geometric framework for nonlinear dimen- Wattenberg, Martin, Fernanda Viégas, and Ian Johnson
sionality reduction,” science 290 (5500), 2319–2323. (2016), “How to use t-sne effectively,” Distill 1 (10), e2.
Tibshirani, Ryan J, et al. (2013), “The lasso problem and Wei, Qianshi, Roger G Melko, and Jeff ZY Chen (2017),
uniqueness,” Electronic Journal of Statistics 7, 1456–1490. “Identifying polymer states by machine learning,” Physical
Tieleman, Tijmen, and Geoffrey Hinton (2009), “Using fast Review E 95 (3), 032504.
weights to improve persistent contrastive divergence,” in Weigt, Martin, Robert A White, Hendrik Szurmant, James A
Proceedings of the 26th Annual International Conference Hoch, and Terence Hwa (2009), “Identification of direct
on Machine Learning (ACM) pp. 1033–1040. residue contacts in protein–protein interaction by message
Tieleman, Tijmen, and Geoffrey Hinton (2012), “Lecture 6.5- passing,” Proceedings of the National Academy of Sciences
rmsprop: Divide the gradient by a running average of its 106 (1), 67–72.
recent magnitude,” COURSERA: Neural networks for ma- Weinstein, Steven (2017), “Learning the einstein-podolsky-
chine learning 4 (2), 26–31. rosen correlations on a restricted boltzmann machine,”
Tishby, Naftali, Fernando C Pereira, and William Bialek arXiv preprint arXiv:1707.03114 .
(2000), “The information bottleneck method,” arXiv Wetzel, Sebastian Johann (2017), “Unsupervised learning of
preprint physics/0004057. phase transitions: from principle component analysis to
Tomczak, P, and H. Puszkarski (2018), “Ferromagnetic res- variational autoencoders,” arXiv preprint arXiv:1703.02435
onance in thin films studied via cross-validation of numer- .
ical solutions of the smit-beljers equation: Application to Wetzel, Sebastian Johann, and Manuel Scherzer (2017),
(ga,mn)as,” Phys. Rev. B 98, 144415. “Machine learning of explicit order parameters: From the
Torgerson, Warren S (1958), “Theory and methods of scaling.” ising model to SU (2) lattice gauge theory,” arXiv preprint
. arXiv:1705.05582 .
Torlai, Giacomo, Guglielmo Mazzola, Juan Carrasquilla, White, Steven R (1992), “Density matrix formulation for
Matthias Troyer, Roger Melko, and Giuseppe Carleo quantum renormalization groups,” Physical review letters
(2017), “Many-body quantum state tomography with neu- 69 (19), 2863.
ral networks,” arXiv preprint arXiv:1703.05334 . White, Tom (2016), “Sampling generative networks:
Torlai, Giacomo, Guglielmo Mazzola, Juan Carrasquilla, Notes on a few effective techniques,” arXiv preprint
Matthias Troyer, Roger Melko, and Giuseppe Carleo arXiv:1609.04468.
(2018), “Neural-network quantum state tomography,” Na- Williams, DRGHR, and Geoffrey Hinton (1986), “Learn-
ture Physics 14 (5), 447. ing representations by back-propagating errors,” Nature
Tramel, Eric W, Marylou Gabrié, Andre Manoel, Francesco 323 (6088), 533–538.
Caltagirone, and Florent Krzakala (2017), “A determin- Wilson, Ashia C, Rebecca Roelofs, Mitchell Stern, Nathan
istic and generalized framework for unsupervised learn- Srebro, and Benjamin Recht (2017), “The marginal value
ing with restricted boltzmann machines,” arXiv preprint of adaptive gradient methods in machine learning,” arXiv
arXiv:1702.03260 . preprint arXiv:1705.08292.
Tubiana, Jérôme, and Rémi Monasson (2017), “Emergence Wilson, Kenneth G, and John Kogut (1974), “The renormal-
of compositional representations in restricted boltzmann ization group and the epsilon expansion,” Physics Reports
machines,” Physical Review Letters 118 (13), 138301. 12 (2), 75–199.
Van Der Maaten, Laurens (2014), “Accelerating t-sne using Witte, RS, and J.S. Witte (2013), Statistics (Wiley).
tree-based algorithms,” The Journal of Machine Learning Wu, Yadong, Pengfei Zhang, Huitao Shen, and Hui Zhai
Research 15 (1), 3221–3245. (2018), “Visualizing neural network developing perturba-
Venderley, Jordan, Vedika Khemani, and Eun-Ah Kim tion theory,” arXiv preprint arXiv:1802.03930 .
(2017), “Machine learning out-of-equilibrium phases of Xie, Junyuan, Ross Girshick, and Ali Farhadi (2016), “Un-
matter,” arXiv preprint arXiv:1711.00020 . supervised deep embedding for clustering analysis,” in In-
Vergassola, Massimo, Emmanuel Villermaux, and Boris I ternational conference on machine learning, pp. 478–487.
Shraiman (2007), Nature 445 (7126), 406. Yang, Tynia, Jinze Liu, Leonard McMillan, and Wei Wang
Vidal, Guifre (2007), “Entanglement renormalization,” Phys- (2006), “A fast approximation to multidimensional scaling,”
ical review letters 99 (22), 220405. in IEEE workshop on Computation Intensive Methods for
Wainwright, Martin J, Michael I Jordan, et al. (2008), Computer Vision.
“Graphical models, exponential families, and variational in- Yang, Xu-Chen, Man-Hong Yung, and Xin Wang
ference,” Foundations and Trends® in Machine Learning (2017), “Neural network designed pulse sequences for
1 (1–2), 1–305. robust control of single-triplet qubits,” arXiv preprint
Wang, Ce, and Hui Zhai (2017), “Unsupervised learning stud- arXiv:1708.00238 .
ies of frustrated classical spin models i: Principle compo- Yedidia, Jonathan (2001), “An idiosyncratic journey beyond
nent analysis,” arXiv preprint arXiv:1706.07977 . mean field theory,” Advanced mean field methods: Theory
Wang, Ce, and Hui Zhai (2018), “Machine learning of frus- and practice , 21–36.
trated classical spin models. ii. kernel principal component Yedidia, Jonathan S, William T Freeman, and Yair Weiss
analysis,” arXiv preprint arXiv:1803.01205 . (2003), “Understanding belief propagation and its general-
Wang, Lei (2017), “Can boltzmann machines discover cluster izations,” Morgan Kaufmann Publishers Inc. San Francisco,
updates?” arXiv preprint arXiv:1702.08586 . CA, USA .
Wang, Yi-Nan, and Zhibai Zhang (2018), “Learning non- Yoshioka, Nobuyuki, Yutaka Akagi, and Hosho Katsura
122
(2017), “Learning disordered topological phases by statisti- learning requires rethinking generalization,” arXiv preprint
cal recovery of symmetry,” arXiv preprint arXiv:1709.05790 arXiv:1611.03530.
. Zhang, Pengfei, Huitao Shen, and Hui Zhai (2017a), “Ma-
You, Yi-Zhuang, Zhao Yang, and Xiao-Liang Qi (2017), chine learning topological invariants with neural networks,”
“Machine learning spatial geometry from entanglement fea- arXiv preprint arXiv:1708.09401 .
tures,” arXiv preprint arXiv:1709.01223 . Zhang, Xiao-Ming, Zi-Wei Cui, Xin Wang, and Man-Hong
Yu, Chao-Hua, Fei Gao, and Qiao-Yan Wen (2017), Yung (2018), “Automatic spin-chain learning to explore the
“Quantum algorithms for ridge regression,” arXiv preprint quantum speed limit,” arXiv preprint arXiv:1802.09248 .
arXiv:1707.09524 . Zhang, Yi, Roger G Melko, and Eun-Ah Kim (2017b), “Ma-
Zdeborová, Lenka, and Florent Krzakala (2016), “Statistical chine learning Z2 quantum spin liquids with quasi-particle
physics of inference: Thresholds and algorithms,” Advances statistics,” arXiv preprint arXiv:1705.01947 .
in Physics 65 (5), 453–552. Zimek, Arthur, Erich Schubert, and Hans-Peter Kriegel
Zeiler, Matthew D (2012), “Adadelta: an adaptive learning (2012), “A survey on unsupervised outlier detection in high-
rate method,” arXiv preprint arXiv:1212.5701. dimensional numerical data,” Statistical Analysis and Data
Zhang, Chengxian, and Xin Wang (2018), “Spin-qubit noise Mining: The ASA Data Science Journal 5 (5), 363–387.
spectroscopy from randomized benchmarking by supervised Zou, Hui, and Trevor Hastie (2005), “Regularization and vari-
learning,” arXiv preprint arXiv:1810.07914 . able selection via the elastic net,” Journal of the Royal Sta-
Zhang, Chiyuan, Samy Bengio, Moritz Hardt, Benjamin tistical Society: Series B (Statistical Methodology) 67 (2),
Recht, and Oriol Vinyals (2016), “Understanding deep 301–320.