Year: 2019

Abstract:
Over-parameterization is ubiquitous nowadays in training neural networks to benefit both optimization in seeking global optima and generalization in reducing prediction error. However, compressive networks are desired in many real world applications and direct training of small networks may be trapped in local optima. In this paper, instead of pruning or distilling over-parameterized models to compressive ones, we propose a new approach based on differential inclusions of inverse scale spaces, that generates a family of models from simple to complex ones by coupling gradient descent and mirror descent to explore model structural sparsity. It has a simple discretization, called the Split Linearized Bregman Iteration (SplitLBI), whose global convergence analysis in deep learning is established that from any initializations, algorithmic iterations converge to a critical point of empirical risks. Experimental evidence shows that SplitLBI may achieve comparable and even better performance than other training algorithms on ResNet-18 in large scale training on ImageNet-2012 dataset etc., while with early stopping it unveils effective subnet architecture with comparable test accuracy to dense models after retraining instead of pruning well-trained ones.

heading: 1 INTRODUCTION
text: The expressive power of deep neural networks comes from the millions of parameters, which are optimized by Stochastic Gradient Descent (SGD) (Bottou, 2010) and variants like Adam (Kingma & Ba, 2015). Remarkably, model over-parameterization helps both optimization and generalization. For optimization, over-parameterization may simplify the landscape of empirical risks toward locating global optima efficiently by gradient descent method (Mei et al., 2018; 2019; Venturi et al., 2018; Allen-Zhu et al., 2018; Du et al., 2018). On the other hand, over-parameterization does not necessarily result in a bad generalization or overfitting (Zhang et al., 2017), especially when some weight-size dependent complexities are controlled (Bartlett, 1997; Bartlett et al., 2017; Golowich et al., 2018; Neyshabur et al., 2019).
However, compressive networks are desired in many real world applications, e.g. robotics, selfdriving cars, and augmented reality. Despite that `1 regularization has been applied to deep learning to enforce the sparsity on weights toward compact, memory efficient networks, it sacrifices some prediction performance (Collins & Kohli, 2014). This is because that the weights learned in neural networks are highly correlated, and `1 regularization on such weights violates the incoherence or irrepresentable conditions needed for sparse model selection (Donoho & Huo, 2001; Tropp, 2004; Zhao & Yu, 2006), leading to spurious selections with poor generalization. On the other hand, `2 regularization is often utilized for correlated weights as some low-pass filtering, sometimes in the form of weight decay (Loshchilov & Hutter, 2019) or early stopping (Yao et al., 2007; Wei et al., 2017). Furthermore, group sparsity regularization (Yuan & Lin, 2006) has also been applied to neural networks, such as finding optimal number of neuron groups (Alvarez & Salzmann, 2016) and exerting good data locality with structured sparsity (Wen et al., 2016; Yoon & Hwang, 2017).
Yet, without the aid of over-parameterization, directly training a compressive model architecture may meet the obstacle of being trapped in local optima in contemporary experience. Alternatively, researchers in practice typically start from training a big model using common task datasets like ImageNet, and then prune or distill such big models to small ones without sacrificing too much of the performance (Jaderberg et al., 2014; Han et al., 2015; Zhu et al., 2017; Zhou et al., 2017; Zhang et al.,
2016; Li et al., 2017; Abbasi-Asl & Yu, 2017; Yang et al., 2018; Arora et al., 2018). In particular, a recent study (Frankle & Carbin, 2019) created the lottery ticket hypothesis based on empirical observations: “dense, randomly-initialized, feed-forward networks contain subnetworks (winning tickets) that – when trained in isolation – reach test accuracy comparable to the original network in a similar number of iterations". How to effectively reduce an over-parameterized model thus becomes the key to compressive deep learning. Yet, Liu et al. (2019) raised a question, is it necessary to fully train a dense, over-parameterized model before finding important structural sparsity?
In this paper, we provide a novel answer by exploiting a dynamic approach to deep learning with structural sparsity. We are able to establish a family of neural networks, from simple to complex, by following regularization paths as solutions of differential inclusions of inverse scale spaces. Our key idea is to design some dynamics that simultaneously exploit over-parameterized models and structural sparsity. To achieve this goal, the original network parameters are lifted to a coupled pair, with one weight set W of parameters following the standard gradient descend to explore the over-parameterized model space, while the other set of parameters learning structure sparsity in an inverse scale space, i.e., structural sparsity set Γ. The large-scale important parameters are learned at a fast speed while the small unimportant ones are learned at a slow speed. The two sets of parameters are coupled in an `2 regularization. The dynamics enjoys a simple discretization, i.e. the Split Linearized Bregman Iteration (SplitLBI), with provable global convergence guarantee shown in this paper. Here, SplitLBI is a natural extension of SGD with structural sparsity exploration: SplitLBI reduces to the standard gradient descent method when the coupling regularization is weak, while it leads to a sparse mirror descent when the coupling is strong.
Critically, SplitLBI enjoys a nice property that important subnet architecture can be rapidly learned via the structural sparsity parameter Γ following the iterative regularization path, without fully training a dense network first. Particularly, the support set of structural sparsity parameter Γ learned in the early stage of this inverse scale space discloses important sparse subnet architectures. Such architectures can be fine-tuned or retrained to achieve comparable test accuracy as the dense, over-parameterized networks. As a result, the structural sparsity parameter Γ may enable us to rapidly find “winning tickets” in early training epochs for the “lottery” of identifying successful subnetworks that bear comparable test accuracy to the dense ones. This point is empirically validated in our experiments.
Historically, the Linearized Bregman Iteration (LBI) was firstly proposed in applied mathematics as iterative regularization paths for image reconstruction and compressed sensing (Osher et al., 2005; Yin et al., 2008), later applied to logistic regression (Shi et al., 2013). The convergence analysis was given for convex problems (Yin et al., 2008; Cai et al., 2009), yet remaining open for non-convex problems met in deep learning. Osher et al. (2016) established statistical model selection consistency for high dimensional linear regression under the same irrepresentable condition as Lasso, later extended to generalized linear models (Huang & Yao, 2018). To relax such conditions, SplitLBI was proposed by Huang et al. (2016) to learn structural sparsity in linear models under weaker conditions than generalized Lasso, that was successfully applied in medical image analysis (Sun et al., 2017) and computer vision (Zhao et al., 2018). In this paper, it is the first time that SplitLBI is exploited to train highly non-convex neural networks with structural sparsity, together with a global convergence analysis based on the Kurdyka-Łojasiewicz framework Łojasiewicz (1963).
Contributions. (1) SplitLBI, as an extension of SGD, is applied to deep learning by exploring both over-parameterized models and structural sparsity in the inverse scale space. (2) Global convergence of SplitLBI in such a nonconvex optimization is established based on the KurdykaŁojasiewicz framework, that the whole iterative sequence converges to a critical point of the empirical loss function from arbitrary initializations. (3) Stochastic variants of SplitLBI demonstrate the comparable and even better performance than other training algorithms on ResNet-18 in large scale training such as ImageNet-2012, among other datasets, together with additional structural sparsity in successful models for interpretability. (4) Structural sparsity parameters in SplitLBI provide important information about subnetwork architecture with comparable or even better accuracies than dense models before and after retraining -- SplitLBI with early stopping can provide fast “winning tickets” without fully training dense, over-parameterized models.

heading: 2 METHODOLOGY
text: The supervised learning task learns a mapping ΦW : X → Y , from input space X to output space Y , with a parameter W such as weights in neural networks, by minimizing certain loss functions on training samples L̂n(W ) = 1n ∑n i=1 `(yi,ΦW (xi)). For example, a neural network of l-layer is defined as ΦW (x) = σl(W lσl−1(W l−1 · · ·σ1(W 1x)), where W = {W i}li=1, σi is the nonlinear activation function of the i-th layer.
Differential Inclusion of Inverse Scale Space. Consider the following dynamics,
Ẇt κ = −∇W L̄ (Wt,Γt) (1a)
V̇t = −∇ΓL̄ (Wt,Γt) (1b) Vt ∈ ∂Ω̄(Γt) (1c)
where V is a sub-gradient of Ω̄(Γ) := Ωλ(Γ) + 12κ‖Γ‖ 2 for some sparsity-enforced regularization Ωλ(Γ) = λΩ1(Γ) (λ ∈ R+) such as Lasso or group Lasso penalties for Ω1(Γ), κ > 0 is a damping parameter such that the solution path is continuous, and the augmented loss function is
L̄ (W,Γ) = L̂n (W ) + 1
2ν ‖W − Γ‖22, (2)
with ν > 0 controlling the gap admitted between W and Γ. Compared to the original loss function L̂n (W ), the L̄ (W,Γ) additionally adopt the variable splitting strategy, by lifting the original neural network parameter W to (W,Γ) with Γ modeling the structural sparsity of W . For simplicity, we assumed L̄ is differentiable with respect to W here, otherwise the gradient in Eq. (1a) is understood as subgradient and the equation becomes an inclusion.
The differential inclusion system (1), called Split Inverse Scale Space (SplitISS), can be understood as a gradient descent flow of Wt in the proximity of Γt and a mirror descent flow (Nemirovski & Yudin, 1983) of Γt associated with a sparsity enforcement penalty Ω̄. In mirror descent flow, gradient descent goes on the dual space consisting of sub-gradients Vt, driving the flow in sparse primal space of Γt. For a large enough ν, it reduces to the gradient descent method for Wt. Yet the solution path of Γt exhibits the following property in the separation of scales: starting at the zero, important parameters of large scale will be learned fast, popping up to be nonzeros early, while unimportant parameters of small scale will be learned slowly, appearing to be nonzeros late. In fact, taking Ωλ(Γ) = ‖Γ‖1 and κ→∞ for simplicity, Vt as the subgradient of Ω̄t, undergoes a gradient descent flow before reaching the `∞-unit box, which implies that Γt = 0 in this stage. The earlier a component in Vt reaches
the `∞-unit box, the earlier a corresponding component in Γt becomes nonzero and rapidly evolves toward a critical point of L̄ under gradient flow. On the other hand, the Wt follows the gradient descent with a standard `2-regularization. Therefore, Wt closely follows dynamics of Γt whose important parameters are selected. Such a property is called as the inverse scale space in applied mathematics (Burger et al., 2006) and recently was shown to achieve statistical model selection consistency in high dimensional linear regression (Osher et al., 2016) and general linear models (Huang & Yao, 2018), with a reduction of bias as κ increases. In this paper, we shall see that the inverse scale space property still holds empirically for the highly nonconvex neural network training via Eq. (1). For example, Fig. 1 shows a LeNet trained on MNIST by the discretized dynamics, where important sparse filters are selected in early epochs while the popular SGD returns dense filters.
Compared with directly enforcing a penalty function such as `1 or `2 regularization
min W R̂n(W ) := L̂n (W ) +Ωλ (W ) , λ ∈ R+. (3)
SplitISS avoids the parameter correlation problem in over-parameterized models. In fact, a necessary and sufficient condition for Lasso or `1-type sparse model selection is the incoherence or irrepresentable conditions (Tropp (2004); Zhao & Yu (2006)) that are violated for highly correlated weight parameters, leading to spurious discoveries. In contrast, Huang et al. (2018) showed that equipped with such a variable splitting where Γ enjoys an orthogonal design where the restricted Hessian of the augmented loss on Γ is orthogonal, the SplitISS can achieve model selection consistency under weaker conditions than generalized Lasso, relaxing the incoherence or irrepresentable conditions when parameters are highly correlated. For weight parameter W , instead of directly being imposed with `1-sparsity, it adopts `2-regularization in the proximity of the sparse path of Γ that admits simultaneously exploring highly correlated parameters in over-parameterized models and sparsity regularization.
Split Linearized Bregman Iterations. SplitISS admits an extremely simple discrete approximation, using the Euler forward discretization of dynamics (1):
Wk+1 = Wk − καk · ∇W L̄ (Wk,Γk) , (4a) Vk+1 = Vk − αk · ∇ΓL̄ (Wk,Γk) , (4b) Γk+1 = κ · ProxΩλ (Vk+1) , (4c)
where V0 = Γ0 = 0, W0 can be small random numbers such as Gaussian distribution in neural networks, for some complex networks it can be initialized as common setting.The proximal map in Eq. (4c) that controls the sparsity of Γ is given by
ProxΩλ(V ) = arg min Γ
{ 1
2 ‖Γ− V ‖22 + Ωλ (Γ)
} , (5)
We shall call such an iterative procedure as Split Linearized Bregman Iteration (SplitLBI), that was firstly coined in Huang et al. (2016) as an iterative regularization path for sparse modeling in high dimensional statistics. In the application to neural networks, the loss becomes highly non-convex, the SplitLBI returns a sequence of sparse models from simple to complex ones whose global convergence condition to be shown below, while solving Eq. (3) at various levels of λ might not be tractable except for over-parameterized models.
The sparsity-enforcement penalty used in convolutional neural networks can be chosen as follows. Our sparsity framework aims at regularizing the groups of weight parameters using group Lasso penalty
(Yuan & Lin, 2006), Ω1(Γ) = ∑ g ‖Γg‖2, where ‖Γg‖2 = √∑|Γg| i=1 (Γ g i )
2 and |Γg| is the number of weights in Γg. Thus Eq. (4c) has a closed form solution Γg = κ ·max (0, 1− 1/‖V g‖2)V g for the g-th filter. We treat convolutional and fully connected layers in different ways.
(1) For a convolutional layer, Γg = Γg(cin, cout, size) denote the convolutional filters where size denotes the kernel size and cin and cout denote the numbers of input channels and output channels, respectively. When we regard each group as each convolutional filter, g = cout; otherwise for weight sparsity, g can be every element in the filter that reduces to the Lasso.
(2) For a fully connected layer, Γ = Γ(cin, cout) where cin and cout denote the numbers of inputs and outputs of the fully connected layer. Each group g corresponds to each element (i, j), and the group Lasso penalty degenerates to the Lasso penalty.

heading: 3 GLOBAL CONVERGENCE OF SPLITLBI FOR NEURAL NETWORKS
text: We present a theorem that guarantees the global convergence of SplitLBI, i.e. from any intialization, the SplitLBI sequence converges to a critical point of L̄. Our treatment extends the block coordinate descent (BCD) studied in Zeng et al. (2019), with a crucial difference being the mirror descent involved in SplitLBI. Instead of the splitting loss in BCD (Zeng et al., 2019), a new Lyapunov function is developed here to meet the Kurdyka-Łojasiewicz property Łojasiewicz (1963). Xue & Xin (2018) studied convergence of variable splitting method for single hidden layer networks with Gaussian inputs.
Let P := (W,Γ). Following Huang & Yao (2018), the SplitLBI algorithm in Eq. (4a-4c) can be rewritten as the following standard Linearized Bregman Iteration,
Pk+1 = arg min P
{ 〈P − Pk, α∇L̄(Pk)〉+BpkΨ (P, Pk) } , (6)
where
Ψ(P ) = Ωλ(Γ) + 1
2κ ‖P‖22 = Ωλ(Γ) +
1
2κ ‖W‖22 +
1
2κ ‖Γ‖22, (7)
pk ∈ ∂Ψ(Pk), and BqΨ is the Bregman divergence associated with convex function Ψ, defined by
BqΨ(P,Q) := Ψ(P )−Ψ(Q)− 〈q, P −Q〉, for some q ∈ ∂Ψ(Q). (8)
Without loss of generality, consider λ = 1 in the sequel. One can establish the global convergence of SplitLBI under the following assumptions.
Assumption 1. Suppose that: (a) L̂n(W ) = 1n ∑n i=1 `(yi,ΦW (xi)) is continuous differentiable and ∇L̂n is Lipschitz continuous with a positive constant Lip; (b)L̂n(W ) has bounded level sets; (c) L̂n(W ) is lower bounded (without loss of generality, we assume that the lower bound is 0); (d) Ω is a proper lower semi-continuous convex function and has locally bounded subgradients, that is, for every compact set S ⊂ Rn, there exists a constant C > 0 such that for all Γ ∈ S and all g ∈ ∂Ω(Γ), there holds ‖g‖ ≤ C; and (e) the Lyapunov function
F (P,G) := αL̄(W,Γ) +Bg̃Ω(Γ, Γ̃), (9)
is a Kurdyka-Łojasiewicz function on any bounded set, whereBg̃Ω(Γ, Γ̃) := Ω(Γ)−Ω(Γ̃)−〈g̃,Γ− Γ̃〉, Γ̃ ∈ ∂Ω∗(g̃), and Ω∗ is the conjugate of Ω defined as
Ω∗(g) := sup U∈Rn {〈U, g〉 − Ω(U)}.
Remark 1. Assumption 1 (a)-(c) are regular in the analysis of nonconvex algorithm (see, Attouch et al. (2013) for instance), while Assumption 1 (d) is also mild including all Lipschitz continuous convex function over a compact set. Some typical examples satisfying Assumption 1(d) are the `1 norm, group `1 norm, and every continuously differentiable penalties. By Eq. (9) and the definition of conjugate, the Lyapunov function F can be rewritten as follows,
F (W,Γ, g) = αL̄(W,Γ) + Ω(Γ) + Ω∗(g)− 〈Γ, g〉. (10)
Now we are ready to present the main theorem.
Theorem 1. [Global Convergence of SplitLBI] Suppose that Assumption 1 holds. Let (Wk,Γk) be the sequence generated by SplitLBI (Eq. (4a-4c)) with a finite initialization. If
0 < αk = α < 2
κ(Lip+ ν−1) ,
then (Wk,Γk) converges to a critical point of L̄ defined in Eq. (2), and {W k} converges to a critical point of L̂n(W ).
Applying to the neural networks, typical examples are summarized in the following corollary.
Corollary 1. Let {Wk,Γk, gk} be a sequence generated by SLBI (16a-16c) for neural network training where (a) ` is any smooth definable loss function, such as the square loss (t2), exponential loss (et), logistic loss log(1 +e−t), and cross-entropy loss; (b) σi is any smooth definable activation, such as linear activation (t), sigmoid ( 11+e−t ), hyperbolic tangent ( et−e−t et+e−t ), and softplus ( 1 c log(1 + e
ct) for some c > 0) as a smooth approximation of ReLU; (c) Ω is the group Lasso. Then the sequence {Wk} converges to a stationary point of L̂n(W ) under the conditions of Theorem 1.
Proofs of Theorem 1 and Corollary 1 are given in Appendix A.

heading: 4 EXPERIMENTS WITH STOCHASTIC SPLITLBI
text: We begin with some stochastic variants of SplitLBI and implementations, followed by four groups of experiments demonstrating the utilities of weight parameter Wt and structural sparsity parameter Γt in prediction, interpretability, and capturing effective sparse subnetworks.
Batch Split LBI. For neural network training with large datasets, stochastic approximation of the gradients in Split LBI over the mini-batch (X,Y)batcht is adopted to update the parameter W ,
∇̃tW = ∇W L̂n (W ) | (X,Y)batcht . (11)
SplitLBI with momentum (Mom). Inspired by the variants of SGD, the momentum term can be also incorporated to the standard Split LBI that leads to the following updates of W by replacing Eq (4a) with,
vt+1 = τvt + ∇̃W L̄ (Wt,Γt) (12a) Wt+1 = Wt − καvt+1 (12b)
where τ is the momentum factor, empirically setting as 0.9 in default. One immediate application of such stochastic algorithms of SplitLBI is to “boost networks", i.e. growing a network from the null to a complex one by sequentially applying our algorithm on subnets with increasing complexities.
SplitLBI with momentum and weight decay (Mom-Wd). The update formulation is,
vt+1 = τvt + ∇̃W L̄ (Wt,Γt) (13) Wt+1 = Wt − καvt+1 − βWt (14)
where β is set as 1e−4.
Implementation. Various algorithms are evaluated over the various backbones – LeNet (LeCun et al., 2015), AlexNet (Krizhevsky et al., 2012), VGG (Simonyan & Zisserman, 2014), and ResNet (He et al., 2016) etc., respectively. For MNIST and Cifar-10, the default hyper-parameters of Split LBI are κ = 1, ν = 10 and αk is set as 0.1, decreased by 1/10 every 30 epochs. In ImageNet-2012, the Split LBI utilizes κ = 1, ν = 1000, and αk is initially set as 0.1, decays 1/10 every 30 epochs. We set λ = 1 in Eq. (5) by default, unless otherwise specified. On MNIST and Cifar-10, the batch size is set as 128; and for all methods, the batch size of ImageNet 2012 is 256. The standard data augmentation implemented in pytorch is applied to Cifar-10 and ImageNet2012 datasets, as He et al. (2016). The weights of all models are initialized as He et al. (2015). In the following experiments, we define sparsity as percentage of non-zero parameters, i.e. the number of non-zero weights dividing the total number of weights in consideration, that equals to one minus the pruning rate of the network. We also have the reproducible source codes 1.

heading: 4.1 IMAGE CLASSIFICATION
text: In SplitLBI, the weight parameter Wt explores over-parameterized models that can achieve the state-of-the-art performance in large scale training such as ImageNet-2012 classification.
Experimental Design. We compare different variants of SGD and Adam in the experiments. By default, the learning rate of competitors is set as 0.1 for SGD and its variant and 0.001 for Adam and its variants, and gradually decreased by 1/10 every 30 epochs. In particular, we have,
1https://anonymous.4open.science/repository/d22bbbc8-50d5-4e60-b2e8-4ded4e93db63/ Split_LBI_code
SGD: (1) Naive SGD: the standard SGD with batch input. (2) SGD with l1 penalty (Lasso). The l1 norm is applied to penalize the weights of SGD by encouraging the sparsity of learned model, with the regularization parameter of the l1 penalty term being set as 1e−3 (3) SGD with momentum (Mom): we utilize momentum 0.9 in SGD. (4) SGD with momentum and weight decay (Mom-Wd): we set the momentum 0.9 and the standard l2 weight decay with the coefficient weight 1e−4. (5) SGD with Nesterov (Nesterov): the SGD uses nesterov momentum 0.9.
Adam: (1) Naive Adam: it refers to the standard version of Adam. We report the results of several recent variants of Adam, including (2) Adabound, (3) Adagrad, (4) Amsgrad, and (5) Radam.
SplitLBI achieves the state-of-the-art performance on ImageNet-2012, etc. Tab. 1 shows the experimental results on ImageNet-2012, Cifar-10, and MNIST of some classical networks -- LeNet, AlexNet and ResNet. Our SplitLBI variants may achieve comparable or even better performance than SGD variants in 100 epochs, indicating the efficacy in learning dense, over-parameterized models.

heading: 4.2 SPLITLBI LEARNS SPARSE FILTERS FOR IMPROVED INTERPRETATION
text: In SplitLBI, the structural sparsity parameter Γt explores important sub-network architectures that contributes significantly to the loss or error reduction in early training stages. Through the `2-coupling, structural sparsity parameter Γt may guide the weight parameter to explore those sparse models in favour of improved interpretabiity. For example, Fig. 1 visualizes some sparse filters learned by SplitLBI of LeNet-5 trained on MNIST (with κ = 10 and weight decay every 40 epochs), in comparison with dense filters learned by SGD. The activation pattern of such sparse filters favours high order global correlations between pixels of input images. To further reveal the insights of learned patterns of SplitLBI, we visualize the first convolutional layer of ResNet-18 on ImageNet-2012 along the training path of our SplitLBI as in Fig. 2. The left figure compares the training and validation accuracy of SplitLBI and SGD. The right figure compares visualizations of the filters learned by SplitLBI and SGD using Springenberg et al. (2014).
Implementation. To be specific, denote the weights of an l-layer network as {W 1,W 2, · · · ,W l}. For the i−th layer weights W i, denote the j−th channel W ij . Then we compute the gradient of the sum of the feature map computed from each filter W ij with respect to the input image (here a snake image). We further conduct the min-max normalization to the gradient image, and generate the final visualization map. The right figure compares the visualized gradient images of first convolutional layer of 64 filters with 7× 7 receptive fields. We visualize the models parameters at 20 (purple), 40 (green), and 60 (black) epochs, respectively, which corresponds to the bounding boxes in the right figure annotated by the corresponding colors, i.e., purple, green, and black. We order the gradient images produced from 64 filters by the descending order of the magnitude (`2-norm) of filters, i.e., images are ordered from the upper left to the bottom right. For comparison, we also provide the visualized gradient from random initialized weights.
Filters learned by ImageNet prefer to non-semantic texture rather than shape and color. The filters of high norms mostly focus on the texture and shape information, while color information is with the filters of small magnitudes. This phenomenon is in accordance with observation of AbbasiAsl & Yu (2017) that filters mainly of color information can be pruned for saving computational cost. Moreover, among the filters of high magnitudes, most of them capture non-semantic textures while few pursue shapes. This shows that the first convolutional layer of ResNet-18 trained on ImageNet learned non-semantic textures rather than shape to do image classification tasks, in accordance with recent studies (Geirhos et al., 2019). How to enhance the semantic shape invariance learning, is arguably a key to improve the robustness of convolutional neural networks.

heading: 4.3 GLOBAL CONVERGENCE AND STRUCTURAL SPARSITY OF SPLITLBI
text: We conduct ablation studies based on Cifar-10 dataset with VGG-16 and ResNet-56 to evaluate (i) global convergence of L̄; and (ii) the structural sparsity learned by Γt via exploring test accuracies of sparse models obtained by projecting Wt onto the support set of Γt (mask), by varying two key hyper-parameters κ and ν.
Implementation. We choose SplitLBI with momentum and weight decay, since it achieves very good performance on large-scale experiments. Specifically, we have two set of experiments, where each experiment is repeated for 5 times: (1) we fix ν = 100 and vary κ = 1, 2, 5, 10, where sparsity of Γt and validation accuracies of sparse models are shown in top row of Fig. 4. Note that we keep κ · αk = 0.1 in Eq (1a), to make comparable learning rate of each variant, and also consistent with SGD. Thus the learning rate αk will be adjusted by different κ values. (2) we fix κ = 1, and validate the results of SplitLBI with ν = 10, 20, 50, 100, 200, 500, 1000, 2000 in the second row of Fig. 4 with the learning rate αk = 0.1. Moreover, rather than using sparse models associated with Γt, Fig. 6 in Appendix shows the validation accuracies of full models learned by Wt.
SplitLBI converges to Critical Point. Figure 3 shows the curves of training loss (L̂n) and accuracies, with each point representing the average and variance bar over 5 times. As shown, both training loss and training accuracy will converge, which validates our theoretical result in Theorem 1. Besides, larger κ brings in slower convergence, which agrees with the analysis the convergence rate is inversely scale to κ in Lemma A.5.
Sparse subnetworks achieve comparable performance to dense models without fine-tuning or retraining. From the experiments above, the sparsity of Γ grows as κ and ν increase. While large κ may cause a small number of important parameters growing rapidly, large ν will decouple Wt and Γt such that the growth of Wt does not affect Γt that may over-sparsify and deteriorate model accuracies.
Thus a moderate choice of κ and ν is preferred in practice. In all cases, one can see that moderate sparse models can achieve comparable predictive power to dense models, even without fine-tuning or retraining. This shows that the structural sparsity parameter Γt can indeed capture important weight parameter Wt through their coupling.

heading: 4.4 SPLITLBI WITH EARLY STOPPING AND RETRAIN FINDS EFFECTIVE SUB-NETWORKS
text: Equipped with early stopping, Γt in early epochs may learn effective subnetworks (i.e. “winning tickets” (Frankle & Carbin, 2019; Liu et al., 2019)) that achieve comparable or even better performance after retraining than existing pruning strategies by SGD.
Experimental Design. We adopt a comparison baseline as the one-shot pruning strategy in Frankle & Carbin (2019), which firstly trains a dense over-parameterized model by SGD for T = 160 epochs and find the sparse structure by pruning weights or filters (Liu et al., 2019), then secondly retrains the structure from the scratch with T epochs from the same initialization as the first step. For SplitLBI, instead of pruning weights/filters from dense models, we directly utilize the structural sparsity Γt at different training epochs to define the subnet architecture, followed by retrain-from-scratch (Fine-tune is shown in Appendix Sec. D with preliminary results). Experiments are conducted on Cifar-10 dataset where we still use VGG–16, ResNet-50, and ResNet-56 as the networks to make direct comparisons to previous works. SplitLBI uses momentum and weight decay with hyperparameters shown in Tab. 10 in Appendix. In particular, we set λ = 0.1, and 0.05 for VGG-16, and ResNet-56 respectively, since ResNet-56 has less parameters than VGG-16. Furthermore, we introduce another variant of our SplitLBI by using Lasso rather than group lasso penalty for Γt to sparsify the weights of convolutional filters; and the corresponding models are denoted as VGG-16 (Lasso) and ResNet-50 (Lasso). Every experiment is repeated for five times and the results are shown in Fig. 5. Note that in different runs of SplitLBI, the sparsity of Γt slightly varies.
Sparse subnets found by early stopping of SplitLBI achieve remarkably good accuracy after retrain from scratch. In Fig.5 (a-b), sparse filters discovered by Γt at different epochs are compared against the methods of Network Slimming (Liu et al., 2017), Soft Filter Pruning (Yang et al., 2018), Scratch-B, and Scratch-E, whose results are reported from Liu et al. (2019). At similar sparsity levels, SplitLBI can achieve comparable or even better accuracy than competitors, even with sparse architecture learned from very early epochs (e.g. t = 20 or 10). Moreover in Fig.5 (c-d), we can draw the same conclusion for the sparse weights of VGG-16 (Lasso) and ResNet-50 (Lasso), against the results reported in Liu et al. (2019). These results shows that the structural sparsity parameter Γt found by early stopping of SplitLBI already discloses important subnetwork architecture that may achieve remarkably good accuracy after retrain from scratch. Therefore, it is not necessary to fully train a dense model to find a successful sparse subnet architecture with comparable performance to the dense ones -- one can early stop SplitLBI properly where the structural parameter Γt unveils “winning tickets” (Frankle & Carbin, 2019).

heading: 5 CONCLUSION
text: In this paper, a parsimonious deep learning method is proposed based on differential inclusions of inverse scale spaces. Implemented by a variable splitting scheme, such a dynamics system can exploit over-parameterized models and structural sparsity simultaneously. Besides, its simple discretization, i.e., the SplitLBI, has a proven global convergence and hence can be employed to train deep networks. We have experimentally shown that it can achieve the state-of-the-art performance on many datasets including ImageNet-2012, with better interpretability than SGD. What’s more, equipped with early stopping, such a structural sparsity can unveil the ”winning tickets” – the architecture of sub-networks which after re-training can achieve comparable and even better accuracy than original dense networks.

heading: A PROOF OF THEOREM 1
text: First of all, we reformulate Eq. (6) into an equivalent form. Without loss of generality, consider Ω = Ω1 in the sequel.
Denote R(P ) := Ω(Γ), then Eq. (6) can be rewritten as,
Pk+1 = ProxκR(Pk + κ(pk − α∇L̄(Pk))), (15a) pk+1 = pk − κ−1(Pk+1 − Pk + κα∇L̄(Pk)), (15b)
where pk = [0, gk]T ∈ ∂R(Pk) and gk ∈ ∂Ω(Γk). Thus SplitLBI is equivalent to the following iterations,
Wk+1 = Wk − κα∇W L̄(Wk,Γk), (16a) Γk+1 = ProxκΩ(Γk + κ(gk − α∇ΓL̄(Wk,Γk))), (16b) gk+1 = gk − κ−1(Γk+1 − Γk + κα · ∇ΓL̄(Wk,Γk)). (16c)
Exploiting the equivalent reformulation (16a-16c), one can establish the global convergence of (Wk,Γk, gk) based on the Kurdyka-Łojasiewicz framework. In this section, the following extended version of Theorem 1 is actually proved.
Theorem 2. [Global Convergence of SplitLBI] Suppose that Assumption 1 holds. Let (Wk,Γk, gk) be the sequence generated by SplitLBI (Eq. (16a-16c)) with a finite initialization. If
0 < αk = α < 2
κ(Lip+ ν−1) ,
then (Wk,Γk, gk) converges to a critical point of F . Moreover, {(Wk,Γk)} converges to a stationary point of L̄ defined in Eq. 2, and {W k} converges to a stationary point of L̂n(W ).
A.1 KURDYKA-ŁOJASIEWICZ PROPERTY
To introduce the definition of the Kurdyka-Łojasiewicz (KL) property, we need some notions and notations from variational analysis, which can be found in Rockafellar & Wets (1998).
The notion of subdifferential plays a central role in the following definitions. For each x ∈ dom(h) := {x ∈ Rp : h(x) < +∞}, the Fréchet subdifferential of h at x, written ∂̂h(x), is the set of vectors v ∈ Rp which satisfy
lim inf y 6=x,y→x h(y)− h(x)− 〈v,y − x〉 ‖x− y‖ ≥ 0.
When x /∈ dom(h), we set ∂̂h(x) = ∅. The limiting-subdifferential (or simply subdifferential) of h introduced in Mordukhovich (2006), written ∂h(x) at x ∈ dom(h), is defined by
∂h(x) := {v ∈ Rp : ∃xk → x, h(xk)→ h(x), vk ∈ ∂̂h(xk)→ v}. (17)
A necessary (but not sufficient) condition for x ∈ Rp to be a minimizer of h is 0 ∈ ∂h(x). A point that satisfies this inclusion is called limiting-critical or simply critical. The distance between a point x to a subset S of Rp, written dist(x,S), is defined by dist(x,S) = inf{‖x− s‖ : s ∈ S}, where ‖ · ‖ represents the Euclidean norm. Let h : Rp → R ∪ {+∞} be an extended-real-valued function (respectively, h : Rp ⇒ Rq be a point-to-set mapping), its graph is defined by
Graph(h) := {(x, y) ∈ Rp × R : y = h(x)}, (resp. Graph(h) := {(x,y) ∈ Rp × Rq : y ∈ h(x)}),
and its domain by dom(h) := {x ∈ Rp : h(x) < +∞} (resp. dom(h) := {x ∈ Rp : h(x) 6= ∅}). When h is a proper function, i.e., when dom(h) 6= ∅, the set of its global minimizers (possibly empty) is denoted by
arg minh := {x ∈ Rp : h(x) = inf h}.
The KL property (Łojasiewicz, 1963; 1993; Kurdyka, 1998; Bolte et al., 2007a;b) plays a central role in the convergence analysis of nonconvex algorithms (Attouch et al., 2013; Wang et al., 2019). The following definition is adopted from Bolte et al. (2007b). Definition 1. [Kurdyka-Łojasiewicz property] A function h is said to have the Kurdyka-Łojasiewicz (KL) property at ū ∈ dom(∂h) := {v ∈ Rn|∂h(v) 6= ∅}, if there exists a constant η ∈ (0,∞), a neighborhood N of ū and a function φ : [0, η) → R+, which is a concave function that is continuous at 0 and satisfies φ(0) = 0, φ ∈ C1((0, η)), i.e., φ is continuous differentiable on (0, η), and φ′(s) > 0 for all s ∈ (0, η), such that for all u ∈ N ∩ {u ∈ Rn|h(ū) < h(u) < h(ū) + η}, the following inequality holds
φ′(h(u)− h(ū)) · dist(0, ∂h(u)) ≥ 1. (18) If h satisfies the KL property at each point of dom(∂h), h is called a KL function.
KL functions include real analytic functions, semialgebraic functions, tame functions defined in some o-minimal structures (Kurdyka, 1998; Bolte et al., 2007b), continuous subanalytic functions (Bolte et al., 2007a) and locally strongly convex functions. In the following, we provide some important examples that satisfy the Kurdyka-Łojasiewicz property. Definition 2. [Real analytic] A function h with domain an open set U ⊂ R and range the set of either all real or complex numbers, is said to be real analytic at u if the function h may be represented by a convergent power series on some interval of positive radius centered at u: h(x) = ∑∞ j=0 αj(x−u)j , for some {αj} ⊂ R. The function is said to be real analytic on V ⊂ U if it is real analytic at each u ∈ V (Krantz & Parks, 2002, Definition 1.1.5). The real analytic function f over Rp for some positive integer p > 1 can be defined similarly.
According to Krantz & Parks (2002), typical real analytic functions include polynomials, exponential functions, and the logarithm, trigonometric and power functions on any open set of their domains. One can verify whether a multivariable real function h(x) on Rp is analytic by checking the analyticity of g(t) := h(x + ty) for any x,y ∈ Rp. Definition 3. [Semialgebraic]
(a) A set D ⊂ Rp is called semialgebraic (Bochnak et al., 1998) if it can be represented as
D = s⋃ i=1 t⋂ j=1 {x ∈ Rp : Pij(x) = 0, Qij(x) > 0} ,
where Pij , Qij are real polynomial functions for 1 ≤ i ≤ s, 1 ≤ j ≤ t.
(b) A function h : Rp → R ∪ {+∞} (resp. a point-to-set mapping h : Rp ⇒ Rq) is called semialgebraic if its graph Graph(h) is semialgebraic.
According to (Łojasiewicz, 1965; Bochnak et al., 1998) and (Shiota, 1997, I.2.9, page 52), the class of semialgebraic sets are stable under the operation of finite union, finite intersection, Cartesian product or complementation. Some typical examples include polynomial functions, the indicator function of a semialgebraic set, and the Euclidean norm (Bochnak et al., 1998, page 26).
A.2 KL PROPERTY IN DEEP LEARNING AND PROOF OF COROLLARY 1
In the following, we consider the deep neural network training problem. Consider a l-layer feedforward neural network including l − 1 hidden layers of the neural network. Particularly, let di be the number of hidden units in the i-th hidden layer for i = 1, . . . , l − 1. Let d0 and dl be the number of units of input and output layers, respectively. Let W i ∈ Rdi×di−1 be the weight matrix between the (i− 1)-th layer and the i-th layer for any i = 1, . . . l2.
2To simplify notations, we regard the input and output layers as the 0-th and the l-th layers, respectively, and absorb the bias of each layer into W i.
According to Theorem 2, one major condition is to verify the introduced Lyapunov function F defined in (9) satisfies the Kurdyka-Łojasiewicz property. For this purpose, we need an extension of semialgebraic set, called the o-minimal structure (see, for instance Coste (1999), van den Dries (1986), Kurdyka (1998), Bolte et al. (2007b)). The following definition is from Bolte et al. (2007b). Definition 4. [o-minimal structure] An o-minimal structure on (R,+, ·) is a sequence of boolean algebras On of “definable” subsets of Rn, such that for each n ∈ N
(i) if A belongs to On, then A× R and R×A belong to On+1;
(ii) if Π : Rn+1 → Rn is the canonical projection onto Rn, then for any A in On+1, the set Π(A) belongs to On;
(iii) On contains the family of algebraic subsets of Rn, that is, every set of the form {x ∈ Rn : p(x) = 0},
where p : Rn → R is a polynomial function.
(iv) the elements of O1 are exactly finite unions of intervals and points.
Based on the definition of o-minimal structure, we can show the definition of the definable function. Definition 5. [Definable function] Given an o-minimal structure O (over (R,+, ·)), a function f : Rn → R is said to be definable in O if its graph belongs to On+1.
According to van den Dries & Miller (1996); Bolte et al. (2007b), there are some important facts of the o-minimal structure, shown as follows.
(i) The collection of semialgebraic sets is an o-minimal structure. Recall the semialgebraic sets are Bollean combinations of sets of the form
{x ∈ Rn : p(x) = 0, q1(x) < 0, . . . , qm(x) < 0}, where p and qi’s are polynomial functions in Rn.
(ii) There exists an o-minimal structure that contains the sets of the form
{(x, t) ∈ [−1, 1]n × R : f(x) = t} where f is real-analytic around [−1, 1]n.
(iii) There exists an o-minimal structure that contains simultaneously the graph of the exponential function R 3 x 7→ exp(x) and all semialgebraic sets.
(iv) The o-minimal structure is stable under the sum, composition, the inf-convolution and several other classical operations of analysis.
The Kurdyka-Łojasiewicz property for the smooth definable function and non-smooth definable function were established in (Kurdyka, 1998, Theorem 1) and (Bolte et al., 2007b, Theorem 11), respectively. Now we are ready to present the proof of Corollary 1.
Proof. [Proof of Corollary 1] To justify this corollary, we only need to verify the associated Lyapunov function F satisfies Kurdyka-Łojasiewicz inequality. In this case and by (10), F can be rewritten as follows
F (W,Γ,G) = α ( L̂n(W,Γ) + 1
2ν ‖W − Γ‖2
) + Ω(Γ) + Ω∗(g)− 〈W, g〉.
Because ` and σi’s are definable by assumptions, then L̂n(W,Γ) are definable as compositions of definable functions. Moreover, according to Krantz & Parks (2002), ‖W − Γ‖2 and 〈W, g〉 are semi-algebraic and thus definable. Since the group Lasso Ω(Γ) = ∑ g ‖Γ‖2 is the composition of `2 and `1 norms, and the conjugate of group Lasso penalty is the maximum of group `2-norm, i.e. Ω∗(Γ) = maxg ‖Γg‖2, where the `2, `1, and `∞ norms are definable, hence the group Lasso and its conjugate are definable as compositions of definable functions. Therefore, F is definable and hence satisfies Kurdyka-Łojasiewicz inequality by (Kurdyka, 1998, Theorem 1).
The verifications of other cases listed in assumptions can be found in the proof of (Zeng et al., 2019, Proposition 1). This finishes the proof of this corollary.
A.3 PROOF OF THEOREM 2
Our analysis is mainly motivated by a recent paper (Benning et al., 2017), as well as the influential work (Attouch et al., 2013). According to (Attouch et al., 2013, Lemma 2.6), there are mainly four ingredients in the analysis, that is, the sufficient descent property, relative error property, continuity property of the generated sequence and the Kurdyka-Łojasiewicz property of the function. More specifically, we first establish the sufficient descent property of the generated sequence via exploiting the Lyapunov function F (see, (9)) in Lemma A.4 in Section A.4, and then show the relative error property of the sequence in Lemma A.5 in Section A.5. The continuity property is guaranteed by the continuity of L̄(W,Γ) and the relation limk→∞BgkΩ (Γk+1,Γk) = 0 established in Lemma 1(i) in Section A.4. Thus, together with the Kurdyka-Łojasiewicz assumption of F , we establish the global convergence of SLBI following by (Attouch et al., 2013, Lemma 2.6).
Let (W̄ , Γ̄, ḡ) be a critical point of F , then the following holds
∂WF (W̄ , Γ̄, ḡ) = α(∇L̂n(W̄ ) + ν−1(W̄ − Γ̄)) = 0, ∂ΓF (W̄ , Γ̄, ḡ) = αν
−1(Γ̄− W̄ ) + ∂Ω(Γ̄)− ḡ 3 0, (19) ∂gF (W̄ , Γ̄, ḡ) = Γ̄− ∂Ω∗(ḡ) 3 0.
By the final inclusion and the convexity of Ω, it implies ḡ ∈ ∂Ω(Γ̄). Plugging this inclusion into the second inclusion yields αν−1(Γ̄− W̄ ) = 0. Together with the first equality imples
∇L̄(W̄ , Γ̄) = 0, ∇L̂n(W̄ ) = 0. This finishes the proof of this theorem.
A.4 SUFFICIENT DESCENT PROPERTY ALONG LYAPUNOV FUNCTION
Let Pk := (Wk,Γk), and Qk := (Pk, gk−1), k ∈ N. In the following, we present the sufficient descent property of Qk along the Lyapunov function F .
Lemma. Suppose that L̂n is continuously differentiable and ∇L̂n is Lipschitz continuous with a constant Lip > 0. Let {Qk} be a sequence generated by SLBI with a finite initialization. If 0 < α < 2κ(Lip+ν−1) , then
F (Qk+1) ≤ F (Qk)− ρ‖Qk+1 −Qk‖22,
where ρ := 1κ − α(Lip+ν−1) 2 .
Proof. By the optimality condition of (15a) and also the inclusion pk = [0, gk]T ∈ ∂R(Pk), there holds
κ(α∇L̄(Pk) + pk+1 − pk) + Pk+1 − Pk = 0, which implies
−〈α∇L̄(Pk), Pk+1 − Pk〉 = κ−1‖Pk+1 − Pk‖22 +D(Γk+1,Γk) (20) where D(Γk+1,Γk) := 〈gk+1 − gk,Γk+1 − Γk〉. Noting that L̄(P ) = L̂n(W ) + 12ν ‖W − Γ‖ 2 2 and by the Lipschitz continuity of ∇L̂n(W ) with a constant Lip > 0 implies∇L̄ is Lipschitz continuous with a constant Lip+ ν−1. This implies
L̄(Pk+1) ≤ L̄(Pk) + 〈∇L̄(Pk), Pk+1 − Pk〉+ Lip+ ν−1
2 ‖Pk+1 − Pk‖22.
Substituting the above inequality into (20) yields
αL̄(Pk+1) +D(Γk+1,Γk) + ρ‖Pk+1 − Pk‖22 ≤ αL̄(Pk). (21) Adding some terms in both sides of the above inequality and after some reformulations implies
αL̄(Pk+1) +BgkΩ (Γk+1,Γk) (22) ≤ αL̄(Pk) +B gk−1 Ω (Γk,Γk−1)− ρ‖Pk+1 − Pk‖ 2 2 − ( D(Γk+1,Γk) +B gk−1 Ω (Γk,Γk−1)−B gk Ω (Γk+1,Γk) ) = αL̄(Pk) +B gk−1 Ω (Γk,Γk−1)− ρ‖Pk+1 − Pk‖ 2 2 −B gk+1 Ω (Γk,Γk−1)−B gk−1 Ω (Γk,Γk−1),
where the final equality holds for D(Γk+1,Γk)−BgkΩ (Γk+1,Γk) = B gk+1 Ω (Γk,Γk−1). That is,
F (Qk+1) ≤ F (Qk)− ρ‖Pk+1 − Pk‖22 −B gk+1 Ω (Γk,Γk−1)−B gk−1 Ω (Γk,Γk−1) (23)
≤ F (Qk)− ρ‖Pk+1 − Pk‖22, (24)
where the final inequality holds for Bgk+1Ω (Γk,Γk−1) ≥ 0 and B gk−1 Ω (Γk,Γk−1) ≥ 0. Thus, we finish the proof of this lemma.
Based on Lemma A.4, we directly obtain the following lemma.
Lemma 1. Suppose that assumptions of Lemma A.4 hold. Suppose further that Assumption 1 (b)-(d) hold. Then
(i) both α{L̄(Pk)} and {F (Qk)} converge to the same finite value, and limk→∞B gk Ω (Γk+1,Γk) = 0.
(ii) the sequence {(Wk,Γk, gk)} is bounded,
(iii) limk→∞ ‖Pk+1 − Pk‖22 = 0 and limk→∞D(Γk+1,Γk) = 0,
(iv) 1K ∑K k=0 ‖Pk+1 − Pk‖22 → 0 at a rate of O(1/K).
Proof. By (21), L̄(Pk) is monotonically decreasing due to D(Γk+1,Γk) ≥ 0. Similarly, by (24), F (Qk) is also monotonically decreasing. By the lower boundedness assumption of L̂n(W ), both L̄(P ) and F (Q) are lower bounded by their definitions, i.e., (2) and (9), respectively. Therefore, both {L̄(Pk)} and {F (Qk)} converge, and it is obvious that limk→∞ F (Qk) ≥ limk→∞ αL̄(Pk). By (23),
B gk−1 Ω (Γk,Γk−1) ≤ F (Qk)− F (Qk+1), k = 1, . . . .
By the convergence of F (Qk) and the nonegativeness of B gk−1 Ω (Γk,Γk−1), there holds
lim k→∞
B gk−1 Ω (Γk,Γk−1) = 0.
By the definition of F (Qk) = αL̄(Pk) +B gk−1 Ω (Γk,Γk−1) and the above equality, it yields
lim k→∞ F (Qk) = lim k→∞
αL̄(Pk).
Since L̂n(W ) has bounded level sets, then Wk is bounded. By the definition of L̄(W,Γ) and the finiteness of L̄(Wk,Γk), Γk is also bounded due to Wk is bounded. The boundedness of gk is due to gk ∈ ∂Ω(Γk), condition (d), and the boundedness of Γk. By (24), summing up (24) over k = 0, 1, . . . ,K yields
K∑ k=0 ( ρ‖Pk+1 − Pk‖2 +D(Γk+1,Γk) ) < αL̄(P0) <∞. (25)
Letting K →∞ and noting that both ‖Pk+1 − Pk‖2 and D(Γk+1,Γk) are nonnegative, thus
lim k→∞ ‖Pk+1 − Pk‖2 = 0, lim k→∞ D(Γk+1,Γk) = 0.
Again by (25),
1
K K∑ k=0 ( ρ‖Pk+1 − Pk‖2 +D(Γk+1,Γk) ) < K−1αL̄(P0),
which implies 1K ∑K k=0 ‖Pk+1 − Pk‖2 → 0 at a rate of O(1/K).
A.5 RELATIVE ERROR PROPERTY
In this subsection, we provide the bound of subgradient by the discrepancy of two successive iterates. By the definition of F (9),
Hk+1 :=  α∇W L̄(Wk+1,Γk+1)α∇ΓL̄(Wk+1,Γk+1) + gk+1 − gk Γk − Γk+1  ∈ ∂F (Qk+1), k ∈ N. (26) Lemma. Under assumptions of Lemma 1, then
‖Hk+1‖ ≤ ρ1‖Qk+1 −Qk‖, for Hk+1 ∈ ∂F (Qk+1), k ∈ N, where ρ1 := 2κ−1 + 1 + α(Lip+ 2ν−1). Moreover, 1K ∑K k=1 ‖Hk‖2 → 0 at a rate of O(1/K).
Proof. Note that
∇W L̄(Wk+1,Γk+1) = (∇W L̄(Wk+1,Γk+1)−∇W L̄(Wk+1,Γk)) (27) + (∇W L̄(Wk+1,Γk)−∇W L̄(Wk,Γk)) +∇W L̄(Wk,Γk).
By the definition of L̄ (see (2)),
‖∇W L̄(Wk+1,Γk+1)−∇W L̄(Wk+1,Γk)‖ = ν−1‖Γk − Γk+1‖, ‖∇W L̄(Wk+1,Γk)−∇W L̄(Wk,Γk)‖ = ‖(∇L̂n(Wk+1)−∇L̂n(Wk)) + ν−1(Wk+1 −Wk)‖
≤ (Lip+ ν−1)‖Wk+1 −Wk‖,
where the last inequality holds for the Lipschitz continuity of ∇L̂n with a constant Lip > 0, and by (16a),
‖∇W L̄(Wk,Γk)‖ = (κα)−1‖Wk+1 −Wk‖.
Substituting the above (in)equalities into (27) yields ‖∇W L̄(Wk+1,Γk+1)‖ ≤ [ (κα)−1 + Lip+ ν−1 ] · ‖Wk+1 −Wk‖+ ν−1‖Γk+1 − Γk‖
Thus, ‖α∇W L̄(Wk+1,Γk+1)‖ ≤ [ κ−1 + α(Lip+ ν−1) ] · ‖Wk+1 −Wk‖+ αν−1‖Γk+1 − Γk‖. (28)
By (16c), it yields
gk+1 − gk = κ−1(Γk − Γk+1)− α∇ΓL̄(Wk,Γk).
Noting that∇ΓL̄(Wk,Γk) = ν−1(Γk −Wk), and after some simplifications yields
‖α∇ΓL̄(Wk+1,Γk+1) + gk+1 − gk‖ = ‖(κ−1 − αν−1) · (Γk − Γk+1) + αν−1(Wk −Wk+1)‖ ≤ αν−1‖Wk −Wk+1‖+ (κ−1 − αν−1)‖Γk − Γk+1‖,
(29)
where the last inequality holds for the triangle inequality and κ−1 > αν−1 by the assumption.
By (28), (29), and the definition of Hk+1 (26), there holds ‖Hk+1‖ ≤ [ κ−1 + α(Lip+ 2ν−1) ] · ‖Wk+1 −Wk‖+ (κ−1 + 1)‖Γk+1 − Γk‖
≤ [ 2κ−1 + 1 + α(Lip+ 2ν−1) ] · ‖Pk+1 − Pk‖ (30)
≤ [ 2κ−1 + 1 + α(Lip+ 2ν−1) ] · ‖Qk+1 −Qk‖.
By (30) and Lemma 1(iv), 1K ∑K k=1 ‖Hk‖2 → 0 at a rate of O(1/K).
This finishes the proof of this lemma.

heading: B SUPPLEMENTARY EXPERIMENTS
text: B.1 ABLATION STUDY OF VGG16 AND RESNET56 ON CIFAR10
To further study the influence of hyperparameters, we record performance of Wt for each epoch t with different combinations of hyperparameters. The experiments is conducted 5 times each, we show the mean in the table, the standard error can be found in the corresponding figure. We perform experiments on Cifar10 and two commonly used network VGG16 and ResNet56.
On κ , we keep ν = 100 and try κ = 1, 2, 5, 10, the validation curves of models Wt are shown in Fig. 6 and Table 2 summarizes the mean accuracies. Table 3 summarizes best validation accuracies achieved at some epochs, together with their sparsity rates. These results show that larger kappa leads to slightly lower validation accuracies, where the numerical results are shown in Table 2 . We can find that κ = 1 achieves the best test accuracy.
On ν , we keep κ = 1 and try ν = 10, 20, 50, 100, 200, 500, 1000, 2000 the validation curve and mean accuracies are show in Fig. 6 and Table 4. Table 5 summarizes best validation accuracies achieved at some epochs, together with their sparsity rates. By carefully tuning ν we can achieve similar or even better results compared to SGD. Different from κ, ν has less effect on the generalization performance. By tuning it carefully, we can even get a sparse model with slightly better performance than SGD trained model.

heading: C COMPUTATIONAL COST OF SPLITLBI
text: We further compare the computational cost of different optimizers: SGD (Mom), SplitLBI (Mom) and Adam (Naive). We test each optimizer on one GPU, and all the experiments are done on one GTX2080. For computational cost, we judge them from two aspects : GPU memory usage and time needed for one batch. The batch size here is 64, experiment is performed on VGG-16 as shown in Table 6.

heading: D FINE-TUNING OF SPARSE SUBNETWORKS
text: We design the experiment on MNIST, inspired by Frankle & Carbin (2019). Here, we explore the subnet obtained by ΓT after T = 100 epochs of training. As in Frankle et al. (2019), we adopt the “rewind” trick: re-loading the subnet mask of Γ100 at different epochs, followed by fine-tuning. In particular, along the training paths, we reload the subnet models at Epoch 0, Epoch 30, 60, 90, and 100, and further fine-tune these models by SplitLBI (Mom-Wd). All the models use the same initialization and hence the subnet model at Epoch 0 gives the retraining with the same random initialization as proposed to find winning tickets of lottery in Frankle & Carbin (2019). We will denote the rewinded fine-tuned model at epoch 0 as (Lottery), and those at epoch 30, 60, 90, and 100, as F-epoch30, F-epoch60, F-epoch90, and F-epoch100, respectively. Three networks are studied here – LeNet-3, Conv-2, and Conv-4. LeNet-3 removes one convolutional layer of LeNet-5; and it is thus less over-parameterized than the other two networks. Conv-2 and Conv-4, as the scaled-down variants of VGG family as done in Frankle & Carbin (2019), have two and four fully-connected layers, respectively, followed by max-pooling after every two convolutional layer.
The whole sparsity for Lenet-3 is 0.055, Conv-2 is 0.0185, and Conv-4 is 0.1378. Detailed sparsity for every layer of the model is shown in Table 7, 8, 9. We find that fc-layers are sparser than conv-layers.
We compare SplitLBI variants to the SGD (Mom-Wd) and SGD (Lottery) (Frankle & Carbin, 2019) in the same structural sparsity and the results are shown in Fig. 7. In this exploratory experiment, one can see that for overparameterized networks – Conv-2 and Conv-4, fine-tuned rewinding subnets – F-epoch30, F-epoch60, F-epoch90, and F-epoch100, can produce better results than the full models; while for the less over-parameterized model LeNet-3, fine-tuned subnets may achieve less yet still comparable performance to the dense models and remarkably better than the retrained sparse subnets from beginning (i.e. SplitLBI/SGD (Lottery)). These phenomena suggest that the subnet architecture disclosed by structural sparsity parameter ΓT is valuable, for fine-tuning sparse models with comparable or even better performance than the dense models of WT .

heading: E RETRAINING OF SPARSE SUBNETS FOUND BY SPLITLBI (LOTTERY)
text: Here we provide more details on the experiments in Fig. 5. Table 10 gives the details on hyperparameter setting. Moreover, Figure 8 provides the sparsity variations during SplitLBI training in Fig. 5.

Review summaries:
This paper studies the problem of finding important structural sparsity of over-parameterized neural networks. Compared to the traditional pruning and distillation based compression methods, the author propose a novel optimization based algorithm called SplitLBI, which can efficiently and effectively find the important subnetwork architecture. In addition, the authors provide the convergence guarantee of the proposed method for solving nonconvex optimization problems under certain conditions. Thorough experiments and ablation studies validate the advantage of the proposed method.