Exploring Structural Sparsity of Deep Networks via Inverse Scale Spaces

TL;DR

This paper introduces a novel differential inclusion-based method called DessiLBI for exploring structural sparsity in deep neural networks, enabling efficient learning of sparse models with competitive performance and early identification of effective subnetworks.

Contribution

It proposes a new inverse scale space approach with a discretization called DessiLBI, providing a systematic way to learn sparse deep network structures with theoretical convergence guarantees.

Findings

Achieves comparable or better performance than existing optimizers.

Unveils 'winning tickets' early in training with comparable accuracy.

Enables efficient network growth with less computational cost.

Abstract

The great success of deep neural networks is built upon their over-parameterization, which smooths the optimization landscape without degrading the generalization ability. Despite the benefits of over-parameterization, a huge amount of parameters makes deep networks cumbersome in daily life applications. Though techniques such as pruning and distillation are developed, they are expensive in fully training a dense network as backward selection methods, and there is still a void on systematically exploring forward selection methods for learning structural sparsity in deep networks. To fill in this gap, this paper proposes a new approach based on differential inclusions of inverse scale spaces, which generate a family of models from simple to complex ones along the dynamics via coupling a pair of parameters, such that over-parameterized deep models and their structural sparsity can be explored simultaneously. This kind of differential inclusion scheme has a simple discretization, dubbed Deep structure splitting Linearized Bregman Iteration (DessiLBI), whose global convergence in learning deep networks could be established under the Kurdyka-Lojasiewicz framework. Experimental evidence shows that our method achieves comparable and even better performance than the competitive optimizers in exploring the sparse structure of several widely used backbones on the benchmark datasets. Remarkably, with early stopping, our method unveils `winning tickets' in early epochs: the effective sparse network structures with comparable test accuracy to fully trained over-parameterized models, that are further transferable to similar alternative tasks. Furthermore, our method is able to grow networks efficiently with adaptive filter configurations, demonstrating a good performance with much less computational cost. Codes and models can be downloaded at {https://github.com/DessiLBI2020/DessiLBI}.

Tables

Table 1. TABLE I: Top-1/Top-5 accuracy(%) on ImageNet-2012. All models are run in 100 epochs. ⋆ : results from the official pytorch website. We use the official pytorch codes to run the competitors. More results on MNIST/CIFAR10, please refer Tab. VI in supplementary.

Dataset ImageNet-2012 Models Variants AlexNet ResNet-18 SGD Naive –/– 60.76/79.18 $l_1$ 46.49/65.45 51.49/72.45 Mom 55.14/78.09 66.98/86.97 Mom-Wd⋆ 56.55/79.09 69.76/89.18 Nesterov -/- 70.19/89.30 Adam Naive –/– 59.66/83.28 DessiLBI Naive 55.06/77.69 65.26/86.57 Mom 56.23/78.48 68.55/87.85 Mom-Wd 57.09/79.86 70.55/89.56

Table 2. TABLE II: The performance of the sparse structure under weight pruning setting found via our method. Here FT n 𝑛 n and RT n 𝑛 n refer to fine-tuning and retraining using the weights and structure after training with n 𝑛 n epoch, respectively. In particular, FT* 200 means training for 200 epochs and using fine-tune for only 10 epochs. The definition of sparsity is the number of non-zero parameters divided by the number of whole parameters. Sparse Acc denotes the accuracy of pruned models.

Method		ResNet50			Method		VGG16
		Acc	Sparse Acc	Sparisty			Acc	Sparse Acc	Sparsity
Iterative-Pruning-A [116]		92.67	93.08	0.40	Iterative-Pruning-A [116]		93.64	93.58	0.40
Iterative-Pruning-A [116]		92.67	85.16	0.10	Iterative-Pruning-A [116]		93.64	91.88	0.10
Iterative-Pruning-B [117]		93.02	92.14	0.40	Iterative-Pruning-B [117]		93.77	93.74	0.40
Iterative-Pruning-B [117]		93.02	73.75	0.10	Iterative-Pruning-B [117]		93.77	91.09	0.10
DessiLBI	FT 60	94.13	91.70	0.09	DessiLBI	FT 60	93.45	93.80	0.05
	RT 60	94.13	90.68	0.09		RT 60	93.45	93.08	0.05
	FT 120	94.13	92.64	0.11		FT 120	93.45	94.18	0.06
	RT 120	94.13	91.35	0.11		RT 120	93.45	93.60	0.06
	FT 200	94.13	92.54	0.11		FT 200	93.45	93.89	0.06
	FT* 200	94.13	90.03	0.11		FT* 200	93.45	91.11	0.06
	RT 200	94.13	91.48	0.11		RT 200	93.45	93.20	0.06

Table 3. TABLE III: The performance of the sparse structure under filter pruning setting found via our method. Here FT n 𝑛 n and RT n 𝑛 n refer to fine-tuning and retraining using the weights and structure after training with n 𝑛 n epoch, respectively. In particular, FT* 200 means training for 200 epochs and using fine-tune for only 10 epochs. Layer Selection contains experiments of pruning layers. MFLOP COUNT indicates the total number of floating-point operations executed in millions. The percentage means the ratio between the current MFLOP count and the MFLOP count of the full model. Sparse Acc denotes the accuracy of pruned model.

Method		ResNet56			Method		VGG16
		Acc	Sparse Acc	MFLOP COUNT			Acc	Sparse Acc	MFLOP COUNT
LEGR [82]		93.90	93.70	59 (47%)	Hinge [123]		93.25	92.91	191 (61%)
AMC [120]		92.80	91.90	63 (50%)	Variational [122]		93.25	93.18	190 (61%)
SFP [27]		93.60	93.40	59 (47%)	L1 [121]		93.25	93.11	200 (64%)
HRank [118]		93.26	93.17	63 (50%)	HRank [118]		93.96	93.43	146 (54%)
DessiLBI	FT 60	93.46	93.30	77 (61%)	DessiLBI	FT 60	93.43	93.00	91 (29%)
	RT 60	93.46	93.31	77 (61%)		RT 60	93.43	93.01	91 (29%)
	FT 120	93.46	93.59	78 (62%)		FT 120	93.43	93.21	97 (31%)
	RT 120	93.46	93.16	78 (62%)		RT 120	93.43	92.92	97 (31%)
	FT 200	93.46	93.94	78 (62%)		FT 200	93.43	93.59	100 (32%)
	FT* 200	93.46	89.97	78 (62%)		FT* 200	93.43	89.76	100 (32%)
	RT 200	93.46	93.15	78 (62%)		RT 200	93.43	93.01	100 (32%)
	Layer Selection	93.73	93.47	56 (45%)		Layer Selection	93.47	94.06	106 (34%)

Table 4. TABLE IV: Filter configuration of model pruned by L1-prune and model learned by our growing method on CIFAR10 dataset. L 1 subscript 𝐿 1 L_{1} -P is the L 1 subscript 𝐿 1 L_{1} pruning. The fully connected layer is altered according to   [ 121 ] .

		Model	$L_1$-P [121]	Ours
VGG16	Output size	Maps	Maps	Maps
Conv_1	32$\times$32	64	32	32
Conv_2	32$\times$32	64	64	96
Conv_3	16$\times$16	128	128	126
Conv_4	16$\times$16	128	128	220
Conv_5	8$\times$8	256	256	242
Conv_6	8$\times$8	256	256	200
Conv_7	8$\times$8	256	256	182
Conv_8	4$\times$4	512	256	261
Conv_9	4$\times$4	512	256	151
Conv_10	4$\times$4	512	256	105
Conv_11	2$\times$2	512	256	88
Conv_12	2$\times$2	512	256	88
Conv_13	2$\times$2	512	256	80
FC1	-	512	512	512
FC2	-	10	10	10
Params	-	15 M	5.40 M	2.96 M
FLOPs	-	313 M	206 M	217 M
Acc.	-	93.25%	93.40%	93.80%

Table 5. TABLE V: Comparison between AutoGrow  [ 124 ] and our growing method.

Dataset	Method	Params.	Acc(%)
CIFAR10	AutoGrow [124]	4.06 M	94.27
	Ours	2.69 M	94.82
CIFAR100	AutoGrow [124]	5.13 M	74.72
	Ours	3.37 M	76.86

Table 6. TABLE VI: Top-1/Top-5 accuracy(%) on ImageNet-2012 and test accuracy on MNIST/CIFAR10. ⋆ : 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.

Dataset MNIST CIFAR10 ImageNet-2012 Models Variants LeNet ResNet-20 AlexNet ResNet-18 SGD Naive 98.87 86.46 –/– 60.76/79.18 $l_1$ 98.52 67.60 46.49/65.45 51.49/72.45 Mom 99.16 89.44 55.14/78.09 66.98/86.97 Mom-Wd⋆ 99.23 90.31 56.55/79.09 69.76/89.18 Nesterov 99.23 90.18 -/- 70.19/89.30 Adam Naive 99.19 89.14 –/– 59.66/83.28 Adabound 99.15 87.89 –/– –/– Adagrad 99.02 88.17 –/– –/– Amsgrad 99.14 88.68 –/– –/– Radam 99.08 88.44 –/– –/– DessiLBI Naive 99.02 89.26 55.06/77.69 65.26/86.57 Mom 99.19 89.72 56.23/78.48 68.55/87.85 Mom-Wd 99.20 89.95 57.09/79.86 70.55/89.56

Table 7. TABLE VII: Results of pruning layers on CIFAR10. MFLOP indicates the total number of floating-point operations executed in millions. No. D-L is short for the number of dropped layers. The percentage means the ratio between the current MFLOP count and the MFLOP COUNT of full model.

	Acc.	Sparse Acc.	Sparsity	MFLOP	No. D-L
ResNet56	93.73	93.47	0.4157	56 (44.89%)	2
VGG16	93.47	94.06	0.1095	106(33.84%)	6

Table 8. TABLE VIII: Computational and Memory Costs.

optimizer	SGD	DessiLBI	Adam
Mean Batch Time	0.0197	0.0221	0.0210
GPU Memory	1161MB	1459MB	1267MB

Table 9. TABLE IX: This table shows the comparison between our methods and the selected growing methods.Here accuracy denotes the testing accuracy, number of parameters means the number of parameters of the obtained model and running time means the training and growing time.

Methods	Accuracy	number of parameter	running time(hour)		Methods	accuracy	number of parameter	running time(hour)
Splitting  [126]	87.17	0.07M	3.26h		Splitting [126]	92.64	0.51M	8.45h
EA-Splitting [127]	88.99	0.12M	3.03h		EA-Splitting [127]	93.19	0.86M	7.38h
Firefly [128]	88.77	0.10M	2.78h		Firefly [128]	92.77	0.77M	6.53h
Nash [125]	89.01	0.13M	4.13h		Nash [125]	92.52	0.75M	11.26h
Ours	88.54	0.11M	0.72h		Ours	93.62	0.82M	2.44h
Splitting [126]	89.47	0.13M	4.80h		Splitting [126]	92.70	0.86M	10.06h
EA-Splitting [127]	91.15	0.23M	4.43h		EA-Splitting [127]	93.46	1.60M	8.90h
Firefly [128]	90.66	0.21M	4.02h		Firefly [128]	93.18	1.47M	7.96h
Nash [125]	90.40	0.23M	6.53h		Nash [125]	93.10	1.17M	13.63h
Ours	91.05	0.20M	1.03h		Ours	93.72	1.41M	4.33h
Splitting [126]	91.55	0.25M	6.55h		Splitting [126]	93.40	1.64M	13.3h
EA-Splitting [127]	91.99	0.45M	5.88h		EA-Splitting [127]	93.56	2.89M	10.55h
Firefly [128]	92.07	0.40M	5.21h		Firefly [128]	93.86	2.74M	9.55h
Nash [125]	91.38	0.44M	8.96h		Nash [125]	93.53	2.11M	16.45h
Ours	92.49	0.40M	1.51h		Ours	94.03	2.00M	5.31h

Table 10. TABLE X: This table shows the detailed structure comparison between our method and Firefly  [ 128 ] . We present the last two model of the two methods.

	Firefly	Ours	Firefly	Ours
Layer1	68	17	73	17
Layer2	122	59	192	59
Layer3	128	102	193	102
Layer4	112	180	161	186
Layer5	142	220	207	308
Layer6	149	222	199	284
Layer7	134	122	180	128
Layer8	136	72	169	72
Layer9	115	48	147	48
Layer10	99	28	118	28
Layer11	73	20	90	20
Layer12	54	18	70	18
Layer13	48	12	69	12
Layer14	37	12	64	12
Layer15	62	12	90	12
Layer16	78	18	80	18
number of parameters	1.47M	1.41M	2.74M	2.00M

Equations

$$V_{k+1}$$

$$W_{k+1}$$

$${\dot{V}}_t$$

$$W_t$$

$${\Phi }_W:\mathcal{X}\to \mathcal{Y}$$

$${\hat{\mathcal{L}}}_n(W)=\frac{1}{n}\sum _{i=1}^n\ell (y_i,{\Phi }_W(x_i)).$$

$${\Phi }_W(x)={\sigma }_l\left(W^l{\sigma }_{l-1}\left(W^{l-1}\cdots {\sigma }_1\left(W^{1}x\right)\right)\right)$$

$$\frac{\dot{W_t}}{\kappa }$$

$$\dot{V_t}$$

$$V_t$$

$$\bar{\Omega }(\Gamma ):={\Omega }_{\lambda }(\Gamma )+\frac{1}{2\kappa }\Vert \Gamma {\Vert }^2$$

$$\bar{\mathcal{L}}\left(W,\Gamma \right)={\hat{\mathcal{L}}}_n\left(W\right)+\frac{1}{2\nu }\Vert W-\Gamma {\Vert }_F^2,$$

$$\underset{W}{\min }{\hat{\mathcal{R}}}_n(W):={\hat{\mathcal{L}}}_n\left(W\right)+$$

$$W_{k+1}=W_k-\kappa {\alpha }_k\cdot {\nabla }_W\bar{\mathcal{L}}\left(W_k,{\Gamma }_k\right),$$

$$V_{k+1}=V_k-{\alpha }_k\cdot {\nabla }_{\Gamma }\bar{\mathcal{L}}\left(W_k,{\Gamma }_k\right),$$

$${\Gamma }_{k+1}=\kappa \cdot {\mathrm{Prox}}_{{\Omega }_{\lambda }}\left(V_{k+1}\right),$$

$${\mathrm{Prox}}_{{\Omega }_{\lambda }}(V)=\arg \underset{\Gamma }{\min }\left\{\frac{1}{2}\Vert \Gamma -V{\Vert }_2^2+{\Omega }_{\lambda }\left(\Gamma \right)\right\},$$

$$P_{k+1}=\arg \underset{P}{\min }\left\{⟨P-P_k,\alpha \nabla \bar{\mathcal{L}}(P_k)⟩+B_{\Psi }^{p_k}(P,P_k)\right\}\text{vspace}-0.1in$$

$$\text{vspace}-0.1in\Psi (P)$$

$$={\Omega }_{\lambda }(\Gamma )+\frac{1}{2\kappa }\Vert W{\Vert }_2^2+\frac{1}{2\kappa }\Vert \Gamma {\Vert }_2^2,$$

$$B_{\Psi }^q(P,Q)$$

$$F(P,\tilde{g}):=\alpha \bar{\mathcal{L}}(W,\Gamma )+B_{\Omega }^{\tilde{g}}(\Gamma ,\tilde{\Gamma }),$$

$${\Omega }^{\ast }(g):=\underset{U\in {\mathbb{R}}^n}{\sup }\{⟨U,g⟩-\Omega (U)\}.$$

$$F(W,\Gamma ,g)=\alpha \bar{\mathcal{L}}(W,\Gamma )+\Omega (\Gamma )+{\Omega }^{\ast }(g)-⟨\Gamma ,g⟩.\text{vspace}-0.1in$$

$$0<{\alpha }_k=\alpha <\frac{2}{\kappa (Lip+{\nu }^{-1})},\text{vspace}-0.1in$$

$${\tilde{\nabla }}_W^t={\nabla }_W\bar{\mathcal{L}}\left(W\right)\mid {}_{(\mathbf{X},\mathbf{Y}{)}_{{\mathrm{batch}}_t}}.$$

$$v_{t+1}$$

$$W_{t+1}$$

$$v_{t+1}$$

$$W_{t+1}$$

$$V_{k+1}=V_k-{\alpha }_k{\beta }_k\odot {\nabla }_{\Gamma }\bar{\mathcal{L}}\left(W_k,{\Gamma }_k\right),$$

$${\Gamma }_{k+1}=\kappa {ϵ}_k\odot {\mathrm{Prox}}_{{\Omega }_{\lambda }}\left(V_{k+1}\right),$$

$${\beta }_{k,j}^i=\min \left(1,\frac{1}{\Vert W_{k,j}^i{\Vert }_{\star }}\right)\cdot \left(1-\frac{\#(\mathrm{supp}({\Gamma }_k^i))}{\#(\mathrm{supp}(W_k^i))}\right),$$

$$M^i=|\mathrm{sign}({\Gamma }^i)|=1_{\mathrm{supp}({\Gamma }^i)}$$

$$M^{i,g}=|\mathrm{sign}({\Gamma }^{i,g})|=1_{\mathrm{supp}({\Gamma }^{i,g})}$$

Full text

1 Introduction

Nowadays deep neural networks have shown great expressive power in many research areas such as image recognition [1], object detection [2], and point cloud estimation [3]. Such power is attributed to an avalanche of network parameters learned by supervision on large-scale datasets, i.e., model over-parameterization. Typically, the total number of parameters is orders of the magnitude higher than the number of training samples. And the over-parameterized neural networks can be trained with the loss functions by Stochastic Gradient Descent (SGD) [4] or modified optimization methods with adaptive stepsize, e.g., Adam [5], accompanied by early stopping.

The over-parameterization can benefit the training process of deep neural networks (DNN), and not necessarily result in a bad generalization or overfitting [6], especially when some weight size dependent complexities are controlled [7, 8, 9, 10]. Particularly, some recent empirical works show that model over-parameterization may help both optimization and generalization of networks, by simplifying the optimization landscape of empirical risks toward locating global optima [11, 12, 13, 14], and improving the generalization ability of deep neural networks for both discriminative [6] and generative models [15].

However, compressive networks are desired in many real world applications, e.g. robotics, self-driving cars, and augmented reality. For instance, the inference of large DNN models typically demands the support of GPUs, which are expensive for many real-world applications. Thus, it is essential to produce compressive networks. For this purpose, the classical way is to employ the norm-based regularization such as $L_{1}$ regularization [16] and enforce the sparsity on weights toward the compact, and memory efficient networks. This type of methods, unfortunately, may cause the decline of expressive power as empirically validated in [17]. This is because that the weights learned in neural networks are highly correlated, and $L_{1}$ regularization on such weights violates the incoherence or irrepresentable condition needed for sparse model selection [18, 19, 20], leading to spurious selections with poor generalization. On the other hand, the general type of regularization such as $L_{2}$ norm typically takes the function of low-pass filtering, sometimes in the form of weight decay [21] or early stopping [22, 23]. Sparsity has not been explicitly enforced on the models in this regularization, which may not produce a compressive model directly. Alternatively, Group Lasso [24] has also been utilized for finding sparse structures in DNN [25], and exerting good data locality with structured sparsity [25].

The difficulty of efficiently training a sparse network without over-parameterization results in the common practice in the community resorting to the backward selection, i.e., starting from training a big model using common task datasets like ImageNet, and then conduct the pruning [26, 27, 28, 29] or distilling [30] such big models to small ones without sacrificing too much of the performance. In particular, the recent Lottery Ticket Hypothesis (LTH) proposed in [31] made the following key empirical observation: dense, randomly-initialized, feed-forward networks contain small, sparse subnetworks, i.e., “winning tickets” structures, capable of being trained to comparable performance as the original network at a similar speed. To find such winning tickets, LTH works in backward selection, relying on the methods of one-shot or iterative pruning, which however, demands expensive computations and rewinding from initializations in  [32, 33].

Is there any alternative approach to find effective subnetworks without fully training a dense network? In this paper, we pursue the methodology in a reverse order, forward selection, a sharp contrast to the backward selection methods above. Particularly, we design some dynamics that starts from simple yet interpretable models, delving into complex models progressively, and simultaneously exploits over-parameterized models and structural sparsity. Our forward selection method enables finding important structural sparsity even before fully training a dense, over-parameterized model, avoiding the expensive computations in backward selection.

To achieve this goal, the Inverse Scale Space (ISS) method in applied mathematics [34, 35] is introduced to training deep neural networks, for the first time up to our knowledge. The “inverse scale space” method, was firstly proposed in [36] with Total-Variation sparsity for image reconstruction. The name comes from the fact that the features in the inverse scale space shown early in small scales are coarse-grained shapes, while fine details appeared later, in a reverse order of wavelet scale space where coarse-grained features appear in large scale spaces. Recently the ISS was shown as sparse regularization paths with statistical model selection consistency in high dimensional linear regression [35] and generalized linear models [37]. Moreover, Huang et al. [38, 39] further improved this by relaxing model selection consistency conditions using variable splitting.

Our inverse scale space dynamics of training neural networks can be described as differential inclusions, where important network parameters are learned at a faster speed than unimportant ones. Specifically, 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 $\Gamma$ learning structure sparsity in an inverse scale space. The two sets of parameters are coupled in an $L_{2}$ regularization. The ISS follows the gradient descent flow when the coupling regularization is weak, while reduces to a sparse mirror descent flow when the coupling is strong. During the training process, the parameters $\Gamma$ plays the role of exploring the sparse structure of the model parameters in inverse scale space, where important structures are learned faster than unimportant ones.

Such differential inclusion dynamics enjoy a simple discretization even in a highly non-convex setting of training deep neural networks, where we call such a discretization as Deep structure splitting Linearized Bregman Iteration (DessiLBI). A proof is provided to guarantee the global convergence of DessiLBI under the Kurdyka-Łojasiewicz framework. DessiLBI is a natural extension of SGD with sparse structure exploration in an inverse scale space. Critically, DessiLBI finds the important structure faster than unimportant ones, which enables a totally new way of exploring and exploiting the compact structure in DNNs. This paper presents the applications of DessiLBI in network sparsification, finding winning tickets, and growing networks. Particularly, we address: (1) how to find sparse network structures directly from our augmented variables $\Gamma$ computed by DessiLBI, where in particular, DessiLBI can help find the winning tickets without the expensive rewinding; and (2) an effective way to grow networks from a simple seed network to complex ones.

Network sparsification: DessiLBI may find sparse network structures that effective subnetworks can be rapidly learned via the structural sparsity parameter along the early iterative dynamics without fully training a dense network first. The support set of structural sparsity parameter learned in the early stage of this inverse scale space discloses important sparse subnetworks, including important weights, filters, and even layers. After obtaining the important sparse structures, both fine-tuning and retraining can be selected as post-processing. The priority of fine-tuning and retraining received wide discussions recently [40, 41], while we conduct extensive experiments to compare their performances. Our experimental results illustrate that the sparse structure found by DessiLBI is relatively robust to post-processing. Sparse structure found along the regularization path shows good performance on several widely-used network structures compared with their dense counterparts. In addition, training with DessiLBI does no harm to, or even enhances the performance of the dense model. As a result, the structural sparsity parameter may enable us to rapidly find sparse structure in early training epochs which saves plenty of training time and computational cost.

Finding winning tickets: DessiLBI also demonstrates new inspiring performance on finding a winning ticket structure in LTH [31]. We conduct several experiments to explore the performance of winning ticket subnetworks found by our methods. These experiments show that our method can find a winning ticket subnetwork at an early stage, while having similar or even better generalization ability if compared against fully trained dense models. Besides, experiments also show that the winning tickets obtained by our method generalize across different natural image datasets, exhibiting transferability as studied in [33].

Growing networks: DessiLBI can result in an elegant way to grow network dynamically. Here we propose to use regularization paths in inverse scale spaces to construct a lite growing method. In detail, we start with a small seed network with only a few filters for each layer. During the exploration of inverse scale spaces, important parameters are selected at an early stage. When the majority of the filters in one layer are selected, we assume that the complexity of this layer should be increased to enhance the model capacity, and more filters will be added to this layer. The early stopping property of DessiLBI will greatly reduce computational cost while maintain the model performance.

Contributions. We highlight the contributions in this paper. (1) The Inverse Scale Space method is, for the first time, applied to explore the structural sparsity of over-parameterized deep networks. DessiLBI can be interpreted as the discretization of solution paths of differential inclusion dynamics for the inverse scale spaces. (2) Global convergence of DessiLBI 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 DessiLBI demonstrate comparable and even better performance than other training algorithms on ResNet-18 in large scale training such as ImageNet-2012, among other datasets, jointly exploring structural sparsity with interpretability. (4) Structural sparsity parameters in DessiLBI provide important information about subnetwork architecture with comparable or even better accuracies than dense models after retraining or finetuning – DessiLBI with early stopping can provide fast winning tickets without fully training dense, over-parameterized models. (5) By using DessiLBI, we present two elegant ways to explore compact model: selecting important structures in the original network and expanding a seed network to ones with sufficient capacity.

Extensions. We explain the extension from our conference paper [42]. (1) Fundamentally, despite the essential idea is still the same as [42], we equip DessiLBI with a new magnitude scaling update strategy, that significantly alleviates the imbalance of magnitude scales across different layers in DNNs, as empirically validated in our experiments. (2) By using the DessiLBI updated from our conference version, we further propose a series of ways to pruning the neural network including weight pruning, filter pruning and our novel layer pruning. (3) We further study the properties of early stopping and transferability of the winning tickets found by DessiLBI. (4) We propose an elegant way to dynamically grow a network via exploring the inverse scale space which needs much less training time and computational cost compared with other methods. (5) Extensive new experiments and ablation studies that are added in addition to our conference version, further reveal the insights and efficacy of our methods.

2 Related Works

Our DessiLBI is built upon the Linearized Bregman Iterations and has a tight relationship to the classical mirror descent algorithm and ADMM. Some other related topics of finding sparse networks are also discussed.

2.1 Mirror Descent Algorithm

Mirror Descent Algorithm (MDA) firstly proposed by [43] to solve constrained convex optimization $L^{\star}:=\min_{W\in K}\mathcal{L}(W)$ ($K$ is convex and compact), can be understood as a generalized projected gradient descent [44] with respect to Bregman distance $B_{\Omega}(u_{0},u_{1}):=\Omega(u_{0})-\Omega(u_{1})-\langle\nabla\Omega(u_{1}),u_{0}-u_{1}\rangle$ induced by a convex and differentiable function $\Omega(\cdot)$,

[TABLE]

where the conjugate function of $\Omega(\cdot)$ is defined as $\Omega^{\star}(V):=\sup_{W}\langle W,V\rangle-\Omega(W)$.

At the $k$ iteration, Equation (1) uses two steps to optimize $W_{k+1}=\arg\min_{v}\langle v,\alpha\mathcal{L}(W_{k})\rangle+B_{\Omega}(v,W_{k})$ [45] : Eq (1a) implements the gradient descent on $V$ that is an element in dual space $V_{k}=\nabla\Omega(W_{k})$; and Eq (1b) projects it back to the primal space. As step size $\alpha\to 0$, MDA has the following limit dynamics as ordinary differential equation (ODE) [43]:

[TABLE]

where $\dot{V}_{t}$ denotes the right derivative of $V_{t}$ at the time $t>0$.

Convergence analysis with rates for convex loss has been well studied. Researchers also extend the analysis to stochastic version [46, 47] and Nesterov acceleration scheme [48, 49]. In deep learning, we have to deal with highly non-convex loss, recent work [50] has established the convergence to global optima for overparameterized under two assumptions: (i) the initial point is close enough to the manifold of global optima; (ii) the $\Omega(\cdot)$ is strongly convex and differentiable.

For non-differentiable $\Omega$ such as the Elastic Net penalty in compressed sensing and high dimensional statistics ($\Omega(W)=\|W\|_{1}+\frac{1}{2\kappa}\|W\|_{F}^{2}$ with damping factor $\kappa>0$), Equation (1) is studied as the Linearized Bregman Iteration (LBI) in applied mathematics [51, 35] that follows a discretized solution path of differential inclusions, to be discussed below. Such solution paths play a role of sparse regularization path where early stopped solutions are often better than the convergent ones when noise is present. In this paper, we investigate a varied form of LBI for the highly non-convex loss in deep learning models, exploiting the sparse paths, and establishing its convergence to a Karush–Kuhn–Tucker (KKT) point for general networks from arbitrary initializations. Furthermore, our method is a natural extension of SGD with sparse structure exploration. It reduces to the standard gradient descent method when the coupling regularization is weak, while reduces to a sparse mirror descent when the coupling is strong.

2.2 *Linearized Bregman Iteration *

Linearized Bregman Iteration (LBI), was proposed in [52, 51] that firstly studies Eq. (1) when $\Omega(W)$ involves $\ell_{1}$ or total variation non-differentiable penalties met in compressed sensing and image denoising. Beyond convergence for convex loss [51, 53], Osher et al. [35] and Huang et al. [37], particularly showed that LBI is a discretization of differential inclusion dynamics whose solutions generate iterative sparse regularization paths, and established the statistical model selection consistency for high-dimensional generalized linear models. Moreover, Huang et al. [38, 39] further improved this by proposing SplitLBI, incorporating into LBI a variable splitting strategy such that the restricted Hessian with respect to augmented variable ($\Gamma$ in Eq. 6) is orthogonal. This can alleviate the multicollinearity problem when the features are highly correlated; and thus can relax the irrepresentable condition,* i.e.*, the necessary condition for Lasso to have model selection consistency [19, 20, 54]. A variety of applications (e.g. [55, 56, 57, 58, 59, 60]) have been found for this algorithm since its inception.

However, existing work on SplitLBI is restricted to convex problems in generalized linear modes. It remains unknown whether the algorithm can exploit the structural sparsity in highly non-convex deep networks. To fill in this gap, in this paper, we propose the Deep structure splitting LBI that simultaneously explores the overparameterized networks and the structural sparsity of parameters in such networks, which enables us to generate an iterative regularization path of deep models whose important sparse architectures are unveiled in early stopping.

2.3 Alternating Direction Method of Multipliers

Alternating Direction Method of Multipliers (ADMM) which also adopted variable splitting strategy, breaks original complex loss into smaller pieces with each one can be easily solved iteratively [61, 62]. Recent works [63] established the convergence result of ADMM in convex, stochastic and non-convex setting, respectively. Wang et al. [64] studied convergence analysis with respect to Bregman distance. Training neural networks by ADMM has been studied in [65, 66]. Recently, Wang et al. [67] established the convergence of ADMM in a very general nonconvex setting, and Franca et al. [68] derived the limit ODE dynamics of ADMM for convergent analysis.

However, one should distinguish the LBI dynamics from ADMM that LBI should be viewed as a discretization of differential inclusion of inverse scale space that generalizes a sparse regularization solution path from simple to complex models where early stopping helps find important sparse models; in a contrast, the ADMM, as an optimization algorithm for a given objective function, focuses on the convergent property of the iterations.

2.4 Early Stopping Regularization of Gradient Method

To optimize the not differentiable target function, subgradient is typically utilized as the generalization of the gradients in classical optimization textbooks [69, 70].

The early stopping is commonly used as a regularization technique to avoid overfitting. It has been studied in mathematics as a regularization method in inverse problems [71]. In statistical machine learning, early stopping has been studied in Boosting as gradient descent method, e.g. $L_{2}$-Boost [72] and Boosting in classification  [73, 74]. In particular $L_{2}$-Boost is generalized by  [22] to gradient descent learning of regression functions in Reproducing Kernel Hilbert Spaces (RKHS) with random designs, showing that early stopping is a polynomial regularization better than Tikhonov or Ridge regularization for avoiding the saturation issue of the latter. A generalization of such kernel boosting algorithms to convex losses is given in [23, 75] using localized Rademacher complexities. Nonetheless, these works do not take sparsity constraints into consideration as high dimensional statistics and our work here. To handle sparse linear regressions, [76], establishes the the inverse scale space approach using differential inclusions, showing the model selection consistency for discovering causal variables by early stopping under the Irrepresentable or Incoherence condition equivalent to that of LASSO [20, 54]. Such results are later extended by [77] to general convex losses including logistic regression and various graphical models. For general structural sparsity where parameters are sparse under a linear transformation, [78, 39] proposes Split LBI and its limit differential inclusions, establishing their model selection consistency under weaker conditions than the incoherence condition of generalized LASSO. This last work lays down a foundation of current exploration of early stopping to find important structure in deep neural networks.

2.5 Finding Compact Networks

Pruning Networks. In real world applications, limited computational resource makes compact models more demanding. In the manner of backward selection, pruning [79] is one of the most direct ways of producing a light model. Network pruning can be roughly categorized as weight and filter pruning, by whether removing some structural parameters such as convolutional filters. For unstructural weight pruning, Han et al. [26] drops small weights of a well-trained dense network. The filter pruning methods take into consideration the network sparsity, memory footprint, and computational cost. Several works [27, 80] attempt to train a network firstly and then prune the network according to specific metrics such as $L_{2}$ norm of weights. Centripetal-SGD [81] groups the weight by their initialization and forces the weights inner one group to have the same values. LEGR [82] studies the scale variation across layers and suggests that using genetic algorithm can find decent affine transformation, making $L_{2}$ normalization more reasonable. In contrast to these backward selection methods that are expensive in both memory and computational cost, forward selection methods have been widely used in traditional statistical machine learning. Examples of such forward selection or parsimonious learning include boosting as functional gradient descent [83, 84] or coordinate descent method to solve LASSO [16], regularization paths associated with $L_{2}$ (Tikhonov regularization or Ridge regression) [85] and/or $L_{1}$ (Lasso) penalties [86], etc. None of these methods has been systematically studied in deep learning. In [40], the authors investigated a greedy forward selection method to sparsify filters. In a contrast, we investigate the inverse scale space method to network sparsification and our DessiLBI can select structural sparsity in weight level, filter level, and even layer level with improved efficiency.

Lottery Ticket Hypothesis (LTH). LTH [31] states that one can find effective subnets in a well-trained dense model, i.e. winning tickets. When training a winning ticket subnetwork in isolation from the initialization of the dense network, its test accuracy can match the dense model test accuracy with at most the same number of iterations. By finding a winning ticket using backward-selection-based algorithms, one can obtain an extremely sparse network with comparable or even better generalization ability.

Frankle et al. [32] studied the LTH in the even deeper neural networks by rewinding, which means retrain from a very early stage of training (0.1% to 7%). LTH also appeared in reinforcement learning and natural language processing [87, 88, 89]. Additionally, winning ticket networks for filter pruning are also studied in several works [41, 42]. The transferability is studied in [33] that the winning tickets can be transferred across datasets as well as optimizers. Recently it found that utilizing gradient information can unveil a winning ticket in a randomly initialized network without modifying the weight values [90]. Furthermore, Malach et al. [91] proved that there exists a subnetwork with similar accuracy in a sufficiently over-parameterized neural network with bounded random weights. More theoretical analysis appeared in [92, 93, 94].

Although LTH has been widely studied, generating winning tickets could be computationally expensive. The one shot pruning and iterative pruning proposed in [31] are the de facto algorithms in winning ticket generation. One shot pruning utilizes a backward selection that firstly trains a network and prunes it according to weight magnitudes. Iterative pruning improves one shot pruning by conducting one shot pruning multiple times with a smaller pruning rate for each time which may cost unbearable GPU hours when finding extremely sparse winning tickets.

In this paper, our novel way to get winning tickets is built upon the forward selection, in the iterative procedure of DessiLBI. Critically, we show that by using early stopping, our DessiLBI can successfully discover the winning ticket structure without fully training a dense network. Empirically, we also show that our winning ticket structure has the nice property of transferability to the new dataset in the same domain of natural images.

Searching Network Structure. Plenty of recent efforts are made on searching for a good sparse network, rather than pruning. For example, Neural Architecture Search (NAS) [95, 96, 97, 98, 99] aims to search better network architectures automatically. The early works [95, 100] of NAS use reinforcement learning or evolutionary algorithms to search the whole network architectures, which need huge amount of computation cost. For example, Zoph et al. [95] spend more than 2000 GPU days in searching a network architecture on CIFAR10 dataset. Many works try to reduce the searching cost. Instead of searching the whole network architectures, Zoph et al. [96] and Zhong et al. [97] search cells and stack the searched cells into the complete networks, which greatly reduces the search space. However, these improved methods still need nearly 100 GPU days [97]. There are other strategies to improve the efficiency of the searching process. ENAS [101] and one-shot NAS [102, 103] share the parameters between child models, which significantly accelerates the searching process. DARTS [98] proposes gradient based methods that just use one or several GPU days for searching. Different from these NAS works, this paper presents a lite network growing method, benefiting from the forward selection in training networks by DessiLBI. In our method, we jointly grow the network structures, and train the network parameters by using our DessiLBI. It thus saves significant computational cost, while still maintains reasonably good performance.

3 Inverse Scale Spaces and DessiLBI

Supervised learning aims to learn a mapping

[TABLE]

from input space $\mathcal{X}$ to output space $\mathcal{Y}$, with a parameter $W$ such as weights in neural networks, by minimizing certain empirical loss function on training samples

[TABLE]

For example, a neural network of $l$-layer is defined as

[TABLE]

where $W=\{W^{i}\}_{i=1}^{l}$, and $\sigma_{i}$ is the nonlinear activation function of the $i$-th layer.

3.1 Differential Inclusion of Inverse Scale Space

Differential Inclusion of Inverse Scale Space. Consider the following dynamics adapted to neural network training,

[TABLE]

where $\dot{W_{t}}$ and $\dot{V_{t}}$ are the right derivatives of $W_{t}$ and $V_{t}$, individually. The $V$ is a sub-gradient of

[TABLE]

where for some sparsity-enforced, often non-differentiable regularization, we have $\Omega_{\lambda}(\Gamma)=\lambda\Omega_{1}(\Gamma)$ ($\lambda\in\mathbb{R}_{+}$) such as Lasso or group Lasso penalties for $\Omega_{1}(\Gamma)$; $\kappa>0$ is a damping parameter such that the solution path is continuous, and the augmented loss function is

[TABLE]

with $\nu>0$ controlling the gap admitted between $W$ and $\Gamma$. Compared to the original loss function $\widehat{\mathcal{L}}_{n}\left(W\right)$, our loss $\bar{\mathcal{L}}\left(W,\Gamma\right)$ additionally uses the variable splitting strategy by lifting the original neural network parameter $W$ to $(W,\Gamma)$ with $\Gamma$ modeling the structural sparsity of $W$. For simplicity, we assume $\bar{\mathcal{L}}$ is differentiable with respect to $W$ here, otherwise the gradient in Eq. (6a) is understood as subgradient and the equation becomes an inclusion.

Differential inclusion system (Eq. 6) is a coupling of gradient descent on $W$ with non-convex loss and mirror descent (LBI) of $\Gamma$ (Eq. 2) with non-differentiable sparse penalty. It may explore dense over-parameterized models $W_{t}$ in the proximity of structural parameter $\Gamma_{t}$ with gradient descent, while $\Gamma_{t}$ records important sparse model structures. Specifically, the solution path of $\Gamma_{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, Equation 7 takes the $\Omega_{\lambda}(\Gamma)=\|\Gamma\|_{1}$ and $\kappa\to\infty$ for simplicity, $V_{t}$ as the subgradient of $\bar{\Omega}_{t}$, undergoes a gradient descent flow before reaching the $\ell_{\infty}$-unit box, which implies that $\Gamma_{t}=0$ in this stage. The earlier a component in $V_{t}$ reaches the $\ell_{\infty}$-unit box, the earlier a corresponding component in $\Gamma_{t}$ becomes nonzero and rapidly evolves toward a critical point of $\bar{\mathcal{L}}$ under gradient flow. On the other hand, the $W_{t}$ follows the gradient descent with a standard $\ell_{2}$-regularization. Therefore, $W_{t}$ closely follows the dynamics of $\Gamma_{t}$ whose important parameters are selected.

Compared with directly enforcing a penalty function such as $\ell_{1}$ or $\ell_{2}$ regularization

[TABLE]

dynamics Eq. (6) can relax the irrepresentable conditions for model selection by Lasso [38], which can be violated for highly correlated weight parameters. The weight $W$, instead of directly being imposed with $\ell_{1}$-sparsity, adopts $\ell_{2}$-regularization in the proximity of the sparse path of $\Gamma$ that admits simultaneously exploring highly correlated parameters in over-parameterized models and sparse regularization.

The key insight lies in that differential inclusion of Eq. (6c) drives the important features in $\Gamma_{t}$ that earlier reaches the $\ell_{\infty}$-unit box to be selected earlier. Hence, the importance of features is related to the “time scale” of dynamic hitting time to the $\ell_{\infty}$ unit box, and such a time scale is inversely proportional to lasso regularization parameter $\lambda=1/t$ [35]. Such a differential inclusion is firstly studied in [34] with Total-Variation (TV) sparsity for image reconstruction, where important features in early dynamics are coarse-grained shapes with fine details appeared later. This is in contrast to wavelet scale space that coarse-grained features appear in large scale spaces, thus named “inverse scale space”. In this paper, we shall see that Eq. (6) inherits such an inverse scale space property empirically even for the highly nonconvex neural network training. Figure 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.

3.2 Deep Structure Splitting LBI

Deep Structure Splitting Linearized Bregman Iteration. Equation (6) admits an extremely simple discrete approximation, using Euler forward discretization of dynamics and called DessiLBI in the sequel:

[TABLE]

where $\alpha_{k}$ is the step size at the $k$ iteration; $V_{0}=\Gamma_{0}=0$, $W_{0}$ can be small random numbers such as Gaussian initialization. Here we add interpretation for $\alpha$ and $\kappa$ in Appendix. F. For some complex networks, it can be initialized as common setting. The proximal map in Eq. (10c) that controls the sparsity of $\Gamma$,

[TABLE]

Such an iterative procedure returns a sequence of sparse networks from simple to complex ones whose global convergence condition to be shown below, while solving Eq. (9) at various levels of $\lambda$ might not be tractable, especially for over-parameterized networks.

Structural Sparsity. Our DessiLBI explores structural sparsity in fully connected and convolutional layers, which can be unified in framework of group lasso penalty, $\Omega_{1}(\Gamma)=\sum_{g}\|\Gamma^{g}\|_{2}$, where $\|\Gamma^{g}\|_{2}=\sqrt{\sum_{i=1}^{\mid\Gamma^{g}\mid}\left(\Gamma_{i}^{g}\right)^{2}}$ and $\left|\Gamma^{g}\right|$ is the number of weights in $\Gamma^{g}$. Thus Eq. (10c) has a closed form solution $\Gamma^{g}=\kappa\cdot\max\left(0,1-1/\|V^{g}\|_{2}\right)V^{g}$. Typically,

• (a)

For a convolutional layer, $\Gamma^{g}=\Gamma^{g}(c_{in},c_{out},\mathtt{size})$ denote the convolutional filters where $\mathtt{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.

• (b)

For a fully connected layer, $\Gamma=\Gamma(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.

4 Global Convergence of DessiLBI

We present a theorem that guarantees the global convergence of DessiLBI, i.e. from any initialization, the DessiLBI sequence converges to a critical point of $\bar{\mathcal{L}}$. Our treatment extends the block coordinate descent (BCD) studied in [105], with a crucial difference being the mirror descent involved in DessiLBI. Instead of the splitting loss in BCD, a new Lyapunov function is developed here to meet the Kurdyka-Łojasiewicz property [106]. [107] studied the convergence of variable splitting method for single hidden layer networks with Gaussian inputs.

Let $P:=(W,\Gamma)$. Following [37], the DessiLBI algorithm in Eq. (10a-10c) can be rewritten as the following standard Linearized Bregman Iteration,

[TABLE]

where

[TABLE]

$p_{k}\in\partial\Psi(P_{k})$, and $B_{\Psi}^{q}$ is the Bregman divergence associated with convex function $\Psi$, defined by

[TABLE]

for some $q\in\partial\Psi(Q)$. Without the loss of generality, consider $\lambda=1$ in the sequel. One can establish the global convergence of DessiLBI under the following assumptions.

Assumption 1**.**

Suppose that:

• (a)

$\widehat{\mathcal{L}}_{n}(W)=\frac{1}{n}\sum_{i=1}^{n}\ell(y_{i},\Phi_{W}(x_{i}))$* is continuous differentiable and $\nabla\widehat{\mathcal{L}}_{n}$ is Lipschitz continuous with a positive constant $Lip$;*

• (b)

$\widehat{\mathcal{L}}_{n}(W)$* has bounded level sets;*

• (c)

$\widehat{\mathcal{L}}_{n}(W)$* is lower bounded (without loss of generality, we assume that the lower bound is [math]);*

• (d)

$\Omega$* is a proper lower semi-continuous convex function and has locally bounded subgradients, that is, for every compact set ${\cal S}\subset\mathbb{R}^{n}$, there exists a constant $C>0$ such that for all $\Gamma\in{\cal S}$ and all $g\in\partial\Omega(\Gamma)$, there holds $\|g\|\leq C$;*

• (e)

the Lyapunov function

[TABLE]

is a Kurdyka-Łojasiewicz function **[105, 108]** on any bounded set, where $B_{\Omega}^{\tilde{g}}(\Gamma,\tilde{\Gamma}):=\Omega(\Gamma)-\Omega(\tilde{\Gamma})-\langle\tilde{g},\Gamma-\tilde{\Gamma}\rangle$, $\tilde{\Gamma}\in\partial\Omega^{*}(\tilde{g})$, and $\Omega^{*}$ is the conjugate of $\Omega$ defined as

[TABLE]

Remark 1**.**

Assumption 1 (a)-(c) are regular in the analysis of nonconvex algorithm (see, [108] 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 $\ell_{1}$ norm, group $\ell_{1}$ norm, and every continuously differentiable penalties. By Eq. (15) and the definition of conjugate, the Lyapunov function $F$ can be rewritten as follows,

[TABLE]

Now we are ready to present the main theorem.

Theorem 1**.**

[Global Convergence of DessiLBI] Suppose that Assumption 1 holds. Let $\{(W_{k},\Gamma_{k})\}$ be the sequence generated by DessiLBI (Eq. (10a-10c)) with a finite initialization. If

[TABLE]

then $\{(W_{k},\Gamma_{k})\}$ converges to a critical point of $\bar{\mathcal{L}}$ defined in Eq. (8), and $\{W^{k}\}$ converges to a critical point of $\widehat{\mathcal{L}}_{n}(W)$.

Applying to the neural networks, typical examples are summarized in the following corollary.

Corollary 1**.**

Let $\{(W_{k},{\Gamma}_{k},g_{k})\}$ be a sequence generated by DessiLBI (26a-26c) for neural network training where (a) $\ell$ is any smooth definable loss function [109], such as the square loss $t^{2}$, exponential loss $e^{t}$, logistic loss $\log(1+e^{-t})$, and cross-entropy loss; (b) $\sigma_{i}$ is any smooth definable activation [109], such as linear activation $t$, sigmoid $\frac{1}{1+e^{-t}}$, hyperbolic tangent $\frac{e^{t}-e^{-t}}{e^{t}+e^{-t}}$, and softplus $\frac{1}{c}\log(1+e^{ct})$ for some $c>0$) as a smooth approximation of ReLU; (c) $\Omega$ is the group Lasso. Then the sequence $\{W_{k}\}$ converges to a stationary point of $\widehat{\mathcal{L}}_{n}(W)$ under the conditions of Theorem 1.

5 Implementation of DessiLBI

The DessiLBI is a natural extension of SGD with exploring sparse structures of network. As referring to the variants of SGD, we further introduce several variants of our DessiLBI, by learning with data batches, using momentum to accelerate the convergence, adding weight decay on the parameters $W$, as well as the magnitude scaling updates.

Batch DessiLBI. To train networks on large datasets, stochastic approximation of the gradients in DessiLBI over the mini-batch $(\mathbf{X},\mathbf{Y})_{{\mathrm{batch_{t}}}}$ is adopted to update the parameter $W$,

[TABLE]

DessiLBI with Momentum (Mom). Inspired by the variants of SGD, the momentum term can be also incorporated to the standard DessiLBI that leads to the following updates of $W$ by replacing Eq (10a) with,

[TABLE]

where $\tau$ is the momentum factor, empirically setting as 0.9.

**DessiLBI with Momentum and Weight Decay (Mom-Wd). ** The update formulation for Eq. (10a) is

[TABLE]

Magnitude Scaling Update Strategy. The imbalance of weight scales across different layers during training may degrade the performance of network sparsification and winning ticket generation. We visualize the distribution of $V$ in Figure 6 and Figure 5. The magnitude scale varies significantly across layers. It adds to the difficulties of finding good structure. On the other hand, for ReLU-based neural networks, the parameter magnitude scales of different layers do not negatively affect the network training, as the positive homogeneity of ReLU [110] guarantees that the scale has no difference to the class prediction. To this end, we introduce a magnitude scaling update strategy to improve our DessiLBI, such that the magnitude scales of $V$ are comparable across different layers. The detailed updating rule is to replace Eq. (10b, 10c) by the following:

[TABLE]

where $\beta_{k}=\left[\beta_{k,j}^{i}\right]$, and $\epsilon_{k}$ are scaling factor tensors. For the $i$-th layer, scaling factors are defined by

[TABLE]

and $\epsilon_{k}=\left[\|W_{k,j}^{i}\|_{\star}\right]$. Here, we introduce the notation of $W_{k,j}^{i}$ in Eq (21) to denote the $j$-th weight/filter of $k$-th iteration and $i$-th layer; for the $j$-th weight or filter in the $i$-th layer, $\|W_{k,j}^{i}\|_{\star}=\left|W_{k,j}^{i}\right|$ or $\|W_{k,j}^{i}\|_{\star}=\|W_{k,j}^{i}\|_{2}$, respectively ; $\#$ denotes the number of elements; and $\mathrm{supp}$ is the support of a set. Note that for $j$-th filter of $i$-th layer, we will expand the size of $\beta_{k,j}^{i}$ to match the size of the filter. To avoid $\beta_{k}=\mathbf{0}$, we will set a minimum value for it.

For the parameter $V$, its gradient is determined by the magnitude of the residue between corresponding parameter $W$ and $\Gamma$. Here, we normalize the gradient for $V$ by the magnitude of $W$, so that the gradient contains more direction information. For instance, one filter may be more important if the updating direction altered in lower frequency. To avoid the gradient explosion, we set the maximum of $\beta_{k}$ to be 1. In addition, we use the ratio between numbers of selected parameters and numbers of whole parameters in one layer as the penalty.

This idea is enlightened by Lipschitz renormalization/scaling for ReLu-activation networks  [110, 111] which is of positive homogeneity, i.e. invariant under a multiplication of a positive constant on the weight matrix such that its norm is renormalized to be close to a unit ball. These works attempt to avoid driving the model weights to infinity via the renomalization. Different from their motivation, we aim to alleviate the scale imbalance problem in neural networks when selecting structure. As mentioned in  [110, 112], relu-based neural networks own the property called positive homogeneity or rescaling invariance. It means that we can modify the weights of relu-based neural networks in an appropriate way that the prediction of the model is not changed. So the scale of magnitude for different layers can vary significantly and this property sets some obstacles for selecting important structures for the scale imbalance. In our work, we attempt to normalize the gradient to $V$ so that the update relies more on the direction instead of magnitude.

6 Applications of DessiLBI

6.1 *Network Sparsification *

In the process of training networks by DessiLBI, we are able to produce the sparse network directly. This enables us to explore the over-parameterization and structural sparsity simultaneously. We take the coupled parameters $\Gamma$ to record the structural sparsity information along the regularization paths. The support set of $\Gamma$ is utilized to conduct the pruning. The key principle is that it will be taken as less important parameters if the convolutional filters have the corresponding coupled parameter $\Gamma=0$ in the inverse scale space. We have three levels for pruning parameters in network, described below. After getting the compact one, some post-processing steps such as fine-tuning or retraining are used to improve the performance of a pruned network. In [41] the authors point out that the performance of pruned network relies on the structure heavily; and retraining them from scratch can get even better results111Potentially this is an arguable point, which is challenged by  [40]. In our work, we will validate this in the experiment.

Pruning Weights. We define $M^{i}$ as the mask for each parameter of the $i$-th layer by using support set $\Gamma^{i}$,

[TABLE]

We directly use the mask to remove some weight parameters by Hadamard (elementwise) product as $\widetilde{W}^{i}=M^{i}\odot W^{i}$. This will enable us to get a sparse network. Typically, such pruning is commonly adopted to prune the weights of fully connected layers.

Pruning Filters. We can remove the convolutional filters by

[TABLE]

where $\mathrm{supp}(\Gamma^{i,g})$ is the support set of filters $\Gamma_{g}$ on the $i$-th layer. Accordingly, the filters are changed as $\widetilde{W}^{i}=M^{i,g}\odot W^{i}$. Generally, the pruning is conducted on the convolutional layers to remove the filters.

Pruning Layers. Benefiting from the magnitude scaling update of DessiLBI, the weight magnitudes of different layers learned by our DessiLBI are comparable or balanced now. This enables us to select the sparse structure at a higher level of whole layers than the level of individual filters or weights. Specifically if $\mathrm{supp}(\Gamma^{i,g})$ or $\mathrm{supp}(\Gamma^{i})$ is the empty set for the $i$-th layer, we remove this layer totally from the network. Note that we will utilize the necessary transition layers, including pooling and concatenation operations, to make the consistent feature dimension after removing layers. Heavy redundancy of deep neural networks results in that several layers sometimes can be removed completely.

6.2 Finding Winning Tickets

As shown in LTH [31], finding winning tickets as subnets in a fully trained over-parameterized network is quite computationally expensive, especially using iterative pruning to unveil extremely sparse networks. In contrast, our DessiLBI empowers a much easier way to generate winning ticket structure. Here a forward selection method is implemented through exploring the Inverse Scale Space. By using the augmented variables of DessiLBI, we can obtain the winning ticket structure of a network in a more efficient way. Particularly, we search for the subnetwork in the inverse scale space computed by DessiLBI: the $\Gamma$ is utilized as the metric for evaluating the importance of weights. That is, the mask $M^{i}$ of each $i$-th layer, derived from the support set of $\Gamma^{i}$ is utilized to define the winning ticket structure, where the winning ticket can be obtained without fully training the network but using early stopping of DessiLBI. In fact, only a few epochs are required to get the winning ticket via DessiLBI in our experiments below. Fine-tuning or retraining after receiving the winning ticket structure exhibit comparable or better performance than the fully trained networks in a similar speed. Hence in a forward selection manner, our DessiLBI offers a more efficient and elegant way to get the winning ticket.

Transferability of winning tickets found in traditional ways has been observed in the empirical experiments in [33]. Here the transferability of winning tickets found by DessiLBI is further studied. In our experiment, for each target dataset, we find winning ticket structures by DessiLBI on different source datasets and show their performance close to the winning tickets generated on the same target dataset.

In a summary, we highlight the two key points in our winning ticket generation.

• (a)

Early stopping. We can define the winning ticket structures by exploiting the structural sparsity parameter computed in the early stage during the training process by DessiLBI, known as the early stopping regularization.

• (b)

Transferrability of winning ticket structure. The winning ticket structure discovered by early stopping of DessiLBI can generalize across a variety of datasets in the natural image domain.

6.3 * Growing Networks by DessiLBI *

By exploiting the inverse scale space by DessiLBI, the structure of our network can be dynamically altered during the training procedure. We present a novel network growing process. Specifically, we start with a simple initialized network and gradually increase its capacity by adding parameters along with the training process by DessiLBI. Starting from very few filters of each convolutional layer, our growing method requires not only efficiently optimizing the parameters of filters, but also adding more filters if the existing filters do not have enough capacity to model the distribution of training data. The key idea is to measure whether the network at the current training iteration has enough capacity.

To this end, we monitor in training, the difference between support set of model parameters $\mathrm{supp}(W^{i})$ and support set $\mathrm{supp}(\Gamma^{i})$ of its augmented parameter $\Gamma^{i}$ on the $i$-th layer of network,

[TABLE]

where $\mathrm{supp}$ is the support set and $\#$ denotes the number of elements in the set. When the ratio becomes close to 1, it means most of the filters in $W^{i}$ are important and they are selected into $\Gamma^{i}$, and thus we should expand the network capacity by adding new filters to the corresponding layer. In our network growing process, we start from a small number of filters (e.g. 4) of convolutional layers, more and more filters tend to be with non-zero values as the training iterates. Every $J$ epochs, we can compute the ratio of Eq. (24) and consider the following two scenarios.

• (a)

If this ratio is larger than a pre-set threshold (denote as $\tau$), we add some new filters to existing filters $W$, randomly initialized as [113], and $\Gamma$ will add corresponding dimensions, initialized as zeros;

• (b)

Otherwise, we will not grow any filter in this epoch.

Then we continue DessiLBI to learn all the weights from training data. This process is repeated until the loss does not change much, or the maximum number of epoch is reached. After stopping the growing process, the network model is trained for extra few epochs (typically 30) with a smaller learning rate to finish the training process.

7 Experiments

Dataset. In this section, four widely used datasets are selected: CIFAR10, CIFAR100, SVHN, and ImageNet. For CIFAR10, it owns 10 categories and each category has 5000 and 1000 images for training and testing respectively. CIFAR100 has 100 classes containing 500 training images and 100 testing images each. SVHN contains real-world number images from 0 to 9. In SVHN, there are 73257 images for training and 26032 images for testing. For CIFAR10, CIFAR100, and SVHN, the images have a spatial size of $32\times 32$. For ImageNet-2012, there are 1.28 million training images and 50k validation images. These images are sampled from 1000 classes.

This section introduces some stochastic variants of DessiLBI, followed by four sets of experiments revealing the key insights and properties of DessiLBI in exploring the structural sparsity of deep networks.

Implementation. Experiments are conducted over various backbones, e.g., LeNet, AlexNet, VGG, and ResNet. For MNIST and CIFAR10, the default hyper-parameters of DessiLBI are $\kappa=1$, $\nu=10$ and $\alpha_{k}$ is set as $0.1$, decreased by 1/10 every 30 epochs. In ImageNet-2012, the DessiLBI utilizes $\kappa=1$, $\nu=1000$, and $\alpha_{k}$ is initially set as 0.1, decays 1/10 every 30 epochs. We set $\lambda=1$ in Eq. (11) by default, unless otherwise specified. On MNIST and CIFAR10, we have batch size as 128; and for all methods, the batch size of ImageNet 2012 is 256. The standard data augmentation implemented in pytorch is applied to CIFAR10 and ImageNet-2012, as [114]. The weights of all models are initialized as [113]. In the experiments, we define the sparsity as percentage of non-zero parameters, i.e., the number of non-zero weights dividing the total number of weights in consideration. Runnable codes can be downloaded222https://github.com/DessiLBI2020/DessiLBI.

7.1 Image Classification Experiments

Settings. By referring the settings in [114], we conduct the experiments on ImageNet-2012. We compare variants of DessiLBI with that 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. (1) Naive SGD: the standard SGD with batch input. (2) SGD with $\mathit{l}_{1}$ penalty (Lasso). The $\mathit{l}_{1}$ norm is applied to penalize the weights of SGD by encouraging the sparsity of the learned model, with the regularization parameter of the $\mathit{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 $\mathit{l}_{2}$ weight decay with the coefficient weight $1e^{-4}$. (5) SGD with Nesterov (Nesterov): the SGD uses nesterov momentum 0.9. (6) Naive Adam: it refers to standard Adam333In Tab. VI of Appendix, we further give more results for Adabound, Adagrad, Amsgrad, and Radam, which, we found, are hard to be trained on ImageNet-2012 in practice..

As a natural extension of SGD with exploring sparse structure, DessiLBI can serve as a network training method. To validate this point, we training different backbones by DessiLBI on the large-scale dataset – ImageNet 2012. The results of image classification are shown in Tab. I. Our DessiLBI variants may achieve comparable or even better performance than SGD variants in 100 epochs, indicating the efficacy in learning dense, over-parameterized models. It verifies that our DessiLBI can explore the structural sparsity while keeping the performance of the dense model.

Learning Sparse Filters for Interpretation. In DessiLBI, the structural sparsity parameter $\Gamma_{t}$ explores important sub-network architectures that contribute significantly to the loss or error reduction in early training stages. Through the $\ell_{2}$-coupling, structural sparsity parameter $\Gamma_{t}$ may guide the weight parameter to explore those sparse models in favour of improved interpretability. Figure 1 visualizes some sparse filters learned by DessiLBI of LeNet-5 trained on MNIST (with $\kappa=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. Further visualization of convolutional filters learned by ImageNet is shown in Figure 8 in Appendix, demonstrating the texture bias in ImageNet training [115].

7.2 Experiments on Network Sparsification

In this section, we conduct experiments on CIFAR10 dataset. Our DessiLBI is utilized for network sparsification. Emprically, with the augmented $\Gamma$, we can find sparse structure in the inverse scale space. In Appendix, Figure 6 and Figure 5, illustrate the distribution of filter norm for $W$ and $V$. It is clear that the magnitude scale of the augmented variable $V$ is balanced and is decoupled with arbitrary value of $W$.

To validate the efficacy of the structure, we use several experiments. By exploiting the magnitude scaling update strategy in computing $\Gamma$ for convolutional filters, and connections in convolutional layers and fully connected layers respectively, we introduce the filter and weight pruning in our experiments. For filter pruning, we aim to prune redundant convolutional filters, while the target of weight pruning is to prune connections in the network.

Common Setting. We follow the setting of LeGr [82], training the network for 200 epochs and fine-tuning the network for 400 epochs. For the training process, the learning rate is decayed by $5\times$ every 60 epochs with an initial learning rate 0.1. In post-processing, we compare fine-tuning and retraining. For fine-tuning, the learning rate is decayed by $5\times$ every 120 epochs. The initial learning rate for fine-tuning is 0.01. For retraining, the setting is the same as training. For all the processes, the weight decay is set as 5e-4.

Pruning Weights. Setting: We conduct experiments on CIFAR10 and select two widely used network structure VGGNet [119] and ResNet [114]. Several methods such Iterative-Pruning-A [116] as are selected for comparison. The training setting is the same as filter pruning. Here, the threshold proximal mapping is altered: 0.002 for ResNet50 and 0.0006 for VGGNet16.

Results. For ResNet50, it is relatively compact, so the extreme sparsity leads to a relative decrease in accuracy. However, our results are often better than our competitors. For VGGNet16, our method can find subnets with sparsity 0.06 and the accuracy is still comparable or even improved. Here fine-tuning is often better than retraining, often with an increase in terms of accuracy. The results are shown in Tab. II. Results for fine-tuning with only 10 epochs are also present, which shows that the pruned network can get decent performance with just a few epochs of fine-tuning. The keep ratio for each layer is visualized in Figure 9 in Appendix. For VGGNet16, most of the weight in the middle of the network can be pruned without hurting the performance. For the input conv layer and output fc layer, a high percent of weights are kept and the total number of weights for them is also smaller than other layers. For ResNet50, we can find an interesting phenomenon, that most layers inside the block can be pruned to a very sparse level. Meanwhile, for the input and output of a block, a relatively high percentage of weights should be kept.

Pruning Filters. Setting: We conduct experiments of filter pruning, and select several methods for comparison. For ResNet56, we pick up LEGR [82], HRank [118], AMC [120] and SFP [27] for comparison. And for VGGNet16, we choose L1 [121], HRank [118] ,Variational [122] and Hinge [123]. Here we set $\nu=500$ and $\kappa=1$. The batch size is 128.

Results. The results are shown in Tab. III. To control the whole sparsity, we set an upper bound for every layer which is 60% for both VGGNet16 and ResNet56. For both of the two networks, we can find that our pruning method achieves good results. For ResNet56, our model can reduce the FLOPs to 62% without hurting the performance. For fine-tuning setting with full (pre-)training, the pruned model can even increase the accuracy by about 0.5. Our early stopping experiments show that our method indeed finds good sparse structure at the early stage. By viewing the comparison between fine-tuning and retraining, we can find that for structure found by our method, fine-tuning slightly outperforms retraining. And both of the post-processing steps show a very decent performance, i.e., our method is robust to the post-processing. Similar results can be observed on VGGNet16 in Tab. III.

Pruning Layers. Setting: To push ahead the performance of our network sparsification, we alter the hyperparameters to get a further sparsification by using pruning layers as in Sec. 6.1. Here we use $\kappa=1$. We select $\nu=1200$ for VGGNet16 and $\nu=1500$ for ResNet56.

Results. The results are shown in Tab. III. For ResNet56, our method reduces the FLOP count to about $45\%$, while the performance is only dropped by $0.26\%$. The detailed structure is shown in Figure 2. Interestingly, the sparsified VGG16 network actually has improved performance over the original VGG16 significantly. Note that although the Flop counts of VGGNet 16 are similar for pruning filters and pruning layers, the network of pruning layers is much sparser than that of filter pruning, with $10.95\%$ of the parameters remained, in contrast to the sparsity for pruning filters at about 42.82%. For VGG16, most of the filters close to the input layer are selected by our network sparsification and much of the pruning occurs near the output layer. By viewing the structure of VGGNet, it is clear that redundancy exists in the layers close to the output layer. It is in accordance with our results. For ResNet56, we drop two layers in the middle of the corresponding blocks. The whole structure shows a dense selection in the beginning and end of channel alternating stage and a sparse selection inside each stage. A more detailed table is shown in Tab. VII in Appendix.

7.3 Experiments on Winning Ticket

In this section, with early stopping, $\Gamma_{t}$ may learn effective subnetworks (i.e. “winning tickets” in [31]) in the inverse scale space. After retraining, these subnetworks achieve comparable or even better performance than the dense network. Our method is comparable to existing pruning strategies with an improved efficiency.

Settings. We find winning tickets on CIFAR10 dataset and two network structures VGGNet16 and ResNet50. When DessiLBI is used for finding winning tickets, we set $\mu=500,\lambda=0.05$ for ResNet50 and $\mu=200,\lambda=0.03$ for VGG16 respectively. $\kappa=1$, $\alpha=0.1$ and training epoch $T=160$ is used in both case. Then, we retrain the winning tickets using SGD with 160 epochs from the end of the second epoch in the first stage. The initial learning rate for SGD is 0.1 and decayed by 10 at the 80 and 120 epoch. To make a fair comparison, in one shot pruning and iterative pruning, we also use SGD with this setting for training and retraining.

Winning tickets in early stopping. The results are shown in Figure 3. DessiLBI gives a sequence of winning tickets with increasing sparsity. All the winning tickets found by our method achieve similar or even better performance than the winning tickets found by one shot pruning and iterative pruning. All the winning tickets outperform the two baseline methods, especially in extremely sparse cases. It is worth noting that DessiLBI with early stopping at a very early epoch (for example, the end of second epoch) is computationally effective and gives extremely sparse winning tickets. Comparing with one shot pruning, we can avoid complete training in the first stage and save up to training 158 epochs. Since iterative pruning consists of 10 iterative pruning, We can save more computation by using early stopped DessiLBI followed by retraining once. The high efficiency of DessiLBI is due to exploring subnets via inverse scale space in which, important subnets come out first.

Transferability of winning tickets. We further explore the generalization of winning tickets found by DessiLBI cross different datasets like [33]. The results are shown in Figure 4. Specifically, we search winning tickets on source datasets and retrain with target datasets.

Similar to the observations in [33], our winning tickets discovered by DessiLBI in early stopping are transferable in the sense that they generalize across a variety of natural image datasets, such as CIFAR10, CIFAR100, and SVHN, often achieving performance close to the winning tickets generated on the same dataset. Moreover, winning tickets of DessiLBI may significantly outperform random tickets in very sparse cases. This shows that our algorithm can find transferable winning tickets.

We try to give some tentative explanations about the transferability. Despite that different datasets are used here, they all come from the same domain, i.e., natural images. In that sense, the images and models on these datasets shall share some common features; and winning tickets found on different datasets may encode the shared inductive bias and benefit the performance of sparse networks trained on another dataset.

7.4 Experiments on Network Growing

In this section, we grow a network to explore compact network structures via Inverse Scale Space. In detail, our method attempts to find a better filter configuration for each layer during training. To verify its efficacy, we conduct several experiments with CIFAR10 dataset.

Settings. We first try to verify our method with a plain network structure which means we do not use some special structure such as residue connection [114]. Here we utilize VGG16 [119] as the backbone. The backbone is initialized with eight filters in each convolutional layer and there are two fully connected layers following the setting of  [121]. The hyper-parameters of DessiLBI here are $\kappa$ = 1, $\nu$ = 100 and the learning rate is set as 0.01. And the batch size is 128. We set the hyper-parameters of our method as $J=10$, $\tau=0.8$. After finishing growing structure, we decrease the learning rate by $10\times$ and continue training for 30 epochs.

Results. Table. IV illustrates the filter configuration of models found by $L_{1}$-prune [121] and our method. We can observe that among the 13 convolutional layers, $L_{1}$-prune only reduces filters in the first convolution layer and the last 6 convolution layers by half, with the total number of final parameters being 5.40 M. By comparison, our method obtains a more compact model with only 2.96M parameters.

A crucial observation is that our growing method prefers to filters of large feature spatial size, which helps our model with a better performance. Such results suggest that our growing method may learn better model structures by learning better filter configurations. Specifically, pruning has the limitation that it only tries to reduce the number of filters in each layer based on the originally designed number, while our method gets rid of this limitation and adapts to better configurations.

Note that our growing method is different from network architecture search (NAS) [95, 98, 99]. NAS searches the overall structure of the new network in the predefined space and spends much more time and computational cost on this searching process. However, our method is much more light and can find good network structure during training. Besides, our model adds little extra memory cost and is more device-friendly.

To further validate our method, AutoGrow [124] is picked for comparison. ResNet [114] is selected as the backbone. For AutoGrow, the results are from the original paper. For our method, we use ResNet56 with initial channel number 8. The results are in Table V. It is clear that the networks grown by our method are more compact while having a better performance than AutoGrow. For more sufficient comparison, we try to compare our approach with several recent method, such as NASH [125], Splitting [126], Energy-aware Splitting [127], Firefly growing [128].Please refer to the Appendix. E for more experimental comparison on our growing algorithm.

8 Conclusion

In this paper, a parsimonious deep learning method is proposed based on the dynamics of inverse scale spaces. Its simple discretization – DessiLBI has a provable global convergence on DNNs and thus is employed to train deep networks. Our DessiLBI can explore the structural sparsity without sacrificing the performance of the dense, over-parameterized models, in favour of interpretability. It helps identify effective sparse structure at the early stage which is verified with several experiments on network sparsification and finding winning ticket subnets which are transferable across different natural image datasets. Particularly, our method can select the structure in the level of layers which is not well studied in the existing literature. Furthermore, the proposed method can be applied to grow a network by simultaneously learning both network structure and parameters.

Appendix to DessiLBI for deep learning: structural sparsity via differential inclusion paths

Appendix A Proof of Theorem 1

First of all, we reformulate Eq. (12) into an equivalent form. Without loss of generality, consider $\Omega=\Omega_{1}$ in the sequel.

Denote $R(P):=\Omega(\Gamma)$, then Eq. (12) can be rewritten as, DessiLBI

[TABLE]

where $p_{k}=[0,g_{k}]^{T}\in\partial R(P_{k})$ and $g_{k}\in\partial\Omega(\Gamma_{k})$. Thus DessiLBI is equivalent to the following iterations,

[TABLE]

Exploiting the equivalent reformulation (26a-26c), one can establish the global convergence of $(W_{k},\Gamma_{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 DessiLBI] Suppose that Assumption 1 holds. Let $\{(W_{k},\Gamma_{k},g_{k})\}$ be the sequence generated by DessiLBI (Eq. (26a-26c)) with a finite initialization. If

[TABLE]

then $\{(W_{k},\Gamma_{k},g_{k})\}$ converges to a critical point of $F$. Moreover, $\{(W_{k},\Gamma_{k})\}$ converges to a stationary point of $\bar{\mathcal{L}}$ defined in Eq. 8, and $\{W^{k}\}$ converges to a stationary point of $\widehat{\mathcal{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 [109].

The notion of subdifferential plays a central role in the following definitions. For each ${\bf x}\in\mathrm{dom}(h):=\{{\bf x}\in\mathbb{R}^{p}:h({\bf x})<+\infty\}$, the Fréchet subdifferential of $h$ at ${\bf x}$, written $\widehat{\partial}h({\bf x)}$, is the set of vectors ${\bf v}\in\mathbb{R}^{p}$ which satisfy

[TABLE]

When ${\bf x}\notin\mathrm{dom}(h),$ we set $\widehat{\partial}h({\bf x})=\varnothing.$ The limiting-subdifferential (or simply subdifferential) of $h$ introduced in [129], written $\partial h({\bf x})$ at ${\bf x}\in\mathrm{dom}(h)$, is defined by

[TABLE]

A necessary (but not sufficient) condition for ${\bf x}\in\mathbb{R}^{p}$ to be a minimizer of $h$ is $\mathbf{0}\in\partial h({\bf x})$. A point that satisfies this inclusion is called limiting-critical or simply critical. The distance between a point ${\bf x}$ to a subset ${\cal S}$ of $\mathbb{R}^{p}$, written $\mathrm{dist}({\bf x},{\cal S})$, is defined by $\mathrm{dist}({\bf x},{\cal S})=\inf\{\|{\bf x}-{\bf s}\|:{\bf s}\in{\cal S}\}$, where $\|\cdot\|$ represents the Euclidean norm.

Let $h:\mathbb{R}^{p}\to\mathbb{R}\cup\{+\infty\}$ be an extended-real-valued function (respectively, $h:\mathbb{R}^{p}\rightrightarrows\mathbb{R}^{q}$ be a point-to-set mapping), its graph is defined by

[TABLE]

and its domain by $\mathrm{dom}(h):=\{{\bf x}\in\mathbb{R}^{p}:h({\bf x})<+\infty\}$ (resp. $\mathrm{dom}(h):=\{{\bf x}\in\mathbb{R}^{p}:h({\bf x})\neq\varnothing\}$). When $h$ is a proper function, i.e., when $\mathrm{dom}(h)\neq\varnothing,$ the set of its global minimizers (possibly empty) is denoted by

[TABLE]

The KL property [106, 130, 131, 132, 133] plays a central role in the convergence analysis of nonconvex algorithms [108, 67]. The following definition is adopted from [133].

Definition 1**.**

[Kurdyka-Łojasiewicz property] A function $h$ is said to have the Kurdyka-Łojasiewicz (KL) property at $\bar{u}\in\mathrm{dom}(\partial h):=\{v\in\mathbb{R}^{n}|\partial h(v)\neq\emptyset\}$, if there exists a constant $\eta\in(0,\infty)$, a neighborhood ${\cal N}$ of $\bar{u}$ and a function $\phi:[0,\eta)\rightarrow\mathbb{R}_{+}$, which is a concave function that is continuous at [math] and satisfies $\phi(0)=0$, $\phi\in{\cal C}^{1}((0,\eta))$, i.e., $\phi$ is continuous differentiable on $(0,\eta)$, and $\phi^{\prime}(s)>0$ for all $s\in(0,\eta)$, such that for all $u\in{\cal N}\cap\{u\in\mathbb{R}^{n}|h(\bar{u})<h(u)<h(\bar{u})+\eta\}$, the following inequality holds

[TABLE]

If $h$ satisfies the KL property at each point of $\mathrm{dom}(\partial h)$, $h$ is called a KL function.

KL functions include real analytic functions, semialgebraic functions, tame functions defined in some o-minimal structures [131, 133], continuous subanalytic functions [132] 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\subset\mathbb{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)=\sum_{j=0}^{\infty}\alpha_{j}(x-u)^{j},$ for some $\{\alpha_{j}\}\subset\mathbb{R}$. The function is said to be real analytic on $V\subset U$ if it is real analytic at each $u\in V$ [134, Definition 1.1.5]. The real analytic function $f$ over $\mathbb{R}^{p}$ for some positive integer $p>1$ can be defined similarly.

According to [134], 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({\bf x)}$ on $\mathbb{R}^{p}$ is analytic by checking the analyticity of $g(t):=h({\bf x}+t{\bf y})$ for any ${\bf x},{\bf y}\in\mathbb{R}^{p}$.

Definition 3**.**

[Semialgebraic]

1. (a)

A set ${\cal D}\subset\mathbb{R}^{p}$ is called semialgebraic **[135]** if it can be represented as

[TABLE]

where $P_{ij},Q_{ij}$ are real polynomial functions for $1\leq i\leq s,1\leq j\leq t.$ 2. (b)

A function $h:\mathbb{R}^{p}\rightarrow\mathbb{R}\cup\{+\infty\}$ (resp. a point-to-set mapping $h:\mathbb{R}^{p}\rightrightarrows\mathbb{R}^{q}$) is called semialgebraic if its graph $\mathrm{Graph}(h)$ is semialgebraic.

According to [136, 135] and [137, 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 [135, 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,\ldots,l-1$. Let $d_{0}$ and $d_{l}$ be the number of units of input and output layers, respectively. Let $W^{i}\in\mathbb{R}^{d_{i}\times d_{i-1}}$ be the weight matrix between the $(i-1)$-th layer and the $i$-th layer for any $i=1,\ldots l$444To simplify notations, we regard the input and output layers as the [math]-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 (15) satisfies the Kurdyka-Łojasiewicz property. For this purpose, we need an extension of semialgebraic set, called the o-minimal structure (see, for instance [138], [139], [131], [133]). The following definition is from [133].

Definition 4**.**

[o-minimal structure] An o-minimal structure on $(\mathbb{R},+,\cdot)$ is a sequence of boolean algebras ${\cal O}_{n}$ of “definable” subsets of $\mathbb{R}^{n}$, such that for each $n\in\mathbb{N}$

1. (i)

if $A$ belongs to ${\cal O}_{n}$, then $A\times\mathbb{R}$ and $\mathbb{R}\times A$ belong to ${\cal O}_{n+1}$; 2. (ii)

if $\Pi:\mathbb{R}^{n+1}\rightarrow\mathbb{R}^{n}$ is the canonical projection onto $\mathbb{R}^{n}$, then for any $A$ in ${\cal O}_{n+1}$, the set $\Pi(A)$ belongs to ${\cal O}_{n}$; 3. (iii)

${\cal O}_{n}$* contains the family of algebraic subsets of $\mathbb{R}^{n}$, that is, every set of the form*

[TABLE]

where $p:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a polynomial function. 4. (iv)

the elements of ${\cal O}_{1}$ 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 ${\cal O}$ (over $(\mathbb{R},+,\cdot)$), a function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is said to be definable in ${\cal O}$ if its graph belongs to ${\cal O}_{n+1}$.

According to [140, 133], there are some important facts of the o-minimal structure, shown as follows.

1. (i)

The collection of semialgebraic sets is an o-minimal structure. Recall the semialgebraic sets are Bollean combinations of sets of the form

[TABLE]

where $p$ and $q_{i}$’s are polynomial functions in $\mathbb{R}^{n}$. 2. (ii)

There exists an o-minimal structure that contains the sets of the form

[TABLE]

where $f$ is real-analytic around $[-1,1]^{n}$. 3. (iii)

There exists an o-minimal structure that contains simultaneously the graph of the exponential function $\mathbb{R}\ni x\mapsto\exp(x)$ and all semialgebraic sets. 4. (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 [131, Theorem 1] and [133, Theorem 14], 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 (16), $F$ can be rewritten as follows

[TABLE]

Because $\ell$ and $\sigma_{i}$’s are definable by assumptions, then $\widehat{\mathcal{L}}_{n}(W,\Gamma)$ are definable as compositions of definable functions. Moreover, according to [134], $\|W-\Gamma\|^{2}$ and $\langle\Gamma,g\rangle$ are semi-algebraic and thus definable. Since the group Lasso $\Omega(\Gamma)=\sum_{g}\|\Gamma\|_{2}$ is the composition of $\ell_{2}$ and $\ell_{1}$ norms, and the conjugate of group Lasso penalty is the maximum of group $\ell_{2}$-norm, i.e. $\Omega^{*}(\Gamma)=\max_{g}\|\Gamma_{g}\|_{2}$, where the $\ell_{2}$, $\ell_{1}$, and $\ell_{\infty}$ 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 [131, Theorem 1].

The verifications of other cases listed in assumptions can be found in the proof of [105, Proposition 1]. This finishes the proof of this corollary. ∎

A.3 Proof of Theorem 2

Our analysis is mainly motivated by a recent paper [141], as well as the influential work [108]. According to Lemma 2.6 in [108], 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, (15)) 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 $\bar{\mathcal{L}}(W,\Gamma)$ and the relation $\lim_{k\rightarrow\infty}B_{\Omega}^{g_{k}}(\Gamma_{k+1},\Gamma_{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 [108, Lemma 2.6].

Let $(\bar{W},\bar{\Gamma},\bar{g})$ be a critical point of $F$, then the following holds

[TABLE]

By the final inclusion and the convexity of $\Omega$, it implies $\bar{g}\in\partial\Omega(\bar{\Gamma})$. Plugging this inclusion into the second inclusion yields $\alpha\nu^{-1}(\bar{\Gamma}-\bar{W})=0$. Together with the first equality imples

[TABLE]

This finishes the proof of this theorem.

A.4 Sufficient Descent Property along Lyapunov Function

Let $P_{k}:=(W_{k},\Gamma_{k})$, and $Q_{k}:=(P_{k},g_{k-1}),k\in\mathbb{N}$. In the following, we present the sufficient descent property of $Q_{k}$ along the Lyapunov function $F$.

Lemma. Suppose that $\widehat{\mathcal{L}}_{n}$ is continuously differentiable and $\nabla\widehat{\mathcal{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<\alpha<\frac{2}{\kappa(Lip+\nu^{-1})}$, then

[TABLE]

where $\rho:=\frac{1}{\kappa}-\frac{\alpha(Lip+\nu^{-1})}{2}$.

Proof.

By the optimality condition of (25a) and also the inclusion $p_{k}=[0,g_{k}]^{T}\in\partial R(P_{k})$, there holds

[TABLE]

which implies

[TABLE]

where

[TABLE]

Noting that $\bar{\mathcal{L}}(P)=\widehat{\mathcal{L}}_{n}(W)+\frac{1}{2\nu}\|W-\Gamma\|_{2}^{2}$ and by the Lipschitz continuity of $\nabla\widehat{\mathcal{L}}_{n}(W)$ with a constant $Lip>0$ implies $\nabla\bar{\mathcal{L}}$ is Lipschitz continuous with a constant $Lip+\nu^{-1}$. This implies

[TABLE]

Substituting the above inequality into (30) yields

[TABLE]

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

[TABLE]

where the final equality holds for $D(\Gamma_{k+1},\Gamma_{k})-B_{\Omega}^{g_{k}}(\Gamma_{k+1},\Gamma_{k})=B_{\Omega}^{g_{k+1}}(\Gamma_{k},\Gamma_{k-1}).$ That is,

[TABLE]

where the final inequality holds for $B_{\Omega}^{g_{k+1}}(\Gamma_{k},\Gamma_{k-1})\geq 0$ and $B_{\Omega}^{g_{k-1}}(\Gamma_{k},\Gamma_{k-1})\geq 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

1. (i)

both $\alpha\{\bar{\mathcal{L}}(P_{k})\}$ and $\{F(Q_{k})\}$ converge to the same finite value, and $\lim_{k\rightarrow\infty}B_{\Omega}^{g_{k}}(\Gamma_{k+1},\Gamma_{k})=0.$ 2. (ii)

the sequence $\{(W_{k},\Gamma_{k},g_{k})\}$ is bounded, 3. (iii)

$\lim_{k\rightarrow\infty}\|P_{k+1}-P_{k}\|_{2}^{2}=0$* and $\lim_{k\rightarrow\infty}D(\Gamma_{k+1},\Gamma_{k})=0,$* 4. (iv)

$\frac{1}{K}\sum_{k=0}^{K}\|P_{k+1}-P_{k}\|_{2}^{2}\rightarrow 0$* at a rate of ${\cal O}(1/K)$.*

Proof.

By (31), $\bar{\mathcal{L}}(P_{k})$ is monotonically decreasing due to $D(\Gamma_{k+1},\Gamma_{k})\geq 0$. Similarly, by (34), $F(Q^{k})$ is also monotonically decreasing. By the lower boundedness assumption of $\widehat{\mathcal{L}}_{n}(W)$, both $\bar{\mathcal{L}}(P)$ and $F(Q)$ are lower bounded by their definitions, i.e., (8) and (15), respectively. Therefore, both $\{\bar{\mathcal{L}}(P_{k})\}$ and $\{F(Q_{k})\}$ converge, and it is obvious that $\lim_{k\rightarrow\infty}F(Q_{k})\geq\lim_{k\rightarrow\infty}\alpha\bar{\mathcal{L}}(P_{k})$. By (33),

[TABLE]

By the convergence of $F(Q_{k})$ and the nonegativeness of $B_{\Omega}^{g_{k-1}}(\Gamma_{k},\Gamma_{k-1})$, there holds

[TABLE]

By the definition of $F(Q_{k})=\alpha\bar{\mathcal{L}}(P_{k})+B_{\Omega}^{g_{k-1}}(\Gamma_{k},\Gamma_{k-1})$ and the above equality, it yields

[TABLE]

Since $\widehat{\mathcal{L}}_{n}(W)$ has bounded level sets, then $W_{k}$ is bounded. By the definition of $\bar{\mathcal{L}}(W,\Gamma)$ and the finiteness of $\bar{\mathcal{L}}(W_{k},\Gamma_{k})$, $\Gamma_{k}$ is also bounded due to $W_{k}$ is bounded. The boundedness of $g_{k}$ is due to $g_{k}\in\partial\Omega(\Gamma_{k})$, condition (d), and the boundedness of $\Gamma_{k}$.

By (34), summing up (34) over $k=0,1,\ldots,K$ yields

[TABLE]

Letting $K\rightarrow\infty$ and noting that both $\|P_{k+1}-P_{k}\|^{2}$ and $D(\Gamma_{k+1},\Gamma_{k})$ are nonnegative, thus

[TABLE]

Again by (35),

[TABLE]

which implies $\frac{1}{K}\sum_{k=0}^{K}\|P_{k+1}-P_{k}\|^{2}\rightarrow 0$ at a rate of ${\cal 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$ (15),

[TABLE]

**Lemma. ** Under assumptions of Lemma 1, then

[TABLE]

where $\rho_{1}:=2\kappa^{-1}+1+\alpha(Lip+2\nu^{-1})$. Moreover, $\frac{1}{K}\sum_{k=1}^{K}\|H_{k}\|^{2}\rightarrow 0$ at a rate of ${\cal O}(1/K)$.

Proof.

Note that

[TABLE]

By the definition of $\bar{\mathcal{L}}$ (see (8)),

[TABLE]

where the last inequality holds for the Lipschitz continuity of $\nabla\widehat{\mathcal{L}}_{n}$ with a constant $Lip>0$, and by (26a),

[TABLE]

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

[TABLE]

Thus,

[TABLE]

By (26c), it yields

[TABLE]

Noting that $\nabla_{\Gamma}\bar{\mathcal{L}}(W_{k},\Gamma_{k})=\nu^{-1}(\Gamma_{k}-W_{k})$, and after some simplifications yields

[TABLE]

where the last inequality holds for the triangle inequality and $\kappa^{-1}>\alpha\nu^{-1}$ by the assumption.

By (41), (A.5), and the definition of $H_{k+1}$ (39), there holds

[TABLE]

By (43) and Lemma 1(iv), $\frac{1}{K}\sum_{k=1}^{K}\|H_{k}\|^{2}\rightarrow 0$ at a rate of ${\cal O}(1/K)$.

This finishes the proof of this lemma. ∎

Appendix B Supplementary Experiments

B.1 Ablation Study on Image 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,

SGD: (1) Naive SGD: the standard SGD with batch input. (2) SGD with $\mathit{l}_{1}$ penalty (Lasso). The $\mathit{l}_{1}$ norm is applied to penalize the weights of SGD by encouraging the sparsity of learned model, with the regularization parameter of the $\mathit{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 $\mathit{l}_{2}$ 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.

The results of image classification are shown in Tab.  VI . It shows the experimental results on ImageNet-2012, CIFAR10, and MNIST of some classical networks – LeNet, AlexNet and ResNet. Our DessiLBI variants may achieve comparable or even better performance than SGD variants in 100 epochs, indicating the efficacy in learning dense, over-parameterized models. The visualization of learned ResNet-18 on ImageNet-2012 is given in Figure 8.

B.2 Layer Selection

To push the performance of our network sparsification, we alter the hyperparameters to get a much higher selection by using pruning layers as in Sec. 6.1. The results are shown in Tab. VII. Here we use $\kappa=1$. We select $\nu=1200$ for VGGNet16 and $\nu=1500$ for ResNet56. For ResNet56, our method drops 2 layers with the correponding FLOP count reduced to 44.89%. By using pruning layers in VGG16, about 90% parameters can be removed. Interestingly, the sparsified VGG16 network actually has the improved performance over the original VGG16. The detailed structure is shown in Figure 2. For VGG16, most of the filters close to the input layer are selected by our network sparsification and much of the pruning exists near the output layer. By viewing the structure of VGGNet, it is clear that the redundancy existed in the layers close to the output layer. It is in accordance with our results. For ResNet56, we drop two layers in the middle of the corresponding block. The whole structure shows an intense selection in the beginning and ending of channel alternating stage and a sparse selection inside each stage.

Appendix C Computational Cost of DessiLBI

We further compare the computational cost of different optimizers: SGD (Mom), DessiLBI (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 VIII.

Appendix D Visualization

D.1 Visualization of Filters

Filters learned by ImageNet prefer to non-semantic texture rather than shape and color as shown in Figure 8. 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 [142] 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 [115]. How to enhance the semantic shape invariance learning, is arguably a key to improve the robustness of convolutional neural networks.

D.2 Visualization of Training and Testing Curves

In

To validate the convergence after our modification empirically, we give a visualization of the training loss curve and testing accuracy curve for ResNet56 [114] and VGG16 [119] in Figure 7. It is clear that adding the scaling strategy still holds the convergence property as shown in the training and testing curve.

D.3 Visualization of Augmented Variable

To validate the efficacy of our modification , we give a visualization of the filter distribution across different layers for ResNet56 [114] in Figure 5 and VGG16 [119] in Figure 6. It can be found that after using the scaling strategy, the filter norms own more stable scale, which can be beneficial for finding the important structure.

D.4 Visualization of Selected Weight

The keep ratio for each layer is visualized in Figure 9. For VGGNet16, most of the weight in the middle of the network can be pruned without hurting the performance. For the input conv layer and output fc layer, a high percent of weights are kept and the whole number of weights for them is also smaller than other layers. For ResNet50, we can find an interesting phenomenon, most layers inside the block can be pruned to a very sparse level. Meanwhile, for the input and output of a block, a relatively high percent of weights should be kept.

Appendix E Comparison with Other Growing Methods

E.1 Performance and Computational Cost Comparison

We compare our growing algorithm against several recent methods. Particularly, we compare NASH [125], Splitting [126], Energy-aware Splitting [127], Firefly growing [128]. To give a complete comparison, we take the both the parameter of the grew model and the time to grow the network into consideration. We conduct experiments on CIFAR10 [144]. VGG19 [119] is selected as the backbone. The initial seed network are set the same as the splitting based growing methods [126, 127, 128]. For VGG19, the initial network is [1, 1, ’M’, 2, 2, ’M’, 4, 4, 4, 4, ’M’, 8, 8, 8, 8, ’M’, 8, 8, 8, 8, ’M’]. Here ’M’ means the max pooling layer. We use the value 2 to denote 2 convolutional kernels for that layer. For these algorithms, the models are grow for 6 times sequentially. The methods of NASH, Splitting, Energy-aware Splitting, Firefly growing are trained on a single NVIDIA 3090ti GPU, as these methods demands relatively large GPU memory. In contrast, our growing algorithm is trained by a single NVIDIA 1080ti GPU.

As shown in Figure. 10, our growing method outperforms the splitting methods significantly along the growing path. In this figure full model means training the VGG19 model with full size. EA-Splitting and Firefly denotes Energy-aware Splitting [127] and Firefly growing [128] respectively. All the results except for our grow algorithm are obtained by run codes from the official open sourced github repo of Firefly growing [128]. Our algorithm is outperformed by EA-splitting at the initial process, but our method has significantly better performance for the following growing process.

For more complete comparison, we illustrate the detailed results and the running time in Table. IX. In this table, we present the running time to get the model which is shown in the unit of hours, the number of parameters of the model in millions and the test accuracy of the model. For these methods, we present the results of growing for 6 times. We can conclude that our grow algorithm can have better performance in this group of comparison experiments with significantly less training cost. In terms of running time, our method use only 5.31h to get a final model with 2.00M parameters and 0.9403 accuracy which outperforms others model obtained by running for more than 9 hours.

E.2 Detailed Model Structure

To further compare our methods with the selected methods, we present the detailed model structure in this section. Here for VGG19, the original version has 16 convolutional layers and 3 fully connected layers. For the growing experiments, we only consider the convolutional layers. For convenience, we only show the comparison between Firefly and our method for the last 2 growing. As shown in Table. X, we can find some significant difference. First of all, our method tends to add more filters to the middle of the network. Besides, our method adds only a few filters to the layers close to the output. For Firefly [128], it prefers more balanced structure compared with ours.

Appendix F Interpretation of $\alpha$ and $\kappa$

In practice, we donot apply the grid search for $\alpha$, and $\kappa$, while the defaulted configuration for them would to make $\alpha$ small enough, and $\kappa$ large enough. Since parameter $\alpha$ is the step size in Euler discretization and $\kappa$ is the damping factor, the basic principle to configure these parameters is, to make $\alpha$ small enough and $\kappa$ large but not violating $\alpha\cdot\kappa<2/(Lip+\nu-1)$, rather than the grid search for them. Particularly, we highlight several points about the hyperparameter $\kappa$ and $\alpha$.

1. damping factor $\kappa>0$ is to make the path continuous while $\kappa\to\infty$ is to approximate the ISS dynamics exponentially; 2) step size $\alpha>0$ is small enough such that $\alpha\cdot\kappa<2/(Lip+\nu-1)$ with convergence guarantee.

To see these points, consider the problem $y=X\beta^{*}+\epsilon$, $\gamma^{*}=D\beta^{*}$. Here $\gamma^{*}$ is sparse and $D$ denotes the transformation. The purpose is to estimate both $\beta^{*}$ and $\gamma^{*}$. Split LBI [39] utilizes variable splitting to combine $L_{2}$-Boost of model parameter and LBI of split structure parameter to find structural sparsity. Specifically, Split LBI [39] derives the following updating equations,

[TABLE]

where the loss function is,

[TABLE]

Here the factor $\nu$ controls the relaxation of $\gamma$ in the neighborhood of $D\beta$ and $prox_{\|\cdot\|_{1}}(\cdot)$ denotes the proximal mapping with $\ell_{1}$-norm. We can rewritten the equation, such that it can be viewed as a discretization of a differential inclusion.

[TABLE]

Here $\alpha$ is the step size. And let step size $\alpha\to 0$, we can get the following dynamics,

[TABLE]

Here $\kappa$ serves as the damping factor. When we let $\kappa\to\infty$, We can get the following dynamics,

[TABLE]

It is called Split ISS [39]. So $\kappa$ is to approximate the Split ISS dynamics exponentially. One can see that $\alpha$ and $\kappa$ have their meanings and we have to set their values accordingly. First, as to the damping factor $\kappa$, we can increase $\kappa$ to approximate ISS dynamics which is piece-wise constant. With larger $\kappa$, the regularization path is expected to be closer to the path of ISS, as shown in Figure. 11 using the simulation example from [76]. Second, the step size $\alpha$ should be small enough. The product of the damping factor $\kappa$ and step size $\alpha$ controls the update size for the weights $\beta$. And they are not totally free. We have to control $\alpha\cdot\kappa$ in an appropriate range ($\alpha\cdot\kappa<2/(Lip+\nu-1)$) for convergence. In Figure. 12, we show an example of fixing $\kappa=1$ using the simulation dataset from [76]. With a small $\alpha=0.05$, we can get a stable path, while setting a large $\alpha=0.515$ causes the divergence of training loss and the oscillation of path. In our practice, we fix $\kappa$ and try to find suitably small $\alpha$ when using Split LBI. When extending Split LBI to train deep networks, we prefer to set $\kappa=1$ and find a suitable $\alpha$ so that the model can be trained in a stable way. We have updated the paper, and include the explanation above into the Appendix of our paper.

Finally we make a note on ISS history. For the Inverse Scale Space (ISS) method [36], it was firstly used for image restoration. And the name comes from the observation that large-scale features are recovered faster than small-scale feature. In  [76], this kind of dynamics is imported to recover sparse signals from noisy measurements.  [76] utilizes differential inclusion to formalize two dynamics, Bregman ISS and Linearized Bregman ISS. And the Linearized Bregman ISS is the damping version of Bregman ISS with damping factor $\kappa$. When $\kappa\to\infty$, it is reduced to Bregman ISS. The path of Linearized Bregman ISS can approximate the ISS dynamic exponentially as $\kappa$ increases. And Linearized Bregman Iteration(LBI) is the discretization of Linearized Bregman ISS with step size $\alpha$.