Topology and Geometry of Half-Rectified Network Optimization

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(\theta)$ of a given model can be expressed in terms of its level sets $\Omega_{\lambda}$, which contain for each energy level $\lambda$ all parameters $\theta$ yielding a loss smaller or equal than $\lambda$. 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(\theta)$ 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.

Topology of Level Sets

Let $P$ be a probability measure on a product space $\mathcal{X}\times\mathcal{Y}$, where we assume $\mathcal{X}$ and $\mathcal{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

We define the level set of $F(\theta)$ as

If $\Omega_{F}(\lambda)$ is connected for all $\lambda$ then every local minima of $F(\theta)$ is a global minima.

Strict local minima implies that $\nabla F(\theta)=0$ and $HF(\theta)\succeq 0$, but avoids degenerate cases where $F$ is constant along a manifold intersecting $\theta$. In that scenario, if $\mathcal{U}_{\theta}$ denotes that manifold, our reasoning immediately implies that if $\Omega_{F}(\lambda)$ are connected, then for all $\epsilon>0$ there exists $\theta^{\prime}$ with $\text{dist}(\theta^{\prime},\mathcal{U}_{\theta})\leq\epsilon$ and $F(\theta^{\prime})<F(\theta)$. In other words, some element at the boundary of $\mathcal{U}_{\theta}$ must be a saddle point. A stronger property that eliminates the risk of gradient descent getting stuck at $\mathcal{U}_{\theta}$ is that all elements at the boundary of $\mathcal{U}_{\theta}$ are saddle points. This can be guaranteed if one can show that there exists a path connecting any $\theta$ to the lowest energy level such that $F$ is strictly decreasing along it.

In particular, $\mathcal{U}_{\theta}$ 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}=\text{diag}(\alpha_{1},\dots,\alpha_{N_{k}})$ with $\alpha_{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 The Linear Case

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

and the ridge regression $\mathcal{R}(\theta)=\|\theta\|^{2}$. This model defines a non-convex (and non-concave) loss $F_{e}(\theta)$. When $\kappa=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 $\kappa>0$ in the case $K=2$.

Let $W_{1},W_{2},\dots,W_{K}$ be weight matrices of sizes $n_{k}\times n_{k+1}$, $k<K$, and let $F_{e}(\theta)$, $F_{o}(\theta)$ denote the risk minimizations using $\Phi$ as in (4). Assume that $n_{j}>\min(n_{1},n_{K})$ for $j=2\dots K-1$. Then $\Omega_{F_{e}}(\lambda)$ (and $\Omega_{F_{o}}$) is connected for all $\lambda$ and all $K$ when $\kappa=0$, and for $\kappa>0$ when $K=2$; and therefore there are no poor local minima in these cases. Moreover, any $\theta$ 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\dots K-1$ and our analysis does not inform about the order of the saddle points, we do not need full rank assumptions on $\Sigma_{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.

3 Half-Rectified Nonlinear Case

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

where $b_{i}$ contains the biases for each layer. For simplicity, we continue to use $W_{i}$ and $x$ in the following.

By adjusting the mixture component $\pi$ we can clearly change the risk at $\theta$ and $\theta^{\prime}$ and make them different, but we conjecture that this preserves the status of local minima of $\theta$ and $\theta^{\prime}$. 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.

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$.

to be the oracle risk using $m$ hidden units with norm $\leq 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\to\infty}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)\sim P$ relative to the nonlinear activation family we have chosen.

Let $\alpha=\cos^{-1}(\langle w_{1},w_{2}\rangle)$ be the angle between unitary vectors $w_{1}$ and $w_{2}$ and let $w_{m}=\frac{w_{1}+w_{2}}{\|w_{1}+w_{2}\|}$ be their unitary bisector. Then

The term $\|\Sigma_{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\|\Sigma_{X}\|$). When $\alpha$ is small, a Taylor expansion of the trigonometric terms reveals that

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.

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 $\alpha$. 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).

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 $\theta^{A}=(W_{1}^{A},W_{2}^{A})\in\mathcal{W}$ and $\theta^{B}=(W_{1}^{B},W_{2}^{B})$ with $F_{o}(\theta^{\{A,B\}})\leq\lambda$, we show that there exists a continuous path $\gamma:\to\mathcal{W}$ connecting $\theta^{A}$ and $\theta^{B}$ such that its oracle risk is uniformly bounded by $\max(\lambda,\epsilon)$, where $\epsilon$ decreases with model overparametrization.

where $C_{1}$ is an absolute constant depending only on $\kappa$ and $P$.

As $m$ increases, the energy gap $\epsilon$ satisfies $\epsilon=O(m^{-\frac{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.

Geometry of Level Sets

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}{\theta^{A,B}}\leq\lambda$. A sensible measure of this ease-of-connectedness is the normalized length of the geodesic connecting one model to the other: $|\gamma_{A,B}(t)|/|\theta_{A}-\theta_{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 $\theta_{i}$, $\theta_{j}$, each with $F_{e}(\theta)$ 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 $\Omega_{F}(\lambda)$–that is, the level sets of the loss function of interest.

The algorithm recursively builds a string of models in the space of weights which continuously connect $\theta_{i}$ to $\theta_{j}$. Models are added and trained until the pairwise linearly interpolated loss, i.e. $\rm{max}_{t}F_{e}(t\theta_{i}\ +\ (1-t)\theta_{j})$ for $t\in(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.

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 $\rm{\mathbf{MaxError}}$ exceeds $L_{0}$ for every new recursive branch as the algorithm progresses, the worst case runtime scales as $O(\rm{exp}(\rm{\mathbf{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(\rm{poly}(\rm{\mathbf{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 $\Phi_{3}$ to $\alpha L_{0}$ for $\alpha<1$ in line $8$ of the algorithm tends to aid convergence without noticeably impacting our numerics. We provide further implementation details in 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.

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\in$ and $f(x)\in$. 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.

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.

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%.

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.

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.

We would like to thank Mark Tygert for pointing out the reference to the $\epsilon$-nets and Kolmogorov capacity, and Martin Arjovsky for spotting several bugs in early version of the results. We would also like to thank Maithra Raghu and Jascha Sohl-Dickstein for enlightening discussions, as well as Yasaman Bahri for helpful feedback on an early version of the manuscript. CDF was supported by the NSF Graduate Research Fellowship under Grant DGE-1106400.

References

Appendix 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.

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 $\sigma_{i}$ and $\sigma_{j}$ where $|\sigma_{i}-\sigma_{j}|<L_{0}$, we aim to follow the evolution of the family of models connecting $\sigma_{i}$ to $\sigma_{j}$. Intuitively, almost every continuous path in the space of random models connecting $\sigma_{i}$ to $\sigma_{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.

(0.) Initialize model string to have two models, $\sigma_{i}$ and $\sigma_{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 $\sigma_{n}$ and $\sigma_{n+1}$ exceeds $L_{0}(t)$, insert a new model at the maximum of the interpolated loss (or halfway) between these two models.

3. Repeat steps (1) and (2) until all models (and interpolated errors) are below a threshold loss $L_{0}(t_{\rm{final}}):=L_{0}$, or until a chosen failure condition (see 3.2).

Appendix B Proofs

Suppose that $\theta_{1}$ is a local minima and $\theta_{2}$ is a global minima, but $F(\theta_{1})>F(\theta_{2})$. If $\lambda=F(\theta_{1})$, then clearly $\theta_{1}$ and $\theta_{2}$ both belong to $\Omega_{F}(\lambda)$. Suppose now that $\Omega_{F}(\lambda)$ is connected. Then we could find a smooth (i.e. continuous and differentiable) path $\gamma(t)$ with $\gamma(0)=\theta_{1}$, $\gamma(1)=\theta_{2}$ and $F(\gamma(t))\leq\lambda=F(\theta_{1})$. But this contradicts the strict local minima status of $\theta_{1}$, and therefore $\Omega_{F}(\lambda)$ cannot be connected $\square$.

B.2 Proof of Proposition 2.2

Let us first consider the case with $\kappa=0$. We proceed by induction over the number of layers $K$. For $K=1$, the loss $F(\theta)$ is convex. Let $\theta^{A}$, $\theta^{B}$ be two arbitrary points in a level set $\Omega_{\lambda}$. Thus $F(\theta^{A})\leq\lambda$ and $F(\theta^{B})\leq\lambda$. By definition of convexity, a linear path is sufficient in that case to connect $\theta^{A}$ and $\theta^{B}$:

$W_{k^{*}-1}(0)=W_{k^{*}-1}^{A}$, $W_{k^{*}-1}(1)=W_{k^{*}-1}^{B}$,

$W_{k^{*}}(0)=W_{k^{*}}^{A}$, $W_{k^{*}}(1)=W_{k^{*}}^{B}$,

$W_{k^{*}-1}(t)$ has the property that $W_{k^{*}-1}(0)=W_{k^{*}-1}^{A}$, $W_{k^{*}-1}(1)=W_{k^{*}-1}^{B}$. Thanks to the fact that $\text{rank}(W_{k^{*}-1}(t))=m$ for all $t\in(0,1)$, there exists $W_{k^{*}}(t)$ such that

Finally, we need to show that the path $W_{k^{*}}(t)$ is continuous and satisfies $W_{k^{*}}(0)=W_{k^{*}}^{A}$, $W_{k^{*}}(1)=W_{k^{*}}^{B}$. Since by construction the paths are continuous in $t\in(0,1)$, it only remains to be shown that

Finally, if either $\text{rank}(W_{k^{*}-1}^{A})<m$ or $\text{rank}(W_{k^{*}-1}^{B})<m$, we denote by $P_{A}$ (resp $P_{B}$) the orthogonal complement of $\text{span}(W_{k^{*}-1}^{A})$ (resp $\text{span}(W_{k^{*}-1}^{B})$), and by $Q_{A}$ (resp $Q_{B}$) the orthogonal complement of $\text{Null}(W_{k^{*}}^{A})$ (resp $\text{Null}(W_{k^{*}}^{B})$). Observe that if either $Q_{A}$ intersects with $P_{A}$ (resp $Q_{B}$ intersects with $P_{B}$), we can shrink $W_{k^{*}}^{A}$ in the intersection with no effect in the loss. We can thus assume without loss of generality that $P_{A}\cap Q_{A}=\emptyset$. In that case, increasing the range of $W_{k*-1}^{A}$ until it has rank $m$ has no effect in the loss either, since the new directions will fall in the kernel of $W_{k^{*}}^{A}$. Therefore, by applying the necessary corrections to $W_{k*-1}^{A}$ and $W_{k*}^{A}$ (resp $W_{k*-1}^{B}$ and $W_{k*}^{B}$) we can reduce ourselves to the previous case.

Finally, let us prove that the result is also true when $K=2$ and $\kappa>0$. We construct the path using the variational properties of atomic norms . When we pick the ridge regression regularization, the corresponding atomic norm is the nuclear norm:

One can verify that we can first consider a path $(\beta^{A}_{1}(s),\beta^{A}_{2}(s))$ from $(W_{1}^{A},W_{2}^{A})$ to $(W_{1}(0),W_{2}(0)$ such that

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

Finally, we verify that the paths we have just created, when applied to $\theta^{A}$ arbitrary and $\theta^{B}=\theta^{*}$ 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 $\theta(t)=(W_{k}(t))_{k\leq K}$ satisfies $F_{o}(\theta(t))<F_{o}(\theta^{A})$. This concludes the proof $\square$.

B.3 Proof of Proposition 2.3

where $Q$ is the orthogonal projection onto the space spanned by $w_{1}$ and $w_{2}$ and $d\bar{P}(x)=d\bar{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.

Now, let $r=\frac{1}{2}\|w_{1}+w_{2}\|=\frac{1+\cos(\alpha)}{2}$ and $d=\frac{w_{2}-w_{1}}{2}$. By construction we have

We conclude by bounding each error term $E_{1}$ and $E_{2}$ separately:

since every point in $B$ by definition has angle greater than $\pi/2-\alpha$ from $w_{m}$. Also,

by direct application of Cauchy-Schwartz. The proof is completed by plugging the bounds from (22) and (23) into (B.3) $\square$.

B.4 Proof of Theorem 2.4

Consider a generic $\alpha$ and $l\leq m$. A path from $\theta^{A}$ to $\theta^{B}$ will be constructed by concatenating the following paths:

from $\theta^{A}$ to $\theta_{lA}$, the best linear predictor using the same first layer as $\theta^{A}$,

from $\theta_{lA}$ to $\theta_{sA}$, the best $(m-l)$-term approximation using perturbed atoms from $\theta^{A}$,

from $\theta_{sA}$ to $\theta^{*}$ the oracle $l$ term approximation,

from $\theta^{*}$ to $\theta_{sB}$, the best $(m-l)$-term approximation using perturbed atoms from $\theta^{B}$,

from $\theta_{sB}$ to $\theta_{lB}$, the best linear predictor using the same first layer as $\theta^{B}$,

The proof will study the increase in the loss along each subpath and aggregate the resulting increase into a common bound.

Subpaths (1) and (6) only involve changing the parameters of the second layer while leaving the first-layer weights fixed, which define a convex loss. Therefore a linear path is sufficient to guarantee that the loss along that path will be upper bounded by $\lambda$ on the first end and $\delta_{W_{1}^{A}}(m,0,m)$ on the other end.

and similarly for $\overline{\beta}_{2}$. By convexity, the augmented linear path $\eta(t)=(1-t)\overline{\beta}_{1}+t\overline{\beta}_{2}$ thus satisfies

Let us now approximate this augmented linear path with a path in terms of first and second layer weights. We consider

Indeed, from Proposition 2.3, and using the fact that

We have just constructed a path from $\theta^{A}$ to $\theta^{B}$, in which all subpaths except (2) and (5) have energy maximized at the extrema due to convexity, given respectively by $\lambda$, $\delta_{W_{A}^{1}}(m,0,m)$, $\delta_{W_{A}^{1}}(m-l,\alpha,m)$, $e(l)$, $\delta_{W_{B}^{1}}(m-l,\alpha,m)$, and $\delta_{W_{B}^{1}}(m,0,m)$. For the two subpaths (2) and (5), (26) shows that it is sufficient to add the corresponding upper bound to the linear subpath, which is of the form $C\alpha+o(\alpha^{2})$ where $C$ is an explicit constant independent of $\theta$. Since $l$ and $\alpha$ are arbitrary, we are free to pick the infimum, which concludes the proof. $\square$

B.5 Proof of Corollary 2.5

From [Lemma 5.2] we verify that the covering number $\mathcal{N}(S^{n-1},\epsilon)$ of the Euclidean unit sphere $S^{n-1}$ satisfies

which means that we can cover the unit sphere with an $\epsilon$-net of size $\mathcal{N}(S^{n-1},\epsilon)$.

Let $0<\eta<n^{-1}(1+n^{-1})^{-1}$, and let us pick, for each $m$, $\epsilon_{m}=m^{\frac{\eta-1}{n}}$. Let us consider its corresponding $\epsilon$-net of size

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_{m}^{-1}\simeq m^{\eta}$ vectors; in other words, we are guaranteed to find amongst our weight vector $W$ a collection $Q_{m}$ of $v_{m}\simeq m^{\eta}$ vectors that are all at an angle at most $2\epsilon_{m}$ apart. Let us now apply Theorem 2.4 by picking $n=v_{m}$ and $\alpha=\epsilon_{m}$. We need to see that the terms involved in the bound all converge to as $m\to\infty$.

The contribution of the oracle error $e(v_{m})-e(m)$ goes to zero as $m\to\infty$ by the fact that $\lim_{m\to\infty}e(m)$ exists (it is a decreasing, positive sequence) and that $v_{m}\to\infty$.

Let us now verify that $\delta(m-v_{m},\epsilon_{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 $\epsilon_{m}$. Indeed, let $w_{k}$ be one of such vectors and let $\beta^{\prime}$ be the solution of

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(\beta_{p})-E(\beta)$, where

since now $\beta_{p}$ is a feasible solution for the pruned first layer.

Let us finally bound $E(\beta_{p})-E(\beta)$.

Since $\angle(w_{k},w_{k-1})\leq\epsilon_{m}$, it results from Proposition 2.3 that

It results that removing $|Q_{m}|$ of such vectors incurs an increase of the loss at most $|Q_{m}|\epsilon_{m}\simeq m^{\eta}m^{\frac{\eta-1}{n}}=m^{\eta+\frac{\eta-1}{n}}$. Since we picked $\eta$ such that $\eta+\frac{\eta-1}{n}<0$, this term converges to zero. The proof is finished. $\square$

Appendix C Cartoon of Algorithm

Appendix 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 $\theta_{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.

Appendix E A Disconnection

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.