arXiv:1711.06959v1 [stat.ML] 19 Nov 2017

BPGrad: Towards Global Optimality in Deep Learning via Branch and Pruning
Ziming Zhang∗
Mitsubishi Electric Research Laboratories
201 Broadway, Cambridge, MA 02139-1955

Yuanwei Wu∗, Guanghui Wang
EECS, The University of Kansas
1450 Jayhawk Blvd., Lawrence, KS 66045

zzhang@merl.com

{y262w558, ghwang}@ku.edu

Abstract
Understanding the global optimality in deep learning
(DL) has been attracting more and more attention recently.
Conventional DL solvers, however, have not been developed
intentionally to seek for such global optimality. In this paper we propose a novel approximation algorithm, BPGrad,
towards optimizing deep models globally via branch and
pruning. Our BPGrad algorithm is based on the assumption
of Lipschitz continuity in DL, and as a result it can adaptively determine the step size for current gradient given the
history of previous updates, wherein theoretically no smaller
steps can achieve the global optimality. We prove that, by repeating such branch-and-pruning procedure, we can locate
the global optimality within finite iterations. Empirically an
efficient solver based on BPGrad for DL is proposed as well,
and it outperforms conventional DL solvers such as Adagrad, Adadelta, RMSProp, and Adam in the tasks of object
recognition, detection, and segmentation.

Figure 1. Illustration of how BPGrad works, where each black dot denotes
the solution at each iterations (i.e. branch), directed dotted lines denote the
current gradients, and red dotted circles denote the regions wherein there
should be no solutions achieving global optimality (i.e. pruning). BPGrad
can automatically estimate the scales of these regions based on the function
evaluation of solutions and the Lipschitz continuity assumption.

Global optimality is always desirable and preferred in optimization. Locating global optimality in deep learning, however, is extremely challenging due to its high non-convexity,
and thus no conventional DL solvers, e.g. stochastic gradient descent (SGD) [2], Adagrad [7], Adadelta [37], RMSProp [33] and Adam [18], is intentionally developed for
this purpose, to our best knowledge. Alternatively different
regularization techniques are applied to smooth the objective
functions in DL so that the solvers can converge to some geometrically wider and flatter regions in the parameter space
where good model solutions may exist [39, 4, 40]. But these
solutions may not necessarily be the global optimum.
Inspired by the techniques in global optimization of nonconvex functions, we propose a novel approximation algorithm, BPGrad, which has the ability of locating global
optimality in DL via branch and pruning (BP). BP [29] is
a well-known algorithm developed for searching for global
solutions for nonconvex optimization problems. Its basic
idea is to effectively and gradually shrink the gap between
the lower and upper bounds of global optimum by efficiently
branching and pruning the parameter space. Fig. 1 illustrates
the optimization procedure in BPGrad.
In order to branch and prune the space we assume that the
objective functions in DL are Lipschitz continuous [8], or can
be approximated by Lipschitz functions. This is motivated
by the facts that (1) Lipschitz continuity provides a natural
way to estimate the lower and upper bounds of the global
optimum (see Sec. 2.3.1) used in BP, and (2) it can also
serve as regularization, if needed, to smoothen the objective
functions so that the returned solutions can generalize well.

1. Introduction
Deep learning (DL) has been demonstrated successfully
in many different research areas such as image classification [20], speech recognition [16] and natural language processing [32]. In general, its empirical success stems mainly
from better network architectures [15], larger mount of training data [6], and better learning algorithms [12].
However, theoretical understanding of DL for its success in applications still remains elusive. Very recently researchers start to understand DL from the perspective of
optimization such as the optimality of the learned models
[13, 14, 36]. It has been proved that under certain (very
restrictive) conditions the critical points learned for the deep
models actually achieve global optimality, even though the
optimization in deep learning is highly nonconvex. These
theoretical results may partially explain why such deep models work well in practice.
∗ Joint

first authors for the paper.

1

In Fig. 2 we illustrate
the functionality of Lipschitz
continuity as regularization,
where the noisy narrower but
deeper valley is smoothed
out, while the wider but shalFigure 2. Illustration of Liplower valley is preserved.
schitz continuity as regularizaSuch regularization behavior
tion (red) to smoothen a function
can prevent algorithms from
(blue).
being stuck in bad local minima. Also this is advocated and demonstrated to be crucial
in order to achieve good generalization of learned DL models in several recent works such as [4]. In this sense, our
BPGrad algorithm/solver essentially aims to locate global
optimality in the smoothed objective functions for DL.
Further BPGrad can generate solutions along the directions of gradients (i.e. branch) based on the estimated regions
wherein no global optimum should exist theoretically (i.e.
pruning), and by repeating such branch-and-pruning procedure BPGrad can locate global optimum. Empirically the
high demand of computation as well as footprint in memory
for running BPGrad inspires us to develop an efficient DL
solver to approximate BPGrad towards global optimization.
Contributions: The main contributions of our work are:
C1. We propose a novel approximation algorithm, BPGrad,
which is intent on locating global optimum in DL. To
our best knowledge, our approach is the first algorithmic attempt towards global optimization in DL.
C2. Theoretically we prove that BPGrad can converge to
global optimality within finite iterations.
C3. Empirically we propose a novel and efficient DL solver
based on BPGrad to reduce the requirement of computation as well as footprint in memory. We provide both
theoretical and empirical justification for our solver towards preserving the theoretical properties of BPGrad.
We demonstrate that our solver outperforms conventional DL solvers in the applications of object recognition, detection, and segmentation.

1.1. Related Work
Global Optimality in DL: The empirical loss minimization
problem in learning deep models is highly dimensional and
nonconvex with potentially numerous local minima and saddle points. Blum and Rivest [1] showed that it is difficult
to find the global optima because in the worst case even
learning a simple 3-node neural network is NP-complete.
In spite of the difficulties in optimizing deep models, researchers have attempted to provide empirical as well as
theoretical justification for the success of these models w.r.t.
global optimality in learning. Zhang et al. [38] empirically
demonstrated that sufficiently over-parametrized networks
trained with stochastic gradient descent can reach global
optimality. Choromanska et al. [5] studied the loss surface

of multilayer networks using spin-glass model and showed
that for many large-size decoupled networks, there exists a
band with many local optima, whose objective values are
small and close to that of a global optimum. Brutzkus and
Globerson [3] showed that gradient descent converges to
the global optimum in polynomial time on a shallow neural
network with one hidden layer and a convolutional structure
and a ReLU activation function. Kawaguchi [17] proved
that the error landscape does not have bad local minima
in the optimization of linear deep neural networks. Yun et
al. [36] extended these results and proposed sufficient and
necessary conditions for a critical point to be a global minimum. Haeffele and Vidal [13] suggested that it is critical
to balance the degrees of positive homogeneity between the
network mapping and the regularization function to prevent
non-optimal local minima in the loss surface of neural networks. Nguyen and Hein [27] argued that almost all local
minima are global optimal in fully connected wide neural
networks, whose number of hidden neurons of one layer is
larger than that of training points. Soudry and Carmon [30]
employed smoothed analysis techniques to provide theoretical guarantee that the highly nonconvex loss functions in
multilayer networks can be easily optimized using local gradient descent updates. Hand and Voroninski [14] provided
theoretical properties for the problem of enforcing priors provided by generative deep neural networks via empirical risk
minimization by establishing the favorable global geometry.
DL Solvers: SGD [2] is the most widely used DL solver
due to its simplicity, whose learning rate (i.e., step size for
gradient) is predefined. In general, SGD suffers from slow
convergence, and thus its learning rate needs to be carefully
tuned. To improve the efficiency of SGD, several DL solvers
with adaptive learning rates have been proposed, including
Adagrad [7], Adadelta [37], RMSProp [33] and Adam [18].
These solvers integrate the advantages from both stochastic and batch methods where small mini-batches are used
to estimate diagonal second-order information heuristically.
These solvers have the capability of escaping saddle points
and often yield faster convergence empirically.
Specifically, Adagrad is well suited for dealing with
sparse data, as it adapts the learning rate to the parameters, performing smaller updates on frequent parameters
and larger updates on infrequent parameters. However, it
suffers from shrinking on the learning rate, which motivates Adadelta, RMSProp and Adam. Adadelta accumulates
squared gradients to be fixed values rather than over time
in Adagrad, RMSProp updates the parameters based on the
rescaled gradients, and Adam does so based on the estimated
mean and variance of the gradients. Very recently, Mukkamala et al. [26] proposed variants of RMSProp and Adagrad
with logarithmic regret bounds.
Convention vs. Ours: Though the properties of global optimality in DL are very attractive, as far as we know, however,

there is no solver developed intentionally to capture such
global optimality so far. To fill this void, we propose our
BPGrad algorithm towards global optimization in DL.
From the optimization perspective, our algorithm shares
similarities with the recent work [25] on global optimization of general Lipschitz functions (not specifically for DL).
In [25] a uniform sampler is utilized to maximize the lower
bound of the maximizer (equivalently minimizing the upper bound of the minimizer) subject to Lipschitz conditions.
Convergence properties w.h.p. are derived. In contrast, our
approach considers estimating both lower and upper bounds
of the global optimum, and employs the gradients as guidance to more effectively sample the parameter space for
pruning. Convergence is proved to show that our algorithm
will terminate within finite iterations.
From the empirical solver perspective, our solver shares
similarities with the recent work [19] on improving SGD
using the feedback from the objective function. Specifically [19] tracks the relative changes in the objective function with a running average, and uses it to adaptively tune
the learning rate in SGD. No theoretical analysis, however, is
provided for justification. In contrast, our solver does use the
feedback from the object function to determine the learning
rate adaptively but based on the rescaled distance between
the feedback and the current lower bound estimation. Theoretical as well as empirical justifications are established.

2. BPGrad Algorithm for Deep Learning
2.1. Key Notation
We denote x ∈ X ⊆ Rd as the parameters in the neural
network, (ω, y) ∈ Ω × Y as a pair of a data sample ω and
its associated label y, φ : Ω × X → Y as the nonconvex
prediction function represented by the network, f as the
objective function for training the network with Lipschitz
constant L ≥ 0, ∇f as the gradient of f over parameters x1 ,
∇f
∇f˜ = k∇f
k2 denotes the normalized gradient (i.e. direction
of the gradient), f ∗ as the global minimum, and k · k2 as the
`2 -norm operator over vectors.
Definition 1 (Lipschitz Continuity [8]). A function f is Lipschitz continuous with Lipschitz constant L on X , if there is
a (necessarily nonnegative) constant L such that
|f (x1 ) − f (x2 )| ≤ Lkx1 − x2 k2 , ∀x1 , x2 ∈ X .

(1)

2.2. Problem Setup
We would like to learn the parameters for a given network
by minimizing the following objective function f :
h
i
min f (x) ≡ E(ω×y)∈Ω×Y L(y, φ(ω, x)) + R(x), (2)
x∈X

where E denotes the expectation over data pairs, L denotes a
loss function (e.g., hinge loss) for measuring the difference
between the ground-truth labels and the predicted labels
given data samples, and R denotes a regularizer over parameters. Particularly we assume that:
F1. f is lower bounded by 0 and upper bounded as well, i.e.
0 ≤ f (x) < +∞, ∀x ∈ X ;
F2. f is differentiable everywhere in the bounded space X ;
F3. f is Lipschitz continuous, or can be approximated by
Lipschitz functions, with constant L ≥ 0.

2.3. Algorithm
2.3.1

Lower & Upper Bound Estimation

Consider the situation where samples x1 , · · · , xt ∈ X exist for evaluation by function f with Lipschitz constant L,
whose global minimum f ∗ is reached by the sample x∗ .
Then based on Eq. 1 and simple algebra, we can obtain
n
o
max f (xi ) − Lkxi − x∗ k2 ≤ f ∗ ≤ min f (xi ).
i=1,··· ,t

i=1,··· ,t

(3)
This provides us a tractable upper bound and an intractable
lower bound, unfortunately, of the global minimum. The
intractability comes from the fact that x∗ is unknown, and
thus makes the lower bound in Eq. 3 unusable empirically.
To address this problem, we propose a novel tractable
estimator, ρ mini=1,··· ,t f (xi ) (0 ≤ ρ < 1). This estimator
intentionally introduces a gap from the upper bound, which
will be shrunk by either decreasing the upper bound or increasing ρ. As proved in Thm. 1 (see Sec. 2.4), when the
parameter space X is fully covered by the samples {xi }, this
estimator will become the lower bound of f ∗ .
In summary, we define our lower and upper bound estimators for the global minimum as ρ mini=1,··· ,t f (xi ) and
mini=1,··· ,t f (xi ), respectively.
2.3.2

Branch & Pruning

Based on our estimators, we propose a novel approximation
algorithm, BPGrad, towards global optimization in DL via
branch and pruning. We show it in Alg. 1 where the predefined constant  ≥ 0 controls the precision of the solution.
Branch: The inner loop in Alg. 1 conducts the branch operation to split the parameter space recursively by sampling.
Towards this goal, we need a mapping between the parameter space and the bounds. Considering the lower bound in
Eq. 3, we propose sampling xt+1 ∈ X based on the previous
samples x1 , · · · , xt ∈ X so that it satisfies
n
o
max f (xi ) − Lkxi − xt+1 k2 ≤ ρ min f (xi ). (4)
i=1,··· ,t

i=1,··· ,t

1 We

assume ∇f 6= 0 w.l.o.g. Empirically we can randomly sample a
non-zero direction for update wherever ∇f = 0.

Note that an equivalent constraint has been used in [25].

Algorithm 1 BPGrad Algorithm for Deep Learning
Input :objective function f with Lipschitz constant L ≥ 0,
precision  ≥ 0
Output :minimizer x∗
Randomly initialize x1 , t ← 1, ρ ← 0;

do
while mini=1,··· ,t f (xi ) ≤ 1−ρ
while ∃xt+1 ∈ X satisfies Eq. 4 do
Compute xt+1 by solving Eq. 5;
t ← t + 1;
end
Increase ρ such that 0 ≤ ρ < 1 still holds;
end
return x∗ = xi∗ where i∗ ∈ arg mini=1,··· ,t f (xi );

2.4. Theoretical Analysis
Theorem 1 (Lower & Upper Bounds). Whenever XR (t) ≡
X holds, the samples generated by Alg. 1 satisfies
ρ min f (xi ) ≤ f ∗ ≤ min f (xi ).
i=1,··· ,t

To improve sampling efficiency for decreasing the objective, we propose sampling along the directions of (stochastic) gradients with small distortion. Though gradients only
encode local structures of (nonconvex) functions in a high dimensional space, they are good indicators for locating local
minima [23, 28]. Specifically, we propose a minimization
problem for generating samples:

 2
min
xt+1 − xt − ηt ∇f˜(xt )
+ γηt2 ,
(5)
xt+1 ∈X ,ηt ≥0
2
n
o
s.t. max f (xi ) − Lkxi − xt+1 k2 ≤ ρ min f (xi ),
i=1,··· ,t

i=1,··· ,t

where γ ≥ 0 is a predefine constant controlling the trade-off
between the distortion and the step size ηt ≥ 0. That is,
under the condition in Eq. 4, the objective in Eq. 5 aims to
generate a sample that has small distortion from an anchor
point, whose step size is small as well due to the locality
property of gradients, along the direction of the gradient.
Note that other reasonable objective functions may also
be utilized here for sampling purpose as long as the condition
in Eq. 4 is satisfied. More efficient sampling objectives will
be investigated in our future work.
Pruning: In fact Eq. 4 specifies that new samples should
be generated outside the union of a set of balls defined by
previous samples. To precisely describe this requirement,
we introduce a new concept of removable solution space in
our work as follows:
Definition 2 (Removable Parameter Space (RPS)). We define the RPS, denoted as XR , as
def
XR (t) = ∪j=1,··· ,t B (xj , rj ) ,

Therefore, when the lower bound estimator is higher than the
global minimum f ∗ , we can safely remove all the points in
RPS without evaluation. However, when it becomes smaller
than f ∗ , we risk missing the global solutions.
To address this issue, we propose the outer loop in Alg. 1
to increase the lower bound for drawing more samples which
may further decrease the upper bound later.

(6)

where B(xj , rj ) = {x | kx − xj k2 < rj , x ∈ X }, ∀j
defines a ball centered at sample xj ∈ X with radius rj =
1
L [f (xj ) − ρ mini=1,··· ,t f (xi )] , ∀j.
RPS specifies a region wherein the function evaluations
of all the points cannot be smaller than the lower bound estimator conditioning on the Lipschitz continuity assumption.

(7)

i=1,··· ,t

Proof. Since f ∗ is the global minimum, it always holds that
f ∗ ≤ mini=1,··· ,t f (xi ). Now suppose that if XR (t) ≡ X
holds, ρ mini=1,··· ,T f (xi ) > f ∗ holds as well. Then there
would exist at least one point (i.e. global minimum) left for
sampling, contradicting the condition of XR (t) ≡ X . We
then complete the proof.
Corollary 1 (Approximation Error Bound). Whenever both

mini=1,··· ,t f (xi ) ≤ 1−ρ
and XR (t) ≡ X hold, it is satisfied that
min f (xi ) − f ∗ ≤ .

(8)

i=1,··· ,t

Theorem 2 (Convergence within Finite Samples). The total
number of samples, T , in Alg. 1 is upper bounded by:

d
2L
VX
T ≤
,
(9)
·
(1 − ρ)fmin
C
d

where VX denotes the volume of the space X , C =

π2
Γ( d
2 +1)

denotes a constant, and fmin = mini=1,··· ,T f (xi ) denotes
the minimum evaluation.
Proof. Given ∀j, ∀t such that 1 ≤ j ≤ t ≤ T − 1, we have


1
kxt+1 − xj k2 ≥
f (xj ) − ρ min f (xi )
(10)
i=1,··· ,t
L
1−ρ
(1 − ρ)fmin
≥
· min f (xi ) ≥
.
i=1,··· ,t
L
L


min
This allows us to generate two balls B xt+1 , (1−ρ)f
2L


(1−ρ)fmin
and B xj ,
so that they have no overlap with
2L
each other. As a result we can generate T balls with radius
min
of (1−ρ)f
and no overlaps, and their accumulated volume
2L
should be no bigger than VX . That is,
VX ≥

T
X
t=1


V B x

(1−ρ)fmin
t,
2L



=C

(1 − ρ)fmin
2L

d
T. (11)

Further using simple algebra we can complete the proof.

Algorithm 2 BPGrad based Solver for Deep Learning
Input :number of evaluations n repeating N times at most,
objective function f with Lipschitz constant L ≥ 0,
momentum 0 ≤ µ ≤ 1
Output :minimizer x∗
t ← 1, v1 ← 0, and randomly initialize x1 ;
for m ← 1 to N do
1
ρ←1− m
;
while t < mn do
f (xt )−ρ mini=1,··· ,t f (xi )
∇f (xt )
vt+1 ← µvt −
· k∇f
L
(xt )k2 ;
xt+1 ← xt + vt+1 ;
t ← t + 1;
end

if mini=1,··· ,t f (xi ) ≤ 1−ρ
holds then Break ;
end
return x∗ = xi∗ where i∗ ∈ arg mini=1,··· ,n f (xi );

3. Approximate DL Solver based on BPGrad
Though the BPGrad algorithm has nice theoretical properties for global optimization, directly applying Alg. 1 to
deep learning will incur the following problems that limit its
empirical usage:
P1. From Thm. 2 we can see that due to the high dimensionality of the parameter space in DL it is impractical
to draw sufficient samples to cover the entire space.
P2. Solving Eq. 5 involves the knowledge of previous samples, which incurs significant amount of both computational and storage burden for deep learning.
P3. Computing f (xt ) and ∇f˜(xt ), ∀xt ∈ X is timeconsuming, especially for large-scale data.
To address problem P1, in practice we manually set the
maximum iterations for both inner and outer loops in Alg. 1.
To address problem P2, we further make some extra assumptions to simplify the branching/sampling procedure
based on Eq. 5 as follows:
A1. Minimizing distortion is much important than minimizing step sizes, i.e. γ  1;
A2. X is sufficiently large where ∃ηt ≥ 0 so that xt+1 =
xt − ηt ∇f˜(xt ) ∈ X \ XR (t) always holds;
A3. ηt ≥ 0 is always sufficiently small for local update.
A4. xt+1 can be sampled only based on xt and ∇f˜(xt ).
By imposing these assumptions upon Eq. 5, we can directly
compute the solution as follows:


1
ηt =
f (xt ) − ρ min f (xi ) .
(12)
i=1,··· ,t
L
To address problem P3, we utilize mini-batches to estimate f (xt ) and ∇f˜(xt ) efficiently in each iteration.

(a) Sampling using Eq. 5

(b) Sampling using Eq. 12

Figure 3. 1D illustration of difference in sampling between (a) using Eq. 5
and (b) using Eq. 12. Here the solid blue lines denote function f , the black
dotted lines denote the sampling paths starting from xt−1 → xt → xt+1 ,
and each big triangle surrounded by blue dotted lines denotes the RPS of
each sample. As we see, (b) suffers from being stuck locally, while (a) can
avoid the locality based on the RPS.

In summary, we list our BPGrad solver in Alg. 2 by modifying Alg. 1 for the sake of fast sampling as well as low
memory footprint in DL, but at the risk of being stuck in
local regions. Fig. 3 illustrates such scenarios in a 1D example. In (b) the sampling method falls into a loop because
it does not consider the history of samples but only current
one. In contrast, the sampling method in (a) is able to keep
generating new samples by avoiding the RPS of previous
samples with more computation and storage, as expected.

3.1. Theoretical Analysis
Theorem 3 (Global Property Preservation). Let xt+1 =
xt − ηt ∇f˜(xt ) where ηt is computed using Eq. 12. Then
xt+1 satisfies Eq. 4 if it holds that
D
E f (x ) − f (x )
i
t
xi − xt , ∇f˜(xt ) ≥
, ∀i = 1, · · · , t, (13)
L
where h·, ·i denotes the inner product between two vectors.
Proof. Based on Eq. 1, Eq. 12, and Eq. 13, we have
kxi − xt+1 k2

D
E 21
= kxi − xt k22 + ηt2 + 2ηt xi − xt , ∇f˜(xt )


1
≥
f (xi ) − ρ min f (xi ) , ∀i = 1, · · · , t, (14)
i=1,··· ,t
L
which is essentially equivalent to Eq. 4 based on algebra.
We then can complete the proof.
Corollary 2. Suppose that a monotonically decreasing sequence {f (xi )}i=1,··· ,t is generated to minimize function f
by sampling using Eq. 12. Then the condition in Eq. 13 can
be rewritten as follows:
D
E
xi − xj , ∇f˜(xj ) ≥ 0, 1 ≤ ∀i < ∀j ≤ t.
(15)
Discussion: Both Thm. 3 and Cor. 2 imply that our solver
prefers sampling the parameter space along a path towards a
single direction, roughly speaking. However, the gradients in
conventional backpropagation have little guarantee to satisfy
Eq. 13 or Eq. 15 due to lack of such constraints in learning.

mnist, L=15, =0

0.3

0

value

0.2
-0.5

0.1
-1
Left-Hand
Right-Hand

0
0

600

1200

1800

2400

-1.5
0

value

t

0.5

500

1000

2400

0
-1
-2
-3

0
0

1800

cifar10, L=50, =0.9

1
Left-Hand
Right-Hand

value

value

1200

# iterations

cifar10, L=50, =0
1

mnist
cifar10

600

# iterations

3.2. Empirical Justification
In this section we discuss the feasibility of the assumptions A1-A4 for reducing computation and storage as well
as preserving the properties towards global optimization in
deep learning. We utilize MatConvNet [35] as our testbed,
and run our solver in Alg. 2 to train the default networks
in MatConvNet for MNIST [21] and CIFAR10 [20], respectively, using the default parameters without explicit mention.
Also we set N = 1, µ = 0, L = 15 for MNIST and L = 50
for CIFAR10 by default. For justification purpose we only
run 4 epochs on each dataset, 600 and 500 iterations per
epoch for MNIST and CIFAR10, respectively. For more
experimental details, please refer to Sec. 4.
Essentially assumption A1 is made to support the other
three to simplify the objective in Eq. 5, and assumption A2
usually holds in deep learning due to its high dimensionality.
Therefore, below we only focus on empirical justification of
assumptions A3 and A4.
0.14
Feasibility of A3: To jus0.12
tify this, we collect ηt ’s
0.1
by running Alg. 2 on both
0.08
datasets, and plot them
0.06
in Fig. 4. Overall these
0.04
0.02
numbers are indeed suffi0
ciently small for local up0
500
1000
1500
2000
2500
# iterations
date based on gradients,
Figure 4. Plots of ηt on MNIST
and ηt decreases with the
and CIFAR10, respectively.
increase of iterations, in
general. This behavior is expected as the objective f is supposed to decrease as well w.r.t. iterations. The value gap
at the beginning on the two datasets is induced mainly by
different L’s.
Feasibility of A4: To justify this, we show some evidences
in Fig. 5, where we plot the left-hand side (LHS) and righthand side (RHS) of Eq. 4 based on xt returned by Alg. 2.
As we see in all the subfigures on the right with µ = 0 the
values on RHS are always no smaller than those on LHS
correspondingly. In contrast, in the remaining subfigures
on the left with µ = 0 (i.e. conventional SGD update) the
values on RHS are always no bigger than those on LHS correspondingly. These observations appear to be robust across
different datasets, and irrelevant to parameter L which determines the radius of balls, i.e. step sizes for gradients. The
momentum parameter µ, which is related to the directions of
gradients for updating models, appear to be the only factor
to make the samples of our solver satisfy Eq. 4. This also

mnist, L=15, =0.9

0.5
Left-Hand
Right-Hand

value

On the other hand, momentum [31] is a well-known technique in deep learning to dampen oscillations in gradients
and accelerate directions of low curvature. Therefore, our
solver in Alg. 2 involves momentum to compensate such
drawbacks in backpropagation for better approximation of
Alg. 1.

1500

2000

# iterations

-4
0

Left-Hand
Right-Hand

500

1000

1500

2000

# iterations

Figure 5. Comparison between LHS and RHS of Eq. 4 based on xt
returned by Alg. 2 using different values for momentum parameter µ.

supports our claims in Thm. 3 and Cor. 2 about the relation
between model update and gradient in order to satisfy Eq. 4.
More evidences have been provided in Sec. 4.1.1. Given
these evidences we hypothesize that assumption A4 may
hold empirically when using sufficiently large values for µ.

4. Experiments
To demonstrate the generalization of our BPGrad solver,
we test it in the applications of object recognition, detection,
and segmentation by training deep convolutional neural networks (CNNs). We utilize MatConvNet as our testbed, and
employ its demo code as well as default network architectures for different tasks. Since our solver has the ability of
determining learning rates adaptively, we compare ours with
another four widely used DL solvers with adaptive learning
rates, namely Adagrad, Adadelta, RMSProp, and Adam. We
tune the parameters in these solvers to achieve their best
performance as we can.

4.1. Object Recognition
4.1.1

MNIST & CIFAR10

The MNIST digital dataset consists of a training set of 60K
images and a test set of 10K images in 10 classes labeled
from 0 to 9, where all images have the resolution of 28 × 28
pixels. The CIFAR-10 dataset consists of a training set
of 50K images and a test set of 10K images in 10 object
classes, where the image resolution is 32 × 32 pixels.
We follow the default implementation to train an individual CNN similar to LeNet-5 [22] on each dataset. For
the details of network architectures please refer to the demo
code. Specifically for all the solvers, we train the networks
for 50 and 100 epochs on MNIST and CIFAR10, respectively, with a mini-batch size 100, weight decay 0.0005, and
momentum 0.9. In addition, we fix the initial weights for
two networks and the feeding order of mini-batches for fair
comparison. The global learning rate is set to 0.001 on
MNIST for Adagrad, RMSProp and Adam. On CIFAR10,

0.06

L=10
L=15
L=20
L=50
L=100
L=500
L=1000

0.04
0.02
0

1

2

3

0.8

0.6

L=10
L=20
L=50
L=100
L=200
L=500
L=1000

0.4

0.2

4

Training objective

Training objective

Training objective

0.1
0.08

6

0.9

5

0.8

cifar10, =0.9

1

1

2

# Epochs

3

4

Validation top1 error

mnist, =0.9

0.12

4

3

2

1

0.6

0.5

0.4
0

# Epochs

0.7

5

10

15

20

0

Epochs

Figure 6. Illustration of robustness of Lipschitz constant L in our solver.

5

10

15

20

Epochs

Figure 8. Comparison on (left) training objectives and (right) validation
0

10

-1

10

-2

10

-3

10

Test top1 error (log scale)

Training objective (log scale)

10

0

10

20

30

40

50

top-1 errors for object recognition using ImageNet ILSVRC2012.

-1

10

-2

10

-3

training
validation

Table 1. Top-1 recognition error (%) on ImageNet ILSVRC2012 dataset.
0

10

Epochs

20

30

40

50

Epochs

2

0.5

0.4

Test top1 error

Training objective

0.45
1.5

1

0.5

0.35
0.3
0.25
0.2
0.15

0
0

20

40

Epochs

60

80

100

Adagrad Adadelta RMSProp Adam BPGrad
49.0
71.6
46.0
70.0
33.0
54.8
76.7
47.2
72.8
44.0

0.1
0

20

40

60

80

100

Epochs

Figure 7. Comparison on (left) training objectives and (right) test top-1
errors for object recognition using (top) MNIST and (bottom) CIFAR10.

the global learning rate is set to 0.001 for RMSProp, but to
0.01 for Adagrad, Adam and Eve [19], and it is reduced to
0.005 and 0.001 at the 31-st and 61-st epoch. Adadelta does
not require the global learning rate.
For our solver, the parameters n and N typically depend
on the numbers of mini-batches and epochs, respectively.
Empirically we find that N = 1 seems to work well, and
thus we use it by default for all the experiments. Accordingly
by default n will be set to the product of the numbers of minibatches and epochs.
Also we find that the parameter L as Lipschitz constant
is quite robust w.r.t. performance, indicating that heavily
tuning this parameter is unnecessary in practice. To demonstrate this, we compare the training objectives of our solver
by varying L in Fig. 6. To highlight the differences, here
we crop and show the results in the first four epochs, but
note that the remaining results have similar behavior. As
we can see on MNIST when L varies from 10 to 100, the
corresponding curves are clustered, similarly on CIFAR10
for L from 50 to 1000. We decide to set L = 15 for MNIST
and L = 50 for CIFAR10, respectively, in our solver.
Next we show the solver comparison results in Fig. 7. To
illustrate the effect of momentum in our solver in terms of
performance, here we plot two variants of our solver with
µ = 0 and µ = 0.9, respectively. As we see our solver with

µ = 0.9 works much better than the counterpart, achieving
lower training objectives as well as lower top-1 error at test
time. This again provides evidence to support the importance
of satisfying Eq. 4 in our solver to search for good solutions
toward global optimality.
Overall, our solver performs best on MNIST and slightly
inferior on CIFAR10 at test time, although in terms of training objective it achieves competitive performance on MNIST
and the best on CIFAR10. We hypothesize that this behavior
comes from the effect of regularization on Lipschitz continuity. However, our solver can decrease the objectives
much faster than all the competitors in the first few epochs.
This observation reflects the superior ability of our solver
in determining adaptive learning rates for gradients. Especially on CIFAR10 we also compare an extra solver Eve
based on our implementation. Eve was proposed in recent
related work [19] that improves Adam with the feedbacks
from the objective function, and tested on CIFAR10 as well.
As we can see, our solver is much more reliable, performing
consistently over epochs.
4.1.2

ImageNet ILSVRC2012 [20]

This dataset contains about 1.28M training images and 50K
validation images among 1000 object classes. Following
the demo code, we train the same AlexNet [20] on it from
the scratch using different solvers. We perform training
for 20 epochs, with a mini-batch size 256, weight decay
0.0005, momentum 0.9, and default learning rates for the
competitors. For our solver we set L = 100.
We show the comparison results in Fig. 8. It is evident
that our solver works the best at both training and test time.
Namely, it converges faster to achieve lower objective as
well as lower top-1 error on validation dataset. In terms of
numbers, ours is 3.2% lower than the second best, RMSProp,
at the 20-th epoch as listed in Table 1.
Based on all the experiments above we conclude that our
solver is suitable to train deep models for object recognition.

Adagrad
RMSProp
Adam
BPGrad

aero
67.5
69.1
68.9
69.4

bike
71.5
75.8
79.9
77.7

bird
60.7
61.5
64.1
66.4

boat
47.1
47.9
56.6
55.1

bottle
28.3
30.2
37.0
37.2

bus
72.7
74.7
77.4
76.1

car
76.7
77.1
77.7
77.7

cat
77.0
79.4
82.5
83.6

chair
34.3
33.2
38.2
38.6

cow
70.2
71.1
71.5
73.8

table
64.0
66.3
64.7
67.4

dog
72.0
74.4
77.6
76.0

horse
74.2
76.3
77.7
81.9

mbike
69.5
69.9
75.0
72.7

persn
64.9
65.1
66.8
66.3

plant
28.8
28.9
30.6
31.0

sheep
57.4
62.9
65.9
64.2

sofa
60.5
62.5
65.1
66.2

train
73.1
73.2
74.4
73.8

tv
61.1
60.8
67.9
64.9

mAP
61.7
63.0
66.0
66.0

10

20

30

40

Table 2. Average precision (AP, %) of object detection on VOC2007 test dataset.
0.45

0.95

0.7

1

0.4

0.6
0.9

0.25
0.2
0.15

Training objective

Training losscls

Training lossbbox

0.3

0.6

0.4

0.2

0.1
0.05
2

4

6

8

10

12

Epochs

0.5
0.4
0.3
0.2
0.1
0

0
0

Validation accuracy

0.8
0.35

0

2

4

6

8

10

12

0.85

0.8

0.75

0.7
0

10

20

30

40

50

0

Epochs

Epochs

50

Epochs

Figure 9. Loss comparison on VOC2007 trainval dataset, including (left)

Figure 10. Segmentation performance comparison using FCN-32s model

the regression loss using bounding boxes and (right) the classification loss.

on VOC2011 training and validation datasets.

4.2. Object Detection
Following Fast RCNN [11] in the demo code, we conduct
the solver comparison on the PASCAL VOC2007 dataset [9]
with 20 object classes using selective search [34] as default
object proposal approach. For all solvers, we train the network for 12 epochs using the 5K images in VOC2007 trainval set and test it using 4.9K images in VOC2007 test set.
We set the weight decay and momentum to 0.0005 and 0.9,
respectively, and use default learning rates for the competitors. We do not compare with Adadelta because we cannot
obtain reasonable performance after heavy parameter tuning.
For our solver we set L = 100 and N = 12.
We show the training comparison in Fig. 9, and test results
in Table 2. Though our training losses are inferior to those of
Adam in this case, our solver works as well as Adam at test
time on average, achieving best AP on 11 out of 20 classes.
This demonstrates the suitability of our solver in training
deep models for object detection.

4.3. Object Segmentation
Following the work [24] for semantic segmentation based
on fully convolutional networks (FCN), we train FCN-32s
with per-pixel multinomial logistic loss and validate it with
the standard metric of mean pixel intersection over union
(IU), pixel accuracy, and mean accuracy. For all the solvers,
we conduct training for 50 epochs with momentum 0 and
weight decay 0.0005 on PASCAL VOC2011 [10] segmentation set. For Adagrad, RMSProp and Adam, we find that the
default parameters are able to achieve the best performance.
For Adadelta, we tune its parameters with  = 10−9 . The
global learning rate for RMSProp is set to 10−5 and 10−4
for both Adagrad and Adam. Adadelta does not require the
global learning rate. For our solver, we set L = 500.
We show the learning curves on training and validation
datasets in Fig. 10, and list the test-time comparison results

Adagrad
Adadelta
RMSProp
Adam
BPGrad

mean IU pixel accuracy mean accuracy
60.8
89.5
77.4
46.6
86.0
54.4
60.5
90.2
71.0
50.9
87.2
66.4
62.4
89.8
79.6

average
75.9
62.3
73.9
68.2
77.3

Table 3. Numerical comparison on semantic segmentation performance
(%) using VOC2011 test dataset at the 50-th epoch.

in Table 3. In this case our solver has very similar learning
behavior as Adagrad, but achieves the best performance at
test time. The smaller fluctuation over epochs on the validation dataset demonstrates again the superior reliability of our
solver, compared with the competitors. Taking these observations into account, we believe that our solver has the ability
of learning robust deep models for object segmentation.

5. Conclusion
In this paper we propose a novel approximation algorithm,
namely BPGrad, towards searching for global optimality in
DL via branch and pruning based on Lipschitz continuity assumption. Our basic idea is to keep generating new samples
from the parameter space (i.e. branch) outside the removable
parameter space (i.e. pruning). Lipschitz continuity not only
provides us a way to estimate the lower and upper bounds
of global optimality, but also serves as regularization to further smooth the objective functions in DL. Theoretically we
prove that under some conditions our BPGrad algorithm can
converge to global optimality within finite iterations. Empirically in order to avoid the high demand of computation as
well as storage for BPGrad in DL, we propose a new efficient
solver. Theoretical and empirical justification on preserving
the properties of BPGrad is provided. We demonstrate the
superiority of our solver to several conventional DL solvers
in object recognition, detection, and segmentation.

References
[1] A. Blum and R. L. Rivest. Training a 3-node neural network is
np-complete. In NIPS, pages 494–501, 1989. 2
[2] L. Bottou, F. E. Curtis, and J. Nocedal. Optimization methods for
large-scale machine learning. arXiv preprint arXiv:1606.04838, 2016.
1, 2
[3] A. Brutzkus and A. Globerson. Globally optimal gradient descent
for a convnet with gaussian inputs. arXiv preprint arXiv:1702.07966,
2017. 2
[4] P. Chaudhari, A. Choromanska, S. Soatto, and Y. LeCun. Entropysgd: Biasing gradient descent into wide valleys. arXiv preprint
arXiv:1611.01838, 2016. 1, 2
[5] A. Choromanska, M. Henaff, M. Mathieu, G. B. Arous, and Y. LeCun.
The loss surfaces of multilayer networks. In AISTATS, pages 192–204,
2015. 2
[6] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei. ImageNet:
A Large-Scale Hierarchical Image Database. In CVPR, 2009. 1
[7] J. Duchi, E. Hazan, and Y. Singer. Adaptive subgradient methods
for online learning and stochastic optimization. JMLR, 12(Jul):2121–
2159, 2011. 1, 2
[8] K. Eriksson, D. Estep, and C. Johnson. Applied Mathematics Body
and Soul: Vol I-III. Springer-Verlag Publishing, 2003. 1, 3
[9] M. Everingham, L. Van Gool, C. K. I. Williams, J. Winn,
and A. Zisserman.
The PASCAL Visual Object Classes
Challenge 2007 (VOC2007) Results.
http://www.pascalnetwork.org/challenges/VOC/voc2007/workshop/index.html.
8
[10] M. Everingham, L. Van Gool, C. K. I. Williams, J. Winn,
and A. Zisserman.
The PASCAL Visual Object Classes
Challenge 2011 (VOC2011) Results.
http://www.pascalnetwork.org/challenges/VOC/voc2011/workshop/index.html.
8
[11] R. Girshick. Fast r-cnn. In CVPR, pages 1440–1448, 2015. 8
[12] P. Goyal, P. Dollár, R. Girshick, P. Noordhuis, L. Wesolowski, A. Kyrola, A. Tulloch, Y. Jia, and K. He. Accurate, large minibatch sgd:
Training imagenet in 1 hour. arXiv preprint arXiv:1706.02677, 2017.
1
[13] B. D. Haeffele and R. Vidal. Global optimality in neural network
training. In CVPR, pages 7331–7339, 2017. 1, 2
[14] P. Hand and V. Voroninski. Global guarantees for enforcing deep
generative priors by empirical risk. arXiv preprint arXiv:1705.07576,
2017. 1, 2
[15] K. He, X. Zhang, S. Ren, and J. Sun. Deep residual learning for image
recognition. In CVPR, pages 770–778, 2016. 1
[16] G. Hinton, L. Deng, D. Yu, G. E. Dahl, A.-r. Mohamed, N. Jaitly,
A. Senior, V. Vanhoucke, P. Nguyen, T. N. Sainath, et al. Deep neural
networks for acoustic modeling in speech recognition: The shared
views of four research groups. IEEE Signal Processing Magazine,
29(6):82–97, 2012. 1
[17] K. Kawaguchi. Deep learning without poor local minima. In NIPS,
pages 586–594, 2016. 2
[18] D. Kingma and J. Ba. Adam: A method for stochastic optimization.
arXiv preprint arXiv:1412.6980, 2014. 1, 2
[19] J. Koushik and H. Hayashi. Improving stochastic gradient descent
with feedback. arXiv preprint arXiv:1611.01505, 2016. 3, 7
[20] A. Krizhevsky, I. Sutskever, and G. E. Hinton. Imagenet classification
with deep convolutional neural networks. In NIPS, pages 1097–1105,
2012. 1, 6, 7
[21] Y. LeCun. The mnist database of handwritten digits. http://yann.
lecun.com/exdb/mnist/, 1998. 6
[22] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner. Gradient-based
learning applied to document recognition. Proceedings of the IEEE,
86(11):2278–2324, 1998. 6

[23] J. D. Lee, M. Simchowitz, M. I. Jordan, and B. Recht. Gradient
descent only converges to minimizers. In COLT, pages 1246–1257,
2016. 4
[24] J. Long, E. Shelhamer, and T. Darrell. Fully convolutional networks
for semantic segmentation. In CVPR, pages 3431–3440, 2015. 8
[25] C. Malherbe and N. Vayatis. Global optimization of lipschitz functions. In ICML, 2017. 3
[26] M. C. Mukkamala and M. Hein. Variants of rmsprop and adagrad
with logarithmic regret bounds. arXiv preprint arXiv:1706.05507,
2017. 2
[27] Q. Nguyen and M. Hein. The loss surface of deep and wide neural
networks. arXiv preprint arXiv:1704.08045, 2017. 2
[28] I. Panageas and G. Piliouras. Gradient descent only converges to
minimizers: Non-isolated critical points and invariant regions. arXiv
preprint arXiv:1605.00405, 2016. 4
[29] D. G. Sotiropoulos and T. N. Grapsa. A branch-and-prune method
for global optimization. In Scientific Computing, Validated Numerics,
Interval Methods, pages 215–226. Springer, 2001. 1
[30] D. Soudry and Y. Carmon. No bad local minima: Data independent training error guarantees for multilayer neural networks. arXiv
preprint arXiv:1605.08361, 2016. 2
[31] I. Sutskever, J. Martens, G. Dahl, and G. Hinton. On the importance
of initialization and momentum in deep learning. In ICML, pages
1139–1147, 2013. 6
[32] I. Sutskever, O. Vinyals, and Q. V. Le. Sequence to sequence learning
with neural networks. In NIPS, pages 3104–3112, 2014. 1
[33] T. Tieleman and G. Hinton. Lecture 6.5—RmsProp: Divide the
gradient by a running average of its recent magnitude. COURSERA:
Neural Networks for Machine Learning, 2012. 1, 2
[34] J. R. Uijlings, K. E. Van De Sande, T. Gevers, and A. W. Smeulders.
Selective search for object recognition. IJCV, 104(2):154–171, 2013.
8
[35] A. Vedaldi and K. Lenc. Matconvnet: Convolutional neural networks
for matlab. In ACM Multimedia, pages 689–692, 2015. 6
[36] C. Yun, S. Sra, and A. Jadbabaie. Global optimality conditions for
deep neural networks. arXiv preprint arXiv:1707.02444, 2017. 1, 2
[37] M. D. Zeiler. Adadelta: an adaptive learning rate method. arXiv
preprint arXiv:1212.5701, 2012. 1, 2
[38] C. Zhang, S. Bengio, M. Hardt, B. Recht, and O. Vinyals. Understanding deep learning requires rethinking generalization. arXiv preprint
arXiv:1611.03530, 2016. 2
[39] S. Zhang, A. E. Choromanska, and Y. LeCun. Deep learning with
elastic averaging sgd. In NIPS, pages 685–693, 2015. 1
[40] Z. Zhang and M. Brand. Convergent block coordinate descent for
training tikhonov regularized deep neural networks. In NIPS, 2017. 1

