Machine Learning from a Continuous Viewpoint

Introduction

We present a continuous formulation of machine learning. As usual this continuous formulation consists of three components: a representation of functions, a loss functional and a training dynamics. For representations of functions, we will discuss the integral transform-based models and the more advanced flow-based models. For the loss functional, we give examples that arise in supervised and unsupervised learning, as well as examples from calculus of variations and partial differential equations (PDEs). For training dynamics, we divide the unknown parameters into two classes: conserved and non-conserved. For non-conserved parameters, we use what is known in the physics literature as the model A dynamics , namely gradient flow in the usual $L^{2}$ metric. For conserved parameters, we use what is known as the model B dynamics , namely the gradient flow in the Wasserstein metric .

In this framework, machine learning becomes a calculus of variations or PDE-like problem, and different numerical algorithms can be used to discretize these continuous models. In particular, two-layer neural network and deep residual neural network (ResNet) models can be recovered, in a scaled form, when the particle method is applied to particular versions of the integral transform-based and flow-based models respectively. New machine learning models and algorithms can also be constructed using this continuous framework. As examples we will discuss a new flow-based random feature model, a new class of transform-based model, smoothed particle methods and spectral methods.

In addition to recovering existing machine learning models and constructing new ones, this continuous framework is also useful for the theoretical understanding of machine learning. We conjecture that the at least for standard supervised learning, the variational problems that arise from minimizing the population and empirical risks are nice variational problems in this continuous formulation, though the precise meaning of this remains to be clarified. Thus, the training models are some versions of the gradient flow of a reasonably nice functional. Hence it is not surprising that stable numerical discretizations of these continuous models perform well. In particular, this continuous viewpoint suggests that over-parametrized models should behave better since they give rise to more accurate discretizations of the continuous gradient flow. The behavior of the training algorithm should follow more closely the behavior of the continuous gradient flow. This viewpoint also suggests that one should expect trouble for very deep fully connected neural network models (which are not ResNets) since they do not have continuum limits. In fact, they suffer from numerical instabilities in the form of exploding gradients .

This work builds upon previous work. In particular, various components of this continuous framework have already appeared in the following set of works.

Continuous differential equation formulation of machine learning .

The work on the integral representations of shallow neural networks .

Mean field analysis of (stochastic) gradient descent for two-layer neural networks and multi-layer fully connected networks .

Also related are the work in . The work presented here is a natural extension of these ideas. However, the present paper is the first that systematically explores the continuous viewpoint.

Philosophically the approach advocated here bears a lot of similarity to that of the PDE-type models in image processing, such as the Mumford-Shah model or the Rudin-Osher-Fatemi model . There the idea is to first formulate a continuous variational problem that presumably represents the “first principle” for the particular image processing task such as image denoising, and then discretize that continuous problem to obtain a specific algorithm. This is in contrast to the more traditional approach in image processing in which different types of filters or algorithms are applied directly to the image, without the need to formulate the underlying mathematical problem first. This latter approach resembles the current practice in machine learning in which different algorithms are applied directly to the given dataset.

Despite its unprecedented successes across a wide spectrum of applications, machine learning still remains to be unsatisfactory as a scientific discipline. The main problem is the lack of fundamental guiding principles for designing machine learning models and algorithms and understanding their performance. Many of the techniques used in practice are still quite ad hoc, and require heavy parameter tuning. Often times the performance of these models is quite fragile and sensitive to the choice of the hyper-parameters.

The situation is reminiscent of what happened during the 1950’s when finite difference and finite element methods were just invented and used to solve PDEs. The performance of the algorithms were found to be sensitive to the particular discretization schemes and particular finite element meshes used. Some seemingly reasonable schemes simply did not run, since they quickly led to overflow on the computer. Some schemes performed reasonably well on coarse grids but blowed up upon refining grids.

Efforts for building a theoretical foundation for these finite difference and finite element methods did not go smoothly either. For example, to explain the overflow phenomenon often encountered in practice, different stability concepts and criteria were proposed . Some were easy to use in practice, but were not robust under perturbations. One such example was the concept of weak stability . Some were more robust but difficult to use in practice. It took a while for the numerical analysis community to finally settle down on the right concepts and criteria . But after all the dusts were settled, what emerged was a solid and reasonably simple picture about the basic concepts and principles behind the designing and understanding of these algorithms .

Our current work is very much motivated by the same objective, namely to develop a reasonably simple and transparent framework for machine learning. However, there is a key difference between machine learning and classical numerical analysis: While classical numerical analysis is mainly concerned with problems in low dimension, machine learning has to face problems in very high dimensions. In fact, our interest is really on machine learning models and algorithms that can overcome “the curse of dimensionality”. While the exact meaning of this terminology requires qualification, the dimensionality issue is certainly among the most important considerations in developing machine learning models and machine learning theory today. In fact, one can roughly divide all machine learning models into two categories: The ones that do suffer from the curse of dimensionality and the ones that do not. We refer to for a discussion on this.

One important class of algorithms that do not suffer from the curse of dimensionality problem are the Monte Carlo algorithms for numerical integration. In this case, one can establish simple dimension-independent error rates. In contrast, grid-based numerical integration methods such as Simpson’s rule do not share this property. Indeed, their performance deteriorates rapidly as the dimensionality goes up.

This example has some important consequences on the formulation that we will present:

We will focus on ways of representing functions as expectations, since there are algorithms for computing expectations with dimension-independent error rates.

For the same reason, particle methods stand out for the training dynamics since they are the analog of Monte Carlo methods for dynamic problems.

There are important aspects of machine learning that cannot be easily formulated at the continuous level. One example is the stochastic gradient descent algorithm.

Representations of functions

We are mainly interested in representations that are potentially effective in high dimensions. Therefore we will focus on the ones that can be expressed as expectations. As an example, instead of the Fourier representation:

where the sum is performed on a regular grid $\{\bm{\omega}_{j}\}_{j=1}^{m}$ in the Fourier space. It is well-known that this kind of grid-based approximations satisfies

where $C(f)$ and $\alpha$ are fixed quantities depending on $f$. The appearance of $1/d$ in the exponent of $m$ signals the curse of dimensionality. In contrast, for (2), by independently sample $\{\bm{\omega}_{j}\}_{j=1}^{m}$ from $\pi$, we obtain an approximation to $f$ with a dimension-independent error rate:

where $\rho(da,d\bm{\omega})=\delta(a-a(\bm{\omega}))da\pi(d\bm{\omega})$.

From an algorithmic viewpoint, (1) is typically associated with non-adaptive discretizations such as the spectral method or the ridglets and curvelets used in signal processing . We will see later that the forms (2) and (5) are closely associated with the random feature model and the two-layer neural network model. In fact one can write

where $\sigma$ is defined by $\sigma(z)=e^{iz}$. This is a two-layer neural network with activation function $\sigma$.

The original ridgelet transform representation is as follows

This representation is better suited for high dimensional situations. More generally, one can also use:

High co-dimensional representation

Ridgelet transforms express function in terms of superpositions of ridge-like structures which are co-dimension one objects. One can also refine this representation, using structures of high co-dimension. For example, the following representation uses co-dimension 2 objects:

where $\bm{w}=(\bm{w}_{1},\bm{w}_{2})$, $\sigma_{1}$ and $\sigma_{2}$ are two nonlinear scalar functions.

More generally, we can consider functions of the form

where $\bm{w}\in\Omega$ and $\rho\in\mathcal{P}(\Omega)$. Note that (6) and (2.1) are both special cases of the representation above.

Compositional structure

The representations discussed above correspond to neural network models with one hidden layer. It can be straightforwardly extended to include more hidden layers using a compositional structure. An example with two hidden layers is given by:

Another way to construct compositional structures is as follows:

2 Flow-based representation

In the flow-based representation, the trial functions are generated by the flow map of a (continuous) dynamical system:

The flow-map at time $1$ is defined as the map: $\bm{x}\rightarrow\bm{z}(1)$. More generally, one can allow a change of dimension between $\bm{x}$ and $\bm{z}$:

The set of functions that can be generated this way depend on how we choose $g$. One natural way is to use the representation discussed above, e.g.

Here $(\pi_{\tau})_{\tau\in}$ is a family of probability distributions parametrized by $\tau$. This gives us the flow:

If we compare the models in (14) and (15) with the model proposed originally in :

we see that (17) is the special case of (15) with $\rho_{\tau}=\delta(\bm{a}-\bm{u}(\tau))\delta(\bm{w}-\bm{w}(\tau))$.

However, as we learn from the work of , the more general representation in (14) and (15) is needed in order to capture the continuum limit of residual neural networks.

The optimization problem

The next step is to formulate the loss function that will be used in order to turn the problem into an optimization problem. In the continuous setting, these optimization problems are calculus of variations problems. Here we will discuss four examples of machine learning tasks.

In the following we will use $\theta$ to denote generically the set of parameters that occur in the representation. For example, for (6), we have $\theta=(a(\cdot),\pi(\cdot))$.

In supervised learning, our objective is to find the best approximation of some target function $f^{*}$ that minimizes the so-called population risk:

Here $\mu$ is a probability distribution. (18) is the $L^{2}$ loss function. Obviously one can also define other loss functions by replacing the square function by some other convex functions with the global minimum at the origin. The general form of loss function is given by:

In reality, we are only given partial information about $f^{*}$ and $\mu$ through a finite sample: $S=\{(\bm{x}_{i},y_{i})\}_{i=1}^{n}$ where $y_{i}=f^{*}(\bm{x}_{i})$ is the label for $\bm{x}_{i}$. Therefore in practice we have to work instead with the “empirical risk”:

2 Dimension reduction

This is the analog of the population risk. Again in practice, one has to work with the empirical risk, defined by:

3 Calculus of variations

where $\varphi^{*}$ is the complex conjugate of $\varphi$. This energy can also be rewritten as

where $\mu_{\varphi}$ is the probability distribution defined by:

$\mathcal{I}$ serves as the analog of the population risk.

In these problems, one typically attempts to compute $\mathcal{I}$ accurately by producing sufficient number of samples from the distribution $\mu_{\varphi}$. This means that one attempts to work directly with the population risk in these problems. In practice, however, there is an issue that the errors in the approximation of $\varphi$ may interact with the errors in the sampling. This issue has not been systematically investigated yet.

4 Nonlinear parabolic PDEs

An important application of machine learning is the numerical solution of high dimensional PDEs . Formulating these PDEs as variational problems is an important step in formulating machine learning based algorithms. In principle, one can always use the “least square” approach, as was done in . But better performance can be achieved if more sophisticated formulations are used.

with the terminal condition $u(T,x)=g(x)$. Among other things, this kinds of PDEs arise in option pricing with default risk or other nonlinear effects taken into account.

It can be shown that this PDE problem is equivalent to the following variational problem

The constraints are backward stochastic differential equations (BSDE) . This was the starting point of the “Deep BSDE method” proposed in .

Analyzing these variational problems is a major task in the mathematical theory of machine learning.

From now on we will focus on the supervised learning problem.

Gradient flows

The third component in machine learning is an algorithm for solving the optimization problem. In this section, we will discuss various gradient flow dynamics for the population or empirical risk. For simplicity we focus on the following loss functional

We first discuss gradient flows using a physics language . The loss functions or functionals defined above serve as the “free energy” of the problem.

To begin with, we need to distinguish conserved and non-conserved “order parameters”. The coefficient $a$ in (6) is non-conserved. The probability distributions $\pi$ or $\rho$ are obviously conserved.

First, let us examine the situation with the representation (6). Let $I=I(a,\pi)$ be the loss functional. Denote by $\frac{\delta I}{\delta a}$ and $\frac{\delta I}{\delta\pi}$ the formal variational derivative of $I$ with respect to $a$ and $\pi$ respectively, under the standard $L^{2}$ metric. The gradient flow for $a$ is simply given by

In the physics literature, this is known as the “model A” dynamics .

The gradient flow for $\pi$ is given by a continuity equation:

where the current $\mathbf{J}$ is given by:

This is known as the “model B” dynamics and $V$ is known as the “chemical potential”.

It is well-known that the model B dynamics is also the gradient flow under the 2-Wasserstein metric .

For flow-based models, the parameters $a$ and $\pi$ are themselves one-parameter families of coefficients or probability distributions respectively: $a=(a_{\tau})_{\tau\in},\pi=(\pi_{\tau})_{\tau\in}$. Given a functional $I,I=I((a_{\tau}),(\pi_{\tau}))$, a natural extension of the gradient flow to $(a_{\tau}),(\pi_{\tau})$ is given by:

Note that the varitional derivatives $\delta I/\delta a_{\tau}$ and $\delta I/\delta\pi_{\tau}$ appeared above are not well-defined, since $I$ is the integral of the influences of $\bm{a}_{\tau}$ and $\pi_{\tau}$ from $\tau=0$ to $\tau=1$. So strictly speaking, these derivatives are infinitesimal quantities. In section 4.3 and 4.4, we will provide rigorous forms of these equations.

The chemical potential for this functional is given by

The model B gradient flow in this case is given by

This is nothing but the “mean field” limit of the gradient descent dynamics for two-layer neural networks .

Example 2: An example of non-conserved parameter

with $\pi$ being fixed. The variational derivative of (28) with respect to $L^{2}(\pi)$ is given by

This is the continuous version of the gradient flow for random feature models.

2 Pontryagin’s maximum principle for flow-based models

Consider a general flow-based model $f(\bm{z};\theta)=\bm{1}^{T}\bm{z}_{1}^{\bm{x}}$ with $\bm{z}_{1}^{\bm{x}}$ given by the following ODE,

Minimizing the risk subject to the dynamics defined by (32) is a control problem where $\{\bm{z}^{\bm{x}}_{\tau}\}$ are the states and the parameters $\theta=\{\theta_{\tau}\}$ serve as the control. Naturally we will borrow concepts from control theory. To simplify the statement of the results, we define the following quantity:

Following the convention in control theory, we call $\bm{z},\bm{p},H$ the state, co-state and the Hamiltonian, respectively.

The Pontryagin’s maximum principle (PMP). This is a necessary condition for the optimal solutions of control problem . In the current case, let $\theta$ be a global minimum of the risk functional. Then it must satisfy

where for each $\bm{x}$, $(\bm{z}^{\bm{x}}_{\tau},\bm{p}^{\bm{x}}_{\tau})$ is given by the Hamiltonian dynamics:

Note that the dynamics of the state $\bm{z}^{\bm{x}}_{\tau}$ is forward in time from $\tau=0$ to $\tau=1$, whereas the dynamics of the co-state $\bm{p}^{\bm{x}}_{\tau}$ is backward in time from $\tau=1$ to $\tau=0$. We refer the reader to for the proof and more discussions.

3 Flow-based random feature model

First a remark about notation. We will use $t$ to denote the time for the gradient flow, and $\tau$ to denote the “time” used to define the flow-based models.

In this case, the Hamiltonian is given by

where $\bm{z}_{\tau}^{\bm{x}},\bm{p}_{\tau}^{\bm{x}}$ satisfies the following Hamiltonian dynamics,

The variational derivative of the loss functional is given by

Obviously, the co-state satisfies the following backward ODE,

Using the definition of $H$, it is easy to verify that $\bm{z}_{\tau}^{x}$ and $\bm{p}^{\bm{x}}_{\tau}$ satisfy the dynamic equations stated above. ∎

The gradient flow of the flow-based random feature model (37) is given by

where $\bm{z}^{\bm{x}}(t)$ and $\bm{p}^{\bm{x}}(t)$ are the state and co-state at time $t$ generated by $\bm{a}(\cdot,\cdot)$ through Eqn. (39).

Note that for each value of $\tau$, there is a gradient flow equation for $\bm{a}_{\tau}$. The coupling between different values of $\tau$’s is through the Eqn. (39).

Assume that $\varphi=\varphi(\bm{z},\bm{w})$ is continuous with respect to $\bm{z},w$, and there is a constant $C$ such that $\max\{|\varphi(\bm{z},\bm{w})|,\|\nabla_{\bm{z}}\varphi(\bm{z},\bm{w})\|\}\leq C$. Moreover, assume that the family $\{\sum_{k=1}^{m}a_{k}\varphi(\bm{z},\bm{w}_{k})\}$ has the universal approximation property, namely any continuous function can be uniformly approximated by functions of the form $\{\sum_{k=1}^{m}a_{k}\varphi(\bm{z},\bm{w}_{k})\}$.

By definition, the following holds for any $\bm{w}\in\Omega$

Therefore, for any $\{\bm{w}_{k}\}_{k=1}^{m}$ we have

From the universal approximation property, we obtain

The assumption implies that $\|\bm{u}(\bm{z},\tau)\|\leq C$ and $\|\nabla_{\bm{z}}\bm{u}(\bm{z},\tau)\|\leq C$. By the Picard-Lindelof theorem, the solution of ODE (51) is unique. Since $\text{rank}(V)=d+1$, the mapping $\bm{x}\to\bm{z}_{\tau}^{\bm{x}}$ is non-degenerate. Therefore, $g(\bm{z}_{\tau}^{\bm{x}})$ can represent any continuous function of $\bm{x}$. Hence, the following holds for any continuous function $h$

The proposition above is concerned with the stationary points of the loss functional. We now turn to the stationary points of the gradient flow.

Assume that $\Omega=\SS^{D}$ and $\pi_{1}$ is absolute continuous with respect to the Lebesgue measure on $\SS^{D}$. Moreover, assume that $\pi_{1}$ has a continuous, positive density on $\bm{w}\in\SS^{D}$.

The dissipation of the gradient flow (50) is given by

We first show that $g(\bm{w},\tau)$ is continuous with respect to $\tau$. Let $(\bm{z}^{x}_{\tau},\bm{p}^{x}_{\tau})$ be the solution of (39). By definition, we have for any $\tau_{1},\tau_{2}\in$,

This implies that $g(\bm{w},\cdot)$ is uniformly bounded. Therefore, we have

Since $g(\bm{w},\tau)$ is continuous with respect to $\bm{w}$ and $\pi_{1}$ has full support, we have for all $\bm{w}\in\SS^{D}$.

4 Gradient flow for the flow-based neural networks

To derive the gradient flow for (28), we first need to define the parameter space appropriately. Denote by $X:=\{\pi:\mapsto\mathcal{P}_{2}(\Omega)\}$, the space of all feasible parameters. For any $\pi^{1},\pi^{2}\in X$, consider the following metric:

where $W_{2}(\cdot,\cdot)$ is the 2-Wasserstein distance. In this case, the Hamiltonian is given by

The gradient flow in the metric space $(X,d)$ for the objective function (28) is given by

and for each $\bm{x}$, $(\bm{z}_{\tau}^{\bm{x}}(t),\bm{p}_{\tau}^{\bm{x}}(t))$ satisfies

For this gradient flow, the energy dissipation relation is given by

Moreover, from the definition of $H$, it is easy to see that

Plugging the above equation into Eqn. (4.4) leads to

Now we turn to the derivation of the gradient flow, defined as the limit of the generalized minimizing movements (GMM) scheme :

The limit of the above is exactly the 2-Wasserstein gradient flow for minimizing $H(\bm{z}_{\tau},\bm{p}_{\tau},\mu)$. This gives us

Lastly, taking expectation with respect to $\bm{x}$, we complete the proof. ∎

To make this argument rigorous, we need to establish the existence and uniqueness of the limit of GMM (72). This is a lengthy but straightforward argument. We will leave the details to interested reader.

Similar results have also been independently obtained in .

Discretizations

There are two kinds of discretization: discretization in the real space for the probability distribution $\mu$ and discretization in the parameter space for the variational problem and the flow. The discretization in the real space is relatively straightforward for a typical supervised or unsupervised learning problem. For problems in reinforcement learning or solving PDEs, this can be more tricky. However, we will skip this issue here and leave it for future work. Instead, we will focus on the discretization in the parameter space.

There are also two levels of discretization: One can either discretize the variational problem for the loss function and use one’s favorite optimization algorithm on the discretized problem, or one can discretize the continuous integral-differential equation for the training dynamics. We will focus on the latter.

Consider the functions admitting the following expectation representation,

The corresponding gradient flow is given by

where $\bm{v}$ is the velocity field given by

Let us consider the simplest particle method discretization of the model (73) and the gradient flow (74). We approximate $\pi$ by

Here $m$ is the number of particles, and $\{\bm{w}_{k}\}_{k=1}^{m}$ are the $m$ particles. In this approximation, the evolution of $\hat{\pi}$ will be completely determined by the $m$ particles.

First, the function represented by $\hat{\pi}$ is given by

where $g$ is a test function. Plugging (76) into the above equation, we get

Therefore, the dynamics of the particles follows

If we taking $\varphi(\bm{x};\bm{w})=a\sigma(\bm{b}^{T}\bm{x})$ with $\bm{w}:=(a,\bm{b})$, the particle method discretization is given by

This is the the (continuous time) gradient descent dynamics for the scaled (i.e. with the factor $1/m$ in front of the expression for $f$) two layer neural network model.

It can be shown that the dynamics described above is exactly the same as that of the GD for scaled two-layer neural networks. In fact, we have:

Given a set of initial data $\{\bm{w}_{k}^{0},k\in[m]\}$. The solution of (74) with initial data $\pi(0)=\frac{1}{m}\sum_{k=1}^{m}\delta_{\bm{w}_{k}^{0}}$ is given by

where $\{\bm{w}_{k}(\cdot),k\in[m]\}$ solves the following systems of ODEs:

In particular, this lemma says that the continuous flow equation (74) also holds for the discrete case with a finite set of neurons.

2 A smoothed particle method

A popular modification of the particle method is the smoothed particle method. Here we illustrate how one can formulate the smoothed particle method for the integral transform-based model (73) and the gradient flow (74). We will consider the special case when $\phi(\bm{x};\bm{w})=a\sigma(\bm{b}^{T}\bm{x})$. Here $\bm{w}=(a,\bm{b})$ and $\sigma(t)=\max(0,t)$ is the ReLU activation function.

Consider a smoothed particle approximation to $\pi_{t}$ This also coincides with the Gaussian mixture approximation suggested by Jianfeng Lu.

where $\phi_{h}$ is the probability density function of $\mathcal{N}(0,h^{2}I)$. The smoothed particle discretization of the flow-based model and the gradient flow is given by

where $\bm{\xi}\sim\mathcal{N}(0,I_{d+1})$. The right hand side of the last equality is the smoothed velocity.

For this to be a practical numerical algorithm, we need a way to evaluate the terms in (5.2) and (83). We will defer this to a future publication. To get some insight about the nature of this smooth particle method, we consider the special case when the data lies on the sphere, i.e. $\|\bm{x}\|=1$.

where in the last equation we have used the assumption that $\|\bm{x}\|=1$. Define a new activation function

where $\phi,\Phi$ are the probability density and cumulative density functions of the standard normal distribution, respectively. Then the discretized model can be rewritten as

From Eqn. (75) and (83), we see that the evolution of particles follows

This is exactly the gradient descent dynamics for the two-layer smoothed ReLU network (86).

It should be noted that the dominate term $t\Phi(t/h)$ in (5.2) is exactly the activation function Gaussian Error Linear Unit (GELU) , which has become quite popular recently .

3 A new algorithm for integral transform-based models

We now view both $a$ and $\pi$ as parameters.

It is tricky to design a particle method for the combined model A and model B dynamics for this problem (29) and (30). Therefore we consider instead the modified ”gradient flow”:

Let $\hat{\pi}_{t}=\frac{1}{m}\sum_{j=1}^{m}\delta(\cdot-\bm{b}_{j}(t))$. For each particle $\bm{w}_{j}$, define two quantities:

Denote by $\hat{\bm{a}}=\{a(\bm{b}_{j}(t),t)\}_{j=1}^{m}$. Then the function represented by $\hat{a}$ and $\hat{\pi}$ is given by

It is now straightforward to derive the dynamics for the particle method:

where $\{\bm{w}_{i}\}$ are randomly drawn from $\pi$.

The results for this experiment are reported in Figure 2. Figure 2 shows that for both target functions, the testing errors decrease nicely in the rate $O(1/t)$ during the training process.

The generalization error

Let $\mathcal{H},\mathcal{H}_{m}$ denote the spaces of functions represented by the continuous and discretized model, respectively. Here $m$ denotes the number of grid points or particles in the discretization. Let $S$ denote the training set and $|S|=n$. Denote by $\hat{f}_{m,n,t}$ the solution generated by the gradient descent dynamics at time $t$, and let

One way to address the generalization problem is to look for an estimate of the following type

where $\|f^{*}\|$ is some norm of the target function.

There are two ways to obtain estimates of the type in (100). One is through the a priori estimates of the discretized gradient descent dynamics. The other is through the a priori estimates of the gradient flow, i.e. the PDEs. In the following, we illustrate both approaches using the random feature model. We recover results proved in with simpler arguments.

Assume that the target function is given by

where $\pi$ is a fixed probability distribution. Assume that $|\varphi(\bm{x};\bm{b})|\leq 1$. The RKHS norm of $f^{*}$ is given by

The particle method discretization is given by

The subtlety of the problem can be appreciated from the work of which shows that the generalization error of this model can be very large in the regime where $m\approx n$.

1 Analyzing the discretized model

We decompose the generalization error into two terms:

Here $I_{1},I_{2}$ are the optimization (training) error and generalization gap, respectively.

The general philosophy is that the generalization gap is bounded by a term of the form $\|f_{m,n,t}\|/\sqrt{n}$. Here $\|\cdot\|$ is some norm determined by the model. For example, for random feature models, this is the RKHS norm. For two-layer neural network models, this is the Barron norm . Therefore to estimate the generalization gap, one needs to derive a priori bounds on these norms.

Since $\hat{\mathcal{R}}_{n}(\bm{a})$ is convex, we have $dJ/dt\leq 0$. So $J(t)\leq J(0)$, i.e.

The first inequality gives a bound on the training error. The second inequality provides a bound for the norm of the parameters.

Using (103) and the Rademacher complexity bound for the generalization gap (see Eqn. (93-95) in ), we have the following i estimates.

For any $\delta\in(0,1)$, with probability $1-\delta$ over the training examples, we have

Let $\mathcal{F}_{C}:=\{f_{m}(\cdot;\bm{a},\mathbf{B}^{0}):\|\bm{a}\|/\sqrt{m}\leq C\}$ and $\mathcal{H}_{C}:=\{(f_{m}(\cdot;\bm{a},\mathbf{B}_{0})-f^{*})^{2}:\|\bm{a}\|/\sqrt{m}\leq C\}$. By Cauchy-Schwarz inequality, $|f_{m}(\bm{x};\bm{a},B_{0})|\leq C$ and $|f^{*}(\bm{x})|\leq\sqrt{\int a(\bm{b})^{2}d\pi(\bm{b})}\leq\|f^{*}\|_{\mathcal{H}}$. Hence, $g(t)=(t-y_{i})^{2}$ is $2(C+\|f^{*}\|_{\mathcal{H}})$-Lipschitz continuous. Then by the contraction property of Rademacher complexity, we have

where the last inequality follows from the fact that $\text{Rad}(\mathcal{F}_{C})\leq\frac{C}{\sqrt{n}}$ . Hence, with probability $1-\delta$ we have for any $\bm{a}$ satisfying $\|\bm{a}\|/\sqrt{m}\leq C$,

The following proposition provides a bound on the approximation error on finite training samples.

In addition, by Hoeffding’s inequality and $|a(\bm{b}_{j}^{0})|^{2}\leq\|f\|_{\infty}^{2}$, we have

Combing Proposition 6 and Proposition 7, we have the following a priori estimates of the generalization error of GD solutions.

For any $\delta\in(0,1)$, assume that $m\geq\log^{2}(n/\delta)$. With probability $1-\delta$, we have

Taking $\bar{\bm{a}}$ be the solution constructed in Proposition 7 and plugging into Eqn. (105), we then have

where $Q_{1}=1+\log(n/\delta)t,Q_{2}=\sqrt{\log(2/\delta)/m}+t\log^{2}(2n/\delta)/m$. Plugging the above estimates into Eqn. (106), we obtain

where in the second inequality we used the fact that $\|f\|_{\mathcal{H}}\leq\|f\|_{\infty}$. Using the definition of $Q_{1},Q_{2}$ and $m\gtrsim\log^{2}(n/\delta)$ gives us that

Moreover, it follows from $m\gtrsim\log^{2}(n/\delta)$ that $I_{2}\lesssim I_{1}$. This completes the proof. ∎

2 Analyzing the continuous model

The approach presented above is the standard approach in machine learning theory. It works since the loss functional is convex in this case. It is difficult to generalize this to more complicated situations due to the lack of convexity. Here we explore an alternative approach by studying the continuous problem. Our hope is that some of the PDE techniques can be leveraged to help our understanding. One such example is found in , which proves a global convergence result for the gradient flow for two-layer neural networks by analyzing the PDE.

We decompose the generalization error as follows,

where $f_{\infty,n,t}$ is the solution given by the gradient flow of the continuous model. The three terms are respectively the discretization error, the generalization gap and the training error (for the continuous problem). The latter two terms require a priori estimates of the gradient flow.

Consider the random feature model, the gradient flow is given by

where the second inequality follows from the convexity of $\mathcal{R}_{n}$ with respect to $a$. It follows that

Since $\bm{a}_{0}=0$ and $\hat{\mathcal{R}}_{n}(a^{*})=0$, we get

Let $c=2\|a^{*}\|_{L^{2}(\pi)}$. Then the function $\hat{f}_{\infty,n,t}$ must lie in $\mathcal{F}_{c}=\{f_{\infty}(\cdot;a):\|a\|_{L^{2}(\pi)}\leq c\}$, with $f_{\infty}(x;a):=\int a(\bm{b})\varphi(\bm{x};\bm{b})d\pi(\bm{b})$. Let $\mathcal{H}_{c}=\{(f_{\infty}(\cdot;a)-f^{*})^{2}:\|a\|_{L^{2}(\pi)}\leq c\}$. By the contraction property of Rademacher complexity, we have

Moreover, $|h|\leq 4c^{2}$ for any $h\in\mathcal{H}_{c}$.

Following Eqn. (116) (117) and using the Rademacher complexity-based bound, we have

The treatment of the discretization error in (115) is more complex. This requires substantial machinery in numerical analysis. We will postpone this to future publications.

An example

In this section, we study a simple $1$-dimensional case of the integral transform-based model proposed in Section 2.1. Specifically, we consider the following conservative gradient flow,

where $w\in[0,2\pi]$, $\rho^{*}$ is a fixed probability distribution that determines the target function:

$\rho_{t}$ obeys the periodic boundary condition. $K$ is given by

It is easy to see that $K$ can be written as $K=K(w-w^{\prime})$ . Hence (125) can be written in a convolutional form,

In the following analysis we consider the case where $K$ is positive definite, i.e.,

holds for any measure $\nu$. This condition is easily satisfied in practice. In addition, we assume that $K$ is three-times differentiable and its derivatives are bounded.

First, we study the situation when $\rho^{*}$ is uniform. In this case, one can prove global convergence of the gradient flow (125). (or free energy of (125)).

Assume $\rho^{*}$ is the uniform distribution. Let $\rho_{t}$ be the solution of (126) initialized from $\rho_{0}$. Assume that $\rho_{0}$ has differentiable density function, then we have $\lim_{t\rightarrow\infty}W_{2}(\rho_{t},\rho^{*})=0$.

By an abuse of notation, we let $\rho_{t}$ and $\rho^{*}$ be the density function of $\rho_{t}$ and $\rho^{*}$, respectively. First, we assume $\rho_{t}$ exists and has differentiable density function. For any probability distribution $\rho$, consider the relative entropy of $\rho$ and $\rho^{*}$,

Let $\hat{K}(k)$, $\hat{\rho}(k)$, $\hat{\rho}^{*}(k)$ be the coefficients of the Fourier series of $K$, $\rho$ and $\rho^{*}$, respectively. The Fourier expansions exist due to the differentiability assumptions on $K$ and $\rho_{t}$. By (129) we have

Therefore, $\mathcal{H}(\rho|\rho^{*})$ is a Lyapunov function for the dynamics (126). Since the set of probability distributions $\rho$ on $[0,2\pi]$ is compact in the space $W_{2}$, any sublevel set of $\mathcal{H}$ is compact in $W_{2}$. Hence, the trajectory $\rho_{t}$ converges to the set where $\frac{d}{dt}\mathcal{H}(\rho_{t}|\rho^{*})=0$. By (130), this set contains only $\rho^{*}$. This proves the statements in the theorem.

We next prove the existence and boundedness of $\partial_{w}\rho_{t}$. The existence and uniqueness of the solution of (126) can be proved in the same way as in . Therefore we only provide the main ideas here.

Taking the partial derivative with respect to $w$ on both sides of (126), and noting that $w\in{R}$, we get

Hence, $\partial_{w}\rho_{t}$ is the solution of the following linear hyperbolic PDE for $u(w,t)$:

By the conditions on $K$, the coefficients of the PDE above are uniformly bounded. Now it follows from standard PDE argument that $u$ is bounded for any finite interval of time $[0,T]$. ∎

2 Local convergence for the general case

The previous global convergence result only holds for the case when $\rho^{*}$ is uniform. The next result shows that as long as $\rho_{0}$ is initialized close to $\rho^{*}$, the gradient flow converges to the global minimum with an $\mathcal{O}(1/t)$ rate.

Assume the conditions of Theorem 9 hold. Furthermore assume that there are constants $C_{0}$, $C_{1}$ and $C^{*}$ such that

hold for any $k\neq 0$. Let $C$ and $t_{0}$ be two constants that satisfy

Let $u_{t}=\rho_{t}-\rho^{*}$. By the conditions we imposed, $\hat{u}_{t}(0)=0$, and

for any $k\neq 0$. From equation (126), the dynamics of $u_{t}$ is

Writing (139) in the Fourier space, we get

and show that this is an invariant set for the dynamics, i.e. trajectories $(\hat{u}_{t}(k))$ initialized inside of $\mathcal{X}$ will not escape from $\mathcal{X}$. To prove this, assume that $(\hat{u}_{t}(\cdot))$ is at the boundary of $\mathcal{X}$, which means there exists a non-empty set $\mathcal{K}$ such that for any $k\in\mathcal{K}$ we have

Then, for any $k\in\mathcal{K}$, by (140) we have

For the second term on the right hand side of (143), we have

where the second inequality holds as a consequence of the condition (135), which implies

Use again condition (135) together with (145), we have

Since (147) holds for any $k\in\mathcal{K}$, the vector field at $(\hat{u}_{t}(k))$ points inside $\mathcal{X}$. Therefore, the trajectory $\{(\hat{u}_{t}(k)):t\geq 0\}$ stays in $\mathcal{X}$ for any $t>0$, which completes the proof. ∎

The theorem above shows local convergence of the gradient descent dynamics with $\mathcal{O}(1/t)$ rate. For simplicity of the proof we assumed that the Fourier coefficients of $K$ decays with an $\mathcal{O}(1/|k|)$ rate. This condition is inessential and can be relaxed, at the expense of a faster decay rate imposed on $\hat{\rho}^{*}(k)$ and $\hat{\rho}_{0}(k)$.

3 Numerical results

A pseudo-spectrum method is implemented to numerically solve the equation (126). Specifically, we consider a $1$-D model

with the feature $\varphi(\bm{x},w)$ given by

where $h$ is the standard deviation, and $x,w\in[0,2\pi]$. It is easy to see that the summation in (149) is finite for any $x$ and $w$, and $\varphi(x,w)$ is $2\pi$-periodic for both $x$ and $w$. A direct calculation gives:

In the experiments, we take $h=1$, and $\rho^{*}$ to be

The target function $f^{*}$ is displayed in the left panel of Figure 3. We see that this simple function contains three components: a mean value, a low frequency part (generated by $sin(w)$ in $\rho^{*}$), and a high frequency part (generated by $\sin(3s)$ in $\rho^{*}$). We take $\rho_{0}$ to be the uniform distribution on $[0,2\pi]$, and solve (126) for $10^{4}$ time units. The error between $f_{\rho_{t}}$ and $f^{*}$ along the path is shown in the right panel of Figure 3. We see that the dynamics proceeds in three different regimes: a nearly flat regime initially, followed by two faster regimes. This is related to the so-called frequency principle discussed next.

It is interesting to study the analog of the empirical risk, defined using the kernel:

where $\{x_{i}\}$ is a set of data samples. We take $n=100$ and sample the $x_{i}$’s from the uniform distribution on $[0,2\pi]$. The results, presented in Figure 4, suggest that the empirical loss converges to , and the $L^{2}$ norm of the density function $\rho_{t}(w)$ stays bounded. As was argued in the previous section, under this circumstance, the generalization error is bounded by $C/\sqrt{n}$.

4 The frequency principle

The frequency principle was suggested by Xu et al in . The idea was that if one uses the gradient descent to train neural network models, then the low frequency part of the target function is recovered before the high frequency component. Here we examine this issue in some detail.

For this purpose it is useful to consider the dynamics in real space, i.e. we study the evolution of the function $f$ in (148). Let $f_{t}$ be the function generated by $\rho_{t}$, we have

Therefore, the dynamics of $f_{t}$ is governed by an integral equation. This fact has important implications. To see this more clearly, let us linearize the kernel in the above equation around $\rho_{0}$, then we get

Figure 5 displays the function $f_{t}$ at different times along the gradient flow path, compared to the target function $f^{*}$. One can see that the low frequency components converge faster than the high frequency components. The is consistent with the frequency principle.

However, one should not expect this simple picture to literally hold in the general case. In Figure 6, we show the results for an example with $h=0.2$ and

In this case, for $k\leq 2/h=5$, the eigenvalue increases with $k$. Thus, we see that the high-frequency part ($k=5$) converges faster than the low frequency part ($k=1$). This is the consequence of the interplay between the frequency components in the target function and the spectrum of $K$. When there is a concentration of energy in the intermediate range of the spectrum for the target function, one should expect the scenario shown in Figure 6 to happen.

Discussions

The continuous viewpoint presented here offers a more abstract way of thinking about machine learning. Instead of thinking about features and neurons, one focuses on the representation of functions, the calculus of variation problem, and the continuous gradient flow. Features and neurons arise as objects used in special discretizations of these continuous problems.

We learn at least two things from this thought process. On one hand we can discuss machine learning without appealing to the idea of neurons, and indeed there are plenty of algorithms and models besides the neural network models. On the other hand, we also see why neural networks, both shallow and deep (ResNet), are inevitable choices: They are the simplest particle method discretization of the simplest continuous gradient flow models (for the integral transform-based and flow-based representations respectively).

One main theme in classical numerical analysis is to come up with design principles for better models and better algorithms. In that spirit, one can suggest the following set of principles for the continuous approach:

The target functions should be represented as expectations in various forms.

The risk functionals should be nice functionals. Even if not convex, they should share many features of convex functionals. A good thing is that if we start from as continuous mode, it is likely that the discretized model will not be plagued by local minima that results of discrete effects.

The different gradient flows are nice flows in the sense that the relevant norms should behave well under the flow. Here the “relevant norm” means the norm associated with the particular representation (e.g. Barron norm for the integral transform-based representation).

The numerical discretization of the flow should be stable over long time intervals.

We suspect that if one follows this set of design principles, the resulting models and algorithms will behave in a rather robust fashion, in contrast to current machine learning models which tend to depend sensitively on the choice of hyper-parameters.

Some of the subtleties in current machine learning algorithms can already be appreciated just by looking at things from a continuous viewpoint. For example, very deep fully connected networks should cause problems since they do not have nice continuum limits .

The work presented here is supported in part by a gift to Princeton University from iFlytek and the ONR grant N00014-13-1-0338. We are grateful to Jianfeng Lu, Stephan Wojtowytsch, Lexing Ying and Shuhai Zhao for helpful discussions.

References