Title: SPLIT LBI FOR DEEP LEARNING: STRUCTURAL SPAR-SITY VIA DIFFERENTIAL INCLUSION PATHS

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.

1. INTRODUCTION:

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

2. METHODOLOGY:

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 ) = 1 n n i=1 (y i , Φ W (x i )). For example, a neural network of l-layer is defined as Φ W (x) = σ l (W l σ l-1 (W l-1 • • • σ 1 (W 1 x)), where W = {W i } l

i=1 , σ i is the nonlinear activation function of the i-th layer.

Differential Inclusion of Inverse Scale Space. Consider the following dynamics,

Ẇt κ = -∇ W L (W t , Γ t ) (1a) Vt = -∇ Γ L (W t , Γ t ) (1b) V t ∈ ∂ Ω(Γ t ) (1c)

where V is a sub-gradient of Ω(Γ) := Ω λ (Γ) + 1 2κ Γ 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

EQUATION

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 W t 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 V t , driving the flow in sparse primal space of Γ t . For a large enough ν, it reduces to the gradient descent method for W t . 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, V t 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 V t 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 W t follows the gradient descent with a standard 2 -regularization. Therefore, W t 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

EQUATION

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):

EQUATION

EQUATION

EQUATION

where V 0 = Γ 0 = 0, W 0 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

EQUATION

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 (c in , c out , size) denote the convolutional filters where size denotes the kernel size and c in and c out denote the numbers of input channels and output channels, respectively. When we regard each group as each convolutional filter, g = c out ; otherwise for weight sparsity, g can be every element in the filter that reduces to the Lasso.

(2) For a fully connected layer, Γ = Γ(c in , c out ) where c in and c out 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.

3. GLOBAL CONVERGENCE OF SPLITLBI FOR NEURAL NETWORKS:

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,

EQUATION

where

EQUATION

p k ∈ ∂Ψ(P k ), and B q Ψ is the Bregman divergence associated with convex function Ψ, defined by

EQUATION

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 ) = 1 n n i=1 (y i , Φ W (x i )) 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 ⊂ R n , there exists a constant C > 0 such that for all Γ ∈ S and all g ∈ ∂Ω(Γ), there holds g ≤ C; and (e) the Lyapunov function

EQUATION

is a Kurdyka-Łojasiewicz function on any bounded set, where B g Ω (Γ, Γ) := Ω(Γ) -Ω( Γ) -g, Γ -Γ , Γ ∈ ∂Ω * (g), and Ω * is the conjugate of Ω defined as 9) and the definition of conjugate, the Lyapunov function F can be rewritten as follows,

Ω * (g) := sup U ∈R n { U, g -Ω(U )}.

EQUATION

Now we are ready to present the main theorem.

Theorem 1. [Global Convergence of SplitLBI] Suppose that Assumption 1 holds. Let (W k , Γ k ) be the sequence generated by SplitLBI (Eq. (4a-4c)) with a finite initialization. If

0 < α k = α < 2 κ(Lip + ν -1 )

, then (W k , Γ 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 {W k , Γ k , g k } 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 (t 2 ), exponential loss (e t ), logistic loss log(1 + e -t ), and cross-entropy loss; (b) σ i is any smooth definable activation, such as linear activation (t), sigmoid (1 1+e -t ), hyperbolic tangent ( e t -e -t e t +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 {W k } 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.

4. EXPERIMENTS WITH STOCHASTIC SPLITLBI:

We begin with some stochastic variants of SplitLBI and implementations, followed by four groups of experiments demonstrating the utilities of weight parameter W t 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 ,

EQUATION

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,

EQUATION

)

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,

EQUATION

)

where β is set as 1e -4 . 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 .

Implementation:



4.1. IMAGE CLASSIFICATION:

In SplitLBI, the weight parameter W t 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, 2) SGD with l 1 penalty (Lasso). The l 1 norm is applied to penalize the weights of SGD by encouraging the sparsity of learned model, with the regularization parameter of the l 1 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 l 2 weight decay with the coefficient weight 1e -4 . (5) SGD with Nesterov (Nesterov): the SGD uses nesterov momentum 0.9. 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 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 i j . Then we compute the gradient of the sum of the feature map computed from each filter W i j 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 Abbasi-Asl & 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.

4.3. GLOBAL CONVERGENCE AND STRUCTURAL SPARSITY OF SPLITLBI:

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 W t 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 W t .

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 W t and Γ t such that the growth of W t 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 W t through their coupling.

4.4. SPLITLBI WITH EARLY STOPPING AND RETRAIN FINDS EFFECTIVE SUB-NETWORKS:

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 

5. CONCLUSION:

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.

APPENDIX TO SplitLBI for deep learning: structural sparsity via differential inclusion paths:



A PROOF OF THEOREM 1:

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,

EQUATION

EQUATION

where

p k = [0, g k ] T ∈ ∂R(P k ) and g k ∈ ∂Ω(Γ k ).

Thus SplitLBI is equivalent to the following iterations,

EQUATION

EQUATION

EQUATION

Exploiting the equivalent reformulation (16a-16c), one can establish the global convergence of (W k , Γ k , g k ) 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 (W k , Γ k , g k ) be the sequence generated by SplitLBI (Eq. (16a-16c)) with a finite initialization. If

0 < α k = α < 2 κ(Lip + ν -1 )

, then (W k , Γ k , g k ) converges to a critical point of F . Moreover, {(W k , Γ 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 ∈ R p : h(x) < +∞}, the Fréchet subdifferential of h at x, written ∂h(x), is the set of vectors v ∈ R p which satisfy lim inf

y =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

EQUATION

A necessary (but not sufficient) condition for x ∈ R p 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 R p , written dist(x, S), is defined by dist(x, S) = inf{ xs : s ∈ S}, where • represents the Euclidean norm.

Let h : R p → R ∪ {+∞} be an extended-real-valued function (respectively, h : R p ⇒ R q be a point-to-set mapping), its graph is defined by

Graph(h) := {(x, y) ∈ R p × R : y = h(x)}, (resp. Graph(h) := {(x, y) ∈ R p × R q : y ∈ h(x)}),

and its domain by dom(h) := {x ∈ R p : h(x) < +∞} (resp. dom(h) := {x ∈ R p : h(x) = ∅}). When h is a proper function, i.e., when dom(h) = ∅, the set of its global minimizers (possibly empty) is denoted by arg min h := {x ∈ R p : h(x) = inf h}.

The KL property (Łojasiewicz, 1963; 1993; Kurdyka, 1998 

) := {v ∈ R n |∂h(v) = ∅}, if there exists a constant η ∈ (0, ∞),

a neighborhood N of ū and a function φ : [0, η) → R + , which is a concave function that is continuous at 0 and satisfies φ(0) = 0, φ ∈ C 1 ((0, η)), i.e., φ is continuous differentiable on (0, η), and φ (s) > 0 for all s ∈ (0, η), such that for all u ∈ N ∩ {u ∈ R n |h(ū) < h(u) < h(ū) + η}, the following inequality holds

EQUATION

)

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 R p 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 R p is analytic by checking the analyticity of g(t) := h(x + ty) for any x, y ∈ R p . Definition 3.

[Semialgebraic] (a) A set D ⊂ R p is called semialgebraic (Bochnak et al., 1998) if it can be represented as

D = s i=1 t j=1 {x ∈ R p : P ij (x) = 0, Q ij (x) > 0} , where P ij , Q ij are real polynomial functions for 1 ≤ i ≤ s, 1 ≤ j ≤ t. (b) A function h : R p → R ∪ {+∞} (resp. a point-to-set mapping h : R p ⇒ R q ) 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 d i be the number of hidden units in the i-th hidden layer for i = 1, . . . , l -1. Let d 0 and d l be the number of units of input and output layers, respectively. Let W i ∈ R di×di-1 be the weight matrix between the (i -1)-th layer and the i-th layer for any i = 1, . . . l2 .

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 (i) The collection of semialgebraic sets is an o-minimal structure. Recall the semialgebraic sets are Bollean combinations of sets of the form

{x ∈ R n : p(x) = 0, q 1 (x) < 0, . . . , q m (x) < 0},

where p and q i 's are polynomial functions in R n . (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 x → 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 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 lim k→∞ B g k Ω (Γ 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 , Γ, ḡ) be a critical point of F , then the following holds

EQUATION

∂ g F ( W , Γ, ḡ) = Γ -∂Ω * (ḡ) 0.

By the final inclusion and the convexity of Ω, it implies ḡ ∈ ∂Ω( Γ). 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

EQUATION

Q k := (P k , g k-1 ), k ∈ N.

In the following, we present the sufficient descent property of Q k 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 {Q k } be a sequence generated by SLBI with a finite initialization. If

0 < α < 2 κ(Lip+ν -1 ) , then F (Q k+1 ) ≤ F (Q k ) -ρ Q k+1 -Q k 2 2 , where ρ := 1 κ -α(Lip+ν -1 ) 2 .

Proof. By the optimality condition of (15a) and also the inclusion p k = [0, g k ] T ∈ ∂R(P k ), there holds κ(α∇ L(P k ) + p k+1 -p k ) + P k+1 -P k = 0, which implies

EQUATION

where

D(Γ k+1 , Γ k ) := g k+1 -g k , Γ k+1 -Γ k .

Noting that

L(P ) = L n (W ) + 1 2ν 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(P k+1 ) ≤ L(P k ) + ∇ L(P k ), P k+1 -P k + Lip + ν -1 2 P k+1 -P k 2 2 .

Substituting the above inequality into (20) yields

EQUATION

Adding some terms in both sides of the above inequality and after some reformulations implies

α L(P k+1 ) + B g k Ω (Γ k+1 , Γ k ) (22) ≤ α L(P k ) + B g k-1 Ω (Γ k , Γ k-1 ) -ρ P k+1 -P k 2 2 -D(Γ k+1 , Γ k ) + B g k-1 Ω (Γ k , Γ k-1 ) -B g k Ω (Γ k+1 , Γ k ) = α L(P k ) + B g k-1 Ω (Γ k , Γ k-1 ) -ρ P k+1 -P k 2 2 -B g k+1 Ω (Γ k , Γ k-1 ) -B g k-1 Ω (Γ k , Γ k-1 ),

where the final equality holds for

EQUATION

where the final inequality holds for B

g k+1 Ω (Γ k , Γ k-1 ) ≥ 0 and B g k-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(P k )} and {F (Q k )} converge to the same finite value, and

lim k→∞ B g k Ω (Γ k+1 , Γ k ) = 0. (ii) the sequence {(W k , Γ k , g k )} is bounded, (iii) lim k→∞ P k+1 -P k 2 2 = 0 and lim k→∞ D(Γ k+1 , Γ k ) = 0, (iv) 1 K K k=0 P k+1 -P k 2 2 → 0 at a rate of O(1/K).

Proof. By (21), L(P k ) is monotonically decreasing due to D(Γ k+1 , Γ k ) ≥ 0. Similarly, by (24), F (Q k ) 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(P k )} and {F (Q k )} converge, and it is obvious that lim k→∞ F (Q k ) ≥ lim k→∞ α L(P k ). By (23),

B g k-1 Ω (Γ k , Γ k-1 ) ≤ F (Q k ) -F (Q k+1 ), k = 1, . . . .

By the convergence of F (Q k ) and the nonegativeness of B g k-1 Ω (Γ k , Γ k-1 ), there holds

lim k→∞ B g k-1 Ω (Γ k , Γ k-1 ) = 0. By the definition of F (Q k ) = α L(P k ) + B g k-1 Ω (Γ k , Γ k-1

) and the above equality, it yields

lim k→∞ F (Q k ) = lim k→∞ α L(P k ).

Since L n (W ) has bounded level sets, then W k is bounded. By the definition of L(W, Γ) and the finiteness of L(W k , Γ k ), Γ k is also bounded due to W k is bounded. The boundedness of g k is due to g k ∈ ∂Ω(Γ k ), condition (d), and the boundedness of Γ k . By (24), summing up (24) over k = 0, 1, . . . , K yields

EQUATION

)

Letting K → ∞ and noting that both P k+1 -P k 2 and D(Γ k+1 , Γ k ) are nonnegative, thus

lim k→∞ P k+1 -P k 2 = 0, lim k→∞ D(Γ k+1 , Γ k ) = 0.

Again by (25),

1 K K k=0 ρ P k+1 -P k 2 + D(Γ k+1 , Γ k ) < K -1 α L(P 0 ), which implies 1 K K k=0 P k+1 -P k 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),

EQUATION

Lemma. Under assumptions of Lemma 1, then

H k+1 ≤ ρ 1 Q k+1 -Q k , for H k+1 ∈ ∂F (Q k+1 ), k ∈ N,

where

ρ 1 := 2κ -1 + 1 + α(Lip + 2ν -1 ). Moreover, 1 K K k=1 H k 2 → 0 at a rate of O(1/K).

Proof. Note that

EQUATION

+ (∇ W L(W k+1 , Γ k ) -∇ W L(W k , Γ k )) + ∇ W L(W k , Γ k ).

By the definition of L (see ( 2)),

∇ W L(W k+1 , Γ k+1 ) -∇ W L(W k+1 , Γ k ) = ν -1 Γ k -Γ k+1 , ∇ W L(W k+1 , Γ k ) -∇ W L(W k , Γ k ) = (∇ L n (W k+1 ) -∇ L n (W k )) + ν -1 (W k+1 -W k ) ≤ (Lip + ν -1 ) W k+1 -W k ,

where the last inequality holds for the Lipschitz continuity of ∇ L n with a constant Lip > 0, and by (16a),

∇ W L(W k , Γ k ) = (κα) -1 W k+1 -W k .

Substituting the above (in)equalities into (27) yields

EQUATION

By (16c), it yields

g k+1 -g k = κ -1 (Γ k -Γ k+1 ) -α∇ Γ L(W k , Γ k ). Noting that ∇ Γ L(W k , Γ k ) = ν -1 (Γ k -W k )

, and after some simplifications yields

EQUATION

where the last inequality holds for the triangle inequality and κ -1 > αν -1 by the assumption.

By (28), (29), and the definition of H k+1 (26), there holds

H k+1 ≤ κ -1 + α(Lip + 2ν -1 ) • W k+1 -W k + (κ -1 + 1) Γ k+1 -Γ k ≤ 2κ -1 + 1 + α(Lip + 2ν -1 ) • P k+1 -P k (30) ≤ 2κ -1 + 1 + α(Lip + 2ν -1 ) • Q k+1 -Q k .

By (30) and Lemma 1(iv),

1 K K k=1 H k 2 → 0 at a rate of O(1/K).

This finishes the proof of this lemma. On κ , we keep ν = 100 and try κ = 1, 2, 5, 10, the validation curves of models W t 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.

C COMPUTATIONAL COST OF SPLITLBI:

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 . 

D FINE-TUNING OF SPARSE SUBNETWORKS:

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 W T .

E RETRAINING OF SPARSE SUBNETS FOUND BY SPLITLBI (LOTTERY):

Here we provide more details on the experiments in Fig. 5 . Table 10 gives the details on hyperparameter setting. Moreover, Figure 8 Group Lasso SLBI 0.1 100 1 0.1 0.9 Yes ResNet-56 Group Lasso SLBI 0.1 100 1 0.05 0.9 Yes VGG-16(Lasso) Lasso SLBI 0.1 500 1 0.05 0.9 Yes ResNet-50(Lasso) Lasso SLBI 0.1 200 1 0.03 0.9 Yes Table 10 : Hyperparameter setting for the experiments in Figure 5 . The vertical blue line shows the epochs that we choose to early stop. We choose the log-scale epochs for achieve larger range of sparsity.

:

https://anonymous.4open.science/repository/d22bbbc8-50d5-4e60-b2e8-4ded4e93db63/ Split_LBI_code

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

and 1 norms, and the conjugate of group Lasso penalty is the maximum of group 2 -norm, i.e. Ω * (Γ) = max g Γ 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.

Figures and Tables:

Figure 1: Visualization of solution path and filter patterns in the third convolutional layer (i.e., conv.c5) of LetNet-5, trained on MNIST. The left figure shows the magnitude changes for each filter of the models trained by SplitLBI and SGD, where x-axis and y-axis indicate the training epochs, and filter magnitudes ( 2-norm), respectively. The SplitLBI path of filters selected in the support of Γ are drawn in blue color, while the red color curves represent the filters that are not important and outside the support of Γ. We visualize the corresponding learned filters by Erhan et al. (2009) at 20 (blue), 40 (green), and 80 (black) epochs, which are shown in the right figure with the corresponding color bounding boxes, i.e., blue, green, and black, respectively. It shows that our SplitLBI enjoys a sparse selection of filters without sacrificing accuracy (see Table1).

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

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, asHe et al. (2016). The weights of all models are initialized asHe et al. (2015).

Figure 2: Visualization of the first convolutional layer filters of ResNet-18 trained on ImageNet-2012. Given the input image and initial weights visualized in the middle, filter response gradients at 20 (purple), 40 (green), and 60 (black) epochs are visualized by Springenberg et al. (2014). The "SLBI-10" ("SLBI-1") in the right figure refers to SplitLBI with κ = 10 and κ = 1, respectively.

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

Figure 3: Loss and training accuracy by different κ and ν validates the global convergence of SplitLBI. The X-axis and Y-axis indicate the training epochs, and sparsity/accuracy. The results are repeated for 5 rounds, by keeping the exactly same initialization for each model. In each round, we use the same initialization for every hyperparameter. For all models, we train for 160 epochs with initial learning rate (lr) of 0.1 and decrease by 0.1 at epoch 80 and 120.

Figure4: Sparsity and validation accuracy by different κ and ν show that moderate sparse models may achieve comparable test accuracies to dense models without fine-tuning. Sparsity is obtained as the percentage of nonzeros in Γ t and sparse model at epoch t is obtained by projection of W t onto the support set of Γ t , i.e. pruning the weights corresponding to zeros in Γ t . The best accuracies achieved are recorded in Tab. 2 and 4 of Appendix for different κ and ν, respectively. X-axis and Y-axis indicate the training epochs, and sparsity/accuracy. The results are repeated for 5 times. Shaded area indicates the variance; and in each round, we keep the exactly same initialization for each model. In each round, we use the same initialization for every hyperparameter. For all the model, we train for 160 epochs with initial learning rate (lr) of 0.1 and decrease by 0.1 at epoch 80 and 120.

Figure 5: SplitLBI with early stopping finds sparse subnets whose test accuracies (stars) after retrain are comparable or even better than the baselines (Network Slimming (reproduced by the released codes from Liu et al. (2019) ) , Soft-Filter Pruning(Tab. 10), Scratch-B(Tab. 10), Scratch-E(Tab. 10), and "Rethinking-Lottery" (Tab. 9a)) as reported in Liu et al. (2019). Sparse filters of VGG-16 and ResNet-56 are show in (a) and (b), while sparse weights of VGG-16 and ResNet-50 are shown in (c) and (d).

Figure6: Validation curves of dense models W t for different κ and ν. For SLBI we find that the model accuracy is robust to the hyperparameters both in terms of convergence rate and generalization ability. Here validation accuracy means the accuracy on test set of Cifar10. The first one is the result for VGG16 ablation study on κ, the second one is the result for ResNet56 ablation study on κ, the third one is the result for VGG16 ablation study on ν and the forth one is the result for ResNet56 ablation study on ν.

Figure 7: Fine-tuning of sparse subnets learned by SplitLBI may achieve comparable or better performance than dense models. F-epochk indicates the fine-tuned model comes from the Epoch k. SplitLBI (Lottery) and SGD (Lottery) use the same sparsity rate for each layer and the same initialization for retrain.

Figure8: Sparsity changing during training process of SplitLBI (Lottery) for VGG and ResNets (corresponding to Fig.5). We calculate the sparsity in every epoch and repeat five times. The black curve represents the mean of the sparsity and shaded area shows the standard deviation of sparsity. The vertical blue line shows the epochs that we choose to early stop. We choose the log-scale epochs for achieve larger range of sparsity.

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;

Table content: <table/>



Table content: <table/>

Top-1/Top-5 accuracy(%) on ImageNet-2012 and test accuracy on MNIST/Cifar-10. : results from the official pytorch website. We use the official pytorch codes to run the competitors. All models are trained by 100 epochs. In this table, we run the experiment by ourselves except for SGD Mom-Wd on ImageNet which is reported in https://pytorch.org/docs/stable/torchvision/models.html.

Table content: <table><tr><td/><td>Dataset</td><td>MNIST</td><td>Cifar-10</td><td colspan="2">ImageNet-2012</td></tr><tr><td>Models</td><td>Variants</td><td>LeNet</td><td>ResNet-20</td><td>AlexNet</td><td>ResNet-18</td></tr><tr><td/><td>Naive</td><td>98.87</td><td>86.46</td><td>-/-</td><td>60.76/79.18</td></tr><tr><td/><td>l1</td><td>98.52</td><td>67.60</td><td>-/-</td><td>-/-</td></tr><tr><td>SGD</td><td>Mom</td><td>99.16</td><td>89.44</td><td>55.14/78.09</td><td>66.98/86.97</td></tr><tr><td/><td>Mom-Wd</td><td>99.23</td><td>90.31</td><td>56.55/79.09</td><td>69.76/89.18</td></tr><tr><td/><td>Nesterov</td><td>99.23</td><td>90.18</td><td>-/-</td><td>70.19/89.30</td></tr><tr><td/><td>Naive</td><td>99.19</td><td>89.14</td><td>-/-</td><td>59.66/83.28</td></tr><tr><td/><td>Adabound</td><td>99.15</td><td>87.89</td><td>-/-</td><td>-/-</td></tr><tr><td>Adam</td><td>Adagrad</td><td>99.02</td><td>88.17</td><td>-/-</td><td>-/-</td></tr><tr><td/><td>Amsgrad</td><td>99.14</td><td>88.68</td><td>-/-</td><td>-/-</td></tr><tr><td/><td>Radam</td><td>99.08</td><td>88.44</td><td>-/-</td><td>-/-</td></tr><tr><td/><td>Naive</td><td>99.02</td><td>89.26</td><td>55.06/77.69</td><td>65.26/86.57</td></tr><tr><td>SplitLBI</td><td>Mom</td><td>99.19</td><td>89.72</td><td>56.23/78.48</td><td>68.55/87.85</td></tr><tr><td/><td>Mom-Wd</td><td>99.20</td><td>89.95</td><td>57.09/79.86</td><td>70.55/89.56</td></tr></table>

; 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 ū ∈ dom(∂h

Table content: <table/>

den Dries (1986),Kurdyka (1998),Bolte et al. (2007b)). The following definition is fromBolte et al. (2007b). Definition 4. [o-minimal structure] An o-minimal structure on (R, +, •) is a sequence of boolean algebras O n of "definable" subsets of R n , such that for each n ∈ N(i) if A belongs to O n , then A × R and R × A belong to O n+1 ;(ii) if Π : R n+1 → R n is the canonical projection onto R n , then for any A in O n+1 , the set Π(A) belongs to O n ;(iii) O n contains the family of algebraic subsets of R n , that is, every set of the form{x ∈ R n : p(x) = 0},where p : R n → R is a polynomial function.(iv) the elements of O 1 are exactly finite unions of intervals and points.

Table content: <table/>

6: Validation curves of dense models W t for different κ and ν. For SLBI we find that the model accuracy is robust to the hyperparameters both in terms of convergence rate and generalization ability. Here validation accuracy means the accuracy on test set of Cifar10. The first one is the result for VGG16 ablation study on κ, the second one is the result for ResNet56 ablation study on κ, the third one is the result for VGG16 ablation study on ν and the forth one is the result for ResNet56 ablation study on ν. This table shows results for different κ, the results are all the best test accuracy. Here we test two widely-used models: VGG16 and ResNet56 on Cifar10. For results in this table, we keep ν = 100. Full means that we use the trained model weights directly, Sparse means the model weights are combined with mask generated by Γ support. Sparse result has no finetuning process, the result is comparable to its Full counterpart. For this experiment, we propose that κ = 1 is a good choice. For all the model, we train for 160 epochs with initial learning rate (lr) of 0. 1 and decrease by 0.1 at epoch 80 and 120.To further study the influence of hyperparameters, we record performance of W t 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.

Table content: <table><tr><td>Type</td><td>Model</td><td colspan="3">κ = 1 κ = 2 κ = 5 κ = 10 SGD</td></tr><tr><td>Full</td><td colspan="2">Vgg16 ResNet56 92.71 92.18 91.50 93.46 93.27 92.77</td><td>92.03 90.92</td><td>93.57 93.08</td></tr><tr><td>Sparse</td><td colspan="2">Vgg16 ResNet56 92.37 91.85 89.48 93.31 93.00 92.36</td><td>76.25 87.02</td><td>--</td></tr><tr><td colspan="3">B SUPPLEMENTARY EXPERIMENTS</td><td/></tr><tr><td colspan="4">B.1 ABLATION STUDY OF VGG16 AND RESNET56 ON CIFAR10</td></tr></table>

Sparsity rate and validation accuracy for different κ at different epochs. Here we pick the test accuracy for specific epoch. In this experiment, we keep ν = 100. We pick epoch 20, 40, 80 and 160 to show the growth of sparsity and sparse model accuracy. Here Sparsity is defined in Sec. 4, and Acc means the test accuracy for sparse model. A sparse model is a model at designated epoch t combined with the mask as the support of Γ t .

Table content: <table><tr><td/><td>Model</td><td/><td>Ep20</td><td/><td>Ep40</td><td/><td/><td>Ep80</td><td colspan="2">Ep160</td></tr><tr><td/><td/><td>Term</td><td>Sparsity</td><td>Acc</td><td>Spasity</td><td>Acc</td><td>Spasity</td><td>Acc</td><td>Spasity</td><td>Acc</td></tr><tr><td/><td/><td>κ = 1</td><td colspan="8">96.62 71.51 96.62 76.92 96.63 77.48 96.63 93.31</td></tr><tr><td/><td>Vgg16</td><td>κ = 2</td><td colspan="8">51.86 72.98 71.99 73.64 75.69 74.54 75.72 93.00</td></tr><tr><td/><td/><td>κ = 5</td><td colspan="8">8.19 10.00 17.64 34.25 29.76 69.92 30.03 92.36</td></tr><tr><td/><td colspan="2">κ = 10</td><td colspan="8">0.85 10.00 6.62 10.00 12.95 38.38 13.26 76.25</td></tr><tr><td/><td/><td>Term</td><td>Sparsity</td><td>Acc</td><td>Spasity</td><td>Acc</td><td>Spasity</td><td>Acc</td><td>Spasity</td><td>Acc</td></tr><tr><td/><td/><td>κ = 1</td><td colspan="8">96.79 73.50 96.87 75.27 96.69 77.47 99.68 92.37</td></tr><tr><td/><td>ResNet56</td><td>κ = 2</td><td colspan="8">76.21 72.85 81.41 74.72 84.17 75.64 84.30 91.85</td></tr><tr><td/><td/><td>κ = 5</td><td colspan="8">36.58 60.43 53.07 76.00 57.48 75.67 57.74 89.48</td></tr><tr><td/><td colspan="2">κ = 10</td><td colspan="8">3.12 10.20 29.43 53.36 41.18 74.56 41.14 87.02</td></tr><tr><td>Type</td><td>Model</td><td>ν = 10</td><td>ν = 20</td><td>ν = 50</td><td>ν = 100</td><td/><td>ν = 200</td><td>ν = 500</td><td>ν = 1000</td><td>ν = 2000</td><td>SGD</td></tr><tr><td>Full</td><td>Vgg16 ResNet56</td><td colspan="3">93.66 93.59 93.57 93.12 92.68 92.78</td><td>93.39 92.45</td><td/><td>93.38 92.95</td><td>93.35 93.11</td><td>93.43 93.16</td><td>93.46 93.31</td><td>93.57 93.08</td></tr><tr><td>Sparse</td><td>Vgg16 ResNet56</td><td colspan="3">93.39 93.42 93.39 92.81 92.19 92.40</td><td>93.23 92.10</td><td/><td>93.21 92.68</td><td>93.01 92.81</td><td>92.68 92.84</td><td>10 88.96</td><td>--</td></tr></table>

Results for different ν, the results are all the best test accuracy. Here we test two widely-used model : VGG16 and ResNet56 on Cifar10. For results in this table, we keep κ = 1. Full means that we use the trained model weights directly, Sparse means the model weights are combined with mask generated by Γ support. Sparse result has no finetuning process, the result is comparable to its Full counterpart. For all the model, we train for 160 epochs with initial learning rate (lr) of 0.1 and decrease by 0.1 at epoch 80 and 120.

Table content: <table><tr><td/><td>Model</td><td>Ep20</td><td/><td>Ep40</td><td/><td>Ep80</td><td/><td>Ep160</td></tr><tr><td/><td>Term</td><td>Sparsity</td><td>Acc</td><td>Spasity</td><td>Acc</td><td>Spasity</td><td>Acc</td><td>Spasity</td><td>Acc</td></tr><tr><td/><td>ν = 10</td><td colspan="8">96.64 71.07 96.64 77.70 96.65 79.46 96.65 93.34</td></tr><tr><td/><td>ν = 20</td><td colspan="8">96.64 69.11 96.64 77.63 96.65 77.08 96.65 93.42</td></tr><tr><td/><td>ν = 50</td><td colspan="8">96.64 74.91 96.65 74.21 96.65 79.15 96.65 93.38</td></tr><tr><td>Vgg16</td><td>ν = 100</td><td colspan="8">96.64 74.82 96.64 73.22 96.64 78.09 96.64 93.23</td></tr><tr><td/><td>ν = 200</td><td colspan="8">91.69 73.67 94.06 74.67 94.15 75.20 94.15 93.21</td></tr><tr><td/><td>ν = 500</td><td colspan="8">18.20 10.00 59.94 67.88 82.03 78.69 82.32 93.01</td></tr><tr><td/><td>ν = 1000</td><td colspan="8">6.43 10.00 17.88 10.00 49.75 61.31 51.21 92.68</td></tr><tr><td/><td>ν = 2000</td><td colspan="8">0.22 10.00 6.89 10.00 18.15 10.00 19.00 10.00</td></tr><tr><td/><td>Term</td><td>Sparsity</td><td>Acc</td><td>Spasity</td><td>Acc</td><td>Spasity</td><td>Acc</td><td>Spasity</td><td>Acc</td></tr><tr><td/><td>ν = 10</td><td colspan="8">99.97 73.37 99.95 71.64 99.74 76.46 99.74 92.81</td></tr><tr><td/><td>ν = 20</td><td colspan="8">99.97 72.58 99.84 74.16 99.69 72.37 99.72 92.19</td></tr><tr><td/><td>ν = 50</td><td colspan="8">99.96 70.72 99.89 73.96 99.79 74.93 99.77 92.40</td></tr><tr><td>ResNet56</td><td>ν = 100</td><td colspan="8">96.31 73.63 96.63 75.79 96.55 72.94 96.57 92.10</td></tr><tr><td/><td>ν = 200</td><td colspan="8">91.98 75.30 94.38 72.13 94.87 73.75 94.88 92.68</td></tr><tr><td/><td>ν = 500</td><td colspan="8">74.44 65.58 90.00 74.12 92.96 71.91 92.99 92.81</td></tr><tr><td/><td>ν = 1000</td><td colspan="8">24.32 10.85 75.68 70.23 88.56 79.67 88.80 92.48</td></tr><tr><td/><td>ν = 2000</td><td colspan="8">0.65 10.00 26.66 13.30 74.98 70.38 75.92 88.95</td></tr></table>

Sparsity rate and validation accuracy for different ν at different epochs. Here we pick the test accuracy for specific epoch. In this experiment, we keep κ = 1. We pick epoch 20, 40, 80 and 160 to show the growth of sparsity and sparse model accuracy. Here Sparsity is defined in Sec. 4 as the percentage of nonzero parameters, and Acc means the test accuracy for sparse model. A sparse model is a model at designated epoch t combined with mask as the support of Γ t .

Table content: <table><tr><td>optimizer</td><td>SGD</td><td>SLBI</td><td>Adam</td></tr><tr><td>Mean Batch Time</td><td>0.0197</td><td>0.0221</td><td>0.0210</td></tr><tr><td>GPU Memory</td><td colspan="3">1161MB 1459MB 1267MB</td></tr></table>

Computational and Memory Costs. ( And GPU memory means the whole memory model weights and the activation cache. )

Table content: <table><tr><td>Layer</td><td>FC1</td><td>FC2</td><td>FC3</td></tr><tr><td>Sparsity</td><td>0.049</td><td colspan="2">0.087 0.398</td></tr><tr><td colspan="4">Number of Weights 235200 30000 1000</td></tr></table>

This table shows the sparsity for every layer of Lenet-3. Here sparsity is defined in Sec. 4, number of weights denotes the total number of parameters in the designated layer. It is interesting that the Γ tends to put lower sparsity on layer with more parameters.

Table content: <table/>

provides the sparsity variations during SplitLBI training in Fig.5. This table shows the sparsity for every layer of Conv-2. Here sparsity is defined in Sec. 4, number of weights denotes the total number of parameters in the designated layer. The sparsity is more significant in fully connected (FC) layers than convolutional layers.

Table content: <table><tr><td>Layer</td><td colspan="2">Conv1 Conv2</td><td>FC1</td><td>FC2</td><td>FC3</td></tr><tr><td>Sparsity</td><td>0.9375</td><td>1</td><td>0.0067</td><td colspan="2">0.0284 0.1551</td></tr><tr><td>Number of Weights</td><td>576</td><td colspan="3">36864 3211264 65536</td><td>2560</td></tr></table>

This table shows the sparsity for every layer of Conv-4. Here sparsity is defined in Sec. 4, number of weights denotes the total number of parameters in the designated layer. Most of the convolutional layers are kept while the FC layers are very sparse.

Table content: <table><tr><td>Accuracy</td><td>98.75 99.00 99.25 99.50 99.75 100.00</td><td colspan="3">The Important Subnet Learned by SplitLBI SGD(Mom-Wd) SGD(Lottery) SplitLBI(Mom-Wd) SplitLBI(Lottery) F-epoch30 F-epoch60 F-epoch90 F-epoch100</td></tr><tr><td/><td>98.50</td><td/><td/></tr><tr><td/><td>98.25</td><td/><td/></tr><tr><td/><td>98.00</td><td>Lenet-3</td><td>Conv-2 Networks</td><td>Conv-4</td></tr></table>

Peer review (iclr2020):
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. The contributions of this paper are as follows: 1. A novel optimization based method for finding important sparse structure of large-scale neural networks using the idea of coupling the learning of weight matrix and sparsity constraints. 2. The convergence guarantee of the proposed method for solving nonconvex optimization is established based on a novel Lyapunov function. 3.Extensive experiments show that the proposed algorithm can: (a) train over-parameterized neural networks to achieve the state-of-the-art performances on image classification tasks; (b) find sparse networks with competitive networks; and (c) early stopping plus retrain from scratch can achieve similar or better performance than existing compression method. I believe that this work is very interesting and will be of interest to the ICLR community. The presentation of this paper is clear and the idea of coupling the gradient descent and mirror descent is very interesting. In addition, the empirical evaluations, including the sensitivity of the parameters and computation and memory costs, of the proposed method are solid. Minor comments: 1.There is no definition of KL function in the main content. 2.Before presenting Corollary 1, please specify the g_k as well as the updates for training neural networks. 3.Is it possible to extend the convergence guarantee to the ReLU activation function since the architectures considered in experiments using such activation. 4.In experiments it is better to present the convergence of the training loss to validate the convergence results. 5.In Figure 2, what is SLBI-10 and SLBI-1? Update: Thanks for your clarification, and I recommend acceptance of this work.