Towards a Mathematical Understanding of Neural Network-Based Machine Learning: what we know and what we don't

Introduction

Neural network-based machine learning is both very powerful and very fragile. On the one hand, it can be used to approximate functions in very high dimensions with the efficiency and accuracy never possible before. This has opened up brand new possibilities in a wide spectrum of different disciplines. On the other hand, it has got the reputation of being somewhat of a “black magic”: Its success depends on lots of tricks, and parameter tuning can be quite an art. The main objective for a mathematical study of machine learning is to

explain the reasons behind the success and the subtleties, and

propose new models that are equally successful but much less fragile.

We are still quite far from completely achieving these goals but it is fair to say that a reasonable big picture is emerging.

The purpose of this article is to review the main achievements towards the first goal and discuss the main remaining puzzles. In the tradition of good old applied mathematics, we will not only give attention to rigorous mathematical results, but also discuss the insight we have gained from careful numerical experiments as well as the analysis of simplified models.

At the moment much less attention has been given to the second goal. One proposal that we should mention is the continuous formulation advocated in . The idea there is to first formulate “well-posed” continuous models of machine learning problems and then discretize to get concrete algorithms. What makes this proposal attractive is the following:

many existing machine learning models and algorithms can be recovered in this way in a scaled form;

there is evidence suggesting that indeed machine learning models obtained this way is more robust with respect to the choice of hyper-parameters than conventional ones (see for example Figure 5 below);

new models and algorithms are borne out naturally in this way. One particularly interesting example is the maximum principle-based training algorithm for ResNet-like models .

However, at this stage one cannot yet make the claim that the continuous formulation is the way to go. For this reason we will postpone a full discussion of this issue to future work.

The model problem of supervised learning which we focus on in this article can be formulated as follows: given a dataset $S=\{({\bm{x}}_{i},y_{i}=f^{*}({\bm{x}}_{i})),i\in[n]\}$, approximate $f^{*}$ as accurately as we can. If $f^{*}$ takes continuous values, this is called a regression problem. If $f^{*}$ takes discrete values, this is called a classification problem.

We will focus on the regression problem. For simplicity, we will neglect the so-called “measurement noise” since it does not change much the big picture that we will describe, even though it does matter for a number of important specific issues. We will assume ${\bm{x}}_{i}\in X:=^{d}$, and we denote by $P$ the distribution of $\{{\bm{x}}_{i}\}$. We also assume for simplicity that $\sup_{{\bm{x}}\in X}|f^{*}({\bm{x}})|\leq 1$.

Obviously this is a problem of function approximation. As such, it can either be regarded as a problem in numerical analysis or a problem in statistics. We will take the former viewpoint since it is more in line with the algorithmic and analysis issues that we will study.

The standard procedure for supervised learning is as follows:

Choose a loss function. Our primary goal is to fit the data. Therefore the most popular choice is the “empirical risk”:

Sometimes one adds some regularization terms.

Choose an optimization algorithm and the hyper-parameters. The most popular choices are gradient descent (GD), stochastic gradient descent (SGD) and advanced optimizers such as Adam , RMSprop .

The overall objective is to minimize the “population risk”, also known as the “generalization error”:

In practice, this is estimated on a finite data set (which is unrelated to any data used to train the model) and called test error, whereas the empirical risk (which is used for training purposes) is called the training error.

2 The main issues of interest

From a mathematical perspective, there are three important issues that we need to study:

Properties of the hypothesis space. In particular, what kind of functions can be approximated efficiently by a particular machine learning model? What can we say about the generalization gap, i.e. the difference between training and test errors.

Properties of the loss function. The loss function defines the variational problem used to find the solution to the machine learning problem. Questions such as the landscape of the variational problem are obviously important. The landscape of neural network models is typically non-convex, and there may exist many saddle points and bad local minima.

Properties of the training algorithm. Two obvious questions are: Can we optimize the loss function using the selected training algorithm? Does the solutions obtained from training generalize well to test data?

The second and third issues are closely related. In the under-parametrized regime (when the size of the training dataset is larger than the number of free parameters in the hypothesis space), this loss function largely determines the solution of the machine learning model. In the opposite situation, the over-parametrized regime, this is no longer true. Indeed it is often the case that there are infinite number of global minimizers of the loss function. Which one is picked depends on the details of the training algorithm.

The most important parameters that we should keep in mind are:

Typically we are interested in the situation when $m,n,t\rightarrow\infty$ and $d\gg 1$.

3 Approximation and estimation errors

Denote by $\hat{f}$ the output of the machine learning (abbreviated ML) model. Let

We can decompose the error $f^{*}-\hat{f}$ into:

$f^{*}-f_{m}$ is the approximation error, due entirely to the choice of the hypothesis space. $f_{m}-\hat{f}$ is the estimation error, the additional error due to the fact that we only have a finite dataset.

To get some basic idea about the approximation error, note that classically when approximating functions using polynomials, piecewise polynomials, or truncated Fourier series, the error typically satisfies

where $H^{\alpha}$ denotes the Sobolev space of order $\alpha$. The appearance of $1/d$ in the exponent of $m$ is a signature of an important phenomenon, the curse of dimensionality (CoD): The number of parameters required to achieve certain accuracy depends exponentially on the dimension. For example, if we want $m^{-\alpha/d}=0.1$, then we need $m=10^{d/\alpha}=10^{d}$, if $\alpha=1$.

At this point, it is useful to recall the one problem that has been extensively studied in high dimension: the problem of evaluating an integral, or more precisely, computing an expectation. Let $g$ be a function defined on $X$. We are interested in computing approximately

where the expectation is taken with respect to the uniform distribution. Typical grid-based quadrature rules, such as the Trapezoidal rule and the Simpson’s rule, all suffer from CoD. The one algorithm that does not suffer from CoD is the Monte Carlo algorithm which works as follows. Let $\{{\bm{x}}_{i}\}_{i=1}^{n}$ be a set of independent, uniformly distributed random variables on $X$. Let

The $O(1/\sqrt{n})$ rate is independent of $d$. It turns out that this rate is almost the best one can hope for.

In practice, $\mbox{Var}(g)$ can be very large in high dimension. Therefore, variance reduction techniques are crucial in order to make Monte Carlo methods truly practical.

Turning now to the estimation error. Our concern is how the approximation produced by the machine learning algorithm behaves away from the training dataset, or in practical terms, whether the training and test errors are close. Shown in Figure 1 is the classical Runge phenomenon for interpolating functions on uniform grids using high order polynomials. One can see that while on the training set, here the grid points, the error of the interpolant is 0, away from the training set, the error can be very large.

It is often easier to study a related quantity, the generalization gap. Consider the solution that minimizes the empirical risk, $\hat{f}=\operatorname{argmin}_{f\in\mathcal{H}_{m}}\hat{\mathcal{R}}_{n}(f)$. The “generalization gap” of $\hat{f}$ is the quantity $|{\mathcal{R}}(\hat{f})-\hat{\mathcal{R}}_{n}(\hat{f})|$. Since it is equal to $|I(g)-I_{n}(g)|$ with $g({\bm{x}})=(\hat{f}({\bm{x}})-f^{*}({\bm{x}}))^{2}$, one might be tempted to conclude that

based on (1). This is NOT necessarily true since $\hat{f}$ is highly correlated with $\{{\bm{x}}_{i}\}$. In fact, controlling this gap is among the most difficult problems in ML.

Studying the correlations of $\hat{f}$ is a rather impossible problem. Therefore to estimate the generalization gap, we resort to the uniform bound:

The RHS of this equation depends heavily on the nature of $\mathcal{H}_{m}$. If we take $\mathcal{H}_{m}$ to be the unit ball in the Lipschitz space, we have

with $g=(h-f^{*})^{2}$. This gives rise to CoD for the size of the dataset (commonly referred to as “sample complexity”). However if we take $\mathcal{H}_{m}$ to be the unit ball in the Barron space, to be defined later, we have

This is the kind of estimates that we should look for.

Assuming that all the functions under consideration are bounded, the problem of estimating the RHS of (2) reduces to the estimation of $\sup_{h\in\mathcal{H}_{m}}|I(h)-I_{n}(h)|$. One way to do this is to use the notion of Rademacher complexity .

Let $\mathcal{F}$ be a set of functions, and $S=({\bm{x}}_{1},{\bm{x}}_{2},...,{\bm{x}}_{n})$ be a set of data points. Then, the Rademacher complexity of $\mathcal{F}$ with respect to $S$ is defined as

where $\{\xi_{i}\}_{i=1}^{n}$ are i.i.d. random variables taking values $\pm 1$ with equal probability.

Roughly speaking, Rademacher complexity quantifies the degree to which functions in the hypothesis space can fit random noise on the given dataset. It bounds the quantity of interest, $\sup_{h\in\mathcal{H}}|I(h)-I_{n}(h)|$, from above and below.

For any $\delta\in(0,1)$, with probability at least $1-\delta$ over the random samples $S=({\bm{x}}_{1},\cdots,{\bm{x}}_{n})$, we have

For a proof, see for example [76, Theorem 26.5]. For this reason, a very important part of theoretical machine learning is to study the Rademacher complexity of a given hypothesis space.

It should be noted that there are other ways of analyzing the generalization gap, such as the stability method .

Preliminary remarks

In this section we set the stage for our discussion by going over some of the classical results as well as recent qualitative studies that are of general interest.

Let us consider the two-layer neural network hypothesis space:

where $\sigma$ is a nonlinear function, the activation function. The Universal Approximation Theorem (UAT) states that under very mild conditions, any continuous function can be uniformly approximated on compact domains by neural network functions.

[24, Theorem 5] If $\sigma$ is sigmoidal, in the sense that $\lim_{z\rightarrow-\infty}\sigma(z)=0,\lim_{z\rightarrow\infty}\sigma(z)=1$, then any function in $C(^{d})$ can be approximated uniformly by two layer neural network functions.

This result can be extended to any activation functions that are not exactly a polynomial .

UAT plays the role of Weierstrass Theorem on the approximation of continuous functions by polynomials. It is obviously an important result, but it is of limited use due to the lack of quantitative information about the error in the approximation. For one thing, the same conclusion can be drawn for polynomial approximation, which we know is of limited use in high dimension.

Many quantitative error estimates have been established since then. But most of these estimates suffer from CoD. The one result that stands out is the estimate proved by Barron :

2 The loss landscape of large neural network models

The landscape of the loss function for large neural networks was studied in using an analogy with high dimensional spherical spin glasses and numerical experiments. The landscape of high dimensional spherical spin glass models has been analyzed in . It was shown that the lowest critical values of the Hamiltonians of these models form a layered structure, ordered by the indices of the critical points. They are located in a well-defined band lower-bounded by the global minimum. The probability of finding critical points outside the band diminishes exponentially with the dimension of the spin-glass model. Choromanska et al used the observation that neural networks with independent weights can be mapped onto spherical spin glass models, and suggested that the same picture should hold qualitatively for large neural network models. They provided numerical evidence to support this suggestion.

This work is among the earliest that suggests even though it is highly non-convex, the landscape of larger neural networks might be simpler than the ones for smaller networks.

3 Over-parametrization, interpolation and implicit regularization

Modern deep learning often works in the over-parametrized regime where the number of free parameters is larger than the size of the training dataset. This is a new territory as far as machine learning theory is concerned. Conventional wisdom would suggest that one should expect overfitting in this regime, i.e. an increase in the generalization gap. Whether overfitting really happens is a problem of great interest in theoretical machine learning. As we will see later, one can show that, overfitting does happen in the highly over-parametrized regime.

An enlightening numerical study of the optimization and generalization properties in this regime was carried out in . Among other things, it was discovered that in this regime, the neural networks are so expressive that they can fit any data, no matter how much noise one adds to the data. Later it was shown by Cooper that the set of global minima with no training error forms a manifold of dimension $m-n$ .

Some of these global minima generalize very poorly. Therefore an important question is how to select the ones that do generalize well. It was suggested in that by tuning the hyper-parameters of the optimization algorithm, one can obtain models with good generalization properties without adding explicit regularization. This means that the training dynamics itself has some implicit regularization mechanism which ensures that bad global minima are not selected during training. Understanding the possible mechanism for this implicit regularization is one of the key issues in understanding modern machine learning.

4 The selection of topics

There are several rather distinct ways for a mathematical analysis of machine learning and particularly neural network models:

The numerical analysis perspective. Basically machine learning problems are viewed as (continuous) function approximation and optimization problems, typically in high dimension.

The harmonic analysis perspective. Deep learning is studied from the viewpoint of building hierarchical wavelets-like transforms. Examples of such an approach can be found in .

The statistical physics perspective. A particularly important tool is the replica trick. For special kinds of problems, this approach allows us to perform asymptotically exact (hence sharp) calculations. It is also useful to study the performance of models and algorithms from the viewpoint of phase transitions. See for example .

The information theory perspective. See for an example of the kinds of results obtained along this line.

The PAC learning theory perspective. This is closely related to the information theory perspective. It studies machine learning from the viewpoint of complexity theory. See for example .

It is not yet clear how the different approaches are connected, even though many results may fall in several categories. In this review, we will only cover the results obtained along the lines of numerical analysis. We encourage the reader to consult the papers referenced above to get a taste of these alternative perspectives.

Within the numerical analysis perspective, supervised machine learning and neural network models are still vast topics. By necessity, this article focusses on a few aspects which we believe to be key problems in machine learning. As a model problem, we focus on $L^{2}$-regression here, where the data is assumed to be of the form $({\bm{x}}_{i},f^{*}({\bm{x}}_{i}))$ without uncertainty in the $y$-direction.

In Section 3, we focus on the function spaces developed for neural network models. Section 4 discusses the (very short) list of results available for the energy landscape of loss functionals in machine learning. Training dynamics for network weights are discussed in Section 5 with a focus on gradient descent. The specific topics are selected for two criteria:

We believe that they have substantial importance for the mathematical understanding of machine learning.

We are reasonably confident that the mathematical models developed so far will over time find their way into the standard language of the theoretical machine learning community.

There are many glaring omissions in this article, among them:

Classification problems. While most benchmark problems for neural networks fall into the framework of classification, we focus on the more well-studied area of regression.

Some other common neural network architectures, such as convolutional neural networks, long short-term memory (LSTM) networks, encoder-decoder networks.

The impact of stochasticity. Many training algorithms and initialization schemes for neural networks use random variables. While toy models with standard Gaussian noise are well understood, the realistic case remains out of reach .

Simplified neural networks such as linear and quadratic networks. While these models are not relevant for applications, the simpler model allows for simpler analysis. Some insight is available for these models which has not been achieved in the non-linear case .

Important “tricks”, such as dropout, batch normalization, layer normalization, gradient clipping, etc. To our knowledge, these remain empirically useful, but mysterious from a mathematical perspective.

Highly data-dependent results. The geometry of data is a large field which we avoid in this context. While many data-distributions seem to be concentrated close to relatively low-dimensional manifolds in a much higher-dimensional ambient space, these ‘data-manifolds’ are in many cases high-dimensional enough that CoD, a central theme in this review, is still an important issue.

The approximation property and the Rademacher complexity of the hypothesis space

The most important issue in classical approximation theory is to identify the function space naturally associated with a particular approximation scheme, e.g. approximation by piecewise polynomials on regular grids. These spaces are typically some Sobolev or Besov spaces, or their variants. They are the natural spaces for the particular approximation scheme, since one can prove matching direct and inverse approximation theorems, namely any function in the space can be approximated using the given approximation scheme with the specified rate of convergence, and conversely any function that can be approximated to the specified order of accuracy belongs to that function space.

Machine learning is just another way to approximate functions, therefore we can ask similar questions, except that our main interest in this case is in high dimension. Any machine learning model hits the ‘curse of dimensionality’ when approximating the class of Lipschitz functions. Nevertheless, many important problems seem to admit accurate approximations by neural networks. Therefore it is important to understand the class of functions that can be well approximated by a particular machine learning model.

There is one important difference from the classical setting. In high dimension, the rate of convergence is limited to the Monte Carlo rate and its variants. There is limited room regarding order of convergence and consequently there is no such thing as the order of the space as is the case for Sobolev spaces.

Ideally, we would like to accomplish the following:

Given a type of hypothesis space $\mathcal{H}_{m}$, say two-layer neural networks, identify the natural function space associated with them (in particular, identify a norm $\|f^{*}\|_{*}$) that satisfies:

Inverse approximation theorem: If a function $f^{*}$ can be approximated efficiently by the functions in $\mathcal{H}_{m}$, as $m\rightarrow\infty$ with some uniform bounds, then $\|f^{*}\|_{*}$ is finite.

Study the generalization gap for this function space. One way of doing this is to study the Rademacher complexity of the set $\mathcal{F}_{Q}=\{f\,:\,\|f\|_{*}\leq Q\}$. Ideally, we would like to have:

If both holds, then a combination gives us, up to logarithmic terms:

It should be noted that what we are really interested in is the quantitative measures of the target function that control the approximation and estimation errors. We call these quantities “norms” but we are not going to insist that they are really norms. In addition, we would like to use one norm to control both the approximation and estimation errors. This way we have one function space that meets both requirements. However, it could very well be the case that we need different quantities to control different errors. See the discussion about residual networks below. This means that we will be content with a generalized version of (3):

We will see that this can indeed be done for the most popular neural network models.

Let $\phi(\cdot;\bm{w})$ be the feature function parametrized by $\bm{w}$, e.g. $\phi({\bm{x}};\bm{w})=\sigma(\bm{w}^{T}{\bm{x}})$. A random feature model is given by

Denote by $\mathcal{H}_{k}$ this RKHS. Then for any $f\in\mathcal{H}_{k}$, there exists $a(\cdot)\in L^{2}(\pi_{0})$ such that

Assume $f^{*}\in\mathcal{H}_{k}$, then there exists $a(\cdot)$, such that

Let $(\bm{w}_{j}^{0})_{j=0}^{\infty}$ be a sequence of i.i.d. random variables drawn from $\pi_{0}$. Let $f^{*}$ be a continuous function on $X$. Assume that there exist constants $C$ and a sequence $(a_{j})_{j=0}^{\infty}$ satisfying $\sup_{j}|a_{j}|\leq C$, such that

To see how these approximation theory results can play out in a realistic ML setting, consider the regularized model:

For any $\delta\in(0,1)$, with probability $1-\delta$, the population risk of the regularized estimator satisfies

These results should be standard. However, they do not seem to be available in the literature. In the appendix, we provide a proof for these results.

It is worth noting that the dependence on $\|f\|_{\infty}$ and $\log(n/\delta)$ can be removed by a more sophisticated analysis . However, to achieve the rate of $O(1/m+1/\sqrt{n})$, one must make an explicit assumption on the decay rate of eigenvalues of the corresponding kernel operator .

2 Two-layer neural network model

The hypothesis space for two-layer neural networks is defined by:

We will focus on the case when the activation function $\sigma$ is ReLU: $\sigma(z)=\max(z,0)$. Many of the results discussed below can be extended to more general activation functions .

Functions in $\mathcal{B}_{p}$ are called Barron functions. As shown in , for the ReLU activation function, we actually have $\|\cdot\|_{\mathcal{B}_{p}}=\|\cdot\|_{\mathcal{B}_{q}}$ for any $1\leq p\leq q\leq\infty$. Hence, we will use $\|\cdot\|_{\mathcal{B}}$ and $\mathcal{B}$ denote the Barron norm and Barron space, respectively.

Barron space and Barron functions are named in reference to the article which was the first to recognize and rigorously establish the advantages of non-linear approximation over linear approximation by considering neural networks with a single hidden layer.

It should be stressed that the Barron norm introduced above is not the same as the one used in , which was based on the Fourier transform (see (11)). To highlight this distinction, we will call the kind of norm in (11) spectral norm.

An important property of the ReLU activation function is the homogeneity property $\sigma(\lambda z)=\lambda\,\sigma(z)$ for all $\lambda>0$. A discussion of the representation for Barron functions with partial attention to homogeneity can be found in . Barron spaces for other activation functions are discussed in .

One natural question is what kind of functions are Barron functions. The following result gives a partial answer.

where $\hat{f}$ is the Fourier transform of $f$, then $f$ can be represented as

where $\sigma(x)=\max(0,x)$. Moreover $\|f\|_{\mathcal{B}}\leq 2\Delta(f)+2\|\nabla f(0)\|_{1}+2f(0)$.

On the other hand, every Barron function is Lipschitz-continuous. An important criterion to establish that certain functions are not in Barron space is the following structure theorem.

As a consequence, distance functions to curved surfaces are not Barron functions.

The claim that the Barron space is the natural space associated with two-layer networks is justified by the following series of results.

We present a brief self-contained proof of the $L^{\infty}$-direct approximation theorem in the appendix. We believe the idea to be standard, but have been unable to locate a good reference for it.

In fact, there exists a constant $C>0$ such that

and for every $\varepsilon>0$ there exists $f\in\mathcal{B}$ such that

see . Further approximation results, including in classical functions spaces, can be found in .

Let $f^{*}$ be a continuous function. Assume there exists a constant $C$ and a sequence of functions $f_{m}\in\mathcal{N}_{C}$ such that

for all ${\bm{x}}\in X$, then there exists a probability distribution $\rho^{*}$ on $\Omega$, such that

for all ${\bm{x}}\in X$ and $\|f^{*}\|_{\mathcal{B}}\leq C$.

Both theorems are proved in . In addition, it just so happens that the Rademacher complexity is also controlled by a Monte Carlo like rate:

Let $\mathcal{F}_{Q}=\{f\in\mathcal{B},\|f\|_{\mathcal{B}}\leq Q\}$. Then we have

In the same way as before, one can now consider the regularized model:

: Assume $f^{*}:X\mapsto\in\mathcal{B}$. There exist constants absolute $C_{0}$, such that for any $\delta>0$, if $\lambda\geq C_{0}$, then with probability at least $1-\delta$ over the choice of the training set, we have

3 Residual networks

We use $\Theta:=\{{\bm{U}}_{1},\dots,{\bm{U}}_{L},{\bm{W}}_{l},\dots,{\bm{W}}_{L},\bm{\alpha}\}$ to denote all the parameters to be learned from data.

This ODE system can be viewed as the limit of the residual network (3.3) (). Consider the following linear ODEs ($p\geq 1$)

Let $f$ be a function that admits the form $f=f_{\bm{\alpha},\{\rho_{t}\}}$ for a pair of ($\bm{\alpha},\{\rho_{t}\}$), then we define

to be the $\mathcal{D}_{p}$ norm of $f$ with respect to the pair ($\bm{\alpha}$, $\{\rho_{t}\}$). Here $|\bm{\alpha}|$ is obtained from $\bm{\alpha}$ by taking element-wise absolute values. We define

to be the $\mathcal{D}_{p}$ norm of $f$, and let $\mathcal{D}_{p}=\{f:\|f\|_{\mathcal{D}_{p}}<\infty\}$ be the flow-induced function space.

Given a family of probability distribution $\{\rho_{t},\ t\in\}$, the “Lipschitz coefficient” of $\{\rho_{t}\}$, denoted by $\textit{Lip}_{\{\rho_{t}\}}$, is defined as the infimum of all the numbers $L$ that satisfies

for any $t,s\in$, where $\|\cdot\|_{1,1}$ is the sum of the absolute values of all the entries in a matrix. The “Lipschitz norm” of $\{\rho_{t}\}$ is defined as

Let $f$ be a function that satisfies $f=f_{\bm{\alpha},\{\rho_{t}\}}$ for a pair of ($\bm{\alpha},\{\rho_{t}\}$), then we define

One can easily see that the “norms” defined here are all non-negative quantities (despite the $-D$ term), even though it is not clear that they are really norms. The following embedding theorem shows that flow-induced function space is larger than Barron space.

Finally, we define a discrete “path norm” for residual networks.

For a residual network defined by (3.3) with parameters $\Theta=\{\bm{\alpha},{\bm{U}}_{l},{\bm{W}}_{l},l=0,1,\cdots,L-1\}$, we define the $l_{1}$ path norm of $\Theta$ to be

With the definitions above, we are ready to state the direct and inverse approximation theorems for the flow-induced function spaces .

there is an $L$-layer residual network $f_{L}(\cdot;\Theta)$ that satisfies

Let $f$ be a function defined on $X$. Assume that there is a sequence of residual networks $\{f_{L}(\cdot;\Theta_{L})\}_{L=1}^{\infty}$ such that $\|f({\bm{x}})-f_{L}({\bm{x}};\Theta_{L})\|\rightarrow 0$ as $L\rightarrow\infty$. Assume further that the parameters in $\{f_{L}(\cdot;\Theta)\}_{L=1}^{\infty}$ are (entry-wise) bounded by $c_{0}$. Then, we have $f\in\mathcal{D}_{\infty}$, and

Moreover, if there exists constant $c_{1}$ such that $\|f_{L}\|_{\mathcal{D}_{1}}\leq c_{1}$ holds for any $L>0$, then we have

The Rademacher complexity estimate is only established for a family of modified flow-induced function norms $\|\cdot\|_{\hat{\mathcal{D}}_{p}}$ (see the factor 2 in the definition below). It is not clear at this stage whether this is only a technical difficulty.

Denote by $\hat{\mathcal{D}}_{p}$ the space of functions with finite $\hat{\mathcal{D}}_{p}$ norm. Then, we have

Let $\hat{\mathcal{D}}_{p}^{Q}=\{f\in\hat{\mathcal{D}}_{p}:\|f\|_{\hat{\mathcal{D}}_{p}}\leq Q\}$, then we have

Next we turn to the generalization error estimates for the regularized estimator. At the moment, for the same reason as above, such estimates have only been proved when the empirical risk is regularized by a weighted path norm

which is the discrete version of (24). This norm assigns larger weights to paths that pass through more non-linearities. Now consider the residual network (3.3) and the regularized empirical risk:

Let $f^{*}:X\to$. Fix any $\lambda\geq 4+2/(3\sqrt{2\log(2d)})$. Assume that $\hat{\Theta}$ is an optimal solution of the regularized model (27). Then for any $\delta\in(0,1)$, with probability at least $1-\delta$ over the random training samples, the population risk satisfies

4 Multi-layer networks: Tree-like function spaces

A neural network with $L$ hidden layers is a function of the form

and the path-norm of the function as the infimum of the path-norm proxies over all weights inducing the function.

Rademacher complexity/generalization gap: Let $\mathcal{F}_{Q}=\{f\in C^{0},\|f\|_{\mathcal{W}^{L}}\leq Q\}$ be the closed ball of radius $Q>0$ in the tree-like function space. Then $\operatorname{Rad}_{S}(\mathcal{F}_{Q})\leq 2^{L+1}Q\sqrt{\frac{2\ln(2d+2)}{n}}$. If a stronger version of the path-norm is controlled, then the dependence on depth can be weakened to $L^{3}$ instead of $2^{L}$ . But remember, here we are thinking of keep $L$ fixed at some finite value and increase the widths of the layers.

In particular, an inverse approximation theorem holds: If $\|f_{m}\|_{\mathcal{W}^{L}}\leq C$ and $f_{m}\to f$ in $L^{2}(P)$, then $f$ is in the tree-like function space of the same depth.

Note that the network has $O(m^{2L-1})$ weights. Part of this stems from the recursive definition, and by rearranging the index set into a tree-like structure, the number of free parameters (but dimension-independent) can be dropped to $O(m^{L+1})$.

It is unclear whether a depth-independent (or less depth-dependent) approximation rate can be expected. Unlike two-layer networks and residual neural networks, layers of a multi-layer network discretize into conditional expectations, whereas the other function classes are naturally expressed as expectations. A larger number of parameters for multi-layer networks is therefore natural.

Consider the regularized loss function: regularized risk functional

In particular, there exists a function $f_{m}$ satisfying the risk estimate. Using a natural cut-off, the constant $\bar{c}$ can be replaced with $\|f^{*}\|_{L^{\infty}}\leq C\,\|f^{*}\|_{\mathcal{W}^{L}}$.

5 Indexed representation and multi-layer spaces

Neural networks used in applications are not tree-like. To capture the structure of the functions represented by practical multi-layer neural networks, introduced an indexed representation of neural network functions.

If the index spaces are finite, this coincides with finite neural networks and the integrals are replaced with finite sums. From this we see that the class of arbitrarily wide neural networks modeled on certain index spaces may not be a vector space. If all weight spaces $\{\Omega_{i}\}$ are sufficiently expressive (e.g. the unit interval with Lebesgue measure), then the set of multi-layer networks modeled on $\Omega_{L},\dots,\Omega_{0}$ becomes a vector space. The key property is that $(0,1)$ can be decomposed into two sets of probability $1/2$, each of which is isomorphic to itself, so adding two neural networks of equal depth is possible.

The space of arbitrarily wide neural networks is a subspace of the tree-like function space of depth $L$ that consists of functions $f$ with a finite path-norm

Since the coefficients of non-consecutive layers do not share an index, this may be a proper subspace for networks with mutiple hidden layers. If $\Omega_{i}=(0,1)$, the space of arbitrarily wide networks with one hidden layer coincides with Barron space.

The space of measurable weight functions which render the path-norm finite is inconveniently large when considering training dynamics. To allow a natural gradient flow structure, we consider the subset of functions with $L^{2}$-weights. This is partially motivated by the observation that the $L^{2}$-norm of weights controls the path-norm.

We sketch the proof for two hidden layers.

Note that the proof is entirely specific to network-like architectures and does not generalize to tree-like structures. We define the measure of complexity of a function (which is not a norm) as

We can equip the class of neural networks modeled on index spaces $\{\Omega_{i}\}$ with a metric which respects the parameter structure.

The space of arbitrarily wide neural networks with $L^{2}$-weights can be metrized with the Hilbert-weight metric

The normalization across layers is required to ensure that functions in which one layer can be chosen identical do not have zero distance by shifting all weight to the one layer.

We refer to these metric spaces (metric vector spaces if $\Omega_{i}=(0,1)$ for all $i\geq 1$) as multi-layer spaces. They are complete metric spaces.

6 Depth separation in multi-layer networks

We can ask how much larger $L$-layer space is compared to $(L-1)$-layer space. A satisfying answer to this question is still outstanding, but partial answers have been found, mostly concerning the differences between networks with one and two hidden layers.

The structure theorem for Barron functions Theorem 10 shows that functions which are non-differentiable on a curved hypersurface are not Barron. In particular, this includes distance functions from hypersurfaces like

It is obvious, however, that $f$ is the composition of two Barron functions and therefore can be represented exactly by a neural network with two hidden layers. This criterion is easy to check in practice and therefore of greater potential impact than mere existence results. But it says nothing about approximation by finite neural networks.

A separation result like this can also be obtained with standard activation functions. If $f$ is a Barron function, then there exists an $f_{m}({\bm{x}})=\sum_{i=1}^{m}a_{i}\,\sigma(\bm{w}_{i}^{T}x)$ such that

In , the authors show that there exists a function $f$ such that

There are some results on functions which can be approximated better with significantly deeper networks than few hidden layers, but a systematic picture is still missing. To the best of our knowledge, there are no results for the separation between $L$ and $L+1$ hidden layers.

7 Tradeoffs between learnability and approximation

This is a typical phenomenon shared by all machine learning models of low complexity. Let ${P}_{d}$ be Lebesgue-measure on the $d$-dimensional unit cube (which we take as the archetype of a truly ‘high-dimensional’ data distribution).

Let $Z$ be a Banach space of functions such that the unit ball $B^{Z}$ in $Z$ satisfies

i.e. the Rademacher complexity on a set of $N$ sample decays at the optimal rate in the number of data points. Then

The Kolmogorov width of $Z$ in the space of Lipschitz functions with respect to the $L^{2}$-metric is low in the sense that

There exists a function $f$ with Lispchitz constant $1$ such that $f(0)=0$, but

This resembles the result of for approximation by ridge functions under a constraint on the number of parameters, whereas here a complexity bound is assumed instead and no specific form of the model is prescribed (and the result thus applies to multi-layer networks as well).

Thus function spaces of low complexity are ‘poor approximators’ for general classes like Lipschitz functions since we need functions of large $Z$-norm to approximate functions to a prescribed level of accuracy. This includes all function spaces discussed in this review, although some spaces are significantly larger than others (e.g. there is a large gap between reproducing kernel Hilbert spaces, Barron space, and tree-like three layer space).

8 A priori vs. a posteriori estimates

The error estimate given above should be compared with a more typical form of estimate in the machine learning literature:

where $\|\hat{\theta}_{n}\|$ is some suitably defined norm. Aside from the fact that (16) gives a bound on the total generalization error and (34) gives a bound on the generalization gap, there is an additional important difference: The right hand side of (16) depends only on the target function $f^{*}$, not the output of the machine learning model. The right hand side of (34) depends only on the output of the machine learning model, not the target function. In accordance with the practice in finite element methods, we call (16) a priori estimates and (34) a posteriori estimates.

How good and how useful are these estimates? A priori estimates discussed here tell us in particular that there exist functions in the hypothesis space for which the generalization error does not suffer from the CoD if the target function lies in the appropriate function space. It is likely that these estimates are nearly optimal in the sense that they are comparable to Monte Carlo error rates (except for multi-layer neural networks, see below). It is possible to improve these estimates, for example using standard tricks for Monte Carlo sampling for the approximation error and local Rademacher complexity for the estimation error . However, these would only improve the exponents in $m$ and $n$ by $O(1/d)$ which diminishes for large $d$.

Regarding the quantitative value of these a priori estimates, the situation is less satisfactory. The first issue is that the values of the norms are not known since the target function is not known. One can estimate these values using the output of the machine learning model, but this does not give us rigorous bounds. An added difficulty is that the norms are defined as an infimum over all possible representations, the output of the machine learning model only gives one representation. But even if we use the exact values of these norms, the bounds given above are still not tight. For one thing, the use of Monte Carlo sampling to control the approximation error does not give a tight bound. This by itself is an interesting issue.

The obvious advantage of the a posteriori bounds is that they can be readily evaluated and give us quantitative bounds for the size of the generalization gap. Unfortunately this has not been borned out in practice: The values of these norms are so enormous that these bounds are almost always vacuous .

In finite element methods, a posteriori estimates are used to help refining the mesh in adaptive methods. Ideally one would like to do the same for machine learning models. However, little has been done in this direction.

Since the a posteriori bounds only controls the generalization gap, not the full generalization error, it misses an important aspect of the whole picture, namely, the approximation error. In fact, by choosing a very strong norm, one can always obtain estimates of the type in (34). However, with such strongly constrained hypothesis space, the approximation error might be huge. This is indeed the case for some of the norm-based a posteriori estimates in the literature. See for examples.

9 What’s not known?

Here is a list of problems that we feel are most pressing.

1. Sharper estimates. There are two obvious places where one should be able to improve the estimates.

In the current analysis, the approximation error is estimated with the help of Monte Carlo sampling. This gives us the typical size of the error for randomly picked parameters. However, in machine learning, we are only interested in the smallest error. This is a clean mathematical problem that has not received attention.

The use of Rademacher complexity to bound the generalization gap neglects the fact that the integrand in the definition of the population risk (as well as the empirical risk) should itself be small, since it is the point-wise error. This should be explored further. We refer to for some results in this direction.

2. The rate for the approximation error for functions in multi-layer spaces is not the same as Monte Carlo. Can this be improved?

More generally, it is not clear whether the multi-layer spaces defined earlier are the right spaces for multi-layer neural networks.

4. Function space for convolutional neural networks that fully explores the benefit of symmetry.

5. Another interesting network structure is the DenseNet . Naturally one is interested in the natural function space associated with DenseNets.

The loss function and the loss landscape

It is a surprise to many that simple gradient descent algorithms work quite well for optimizing the loss function in common machine learning models. To put things into perspective, no one would dream of using the same kind of algorithms for protein folding – the energy landscape for am typical protein is so complicated with lots of local minima that gradient descent algorithms will not go very far. The fact that they seem to work well for machine learning models strongly suggests that the landscapes for the loss functions are qualitatively different. An important question is to quantify exactly how the loss landscape looks like. Unfortunately, theoretical results on this important problem is still quite scattered and there is not yet a general picture that has emerged. But generally speaking, the current understanding is that while it is possible to find arbitrarily bad examples for finite sized neural networks, their landscape simplifies as the size increases.

To begin with, the loss function is non-convex and it is easy to cook up models for which the loss landscape has bad local minima (see for example ). Moreover, it has been suggested that for small size two-layer ReLU networks with teacher networks as the target function, nearly all target networks lead to spurious local minima, and the probability of hitting such local minima is quite high . It has also been suggested the over-parametrization helps to avoid these bad spurious local minima.

The loss landscape of large networks can be very complicated. presented some amusing numerical results in which the authors demonstrated that one can find arbitrarily complex patterns near the global minima of the loss function. Some theoretical results along this direction were proved in . Roughly speaking, it was shown that for any $\varepsilon>0$, every low-dimensional pattern can be found in a loss surface of a sufficiently deep neural network, and within the pattern there exists a point whose loss is within $\varepsilon$ of the global minimum .

On the positive side, a lot is known for linear and quadratic neural network models. For linear neural network models, it has been shown that : (1) every local minimum is a global minimum, (2) every critical point that is not a global minimum is a saddle point, (3) for networks with more than three layers there exist “bad” saddle points where the Hessian has no negative eigenvalue and (4) there are no such bad saddle points for networks with three layers.

Similar results have been obtained for over-parametrized two-layer neural network models with quadratic activation . In this case it has been shown under various conditions that (1) all local minima are global and (2) all saddle points are strict, namely there are directions of strictly negative curvature. Another interesting work for the case of quadratic activation function is . In the case when the target function is a “single neuron”, gives an asymptotically exact (as $d\rightarrow\infty$) characterization of the number of training data samples needed for the global minima of the empirical loss to give rise to a unique function, namely the target function.

For networks where half of the weights in every layer can be set to zero if the remaining weights are rescaled appropriately, the set of global minimizers is connected .

We also mention the interesting empirical work reported in . Among other things, it was demonstrated that adding skip connection has a drastic effect on smoothing the loss landscape.

We still lack a good mathematical tool to describe the landscape of the loss function for large neural networks. In particular, are there local minima and how large is the basin of attraction of these local minima if they do exist?

For fixed, finite dimensional gradient flows, knowledge about the landscape allows us to draw conclusions about the qualitative behavior of the gradient descent dynamics independent of the detailed dynamics. In machine learning, the dimensionality of the loss function is $m$, the number of free parameters, and we are interested in the limit as $m$ goes to infinity. So it is tempting to ask about the landscape of the limiting (infinite dimensional) problem. It is not clear whether this can be formulated as a well-posed mathematical problem.

The training process: convergence and implicit regularization

The results of Section 3 tell us that good solutions do exist in the hypothesis space. The amazing thing is that simple minded gradient descent algorithms are able to find them, even though one might have to be pretty good at parameter tuning. In comparison, one would never dream of using gradient descent to perform protein folding, since the landscape of protein folding is so complicated with lots of bad local minima.

The basic questions about the training process are:

Optimization: Does the training process converge to a good solution? How fast?

Generalization: Does the solution selected by the training process generalize well? In particular, is there such thing as “implicit regularization”? What is the mechanism for such implicit regularization?

At the moment, we are still quite far from being able to answering these questions completely, but an intuitive picture has started to emerge.

We will mostly focus on the gradient descent (GD) training dynamics. But we will touch upon some important qualitative features of other training algorithms such as stochastic gradient descent (SGD) and Adam.

“Mean-field” is a notion in statistical physics that describes a particular form of interaction between particles. In the mean-field situation, particles interact with each other only through a mean-field which every particle contributes to more or less equally. The most elegant mean-field picture in machine learning is found in the case of two-layer neural networks: If one views the neurons as interacting particles, then these particles only interact with each other through the function represented by the neural network, the mean-field in this case. This observation was first made in . By taking the hydrodynamic limit for the gradient flow of finite neuron systems, these authors obtained a continuous integral differential equation that describes the evolution of the probability measure for the weights associated with the neurons.

then the GD dynamics (35) can be expressed equivalently as:

(36) is the mean-field equation that describes the evolution of the probability distribution for the weights associated with each neuron. The lemma above simply states that (36) is satisfied for finite neuron systems.

It is well-known that (36) is the gradient flow of $\hat{\mathcal{R}}_{n}$ under the Wasserstein metric. This brings the hope that the mathematical tools developed in the theory of optimal transport can be brought to bear for the analysis of (36) . In particular, we would like to use these tools to study the qualitative behavior of the solutions of (36) as $t\rightarrow\infty$. Unfortunately the most straightforward application of the results from optimal transport theory requires that the risk functional be displacement convex , a property that rarely holds in machine learning (see however the example in ). As a result, less than expected has been obtained using optimal transport theory.

The one important result, due originally to Chizat and Bach , is the following. We will state the result for the population risk.

Let $\{\rho_{t}\}$ be a solution of the Wasserstein gradient flow such that

$\rho_{0}$ is a probability distribution on the cone $\Theta:=\{|a|^{2}\leq|w|^{2}\}$.

Every open cone in $\Theta$ has positive measure with respect to $\rho_{0}$.

The velocity potentials $\frac{\delta\mathcal{R}}{\delta\rho}(\rho_{t},\cdot)$ converge to a unique limit as $t\to\infty$.

$\mathcal{R}(\rho_{t})$ decays to the global infimum value as $t\to\infty$.

If either condition is met, the unique limit of $\mathcal{R}(\rho_{t})$ is zero. If $\rho_{t}$ also converges in the Wasserstein metric, then the limit $\rho_{\infty}$ is a minimizer.

Intuitively, the theorem is slightly stronger than the statement that $\mathcal{R}(\rho_{t})$ converges to its infimum value if and only if its derivative converges to zero in a suitable sense (if we approach from a flow of omni-directional distributions). The theorem is more a statement about a training algorithm than the energy landscape, and specific PDE arguments are used in its proof. A few remarks are in order:

There are further technical conditions for the theorem to hold.

Convergence of subsequences of $\frac{\delta\mathcal{R}}{\delta\rho}(\rho_{t},\cdot)$ is guaranteed by compactness. The question of whether they converge to a unique limit can be asked independently of the initial distribution and therefore may be more approachable by standard means.

The first assumption on $\rho_{0}$ is a smoothness assumption needed for the existence of the gradient flow.

The second assumption on $\rho_{0}$ is called omni-directionality. It ensures that $\rho$ can shift mass in any direction which reduces risk. The cone formulation is useful due to the homogeneity of ReLU.

This is almost the only non-trivial rigorous result known on the global convergence of gradient flows in the nonlinear regime. In addition, it reveals the fact that having full support for the probability distribution (or something to that effect) is an important property that helps for the global convergence.

The result is insensitive to whether a minimizer of the risk functional exists (i.e. whether the target function is in Barron space). If the target function is not in Barron space, convergence may be very slow since the Barron norm increases sub-linearly during gradient descent training.

[87, Lemma 3.3] If $\{\rho_{t}\}$ evolves by the Wasserstein-gradient flow of $\mathcal{R}$, then

In high dimensions, even reasonably smooth functions are difficult to approximate with functions of low Barron norm in the sense of Theorem 34. Thus a “dynamic curse of dimensionality” may affect gradient descent training if the target function is not in Barron space.

Consider population and empirical risk expressed by the functionals

where the points $\{{\bm{x}}_{i}\}$ are i.i.d samples from the uniform distribution on $X$. There exists $f^{*}$ with Lipschitz constant and $L^{\infty}$-norm bounded by $1$ such that the parameter measures $\{\rho_{t}\}$ defined by the $2$-Wasserstein gradient flow of either $\hat{\mathcal{R}}_{n}$ or $\mathcal{R}$ satisfy

2 Two-layer neural networks with conventional scaling

In practice, people often use the conventional scaling (instead of the mean-field scaling) which takes the form:

With this “scaling”, the one case where a lot is known is the so-called highly over-parametrized regime. There is both good and bad news in this regime. The good news is that one can prove exponential convergence to global minima of the empirical risk.

Let $\lambda_{n}=\lambda_{\min}(K)$ and assume $\beta=0$. For any $\delta\in(0,1)$, assume that $m\gtrsim n^{2}\lambda_{n}^{-4}\delta^{-1}\ln(n^{2}\delta^{-1})$. Then with probability at least $1-6\delta$ we have

Now the bad news: the generalization property of the converged solution is no better than that of the associated random feature model, defined by freezing $\{\bm{b}_{j}\}=\{\bm{b}_{j}(0)\}$, and only training $\{a_{i}\}$.

The first piece of insight that the underlying dynamics in this regime is effectively linear is given in . termed the effective kernel the “neural tangent kernel”. Later it was proved rigorously that in this regime, the entire GD path for the two-layer neural network model is uniformly close to that of the associated random feature model .

In particular, there is no “implicit regularization” in this regime.

Even though convergence of the training process can be proved, overall this is a disappointing result. At the theoretical level, it does not shed any light on possible implicit regularization. At the practical level, it tells us that high over-parametrization is indeed a bad thing for these models.

What happens in practice? Do less over-parametrized regimes exist for which implicit regularization does actually happen? Some insight has been gained from the numerical study in .

Let us look at a simple example: the single neuron target function: $f_{1}^{*}({\bm{x}})=\sigma(\bm{b}^{*}\cdot{\bm{x}}),\quad\bm{b}^{*}=\bm{e}_{1}$. This is admittedly a very simple target function. Nevertheless, the training dynamics for this function is quite representative of target functions which can be accurately approximated by a small number of neurons (i.e. effectively “over-parametrized”).

Figure 3 shows some results from . First note that the bottom two figures represent exactly the kind of results stated in Theorem 40. The top two figures suggest that the training dynamics displays two phases. In the first phase, the training dynamics follows closely that of the associated random feature model. This phase quickly saturates. The training dynamics for the neural network model is able to reduce the training (and test) error further in the second phase through a quenching-activation process: Most neurons are quenched in the sense that their outer layer coefficients $a_{j}$ become very small, with the exception of only a few activated ones.

This can also be seen from the $m-n$ hyper-parameter space. Shown in Figure 4 are the heat maps of the test errors under the conventional and mean-field scaling, respectively. We see that the test error does not change much as $m$ changes for the mean-field scaling. In contrast, there is a clear “phase transition” in the heat map for the conventional scaling when the training dynamics undergoes a change from the neural network-like to a random feature-like behavior. This is further demonstrated in Figure 5: One can see that the performance of the neural network models indeed becomes very close to that of the random feature model as the network width $m$ increases. At the same time, the path norm also undergoes a sudden increase from the mildly over-parameterized/under-parameterized regime ($m\approx n/(d+1)$) to the highly over-parameterized regime ($m\approx n$). This may provide an explanation for the increase of the test error.

The existence of different kinds of behavior with drastically different generalization properties is one of the reasons behind the fragility of deep learning: If the network architecture happens to fall into the random feature-like regime, the performance will deteriorate.

3 Other convergence results for the training of neural network models

Convergence of GD for linear neural networks was proved in . analyzes the training of neural networks with quadratic activation function. It is proved that a modified GD with small initialization converges to the ground truth in polynomial time if the sample size $n\geq O(dr^{5})$ where $d,r$ are the input dimension and the rank of the target weight matrix, respectively. considers the learning of single neuron with SGD. It was proved that as long as the sample size is large enough, SGD converges to the ground truth exponentially fast if the model also consists of a single neuron. Similar analysis for the population risk was presented in .

4 Double descent and slow deterioration for the random feature model

At the continuous (in time) level, the training dynamics for the random feature model is a linear system of ODEs defined by the Gram matrix. It turns out that the generalization properties for this linear model is surprisingly complex, as we now discuss.

Consider a random feature model with features $\{\phi(\cdot;\bm{w})\}$ and probability distribution $\pi$ over the feature vectors $\bm{w}$. Let $\{\bm{w}_{1},\bm{w}_{2},...,\bm{w}_{m}\}$ be the random feature vectors sampled from $\pi$. Let $\Phi$ be an $n\times m$ matrix with $\Phi_{ij}=\phi({\bm{x}}_{i};\bm{w}_{j})$ where $\{{\bm{x}}_{1},{\bm{x}}_{2},...,{\bm{x}}_{n}\}$ is the training data, and let

where $\bm{a}=(a_{1},b_{2},...,a_{m})^{T}$are the parameters. To find $\bm{a}$, GD is used to optimize the following least squares objective function,

starting from the origin, where ${\bm{y}}$ is a vector containing the values of the target function at $\{{\bm{x}}_{i},\ i=1,2,...,n\}$. The dynamics of $\bm{a}$ is then given by

Let $\Phi=U\Sigma V^{T}$ be the singular value decomposition of $\Phi$, where

with $\{\lambda_{i}\}$ being the singular values of $\Phi$, in descending order. Then the GD solution of (40) at time $t\geq 0$ is given by

With this solution, we can conveniently compute the training and test error at any time $t$. Specifically, let $B=[\bm{w}_{1},...,\bm{w}_{m}]$, and with an abuse of notation let $\phi({\bm{x}};B)=(\phi({\bm{x}};\bm{w}_{1}),...,\phi({\bm{x}};\bm{w}_{m}))^{T}$, the prediction function at time $t$ is given by

Shown in the left figure of Figure 6 is the test error of the minimum norm solution (which is the limit of the GD path when initialized at 0) of the random feature model as a function of the size of the model $m$ for $n=500$ for the MNIST dataset . Also shown is the smallest eigenvalue of the Gram matrix. One can see that the test error peaks at $m=n$ where the smallest eigenvalue of the Gram matrix becomes exceedingly small. This is the same as the “double descent” phenomenon reported in . But as can be seen from the figure (and discussed in more detail below), it is really a resonance kind of phenomenon caused by the appearance of very small eigenvalues of the Gram matrix when $m=n$.

Shown in right is the test error of the solution when the gradient descent training is stopped after different number of steps. One curious thing is that when training is stopped at moderately large values of training steps, the resonance behavior does not show up much. Only when training is continued to very large number of steps does resonance show up.

Intuitively this is easy to understand. The large test error is caused by the small eigenvalues, because each term in (43) convergences to $1/\lambda_{i}$, and this limit is large when the corresponding eigenvalue is small. However, still by (43), the contribution of the eigenvalues enter through terms like $e^{-\lambda^{2}t}$, and therefore shows up very slowly in time.

Some rigorous analysis of this can be found in and will not be repeated here. Instead we show in Figure 7 an example of the detailed dynamics of the training process. One can see that the training dynamics can be divided into three regimes. In the first regime, the test error decreases from an $O(1)$ initial value to something small. In the second regime, which spans several decades, the test error remains small. It eventually becomes big in the third regime when the small eigenvalues start to contribute significantly. This slow deterioration phenomenon may also happen for more complicated models, like deep neural networks, considering that linear approximation of the loss function is effective near the global minimum, and hence in this regime the GD dynamics is similar to that of a linear model (e.g. random feature). Though, the conjecture needs to be confirmed by further theoretical and numerical study.

One important consequence of this resonance phenomenon is that it affects the training of two-layer neural networks under the conventional scaling. For example in Figure 4, the test error of neural networks under the conventional scaling peaks around $m=n$, with $m$ being the width of the network. Though at $m=n$ the neural network has more parameters than $n$, it still suffers from this resonance phenomenon since as we saw before, in the early phase of the training, the GD path for the neural network model follows closely that of the corresponding random feature model.

5 Global minima selection

In the over-parametrized regime, there usually exist many global minima. Different optimization algorithms may select different global minima. For example, it has been widely observed that SGD tends to pick “flatter solutions” than GD . A natural question is which global minima are selected by a particular optimization algorithm.

An interesting illustration of this is shown in Figure 8. Here GD was used to train the FashionMNIST dataset to near completion, and was suddenly replaced by SGD. Instead of finishing the last step in the previous training process, the subsequent SGD path escapes from the vicinity of the minimum that GD was converging to, and eventually converges to a different global minimum, with a slightly smaller test error .

This phenomenon can be well-explained by considering the dynamic stability of the optimizers, as was done in . It was found that the set of dynamically stable global minima is different for different training algorithm. In the example above, the global minimum that GD was converging to was unstable for SGD.

The gist of this phenomenon can be understood from the following simple one-dimensional optimization problem,

with $a_{i}\geq 0\,\,\forall i\in[n]$. The minimum is at $x=0$. For GD with learning rate $\eta$ to be stable, the following has to hold:

where $a=\sum_{i=1}^{n}a_{i}/n$, $s=\sqrt{\sum_{i=1}^{n}a_{i}^{2}/n-a^{2}}$. Therefore, for SGD to be stable at $x=0$, we not only need $|1-\eta a|\leq 1$, but also $(1-\eta a)^{2}+\eta^{2}s^{2}\leq 1$. In particular, SGD can only converge with the additional requirement that $s\leq 2/\eta$.

The quantity $a$ is called sharpness, and $s$ is called non-uniformity in . The above simple argument suggests that the global minima selected by SGD tend to be more uniform than the ones selected by GD.

This sharpness-non-uniformity criterion can be extended to multi-dimension. It turns out that this theoretical prediction is confirmed quite well by practical results. Figure 9 shows the sharpness and non-uniformity results of SGD solutions for a VGG-type network for the FashionMNIST dataset. We see that the predicted upper bound for the non-uniformity is both valid and quite sharp.

6 Qualitative properties of adaptive gradient algorithms

Adaptive gradient algorithms are a family of optimization algorithms widely used for training neural network models. These algorithms use a coordinate-wise scaling of the update direction (gradient or gradient with momentum) according to the history of the gradients. Two most popular adaptive gradient algorithms are RMSprop and Adam, whose update rules are

In (47) and (48), $\epsilon$ is a small constant used to avoid division by , usually taken to be $10^{-8}$. It is added to each component of the vector in the denominators of these equations. The division should also be understood as being component-wise. Here we will focus on Adam. For further discussion, we refer to .

One important tool for understanding these adaptive gradient algorithms is their continuous limit. Different continuous limits can be obtained from different ways of taking limits. If we let $\eta$ tends to while keeping $\alpha$ and $\beta$ fixed, the limiting dynamics will be

This becomes signGD when $\epsilon=0$. On the other hand, if we let $\alpha=1-a\eta$ and $\beta=1-b\eta$ and take $\eta\rightarrow 0$, then we obtain the limiting dynamics

In practice the loss curves of Adam can be very complicated. The left panel of Figure 11 shows an example. Three obvious features can be observed from these curves:

Fast initial convergence: the loss curve decreases very fast, sometimes even super-linearly, at the early stage of the training.

Small oscillations: The fast initial convergence is followed by oscillations around the minimum.

Large spikes: spikes are sudden increase of the value of the loss. They are followed by an oscillating recovery. Different from small oscillations, spikes make the loss much larger and the interval between two spikes is longer.

The fast initial convergence can be partly explained by the convergence property of signGD, which attains global minimum in finite time for strongly convex objective functions. Specifically, we have the following proposition .

Assume that the objective function satisfies the Polyak-Lojasiewicz (PL) condition: $\|\nabla f({\bm{x}})\|_{2}^{2}\geq\mu f({\bm{x}})$, for some positive constant $\mu$. Define continuous signGD dynamics,

The small oscillations and spikes may potentially be explained by linearization around the stationary point, though in this case, linearization is quite tricky due to the singular nature of the stationary point .

The performance of Adam depends sensitively on the values of $\alpha$ and $\beta$. This has also been studied in . Recall that $\alpha=1-a\eta$ and $\beta=1-b\eta$. Focusing on the region where $a$ and $b$ are not too large, three regimes with different behavior patterns were observed in the hyper-parameter space of $(a,b)$:

The spike regime: This happens when $b$ is sufficiently larger than $a$. In this regime large spikes appear in the loss curve, which makes the optimization process unstable.

The oscillation regime: This happens when $a$ and $b$ have similar magnitude (or in the same order). In this regime the loss curve exhibits fast and small oscillations. Small loss and stable loss curve can be achieved.

The divergence regime: This happens when $a$ is sufficiently larger than $b$. In this regime the loss curve is unstable and usually diverges after some period of training. This regime should be avoided in practice since the training loss stays large.

In Figure 10 we show one typical loss curve for each regime for a typical neural network model.

In Figure 11, we study Adam on a neural network model and show the final average training loss for different values of $a$ and $b$. One can see that Adam can achieve very small loss in the oscillation regime. The algorithm is not stable in the spike regime and may blow up in the divergence regime. These observations suggest that in practice one should take $a\approx b$ with small values of $a$ and $b$. Experiments on more complicated models, such as ResNet18, also support these conclusions .

7 Exploding and vanishing gradients for multi-layer neural networks

Exploding and vanishing gradients is one of the main obstacles for training of multi (many)-layer neural networks. Intuitively it is easy to see why this might be an issue. The gradient of the loss function with respect to the parameters involve a product of many matrices. Such a product can easily grow or diminish very fast as the number of products, namely the layers, increases.

Consider the multi-layer neural network with depth $L$:

To make quantitative statements, we consider the case when the weights are i.i.d. random variables. This is typically the case when the neural networks are initialized. This problem has been studied in . Below is a brief summary of the results obtained in these papers.

$\mu^{(l)}$, $\nu^{(l)}$ are symmetric around 0 for every $1\leq l\leq L$;

the variance of $\mu^{(l)}$ is $2/n_{l-1}$;

We consider the random network obtained by:

for any $i=1,2,\cdots,m_{l}$ and $j=1,2,\cdots,m_{l-1}$, i.e. the weights and biases at layer $l$ are drawn independently from $\mu^{(l)}$, $\nu^{(l)}$ respectively. Let

In contrast, the fourth moment of $Z_{p,q}$ is exponential in $\sum_{l}\frac{1}{m_{l}}$:

where $C_{\mu}>0$ is a constant only related to the (fourth) moment of $\mu=\{\mu^{(l)}\}_{l=1}^{L}$.

where $C_{\mu,K}>0$ is a constant depending on $K$ and the (first $2K$) moments of $\mu$.

Obviously, by Theorem 42, to avoid the exploding and vanishing gradient problem, we want the quantity $\sum_{l=1}^{L-1}\frac{1}{m_{l}}$ to be small. This is also borne out from numerical experiments .

8 What’s not known?

There are a lot that we don’t know about training dynamics. Perhaps the most elegant mathematical question is the global convergence of the mean-field equation for two-layer neural networks.

if $\rho_{0}$ is a smooth distribution with full support, then the dynamics described by this flow should converge and the population risk $\mathcal{R}$ should converge to 0;

if the target function $f^{*}$ lies in the Barron space, then the Barron norm of the output function stays uniformly bounded.

Similar statements should also hold for the empirical risk.

Another core question is when two-layer neural networks can be trained efficiently The authors would like to thank Jason Lee for helpful conversations on the topic.. The work shows that there exist target functions with Barron norms of size $poly(d)$, such that the training is exponentially slow in the statistical-query (SQ) setting . A concrete example is $f^{*}({\bm{x}})=\sin(dx_{1})$ with ${\bm{x}}\sim\mathcal{N}(0,I_{d})$, for which it is easy to verify that the Barron norm of $f^{*}$ is $poly(d)$. However, shows that the gradients of the corresponding neural networks are exponentially small, i.e. $O(e^{-d})$. Thus even for a moderate $d$, it is impossible to evaluate the gradients accurately on a finite-precision machine due to the floating-point error. Hence gradient-based optimizers are unlikely to succeed. This suggests that the Barron space is very likely too large for studying the training of two-layer neural networks. It is an open problem to identify the right function space, such that the functions can be learned in polynomial time by two-layer neural networks.

There are (at least) two possibilities for why the training becomes slow for certain (Barron) target functions in high dimension:

The training is slow in the continuous model due to the large parameter space or

The training is fast in the continuous model with dimension-independent rates, but the discretization becomes more difficult in high dimension.

It is unknown which of these explanations applies. Under strong conditions, it is known that parameters which are initialized according to a law $\pi_{0}$ and trained by gradient descent perform better at time $t>0$ than parameters which are drawn from the distribution $\pi_{t}$ given by the Wasserstein gradient flow starting at $\pi_{0}$ . This might suggest that also gradient flows with continuous initial condition may not reduce risk at a dimension-independent rate.

One can ask similar questions for multi-layer neural networks and residual neural networks. However, it is much more natural to formulate them using the continuous formulation. We will postpone this to a separate article.

Even less is known for the training dynamics of neural network models under the conventional scaling, except for the high over-parametrized regime. This issue is all the more important since conventional scaling is the overwhelming scaling used in practice.

In particular, since the neural networks used in practice are often over-parametrized, and they seem to perform much better than random feature models, some form of implicit regularization must be at work. Identifying the presence and the mechanism of such implicit regularization is a very important question for understanding the success of neural network models.

We have restricted our discussion to gradient descent training. What about stochastic gradient descent and other training algorithms?

Concluding remarks

Very briefly, let us summarize the main theoretical results that have been established so far.

Approximation/generalization properties of hypothesis space:

Function spaces and quantitative measures for the approximation properties for various machine learning models. The random feature model is naturally associated with the corresponding RKHS. In the same way, the Barron norm is identified as the natural measure associated with two-layer neural network models (note that this is different from the spectral norm defined by Barron). For residual networks, the corresponding quantities are defined for the flow-induced spaces. For multi-layer networks, a promising candidate is provided by the multi-layer norm.

Generalization error estimates of regularized model. Dimension-independent error rates have been established for these models. Except for multi-layer neural networks, these error rates are comparable to Monte Carlo. They provide a way to compare these different machine learning models and serve as a benchmark for studying implicit regularization.

Training dynamics for highly over-parametrized neural network models:

Exponential convergence for the empirical risk.

Their generalization properties are no better than the corresponding random feature model or kernel method.

Mean-field training dynamics for two-layer neural networks:

If the initial distribution has full support and the GD path converges, then it must converge to a global minimum.

A lot has also been learned from careful numerical experiments and partial analytical arguments, such as:

Over-parametrized networks may be able to interpolate any training data.

The “double descent” and “slow deterioration” phenomenon for the random feature model and their effect on the corresponding neural network model.

The qualitative behavior of adaptive optimization algorithms.

The global minima selection mechanism for different optimization algorithms.

The phase transition of the generalization properties of two-layer neural networks under the conventional scaling.

We have mentioned many open problems throughout this article. Besides the rigorous mathematical results that are called for, we feel that carefully designed numerical experiments should also be encouraged. In particular, they might give some insight on the difference between neural network models of different depth (for example, two-layer and three-layer neural neworks), and the difference between scaled and unscaled residual network models.

One very important issue that we did not discuss much is the continuous formulation of machine learning. We feel this issue deserves a separate article when the time is ripe.

Acknowledgement:. This work is supported in part by a gift to the Princeton University from iFlytek.

References

Appendix A Proofs for Section 3.1

Since $\text{supp }(\rho_{m})\subseteq K:=[-C,C]\times\Omega$, the sequence of probability measures $(\rho_{m})$ is tight. By Prokhorov’s theorem, there exists a subsequence $(\rho_{m_{k}})$ and a probability measure $\rho^{*}\in\mathcal{P}(K)$ such that $\rho_{m_{k}}$ converges weakly to $\rho^{*}$. Due to that $g(a,\bm{w};{\bm{x}})=a\phi({\bm{x}};\bm{w})$ is bounded and continuous with respect to $(a,\bm{w})$ for any ${\bm{x}}\in^{d}$, we have

Denote by $a^{*}(\bm{w})=\int a\rho^{*}(a|\bm{w})da$ the conditional expectation of $a$ given $\bm{w}$. Then we have

By the strong law of large numbers, the last equality holds with probability $1$. Taking $g(\bm{w})=a^{*}(\bm{w})\phi({\bm{x}};\bm{w})$, we obtain

To prove Theorem 7, we first need the following result ([33, Proposition 7]).

Proof of Theorem 7

By inserting Eqn. (61), we complete the proof. ∎

Appendix B Proofs for Section 3.2

where $\xi_{i}=\pm 1$ with probability $1/2$ independently of $\xi_{j}$ are Rademacher variables. Now we bound

by using the Rademacher complexity of the unit ball in Hilbert spaces [76, Lemma 26.10]. Applying the same argument to $-\left[\frac{1}{m}\sum_{i=1}^{m}a_{i}\sigma(\bm{w}_{i}^{T}{\bm{x}})-f({\bm{x}})\right]$, we find that

In particular, there exists weights $(a_{i},\bm{w}_{i})_{i=1}^{m}$ such that the inequality is true.