Year: 2020

Abstract:
Sparsifying deep neural networks is of paramount interest in many areas, especially when those networks have to be implemented on low-memory devices. In this article, we propose a new formulation of the problem of generating sparse weights for a neural network. By leveraging the properties of standard nonlinear activation functions, we show that the problem is equivalent to an approximate subdifferential inclusion problem. The accuracy of the approximation controls the sparsity. We show that the proposed approach is valid for a broad class of activation functions (ReLU, sigmoid, softmax). We propose an iterative optimization algorithm to induce sparsity whose convergence is guaranteed. Because of the algorithm flexibility, the sparsity can be ensured from partial training data in a minibatch manner. To demonstrate the effectiveness of our method, we perform experiments on various networks in different applicative contexts: image classification, speech recognition, natural language processing, and time-series forecasting.

heading: 1 INTRODUCTION
text: Deep neural networks have evolved to the state-of-the-art techniques in a wide array of applications: computer vision (Simonyan & Zisserman, 2015; He et al., 2016; Huang et al., 2017), automatic speech recognition (Hannun et al., 2014; Dong et al., 2018; Li et al., 2019; Watanabe et al., 2018; Hayashi et al., 2019; Inaguma et al., 2020), natural language processing (Turc et al., 2019; Radford et al., 2019; Dai et al., 2019b; Brown et al., 2020), and time series forecasting (Oreshkin et al., 2020). While their performance in various applications has matched and often exceeded human capabilities, neural networks may remain difficult to apply in real-world scenarios. Deep neural networks leverage the power of Graphical Processing Units (GPUs), which are power-hungry. Using GPUs to make billions of predictions per day, thus comes with a substantial energy cost. In addition, despite their quite fast response time, deep neural networks are not yet suitable for most real-time applications where memory-limited low-cost architectures need to be used. For all those reasons, compression and efficiency have become a topic of high interest in the deep learning community.
Sparsity in DNNs has been an active research topic generating numerous approaches. DNNs achieving the state-of-the-art in a given problem usually have a large number of layers with non-uniform parameter distribution across layers. Most sparsification methods are based on a global approach, which may result in a sub-optimal compression for a reduced accuracy. This may occur because layers with a smaller number of parameters may remain dense, although they may contribute more in terms of computational complexity (e.g., for convolutional layers). Some methods, also known as magnitude pruning, use a hard or soft-thresholding to remove less significant parameters. Soft thresholding techniques achieve a good sparsity-accuracy trade-off at the cost of additional parameters and increased computation time during training. Searching for a hardware efficient network is another area that has been proven quite useful, but it requires a huge amount of computational resources. Convex optimization techniques such as those used in (Aghasi et al., 2017) often rely upon fixed point iterations that make use of the proximity operator (Moreau, 1962). The related concepts are fundamental for tackling nonlinear problems and have recently come into play in the analysis of neural networks (Combettes & Pesquet, 2020a) and nonlinear systems (Combettes & Woodstock, 2020).
This paper shows that the properties of nonlinear activation functions can be utilized to identify highly sparse subnetworks. We show that the sparsification of a network can be formulated as an approximate
subdifferential inclusion problem. We provide an iterative algorithm called subdifferential inclusion for sparsity (SIS) that uses partial training data to identify a sparse subnetwork while maintaining good accuracy. SIS makes even small-parameter layers sparse, resulting in models with significantly lower inference FLOPs than the baselines. For example, SIS for 90% sparse MobileNetV3 on ImageNet-1K achieves 66.07% top-1 accuracy with 33% fewer inference FLOPs than its dense counterpart and thus provides better results than the state-of-the-art method RigL. For non-convolutional networks like Transformer-XL trained on WikiText-103, SIS is able to achieve 70% sparsity while maintaining 21.1 perplexity score. We evaluate our approach across four domains and show that our compressed networks can achieve competitive accuracy for potential use on commodity hardware and edge devices.

heading: 2 RELATED WORK

heading: 2.1 INDUCING SPARSITY POST TRAINING
text: Methods inducing sparsity after a dense network is trained involve several pruning and fine-tuning cycles till desired sparsity and accuracy are reached (Mozer & Smolensky, 1989; LeCun et al., 1990; Hassibi et al., 1993; Han et al., 2015; Molchanov et al., 2017; Guo et al., 2016; Park et al., 2020). (Renda et al., 2020) proposed weight rewinding technique instead of vanilla fine-tuning post-pruning. Net-Trim algorithm (Aghasi et al., 2017) removes connections at each layer of a trained network by convex programming. The proposed method works for networks using rectified linear units (ReLUs). Lowering rank of parameter tensors (Jaderberg et al., 2014; vahid et al., 2020; Lu et al., 2016), removing channels, filters and inducing group sparsity (Wen et al., 2016; Li et al., 2017; Luo et al., 2017; Gordon et al., 2018; Yu et al., 2019; Liebenwein et al., 2020) are some methods that take network structure into account. All these methods rely on pruning and fine-tuning cycle(s) often from full training data.

heading: 2.2 INDUCING SPARSITY DURING TRAINING
text: Another popular approach has been to induce sparsity during training. This can be achieved by modifying the loss function to consider sparsity as part of the optimization (Chauvin, 1989; CarreiraPerpiñán & Idelbayev, 2018; Ullrich et al., 2017; Neklyudov et al., 2017). Dynamically pruning during training (Zhu & Gupta, 2018; Bellec et al., 2018; Mocanu et al., 2018; Dai et al., 2019a; Lin et al., 2020b) by observing network flow. (Mostafa & Wang, 2019; Dettmers & Zettlemoyer, 2020; Evci et al., 2020) computes weight magnitude and reallocates weights at every step. Bayesian priors (Louizos et al., 2017), L0, L1 regularization (Louizos et al., 2018), and variational dropout (Molchanov et al., 2017) get accuracy comparable to (Zhu & Gupta, 2018) but at a cost of 2×memory and 4× computations during training. (Liu et al., 2019; Savarese et al., 2020; Kusupati et al., 2020; Lee, 2019; Xiao et al., 2019; Azarian et al., 2020) have proposed learnable sparsity methods through training of the sparse masks and weights simultaneously with minimal heuristics. Although these methods are cheaper than pruning after training, they need at least the same computational effort as training a dense network to find a sparse sub-network. This makes them expensive when compressing big networks where the number of parameters ranges from hundreds of millions to billions (Dai et al., 2019b; Li et al., 2019; Brown et al., 2020).

heading: 2.3 TRAINING SPARSELY INITIALIZED NETWORKS
text: (Frankle & Carbin, 2019) showed that it is possible to find sparse sub-networks that, when trained from scratch, were able to match or even outperform their dense counterparts. (Lee et al., 2019) presented SNIP, a method to estimate, at initialization, the importance that each weight could have later during training. In (Lee et al., 2020) the authors perform a theoretical study of pruning at initialization from a signal propagation perspective, focusing on the initialization scheme. Recently, (Wang et al., 2020) proposed GraSP, a different method based on the gradient norm after pruning, and showed a significant improvement for moderate levels of sparsity. (Ye et al., 2020) starts with a small subnetwork and progressively grow it to a subnetwork that is as accurate as its dense counterpart. (Tanaka et al., 2020) proposes SynFlow that avoids flow collapse of a pruned network during training. (Jorge et al., 2020) proposed FORCE, an iterative pruning method that progressively removes a small number of weights. This method is able to achieve extreme sparsity at little accuracy expense. These
methods are not usable for big pre-trained networks and are expensive as multiple training rounds are required for different sparse models depending on deployment scenarios (computing devices).

heading: 2.4 EFFICIENT NEURAL ARCHITECTURE SEARCH
text: Hardware-aware NAS methods (Zoph et al., 2018; Real et al., 2019; Cai et al., 2018; Wu et al., 2019; Tan et al., 2019; Cai et al., 2019; Howard et al., 2019) directly incorporate the hardware feedback into efficient neural architecture search. (Cai et al., 2020) proposes to learn a single network composing of a large number of subnetworks from which a hardware aware subnetwork can be extracted in linear time. (Lin et al., 2020a) proposes a similar approach wherein they identify subnetworks that can be run efficiently on microcontrollers (MCUs).
Our proposed algorithm applies to possibly large pre-trained networks. In contrast with methods presented in Section 2.1, ours can use a small amount of training data during pruning and fewer epochs during fine-tuning. As we will see in the next section, a key feature of our approach is that it is based on a fine analysis of the mathematical properties of activation functions, so allowing the use of powerful convex optimization tools that offer sound convergence guarantees.

heading: 3 PROPOSED METHOD

heading: 3.1 VARIATIONAL PRINCIPLES
text: A basic neural network layer can be described by the relation:
y = R(Wx+ b) (1)
where x ∈ RM is the input, y ∈ RN the output, W ∈ RN×M is the weight matrix, b ∈ RN the bias vector, and R is a nonlinear activation operator from RN to RN . A key observation is that most of the activation operators currently used in neural networks are proximity operators of convex functions (Combettes & Pesquet, 2020a;b). We will therefore assume that there exists a proper lower-semicontinuous convex function f from RN to R ∪ {+∞} such that R = proxf . We recall that f is a proper lower-semicontinuous convex function if the area overs its graph, its epigraph{
(y, ξ) ∈ RN × R ∣∣ f(y) 6 ξ}, is a nonempty closed convex set. For such a function the proximity
operator of f at z ∈ RN (Moreau, 1962) is the unique point defined as
proxf (z) = argmin p∈RN
1 2 ‖z − p‖2 + f(p). (2)
It follows from standard subdifferential calculus that Eq. (1) can be re-expressed as the following inclusion relation:
Wx+ b− y ∈ ∂f(y), (3)
where ∂f(y) is the Moreau subdifferential of f at y defined as ∂f(y) = { t ∈ RN ∣∣ (∀z ∈ RN )f(z) > f(y) + 〈t | z − y〉}. (4) The subdifferential constitutes a useful extension of the notion of differential, which is applicable to nonsmooth functions. The set ∂f(y) is closed and convex and, if y satisfies Eq. (1), it is nonempty. The distance to this set of a point z ∈ RN is given by
d∂f(y)(z) = inf t∈∂f(y)
‖z − t‖. (5)
We thus see that the subdifferential inclusion in Eq. (3) is also equivalent to
d∂f(y)(Wx+ b− y) = 0. (6)
Therefore, a suitable accuracy measure for approximated values of the layer parameters (W, b) is d∂f(y)(Wx+ b− y).

heading: 3.2 OPTIMIZATION PROBLEM
text: Compressing a network consists of a sparsification of its parameters while keeping a satisfactory accuracy. Assume that, for a given layer, a training sequence of input/output pairs is available which results from a forward pass performed on the original network for some input dataset of length K. The training sequence is split in J minibatches of size T so that K = JT . The j-th minibatch with j ∈ {1, . . . , J} is denoted by (xj,t, yj,t)16t6T . In order to compress the network, we propose to solve the following constrained optimization problem.
Problem 1 We want to minimize (W,b)∈C g(W, b) (7)
with
C = { (W, b) ∈ RN×M × RN ∣∣ (∀j ∈ {1, . . . , J}) T∑
t=1
d2∂f(yj,t)(Wxj,t + b− yj,t) 6 Tη } , (8)
where g is a sparsity measure defined on RN×M × RN and η ∈ [0,+∞[ is some accuracy tolerance. Since, for every j ∈ {1, . . . , J}, the function (W, b) 7→ ∑T t=1 d 2 ∂f(yj,t)
(Wxj,t + b − yj,t) is continuous and convex, C is a closed and convex subset of RN×M × RN . In addition, this set is nonempty when there exist W ∈ RN×M and b ∈ RN such that, for every j ∈ {1, . . . , J} and t ∈ {1, . . . , T}, d2∂f(yj,t)(Wxj,t + b− yj,t) = 0. As we have seen in Section 3.1, this condition is satisfied when (W, b) are the parameters of the uncompressed layer. Often, the sparsity of the weight matrix is the determining factor whereas the bias vector represents a small number of parameters, so that we can make the following assumption.
Assumption 2 For every W ∈ RN×M and b ∈ RN , g(W, b) = h(W ) where h is a function from RN×M to R ∪ {+∞}, which is lower-semicontinuous, convex, and coercive (i.e. lim‖W‖F→+∞ h(W ) = +∞). In addition, there exists (W, b) ∈ C such that h(W ) < +∞ and there exists (j∗, t∗) ∈ {1, . . . , J} × {1, . . . , T} such that yj∗,t∗ lies in the interior of the range of R.
Under this assumption, the existence of a solution to Problem 1 is guaranteed (see Appendix A). A standard choice for such a function is the `1-norm of the matrix elements, h = ‖ · ‖1, but other convex sparsity measures could also be easily incorporated within this framework, e.g. group sparsity measures. Another point worth being noticed is that constraints other than (8) could be imposed. For example, one could make the following alternative choice for the constraint set
C = { (W, b) ∈ RN×M × RN ∣∣ sup j∈{1,...,J},t∈{1,...,T} d∂f(yj,t)(Wxj,t + b− yj,t) 6 √ η } . (9)
Although the resulting optimization problem could be tackled by the same kind of algorithm as the one we will propose, Constraint (8) leads to a simpler implementation.

heading: 3.3 OPTIMIZATION ALGORITHM
text: A standard proximal method for solving Problem 1 is the Douglas-Rachford algorithm (Lions & Mercier, 1979; Combettes & Pesquet, 2007). This algorithm alternates between a proximal step aiming at sparsifying the weight matrix and a projection step allowing a given accuracy to be reached. Assume that a solution to Problem 1 exists. Then, this algorithm reads as shown on the top of the next page. The Douglas-Rachford algorithm uses parameters γ ∈ ]0,+∞[ and (λn)n∈N in ]0, 2[ such that∑ n∈N λn(2− λn) = +∞. Throughout this article, projS denotes the projection onto a nonempty closed convex set S. Under these conditions, the sequence (Wn, bn)n∈N generated by Algorithm 1 is guaranteed to converge to a solution to Problem 1 if there exists (W, b) ∈ C such W is a point in the relative interior of the domain of h Combettes & Pesquet (2007) (see illustrations in Appendix E).
Algorithm 1: Douglas-Rachford algorithm for network compression
Initialize :Ŵ0 ∈ RN×M and b0 ∈ RN for n = 0, 1, . . . do
Wn = proxγh(Ŵn)
(W̃n, b̃n) = projC(2Wn − Ŵn, bn) Ŵn+1 = Ŵn + λn(W̃n −Wn) bn+1 = bn + λn(̃bn − bn).
The proximity operator of function γh has a closed-form for standard choices of sparsity measures1. For example, when h = ‖ · ‖1, this operator reduces to a soft-thresholding (with threshold value γ) of the input matrix elements. In turn, since the convex set C has an intricate form, an explicit expression of projC does not exist. Finding an efficient method for computing this projection for large datasets thus constitutes the main challenge in the use of the above Douglas-Rachford strategy, which we will discuss in the next section.

heading: 3.4 COMPUTATION OF THE PROJECTION ONTO THE CONSTRAINT SET
text: For every mini-batch index j ∈ {1, . . . , J}, let us define the following convex function:
(∀(W, b) ∈ RN×M × RN ) cj(W, b) = T∑ t=1 d2∂f(yj,t)(Wxj,t + b− yj,t)− Tη. (10)
Note that, for every j ∈ {1, . . . , J}, function cj is differentiable and its gradient at (W, b) ∈ RN×M × RN is given by
∇cj(W, b) = (∇Wcj(W, b),∇bcj(W, b)), (11)
where
∇Wcj(W, b) = 2 T∑ t=1 ej,tx > j,t, ∇bcj(W, b) = 2 T∑ t=1 ej,t (12)
with (∀t ∈ {1, . . . , T}) ej,t = Wxj,t + b− yj,t − proj∂f(yj,t)(Wxj,t + b− yj,t). (13)
A pair of weight/bias parameters belongs to C if and only if it lies in the intersection of the 0-lower level sets of the functions (cj)16j6J . To compute the projection of some (W, b) ∈ RN×M × RN onto this intersection, we use Algorithm 2 (‖ · ‖F denotes here the Frobenius norm). This iterative algorithm has the advantage of proceeding in a minibatch manner. It allows us to choose the mini-batch index jn at iteration n in a quasi-cyclic manner. The simplest rule is to activate each minibatch once within J successive iterations of the algorithm so that they correspond to an epoch. The proposed algorithm belongs to the family of block-iterative outer approximation schemes for solving constrained quadratic problems, which was introduced in (Combettes, 2003). The convergence of the sequence (Wn, bn)n∈N generated by Algorithm 2 to projC(W, b) is thus guaranteed. One of the main features of the algorithm is that it does not require to perform any projection onto the 0-lower level sets of the functions cj , which would be intractable due to their expressions. Instead, these projections are implicitly replaced by subgradient projections, which are much easier to compute in our context.

heading: 3.5 DEALING WITH VARIOUS NONLINEARITIES
text: For any choice of activation operator R, we have to calculate the projection onto ∂f(y) for every vector y satisfying Eq. (1). This projection is indeed required in the computation of the gradients of functions (cj)16j6J , as shown by Eq. (13). Two properties may facilitate this calculation. First, if f is differentiable at y, then ∂f(y) reduces to a singleton containing the gradient ∇f(y) of
1http://proximity-operator.net
Algorithm 2: Minibatch algorithm for computing projC(W, b) Initialize :W0 = W and b0 = b for n = 0, 1, . . . do
Select a batch of index jn ∈ {1, . . . , J} if cjn(Wn, bn) > 0 then
Compute∇Wcjn(Wn, bn) and ∇bcjn(Wn, bn) by using Eq. (12) and Eq. (13) δWn =
cjn (Wn,bn)∇Wcjn (Wn,bn) ‖∇Wcjn,n(Wn,bn)‖2F+‖∇bcjn (Wn,bn)‖2
δbn = cjn (Wn,bn)∇bcjn (Wn,bn)
‖∇Wcjn,n(Wn,bn)‖2F+‖∇bcjn (Wn,bn)‖2
πn = tr((W0 −Wn)>δWn) + (b0 − bn)>δbn µn = ‖W0 −Wn‖2F + ‖b0 − bn‖2 νn = ‖δWn‖2F + ‖δbn‖2 ζn = µnνn − π2n if ζn = 0 and πn > 0 then
Wn+1 = Wn − δWn bn+1 = bn − δbn
else if ζn > 0 and πnνn > ζn then Wn+1 = W0 − (1 + πnνn )δWn bn+1 = b0 − (1 + πnνn )δbn else Wn+1 = Wn +
νn ζn (πn(W0 −Wn)− µnδWn) bn+1 = bn +
νn ζn (πn(b0 − bn)− µnδbn)
else Wn+1 = Wn bn+1 = bn
f at y, so that, for every z ∈ RN , proj∂f(y)(z) = ∇f(y). Second, R is often separable, i.e. consists of the application of a scalar activation function ρ : R→ R to each component of its input argument. According to our assumptions, there thus exists a proper lower-semicontinuous convex function ϕ from R to R ∪ {+∞} such that ρ = proxϕ and, for every z = (ζ(k))16k6N ∈ RN , f(z) = ∑N k=1 ϕ(ζ
(k)). This implies that, for every z = (ζ(k))16k6N ∈ RN , proj∂f(y)(z) = (proj∂ϕ(υ(k))(ζ
(k)))16k6N , where the components of y are denoted by (υ(k))16k6N . Based on these properties, a list of standard activation functions ρ is given in Table 1, for which we provide the associated expressions of the projection onto ∂ϕ. The calculations are detailed in Appendix B.
An example of non-separable activation operator frequently employed in neural network architectures is the softmax operation defined as: (∀z = (ζ(k))16k6N ∈ RN ) R(z) = (
exp(ζ(k))∑N j=1 exp(ζ (j)) ) 16k6N .
It is shown in Appendix C that, for every y = (υ(k))16k6N in the range of R,
(∀z ∈ RN ) proj∂f(y)(z) = Q(y) + 1>(z −Q(y))
N 1, (14)
where 1 = [1, . . . , 1]> ∈ RN and Q(y) = (ln υ(k) + 1− υ(k))16k6N .

heading: 3.6 SIS ON MULTI-LAYERED NETWORKS
text: Algorithm 3 describes how we make use of SIS for a multi-layered neural network. We use a pretrained network and part of the training sequence to extract layer-wise input-output features. Then we apply SIS on each individual layer l by passing η, layer parameters (W (l), b(l)) and extracted input-output features (Y (l−1), Y (l)) to Algorithm 1. The benefit of applying SIS to each layer independently is that we can apply SIS on all the layers of a network in parallel. This reduces the time required to process the whole network and compute resources are optimally utilized.
Algorithm 3: Parallel SIS for multi-layered network Input: input sequence X ∈ RM×K , compression parameter η > 0, weight matrices W (1), . . . ,W (L), and bias vectors b(1), . . . , b(L) Y (0) ← X for l = 1, . . . , L do
Y (l) = R`(W (l)>Y (l−1) + b(l))
Ŵ (l), b̂(l) ← SIS(η,W (l), b(l), Y (l), Y (l−1)) Output: Ŵ (1), . . . , Ŵ (L) and b̂(1), . . . , b̂(L)

heading: 4 EXPERIMENTS
text: In this section, we conduct various experiments to validate the effectiveness of SIS in terms of the test accuracy vs. sparsity and inference time FLOPs vs. sparsity by comparing against RigL (Evci et al., 2020). We also include SNIP (Lee et al., 2019), GraSP (Wang et al., 2020), SynFlow (Tanaka et al., 2020), STR (Kusupati et al., 2020), and FORCE (Jorge et al., 2020). These methods start training from a sparse network and have some limitations when compared to methods that prune a pretrained network Blalock et al. (2020); Gale et al. (2019). For a fair comparison we also include LRR (Renda et al., 2020) which uses a pretrained network and multiple rounds of pruning and retraining by leveraging learning rate rewinding. The experimental setup is described in Appendix D.

heading: 4.1 MODERN CONVNETS ON CIFAR AND IMAGENET
text: We compare SIS with competitive baselines on CIFAR-10/100 for three different sparsity regimes 90%, 95%, 98%, and the results are listed in Table 2. It can be observed that LRR, RigL and SIS are able to maintain high accuracy with increasing sparsity. LRR performs better than both RigL and SIS for VGG19 on CIFAR-10 at 90% and 95% sparsity. When compared to SNIP, our method achieves impressive performance for VGG19 on CIFAR-100 (58.46 -> 71.17). In the case of ResNet50, SIS outperforms all the other methods for CIFAR-10/100 except for CIFAR-100 at 90%.
Due to its small size and controlled nature, CIFAR-10/100 may not appear sufficient to draw solid conclusions. We thus conduct further experiments on ImageNet using ResNet50 and MobileNets. Table 3 shows that, in the case of ResNet50, LRR performs marginally better than SIS at 60% sparsity.
At 80%, 90%, and 96.5% sparsity SIS outperforms all other methods. For all sparsity regimes, SIS achieves least inference FLOPs. RigL achieves these results in same training time as SIS with less training time FLOPs. This may be related to the fact that SIS can achieve better compression in the last layer before SoftMax. MobileNets are compact architectures designed specifically for resource-constrained devices. Table 4 shows results for RigL and SIS on MobileNets. We observe that SIS outperforms all MobileNet versions at 75% sparsity level. For a 90% sparsity level, SIS outperforms RigL for MobileNet V1 and V3 whereas, for MobileNetV2, RigL performs slightly better than SIS at 90% sparsity level. In all the cases, we can see that the resulting SIS sparse network uses fewer FLOPs than RigL. A possible explanation for this fact is that SIS leverages activation function properties during the sparsification process.

heading: 4.2 SEQUENTIAL TASKS
text: Jasper on LibriSpeech. Jasper is a speech recognition model that uses 1D convolutions. The trained network is a 333 million parameter model and has a word error rate (WER) of 12.2 on the test set. We apply SIS on this network and compare it with RigL and SNIP in terms of sparsity. Table 5 reports WER and inference FLOPs for all three methods. SIS marginally performs better than SIS on this task in terms of WER and FLOPs for 70% sparsity. The main advantage of our approach lies in the fact that we can use a single pre-trained Jasper network and achieve different sparsity level for different types of deployment scenarios with less computational resources than RigL.
Transformer-XL on WikiText-103. Transformer-XL is a language model with 246 million parameters. The trained network on WikiText-103 has a perplexity score (PPL) of 18.6. In Table 5, we see that SIS performs better than SNIP and RigL in terms of PPL and has 68% fewer inference FLOPs. This is due to the fact that large language models can be efficiently trained and then compressed easily, but training a sparse sub-network from scratch is hard (Li et al., 2020), as is the case with SNIP and RigL. SNIP uses one-shot pruning to obtain a random sparse sub-network, whereas RigL is able to change its structure during training, which allows it to perform better than SNIP.
N-BEATS on M4. N-BEATS is a very deep residual fully-connected network to perform forecasting in univariate time-series problems. It is a 14 million parameter network. The Symmetric Mean Absolute Percentage Error (SMAPE) of the dense network on the M4 dataset is 8.3%. We apply SIS on this network and compare its performance with respect to RigL and SIS. SIS performs better than both methods and results in 65% fewer inference FLOPs.

heading: 5 CONCLUSION
text: In this article, we have proposed a novel method for sparsifying neural networks. The compression problem for each layer has been recast as the minimization of a sparsity measure under accuracy constraints. This constrained optimization problem has been solved by means of advanced convex optimization tools. The resulting SIS algorithm is i) reliable in terms of iteration convergence guarantees, ii) applicable to a wide range of activation operators, and iii) able to deal with large datasets split into mini-batches. Our numerical tests demonstrate that the approach is not only appealing from a theoretical viewpoint but also practically efficient.

heading: A EXISTENCE OF A SOLUTION TO PROBLEM 1 UNDER ASSUMPTION 2
text: Under Assumption 2, Problem 1 is equivalent to
minimize (W,b)∈C h(W ) (15)
with C = { (W, b) ∈ RN×M × RN ∣∣ max j∈{1,...,J} cj(W, b) 6 0 } , (16)
where the functions (cj)16j6J are defined in Eq. (10). These functions being convex, Φ = maxj∈{1,...,J} cj is convex (Bauschke & Combettes, 2019, Proposition 8.16). We deduce that Ψ = infb∈RN Φ(·, b) is also a convex function (Bauschke & Combettes, 2019, Proposition 8.35). Since Φ > −ηT , Ψ is finite valued. It is thus continuous on RN×M (Bauschke & Combettes, 2019, Corollary 8.40). Let us now consider the problem:
minimize W∈lev60 Ψ h(W ) (17)
where lev60 Ψ is the 0-lower level set of Ψ defined as lev60 Ψ = { W ∈ RN×M ∣∣ Ψ(W ) 6 0}, (18) Ψ being both convex and continuous, lev60 Ψ is closed and convex. According to Assumption 2, there exists (W, b) ∈ RN×M × RN such that h(W ) < +∞ and Φ(W, b) 6 0, which implies that Ψ(W ) 6 0. This shows that lev60 Ψ has a nonempty intersection with the domain of h. By invoking now the coercivity property of h, the existence of a solution Ŵ to Problem (17) is guaranteed by standard convex analysis results (Bauschke & Combettes, 2019, Theorem 11.10).
To show that (Ŵ , b̂) is a solution to (15), it is sufficient to show that there exists b̂ ∈ RN such that Φ(Ŵ , b̂) = Ψ(Ŵ ). This is equivalent to prove that there exists a solution b̂ to the problem:
minimize b∈RN Φ(Ŵ , b). (19)
We know that Φ(Ŵ , ·) is a continuous function. In addition, we have assumed that there exists (j∗, t∗) ∈ {1, . . . , J} × {1, . . . , T} such that yj∗,t∗ is an interior point of R(RN ), which is also equal to the domain of ∂f and a subset of the domain of f . Since f is continuous on the interior of its domain, ∂f(yj∗,t∗) is bounded (Bauschke & Combettes, 2019, Proposition 16.17(ii)). Then d∂f(yj∗,t∗ ) is coercive, hence cj∗(Ŵ , ·) is coercive, and so is Φ(Ŵ , ·) > cj∗(Ŵ , ·). The existence of b̂ thus follows from the Weierstrass theorem.

heading: B RESULTS IN TABLE 1
text: The results are derived from the expression of the convex function ϕ associated with each activation function ρ (Combettes & Pesquet, 2020a, Section 2.1) (Combettes & Pesquet, 2020b, Section 3.2).
Sigmoid
(∀ζ ∈ R) ϕ(ζ) =  (ζ + 1/2) ln(ζ + 1/2) + (1/2− ζ) ln(1/2− ζ)− 1 2 (ζ2 + 1/4) if |ζ| < 1/2 −1/4 if |ζ| = 1/2 +∞ if |ζ| > 1/2.
(20)
The range of the Sigmoid function is ]− 1/2, 1/2[ and the above function is differentiable on this interval and its derivative at every υ ∈]− 1/2, 1/2[ is
ϕ′(υ) = ln(υ + 1/2)− ln(υ − 1/2)− υ. (21) We deduce that, for every ζ ∈ R, proj∂ϕ(υ)(ζ) = ϕ′(υ).
Arctangent
(∀ζ ∈ R) ϕ(ζ) =
{ − 2π ln ( cos ( πζ 2 )) − 12ζ
2, if |ζ| < 1 +∞, if |ζ| > 1.
(22)
By proceeding for this function similarly to the Sigmoid function, we have, for every υ ∈ ρ(R) = ]− 1, 1[,
(∀ζ ∈ R) proj∂ϕ(υ)(ζ) = ϕ′(υ) = tan(πυ/2)− υ. (23)
ReLU
(∀ζ ∈ R) ϕ(ζ) = {
0 if ζ > 0 +∞ otherwise. (24)
For every υ ∈ ρ(R) = [0,+∞[, we have
∂ϕ(υ) = { {0} if υ > 0 ]−∞, 0] if υ = 0. (25)
We deduce that
(∀ζ ∈ R) proj∂ϕ(υ)(ζ) = {
0 if υ > 0 or ζ > 0 ζ otherwise.
(26)
Leaky ReLU
(∀ζ ∈ R) ϕ(ζ) = {
0, if ζ > 0 (1/α− 1)ζ2/2 if ζ 6 0. (27)
Since this function is differentiable on R, for every υ ∈ R, (∀ζ ∈ R) proj∂ϕ(υ)(ζ) = ϕ′(υ) = {
0 if υ > 0 (1/α− 1)υ otherwise. (28)
Capped ReLU
(∀ζ ∈ R) ϕ(ζ) = {
0 if ζ ∈ [0, α] +∞ otherwise. (29)
We have thus for every υ ∈ [0, α],
∂ϕ(υ) =  {0} if υ ∈]0, α[ ]−∞, 0] if υ = 0 [0,+∞[ if υ = α.
(30)
This leads to
(∀ζ ∈ R) proj∂ϕ(υ)(ζ) =  ζ if (υ = 0 and ζ < 0)
or (υ = α and ζ > 0) 0 otherwise.
(31)
ELU
(∀ζ ∈ R) ϕ(ζ) =  0 if ζ > 0; (ζ + α) ln ( ζ+α α ) − ζ − ζ 2 2 , if − α < ζ < 0 α− α 2
2 , if ζ = −α +∞, if ζ < −α.
(32)
This function being differentiable on ρ(R) =]− α,+∞[, we have for every υ ∈]− α,+∞[,
(∀ζ ∈ R) proj∂ϕ(υ)(ζ) = ϕ′(υ) =
{ 0 if υ > 0
ln ( υ+α α ) − υ otherwise. (33)
QuadReLU Unlike the previous ones, this function does not seem to have been investigated before. It can be seen as a surrogate to the hard swish activation function, which is not a proximal activation function. Let us define
(∀ζ ∈ R) ϕ(ζ) =  +∞ if ζ < 0 − ζ 2 2 + 4 3 √ αζ3/2 − αζ if ζ ∈ [0, α] ζ2
2 − αζ + α2 3 if ζ > α. (34)
ϕ is a lower-semicontinuous convex function whose subdifferential is
(∀υ ∈ [0,+∞[) ∂ϕ(υ) =  ]−∞,−α] if υ = 0 {−υ + 2 √ αυ − α} if υ ∈]0, α]
{υ − α} if υ > α. (35)
From the definition of the proximity operator, for every (υ, ζ) ∈ R2, we have υ = proxϕ(ζ) if and only if
ζ ∈ υ + ∂ϕ(υ) ⇔  ζ ∈]−∞,−α] if υ = 0 ζ = 2 √ αυ − α if υ ∈]0, α]
ζ = 2υ − α if υ > α.
⇔ υ =  0 if ζ ∈]−∞,−α] (ζ+α)2 4α if ζ ∈]− α, α] ζ+α
2 if ζ > α. (36)
This shows that proxϕ(ζ) = (4α) −1(ζ + α) ReLU2α(ζ + α). In addition, for every υ ∈ [0,+∞[, it follows from Eq. (35) that the projection onto ∂f(υ) is
(∀ζ ∈ R) proj∂f(υ)(ζ) =  υ if υ = 0 and ζ 6 −α −υ + 2 √ αυ − α if υ ∈]0, α]
or (υ = 0 and ζ > −α) υ − α if υ > α.
(37)

heading: C SOFTMAX ACTIVATION
text: Let C denote the closed hypercube [0, 1]N , let V be the vector hyperplane defined as
V = { z = (ζ(k))16k6N ∈ RN ∣∣ N∑ k=1 ζ(k) = 0 } , (38)
and let A be the affine hyperplane defined as
A = { z = (ζ(k))16k6N ∈ RN ∣∣ N∑ k=1 ζ(k) = 1 } = V + u, (39)
where u = [1, . . . , 1]> = 1/N ∈ RN . If R is the Softmax activation operator, the convex function f such that proxf = R is (Combettes & Pesquet, 2020a, Example 2.23):
(∀z = (ζ(k))16k6N ) f(z) =
{∑N i=1 ϕ(ζ
(k)) if z ∈ C ∩A +∞ otherwise,
(40)
where
(∀ζ ∈ [0,+∞[) ϕ(ζ) = ζ ln ζ − ζ 2
2 (41)
(with the convention 0 ln 0 = 0). The latter function can be extended on R, say by a quadratic function on ]−∞, 0[, yielding a convex function ϕ̃ which is differentiable on R. We have then
(∀z = (ζ(k))16k6N ) ∈ RN ) f(z) = N∑ i=1 ϕ̃(ζ(k)) + ιC∩A(z), (42)
where ιC∩A denotes the indicator function of the intersection of C and A (equal to 0 on this set and +∞ elsewhere). It then follows from standard subdifferential calculus rules that, for every y = (υ(k))16k6N ∈ RN ,
∂f(y) = (ϕ̃′(υ(k)))16k6N +NC(y) +NA(y), (43)
where ϕ̃′ is the derivative of ϕ̃ and ND denotes the normal cone to a nonempty closed convex set D, which is defined as ND(y) = { t ∈ RN ∣∣ (∀z ∈ D) 〈t | z − y〉 6 0}. (44) Thus NA(y) = NV (y) is the orthogonal space V ⊥ of V . Let us now assume that y ∈ R(RN ) ⊂]0, 1[N . Then, since y is an interior point of C, NC(y) = {0}. We then deduce from Eq. (43) that
∂f(y) = Q(y) + V ⊥, (45)
where Q(y) = (ϕ′(υ(k)))16k6N = (ln υ(k) + 1− υ(k))16k6N . It follows that, for every z ∈ RN , proj∂f(y)(z) = Q(y) + projV ⊥(z −Q(y)). (46)
By using the expression of the projection projV = Id−projV ⊥ onto hyperplane V , we finally obtain
proj∂f(y)(z) = Q(y) + 1>(z −Q(y))
N 1. (47)

heading: D EXPERIMENTAL SETUP
text: PyTorch is employed to implement our method. We use and extend SNIP and RigL code available here2, LRR3, GraSP4, SynFlow5, STR6, and FORCE7. In order to manage our experiments we use Polyaxon8 on a Kubernetes9 cluster and use five computing nodes with eight V100 GPUs each. Floating point operations per second (FLOPs) is calculated as equal to one multiply-add accumulator using the code10.
SIS has the following parameters: number of iterations of Algorithm 1, number of iterations of Algorithm 2, step size parameter γ in Algorithm 1, constraint bound parameter η used to control the sparsity, and relaxation parameter λn ≡ λ of Algorithm 1. In our experiments, the maximum numbers of iterations of Algorithm 1 and Algorithm 2 are set to 2000 and 1000, respectively. λ is set to 1.5 and γ is set to 0.1 for all the SIS experiments. η value depends on the network and dataset. With few experiments, we search for a good η value that gives suitable sparsity and accuracy.
VGG19 and ResNet50 on CIFAR-10/100. We train VGG19 on CIFAR-10 for 160 epochs with a batch size of 128, learning rate of 0.1 and weight decay of 5× 10−4 applied at epochs 81 and 122. A momentum of 0.9 is used with stochastic gradient descent (SGD). We make use of 1000 images per training class when using SIS. We fine-tune the identified sparse subnetwork for 10 epochs at a learning rate of 10−3. For CIFAR-100 we keep the same training hyperparameters as for CIFAR-10. When applying SIS to the dense network, we use 300 images per class from the training samples. We fine-tune the identified sparse subnetwork for 40 epochs on the training set with a learning rate of 10−3. ResNet50 employs the same hyperparameters as VGG19, except the weight decay that we set to 10−4. When applying SIS to train dense ResNet50, we use the same partial training set and the same hyperparameters during fine-tuning. In case of VGG19 for CIFAR-10 and CIFAR-100, we found that η values in range (1.5, 2) works best for sparsity range (90%, 98%). In case of ResNet50, η values in range (1, 2) is used.
2https://github.com/google-research/rigl 3https://github.com/lottery-ticket/rewinding-iclr20-public/tree/master/
vision/gpu-src/official 4https://github.com/alecwangcq/GraSP 5https://github.com/ganguli-lab/Synaptic-Flow 6https://github.com/RAIVNLab/STR 7https://github.com/naver/force 8https://github.com/polyaxon/polyaxon 9https://kubernetes.io/
10https://github.com/Lyken17/pytorch-OpCounter
ResNet50 on ImageNet We use the weights of ResNet50 pre-trained on ImageNet available at PyTorch hub11. When applying SIS to the dense pre-trained network we use 20% samples per class from the training set. We fine-tune the identified sparse subnetwork for 40 epochs on the training set with a learning rate of 10−4. We use different η values in range (0.7, 1.5) for sparsity range (60%, 90%). We found that η = 2.3 achieves 96.5% sparsity.
MobileNets on ImageNet We use MobileNetV1 dense pre-trained model from here12 and MobileNetV2 from PyTorch hub13. In case of MobileNetV3, we replace the hard swish activation function used in the original paper Howard et al. (2019) with our QuadReLU function (see the last row of Table 1). We use hyperparameters provided in the original paper to train MobileNetV3. When applying SIS to the dense pre-trained MobileNets, we use 20% samples per class from the training set. We fine-tune the identified sparse subnetwork for 30 epochs on the training set with a learning rate of 10−4. For MobileNets, we search η values in range (0.6, 1.75) for sparsity range (75, 90).
Jasper on LibriSpeech A BxR Jasper network has B blocks, each consisting of R repeating sub-blocks. Each sub-block consists of 1D-Convolution, Batch Normalization, ReLU activation, and Dropout. The kernel size of convolutions increases with depth. The network has one convolution block at the beginning and three at the end. We train a network of 13 encoding blocks and one decoding block, having 54 1D-Convolution layers on the LibriSpeech dataset. The total number of parameters in our trained network is 333 million. Jasper network is trained on train-clean-100, train-clean-360, and train-other-500 splits of the LibriSpeech dataset (Panayotov et al., 2015). The training configuration can be found here14. We use train-clean-100 when using SIS. We fine-tune the identified sparse sub-network on the completed training set for ten epochs with a learning rate of 10−4. We use η values in range (0.6, 1.75) for sparsity range (70, 90).
Transformer-XL on WikiText-103 We train the Transformer-XL network (Dai et al., 2019b) on the base version of WikiText-103 (Merity et al., 2017). We use the training configuration available here15. We use 10% of the training set articles when using SIS. We use η values in range (0.5, 0.75) for sparsity range (40, 70).
N-BEATS on M4 We train the interpretable architecture network of N-BEATS on the M4 dataset. The trained network has six residual blocks. Each block consists of four fully-connected layers and two linear projection layers. With 24 fully-connected layers, this network has 14 million trainable parameters. To compare different methods, we only train a single network on a 48-hour window instead of 180 networks on different timescales. We use the training configuration available here16. The training set has 50K time-series samples. We use 10K training samples to generate a sparse sub-network using SIS. We use η values in range (0.75, 1.5) for sparsity range (70, 90).

heading: E EMPIRICAL CONVERGENCE ANALYSIS
text: We illustrate the convergence of our method on LeNet-FCN trained on MNIST. LeNet-FCN is a fully-connected network having four layers with 784-300-1000-300-10 nodes (two 300 nodes and one 1000 node hidden layers). Figure 1 shows the convergence of SIS when applied to dense LeNet-FCN. We observe that the convergence is smooth and SIS finds a global solution for the first (ReLU activated) and last (softmax) layer cases. This fact is in agreement with our theoretical claims. SIS attains a sparsity of 99.21% at an error of 1.86%. The trained dense network has an error of 1.65%. This result is obtained at η = 2.
11https://pytorch.org/hub/pytorch_vision_resnet/ 12https://github.com/RAIVNLab/STR 13https://pytorch.org/hub/pytorch_vision_mobilenet_v2/ 14https://github.com/NVIDIA/DeepLearningExamples/blob/master/PyTorch/
SpeechRecognition/Jasper/configs/jasper10x5dr_sp_offline_specaugment.toml 15https://github.com/NVIDIA/DeepLearningExamples/blob/master/PyTorch/ LanguageModeling/Transformer-XL/pytorch/wt103_base.yaml 16https://github.com/ElementAI/N-BEATS/blob/master/experiments/m4/ interpretable.gin
The η parameter in our algorithm controls the accuracy tolerance. The higher, the more tolerant we are on the loss of precision and the sparser the network is. Thus, this parameter also controls the network sparsity. The choice of this parameter should be the result of an accuracy-sparsity trade-off. This is illustrated in Figure 2.

Review summaries:
The paper propose a network compression algorithm by exploiting a reformulation of activation function as proximity operator. The latter is an optimization problem whose optimality condition reveals constraints on the weight matrix W of the neural net. The main idea is then to "biasedly" select W as a minimizer of a sparsity inducing penalties under a relaxation of the previous optimality conditions. The authors provide details on solving such problem as well as numerical experiments that leads to similar results than competitors.

In this paper the authors propose a new model compression method based on subdifferential inclusion. The key idea is to make the outputs of the neurons in the sparse and dense networks at the same input close enough. They rewrite the activation function as the proximity operator of a proper convex function and finally formulate the compression problem into a constraint minimization problem using the technique of subdifferential inclusion. They conduct a series of experiments to evaluate the performance of their proposed methods.