Abstract
Deep learning has arguably achieved tremendous success in recent years. Insimple words, deep learning uses the composition of many nonlinear functions tomodel the complex dependency between input features and labels. While neuralnetworks have a long history, recent advances have greatly improved theirperformance in computer vision, natural language processing, etc. From thestatistical and scientific perspective, it is natural to ask: What is deeplearning? What are the new characteristics of deep learning, compared withclassical methods? What are the theoretical foundations of deep learning? Toanswer these questions, we introduce common neural network models (e.g.,convolutional neural nets, recurrent neural nets, generative adversarial nets)and training techniques (e.g., stochastic gradient descent, dropout, batchnormalization) from a statistical point of view. Along the way, we highlightnew characteristics of deep learning (including depth and overparametrization)and explain their practical and theoretical benefits. We also sample recentresults on theories of deep learning, many of which are only suggestive. Whilea complete understanding of deep learning remains elusive, we hope that ourperspectives and discussions serve as a stimulus for new statistical research.
Quick Read (beta)
A Selective Overview of Deep Learning^{0}^{0}footnotetext: Author names are sorted alphabetically.
Abstract
Deep learning has arguably achieved tremendous success in recent years. In simple words, deep learning uses the composition of many nonlinear functions to model the complex dependency between input features and labels. While neural networks have a long history, recent advances have greatly improved their performance in computer vision, natural language processing, etc. From the statistical and scientific perspective, it is natural to ask: What is deep learning? What are the new characteristics of deep learning, compared with classical methods? What are the theoretical foundations of deep learning?
To answer these questions, we introduce common neural network models (e.g., convolutional neural nets, recurrent neural nets, generative adversarial nets) and training techniques (e.g., stochastic gradient descent, dropout, batch normalization) from a statistical point of view. Along the way, we highlight new characteristics of deep learning (including depth and overparametrization) and explain their practical and theoretical benefits. We also sample recent results on theories of deep learning, many of which are only suggestive. While a complete understanding of deep learning remains elusive, we hope that our perspectives and discussions serve as a stimulus for new statistical research.
Keywords: neural networks, overparametrization, stochastic gradient descent, approximation theory, generalization error.
Contents
1 Introduction
Modern machine learning and statistics deal with the problem of learning from data: given a training dataset ${\{({y}_{i},{\bm{x}}_{i})\}}_{1\le i\le n}$ where ${\bm{x}}_{i}\in {\mathbb{R}}^{d}$ is the input and ${y}_{i}\in \mathbb{R}$ is the output^{1}^{1} 1 When the label $y$ is given, this problem is often known as supervised learning. We mainly focus on this paradigm throughout this paper and remark sparingly on its counterpart, unsupervised learning, where $y$ is not given., one seeks a function $f:{\mathbb{R}}^{d}\mapsto \mathbb{R}$ from a certain function class $\mathcal{F}$ that has good prediction performance on test data. This problem is of fundamental significance and finds applications in numerous scenarios. For instance, in image recognition, the input $\bm{x}$ (reps. the output $y$) corresponds to the raw image (reps. its category) and the goal is to find a mapping $f(\cdot )$ that can classify future images accurately. Decades of research efforts in statistical machine learning have been devoted to developing methods to find $f(\cdot )$ efficiently with provable guarantees. Prominent examples include linear classifiers (e.g., linear$$/$$logistic regression, linear discriminant analysis), kernel methods (e.g., support vector machines), treebased methods (e.g., decision trees, random forests), nonparametric regression (e.g., nearest neighbors, local kernel smoothing), etc. Roughly speaking, each aforementioned method corresponds to a different function class $\mathcal{F}$ from which the final classifier $f(\cdot )$ is chosen.
Deep learning [70], in its simplest form, proposes the following compositional function class:
$$\left\{f(\bm{x};\bm{\theta})={\mathbf{W}}_{L}{\bm{\sigma}}_{L}({\mathbf{W}}_{L1}\mathrm{\cdots}{\bm{\sigma}}_{2}({\mathbf{W}}_{2}{\bm{\sigma}}_{1}({\mathbf{W}}_{1}\bm{x})))\right\bm{\theta}=\{{\mathbf{W}}_{1},\mathrm{\dots},{\mathbf{W}}_{L}\}\}.$$  (1) 
Here, for each $1\le l\le L$, ${\bm{\sigma}}_{\mathrm{\ell}}(\cdot )$ is some nonlinear function, and $\bm{\theta}=\{{\mathbf{W}}_{1},\mathrm{\dots},{\mathbf{W}}_{L}\}$ consists of matrices with appropriate sizes. Though simple, deep learning has made significant progress towards addressing the problem of learning from data over the past decade. Specifically, it has performed close to or better than humans in various important tasks in artificial intelligence, including image recognition [50], game playing [113], and machine translation [131]. Owing to its great promise, the impact of deep learning is also growing rapidly in areas beyond artificial intelligence; examples include statistics [15, 110, 76, 103, 41], applied mathematics [129, 22], clinical research [28], etc.
Model  Year  # Layers  # Params  Top5 error 

Shallow  $$  —  —  $>25\%$ 
AlexNet  $2012$  $8$  $61$M  $16.4\%$ 
VGG19  $2014$  $19$  $144$M  $7.3\%$ 
GoogleNet  $2014$  $22$  $7$M  $6.7\%$ 
ResNet$152$  $2015$  $152$  $60$M  $3.6\%$ 
To get a better idea of the success of deep learning, let us take the ImageNet Challenge [106] (also known as ILSVRC) as an example. In the classification task, one is given a training dataset consisting of 1.2 million color images with $1000$ categories, and the goal is to classify images based on the input pixels. The performance of a classifier is then evaluated on a test dataset of 100 thousand images, and in the end the top5 error^{2}^{2} 2 The algorithm makes an error if the true label is not contained in the $5$ predictions made by the algorithm. is reported. Table 1 highlights a few popular models and their corresponding performance. As can be seen, deep learning models (the second to the last rows) have a clear edge over shallow models (the first row) that fit linear models$$/$$treebased models on handcrafted features. This significant improvement raises a foundational question:

Why is deep learning better than classical methods on tasks like image recognition?
1.1 Intriguing new characteristics of deep learning
It is widely acknowledged that two indispensable factors contribute to the success of deep learning, namely (1) huge datasets that often contain millions of samples and (2) immense computing power resulting from clusters of graphics processing units (GPUs). Admittedly, these resources are only recently available: the latter allows to train larger neural networks which reduces biases and the former enables variance reduction. However, these two alone are not sufficient to explain the mystery of deep learning due to some of its “dreadful” characteristics: (1) overparametrization: the number of parameters in stateoftheart deep learning models is often much larger than the sample size (see Table 1), which gives them the potential to overfit the training data, and (2) nonconvexity: even with the help of GPUs, training deep learning models is still NPhard [8] in the worst case due to the highly nonconvex loss function to minimize. In reality, these characteristics are far from nightmares. This sharp difference motivates us to take a closer look at the salient features of deep learning, which we single out a few below.
1.1.1 Depth
Deep learning expresses complicated nonlinearity through composing many nonlinear functions; see (1). The rationale for this multilayer structure is that, in many realworld datasets such as images, there are different levels of features and lowerlevel features are building blocks of higherlevel ones. See [133] for a visualization of trained features of convolutional neural nets; here in Figure 1, we sample and visualize weights from a pretrained AlexNet model. This intuition is also supported by empirical results from physiology and neuroscience [56, 2]. The use of function composition marks a sharp difference from traditional statistical methods such as projection pursuit models [38] and multiindex models [73, 27]. It is often observed that depth helps efficiently extract features that are representative of a dataset. In comparison, increasing width (e.g., number of basis functions) in a shallow model leads to less improvement. This suggests that deep learning models excel at representing a very different function space that is suitable for complex datasets.
1.1.2 Algorithmic regularization
The statistical performance of neural networks (e.g., test accuracy) depends heavily on the particular optimization algorithms used for training [130]. This is very different from many classical statistical problems, where the related optimization problems are less complicated. For instance, when the associated optimization problem has a relatively simple structure (e.g., convex objective functions, linear constraints), the solution to the optimization problem can often be unambiguously computed and analyzed. However, in deep neural networks, due to overparametrization, there are usually many local minima with different statistical performance [72]. Nevertheless, common practice runs stochastic gradient descent with random initialization and finds model parameters with very good prediction accuracy.
1.1.3 Implicit prior learning
It is well observed that deep neural networks trained with only the raw inputs (e.g., pixels of images) can provide a useful representation of the data. This means that after training, the units of deep neural networks can represent features such as edges, corners, wheels, eyes, etc.; see [133]. Importantly, the training process is automatic in the sense that no human knowledge is involved (other than hyperparameter tuning). This is very different from traditional methods, where algorithms are designed after structural assumptions are posited. It is likely that training an overparametrized model efficiently learns and incorporates the prior distribution $p(\bm{x})$ of the input, even though deep learning models are themselves discriminative models. With automatic representation of the prior distribution, deep learning typically performs well on similar datasets (but not very different ones) via transfer learning.
(a) MNIST images  (b) training and test accuracies 
1.2 Towards theory of deep learning
Despite the empirical success, theoretical support for deep learning is still in its infancy. Setting the stage, for any classifier $f$, denote by $\mathbb{E}(f)$ the expected risk on fresh sample (a.k.a. test error, prediction error or generalization error), and by ${\mathbb{E}}_{n}(f)$ the empirical risk$$/$$training error averaged over a training dataset. Arguably, the key theoretical question in deep learning is
why is $\mathrm{E}\mathit{}\mathrm{(}{\widehat{f}}_{n}\mathrm{)}$ small, where ${\widehat{f}}_{n}$ is the classifier returned by the training algorithm?
We follow the conventional approximationestimation decomposition (sometimes, also biasvariance tradeoff) to decompose the term $\mathbb{E}({\widehat{f}}_{n})$ into two parts. Let $\mathcal{F}$ be the function space expressible by a family of neural nets. Define ${f}^{*}={\text{argmin}}_{f}\mathbb{E}(f)$ to be the best possible classifier and ${f}_{\mathcal{F}}^{*}={\text{argmin}}_{f\in \mathcal{F}}\mathbb{E}(f)$ to be the best classifier in $\mathcal{F}$. Then, we can decompose the excess error $\mathcal{E}\triangleq \mathbb{E}({\widehat{f}}_{n})\mathbb{E}({f}^{*})$ into two parts:
$$\mathcal{E}=\underset{\text{approximation error}}{\underset{\u23df}{\mathbb{E}({f}_{\mathcal{F}}^{*})\mathbb{E}({f}^{*})}}+\underset{\text{estimation error}}{\underset{\u23df}{\mathbb{E}({\widehat{f}}_{n})\mathbb{E}({f}_{\mathcal{F}}^{*})}}.$$  (2) 
Both errors can be small for deep learning (cf. Figure 2), which we explain below.

•
The approximation error is determined by the function class $\mathcal{F}$. Intuitively, the larger the class, the smaller the approximation error. Deep learning models use many layers of nonlinear functions (Figure 3)that can drive this error small. Indeed, in Section 5, we provide recent theoretical progress of its representation power. For example, deep models allow efficient representation of interactions among variable while shallow models cannot.

•
The estimation error reflects the generalization power, which is influenced by both the complexity of the function class $\mathcal{F}$ and the properties of the training algorithms. Interestingly, for overparametrized deep neural nets, stochastic gradient descent typically results in a nearzero training error (i.e., ${\mathbb{E}}_{n}({\widehat{f}}_{n})\approx 0$; see e.g. left panel of Figure 2). Moreover, its generalization error $\mathbb{E}({\widehat{f}}_{n})$ remains small or moderate. This “counterintuitive” behavior suggests that for overparametrized models, gradientbased algorithms enjoy benign statistical properties; we shall see in Section 7 that gradient descent enjoys implicit regularization in the overparametrized regime even without explicit regularization (e.g., ${\mathrm{\ell}}_{2}$ regularization).
The above two points lead to the following heuristic explanation of the success of deep learning models. The large depth of deep neural nets and heavy overparametrization lead to small or zero training errors, even when running simple algorithms with moderate number of iterations. In addition, these simple algorithms with moderate number of steps do not explore the entire function space and thus have limited complexities, which results in small generalization error with a large sample size. Thus, by combining the two aspects, it explains heuristically that the test error is also small.
1.3 Roadmap of the paper
We first introduce basic deep learning models in Sections 2–4, and then examine their representation power via the lens of approximation theory in Section 5. Section 6 is devoted to training algorithms and their ability of driving the training error small. Then we sample recent theoretical progress towards demystifying the generalization power of deep learning in Section 7. Along the way, we provide our own perspectives, and at the end we identify a few interesting questions for future research in Section 8. The goal of this paper is to present suggestive methods and results, rather than giving conclusive arguments (which is currently unlikely) or a comprehensive survey. We hope that our discussion serves as a stimulus for new statistics research.
2 Feedforward neural networks
Before introducing the vanilla feedforward neural nets, let us set up necessary notations for the rest of this section. We focus primarily on classification problems, as regression problems can be addressed similarly. Given the training dataset ${\{({y}_{i},{\bm{x}}_{i})\}}_{1\le i\le n}$ where ${y}_{i}\in [K]\triangleq \{1,2,\mathrm{\dots},K\}$ and ${\bm{x}}_{i}\in {\mathbb{R}}^{d}$ are independent across $i\in [n]$, supervised learning aims at finding a (possibly random) function $\widehat{f}(\bm{x})$ that predicts the outcome $y$ for a new input $\bm{x}$, assuming $(y,\bm{x})$ follows the same distribution as $({y}_{i},{\bm{x}}_{i})$. In the terminology of machine learning, the input ${\bm{x}}_{i}$ is often called the feature, the output ${y}_{i}$ called the label, and the pair $({y}_{i},{\bm{x}}_{i})$ is an example. The function $\widehat{f}$ is called the classifier, and estimation of $\widehat{f}$ is training or learning. The performance of $\widehat{f}$ is evaluated through the prediction error $\mathbb{P}(y\ne \widehat{f}(\bm{x}))$, which can be often estimated from a separate test dataset.
As with classical statistical estimation, for each $k\in [K]$, a classifier approximates the conditional probability $\mathbb{P}(y=k\bm{x})$ using a function ${f}_{k}(\bm{x};{\bm{\theta}}_{k})$ parametrized by ${\bm{\theta}}_{k}$. Then the category with the highest probability is predicted. Thus, learning is essentially estimating the parameters ${\bm{\theta}}_{k}$. In statistics, one of the most popular methods is (multinomial) logistic regression, which stipulates a specific form for the functions ${f}_{k}(\bm{x};{\bm{\theta}}_{k})$: let ${z}_{k}={\bm{x}}^{\top}{\bm{\beta}}_{k}+{\alpha}_{k}$ and ${f}_{k}(\bm{x};{\bm{\theta}}_{k})={Z}^{1}\mathrm{exp}({z}_{k})$ where $Z={\sum}_{k=1}^{K}\mathrm{exp}({z}_{k})$ is a normalization factor to make ${\{{f}_{k}(\bm{x};{\bm{\theta}}_{k})\}}_{1\le k\le K}$ a valid probability distribution. It is clear that logistic regression induces linear decision boundaries in ${\mathbb{R}}^{d}$, and hence it is restrictive in modeling nonlinear dependency between $y$ and $\bm{x}$. The deep neural networks we introduce below provide a flexible framework for modeling nonlinearity in a fairly general way.
2.1 Model setup
From the high level, deep neural networks (DNNs) use composition of a series of simple nonlinear functions to model nonlinearity
$${\bm{h}}^{(L)}={\mathbf{g}}^{(L)}\circ {\mathbf{g}}^{(L1)}\circ \mathrm{\dots}\circ {\mathbf{g}}^{(1)}(\bm{x}),$$ 
where $\circ $ denotes composition of two functions and $L$ is the number of hidden layers, and is usually called depth of a NN model. Letting ${\bm{h}}^{(0)}\triangleq \bm{x}$, one can recursively define ${\bm{h}}^{(l)}={\mathbf{g}}^{(l)}\left({\bm{h}}^{(l1)}\right)$ for all $\mathrm{\ell}=1,2,\mathrm{\dots},L$. The feedforward neural networks, also called the multilayer perceptrons (MLPs), are neural nets with a specific choice of ${\mathbf{g}}^{(l)}$: for $\mathrm{\ell}=1,\mathrm{\dots},L$, define
$${\bm{h}}^{(\mathrm{\ell})}={\mathbf{g}}^{(l)}\left({\bm{h}}^{(l1)}\right)\triangleq \bm{\sigma}\left({\mathbf{W}}^{(\mathrm{\ell})}{\bm{h}}^{(\mathrm{\ell}1)}+{\bm{b}}^{(\mathrm{\ell})}\right),$$  (3) 
where ${\mathbf{W}}^{(l)}$ and ${\bm{b}}^{(l)}$ are the weight matrix and the bias$$/$$intercept, respectively, associated with the $l$th layer, and $\bm{\sigma}(\cdot )$ is usually a simple given (known) nonlinear function called the activation function. In words, in each layer $\mathrm{\ell}$, the input vector ${\bm{h}}^{(\mathrm{\ell}1)}$ goes through an affine transformation first and then passes through a fixed nonlinear function $\bm{\sigma}(\cdot )$. See Figure 3 for an illustration of a simple MLP with two hidden layers. The activation function $\bm{\sigma}(\cdot )$ is usually applied elementwise, and a popular choice is the ReLU (Rectified Linear Unit) function:
$${[\bm{\sigma}(\bm{z})]}_{j}=\mathrm{max}\{{z}_{j},0\}.$$  (4) 
Other choices of activation functions include leaky ReLU, $\mathrm{tanh}$ function [79] and the classical sigmoid function ${(1+{e}^{z})}^{1}$, which is less used now.
Given an output ${\bm{h}}^{(L)}$ from the final hidden layer and a label $y$, we can define a loss function to minimize. A common loss function for classification problems is the multinomial logistic loss. Using the terminology of deep learning, we say that ${\bm{h}}^{(L)}$ goes through an affine transformation and then the softmax function:
$${f}_{k}(\bm{x};\bm{\theta})\triangleq \frac{\mathrm{exp}({z}_{k})}{{\sum}_{k}\mathrm{exp}({z}_{k})},\forall k\in [K],\text{where}\bm{z}={\mathbf{W}}^{(L+1)}{\bm{h}}^{(L)}+{\bm{b}}^{(L+1)}\in {\mathbb{R}}^{K}.$$ 
Then the loss is defined to be the crossentropy between the label $y$ (in the form of an indicator vector) and the score vector ${({f}_{1}(\bm{x};\bm{\theta}),\mathrm{\dots},{f}_{K}(\bm{x};\bm{\theta}))}^{\top}$, which is exactly the negative loglikelihood of the multinomial logistic regression model:
$$\mathcal{L}(\bm{f}(\bm{x};\bm{\theta}),y)=\sum _{k=1}^{K}\mathrm{\U0001d7d9}\{y=k\}\mathrm{log}{p}_{k},$$  (5) 
where $\bm{\theta}\triangleq \{{\mathbf{W}}^{(\mathrm{\ell})},{\bm{b}}^{(\mathrm{\ell})}:1\le \mathrm{\ell}\le L+1\}$. As a final remark, the number of parameters scales with both the depth $L$ and the width (i.e., the dimensionality of ${\mathbf{W}}^{(\mathrm{\ell})}$), and hence it can be quite large for deep neural nets.
2.2 Backpropagation in computational graphs
Training neural networks follows the empirical risk minimization paradigm that minimizes the loss (e.g., (5)) over all the training data. This minimization is usually done via stochastic gradient descent (SGD). In a way similar to gradient descent, SGD starts from a certain initial value ${\bm{\theta}}^{0}$ and then iteratively updates the parameters ${\bm{\theta}}^{t}$ by moving it in the direction of the negative gradient. The difference is that, in each update, a small subsample $\mathcal{B}\subset [n]$ called a minibatch—which is typically of size 32–512—is randomly drawn and the gradient calculation is only on $\mathcal{B}$ instead of the full batch $[n]$. This saves considerably the computational cost in calculation of gradient. By the law of large numbers, this stochastic gradient should be close to the full sample one, albeit with some random fluctuations. A pass of the whole training set is called an epoch. Usually, after several or tens of epochs, the error on a validation set levels off and training is complete. See Section 6 for more details and variants on training algorithms.
The key to the above training procedure, namely SGD, is the calculation of the gradient $\nabla {\mathrm{\ell}}_{\mathcal{B}}(\bm{\theta})$, where
$${\mathrm{\ell}}_{\mathcal{B}}(\bm{\theta})\triangleq {\mathcal{B}}^{1}\sum _{i\in \mathcal{B}}\mathcal{L}(\bm{f}({\bm{x}}_{i};\bm{\theta}),{y}_{i}).$$  (6) 
Gradient computation, however, is in general nontrivial for complex models, and it is susceptible to numerical instability for a model with large depth. Here, we introduce an efficient approach, namely backpropagation, for computing gradients in neural networks.
Backpropagation [105] is a direct application of the chain rule in networks. As the name suggests, the calculation is performed in a backward fashion: one first computes $\partial {\mathrm{\ell}}_{\mathcal{B}}/\partial {\bm{h}}^{(L)}$, then $\partial {\mathrm{\ell}}_{\mathcal{B}}/\partial {\bm{h}}^{(L1)}$, $\mathrm{\dots}$, and finally $\partial {\mathrm{\ell}}_{\mathcal{B}}/\partial {\bm{h}}^{(1)}$. For example, in the case of the ReLU activation function^{3}^{3} 3 The issue of nondifferentiability at the origin is often ignored in implementation., we have the following recursive$$/$$backward relation
$$\frac{\partial {\mathrm{\ell}}_{\mathcal{B}}}{\partial {\bm{h}}^{(\mathrm{\ell}1)}}=\frac{\partial {\bm{h}}^{(\mathrm{\ell})}}{\partial {\bm{h}}^{(\mathrm{\ell}1)}}\cdot \frac{\partial {\mathrm{\ell}}_{\mathcal{B}}}{\partial {\bm{h}}^{(\mathrm{\ell})}}={({\mathbf{W}}^{(\mathrm{\ell})})}^{\top}\mathrm{\U0001d5bd\U0001d5c2\U0001d5ba\U0001d5c0}\left(\mathrm{\U0001d7d9}\{{\mathbf{W}}^{(\mathrm{\ell})}{\bm{h}}^{(\mathrm{\ell}1)}+{\bm{b}}^{(\mathrm{\ell})}\ge \mathrm{\U0001d7ce}\}\right)\frac{\partial {\mathrm{\ell}}_{\mathcal{B}}}{\partial {\bm{h}}^{(\mathrm{\ell})}}$$  (7) 
where $\mathrm{\U0001d5bd\U0001d5c2\U0001d5ba\U0001d5c0}(\cdot )$ denotes a diagonal matrix with elements given by the argument. Note that the calculation of $\partial {\mathrm{\ell}}_{\mathcal{B}}/\partial {\bm{h}}^{(\mathrm{\ell}1)}$ depends on $\partial {\mathrm{\ell}}_{\mathcal{B}}/\partial {\bm{h}}^{(\mathrm{\ell})}$, which is the partial derivatives from the next layer. In this way, the derivatives are “backpropagated” from the last layer to the first layer. These derivatives $\{\partial {\mathrm{\ell}}_{\mathcal{B}}/\partial {\bm{h}}^{(\mathrm{\ell})}\}$ are then used to update the parameters. For instance, the gradient update for ${\mathbf{W}}^{(\mathrm{\ell})}$ is given by
$${\mathbf{W}}^{(\mathrm{\ell})}\leftarrow {\mathbf{W}}^{(\mathrm{\ell})}\eta \frac{\partial {\mathrm{\ell}}_{\mathcal{B}}}{\partial {\mathbf{W}}^{(\mathrm{\ell})}},\text{where}\mathit{\hspace{1em}}\frac{\partial {\mathrm{\ell}}_{\mathcal{B}}}{\partial {W}_{jm}^{(\mathrm{\ell})}}=\frac{\partial {\mathrm{\ell}}_{\mathcal{B}}}{\partial {h}_{j}^{(\mathrm{\ell})}}\cdot {\sigma}^{\prime}\cdot {h}_{m}^{(\mathrm{\ell}1)},$$  (8) 
where ${\sigma}^{\prime}=1$ if the $j$th element of ${\mathbf{W}}^{(\mathrm{\ell})}{\bm{h}}^{(\mathrm{\ell}1)}+{\bm{b}}^{(\mathrm{\ell})}$ is nonnegative, and ${\sigma}^{\prime}=0$ otherwise. The step size $\eta >0$, also called the learning rate, controls how much parameters are changed in a single update.
A more general way to think about neural network models and training is to consider computational graphs. Computational graphs are directed acyclic graphs that represent functional relations between variables. They are very convenient and flexible to represent function composition, and moreover, they also allow an efficient way of computing gradients. Consider an MLP with a single hidden layer and an ${\mathrm{\ell}}_{2}$ regularization:
$${\mathrm{\ell}}_{\mathcal{B}}^{\lambda}(\bm{\theta})={\mathrm{\ell}}_{\mathcal{B}}(\bm{\theta})+{r}_{\lambda}(\bm{\theta})={\mathrm{\ell}}_{\mathcal{B}}(\bm{\theta})+\lambda \left(\sum _{j,{j}^{\prime}}{\left({W}_{j,{j}^{\prime}}^{(1)}\right)}^{2}+\sum _{j,{j}^{\prime}}{\left({W}_{j,{j}^{\prime}}^{(2)}\right)}^{2}\right),$$  (9) 
where ${\mathrm{\ell}}_{\mathcal{B}}(\bm{\theta})$ is the same as (6), and $\lambda \ge 0$ is a tuning parameter. A similar example is considered in [45]. The corresponding computational graph is shown in Figure 4. Each node represents a function (inside a circle), which is associated with an output of that function (outside a circle). For example, we view the term ${\mathrm{\ell}}_{\mathcal{B}}(\bm{\theta})$ as a result of $4$ compositions: first the input data $\bm{x}$ multiplies the weight matrix ${\mathbf{W}}^{(1)}$ resulting in ${\bm{u}}^{(1)}$, then it goes through the ReLU activation function relu resulting in ${\bm{h}}^{(1)}$, then it multiplies another weight matrix ${\mathbf{W}}^{(2)}$ leading to $\bm{p}$, and finally it produces the crossentropy with label $y$ as in (5). The regularization term is incorporated in the graph similarly.
A forward pass is complete when all nodes are evaluated starting from the input $\bm{x}$. A backward pass then calculates the gradients of ${\mathrm{\ell}}_{\mathcal{B}}^{\lambda}$ with respect to all other nodes in the reverse direction. Due to the chain rule, the gradient calculation for a variable (say, $\partial {\mathrm{\ell}}_{\mathcal{B}}/\partial {\bm{u}}^{(1)}$) is simple: it only depends on the gradient value of the variables ($\partial {\mathrm{\ell}}_{\mathcal{B}}/\partial \bm{h}$) the current node points to, and the function derivative evaluated at the current variable value (${\bm{\sigma}}^{\prime}({\bm{u}}^{(1)})$). Thus, in each iteration, a computation graph only needs to (1) calculate and store the function evaluations at each node in the forward pass, and then (2) calculate all derivatives in the backward pass.
3 Popular models
Moving beyond vanilla feedforward neural networks, we introduce two other popular deep learning models, namely, the convolutional neural networks (CNNs) and the recurrent neural networks (RNNs). One important characteristic shared by the two models is weight sharing, that is some model parameters are identical across locations in CNNs or across time in RNNs. This is related to the notion of translational invariance in CNNs and stationarity in RNNs. At the end of this section, we introduce a modular thinking for constructing more flexible neural nets.
3.1 Convolutional neural networks
The convolutional neural network (CNN) [71, 40] is a special type of feedforward neural networks that is tailored for image processing. More generally, it is suitable for analyzing data with salient spatial structures. In this subsection, we focus on image classification using CNNs, where the raw input (image pixels) and features of each hidden layer are represented by a 3D tensor $\bm{X}\in {\mathbb{R}}^{{d}_{1}\times {d}_{2}\times {d}_{3}}$. Here, the first two dimensions ${d}_{1},{d}_{2}$ of $\bm{X}$ indicate spatial coordinates of an image while the third ${d}_{3}$ indicates the number of channels. For instance, ${d}_{3}$ is $3$ for the raw inputs due to the red, green and blue channels, and ${d}_{3}$ can be much larger (say, 256) for hidden layers. Each channel is also called a feature map, because each feature map is specialized to detect the same feature at different locations of the input, which we will soon explain. We now introduce two building blocks of CNNs, namely the convolutional layer and the pooling layer.

1.
Convolutional layer (CONV). A convolutional layer has the same functionality as described in (3), where the input feature $\bm{X}\in {\mathbb{R}}^{{d}_{1}\times {d}_{2}\times {d}_{3}}$ goes through an affine transformation first and then an elementwise nonlinear activation. The difference lies in the specific form of the affine transformation. A convolutional layer uses a number of filters to extract local features from the previous input. More precisely, each filter is represented by a 3D tensor ${\bm{F}}_{k}\in {\mathbb{R}}^{w\times w\times {d}_{3}}$ ($1\le k\le {\stackrel{~}{d}}_{3}$), where $w$ is the size of the filter (typically 3 or 5) and ${\stackrel{~}{d}}_{3}$ denotes the total number of filters. Note that the third dimension ${d}_{3}$ of ${\bm{F}}_{k}$ is equal to that of the input feature $\bm{X}$. For this reason, one usually says that the filter has size $w\times w$, while suppressing the third dimension ${d}_{3}$. Each filter ${\bm{F}}_{k}$ then convolves with the input feature $\bm{X}$ to obtain one single feature map ${\bm{O}}^{k}\in {\mathbb{R}}^{({d}_{1}w+1)\times ({d}_{1}w+1)}$, where^{4}^{4} 4 To simplify notation, we omit the bias/intercept term associated with each filter.
$${O}_{ij}^{k}=\u27e8{\left[\bm{X}\right]}_{ij},{\bm{F}}_{k}\u27e9=\sum _{{i}^{\prime}=1}^{w}\sum _{{j}^{\prime}=1}^{w}\sum _{l=1}^{{d}_{3}}{[\bm{X}]}_{i+{i}^{\prime}1,j+{j}^{\prime}1,l}{[{\bm{F}}_{k}]}_{{i}^{\prime},{j}^{\prime},l}.$$ (10) Here ${[\bm{X}]}_{ij}\in {\mathbb{R}}^{w\times w\times {d}_{3}}$ is a small “patch” of $\bm{X}$ starting at location $(i,j)$. See Figure 5 for an illustration of the convolution operation. If we view the 3D tensors ${[\bm{X}]}_{ij}$ and ${\bm{F}}_{k}$ as vectors, then each filter essentially computes their inner product with a part of $\bm{X}$ indexed by $i,j$ (which can be also viewed as convolution, as its name suggests). One then pack the resulted feature maps $\{{\bm{O}}^{k}\}$ into a 3D tensor $\bm{O}$ with size $({d}_{1}w+1)\times ({d}_{1}w+1)\times {\stackrel{~}{d}}_{3}$, where
$${[\bm{O}]}_{ijk}={[{\bm{O}}^{k}]}_{ij}.$$ (11) The outputs of convolutional layers are then followed by nonlinear activation functions. In the ReLU case, we have
$${\stackrel{~}{X}}_{ijk}=\sigma ({O}_{ijk}),\forall i\in [{d}_{1}w+1],j\in [{d}_{2}w+1],k\in [{\stackrel{~}{d}}_{3}].$$ (12) The convolution operation (10) and the ReLU activation (12) work together to extract features $\stackrel{~}{\bm{X}}$ from the input $\bm{X}$. Different from feedforward neural nets, the filters ${\bm{F}}_{k}$ are shared across all locations $(i,j)$. A patch ${[\bm{X}]}_{ij}$ of an input responds strongly (that is, producing a large value) to a filter ${\bm{F}}_{k}$ if they are positively correlated. Therefore intuitively, each filter ${\bm{F}}_{k}$ serves to extract features similar to ${\bm{F}}_{k}$.
As a side note, after the convolution (10), the spatial size ${d}_{1}\times {d}_{2}$ of the input $\bm{X}$ shrinks to $({d}_{1}w+1)\times ({d}_{2}w+1)$ of $\stackrel{~}{\bm{X}}$. However one may want the spatial size unchanged. This can be achieved via padding, where one appends zeros to the margins of the input $\bm{X}$ to enlarge the spatial size to $({d}_{1}+w1)\times ({d}_{2}+w1)$. In addition, a stride in the convolutional layer determines the gap ${i}^{\prime}i$ and ${j}^{\prime}j$ between two patches ${\bm{X}}_{ij}$ and ${\bm{X}}_{{i}^{\prime}{j}^{\prime}}$: in (10) the stride is $1$, and a larger stride would lead to feature maps with smaller sizes.

2.
Pooling layer (POOL). A pooling layer aggregates the information of nearby features into a single one. This downsampling operation reduces the size of the features for subsequent layers and saves computation. One common form of the pooling layer is composed of the $2\times 2$ maxpooling filter. It computes $\mathrm{max}\{{X}_{i,j,k},{X}_{i+1,j,k},{X}_{i,j+1,k},{X}_{i+1,j+1,k}\}$, that is, the maximum of the $2\times 2$ neighborhood in the spatial coordinates; see Figure 6 for an illustration. Note that the pooling operation is done separately for each feature map $k$. As a consequence, a $2\times 2$ maxpooling filter acting on $\bm{X}\in {\mathbb{R}}^{{d}_{1}\times {d}_{2}\times {d}_{3}}$ will result in an output of size ${d}_{1}/2\times {d}_{2}/2\times {d}_{3}$. In addition, the pooling layer does not involve any parameters to optimize. Pooling layers serve to reduce redundancy since a small neighborhood around a location $(i,j)$ in a feature map is likely to contain the same information.
In addition, we also use fullyconnected layers as building blocks, which we have already seen in Section 2. Each fullyconnected layer treats input tensor $\bm{X}$ as a vector $\mathrm{Vec}(\bm{X})$, and computes $\stackrel{~}{\bm{X}}=\bm{\sigma}(\mathbf{W}\mathrm{Vec}(\bm{X}))$. A fullyconnected layer does not use weight sharing and is often used in the last few layers of a CNN. As an example, Figure 7 depicts the wellknown LeNet 5 [71], which is composed of two sets of CONVPOOL layers and three fullyconnected layers.
3.2 Recurrent neural networks
Recurrent neural nets (RNNs) are another family of powerful models, which are designed to process time series data and other sequence data. RNNs have successful applications in speech recognition [107], machine translation [131], genome sequencing [21], etc. The structure of an RNN naturally forms a computational graph, and can be easily combined with other structures such as CNNs to build large computational graph models for complex tasks. Here we introduce vanilla RNNs and improved variants such as long shortterm memory (LSTM).
(a) Onetomany  (b) Manytoone  (c) Manytomany 
3.2.1 Vanilla RNNs
Suppose we have general time series inputs ${\bm{x}}_{1},{\bm{x}}_{2},\mathrm{\dots},{\bm{x}}_{T}$. A vanilla RNN models the “hidden state” at time $t$ by a vector ${\bm{h}}_{t}$, which is subject to the recursive formula
$${\bm{h}}_{t}={\bm{f}}_{\bm{\theta}}({\bm{h}}_{t1},{\bm{x}}_{t}).$$  (13) 
Here, ${f}_{\bm{\theta}}$ is generally a nonlinear function parametrized by $\bm{\theta}$. Concretely, a vanilla RNN with one hidden layer has the following form^{5}^{5} 5 Similar to the activation function $\bm{\sigma}(\cdot )$, the function $\text{\mathbf{t}\mathbf{a}\mathbf{n}\mathbf{h}}(\cdot )$ means elementwise operations.
${\bm{h}}_{t}$  $=\text{\mathbf{t}\mathbf{a}\mathbf{n}\mathbf{h}}\left({\mathbf{W}}_{hh}{\bm{h}}_{t1}+{\mathbf{W}}_{xh}{\bm{x}}_{t}+{\bm{b}}_{\bm{h}}\right),\text{where}\mathrm{tanh}(a)=\frac{{e}^{2a}1}{{e}^{2a}+1},$  
${\bm{z}}_{t}$  $=\bm{\sigma}\left({\mathbf{W}}_{hy}{\bm{h}}_{t}+{\bm{b}}_{\bm{z}}\right),$ 
where ${\mathbf{W}}_{hh},{\mathbf{W}}_{xh},{\mathbf{W}}_{hy}$ are trainable weight matrices, ${\bm{b}}_{\bm{h}},{\bm{b}}_{\bm{z}}$ are trainable bias vectors, and ${\bm{z}}_{t}$ is the output at time $t$. Like many classical time series models, those parameters are shared across time. Note that in different applications, we may have different input/output settings (cf. Figure 8). Examples include

•
Onetomany: a single input with multiple outputs; see Figure 8(a). A typical application is image captioning, where the input is an image and outputs are a series of words.

•
Manytoone: multiple inputs with a single output; see Figure 8(b). One application is text sentiment classification, where the input is a series of words in a sentence and the output is a label (e.g., positive vs. negative).

•
Manytomany: multiple inputs and outputs; see Figure 8(c). This is adopted in machine translation, where inputs are words of a source language (say Chinese) and outputs are words of a target language (say English).
As the case with feedforward neural nets, we minimize a loss function using backpropagation, where the loss is typically
$${\mathrm{\ell}}_{\mathcal{T}}(\bm{\theta})=\sum _{t\in \mathcal{T}}\mathcal{L}({y}_{t},{\bm{z}}_{t})=\sum _{t\in \mathcal{T}}\sum _{k=1}^{K}\mathrm{\U0001d7d9}\{{y}_{t}=k\}\mathrm{log}\left(\frac{\mathrm{exp}({[{\bm{z}}_{t}]}_{k})}{{\sum}_{k}\mathrm{exp}({[{\bm{z}}_{t}]}_{k})}\right),$$ 
where $K$ is the number of categories for classification (e.g., size of the vocabulary in machine translation), and $\mathcal{T}\subset [T]$ is the length of the output sequence. During the training, the gradients $\partial {\mathrm{\ell}}_{\mathcal{T}}/\partial {\bm{h}}_{t}$ are computed in the reverse time order (from $T$ to $t$). For this reason, the training process is often called backpropagation through time.
One notable drawback of vanilla RNNs is that, they have difficulty in capturing longrange dependencies in sequence data when the length of the sequence is large. This is sometimes due to the phenomenon of exploding$$/$$vanishing gradients. Take Figure 8(c) as an example. Computing $\partial {\mathrm{\ell}}_{\mathcal{T}}/\partial {\bm{h}}_{1}$ involves the product ${\prod}_{t=1}^{3}(\partial {\bm{h}}_{t+1}/\partial {\bm{h}}_{t})$ by the chain rule. However, if the sequence is long, the product will be the multiplication of many Jacobian matrices, which usually results in exponentially large or small singular values. To alleviate this issue, in practice, the forward pass and backward pass are implemented in a shorter sliding window $\{{t}_{1},{t}_{1}+1,\mathrm{\dots},{t}_{2}\}$, instead of the full sequence $\{1,2,\mathrm{\dots},T\}$. Though effective in some cases, this technique alone does not fully address the issue of longterm dependency.
3.2.2 GRUs and LSTM
There are two improved variants that alleviate the above issue: gated recurrent units (GRUs) [26] and long shortterm memory (LSTM) [54].

•
A GRU refines the recursive formula (13) by introducing gates, which are vectors of the same length as ${\bm{h}}_{t}$. The gates, which take values in $[0,1]$ elementwise, multiply with ${\bm{h}}_{t1}$ elementwise and determine how much they keep the old hidden states.

•
An LSTM similarly uses gates in the recursive formula. In addition to ${\bm{h}}_{t}$, an LSTM maintains a cell state, which takes values in $\mathbb{R}$ elementwise and are analogous to counters.
Here we only discuss LSTM in detail. Denote by $\odot $ the elementwise multiplication. We have a recursive formula in replace of (13):
$\left(\begin{array}{c}\hfill {\bm{i}}_{t}\hfill \\ \hfill {\bm{f}}_{t}\hfill \\ \hfill {\bm{o}}_{t}\hfill \\ \hfill {\bm{g}}_{t}\hfill \end{array}\right)$  $=\left(\begin{array}{c}\hfill \bm{\sigma}\hfill \\ \hfill \bm{\sigma}\hfill \\ \hfill \bm{\sigma}\hfill \\ \hfill \text{\mathbf{t}\mathbf{a}\mathbf{n}\mathbf{h}}\hfill \end{array}\right)\mathbf{W}\left(\begin{array}{c}\hfill {\bm{h}}_{t1}\hfill \\ \hfill {\bm{x}}_{t}\hfill \\ \hfill 1\hfill \end{array}\right),$  
${\bm{c}}_{t}$  $={\bm{f}}_{t}\odot {\bm{c}}_{t1}+{\bm{i}}_{t}\odot {\bm{g}}_{t},$  
${\bm{h}}_{t}$  $={\bm{o}}_{t}\odot \text{\mathbf{t}\mathbf{a}\mathbf{n}\mathbf{h}}({\bm{c}}_{t}),$ 
where $\mathbf{W}$ is a big weight matrix with appropriate dimensions. The cell state vector ${\bm{c}}_{t}$ carries information of the sequence (e.g., singular/plural form in a sentence). The forget gate ${\bm{f}}_{t}$ determines how much the values of ${\bm{c}}_{t1}$ are kept for time $t$, the input gate ${\bm{i}}_{t}$ controls the amount of update to the cell state, and the output gate ${\bm{o}}_{t}$ gives how much ${\bm{c}}_{t}$ reveals to ${\bm{h}}_{t}$. Ideally, the elements of these gates have nearly binary values. For example, an element of ${\bm{f}}_{t}$ being close to $1$ may suggest the presence of a feature in the sequence data. Similar to the skip connections in residual nets, the cell state ${\bm{c}}_{t}$ has an additive recursive formula, which helps backpropagation and thus captures longrange dependencies.
3.2.3 Multilayer RNNs
Multilayer RNNs are generalization of the onehiddenlayer RNN discussed above. Figure 9 shows a vanilla RNN with two hidden layers. In place of (13), the recursive formula for an RNN with $L$ hidden layers now reads
$${\bm{h}}_{t}^{\mathrm{\ell}}=\text{\mathbf{t}\mathbf{a}\mathbf{n}\mathbf{h}}\left[{\mathbf{W}}^{\mathrm{\ell}}\left(\begin{array}{c}\hfill {\bm{h}}_{t}^{\mathrm{\ell}1}\hfill \\ \hfill {\bm{h}}_{t1}^{\mathrm{\ell}}\hfill \\ \hfill 1\hfill \end{array}\right)\right],\text{for all}\mathrm{\ell}\in [L],{\bm{h}}_{t}^{0}\triangleq {\bm{x}}_{t}.$$ 
Note that a multilayer RNN has two dimensions: the sequence length $T$ and depth $L$. Two special cases are the feedforward neural nets (where $T=1$) introduced in Section 2, and RNNs with one hidden layer (where $L=1$). Multilayer RNNs usually do not have very large depth (e.g., $2$–$5$), since $T$ is already very large.
Finally, we remark that CNNs, RNNs, and other neural nets can be easily combined to tackle tasks that involve different sources of input data. For example, in image captioning, the images are first processed through a CNN, and then the highlevel features are fed into an RNN as inputs. Theses neural nets combined together form a large computational graph, so they can be trained using backpropagation. This generic training method provides much flexibility in various applications.
3.3 Modules
Deep neural nets are essentially composition of many nonlinear functions. A component function may be designed to have specific properties in a given task, and it can be itself resulted from composing a few simpler functions. In LSTM, we have seen that the building block consists of several intermediate variables, including cell states and forget gates that can capture longterm dependency and alleviate numerical issues.
This leads to the idea of designing modules for building more complex neural net models. Desirable modules usually have low computational costs, alleviate numerical issues in training, and lead to good statistical accuracy. Since modules and the resulting neural net models form computational graphs, training follows the same principle briefly described in Section 2.
Here, we use the examples of Inception and skip connections to illustrate the ideas behind modules. Figure 10(a) is an example of “Inception” modules used in GoogleNet [122]. As before, all the convolutional layers are followed by the ReLU activation function. The concatenation of information from filters with different sizes give the model great flexibility to capture spatial information. Note that $1\times 1$ filters is an $1\times 1\times {d}_{3}$ tensor (where ${d}_{3}$ is the number of feature maps), so its convolutional operation does not interact with other spatial coordinates, only serving to aggregate information from different feature maps at the same coordinate. This reduces the number of parameters and speeds up the computation. Similar ideas appear in other work [78, 57].
(a) “Inception” module  (b) Skip connections 
Another module, usually called skip connections, is widely used to alleviate numerical issues in very deep neural nets, with additional benefits in optimization efficiency and statistical accuracy. Training very deep neural nets are generally more difficult, but the introduction of skip connections in residual networks [50, 51] has greatly eased the task.
The high level idea of skip connections is to add an identity map to an existing nonlinear function. Let $\mathbf{F}(\bm{x})$ be an arbitrary nonlinear function represented by a (fragment of) neural net, then the idea of skip connections is simply replacing $\mathbf{F}(\bm{x})$ with $\bm{x}+\mathbf{F}(\bm{x})$. Figure 10(b) shows a wellknown structure from residual networks [50]—for every two layers, an identity map is added:
$$\bm{x}\u27fc\bm{\sigma}(\bm{x}+\mathbf{F}(\bm{x}))=\bm{\sigma}(\bm{x}+{\mathbf{W}}^{\prime}\bm{\sigma}(\mathbf{W}\bm{x}+\bm{b})+{\bm{b}}^{\prime}),$$  (14) 
where $\bm{x}$ can be hidden nodes from any layer and $\mathbf{W},{\mathbf{W}}^{\prime},\bm{b},{\bm{b}}^{\prime}$ are corresponding parameters. By repeating (namely composing) this structure throughout all layers, [50, 51] are able to train neural nets with hundreds of layers easily, which overcomes wellobserved training difficulties in deep neural nets. Moreover, deep residual networks also improve statistical accuracy, as the classification error on ImageNet challenge was reduced by $46\%$ from 2014 to 2015. As a side note, skip connections can be used flexibly. They are not restricted to the form in (14), and can be used between any pair of layers $\mathrm{\ell},{\mathrm{\ell}}^{\prime}$ [55].
4 Deep unsupervised learning
In supervised learning, given labelled training set $\{({y}_{i},{\mathbf{x}}_{i})\}$, we focus on discriminative models, which essentially represents $\mathbb{P}(y\bm{x})$ by a deep neural net $f(\bm{x};\bm{\theta})$ with parameters $\bm{\theta}$. Unsupervised learning, in contrast, aims at extracting information from unlabeled data $\{{\bm{x}}_{i}\}$, where the labels $\{{y}_{i}\}$ are absent. In regard to this information, it can be a lowdimensional embedding of the data $\{{\bm{x}}_{i}\}$ or a generative model with latent variables to approximate the distribution ${\mathbb{P}}_{\mathbf{X}}(\bm{x})$. To achieve these goals, we introduce two popular unsupervised deep leaning models, namely, autoencoders and generative adversarial networks (GANs). The first one can be viewed as a dimension reduction technique, and the second as a density estimation method. DNNs are the key elements for both of these two models.
4.1 Autoencoders
Recall that in dimension reduction, the goal is to reduce the dimensionality of the data and at the same time preserve its salient features. In particular, in principal component analysis (PCA), the goal is to embed the data ${\{{\bm{x}}_{i}\}}_{1\le i\le n}$ into a lowdimensional space via a linear function $\bm{f}$ such that maximum variance can be explained. Equivalently, we want to find linear functions $\bm{f}:{\mathbb{R}}^{d}\to {\mathbb{R}}^{k}$ and $\bm{g}:{\mathbb{R}}^{k}\to {\mathbb{R}}^{d}$ ($k\le d$) such that the difference between ${\bm{x}}_{i}$ and $\bm{g}(\bm{f}({\bm{x}}_{i}))$ is minimized. Formally, we let
$$\bm{f}\left(\bm{x}\right)={\bm{W}}_{f}\bm{x}\triangleq \bm{h}\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}\bm{g}\left(\bm{h}\right)={\bm{W}}_{g}\bm{h},\text{where}\mathit{\hspace{1em}}{\bm{W}}_{f}\in {\mathbb{R}}^{k\times d}\text{and}{\bm{W}}_{g}\in {\mathbb{R}}^{d\times k}.$$ 
Here, for simplicity, we assume that the intercept/bias terms for $\bm{f}$ and $\bm{g}$ are zero. Then, PCA amounts to minimizing the quadratic loss function
$${\text{minimize}}_{{\bm{W}}_{f},{\bm{W}}_{g}}\mathit{\hspace{1em}\hspace{1em}}\frac{1}{n}\sum _{i=1}^{n}{\parallel {\bm{x}}_{i}{\bm{W}}_{f}{\bm{W}}_{g}{\bm{x}}_{i}\parallel}_{2}^{2}.$$  (15) 
It is the same as minimizing ${\parallel \mathbf{X}\mathrm{\mathbf{W}\mathbf{X}}\parallel}_{\mathrm{F}}^{2}$ subject to $\mathrm{rank}(\mathbf{W})\le k$, where $\mathbf{X}\in {\mathbb{R}}^{p\times n}$ is the design matrix. The solution is given by the singular value decomposition of $\mathbf{X}$ [44, Thm. 2.4.8], which is exactly what PCA does. It turns out that PCA is a special case of autoencoders, which is often known as the undercomplete linear autoencoder.
More broadly, autoencoders are neural network models for (nonlinear) dimension reduction, which generalize PCA. An autoencoder has two key components, namely, the encoder function $\bm{f}(\cdot )$, which maps the input $\bm{x}\in {\mathbb{R}}^{d}$ to a hidden code/representation $\bm{h}\triangleq \bm{f}(\bm{x})\in {\mathbb{R}}^{k}$, and the decoder function $\bm{g}(\cdot )$, which maps the hidden representation $\bm{h}$ to a point $\bm{g}(\bm{h})\in {\mathbb{R}}^{d}$. Both functions can be multilayer neural networks as (3). See Figure 11 for an illustration of autoencoders. Let $\mathcal{L}({\bm{x}}_{1},{\bm{x}}_{2})$ be a loss function that measures the difference between ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$ in ${\mathbb{R}}^{d}$. Similar to PCA, an autoencoder is used to find the encoder $\bm{f}$ and decoder $\bm{g}$ such that $\mathcal{L}(\bm{x},\bm{g}(\bm{f}(\bm{x})))$ is as small as possible. Mathematically, this amounts to solving the following minimization problem
$${\text{minimize}}_{\bm{f},\bm{g}}\mathit{\hspace{1em}}\frac{1}{n}\sum _{i=1}^{n}\mathcal{L}({\bm{x}}_{i},\bm{g}\left({\bm{h}}_{i}\right))\mathit{\hspace{1em}}\text{with}\mathit{\hspace{1em}}{\bm{h}}_{i}=\bm{f}\left({\bm{x}}_{i}\right),\text{for all}i\in [n].$$  (16) 
One needs to make structural assumptions on the functions $\bm{f}$ and $\mathbf{g}$ in order to find useful representations of the data, which leads to different types of autoencoders. Indeed, if no assumption is made, choosing $\bm{f}$ and $\bm{g}$ to be identity functions clearly minimizes the above optimization problem. To avoid this trivial solution, one natural way is to require that the encoder $f$ maps the data onto a space with a smaller dimension, i.e., $$. This is the undercomplete autoencoder that includes PCA as a special case. There are other structured autoencoders which add desired properties to the model such as sparsity or robustness, mainly through regularization terms. Below we present two other common types of autoencoders.

•
Sparse autoencoders. One may believe that the dimension $k$ of the hidden code ${\bm{h}}_{i}$ is larger than the input dimension $d$, and that ${\bm{h}}_{i}$ admits a sparse representation. As with LASSO [125] or SCAD [36], one may add a regularization term to the reconstruction loss $\mathcal{L}$ in (16) to encourage sparsity [98]. A sparse autoencoder solves
$${\text{min}}_{\bm{f},\bm{g}}\underset{\text{loss}}{\underset{\u23df}{\frac{1}{n}\sum _{i=1}^{n}\mathcal{L}({\bm{x}}_{i},\bm{g}\left({\bm{h}}_{i}\right))}}+\underset{\text{regularizer}}{\underset{\u23df}{\lambda {\parallel {\bm{h}}_{i}\parallel}_{1}}}\mathit{\hspace{1em}}\text{with}\mathit{\hspace{1em}}{\bm{h}}_{i}=\bm{f}\left({\bm{x}}_{i}\right),\text{for all}i\in [n].$$ This is similar to dictionary learning, where one aims at finding a sparse representation of input data on an overcomplete basis. Due to the imposed sparsity, the model can potentially learn useful features of the data.

•
Denoising autoencoders. One may hope that the model is robust to noise in the data: even if the input data ${\bm{x}}_{i}$ are corrupted by small noise ${\bm{\xi}}_{i}$ or miss some components (the noise level or the missing probability is typically small), an ideal autoencoder should faithfully recover the original data. A denoising autoencoder [127] achieves this robustness by explicitly building a noisy data ${\stackrel{~}{\bm{x}}}_{i}={\bm{x}}_{i}+{\bm{\xi}}_{i}$ as the new input, and then solves an optimization problem similar to (16) where $\mathcal{L}({\bm{x}}_{i},\bm{g}\left({\bm{h}}_{i}\right))$ is replaced by $\mathcal{L}({\bm{x}}_{i},\bm{g}\left(\bm{f}({\stackrel{~}{\bm{x}}}_{i})\right))$. A denoising autoencoder encourages the encoder/decoder to be stable in the neighborhood of an input, which is generally a good statistical property. An alternative way could be constraining $f$ and $g$ in the optimization problem, but that would be very difficult to optimize. Instead, sampling by adding small perturbations in the input provides a simple implementation. We shall see similar ideas in Section 6.3.3.
4.2 Generative adversarial networks
Given unlabeled data ${\{{\bm{x}}_{i}\}}_{1\le i\le n}$, density estimation aims to estimate the underlying probability density function ${\mathbb{P}}_{\bm{X}}$ from which the data is generated. Both parametric and nonparametric estimators [114] have been proposed and studied under various assumptions on the underlying distribution. Different from these classical density estimators, where the density function is explicitly defined in relatively low dimension, generative adversarial networks (GANs) [46] can be categorized as an implicit density estimator in much higher dimension. The reasons are twofold: (1) GANs put more emphasis on sampling from the distribution ${\mathbb{P}}_{\bm{X}}$ than estimation; (2) GANs define the density estimation implicitly through a source distribution ${\mathbb{P}}_{\bm{Z}}$ and a generator function $g(\cdot )$, which is usually a deep neural network. We introduce GANs from the perspective of sampling from ${\mathbb{P}}_{\bm{X}}$ and later we will generalize the vanilla GANs using its relation to density estimators.
4.2.1 Sampling view of GANs
Suppose the data ${\{{\bm{x}}_{i}\}}_{1\le i\le n}$ at hand are all real images, and we want to generate new natural images. With this goal in mind, GAN models a zerosum game between two players, namely, the generator $\mathcal{G}$ and the discriminator $\mathcal{D}$. The generator $\mathcal{G}$ tries to generate fake images akin to the true images ${\{{\bm{x}}_{i}\}}_{1\le i\le n}$ while the discriminator $\mathcal{D}$ aims at differentiating the fake ones from the true ones. Intuitively, one hopes to learn a generator $\mathcal{G}$ to generate images where the best discriminator $\mathcal{D}$ cannot distinguish. Therefore the payoff is higher for the generator $\mathcal{G}$ if the probability of the discriminator $\mathcal{D}$ getting wrong is higher, and correspondingly the payoff for the discriminator correlates positively with its ability to tell wrong from truth.
Mathematically, the generator $\mathcal{G}$ consists of two components, an source distribution ${\mathbb{P}}_{\bm{Z}}$ (usually a standard multivariate Gaussian distribution with hundreds of dimensions) and a function $\mathbf{g}(\cdot )$ which maps a sample $\bm{z}$ from ${\mathbb{P}}_{\bm{Z}}$ to a point $\mathbf{g}(\bm{z})$ living in the same space as $\bm{x}$. For generating images, $\mathbf{g}(\bm{z})$ would be a 3D tensor. Here $\mathbf{g}(\bm{z})$ is the fake sample generated from $\mathcal{G}$. Similarly the discriminator $\mathcal{D}$ is composed of one function which takes an image $\bm{x}$ (real or fake) and return a number $d(\bm{x})\in [0,1]$, the probability of $\bm{x}$ being a real sample from ${\mathbb{P}}_{\bm{X}}$ or not. Oftentimes, both the generating function $\mathbf{g}(\cdot )$ and the discriminating function $d(\cdot )$ are realized by deep neural networks, e.g., CNNs introduced in Section 3.1. See Figure 12 for an illustration for GANs. Denote ${\bm{\theta}}_{\mathcal{G}}$ and ${\bm{\theta}}_{\mathcal{D}}$ the parameters in $\mathbf{g}(\cdot )$ and $d(\cdot )$, respectively. Then GAN tries to solve the following minmax problem:
$$\underset{{\bm{\theta}}_{\mathcal{G}}}{\mathrm{min}}\underset{{\bm{\theta}}_{\mathcal{D}}}{\mathrm{max}}\mathit{\hspace{1em}\hspace{1em}}{\mathbb{E}}_{\bm{x}\sim {\mathbb{P}}_{\bm{X}}}\left[\mathrm{log}\left(d\left(\bm{x}\right)\right)\right]+{\mathbb{E}}_{\bm{z}\sim {\mathbb{P}}_{\bm{Z}}}\left[\mathrm{log}\left(1d\left(\mathbf{g}\left(\bm{z}\right)\right)\right)\right].$$  (17) 
Recall that $d(\bm{x})$ models the belief / probability that the discriminator thinks that $\bm{x}$ is a true sample. Fix the parameters ${\bm{\theta}}_{\mathcal{G}}$ and hence the generator $\mathcal{G}$ and consider the inner maximization problem. We can see that the goal of the discriminator is to maximize its ability of differentiation. Similarly, if we fix ${\bm{\theta}}_{\mathcal{D}}$ (and hence the discriminator), the generator tries to generate more realistic images $\mathbf{g}(\bm{z})$ to fool the discriminator.
4.2.2 Density estimation view of GANs
Let us now take a densityestimation view of GANs. Fixing the source distribution ${\mathbb{P}}_{\bm{Z}}$, any generator $\mathcal{G}$ induces a distribution ${\mathbb{P}}_{\mathcal{G}}$ over the space of images. Removing the restrictions on $d(\cdot )$, one can then rewrite (17) as
$$\underset{{\mathbb{P}}_{\mathcal{G}}}{\mathrm{min}}\underset{d(\cdot )}{\mathrm{max}}\mathit{\hspace{1em}\hspace{1em}}{\mathbb{E}}_{\bm{x}\sim {\mathbb{P}}_{\bm{X}}}\left[\mathrm{log}\left(d\left(\bm{x}\right)\right)\right]+{\mathbb{E}}_{\bm{x}\sim {\mathbb{P}}_{\mathcal{G}}}\left[\mathrm{log}\left(1d\left(\bm{x}\right)\right)\right].$$  (18) 
Observe that the inner maximization problem is solved by the likelihood ratio, i.e.
$${d}^{*}\left(\bm{x}\right)=\frac{{\mathbb{P}}_{\bm{X}}\left(\bm{x}\right)}{{\mathbb{P}}_{\bm{X}}\left(\bm{x}\right)+{\mathbb{P}}_{\mathcal{G}}\left(\bm{x}\right)}.$$ 
As a result, (18) can be simplified as
$$\underset{{\mathbb{P}}_{\mathcal{G}}}{\mathrm{min}}\mathit{\hspace{1em}\hspace{1em}}\text{JS}\left({\mathbb{P}}_{\bm{X}}\parallel {\mathbb{P}}_{\mathcal{G}}\right),$$  (19) 
where $\text{JS}(\cdot \parallel \cdot )$ denotes the Jensen–Shannon divergence between two distributions
$$\text{JS}\left({\mathbb{P}}_{\bm{X}}\parallel {\mathbb{P}}_{\mathcal{G}}\right)=\frac{1}{2}\text{KL}\left({\mathbb{P}}_{\bm{X}}\parallel \frac{{\mathbb{P}}_{\bm{X}}+{\mathbb{P}}_{\mathcal{G}}}{2}\right)+\frac{1}{2}\text{KL}\left({\mathbb{P}}_{\mathcal{G}}\parallel \frac{{\mathbb{P}}_{\bm{X}}+{\mathbb{P}}_{\mathcal{G}}}{2}\right).$$ 
In words, the vanilla GAN (17) seeks a density ${\mathbb{P}}_{\mathcal{G}}$ that is closest to ${\mathbb{P}}_{\bm{X}}$ in terms of the Jensen–Shannon divergence. This view allows to generalize GANs to other variants, by changing the distance metric. Examples include fGAN [90], Wasserstein GAN (WGAN) [6], MMD GAN [75], etc. We single out the Wasserstein GAN (WGAN) [6] to introduce due to its popularity. As the name suggests, it minimizes the Wasserstein distance between ${\mathbb{P}}_{\bm{X}}$ and ${\mathbb{P}}_{\mathcal{G}}$:
$$\underset{{\bm{\theta}}_{\mathcal{G}}}{\mathrm{min}}\mathit{\hspace{1em}}\text{WS}\left({\mathbb{P}}_{\bm{X}}\parallel {\mathbb{P}}_{\mathcal{G}}\right)=\underset{{\bm{\theta}}_{\mathcal{G}}}{\mathrm{min}}\underset{f:f\text{1Lipschitz}}{sup}{\mathbb{E}}_{\bm{x}\sim {\mathbb{P}}_{\bm{X}}}\left[f\left(\bm{x}\right)\right]{\mathbb{E}}_{\bm{x}\sim {\mathbb{P}}_{\mathcal{G}}}\left[f\left(\bm{x}\right)\right],$$  (20) 
where $f(\cdot )$ is taken over all Lipschitz functions with coefficient 1. Comparing WGAN (20) with the original formulation of GAN (17), one finds that the Lipschitz function $f$ in (20) corresponds to the discriminator $\mathcal{D}$ in (17) in the sense that they share similar objectives to differentiate the true distribution ${\mathbb{P}}_{\bm{X}}$ from the fake one ${\mathbb{P}}_{\mathcal{G}}$. In the end, we would like to mention that GANs are more difficult to train than supervised deep learning models such as CNNs [109]. Apart from the training difficulty, how to evaluate GANs objectively and effectively is an ongoing research.
5 Representation power: approximation theory
Having seen the building blocks of deep learning models in the previous sections, it is natural to ask: what is the benefits of composing multiple layers of nonlinear functions. In this section, we address this question from a approximation theoretical point of view. Mathematically, letting $\mathscr{H}$ be the space of functions representable by neural nets (NNs), how well can a function $f$ (with certain properties) be approximated by functions in $\mathscr{H}$. We first revisit universal approximation theories, which are mostly developed for shallow neural nets (neural nets with a single hidden layer), and then provide recent results that demonstrate the benefits of depth in neural nets. Other notable works include KolmogorovArnold superposition theorem [7, 119], and circuit complexity for neural nets [91].
5.1 Universal approximation theory for shallow NNs
The universal approximation theories study the approximation of $f$ in a space $\mathcal{F}$ by a function represented by a onehiddenlayer neural net
$$g(\bm{x})=\sum _{j=1}^{N}{c}_{j}{\sigma}_{*}({\bm{w}}_{j}^{\top}\bm{x}{b}_{j}),$$  (21) 
where ${\sigma}_{*}:\mathbb{R}\to \mathbb{R}$ is certain activation function and $N$ is the number of hidden units in the neural net. For different space $\mathcal{F}$ and activation function ${\sigma}_{*}$, there are upper bounds and lower bounds on the approximation error $\parallel fg\parallel $. See [93] for a comprehensive overview. Here we present representative results.
First, as $N\to \mathrm{\infty}$, any continuous function $f$ can be approximated by some $g$ under mild conditions. Loosely speaking, this is because each component ${\sigma}_{*}({\bm{w}}_{j}^{\top}\bm{x}{b}_{j})$ behaves like a basis function and functions in a suitable space $\mathcal{F}$ admits a basis expansion. Given the above heuristics, the next natural question is: what is the rate of approximation for a finite $N$?
Let us restrict the domain of $\bm{x}$ to a unit ball ${B}^{d}$ in ${\mathbb{R}}^{d}$. For $p\in [1,\mathrm{\infty})$ and integer $m\ge 1$, consider the ${L}^{p}$ space and the Sobolev space with standard norms
${\parallel f\parallel}_{p}={\left[{\displaystyle {\int}_{{B}^{n}}}{g(\bm{x})}^{p}\mathit{d}\bm{x}\right]}^{1/p},{\parallel f\parallel}_{m,p}={\left[{\displaystyle \sum _{0\le \bm{k}\le m}}{\parallel {D}^{\bm{k}}f\parallel}_{p}^{p}\right]}^{1/p},$ 
where ${D}^{\bm{k}}f$ denotes partial derivatives indexed by $\bm{k}\in {\mathbb{Z}}_{+}^{d}$. Let $\mathcal{F}\triangleq {\mathcal{F}}_{p}^{m}$ be the space of functions $f$ in the Sobolev space with ${\parallel f\parallel}_{m,p}\le 1$. Note that functions in $\mathcal{F}$ have bounded derivatives up to $m$th order, and that smoothness of functions is controlled by $m$ (larger $m$ means smoother). Denote by ${\mathscr{H}}_{N}$ the space of functions with the form (21). The following general upper bound is due to [85].
Theorem 1 (Theorem 2.1 in [85]).
Assume ${\sigma}_{\mathrm{*}}\mathrm{:}\mathrm{R}\mathrm{\to}\mathrm{R}$ is such that ${\sigma}_{\mathrm{*}}$ has arbitrary order derivatives in an open interval $I$, and that ${\sigma}_{\mathrm{*}}$ is not a polynomial on $I$. Then, for any $p\mathrm{\in}\mathrm{[}\mathrm{1}\mathrm{,}\mathrm{\infty}\mathrm{)}$, $d\mathrm{\ge}\mathrm{2}$, and integer $m\mathrm{\ge}\mathrm{1}$,
$$\underset{f\in {\mathcal{F}}_{p}^{m}}{sup}\underset{g\in {\mathscr{H}}_{N}}{inf}{\parallel fg\parallel}_{p}\le {C}_{d,m,p}{N}^{m/d},$$ 
where ${C}_{d\mathrm{,}m\mathrm{,}p}$ is independent of $N$, the number of hidden units.
In the above theorem, the condition on ${\sigma}_{*}(\cdot )$ is mainly technical. This upper bound is useful when the dimension $d$ is not large. It clearly implies that the onehiddenlayer neural net is able to approximate any smooth function with enough hidden units. However, it is unclear how to find a good approximator $g$; nor do we have control over the magnitude of the parameters (huge weights are impractical). While increasing the number of hidden units $N$ leads to better approximation, the exponent $m/d$ suggests the presence of the curse of dimensionality. The following (nearly) matching lower bound is stated in [80].
Theorem 2 (Theorem 5 in [80]).
Let $p\mathrm{\ge}\mathrm{1}$, $m\mathrm{\ge}\mathrm{1}$ and $N\mathrm{\ge}\mathrm{2}$. If the activation function is the standard sigmoid function $\sigma \mathit{}\mathrm{(}t\mathrm{)}\mathrm{=}{\mathrm{(}\mathrm{1}\mathrm{+}{e}^{\mathrm{}t}\mathrm{)}}^{\mathrm{}\mathrm{1}}$, then
$$\underset{f\in {\mathcal{F}}_{p}^{m}}{sup}\underset{g\in {\mathscr{H}}_{N}}{inf}{\parallel fg\parallel}_{p}\ge {C}_{d,m,p}^{\prime}{(N\mathrm{log}N)}^{m/d},$$  (22) 
where ${C}_{d\mathrm{,}m\mathrm{,}p}^{\mathrm{\prime}}$ is independent of $N$.
Results for other activation functions are also obtained by [80]. Moreover, the term $\mathrm{log}N$ can be removed if we assume an additional continuity condition [85].
For the natural space ${\mathcal{F}}_{p}^{m}$ of smooth functions, the exponential dependence on $d$ in the upper and lower bounds may look unappealing. However, [12] showed that for a different function space, there is a good dimensionfree approximation by the neural nets. Suppose that a function $f:{\mathbb{R}}^{d}\mapsto \mathbb{R}$ has a Fourier representation
$$f(\bm{x})={\int}_{{\mathbb{R}}^{d}}{e}^{i\u27e8\bm{\omega},\bm{x}\u27e9}\stackrel{~}{f}(\bm{\omega})\mathit{d}\bm{\omega},$$  (23) 
where $\stackrel{~}{f}(\bm{\omega})\in \u2102$. Assume that $f(\mathrm{\U0001d7ce})=0$ and that the following quantity is finite
$${C}_{f}={\int}_{{\mathbb{R}}^{d}}{\parallel \bm{\omega}\parallel}_{2}\stackrel{~}{f}(\bm{\omega})\mathit{d}\bm{\omega}.$$  (24) 
[12] uncovers the following dimensionfree approximation guarantee.
Theorem 3 (Proposition 1 in [12]).
Fix a $C\mathrm{>}\mathrm{0}$ and an arbitrary probability measure $\mu $ on the unit ball ${B}^{d}$ in ${\mathrm{R}}^{d}$. For every function $f$ with ${C}_{f}\mathrm{\le}C$ and every $N\mathrm{\ge}\mathrm{1}$, there exists some $g\mathrm{\in}{\mathrm{H}}_{N}$ such that
$${\left[{\int}_{{B}^{d}}{(f(\bm{x})g(\bm{x}))}^{2}\mu (d\bm{x})\right]}^{1/2}\le \frac{2C}{\sqrt{N}}.$$ 
Moreover, the coefficients of $g$ may be restricted to satisfy ${\mathrm{\sum}}_{j\mathrm{=}\mathrm{1}}^{N}\mathrm{}{c}_{j}\mathrm{}\mathrm{\le}\mathrm{2}\mathit{}C$.
The upper bound is now independent of the dimension $d$. However, ${C}_{f}$ may implicitly depend on $d$, as the formula in (24) involves an integration over ${\mathbb{R}}^{d}$ (so for some functions ${C}_{f}$ may depend exponentially on $d$). Nevertheless, this theorem does characterize an interesting function space with an improved upper bound. Details of the function space are discussed by [12]. This theorem can be generalized; see [81] for an example.
To help understand why a dimensionalityfree approximation holds, let us appeal to a heuristic argument given by Monte Carlo simulations. It is wellknown that Monte Carlo approximation errors are independent of dimensionality in evaluation of highdimensional integrals. Let us generate ${\{{\bm{\omega}}_{j}\}}_{1\le j\le N}$ randomly from a given density $p(\cdot )$ in ${\mathbb{R}}^{d}$. Consider the approximation to (23) by
$${g}_{N}(\bm{x})=\frac{1}{N}\sum _{j=1}^{N}{c}_{j}{e}^{i\u27e8{\bm{\omega}}_{j},\bm{x}\u27e9},{c}_{j}=\frac{\stackrel{~}{f}({\bm{\omega}}_{j})}{p({\bm{\omega}}_{j})}.$$  (25) 
Then, ${g}_{N}(\bm{x})$ is a onehiddenlayer neural network with $N$ units and the sinusoid activation function. Note that $\mathbb{E}{g}_{N}(\bm{x})=f(\bm{x})$, where the expectation is taken with respect to randomness $\{{\bm{\omega}}_{j}\}$. Now, by independence, we have
$$\mathbb{E}{({g}_{N}(\bm{x})f(\bm{x}))}^{2}=\frac{1}{N}\mathrm{Var}({c}_{j}{e}^{i\u27e8{\bm{\omega}}_{j},\bm{x}\u27e9})\le \frac{1}{N}\mathbb{E}{c}_{j}^{2},$$ 
if $$. Therefore, the rate is independent of the dimensionality $d$, though the constant can be.
5.2 Approximation theory for multilayer NNs
The approximation theory for multilayer neural nets is less understood compared with neural nets with one hidden layer. Driven by the success of deep learning, there are many recent papers focusing on expressivity of deep neural nets. As studied by [124, 35, 84, 94, 15, 110, 77, 102], deep neural nets excel at representing composition of functions. This is perhaps not surprising, since deep neural nets are themselves defined by composing layers of functions. Nevertheless, it points to a new territory rarely studied in statistics before. Below we present a result based on [77, 102].
Suppose that the inputs $\bm{x}$ have a bounded domain ${[1,1]}^{d}$ for simplicity. As before, let ${\sigma}_{*}:\mathbb{R}\to \mathbb{R}$ be a generic function, and ${\bm{\sigma}}_{*}={({\sigma}_{*},\mathrm{\cdots},{\sigma}_{*})}^{\top}$ be elementwise application of ${\sigma}_{*}$. Consider a neural net which is similar to (3) but with scaler output: $g(\bm{x})={\mathbf{W}}_{\mathrm{\ell}}{\bm{\sigma}}_{*}(\mathrm{\cdots}{\bm{\sigma}}_{*}({\mathbf{W}}_{2}{\bm{\sigma}}_{*}({\mathbf{W}}_{1}\bm{x}))\mathrm{\cdots})$. A unit or neuron refers to an element of vectors ${\bm{\sigma}}_{*}({\mathbf{W}}_{k}\mathrm{\cdots}{\bm{\sigma}}_{*}({\mathbf{W}}_{2}{\bm{\sigma}}_{*}({\mathbf{W}}_{1}\bm{x}))\mathrm{\cdots})$ for any $k=1,\mathrm{\dots},\mathrm{\ell}1$. For a multivariate polynomial $p$, define ${m}_{k}(p)$ to be the smallest integer such that, for any $\u03f5>0$, there exists a neural net $g(\bm{x})$ satisfying $$, with $k$ hidden layers (i.e., $\mathrm{\ell}=k+1$) and no more than ${m}_{k}(p)$ neurons in total. Essentially, ${m}_{k}(p)$ is the minimum number of neurons required to approximate $p$ arbitrarily well.
Theorem 4 (Theorem 4.1 in [102]).
Let $p\mathit{}\mathrm{(}\mathbf{x}\mathrm{)}$ be a monomial ${x}_{\mathrm{1}}^{{r}_{\mathrm{1}}}\mathit{}{x}_{\mathrm{2}}^{{r}_{\mathrm{2}}}\mathit{}\mathrm{\cdots}\mathit{}{x}_{d}^{{r}_{d}}$ with $q\mathrm{=}{\mathrm{\sum}}_{j\mathrm{=}\mathrm{1}}^{d}{r}_{j}$. Suppose that ${\sigma}_{\mathrm{*}}$ has derivatives of order $\mathrm{2}\mathit{}q$ at the origin, and that they are nonzero. Then,
(i) ${m}_{\mathrm{1}}\mathit{}\mathrm{(}p\mathrm{)}\mathrm{=}{\mathrm{\prod}}_{j\mathrm{=}\mathrm{1}}^{d}\mathrm{(}{r}_{j}\mathrm{+}\mathrm{1}\mathrm{)}$;
(ii) ${\mathrm{min}}_{k}\mathit{}{m}_{k}\mathit{}\mathrm{(}p\mathrm{)}\mathrm{\le}{\mathrm{\sum}}_{j\mathrm{=}\mathrm{1}}^{d}\mathrm{\left(}\mathrm{7}\mathit{}\mathrm{\lceil}{\mathrm{log}}_{\mathrm{2}}\mathit{}\mathrm{(}{r}_{j}\mathrm{)}\mathrm{\rceil}\mathrm{+}\mathrm{4}\mathrm{\right)}$.
This theorem reveals a sharp distinction between shallow networks (one hidden layer) and deep networks. To represent a monomial function, a shallow network requires exponentially many neurons in terms of the dimension $d$, whereas linearly many neurons suffice for a deep network (with bounded ${r}_{j}$). The exponential dependence on $d$, as shown in Theorem 4(i), is resonant with the curse of dimensionality widely seen in many fields; see [30]. One may ask: how does depth help? Depth circumvents this issue, at least for certain functions, by allowing us to represent function composition efficiently. Indeed, Theorem 4(ii) offers a nice result with clear intuitions: it is known that the product of two scalar inputs can be represented using $4$ neurons [77], so by composing multiple products, we can express monomials with $O(d)$ neurons.
Recent advances in nonparametric regressions also support the idea that deep neural nets excel at representing composition of functions [15, 110]. In particular, [15] considered the nonparametric regression setting where we want to estimate a function ${\widehat{f}}_{n}(\bm{x})$ from i.i.d. data ${\mathcal{D}}_{n}={\{({y}_{i},{\bm{x}}_{i})\}}_{1\le i\le n}$. If the true regression function $f(\bm{x})$ has certain hierarchical structure with intrinsic dimensionality^{6}^{6} 6 Roughly speaking, the true regression function can be represented by a tree where each node has at most ${d}^{*}$ children. See [15] for the precise definition. ${d}^{*}$, then the error
$${\mathbb{E}}_{{\mathcal{D}}_{n}}{\mathbb{E}}_{\bm{x}}{\left{\widehat{f}}_{n}(\bm{x})f(\bm{x})\right}^{2}$$ 
has an optimal minimax convergence rate $O({n}^{\frac{2q}{2q+{d}^{*}}})$, rather than the usual rate $O({n}^{\frac{2q}{2q+d}})$ that depends on the ambient dimension $d$. Here $q$ is the smoothness parameter. This provides another justification for deep neural nets: if data are truly hierarchical, then the quality of approximators by deep neural nets depends on the intrinsic dimensionality, which avoids the curse of dimensionality.
We point out that the approximation theory for deep learning is far from complete. For example, in Theorem 4, the condition on ${\sigma}_{*}$ excludes the widely used ReLU activation function, there are no constraints on the magnitude of the weights (so they can be unreasonably large).
6 Training deep neural nets
The existence of a good function approximator in the NN function class does not explain why in practice we can easily find them. In this section, we introduce standard methods, namely stochastic gradient descent (SGD) and its variants, to train deep neural networks (or to find such a good approximator). As with many statistical machine learning tasks, training DNNs follows the empirical risk minimization (ERM) paradigm which solves the following optimization problem
$${\text{minimize}}_{\bm{\theta}\in {\mathbb{R}}^{p}}\mathit{\hspace{1em}\hspace{1em}}{\mathrm{\ell}}_{n}\left(\bm{\theta}\right)\triangleq \frac{1}{n}\sum _{i=1}^{n}\mathcal{L}(f({\bm{x}}_{i};\bm{\theta}),{y}_{i}).$$  (26) 
Here $\mathcal{L}(f({\bm{x}}_{i};\bm{\theta}),{y}_{i})$ measures the discrepancy between the prediction $f({\bm{x}}_{i};\bm{\theta})$ of the neural network and the true label ${y}_{i}$. Correspondingly, denote by $\mathrm{\ell}(\bm{\theta})\triangleq {\mathbb{E}}_{(\bm{x},y)\sim \mathcal{D}}[\mathcal{L}(f(\bm{x};\bm{\theta}),\bm{y})]$ the outofsample error, where $\mathcal{D}$ is the joint distribution over $(y,\bm{x})$. Solving ERM (26) for deep neural nets faces various challenges that roughly fall into the following three categories.

•
Scalability and nonconvexity. Both the sample size $n$ and the number of parameters $p$ can be huge for modern deep learning applications, as we have seen in Table 1. Many optimization algorithms are not practical due to the computational costs and memory constraints. What is worse, the empirical loss function ${\mathrm{\ell}}_{n}(\bm{\theta})$ in deep learning is often nonconvex. It is a priori not clear whether an optimization algorithm can drive the empirical loss (26) small.

•
Numerical stability. With a large number of layers in DNNs, the magnitudes of the hidden nodes can be drastically different, which may result in the “exploding gradients” or “vanishing gradients” issue during the training process. This is because the recursive relations across layers often lead to exponentially increasing$$/$$decreasing values in both forward passes and backward passes.

•
Generalization performance. Our ultimate goal is to find a parameter $\widehat{\bm{\theta}}$ such that the outofsample error $\mathrm{\ell}(\widehat{\bm{\theta}})$ is small. However, in the overparametrized regime where $p$ is much larger than $n$, the underlying neural network has the potential to fit the training data perfectly while performing poorly on the test data. To avoid this overfitting issue, proper regularization, whether explicit or implicit, is needed in the training process for the neural nets to generalize.
In the following three subsections, we discuss practical solutions$$/$$proposals to address these challenges.
6.1 Stochastic gradient descent
Stochastic gradient descent (SGD) [100] is by far the most popular optimization algorithm to solve ERM (26) for largescale problems. It has the following simple update rule:
$${\bm{\theta}}^{t+1}={\bm{\theta}}^{t}{\eta}_{t}G({\bm{\theta}}^{t})\mathit{\hspace{1em}\hspace{1em}}\text{with}\mathit{\hspace{1em}\hspace{1em}}G\left({\bm{\theta}}^{t}\right)=\nabla \mathcal{L}(f({\bm{x}}_{{i}_{t}};{\bm{\theta}}^{t}),{y}_{{i}_{t}})$$  (27) 
for $t=0,1,2,\mathrm{\dots}$, where ${\eta}_{t}>0$ is the step size (or learning rate), ${\bm{\theta}}^{0}\in {\mathbb{R}}^{p}$ is an initial point and ${i}_{t}$ is chosen randomly from $\{1,2,\mathrm{\cdots},n\}$. It is easy to verify that $G({\bm{\theta}}^{t})$ is an unbiased estimate of $\nabla {\mathrm{\ell}}_{n}({\bm{\theta}}^{t})$. The advantage of SGD is clear: compared with gradient descent, which goes over the entire dataset in every update, SGD uses a single example in each update and hence is considerably more efficient in terms of both computation and memory (especially in the first few iterations).
Apart from practical benefits of SGD, how well does SGD perform theoretically in terms of minimizing ${\mathrm{\ell}}_{n}(\bm{\theta})$? We begin with the convex case, i.e., the case where the loss function is convex w.r.t. $\bm{\theta}$. It is well understood in literature that with proper choices of the step sizes $\{{\eta}_{t}\}$, SGD is guaranteed to achieve both consistency and asymptotic normality.

•
Consistency. If $\mathrm{\ell}(\bm{\theta})$ is a strongly convex function^{7}^{7} 7 For results on consistency and asymptotic normality, we consider the case where in each step of SGD, the stochastic gradient is computed using a fresh sample $(y,\bm{x})$ from $\mathcal{D}$. This allows to view SGD as an optimization algorithm to minimize the population loss $\mathrm{\ell}(\bm{\theta})$., then under some mild conditions^{8}^{8} 8 One example of such condition can be constraining the second moment of the gradients: $\mathbb{E}\left[{\parallel \nabla \mathcal{L}({\bm{x}}_{i},{y}_{i};{\bm{\theta}}^{t})\parallel}_{2}^{2}\right]\le {C}_{1}+{C}_{2}{\parallel {\bm{\theta}}^{t}{\bm{\theta}}^{*}\parallel}_{2}^{2}$ for some ${C}_{1},{C}_{2}>0$. See [16] for details., learning rates that satisfy
$$ (28) guarantee almost sure convergence to the unique minimizer ${\bm{\theta}}^{*}\triangleq {\text{argmin}}_{\bm{\theta}}\mathrm{\ell}(\bm{\theta})$, i.e., ${\bm{\theta}}^{t}\stackrel{\text{a.s.}}{\to}{\bm{\theta}}^{*}$ as $t\to \mathrm{\infty}$ [100, 64, 16, 69]. The requirements in (28) can be viewed from the perspective of biasvariance tradeoff: the first condition ensures that the iterates can reach the minimizer (controlled bias), and the second ensures that stochasticity does not prevent convergence (controlled variance).

•
Asymptotic normality. It is proved by [97] that for robust linear regression with fixed dimension $p$, under the choice ${\eta}_{t}={t}^{1}$, $\sqrt{t}({\bm{\theta}}^{t}{\bm{\theta}}^{*})$ is asymptotically normal under some regularity conditions (but ${\bm{\theta}}^{t}$ is not asymptotically efficient in general). Moreover, by averaging the iterates of SGD, [96] proved that even with a larger step size ${\eta}_{t}\propto {t}^{\alpha},\alpha \in (1/2,1)$, the averaged iterate ${\overline{\bm{\theta}}}^{t}={t}^{1}{\sum}_{s=1}^{t}{\bm{\theta}}^{s}$ is asymptotic efficient for robust linear regression. These strong results show that SGD with averaging performs as well as the MLE asymptotically, in addition to its computational efficiency.
These classical results, however, fail to explain the effectiveness of SGD when dealing with nonconvex loss functions in deep learning. Admittedly, finding global minima of nonconvex functions is computationally infeasible in the worst case. Nevertheless, recent work [4, 32] bypasses the worst case scenario by focusing on losses incurred by overparametrized deep learning models. In particular, they show that (stochastic) gradient descent converges linearly towards the global minimizer of ${\mathrm{\ell}}_{n}(\bm{\theta})$ as long as the neural network is sufficiently overparametrized. This phenomenon is formalized below.
Theorem 5 (Theorem 2 in [4]).
Let ${\mathrm{\{}\mathrm{(}{y}_{i}\mathrm{,}{\mathbf{x}}_{i}\mathrm{)}\mathrm{\}}}_{\mathrm{1}\mathrm{\le}i\mathrm{\le}n}$ be a training set satisfying ${\mathrm{min}}_{i\mathrm{,}j\mathrm{:}i\mathrm{\ne}j}\mathit{}{\mathrm{\parallel}{\mathbf{x}}_{i}\mathrm{}{\mathbf{x}}_{j}\mathrm{\parallel}}_{\mathrm{2}}\mathrm{\ge}\delta \mathrm{>}\mathrm{0}$. Consider fitting the data using a feedforward neural network (1) with ReLU activations. Denote by $L$ (resp. $W$) the depth (resp. width) of the network. Suppose that the neural network is sufficiently overparametrized, i.e.,
$$W\gg \mathrm{\U0001d5c9\U0001d5c8\U0001d5c5\U0001d5d2}(n,L,\frac{1}{\delta}),$$  (29) 
where $\mathrm{poly}$ means a polynomial function. Then with high probability, running SGD (27) with certain random initialization and properly chosen step sizes yields ${\mathrm{\ell}}_{n}\mathit{}\mathrm{(}{\mathbf{\theta}}^{t}\mathrm{)}\mathrm{\le}\epsilon $ in $t\mathrm{\asymp}\mathrm{log}\mathit{}\frac{\mathrm{1}}{\epsilon}$ iterations.
Two notable features are worth mentioning: (1) first, the network under consideration is sufficiently overparametrized (cf. (29)) in which the number of parameters is much larger than the number of samples, and (2) one needs to initialize the weight matrices to be in nearisometry such that the magnitudes of the hidden nodes do not blow up or vanish. In a nutshell, overparametrization and random initialization together ensure that the loss function (26) has a benign landscape^{9}^{9} 9 In [4], the loss function ${\mathrm{\ell}}_{n}(\bm{\theta})$ satisfies the PL condition. around the initial point, which in turn implies fast convergence of SGD iterates.
There are certainly other challenges for vanilla SGD to train deep neural nets: (1) training algorithms are often implemented in GPUs, and therefore it is important to tailor the algorithm to the infrastructure, (2) the vanilla SGD might converge very slowly for deep neural networks, albeit good theoretical guarantees for wellbehaved problems, and (3) the learning rates $\{{\eta}_{t}\}$ can be difficult to tune in practice. To address the aforementioned challenges, three important variants of SGD, namely minibatch SGD, momentumbased SGD, and SGD with adaptive learning rates are introduced.
6.1.1 Minibatch SGD
Modern computational infrastructures (e.g., GPUs) can evaluate the gradient on a number (say 64) of examples as efficiently as evaluating that on a single example. To utilize this advantage, minibatch SGD with batch size $K\ge 1$ forms the stochastic gradient through $K$ random samples:
$${\bm{\theta}}^{t+1}={\bm{\theta}}^{t}{\eta}_{t}G({\bm{\theta}}^{t})\mathit{\hspace{1em}\hspace{1em}}\text{with}\mathit{\hspace{1em}\hspace{1em}}G({\bm{\theta}}^{t})=\frac{1}{K}\sum _{k=1}^{K}\nabla \mathcal{L}(f({\bm{x}}_{{i}_{t}^{k}};{\bm{\theta}}^{t}),{y}_{{i}_{t}^{k}}),$$  (30) 
where for each $1\le k\le K$, ${i}_{t}^{k}$ is sampled uniformly from $\{1,2,\mathrm{\cdots},n\}$. Minibatch SGD, which is an “interpolation” between gradient descent and stochastic gradient descent, achieves the best of both worlds: (1) using $1\ll K\ll n$ samples to estimate the gradient, one effectively reduces the variance and hence accelerates the convergence, and (2) by taking the batch size $K$ appropriately (say 64 or 128), the stochastic gradient $G({\bm{\theta}}^{t})$ can be efficiently computed using the matrix computation toolboxes on GPUs.
6.1.2 Momentumbased SGD
While minibatch SGD forms the foundation of training neural networks, it can sometimes be slow to converge due to its oscillation behavior [121]. Optimization community has long investigated how to accelerate the convergence of gradient descent, which results in a beautiful technique called momentum methods [95, 88]. Similar to gradient descent with moment, momentumbased SGD, instead of moving the iterate ${\bm{\theta}}^{t}$ in the direction of the current stochastic gradient $G({\bm{\theta}}^{t})$, smooth the past (stochastic) gradients $\{G({\bm{\theta}}^{t})\}$ to stabilize the update directions. Mathematically, let ${\bm{v}}^{t}\in {\mathbb{R}}^{p}$ be the direction of update in the $t$th iteration, i.e.,
$${\bm{\theta}}^{t+1}={\bm{\theta}}^{t}{\eta}_{t}{\bm{v}}^{t}.$$ 
Here ${\bm{v}}^{0}=G({\bm{\theta}}^{0})$ and for $t=1,2,\mathrm{\cdots}$
$${\bm{v}}^{t}=\rho {\bm{v}}^{t1}+G({\bm{\theta}}^{t})$$  (31) 
with $$. A typical choice of $\rho $ is 0.9. Notice that $\rho =0$ recovers the minibatch SGD (30), where no past information of gradients is used. A simple unrolling of ${\bm{v}}^{t}$ reveals that ${\bm{v}}^{t}$ is actually an exponential averaging of the past gradients, i.e., ${\bm{v}}^{t}={\sum}_{j=0}^{t}{\rho}^{tj}G({\bm{\theta}}^{j}).$ Compared with vanilla minibatch SGD, the inclusion of the momentum “smoothes” the oscillation direction and accumulates the persistent descent direction. We want to emphasize that theoretical justifications of momentum in the stochastic setting is not fully understood [63, 60].
6.1.3 SGD with adaptive learning rates
In optimization, preconditioning is often used to accelerate firstorder optimization algorithms. In principle, one can apply this to SGD, which yields the following update rule:
$${\bm{\theta}}^{t+1}={\bm{\theta}}^{t}{\eta}_{t}{\bm{P}}_{t}^{1}G({\bm{\theta}}^{t})$$  (32) 
with ${\bm{P}}_{t}\in {\mathbb{R}}^{p\times p}$ being a preconditioner at the $t$th step. Newton’s method can be viewed as one type of preconditioning where ${\bm{P}}_{t}={\nabla}^{2}\mathrm{\ell}({\bm{\theta}}^{t})$. The advantages of preconditioning are twofold: first, a good preconditioner reduces the condition number by changing the local geometry to be more homogeneous, which is amenable to fast convergence; second, a good preconditioner frees practitioners from laboring tuning of the step sizes, as is the case with Newton’s method. AdaGrad, an adaptive gradient method proposed by [33], builds a preconditioner ${\bm{P}}_{t}$ based on information of the past gradients:
$${\bm{P}}_{t}={\left\{\mathrm{\U0001d5bd\U0001d5c2\U0001d5ba\U0001d5c0}\left(\sum _{j=0}^{t}G\left({\bm{\theta}}^{t}\right)G{\left({\bm{\theta}}^{t}\right)}^{\top}\right)\right\}}^{1/2}.$$  (33) 
Since we only require the diagonal part, this preconditioner (and its inverse) can be efficiently computed in practice. In addition, investigating (32) and (33), one can see that AdaGrad adapts to the importance of each coordinate of the parameters by setting smaller learning rates for frequent features, whereas larger learning rates for those infrequent ones. In practice, one adds a small quantity $\delta >0$ (say ${10}^{8}$) to the diagonal entries to avoid singularity (numerical underflow). A notable drawback of AdaGrad is that the effective learning rate vanishes quickly along the learning process. This is because the historical sum of the gradients can only increase with time. RMSProp [52] is a popular remedy for this problem which incorporates the idea of exponential averaging:
$${\bm{P}}_{t}={\left\{\mathrm{\U0001d5bd\U0001d5c2\U0001d5ba\U0001d5c0}\left(\rho {\bm{P}}_{t1}+(1\rho )G\left({\bm{\theta}}^{t}\right)G{\left({\bm{\theta}}^{t}\right)}^{\top}\right)\right\}}^{1/2}.$$  (34) 
Again, the decaying parameter $\rho $ is usually set to be $0.9$. Later, Adam [65, 99] combines the momentum method and adaptive learning rate and becomes the default training algorithms in many deep learning applications.
6.2 Easing numerical instability
For very deep neural networks or RNNs with long dependencies, training difficulties often arise when the values of nodes have different magnitudes or when the gradients “vanish” or “explode” during backpropagation. Here we discuss three partial solutions to alleviate this problem.
6.2.1 ReLU activation function
One useful characteristic of the ReLU function is that its derivative is either $0$ or $1$, and the derivative remains $1$ even for a large input. This is in sharp contrast with the standard sigmoid function ${(1+{e}^{t})}^{1}$ which results in a very small derivative when inputs have large magnitude. The consequence of small derivatives across many layers is that gradients tend to be “killed”, which means that gradients become approximately zero in deep nets.
6.2.2 Skip connections
We have introduced skip connections in Section 3.3. Why are skip connections helpful for reducing numerical instability? This structure does not introduce a larger function space, since the identity map can be also represented with ReLU activations: $\bm{x}=\bm{\sigma}(\bm{x})\bm{\sigma}(\bm{x})$.
One explanation is that skip connections bring ease to the training$$/$$optimization process. Suppose that we have a general nonlinear function $\mathbf{F}({\bm{x}}_{\mathrm{\ell}};{\bm{\theta}}_{\mathrm{\ell}})$. With a skip connection, we represent the map as ${\bm{x}}_{\mathrm{\ell}+1}={\bm{x}}_{\mathrm{\ell}}+\mathbf{F}({\bm{x}}_{\mathrm{\ell}};{\bm{\theta}}_{\mathrm{\ell}})$ instead. Now the gradient $\partial {\bm{x}}_{\mathrm{\ell}+1}/\partial {\bm{x}}_{\mathrm{\ell}}$ becomes
$$\frac{\partial {\bm{x}}_{\mathrm{\ell}+1}}{\partial {\bm{x}}_{\mathrm{\ell}}}=\mathbf{I}+\frac{\partial \mathbf{F}({\bm{x}}_{\mathrm{\ell}};{\bm{\theta}}_{\mathrm{\ell}})}{\partial {\bm{x}}_{\mathrm{\ell}}}\mathit{\hspace{1em}\hspace{1em}}\text{instead of}\mathit{\hspace{1em}\hspace{1em}}\frac{\partial \mathbf{F}({\bm{x}}_{\mathrm{\ell}};{\bm{\theta}}_{\mathrm{\ell}})}{\partial {\bm{x}}_{\mathrm{\ell}}},$$  (35) 
where $\mathbf{I}$ is an identity matrix. By the chain rule, gradient update requires computing products of many components, e.g., $\frac{\partial {\bm{x}}_{L}}{\partial {\bm{x}}_{1}}={\prod}_{\mathrm{\ell}=1}^{L1}\frac{\partial {\bm{x}}_{\mathrm{\ell}+1}}{\partial {\bm{x}}_{\mathrm{\ell}}}$, so it is desirable to keep the spectra (singular values) of each component $\frac{\partial {\bm{x}}_{\mathrm{\ell}+1}}{\partial {\bm{x}}_{\mathrm{\ell}}}$ close to $1$. In neural nets, with skip connections, this is easily achieved if the parameters have small values; otherwise, this may not be achievable even with careful initialization and tuning. Notably, training neural nets with hundreds of layers is possible with the help of skip connections.
6.2.3 Batch normalization
Recall that in regression analysis, one often standardizes the design matrix so that the features have zero mean and unit variance. Batch normalization extends this standardization procedure from the input layer to all the hidden layers. Mathematically, fix a minibatch of input data ${\{({\bm{x}}_{i},{y}_{i})\}}_{i\in \mathcal{B}}$, where $\mathcal{B}\subset [n]$. Let ${\bm{h}}_{i}^{(\mathrm{\ell})}$ be the feature of the $i$th example in the $\mathrm{\ell}$th layer ($\mathrm{\ell}=0$ corresponds to the input ${\bm{x}}_{i}$). The batch normalization layer computes the normalized version of ${\bm{h}}_{i}^{(\mathrm{\ell})}$ via the following steps:
$\bm{\mu}$  $\triangleq {\displaystyle \frac{1}{\left\mathcal{B}\right}}{\displaystyle \sum _{i\in \mathcal{B}}}{\bm{h}}_{i}^{(\mathrm{\ell})},{\bm{\sigma}}^{2}\triangleq {\displaystyle \frac{1}{\left\mathcal{B}\right}}{\displaystyle \sum _{i\in \mathcal{B}}}{\left({\bm{h}}_{i}^{(\mathrm{\ell})}\bm{\mu}\right)}^{2}\mathit{\hspace{1em}\hspace{1em}}\text{and}\mathit{\hspace{1em}\hspace{1em}}{\bm{h}}_{i,\text{norm}}^{(l)}\triangleq {\displaystyle \frac{{\bm{h}}_{i}^{(\mathrm{\ell})}\bm{\mu}}{\bm{\sigma}}}.$ 
Here all the operations are elementwise. In words, batch normalization computes the zscore for each feature over the minibatch $\mathcal{B}$ and use that as inputs to subsequent layers. To make it more versatile, a typical batch normalization layer has two additional learnable parameters ${\bm{\gamma}}^{(\mathrm{\ell})}$ and ${\bm{\beta}}^{(\mathrm{\ell})}$ such that
$${\bm{h}}_{i,\text{new}}^{(l)}={\bm{\gamma}}^{(l)}\odot {\bm{h}}_{i,\text{norm}}^{(l)}+{\bm{\beta}}^{(l)}.$$ 
Again $\odot $ denotes the elementwise multiplication. As can be seen, ${\bm{\gamma}}^{(\mathrm{\ell})}$ and ${\bm{\beta}}^{(\mathrm{\ell})}$ set the new feature ${\bm{h}}_{i\text{new}}^{(l)}$ to have mean ${\bm{\beta}}^{(\mathrm{\ell})}$ and standard deviation ${\bm{\gamma}}^{(\mathrm{\ell})}$. The introduction of batch normalization makes the training of neural networks much easier and smoother. More importantly, it allows the neural nets to perform well over a large family of hyperparameters including the number of layers, the number of hidden units, etc. At test time, the batch normalization layer needs more care. For brevity we omit the details and refer to [58].
6.3 Regularization techniques
So far we have focused on training techniques to drive the empirical loss (26) small efficiently. Here we proceed to discuss common practice to improve the generalization power of trained neural nets.
6.3.1 Weight decay
One natural regularization idea is to add an ${\mathrm{\ell}}_{2}$ penalty to the loss function. This regularization technique is known as the weight decay in deep learning. We have seen one example in (9). For general deep neural nets, the loss to optimize is ${\mathrm{\ell}}_{n}^{\lambda}(\bm{\theta})={\mathrm{\ell}}_{n}(\bm{\theta})+{r}_{\lambda}(\bm{\theta})$ where
$${r}_{\lambda}(\bm{\theta})=\lambda \sum _{\mathrm{\ell}=1}^{L}\sum _{j,{j}^{\prime}}{\left[{W}_{j,{j}^{\prime}}^{(\mathrm{\ell})}\right]}^{2}.$$ 
Note that the bias (intercept) terms are not penalized. If ${\mathrm{\ell}}_{n}(\bm{\theta})$ is a least square loss, then regularization with weight decay gives precisely ridge regression. The penalty ${r}_{\lambda}(\bm{\theta})$ is a smooth function and thus it can be also implemented efficiently with backpropagation.
6.3.2 Dropout
Dropout, introduced by [53], prevents overfitting by randomly dropping out subsets of features during training. Take the $l$th layer of the feedforward neural network as an example. Instead of propagating all the features in ${\bm{h}}^{(\mathrm{\ell})}$ for later computations, dropout randomly omits some of its entries by
$${\bm{h}}_{\text{drop}}^{(\mathrm{\ell})}={\bm{h}}^{(\mathrm{\ell})}\odot {\mathrm{\U0001d5c6\U0001d5ba\U0001d5cc\U0001d5c4}}^{\mathrm{\ell}},$$ 
where $\odot $ denotes elementwise multiplication as before, and ${\mathrm{\U0001d5c6\U0001d5ba\U0001d5cc\U0001d5c4}}^{\mathrm{\ell}}$ is a vector of Bernoulli variables with success probability $p$. It is sometimes useful to rescale the features ${\bm{h}}_{\text{inv drop}}^{(\mathrm{\ell})}={\bm{h}}_{\text{drop}}^{(\mathrm{\ell})}/p$, which is called inverted dropout. During training, ${\mathrm{\U0001d5c6\U0001d5ba\U0001d5cc\U0001d5c4}}^{\mathrm{\ell}}$ are i.i.d. vectors across minibatches and layers. However, when testing on fresh samples, dropout is disabled and the original features ${\bm{h}}^{(\mathrm{\ell})}$ are used to compute the output label $y$. It has been nicely shown by [128] that for generalized linear models, dropout serves as adaptive regularization. In the simplest case of linear regression, it is equivalent to ${\mathrm{\ell}}_{2}$ regularization. Another possible way to understand the regularization effect of dropout is through the lens of bagging [45]. Since different minibatches has different masks, dropout can be viewed as training a large ensemble of classifiers at the same time, with a further constraint that the parameters are shared. Theoretical justification remains elusive.
6.3.3 Data augmentation
Data augmentation is a technique of enlarging the dataset when we have knowledge about invariance structure of data. It implicitly increases the sample size and usually regularizes the model effectively. For example, in image classification, we have strong prior knowledge about what invariance properties a good classifier should possess. The label of an image should not be affected by translation, rotation, flipping, and even crops of the image. Hence one can augment the dataset by randomly translating, rotating and cropping the images in the original dataset.
Formally, during training we want to minimize the loss ${\mathrm{\ell}}_{n}(\bm{\theta})={\sum}_{i}\mathcal{L}(f({\bm{x}}_{i};\bm{\theta}),{y}_{i})$ w.r.t. parameters $\bm{\theta}$, and we know a priori that certain transformation $T\in \mathcal{T}$ where $T:{\mathbb{R}}^{d}\to {\mathbb{R}}^{d}$ (e.g., affine transformation) should not change the category$$/$$label of a training sample. In principle, if computation costs were not a consideration, we could convert this knowledge to a constraint ${f}_{\bm{\theta}}(T{\bm{x}}_{i})={f}_{\bm{\theta}}({\bm{x}}_{i}),\forall T\in \mathcal{T}$ in the minimization formulation. Instead of solving a constrained optimization problem, data augmentation enlarges the training dataset by sampling $T\in \mathcal{T}$ and generating new data $\{(T{\bm{x}}_{i},{y}_{i})\}$. In this sense, data augmentation induces invariance properties through sampling, which results in a much bigger dataset than the original one.
7 Generalization power
Section 6 has focused on the insample$$/$$ training error obtained via SGD, but this alone does not guarantee good performance with respect to the outofsample$$/$$test error. The gap between the insample error and the outofsample error, namely the generalization gap, has been the focus of statistical learning theory since its birth; see [111] for an excellent introduction to this topic.
While understanding the generalization power of deep neural nets is difficult, we sample recent endeavors in this section. From a high level point of view, these approaches can be divided into two categories, namely algorithmindependent controls and algorithmdependent controls. More specifically, algorithmindependent controls focus solely on bounding the complexity of the function class represented by certain deep neural networks. In contrast, algorithmdependent controls take into account the algorithm (e.g., SGD) used to train the neural network.
7.1 Algorithmindependent controls: uniform convergence
The key to algorithmindependent controls is the notion of complexity of the function class parametrized by certain neural networks. Informally, as long as the complexity is not too large, the generalization gap of any function in the function class is wellcontrolled. However, the standard complexity measure (e.g., VC dimension [126]) is at least proportional to the number of weights in a neural network [5, 111], which fails to explain the practical success of deep learning. The caveat here is that the function class under consideration is all the functions realized by certain neural networks, with no restrictions on the size of the weights at all. On the other hand, for the class of linear functions with bounded norm, i.e., $\{\bm{x}\mapsto {\bm{w}}^{\top}\bm{x}{\parallel \bm{w}\parallel}_{2}\le M\}$, it is well understood that the complexity of this function class (measured in terms of the empirical Rademacher complexity) with respect to a random sample ${\{{\bm{x}}_{i}\}}_{1\le i\le n}$ is upper bounded by ${\mathrm{max}}_{i}{\parallel {\bm{x}}_{i}\parallel}_{2}M/\sqrt{n}$, which is independent of the number of parameters in $\bm{w}$. This motivates researchers to investigate the complexity of normcontrolled deep neural networks^{10}^{10} 10 Such attempts have been made in the seminal work [13]. [89, 14, 43, 74]. Setting the stage, we introduce a few necessary notations and facts. The key object under study is the function class parametrized by the following fullyconnected neural network with depth $L$:
$${\mathcal{F}}_{L}\triangleq \left\{\bm{x}\mapsto {\bm{W}}_{L}\bm{\sigma}\left({\bm{W}}_{L1}\bm{\sigma}\left(\mathrm{\cdots}{\bm{W}}_{2}\bm{\sigma}\left({\bm{W}}_{1}\bm{x}\right)\right)\right)\right({\bm{W}}_{1},\mathrm{\cdots},{\bm{W}}_{L})\in \mathcal{W}\}.$$  (36) 
Here $({\bm{W}}_{1},{\bm{W}}_{2},\mathrm{\cdots},{\bm{W}}_{L})\in \mathcal{W}$ represents a certain constraint on the parameters. For instance, one can restrict the Frobenius norm of each parameter ${\bm{W}}_{l}$ through the constraint ${\parallel {\bm{W}}_{l}\parallel}_{\mathrm{F}}\le {M}_{\mathrm{F}}(l)$, where ${M}_{\mathrm{F}}(l)$ is some positive quantity. With regard to the complexity measure, it is standard to use Rademacher complexity to control the capacity of the function class of interest.
Definition 1 (Empirical Rademacher complexity).
The empirical Rademacher complexity of a function class $\mathrm{F}$ w.r.t. a dataset $S\mathrm{\triangleq}{\mathrm{\{}{\mathbf{x}}_{i}\mathrm{\}}}_{\mathrm{1}\mathrm{\le}i\mathrm{\le}n}$ is defined as
$${\mathcal{R}}_{S}\left(\mathcal{F}\right)={\mathbb{E}}_{\bm{\epsilon}}\left[\underset{f\in \mathcal{F}}{sup}\frac{1}{n}\sum _{i=1}^{n}{\epsilon}_{i}f\left({\bm{x}}_{i}\right)\right],$$  (37) 
where $\mathbf{\epsilon}\mathrm{\triangleq}\mathrm{(}{\epsilon}_{\mathrm{1}}\mathrm{,}{\epsilon}_{\mathrm{2}}\mathrm{,}\mathrm{\cdots}\mathrm{,}{\epsilon}_{n}\mathrm{)}$ is composed of i.i.d. Rademacher random variables, i.e., $\mathrm{P}\mathrm{(}{\epsilon}_{i}\mathrm{=}\mathrm{1}\mathrm{)}\mathrm{=}\mathrm{P}\mathrm{(}{\epsilon}_{i}\mathrm{=}\mathrm{}\mathrm{1}\mathrm{)}\mathrm{=}\mathrm{1}\mathrm{/}\mathrm{2}$.
In words, Rademacher complexity measures the ability of the function class to fit the random noise represented by $\bm{\epsilon}$. Intuitively, a function class with a larger Rademacher complexity is more prone to overfitting. We now formalize the connection between the empirical Rademacher complexity and the outofsample error; see Chapter 24 in [111].
Theorem 6.
Assume that for all $f\mathrm{\in}\mathrm{F}$ and all $\mathrm{(}y\mathrm{,}\mathbf{x}\mathrm{)}$ we have $\mathrm{}\mathrm{L}\mathit{}\mathrm{(}f\mathit{}\mathrm{(}\mathbf{x}\mathrm{)}\mathrm{,}y\mathrm{)}\mathrm{}\mathrm{\le}\mathrm{1}$. In addition, assume that for any fixed $y$, the univariate function $\mathrm{L}\mathit{}\mathrm{(}\mathrm{\cdot}\mathrm{,}y\mathrm{)}$ is Lipschitz with constant 1. Then with probability at least $\mathrm{1}\mathrm{}\delta $ over the sample $S\mathrm{\triangleq}{\mathrm{\{}\mathrm{(}{y}_{i}\mathrm{,}{\mathbf{x}}_{i}\mathrm{)}\mathrm{\}}}_{\mathrm{1}\mathrm{\le}i\mathrm{\le}n}\mathit{}\stackrel{\mathrm{i}\mathrm{.}\mathrm{i}\mathrm{.}\mathrm{d}\mathrm{.}}{\mathrm{\sim}}\mathit{}\mathrm{D}$, one has for all $f\mathrm{\in}\mathrm{F}$
$$\underset{\mathrm{out}\mathit{\text{}}\mathrm{of}\mathit{\text{}}\mathrm{sample}\mathrm{error}}{\underset{\u23df}{{\mathbb{E}}_{(y,\bm{x})\sim \mathcal{D}}\left[\mathcal{L}(f(\bm{x}),y)\right]}}\le \underset{\mathrm{in}\mathit{\text{}}\mathrm{sample}\mathrm{error}}{\underset{\u23df}{\frac{1}{n}\sum _{i=1}^{n}\mathcal{L}(f({\bm{x}}_{i}),{y}_{i})}}+2{\mathcal{R}}_{S}\left(\mathcal{F}\right)+4\sqrt{\frac{\mathrm{log}\left(4/\delta \right)}{n}}.$$ 
In English, the generalization gap of any function $f$ that lies in $\mathcal{F}$ is wellcontrolled as long as the Rademacher complexity of $\mathcal{F}$ is not too large. With this connection in place, we single out the following complexity bound.
Theorem 7 (Theorem 1 in [43]).
Consider the function class ${\mathrm{F}}_{L}$ in (36), where each parameter ${\mathbf{W}}_{l}$ has Frobenius norm at most ${M}_{\mathrm{F}}\mathit{}\mathrm{(}l\mathrm{)}$. Further suppose that the elementwise activation function $\sigma \mathit{}\mathrm{(}\mathrm{\cdot}\mathrm{)}$ is $\mathrm{1}$Lipschitz and positivehomogeneous (i.e., $\sigma \mathit{}\mathrm{(}c\mathrm{\cdot}x\mathrm{)}\mathrm{=}c\mathit{}\sigma \mathit{}\mathrm{(}x\mathrm{)}$ for all $c\mathrm{\ge}\mathrm{0}$). Then the empirical Rademacher complexity (37) w.r.t. $S\mathrm{\triangleq}{\mathrm{\{}{\mathbf{x}}_{i}\mathrm{\}}}_{\mathrm{1}\mathrm{\le}i\mathrm{\le}n}$ satisfies
$${\mathcal{R}}_{S}\left({\mathcal{F}}_{L}\right)\le \underset{i}{\mathrm{max}}{\parallel {\bm{x}}_{i}\parallel}_{2}\cdot \frac{4\sqrt{L}{\prod}_{l=1}^{L}{M}_{\mathrm{F}}(l)}{\sqrt{n}}.$$  (38) 
The upper bound of the empirical Rademacher complexity (38) is in a similar vein to that of linear functions with bounded norm, i.e., ${\mathrm{max}}_{i}{\parallel {\bm{x}}_{i}\parallel}_{2}M/\sqrt{n}$, where $\sqrt{L}{\prod}_{l=1}^{L}{M}_{\mathrm{F}}(l)$ plays the role of $M$ in the latter case. Moreover, ignoring the term $\sqrt{L}$, the upper bound (38) does not depend on the size of the network in an explicit way if ${M}_{F}(l)$ sharply concentrates around $1$. This reveals that the capacity of the neural network is wellcontrolled, regardless of the number of parameters, as long as the Frobenius norm of the parameters is bounded. Extensions to other norm constraints, e.g., spectral norm constraints, path norm constraints have been considered by [89, 14, 74, 67, 34]. This line of work improves upon traditional capacity analysis of neural networks in the overparametrized setting, because the upper bounds derived are often sizeindependent. Having said this, two important remarks are in order: (1) the upper bounds (e.g., ${\prod}_{l=1}^{L}{M}_{\mathrm{F}}(l)$) involve implicit dependence on the size of the weight matrix and the depth of the neural network, which is hard to characterize; (2) the upper bound on the Rademacher complexity offers a uniform bound over all functions in the function class, which is a pure statistical result. However, it stays silent about how and why standard training algorithms like SGD can obtain a function whose parameters have small norms.
7.2 Algorithmdependent controls
In this subsection, we bring computational thinking into statistics and investigate the role of algorithms in the generalization power of deep learning. The consideration of algorithms is quite natural and well motivated: (1) local/global minima reached by different algorithms can exhibit totally different generalization behaviors due to extreme nonconvexity, which marks a huge difference from traditional models, (2) the effective capacity of neural nets is possibly not large, since a particular algorithm does not explore the entire parameter space.
These demonstrate the fact that on top of the complexity of the function class, the inherent property of the algorithm we use plays an important role in the generalization ability of deep learning. In what follows, we survey three different ways to obtain upper bounds on the generalization errors by exploiting properties of the algorithms.
7.2.1 Mean field view of neural nets
As we have emphasized, modern deep learning models are highly overparametrized. A line of work [83, 116, 104, 25, 82, 61] approximates the ensemble of weights by an asymptotic limit as the number of hidden units tends to infinity, so that the dynamics of SGD can be studied via certain partial different equations.
More specifically, let $\widehat{f}(\bm{x};\bm{\theta})={N}^{1}{\sum}_{i=1}^{N}\sigma ({\bm{\theta}}_{i}^{\top}\bm{x})$ be a function given by a onehiddenlayer neural net with $N$ hidden units, where $\sigma (\cdot )$ is the ReLU activation function and parameters $\bm{\theta}\triangleq {[{\bm{\theta}}_{1},\mathrm{\dots},{\bm{\theta}}_{N}]}^{\top}\in {\mathbb{R}}^{N\times d}$ are suitably randomly initialized. Consider the regression setting where we want to minimize the population risk ${R}_{N}(\bm{\theta})=\mathbb{E}[{(y\widehat{f}(\bm{x};\bm{\theta}))}^{2}]$ over parameters $\bm{\theta}$. A key observation is that this population risk depends on the parameters $\bm{\theta}$ only through its empirical distribution, i.e., ${\widehat{\rho}}^{(N)}={N}^{1}{\sum}_{i=1}^{N}{\delta}_{{\bm{\theta}}_{i}}$ where ${\delta}_{{\bm{\theta}}_{i}}$ is a point mass at ${\bm{\theta}}_{i}$. This motivates us to view express ${R}_{N}(\bm{\theta})$ equivalently as $R({\widehat{\rho}}^{(N)})$, where $R(\cdot )$ is a functional that maps distributions to real numbers. Running SGD for ${R}_{N}(\cdot )$—in a suitable scaling limit—results in a gradient flow on the space of distributions endowed with the Wasserstein metric that minimizes $R(\cdot )$. It turns out that the empirical distribution ${\widehat{\rho}}_{k}^{(N)}$ of the parameters after $k$ steps of SGD is well approximated by the gradient follow, as long as the the neural net is overparametrized (i.e., $N\gg d$) and the number of steps is not too large. In particular, [83] have shown that under certain regularity conditions,
$$\underset{k\in [0,T/\epsilon ]\cap \mathbb{N}}{sup}\leftR({\widehat{\rho}}^{(N)})R\left({\rho}_{k\epsilon}\right)\right\lesssim {e}^{T}\sqrt{\frac{1}{N}\vee \epsilon}\cdot \sqrt{d+\mathrm{log}\frac{N}{\epsilon}},$$ 
where $\epsilon >0$ is an proxy for the step size of SGD and ${\rho}_{k\epsilon}$ is the distribution of the gradient flow at time $k\epsilon $. In words, the outofsample error under ${\bm{\theta}}^{k}$ generated by SGD is wellapproximated by that of ${\rho}_{k\epsilon}$. Viewing the optimization problem from the distributional aspect greatly simplifies the problem conceptually, as the complicated optimization problem is now passed into its limit version—for this reason, this analytical approach is called the mean field perspective. In particular, [83] further demonstrated that in some simple settings, the outofsample error $R({\rho}_{k\epsilon})$ of the distributional limit can be fully characterized. Nevertheless, how well does $R({\rho}_{k\epsilon})$ perform and how fast it converges remain largely open for general problems.
7.2.2 Stability
A second way to understand the generalization ability of deep learning is through the stability of SGD. An algorithm is considered stable if a slight change of the input does not alter the output much. It has long been observed that a stable algorithm has a small generalization gap; examples include $k$ nearest neighbors [101, 29], bagging [18, 19], etc. The precise connection between stability and generalization gap is stated by [17, 112]. In what follows, we formalize the idea of stability and its connection with the generalization gap. Let $\mathcal{A}$ denote an algorithm (possibly randomized) which takes a sample $S\triangleq {\{({y}_{i},{\bm{x}}_{i})\}}_{1\le i\le n}$ of size $n$ and returns an estimated parameter $\widehat{\bm{\theta}}\triangleq \mathcal{A}(S)$. Following [49], we have the following definition for stability.
Definition 2.
An algorithm (possibly randomized) $\mathrm{A}$ is $\epsilon $uniformly stable with respect to the loss function $\mathrm{L}\mathit{}\mathrm{(}\mathrm{\cdot}\mathrm{,}\mathrm{\cdot}\mathrm{)}$ if for all datasets $S\mathrm{,}{S}^{\mathrm{\prime}}$ of size $n$ which differ in at most one example, one has
$$\underset{\bm{x},y}{sup}{\mathbb{E}}_{\mathcal{A}}\left[\mathcal{L}(f(\bm{x};\mathcal{A}\left(S\right)),y)\mathcal{L}(f(\bm{x};\mathcal{A}\left({S}^{\prime}\right)),y)\right]\le \epsilon .$$ 
Here the expectation is taken w.r.t. the randomness in the algorithm $\mathrm{A}$ and $\epsilon $ might depend on $n$. The loss function $\mathrm{L}\mathit{}\mathrm{(}\mathrm{\cdot}\mathrm{,}\mathrm{\cdot}\mathrm{)}$ takes an example (say $\mathrm{(}\mathbf{x}\mathrm{,}y\mathrm{)}$) and the estimated parameter (say $\mathrm{A}\mathit{}\mathrm{(}S\mathrm{)}$) as inputs and outputs a real value.
Surprisingly, an $\epsilon $uniformly stable algorithm incurs small generalization gap in expectation, which is stated in the following lemma.
Lemma 1 (Theorem 2.2 in [49]).
Let $\mathrm{A}$ be $\epsilon $uniformly stable. Then the expected generalization gap is no larger than $\epsilon $, i.e.,
$$\left{\mathbb{E}}_{\mathcal{A},S}\left[\frac{1}{n}\sum _{i=1}^{n}\mathcal{L}(f({\bm{x}}_{i};\mathcal{A}\left(S\right)),{y}_{i}){\mathbb{E}}_{(\bm{x},y)\sim \mathcal{D}}\left[\mathcal{L}(f(\bm{x};\mathcal{A}\left(S\right)),y)\right]\right]\right\le \epsilon .$$  (39) 
With Lemma 1 in hand, it suffices to prove stability bound on specific algorithms. It turns out that SGD introduced in Section 6 is uniformly stable when solving smooth nonconvex functions.
Theorem 8 (Theorem 3.12 in [49]).
Assume that for any fixed $\mathrm{(}y\mathrm{,}\mathbf{x}\mathrm{)}$, the loss function $\mathrm{L}\mathit{}\mathrm{(}f\mathit{}\mathrm{(}x\mathrm{;}\mathbf{\theta}\mathrm{)}\mathrm{,}y\mathrm{)}$, viewed as a function of $\mathbf{\theta}$, is $L$Lipschitz and $\beta $smooth. Consider running SGD on the empirical loss function with decaying step size ${\alpha}_{t}\mathrm{\le}c\mathrm{/}t$, where $c$ is some small absolute constant. Then SGD is uniformly stable with
$$\epsilon \lesssim \frac{{T}^{1\frac{1}{\beta c+1}}}{n},$$ 
where we have ignored the dependency on $\beta \mathrm{,}c$ and $L$.
Theorem 8 reveals that SGD operating on nonconvex loss functions is indeed uniformly stable as long as the number of steps $T$ is not large compared with $n$. This together with Lemma 1 demonstrates the generalization ability of SGD in expectation. Nevertheless, two important limitations are worth mentioning. First, Lemma 1 provides an upper bound on the outofsample error in expectation, but ideally, instead of an onaverage guarantee under ${\mathbb{E}}_{\mathcal{A},S}$, we would like to have a high probability guarantee as in the convex case [37]. Second, controlling the generalization gap alone is not enough to achieve a small outofsample error, since it is unclear whether SGD can achieve a small training error within $T$ steps.
7.2.3 Implicit regularization
In the presence of overparametrization (number of parameters larger than the sample size), conventional wisdom informs us that we should apply some regularization techniques (e.g., ${\mathrm{\ell}}_{1}/{\mathrm{\ell}}_{2}$ regularization) so that the model will not overfit the data. However, in practice, neural networks without explicit regularization generalize well. This phenomenon motivates researchers to look at the regularization effects introduced by training algorithms (e.g., SGD) in this overparametrized regime. While there might exits multiple, if not infinite global minima of the empirical loss (26), it is possible that practical algorithms tend to converge to solutions with better generalization powers.
Take the underdetermined linear system $\bm{X}\bm{\theta}=\bm{y}$ as a starting point. Here $\bm{X}\in {\mathbb{R}}^{n\times p}$ and $\bm{\theta}\in {\mathbb{R}}^{p}$ with $p$ much larger than $n$. Running gradient descent on the loss $\frac{1}{2}{\parallel \bm{X}\bm{\theta}\bm{y}\parallel}_{2}^{2}$ from the origin (i.e., ${\bm{\theta}}^{0}=\mathrm{\U0001d7ce}$) results in the solution with the minimum Euclidean norm, that is GD converges to
$\underset{\bm{\theta}\in {\mathbb{R}}^{p}}{\mathrm{min}}$  $\mathrm{\hspace{1em}}{\parallel \bm{\theta}\parallel}_{2}\mathit{\hspace{1em}\hspace{1em}}\text{subject to}\mathit{\hspace{1em}}\bm{X}\bm{\theta}=\bm{y}.$ 
In words, without any ${\mathrm{\ell}}_{2}$ regularization in the loss function, gradient descent automatically finds the solution with the least ${\mathrm{\ell}}_{2}$ norm. This phenomenon, often called as implicit regularization, not only has been empirically observed in training neural networks, but also has been theoretically understood in some simplified cases, e.g., logistic regression with separable data. In logistic regression, given a training set ${\{({y}_{i},{\bm{x}}_{i})\}}_{1\le i\le n}$ with ${\bm{x}}_{i}\in {\mathbb{R}}^{p}$ and ${y}_{i}\in \{1,1\}$, one aims to fit a logistic regression model by solving the following program:
$$\underset{\bm{\theta}\in {\mathbb{R}}^{p}}{\mathrm{min}}\mathit{\hspace{1em}\hspace{1em}}\frac{1}{n}\sum _{i=1}^{n}\mathrm{\ell}\left({y}_{i}{\bm{x}}_{i}^{\top}{\bm{\theta}}^{t}\right).$$  (40) 
Here, $\mathrm{\ell}(u)\triangleq \mathrm{log}(1+{e}^{u})$ denotes the logistic loss. Further assume that the data is separable, i.e., there exists ${\bm{\theta}}^{*}\in {\mathbb{R}}^{p}$ such that ${y}_{i}{\bm{\theta}}^{*\top}{\bm{x}}_{i}>0$ for all $i$. Under this condition, the loss function (40) can be arbitrarily close to zero for certain $\bm{\theta}$ with ${\parallel \bm{\theta}\parallel}_{2}\to \mathrm{\infty}$. What happens when we minimize (40) using gradient descent? [118] uncovers a striking phenomenon.
Theorem 9 (Theorem 3 in [118]).
Consider the logistic regression (40) with separable data. If we run GD
$${\bm{\theta}}^{t+1}={\bm{\theta}}^{t}\eta \frac{1}{n}\sum _{i=1}^{n}{y}_{i}{\bm{x}}_{i}{\mathrm{\ell}}^{\prime}\left({y}_{i}{\bm{x}}_{i}^{\top}{\bm{\theta}}^{t}\right)$$ 
from any initialization ${\mathbf{\theta}}^{\mathrm{0}}$ with appropriate step size $\eta \mathrm{>}\mathrm{0}$, then normalized ${\mathbf{\theta}}^{t}$ converges to a solution with the maximum ${\mathrm{\ell}}_{\mathrm{2}}$ margin. That is,
$$\underset{t\to \mathrm{\infty}}{lim}\frac{{\bm{\theta}}^{t}}{{\parallel {\bm{\theta}}^{t}\parallel}_{2}}=\widehat{\bm{\theta}},$$  (41) 
where $\widehat{\mathbf{\theta}}$ is the solution to the hard margin support vector machine:
$$\widehat{\bm{\theta}}\triangleq \mathrm{arg}\underset{\bm{\theta}\in {\mathbb{R}}^{p}}{\mathrm{min}}{\parallel \bm{\theta}\parallel}_{2},\mathit{\text{subject to}}\mathit{\hspace{1em}}{y}_{i}{\bm{x}}_{i}^{\top}\bm{\theta}\ge 1\mathit{\hspace{1em}}\mathit{\text{for all}}1\le i\le n.$$  (42) 
The above theorem reveals that gradient descent, when solving logistic regression with separable data, implicitly regularizes the iterates towards the ${\mathrm{\ell}}_{2}$ max margin vector (cf. (41)), without any explicit regularization as in (42). Similar results have been obtained by [62]. In addition, [47] studied algorithms other than gradient descent and showed that coordinate descent produces a solution with the maximum ${\mathrm{\ell}}_{1}$ margin.
Moving beyond logistic regression, which can be viewed as a onelayer neural net, the theoretical understanding of implicit regularization in deeper neural networks is still limited; see [48] for an illustration in deep linear convolutional neural networks.
8 Discussion
Due to space limitations, we have omitted several important deep learning models; notable examples include deep reinforcement learning [86], deep probabilistic graphical models [108], variational autoencoders [66], transfer learning [132], etc. Apart from the modeling aspect, interesting theories on generative adversarial networks [10, 11], recurrent neural networks [3], connections with kernel methods [59, 9] are also emerging. We have also omitted the inverseproblem view of deep learning where the data are assumed to be generated from a certain neural net and the goal is to recover the weights in the NN with as few examples as possible. Various algorithms (e.g., GD with spectral initialization) have been shown to recover the weights successfully in some simplified settings [134, 117, 42, 87, 23, 39].
In the end, we identify a few important directions for future research.

•
New characterization of data distributions. The success of deep learning relies on its power of efficiently representing complex functions relevant to real data. Comparatively, classical methods often have optimal guarantee if a problem has a certain known structure, such as smoothness, sparsity, and lowrankness [120, 31, 20, 24], but they are insufficient for complex data such as images. How to characterize the highdimensional real data that can free us from known barriers, such as the curse of dimensionality is an interesting open question?

•
Understanding various computational algorithms for deep learning. As we have emphasized throughout this survey, computational algorithms (e.g., variants of SGD) play a vital role in the success of deep learning. They allow fast training of deep neural nets and probably contribute towards the good generalization behavior of deep learning in practice. Understanding these computational algorithms and devising better ones are crucial components in understanding deep learning.

•
Robustness. It has been well documented that DNNs are sensitive to small adversarial perturbations that are indistinguishable to humans [123]. This raises serious safety issues once if deploy deep learning models in applications such as selfdriving cars, healthcare, etc. It is therefore crucial to refine current training practice to enhance robustness in a principled way [115].

•
Low SNRs. Arguably, for image data and audio data where the signaltonoise ratio (SNR) is high, deep learning has achieved great success. In many other statistical problems, the SNR may be very low. For example, in financial applications, the firm characteristic and covariates may only explain a small part of the financial returns; in healthcare systems, the uncertainty of an illness may not be predicted well from a patient’s medical history. How to adapt deep learning models to excel at such tasks is an interesting direction to pursue?
Acknowledgements
We thank Ruying Bao, Yuxin Chen, Chenxi Liu, Weijie Su, Qingcan Wang and Pengkun Yang for helpful comments and discussions.
References
 [1] Martín Abadi and et. al. TensorFlow: Largescale machine learning on heterogeneous systems, 2015. Software available from tensorflow.org.
 [2] Reza AbbasiAsl, Yuansi Chen, Adam Bloniarz, Michael Oliver, Ben DB Willmore, Jack L Gallant, and Bin Yu. The deeptune framework for modeling and characterizing neurons in visual cortex area v4. bioRxiv, page 465534, 2018.
 [3] Zeyuan AllenZhu and Yuanzhi Li. Can SGD Learn Recurrent Neural Networks with Provable Generalization? ArXiv eprints, abs/1902.01028, 2019.
 [4] Zeyuan AllenZhu, Yuanzhi Li, and Zhao Song. A convergence theory for deep learning via overparameterization. arXiv preprint arXiv:1811.03962, 2018.
 [5] Martin Anthony and Peter L Bartlett. Neural network learning: Theoretical foundations. cambridge university press, 2009.
 [6] Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. 70:214–223, 06–11 Aug 2017.
 [7] Vladimir I Arnold. On functions of three variables. Collected Works: Representations of Functions, Celestial Mechanics and KAM Theory, 1957–1965, pages 5–8, 2009.
 [8] Sanjeev Arora and Boaz Barak. Computational complexity: a modern approach. Cambridge University Press, 2009.
 [9] Sanjeev Arora, Simon S Du, Wei Hu, Zhiyuan Li, and Ruosong Wang. Finegrained analysis of optimization and generalization for overparameterized twolayer neural networks. arXiv preprint arXiv:1901.08584, 2019.
 [10] Sanjeev Arora, Rong Ge, Yingyu Liang, Tengyu Ma, and Yi Zhang. Generalization and equilibrium in generative adversarial nets (gans). In Proceedings of the 34th International Conference on Machine LearningVolume 70, pages 224–232. JMLR. org, 2017.
 [11] Yu Bai, Tengyu Ma, and Andrej Risteski. Approximability of discriminators implies diversity in gans. arXiv preprint arXiv:1806.10586, 2018.
 [12] Andrew R Barron. Universal approximation bounds for superpositions of a sigmoidal function. IEEE Transactions on Information theory, 39(3):930–945, 1993.
 [13] Peter L Bartlett. The sample complexity of pattern classification with neural networks: the size of the weights is more important than the size of the network. IEEE transactions on Information Theory, 44(2):525–536, 1998.
 [14] Peter L Bartlett, Dylan J Foster, and Matus J Telgarsky. Spectrallynormalized margin bounds for neural networks. In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 6240–6249. Curran Associates, Inc., 2017.
 [15] Benedikt Bauer and Michael Kohler. On deep learning as a remedy for the curse of dimensionality in nonparametric regression. Technical report, Technical report, 2017.
 [16] Léon Bottou. Online learning and stochastic approximations. Online learning in neural networks, 17(9):142, 1998.
 [17] Olivier Bousquet and André Elisseeff. Stability and generalization. Journal of machine learning research, 2(Mar):499–526, 2002.
 [18] Leo Breiman. Bagging predictors. Machine learning, 24(2):123–140, 1996.
 [19] Leo Breiman et al. Heuristics of instability and stabilization in model selection. The annals of statistics, 24(6):2350–2383, 1996.
 [20] Emmanuel J Candès and Terence Tao. The power of convex relaxation: Nearoptimal matrix completion. arXiv preprint arXiv:0903.1476, 2009.
 [21] Chensi Cao, Feng Liu, Hai Tan, Deshou Song, Wenjie Shu, Weizhong Li, Yiming Zhou, Xiaochen Bo, and Zhi Xie. Deep learning and its applications in biomedicine. Genomics, proteomics & bioinformatics, 16(1):17–32, 2018.
 [22] Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. arXiv preprint arXiv:1806.07366, 2018.
 [23] Yuxin Chen, Yuejie Chi, Jianqing Fan, and Cong Ma. Gradient descent with random initialization: Fast global convergence for nonconvex phase retrieval. Mathematical Programming, pages 1–33, 2019.
 [24] Yuxin Chen, Yuejie Chi, Jianqing Fan, Cong Ma, and Yuling Yan. Noisy matrix completion: Understanding statistical guarantees for convex relaxation via nonconvex optimization. arXiv preprint arXiv:1902.07698, 2019.
 [25] Lenaic Chizat and Francis Bach. On the global convergence of gradient descent for overparameterized models using optimal transport. In Advances in neural information processing systems, pages 3040–3050, 2018.
 [26] Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio. Learning phrase representations using rnn encoderdecoder for statistical machine translation. arXiv preprint arXiv:1406.1078, 2014.
 [27] R Dennis Cook et al. Fisher lecture: Dimension reduction in regression. Statistical Science, 22(1):1–26, 2007.
 [28] Jeffrey De Fauw, Joseph R Ledsam, Bernardino RomeraParedes, Stanislav Nikolov, Nenad Tomasev, Sam Blackwell, Harry Askham, Xavier Glorot, Brendan O’Donoghue, Daniel Visentin, et al. Clinically applicable deep learning for diagnosis and referral in retinal disease. Nature medicine, 24(9):1342, 2018.
 [29] Luc Devroye and Terry Wagner. Distributionfree performance bounds for potential function rules. IEEE Transactions on Information Theory, 25(5):601–604, 1979.
 [30] David L Donoho. Highdimensional data analysis: The curses and blessings of dimensionality. AMS math challenges lecture, 1(2000):32, 2000.
 [31] David L Donoho and Jain M Johnstone. Ideal spatial adaptation by wavelet shrinkage. biometrika, 81(3):425–455, 1994.
 [32] Simon S Du, Jason D Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai. Gradient descent finds global minima of deep neural networks. arXiv preprint arXiv:1811.03804, 2018.
 [33] John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of Machine Learning Research, 12(Jul):2121–2159, 2011.
 [34] Weinan E, Chao Ma, and Qingcan Wang. A priori estimates of the population risk for residual networks. arXiv preprint arXiv:1903.02154, 2019.
 [35] Ronen Eldan and Ohad Shamir. The power of depth for feedforward neural networks. In Conference on Learning Theory, pages 907–940, 2016.
 [36] Jianqing Fan and Runze Li. Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American statistical Association, 96(456):1348–1360, 2001.
 [37] Vitaly Feldman and Jan Vondrak. High probability generalization bounds for uniformly stable algorithms with nearly optimal rate. arXiv preprint arXiv:1902.10710, 2019.
 [38] Jerome H Friedman and Werner Stuetzle. Projection pursuit regression. Journal of the American statistical Association, 76(376):817–823, 1981.
 [39] Haoyu Fu, Yuejie Chi, and Yingbin Liang. Local geometry of onehiddenlayer neural networks for logistic regression. arXiv preprint arXiv:1802.06463, 2018.
 [40] Kunihiko Fukushima and Sei Miyake. Neocognitron: A selforganizing neural network model for a mechanism of visual pattern recognition. In Competition and cooperation in neural nets, pages 267–285. Springer, 1982.
 [41] Chao Gao, Jiyi Liu, Yuan Yao, and Weizhi Zhu. Robust estimation and generative adversarial nets. arXiv preprint arXiv:1810.02030, 2018.
 [42] Surbhi Goel, Adam Klivans, and Raghu Meka. Learning one convolutional layer with overlapping patches. arXiv preprint arXiv:1802.02547, 2018.
 [43] Noah Golowich, Alexander Rakhlin, and Ohad Shamir. Sizeindependent sample complexity of neural networks. arXiv preprint arXiv:1712.06541, 2017.
 [44] Gene H Golub and Charles F Van Loan. Matrix computations. JHU Press, 4 edition, 2013.
 [45] Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning. MIT Press, 2016.
 [46] Ian Goodfellow, Jean PougetAbadie, Mehdi Mirza, Bing Xu, David WardeFarley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. In Advances in neural information processing systems, pages 2672–2680, 2014.
 [47] Suriya Gunasekar, Jason Lee, Daniel Soudry, and Nathan Srebro. Characterizing implicit bias in terms of optimization geometry. arXiv preprint arXiv:1802.08246, 2018.
 [48] Suriya Gunasekar, Jason D Lee, Daniel Soudry, and Nati Srebro. Implicit bias of gradient descent on linear convolutional networks. In Advances in Neural Information Processing Systems, pages 9482–9491, 2018.
 [49] Moritz Hardt, Benjamin Recht, and Yoram Singer. Train faster, generalize better: Stability of stochastic gradient descent. arXiv preprint arXiv:1509.01240, 2015.
 [50] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
 [51] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Identity mappings in deep residual networks. In European conference on computer vision, pages 630–645. Springer, 2016.
 [52] Geoffrey Hinton, Nitish Srivastava, and Kevin Swersky. Neural networks for machine learning lecture 6a overview of minibatch gradient descent. 2012.
 [53] Geoffrey E Hinton, Nitish Srivastava, Alex Krizhevsky, Ilya Sutskever, and Ruslan R Salakhutdinov. Improving neural networks by preventing coadaptation of feature detectors. arXiv preprint arXiv:1207.0580, 2012.
 [54] Sepp Hochreiter and Jürgen Schmidhuber. Long shortterm memory. Neural computation, 9(8):1735–1780, 1997.
 [55] Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger. Densely connected convolutional networks. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4700–4708, 2017.
 [56] David H Hubel and Torsten N Wiesel. Receptive fields, binocular interaction and functional architecture in the cat’s visual cortex. The Journal of physiology, 160(1):106–154, 1962.
 [57] Forrest N Iandola, Song Han, Matthew W Moskewicz, Khalid Ashraf, William J Dally, and Kurt Keutzer. Squeezenet: Alexnetlevel accuracy with 50x fewer parameters and< 0.5 mb model size. arXiv preprint arXiv:1602.07360, 2016.
 [58] Sergey Ioffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. arXiv preprint arXiv:1502.03167, 2015.
 [59] Arthur Jacot, Franck Gabriel, and Clément Hongler. Neural tangent kernel: Convergence and generalization in neural networks. In Advances in neural information processing systems, pages 8580–8589, 2018.
 [60] Prateek Jain, Sham M Kakade, Rahul Kidambi, Praneeth Netrapalli, and Aaron Sidford. Accelerating stochastic gradient descent. arXiv preprint arXiv:1704.08227, 2017.
 [61] Adel Javanmard, Marco Mondelli, and Andrea Montanari. Analysis of a twolayer neural network via displacement convexity. arXiv preprint arXiv:1901.01375, 2019.
 [62] Ziwei Ji and Matus Telgarsky. Risk and parameter convergence of logistic regression. arXiv preprint arXiv:1803.07300, 2018.
 [63] Rahul Kidambi, Praneeth Netrapalli, Prateek Jain, and Sham Kakade. On the insufficiency of existing momentum schemes for stochastic optimization. In 2018 Information Theory and Applications Workshop (ITA), pages 1–9. IEEE, 2018.
 [64] Jack Kiefer, Jacob Wolfowitz, et al. Stochastic estimation of the maximum of a regression function. The Annals of Mathematical Statistics, 23(3):462–466, 1952.
 [65] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
 [66] Diederik P Kingma and Max Welling. Autoencoding variational bayes. arXiv preprint arXiv:1312.6114, 2013.
 [67] Jason M Klusowski and Andrew R Barron. Risk bounds for highdimensional ridge function combinations including neural networks. arXiv preprint arXiv:1607.01434, 2016.
 [68] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton. Imagenet classification with deep convolutional neural networks. In Advances in neural information processing systems, pages 1097–1105, 2012.
 [69] Harold Kushner and G George Yin. Stochastic approximation and recursive algorithms and applications, volume 35. Springer Science & Business Media, 2003.
 [70] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton. Deep learning. nature, 521(7553):436, 2015.
 [71] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradientbased learning applied to document recognition. Proceedings of the IEEE, 86(11):2278–2324, 1998.
 [72] Hao Li, Zheng Xu, Gavin Taylor, Christoph Studer, and Tom Goldstein. Visualizing the loss landscape of neural nets. In Advances in Neural Information Processing Systems, pages 6391–6401, 2018.
 [73] KerChau Li. Sliced inverse regression for dimension reduction. Journal of the American Statistical Association, 86(414):316–327, 1991.
 [74] Xingguo Li, Junwei Lu, Zhaoran Wang, Jarvis Haupt, and Tuo Zhao. On tighter generalization bound for deep neural networks: Cnns, resnets, and beyond. arXiv preprint arXiv:1806.05159, 2018.
 [75] Yujia Li, Kevin Swersky, and Rich Zemel. Generative moment matching networks. In International Conference on Machine Learning, pages 1718–1727, 2015.
 [76] Tengyuan Liang. How well can generative adversarial networks (gan) learn densities: A nonparametric view. arXiv preprint arXiv:1712.08244, 2017.
 [77] Henry W Lin, Max Tegmark, and David Rolnick. Why does deep and cheap learning work so well? Journal of Statistical Physics, 168(6):1223–1247, 2017.
 [78] Min Lin, Qiang Chen, and Shuicheng Yan. Network in network. arXiv preprint arXiv:1312.4400, 2013.
 [79] Andrew L Maas, Awni Y Hannun, and Andrew Y Ng. Rectifier nonlinearities improve neural network acoustic models. In Proc. icml, volume 30, page 3, 2013.
 [80] VE Maiorov and Ron Meir. On the near optimality of the stochastic approximation of smooth functions by neural networks. Advances in Computational Mathematics, 13(1):79–103, 2000.
 [81] Yuly Makovoz. Random approximants and neural networks. Journal of Approximation Theory, 85(1):98–109, 1996.
 [82] Song Mei, Theodor Misiakiewicz, and Andrea Montanari. Meanfield theory of twolayers neural networks: dimensionfree bounds and kernel limit. arXiv preprint arXiv:1902.06015, 2019.
 [83] Song Mei, Andrea Montanari, and PhanMinh Nguyen. A mean field view of the landscape of twolayer neural networks. Proceedings of the National Academy of Sciences, 115(33):E7665–E7671, 2018.
 [84] Hrushikesh Mhaskar, Qianli Liao, and Tomaso Poggio. Learning functions: when is deep better than shallow. arXiv preprint arXiv:1603.00988, 2016.
 [85] Hrushikesh N Mhaskar. Neural networks for optimal approximation of smooth and analytic functions. Neural computation, 8(1):164–177, 1996.
 [86] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A Rusu, Joel Veness, Marc G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg Ostrovski, et al. Humanlevel control through deep reinforcement learning. Nature, 518(7540):529, 2015.
 [87] Marco Mondelli and Andrea Montanari. On the connection between learning twolayers neural networks and tensor decomposition. arXiv preprint arXiv:1802.07301, 2018.
 [88] Yurii E Nesterov. A method for solving the convex programming problem with convergence rate o (1/k^ 2). In Dokl. Akad. Nauk SSSR, volume 269, pages 543–547, 1983.
 [89] Behnam Neyshabur, Ryota Tomioka, and Nathan Srebro. Normbased capacity control in neural networks. In Conference on Learning Theory, pages 1376–1401, 2015.
 [90] Sebastian Nowozin, Botond Cseke, and Ryota Tomioka. fgan: Training generative neural samplers using variational divergence minimization. In Advances in Neural Information Processing Systems, pages 271–279, 2016.
 [91] Ian Parberry. Circuit complexity and neural networks. MIT press, 1994.
 [92] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in pytorch. 2017.
 [93] Allan Pinkus. Approximation theory of the mlp model in neural networks. Acta numerica, 8:143–195, 1999.
 [94] Tomaso Poggio, Hrushikesh Mhaskar, Lorenzo Rosasco, Brando Miranda, and Qianli Liao. Why and when can deepbut not shallownetworks avoid the curse of dimensionality: a review. International Journal of Automation and Computing, 14(5):503–519, 2017.
 [95] Boris T Polyak. Some methods of speeding up the convergence of iteration methods. USSR Computational Mathematics and Mathematical Physics, 4(5):1–17, 1964.
 [96] Boris T Polyak and Anatoli B Juditsky. Acceleration of stochastic approximation by averaging. SIAM Journal on Control and Optimization, 30(4):838–855, 1992.
 [97] Boris Teodorovich Polyak and Yakov Zalmanovich Tsypkin. Adaptive estimation algorithms: convergence, optimality, stability. Avtomatika i Telemekhanika, (3):71–84, 1979.
 [98] Christopher Poultney, Sumit Chopra, Yann L Cun, et al. Efficient learning of sparse representations with an energybased model. In Advances in neural information processing systems, pages 1137–1144, 2007.
 [99] Sashank J Reddi, Satyen Kale, and Sanjiv Kumar. On the convergence of adam and beyond. 2018.
 [100] Herbert Robbins and Sutton Monro. A stochastic approximation method. The Annals of Mathematical Statistics, 22(3):400–407, 1951.
 [101] William H Rogers and Terry J Wagner. A finite sample distributionfree performance bound for local discrimination rules. The Annals of Statistics, pages 506–514, 1978.
 [102] David Rolnick and Max Tegmark. The power of deeper networks for expressing natural functions. arXiv preprint arXiv:1705.05502, 2017.
 [103] Yaniv Romano, Matteo Sesia, and Emmanuel J Candès. Deep knockoffs. arXiv preprint arXiv:1811.06687, 2018.
 [104] Grant M Rotskoff and Eric VandenEijnden. Neural networks as interacting particle systems: Asymptotic convexity of the loss landscape and universal scaling of the approximation error. arXiv preprint arXiv:1805.00915, 2018.
 [105] David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams. Learning internal representations by error propagation. Technical report, California Univ San Diego La Jolla Inst for Cognitive Science, 1985.
 [106] Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, Alexander C. Berg, and Li FeiFei. ImageNet Large Scale Visual Recognition Challenge. International Journal of Computer Vision (IJCV), 115(3):211–252, 2015.
 [107] Haşim Sak, Andrew Senior, and Françoise Beaufays. Long shortterm memory recurrent neural network architectures for large scale acoustic modeling. In Fifteenth annual conference of the international speech communication association, 2014.
 [108] Ruslan Salakhutdinov and Geoffrey Hinton. Deep boltzmann machines. In Artificial intelligence and statistics, pages 448–455, 2009.
 [109] Tim Salimans, Ian Goodfellow, Wojciech Zaremba, Vicki Cheung, Alec Radford, and Xi Chen. Improved techniques for training gans. In Advances in Neural Information Processing Systems, pages 2234–2242, 2016.
 [110] Johannes SchmidtHieber. Nonparametric regression using deep neural networks with relu activation function. arXiv preprint arXiv:1708.06633, 2017.
 [111] Shai ShalevShwartz and Shai BenDavid. Understanding machine learning: From theory to algorithms. Cambridge university press, 2014.
 [112] Shai ShalevShwartz, Ohad Shamir, Nathan Srebro, and Karthik Sridharan. Learnability, stability and uniform convergence. Journal of Machine Learning Research, 11(Oct):2635–2670, 2010.
 [113] David Silver, Julian Schrittwieser, Karen Simonyan, Ioannis Antonoglou, Aja Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton, et al. Mastering the game of go without human knowledge. Nature, 550(7676):354, 2017.
 [114] Bernard W Silverman. Density estimation for statistics and data analysis. Chapman & Hall, CRC, 1998.
 [115] Chandan Singh, W James Murdoch, and Bin Yu. Hierarchical interpretations for neural network predictions. arXiv preprint arXiv:1806.05337, 2018.
 [116] Justin Sirignano and Konstantinos Spiliopoulos. Mean field analysis of neural networks. arXiv preprint arXiv:1805.01053, 2018.
 [117] Mahdi Soltanolkotabi. Learning relus via gradient descent. In Advances in Neural Information Processing Systems, pages 2007–2017, 2017.
 [118] Daniel Soudry, Elad Hoffer, Mor Shpigel Nacson, Suriya Gunasekar, and Nathan Srebro. The implicit bias of gradient descent on separable data. The Journal of Machine Learning Research, 19(1):2822–2878, 2018.
 [119] David A Sprecher. On the structure of continuous functions of several variables. Transactions of the American Mathematical Society, 115:340–355, 1965.
 [120] Charles J Stone. Optimal global rates of convergence for nonparametric regression. The annals of statistics, pages 1040–1053, 1982.
 [121] Ilya Sutskever, James Martens, George Dahl, and Geoffrey Hinton. On the importance of initialization and momentum in deep learning. In International conference on machine learning, pages 1139–1147, 2013.
 [122] Christian Szegedy, Wei Liu, Yangqing Jia, Pierre Sermanet, Scott Reed, Dragomir Anguelov, Dumitru Erhan, Vincent Vanhoucke, and Andrew Rabinovich. Going deeper with convolutions. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1–9, 2015.
 [123] Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus. Intriguing properties of neural networks. arXiv preprint arXiv:1312.6199, 2013.
 [124] Matus Telgarsky. Benefits of depth in neural networks. arXiv preprint arXiv:1602.04485, 2016.
 [125] Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288, 1996.
 [126] VN Vapnik and A Ya Chervonenkis. On the uniform convergence of relative frequencies of events to their probabilities. Theory of Probability & Its Applications, 16(2):264–280, 1971.
 [127] Pascal Vincent, Hugo Larochelle, Yoshua Bengio, and PierreAntoine Manzagol. Extracting and composing robust features with denoising autoencoders. In Proceedings of the 25th international conference on Machine learning, pages 1096–1103. ACM, 2008.
 [128] Stefan Wager, Sida Wang, and Percy S Liang. Dropout training as adaptive regularization. In Advances in neural information processing systems, pages 351–359, 2013.
 [129] E Weinan, Jiequn Han, and Arnulf Jentzen. Deep learningbased numerical methods for highdimensional parabolic partial differential equations and backward stochastic differential equations. Communications in Mathematics and Statistics, 5(4):349–380, 2017.
 [130] Ashia C Wilson, Rebecca Roelofs, Mitchell Stern, Nati Srebro, and Benjamin Recht. The marginal value of adaptive gradient methods in machine learning. In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 4148–4158. Curran Associates, Inc., 2017.
 [131] Yonghui Wu, Mike Schuster, Zhifeng Chen, Quoc V Le, Mohammad Norouzi, Wolfgang Macherey, Maxim Krikun, Yuan Cao, Qin Gao, Klaus Macherey, et al. Google’s neural machine translation system: Bridging the gap between human and machine translation. arXiv preprint arXiv:1609.08144, 2016.
 [132] Jason Yosinski, Jeff Clune, Yoshua Bengio, and Hod Lipson. How transferable are features in deep neural networks? In Advances in neural information processing systems, pages 3320–3328, 2014.
 [133] Jason Yosinski, Jeff Clune, Anh Nguyen, Thomas Fuchs, and Hod Lipson. Understanding neural networks through deep visualization. arXiv preprint arXiv:1506.06579, 2015.
 [134] Kai Zhong, Zhao Song, Prateek Jain, Peter L Bartlett, and Inderjit S Dhillon. Recovery guarantees for onehiddenlayer neural networks. In Proceedings of the 34th International Conference on Machine LearningVolume 70, pages 4140–4149. JMLR. org, 2017.