Title: TOPOLOGY AND GEOMETRY OF HALF-RECTIFIED NETWORK OPTIMIZATION

Abstract: The loss surface of deep neural networks has recently attracted interest in the optimization and machine learning communities as a prime example of highdimensional non-convex problem. Some insights were recently gained using spin glass models and mean-field approximations, but at the expense of strongly simplifying the nonlinear nature of the model. In this work, we do not make any such assumption and study conditions on the data distribution and model architecture that prevent the existence of bad local minima. Our theoretical work quantifies and formalizes two important folklore facts: (i) the landscape of deep linear networks has a radically different topology from that of deep half-rectified ones, and (ii) that the energy landscape in the non-linear case is fundamentally controlled by the interplay between the smoothness of the data distribution and model over-parametrization. Our main theoretical contribution is to prove that half-rectified single layer networks are asymptotically connected, and we provide explicit bounds that reveal the aforementioned interplay. The conditioning of gradient descent is the next challenge we address. We study this question through the geometry of the level sets, and we introduce an algorithm to efficiently estimate the regularity of such sets on large-scale networks. Our empirical results show that these level sets remain connected throughout all the learning phase, suggesting a near convex behavior, but they become exponentially more curvy as the energy level decays, in accordance to what is observed in practice with very low curvature attractors.

1. INTRODUCTION:

Optimization is a critical component in deep learning, governing its success in different areas of computer vision, speech processing and natural language processing. The prevalent optimization strategy is Stochastic Gradient Descent, invented by Robbins and Munro in the 50s. The empirical performance of SGD on these models is better than one could expect in generic, arbitrary non-convex loss surfaces, often aided by modifications yielding significant speedups Duchi et al. (2011) ; Hinton et al. (2012) ; Ioffe & Szegedy (2015) ; Kingma & Ba (2014) . This raises a number of theoretical questions as to why neural network optimization does not suffer in practice from poor local minima.

The loss surface of deep neural networks has recently attracted interest in the optimization and machine learning communities as a paradigmatic example of a hard, high-dimensional, non-convex problem. Recent work has explored models from statistical physics such as spin glasses Choromanska et al. (2015) , in order to understand the macroscopic properties of the system, but at the expense of strongly simplifying the nonlinear nature of the model. Other authors have advocated that the real danger in high-dimensional setups are saddle points rather than poor local minima Dauphin et al. (2014) , although recent results rigorously establish that gradient descent does not get stuck on saddle points Lee et al. (2016) but merely slowed down. Other notable recent contributions are Kawaguchi (2016) , which further develops the spin-glass connection from Choromanska et al. (2015) and resolves the linear case by showing that no poor local minima exist; Sagun et al. (2014) which also discusses the impact of stochastic vs plain gradient, Soudry & Carmon (2016) , that studies Empirical Risk Minimization for piecewise multilayer neural networks under overparametrization (which needs to grow with the amount of available data), and Goodfellow et al. (2014) , which provided insightful intuitions on the loss surface of large deep learning models and partly motivated our work. Additionally, the work Safran & Shamir (2015) studies some topological properties of homogeneous nonlinear networks and shows how overparametrization acts upon these properties, and the pioneering Shamir (2016) studied the distribution-specific hardness of optimizing non-convex objectives. Lastly, several papers submitted concurrently and independently of this one deserve note, particularly Swirszcz et al. (2016) which analyzes the explicit criteria under which sigmoid-based neural networks become trapped by poor local minima, as well as Tian (2017) , which offers a complementary study of two layer ReLU based networks, and their learning dynamics.

In this work, we do not make any linearity assumption and study conditions on the data distribution and model architecture that prevent the existence of bad local minima. The loss surface F (θ) of a given model can be expressed in terms of its level sets Ω λ , which contain for each energy level λ all parameters θ yielding a loss smaller or equal than λ. A first question we address concerns the topology of these level sets, i.e. under which conditions they are connected. Connected level sets imply that one can always find a descent direction at each energy level, and therefore that no poor local minima can exist. In absence of nonlinearities, deep (linear) networks have connected level sets Kawaguchi (2016) . We first generalize this result to include ridge regression (in the two layer case) and provide an alternative, more direct proof of the general case. We then move to the half-rectified case and show that the topology is intrinsically different and clearly dependent on the interplay between data distribution and model architecture. Our main theoretical contribution is to prove that half-rectified single layer networks are asymptotically connected, and we provide explicit bounds that reveal the aforementioned interplay.

Beyond the question of whether the loss contains poor local minima or not, the immediate follow-up question that determines the convergence of algorithms in practice is the local conditioning of the loss surface. It is thus related not to the topology but to the shape or geometry of the level sets. As the energy level decays, one expects the level sets to exhibit more complex irregular structures, which correspond to regions where F (θ) has small curvature. In order to verify this intuition, we introduce an efficient algorithm to estimate the geometric regularity of these level sets by approximating geodesics of each level set starting at two random boundary points. Our algorithm uses dynamic programming and can be efficiently deployed to study mid-scale CNN architectures on MNIST, CIFAR-10 and RNN models on Penn Treebank next word prediction. Our empirical results show that these models have a nearly convex behavior up until their lowest test errors, with a single connected component that becomes more elongated as the energy decays. The rest of the paper is structured as follows. Section 2 presents our theoretical results on the topological connectedness of multilayer networks. Section 3 presents our path discovery algorithm and Section 4 covers the numerical experiments.

2. TOPOLOGY OF LEVEL SETS:

Let P be a probability measure on a product space X × Y, where we assume X and Y are Euclidean vector spaces for simplicity. Let {(x i , y i )} i be an iid sample of size L drawn from P defining the training set. We consider the classic empirical risk minimization of the form

EQUATION

where Φ(x; θ) encapsulates the feature representation that uses parameters θ ∈ R S and R(θ) is a regularization term. In a deep neural network, θ contains the weights and biases used in all layers. For convenience, in our analysis we will also use the oracle risk minimization:

EQUATION

Our setup considers the case where R consists on either 1 or 2 norms, as we shall describe below. They correspond to well-known sparse and ridge regularization respectively.

2.1. POOR LOCAL MINIMA CHARACTERIZATION FROM TOPOLOGICAL CONNECTEDNESS:

We define the level set of F (θ) as

EQUATION

The first question we study is the structure of critical points of F e (θ) and F o (θ) when Φ is a multilayer neural network. For simplicity, we consider first a strict notion of local minima: θ ∈ R S is a strict local minima of F if there is > 0 with F (θ ) > F (θ) for all θ ∈ B(θ, ) and θ = θ.

In particular, we are interested to know whether F e has local minima which are not global minima. This question is answered by knowing whether Ω F (λ) is connected at each energy level λ: Proposition 2.1. If Ω F (λ) is connected for all λ then every local minima of F (θ) is a global minima.

Strict local minima implies that ∇F (θ) = 0 and HF (θ) 0, but avoids degenerate cases where F is constant along a manifold intersecting θ. In that scenario, if U θ denotes that manifold, our reasoning immediately implies that if Ω F (λ) are connected, then for all > 0 there exists θ with dist(θ , U θ ) ≤ and F (θ ) < F (θ). In other words, some element at the boundary of U θ must be a saddle point. A stronger property that eliminates the risk of gradient descent getting stuck at U θ is that all elements at the boundary of U θ are saddle points. This can be guaranteed if one can show that there exists a path connecting any θ to the lowest energy level such that F is strictly decreasing along it.

Such degenerate cases arise in deep linear networks in absence of regularization. If θ = (W 1 , . . . , W K ) denotes any parameter value, with N 1 , . . . N K denoting the hidden layer sizes, and

F k ∈ GL + N k (R)

are arbitrary elements of the general linear group of invertible N k × N k matrices with positive determinant, then

U θ = {W 1 F -1 1 , F 1 W 2 F -1 2 , . . . , F K W K ; F k ∈ GL + N k (R)} .

In particular, U θ has a Lie Group structure. In the half-rectified nonlinear case, the general linear group is replaced by the Lie group of homogeneous invertible matrices

F k = diag(α 1 , . . . , α N k ) with α j > 0.

This proposition shows that a sufficient condition to prevent the existence of poor local minima is having connected level sets, but this condition is not necessary: one can have isolated local minima lying at the same energy level. This can be the case in systems that are defined up to a discrete symmetry group, such as multilayer neural networks. However, as we shall see next, this case puts the system in a brittle position, since one needs to be able to account for all the local minima (and there can be exponentially many of them as the parameter dimensionality increases) and verify that their energy is indeed equal.

2.2. THE LINEAR CASE:

We first consider the particularly simple case where F is a multilayer network defined by

EQUATION

)

and the ridge regression R(θ) = θ 2 . This model defines a non-convex (and non-concave) loss F e (θ). When κ = 0, it has been shown in Saxe et al. (2013) and Kawaguchi (2016) that in this case, every local minima is a global minima. We provide here an alternative proof of that result that uses a somewhat simpler argument and allows for κ > 0 in the case K = 2. Proposition 2.2. Let W 1 , W 2 , . . . , W K be weight matrices of sizes n k × n k+1 , k < K, and let F e (θ), F o (θ) denote the risk minimizations using Φ as in (4). Assume that n j ≥ min(n 1 , n K ) for j = 2 . . . K -1. Then Ω Fe (λ) (and Ω Fo ) is connected for all λ and all K when κ = 0, and for κ > 0 when K = 2; and therefore there are no poor local minima in these cases. Moreover, any θ can be connected to the lowest energy level with a strictly decreasing path.

Let us highlight that this result is slightly complementary than that of Kawaguchi (2016) , Theorem 2.3. Whereas we require n j ≥ min(n 1 , n K ) for j = 2 . . . K -1 and our analysis does not inform about the order of the saddle points, we do not need full rank assumptions on Σ X nor the weights W k .

This result does also highlight a certain mismatch between the picture of having no poor local minima and generalization error. Incorporating regularization drastically changes the topology, and the fact that we are able to show connectedness only in the two-layer case with ridge regression is profound; we conjecture that extending it to deeper models requires a different regularization, perhaps using more general atomic norms Bach (2013) . But we now move our interest to the nonlinear case, which is more relevant to our purposes.

2.3. HALF-RECTIFIED NONLINEAR CASE:

We now study the setting given by

EQUATION

where ρ(z) = max(0, z). The biases can be implemented by replacing the input vector x with x = (x, 1) and by rebranding each parameter matrix as

W i = W i b i 0 1 ,

where b i contains the biases for each layer. For simplicity, we continue to use W i and x in the following.

2.3.1. NONLINEAR MODELS ARE GENERALLY DISCONNECTED:

One may wonder whether the same phenomena of global connectedness also holds in the halfrectified case. A simple motivating counterexample shows that this is not the case in general. Consider a simple setup with X ∈ R 2 drawn from a mixture of two Gaussians N -1 and N 1 , and let Y = (X -µ Z ) • Z , where Z is the (hidden) mixture component taking {1, -1} values. Let Ŷ = Φ(X; {W 1 , W 2 }) be a single-hidden layer ReLU network, with two hidden units. Let θ A be a configuration that bisects the two mixture components, and let θ B the same configuration, but swapping the bisectrices. One can verify that they can both achieve arbitrarily small risk by letting the covariance of the mixture components go to 0. However, any path that connects θ A to θ B must necessarily pass through a point in which W 1 has rank 1, which leads to an estimator with risk at least 1/2.

In fact, it is easy to see that this counter-example can be extended to any generic half-rectified architecture, if one is allowed to adversarially design a data distribution. For any given Φ(X; θ) with arbitrary architecture and current parameters θ = (W i ), let P θ = {A 1 , . . . , A S } be the underlying tessellation of the input space given by our current choice of parameters; that is, Φ(X; θ) is piece-wise linear and P θ contains those pieces. Now let X be any arbitrary distribution with density p(x) > 0 for all x ∈ R n , for example a Gaussian, and let Y | X d = Φ(X; θ) . Since Φ is invariant under a subgroup of permutations θ σ of its hidden layers, it is easy to see that one can find two parameter values θ A = θ and θ B = θ σ such that F o (θ A ) = F o (θ B ) = 0, but any continuous path γ(t) from θ A to θ B will have a different tessellation and therefore won't satisfy F o (γ(t)) = 0. Moreover, one can build on this counter-example to show that not only the level sets are disconnected, but also that there exist poor local minima. Let θ be a different set of parameters, and Y | X d = Φ(X; θ ) be a different target distribution. Now consider the data distribution given by the mixture

X | p(x) , z ∼ Bernoulli(π) , Y | X, z d = zΦ(X; θ) + (1 -z)Φ(X; θ ) .

By adjusting the mixture component π we can clearly change the risk at θ and θ and make them different, but we conjecture that this preserves the status of local minima of θ and θ . Appendix E constructs a counter-example numerically.

This illustrates an intrinsic difficulty in the optimization landscape if one is after universal guarantees that do not depend upon the data distribution. This difficulty is non-existent in the linear case and not easy to exploit in mean-field approaches such as Choromanska et al. (2015) , and shows that in general we should not expect to obtain connected level sets. However, connectedness can be recovered if one is willing to accept a small increase of energy and make some assumptions on the complexity of the regression task. Our main result shows that the amount by which the energy is allowed to increase is upper bounded by a quantity that trades-off model overparametrization and smoothness in the data distribution.

For that purpose, we start with a characterization of the oracle loss, and for simplicity let us assume Y ∈ R and let us first consider the case with a single hidden layer and 1 regularization: R(θ) = θ 1 .

2.3.2. PRELIMINARIES:

Before proving our main result, we need to introduce preliminary notation and results. We first describe the case with a single hidden layer of size m.

We define

e(m) = min W1∈R m×n , W1(i) 2 ≤1,W2∈R m E{|Φ(X; θ) -Y | 2 } + κ W 2 1 . (6)

to be the oracle risk using m hidden units with norm ≤ 1 and using sparse regression. It is a well known result by Hornik and Cybenko that a single hidden layer is a universal approximator under very mild assumptions, i.e. lim m→∞ e(m) = 0. This result merely states that our statistical setup is consistent, and it should not be surprising to the reader familiar with classic approximation theory. A more interesting question is the rate at which e(m) decays, which depends on the smoothness of the joint density (X, Y ) ∼ P relative to the nonlinear activation family we have chosen.

For convenience, we redefine W = W 1 and β = W 2 and Z(W ) = max(0, W X). We also write z(w) = max(0, w, X ) where (X, Y ) ∼ P and w ∈ R N is any deterministic vector. Let Σ X = E P XX T ∈ R N ×N be the covariance operator of the random input X. We assume Σ X < ∞.

A fundamental property that will be essential to our analysis is that, despite the fact that Z is nonlinear, the quantity

[w 1 , w 2 ] Z := E P {z(w 1 )z(w 2 )} is locally equivalent to the linear metric w 1 , w 2 X = E P {w T 1 XX T w 2 } = w 1 , Σ X w 2

, and that the linearization error decreases with the angle between w 1 and w 2 . Without loss of generality, we assume here that w 1 = w 2 = 1, and we write w 2 Z = E{|z(w)| 2 }. Proposition 2.3. Let α = cos -1 ( w 1 , w 2 ) be the angle between unitary vectors w 1 and w 2 and let w m = w1+w2 w1+w2 be their unitary bisector. Then

EQUATION

The term Σ X is overly pessimistic: we can replace it by the energy of X projected into the subspace spanned by w 1 and w 2 (which is bounded by 2 Σ X ). When α is small, a Taylor expansion of the trigonometric terms reveals that

2 3 Σ X w 1 , w 2 = 2 3 Σ X cos α = 2 3 Σ X (1 - α 2 2 + O(α 4 )) ≤ (1 -α 2 /4) w m 2 Z -Σ X (α 2 /4 + α 2 ) + O(α 4 ) ≤ [w 1 , w 2 ] Z + O(α 4 ) ,

and similarly

[w 1 , w 2 ] Z ≤ w 1 , w 2 w m 2 Z ≤ Σ X w 1 , w 2 .

The local behavior of parameters w 1 , w 2 on our regression problem is thus equivalent to that of having a linear layer, provided w 1 and w 2 are sufficiently close to each other. This result can be seen as a spoiler of what is coming: increasing the hidden layer dimensionality m will increase the chances to encounter pairs of vectors w 1 , w 2 with small angle; and with it some hope of approximating the previous linear behavior thanks to the small linearization error.

In order to control the connectedness, we need a last definition. Given a hidden layer of size m with current parameters W ∈ R n×m , we define a "robust compressibility" factor as

EQUATION

This quantity thus measures how easily one can compress the current hidden layer representation, by keeping only a subset of l its units, but allowing these units to move by a small amount controlled by α. It is a form of n-width similar to Kolmogorov width Donoho (2006) and is also related to robust sparse coding from Tang et al. ( 2013); Ekanadham et al. (2011) .

2.3.3. MAIN RESULT:

Our main result considers now a non-asymptotic scenario given by some fixed size m of the hidden layer. Given two parameter values

θ A = (W A 1 , W A 2 ) ∈ W and θ B = (W B 1 , W B 2 ) with F o (θ {A,B}

) ≤ λ, we show that there exists a continuous path γ : [0, 1] → W connecting θ A and θ B such that its oracle risk is uniformly bounded by max(λ, ), where decreases with model overparametrization.

Theorem 2.4. For any θ A , θ B ∈ W and λ ∈ R satisfying F o (θ {A,B} ) ≤ λ, there exists a continuous path γ :

EQUATION

EQUATION

where C 1 is an absolute constant depending only on κ and P .

Some remarks are in order. First, our regularization term is currently a mix between 2 norm constraints on the first layer and 1 norm constraints on the second layer. We believe this is an artifact of our proof technique, and we conjecture that more general regularizations yield similar results. Next, this result uses the data distribution through the oracle bound e(m) and the covariance term. The extension to empirical risk is accomplished by replacing the probability measure P by the empirical measure P = 1 L l δ ((x, y) -(x l , y l )). However, our asymptotic analysis has to be carefully reexamined to take into account and avoid the trivial regime when M outgrows L. A consequence of Theorem 2.4 is that as m increases, the model becomes asymptotically connected, as proven in the following corollary.

Corollary 2.5. As m increases, the energy gap satisfies = O(m -1 n ) and therefore the level sets become connected at all energy levels. This is consistent with the overparametrization results from Safran & Shamir (2015) ; Shamir (2016) and the general common knowledge amongst deep learning practitioners. Our next sections explore this question, and refine it by considering not only topological properties but also some rough geometrical measure of the level sets.

3. GEOMETRY OF LEVEL SETS:



3.1. THE GREEDY ALGORITHM:

The intuition behind our main result is that, for smooth enough loss functions and for sufficient overparameterization, it should be "easy" to connect two equally powerful models-i.e., two models with F o θ A,B ≤ λ. A sensible measure of this ease-of-connectedness is the normalized length of the geodesic connecting one model to the other: |γ A,B (t)|/|θ A -θ B |. This length represents approximately how far of an excursion one must make in the space of models relative to the euclidean distance between a pair of models. Thus, convex models have a geodesic length of 1, because the geodesic is simply linear interpolation between models, while more non-convex models have geodesic lengths strictly larger than 1.

Because calculating the exact geodesic is difficult, we approximate the geodesic paths via a dynamic programming approach we call Dynamic String Sampling. We comment on alternative algorithms in Appendix A.

For a pair of models with network parameters θ i , θ j , each with F e (θ) below a threshold L 0 , we aim to efficienly generate paths in the space of weights where the empirical loss along the path remains below L 0 . These paths are continuous curves belonging to Ω F (λ)-that is, the level sets of the loss function of interest.

Algorithm 1 Greedy Dynamic String Sampling 1: L0 ← Threshold below which path will be found 2: Φ1 ← randomly initialize θ1, train Φ(xi θ1) to L0 3: Φ2 ← randomly initialize θ2, train Φ(xi θ2) to L0 4: BeadList ←(Φ1, Φ2) 5: Depth ← 0 6: procedure FINDCONNECTION(Φ1, Φ2) 7:

t * ← t such that dγ(θ 1 ,θ 2 ,t) dt t

= 0 OR t = 0.5

8::

Φ3 ← train Φ(xi; t * θ1 + (1 -t * )θ2) to L0

9::

BeadList ← insert(Φ3, after Φ1, BeadList) 10:

M axError1 ← maxt(Fe(tθ3 + (1 -t)θ1))

11:

M axError2 ← maxt(Fe(tθ2 + (1 -t)θ3))

12:

if M axError1 > L0 then return FindConnection(Φ1, Φ3)

13:

if M axError2 > L0 then return FindConnection(Φ3, Φ2)

14:

Depth ← Depth+1

The algorithm recursively builds a string of models in the space of weights which continuously connect θ i to θ j . Models are added and trained until the pairwise linearly interpolated loss, i.e. max t F e (tθ i + (1 -t)θ j ) for t ∈ (0, 1), is below the threshold, L 0 , for every pair of neighboring models on the string. We provide a cartoon of the algorithm in Appendix C.

3.2. FAILURE CONDITIONS AND PRACTICALITIES:

While the algorithm presented will faithfully certify two models are connected if the algorithm converges, it is worth emphasizing that the algorithm does not guarantee that two models are disconnected if the algorithm fails to converge. In general, the problem of determining if two models are connected can be made arbitrarily difficult by choice of a particularly pathological geometry for the loss function, so we are constrained to heuristic arguments for determining when to stop running the algorithm. Thankfully, in practice, loss function geometries for problems of interest are not intractably difficult to explore. We comment more on diagnosing disconnections more carefully in Appendix E.

Further, if the MaxError exceeds L 0 for every new recursive branch as the algorithm progresses, the worst case runtime scales as O(exp(Depth)). Empirically, we find that the number of new models added at each depth does grow, but eventually saturates, and falls for a wide variety of models and architectures, so that the typical runtime is closer to O(poly(Depth))-at least up until a critical value of L 0 .

To aid convergence, either of the choices in line 7 of the algorithm works in practice-choosing t * at a local maximum can provide a modest increase in algorithm runtime, but can be unstable if the the calculated interpolated loss is particularly flat or noisy. t * = .5 is more stable, but slower. Finally, we find that training Φ 3 to αL 0 for α < 1 in line 8 of the algorithm tends to aid convergence without noticeably impacting our numerics. We provide further implementation details in 4.

4. NUMERICAL EXPERIMENTS:

For our numerical experiments, we calculated normalized geodesic lengths for a variety of regression and classification tasks. In practice, this involved training a pair of randomly initialized models to the desired test loss value/accuracy/perplexity, and then attempting to connect that pair of models via the Dynamic String Sampling algorithm. We also tabulated the average number of "beads", or the number intermediate models needed by the algorithm to connect two initial models. For all of the below experiments, the reported losses and accuracies are on a restricted test set. For more complete architecture and implementation details, see our GitHub page.

The results are broadly organized by increasing model complexity and task difficulty, from easiest to hardest. Throughout, and remarkably, we were able to easily connect models for every dataset and architecture investigated except the one explicitly constructed counterexample discussed in Appendix E.1. Qualitatively, all of the models exhibit a transition from a highly convex regime at high loss to a non-convex regime at low loss, as demonstrated by the growth of the normalized length as well as the monotonic increase in the number of required "beads" to form a low-loss connection.

Published as a conference paper at ICLR 2017

4.1. POLYNOMIAL REGRESSION:

We studied a 1-4-4-1 fully connected multilayer perceptron style architecture with sigmoid nonlinearities and RMSProp/ADAM optimization. For ease-of-analysis, we restricted the training and test data to be strictly contained in the interval x ∈ [0, 1] and f (x) ∈ [0, 1]. The number of required beads, and thus the runtime of the algorithm, grew approximately as a power-law, as demonstrated in Table 1 Fig. 1 . We also provide a visualization of a representative connecting path between two models of equivalent power in Appendix D.

0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14

L0 1.2 1.4 1.6 1.8 2.0 Normalized Length (1a) 2 -10 2 -9 2 -8 2 -7 2 -6 2 -5 2 -4 2 -3 L0 2 1 2 2 2 3 2 4 2 5 2 6 2 7

Number of Beads (1b) 0.00 0.05 0.10 0.15 0.20 0.25 0.30

L0 1.0 1.5 2.0 2.5 3.0 Normalized Length (2a) 2 -7 2 -6 2 -5 2 -4 2 -3 2 -2 L0 2 1 2 2 2 3 2 4 2 5 2 6

Number of Beads The cubic regression task exhibits an interesting feature around L 0 = .15 in Table 1 Fig. 2 , where the normalized length spikes, but the number of required beads remains low. Up until this point, the cubic model is strongly convex, so this first spike seems to indicate the onset of non-convex behavior and a concomitant radical change in the geometry of the loss surface for lower loss.

4.2. CONVOLUTIONAL NEURAL NETWORKS:

To test the algorithm on larger architectures, we ran it on the MNIST hand written digit recognition task as well as the CIFAR10 image recognition task, indicated in Table 1 , Figs. 3 and 4 . Again, the data exhibits strong qualitative similarity with the previous models: normalized length remains low until a threshold loss value, after which it grows approximately as a power law. Interestingly, the MNIST dataset exhibits very low normalized length, even for models nearly at the state of the art in classification power, in agreement with the folk-understanding that MNIST is highly convex and/or "easy". The CIFAR10 dataset, however, exhibits large non-convexity, even at the modest test accuracy of 80%.

4.3. RECURRENT NEURAL NETWORKS:

To gauge the generalizability of our algorithm, we also applied it to an LSTM architecture for solving the next word prediction task on the PTB dataset, depicted in Table 1 Fig. 5 . Noteably, even for a radically different architecture, loss function, and data set, the normalized lengths produced by the DSS algorithm recapitulate the same qualitative features seen in the above datasets-i.e., models can be easily connected at high perplexity, and the normalized length grows at lower and lower perplexity after a threshold value, indicating an onset of increased non-convexity of the loss surface.

5. DISCUSSION:

We have addressed the problem of characterizing the loss surface of neural networks from the perspective of gradient descent algorithms. We explored two angles -topological and geometrical aspects -that build on top of each other.

On the one hand, we have presented new theoretical results that quantify the amount of uphill climbing that is required in order to progress to lower energy configurations in single hidden-layer ReLU networks, and proved that this amount converges to zero with overparametrization under mild conditions. On the other hand, we have introduced a dynamic programming algorithm that efficiently approximates geodesics within each level set, providing a tool that not only verifies the connectedness of level sets, but also estimates the geometric regularity of these sets. Thanks to this information, we can quantify how 'non-convex' an optimization problem is, and verify that the optimization of quintessential deep learning tasks -CIFAR-10 and MNIST classification using CNNs, and next word prediction using LSTMs -behaves in a nearly convex fashion up until they reach high accuracy levels.

That said, there are some limitations to our framework. In particular, we do not address saddle-point issues that can greatly affect the actual convergence of gradient descent methods. There are also a number of open questions; amongst those, in the near future we shall concentrate on:

• Extending Theorem 2.4 to the multilayer case. We believe this is within reach, since the main analytic tool we use is that small changes in the parameters result in small changes in the covariance structure of the features. That remains the case in the multilayer case. • Empirical versus Oracle Risk. A big limitation of our theory is that right now it does not inform us on the differences between optimizing the empirical risk versus the oracle risk.

Understanding the impact of generalization error and stochastic gradient in the ability to do small uphill climbs is an open line of research. • Influence of symmetry groups. Under appropriate conditions, the presence of discrete symmetry groups does not prevent the loss from being connected, but at the expense of increasing the capacity. An important open question is whether one can improve the asymptotic properties by relaxing connectedness to being connected up to discrete symmetry. • Improving numerics with Hyperplane method. Our current numerical experiments employ a greedy (albeit faster) algorithm to discover connected components and estimate geodesics. We plan to perform experiments using the less greedy algorithm described in Appendix A. 

A CONSTRAINED DYNAMIC STRING SAMPLING:

While the algorithm presented in Sec. 3.1 is fast for sufficiently smooth families of loss surfaces with few saddle points, here we present a slightly modified version which, while slower, provides more control over the convergence of the string. We did not use the algorithm presented in this section for our numerical studies.

Instead of training intermediate models via full SGD to a desired accuracy as in step 8 of the algorithm, intermediate models are be subject to a constraint that ensures they are "close" to the neighboring models on the string. Specifically, intermediate models are constrained to the unique hyperplane in weightspace equidistant from its two neighbors. This can be further modified by additional regularization terms to control the "springy-ness" of the string. These heuristics could be chosen to try to more faithfully sample the geodesic between two models.

In practice, for a given model on the string, θ i , these two regularizations augment the standard loss by:

F (θ) = F (θ) + ζ( θ i-1 -θ i + θ i+1 -θ i ) + κ (θi-1-θi+1)/2 (θi-1-θi+1)/2 • (θi-(θi-1-θi+1)/2) (θi-(θi-1-θi+1)/2)

. The ζ regularization term controls the "springy-ness" of the weightstring, and the κ regularization term controls how far off the hyperplane a new model can deviate.

Because adapting DSS to use this constraint is straightforward, here we will describe an alternative "breadth-first" approach wherein models are trained in parallel until convergence. This alternative approach has the advantage that it will indicate a disconnection between two models "sooner" in training. The precise geometry of the loss surface will dictate which approach to use in practice.

Given two random models σ i and σ j where |σ i -σ j | < L 0 , we aim to follow the evolution of the family of models connecting σ i to σ j . Intuitively, almost every continuous path in the space of random models connecting σ i to σ j has, on average, the same (high) loss. For simplicity, we choose to initialize the string to the linear segment interpolating between these two models. If this entire segment is evolved via gradient descent, the segment will either evolve into a string which is entirely contained in a basin of the loss surface, or some number of points will become fixed at a higher loss. These fixed points are difficult to detect directly, but will be indirectly detected by the persistence of a large interpolated loss between two adjacent models on the string.

The algorithm proceeds as follows:

(0.) Initialize model string to have two models, σ i and σ j .

1. Begin training all models to the desired loss, keeping the instantaneous loss, L 0 (t), of all models being trained approximately constant.

2. If the pairwise interpolated loss between σ n and σ n+1 exceeds L 0 (t), insert a new model at the maximum of the interpolated loss (or halfway) between these two models. 2) until all models (and interpolated errors) are below a threshold loss L 0 (t final ) := L 0 , or until a chosen failure condition (see 3.2).

3.. Repeat steps (1) and (:



B PROOFS B.1 PROOF OF PROPOSITION 2.1:

Suppose that θ 1 is a local minima and θ 2 is a global minima, but F (θ 1 ) > F (θ 2 ). If λ = F (θ 1 ), then clearly θ 1 and θ 2 both belong to Ω F (λ). Suppose now that Ω F (λ) is connected. Then we Published as a conference paper at ICLR 2017 could find a smooth (i.e. continuous and differentiable) path γ(t) with γ(0) = θ 1 , γ(1) = θ 2 and F (γ(t)) ≤ λ = F (θ 1 ). But this contradicts the strict local minima status of θ 1 , and therefore Ω F (λ) cannot be connected .

B.2 PROOF OF PROPOSITION 2.2:

Let us first consider the case with κ = 0. We proceed by induction over the number of layers K. For K = 1, the loss F (θ) is convex. Let θ A , θ B be two arbitrary points in a level set Ω λ . Thus F (θ A ) ≤ λ and F (θ B ) ≤ λ. By definition of convexity, a linear path is sufficient in that case to connect θ A and θ B :

F ((1 -t)θ A + tθ B ) ≤ (1 -t)F (θ A ) + tF (θ B ) ≤ λ . Suppose the result is true for K -1. Let θ A = (W A 1 , . . . , W A K ) and θ B = (W B 1 , . . . , W B K ) with F (θ A ) ≤ λ, F (θ B ) ≤ λ. Since n j ≥ min(n 1 , n K ) for j = 2 . . . K -1, we can find k * = {1, K -1} such that n k * ≥ min(n k * -1 , n k * +1 ). For each W 1 , . . . , W K , we denote Wj = W j for j = k * , k * -1 and Wk * = W k * -1 W k * .

By induction hypothesis, the loss expressed in terms of θ = ( W1 , . . . , WK-1 ) is connected between θA and θB . Let Wk * (t) the corresponding linear path projected in the layer k * . We need to produce a path in the variables

W k * -1 (t), W k * (t) such that: i W k * -1 (0) = W A k * -1 , W k * -1 (1) = W B k * -1 , ii W k * (0) = W A k * , W k * (1) = W B k * , iii W k * (t)W k * -1 (t) = Wk * -1 (t) for t ∈ (0, 1).

We construct it as follows. Let

W k * (t) = tW B k * + (1 -t)W A k * + t(1 -t)V , W k * -1 (t) = W k * (t) † Wk * -1 (t) ,

where

W k * (t) † = (W k * (t) T W k * (t)) -1 W k * (t)

T denotes the pseudoinverse and V is a n k * -1 × n k * matrix drawn from a iid distribution. Conditions (i) and (ii) are immediate from the definition, and condition (iii) results from the fact that

W k * (t)W k * (t) † = I N * k , since W * k (t)

has full rank for all t ∈ (0, 1). Finally, let us prove that the result is also true when K = 2 and κ > 0. We construct the path using the variational properties of atomic norms Bach (2013) . When we pick the ridge regression regularization, the corresponding atomic norm is the nuclear norm:

X * = min U V T =X 1 2 ( U 2 + V 2 ) .

The path is constructed by exploiting the convexity of the variational norm

EQUATION

2 ) .

From ( 12) it results that the loss F o (W 1 , W 2 ) can be minored by another loss expressed in terms of W of the form E{|Y -W X| 2 } + 2κ W * , which is convex with respect to W . Thus a linear path in W from W A to W B is guaranteed to be below F o (θ {A,B} ). Let us define

∀ t , W 1 (t), W 2 (t) = arg min U V T = W (t) ( U 2 + V 2 ) .

One can verify that we can first consider a path

(β A 1 (s), β A 2 (s)) from (W A 1 , W A 2 ) to (W 1 (0), W 2 (0) such that ∀ s β 1 (s)β 2 (s) = W A and β 1 (s) 2 + β 2 (s) 2 decreases ,

and similarly for (W B 1 , W B 2 ) to (W 1 (1), W 2 (1). The path (β A {1,2} (s), W {1,2} (t), β B {1,2} (s)) satisfies (i-iii) by definition. We also verify that

W 1 (t) 2 + W 2 (t) 2 = 2 W (t) * ≤ 2(1 -t) W (0) * + 2t W (1) * ≤ (1 -t)( W 2 1 (0) + W 2 2 (0)) + t( W 2 1 (1) + W 2 2 (1)

) . Finally, we verify that the paths we have just created, when applied to θ A arbitrary and θ B = θ * a global minimum, are strictly decreasing, again by induction. For K = 1, this is again an immediate consequence of convexity. For K > 1, our inductive construction guarantees that for any 0 < t < 1, the path θ(t) = (W k (t)) k≤K satisfies F o (θ(t)) < F o (θ A ). This concludes the proof .

B.3 PROOF OF PROPOSITION 2.3 Let:

A(w 1 , w 2 ) = {x ∈ R n ; x, w 1 ≥ 0 , x, w 2 ≥ 0} . By definition, we have w 1 , w 2 Z = E{max(0, X, w 1 ) max(0, X, w 2 )} (13)

=:

A(w1,w2)

x, w 1 x, w 2 dP (x) ,

EQUATION

where Q is the orthogonal projection onto the space spanned by w 1 and w 2 and d P (x) = d P (x 1 , x 2 ) is the marginal density on that subspace. Since this projection does not interfere with the rest of the proof, we abuse notation by dropping the Q and still referring to dP (x) as the probability density. 

EQUATION

We conclude by bounding each error term E 1 and E 2 separately:

EQUATION

since every point in B by definition has angle greater than π/2 -α from w m . Also,

EQUATION

by direct application of Cauchy-Schwartz. The proof is completed by plugging the bounds from ( 18) and ( 19) into (17) .

We have that

EQUATION

Finally, we verify that

E{|Y -η2(t)Z(η1(t))| 2 } -E{|Y -(1 -t)β1Z(W A ) -tγAZ( WA)| 2 } ≤ (22) ≤ 4α max(E|Y | 2 , E|Y 2 |) ΣX (κ -1/2 + α E|Y 2 |κ -1 ) + o(α 2 ) .

Indeed, from Proposition 2.3, and using the fact that

∀ i ≤ M, t ∈ [0, 1] , ∠((1 -t)w A i + t wA i ; w A i ) ≤ α , ∠((1 -t)w A i + t wA i ; wA i ) ≤ α we can write (1 -t)β 1,i z(w A i ) + tγ A,i z( wA i ) d = η 2 (t) i z(η 1 (t) i ) + n i , with E{|n i | 2 } ≤ 4|η 2 (t) i | 2 Σ X α 2 + O(α 4 ) and E|n i | ≤ 2|η 2 (t) i |α Σ X using concavity of the moments. Thus E{|Y -η 2 (t)Z(η 1 (t))| 2 } -E{|Y -(1 -t)β 1 Z(W A ) -tγ A Z( WA )| 2 } ≤ 2E i (Y -η 2 (t)Z(η 1 (t)))n i + E | i n i | 2 ≤ 4 α E|Y 2 | Σ X η 2 + α 2 ( η 2 1 ) 2 Σ X ≤ 4α max(1, E|Y 2 |) Σ X ( η 2 1 + α η 2 2 1 ) + o(α 2 ) ≤ 4α max( E|Y 2 |, E|Y 2 |) Σ X (κ -1 + α E|Y 2 |κ -2 ) + o(α 2 ) , which proves (22).

We have just constructed a path from θ A to θ B , in which all subpaths except (2) and ( 5) have energy maximized at the extrema due to convexity, given respectively by λ, δ W 1 A (m, 0, m), δ W 1 A (ml, α, m), e(l), δ W 1 B (m -l, α, m), and δ W 1 B (m, 0, m). For the two subpaths (2) and ( 5), ( 22) shows that it is sufficient to add the corresponding upper bound to the linear subpath, which is of the form Cα + o(α 2 ) where C is an explicit constant independent of θ. Since l and α are arbitrary, we are free to pick the infimum, which concludes the proof.

B.5 PROOF OF COROLLARY 2.5:

Let us consider a generic first layer weight matrix W ∈ R n×m . Without loss of generality, we can assume that w k = 1 for all k, since increasing the norm of w k within the unit ball has no penalty in the loss, and we can compensate this scaling in the second layer thanks to the homogeneity of the half-rectification. Since this results in an attenuation of these second layer weights, they too are guaranteed not to increase the loss.

From Vershynin (2010) [Lemma 5.2] we verify that the covering number N (S n-1 , ) of the Euclidean unit sphere S n-1 satisfies

N (S n-1 , ) ≤ 1 + 2 n ,

which means that we can cover the unit sphere with an -net of size N (S n-1 , ).

Let 0 < η < n -1 (1 + n -1 ) -1 , and let us pick, for each m, m = m η-1 n . Let us consider its corresponding -net of size

u m = N (S n-1 , m ) 1 + 2 m n m 1-η .

Since we have m vectors in the unit sphere, it results from the pigeonhole principle that at least one element of the net will be associated with at least v m = mu -1 m m η vectors; in other words, we are guaranteed to find amongst our weight vector W a collection Q m of v m m η vectors that are all at an angle at most 2 m apart. Let us now apply Theorem 2.4 by picking n = v m and α = m . We need to see that the terms involved in the bound all converge to 0 as m → ∞.

The contribution of the oracle error e(v m ) -e(m) goes to zero as m → ∞ by the fact that lim m→∞ e(m) exists (it is a decreasing, positive sequence) and that v m → ∞.

Let us now verify that δ(m -v m , m , m) also converges to zero. We are going to prune the first layer by removing one by one the vectors in Q m . Removing one of these vectors at a time incurs in an error of the order of m . Indeed, let w k be one of such vectors and let β be the solution of

min β E(β ) = min β =(β f ;β k )∈R k E{|Y -β T f Z(W -k ) -β k z(w k )| 2 } + κ( β f 1 + |β k |) ,

where W -k is a shorthand for the matrix containing the rest of the vectors that have not been discarded yet. Removing the vector w k from the first layer increases the loss by a factor that is upper bounded by E(β p ) -E(β), where

(β p ) j = β j for j < k -1 , β k-1 + β k otherwise.

, since now β p is a feasible solution for the pruned first layer.

Let us finally bound E(β p ) -E(β).

Since ∠(w k , w k-1 ) ≤ m , it results from Proposition 2.3 that

z(w k ) d = z(w k-1 ) + n , with E{|n| 2 } ≤ Cα 2 for some constant C independent of m. By redefining p 1 = Y -β T p Z(W -k )- 1 2 n and p 2 = 1 2 n, we have E{|Y -β T p Z(W -k )| 2 } -E{|Y -β T Z(W -k ) -β k z(w k )| 2 } = E{|p 1 + p 2 | 2 } -E{|p 1 -p 2 | 2 } = 4E{|p 1 p 2 |} ≤ E Y -β T p Z(W -k ) - 1 2 n 2 E{|n| 2 } ≤ (C + α)α m ,

where C only depends on E{|Y | 2 }. We also verify that β p 1 ≤ β 1 .

It results that removing |Q m | of such vectors incurs an increase of the loss at most

|Q m | m m η m η-1 n = m η+ η-1 n .

Since we picked η such that η + η-1 n < 0, this term converges to zero. The proof is finished.

C CARTOON OF ALGORITHM:

Refer to Fig. 2 .

D VISUALIZATION OF CONNECTION:

Because the weight matrices are anywhere from high to extremely high dimensional, for the purposes of visualization we projected the models on the connecting path into a three dimensionsal subspace. Snapshots of the algorithm in progress for the quadratic regression task are indicated in Fig. 3 . This was done by vectorizing all of the weight matrices for all the beads for a given connecting path, and then performing principal component analysis to find the three highest weight projections for the collection of models that define the endpoints of segments for a connecting path-i.e., the The interpolated model Θ(θ i , θ j , t * ) is added and labeled θ i,j . d) : Stochastic gradient descent is performed on the interpolated model until its loss is below αL 0 . e) : New interpolated loss curves are calculated between the models, pairwise on a chain. f ) : As in step c), a new model is inserted at the maxima of the interpolated loss curve between θ i and θ i,j . g) : As in step d), gradient descent is performed until the model has low enough loss. θ i discussed in the algorithm. We then projected the connecting string of models onto these three directions.

The color of the strings was chosen to be representative of the test loss under a log mapping, so that extremely high test loss mapped to red, whereas test loss near the threshold mapped to blue. An animation of the connecting path can be seen on our Github page.

Finally, projections onto pairs of principal components are indicated by the black curves.

E A DISCONNECTION E.1 A DISCONNECTION:

As a sanity check for the algorithm, we also applied it to a problem for which we know that it is not possible to connect models of equivalent power by the arguments of section 2.3.1. The input data is 3 points in R 2 , and the task is to permute the datapoints, i.e. map {x 1 , x 2 , x 3 } → {x 2 , x 3 , x 1 }. This map requires at least 12 parameters in general for the three linear maps which take x i → x j for i, j ∈ {{1, 2}, {2, 3}, {3, 1}}. Our archticture was a 2-3-2 fully connected neural network with a single relu nonlinearity after the hidden layer-a model which clearly has 12 free parameters by construction. The two models we tried to connect were a single model, θ, and a copy of θ with the first two neurons in the hidden layer permuted, θσ . The algorithm fails to converge when initialized with these two models. We provide a visualization of the string of models produced by the algorithm in Fig. 4 .

In general, a persistent high interpolated loss between two neighboring beads on the string of models could arise from either a slowly converging, connected pair of models or from a truly disconnected pair of models. "Proving" a disconnection at the level of numerical experiments is intractable in general, but a collection of negative results-i.e., failures to converge-are highly suggestive of a true disconnection. 

:



Figures and Tables:

Figure 1: (Column a) Average normalized geodesic length and (Column b) average number of beads versus loss. (1) A quadratic regression task. (2) A cubic regression task. (3) A convnet for MNIST. (4) A convnet inspired by Krizhevsky for CIFAR10. (5) A RNN inspired by Zaremba for PTB next word prediction.

By construction we have w 1 = rw m -d , w 2 = rw m + d , and thus x, w 1 x, w 2 = r 2 | x, w m | 2 -| x, d | 2 . (16) By denoting C(w m ) = {x ∈ R n ; x, w m ≥ 0}, observe that A(w 1 , w 2 ) ⊆ C(w m ). Let us denote by B = C(w m ) \ A(w 1 , w 2 ) the disjoint complement. It results that

Figure 2: A cartoon of the algorithm. a) : The initial two models with approximately the same loss, L 0 . b) : The interpolated loss curve, in red, and its global maximum, occuring at t = t * . c) :The interpolated model Θ(θ i , θ j , t * ) is added and labeled θ i,j . d) : Stochastic gradient descent is performed on the interpolated model until its loss is below αL 0 . e) : New interpolated loss curves are calculated between the models, pairwise on a chain. f ) : As in step c), a new model is inserted at the maxima of the interpolated loss curve between θ i and θ i,j . g) : As in step d), gradient descent is performed until the model has low enough loss.

Figure 3: Snapshots of Dynamic String Sampling in action for the quadratic regression task. The string's coordinates are its projections onto the three most important principal axes of the fully converged string. (Top Left) One step into the algorithm, note the high loss between all of the vertices of the path. (Top Right) An intermediate step of the algorithm. Portions of the string have converged, but there are still regions with high interpolated loss. (Bottom Left) Near the end of the algorithm. Almost the entire string has converged to low loss. (Bottom Right) The algorithm has finished. A continuous path between the models has been found with low loss.

Figure 4: These three figures are projections of the components of the 12-dimensional weight matrices which comprise the models on the string produced by the DSS algorithm. The axes are the principal components of the weight matrices, and the colors indicate test error for the model. For more details on the figure generation, see Appendix D. (Left) The string of models after 1 step. Note the high error at all points except the middle and the endpoints. (Middle) An intermediate stage of the algorithm. Part of the string has converged, but a persistent high-error segment still exists. (Right) Even after running for many steps, the error persists, and the algorithm does not converge.

Gongguo Tang, Badri Narayan Bhaskar, Parikshit Shah, and Benjamin Recht. Compressed sensing off the grid. IEEE Transactions on Information Theory, 59(11):7465-7490, 2013. Yuandong Tian. Symmetry-breaking convergence analysis of certain two-layered neural networks with relu nonlinearity. ICLR Workshop 2017, 2017.

Table content: <table><tr><td>Roman Vershynin. Introduction to the non-asymptotic analysis of random matrices. arXiv preprint</td></tr><tr><td>arXiv:1011.3027, 2010.</td></tr></table>

Peer review (iclr2017):
This paper studies the energy landscape of the loss function in neural networks. It is generally clearly written and nicely provides intuitions for the results. One main contribution is to show that the level sets of the loss becomes connected as the network is increasingly overparameterized. It also quantifies, in a way, the degree of disconnectedness possible in terms of the increase in loss that one must allow to find a connected path. It would seem that this might have some implications for the likelihood of escaping local minima with stochastic gradient descent. The paper also presents a simple algorithm for finding geodesic paths between two networks such that the loss is decreasing along the path. Using this they show that the loss seems to become more nonconvex when the loss is smaller. This is also quite interesting. The work does have some significant limitations, which is not surprising given the difficulty of fully analyzing the network loss function. However, the authors are quite clear about these limitations, which especially include not yet analyzing deep networks and analyzing only the oracle loss, and not the empirical loss. I would have also appreciated a little more practical discussion of the bound in Theorem 2.4. It is hard to tell whether this bound is tight enough to be practically relevant.