Training Structured Neural Networks Through Manifold Identification and Vari- ance Reduction
Year: 2022

Abstract:
This paper proposes an algorithm, RMDA, for training neural networks (NNs) with a regularization term for promoting desired structures. RMDA does not incur computation additional to proximal SGD with momentum, and achieves variance reduction without requiring the objective function to be of the finite-sum form. Through the tool of manifold identification from nonlinear optimization, we prove that after a finite number of iterations, all iterates of RMDA possess a desired structure identical to that induced by the regularizer at the stationary point of asymptotic convergence, even in the presence of engineering tricks like data augmentation that complicate the training process. Experiments on training NNs with structured sparsity confirm that variance reduction is necessary for such an identification, and show that RMDA thus significantly outperforms existing methods for this task. For unstructured sparsity, RMDA also outperforms a state-of-the-art pruning method, validating the benefits of training structured NNs through regularization. Implementation of RMDA is available at https://www.github.com/zihsyuan1214/rmda.

heading: 1 Introduction
text: Training neural networks (NNs) with regularization to obtain a certain desired structure such as structured sparsity or discrete-valued parameters is a problem of increasing interest. Existing approaches either use stochastic subgradients of the regularized objective (Wen et al., 2016; 2018) or combine popular stochastic gradient algorithms for NNs, like SGD with momentum (MSGD) or Adam (Kingma & Ba, 2015), with the proximal operator associated with the regularizer to conduct proximal stochastic gradient updates to obtain a model with preferred structures (Bai et al., 2018; Yang et al., 2019; Yun et al., 2020; Deleu & Bengio, 2021). Such methods come with proven convergence for certain measures of first-order optimality and have shown some empirical success in applications. However, we notice that an essential theoretical support lacking in existing methods is the guarantee for the output iterate to possess the same structure as that at the point of convergence. More specifically, often the imposed regularization is only known to induce a desired structure exactly at optimal or stationary points of the underlying optimization problem (see for example, Zhao & Yu, 2006), but training algorithms are only able to generate iterates asymptotically converging to a stationary point. Without further theoretical guarantees, it is unknown whether the output iterate, which is just an approximation of the stationary point, still has the same structure. For example, let us assume that sparsity is desired, the point of convergence is x∗ = (1, 0, 0), and two algorithms respectively produce iterates {yt = (1, t−1, t−1)} and {zt = (1 + t−1, 0, 0)}. Clearly, both iterate sequences converge to x∗, but only zt has the same desired structure as its limit point x∗, while yt is not useful for sparsity despite that the point of convergence is. This work aims at filling this gap to propose an algorithm for training structured NNs that can provably make all its iterates after
a finite number of iterations possess the desired structure of the stationary point to which the iterates converge. We term the structure at a stationary point a stationary structure, and it should be understood that for multiple stationary points, each might correspond to a different stationary structure, and we aim at identifying the one at the limit point of the iterates of an algorithm, instead of selecting the optimal one among all stationary structures. Although finding the structure at an inferior stationary point might seem not very meaningful, another reason for studying this identification property is that for the same point of convergence, the structure at the limit point is the most preferable one. Consider the same example above, we note that for any sequence {xt} converging to x∗, xt1 6= 0 for all t large enough, for otherwise xt does not converge to x∗. Therefore, xt cannot be sparser than x∗ if xt → x∗.1 Identifying the structure of the point of convergence thus also amounts to finding the locally most ideal structure under the same convergence premise. It is well-known in the literature of nonlinear optimization that generating iterates consistently possessing the structure at the stationary point of convergence is possible if all points with the same structure near the stationary point can be presented locally as a manifold along which the regularizer is smooth. This manifold is often termed as the active manifold (relative to the given stationary point), and the task of generating iterates staying in the active manifold relative to the point of convergence after finite iterations is called manifold identification (Lewis, 2002; Hare & Lewis, 2004; Lewis & Zhang, 2013). To identify the active manifold of a stationary point, we need the regularizer to be partly smooth (Lewis, 2002; Hare & Lewis, 2004) at that point, roughly meaning that the regularizer is smooth along the active manifold around the point, while the change in its value is drastic along directions leaving the manifold. A more technical definition will be given in Section 3. Fortunately, most regularizers used in machine learning are partly smooth, so stationary structure identification is possible, and various deterministic algorithms are known to achieve so (Hare & Lewis, 2007; Hare, 2011; Wright, 2012; Liang et al., 2017a;b; Li et al., 2020; Lee, 2020; Bareilles et al., 2020). On the other hand, for stochastic gradient-related methods to identify a stationary structure, existing theory suggests that the variance of the gradient estimation needs to vanish as the iterates approach a stationary point (Poon et al., 2018), and indeed, it is observed empirically that proximal stochastic gradient descent (SGD) is incapable of manifold identification due to the presence of the variance in the gradient estimation (Lee & Wright, 2012; Sun et al., 2019).2 Poon et al. (2018) showed that variance-reduction methods such as SVRG (Johnson & Zhang, 2013; Xiao & Zhang, 2014) and SAGA (Defazio et al., 2014) that utilize the finite-sum structure of empirical risk minimization to drive the variance of their gradient estimators to zero are suitable for this task. Unfortunately, with the standard practice of data augmentation in deep learning, training of deep learning models with a regularizer should be treated as the following stochastic optimization problem that minimizes the expected loss over a distribution, instead of the commonly seen finite-sum form:
min W∈E F (W ) := Eξ∼D [fξ (W )] + ψ (W ) , (1)
where E is a Euclidean space with inner product 〈·, ·〉 and the associated norm ‖·‖, D is a distribution over a space Ω, fξ is differentiable almost everywhere for all ξ ∈ Ω, and ψ(W ) is a regularizer that might be nondifferentiable. We will also use the notation f(W ) := Eξ∼D[fξ(W )]. Without a finite-sum structure in (1), Defazio & Bottou (2019) pointed out that classical variance-reduction methods are ineffective for deep learning, and one major reason is the necessity of periodically evaluating ∇f(W ) (or at least using a large batch from D to get a precise approximation of it) in variance-reduction methods is intractable, hence manifold identification and therefore finding the stationary structure becomes an extremely tough task for deep learning. Although recently there are efforts in developing variancereduction methods for (1) inspired by online problems (Wang et al., 2019; Nguyen et al., 2021; Pham et al., 2020; Cutkosky & Orabona, 2019; Tran-Dinh et al., 2019), these methods all have multiple hyperparameters to tune and incur computational cost at least twice or
1See a more detailed discussion in Appendix B.1. 2An exception is the interpolation case, in which the variance of plain SGD vanishes asymptot-
ically. But data augmentation often fails this interpolation condition.
thrice to that of (proximal) SGD. As the training of deep learning models is time- and resource-consuming, these drawbacks make such methods less ideal for deep learning. To tackle these difficulties, we extend the recently proposed modernized dual averaging framework (Jelassi & Defazio, 2020) to the regularized setting by incorporating proximal operations, and obtain a new algorithm RMDA (Regularized Modernized Dual Averaging) for (1). The proposed algorithm provably achieves variance reduction beyond finite-sum problems without any cost or hard-to-tune hyperparameters additional to those of proximal momentum SGD (proxMSGD), and we provide theoretical guarantees for its convergence and ability for manifold identification. The key difference between RMDA and the original regularized dual averaging (RDA) of Xiao (2010) is that RMDA incorporates momentum and can achieve better performance for deep learning in terms of the generalization ability, and the new algorithm requires nontrivial proofs for its guarantees. We further conduct experiments on training deep learning models with a regularizer for structured-sparsity to demonstrate the ability of RMDA to identify the stationary structure without sacrificing the prediction accuracy. When the desired structure is (unstructured) sparsity, a popular approach is pruning that trims a given dense model to a specified level, and works like (Gale et al., 2019; Blalock et al., 2020; Evci et al., 2020; Verma & Pesquet, 2021) have shown promising results. However, as a post-processing approach, pruning is essentially different from structured training considered in this work, because pruning is mainly used when a model is available, while structured training combines training and structure inducing in one procedure to potentially reduce the computational cost and memory footprint when resources are scarce. We will also show in our experiment that RMDA can achieve better performance than a state-of-the-art pruning method, suggesting that structured training indeed has its merits for obtaining sparse NNs. The main contributions of this work are summarized as follows. • Principled analysis: We use the theory of manifold identification from nonlinear opti-
mization to provide a unified way towards better understanding of algorithms for training structured neural networks.
• Variance reduction beyond finite-sum with low cost: RMDA achieves variance reduction for problems that consist of an infinite-sum term plus a regularizer (see Lemma 2) while incorporating momentum to improve the generalization performance. Its spatial and computational cost is almost the same as proxMSGD, and there is no additional hyperparameters to tune, making RMDA suitable for large-scale deep learning.
• Structure identification: With the help of variance reduction, our theory shows that under suitable conditions, after a finite number of iterations, iterates of RMDA stay in the active manifold of its limit point.
• Superior empirical performance: Experiments on neural networks with structured sparsity exemplify that RMDA can identify a stationary structure without reducing the validation accuracy, thus outperforming existing methods by achieving higher group sparsity. Another experiment on unstructured sparsity also shows RMDA outperforms a state-of-the-art pruning method.
After this work is finished, we found a very recent paper Kungurtsev & Shikhman (2021) that proposed the same algorithm (with slightly differences in the parameters setting in Line 5 of Algorithm 1) and analyzed the expected convergence of (1) under a specific scheduling of ct = st+1α−1t+1 when both terms are convex. In conrast, our work focuses on nonconvex deep learning problems, and especially on the manifold identification aspect.

heading: 2 Algorithm
text: Details of the proposed RMDA are in Algorithm 1. At the t-th iteration with the iterate W t−1, we draw an independent sample ξt ∼ D to compute the stochastic gradient ∇fξt(W t−1), decide a learning rate ηt, and update the weighted sum Vt of previous stochastic gradients using ηt and the scaling factor βt := √ t:
V0 := 0, Vt := ∑t
k=1 ηkβk∇fξk(W k−1) = Vt−1 + ηtβt∇fξt(W t−1), ∀t > 0.
Algorithm 1: RMDA (W 0, T, η(·), c(·)) input : Initial point W 0, learning rate schedule η(·), momentum schedule c(·),
number of epochs T 1 V0 ← 0, α0 ← 0 2 for t = 1, . . . , T do 3 βt ← √ t, st ← η(t)βt, αt ← αt−1 + st 4 Sample ξt ∼ D and compute V t ← V t−1 + st∇fξt(W t−1) 5 W̃ t ← arg minW 〈V t, W 〉+ βt2
∥∥W −W 0∥∥2 + αtψ(W ) // (2) 6 W t ← (1− c(t))W t−1 + c(t)W̃ t
output: The final model WT
The tentative iterate W̃ t is then obtained by the proximal operation associated with ψ: W̃ t = proxαtβ−1t ψ ( W 0 − β−1t V t ) , αt := ∑t k=1 βkηk, (2)
where for any function g, proxg(x) := arg miny ‖y − x‖ 2 /2 + g(y) is its proximal operator. The iterate is then updated along the direction W̃ t −W t−1 with a factor of ct ∈ [0, 1]: W t = (1− ct)W t−1 + ctW̃ t = W t−1 + ct ( W̃ t −W t−1 ) . (3)
When ψ ≡ 0, RMDA reduces to the modernized dual averaging algorithm of Jelassi & Defazio (2020), in which case it has been shown that mixing W t−1 and W̃ t in (3) equals to introducing momentum (Jelassi & Defazio, 2020; Tao et al., 2018). We found that this introduction of momentum greatly improves the performance of RMDA and is therefore essential for applying it on deep learning problems.

heading: 3 Analysis
text: We provide theoretical analysis of the proposed RMDA in this section. Our analysis shows variance reduction in RMDA and stationarity of the limit point of its iterates, but all of them revolves around our main purpose of identification of a stationary structure within a finite number of iterations. The key tools for this end are partial smoothness and manifold identification (Hare & Lewis, 2004; Lewis, 2002). Our result is the currently missing cornerstone for those proximal algorithms applied to deep learning problems for identifying desired structures. In fact, it is actually well-known in convex optimization that those algorithms based on plain proximal stochastic gradient without variance reduction are unable to identify the active manifold, and the structure of the iterates oscillates due to the variance in the gradient estimation; see, for example, experiments and discussions in Lee & Wright (2012); Sun et al. (2019). Our work is therefore the first one to provide justification for solving the regularized optimization problem in deep learning to really identify a desired structure induced by the regularizer. Throughout, ∇fξ denotes the gradient of fξ, ∂ψ is the (regular) subdifferential of ψ, and relint(C) means the relative interior of the set C. We start from introducing the notion of partial smoothness. Definition 1. A function ψ is partly smooth at a point W ∗ relative to a set MW∗ 3W ∗ if 1. Around W ∗, MW∗ is a C2-manifold and ψ|MW∗ is C2. 2. ψ is regular (finite with the Fréchet subdifferential coincides with the limiting Fréchet
subdifferential) at all points W ∈MW∗ around W ∗ with ∂ψ(W ) 6= ∅. 3. The affine span of ∂ψ(W ∗) is a translate of the normal space to MW∗ at W ∗. 4. ∂ψ is continuous at W ∗ relative to MW∗ .
We often call MW∗ the active manifold at W ∗. Another concept required for manifold identification is prox-regularity (Poliquin & Rockafellar, 1996). Definition 2. A function ψ is prox-regular at W ∗ for V ∗ ∈ ∂ψ(W ∗) if ψ is finite at W ∗, locally lower semi-continuous around W ∗, and there is ρ > 0 such that ψ(W1) ≥ ψ(W2) + 〈V, W1−W2〉− ρ2‖W1 −W2‖
2 whenever W1,W2 are close to W ∗ with ψ(W2) near ψ(W ∗) and V ∈ ∂ψ(W2) near V ∗. ψ is prox-regular at W ∗ if it is so for all V ∈ ∂ψ(W ∗).
To broaden the applicable range, a function ψ prox-regular at some W ∗ is often also assumed to be subdifferentially continuous (Poliquin & Rockafellar, 1996) there, meaning that if W t → W ∗, ψ(W t) → ψ(W ∗) holds when there are V ∗ ∈ ∂ψ(W ∗) and a sequence {V t} such that V t ∈ ∂ψ(W t) and V t → V ∗. Notably, all convex and weakly-convex (Nurminskii, 1973) functions are regular, prox-regular, and subdifferentially continuous in their domain.

heading: 3.1 Theoretical Results
text: When the problem is convex, convergence guarantees for Algorithm 1 under two specific specific schemes are known. First, when ct ≡ 1, RMDA reduces to the classical RDA, and convergence to a global optimum (of W t = W̃ t in this case) on convex problems has been proven by Lee & Wright (2012); Duchi & Ruan (2021), with convergence rates of the expected objective or the regret given by Xiao (2010); Lee & Wright (2012). Second, when ct = st+1α−1t+1 and (βt, αt) in Line 5 of Algorithm 1 are replaced by (βt+1, αt+1), convergence is recently analyzed by Kungurtsev & Shikhman (2021). In our analysis below, we do not assume convexity of either term. We show that if {W̃ t} converges to a point W ∗ (which could be a non-stationary one), {W t} also converges to W ∗. Lemma 1. Consider Algorithm 1 with {ct} satisfying ∑ ct = ∞. If {W̃ t} converges to a point W ∗, {W t} also converges to W ∗.
We then show that if {W̃ t} converges to a point, almost surely this point of convergence is stationary. This requires the following lemma for variance reduction of RMDA, meaning that the variance of using Vt to estimate ∇f(W t−1) reduces to zero, as α−1t Vt converges to ∇f(W t−1) almost surely, and this result could be of its own interest. The first claim below uses a classical result in stochastic optimization that can be found at, for example, (Gupal, 1979, Theorem 4.1, Chapter 2.4), but the second one is, to our knowledge, new. Lemma 2. Consider Algorithm 1. Assume for any ξ ∼ D, fξ is L-Lipschitz-continuouslydifferentiable almost surely for some L, so f is also L-Lipschitz-continuously-differentiable, and there is C ≥ 0 such that Eξt∼D
∥∥∇fξt (W t−1)∥∥2 ≤ C for all t. If {ηt} satisfies∑ βtηtα −1 t =∞, ∑( βtηtα −1 t )2 <∞,
∥∥W t+1 −W t∥∥ (βtηtα−1t )−1 a.s.−−→ 0, (4) then α−1t V t −→ ∇f(W t−1) with probability one. Moreover, if {W t} lies in a bounded set, we get E
∥∥α−1t V t −∇f (W t−1)∥∥2 → 0 even if the second condition in (4) is replaced by a weaker condition of βtηtα−1t → 0.
In general, the last condition in (4) requires some regularity conditions in F to control the change speed of W t. One possibility is when ψ is the indicator function of a convex set, βtηt ∝ tp for t ∈ (1/2, 1) will satisfy this condition. However, in other settings for ηt, even when F and ψ are both convex, existing analyses for the classical RDA such that ct ≡ 1 in Algorithm 1 still need an additional local error bound assumption to control the change of W t+1 −W t. Hence, to stay focused on our main message, we take this assumption for granted, and leave finding suitable sufficient conditions for it as future work. With the help of Lemmas 1 and 2, we can now show the stationarity result for the limit point of the iterates. The assumption of βtα−1t approaching 0 below is classical in analyses of dual averaging in order to gradually remove the influence of the term
∥∥W −W 0∥∥2. Theorem 1. Consider Algorithm 1 with the conditions in Lemmas 1 and 2 hold, and assume the set of stationary points Z := {W | 0 ∈ ∂F (W )} is nonempty and βtα−1t → 0. For any given W 0, consider the event that {W̃ t} converges to a point W ∗ (each event corresponds to a different W ∗), then if ∂ψ is outer semicontinuous at W ∗, and this event has a nonzero probability, W ∗ ∈ Z, or equivalently, W ∗ is a stationary point, with probability one conditional on this event.
Finally, with Lemmas 1 and 2 and Theorem 1, we prove the main result that the active manifold of the limit point is identified in finite iterations of RMDA under nondegeneracy.
Theorem 2. Consider Algorithm 1 with the conditions in Theorem 1 satisfied. Consider the event of {W̃ t} converging to a certain point W ∗ as in Theorem 1, if the probability of this event is nonzero; ψ is prox-regular and subdifferentially continuous at W ∗ and partly smooth at W ∗ relative to the active C2 manifoldM; ∂ψ is outer semicontinuous at W ∗; and the nondegeneracy condition
−∇f (W ∗) ∈ relint ∂ψ (W ∗) (5) holds at W ∗, then conditional on this event, almost surely there is T0 ≥ 0 such that
W̃ t ∈M, ∀t ≥ T0. (6)

heading: In other words, the active manifold at W ∗ is identified by the iterates of Algorithm 1 after a finite number of iterations almost surely.
text: As mentioned in Section 1, an important reason for studying manifold identification is to get the lowest-dimensional manifold representing the structure of the limit point, which often corresponds to a preferred property for the application, like the highest sparsity, lowest rank, or lowest VC dimension locally. See an illustrated example in Appendix B.1.

heading: 4 Applications in Deep Learning
text: We discuss two popular schemes of training structured deep learning models achieved through regularization to demonstrate the applications of RMDA. More technical details for applying our theory to the regularizers in these applications are in Appendix B.

heading: 4.1 Structured Sparsity
text: As modern deep NN models are often gigantic, it is sometimes desirable to trim the model to a smaller one when only limited resources are available. In this case, zeroing out redundant parameters during training at the group level is shown to be useful (Zhou et al., 2016), and one can utilize regularizers promoting structured sparsity for this purpose. The most famous regularizer of this kind is the group-LASSO norm (Yuan & Lin, 2006; Friedman et al., 2010). Given λ ≥ 0 and a collection G of index sets {Ig} of the variable W , this convex regularizer is defined as
ψ(W ) := λ ∑|G| g=1 wg ∥∥WIg∥∥, (7)
with wg > 0 being the pre-specified weight for Ig. For any W ∗, let GW∗ ⊆ G be the index set such that W ∗Ij = 0 for all j ∈ GW∗ , the group-LASSO norm is partly smooth around W ∗ relative to the manifold MW∗ := {W |WIi = 0,∀i ∈ GW∗}, so our theory applies. In order to promote structured sparsity, we need to carefully design the grouping. Fortunately, in NNs, the parameters can be grouped naturally (Wen et al., 2016). For any fully-connected layer, let W ∈ Rout×in be the matrix representation of the associated parameters, where out is the number of output neurons and in is that of input neurons, we can consider the column-wise groups, defined as W:,j for all j, and the row-wise groups of the form Wi,:. For a convolutional layer with W ∈ Rfilter×channel×height×width being the tensor form of the corresponding parameters, we can consider channel-wise, filter-wise, and kernel-wise groups, defined respectively as W:,j,:,:, Wi,:,:,: and Wi,j,:,:.

heading: 4.2 Binary Neural Networks
text: Making the parameters of an NN binary integers is another way to obtain a more compact model during training and deployment (Hubara et al., 2016), but discrete optimization is hard to scale-up. Using a vector representation w ∈ Rm of the variables, Hou et al. (2017) thus proposed to use the indicator function of { w | wIi = αibIi , αi > 0, bIi ∈ {±1}|Ii| } to induce the entries of w to be binary without resorting to discrete optimization tools, where each Ii enumerates all parameters in the i-th layer. Yang et al. (2019) later proposed to use minα∈[0,1]m ∑m i=1 ( αi(wi + 1)2 + (1− αi)(wi − 1)2 ) as the regularizer and to include α as
a variable to train. At any α∗ with I0 := {i | α∗i = 0} and I1 := {i | α∗i = 1}, the objective is partly smooth relative to the manifold {(W,α) | αI0 = 0, αI1 = 1}. Extension to discrete NNs beyond the binary ones is possible, and Bai et al. (2018) have proposed regularizers with closed-form proximal operators for it.

heading: 5 Experiments
text: We use the structured sparsity application in Section 4.1 to empirically exemplify the ability of RMDA to find desired structures in the trained NNs. RMDA and the following methods for structured sparsity in deep learning are compared using PyTorch (Paszke et al., 2019). • ProxSGD (Yang et al., 2019): A simple proxMSGD algorithm. To obtain group sparsity,
we skip the interpolating step in Yang et al. (2019). • ProxSSI (Deleu & Bengio, 2021): This is a special case of the adaptive proximal SGD
framework of Yun et al. (2020) that uses the Newton-Raphson algorithm to approximately solve the subproblem. We directly use the package released by the authors.
We exclude the algorithm of Wen et al. (2016) because their method is shown to be worse than ProxSSI by Deleu & Bengio (2021). To compare these algorithms, we examine both the validation accuracy and the group sparsity level of their trained models. We compute the group sparsity as the percentage of groups whose elements are all zero, so the reported group sparsity is zero when there is no group with a zero norm, and is one when the whole model is zero. For all methods above, we use (7) with column-wise and channel-wise groupings in the regularization for training, but adopt the kernel-wise grouping in their group sparsity evaluation. Throughout the experiments, we always use multi-step learning rate scheduling that decays the learning rate by a constant factor every time the epoch count reaches a pre-specified threshold. For all methods, we conduct grid searches to find the best hyperparameters. All results shown in tables are the mean and standard deviation of three independent runs with the same hyperparameters, while figures use one representative run for better visualization. In convex optimization, a popular way to improve the practical convergence behavior for momentum-based methods is restarting that periodically reset the momentum to zero (O’donoghue & Candes, 2015). Following this idea, we introduce a restart heuristic to RMDA. At each round, we use the output of Algorithm 1 from the previous round as the new input to the same algorithm, and continue using the scheduling η and c without resetting them. For ψ ≡ 0, Jelassi & Defazio (2020) suggested to increase ct proportional to the decrease of ηt until reaching ct = 1. We adopt the same setting for ct and ηt and restart RMDA whenever ηt changes. As shown in Section 3 that W̃ t finds the active manifold, increasing ct to 1 also accords with our interest in identifying the stationary structure.
5.1 Correctness of Identified Structure Using Synthetic Data
Our first step is to numerically verify that RMDA can indeed identify the stationary structure desired. To exactly find a stationary point and its structure a priori, we consider synthetic problems. We first decide a ground truth model W that is structured sparse, generate random data points that can be well separated by W , and then decide their labels using W . The generated data are then taken as our training data. We consider a linear logistic regression model and a small NN that has one fully-connected layer and one convolutional
layer. To ensure convergence to the ground truth, for logistic regression we generate more data points than the problem dimension to ensure the problem is strongly convex so that there is only one stationary/optimal point, and for the small NN, we initialize all algorithms close enough to the ground truth. We report in Fig. 1 training error rates (as an indicator for the proximity to the ground truth) and percentages of the optimal group sparsity pattern of the ground truth identified. Clearly, although all methods converge to the ground truth, only RMDA identifies the correct structure of it, and other methods without guarantees for manifold identification fail.

heading: 5.2 Linear Model and Fully-Connected Neural Network
text: We turn to real-world data and begin with two relatively simple problems. The first is logistic regression with the MNIST dataset (LeCun et al., 1998), and the second is training a fully-connected NN on the FashionMNIST data (Xiao et al., 2017). From Figs. 2a and 2b, we see that similar to the previous experiment, the group sparsity level of RMDA is stable in the last epochs, while that of ProxSGD and ProxSSI oscillates below. This suggests that RMDA is the only method that identifies the structured sparsity at its limit point, and other methods with no variance reduction fail. Moreover, Table 1a shows that manifold identification of RMDA is achieved with no sacrifice of the validation accuracy.

heading: 5.3 Convolutional Neural Network
text: We now move to CNNs to train a version of LeNet5 on MNIST and a modified VGG16 (Simonyan & Zisserman, 2015) with the CIFAR10 dataset (Krizhevsky, 2009). For VGG16, we follow the practice in modern vision tasks to apply data augmentation through random cropping and horizontal flipping so that the training problem is no longer a finite-sum one. The group sparsity results are in Figs. 2c and 3. Notably, for VGG16, RMDA succeeds in finding the structured sparsity pattern despite the presence of data augmentation, showing that RMDA can indeed overcome the difficulty from the infinite-sum setting of modern deep learning tasks, but other methods all have group sparsity levels fluctuating way below that of RMDA. For this more complicated task, we include a dense baseline of MSGD with no sparsity-inducing regularizer in our comparison in Table 1b. Remarkably, RMDA beats ProxSGD and ProxSSI in both structured sparsity and accuracy, and its accuracy is close to that of MSGD. We also report that in the VGG16 task, on our NVIDIA RTX 8000 GPU, MSGD, ProxSGD, and RMDA all have an average per-epoch cost of around 39 seconds, while ProxSSI needs 189 seconds per epoch. RMDA is thus also more suitable for large-scale structured deep learning in terms of practical efficiency. To test RMDA on more models, we further train a residual NN (ResNet50) with CIFAR100 and the same data augmentation above. This is a more difficult task and ProxSSI crashes after 1, 170 epochs, so we report its best result before crashing. Table 1c shows that RMDA still outperforms others in both metrics, and even attains the same accuracy as MSGD.

heading: 5.4 Comparison with Pruning
text: We compare RMDA with a state-of-the-art pruning method RigL (Evci et al., 2020). As pruning focuses on unstructured sparsity, we use RMDA with ψ(W ) = λ‖W‖1 to have a fair comparison, and tune λ to achieve a pre-specified sparsity level. We run RigL with 1000 epochs, as its performance at the default 500 epochs was unstable, and let RMDA use the same number of epochs. Results of 98% sparsity in Table 2 show that RMDA consistently outdoes RigL, indicating regularized training could be a promising alternative to pruning.

heading: 6 Conclusions
text: In this work, we proposed and analyzed a new algorithm, RMDA, for efficiently training structured neural networks with state-of-the-art performance. Even in the presence of data augmentation, RMDA can still achieve variance reduction and provably identify the desired structure at a stationary point using the tools of manifold identification. Experiments show that existing algorithms for the same purpose fail to find a stable stationary structure, while RMDA achieves so with no accuracy drop nor additional time cost.

heading: Acknowledgements
text: This work was supported in part by MOST of R.O.C. grant 109-2222-E-001-003-MY3, and the AWS Cloud Credits for Research program of Amazon Inc.

heading: Appendices

heading: Table of Contents

heading: A Proofs 13
text: A.1 Proof of Lemma 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 A.2 Proof of Lemma 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 A.3 Proof of Theorem 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 A.4 Proof of Theorem 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

heading: B Additional Discussions on Applications 18
text: B.1 Structured Sparsity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 B.2 Binary Neural Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

heading: C Experiment Setting Details 20

heading: D More Details of the Experiments 21
text: E Other Regularizers for Possibly Better Group Sparsity and Generalization 21

heading: A Proofs

heading: A.1 Proof of Lemma 1
text: Proof. Using (3), the distance between W t and W ∗ can be upper bounded through the triangle inequality:∥∥W t −W ∗∥∥ = ∥∥(1− ct) (W t−1 −W ∗)+ ct (W̃ t −W ∗)∥∥
≤ ct ∥∥W̃ t −W ∗∥∥+ (1− ct)∥∥W t−1 −W ∗∥∥. (8)
For any event such that W̃ t →W ∗, for any > 0, we can find T ≥ 0 such that∥∥W̃ t −W ∗∥∥ ≤ , ∀t ≥ T . Let δt := ‖W t −W ∗‖, we see from the above inequality and (8) that
δt ≤ (1− ct) δt−1 + ct , ∀t ≥ T .
By deducting from both sides, we get that
(δt − ) ≤ (1− ct) (δt−1 − ) , ∀t ≥ T . Since ∑ ct =∞, we further deduce that
lim t→∞ (δt − ) ≤ ∞∏ t=T (1− ct) (δT −1 − )
≤ ∞∏ t=T exp (−ct) (δT −1 − )
= exp ( − ∞∑ t=T ct ) (δT −1 − ) = 0,
where in the first inequality we used the fact that exp(x) ≥ 1 +x for all real number x. The result above then implies that
lim t→∞ δt ≤ .
As is arbitrary and δt ≥ 0 from the definition, we conclude that limt→∞ δt = 0, which is equivalent to that W t →W ∗.

heading: A.2 Proof of Lemma 2
text: Proof. We observe that
α−1t V t = t∑ k=1 ηkβk αt ∇fξk ( W k−1 ) = αt−1
αt α−1t−1V t−1 + αt − αt−1 αt
∇fξt ( W t−1 ) = (
1− βtηt αt
) α−1t−1V
t−1 + βtηt αt ∇fξt
( W t−1 ) .
From that f is L-Lipschitz-continuously differentiable, we have that∥∥Eξt+1∼D [∇fξt+1 (W t)]− Eξt∼D [∇fξt (W t−1)]∥∥ = ∥∥∇f (W t)− f (W t−1)∥∥ ≤ L
∥∥W t −W t−1∥∥. (9) Therefore, (4) and (9) imply that
0 ≤ ∥∥Eξt+1∼D [∇fξt+1 (W t)]− Eξt∼D [∇fξt (W t−1)]∥∥
βtηtα −1 t
≤ L ∥∥W t −W t−1∥∥
βtηtα −1 t
a.s.−−→ 0,
which together with the sandwich lemma shows that∥∥Eξt+1∼D [∇fξt+1 (W t)]− Eξt∼D [∇fξt (W t−1)]∥∥ βtηtα −1 t a.s.−−→ 0. (10)
Therefore, the first two conditions of (4) together with (10) and the bounded variance assumption satisfy the requirements of (Gupal, 1979, Chapter 2.4, Theorem 4.1), so the conclusion of almost sure convergence hold.
For the convergence in L2 part, we first define mt := α−1t Vt and τt := βtηtα−1t for notational ease. Consider ∥∥mt+1 −∇F (W t)∥∥2, we have from the update rule in Algorithm 1 that∥∥mt+1 −∇F (W t)∥∥2 = ∥∥(1− τt)mt + τt∇fξt+1(W t)−∇F (W t)∥∥2
= ∥∥(1− τt) (mt −∇F (W t))+ τt (∇fξt+1(W t)−∇F (W t))∥∥2
= (1− τt)2 ∥∥mt −∇F (W t)∥∥2 + τ2t ∥∥∇fξt+1(W t)−∇F (W t)∥∥2
+ 2τt(1− τt)〈mt −∇F (W t), ∇fξt+1(W t)−∇F (W t)〉 = (1− τt)2 ∥∥(mt −∇F (W t−1))+ (∇F (W t−1)−∇F (W t))∥∥2 (11)
+ τ2t ∥∥∇fξt+1(W t)−∇F (W t)∥∥2 + 2τt(1− τt)〈mt −∇F (W t), ∇fξt+1(W t)−∇F (W t)〉.
Let {Ft}t≥0 denote the natural filtration of {(mt,W t)}t≥0. Namely, Ft records the information of W 0, {ci}t−1i=0, {ηi}t−1i=0, and {ξi}ti=1. By defining Ut := ∥∥mt −∇F (W t−1)∥∥2 and taking expectation over (11) conditional on Ft, we obtain from E [ ∇fξt+1(W t) | Ft ] = ∇F (W t) that
E [Ut+1 | Ft] = (1− τt)2 ∥∥(mt −∇F (W t−1))+ (∇F (W t−1)−∇F (W t))∥∥2
+ τ2t E [∥∥∇fξt(W t)−∇F (W t)∥∥2 | Ft] . (12)
From the last condition in (4) and the Lipschitz continuity of∇F , there are random variables { t} and {ut} such that ‖ut‖ = 1, t ≥ 0, and ∇F (W t−1)−∇F (W t) = τt tut for all t > 0, with t ↓ 0 almost surely. We thus obtain that∥∥mt −∇F (W t−1) +∇F (W t−1)−∇F (W t)∥∥2
= ∥∥mt −∇F (W t−1) + τt tut∥∥2
= (1 + τt)2 ∥∥∥∥ 11 + τt (mt −∇F (W t−1))+ τt1 + τt tut ∥∥∥∥2 ≤ (1 + τt)2 ( 1
1 + τt Ut + τt 1 + τt
t 2 ) , (13)
where we used Jensen’s inequality and the convexity of ‖·‖2 in the last inequality. By substituting (13) back into (12), we obtain
E [Ut+1|Ft] ≤ (1− τt)2(1 + τt)Ut + (1− τt)2(1 + τt)τt t2 + τt2E [∥∥∇fξt(W t)−∇F (W t)∥∥2 | Ft]
≤ (1− τt)(Ut + τt t2) + τt2E [∥∥∇fξt(W t)−∇F (W t)∥∥2 | Ft]
≤ (1− τt)Ut + τt t2 + τt2E [∥∥∇fξt(W t)−∇F (W t)∥∥2 | Ft] . (14)
For the last term in (14), we notice that E [∥∥∇fξt(W t)−∇F (W t)∥∥2 | Ft] ≤ 2(E [∥∥∇fξt(W t)∥∥2]+ ∥∥∇F (W t)∥∥2)
≤ 2 ( C + ∥∥∇F (W t)∥∥2) , (15) where the last inequality is from the bounded variance assumption. Since by assumption the {W t} lies in a bounded set K, we have that for any point W ∗ ∈ K, W t −W ∗ is upper bounded, and thus ‖∇F (W t)−∇F (W ∗)‖ is also bounded, implying that ‖∇F (W t)‖2 ≤ C2 for some C2 ≥ 0. Therefore, (15) further leads to
E [∥∥∇fξt(W t)−∇F (W t)∥∥2 | Ft] ≤ C3 (16)
for some C3 ≥ 0. Now we further take expectation on (14) and apply (16) to obtain
EUt+1 ≤ (1− τt)EUt + τt t2 + τt2C3 = (1− τt)EUt + τt ( 2t + τtC3 ) . (17)
Note that the third implies t ↓ 0, so this together with the second condition that τt ↓ 0 means 2t+τtC3 ↓ 0 as well, and thus for any δ > 0, we can find Tδ ≥ 0 such that 2t+τtC3 ≤ δ for all t ≥ Tδ. Thus, (17) further leads to
EUt+1 − δ ≤ (1− τt)EUt + τtδ − δ = (1− τt) (EUt − δ) ,∀t ≥ Tδ. (18) This implies that (EUt − δ) becomes a decreasing sequence starting from t ≥ Tδ, and since Ut ≥ 0, this sequence is lower bounded by −δ, and hence it converges to a certain value. By recursion of (18), we have that
EUt − δ ≤ t∏
i=Tδ
(1− τi) (EUTδ − δ) ,
and from the well-known inequality (1 + x) ≤ expx for all x ∈ R, the above result leads to EUt − δ ≤ exp ( − ∑ i = Tδtτi ) (EUTδ − δ) .
By letting t approach infinity and noting that the first condition of (4) indicates ∞∑ t=k τt =∞
for any k ≥ 0, we see that
−δ ≤ lim t→∞ EUt − δ ≤ exp ( − ∞∑ i=Tδ τi ) (EUTδ − δ) = 0. (19)
As δ is arbitrary, by taking δ ↓ 0 in (19) and noting the nonnegativity of Ut, we conclude that limEUt = 0, as desired. This proves the last result in Lemma 2.

heading: A.3 Proof of Theorem 1
text: Proof. Using Lemma 2, we can view α−1t V t as ∇f(W t) plus some noise that asymptotically decreases to zero with probability one:
α−1t Vt = ∇f(W t) + t, ‖ t‖ a.s.−−→ 0. (20)
We use (20) to rewrite the optimality condition of (2) as (also see Line 5 of Algorithm 1) − ( ∇f ( W t ) + t + βtα−1t ( W̃ t −W 0 )) ∈ ∂ψ ( W̃ t ) . (21)
Now we consider ∂F (W̃ t). Clearly from (21), we have that ∇f ( W̃ t ) −∇f ( W t ) − t − βtα−1t ( W̃ t −W 0 ) ∈ ∂∇f ( W̃ t ) + ψ ( W̃ t ) = ∂F ( W̃ t ) . (22)
Now we consider the said event that W̃ t → W ∗ for a certain W ∗, and let us define this event as A ⊆ Ω. From Lemma 1, we know that W t → W ∗ as well under A. Let us define B ⊆ Ω as the event of t → 0, then we know that since P (A) > 0 and P (B) = 1, where P is the probability function for events in Ω, P (A ∩ B) = P (A). Therefore, conditional on the event of A, we have that t
a.s.−−→ 0 still holds. Now we consider any realization of A∩B. For the right-hand side of (22), as W̃ t is convergent and βtα−1t decreases to zero, by letting t approach infinity, we have that
lim t→∞
t + βtα−1t ( W̃ t −W 0 ) = 0 + 0 ( W ∗ −W 0 ) = 0.
By the Lipschitz continuity of ∇f , we have from (3) and (4) that 0 ≤ ∥∥∇f (W̃ t)−∇f (W t)∥∥ ≤ L∥∥W t − W̃ t∥∥.
As {W t} and {W̃ t} converge to the same point, we see that ∥∥W t − W̃ t∥∥→ 0, so ∇f (W̃ t)−
∇f(W t) also approaches zero. Hence, the limit of the right-hand side of (22) is lim t→∞ ∇f ( W̃ t ) − ( ∇f ( W t ) + t + βtα−1t ( W̃ t −W 0 )) = 0. (23)
On the other hand, for the left-hand side of (22), the outer semicontinuity of ∂ψ at W ∗ and the continuity of ∇f show that
lim t→∞
∇f(W̃ t) + ∂ψ(W̃ t) ⊆ ∂∇f(W ∗) + ψ (W ∗) = ∂F (W ∗). (24)
Substituting (23) and (24) back into (22) then proves that 0 ∈ ∂F (W ∗) and thus W ∗ ∈ Z.

heading: A.4 Proof of Theorem 2
text: Proof. Our discussion in this proof are all under the event that W̃ t →W ∗. From the argument in Appendix A.3, we can view α−1t V t as ∇f(W t) plus some noise that asymptotically decreases to zero with probability one as shown in (20). From Lemma 1, we know that W t →W ∗. From (21), there is U t ∈ ∂ψ ( W̃ t )
such that U t = −α−1t V t + α−1t βt ( W̃ t −W 0 ) . (25)
Moreover, we define γt := W t − W̃ t. (26) By combining (25)–(26) with (20), we obtain min
Y ∈∂F (W̃ t) ‖Y ‖
≤ ∥∥∇f (W̃ t)+ U t∥∥
= ∥∥∇f (W̃ t)−∇f (W t)− t − α−1t βt (W̃ t −W 0)∥∥
≤ ∥∥∇f (W̃ t)−∇f (W t)∥∥+ ‖ t‖+ α−1t βt∥∥W̃ t −W 0∥∥ ≤L‖γt‖+ ‖ t‖+ α−1t βt
(∥∥W ∗ − W̃ t∥∥+ ∥∥W 0 −W ∗∥∥) , (27) where we used the Lipschitz continuity of ∇f and the triangle inequality in the last inequality. We now separately bound the terms in (27). From that W t → W ∗ and W̃ t → W ∗, it is straightforward that ‖γt‖ → 0. The second term decreases to zero almost surely according to (20) and the argument in Appendix A.3. For the last term, since α−1t βt → 0, and∥∥W̃ t −W ∗∥∥→ 0, we know that
α−1t βt ∥∥W 0 −W ∗∥∥→ 0, α−1t βt∥∥W̃ t −W ∗∥∥→ 0.
Therefore, we conclude from the above argument and (27) that min
Y ∈∂F (W̃ t) ‖Y ‖ a.s.−−→ 0.
As f is smooth with probability one, we know that if ψ is partly smooth at W ∗ relative toM, then so is F = f + ψ with probability one. Moreover, Lipschitz-continuously differentiable functions are always prox-regular, and the sum of two prox-regular functions is still proxregular, so F is also prox-regular at W ∗ with probability one. Following the argument identical to that in Appendix A.3, we know that these probability one events are still probability one conditional on the event of W̃ t →W ∗ as this event has a nonzero probability. As W̃ t →W ∗ and ∇f(W̃ t) +U t a.s.−−→ 0 ∈ ∂F (W ∗) (the inclusion is from (5)), we have from the subdifferential continuity of ψ and the smoothness of f that F (W̃ t) a.s.−−→ F (W ∗). Since we also have W̃ t →W ∗ and minY ∈∂F (W̃ t) ‖Y ‖
a.s.−−→ 0, clearly( W̃ t, F
( W t ) , min Y ∈∂F(W̃ t) ‖Y ‖
) a.s.−−→ (W ∗, F (W ∗), 0) . (28)
Therefore, (28) and (5) together with the assumptions on ψ at W ∗ imply that with probability one, all conditions of Lemma 1 of Lee (2020) are satisfied, so from it, (6) holds almost surely, conditional on the event of W̃ t →W ∗.

heading: B Additional Discussions on Applications
text: We now discuss in more technical details the applications in Section 4.1, especially regarding how the regularizers satisfy the properties required by our theory.

heading: B.1 Structured Sparsity
text: We start our discussion with the simple `1 norm as the warm-up for the group-LASSO norm. It is clear that ‖W‖1 is a convex function that is finite everywhere, so it is prox-regular, subdifferentially continuous, and regular everywhere, hence we just need to discuss about the remaining parts in Definition 1. Consider a problem with dimension n > 0. Note that
‖x‖1 = n∑ i=1 |xi|,
and the absolute value is smooth everywhere except the point of origin. Therefore, it is clear that ‖x‖1 is locally smooth if xi 6= 0 for all i. For any point x∗, when there is an index set I such that x∗i = 0 for all i ∈ I and x∗i 6= 0 for i /∈ I, we see that the part of the norm corresponds to IC (the complement of I):∑
i∈IC |x∗i |
is locally smooth around x∗. Without loss of generality, we assume that I = {1, 2, . . . , k} for some k ≥ 0, then the subdifferential of ‖x‖1 at x∗ is the set
{sgn(x1)} × · · · × {sgn(xk)} × [−1, 1]n−k, (29) and clearly if we move from x∗ along any direction y := (y1, . . . , yk, 0, . . . , 0) with a small step, the function value changes smoothly as it is a linear function, satisfying the first condition of Definition 1. Along the same direction y with a small enough step, the set of subdifferential remains the same, so the continuity of subdifferential requirement holds. We can also observe from the above argument that the manifold should be Mx∗ = {x | xi = 0,∀i ∈ I}, and clearly it is a subspace of Rn with its normal space at x∗ being N := {y | 〈x∗, y〉 = 0} = {y | yi = 0,∀i ∈ IC}, which is clearly the affine span of (29) with the translation being (sgn(x1)× · · ·× sgn(xk), 0, . . . , 0). Moreover, indeed the manifolds are low dimensional ones, and for iterates approaching x∗, staying in this active manifold means that the (unstructured) sparsity of the iterates is the same as the limit point x∗. We also provide a graphical illustration of ‖x‖1 with n = 2 in Fig. 4. We can observe that for any x with x1 6= 0 and x2 6= 0, the function is smooth locally around any point, meaning that ‖x‖1 is partly smooth relative to the whole space at x (so actually smooth locally around x). For x with x1 = 0, the function value corresponds to the sharp valley in the graph, and we can see that the function is smooth along the valley, and this valley corresponds to the one-dimensional manifold {x | x1 = 0} for partial smoothness. Next, we use the same graph to illustrate the importance of manifold identification. Consider that the red point x∗ = (0, 1.5) is the limit point of the iterates of a certain algorithm, and the yellow points and black points are two sequences that both converge to x∗. If the iterates of the algorithm are the black points, then clearly except for the limit point itself, all iterates are nonsparse, and thus the final output of the algorithm is also nonsparse unless we can get to exactly the limit point within finite iterations (which is usually impossible for iterative methods). On the other hand, if the iterates are the yellow points, this is the case that the manifold is identified, because all points sit in the valley and enjoy the same sparsity pattern as the limit point x∗. This is why we concern about manifold identification when we solve regularized optimization problems. From this example, we can also see an explanation for why our algorithm with the property of manifold identification performs better than other methods without such a property. Consider a Euclidean space any point x∗ with an index set I such that x∗I = 0 and |I| > 0. This means that x∗ has at least one coordinate being zero, namely x∗ contains sparsity. Now let
0 := min i∈IC |x∗i |,
then from the definition of I, 0 > 0. Fro any sequence {xt} converging to x∗, for any ∈ (0, 0), we can find T ≥ 0 such that∥∥xt − x∗∥∥2 ≤ , ∀t ≥ T . Therefore, for any i /∈ I, we must have that xti 6= 0 for all t ≥ T . Otherwise, ‖xt − x∗‖2 ≥ 0, but 0 > ≥ ‖xt − x∗‖2, leading to a contradiction. On the other hand, for any i ∈ I, we can have xti 6= 0 for all t without violating the convergence. That being said, for any sequence converging to x∗, eventually the iterates cannot be sparser than x∗, so the sparsity level of x∗, or of its active manifold, is the local upper bound for the sparsity level of points converging to x∗. Therefore, if iterates of two algorithms converge to the same limit point, the one with a proven manifold identification ability clearly will produce a higher sparsity level. Similar to our example here, in applications other than sparsity, iterates converging to a limit point dwell on super-manifolds of the active manifold, and the active manifold is the minimum one that locally describes points with the same structure as the limit point, and thus identifying this manifold is equivalent to finding the locally most ideal structure of the application. Now back to the sparsity case. One possible concern is the case that the limit point is (0, 0) in the two-dimension example. In this case, the manifold is the 0-dimensional subspace {0}. If this is the case and manifold identification can be ensured, it means that limit point itself can be found within finite iterations. This case is known as the weak sharp minima (Burke & Ferris, 1993) in nonlinear optimization, and its associated finite termination property is also well-studied. For this example, We also see that ‖x‖1 is partly smooth at any point x∗, but the manifold differs with x∗. This is a specific benign example, and in other cases, partial smoothness might happen only locally at some points of interest instead of everywhere. Next, we further extend our argument above to the case of (7). This can be viewed as the `1 norm for each group and we can easily obtain similar results. Again, since the group-LASSO norm is also convex and finite everywhere, prox-regularity, regularity, and subdifferential continuity are not issues at all. For the other properties, we consider one group first, then the group-LASSO norm reduces to the `2 norm. Clearly, ‖x‖2 is smooth locally if x 6= 0, with the gradient being x/‖x‖2, but it is nonsmooth at the point x = 0, where the subdifferential is the unit ball. This is very similar to the absolute value, whose subdifferential at 0 is the interval [−1, 1]. Thus, we can directly apply similar arguments above, and conclude that for any W ∗, (7) is partly smooth at W ∗ with respect to the manifold MW∗ = {W | WIg = 0,∀g : W ∗Ig = 0}, which is again a lower-dimensional subspace. Therefore, the manifold of defining the partial smoothness for the group-LASSO norm exactly corresponds to its structured sparsity pattern.

heading: B.2 Binary Neural Networks
text: We continue to consider the binary neural network problem. For easier description, for the Euclidean space E we consider, we will use a vectorized representation for W,A ∈ E such that the elements are enumerated as W1, . . . ,Wn and α1, . . . , αn. The corresponding optimization problem can therefore be written as
min W,A∈E Eξ∼D [fξ (W )] + λ n∑ i=1 ( αi (wi + 1)2 + (1− αi) (wi − 1)2 + δ[0,1] (αi) ) , (30)
where given any set C, δC is the indicator function of C, defined as δC(x) = {
0 if x ∈ C, ∞ else.
We see that except for the indicator function part, the objective is smooth, so the real partly smooth term that we treat as the regularizer is
Φ(α) := n∑ i=1 δ[0,1](αi).
We note that for αi ∈ (0, 1), the value of δ[0,1](αi) remains a constant zero in a neighborhood of αi, and for αi /∈ [0, 1], the indicator function is also constantly infinite within a neighborhood. Thus, the point of nonsmoothness, happens only at αi ∈ {0, 1}, and similar to our discussion in the previous subsection, Φ is partly smooth along directions that we fix those αi at the boundary (namely, being either 0 or 1) unchanged. The identified manifold therefore corresponds to the entries of αi that are fixed at 0 or 1, and this can serve as the indicator for the desired binary pattern in this task.

heading: C Experiment Setting Details
text: For the weights wg of each group in (7), we consider the following setting. Let Li be a collection of index sets that belong to the i-th layer in the network, and denote
NLi := ∑ j∈Li |Ij |
the number of parameters in this layer, we set wg = √ NLi for all g ∈ Li for all i in our experiments. All ProxSSI parameter settings expect the regularization weight and the learning rate schedule follow the default values in their package. The detailed model and algorithm settings as well as statistics of the synthetic data of Section 5.1 are shown in Tables 3 and 4. Tables 5 to 7 provide the detailed setting of Section 5.2 while Tables 8 to 12 are for Section 5.3. For the modified VGG16 model, we follow Deleu & Bengio (2021) to eliminate all fully-connected layers except the output layer, and add one batch-norm layer (Ioffe & Szegedy, 2015) after each convolutional layer to simulate modern CNNs like those proposed in He et al. (2016); Huang et al. (2017). For ResNet50 in the structured sparsity experiment in Section 5.3, our version of ResNet50 is the one constructed by the publicly available script at https://github.com/weiaicunzai/ pytorch-cifar100. In the unstructured sparsity experiment presented in Section 5.4, for better comparison with existing works in the literature of pruning, we adopt the version of ResNet50 used by Sundar & Dwaraknath (2021).3 Table 13 provides detailed settings of Section 5.4. For RigL, we use the PyTorch implementation of Sundar & Dwaraknath (2021).
3https://github.com/varun19299/rigl-reproducibility.

heading: RMDA

heading: ProxSSI

heading: ProxSGD

heading: D More Details of the Experiments
text: In this section, we provide more details of the results of the experiments we conducted in the main text. In particular, in Fig. 5, we present the change of validation accuracies and group sparsity levels with epochs for the group sparsity tasks in Sections 5.2 and 5.3. We then present in Fig. 6 validation accuracies and unstructured sparsity level versus epochs for the task in Section 5.4. We note that althuogh it takes more epochs for RMDA to fully stabilize in terms of manifold identification, the sparsity level usually only changes in a very limited range once (sometimes even before) the validation accuracy becomes steady, meaning that we do not need run the algorithm for an unreasonably long time to obtain satisfactory results.

heading: E Other Regularizers for Possibly Better Group Sparsity and Generalization
text: A downside of (7) is that it pushes all groups toward zeros and thus introduces bias in the final model. For its remedy, minimax concave penalty (MCP, Zhang, 2010) is then proposed to penalize only the groups whose norm is smaller than a user-specified threshold. More precisely, given hyperparameters λ ≥ 0, c ≥ 1, the one-dimensional MCP is defined by
MCP(w) := { λ|w| − w 2
2c if|w| < cλ, cλ2
2 if|w| ≥ cλ.
Given the grouping G, MCP is then defined as (Breheny & Huang, 2009)
MCP ( WIg ) := { λ|Ig| ∥∥WIg∥∥2 − ‖WIg‖22c if ∥∥WIg∥∥ < cλ, cλ2|Ig|2
2 if ∥∥WIg∥∥ ≥ cλ, ψ(W ) := |G|∑ g=1 MCP ( WIg ) .
It is shown in Deleu & Bengio (2021) that group MCP regularization may simultaneously provide higher group sparsity and better validation accuracy than the group LASSO norm in vision and language tasks. Another possibility to enhance sparsity is to add another
`1-norm or entry-wise MCP regularization to the group-level regularizer. The major drawback of these approaches is the requirement of additional hyperparameters, and we prefer simpler approaches over those with more hyperparameters, as hyperparameter tuning in the latter can be troublesome for users with limited computational resources, and using a simpler setting can also help us to focus on the comparison of the algorithms themselves. The experiment in this subsection is therefore only for illustrating that these more complicated regularizers can be combined with RMDA if the user wishes, and such regularizers might lead to better results. Therefore, we train a simple CNN (LeNet5) on the MNIST dataset with such regularizers using RMDA and display the respective performance of various regularization schemes in Fig. 7. In this figure, group LASSO is abbreviated as GLASSO, `1-norm plus a group LASSO norm is abbreviated as L1GLASSO, group MCP is abbreviated as GMCP, and element-wise MCP plus group MCP is abbreviated as L1+GMCP. Our results exemplify that different regularization schemes might have different benefits on one of the criterion with proper hyperparameters tuning. The detailed numbers are reported in Table 14 and the experiment settings can be found in Table 15.

heading: RMDA

heading: ProxSSI

heading: ProxSGD

heading: RMDA

heading: ProxSSI

heading: ProxSGD

heading: RMDA

heading: ProxSSI

heading: ProxSGD

heading: RMDA

heading: ProxSSI

heading: ProxSGD

heading: RMDA

heading: ProxSSI

heading: ProxSGD

heading: RigL

heading: RMDA

Review summaries:
This paper proposes an algorithm for training neural networks with a regularization term for promoting desired structures in the model. The paper (claims to) prove that the proposed method can correctly identify within finite steps the underlying structure. The simulations then show that the method outperforms existing packages in identifying structured sparsity without compromising prediction accuracy.

In this paper, the authors propose Regularized Modernized Dual Averaging (RDMA) algorithm to train structured neural networks. This is a proximal method that achieves variance reduction without any extra cost per iteration. The algorithm does not require any extra hyperparameters than what stochastic gradient descent requires. Authors theoretically prove that structure identification is guaranteed after a finite number of iterations.  Experiments have been performed to obtain group sparse neural networks. Networks and datasets used in the experiments are simple logistic regression network on MNIST, 7-layer fully-connected network on FashionMNIST, LeNet5 on MNIST, and VGG16 on CIFAR10. The method has been compared with SGD, ProxSGD, and ProxSSI. 

Regularization is generally used to impose desired structure on NN during training. The authors developed RMDA (Regularized Modernized Dual Averaging) which uses the weighted average of the previous SG to compute the tentative update via proximal operation associated with the regularizer. Then the parameters of the model are updated in the direction of this tentative update with a pre-defined factor. The authors theoretically analyzed the performance of their algorithm and showed that under some assumptions, RMDA can identify the structure of model in finite iterations.