Fit without fear: remarkable mathematical phenomena of deep learning through the prism of interpolation

Preface

In recent years we have witnessed triumphs of Machine Learning in practical challenges from machine translation to playing chess to protein folding. These successes rely on advances in designing and training complex neural network architectures and on availability of extensive datasets. Yet, while it is easy to be optimistic about the potential of deep learning for our technology and science, we may still underestimate the power of fundamental mathematical and scientific principles that can be learned from its empirical successes.

In what follows, I will attempt to assemble some pieces of the remarkable mathematical mosaic that is starting to emerge from the practice of deep learning. This is an effort to capture parts of an evolving and still elusive picture with many of the key pieces still missing. The discussion will be largely informal, aiming to build mathematical concepts and intuitions around empirically observed phenomena. Given the fluid state of the subject and our incomplete understanding, it is necessarily a subjective, somewhat impressionistic and, to a degree, conjectural view, reflecting my understanding and perspective. It should not be taken as a definitive description of the subject as it stands now. Instead, it is written with the aspiration of informing and intriguing a mathematically minded reader and encouraging deeper and more detailed research.

Introduction

In the last decade theoretical machine learning faced a crisis. Deep learning, based on training complex neural architectures, has become state-of-the-art for many practical problems, from computer vision to playing the game of Go to Natural Language Processing and even for basic scientific problems, such as, recently, predicting protein folding . Yet, the mathematical theory of statistical learning extensively developed in the 1990’s and 2000’s struggled to provide a convincing explanation for its successes, let alone help in designing new algorithms or providing guidance in improving neural architectures. This disconnect resulted in significant tensions between theory and practice. The practice of machine learning was compared to “alchemy”, a pre-scientific pursuit, proceeding by pure practical intuition and lacking firm foundations . On the other hand, a counter-charge of practical irrelevance, “looking for lost keys under a lamp post, because that’s where the light is” was leveled against the mathematical theory of learning.

In what follows, I will start by outlining some of the reasons why classical theory failed to account for the practice of “modern” machine learning. I will proceed to discuss an emerging mathematical understanding of the observed phenomena, an understanding which points toward a reconciliation between theory and practice.

The key themes of this discussion are based on the notions of interpolation and over-parameterization, and the idea of a separation between the two regimes:

The classical setting can be characterized by limited model complexity, which does not allow arbitrary data to be fit exactly. The goal is to understand the properties of the (typically unique) classifier with the smallest loss. The standard tools include Uniform Laws of Large Numbers resulting in “what you see is what you get” (WYSIWYG) bounds, where the fit of classifiers on the training data is predictive of their generalization to unseen data. Non-convex optimization problems encountered in this setting typically have multiple isolated local minima, and the optimization landscape is locally convex around each minimum.

Over-parameterized setting deals with rich model classes, where there are generically manifolds of potential interpolating predictors that fit the data exactly. As we will discuss, some but not all of those predictors exhibit strong generalization to unseen data. Thus, the statistical question is understanding the nature of the inductive bias – the properties that make some solutions preferable to others despite all of them fitting the training data equally well. In interpolating regimes, non-linear optimization problems generically have manifolds of global minima. Optimization is always non-convex, even locally, yet it can often be shown to satisfy the so-called Polyak - Łojasiewicz (PL) condition guaranteeing convergence of gradient-based optimization methods.

As we will see, interpolation, the idea of fitting the training data exactly, and its sibling over-parameterization, having sufficiently many parameters to satisfy the constraints corresponding to fitting the data, taken together provide a perspective on some of the more surprising aspects of neural networks and other inferential problems. It is interesting to point out that interpolating noisy data is a deeply uncomfortable and counter-intuitive concept to statistics, both theoretical and applied, as it is traditionally concerned with over-fitting the data. For example, in a book on non-parametric statistics (page 21) the authors dismiss a certain procedure on the grounds that it “may lead to a function which interpolates the data and hence is not a reasonable estimate”. Similarly, a popular reference (page 194) suggests that “a model with zero training error is overfit to the training data and will typically generalize poorly”.

Likewise, over-parameterization is alien to optimization theory, which is traditionally more interested in convex problems with unique solutions or non-convex problems with locally unique solutions. In contrast, as we discuss in Section 4, over-parameterized optimization problems are in essence never convex nor have unique solutions, even locally. Instead, the solution chosen by the algorithm depends on the specifics of the optimization process.

To avoid confusion, it is important to emphasize that interpolation is not necessary for good generalization. In certain models (e.g., ), introducing some regularization is provably preferable to fitting the data exactly. In practice, early stopping is typically used for training neural networks. It prevents the optimization process from full convergence and acts as a type of regularization . What is remarkable is that interpolating predictors often provide strong generalization performance, comparable to the best possible predictors. Furthermore, the best practice of modern deep learning is arguably much closer to interpolation than to the classical regimes (when training and testing losses match). For example in his 2017 tutorial on deep learning Ruslan Salakhutdinov stated that “The best way to solve the problem from practical standpoint is you build a very big system $\ldots$ basically you want to make sure you hit the zero training error”. While more tuning is typically needed for best performance, these “overfitted” systems already work well . Indeed, it appears that the largest technologically feasible networks are consistently preferable for best performance. For example, in 2016 the largest neural networks had fewer than $10^{9}$ trainable parameters , the current (2021) state-of-the-art Switch Transformers have over $10^{12}$ weights, over three orders of magnitude growth in under five years!

Just as a literal physical prism separates colors mixed within a ray of light, the figurative prism of interpolation helps to disentangle a blend of properties within the complex picture of modern Machine Learning. While significant parts are still hazy or missing and precise analyses are only being developed, many important pieces are starting to fall in place.

The problem of generalization

Here $l(f({\mathbf{x}}),y)={\bf 1}_{f({\mathbf{x}})\neq y}$ is the Kronecker delta function called $0-1$ loss function. The expected loss of the Bayes optimal classifier $f^{*}$ it called the Bayes loss or Bayes risk.

We note that $0-1$ loss function can be problematic due to its discontinuous nature, and is entirely unsuitable for regression, where the square loss $l(f({\mathbf{x}}),y)=(f({\mathbf{x}})-y)^{2}$ is typically used. For the square loss, the optimal predictor $f^{*}$ is called the regression function.

In what follows, we will simply denote a general loss by $l(f({\mathbf{x}}),y)$, specifying its exact form when needed.

2 The framework of empirical and structural risk Minimization

Even in that formulation the problem is still under-defined as infinitely many different functions minimize the empirical risk. Yet, it can be made well-posed by restricting the space of candidate functions ${\mathcal{H}}$ to make the solution unique. Thus, we obtain the following formulation of the Empirical Risk Minimization (ERM):

Solving this optimization problem is called “training”. Of course, $f_{\rm emp}$ is only useful to the degree it approximates $f^{*}$. While superficially the predictors $f^{*}$ and $f_{\rm emp}$ appear to be defined similarly, their mathematical relationship is subtle due, in particular, to the choice of the space ${\mathcal{H}}$, the “structural part” of the empirical risk minimization.

According to the discussion in , “the theory of induction” based on the Structural Risk Minimization must meet two mathematical requirements:

The theory of induction is based on the Uniform Law of Large Numbers.

Effective methods of inference must include Capacity Control.

A uniform law of large numbers (ULLN) indicates that for any hypothesis in ${\mathcal{H}}$, the loss on the training data is predictive of the expected (future) loss:

We generally expect that ${\mathcal{R}}(f)\geq{\mathcal{R}}_{\rm emp}(f)$, which allows ULNN to be written as a one-sided inequality, typically of the formThis is the most representative bound, rates faster and slower than $\sqrt{n}$ are also found in the literature. The exact dependence on $n$ does not change our discussion here.

Here ${\mathop{cap}({\mathcal{H}})}$ is a measure of the capacity of the space ${\mathcal{H}}$, such as its Vapnik-Chervonenkis (VC) dimension or the covering number (see ), and $O^{*}$ can contain logarithmic terms and other terms of lower order. The inequality above holds with high probability over the choice of the data sample.

Eq. 2 is a mathematical instantiation of the ULLN condition and directly implies

This guarantees that the true risk of $f_{\rm emp}$ is nearly optimal for any function in ${\mathcal{H}}$, as long as ${\mathop{cap}({\mathcal{H}})}\ll n$.

The structural condition CC is needed to ensure that ${\mathcal{H}}$ also contains functions that approximate $f^{*}$. Combining CC and ULLN and applying the triangle inequality, yields a guarantee that ${\mathcal{R}}_{\rm emp}(f_{\rm emp})$ approximates ${\mathcal{R}}(f^{*})$ and the goal of generalization is achieved.

It is important to point out that the properties ULLN and CC are in tension to each other. If the class ${\mathcal{H}}$ is too small, no $f\in{\mathcal{H}}$ will generally be able to adequately approximate $f^{*}$. In contrast, if ${\mathcal{H}}$ is too large, so that ${\mathop{cap}({\mathcal{H}})}$ is comparable to $n$, the capacity term is large and there is no guarantee that ${\mathcal{R}}_{\rm emp}(f_{\rm emp})$ will be close to the expected risk ${\mathcal{R}}(f_{\rm emp})$. In that case the bound becomes tautological (such as the trivial bound that the classification risk is bounded by $1$ from above).

Hence the prescriptive aspect of Structural Risk Minimization according to Vapnik is to enlarge ${\mathcal{H}}$ until we find the sweet spot, a point where the empirical risk and the capacity term are balanced. This is represented by Fig. 1 (cf. , Fig. 6.2).

This view, closely related to the “bias-variance dilemma” in statistics , had become the dominant paradigm in supervised machine learning, encouraging a rich and increasingly sophisticated line of mathematical research uniform laws of large numbers and concentration inequalities.

3 Margins theory and data-dependent explanations.

Yet, even in the 1990’s it had become clear that successes of Adaboost and neural networks were difficult to explain from the SRM or bias-variance trade-off paradigms. Leo Breiman, a prominent statistician, in his note from 1995 posed the question “Why don’t heavily parameterized neural networks overfit the data?”. In particular, it was observed that increasing complexity of classifiers (capacity of ${\mathcal{H}}$) in boosting did not necessarily lead to the expected drop of performance due to over-fitting. Why did the powerful mathematical formalism of uniform laws of large numbers fail to explain the observed evidenceThis question appears as a refrain throughout the history of Machine Learning and, perhaps, other domains.?

An elegant explanation known as the margins theory, was proposed in . It is based on a more careful examination of the bound in Eq. 2, which identifies a serious underlying issue. We observe that the bound applies to any function $f\in{\mathcal{H}}$. Yet, in the learning context, we are not at all concerned with all functions, only with those that are plausible predictors. Indeed, it is a priori clear that the vast majority of predictors in standard function classes (linear functions, for example), are terrible predictors with performance no better than chance. Whether their empirical risk matches the true risk may be of importance to the theory of empirical processes or to functional analysis, but is of little concern to a “theory of induction”. The plausible candidate functions, those that are in an appropriate sense close to $f^{*}$, form a much narrower subset of ${\mathcal{H}}$. Of course, “closeness” needs to be carefully defined to be empirically observable without the exact prior knowledge of $f^{*}$.

To give an important special case, suppose we believe that our data are separable, so that ${\mathcal{R}}(f^{*})=0$. We can then concentrate our analysis on the subset of the hypothesis set ${\mathcal{H}}$ with small empirical loss

Indeed, since ${\mathcal{R}}(f^{*})=0$, ${\mathcal{R}}_{\rm emp}(f^{*})=0$ and hence $f^{*}\in{\mathcal{H}}_{\epsilon}$.

The capacity $\mathop{cap}({\mathcal{H}}_{\epsilon})$ will generally be far smaller than ${\mathop{cap}({\mathcal{H}})}$ and we thus hope for a tighter bound. It is important to note that the capacity $\mathop{cap}({\mathcal{H}}_{\epsilon})$ is a data-dependent quantity as ${\mathcal{H}}_{\epsilon}$ is defined in terms of the training data. Thus we aim to replace Eq. 2 with a data-dependent bound:

where class capacity ${\mathop{cap}({\mathcal{H}},X)}$ depends both on the hypothesis class ${\mathcal{H}}$ and the training data ${\mathcal{X}}$.

This important insight underlies the margins theory , introduced specifically to address the apparent lack of over-fitting in boosting. The idea of data-dependent margin bounds has led to a line of increasingly sophisticated mathematical work on understanding data-dependent function space complexity with notions such as Rademacher Complexity . Yet, we note that as an explanation for the effectiveness of Adaboost, the margins theory had not been universally accepted (see, e.g., for an interesting discussion).

4 What you see is not what you get

It is important to note that the generalization bounds mentioned above, even the data-dependent bounds such as Eq. 3, are “what you see is what you get” (WYSIWYG): the empirical risk that you see in training approximates and bounds the true risk that you expect on unseen data, with the capacity term providing an upper bound on the difference between expected and empirical risk.

Yet, it had gradually become clear (e.g., ) that in modern ML, training risk and the true risk were often dramatically different and lacked any obvious connection. In an influential paper the authors demonstrate empirical evidence showing that neural networks trained to have zero classification risk in training do not suffer from significant over-fitting. The authors argue that these and similar observations are incompatible with the existing learning theory and “require rethinking generalization”. Yet, their argument does not fully rule out explanations based on data-dependent bounds such as those in which can produce nontrivial bounds for interpolating predictors if the true Bayes risk is also small.

A further empirical analysis in made such explanations implausible, if not outright impossible. The experiments used a popular class of algorithms known as kernel machines, which are mathematically predictors of the form

Here $K({\mathbf{x}},{\mathbf{z}})$ is a positive definite kernel function (see, e.g., for a review), such as the commonly used Gaussian kernel $K({\mathbf{x}},{\mathbf{z}})=e^{-\frac{\|{\mathbf{x}}-{\mathbf{z}}\|^{2}}{2}}$ or the Laplace kernel $K({\mathbf{x}},{\mathbf{z}})=e^{-\|{\mathbf{x}}-{\mathbf{z}}\|}$. It turns out that there is a unique predictor $f_{\rm ker}$ of that form which interpolates the data:

The coefficients $\alpha_{i}$ can be found analytically, by matrix inversion ${\boldsymbol{\alpha}}=K^{-1}{\mathbf{y}}$. Here $K$ is the kernel matrix $K_{ij}=K({\mathbf{x}}_{i},{\mathbf{x}}_{j})$, and ${\mathbf{y}}$ is the vector containing the labels $y_{i}$.

Consider now a probability distribution $P$, “corrupted” by label noise. Specifically (for a two-class problem) with probability $q$ the label for any ${\mathbf{x}}$ is assigned from $\{-1,1\}$ with equal probability, and with probability $1-q$ it is chosen according to the original distribution $P$. Note that $P_{q}$ can be easily constructed synthetically by randomizing the labels on the $q$ fraction of the training and test sets respectively.

It can be seen that the Bayes optimal classifier for the corrupted distribution $P_{q}$ coincides with the Bayes optimal $f^{*}_{P}$ for the original distribution:

Furthermore, it is easy to check that the $0-1$ loss of the Bayes optimal predictor $f^{*}_{P}$ computed with respect to $P_{q}$ (denoted by ${\mathcal{R}}_{P_{q}}$) is bounded from below by the noise level:

It was empirically shown in that interpolating kernel machines $f_{\rm ker,q}$ (see Eq. 4) with common Laplace and Gaussian kernels, trained to interpolate $q$-corrupted data, generalizes nearly optimally (approaches the Bayes risk) to the similarly corrupted test data. An example of that is shown inFor a ten-class problem in panel (b), which makes the point even stronger. For simplicity, we only discuss a two-class analysis here. Fig. 2. In particular, we see that the Laplace kernel tracks the optimal Bayes error very closely, even when as much as $80\%$ of the data are corrupted (i.e., $q=0.8$).

Why is it surprising from the WYISWYG bound point of view? For simplicity, suppose $P$ is deterministic (${\mathcal{R}}({f^{*}_{P}})=0$), which is essentially the case [FOOTNOTE MOVED] in Fig. 2, Panel (b). In that case (for a two-class problem), ${\mathcal{R}}_{P_{q}}({f^{*}_{P}})=\frac{q}{2}$.

On the other hand ${\mathcal{R}}_{\rm emp}(f_{\rm ker,q})=0$ and hence for the left-hand side in Eq. 3 we have

To explain good empirical performance of $f_{\rm ker,q}$, a bound like Eq. 3 needs to be both correct and nontrivial. Since the left hand side is at least $\frac{q}{2}$ and observing that ${\mathcal{R}}_{P_{q}}(f_{\rm ker,q})$ is upper bounded by the loss of a random guess, which is $1/2$ for a two-class problem, we must have

Note that such a bound would require the multiplicative coefficient in $O^{*}$ to be tight within a multiplicative factor $1/q$ (which is $1.25$ for $q=0.8$). No such general bounds are known. In fact, typical bounds include logarithmic factors and other multipliers making really tight estimates impossible. More conceptually, it is hard to see how such a bound can exist, as the capacity term would need to “magically” knowThis applies to the usual capacity definitions based on norms, covering numbers and similar mathematical objects. In principle, it may be possible to “cheat” by letting capacity depend on complex manipulations with the data, e.g., cross-validation. This requires a different type of analysis (see for some recent attempts) and raises the question of what may be considered a useful generalization bound. We leave that discussion for another time. about the level of noise $q$ in the probability distribution. Indeed, a strict mathematical proof of incompatibility of generalization with uniform bounds was recently given in under certain specific settings. The consequent work proved that no good bounds can exist for a broad range of models.

Thus we see that strong generalization performance of classifiers that interpolate noisy data is incompatible with WYSIWYG bounds, independently of the nature of the capacity term.

5 Giving up on WYSIWYG, keeping theoretical guarantees

So can we provide statistical guarantees for classifiers that interpolate noisy data?

Until very recently there had not been many. In fact, the only common interpolating algorithm with statistical guarantees for noisy data is the well-known 1-NN ruleIn the last two or three years there has been significant progress on interpolating guarantees for classical algorithms like linear regression and kernel methods (see the discussion and references below). However, traditionally analyses nearly always used regularization which precludes interpolation.. Below we will go over a sequence of three progressively more statistically powerful nearest neighbor-like interpolating predictors, starting with the classical 1-NN rule, and going to simplicial interpolation and then to general weighted nearest neighbor/Nadaraya-Watson schemes with singular kernels.

Given an input ${\mathbf{x}}$, ${\rm 1{\text{-}}NN}({\mathbf{x}})$ outputs the label for the closest (in Euclidean or another appropriate distance) training example.

While the 1-NN rule is among the simplest and most classical prediction rules both for classification and regression, it has several striking aspects which are not usually emphasized in standard treatments:

It is an interpolating classifier, i.e., ${\mathcal{R}}_{\rm emp}({\rm 1{\text{-}}NN})=0$.

Despite “over-fitting”, classical analysis in shows that the classification risk of ${\cal R}({\rm 1{\text{-}}NN})$ is (asymptotically as $n\to\infty$) bounded from above by $2{\cdot}{\mathcal{R}}(f^{*})$, where $f^{*}$ is the Bayes optimal classifier defined by Eq. 1.

Not surprisingly, given that it is an interpolating classifier, there no ERM-style analysis of 1-NN.

It seems plausible that the remarkable interpolating nature of 1-NN had been written off by the statistical learning community as an aberration due to its high excess riskRecall that the excess risk of a classifier $f$ is the difference between the risk of the classifier and the risk of the optimal predictor ${\mathcal{R}}(f)-{\mathcal{R}}(f^{*})$.. As we have seen, the risk of 1-NN can be a factor of two worse than the risk of the optimal classifier. The standard prescription for improving performance is to use k-NN, an average of $k$ nearest neighbors, which no longer interpolates. As $k$ increases (assuming $n$ is large enough), the excess risk decreases as does the difference between the empirical and expected risks. Thus, for large $k$ (but still much smaller than $n$) we have, seemingly in line with the standard ERM-type bounds,

It is perhaps ironic that an outlier feature of $1$-NN rule, shared with no other common methods in the classical statistics literature (except for the relatively unknown work ), may be one of the cues to understanding modern deep learning.

5.2 Geometry of simplicial interpolation and the blessing of dimensionality

Yet, a modification of 1-NN different from k-NN maintains its interpolating property while achieving near-optimal excess risk, at least in when the dimension is high. The algorithm is simplicial interpolation analyzed statistically in . Consider a triangulation of the data, ${\mathbf{x}}_{1},\ldots,{\mathbf{x}}_{n}$, that is a partition of the convex hull of the data into a set of $d$-dimensional simplices so that:

Vertices of each simplex are data points.

For any data point ${\mathbf{x}}_{i}$ and simplex $s$, ${\mathbf{x}}_{i}$ is either a vertex of $s$ or does not belong to $s$.

The exact choice of the triangulation turns out to be unimportant as long as the size of each simplex is small enough. This is guaranteed by, for example, the well-known Delaunay triangulation.

Given a multi-dimensional triangulation, we define $f_{\rm simp}(x)$, the simplicial interpolant, to be a function which is linear within each simplex and such that $f_{\rm simp}(x_{i})=y_{i}$. It is not hard to check that $f_{\rm simp}$ exists and is unique.

It is worth noting that in one dimension simplicial interpolation based on the Delaunay triangulation is equivalent to 1-NN for classification. Yet, when the dimension $d$ is high enough, simplicial interpolation is nearly optimal both for classification and regression. Specifically, it is was shown in Theorem 3.4 in (Theorem 3.4) that simplicial interpolation benefits from a blessing of dimensionality. For large $d$, the excess risk of $f_{\rm simp}$ decreases with dimension:

Analogous results hold for regression, where the excess risk is similarly the difference between the loss of a predictor and the loss of the (optimal) regression function. Furthermore, for classification, under additional conditions $\sqrt{d}$ can be replaced by $e^{d}$ in the denominator.

Why does this happen? How can an interpolating function be nearly optimal despite the fact that it fits noisy data and why does increasing dimension help?

Suppose also that the probability distribution is uniform on the simplex (the convex hull of ${\mathbf{x}}_{1},\ldots{\mathbf{x}}_{d+1}$) and the “correct” labels are identically $1$. As our training data, we are given $({\mathbf{x}}_{i},y_{i})$, where $y_{i}=1$, except for the one vertex, which is “corrupted by noise”, so that $y_{d+1}=-1$. It is easy to verify that

We see that $f_{\rm simp}$ coincides with $f^{*}\equiv 1$ in the simplex except for the set $s_{1/2}=\{{\mathbf{x}}:\sum_{i=1}^{d}{\mathbf{x}}_{i}\leq{1/2}\}$, which is equal to the simplex $\frac{1}{2}s_{d}$ and thus

We see that the interpolating predictor $f_{\rm simp}$ is different from the optimal, but the difference is highly localized around the “noisy” vertex, while at most points within $s_{d}$ their predictions coincide. This is illustrated geometrically in Fig. 3. The reasons for the blessing of dimensionality also become clear, as small neighborhoods in high dimension have smaller volume relative to the total space. Thus, there is more freedom and flexibility for the noisy points to be localized.

5.3 Optimality of k-NN with singular weighting schemes

While simplicial interpolation improves on ${\rm 1{\text{-}}NN}$ in terms of the excess loss, it is still not consistent. In high dimension $f_{\rm simp}$ is near $f^{*}$ but does not converge to $f^{*}$ as $n\to\infty$. Traditionally, consistency and rates of convergence have been a central object of statistical investigation. The first result in this direction is , which showed statistical consistency of a certain kernel regression scheme, closely related to Shepard’s inverse distance interpolation .

It turns out that a similar interpolation scheme based on weighted $k$-NN can be shown to be consistent for both regression and classification and indeed to be optimal in a certain statistical sense (see for convergence rates for regression and classification and the follow-up work for optimal rates for regression). The scheme can be viewed as a type of Nadaraya-Watson predictor. It can be described as follows. Let $K({\mathbf{x}},{\mathbf{z}})$ be a singular kernel, such as

with an appropriate choice of $\alpha$. Consider the weighted nearest neighbor predictor

Here the sum is taken over the $k$ nearest neighbors of ${\mathbf{x}}$, ${\mathbf{x}}_{(1)},\ldots,{\mathbf{x}}_{(k)}$. While the kernel $K({\mathbf{x}},{\mathbf{x}}_{(i)})$ is infinite at $x={\mathbf{x}}_{i}$, it is not hard to see that $f_{\rm sing}(x)$ involves a ratio that can be defined everywhere due to the cancellations between the singularities in the numerator and the denominator. It is, furthermore, a continuous function of ${\mathbf{x}}$. Note that for classification it suffices to simply take the sign of the numerator $\sum_{i=1}^{k}K({\mathbf{x}},{\mathbf{x}}_{(i)})y_{(i)}$ as the denominator is positive.

To better understand how such an unusual scheme can be consistent for regression, consider an example shown in Fig. 4 for one-dimensional data sampled from a noisy linear model: $y=x+\epsilon$, where $\epsilon$ is normally distributed noise. Since the predictor $f_{\rm sing}(x)$ fits the noisy data exactly, it is far from optimal on the majority of data points. Yet, the prediction is close to optimal for most points in the interval $$! In general, as$n\to\infty$, the fraction of those points tends to$1$.

We will discuss this phenomenon further in connection to adversarial examples in deep learning in Section 5.2.

6 Inductive biases and the Occam’s razor

The realization that, contrary to deeply ingrained statistical intuitions, fitting noisy training data exactly does not necessarily result in poor generalization, inevitably leads to quest for a new framework for a “theory of induction”, a paradigm not reliant on uniform laws of large numbers and not requiring empirical risk to approximate the true risk.

While, as we have seen, interpolating classifiers can be statistically near-optimal or optimal, the predictors discussed above appear to be different from those widely used in ML practice. Simplicial interpolation, weighted nearest neighbor or Nadaraya-Watson schemes do not require training and can be termed direct methods. In contrast, common practical algorithms from linear regression to kernel machines to neural networks are “inverse methods” based on optimization. These algorithms typically rely on algorithmic empirical risk minimization, where a loss function ${\mathcal{R}}_{\rm emp}(f_{\mathbf{w}})$ is minimized via a specific algorithm, such as stochastic gradient descent (SGD) on the weight vector ${\mathbf{w}}$. Note that there is a crucial and sometimes overlooked difference between the empirical risk minimization as an algorithmic process and the Vapnik’s ERM paradigm for generalization, which is algorithm-independent. This distinction becomes important in over-parameterized regimes, where the hypothesis space ${\mathcal{H}}$ is rich enough to fit any data setAssuming that ${\mathbf{x}}_{i}\neq{\mathbf{x}}_{j}$, when $i\neq j$. of cardinality $n$. The key insight is to separate “classical” under-parameterized regimes where there is typically no $f\in{\mathcal{H}}$, such that ${\mathcal{R}}(f)=0$ and “modern” over-parameterized settings where there is a (typically large) set ${\mathcal{S}}$ of predictors that interpolate the training data

First observe that an interpolating learning algorithm ${\mathcal{A}}$ selects a specific predictor $f_{\mathcal{A}}\in{\mathcal{S}}$. Thus we are faced with the issue of the inductive bias: why do solutions, such as those obtained by neural networks and kernel machines, generalize, while other possible solutions do notThe existence of non-generalizing solutions is immediately clear by considering over-parameterized linear predictors. Many linear functions fit the data – most of them generalize poorly.. Notice that this question cannot be answered through the training data alone, as any $f\in{\mathcal{S}}$ fits data equally wellWe note that inductive biases are present in any inverse problem. Interpolation simply isolates this issue.. While no conclusive recipe for selecting the optimal $f\in{\mathcal{S}}$ yet exists, it can be posited that an appropriate notion of functional smoothness plays a key role in that choice. As argued in , the idea of maximizing functional smoothness subject to interpolating the data represents a very pure form of the Occam’s razor (cf. ). Usually stated as

Entities should not be multiplied beyond necessity,

the Occam’s razor implies that the simplest explanation consistent with the evidence should be preferred. In this case fitting the data corresponds to consistency with evidence, while the smoothest function is “simplest”. To summarize, the “maximum smoothness” guiding principle can be formulated as:

Select the smoothest function, according to some notion of functional smoothness, among those that fit the data perfectly.

We note that kernel machines described above (see Eq. 4) fit this paradigm precisely. Indeed, for every positive definite kernel function $K({\mathbf{x}},{\mathbf{z}})$, there exists a Reproducing Kernel Hilbert Space ( functional spaces, closely related to Sobolev spaces, see ) ${\mathcal{H}}_{K}$, with norm $\|\cdot\|_{{\mathcal{H}}_{K}}$ such that

We proceed to discuss how this idea may apply to training more complex variably parameterized models including neural networks.

7 The Double Descent phenomenon

A hint toward a possible theory of induction is provided by the

double descent generalization curve (shown in Fig. 5), a pattern proposed in as a replacement for the classical U-shaped generalization curve (Fig. 1).

When the capacity of a hypothesis class ${\mathcal{H}}$ is below the interpolation threshold, not enough to fit arbitrary data, learned predictors follow the classical U-curve from Figure 1. The shape of the generalization curve undergoes a qualitative change when the capacity of ${\mathcal{H}}$ passes the interpolation threshold, i.e., becomes large enough to interpolate the data. Although predictors at the interpolation threshold typically have high risk, further increasing the number of parameters (capacity of ${\mathcal{H}}$) leads to improved generalization. The double descent pattern has been empirically demonstrated for a broad range of datasets and algorithms, including modern deep neural networks and observed earlier for linear models . The “modern” regime of the curve, the phenomenon that large number of parameters often do not lead to over-fitting has historically been observed in boosting and random forests, including interpolating random forests as well as in neural networks .

Why should predictors from richer classes perform better given that they all fit data equally well? Considering an inductive bias based on smoothness provides an explanation for this seemingly counter-intuitive phenomenon as larger spaces contain will generally contain “better” functions. Indeed, consider a hypothesis space ${\mathcal{H}}_{1}$ and a larger space ${\mathcal{H}}_{2},{\mathcal{H}}_{1}\subset{\mathcal{H}}_{2}$. The corresponding subspaces of interpolating predictors, ${\mathcal{S}}_{1}\subset{\mathcal{H}}_{1}$ and ${\mathcal{S}}_{2}\subset H_{2}$, are also related by inclusion: ${\mathcal{S}}_{1}\subset{\mathcal{S}}_{2}$. Thus, if $\|\cdot\|_{s}$ is a functional norm, or more generally, any functional, we see that

Assuming that $\|\cdot\|_{s}$ is the “right” inductive bias, measuring smoothness (e.g., a Sobolev norm), we expect the minimum norm predictor from ${\mathcal{H}}_{2}$, $f_{{\mathcal{H}}_{2}}=\arg\min_{f\in S_{2}}\|f\|_{s}$ to be superior to that from ${\mathcal{H}}_{1}$, $f_{{\mathcal{H}}_{1}}=\arg\min_{f\in S_{1}}\|f\|_{s}$.

A visual illustration for double descent and its connection to smoothness is provided in Fig. 6 within the random ReLU family of models in one dimension. A very similar Random Fourier Feature family is described in more mathematical detail below.The Random ReLU family consists of piecewise linear functions of the form $f({\mathbf{w}},x)=\sum_{k}w_{k}\min(v_{k}x+b_{k},0)$ where $v_{k},b_{k}$ are fixed random values. While it is quite similar to RFF, it produces better visualizations in one dimension. The left panel shows what may be considered a good fit for a model with a small number of parameters. The middle panel, with the number of parameters slightly larger than the minimum necessary to fit the data, shows textbook over-fitting. However increasing the number of parameters further results in a far more reasonably looking curve. While this curve is still piece-wise linear due to the nature of the model, it appears completely smooth. Increasing the number of parameters to infinity will indeed yield a differentiable function (a type of spline), although the difference between $3000$ and infinitely many parameters is not visually perceptible. As discussed above, over-fitting appears in a range of models around the interpolation threshold which are complex but yet not complex enough to allow smooth structure to emerge. Furthermore, low complexity parametric models and non-parametric (as the number of parameters approaches infinity) models coexist within the same family on different sides of the interpolation threshold.

Given data $\{{\mathbf{x}}_{i},y_{i}\},i=1,\ldots,n$, we can fit $f_{m}\in{\mathcal{H}}_{m}$ by linear regression on the coefficients ${\mathbf{w}}$. In the overparameterized regime linear regression is given by minimizing the norm under the interpolation constraintsAs opposed to the under-parameterized setting when linear regression it simply minimizes the empirical loss over the class of linear predictors.:

We see that increasing the number of parameters $m$ expands the space of interpolating classifiers in ${\mathcal{H}}_{m}$ and allows to obtain progressively better approximations of the ultimate functional smoothness minimizer $f_{\rm ker}$. Thus adding parameters in the over-parameterized setting leads to solutions with smaller norm, in contrast to under-parameterized classical world when more parameters imply norm increase. The norm of the weight vector $\|{\mathbf{w}}\|$ asymptotes to the true functional norm of the solution $f_{\rm ker}$ as $m\to\infty$. This is verified experimentally in Fig. 7. We see that the generalization curves for both 0-1 loss and the square loss follow the double descent curve with the peak at the interpolation threshold. The norm of the corresponding classifier increases monotonically up to the interpolation peak and decreases beyond that. It asymptotes to the norm of the kernel machine which can be computed using the following explicit formula for a function written in the form of Eq. 4) (where $K$ is the kernel matrix):

8 When do minimum norm predictors generalize?

As we have discussed above, considerations of smoothness and simplicity suggest that minimum norm solutions may have favorable generalization properties. This turns out to be true even when the norm does not have a clear interpretation as a smoothness functional. Indeed, consider an ostensibly simple classical regression setup, where data satisfy a linear relation corrupted by noise $\epsilon_{i}$

In the over-parameterized setting, when $d>n$, least square regression yields a minimum norm interpolator given by $y({\mathbf{x}})=\langle{\boldsymbol{\beta}}_{\text{int}},{\mathbf{x}}\rangle$, where

${\boldsymbol{\beta}}_{\text{int}}$ can be written explicitly as

where $\bf X$ is the data matrix, ${\mathbf{y}}$ is the vector of labels and ${\bf X}^{\dagger}$ is the Moore-Penrose (pseudo-)inverseIf ${\mathbf{X}}{\mathbf{X}}^{T}$ is invertible, as is usually the case in over-parameterized settings, ${\mathbf{X}}^{\dagger}={\mathbf{X}}^{T}({\mathbf{X}}{\mathbf{X}}^{T})^{-1}$. In contrast, if ${\mathbf{X}}^{T}{\mathbf{X}}$ is invertible (under the classical under-parameterized setting), ${\mathbf{X}}^{\dagger}=({\mathbf{X}}^{T}{\mathbf{X}})^{-1}{\mathbf{X}}^{T}$. Note that both ${\mathbf{X}}{\mathbf{X}}^{T}$ and ${\mathbf{X}}^{T}{\mathbf{X}}$ matrices cannot be invertible unless $X$ is a square matrix, which occurs at the interpolation threshold.. Linear regression for models of the type in Eq. 8 is no doubt the oldestOriginally introduced by Gauss and, possibly later, Legendre! See . and best studied family of statistical methods. Yet, strikingly, predictors such as those in Eq. 9, have historically been mostly overlooked, at least for noisy data. Indeed, a classical prescription is to regularize the predictor by, e.g., adding a “ridge” $\lambda I$ to obtain a non-interpolating predictor. The reluctance to overfit inhibited exploration of a range of settings where $y({\mathbf{x}})=\langle{\boldsymbol{\beta}}_{\text{int}},{\mathbf{x}}\rangle$ provided optimal or near-optimal predictions. Very recently, these “harmless interpolation” or “benign over-fitting” regimes have become a very active direction of research, a development inspired by efforts to understand deep learning. In particular, the work provided a spectral characterization of models exhibiting this behavior. In addition to the aforementioned papers, some of the first work toward understanding “benign overfitting” and double descent under various linear settings include . Importantly, they demonstrate that when the number of parameters varies, even for linear models over-parametrized predictors are sometimes preferable to any “classical” under-parameterized model.

Notably, even in cases when the norm clearly corresponds to measures of functional smoothness, such as the cases of RKHS or, closely related random feature maps, the analyses of interpolation for noisy data are subtle and have only recently started to appear, e.g., . For a far more detailed overview of the progress on interpolation in linear regression and kernel methods see the parallel Acta Numerica paper .

9 Alignment of generalization and optimization in linear and kernel models

While over-parameterized models have manifolds of interpolating solutions, minimum norm solutions, as we have discussed, have special properties which may be conducive to generalization. For over-parameterized linear and kernel modelsKernel models are linear from the optimization point of view as they can be viewed as a fixed feature map followed by a linear method. Thus we will not make a distinction in the optimization context. They are, however, non-linear functions of the input space. there is a beautiful alignment of optimization and minimum norm interpolation: gradient descent GD or Stochastic Gradient Descent (SGD) initialized at the origin can be guaranteed to converge to ${\boldsymbol{\beta}}_{\text{int}}$ defined in Eq. 9. To see why this is the case we make the following observations:

${\boldsymbol{\beta}}_{\text{int}}\in{\mathcal{T}}$, where ${\mathcal{T}}=\mathop{Span}{\{x_{1},\ldots,x_{n}\}}$ is the span of the training examples (or their feature embeddings in the kernel case). To see that, verify that if ${\boldsymbol{\beta}}_{\text{int}}\notin{\mathcal{T}}$, orthogonal projection of ${\boldsymbol{\beta}}_{\text{int}}$ onto ${\mathcal{T}}$ is an interpolating predictor with even smaller norm, a contradiction to the definition of ${\boldsymbol{\beta}}_{\text{int}}$.

The (affine) subspace of interpolating predictors ${\mathcal{S}}$ (Eq. 6) is orthogonal to ${\mathcal{T}}$ and hence $\{{\boldsymbol{\beta}}_{\text{int}}\}={\mathcal{S}}\cap{\mathcal{T}}$.

These two points together are in fact a version of the Representer theorem briefly discussed in Sec. 3.7.

Consider now gradient descent for linear regression initialized at within the span of training examples ${\boldsymbol{\beta}}_{0}\in{\mathcal{T}}$. Typically, we simply choose ${\boldsymbol{\beta}}_{0}=0$ as the origin has the notable property of belonging to the span of any vectors. It can be easily verified that the gradient of the loss function at any point is also in the span of the training examples and thus the whole optimization path lies within ${\mathcal{T}}$. As the gradient descent converges to a minimizer of the loss function, and ${\mathcal{T}}$ is a closed set, GD must converge to the minimum norm solution ${\boldsymbol{\beta}}_{\text{int}}$. Remarkably, in the over-parameterized settings convergence to ${\boldsymbol{\beta}}_{\text{int}}$ is true for SGD, even with a fixed learning rate (see Sec. 4.4). In contrast, under-parameterized SGD with a fixed learning rate does not converge at all.

10 Is deep learning kernel learning? Transition to linearity in wide neural networks.

But how do these ideas apply to deep neural networks? Why are complicated non-linear systems with large numbers of parameters able to generalize to unseen data?

It is important to recognize that generalization in large neural networks is a robust pattern that holds across multiple dimensions of architectures, optimization methods and datasetsWhile details such as selection of activation functions, initialization methods, connectivity patterns or many specific parameters of training (annealing schedules, momentum, batch normalization, dropout, the list goes on ad infinitum), matter for state-of-the-art performance, they are almost irrelevant if the goal is to simply obtain passable generalization. . As such, the ability of neural networks to generalize to unseen data reflects a fundamental interaction between the mathematical structures underlying neural function spaces, algorithms and the nature of our data. It can be likened to the gravitational force holding the Solar System, not a momentary alignment of the planets.

This point of view implies that understanding generalization in complex neural networks has to involve a general principle, relating them to more tractable mathematical objects. A prominent candidate for such an object are kernel machines and their corresponding Reproducing Kernel Hilbert Spaces. As we discussed above, Random Fourier Features-based networks, a rather specialized type of neural architectures, approximate Gaussian kernel machines. Perhaps general neural networks can also be tied to kernel machines? Strikingly, it turns out to be the case indeed, at least for some classes of neural networks.

The surprising and singular finding of is that for a range of infinitely wide neural network architectures with linear output layer, $\phi_{\mathbf{w}}({{\mathbf{x}}})$ is independent of ${\mathbf{w}}$ in a ball around a random “initialization” point ${\mathbf{w}}_{0}$. That can be shown to be equivalent to the linearity of $f({\mathbf{w}},{\mathbf{x}})$ in ${\mathbf{w}}$ (and hence transition to linearity in the limit of infinite width):

We see that the deviation from the linearity is bounded by the spectral norm of the Hessian:

A general (feed-forward) neural network with $L$ hidden layers and a linear output layer is a function defined recursively as:

The parameter vector ${\mathbf{w}}$ is obtained by concatenation of all weight vectors ${\mathbf{w}}=({\mathbf{w}}^{(1)},\ldots,{\mathbf{w}}^{(L)},{\mathbf{v}})$ and the activation functions $\phi_{l}$ are usually applied coordinate-wise. It turns out these, seemingly complex, non-linear systems exhibit transition to linearity under quite general conditions (see ), given appropriate random initialization ${\mathbf{w}}_{0}$. Specifically, it can be shown that for a ball $\cal B$ of fixed radius around the initialization ${\mathbf{w}}_{0}$ the spectral norm of the Hessian satisfies

It is important to emphasize that linearity is a true emerging property of large systems and does not come from the scaling of the function value with the increasing width $m$. Indeed, for any $m$ the value of the function at initialization and its gradient are all of order $1$, $f({\mathbf{w}},x)=\Omega(1)$, $\nabla f({\mathbf{w}},x)=\Omega(1)$.

For simplicity, assume that $v_{i}\in\{-1,1\}$ are fixed and $w_{i}$ are trainable parameters. It is easy to see that in this case the Hessian $H({\mathbf{w}})$ is a diagonal matrix with

Thus, we see that the structure of the Hessian matrix forces its spectral norm to be a factor of $\sqrt{m}$ smaller compared to the gradient. If (following a common practice) $w_{i}$ are sampled iid from the standard normal distribution

If, furthermore, the second layer weights $v_{i}$ are sampled with expected value zero, $f({\mathbf{w}},x)=O(1)$. Note that to ensure the transition to linearity we need for the scaling in Eq. 15 to hold in ball of radius $O(1)$ around ${\mathbf{w}}$ (rather than just at the point ${\mathbf{w}}$), which, in this case, can be easily verified.

The example above illustrates how the transition to linearity is the result of the structural properties of the network (in this case the Hessian is a diagonal matrix) and the difference between the $2$-norm ind $\infty$-norm in a high-dimensional space. For general deep networks the Hessian is no longer diagonal, and the argument is more involved, yet there is a similar structural difference between the gradient and the Hessian related to different scaling of the $2$ and $\infty$ norms with dimension.

Furthermore, transition to linearity is not simply a property of large systems. Indeed, adding a non-linearity at the output layer, i.e., defining

where $f({\mathbf{w}},x)$ is defined by Eq. 13 and $\phi$ is any smooth function with non-zero second derivative breaks the transition to linearity independently of the width $m$ and the function $\phi$. To see that, observe that the Hessian of $g$, $H_{g}$ can be written, in terms of the gradient and Hessian of $f$, ($\nabla_{\mathbf{w}}f$ and $H({\mathbf{w}})$, respectively) as

We see that the second term in Eq. 16 is of the order $\|\nabla_{\mathbf{w}}f\|^{2}=\Omega(1)$ and does not scale with $m$. Thus the transition to linearity does not occur and the tangent kernel does not become constant in a ball of a fixed radius even as the width of the network tends to infinity. Interestingly, introducing even a single narrow “bottleneck” layer has the same effect even if the activation functions in that layer are linear (as long as some activation functions in at least one of the deeper layers are non-linear).

As we will discuss later in Section 4, the transition to linearity is not needed for optimization, which makes this phenomenon even more intriguing. Indeed, it is possible to imagine a world where the transition to linearity phenomenon does not exist, yet neural networks can still be optimized using the usual gradient-based methods.

It is thus even more fascinating that a large class of very complex functions turn out to be linear in parameters and the corresponding complex learning algorithms are simply training kernel machines. In my view this adds significantly to the evidence that understanding kernel learning is a key to deep learning as we argued in . Some important caveats are in order. While it is arguable that deep learning may be equivalent to kernel learning in some interesting and practical regimes, the jury is still out on the question of whether this point of view can provide a conclusive understanding of generalization in neural networks. Indeed a considerable amount of recent theoretical work has been aimed at trying to understand regimes (sometimes called the “rich regimes”, e.g., ) where the transition to linearity does not happen and the system is non-linear throughout the training process. Other work (going back to ) argues that there are theoretical barriers separating function classes learnable by neural networks and kernel machines . Whether these analyses are relevant for explaining empirically observed behaviours of deep networks still requires further exploration.

Please also see some discussion of these issues in Section 6.2.

The wonders of optimization

The success of deep learning has heavily relied on the remarkable effectiveness of gradient-based optimization methods, such as stochastic gradient descent (SGD), applied to large non-linear neural networks. Classically, finding global minima in non-convex problems, such as these, has been considered intractable and yet, in practice, neural networks can be reliably trained.

Over-parameterization and interpolation provide a distinct perspective on optimization. Under-parameterized problems are typically locally convex around their local minima. In contrast, over-parameterized non-linear optimization landscapes are generically non-convex, even locally. Instead, as we will argue, throughout most (but not all) of the parameter space they satisfy the Polyak - Łojasiewicz condition, which guarantees both existence of global minima within any sufficiently large ball and convergence of gradient methods, including GD and SGD.

Finally, as we discuss in Sec. 4.4, interpolation sheds light on a separate empirically observed phenomenon, the striking effectiveness of mini-batch SGD (ubiquitous in applications) in comparison to the standard gradient descent.

Mathematically, interpolation corresponds to identifying ${\mathbf{w}}$ so that

This is a system of $n$ equations with $M$ variables. Aggregating these equations into a single map,

and setting ${\mathbf{y}}=(y_{1},\ldots,y_{n})$, we can write that ${\mathbf{w}}$ is a solution for a single equation

When can such a system be solved? The question posed in such generality initially appears to be absurd. A special case, that of solving systems of polynomial equations, is at the core of algebraic geometry, a deep and intricate mathematical field. And yet, we can often easily train non-linear neural networks to fit arbitrary data . Furthermore, practical neural networks are typically trained using simple first order gradient-based methods, such as stochastic gradient descent (SGD).

The idea of over-parameterization has recently emerged as an explanation for this phenomenon based on the intuition that a system with more variables than equations can generically be solved. We first observe that solving Eq. 18 (assuming a solution exists) is equivalent to minimizing the loss function

This is a non-linear least squares problem, which is well-studied under classical under-parameterized settings (see , Chapter 10). What property of the over-parameterized optimization landscape allows for effective optimization by gradient descent (GD) or its variants? It is instructive to consider a simple example in Fig. 8 (from ). The left panel corresponds to the classical regime with many isolated local minima. We see that for such a landscape there is little hope that a local method, such as GD can reach a global optimum. Instead we expect it to converge to a local minimum close to the initialization point. Note that in a neighborhood of a local minimizer the function is convex and classical convergence analyses apply.

A key insight is that landscapes of over-parameterized systems look very differently, like the right panel in Fig 8(b). We see that there every local minimum is global and the manifold of minimizers ${\mathcal{S}}$ has positive dimension. It is important to observe that such a landscape is incompatible with convexity even locally. Indeed, consider an arbitrary point $s\in{\mathcal{S}}$ inside the insert in Fig 8(b).

Thus, one of the key lessons of deep learning in optimization: Convexity, even locally, cannot be the basis of analysis for over-parameterized systems.

But what mathematical property encapsulates ability to optimize by gradient descent for landscapes, such as in Fig. 8. It turns out that a simple condition proposed in 1963 by Polyak is sufficient for efficient minimization by gradient descent. This PL-condition (for Polyak and also Łojasiewicz, who independently analyzed a more general version of the condition in a different context ) is a simple first order inequality applicable to a broad range of optimization problems .

We say that ${\mathcal{L}}({\mathbf{w}})$ is $\mu$-PL, if the following holds:

Here ${\mathbf{w}}^{*}$ is a global minimizer and $\mu>0$ is a fixed real number. The original Polyak’s work showed that PL condition within a sufficiently large ball (with radius $O(1/\mu)$) implied convergence of gradient descent.

It is important to notice that, unlike convexity, PL-condition is compatible with curved manifolds of minimizers. However, in this formulation, the condition is non-local. While convexity can be verified point-wise by making sure that the Hessian of ${\mathcal{L}}$ is positive semi-definite, the PL condition requires ”oracle” knowledge of ${\mathcal{L}}({\mathbf{w}}^{*})$. This lack of point-wise verifiability is perhaps the reason PL-condition has not been used more widely in the optimization literature.

However simply removing the $L({\mathbf{w}}^{*})$ from Eq. 19 addresses this issue in over-parameterized settings! Consider the following modification called PL* in and local PL in .

It turns out that PL* condition in a ball of sufficiently large radius implies both existence of an interpolating solution within that ball and exponential convergence of gradient descent and, indeed, stochastic gradient descent.

It is interesting to note that PL* is not a useful concept in under-parameterized settings – generically, there is no solution to $F({\mathbf{w}})={\mathbf{y}}$ and thus the condition cannot be satisfied along the whole optimization path. On the other hand, the condition is remarkably flexible – it naturally extends to Riemannian manifolds (we only need the gradient to be defined) and is invariant under non-degenerate coordinate transformations.

2 Condition numbers of nonlinear systems

Why do over-parameterized systems satisfy the PL* condition? The reason is closely related to the Tangent Kernel discussed in Section 3.10. Consider the tangent kernel of the map $F({\mathbf{w}})$ defined as $n\times n$ matrix valued function

where $DF$ is the differential of the map $F$. It can be shown for the square loss ${\mathcal{L}}({\mathbf{w}})$ satisfies the PL*- condition with $\mu=\lambda_{\text{m}in}(K)$. Note that the rank of $K$ is less or equal to $M$. Hence, if the system is under-parameterized, i.e., $M<n$, $\lambda_{\text{m}in}(K)({\mathbf{w}})\equiv 0$ and the corresponding PL* condition is always trivial.

But how can an analytic condition, like a lower bound on the smallest eigenvalue of the tangent kernel, be verified for models such as neural networks?

3 Controlling PL* condition of neural networks

As discussed above and graphically illustrated in Fig. 9, we expect over-parameterized systems to satisfy the PL* condition over most of the parameter space. Yet, explicitly controlling $\mu=\lambda_{\text{m}in}(K)$ in a ball of a certain radius can be subtle. We can identify two techniques which help establish such control for neural networks and other systems. The first one, the Hessian control, uses the fact that near-linear systems are well-conditioned in a ball, provided they are well-conditioned at the origin. The second, transformation control, is based on the observation that well-conditioned systems stay such under composition with “benign” transformations. Combining these techniques can be used to prove convergence of randomly initialized wide neural networks.

Transition to linearity, discussed in Section 3.10, provides a powerful (if somewhat crude) tool for controlling $\lambda_{\text{m}in}(K)$ for wide networks. The key observation is that $K({\mathbf{w}})$ is closely related to the first derivative of $F$ at ${\mathbf{w}}$. Thus the change of $K({\mathbf{w}})$ from the initialization $K({\mathbf{w}}_{0})$ can be bounded in terms of the norm of the Hessian $H$, the second derivative of $F$ using, essentially, the mean value theorem. We can bound the operator norm to get the following inequality (see ):

where ${\mathcal{B}}_{R}$ is a ball of radius $R$ around ${\mathbf{w}}_{0}$.

Using standard eigenvalue perturbation bounds we have

Recall (Eq. 12) that for networks of width $m$ with linear last layer $\|H\|=O(1/\sqrt{m})$. On the other hand, it can be shown (e.g., and for shallow and deep networks respectively) that $\lambda_{\text{m}in}(K)({\mathbf{w}}_{0})=O(1)$ and is essentially independent of the width. Hence Eq. 21 guarantees that given any fixed radius $R$, for a sufficiently wide network $\lambda_{\text{m}in}(K)({\mathbf{w}})$ is separated from zero in the ball ${\mathcal{B}}_{R}$. Thus the loss function satisfies the PL* condition in ${\mathcal{B}}_{R}$. As we discussed above, this guarantees the existence of global minima of the loss function and convergence of gradient descent for wide neural networks with linear output layer.

3.2 Transformation control

Another way to control the condition number of a system is by representing it as a composition of two or more well-conditioned maps.

Informally, due to the chain rule, if $F$ is well conditioned, so is $\phi\circ F\circ\psi({\mathbf{w}})$, where

are maps with non-degenerate Jacobian matrices.

In particular, combining Hessian control with transformation control, can be used to prove convergence for wide neural networks with non-linear last layer .

4 Efficient optimization by SGD

We have seen that over-parameterization helps explain why Gradient Descent can reach global minima even for highly non-convex optimization landscapes. Yet, in practice, GD is rarely used. Instead, mini-batch stochastic methods, such as SGD or Adam are employed almost exclusively. In its simplest form, mini-batch SGD uses the following update rule:

Here $\{({\mathbf{x}}_{i_{1}},y_{i_{1}}),\ldots,({\mathbf{x}}_{i_{m}},y_{i_{m}})\}$ is a mini-batch, a subset of the training data of size $m$, chosen at random or sequentially and $\eta>0$ is the learning rate.

At a first glance, from a classical point of view, it appears that GD should be preferable to SGD. In a standard convex setting GD converges at an exponential (referred as linear in the optimization literature) rate, where the loss function decreases exponentially with the number of iterations. In contrast, while SGD requires a factor of $\frac{n}{m}$ less computation than GD per iteration, it converges at a far slower sublinear rate (see for a review), with the loss function decreasing proportionally to the inverse of the number of iterations. Variance reduction techniques can close the gap theoretically but are rarely used in practice.

As it turns out, interpolation can explain the surprising effectiveness of plain SGD compared to GD and other non-stochastic methodsNote that the analysis is for the convex interpolated setting. While bounds for convergence under the PL* condition are available , they do not appear to be tight in terms of the step size and hence do not show an unambiguous advantage over GD. However, empirical evidence suggests that analogous results indeed hold in practice for neural networks.

The key observation is that in the interpolated regime SGD with fixed step size converges exponentially fast for convex loss functions. The results showing exponential convergence of SGD when the optimal solution minimizes the loss function at each point go back to the Kaczmarz method for quadratic functions, more recently analyzed in . For the general convex case, it was first shown in . The rate was later improved in .

Intuitively, exponential convergence of SGD under interpolation is due to what may be termed “automatic variance reduction”(). As we approach interpolation, the loss at every data point nears zero, and the variance due to mini-batch selection decreases accordingly. In contrast, under classical under-parameterized settings, it is impossible to satisfy all of the constraints at once, and the mini-batch variance converges to a non-zero constant. Thus SGD will not converge without additional algorithmic ingredients, such as averaging or reducing the learning rate. However, exponential convergence on its own is not enough to explain the apparent empirical superiority of SGD. An analysis in , identifies interpolation as the key to efficiency of SGD in modern ML, and provides a sharp computational characterization of the advantage in the convex case. As the mini-batch size $m$ grows, there are two distinct regimes, separated by the critical value $m^{*}$:

Linear scaling: One SGD iteration with mini-batch of size $m\leq m^{*}$ is equivalent to $m$ iterations of mini-batch of size one up to a multiplicative constant close to $1$.

(saturation) One SGD iterations with a mini-batch of size $m>m^{*}$ is as effective (up to a small multiplicative constant) as one iteration of SGD with mini-batch $m^{*}$ or as one iteration of full gradient descent.

For the quadratic model, $m^{*}=\frac{\max_{i=1}^{n}\{{\|{\mathbf{x}}_{i}\|}^{2}\}}{\lambda_{max}(H)}\leq\frac{{\operatorname{tr}}(H)}{\lambda_{max}(H)}$ , where $H$ is the Hessian of the loss function and $\lambda_{max}$ is its largest eigenvalue. This dependence is graphically represented in Fig. 10 from .

Thus, we see that the computational savings of SGD with mini-batch size smaller than the critical size $m^{*}$ over GD are of the order $\frac{n}{m^{*}}\approx n\frac{\lambda_{max}(H)}{{\operatorname{tr}}(H)}$. In practice, at least for kernel methods $m^{*}$ appears to be a small number, less than $100$ . It is important to note that $m^{*}$ is essentially independent on $n$ – we expect it to converge to a constant as $n\to\infty$. Thus, small (below the critical batch size) mini-batch SGD, has $O(n)$ computational advantage over GD.

Odds and ends

The attentive reader will note that most of our optimization discussions (as well as in much of the literature) involved the square loss. While training using the square loss is standard for regression tasks, it is rarely employed for classification, where the cross-entropy loss function is the standard choice for training. For two class problems with labels $y_{i}\in\{1,-1\}$ the cross-entropy (logistic) loss function is defined as

A striking aspect of cross-entropy is that in order to achieve zero loss we need to have $y_{i}f({\mathbf{x}}_{i})=\infty$. Thus, interpolation only occurs at infinity and any optimization procedure would eventually escape from a ball of any fixed radius. This presents difficulties for optimization analysis as it is typically harder to apply at infinity. Furthermore, since the norm of the solution vector is infinite, there can be no transition to linearity on any domain that includes the whole optimization path, no matter how wide our network is and how tightly we control the Hessian norm (see Section 3.10). Finally, analyses of cross-entropy in the linear case suggest that convergence is much slower than for the square loss and thus we are unlikely to approach interpolation in practice.

Thus the use of the cross-entropy loss leads us away from interpolating solutions and toward more complex mathematical analyses. Does the prism of interpolation fail us at this junction?

The accepted justification of the cross-entropy loss for classification is that it is a better “surrogate” for the 0-1 classification loss than the square loss (e.g., , Section 8.1.2). There is little theoretical analysis supporting this point of view. To the contrary, very recent theoretical works prove that in certain over-parameterized regimes, training using the square loss for classification is at least as good or better than using other loss functions. Furthermore, extensive empirical evaluations conducted in show that modern neural architectures trained with the square loss slightly outperform same architectures trained with the cross-entropy loss on the majority of tasks across several application domains including Natural Language Processing, Speech Recognition and Computer Vision.

A curious historical parallel is that current reliance on cross-entropy loss in classification reminiscent of the predominance of the hinge loss in the era of the Support Vector Machines (SVM). At the time, the prevailing intuition had been that the hinge loss was preferable to the square loss for training classifiers. Yet, the empirical evidence had been decidedly mixed. In his remarkable 2002 thesis , Ryan Rifkin conducted an extensive empirical evaluation and concluded that “the performance of the RLSC [square loss] is essentially equivalent to that of the SVM [hinge loss] across a wide range of problems, and the choice between the two should be based on computational tractability considerations”.

We see that interpolation as a guiding principle points us in a right direction yet again. Furthermore, by suggesting the square loss for classification, it reveals shortcomings of theoretical intuitions and the pitfalls of excessive belief in empirical best practices.

2 Interpolation and adversarial examples

A remarkable feature of modern neural networks is existence of adversarial examples. It was observed in that by adding a small, visually imperceptible, perturbation of the pixels, an image correctly classified as “dog” can be moved to class “ostrich” or to some other obviously visually incorrect class. Far from being an isolated curiosity, this turned out to be a robust and ubiquitous property among different neural architectures. Indeed, modifying a single, carefully selected, pixel is frequently enough to coax a neural net into misclassifying an image .

The full implications and mechanisms for the emergence of adversarial examples are not yet fully understood and are an active area of research. Among other things, the existence and pervasiveness of adversarial examples points to the limitations of the standard iid models as these data are not sampled from the same distribution as the training set. Yet, it can be proved mathematically that adversarial examples are unavoidable for interpolating classifiers in the presence of label noise (Theorem 5.1). Specifically, suppose $f_{\rm int}$ is an interpolating classifier and let ${\mathbf{x}}$ be an arbitrary point. Assume that $f_{\rm int}({\mathbf{x}})=y$ is a correct prediction. Given a sufficiently large dataset, there will be at least one ”noisy” point ${\mathbf{x}}_{i},y_{i},$, such as $f^{*}({\mathbf{x}}_{i})\neq y_{i}$, in a small neighborhood of ${\mathbf{x}}$ and thus a small perturbation of ${\mathbf{x}}$ can be used to flip the label.

If, furthermore, $f_{\rm int}$ is a consistent classifier, such as predictors discussed in Section 3.5.3, it will approach the optimal predictor $f^{*}$ as the data size grows.

Specifically, consider the set where predictions of $f_{\rm int}$ differ from the optimal classification

where $\mu$ is marginal probability measure of the data distribution. On the other hand, as $n\to\infty$, $S_{n}$ becomes a dense subset of the data domain. This can be thought of as a raisin breadAny similarity to the “plum pudding” model of the atom due to J.J.Thompson is purely coincidental.. The are the incorrect classification basins around each misclassified example, i.e., the areas where the output of $f_{\rm int}$ differs from $f^{*}$. While the seeds permeate the bread, they occupy negligible volume inside.

This picture is indeed consistent with the extensive empirical evidence for neural networks. A random perturbation avoids adversarial “raisins” , yet they are easy to find by targeted optimization methods such as PCG . I should point out that there are also other explanations for adversarial examples . It seems plausible that several mathematical effects combine to produce adversarial examples.

Summary and thoughts

We proceed to summarize the key points of this article and conclude with a discussion of machine learning and some key questions still unresolved.

The sharp contrast between the “classical” and “modern” regimes in machine learning, separated by the interpolation threshold, in various contexts, has been a central aspect of the discussion in this paper. A concise summary of some of these differences in a single table is given below.

2 Through a glass darkly

In conclusion, it may be worthwhile to discuss some of the many missing or nebulous mathematical pieces in the gradually coalescing jigsaw puzzle of deep learning.

To my mind, the most puzzling question of machine learning is why inverse methods, requiring optimization or inversion, generally perform better than direct methods such as nearest neighbors. For example, a kernel machine with a positive definite kernel $K({\mathbf{x}},{\mathbf{z}})$, appears to perform consistently and measurably better than a Nadaraya-Watson (NW) classifier using the same kernel (or the same family of kernels), despite the fact that both have the same functional form

The difference is that for a kernel machine ${\boldsymbol{\alpha}}=(K)^{-1}{\mathbf{y}}$, which requires a kernel matrix inversionRegularization, e.g., ${\boldsymbol{\alpha}}=(K+\lambda I)^{-1}{\mathbf{y}}$ does not change the nature of the problem., while NW (for classification) simply puts ${\boldsymbol{\alpha}}={\mathbf{y}}$.

The advantage of inverse methods appears to be a broad empirical pattern, manifested, in particular, by successes of neural networks. Indeed, were it not the case that inverse methods performed significantly better, the Machine Learning landscape would look quite different – there would be far less need for optimization techniques and, likely, less dependence on the availability of computational resources. I am not aware of any compelling theoretical analyses to explain this remarkable empirical difference.

A related question is that of the inductive bias. In over-parameterized settings, optimization methods, such as commonly used SGD and Adam , select a specific point ${\mathbf{w}}^{*}$ in the set of parameters ${\mathcal{S}}$ corresponding to interpolating solutions. In fact, given that ${\mathbf{w}}^{*}$ depends on the initialization typically chosen randomly, e.g., from a normal distribution, we should view ${\mathbf{w}}^{*}$ as sampled from some induced probability distribution $\mu$ on the subset of ${\mathcal{S}}$ reachable by optimization.

Why do parameters sampled from $\mu$ consistently generalize to data unseen in training?

While there is significant recent work on this topic, including a number of papers cited here, and the picture is becoming clearer for linear and kernel methods, we are still far from a thorough theoretical understanding of this alignment in general deep learning. Note that interpolation is particularly helpful in addressing this question as it removes the extra complication of analyzing the trade-off between the inductive bias and the empirical loss.

In this paper we have concentrated on interpolation as it provides insights into the phenomena of deep learning. Yet, in practice, at least for neural networks, precise interpolation is rarely used. Instead, iterative optimization algorithms are typically stopped when the validation loss stops decreasing or according to some other early stopping criterion.

This is done both for computational reasons, as running SGD-type algorithms to numerical convergence is typically impractical and unnecessary, but also to improve generalization, as early stopping can be viewed as a type of regularization (e.g., ) or label denoising that can improve test performance.

For kernel machines with standard Laplacian and Gaussian kernels, a setting where both early stopping and exact solutions can be readily computed, early stopping seems to provide at best a modest improvement to generalization performance . Yet, even for kernel machines, computational efficiency of training on larger datasets seems to require iterative methods similar to SGD , thus making early stopping a computational necessity.

Despite extensive experimental work, the computational and statistical trade-offs of early stopping in the non-convex over-parameterized regimes remain murky.

A remarkable recent theoretical deep learning discovery (discussed in Section 3.10) is that in certain regimes very wide neural networks are equivalent to kernel machines. At this point much of the theoretical discussion centers on understanding the “rich regimes” (e.g., )), often identified with “feature learning”, i.e., learning representations from data. In these regimes, tangent kernels change during the training, hence neural networks are not approximated by kernel machines, i.e., a feature map followed by a linear method. The prevalent view among both theoreticians and practitioners, is that success of neural networks cannot be explained by kernel methods as kernel predictors. Yet kernel change during training does not logically imply useful learning and may be an extraneous side effect. Thus the the question of equivalence remains open. Recent, more sophisticated, kernel machines show performance much closer to the state-of-the-art on certain tasks but have not yet closed the gap with neural networks.

Without going into a detailed analysis of the arguments (unlikely to be fruitful in any case, as performance of networks has not been conclusively matched by kernels, nor is there a convincing “smoking gun” argument why it cannot be), it is worth outlining three possibilities:

Neural network performance has elements which cannot be replicated by kernel machines (linear optimization problems).

Neural networks can be approximated by data-dependent kernels, where the kernel function and the Reproducing Kernel Hilbert Space depend on the training data (e.g., on unlabeled training data like “warped RKHS” ).

Neural networks in practical settings can be effectively approximated by kernels, such as Neural Tangent Kernels corresponding to infinitely wide networks .

I am hopeful that in the near future some clarity on these points will be achieved.

Last and, possibly, least, we would be remiss to ignore the question of depth in a paper with deep in its title. Yet, while many analyses in this paper are applicable to multi-layered networks, it is the width that drives most of the observed phenomena and intuitions. Despite recent efforts, the importance of depth is still not well-understood. Properties of deep architectures point to the limitations of simple parameter counting – increasing the depth of an architecture appears to have very different effects from increasing the width, even if the total number of trainable parameters is the same. In particular, while wider networks are generally observed to perform better than more narrow architectures (, even with optimal early stopping ), the same is not true with respect to the depth, where very deep architectures can be inferior . One line of inquiry is interpreting depth recursively. Indeed, in certain settings increasing the depth manifests similarly to iterating a map given by a shallow network . Furthermore, fixed points of such iterations have been proposed as an alternative to deep networks with some success . More weight for this point of view is provided by the fact that tangent kernels of infinitely wide architectures satisfy a recursive relationship with respect to their depth .

Acknowledgements

A version of this work will appear in Acta Numerica. I would like to thank Acta Numerica for the invitation to write this article and its careful editing. I thank Daniel Hsu, Chaoyue Liu, Adityanarayanan Radhakrishnan, Steven Wright and Libin Zhu for reading the draft and providing numerous helpful suggestions and corrections. I am especially grateful to Daniel Hsu and Steven Wright for insightful comments which helped clarify exposition of key concepts. The perspective outlined here has been influenced and informed by many illuminating discussions with collaborators, colleagues, and students. Many of these discussions occurred in spring 2017 and summer 2019 during two excellent programs on foundations of deep learning at the Simons Institute for the Theory of Computing at Berkeley. I thank it for the hospitality. Finally, I thank the National Science Foundation and the Simons Foundation for financial support.

References