Deep Relaxation: partial differential equations for optimizing deep neural networks

TL;DR

This paper links non-convex neural network training methods to nonlinear PDEs, showing that relaxation techniques can improve optimization performance and provide insights into the energy landscape of deep learning.

Contribution

It introduces a PDE-based relaxation approach for neural network optimization, offering a novel theoretical framework and demonstrating improved performance over traditional methods.

Findings

Relaxation techniques interpreted as viscous Hamilton-Jacobi PDE solutions.

The PDE approach outperforms stochastic gradient descent in expectation.

Analysis of the energy landscape confirms empirical observations.

Abstract

In this paper we establish a connection between non-convex optimization methods for training deep neural networks and nonlinear partial differential equations (PDEs). Relaxation techniques arising in statistical physics which have already been used successfully in this context are reinterpreted as solutions of a viscous Hamilton-Jacobi PDE. Using a stochastic control interpretation allows we prove that the modified algorithm performs better in expectation that stochastic gradient descent. Well-known PDE regularity results allow us to analyze the geometry of the relaxed energy landscape, confirming empirical evidence. The PDE is derived from a stochastic homogenization problem, which arises in the implementation of the algorithm. The algorithms scale well in practice and can effectively tackle the high dimensionality of modern neural networks.

Figures

Tables

Table 1. Table 1 . Summary of experimental results: Validation error (%) @ Effective epochs

Model	Entropy-SGD	HEAT	HJ	SGD
$\mathrm{mnistfc}$	$1.08\pm 0.02@120$	$1.13\pm 0.02@200$	$1.17\pm 0.04@200$	$1.10\pm 0.01@194$
$\mathrm{LeNet}$	$0.5\pm 0.01@80$	$0.59\pm 0.02@75$	$0.5\pm 0.01@70$	$0.5\pm 0.02@67$
All-CNN	$7.96\pm 0.05@160$	$9.04\pm 0.04@150$	$7.89\pm 0.07@145$	$7.94\pm 0.06@195$

Equations

$$\frac{\partial u}{\partial t}$$

$$f_{\gamma }(x)=u(x,\gamma )$$

$$dx(t)=-\nabla f(x(t))dt+{\beta }^{-1/2}dW(t)$$

$$dy(s)=-\left({\gamma }^{-1}y(s)+\nabla f(x(t)-y(s))\right)ds+{\beta }^{-1/2}dW(s)$$

$$dx(t)=-\int \nabla f\left(x(t)-y\right)\rho (dy)$$

$$dx(t)=-\nabla u(x(t),\gamma )dt$$

$$dx(t)=-\nabla u(x(t),T-t)dt$$

$$y(x;\xi )={\sigma }^{\prime }\left(x^p\sigma \left(x^{p-2}\dots \sigma \left(x^1\xi \right)\right)\dots \right);$$

$$f(x):=\frac{1}{N}\sum _{i=1}^N{f}_i(x);$$

$$f_i(x):=1_{\left\{y(x;{\xi }^i)\ne y^i\right\}}$$

$$x^{\ast }=\underset{x}{\mathrm{argmin}}f(x)$$

$$x^{k+1}=x^k-{\eta }_k\nabla f_{i_k}(x^{k-1})\mathrm{where}{\eta }_k>0,k\in \mathbb{N}$$

$$\mathbb{E}\left[f(x^k)\right]-f(x^{\ast })=\mathcal{O}(k^{-1/2});$$

$$\mathbb{E}\left[\nabla f_{\mathrm{mb}}(x)\right]$$

$$\mathbb{E}\left[\Vert \nabla f_{\mathrm{mb}}(x)-\nabla f(x){\Vert }^2\right]$$

$$dx(t)=-\nabla f(x(t))dt+{\beta }^{-1/2}dW(t).$$

$$\mathcal{L}\varphi =-\nabla f\cdot \nabla \varphi +\frac{{\beta }^{-1}}{2}\Delta \varphi .$$

$$u(x,t)=\mathbb{E}\Big{[}V(x(T))\leavevmode\nobreak\ \Big{|}\leavevmode\nobreak\ x(t)=x\Big{]}$$

$$\frac{\partial u}{\partial t}$$

$$\frac{\partial}{\partial t}\ \rho(x,t)=\nabla\cdot\Big{(}\nabla f(x)\ \rho(x,t)\Big{)}\ +\frac{\beta^{-1}}{2}\ \Delta\rho(x,t)$$

$${\rho }^{\infty }(x;\beta )=Z(\beta {)}^{-1}{e}^{-\beta f(x)},$$

$$dx(t)$$

$$dp(t)$$

$$J(\rho )=\int V(x)\rho dx+\frac{{\beta }^{-1}}{2}\int \rho \log \rho dx;$$

$$d_{W_2}\left(\rho (x,t),{\rho }^{\infty }\right)\le d_{W_2}\left(\rho (x,0),{\rho }^{\infty }\right)e^{-\lambda t}.$$

$$\mathbb{E}_{x_{1}}\left[\tau_{x_{2}}\right]\ \propto\frac{1}{\left|f^{{}^{\prime\prime}}(x_{3})\ f^{{}^{\prime\prime}}(x_{1})\right|^{1/2}}\exp\Big{(}\beta(f(x_{3})-f(x_{1}))\Big{)};$$

$$f_{\gamma}(x):=u(x,\gamma)=-\frac{1}{\beta}\log\Big{(}G_{\beta^{-1}\gamma}*\ \exp\left(-\beta f(x)\right)\Big{)};$$

$$\nabla u(x,t)={\int }_{{\mathbb{R}}^n}\frac{y-x}{t}{\rho }_1^{\infty }(dy;x)={\int }_{{\mathbb{R}}^n}\nabla f(x-y){\rho }_2^{\infty }(dy;x)$$

$${\rho }_1^{\infty }(y;x)=Z_1^{-1}\exp \left(-\beta f(y)-\frac{\beta }{2t}{\left|x-y\right|}^2\right){\rho }_2^{\infty }(y;x)=Z_2^{-1}\exp \left(-\beta f(x-y)-\frac{\beta }{2t}{\left|y\right|}^2\right)$$

$$v_t={\beta }^{-1}\Delta v$$

$$v_t=-\beta v{u}_t,$$

Full text

Deep Relaxation: partial differential equations

for optimizing deep neural networks

Pratik Chaudhari1∗

, 

Adam Oberman2∗

, 

Stanley Osher3

, 

Stefano Soatto1

 and 

Guillaume Carlier4

∗ Joint first authors

1 Computer Science Department, University of California, Los Angeles.

2 Department of Mathematics and Statistics, McGill University, Montreal.

3 Department of Mathematics & Institute for Pure and Applied Mathematics, University of California, Los Angeles.

4 CEREMADE, Université Paris IX Dauphine.

Email: pratikac@ucla.edu, adam.oberman@mcgill.ca, sjo@math.ucla.edu, soatto@ucla.edu, carlier@ceremade.dauphine.fr

Abstract: In this paper we establish a connection between non-convex optimization methods for training deep neural networks and nonlinear partial differential equations (PDEs). Relaxation techniques arising in statistical physics which have already been used successfully in this context are reinterpreted as solutions of a viscous Hamilton-Jacobi PDE. Using a stochastic control interpretation allows we prove that the modified algorithm performs better in expectation that stochastic gradient descent. Well-known PDE regularity results allow us to analyze the geometry of the relaxed energy landscape, confirming empirical evidence. The PDE is derived from a stochastic homogenization problem, which arises in the implementation of the algorithm. The algorithms scale well in practice and can effectively tackle the high dimensionality of modern neural networks.

Keywords: Non-convex optimization, deep learning, neural networks, partial differential equations, stochastic gradient descent, local entropy, regularization, smoothing, viscous Burgers, local convexity, optimal control, proximal, inf-convolution

1. Introduction

1.1. Overview of the main results

Deep neural networks have achieved remarkable success in a number of applied domains from visual recognition and speech, to natural language processing and robotics (LeCun et al., , 2015). Despite many attempts, an understanding of the roots of this success remains elusive. Neural networks are a parametric class of functions whose parameters are found by minimizing a non-convex loss function, typically via stochastic gradient descent (SGD), using a variety of regularization techniques. In this paper, we study a modification of the SGD algorithm which leads to improvements in the training time and generalization error.

Baldassi et al., (2015) recently introduced the local entropy loss function, $f_{\gamma}(x)$, which is a modification of a given loss function $f(x)$, defined for $x\in\mathbb{R}^{n}$. Chaudhari et al., (2016) showed that the local entropy function is effective in training deep networks. Our first result, in Sec. 3, shows that $f_{\gamma}(x)$ is the solution of a nonlinear partial differential equation (PDE). Define the viscous Hamilton-Jacobi PDE

[TABLE]

with $u(x,0)=f(x)$. Here we write $\Delta u=\sum_{i=1}^{n}\ \frac{\partial^{2}}{\partial x_{i}^{2}}u$, for the Laplacian. We show that

[TABLE]

Using standard semiconcavity estimates in Sec 6, the PDE interpretation immediately leads to regularity results for the solutions that are consistent with empirical observations.

Algorithmically, our starting point is the continuous time stochastic gradient descent equation,

[TABLE]

for $t\geq 0$ where $W(t)\in\mathbb{R}^{n}$ is the $n$-dimensional Wiener process and the coefficient $\beta^{-1}$ is the variance of the gradient term $\nabla f(x)$ (see Section 2.2 below). Chaudhari et al., (2016) compute the gradient of local entropy $\nabla f_{\gamma}(x)$ by first computing the invariant measure, $\rho(y)$, of the auxiliary stochastic differential equation (SDE)

[TABLE]

then by updating $x(t)$ by averaging against this invariant measure

[TABLE]

The effectiveness of the algorithm may be partially explained by the fact that for small value of $\gamma$, that solutions of Fokker-Planck equation associated with the SDE (1) converge exponentially quickly to the invariant measure $\rho$. Tools from optimal transportation (Santambrogio, , 2015) or functional inequalities (Bakry and Émery, , 1985) imply that there is a threshold parameter related to the semiconvexity of $f(x)$ for which the auxiliary SGD converges exponentially. See Sec. 2.3 and Remark 6. Beyond this threshold, the auxiliary problem can be slow to converge.

In Sec. 4 we prove using homogenization of SDEs (Pavliotis and Stuart, , 2008) that (2) recovers the gradient of the local entropy function, $u(x,\gamma)$. In other words, in the limit of fast dynamics, (2) is equivalent to

[TABLE]

where $u(x,\gamma)$ is the solution of (viscous-HJ). Thus while smoothing the loss function using PDEs directly is not practical, due to the curse of dimensionality, the auxliary SDE (1) accomplishes smoothing of the loss function $f(x)$ using exponentially convergent local dynamics.

The Elastic-SGD algorithm (Zhang et al., , 2015) is an influential algorithm which allowed training of Deep Neural Networks to be accomplished in parallel. We also prove, in Sec. 4.3, that the Elastic-SGD algorithm (Zhang et al., , 2015) is equivalent to regularization by local entropy. The equivalence is a surprising result, since the former runs by computing temporal averages on a single processor while the latter computes an average of coupled independent realizations of SGD over multiple processors. Once the homogenization framework is established for the local entropy algorithm, the result follows naturally from ergodicity, which allows us to replace temporal averages in (1) by the corresponding spatial averages used in Elastic-SGD.

Rather than fix one value of $\gamma$ for the evolution, scope the variable. Scoping has the advantage of more smoothing earlier in the algorithm, while recovering the true gradients (and minima) of the original loss function $f(x)$ towards the end. Scoping is equivalent to replacing (3) by

[TABLE]

where $T$ is a termination time for the algorithm. We still fix $\gamma(t)$ in the inner loop represented by (1). Scoping, which was previously considered to be an effective heuristic, is rigorously shown to be effective in Sec. 5. Additionally, scoping corresponds to nonlinear forward-backward equations, which appear in Mean Field Games, see Remark 9.

It is unusual in stochastic nonconvex optimization to obtain a proof of the superiority of a particular algorithm compared to standard stochastic gradient descent. We obtain such a result in Theorem 12 of Sec. 5 Sec. 5 where we prove (for slightly modified dynamics) an improvement in the expected value of the loss function using the local entropy dynamics as compared to that of SGD. This result is obtained using well-known techniques from stochastic control theory (Fleming and Soner, , 2006) and the comparison principle for viscosity solutions. Again, while these techniques are well-established in nonlinear PDEs, we are not aware of their application in this context.

The effectiveness of the PDE regularization is illustrated Figure 1.

2. Background

2.1. Deep neural networks

A deep network is a nested composition of linear functions of an input datum $\xi\in\mathbb{R}^{d}$ and a nonlinear function $\sigma$ at each level. The output of a neural network, denoted by $y(x;\ \xi)$, can be written as

[TABLE]

where $x^{1},\ldots,x^{p-1}\in\mathbb{R}^{d\times d}$ and $x^{p}\in\mathbb{R}^{d\times K}$ are the parameters or “weights” of the network. We will let $x\in\mathbb{R}^{n}$ denote their concatenation after rearranging each one of them as a vector. The functions $\sigma(\cdot)$ are applied point-wise to each element of their argument. For instance, in a network with rectified linear units (ReLUs), we have $\sigma(z)=\max(0,z)$ for $z\in\mathbb{R}^{d}$. Typically, the last non-linearity is set to $\sigma^{\prime}(z)={\bf 1}_{\left\{z\ \geq\ 0\right\}}$. A neural network is thus a function with $n$ parameters that produces an output $y(x;\ \xi)\in\left\{1,\ldots,K\right\}$ for each input $\xi$.

In supervised learning or “training” of a deep network for the task of image classification, given a finite sample of inputs and outputs $D=\{\xi^{i},y^{i})\}_{i=1}^{N}$ (dataset), one wishes to minimize the empirical loss111The empirical loss is a sample approximation of the expected loss, $\mathbb{E}_{x\sim P}f(x)$, which cannot be computed since the data distribution $P$ is unknown. The extent to which the empirical loss (or a regularized version thereof) approximates the expected loss relates to generalization, i.e., the value of the loss function on (“test” or “validation”) data drawn from $P$ but not part of the training set $D$.

[TABLE]

where $f_{i}(x)$ is a loss function that measures the discrepancy between the predictions $y(x,\xi^{i})$ for the sample $\xi^{i}$ and the true label $y^{i}$. For instance, the zero-one loss is

[TABLE]

whereas for regression, the loss could be a quadratic penalty $f_{i}(x):=\frac{1}{2}\lVert y(x;\ \xi^{i})-y^{i}\rVert^{2}_{2}$. Regardless of the choice of the loss function, the overall objective $f(x)$ in deep learning is a non-convex function of its argument $x$.

In the sequel, we will dispense with the specifics of the loss function and consider training a deep network as a generic non-convex optimization problem given by

[TABLE]

However, the particular relaxation and regularization schemes we describe are developed for, and tested on, classification tasks using deep networks.

2.2. Stochastic gradient descent (SGD)

First-order methods are the default choice for optimizing (6), since the dimensionality of $x$ is easily in the millions. Furthermore, typically, one cannot compute the entire gradient $N^{-1}\ \sum_{i=1}^{N}\nabla f_{i}(x)$ at each iteration (update) because $N$ can also be in the millions for modern datasets.222For example, the ImageNet dataset (Krizhevsky et al., , 2012) has $N=1.25$ million RGB images of size $224\times 224$ (i.e. $d\approx 10^{5}$) and $K=1000$ distinct classes. A typical model, e.g. the Inception network (Szegedy et al., , 2015) used for classification on this dataset has about $N=10$ million parameters and is trained by running (7) for $k\approx 10^{5}$ updates; this takes roughly $100$ hours with $8$ graphics processing units (GPUs). Stochastic gradient descent (Robbins and Monro, , 1951) consists of performing partial computations of the form

[TABLE]

where $i_{k}$ is sampled uniformly at random from the set $\left\{1,\ldots,N\right\}$, at each iteration, starting from a random initial condition $x^{0}$. The stochastic nature of SGD arises from the approximation of the gradient using only a subset of data points. One can also average the gradient over a set of randomly chosen samples, called a “mini-batch”. We denote the gradient on such a mini-batch by $\nabla f_{\mathrm{mb}}(x)$.

If each $f_{i}(x)$ is convex and the gradient $\nabla f_{i}(x)$ is uniformly Lipschitz, SGD converges at a rate

[TABLE]

this can be improved to $\mathcal{O}(k^{-1})$ using additional assumptions such as strong convexity (Nemirovski et al., , 2009) or at the cost of memory, e.g., SAG (Schmidt et al., , 2013), SAGA (Defazio et al., , 2014) among others. One can also obtain convergence rates that scale as $\mathcal{O}(k^{-1/2})$ for non-convex loss functions if the stochastic gradient is unbiased with bounded variance (Ghadimi and Lan, , 2013).

2.3. SGD in continuous time

Developments in this paper hinge on interpreting SGD as a continuous-time stochastic differential equation. As is customary (Ghadimi and Lan, , 2013), we assume that the stochastic gradient $\nabla f_{\mathrm{mb}}(x^{k-1})$ in (7) is unbiased with respect to the full-gradient and has bounded variance, i.e., for all $x\in\mathbb{R}^{n}$,

[TABLE]

for some $\beta_{\mathrm{mb}}^{-1}\geq 0$ We use the notation to emphasize that the noise is coming from the minibatch gradients. In the case where we compute full gradients, this coefficient is zero. For the sequel, we forget the source of the noise, and simply write $\beta^{-1}$. The discrete-time dynamics in (7) can then be modelled by the stochastic differential equation (SGD)

[TABLE]

The generator, $\mathcal{L}$, corresponding to (SGD) is defined for smooth functions $\phi$ as

[TABLE]

The adjoint operator $\mathcal{L}^{*}$ is given by $\mathcal{L}^{*}\rho=\nabla\cdot(\nabla f\rho)+\frac{\beta^{-1}}{2}\Delta\ \rho$. Given the function $V(x):X\to\mathbb{R}$, define

[TABLE]

to be the expected value of $V$ at the final time for the path (SGD) with initial data $x(t)=x$. By taking expectations in the Itô formula, it can be established that $u$ satisfies the backward Kolmogorov equation

[TABLE]

along with the terminal value $u(x,T)=V(x)$. Furthermore, $\rho(x,t)$, the probability density of $x(t)$ at time $t$, satisfies the Fokker-Planck equation (Risken, , 1984)

[TABLE]

along with the initial condition $\rho(x,0)=\rho_{0}(x)$, which represents the probability density of the initial distribution. With mild assumptions on $f(x)$, and even when $f$ is non-convex, $\rho(x,t)$ converges to the unique stationary solution of (FP) as $t\to\infty$ (Pavliotis, , 2014, Section 4.5). Elementary calculations show that the stationary solution is the Gibbs distribution

[TABLE]

for a normalizing constant $Z(\beta)$. The notation for the parameter $\beta>0$ comes from physics, where it corresponds to the inverse temperature. From (10), as $\beta\to\infty$, the Gibbs distribution $\rho^{\infty}$ concentrates on the global minimizers of $f(x)$. In practice, “momentum” is often used to accelerate convergence to the invariant distribution; this corresponds to running the Langevin dynamics (Pavliotis, , 2014, Chapter 6) given by

[TABLE]

In the overdamped limit $c\to\infty$ we recover (SGD) and in the underdamped limit $c=0$ we recover Hamiltonian dynamics. However, while momentum is used in practice, there are more theoretical results available for the overdamped limit.

The celebrated paper (Jordan et al., , 1998) interpreted the Fokker-Planck equation as a gradient descent in the Wasserstein metric $d_{W_{2}}$ (Santambrogio, , 2015) of the energy functional

[TABLE]

If $V(x)$ is $\lambda$-convex function $f(x)$, i.e., $\nabla^{2}V(x)+\lambda I$, is positive definite for all $x$, then the solution $\rho(x,t)$ converges exponentially in the Wasserstein metric with rate $\lambda$ to $\rho^{\infty}$ (Carrillo et al., , 2006),

[TABLE]

Let us note that a less explicit, but more general, criterion for speed of convergence in the weighted $L^{2}$ norm is the Bakry-Emery criterion (Bakry and Émery, , 1985), which uses PDE-based estimates of the spectral gap of the potential function.

2.4. Metastability and simulated annealing

Under mild assumptions for non-convex functions $f(x)$, the Gibbs distribution $\rho^{\infty}$ is still the unique steady solution of (FP). However, convergence of $\rho(x,t)$ can take an exponentially long time. Such dynamics are said to exhibit metastability: there can be multiple quasi-stationary measures which are stable on time scales of order one. Kramers’ formula (Kramers, , 1940) for Brownian motion in a double-well potential is the simplest example of such a phenomenon: if $f(x)$ is a double-well with two local minima at locations $x_{1},x_{2}\in\mathbb{R}$ with a saddle point $x_{3}\in\mathbb{R}$ connecting them, we have

[TABLE]

where $\tau_{x_{2}}$ is the time required to transition from $x_{1}$ to $x_{2}$ under the dynamics in (SGD). The one dimensional example can be generalized to the higher dimensional case (Bovier and den Hollander, , 2006). Observe that there are two terms that contribute to the slow time scale: (i) the denominator involves the Hessian at a saddle point $x_{3}$, and (ii) the exponential dependence on the difference of the energies $f(x_{1})$ and $f(x_{2})$. In the higher dimensional case the first term can also go to infinity if the spectral gap goes to zero.

Simulated annealing (Kushner, , 1987; Chiang et al., , 1987) is a popular technique which aims to evade metastability and locate the global minimum of a general non-convex $f(x)$. The key idea is to control the variance of the Wiener process, in particular, decrease it exponentially slowly by setting $\beta=\beta(t)=\log\ (t+1)$ in (7). There are theoretical results that prove that simulated annealing converges to the global minimum (Geman and Hwang, , 1986); however the time required is still exponential.

2.5. SGD heuristics in deep learning

State-of-the-art results in deep learning leverage upon a number of heuristics. Variance reduction of the stochastic gradients of Sec. 2.2 leads to improved theoretical convergence rates for convex loss functions (Schmidt et al., , 2013; Defazio et al., , 2014) and is implemented in algorithms such as AdaGrad (Duchi et al., , 2011), Adam (Kingma and Ba, , 2014) etc.

Adapting the step-size, i.e., changing $\eta_{k}$ in (7) with time $k$, is common practice. But while classical convergence analysis suggests step-size reduction of the form $\mathcal{O}(k^{-1})$ (Robbins and Monro, , 1951; Schmidt et al., , 2013), staircase-like drops are more commonly used. Merely adapting the step-size is not equivalent to simulated annealing, which explicitly modulates the diffusion term. This can be achieved by modulating the mini-batch size $i_{k}$ in (7). Note that as the mini-batch size goes to $N$ and the step-size $\eta_{k}$ is \mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.6}{\scriptscriptstyle\mathcal{O}}}(N), the stochastic term in (7) vanishes (Li et al., , 2015).

2.6. Regularization of the loss function

A number of prototypical models of deep neural networks have been used in the literature to analyze characteristics of the energy landscape. If one forgoes the nonlinearities, a deep network can be seen as a generalization of principal component analysis (PCA), i.e., as a matrix factorization problem; this is explored by Baldi and Hornik, (1989); Haeffele and Vidal, (2015); Soudry and Carmon, (2016); Saxe et al., (2014), among others. Based on empirical results such as Dauphin et al., (2014), authors in Bray and Dean, (2007); Fyodorov and Williams, (2007); Choromanska et al., (2015) and Chaudhari and Soatto, (2015) have modeled the energy landscape of deep neural networks as a high-dimensional Gaussian random field. In spite of these analyses being drawn from vastly diverse areas of machine learning, they suggest that, in general, the energy landscape of deep networks is highly non-convex and rugged.

Smoothing is an effective way to improve the performance of optimization algorithms on such a landscape and it is our primary focus in this paper. This can be done by convolution with a kernel (Chen, , 2012); for specific architectures this can be done analytically (Mobahi, , 2016) or by averaging the gradient over random, local perturbations of the parameters (Gulcehre et al., , 2016).

In the next section, we present a technique for smoothing that has shown good empirical performance. We compare and contrast the above mentioned smoothing techniques in a unified mathematical framework in this paper.

3. PDE interpretation of local entropy

Local entropy is a modified loss function first introduced as a means for studying the energy landscape of the discrete perceptron, i.e., a “shallow” neural network with one layer and discrete parameters, e.g., $x\in\left\{-1,1\right\}^{n}$ (Baldassi et al., 2016b, ). An analysis of local entropy using the replica method (Mezard et al., , 1987) predicts dense clusters of solutions to the optimization problem (6). Parameters that lie within such clusters yield better error on samples outside the training set $D$, i.e. they have improved generalization error (Baldassi et al., , 2015).

Recent papers by Baldassi et al., 2016a and Chaudhari et al., (2016) have suggested algorithmic methods to exploit the above observations. The later extended the notion of local entropy to continuous variables by replacing the loss function $f(x)$ with

[TABLE]

where $G_{\gamma}(x)={(2\pi\gamma)^{-d/2}}\exp{-\frac{\left|x\right|^{2}}{2\gamma}}$ is the heat kernel.

3.1. Derivation of the viscous Hamilton-Jacobi PDE

The Cole-Hopf transformation (Evans, , 1998, Section 4.4.1) is a classical tool used in PDEs. It relates solutions of the heat equation to those of the viscous Hamilton-Jacobi equation. For the convenience of the reader, we restate it here.

Lemma 1**.**

The local entropy function $f_{\gamma}(x)=u(x,\gamma)$ defined by (11) is the solution of the initial value problem for the viscous Hamilton-Jacobi equation (viscous-HJ) with initial values $u(x,0)=f(x)$. Moreover, the gradient is given by

[TABLE]

where $\rho^{\infty}_{i}(y;\ x)$ are probability distributions given by

[TABLE]

and $Z_{i}=Z_{i}(x)$ are normalizing constants for $i=1,2$.

Proof.

Define $u(x,t)=-{\beta^{-1}}\log v(x,t)$. From (11), $v=\exp(-\beta u)$ solves the heat equation

[TABLE]

with initial data $v(x,0)=\exp(-\beta f(x))$. Taking partial derivatives gives

[TABLE]

Combining these expressions results in (viscous-HJ). Differentiating $v(x,t)=\exp(-\beta u(x,t))$ using the convolution in (11), gives up to positive constants which can be absorbed into the density,

[TABLE]

Evaluating the last or second to last expression for the integral leads to the corresponding parts of (12) and (13).

3.2. Hopf-Lax formula for the Hamilton-Jacobi equation and dynamics for the gradient

In addition to the connection with (viscous-HJ) provided by Lemma 1, we can also explore the non-viscous Hamilton-Jacobi equation. The Hamiliton-Jacobi equation corresponds to the limiting case of (viscous-HJ) as the viscosity term $\beta^{-1}\to 0$.

There are several reasons for studying this equation. It has a simple, explicit formula for the gradient. In addition, semi-concavity of the solution follows directly from the Hopf-Lax formula. Moreover, under certain convexity conditions, the gradient of the solution can be computed as an exponentially convergent gradient flow. The deterministic dynamics for the gradient are a special case of the stochastic dynamics for the viscous-HJ equation which are discussed in another section. Let us therefore consider the Hamilton-Jacobi equation

[TABLE]

In the following lemma, we apply the well-known Hopf-Lax formula (Evans, , 1998) for the solution of the HJ equation. It is also called the $\inf$-convolution of the functions $f(x)$ and $\frac{1}{2t}\left|x\right|^{2}$ (Cannarsa and Sinestrari, , 2004). It is closely related to the proximal operator (Moreau, , 1965; Rockafellar, , 1976).

Lemma 2**.**

Let $u(x,t)$ be the viscosity solution of (HJ) with $u(x,0)=f(x)$. Then

[TABLE]

Define the proximal operator

[TABLE]

If $y^{*}=\textrm{prox}_{tf}(x)$ is a singleton, then $\nabla_{x}u(x,t)$ exists, and

[TABLE]

Proof.

The Hopf-Lax formula for the solution of the HJ equation can be found in  (Evans, , 1998). Danskin’s theorem (Bertsekas et al., , 2003, Prop. 4.5.1) allows us to pass a gradient through an infimum, by applying it at the argminimum. The first equality is a direct application of Danskin’s Theorem to (HL). For the second equality, rewrite (HL) with $z=(x-y)/t$ and the drop the prime, to obtain,

[TABLE]

Applying Danskin’s theorem to the last expresion gives

[TABLE]

Next we give a lemma which shows that under an auxiliary convexity condition, we can find $\nabla(x,t)$ using exponentially convergent gradient dynamics.

Lemma 3** (Dynamics for HJ).**

For a fixed $x$, define

[TABLE]

Suppose that $t>0$ is chosen so that $h(x,y;t)$ is $\lambda$-convex as a function of $y$ (meaning that $D^{2}_{y}\ h(x,y,t)-\lambda I$ is positive definite). The gradient $p=\nabla_{x}\ u(x,t)$ is then the unique steady solution of the dynamics

[TABLE]

Moreover, the convergence to $p$ is exponential, i.e.,

[TABLE]

Proof.

Note that (16) is the gradient descent dynamics on $h$ which is $\lambda$-convex by assumption. Thus by standard techniques from ordinary differential equation (ODE) (see, for example (Santambrogio, , 2015)) the dynamics are contraction, and the convergence is exponential.

3.3. The invariant measure for a quadratic function

Local entropy computes a non-linear average of the gradient in the neighborhood of $x$ by weighing the gradients according to the steady-state measure $\rho^{\infty}(y;\ x)$.

It is useful to compute $\nabla f_{\gamma}(x)$ when $f(x)$ is quadratic, since this gives an estimate of the quasi-stable invariant measure near a local minimum of $f$. In this case, the invariant measure $\rho^{\infty}(y;\ x)$ is a Gaussian distribution.

Lemma 4**.**

Suppose $f(x)$ is quadratic, with $p=\nabla\ f(x)$ and $Q=\nabla^{2}f(x)$, then the invariant measure $\rho^{\infty}_{2}(y;\ x)$ of (13) is a normal distribution with mean $\mu$ and covariance $\Sigma$ given by

[TABLE]

In particular, for $\gamma>\left|Q\right|_{2}$, we have

[TABLE]

Proof.

Without loss of generality, a quadratic $f(x)$ centered at $x=0$ can be written as

[TABLE]

by completing the square. The distribution $\exp\left(-f(y)-\frac{1}{2\gamma}\ \left|y\right|^{2}\right)$ is thus a normal distribution with mean $\mu=x-\Sigma\ p$ and variance $\Sigma=\left(Q+\gamma\ I\right)^{-1}$ and an appropriate normalization constant. The approximation $\Sigma=\gamma^{-1}\ I-\gamma^{-2}Q$ which avoids inverting the Hessian follows from the Neumann series for the matrix inverse; it converges for $\gamma>\left|Q\right|_{2}$ and is accurate up to $\mathcal{O}\left(\gamma^{-3}\right)$.

4. Derivation of local entropy via homogenization of SDEs

This section introduces the technique of homogenization for stochastic differential equations and gives a mathematical treatment to the result of Lemma 1 which was missing in Chaudhari et al., (2016). In addition to this, homogenization allows us to rigorously show that an algorithm called Elastic-SGD (Zhang et al., , 2015) that was heuristically connected to a distributed version of local entropy in Baldassi et al., 2016a is indeed equivalent to local entropy.

Homogenization is a technique used to analyze dynamical systems with multiple time-scales that have a few fast variables that may be coupled with other variables which evolve slowly. Computing averages over the fast variables allows us to obtain averaged equations for the slow variables in the limit that the times scales separate. We refer the reader to Pavliotis and Stuart, (2008, Chap. 10, 17) or E, (2011) for details.

4.1. Background on homogenization

Consider the system of SDEs given by

[TABLE]

where $h,g$ are sufficiently smooth functions, $W(s)\in\mathbb{R}^{n}$ is the standard Wiener process and $\epsilon>0$ is the homogenization parameter which introduces a fast time-scale for the dynamics of $y(s)$. Let us define the generator for the second equation to be

[TABLE]

We can assume that, for fixed $x$, the fast-dynamics, $y(s)$, has a unique invariant probability measure denoted by $\rho^{\infty}(y;\ x)$ and the operator $\mathcal{L}_{0}$ has a one-dimensional null-space characterized by

[TABLE]

here $1(y)$ are all constant functions in $y$ and $\mathcal{L}_{0}^{*}$ is the adjoint of $\mathcal{L}_{0}$. In that case, it follows that in the limit $\epsilon\to 0$, the dynamics for $x(s)$ in (17) converge to

[TABLE]

where the homogenized vector field for $X$ is defined as the average agains the invariant measure. Moreover, by ergodicity, in the limit, the spatial average is equal to a long term average over the dynamics independent of the initial value $y(0)$.

[TABLE]

4.2. Derivation of local entropy via homogenization of SDEs

Recall (12) of Lemma 1

[TABLE]

Hence, let us consider the following system of SDEs

[TABLE]

Write

[TABLE]

The Fokker-Planck equation for the density of $y(s)$ is given by

[TABLE]

The invariant measure for this Fokker-Planck equation is thus

[TABLE]

which agrees with (13) in Lemma 1.

Theorem 5**.**

As $\epsilon\to 0$, the system (Entropy-SGD) converges to the homogenized dynamics given by

[TABLE]

Moreover, $-\nabla f_{\gamma}(x)=-\gamma^{-1}\ \left(x-\left\langle y\right\rangle\right)$ where

[TABLE]

where $y(s)$ is the solution of the second equation in (Entropy-SGD) for fixed $x$.

Proof.

The result follows immediately from the convergence result in Sec. 4.1 for the system (17). The homogenized vector field is

[TABLE]

which is equivalent to $\nabla f_{\gamma}(X)$ from Lemma 1.

By Lemma 1, the following dynamics also converge to the gradient descent dynamics for $f_{\gamma}$.

[TABLE]

Remark 6** (Exponentially fast convergence).**

Theorem 5 relies upon the existence of an ergodic, invariant measure $\rho^{\infty}_{1}(y;x)$. Convergence to such a measure is exponentially fast if the underlying function, $H(x,y;\gamma)$, is convex in $y$; in our case, this happens if $\nabla^{2}f(x)+\gamma^{-1}I$ is positive definite for all $x$.

4.3. Elastic-SGD as local entropy

Elastic-SGD is an algorithm introduced by Zhang et al., (2015) for distributed training of deep neural networks and aims to minimize the communication overhead between a set of workers that together optimize replicated copies of the original function $f(x)$. For $n_{p}>0$ distinct workers, Let $y=(y_{1},\dots,y_{n_{p}})$, and let

[TABLE]

be the average value. Define the modified loss function

[TABLE]

The formulation in (20) lends itself easily to a distributed optimization where worker $i$ performs the following update, which depends only on the average of the other workers

[TABLE]

This algorithm was later shown to be connected to the local entropy loss $f_{\gamma}(x)$ by Baldassi et al., 2016a using arguments from replica theory. The results from the previous section can be modified to prove that the dynamics above converges to gradient descent of the local entropy. Consider the system

[TABLE]

along with (21). As in Sec. 4.2, define $H(x,y;\gamma)=\sum_{i=1}^{n_{p}}f(y_{i})+\frac{1}{2\gamma}\left|x-\bar{y}\right|^{2}$. The invariant measure for the $y(s)$ dynamics corresponding to the $\epsilon$ scaling of (21) is given by

[TABLE]

and the homogenized vector field is

[TABLE]

where $\left\langle\bar{y}\right\rangle=\int\bar{y}\ \rho^{\infty}(d\bar{y};\ X)$. By ergodicity, we can replace the average over the invariant measure with a combination of the temporal average and the average over the workers

[TABLE]

Thus the same argument as in Theorem 5 applies, and we can conclude that the dynamics also converge to gradient descent for local entropy.

Remark 7** (Variational Interpretation of the distributed algorithm).**

An alternate approach to understanding the distributed version of the algorithm, which we hope to study in the future, is modeling it as the gradient flow in the Wasserstein metric of a non-local functional. This functional is obtained by including an additional term in the Fokker-Planck functional,

[TABLE]

see Santambrogio, (2015).

4.4. Heat equation versus the viscous Hamilton-Jacobi equation

Lemma 1 showed that local entropy is the solution to the viscous Hamilton-Jacobi equation (viscous-HJ). The homogenization dynamics in (Entropy-SGD) can thus be interpreted as a way of smoothing of the original loss $f(x)$ to aid optimization. Other partial differential equations could also be used to achieve a similar effect. For instance, the following dynamics that corresponds to gradient descent on the solution of the the heat equation:

[TABLE]

Using a very similar analysis as above, the invariant distribution for the fast variables is $\rho^{\infty}(dy;\ X)=G_{\beta^{-1}\gamma}(y)$. In other words, the fast-dynamics for $y$ is decoupled from $x$ and $y(s)$ converges to a Gaussian distribution with mean zero and variance $\beta^{-1}\gamma\ I$. The homogenized dynamics is given by

[TABLE]

where $v(x,\gamma)$ is the solution of the heat equation $v_{t}=\frac{\beta^{-1}}{2}\Delta v$ with initial data $v(x,0)=f(x)$. The dynamics corresponds to a Gaussian averaging of the gradients.

Remark 8** (Comparison between local entropy and heat equation dynamics).**

The extra term in the $y(s)$ dynamics for (Entropy-SGD-2) with respect to (HEAT) is exactly the distinction between the smoothing performed by local entropy and the smoothing performed by heat equation. The latter is common in the deep learning literature under different forms, e.g., Gulcehre et al., (2016) and Mobahi, (2016). The former version however has much better empirical performance (cf. experiments in Sec. 8 and Chaudhari et al., (2016)) as well as improved convergence in theory (cf. Theorem 12). Moreover, at a critical point, the gradient term vanishes, and the operators (HEAT) and (Entropy-SGD-2) coincide.

5. Stochastic control interpretation

The interpretation of the local entropy function $f_{\gamma}(x)$ as the solution of a viscous Hamiltonian-Jacobi equation provided by Lemma 1 allows us to interpret the gradient descent dynamics as an optimal control problem. While the interpretation does not have immediate algorithmic implications, it allows us to prove that the expected value of a minimum is improved by the local entropy dynamics compared to the stochastic gradient descent dynamics. Controlled stochastic differential equations (Fleming and Soner, , 2006; Fleming and Rishel, , 2012) are generalizations of SDEs discussed earlier in Sec. 2.3.

5.1. Stochastic control for a variant of local entropy

Consider the following controlled SDE

[TABLE]

along with $x(t)=x$. For fixed controls $\alpha$, the generator for (CSGD) is given by

[TABLE]

Define the cost functional for a given solution $x(\cdot)$ of (CSGD) corresponding to the control $\alpha(\cdot)$,

[TABLE]

Here the terminal cost is a given function $V:\mathbb{R}^{n}\to\mathbb{R}$ and we use the prototypical quadratic running cost. Define the value function

[TABLE]

to be the minimum expected cost, over admissible controls, and over paths which start at $x(t)=x$. From the definition, the value function satisfies

[TABLE]

The dynamic programming principle expresses the value function as the unique viscosity solution of a Hamilton-Jacobi-Bellman (HJB) PDE (Fleming and Soner, , 2006; Fleming and Rishel, , 2012), $u_{t}=H(\nabla u,\ D^{2}u)$ where the operator $H$ is defined by

[TABLE]

The minimum is achieved at $\alpha=p$, which results in $H(p,Q)=-\nabla f\cdot p-\frac{1}{2}|p|^{2}+\frac{\beta^{-1}}{2}\ \mathrm{tr}\ Q$. The resulting PDE for the value function is

[TABLE]

along with the terminal values $u(x,T)=V(x)$. We make note that in this case, the optimal control is equal to the gradient of the solution

[TABLE]

Remark 9** (Forward-backward equations).**

Note that the PDE (HJB) is a terminal value problem in backwards time. This equation is well-posed, and by relabelling time, we can obtain an initial value problem in forward time. The corresponding forward equation for the evolution is of a density under the dynamics (CSGD), involving the gradient of the value function. Together, these PDEs are the forward-backward system

[TABLE]

for $0\ \leq\ s\ \leq\ T$ along with the terminal and the initial data $u(x,T)=V(x)$, $\rho(x,0)=\rho_{0}(x)$.

Remark 10** (Stochastic control for (viscous-HJ)).**

We can also obtain a stochastic control interpretation for (viscous-HJ) by dropping the $\nabla f(x(s))$ term in (CSGD). Then, by a similar argument, (HJB) results in (viscous-HJ). Thus we obtain the interpretation of solutions of (viscous-HJ) as the value function of an optimal control problem, minimizing the cost function (24) but with the modified dynamics. The reason for including the $\nabla f$ term in (CSGD) is that it allows us to prove the comparison principle in the next subsection.

Remark 11** (Mean field games).**

More generally, we can consider stochastic control problems where the control and the cost depend on the density of other players. In the case of “small” players, this can lead to mean field games (Lasry and Lions, , 2007; Huang et al., , 2006). In the distributed setting is it natural to consider couplings (control and costs) which depend on the density (or moments of the density). In this context, it may be possible to use mean field game theory to prove an analogue of Theorem 12 and obtain an estimate of the improvement in the terminal reward.

5.2. Improvement in the value function

The optimal control interpretation above allows us to provide the following theorem for the improvement in the value function obtained by the dynamics (CSGD) using the optimal control (25), compared to stochastic gradient descent.

Theorem 12**.**

Let $x_{\mathrm{csgd}}(s)$ and $x_{\mathrm{sgd}}(s)$ be solutions of (CSGD) and (SGD), respectively, with the same initial data $x_{\mathrm{csgd}}(0)=x_{\mathrm{sgd}}(0)=x_{0}$. Fix a time $t\geq 0$ and a terminal function, $V(x)$. Then

[TABLE]

The optimal control is given by $\alpha(x,t)=\nabla u(x,t)$, where $u(x,t)$ is the solution of (HJB) along with terminal data $u(x,T)=V(x)$.

Proof.

Consider (SGD) and let $\mathcal{L}=-\nabla f(x)\cdot\nabla+\frac{1}{2}\ \Delta$ to be the generator, as in (8). The expected value function $v(x,t)=\mathbb{E}\left[V(x(t))\right]$ is the solution of the PDE $v_{t}=\mathcal{L}v$. Let $u(x,t)$ be the solution of (HJB). From the PDE (HJB), we have

[TABLE]

We also have $u(x,T)=v(x,T)=V(x)$. Use the maximum principle (Evans, , 1998) to conclude

[TABLE]

Use the definition to get

[TABLE]

the expectation is taken over the paths of (SGD). Similarly,

[TABLE]

over paths of (CSGD) with $x(t)=x$. Here $\alpha(\cdot)$ is the optimal control. The interpretation along with (27) gives the desired result. Note that the second term on the right hand side in Theorem 12 above can also be written as

[TABLE]

6. Regularization, widening and semi-concavity

The PDE approaches discussed in Sec. 3 result in a smoother loss function. We exploit the fact that our regularized loss function is the solution of the viscous HJ equation. This PDE allows us to apply standard semi-concavity estimates from PDE theory (Evans, , 1998; Cannarsa and Sinestrari, , 2004) and quantify the amount of smoothing. Note that these estimates do not depend on the coefficient of viscosity so they apply to the HJ equation as well. Indeed, semi-concavity estimates apply to solutions of PDEs in a more general setting which includes (viscous-HJ) and (HJ).

We begin with special examples which illustrate the widening of local minima. In particular, we study the limiting case of the HJ equation, which has the property that local minima are preserved for short times. We then prove semiconcavity and establish a more sensitive form of semiconcavity using the harmonic mean. Finally, we explore the relationship between wider robust minima and the local entropy.

6.1. Widening for solutions of the Hamilton-Jacobi PDE

A function $f(x)$ is semi-concave with a constant $C$ if $f(x)-\frac{C}{2}\ \left|x\right|^{2}$ is concave, refer to Cannarsa and Sinestrari, (2004). Semiconcavity is a way of measuring the width of a local minimum: when a function is semiconcave with constant $C$, at a local minimum, no eigenvalue of the Hessian can be larger than $C$. The semiconcavity estimates below establish that near high curvature local minima widen faster than flat ones, when we evolve $f$ by (viscous-HJ).

It is illustrating to consider the following example.

Example 13**.**

Let $u(x,t)$ be the solution of  (viscous-HJ) with $u(x,0)=c\left|x\right|^{2}/2$. Then

[TABLE]

In particular, in the non-viscous case, $u(x,t)=\frac{\left|x\right|^{2}}{2(t+c^{-1})}$.

Proof.

This follows from taking derivatives:

[TABLE]

In the previous example, the semiconcavity constant of $u(x,t)$ is $C(t)=1/(c^{-1}+t)$ where $c$ measures the curvature of the initial data. This shows that for very large $c$, the improvement is very fast, for example $c=10^{8}$ leads to $C(t=.01)\approx 10$, etc. In the sequel, we show how this result applies in general for both the viscous and the non-viscous cases. In this respect, for short times, solutions of (viscous-HJ) and (HJ) widen faster than solutions of the heat equation.

Example 14** (Slower rate for the heat equation).**

In this example we show that the heat equation can result in a very slow improvement of the semiconcavity. Let $v_{1}(x)$ be the first eigenfunction of the heat equation, with corresponding eigenvalue $\lambda_{1}$, (for example $v_{1}(x)=\sin(x)$). In many cases $\lambda_{1}$ is of order 1. Then, since $v(x,t)=\exp(-\lambda_{1}t)\ v_{1}(x)$ is a solution of the heat equation, the seminconcavity constant for $v(x,t)$ is $C(t)=c_{0}\exp(-\lambda_{1}t)$, where $c_{0}$ is the constant for $v_{1}(x)$. For short time, $C(t)\approx c_{0}(1-\lambda_{1}t)$ which corresponds to a much slower decrease in $C(t)$.

Another illustrative example is to study the widening of convex regions for the non-viscous Hamilton-Jacobi equation (HJ). As can be seen in Figure 1, solutions of the Hamilton-Jacobi have additional remarkable properties with respect to local minima. In particular, for small values of $t$ (which depend on the derivatives of the initial function) local minima are preserved. To motivate these properties, consider the following examples. We begin with the Burgers equation.

Example 15** (Burgers equation in one dimension).**

There is a well-known connection between Hamilton-Jacobi equations in one dimension and the Burgers equation (Evans, , 1998) which is the prototypical one-dimensional conservation law. Consider the solution $u(x,t)$ of the non-viscous HJ equation (HJ) with initial value $f(x)$ where $x\in\mathbb{R}$. Let $p(x,t)=u_{x}(x,t)$. Then $p$ satisfies the Burgers equation

[TABLE]

We can solve Burgers equation using the method of characteristics, for a short time, depending on the regularity of $f(x)$. Set $p(x,t)=u_{x}(x,t)$. For short time, we can also use the fact that $p(x,t)$ is the fixed point of (16), which leads to

[TABLE]

In particular, the solutions of this equation remain classical until the characteristics intersect, and this time, $t^{*}$, is determined by

[TABLE]

which recovers the convexity condition of Lemma 3.

Example 16** (Widening of convex regions in one dimension).**

In this example, we show that the widening of convex regions for one dimensional displayed in Figure 1 occurs in general. Let $x$ be a local minimum of the smooth function $f(x)$, and let $0\leq t<t^{*}$, where the critical $t^{*}$ is defined by (29). Define the convexity interval as:

[TABLE]

Let $I(x,0)=[x_{0},\ x_{1}]$, so that $f^{\prime\prime}(x)\geq 0$ on $[x_{0},\ x_{1}]$ and $f^{\prime\prime}(x_{0})=f^{\prime\prime}(x_{1})=0$. Then $I(x,t)$ contains $I(x,0)$ for $0\leq t\leq t^{*}$. Differentiating (28), and the solving for $p_{x}(x,t)$ leads to

[TABLE]

Since $t\leq t^{*}$, the last equation implies that $p_{x}$ has the same sign as $f^{\prime\prime}(x-tp)$. Recall now that $x_{1}$ is an inflection point of $f(x)$. Choose $x$ such that $x_{1}=x-tp(x,t)$, in other words, $u_{xx}(x,t)=0$. Note that $x>x_{1}$, since $p(x,t)>0$. Thus the inflection point moves to the right of $x_{1}$. Similarly the inflection point $x_{0}$ moves to the left. So we have established that the interval is larger.

Remark 17**.**

In the higher dimensional case, well-known properties of the inf-convolution can be used to establish the following facts. If $f(x)$ is convex, $u(x,t)$ given by (HL) is also convex; thus, any minimizer of $f$ also minimizes $u(x,t)$. Moreover, for bounded $f(x)$, using the fact that $\left|y^{*}(x)-x\right|=\mathcal{O}(\sqrt{t})$, one can obtain a local version of the same result, viz., for short times $t$, a local minimum persists. Similar to the example in Figure 1, local minima in the solution vanish for longer times.

6.2. Semiconcavity estimates

In this section we establish semiconcavity estimates. These are well known, and are included for the convenience of the reader.

Lemma 18**.**

Suppose $u(x,t)$ is the solution of (viscous-HJ), and let $\beta^{-1}\geq 0$. If

[TABLE]

then

[TABLE]

Proof.

1. First consider the nonviscous case $\beta^{-1}=0$. Define $g(x,y)=f(y)+\frac{1}{2t}\ \left|x-y\right|^{2}$ to be a function of $x$ parametrized by $y$. Each $g(x,y)$, considered as a function of $x$ is semi-concave with constant $1/t$. Then $u(x,t)$, the solution of (HJ), which is given by the Hopf-Lax formula (HL) is expressed as the infimum of such functions. Thus $u(x,t)$ is also semi-concave with the same constant, $1/t$.

2. Next, consider the viscous case $\beta^{-1}>0$. Let $w=u_{x_{k}x_{k}}$, differentiating (viscous-HJ) twice gives

[TABLE]

Using $u_{x_{k}x_{k}}^{2}=w^{2}$, we have

[TABLE]

Note that $w(x,0)\leq C_{k}$ and the function

[TABLE]

is a spatially independent super-solution of the preceding equation with $w(x,0)\leq\phi(x,0)$. By the comparison principle applied to $w$ and $\phi$, we have

[TABLE]

for all $x$ and for all $t\geq 0$, which gives the first result. Now set $v=\Delta u$ and sum (30) over all $k$ to obtain

[TABLE]

Thus

[TABLE]

by the Cauchy-Schwartz inequality. Now apply the comparison principle to $v^{\prime}(x,t)=\left(C_{\mathrm{Lap}}^{-1}+t/n\right)^{-1}$ to obtain the second result.

6.3. Estimates on the Harmonic mean of the spectrum

Our semi-concavity estimates in Sec. 6.2 gave bounds on $\sup_{x}\ u_{x_{k}x_{k}}(x,t)$ and the Laplacian $\sup_{x}\ \Delta u(x,t)$. In this section, we extend the above estimates to characterize the eigenvalue spectrum more precisely. Our approach is motivated by experimental evidence in Chaudhari et al., (2016, Figure 1) and Sagun et al., (2016). The authors found that the Hessian of a typical deep network at a local minimum discovered by SGD has a very large proportion of its eigenvalues that are close to zero. This trend was observed for a variety of neural network architectures, datasets and dimensionality. For such a Hessian, instead of bounding the largest eigenvalue or the trace (which is also the Laplacian), we obtain estimates of the harmonic mean of the eigenvalues. This effectively discounts large eigenvalues and we get an improved estimate of the ones close to zero, that dictate the width of a minimum. The harmonic mean of a vector $x\in\mathbb{R}^{n}$ is

[TABLE]

The harmonic mean is more sensitive to small values as seen by

[TABLE]

which does not hold for the arithmetic mean. To give an example, the eigenspectrum of a typical network (Chaudhari et al., , 2016, Figure 1) has an arithmetic mean of $0.0029$ while its harmonic mean is $\approx 10^{-10}$ which better reflects the large proportion of almost-zero eigenvalues.

Lemma 19**.**

If $\Lambda$ is the vector of eigenvalues of the symmetric positive definite matrix $A\in\mathbb{R}^{n\times n}$ and $D=(a_{11},\ldots,a_{nn})$ is its diagonal vector,

[TABLE]

Proof.

Majorization is a partial order on vectors with non-negative components (Marshall et al., , 1979, Chap. 3) and for two vectors $x,y\in\mathbb{R}^{n}$, we say that $x$ is majorized by $y$ and write $x\preceq y$ if $x=Sy$ for some doubly stochastic matrix $S$. The diagonal of a symmetric, positive semi-definite matrix is majorized by the vector of its eigenvalues (Schur, , 1923). For any Schur-concave function $f(x)$, if $x\preceq y$ we have

[TABLE]

from Marshall et al., (1979, Chap. 9). The harmonic mean is also a Schur-concave function which can be checked using the Schur-Ostrowski criterion

[TABLE]

and we therefore have

[TABLE]

Lemma 20**.**

If $u(x,t)$ is a solution of (viscous-HJ) and $C\in\mathbb{R}^{n}$ is such that $\sup_{x}\ u_{x_{k}x_{k}}(x,t)\leq C_{k}$ for all $k\leq n$, then at a local minimum $x^{*}$ of $u(x,t)$,

[TABLE]

where $\Lambda$ is the vector of eigenvalues of $\nabla^{2}u(x^{*},t)$.

Proof.

For a local minimum $x^{*}$, we have $\mathrm{HM}(L)\leq\mathrm{HM}(D)$ from the previous lemma; here $D$ is the diagonal of the Hessian $\nabla^{2}u(x_{\mathrm{min}},t)$. From Lemma 18 we have

[TABLE]

The result follows by observing that

[TABLE]

7. Algorithmic details

In this section, we compare and contrast the various algorithms in this paper and provide implementation details that are used for the empirical validation are provided in Sec. 8.

7.1. Discretization of (Entropy-SGD)

The time derivative in the equation for $y(s)$ in (Entropy-SGD) is discretized using the Euler-Maruyama method with time step (learning rate) $\eta_{y}$ as

[TABLE]

where $\epsilon_{k}$ are zero-mean, unit variance Gaussian random variables.

In practice, we do not have access to the full gradient $\nabla f(x)$ for a deep network and instead, only have the gradient over a mini-batch $\nabla f_{\mathrm{mb}}(x)$ (cf. Sec. 2.2). In this case, there are two potential sources of noise: the coefficient $\beta^{-1}_{\mathrm{mb}}$ arising from the mini-batch gradient, and any amount of extrinsic noise, with coefficient $\beta^{-1}_{\mathrm{ex}}$. Combining these two sources leads to the equation

[TABLE]

which corresponds to (31) with

[TABLE]

We initialize $y^{k}=x^{k}$ if $k/L$ is an integer. The number of time steps taken for the $y$ variable is set to be $L$. This corresponds to $\epsilon=1/L$ in (17) and to $T=1$ in (19). This results in a discretization of the equation for $x$ in (Entropy-SGD) given by

[TABLE]

Since we may be far from the limit $\epsilon\to 0$, $T\to\infty$, instead of using the linear averaging for $\left\langle y\right\rangle$ in (19), we will use a forward looking average. Such an averaging gives more weight to the later steps which is beneficial for large values of $\gamma$ when the invariant measure $\rho^{\infty}(y;\ x)$ in (19) does not converge quickly. The two algorithms we describe will differ in the choices of $L$ and $\beta^{-1}_{\mathrm{ex}}$ and the definition of $\left\langle y\right\rangle$. For Entropy-SGD, we set $L=20$, $\beta^{-1}_{\mathrm{ex}}=10^{8}$, and set $\left\langle y\right\rangle$ to be

[TABLE]

$\left\langle y\right\rangle^{k}$ is re-initialized to $x^{k}$ if $k/L$ is an integer. The parameter $\alpha$ is used to perform an exponential averaging of the iterates $y^{k}$. We set $\alpha=0.75$ for experiments in Sec. 8 using Entropy-SGD. This is equivalent to the original implementation of Chaudhari et al., (2016).

The step-size for the $y^{k}$ dynamics is fixed to $\eta_{y}=0.1$. It is a common practice to perform step-size annealing while optimizing deep networks (cf. Sec. 2.5), i.e., we set $\eta=0.1$ and reduce it by a factor of $5$ after every $3$ epochs (iterations on the entire dataset).

7.2. Non-viscous Hamilton-Jacobi (HJ)

For the other algorithm which we denote as HJ, we set $L=5$ and $\beta^{-1}_{\mathrm{ex}}=0$, and set $\left\langle y\right\rangle^{k}=y^{k}$ in (32) and (33), i.e., no averaging is performed and we simply take the last iterate of the $y^{k}$ updates.

Remark 21**.**

We can also construct an equivalent system of updates for HJ by discretizing (Entropy-SGD-2) and again setting $\beta^{-1}_{\mathrm{ex}}=0$; this exploits the two different ways of computing the gradient in Lemma 1 and gives

[TABLE]

and we initialize $y^{k}=0$ if $k/L$ is an integer.

7.3. Heat equation:

We perform the update corresponding to (23)

[TABLE]

where $\epsilon^{i}$ are Gaussian random variables with zero mean and variance $\gamma\ I$. Note that we have implemented the convolution in (23) as an average over Gaussian perturbations of the parameters $x$. We again set $L=20$ for our experiments in Sec. 8.

7.4. Momentum

Momentum is a popular technique for accelerating the convergence of SGD (Nesterov, , 1983). We employ this in our experiments by maintaining an auxiliary variable $z^{k}$ which intuitively corresponds to the velocity of $x^{k}$. The $x^{k}$ update in (33) is modified to be

[TABLE]

if $k/L$ is an integer; $x^{k}$ and $z^{k}$ are left unchanged otherwise. The other updates for $x^{k}$ in (34) and (35) are modified similarly. The momentum parameter $\delta=0.9$ is fixed for all experiments.

7.5. Tuning $\gamma$

Theorem 12 gives a rigorous justification to the technique of ‘scoping” where we set $\gamma$ to a large value and reduce it as training progresses, indeed, as Figure 1 shows, a larger $\gamma$ leads to a smoother loss. Scoping reduces the smoothing effect and as $\gamma\to 0$, we have $f_{\gamma}(x)\to f(x)$. For the updates in Sec. 7.1, 7.2 and 7.3, we update $\gamma$ every time $k/L$ is an integer using

[TABLE]

we pick $\gamma_{0}\in[10^{-4},\ 10^{-1}]$ so as to obtain the best validation error while $\gamma_{1}$ is fixed to $10^{-3}$ for all experiments. We have obtained similar results by normalizing $\gamma_{0}$ by the dimensionality of $x$, i.e., if we set $\gamma_{0}:=\gamma_{0}/n$ this parameter becomes insensitive to different datasets and neural networks.

8. Empirical validation

We now discuss experimental results on deep neural networks that demonstrate that the PDE methods considered in this article achieve good regularization, aid optimization and lead to improved classification on modern datasets.

8.1. Setup for deep networks

In machine learning, one typically splits the data into three sets: (i) the training set $D$, used to compute the optimization objective (empirical loss $f$), (ii) the validation set, used to tune parameters of the optimization process which are not a part of the function $f$, e.g., mini-batch size, step-size, number of iterations etc., and, (iii) a test set, used to quantify generalization as measured by the empirical loss on previously-unseen (sequestered) data. However, for the benchmark datasets considered here, it is customary in the literature to not use a separate test set. Instead, one typically reports test error on the validation set itself; for the sake of enabling direct comparison of numerical values, we follow the same practice here.

We use two standard computer vision datasets for the task of image classification in our experiments. The MNIST dataset (LeCun et al., , 1998) contains $70,000$ gray-scale images of size $28\times 28$, each image portraying a hand-written digit between [math] to $9$. This dataset is split into $60,000$ training images and $10,000$ images in the validation set. The CIFAR-10 dataset (Krizhevsky, , 2009) contains $60,000$ RGB images of size $32\times 32$ of $10$ objects (aeroplane, automobile, bird, cat, deer, dog, frog, horse, ship and truck). The images are split as $50,000$ training images and $10,000$ validation images. Figure 3 and Figure 3 show a few example images from these datasets.

8.1.1. Pre-processing

We do not perform any pre-processing for the MNIST dataset. For CIFAR-10, we perform a global contrast normalization (Coates et al., , 2010) followed by a ZCA whitening transform (sometimes also called the “Mahalanobis transformation”) (Krizhevsky, , 2009). We use ZCA as opposed to principal component analysis (PCA) for data whitening because the former preserves visual characteristics of natural images unlike the latter; this is beneficial for convolutional neural networks. We do not perform any data augmentation, i.e., transformation of input images by mirror-flips, translations, or other transformations.

8.1.2. Loss function

We defined the zero-one loss $f_{i}(x)$ for a binary classification problem in Sec. 2.1. For classification tasks with $K$ classes, it is common practice to construct a network with $K$ distinct outputs. The output $y(x;\ \xi^{i})$ is thus a vector of length $K$ that is further normalized to sum up to $1$ using the $\mathrm{softmax}$ operation defined in Sec. 8.2. The cross-entropy loss typically used for optimizing deep networks is defined as:

[TABLE]

We use the cross-entropy loss for all our experiments.

8.1.3. Training procedure

We run the following algorithms for each neural network, their algorithmic details are given in Sec. 7. The ordering below corresponds to the ordering in the figures.

• •

Entropy-SGD: the algorithm of Chaudhari et al., (2016) described in Sec. 7.1,

• •

HEAT: smoothing by the heat equation described in Sec. 7.3,

• •

HJ: updates for the non-viscous Hamilton-Jacobi equation described in Sec. 7.2,

• •

SGD: stochastic gradient descent given by (7),

All the above algorithms are stochastic in nature; in particular, they sample a mini-batch of images randomly from the training set at each iteration before computing the gradient. A significant difference between Entropy-SGD and HJ is that the former adds additional noise, while the latter only uses the intrinsic noise from the mini-batch gradients. We therefore report and plot the mean and standard deviation across $6$ distinct random seeds for each experiment. The mini-batch size ($\left|i_{t}\right|$ in (7)) is fixed at $128$. An “epoch” is defined to be one pass over the entire dataset, for instance, it consists of $\lceil 60,000/128\rceil=469$ iterations of (7) for MNIST with each mini-batch consisting of $128$ different images.

We use SGD as a baseline for comparing the performance of the above algorithms in all our experiments.

Remark 22** (Number of gradient evaluations).**

A deep network evaluates the average gradient over an entire mini-batch of $128$ images in one single back-propagation operation (Rumelhart et al., , 1988). We define the number of back-props per weight update by $L$. SGD uses one back-prop per update (i.e., $L=1$) while we use $L=20$ back-props for Entropy-SGD (cf. Sec. 7.1) and the heat equation (cf. Sec. 7.3) and $L=5$ back-props for the non-viscous HJ equation (cf. Sec. 7.2). For each of the algorithms, we plot the training loss or the validation error against the number of “effective epochs”, i.e., the number of epochs is multiplied by $L$ which gives us a uniform scale to measure the computational efficiency of various algorithms. This is a direct measure of the wall-clock time and is agnostic to the underlying hardware and software implementation.

8.2. MNIST

8.2.1. Fully-connected network

Consider a “fully-connected” network on MNIST

[TABLE]

The $\mathrm{input}$ layer reshapes each MNIST image as a vector of size $784$. The notation $\mathrm{fc}_{d}$ denotes a “fully connected” dense matrix $x^{k}$ in (4) with $d$ output dimensions, followed by the ReLU nonlinearity (cf. Sec. 2.1) and an operation known as batch normalization (Ioffe and Szegedy, , 2015) which whitens the output of each layer by subtracting the mean across a mini-batch and dividing by the standard deviation. The notation $\mathrm{drop}_{p}$ denotes the dropout layer (Srivastava et al., , 2014) which randomly sets a fraction $p$ of the weights to zero at each iteration; this is a popular regularization technique in deep learning, see Kingma et al., (2015); Achille and Soatto, (2016) for a Bayesian perspective. For $10$ classes in the classification problem, we create an output vector of length $10$ in the last $\mathrm{fc}$ layer followed by $\mathrm{softmax}$ which picks the largest element of this vector, which is interpreted as the output (or “prediction”) by this network. The $\mathrm{mnistfc}$ network has $n=1.86$ million parameters.

As Figure 4(a) and Table 1 show, the final validation error for all algorithms is quite similar. The convergence rate varies, in particular, Entropy-SGD converges fastest in this case. Note that $\mathrm{mnistfc}$ is a small network and the difference in the performance of the above algorithms, e.g., $1.08\ \%$ for Entropy-SGD versus $1.17\ \%$ for HJ, is minor.

8.2.2. LeNet

Our second network for MNIST is a convolutional neural network (CNN) denoted as follows:

[TABLE]

The notation $\mathrm{conv}_{c,k,m}$ denotes a 2D convolutional layer with $c$ output channels, each of which is the sum of a channel-wise convolution operation on the input using a learnable kernel of size $k\times k$. It further adds ReLU nonlinearity, batch normalization and an operation known as max-pooling which down-samples the image by picking the maximum over a patch of size $m\times m$. Convolutional networks typically perform much better than fully-connected ones on image classification task despite having fewer parameters, $\mathrm{LeNet}$ has only $n=131,220$ parameters.

The results for $\mathrm{LeNet}$ are described in Figure 4(b) and Table 1. This is a convolutional neural network and performs better than $\mathrm{mnistfc}$ on the same dataset. The final validation error is very similar for all algorithms at $0.50\ \%$ with the exception of the heat equation which only reaches $0.59\ \%$. We can also see that the other algorithms converge quickly, in about half the number of effective epochs as that of $\mathrm{mnistfc}$.

8.3. CIFAR

The CIFAR-10 dataset is more complex than MNIST and fully-connected networks typically perform poorly. We will hence employ a convolutional network for this dataset. We use the All-CNN-C architecture introduced by Springenberg et al., (2014) and add batch-normalization:

[TABLE]

where

[TABLE]

The final convolutional layer in the above $\mathrm{block}$ denoted by $\mathrm{conv}^{*}$ is different from others; while they perform convolution at every pixel otherwise known as a “stride” of $1$ pixel, $\mathrm{conv}^{*}$ on the other hand uses a stride of $2$; this results in the image being downsampled by a factor of two. Note that $\mathrm{conv}_{c,k,m}$ with $m=1$ does not result in any downsampling. Max-pooling usually results in a drastic reduction of the image size and by replacing it with a strided convolution, All-CNN achieves improved performance with much fewer parameters than many other networks on CIFAR-10. The final layer denoted as mean-pool takes an input of size $10\times 8\times 8$ and computes the spatial average to obtain an output vector of size $10$. This network has $n=1.67$ million parameters.

Figure 5(a) and 5(b) show the training loss and validation error for the All-CNN network on the CIFAR-10 dataset. The Hamilton-Jacobi equation (HJ) obtains a validation error of $7.89\ \%$ in $145$ epochs and thus performs best among the algorithms tested here; it also has the lowest training cross-entropy loss of $0.046$. Note that both HJ and Entropy-SGD converge faster than SGD. The heat equation performs again poorly on this dataset and has a much higher validation error than others ($9.04\%$). This is consistent with our discussion in Sec. 4.4 which suggests that the viscous or non-viscous HJ equations result in better smoothing than the heat equation.

9. Discussion

Our results apply nonlinear PDEs, stochastic optimal control, and stochastic homogenization to the analysis of two recent and effective algorithms for the optimization of neural networks. Our analysis also contributes to devising improved algorithms.

We replaced the standard stochastic gradient descent (SGD) algorithm for the function $f(x)$, with SGD on the two variable function $H(x,y;\gamma)=f(y)+\gamma^{-1}|x-y|^{2}/2$, along with a small parameter $\epsilon$ (Entropy-SGD). Using the new system, in the limit $\epsilon\to 0$, we can provably, and quantitatively, improve the expected value of the original loss function. The effectiveness of our algorithm comes from the connection, via homogenization, of system of SDEs to the gradient of the function $u(x,t)$ which is the solution of the (viscous-HJ) PDE with initial data $f(x)$. The function $H$ is more convex in $y$ that the original function $f$. The convexity of $H$ in $y$ is related to the exponentially fast convergence of the dynamics, a connection explained by the celebrated gradient flow interpretation of Fokker-Planck equation Jordan et al., (1998). Ergodicity of the dynamics for (Entropy-SGD) makes clear that the local entropy SGD is equivalent to the influential distributed algorithm, Elastic-SGD Zhang et al., (2015). These insights may lead to improved distributed algorithms, which are presently of great interest.

On a practical level, the large number of hyper-parameters involved in training deep neural networks is very costly, in terms of both human effort and computer time. Our analysis lead to better understanding of the parameters involved in the algorithms, and provides an insight into the choice of hyper-parameters for these methods. In particular (i) the parameter $\gamma$ is now understood as the time $t$, in the PDE (viscous-HJ) (ii) scoping of $\gamma$, which was seen as a heuristic, is now rigorously justified, (iii) we set the extrinsic noise parameter $\beta^{-1}_{\mathrm{ex}}$ to zero, resulting in a simpler, more effective algorithm, and (iv) we now understand that below a critical value of $\gamma$ (which is determined by the requirement that $H$ be convex in $y$), the convergence of the dynamics in (Entropy-SGD) to the invariant measure is exponentially fast.

Conceptually, while simulated annealing and related methods work by modulating the level of noise in the dynamics, our algorithm works by modulating the smoothness of the underlying loss function. Moreover, while most algorithms used in deep learning derive their motivation from the literature on convex optimization, the algorithms we have presented here are specialized to non-convex loss functions and have been shown to perform well on these problems both in theory and in practice.

10. Acknowledgments

AO is supported by a grant from the Simons Foundation (395980); PC and SS by ONR N000141712072, AFOSR FA95501510229 and ARO W911NF151056466731CS; SO by ONR N000141410683, N000141210838, N000141712162 and DOE DE-SC00183838. AO would like to thank the hospitality of the UCLA mathematics department where this work was completed.