{'fraction_non_alphanumeric': 0.04752316857899382, 'fraction_numerical': 0.018948587819947044, 'mean_word_length': 4.432574168414744, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Parsimonious Deep Learning : A Differential Inclusion Approach with Global Convergence

Over-parameterization is ubiquitous nowadays in training neural networks to benefit both optimization in seeking global optima and generalization in reducing prediction error. However, compressive networks are desired in many real world applications and direct training of small networks may be trapped in local optima. In this paper, instead of pruning or distilling an over-parameterized model to compressive ones, we propose a parsimonious learning approach based on differential inclusions of inverse scale spaces, that generates a family of models from simple to complex ones with a better efficiency and interpretability than stochastic gradient descent in exploring the model space. It enjoys a simple discretization, the Split Linearized Bregman Iterations, with provable global convergence that from any initializations, algorithmic iterations converge to a critical point of empirical risks. One may exploit the proposed method to boost the complexity of neural networks progressively. Numerical experiments with MNIST, Cifar-10/100, and ImageNet are conducted to show the method is promising in training large scale models with a favorite interpretability.


Introduction
The expressive power of Deep Convolutional Neural Networks (DNNs) comes from the millions of parameters, which are optimized by Stochastic Gradient Descent (SGD) algorithms [Bottou, 2010] and variants like Adam [Kingma and Ba, 2015]. Remarkably, model over-parameterization helps both optimization and generalization. For optimization, overparameterization may simplify the landscape of empirical risks toward locating global optima efficiently by gradient descent method [Mei et al., 2018, Venturi et al., 2018, Allen-Zhu et al., 2018, Du et al., 2018. On the other hand, over-parameterization does not necessarily result in a bad generalization or overfitting [Zhang et al., 2017], especially when some weightsize dependent complexities are controlled [Bartlett, 1997, Bartlett et al., 2017, Golowich et al., 2018, Neyshabur et al., 2019.
However, compressive networks are desired in many real world applications, e.g. robotics, self-driving cars, and augmented reality. Despite that 1 regularization has been applied to deep learning to enforce the sparsity on the weights toward compact, memory efficient networks, it sacrifices some prediction performance [Collins and Kohli, 2014]. This is because that the weights learned in neural networks are highly correlated, and 1 regularization on such weights violates the incoherence or irrepresentable conditions needed for sparse model selection [Donoho and Huo, 2001, Tropp, 2004, Zhao and Yu, 2006, leading to spurious selections with poor generalization. On the other hand, 2 regularization is often utilized for correlated weights as some low-pass filtering, sometimes in the form of weight decay [Loshchilov and Hutter, 2019] or early stopping [Yao et al., 2007, Wei et al., 2017. Furthermore, group sparsity regularization [Yuan and Lin, 2006] has also been applied to neural networks, such as finding optimal number of neuron groups [Alvarez and Salzmann, 2016] and exerting good data locality with structured sparsity [Wen et al., 2016, Yoon andHwang, 2017].
Yet, without the aid of over-parameterization, directly training a compressive model architecture may meet the obstacle of being trapped in local optima in contemporary experience. Alternatively, researchers in practice typically start from training a big model using common task datasets like ImageNet, and then prune or distill such big models to small ones without sacrificing too much of the performance [Jaderberg et al., 2014, Han et al., 2015, Zhu et al., 2017, Zhou et al., 2017, Li et al., 2017, Abbasi-Asl and Yu, 2017, Yang et al., 2018, Arora et al., 2018. In particular, a recent study [Frankle and Carbin, 2019] created the lottery ticket hypothesis based on empirical observations: "dense, randomly-initialized, feed-forward networks contain subnetworks (winning tickets) that -when trained in isolationreach test accuracy comparable to the original network in a similar number of iterations". How to effectively reduce an over-parameterized model thus becomes the key to compressive deep learning. Such a backward pruning approach is in a sharp contrast to traditional parsimonious or sparse learning approach where one starts from simple yet interpretable models, delving into complex models progressively. Examples of such a parsimonious learning include boosting as functional gradient descent or coordinate descent method, regularization paths associated with 2 (Tikhonov regularization or Ridge regression) and/or 1 (Lasso) penalties, etc. None of these has been systematically studied in deep learning.
In this paper, we aim to fill in this gap by exploiting a dynamic approach to forward parsimonious learning. We are able to establish a family of neural networks, from simple to complex, by following regularization paths as solutions of differential inclusions of inverse scale spaces. Our key idea is to design some dynamics that simultaneously exploit overparameterized models and structural sparsities. To achieve this goal, the original network parameters are lifted to a coupled pair, with one set of parameters following the standard gradient descend to explore the over-parameterized model space, while the other set of parameters learning structure sparsity in an inverse scale space, that the large-scale important parameters are learned at a fast speed while the small unimportant ones are learned at a slow speed. The two sets of parameters are coupled in an 2 regularization. The dynamics enjoys a simple discretization, the Split Linearized Bregman Iteration, with provable global convergence guarantee. Equipped with such a method, one may efficiently boost or grow a network with different sparse models, in a more efficient way than SGD to explore the model space and in favor of better interpretability. For example, Figure 1 exhibits the performance on training ResNet-18 on ImageNet, using SGD, Adam, and our algorithms, where a proper  Adam, and SplitLBI with κ = 1, 10, respectively. For the choice of κ = 1, SplitLBI enjoys more efficient learning curves than SGD and Adam noticeably in early epochs. Details about the experiments can be found in Section 5. choice of parameters leads to more efficient training than SGD and Adam in early epochs with a comparable accuracy in the end. Moreover, Figure 3 further demonstrates that the sparse convolutional filters trained on ImageNet, focus more on non-semantic textures than shapes and colors, in support of alternative studies [Abbasi-Asl andYu, 2017, Geirhos et al., 2019]. How to enhance the semantic shape invariance learning, is arguably a key to improve the robustness of convolutional neural networks.
The organization of the paper is as follows. Section 2 introduces the methodology; Section 3 establishes the global convergence theory; Section 4 discusses stochastic approximation algorithms and boosting; experimental results are presented in Section 5. Proofs are collected in Appendix A.

Methodology
The supervised learning task learns a mapping Φ W : X → Y, from input space X to output space Y, with a parameter W such as weights in neural networks, by minimizing certain loss functions on training samples L n (W ) = 1 n n i=1 (y i , Φ W (x i )). For example, a neural network of l-layer is defined as , σ i is the nonlinear activation function of the i-th layer and the loss function is usually the cross-entropy.
Differential Inclusion of Inverse Scale Space. Consider the following dynamics, whereΩ(Γ) = Ω(Γ) + 1 2κ Γ 2 for some sparsity-enforced regularization Ω(Γ) such as Lasso or grouped Lasso penalties, κ > 0 as a damping parameter such that the solution path is continuous, and the augmented loss function is with ν > 0 controlling the gap admitted between W and Γ. Compared to the original loss function L n (W ), theL (W, Γ) additionally adopt the variable splitting strategy, by lifting the original neural network parameter W to (W, Γ) with Γ to model the structural sparsity of W . For simplicity, we assumedL is differentiable with respect to W here, otherwise the gradient in (1a) is understood as subgradient and the equation becomes an inclusion. The differential inclusion system (1) is called as the Split Linearized Bregman Inverse Scale Space (SplitISS) for the following reasons. It can be understood as a gradient descent flow of W t in the proximity of Γ t and a mirror descent [Nemirovski and Yudin, 1983] flow of Γ t associated with a sparsity enforcement penaltyΩ. For a large enough ν, it reduces to the gradient descent method for W t . Yet the solution path of Γ t exhibits the following selection property in the separation of scales: starting at the zero, important parameters of large scale will be learned fast, popping up to be nonzeros early, while unimportant parameters of small scale will be learned slowly, appearing to be nonzeros late. In fact, taking Ω(Γ) = Γ 1 and κ → ∞ for simplicity, V t undergoes a gradient descent flow before reaching the ∞ -unit box, which implies that Γ t = 0 in this stage and W t follows the gradient descent with a standard 2 -regularization; the earlier a component in V t reaching the ∞ -unit box, the earlier a corresponding component in Γ t becomes nonzero and rapidly evolves toward a critical point ofL under gradient flow, which drives W t closely following dynamics of Γ t whose important parameters are selected. Such a property is called as the inverse scale space in applied mathematics [Burger et al., 2006] and recently was shown to achieve statistical model selection consistency in high dimensional linear regression [Osher et al., 2016] and general linear models , with a reduction of bias as κ increases. In this paper, we shall see the inverse scale space property still holds for neural networks trained by Eq. (1), e.g. Figure 2 shows a LeNet trained on MNIST by the discretized dynamics, where important sparse filters are selected in early epochs while the popular stochastic gradient descent method returns dense filters.
Compared with directly enforcing a penalty function Ω (W ) (e.g. 1 or 2 regularization) the SplitISS avoids the parameter correlation problem in over-parameterized models. In fact, a necessary and sufficient condition for Lasso or 1 -type sparse model selection is the incoherence or irrepresentable conditions ( [Tropp, 2004, Zhao andYu, 2006]) that are violated for highly correlated weight parameters, leading to spurious discoveries. In contrast, Huang et al.  showed that equipped with such a variable splitting where Γ enjoys an orthogonal design, the SplitISS can achieve model selection consistency under weaker conditions than generalized Lasso, relaxing the incoherence or irrepresentable conditions when parameters are highly correlated. For weight parameter W , instead of directly being imposed with 1 -sparsity, it adopts 2 -regularization in the proximity of the sparse path of Γ that admits simultaneously exploring highly correlated parameters in over-parameterized models and sparsity regularizations. x-axis and y-axis indicate the training epochs, and filter magnitudes ( 2 -norm), respectively. The SplitLBI path of filters selected in the support of Γ are drawn in blue color, while the red color curves represent the filters that are not important and outside the support of Γ. We visualize the corresponding learned filters by Erhan et al. [2009] at 20 (blue), 40 (green), and 80 (black) epochs, which are shown in the right figure with the corresponding color bounding boxes, i.e., blue, green, and black, respectively. It shows that our SplitLBI enjoys a sparse selection of filters without sacrificing accuracy (see Table 1).
Split Linearized Bregman Iterations. The SplitISS admits an extremely simple discrete approximation, using the Euler forward discretization of the dynamics above: where V 0 = Γ 0 = 0, W 0 can be small random numbers such as Gaussian distribution in neural networks, the proximal map in Eq (4c) that controls the sparsity of Γ is given by We shall call such an iterative procedure as Split Linearized Bregman Iteration (SplitLBI), that was firstly proposed in [Huang et al., 2016] as an iterative regularization path for sparse modeling in high dimensional statistics with an improved model selection consistency than generalized Lasso, and has been later used in medical image analysis [Sun et al., 2017] and computer vision [Zhao et al., 2018]. In the application to neural networks, the loss becomes highly non-convex, the SplitLBI returns a sequence of sparse models from simple to complex ones with a global convergence guarantee shown below, while solving Eq (3) at various levels of λ might not be tractable except for over-parameterized models. The sparsity-enforcement penalty Ω(·) used in convolutional neural networks can be chosen as follows. Our choice in this paper aims at regularizing the groups of weight parameters using the group Lasso penalty [Yuan and Lin, 2006] defined by, where |Γ g | is the number of weights in Γ g . We treat convolutional layers and fully connected layers in different ways as follows.
(a) For a convolutional layer, Γ g = Γ g (c in , c out , size) denote the convolutional filters where size denotes the kernel size and c in and c out denote the numbers of input channels and output channels, respectively. We select g = c out and regard each group as each convolutional filter.
(b) For a fully connected layer, Γ = Γ(c in , c out ) where c in and c out denote the numbers of inputs and outputs of the fully connected layer. Each group g corresponds to each element (i, j), and the group Lasso penalty degenerates to the Lasso penalty.

Global Convergence of SplitLBI for Neural Networks
We present a theorem that guarantees the global convergence of SplitLBI, i.e. from any intialization, the SplitLBI sequence converges to a critical point ofL.
Assumption 1. Suppose that: (a) L n (W ) is continuous differentiable and ∇ L n is Lipschitz continuous with a positive constant Lip; (b) L n (W ) has bounded level sets; (c) L n (W ) is lower bounded (without loss of generality, we assume that the lower bound is 0); (d) Ω is a proper lower semi-continuous convex function and has locally bounded subgradients, that is, for every compact set S ⊂ R n , there exists a constant C > 0 such that for all Γ ∈ S and all g ∈ ∂Ω(Γ), there holds g ≤ C; (e) the Lyapunov function is a Kurdyka-Lojasiewicz function on any bounded set, where and Ω * is the conjugate of Ω defined as Remark 1. Assumption 1 (a)-(c) are regular in the analysis of nonconvex algorithm (see, Attouch et al. [2013] for instance), while Assumption 1 (d) is also mild including all Lipschitz continuous convex function over a compact set. Some typical examples satisfying Assumption 1(d) are the 1 norm, group 1 norm, and every continuously differentiable penalties. By Eq (11) and the definition of conjugate, the Lyapunov function F can be rewritten as follows, Now we are ready to present the main theorem.
Theorem 1 (Global Convergence of SplitLBI). Suppose that Assumption 1 holds. Let (W k , Γ k , g k ) be the sequence generated by SplitLBI (Eq (10a)-(10c)) with a finite initialization. If , converges to a stationary point ofL defined in Eq (2), and {W k } converges to a stationary point of L n (W ).
Applying to the neural network training settings, typical examples satisfying the assumption of this theorem are summarized in the following corollary. Corollary 1. Let {W k , Γ k , g k } be a sequence generated by SLBI (10a)-(10c) for the deep neural network training problem where (a) is any smooth definable loss function, such as the square loss t 2 , exponential loss e t , logistic loss log(1 + e −t ), and cross-entropy loss; (b) σ i is any smooth definable activation, such as linear activation, sigmoid 1 1+e −t , hyperbolic tangent e t −e −t e t +e −t , and softplus ( 1 c log(1 + e ct ) for some c > 0) as a smooth approximation of ReLU; (c) Ω is the group Lasso penalty (i.e., Ω(Γ) = g Γ g 2 ).
Moreover, suppose that the empirical risk L n (W ) has lower bounded level set and its gradient is Lipschitz continuous with a constant Lip > 0. Then the sequence {W k } converges to a stationary point of Proofs of Theorem 1 and Corollary 1 are given in Appendix A.

Stochastic Approximation of SplitLBI and Applications
For neural network training with large datasets, stochastic approximation of the gradients in Split LBI over the mini-batch (X, Y) batcht is adopted to update the parameter W , Moreover, inspired by the variants of SGD, the momentum term can be also incorporated to the standard Split LBI that leads to the following updates of W by replacing Eq ((4a)) with, where τ is the momentum factor, empirically setting as 0.9 in default. One immediate application of such stochastic algorithms of SplitLBI is to "boost networks", i.e. growing a network from the null to a complex one by sequentially applying our algorithm on subnets with increasing complexities.
Boosting Networks. We propose the task of boosting network by SplitLBI, i.e., learning to expand the convolutional filters and network parameters during training. Remarkably, our boosting network is very different from previous tasks, including AutoML [Wong et al., 2018], or life-long learning , Thrun and Mitchell, 1995, Pentina and Lampert, 2015, Li and Hoiem, 2016, knowledge distill [Hinton et al., 2014]. Starting from very few filters of each convolutional layer, algorithms of boosting network require 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. Note that the existing work mentioned above does not allow additional expanding or fine-tuning algorithms which are very computational expensive in practice.
To apply the SplitLBI to boost networks, we define a projection of the convolutional filters W on the support set of Γ, The basic idea is to monitor the gap between W and its projection W along the iterations: when the gap becomes small, we are going to expand the network by adding new filters. We give the description of the procedure of expanding a network by SplitLBI. Starting from a small number of filters (e.g. 2) of convolutional layer, more and more filters tend to be with non-zero value as the algorithm iterates. At each step, we can compute the ratio of W k to be non-zero. If this ratio passes a pre-set threshold (empirically setting as 0.8), we add m (e.g. 2) new convolutional filters into W ; and the newly added parameters are randomly initialized. Then we continue learning the network from training data. This process is repeated until model converged or maximum epochs reached.

Experiments
The proposed SplitLBI is evaluated on MNIST, Cifar-10/100 and ImageNet-2012 datasets, which are widely utilized as the benchmarks in deep learning. Particularly, Sec. 5.1 compares SplitLBI against SGD and Adam in image classification. Furthermore, Sec. 5.2 shows that SplitLBI can efficiently boost a network from null to complex ones from data.

Image Classification
In the supervised setting, SGD, Adam and SplitLBI are adopted in optimizing DNNs on these datasets. The results are shown in Table 1. To train ResNet-18 on ImageNet-2012 dataset, the model trained by SplitLBI (κ = 1) can achieve the Top-1 accuracy of 67.76 by 40 epochs, which is 2% lower than the best ResNet-18 results report in  optimized by 120 epochs. Still this indicates that the SplitLBI may have a more efficient training than SGD. Note that, (1) In Fig. 1, we qualitatively compare our SplitLBI of κ = 1, 10 against SGD and Adam. Remarkably, at the early stage, our SplitLBI of κ = 1 has a much better efficiency than Adam and SGD within the same number of training epochs; (2) our SLBI requires the ResNet-18@ImageNet (Conv1) extra 7% − 8% more GPU usage than that of SGD, in order to save the parameters of V k and Γ k . Comparing with SGD, the SplitLBI adds little extra memory cost to store V and Γ. The total running time of each epoch is comparable to that of SGD.
To interpret the learned filters by SplitLBI with structural sparsity, we visualize a LeNet-5 trained on MNIST and a ResNet-18 on ImageNet.
Visualization of LeNet-5 on MNIST. We visualize the solution paths and filter patterns of the third convolutional layer (i.e., conv.c5) of LeNet-5 trained by SplitLBI and SGD trained on MNIST as in Fig. 2. X-axis and Y-axis indicate the training epochs, and magnitudes of each filter respectively. In particular, the left two figures compares the magnitude value changes of each filter trained by SplitLBI and SGD, individually. The filters that are not in the support of Γ, are drawn by red color. In the right figure, the model parameters learned at 20 (blue), 40 (green), and 80 (black) epochs are visualized as Erhan et al. [2009], with the corresponding color bounding boxes, i.e., blue, green, and black, respectively. The conv.c5 has 120 convolutional filters, which are ordered by the magnitude of filters. We can observe that the filters of high magnitude learned by our SplitLBI are much sparser and have very much clearer patterns that those trained by SGD, showing that SplitLBI encourages structural sparsity of filters.
Visualization of ResNet-18 on ImageNet-2012. To further reveal the insights of learned patterns of SplitLBI, we visualize the first convolutional layer of ResNet-18 on ImageNet-2012 in Fig. 3. The left figure compares the training and validation accuracy of SplitLBI (κ = 1, 10) and SGD. The right figure compares visualizations of the filters learned by SplitLBI and SGD using Springenberg et al. [2014]. To be specific, denote the weights of an l-layer network as {W 1 , W 2 , · · · , W l }. For the i−th layer weights W i , denote the j−th channel W i j . Then we compute the gradient of the sum of the feature map computed from each filter W i j with respect to the input image (here a snake image). We further conduct the min-max normalization to  Table 2: Results of boosting networks on MNIST, Cifar-10, and Cifar-100. All models are trained by SplitLBI in 1000 epochs (κ = 1, ν = 1000).
the gradient image, and generate the final visualization map. The right figure compares the visualized gradient images of first convolutional layer of 64 filters with 7×7 receptive fields. We visualize the models parameters at 20 (purple), 40 (green), and 60 (black) epochs, respectively, which corresponds to the bounding boxes in the right figure annotated by the corresponding colors, i.e., purple, green, and black. We order the gradient images produced from 64 filters by the descending order of the magnitude ( 2 -norm) of filters, i.e., images are ordered from the upper left to the bottom right. For comparison, we also provide the visualized gradient from random initialized weights. We find that the filters with high magnitudes mostly focus on the texture and shape information, while color information is kept in the filters with small magnitudes. This phenomenon is in accordance with observation of Abbasi-Asl and Yu [2017] that filters mainly contain color information can be pruned for saving computational cost. Moreover, among the filters of high magnitudes, most of them capture non-semantic textures while few pursue shapes. This shows that the first convolutional layer of ResNet-18 trained on ImageNet learned non-semantic textures rather than shape to do image classification tasks, in accordance with recent studies [Geirhos et al., 2019].

Boosting a Network
In this experiment, we explore the performance that SplitLBI boosts filters of a network with a fixed number of layers. In each individual layer, we start from a small initial number of filters (Min), and gradually grow to a maximum number of filters (Max), adaptively learned by the SplitLBI-Boosting described in Sec. 4. We compare the results of directly using small number of filters, maximum number of filters (Max) and our boosting method in Table 2. Note that, the maximum number of filters are computed by our boosting algorithm. We do not know it beforehand. In addition, we also include the networks with over-parameterized (O.P.) convolutional filters, which is twice the maximum number of filters of each corresponding layer. By default, all models are trained to 1000 epochs. We employ MNIST, Cifar-10 and Cifar-100 as the evaluation benchmarks here. As shown in Table 2, our SplitLBI can successfully boost the filters of networks from Min to Max. In Fig. 4, Fig. 5, Fig. 6, and 7, we give the full training paths of Min (Minimum Model), Max (Maximum Model) and Boosting (Grow) as compared in Table 2. Each black line indicates filters are added. In our experiments, two filters (m = 2) are added each time. To the best of our knowledge, this is the first algorithm that can simultaneously learn the network structures and parameters from data. It is desired that the Max and Boosting hit the same performance with enough training epochs; and impressively, after the 1000 epochs, Max and Boosting can indeed achieve almost comparable results on all of the two layer networks. This further indicates the efficacy of proposed SplitLBI. On Cifar-100 with a three convolutional layer network, the results of boosting are inferior to that of Max. A further inspection on the full training paths of Min, Max and Boosting (Fig. 7) shows that the boosting may require more training epochs in such a relative complex structures. Thus empirically, we fix the Boosted network in Table 2, continue training another 500 epochs, and we find the trained model indeed can report almost the same performance as the Max.

Conclusion
In this paper, a parsimonious deep learning method is proposed based on differential inclusions of inverse scale spaces. Its simple discretization -SplitLBI has a provable global convergence on DNNs, and thus is employed to train deep networks. Furthermore, the proposed method can be applied to boost a network by simultaneously learning both network structure and parameters. The extensive experiments and visualization show that the proposed method may not only have better efficiency than SGD, but also capture the structural sparsity of convolutional filters for interpretability.

Kurdyka-Lojasiewicz Property
To introduce the definition of the Kurdyka-Lojasiewicz (KL) property, we need some notions and notations from variational analysis, which can be found in Rockafellar and Wets [1998].
The notion of subdifferential plays a central role in the following definitions. For each x ∈ dom(h) := {x ∈ R p : h(x) < +∞}, the Fréchet subdifferential of h at x, written ∂h(x), is the set of vectors v ∈ R p which satisfy lim inf When x / ∈ dom(h), we set ∂h(x) = ∅. The limiting-subdifferential (or simply subdifferential ) of h introduced in Mordukhovich [2006], written ∂h(x) at x ∈ dom(h), is defined by A necessary (but not sufficient) condition for x ∈ R p to be a minimizer of h is 0 ∈ ∂h(x). A point that satisfies this inclusion is called limiting-critical or simply critical. The distance between a point x to a subset S of R p , written dist(x, S), is defined by dist(x, S) = inf{ x−s : s ∈ S}, where · represents the Euclidean norm.
Let h : R p → R ∪ {+∞} be an extended-real-valued function (respectively, h : R p ⇒ R q be a point-to-set mapping), its graph is defined by and its domain by dom(h) := {x ∈ R p : h(x) < +∞} (resp. dom(h) := {x ∈ R p : h(x) = ∅}). When h is a proper function, i.e., when dom(h) = ∅, the set of its global minimizers (possibly empty) is denoted by arg min h := {x ∈ R p : h(x) = inf h}.
KL functions include real analytic functions, semialgebraic functions, tame functions defined in some o-minimal structures [Kurdyka, 1998, Bolte et al., 2007b, continuous subanalytic functions [Bolte et al., 2007a] and locally strongly convex functions. In the following, we provide some important examples that satisfy the Kurdyka-Lojasiewicz property.
Definition 2 (Real analytic). A function h with domain an open set U ⊂ R and range the set of either all real or complex numbers, is said to be real analytic at u if the function h may be represented by a convergent power series on some interval of positive radius centered at u: h(x) = ∞ j=0 α j (x − u) j , for some {α j } ⊂ R. The function is said to be real analytic on V ⊂ U if it is real analytic at each u ∈ V [Krantz and Parks, 2002, Definition 1.1.5]. The real analytic function f over R p for some positive integer p > 1 can be defined similarly.
According to Krantz and Parks [2002], typical real analytic functions include polynomials, exponential functions, and the logarithm, trigonometric and power functions on any open set of their domains. One can verify whether a multivariable real function h(x) on R p is analytic by checking the analyticity of g(t) := h(x + ty) for any x, y ∈ R p .

Definition 3 (Semialgebraic).
(a) A set D ⊂ R p is called semialgebraic [Bochnak et al., 1998] if it can be represented as where P ij , Q ij are real polynomial functions for 1 ≤ i ≤ s, 1 ≤ j ≤ t.
According to Lojasiewicz [1965], Bochnak et al. [1998] and Shiota [1997, I.2.9, page 52], the class of semialgebraic sets are stable under the operation of finite union, finite intersection, Cartesian product or complementation. Some typical examples include polynomial functions, the indicator function of a semialgebraic set, and the Euclidean norm [Bochnak et al., 1998, page 26].

A.2 KL Property in Deep Learning and Proof of Corollary 1
In the following, we consider the deep neural network training problem. Consider a l-layer feedforward neural network including l − 1 hidden layers of the neural network. Particularly, let d i be the number of hidden units in the i-th hidden layer for i = 1, . . . , l − 1. Let d 0 and d l be the number of units of input and output layers, respectively. Let W i ∈ R d i ×d i−1 be the weight matrix between the (i − 1)-th layer and the i-th layer for any i = 1, . . . l 1 .
According to Theorem 1, one major condition is to verify the introduced Lyapunov function F defined in (11) satisfies the Kurdyka-Lojasiewicz property. For this purpose, we need an extension of semialgebraic set, called the o-minimal structure (see, for instance Coste [1999], van den Dries [1986], Kurdyka [1998], Bolte et al. [2007b]). The following definition is from Bolte et al. [2007b].
Definition 4 (o-minimal structure). An o-minimal structure on (R, +, ·) is a sequence of boolean algebras O n of "definable" subsets of R n , such that for each n ∈ N (i) if A belongs to O n , then A × R and R × A belong to O n+1 ; (ii) if Π : R n+1 → R n is the canonical projection onto R n , then for any A in O n+1 , the set Π(A) belongs to O n ; (iii) O n contains the family of algebraic subsets of R n , that is, every set of the form where p : R n → R is a polynomial function.
(iv) the elements of 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 O (over (R, +, ·)), a function f : R n → R is said to be definable in O if its graph belongs to O n+1 .
According to van den Dries and Miller [1996], Bolte et al. [2007b], there are some important facts of the o-minimal structure, shown as follows.
(i) The collection of semialgebraic sets is an o-minimal structure. Recall the semialgebraic sets are Bollean combinations of sets of the form {x ∈ R n : p(x) = 0, q 1 (x) < 0, . . . , q m (x) < 0}, where p and q i 's are polynomial functions in R n .
(ii) There exists an o-minimal structure that contains the sets of the form where f is real-analytic around [−1, 1] n .
(iii) There exists an o-minimal structure that contains simultaneously the graph of the exponential function R x → exp(x) and all semialgebraic sets.
(iv) The o-minimal structure is stable under the sum, composition, the inf-convolution and several other classical operations of analysis.
The Kurdyka-Lojasiewicz property for the smooth definable function and non-smooth definable function were established in [Kurdyka, 1998, Theorem 1] and [Bolte et al., 2007b, Theorem 11], respectively. Now we are ready to present the proof of Corollary 1.
Because and σ i 's are definable by assumptions, then L n (W, Γ) are definable as compositions of definable functions. Moreover, according to Krantz and Parks [2002], W − Γ 2 and W, g are semi-algebraic and thus definable. Since the group Lasso Ω(Γ) = g Γ 2 is the composition of 2 and 1 norms, and the conjugate of group Lasso penalty is the maximum of group 2 -norm, i.e. Ω * (Γ) = max g Γ g 2 , where the 2 , 1 , and ∞ norms are definable, hence the group Lasso and its conjugate are definable as compositions of definable functions. Therefore, F is definable and hence satisfies Kurdyka-Lojasiewicz inequality by [Kurdyka, 1998, Theorem 1]. The verifications of these specific cases listed in assumptions can be found in the proof of [Zeng et al., 2019, Proposition 1]. This finishes the proof of this corollary.

A.3 Proof of Theorem 1
Our analysis is mainly motivated by a recent paper [Benning et al., 2017], as well as the influential work [Attouch et al., 2013]. According to [Attouch et al., 2013, Lemma 2.6], there are mainly four ingredients in the analysis, that is, the sufficient descent property, relative error property, continuity property of the generated sequence and the Kurdyka-Lojasiewicz property of the function. More specifically, we first establish the sufficient descent property of the generated sequence via exploiting the Lyapunov function F (see, (11)) 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 ofL(W, Γ) and the relation lim k→∞ B g k Ω (Γ k+1 , Γ k ) = 0 established in Lemma 1(i) in Section A.4. Thus, together with the Kurdyka-Lojasiewicz assumption of F , we establish the global convergence of SLBI following by [Attouch et al., 2013, Lemma 2.6].
Let (W ,Γ,ḡ) be a critical point of F , then the following holds By the final inclusion and the convexity of Ω, it impliesḡ ∈ ∂Ω(Γ). Plugging this inclusion into the second inclusion yields αν −1 (Γ −W ) = 0. Together with the first equality imples This finishes the proof of this theorem.

A.4 Sufficient Descent Property along Lyapunov Function
Let P k := (W k , Γ k ), and Q k := (P k , g k−1 ), k ∈ N. In the following, we present the sufficient descent property of Q k along the Lyapunov function F . Lemma. Suppose that L n is continuously differentiable and ∇ L n is Lipschitz continuous with a constant Lip > 0. Let {Q k } be a sequence generated by SLBI with a finite initialization. If 0 < α < 2 κ(Lip+ν −1 ) , then Proof. By the optimality condition of (9a) and also the inclusion p k = [0, g k ] T ∈ ∂R(P k ), there holds where Noting thatL(P ) = L n (W ) + 1 2ν W − Γ 2 2 and by the Lipschitz continuity of ∇ L n (W ) with a constant Lip > 0 implies ∇L is Lipschitz continuous with a constant Lip + ν −1 . This implies L(P k+1 ) ≤L(P k ) + ∇L(P k ), P k+1 − P k + Lip + ν −1 2 P k+1 − P k 2 2 .
Substituting the above inequality into (19) yields Adding some terms in both sides of the above inequality and after some reformulations implies where the final inequality holds for B Thus, we finish the proof of this lemma.
Based on Lemma A.4, we directly obtain the following lemma.
Proof. Note that By the definition ofL (see (2)), where the last inequality holds for the Lipschitz continuity of ∇ L n with a constant Lip > 0, and by (10a), Substituting the above (in)equalities into (26) yields By (10c), it yields Noting that ∇ ΓL (W k , Γ k ) = ν −1 (Γ k − W k ), and after some simplifications yields where the last inequality holds for the triangle inequality and κ −1 > αν −1 by the assumption. By (27), (28), and the definition of H k+1 (25), there holds By (29) and Lemma 1(iv), 1 K K k=1 H k 2 → 0 at a rate of O(1/K). This finishes the proof of this lemma.