AN EMPIRICAL ANALYSIS OF DEEP NETWORK LOSS SURFACES
The training of deep neural networks is a high-dimension optimization problem with respect to the loss function of a model. Unfortunately, these functions are of high dimension and non-convex and hence difficult to characterize. In this paper, we empirically investigate the geometry of the loss functions for state-of-the-art networks with multiple stochastic optimization methods. We do this through several experiments that are visualized on polygons to understand how and when these stochastic optimization methods find local minima.
1 INTRODUCTION
Deep neural networks are trained by optimizing an extremely high-dimensional loss function with respect to the weights of the network’s linear layers. The objective function minimized is some measure of the error of the network’s predictions based on these weights compared to training data. This loss function is non-convex and has many local minima. These loss functions are usually minimized using first-order gradient descent (Robbins & Monro, 1951; Polyak, 1964) algorithms such as stochastic gradient descent (SGD) (Bottou, 1991). The success of deep learning critically depends on how well we can minimize this loss function, both in terms of the quality of the local minima found and the time to find them. Understanding the geometry of this loss function and how well optimization algorithms can find good local minima is thus of vital importance.
Several works have theoretically analyzed and characterized the geometry of deep network loss functions. However, to make these analyses tractible, they have relied on simplifications of the network structures, including that the networks are linear (Saxe et al., 2014), or assuming the path and variable independence of the neural networks (Choromanska et al., 2015). Orthogonally, the performance of various gradient descent algorithms has been theoretically characterized (Nesterov, 1983). Again, these analyses make simplifying assumptions, in particular that the loss function is strictly convex, i.e. there is only a single local minimum.
In this work, we empirically investigated the geometry of the real loss functions for state-of-the-art networks and data sets. In addition, we investigated how popular optimization algorithms interact with these real loss surfaces. To do this, we plotted low-dimensional projections of the loss function in subspaces chosen to investigate properties of the local minima selected by different algorithms. We chose these subspaces to address the following questions:
• What types of changes to the optimization procedure result in different local minima? • Do different optimization algorithms find qualitatively different types of local minima?
2 RELATED WORK
2.1 LOSS SURFACES
There have been several attempts to understand the loss surfaces of deep neural networks. Some have studied the critical points of the deep linear neural networks (Baldi, 1989; Baldi & Hornik,
∗Work done during an internship at Janelia Research Campus
1989; Baldi & Lu, 2012). Others further investigated the learning dynamics of the deep linear neural networks (Saxe et al., 2014). More recently, several others have attempted to study the loss surfaces of deep non-linear neural networks (Choromanska et al., 2015; Kawaguchi, 2016; Soudry & Carmon, 2016).
One approach is to analogize the states of neurons as the magnetics dipoles used in spherical spinglass Ising models from statistical physics (Parisi, 2016; Fyodorov & Williams, 2007; Bray & Dean, 2007). Choromanska et al. (2015) attempted to understand the loss function of neural networks through studying the random Gaussian error functions of Ising models. Recent results (Kawaguchi, 2016; Soudry & Carmon, 2016) have provided cursory evidence in agreement with the theory provided by Choromanska et al. (2015) in that they found that that there are no “poor” local minima in neural networks still with strong assumptions.
There is some potential disconnect between these theoretical results and what is found in practice due to several strong assumptions such as the activation of the hidden units and output being independent of the previous hidden units and input data. The work of Dauphin et al. (2014) empirically investigated properties of the critical points of neural network loss functions and demonstrated that their critical points behave similarly to the critical points of random Gaussian error functions in high dimensional space. We will expose further evidence along this trajectory.
2.2 OPTIMIZATION
In practice, the local minima of deep network loss functions are for the most part decent. This implies that we probably do not need to take many precautions to avoid bad local minima in practice. If all local minima are decent, then the task of finding a decent local minimum quickly is reduced to the task of finding any local minimum quickly. From an optimization perspective this implies that solely focusing on designing fast methods are of key importance for training deep networks.
In the literature the common method for measuring performance of optimization methods is to analyze them on nice convex quadratic functions (Polyak, 1964; Broyden, 1970; Nesterov, 1983; Martens, 2010; Erdogdu & Montanari, 2015) even though the problems are applied to non-convex problems. For non-convex problems, however, if two methods converge to different local minima their performance will be dictated on how those methods solve those two convex subproblems. It is challenging to show that one method will beat another without knowledge of the sort of convex subproblems, which is generally not known apriori. What we will explore is whether indeed are some characteristics that can found experimentally. If so, perhaps one could validate where these analytical results are valid or even improve methods for training neural networks.
2.2.1 LEARNING PHASES
One of the interesting empirical observation is that we often observe is that the incremental improvement of optimization methods decreases rapidly even in non-convex problems. This behavior has been discussed as a “transient” phase followed by a “minimization” phase (Sutskever et al., 2013)
where the former finds the neighborhood of a decent local minima and the latter finds the local minima within that neighborhood. The existence of these phases implies that if certain methods are better at different phases one could create novel methods that schedule when to apply each method.
3 EXPERIMENTAL SETUP AND TOOLS
3.1 NETWORK ARCHITECTURES AND DATA SETS
We conducted experiments on three state-of-the-art neural network architectures. Network-inNetwork (NIN) (Lin et al., 2014) and the VGG(Simonyan & Zisserman, 2015) network are feedforward convolutional networks developed for image classification, and have excellent performance on the Imagenet (Russakovsky et al., 2014) and CIFAR10 (Krizhevsky, 2009) data sets. The long short-term memory network (LSTM) (Hochreiter & Schmidhuber, 1997) is a recurrent neural network that has been successful in tasks that take variable-length sequences as input and/or produce variable-length sequences as output, such as speech recognition and image caption generation. These are large networks currently used in many machine vision and learning tasks, and the loss functions minimized by each are highly non-convex.
All results using the feed-forward convolutional networks (NIN and VGG) are on the CIFAR10 image classification data set, while the LSTM was tested on the Penn Treebank next-word prediction data set.
3.2 OPTIMIZATION METHODS
We analyzed the performance of five popular gradient-descent optimization methods for these learning frameworks: Stochastic gradient descent (SGD) (Robbins & Monro, 1951), stochastic gradient descent with momentum (SGDM), RMSprop (Tieleman & Hinton, 2012), Adadelta (Zeiler et al., 2011), and ADAM (Kingma & Ba, 2014). These are all first-order gradient descent algorithms that estimate the gradients based on randomly-grouped minibatches of training examples. One of the major differences between these algorithms is how they select the weight-update step-size at each iteration, with SGD and SGDM using fixed schedules, and RMSprop, Adadelta, and ADAM using adaptive, per-parameter step-sizes. Details are provided in Section A.2.
In addition to these five existing optimization methods, we compare to a new gradient descent method we developed based on the family of Runge Kutta integrators. In our experiments, we tested a second-order Runge-Kutta integrator in combination with SGD (RK2) and in combination with ADAM (ADAM&RK2). Details are provided in Section A.3).
3.3 ANALYSIS METHODS
Several of our empirical analyses are based on the technique of Goodfellow et al. (Goodfellow et al., 2015). They visualize the loss function by projecting it down to one carefully chosen dimension. They plot the value of the loss function along a set of samples along this dimension. The projection space is chosen based on important weight configurations, thus they plot the value of the loss function at linear interpolations between two weight configurations. They perform two such analyses: one in which they interpolate between the initialization weights and the final learned weights, and one in which they interpolate between two sets of final weights, each learned from different initializations.
In this work, we use a similar visualization technique, but choose different low-dimensional subspaces for the projection of the loss function. These subspaces are based on the initial weights as well as the final weights learned using the different optimization algorithms and combinations of them, and are chosen to answer a variety of questions about the loss function and how the different optimization algorithms interact with this loss function. In contrast, Goodfellow et al. only looked at SGDM. In addition, we explore the use of two-dimensional projections of the loss function, allowing us to better visualize the space between local minima. We do this via barycentric and bilinar interpolation for triplets and quartets of points respectively (details in Section A.1).
We refer to the critical points found using these variants of SGD, for which the gradient is approximately 0, as local minima. Our evidence that these are local minima as opposed to saddle points is
similar to that presented in Goodfellow et al. (Goodfellow et al., 2015). If we interpolate beyond the critical point, in this one-dimensional projection, the loss increases (Fig. 10).
3.4 TECHNICAL DETAILS
We used the VGG and NIN implementations from https://github.com/szagoruyko/cifar.torch.git.
The batch size was set to 128 and the number of epochs was set to 200. The learning rate was chosen from the discrete range between [0.2, 0.1, 0.05, 0.01] for SGD and [0.002, 0.001, 0.0005, 0.0001] for adaptive learning methods. We doubled the learning rates when we ran our augmented versions with Runge-Kutta because they required two stochastic gradient computations per epoch. We used batchnormalization and dropout to regularize our networks. All experiments were run on a 6-core Intel(R) Xeon(R) CPU @ 2.40GHz with a TITAN X.
4 EXPERIMENTAL RESULTS
4.1 DIFFERENT OPTIMIZATION METHODS FIND DIFFERENT LOCAL MINIMA
We trained the neural networks described above using each optimization method starting from the same initial weights and with the same minibatching. We computed the value of the loss function for weight vectors interpolated between the initial weights, the final weights for one algorithm, and the final weights for a second algorithm for several pairings of algorithms. The results are shown in the lower triangle of Table 1.
For every pair of optimization algorithms, we observe that the training loss between the final weights for different algorithms shows a sharp increase along the interpolated path. This suggests that each optimization algorithm found a different critical point, despite starting at the same initialization. We investigated the space between other triplets and quadruples of weight vectors (Figure 2 and 3), and even in these projections of the loss function, we still see that the local minima returned by different algorithms are separated by high loss weight parameters.
Deep networks are overparameterized. For example, if we switch all corresponding weights for a pair of nodes in our network, we will obtain effectively the same network, with both the original and permuted networks outputting the same prediction for a given input. To ensure that the weight vectors returned by the different algorithms were functionally different, we compared the outpts of the networks on each example in a validation data set:
dist(θ1, θ2) = √√√√ 1 Ntest Ntest∑ i=1 ‖F (xi, θ1)− F (xi, θ2)‖2,
where θ1 and θ2 are the weights learned by two different optimization algorithms, xi is the input for a validation example, and F (x, θ) is the output of the network for weights θ on input x.
We found that, for all pairs of algorithms, the average distance between the outputs of the networks (Equation 4.1) was approximately 0.16, corresponding to a label disagreement of about 8% (upper triangle of Table 1). Given the generalization error of these networks (approximately 11%, Figure 4), the maximum disagreement we could see was 22%. Thus, these networks disagreed on a large fraction of these test examples – over 13 rd. Thus, the local minima found by different algorithms correspond to effectively different networks, not trivial reparameterizations of the same one.
4.2 DIFFERENT OPTIMIZATION ALGORITHMS FIND DIFFERENT TYPES OF LOCAL MINIMA
Next, we investigated whether the local minima found by the different optimization algorithms had distinguishing properties. To do this, we trained the networks with each optimization algorithm using different initial parameters. We then compared differences between runs of the same algorithm but different initializations to differences between different algorithms.
As shown in Figure 4(a), in terms of training accuracy, we do see some stereotypy for the optima found by different algorithms, with SGD finding local minima with the lowest training accuracy and ADAM, Rmsprop, and Adadelta finding local minima with the highest training accuracy. However, this could be attributed to SGD’s asymtotically slow convergence near local minima due to the gradient diminishing near extrema. Despite this limitation, Figure 4(b) shows that the generalization accuracy of these different local minima on validation data was not significantly different between algorithms. We also did not see a relationship between the weight initialization and the validation accuracy. Thus, while these algorithms fall into different local minima, they are not different in terms of their final quality.
We visualized the loss surface around each of the local minima for the multiple runs. To do this, we plotted the value of the loss function between the initial and final weights for each algorithm (Figure 5(a,c)) for each run of the algorithm from a different initialization. In addition, we plotted
the value of the loss function between the final weights for selected pairs of algorithms for each run (Figure 5(b,d)). We see that the surfaces look strikingly similar for different runs of the same algorithm, but characteristically different for different algorithms. Thus, we found evidence that the different algorithms land in qualitatively different types of local minima.
In particular, we see in Figure 5(a,c) that the size of the basins around the local minima found by ADAM and ADAM&RK2 are larger than those found by SGD and RK2, i.e. that the training loss is small for a wider range of α values. This is a relative measure, and the magnitude of the change in the weight vector is ∆α‖θ1 − θ0‖ for a change of size ∆α, where θ0 is the initial weight vector θ1 is the result found by a given optimization algorithm. In Figure 6, we repeat this analysis, instead showing the loss as a function of the absolute distance in parameter space:
θ(λ) = θ1 + λ θ0 − θ1 ‖θ0 − θ1‖
(1)
We again see that the size of the basin around the local minima varies by optimization algorithm. Note that we evaluate the loss for weight vectors beyond the initial configuration, which had a loss of 2.4.
4.3 ANALYZING LEARNING AFTER “TRANSIENT PERIOD”
Recall that, during optimization, it has been observed that there is a short “transient” phase when the loss decreases rapidly and a “minimization” phase in which the loss decreases slowly (Section 2.2.1 and Figure 1). In this set of experiments, we investigated the effects of switching from one type of optimization method to another at various points during training, in particular at late stages of training when it is thought that a local minimum has been chosen and is only being localized. We switched from one optimization method to another 25%, 50%, and 75% of the way through training. The results are plotted in Figure 7d. We emphasize that we are not switching methods to improve performance, but rather to investigate the shape of the loss function in regions explored during the “minimization” phase of optimization.
We found that, regardless of how late we switch optimization algorithms, as shown in the rightmost column of Figure 7, the local minima found were all different. This directly disagrees with the notion that the local minimum has effectively been chosen before the “minimization” phase, but instead that which local minimum is found is still in flux this late in optimization. It appears that this switch from one local minimum to another happens almost immediately after the optimization method switches, with the training accuracy jumping to the characteristic accuracy for the given method within a few epochs (Figure 7, left column). Interestingly, we also see the distance between the initial and current weight vectors changes drastically after switching from one optimization
method to another, and that this distance is characteristic per algorithm (Figure 7, middle column). While distance increases with training epoch for any single optimization method, it actually starts to decrease when switching from ADAM to SGD.
4.4 EFFECTS OF BATCH-NORMALIZATION
To understand how batch normalization affects the types of local minima found, we performed a set of experiments comparing loss surfaces near local minima found with and without batch normal-
ization for each of the optimization methods. We visualized the surface near these local minima by interpolating between the initial weights and the final weights as well as between pairs of final weights found with different algorithms.
We observed clear qualitative differences between optimization with (Figure 5) and without (Figure 8) batch normalization. We see that, without batch normalization, the quality of local minimum found by a given algorithm is much more dependent on the initialization. In addition, the surfaces between different local minima are more complex in appearance: with batch normalization we see sharp unimodal jumps in performance but without batch normalization we obtain wide bumpy shapes that aren’t necessarily unimodal.
The neural networks are typically initialized with very small parameter values (Glorot & Bengio, 2010; He et al., 2015). Instead, we trained NIN with exotic intializations such as initial parameters drawn from N (−10.0, 0.01) or N (−1.0, 1.0) and observe the loss surface behaviours. The details of results are discussed in Appendix A.5.
5 CONCLUSIONS
In this work, we performed a series of empirical analyses to understand the geometry of the loss functions corresponding to deep neural networks, and how different optimization methods minimize this loss to answer the two questions posed in the introduction.
What types of changes to the optimization procedure result in different local minima?
We found that every type of change to the optimization procedure we tested resulted in a different local minimum. Different local minima were found using the different optimization algorithms from the same initialization (Section 4.1). Even switching the optimization algorithm to another very late in optimization – during the slow “mimimization” portion of learning – resulted in a different local minimum (Section 4.3). The quality of the local minima found, in terms of training and generalization error, is similar. These different local minima were not equivalent, and made mistakes on different test examples (Section 4.1). Thus, they were not trivially different local minima, as would occur if nodes in internal layers of the network were permuted. We observed that the quality of these local minima was only consistently good when we used batch normalization for regularization. Without batch normalization, the quality of the critical points found depended on the initialization, and some solutions found were not as good as others. Our observations are in contrast to the conclusions of Goodfellow et al., i.e. that local minima are not a problem in deep learning because, in the region of the loss function explored by SGD algorithms, the loss function is well-behaved (Goodfellow et al., 2015). Instead, our observations are more consistent with the explanation that the local minima found by popular SGD optimization methods are almost all good (Choromanska et al., 2015; Kawaguchi, 2016; Soudry & Carmon, 2016).
Do different optimization algorithms find qualitatively different types of local minima?
Interestingly, we found that, while the local minima found by the same optimization algorithm from different initializations were different, the shape of the loss function around these local minima was strikingly similar, and was a characteristic of the optimization algorithm. In particular, we found that the size of the basin around ADAM-based optimization was larger than that around vanilla SGD (Section 4.2). A large basin is related to a large margin, as small changes in the weight vector will not affect the training error, and perhaps could have some implications for generalization error. In our experiments, however, we did not observe better generalization error for ADAM than SGD. Questions for potential future research are why the shapes of the loss functions around different local minima found by the same algorithm are so similar, and what the practical implications of this are.
A SUPPLEMANTARY MATERIALS
A.1 3D VISUALIZATION
Goodfellow et al. (2015) introduced the idea of visualizing 1D subspace of the loss surface between the parameters. Here, we propose to visualize loss surface in 3D space through interpolating over three and four vertices.
Linear Interpolation Given two parameters θ1 and θ2, θi = αθ1 + (1− α)θ2, ∀α ∈ [0, 1]. (2)
Bilinear Interpolation Given four parameters θ0, θ1, θ2, and θ3, φi = αθ1 + (1− α)θ2 (3) ϕi = αθ3 + (1− α)θ4 (4) θj = βφi + (1− β)ϕi (5)
for all α ∈ [0, 1] and β ∈ [0, 1].
Barycentric Interpolation Given four parameters θ0, θ1, and θ2, let d1 = θ1 − θ0 and d2 = θ2 − θ0. Then, the formulation of the interpolation is
φi = αd1 + θ0 (6) ϕi = αd2 + θ0 (7) θj = βφi + (1− β)ϕi (8)
for all α ∈ [0, 1] and β ∈ [0, 1].
A.2 OPTIMIZATION METHODS
A.2.1 STOCHASTIC GRADIENT DESCENT
In many deep learning applications both the number of parameters and quantity of input data points can be quite large. This makes the full evaluation of U(θ) be prohibitively expensive. A standard technique for aleviating computational loadis to apply an stochastic approximation to the gradient Robbins & Monro (1951). More precisely, one approximates U by a subset of n data points, denoted by {σj}Nj=1 at each timestep:
Un(θ) = 1
n n∑ j=1 `(θ,xσj ) ' 1 N N∑ i=1 `(θ,xi) = U(θ) (9)
Of course this approximation also carries over to the gradient, which is of vital importance to optimization techniques:
∇Un(θ) = 1 n n∑ j=1 ∇`(θ,xσj ) ' ∇U(θ) (10)
This method is what is commonly called Stochastic Gradient Descent or SGD. So long as the data is distributed nicely the approximation error of Un should be sufficiently small such that not only will SGD still behave like normal GD , but it’s wall clock time for to converge should be significantly lower as well.
Usually one uses the stochastic gradient rather than the true gradient, but the inherent noisiness must be kept in mind. In what follows we will always mean the stochastic gradient.
A.2.2 MOMENTUM
In order to aleviate both noise in the input data as well as noise from stochasticity used in computing quantities one often maintains history of previous evaluations. In order to only require one extra variable one usually stores variables of the form
E[F ]t = αFt + βE[F ]t−1. (11)
where Ft is some value changing over time and E[F ]t is the averaged quantity.
An easy scheme to apply this method to is to compute a rolling weighted average of gradients such as E[g]t = (1− α)gt + αE[g]t−1 but there will be other uses in the future.
A.2.3 PERTINENT METHODS
With the aforementioned tools there are a variety of methods that can be constructed. We choose to view these algorithms as implementations of Explicit Euler on a variety of different vector fields to remove the ambiguity between η and gt. We therefore can define a method by the vector field Xt that explicit Euler is applied to with a single η that is never changed.
SGD with Momentum (SGDM) By simply applying momentum to gt one obtains this stabilized stochastic version of gradient descent:
Xt = −E[g]t. (12) This is the most fundamental method that is used in practice and the basis for everything that follows.
Adagrad Adagrad rescalesXt by summing up the sqaures of all previous gradients in a coefficientwise fashion:
Xt = − gt√∑t
i=1 g 2 i +
. (13)
Here is simply set to some small positive value to prevent division-by-zero. In the future we will neglect this term in denominators because it is always necessary.
The concept is to accentuate variations in gt, but because the denominator is monotonically nondecreasing over time this method is doomed to retard its own progress over time. The denominator can also be seen as a form of momentum where α and β are both set to 1.
Rmsprop A simple generalization of ADAGrad is to simply allow for α and β to be changed from 1. In particular one usually chooses a β less than 1, and presumably α = 1− β. Thus one arrives at a method where the effects of the distance history are diminished:
Xt =− gt√ E[g2]t . (14)
Adadelta Adadelta adds another term to RMSprop in order to guarantee that the magnitude of X is balanced with gt Zeiler et al. (2011). More precisely it maintains
Xt√ E[X 2t ] = − gt√ E[g2t ]
(15)
which results in the following vector field: Xt =− √ E[X 2t ]√ E[g2t ] gt. (16)
and η is set to 1.
ADAM By applying momentum to both gt and g2t one arrives at what is called ADAM. This is often considered a combination of SGDM + RMSprop,
Xt = ct E[g]t√ E[g2]t . (17)
ct =
√ 1−βt2
1−βt1 is the initialization bias correction term with β1, β2 ∈ [0, 1) being the β parameters
used in momentum for g and g2 respectively. Initialization bias is caused by the history of the momentum variable being initially set to zero.
A.3 RUNGE KUTTA
Runge-Kutta methods Butcher (1963) are a broad class of numerical integrators categorized by their truncation error. Because the ordinary differential equations Runge-Kutta methods solve generalize gradient descent, our augmentation is quite straightforward. Although our method applies to all explicit Runge-Kutta methods we will only describe second order methods for simplicity.
The general form of second-order explicit Runge-Kutta on a time-independent vector field is
θt+1 = θt + (a1k1 + a2k2)h (18) k1 = X (θt) (19) k2 = X (θt + q1hk1) (20)
where a1, a2, and q1 are parameters that define a given Runge-Kutta method. Table 3 refers to the parameters used for the different Runge-Kutta variants we use in our experiments.
A.3.1 AUGMENTING OPTIMIZATION WITH RUNGE KUTTA
For a given timestep, explicit integrators can be seen as a morphism over vector fields X → X̄ h. For a gradient gt = ∇θU we can solve a modified RK2 gradient ḡt in the following fashion:
θt+1 = θt + ḡth = Advect rk2 g (θ, h) (21)
rearranged with respect to ḡt
ḡt = Advectrk2g (θ, h)− θt
h (22)
= θt + (a1k1 + a2k2)h− θt
h (23)
= (a1k1 + a2k2). (24) If we simply substitute the gradient gt with ḡt one obtains an RK2-augmented optimization technique.
A.4 EXPERIMENTS WITH RUNGE-KUTTA INTEGRATOR
The results in Figure 9 illustrates that, with the exception of the Midpoint method, stochastic RungeKutta methods outperform SGD. “SGD x2” is the stochastic gradient descent with twice of the learning rate of “SGD”. From the figure, we observe that the Runge-Kutta methods perform even better with half the number of gradient computed by SGD. The reason is because SGD has the accumulated truncated error of O(h) while second-order Runge-Kutta methods have the accumulated truncated error of O(h2).
Unfortunately, ADAM outperforms ADAM+RK2 methods. We speculate that this is because the way how ADAM’s renormalization of input gradients in conjunction with momentum eliminates the value added by using our RK-based descent directions.
A.5 EFFECTS OF BATCH-NORMALIZATION AND EXTREME INITIALIZATIONS
The neural networks are typically initialized with very small parameter values (Glorot & Bengio, 2010; He et al., 2015). Instead, we trained NIN with exotic intializations such as initial parameters drawn from N (−10.0, 0.01) or N (−1.0, 1.0) and observe the loss surface behaviours. The results are shown in Figure 11. We can see that NIN without BN does not train at all with any of these initializations. Swirszcz et al. (2016) mentioned that bad performance of neural networks trained with these initializations are due to finding a bad local minima. However, we see that loss surface region around these initializations are plateau 1 rather than a bad local minima as shown in Figure 11b. On
1We used same initializations as (Swirszcz et al., 2016) but we trained different neural networks with SGD on a different dataset. We used NIN and CIFAR10 and Swirszcz et al. (2016) used smaller neural network and MNIST.
the other hand, NIN with BN does train slowly over time but finds a local minima. This implies that BN redeems the ill-posed loss surface (plateau region). Nevertheless, the local minima it found was not good as when the parameters were initialized with small values. However, it is not totally clear whether this is due to difficulty of training or due to falling in a bad local minima.
A.6 SWITCHING OPTIMIZATION METHODS
The paper is dedicated to better understanding the optimization landscape in deep learning, in particular when explored with different optimization algorithms, and thus it also characterizes the behavior of these algorithms. It heavily re-uses the approach of Goodfellow et al. (2015). I find it hard to understand the contributions of the paper, for example: is it surprising that different algorithms reach different solutions when starting from the same initialization? It would be useful if the authors build such basic intuition in the paper. I also did not receive a clear answer to the question I posed to reviewers regarding clarifying how does the findings of the paper can contribute to future works on optimization in deep learning. And this is what I find fundamentally missing. So for example, there are probably plenty of ways to modify approach of Goodfellow et al. (2015), and similar works, and come up with interesting visualization methods for deep learning - but the question is: how is this helpful in terms of designing better algorithms, gaining more intuition how the optimization surface looks like in general, etc.? This is an interesting paper, though I am fairly confident it is a better fit for the journal than this conference. 
It would be interesting and instructive, even for sanity check, to plot the eigenspectra of the solutions recovered by the algorithms to see the order of critical points recovered.