Year: 2019

Abstract:
With the increase of data volume and data dimension, sparse representation learning attracts more and more attention. For high-dimensional data, randomized block coordinate descent methods perform well because they do not need to calculate the gradient along the whole dimension. Existing hard thresholding algorithms evaluate gradients followed by a hard thresholding operation to update the model parameter, which leads to slow convergence. To address this issue, we propose a novel hard thresholding algorithm, called Semi-stochastic Block Coordinate Descent Hard Thresholding Pursuit (SBCD-HTP). Moreover, we present its sparse and asynchronous parallel variants. We theoretically analyze the convergence properties of our algorithms, which show that they have a significantly lower hard thresholding complexity than existing algorithms. Our empirical evaluations on real-world datasets and face recognition tasks demonstrate the superior performance of our algorithms for sparsity-constrained optimization problems.

heading: 1 INTRODUCTION
text: In modern high-dimensional data analytics, where the variable dimension can be equal or even larger than the number of samples, sparse representation learning has become a mainstream method to explore the potential real model of the problem and provides statistically reliable results. It has been applied to many diverse domains such as high-dimensional statistics (Bühlmann & Van De Geer, 2011), signal processing (Lai et al., 2013), and computer vision (Wright et al., 2008). Many sparse learning algorithms such as `1-norm convex relaxation methods (Negahban et al., 2009; Van de Geer et al., 2008) have been proposed in the past few decades. In order to get a smaller estimation error, `0-norm constrained algorithms are more prominent than `1-norm convex relaxation algorithms (Zhang et al., 2010). In this paper, we mainly focus on the following sparsity-constrained optimization problem:
min w
F (w) := 1
n n∑ i=1 fi(w), s.t., ‖w‖0 ≤ s, (1)
where F (w) is a finite-sum convex and smooth function, each function fi(w) is associated with the i-th sample, ‖w‖0 is the number of nonzero entries in variable w, and s represents the sparsity level. The goal of sparse representation learning is to recover w∗ based on the given data. Such a formulation encapsulates plentiful important problems, including sparse graphical model learning (Zhou et al., 2019), sparse linear/logistic regression (Blumensath & Davies, 2009; Foucart, 2011; Pati et al., 1993; Bahmani et al., 2013), and low-rank regression (Rohde et al., 2011).
However, Problem (1) is in general NP hard, which is caused by the non-convexity of the sparsity constraint. It makes us to obtain an approximate solution to Problem (1). One of the most widely used methods for solving this problem is the hard thresholding based gradient descent method. In recent years, sparse representation learning is becoming more common for large-scale and highdimensional data. However, deterministic gradient descent hard thresholding algorithms such as fast gradient descent hard thresholding (FG-HT) (Jain et al., 2014; Yuan et al., 2014) have to compute the gradients of all n component functions, which leads to a huge computing overhead. To address this issue, many stochastic hard thresholding methods have been proposed. For instance, Nguyen et al. (2017) proposed a stochastic gradient descent hard thresholding (SG-HT) algorithm, whose
gradient and hard thresholding complexities are both O(κs̃ log( 1 )). However, due to the stochastic sampling, SG-HT can only attain a sub-optimal estimation bound, as shown in Table 1, which is inferior to those of deterministic gradient methods such as FG-HT. Another limitation of SG-HT is that it requires that the restricted condition number κs̃ should not be larger than 4/3, which makes SG-HT hard for solving high-dimensional representation learning problems.
Recently, many stochastic variance reduced methods (e.g., SAG (Roux et al., 2012), SVRG (Johnson & Zhang, 2013)) and their accelerated variants such as (Defazio, 2016; Allen-Zhu, 2018) have been proposed to accelerate stochastic gradient methods for convex optimization. All these methods enjoy low per-iteration complexities comparable with stochastic gradient descent (SGD), and they also can attain improved convergence rates. By incorporating the variance reduction technique into sparsityconstrained optimization domain, Li et al. (2016b) proposed a stochastic variance reduced gradient hard thresholding (SVRG-HT) algorithm, which can converge more quickly and obtain a smaller estimation error than SG-HT. Moreover, SVRG-HT allows an arbitrarily large condition number in its theoretical analysis similar to FG-HT (Yuan et al., 2014). The gradient oracle and hard thresholding complexities of SVRG-HT are O((n+κs̃) log(1 )) and O(κs̃ log( 1 )), respectively. Nevertheless, it can not leverage the coordinate block to accelerate convergence. Chen & Gu (2016) proposed an accelerated stochastic block coordinate gradient descent hard thresholding (ASBCDHT) algorithm. The gradient oracle complexity of ASBCDHT is O((n+ κs̃|B|k ) log( 1 )), which is superior to SVRG-HT. Moreover, ASBCDHT also has the hard thresholding complexity, O(κs̃ log( 1 )). However, the hard thresholding complexity of ASBCDHT still scales linearly with κs̃, which is usually expensive for real-world sparse learning problems. For large-scale and high-dimensional data, Li et al. (2016a) also proposed the asynchronous parallel variant of SVRG-HT (called ASVRG-HT) by utilizing multicore architectures. Although it makes each processor to evaluate a stochastic gradient update on a global parameter stored in a shared memory in an asynchronous and lock-free mode, ASVRG-HT attains the similar gradient and hard thresholding complexities as its general version (i.e., SVRG-HT). Therefore, all the algorithms mentioned above have a linearly κs̃-dependent hard thresholding complexity. This motivates us to address the following key issue:
Can we design such algorithms that have a κs̃-independent hard thresholding complexity and a lower gradient oracle complexity?
To answer the above problem, we propose an efficient Semi-stochastic Block Coordinate Descent Hard Thresholding Pursuit (SBCD-HTP) algorithm and its sparse Asynchronous variant (ASBCDHTP). The oracle complexities and statistical estimation error of the proposed algorithms and other hard thresholding methods are summarized in Table 1. We highlight several theoretical advantages of the proposed algorithms over the state-of-the-art methods as follows:
• SBCD-HTP and ASBCD-HTP substantially improve the restricted condition on the sparsity level s, i.e., s=Ω(κŝs∗), while most related algorithms such as SVRG-HT require s=Ω(κ2s̃s ∗).
• The statistical estimation error of both our algorithms is better than those of FG-HT, SVRG-HT and ASBCDHT, as our algorithms have a smaller sparsity level s=Ω(κŝs∗) and thus have a smaller cardinality, i.e., ŝ= 2s+s∗. This makes the restricted condition number κŝ of SBCD-HTP smaller than other algorithms. Moreover, the statistical error of our algorithms is also smaller than that of SG-HT (i.e., 1n ∑n i=1‖∇fi(w∗)‖) due to the large magnitude of ‖∇fi(w∗)‖ (Zhou et al., 2018b).
• Both SBCD-HTP and ASBCD-HTP enjoy a κs̃-independent hard thresholding complexity, which is significantly better than the state-of-the-art hard thresholding algorithms. That is, the hard thresholding complexity of our algorithms is κs̃ times lower than that of FG-HT, SG-HT, SVRG-HT and ASBCDHT with s̃ = 2Ω(κ2s̃s
∗) + s∗. Since both our algorithms have a significantly lower hard thresholding complexity, they are more suitable for handling large-scale sparse representation learning problems, especially for high-dimensional data.
• For sparsity-constrained problems, the gradient oracle complexities of both SBCD-HTP and ASBCD-HTP (or S2BCD-HTP) are much lower than those of the state-of-the-art hard thresholding algorithms, and that of S2BCD-HTP is similar to ASBCD-HTP under conditions on the delay. Actually, the gradient oracle complexity of S2BCD-HTP is better than SBCD-HTP when we deal with sparse data sets, since S2BCD-HTP only evaluates the common nonzero point of the random sample and the sampled block. Furthermore, their gradient oracle complexities are significantly lower than those of ASBCDHT and SVRG-HT, and much lower than that of FG-HT, since they require a smaller sparsity level s= Ω(κŝs∗) and also have a smaller value of κŝ with ŝ= 2s+s∗ compared with the restricted condition number κs̃ used in other hard thresholding algorithms (Zhou et al., 2018b).

heading: 2 RELATED WORK
text: In this section, we briefly discuss the relevant research to our work, which beyonds the sparsityconstrained optimization domain. Traditional gradient descent is computational expensive at each iteration, and stochastic gradient descent has a low per-iteration complexity while obtaining a large variance in estimating. Thus, many stochastic variance reduced algorithms (Defazio et al., 2014; Roux et al., 2012; Shalev-Shwartz & Zhang, 2013; Shang et al., 2018; Johnson & Zhang, 2013) and their variants (Schmidt et al., 2017; Konečnỳ et al., 2015; Xiao & Zhang, 2014; Hannah et al., 2018) have been proposed. The stochastic variance reduced methods are very promising for many machine learning problems including sparse learning problems.
In contrast to gradient descent methods, coordinate descent algorithms have received increasing attention due to their successful applications in high dimensional problems (Breheny & Huang, 2011; Friedman et al., 2007). Among them, randomized block coordinate descent (RBCD) updates a block of coordinate with respect to the entire training instances. The per-iteration cost is significantly lower than gradient descent, while it still has a relationship with n component functions. Some stochastic block coordinate descent (SBCD) algorithms such as (Dang & Lan, 2015; Xu & Yin, 2015; Konečnỳ et al., 2017) were proposed. Such algorithms compute the stochastic gradient with respect to a random sample restricted to one randomized coordinate block. Therefore, these algorithms randomly sample a block of features and data instances at each iteration. However, they can only achieve a sublinear convergence rate (Zhang & Gu, 2016). Recently, some mini-batch randomized block coordinate descent algorithms such as (Zhao et al., 2014; Wang & Banerjee, 2014) were proposed to accelerate the convergence of stochastic block coordinate gradient descent. Our work studies over the above research by considering a sparsity-constrained optimization problem. The proposed algorithms enjoy a lower computational complexity in both gradient evaluation and hard thresholding computation while obtaining a linear convergence performance. Moreover, its sparse and asynchronous variants are also proposed to deal with sparse high-dimensional data.

heading: 3 PRELIMINARIES
text: Throughout this paper, we use w∗ to denote the optimal solution of Problem (1), and its optimal sparsity level is s∗ with ‖w∗‖0≤ s∗. ‖w‖ is the Euclidean norm for a vector w∈Rd and ‖w‖∞ is the largest absolute entry in w. The hard thresholding operation HT (w, s) preserves the s largest entries of w in magnitude for vector w and the rest entries are set to zero. We use supp(w) to denote the support set of w, i.e., the indices of its nonzero entries. Given an index set I, we define IC as the complement set of I, and vI ∈ Rd, where [vI ]j=vj if j∈I and [vI ]j=0 if j /∈I. Given an integer
Algorithm 1 Semi-stochastic Block Coordinate Descent Hard Thresholding Pursuit (SBCD-HTP) Input: Number of outer-loops R, number of inner-loops m, step size η, sparsity level s. Initialize: w̃0.
1: for r = 0, 1, . . . , R− 1 do 2: w0 = w̃ = w̃r, ∇F (w̃) = 1n ∑n i=1∇fi(w̃), G̃ = supp(w̃); 3: for t = 0, 1, . . . ,m− 1 do 4: Randomly sample a mini-batch B from [n] uniformly; 5: Randomly sample a block jt from [k] uniformly, and S = G̃ ∪ Gjt ; 6: ∇Sg(wt) = 1|B| ∑ it∈B∇Sfit(w
t)−∇Sfit(w̃) +∇SF (w̃); 7: wt+1 = wt − η∇Sg(wt); 8: end for 9: w̃r+1 = HT (wm, s);
10: end for Output: w̃R.
n ≥ 1, we define [n] = {1, · · · , n}. For a set B, we denote its cardinality by |B|. Moreover, we use the common notations of Ω(·) and O(·) to characterize the asymptotics of two real sequences. Throughout the analysis, we make two important assumptions on the objective function, which are commonly used in the analysis of hard thresholding algorithms (Li et al., 2016b; Shen & Li, 2017; Chen & Gu, 2016; Gao & Huang, 2018). In high-dimensional sparse learning, the per-iteration hard thresholding operation can be time-consuming or even more expensive than gradient computation. Thus, we also take the complexity of hard thresholding into our consideration.
Assumption 1 (Restricted Strong Convexity) A differentiable function F (·) is restricted ρ−ŝ - strongly convex at sparsity level ŝ, if there exists a uniform constant ρ−ŝ > 0 such that for any w,w′ ∈ Rd with ‖w − w′‖0 ≤ ŝ, we have
F (w)− F (w′)− 〈∇F (w′), w − w′〉 ≥ ρ−ŝ 2 ‖w − w′‖2. (2)
Assumption 2 (Restricted Strong Smoothness) For any i ∈ [n], the differentiable function fi(·) is restricted ρ+ŝ -strongly smooth at sparsity level ŝ, if there exists a uniform constant ρ + ŝ > 0 such that for any w,w′ ∈ Rd with ‖w − w′‖0 ≤ ŝ, we have
fi(w)− fi(w′)− 〈∇fi(w′), w − w′〉 ≤ ρ+ŝ 2 ‖w − w′‖2. (3)
Moreover, the restricted condition number is defined as: κŝ = ρ+ŝ /ρ − ŝ . Since the complexity of one hard thresholding operation is linear with d, which is similar to a computational cost of a stochastic gradient, we define the hard thresholding complexity as a metric for hard thresholding algorithms.
Definition 1 (Hard Thresholding Complexity) In a hard thresholding operation, a vector w is fed into HT (·, s) and then the output HT (w, s) is obtained.
Both gradient oracle and hard thresholding complexities can more comprehensively reflect the overall computational cost of first-order hard thresholding algorithms. Therefore, their per-iteration cost is dominated by both the gradient evaluation and hard thresholding operations.
We also define the optimization error and statistical error to clearly understand our theoretical results.
Definition 2 (Optimization Error and Statistical Error) The optimization error is the difference which will decrease with the increasing of iterations. A statistical error is the (unknown) difference between the retained value and the true value.
This means that the optimization error will decrease to zero when the number of iteration is quite large. However, a statistical error will still work during the optimized process. Therefore, the larger the statistical error, the more difference to the true value.

heading: 4 SEMI-STOCHASTIC BLOCK COORDINATE DESCENT HARD THRESHOLDING PURSUIT
text: In this section, we propose a novel Semi-stochastic Block Coordinate Descent Hard Thresholding Pursuit (SBCD-HTP) algorithm, and also theoretically analyze its convergence properties.

heading: 4.1 OUR SBCD-HTP ALGORITHM
text: The proposed SBCD-HTP algorithm is summarized in Algorithm 1. In each outer-loop, we select a snapshot point w̃ and compute its full gradient ∇F (w̃). In each inner-loop, we uniformly and randomly select a mini-batch samples B and a block of coordinates Gjt , where jt ∈ {1, 2, . . . , k}, {G1 . . .Gk} is a partition of all the d coordinates and is divided uniformly at random. We usually set k=10 in our experiments. Note that, different from ASBCDHT, in line 5 we set S to be a union set of the support set of snapshot point w̃ and Gjt . In this way, the coordinate of the support set S which includes the possibly optimal coordinates can be updated, and thus this can reduce the error caused by randomness. Moreover, SBCD-HTP enjoys the dual advantages of randomized block coordinate descent and semi-stochastic gradient descent to optimize the objective as shown in line 6. Note that we perform the hard thresholding operation in only the outer-loop, which is different from all existing hard thresholding algorithms. Performing hard thresholding operations with high frequency or prematurely in inner-loop will make the gradient information lost, which is the main reason for slow convergence and low accuracy of existing hard thresholding sparsity-constrained algorithms. Our hard thresholding algorithm not only guarantees the sparsity of model parameters, but also uses the full dimensional information without the interference of hard thresholding operations to update model parameter wt+1. Our theoretical analysis and experimental results show that SBCD-HTP has a faster convergence rate and more accurate results.
We have to clarify the difference between the proposed SBCD-HTP and ASBCDHT in Chen & Gu (2016) for clear improvements. There are mainly two changes: a) we design a bigger support set with the union of the sampled block Gjt with the one of the current support of the snapshot point w̃, since the current support of snapshot point contain the most related optimized information, while that of Chen & Gu (2016) only use the sampled block. b) we move the hard thresholding operation outside of the inner-loop, since the truncated function (hard thresholding) may loss many information and bringing higher complexity. Moreover, the asynchronous variant of SBCD-HTP is also proposed in this paper to optimize the sparse problem 1 in the sparse data. This is the first to design this type of algorithm in sparse learning field, and the theoretical proof is also provided.

heading: 4.2 CONVERGENCE ANALYSIS
text: We first analyze the convergence behavior of SBCD-HTP. The main result is summarized in the following theorem, whose proof is provided in Section C.
Theorem 1 Suppose F (w) is ρ−ŝ -strongly convex and each function fi(w) is ρ + ŝ -strongly smooth with parameter ŝ= 2s+s∗. Let κŝ = ρ+ ŝ
ρ− ŝ and the sparsity level s ≥ Ω(κŝs∗). In addition, assume that the learning rate η≤ 1
48ρ+ ŝ
, the number of inner-loops m≥ 1800κŝ, and the number of blocks
k=10, then the convergence rate γ= k 2ηmρ−
ŝ (1−ω) +
8ρ+ ŝ η(n−|B|)
|B|(1−ω)(n−1) ≤ 1 2 . For the sparsity-constrained
problem (1), the output of Algorithm 1 satisfies
E[F (w̃r)− F (w∗)] ≤ (1 2 )rE[F (w̃0)− F (w∗)] + η 1− ω ‖∇ĨF (w ∗)‖2, (4)
where Ĩ = supp(HT (∇F (w∗), 2s)) ∪ supp(w∗) and ω = 8ρ+ŝ η(1 + n−|B| |B|(n−1) ).
Theorem 1 shows that SBCD-HTP enjoys a linear convergence rate for solving sparsity-constrained problem (1). In order to ensure that the linear converging term in Eq. (4) satisfies ( 12 )
rE[F (w̃0) − F (w∗)] ≤ , we can obtain the following corollary.
Corollary 1 Suppose the conditions in Theorem 1 hold. To achieve ( 12 ) rE[F (w̃0)−F (w∗)] ≤ , the gradient oracle complexity of SBCD-HTP is O((n+ κŝ|B|k ) log( 1 )), and its hard thresholding complexity is O(log( 1 )).
Remark 1 This corollary is inferred under the high-dimensional condition. Actually, the worst condition of the gradient oracle complexity of SBCD-HTP is O((n + κŝ|B| d k+s d ) log( 1 )). Due to the small sparsity level s required in our Theorem 1 and the overlaps between Gjt and G̃, we can obtain the gradient oracle complexity O((n+ κŝ|B|k ) log( 1 )) for SBCD-HTP.
Table 1 summarizes the properties of some hard thresholding methods. Compared with the gradient oracle complexities of FG-HT (Yuan et al., 2014) and SVRG-HT (Li et al., 2016b) (i.e., O(nκs̃ log( 1 )) andO((n+κs̃) log( 1 )), respectively), SBCD-HTP has a much lower oracle complexity. For ASBCDHT (Chen & Gu, 2016), whose gradient oracle complexity isO((n+κs̃|B|k ) log( 1 )), the gradient oracle complexity of SBCD-HTP is O((n+ κŝ|B|k ) log( 1 )) and is much lower than ASBCDHT, as κŝ of SBCD-HTP is usually comparable to or even smaller than those of others (note that s in ŝ=2s+s∗ is required to be smaller than others as shown below). Moreover, SBCD-HTP allows s= Ω(κŝs ∗), which is considerably superior to the condition of s= Ω(κ2s̃s ∗) required in other hard thresholding algorithms such as (Zhou et al., 2018b). In particular, with the same parameter settings of |B| and k, the gradient oracle complexity of our algorithm outperforms ASBCDHT. SBCD-HTP also attains a lower hard thresholding complexity than the state-of-the-art hard thresholding methods. That is, the hard thresholding complexity of SBCD-HTP is O(log( 1 )), which is κs̃-independent and is significantly lower than those of other methods. This is a significant improvement on hard thresholding complexity. Theorem 1 immediately implies the following results.
Corollary 2 Suppose the conditions in Theorem 1 hold. The output of Algorithm 1 satisfies
E‖w̃r − w∗‖ ≤
√ 2( 12 )
r[F (w̃0)− F (w∗)] ρ−ŝ +
( 2
ρ−ŝ +
2η
1− ω
)√ ŝ‖∇F(w∗)‖∞. (5)
The right-hand side of Eq. (5) consists of two terms. The first term is the optimization error, which approaches zero with the increase of r. The second term corresponds to the statistical error, which is proportional to √ ŝ‖∇F (w∗)‖∞. One can also observe that the statistical error bound is better
than those of FG-HT and ASBCDHT, since ŝ in the statistical error O( √ ŝ‖∇F (w∗)‖∞) has a smaller cardinality (i.e., 2Ω(κŝs∗)+s∗), while existing algorithms require a larger sparsity level s′ = Ω(κ2s̃s
∗) and have a larger cardinality. Therefore, compared with other existing algorithms, SBCD-HTP enjoys a smaller statistical error, especially for the problems with a large restricted condition number. It is usually better than the error bound O( √ s̃‖∇F (w∗)‖∞+‖∇ĨF (w∗)‖) with s̃=2Ω(κ2s̃s ∗)+s∗ in SVRG-HT. Moreover, the error bound of SBCD-HTP is much smaller than that
of SG-HT (Nguyen et al., 2017) (i.e., O( 1n ∑n i=1‖∇fi(w∗)‖)), since the magnitude of the individual gradient norm ‖∇fi(w∗)‖ can still be relatively large (Zhou et al., 2018b).

heading: 5 ASYNCHRONOUS VARIANT OF SBCD-HTP
text: In this section, we propose the serial sparse and asynchronous parallel variants of SBCD-HTP for high-dimensional sparse data sets. Our Asynchronous Sparse Semi-stochastic Block Coordinate Descent Hard Thresholding Pursuit (ASBCD-HTP) algorithm is summarized in Algorithm 2. The main difference between ASBCD-HTP and SBCD-HTP is sparse approximate gradients as in (Mania et al., 2017; Zhou et al., 2018a). In order to perform fully sparse updates, we use a diagonal matrix D to re-weigh the dense fully gradient∇F (w̃). The entries of D are the inverse probabilities of the corresponding coordinates belonging to a uniformly sampled support Tit of sample it. Let Pit be the projection matrix for the support Tit , we can get Dit since Dit = PitD. Then we can find that Eit[Dit∇F (w̃)]=∇F (w̃). In this setting, we can compute the full gradient ∇F (w̃) also in a parallel way. In order to reduce thread interference and fully utilize the sparsity of datasets, we restrict the sparse gradient ∇g([ŵ]Tit) to the intersection set of S and Tit and obtain ∇Sg([ŵ]Tit). Specifically, each thread of ASBCD-HTP computes sparse approximate gradient independently and then atomically write this gradient to the shared variable w̃r. Note that Sparse variant of SBCDHTP (S2BCD-HTP) can be viewed as the single-thread version of ASBCD-HTP, which has linear convergence and the same oracle complexity as ASBCD-HTP. The detailed description and analysis of S2BCD-HTP is listed in Section D. Next we give the main theoretical result of ASBCD-HTP.
Algorithm 2 Asynchronous Sparse SBCD-HTP (ASBCD-HTP) Input: The number of outer-loops R, number of inner-loops m, step size η, sparsity level s. Initialize: w̃0.
1: for r = 0, 1, . . . , R− 1 do 2: w0 = w̃ = w̃r, G̃ = supp(w̃), and ∇F (w̃) = 1n ∑n i=1∇fi(w̃); //computed in parallel 3: t = 0; //inner loop counter 4: compute the while loop in parallel 5: while t < m do 6: Randomly sample it from [n] uniformly, and randomly sample jt from [k] uniformly; 7: Tit :=support of sample it, and t = t+ 1; //atomic increase counter t 8: S = G̃ ∪ Gjt , [ŵ]Tit := inconsistent read of shared variable [w̃
r]Tit ; 9: ∇Sg([ŵ]Tit ) = ∇Sfit([ŵ]Tit )−∇Sfit([w̃]Tit ) +Dit∇SF (w̃);
10: [w̃r]S⋂Tit = [w̃r]S⋂Tit − η∇Sg([ŵ]Tit ); //atomic write 11: end while 12: w̃r+1 = HT (w̃r, s); 13: end for Output: w̃R.
Theorem 2 Suppose F (w) is ρ−ŝ -strongly convex and each function fi(w) is ρ + ŝ -strongly smooth with parameter ŝ = 2s+s∗. Let κŝ = ρ+ ŝ
ρ− ŝ and the sparsity level s ≥ Ω(κŝs∗). We assume that the step size η = 1
60ρ+ ŝ
, the number of inner-loops m ≥ 120κŝ, the number of blocks k = 10 and
τ ≤min{ 3 5 √ ∆ , 2κŝ,
√ 2κŝ√
∆ } (i.e., the linear speedup condition). Then for the sparsity-constrained
problem (1), the gradient oracle and hard thresholding complexities of Algorithm 2 are
O ((n+ κŝ/k) log(1/ )) and O (log(1/ )) , (6)
where τ denotes the maximum number of concurrent threads (Mania et al., 2017) and ∆ = maxj=1···d pj , which is a indicator to measure the sparsity of datasets (Leblond et al., 2017).
Theorem 2 implies that ASBCD-HTP (or S2BCD-HTP) can still obtain similar gradient oracle complexity and hard thresholding complexity as SBCD-HTP while obtaining a linear speedup ratio. Compared with ASVRG-HT (Li et al., 2016a), ASBCD-HTP has a much lower gradient complexity, i.e.,O((n+κŝk ) log 1 ) for ASBCD-HTP vs. O((n+κs̃) log 1 ) for ASVRG-HT. For hard thresholding complexity, ASBCD-HTP is κs̃ times lower than ASVRG-HT, which is a significant improvement in the asynchronous setting. Since ASVRG-HT has to compute hard thresholding operations at each inner-loop, the practical performance may be worse, especially for high-dimensional data. We also show this improvement by our empirical evaluations.

heading: 6 EXPERIMENTS
text: In this section, we evaluate the performance of SBCD-HTP, S2BCD-HTP and ASBCD-HTP for solving sparse linear regression and sparse logistic regression on real-world datasets. All the dense algorithms were implemented in MATLAB, while the sparse and asynchronous algorithms were implemented in C++ and executed through MATLAB interface for a fair comparison. We also test the performance for face recognition tasks. The detailed information is described in Section F. Since there is no ground truth on real-world datasets, we run all these baselines sufficiently long until ‖wr−wr+1‖/‖wr‖ ≤ 10−6, and then we use the minimum F (wr) as the approximate optimal value F ∗ for sub-optimality estimation in our experimental results Zhou et al. (2018b).

heading: 6.1 DENSE ALGORITHMS
text: We ran the experiments of all the dense algorithms on a PC with an Intel i7-7700 CPU and 32GB RAM. We compare SBCD-HTP with several state-of-the-art sparsity-constrained algorithms, including FG-HT (Yuan et al., 2014), SG-HT (Nguyen et al., 2017), SVRG-HT (Li et al., 2016b), ASBCDHT (Chen & Gu, 2016) and FNHTP (Chen & Gu, 2017). The detailed descriptions of the
algorithms and datasets are provided in Section F. We set the sparsity level s= 200 for real-word datasets. For the number of inner-loops, we set m = 2n for SBCD-HTP and ASBCDHT, m = n for SVRG-HT as suggested in (Li et al., 2016b; Chen & Gu, 2016) and 1000 (sparse linear regression)/3000 (sparse logistic regression) for FNHTP. We also set the batch-size of FNHTP to s, and that of ASBCDHT and SBCD-HTP to 5. For the block coordinate descent methods, we set the number of coordinate blocks k to 10. All the algorithms are tuned to their best performance.
We report the performance of all the algorithms in terms of both effective passes and CPU time. More experimental results are presented in Section F. Figures 1 and 2 show that for both sparse linear regression and sparse logistic regression problems, our algorithm converges significantly faster than the baseline algorithms in terms of both effective passes and CPU time. The main reason is that we use all the information of gradients, while the related algorithms use less information caused by too frequent and premature hard thresholding operation in each inner-loop.

heading: 6.2 ASYNCHRONOUS ALGORITHMS
text: We ran the experiments of all the asynchronous algorithms on a PC with an Intel Xeon(R) Gold 5120 CPU and 64GB RAM. We compare our sparse and asynchronous algorithms with several stateof-the-art algorithms, including ASG-HT, ASVRG-HT and A2SBCD-HT (asynchronous variant of ABSCDHT). For the algorithms, we set the number of threads to 20.
We only report the performance of all the algorithms in terms of CPU time since the performance in terms of effective passes is similar to the dense case. More experimental results are reported in Section F. Figure 3 shows that ASBCD-HTP and S2BCD-HTP significantly outperform than other algorithms. Moreover, ASBCD-HTP has a faster convergence performance than S2BCD-HTP because of asynchronous parallel acceleration. There are three main reasons for the advantages of our ASBCD-HTP: a) We get the set S from the support set of snapshot point w̃ and Gjt . In this way, the coordinate of the support set which appears as the possible optimal solution can be updated, and thus reducing the error caused by randomness. b) Because we put the hard thresholding operation in each outer-loop, the computation cost per thread of our algorithm at each inner-loop is less, while each thread of the baseline algorithms has a hard thresholding operation, which is expensive especially for high-dimensional datasets. c) Our algorithm updates w̃r with all information of gradients, while the baseline algorithms update w̃r with less information caused by too frequent and premature hard thresholding operation in each inner-loop. We also evaluate the improvement of asynchronous parallel by running the same passes with different numbers of threads. We calculate the speed-up
ratio based on the running time of a single thread. From Figure 4, we can see that our ASBCD-HTP is accelerated by a nearly linear ratio.

heading: 7 CONCLUSION
text: In this paper, we proposed an efficient semi-stochastic block coordinate descent hard thresholding pursuit (SBCD-HTP) method for sparse representation learning. We proved that SBCD-HTP attains the gradient oracle complexity of O((n+ κŝ|B|k ) log( 1 )), and the hard thresholding complexity of O(log( 1 )), respectively, which is an improvement over the existing algorithms. Moreover, we also presented the sparse and asynchronous parallel variants of SBCD-HTP. As far as we know, our algorithms are the first to use the full information of gradient to update. This is an improvement on sparse representation learning because it can be easily extended to other sparse learning algorithms.

heading: APPENDIX

heading: A KEY LEMMAS
text: Lemma 1 ((Jain et al., 2014)) For any index set I , any w ∈ RI and HT (·, s) : Rd → Rd be the hard thresholding operator, which keeps the largest s entries (in magnitude) and sets the other entries equal to zero. Then for any w∗ ∈ RI such that ‖w∗‖0 ≤ s∗, we have
‖HT (w, s)− w‖2 ≤ |I| − s |I| − s∗ ‖w − w∗‖2.
Lemma 2 ((Li et al., 2016b)) Let w∗ ∈ Rd be the optimal sparse vector such that ‖w∗‖0 ≤ s∗ and HT (·, s) : Rd → Rd be the hard thresholding operator, which keeps the largest s entries (in magnitude) and sets the other entries equal to zero. Given s > s∗ for any vector w ∈ Rd, we have
‖HT (w, s)− w∗‖2 ≤ (1 + 2 √ s∗√
s− s∗ )‖w − w∗‖2.
Lemma 3 ((Nesterov, 2004)) For any given vector w,w′ ∈ Rd, suppose F (·) also satisfies RSC and RSS conditions, then the following inequality holds
‖∇F (w)−∇F (w′)‖2 ≤ 2ρ+ŝ [F (w)− F (w ′) + 〈∇F (w′), w′ − w〉].
Lemma 4 Suppose F (w) is ρ−ŝ -strongly convex and each function fi(w) is ρ + ŝ -strongly smooth with parameter ŝ = 2s + s∗. I = supp(w∗) ∪ supp(w̃r) ∪ supp(w̃r+1), |B| is the batch size. For any wt ∈ Rd and the sample it ∈ B, denote νit = (∇Ifit(wt)−∇Ifit(w∗))−(∇IF (wt)−∇IF (w∗)), we have
E ∥∥∥∥∥ 1|B|∑ it∈B νit ∥∥∥∥∥ 2 = n− |B| |B|(n− 1) E‖∇Ifi(wt)−∇Ifi(w∗) +∇IF (w∗)−∇IF (wt)‖2. (7)

heading: Proof. See Section B.1 for the proof of Lemma 4.
text: Lemma 5 Suppose F (w) is ρ−ŝ -strongly convex and each function fi(w) is ρ + ŝ -strongly smooth with parameter ŝ = 2s + s∗. Let w∗ ∈ Rd be the optimal sparse vector with ‖w∗‖0 ≤ s∗. I = supp(w∗) ∪ supp(w̃r) ∪ supp(w̃r+1), δ > 1 is a uniform constant factor, and |B| is the batch size. For any wt, w̃ ∈ Rd and the sample it, denote∇Sg(wt) = 1|B| ∑ it∈B∇Sfit(w
t)−∇Sfit(w̃)+ ∇SF (w̃), then we can bound E‖∇SgI(wt)‖2 as follows:
E‖∇SgI(wt)‖2 ≤ 1
k (16δρ+ŝ (1 + n− |B| |B|(n− 1) )E[F (wt)− F (w∗)] + 4‖∇IF (w∗)‖2
+ 16δρ+ŝ n− |B| |B|(n− 1) E[F (w̃)− F (w∗)]),
Proof. See Section B.2 for the detailed proof of Lemma 5.
Lemma 6 Suppose F (w) is ρ−ŝ -strongly convex and each function fi(w) is ρ + ŝ -strongly smooth with parameter ŝ = 2s + s∗. Let w∗ ∈ Rd be the optimal sparse vector with ‖w∗‖0 ≤ s∗. I = supp(w∗)∪supp(w̃r)∪supp(w̃r+1), δ > 1 is a uniform constant factor, andDm = maxj=1,··· ,d 1pj . For any wt, w̃ ∈ Rd and the sample it, denote ∇Sg(wt) = ∇Sfit(wt)−∇Sfit(w̃) +Dit∇SF (w̃), we can bound E‖∇SgI(wt)‖2 as follows:
E‖∇SgI(wt)‖2 ≤ 1
k (12δρ+ŝ [F (w t)− F (w∗)] + (3 + 12(D2m − 1))‖∇IF (w∗)‖2
+ (24δρ+ŝ + 48δρ + ŝ (D 2 m − 1))[F (w̃)− F (w∗)]),
where pj is the probability that dimension j belonging to the support set of randomly sampling sample it.
Proof. See Section B.3 for the proof of Lemma 6.
This lemma is different from other variance reduced results for sparsity-constrained problem (Mania et al., 2017; Johnson & Zhang, 2013; Chen & Gu, 2016; Shang et al., 2018). Since this variance reduced approximate gradients introduce an approximate full gradient to make fully sparse update in cardinality constraint problem.

heading: B PROOFS OF KEY LEMMAS
text: In this section, we prove three key lemmas.

heading: B.1 PROOF OF LEMMA 4
text: Proof. Let νit = (∇Ifit(wt)−∇Ifit(w∗))− (∇IF (wt)−∇IF (w∗)). We have
E ∥∥∥∥∥ 1|B|∑ it∈B νit ∥∥∥∥∥ 2 = 1 |B|2 E ∑ it,i′t∈B ν>it νi′t
= 1 |B|2 E ∑
it 6=i′t∈B
ν>it νi′t + 1
|B| E‖νi‖2
= |B| − 1 |B|n(n− 1) ∑ i6=i′ ν>i νi′ + 1 |B| E‖νi‖2
= |B| − 1 |B|n(n− 1) ∑ i,i′ ν>i νi′ − |B| − 1 |B|(n− 1) E‖νi‖2 + 1 |B| E‖νi‖2
= n− |B| |B|(n− 1) E‖νi‖2.
(11)
This proof is similar to that of (Zheng & Kwok, 2016).

heading: B.2 PROOF OF LEMMA 5
text: Proof. It is obvious that the stochastic variance reduced gradient satisfies
E[∇g(wt)] = E
[ 1
|B| ∑ it∈B (∇fit(wt)−∇fit(w̃)) + µ
] = ∇F (wt). (12)
Thus,∇g(wt) is an unbiased estimator of the full gradient∇F (wt). Now, we bound E‖∇gI(wt)‖2. For any i ∈ {1, 2, · · · , n}, we define the following function
hi(w) = fi(w)− fi(w∗)− 〈∇fi(w∗), w − w∗〉. (13) It is easy to find that ∇hi(w∗) = 0, which implies that hi(w∗) = minw hi(w). For any vector w, we have
0 = hi(w ∗) ≤ min η hi(w − η∇Ihi(w))
≤ min η hi(w)− η〈∇hi(w),∇Ihi(w)〉+ η2ρ+ŝ 2 ‖∇Ihi(w)‖2
= min η hi(w)− η‖∇Ihi(w)‖2 + η2ρ+ŝ 2 ‖∇Ihi(w)‖2
= hi(w)− 1
2ρ+ŝ ‖∇Ihi(w)‖2,
(14)
where the second inequality follows from the RSS condition, and the second equality holds due to the fact that 〈∇hi(w),∇Ihi(w)〉 = ‖∇Ihi(w)‖2, and the last equality holds due to the fact that η = 1
ρ+ ŝ
minimizes the function. Then, we have
‖∇Ifi(w)−∇Ifi(w∗)‖2 ≤ 2ρ+ŝ [fi(w)− fi(w ∗)− 〈∇Ifi(w∗), w − w∗〉]. (15)
Since the sampling i is chosen uniformly from {1, 2, · · · , n}, we have
E‖∇Ifi(w)−∇Ifi(w∗)‖2
= E‖∇IF (w)−∇IF (w∗)‖2
= 1
n n∑ i=1 ‖∇Ifi(w)−∇Ifi(w∗)‖2
≤ 2ρ+ŝ [F (w)−F(w ∗)− 〈∇IF (w∗), w − w∗〉] ≤ 2ρ+ŝ [F (w)− F (w ∗) + |〈∇IF (w∗), w − w∗〉|] ≤ 4δρ+ŝ [F (w)− F (w ∗)],
(16)
where the last inequality follows from the restricted strong convexity of F (w). And we use a simple underlying consensus that there must exist a uniform factor δ > 1 which makes the following inequality true:
δ(F (w)− F (w∗)) ≥ |〈∇IF (w∗), w − w∗〉|+ ρ−ŝ 2 ‖w − w∗‖2.
This property is called the extended restricted strong convexity used in the following proof. Therefore, we have
E‖∇gI(wt)‖2 = E ∥∥∥∥∥ 1|B|∑ it∈B ( ∇Ifit(wt)−∇Ifit(w̃) ) +∇IF (w̃) ∥∥∥∥∥ 2
= E ∥∥∥∥∥ ( 1 |B| ∑ it∈B ∇Ifit(wt)−∇IF (wt)−∇Ifit(w∗) +∇IF (w∗) ) +∇IF (wt)−∇IF (w∗) +∇IF (w∗)
−
( 1
|B| ∑ it∈B ∇Ifit(w̃)−∇IF (w̃)−∇Ifit(w∗) +∇IF (w∗) )∥∥∥∥∥ 2
≤ 4E ∥∥∥∥∥ 1|B|∑ it∈B ∇Ifit(wt)−∇IF (wt)−∇Ifit(w∗) +∇IF (w∗) ∥∥∥∥∥ 2
+ 4E ∥∥∥∥∥ 1|B|∑ it∈B ∇Ifit(w̃)−∇IF (w̃)−∇Ifit(w∗) +∇IF (w∗) ∥∥∥∥∥ 2 + 4E‖∇IF (wt)−∇IF (w∗)‖2 + 4‖∇IF (w∗)‖2
= 4 n− |B| |B|(n− 1) E‖∇Ifi(wt)−∇IF (wt)−∇Ifi(w∗) +∇IF (w∗)‖2
+ 4 n− |B| |B|(n− 1) E‖∇Ifi(w̃)−∇IF (w̃)−∇Ifi(w∗) +∇IF (w∗)‖2
+ 4E‖∇IF (wt)−∇IF (w∗)‖2 + 4‖∇IF (w∗)‖2
≤ (4 + 4 n− |B| |B|(n− 1) )E‖∇Ifi(wt)−∇Ifi(w∗)‖2 + 4‖∇IF (w∗)‖2
+ 4 n− |B| |B|(n− 1) E‖∇Ifi(w̃)−∇Ifi(w∗)‖2
≤ 16δρ+ŝ (1 + n− |B| |B|(n− 1) )E[F (wt)− F (w∗)] + 4‖∇IF (w∗)‖2
+ 16δρ+ŝ n− |B| |B|(n− 1) E[F (w̃)− F (w∗)]
(17)
where the first inequality follows from ‖a+ b+ c+ d‖2≤4‖a‖2+4‖b‖2+4‖c‖2+4‖d‖2, the third equality follows from Lemma 4, the second inequality follows from E‖x−Ex‖2 ≤ E‖x‖2 and (16),
and the last inequality follows from (16). Due to the fact that E‖∇SgI(wt)‖2 ≤ 1/kE‖∇gI(wt)‖2 where I = supp(w∗) ∪ supp(w̃r+1) ∪ supp(w̃r) and S = Gjt ∪ supp(w̃r), we have
E‖∇SgI(wt)‖2 ≤ 1
k (16δρ+ŝ (1 + n− |B| |B|(n− 1) )E[F (wt)− F (w∗)] + 4‖∇IF (w∗)‖2
+ 16δρ+ŝ n− |B| |B|(n− 1)
E[F (w̃)− F (w∗)]), (18)
where k is the number of blocks. Note that we take expectation with respect to randomized block in the above equation. This completes the proof.

heading: B.3 PROOF OF LEMMA 6
text: Proof. By the definition∇g(wt) = ∇fit(wt)−∇fit(w̃)+Dit∇F (w̃), we have E[Dit∇F (w̃)] = ∇F (w̃). Therefore, ∇g(wt) is an unbiased estimator of∇F (wt). Based on (16), we have the following result:
E‖∇gI(wt)‖2 = E‖∇Ifit(wt)−∇Ifij (w̃) +Dit∇IF (w̃)‖2
≤ 3E‖[∇Ifit(w̃)−∇Ifit(w∗)]−∇IF (w̃) +∇IF (w∗) +∇IF (w̃)−Dit∇IF (w̃)‖2
+ 3E‖∇Ifit(wt)−∇Ifit(w∗)‖2 + 3‖∇IF (w∗)‖2
≤ 6E‖[∇Ifit(w̃)−∇Ifit(w∗)]−∇IF (w̃) +∇IF (w∗)‖2
+ 6E‖∇IF (w̃)−D∇IF (w̃)‖2
+ 3E‖∇Ifit(wt)−∇Ifit(w∗)‖2 + 3‖∇IF (w∗)‖2
≤ 12δρ+ŝ [F (w t)− F (w∗)] + (3 + 12(D2m − 1))‖∇IF (w∗)‖2
+ (24δρ+ŝ + 48δρ + ŝ (D 2 m − 1))[F (w̃)− F (w∗)],
(19)
where the first inequality follows from ‖a+b+c‖2 ≤ 3‖a‖2 +3‖b‖2 +3‖c‖2, the second inequality follows from ‖a+ b‖2 ≤ 2‖a‖2 + 2‖b‖2, the third inequality follows from E‖w− Ew‖2 ≤ E‖w‖2 with E[∇Ifit(w̃) − ∇Ifit(w∗)] = ∇IF (w̃) − ∇IF (w∗), and E‖∇Ifi(wt) − ∇Ifi(w∗)‖2 ≤ 4δρ+ŝ [F (w
t) − F (w∗)] with δ > 1, which is also satisfied for w̃ with δ > 1. Since Dm = maxv=1,··· ,d
1 pv , we bound the ‖∇IF (w̃)−D∇IF (w̃)‖2 term as follows:
‖(D − I)∇IF (w̃)‖2
= d∑ v=1 ( 1 p2v − 1)[∇IF (w̃)]2v
≤ (D2m − 1)‖∇IF (w̃)‖2 = (D2m − 1)‖∇IF (w̃)−∇IF (w∗) +∇IF (w∗)‖2 ≤ 2(D2m − 1)‖∇IF (w̃)−∇IF (w∗)‖2 + 2(D2m − 1)‖∇IF (w∗)‖2 ≤ 8δρ+ŝ (D 2 m − 1)[F (w̃)− F (w∗)] + 2(D2m − 1)‖∇IF (w∗)‖2.
(20)
Due to the fact that E‖∇SgI(wt)‖2 ≤ 1kE‖∇gI(w t)‖2, we have
E‖∇SgI(wt)‖2 ≤ 1
k (12δρ+ŝ [F (w t)− F (w∗)] + (3 + 12(D2m − 1))‖∇IF (w∗)‖2
+ (24δρ+ŝ + 48δρ + ŝ (D 2 m − 1))[F (w̃)− F (w∗)]),
(21)
where k is the number of blocks. Note that we take expectation with respect to randomized block in the above equation. This completes the proof.

heading: C PROOFS OF THEOREM 1 AND COROLLARIES 1 AND 2
text: In this section, we give the detailed proofs of Theorem 1 and Corollaries 1 and 2, which can guarantee the convergence properties of our SBCD-HTP.

heading: C.1 PROOF OF THEOREM 1
text: Proof. In this subsection, we provide the proof of Theorem 1. Let wt+1 = wt − η∇SgI(wt) and I = I∗∪Ir∪Ir+1, where I∗ = supp(w∗), Ir = supp(w̃r) and Ir+1 = supp(w̃r+1). Conditioning on wt, we have the following expectation
E‖wt+1 − w∗‖2
= E‖wt − η∇SgI(wt)− w∗‖2
= E‖wt − w∗‖2 + η2E‖∇SgI(wt)‖2 − 2η〈wt − w∗,E∇SgI(wt)〉
≤ E‖wt − w∗‖2 + η2E‖∇SgI(wt)‖2 + 2η
k 〈w∗ − wt,∇IF (wt)〉
≤ E‖wt − w∗‖2 + η2E‖∇SgI(wt)‖2 + 2ηδ
k [F (w∗)− F (wt)]
≤ E‖wt − w∗‖2 + 1 k (16δρ+ŝ η 2(1 + n− |B| |B|(n− 1) )E[F (wt)− F (w∗)]
+ 16δρ+ŝ η 2 n− |B| |B|(n− 1) E[F (w̃)− F (w∗)]− 2ηδ[F (wt)− F (w∗)]) + 4η 2 k ‖∇IF (w∗)‖2
= E‖wt − w∗‖2 + 1 k ((16δρ+ŝ η 2(1 + n− |B| |B|(n− 1) )− 2ηδ)E[F (wt)− F (w∗)]
+ 16δρ+ŝ η 2 n− |B| |B|(n− 1) E[F (w̃)− F (w∗)]) + 4η 2 k ‖∇IF (w∗)‖2
(22)
where the first inequality follows from the expectation with respect to randomized block, the second inequality follows from our extended restricted strong convexity in the proof of Lemma 5, and the third inequality holds by Lemma 5. Notice that w̃ = w0 = w̃r−1. Summing (22) over t = 0, · · · ,m− 1 and taking expectation with respect to all randomness, we have
E‖wm − w∗‖2 ≤ E‖w̃r−1 − w∗‖2 + 1 k ((16δρ+ŝ η 2m(1 + n− |B| |B|(n− 1) )− 2mηδ)E[F (wm)− F (w∗)]
+ 16δρ+ŝ η 2m n− |B| |B|(n− 1) E[F (w̃r−1)− F (w∗)]) + 4η 2m k ‖∇IF (w∗)‖2
(23) Notice that the convergence property used here is similar to the proofs of the convergence results for many hard thresholding algorithms (Chen & Gu, 2016; 2017; Li et al., 2016b).
Moreover, we can obtain
E‖wm − w∗‖2 ≤ E‖w̃r−1 − w∗‖2 + 1 k ((16δρ+ŝ η 2m(1 + n− |B| |B|(n− 1) )− 2mηδ)E[F (wm)− F (w∗)]
+ 16δρ+ŝ η 2m n− |B| |B|(n− 1) E[F (w̃r−1)− F (w∗)]) + 4η 2m k ‖∇IF (w∗)‖2
= E‖w̃r−1 − w∗‖2 + 1 k ((16δρ+ŝ η 2m(1 + n− |B| |B|(n− 1) )− 2mηδ)E[F (w̃r)− F (w∗)]
− (16δρ+ŝ η 2m(1 + n− |B| |B|(n− 1) )− 2mηδ)E[F (w̃r)− F (wm)]
+ 16δρ+ŝ η 2m n− |B| |B|(n− 1) E[F (w̃r−1)− F (w∗)]) + 4η 2m k ‖∇IF (w∗)‖2
≤ E‖w̃r−1 − w∗‖2 + 1 k ((16δρ+ŝ η 2m(1 + n− |B| |B|(n− 1) )− 2mηδ)E[F (w̃r)− F (w∗)]
− (16δρ+ŝ η 2m(1 + n− |B| |B|(n− 1) )− 2mηδ) ρ+ŝ (d− s) 2σ(d− s∗) E‖wm − w∗‖2
+ 16δρ+ŝ η 2m n− |B| |B|(n− 1) E[F (w̃r−1)− F (w∗)]) + 4η 2m k ‖∇IF (w∗)‖2
(24)
where the second inequality holds due to the fact that we have the following results: When F (w̃r)− F (wm) > 0, by using Lemma 3 on wm and w̃r, we have
〈∇F (wm), w̃r − wm〉 ≤ F (w̃r)− F (wm)− 1 2ρ+ŝ ‖∇F (w̃r)−∇F (wm)‖2
≤ (1− σ)[F (w̃r)− F (wm)], (25)
where the last inequality holds due to the fact that we assume that there exists a constant factorσ>0 making this inequality true. Then by using the RSS condition, we have
F (w̃r)− F (wm) ≤ 〈∇F (wm), w̃r − wm〉+ ρ+ŝ 2 ‖w̃r − wm‖2
≤ (1− σ)[F (w̃r)− F (wm)] + ρ+ŝ 2 ‖w̃r − wm‖2.
(26)
After simplification, we have F (w̃r)− F (wm) ≤ ρ + ŝ
2σ ‖w̃ r − wm‖2. Using Lemma 1, we have
‖w̃r − wm‖2 = ‖Hk(wm)− wm‖2 ≤ |I| − s |I| − s∗ ‖wm − w∗‖2, (27)
and
F (w̃r)− F (wm) ≤ ρ+ŝ (d− s) 2σ(d− s∗) ‖wm − w∗‖2. (28)
When F (w̃r)−F (wm) < 0, we can omit−(16δρ+ŝ η 2m(1+ n−|B||B|(n−1) )−2mηδ)E[F (w̃ r)−F (wm)] directly since 16δρ+ŝ η
2m(1+ n−|B||B|(n−1) )−2mηδ is less than zero. Because this is the simple situation, we only analyze the complex one in our whole proof. Therefore, we have the following inequality:
(1 + 1
k (16δρ+ŝ η 2m(1 + n− |B| |B|(n− 1) )− 2mηδ) ρ+ŝ (d− s) 2σ(d− s∗) )E‖wm−w∗‖2
≤ E‖w̃r−1 − w∗‖2 + 1 k ((16δρ+ŝ η 2m(1 + n− |B| |B|(n− 1) )− 2mηδ)E[F (w̃r)− F (w∗)]
+ 16δρ+ŝ η 2m n− |B| |B|(n− 1) E[F (w̃r−1)− F (w∗)]) + 4η 2m k ‖∇IF (w∗)‖2
(29)
Since w̃r = Hk(wm), i.e., w̃r is the best s-sparse approximation of wm, then we have the following result due to Lemma 2,
‖w̃r − w∗‖2 ≤ (1 + 2 √ s∗√
s− s∗ )‖wm − w∗‖2. (30)
Let α = 1 + 2 √ s∗√
s−s∗ , and by combining (29) and (30), we have
(1 + 1
k (16δρ+ŝ η 2m(1 + n− |B| |B|(n− 1) )− 2mηδ) ρ+ŝ (d− s) 2σ(d− s∗) )E‖w̃r−w∗‖2
≤ αE‖w̃r−1 − w∗‖2 + α k ((16δρ+ŝ η 2m(1 + n− |B| |B|(n− 1) )− 2mηδ)E[F (w̃r)− F (w∗)]
+ 16δρ+ŝ η 2m n− |B| |B|(n− 1) E[F (w̃r−1)− F (w∗)]) + 4η 2mα k ‖∇IF (w∗)‖2
≤ 2α ρ−ŝ E[F (w̃r−1)− F (w∗)] + α k ((16δρ+ŝ η 2m(1 + n− |B| |B|(n− 1) )− 2mηδ)E[F (w̃r)− F (w∗)]
+ 16δρ+ŝ η 2m n− |B| |B|(n− 1) E[F (w̃r−1)− F (w∗)]) + 4η 2mα k ‖∇ĨF (w ∗)‖2
(31)
where the last inequality follows from the RSC condition and the definition of Ĩ, i.e., Ĩ = supp(HT (∇F (w∗), 2s)) ∪ supp(w∗).
Through proper simplification, we have the following result:
α k (2mηδ − 16δρ+ŝ η 2m(1 + n− |B| |B|(n− 1) ))E[F (w̃r)− F (w∗)]
≤ ( 2α ρ−ŝ + 16δρ+ŝ η 2mα n− |B| |B|k(n− 1) )E[F (w̃r−1)− F (w∗)] + 4η 2mα k ‖∇ĨF (w
∗)‖2 (32)
where the first term of the LHS of (31) is omitted due to the fact that the coefficient of E‖w̃r−w∗‖2 is larger than zero, when 8ρ+ŝ mη 2(1 + n−|B||B|(n−1) ) − mη +
σk(d−s∗) δρ+ ŝ (d−s) > 0, which can be naturally
satisfied under our parameter setting with s ≥ κŝs∗. Following from (32), we have
E[F (w̃r)− F (w∗)] ≤ βE[F (w̃r−1)− F (w∗)] + 2η δ(1− ω) ‖∇ĨF (w ∗)‖2 (33)
where β = k ηmρ−
ŝ δ(1−ω) +
8ρ+ ŝ η(n−|B|)
|B|(1−ω)(n−1) and ω = 8ρ + ŝ η(1+ n−|B| |B|(n−1) ). By applying (33) recursively,
then we have the following desired bound with β ≤ 12 < 1:
E[F (w̃r)− F (w∗)] ≤ (1 2 )rE[F (w̃0)− F (w∗)] + η 1− ω ‖∇ĨF (w ∗)‖2. (34)
This completes the proof.

heading: C.2 PROOF OF COROLLARY 1
text: Proof. We then give the proof for the statistical estimation analysis. Following from the RSC condition of F (·), we have
F (w∗) ≤ F (w̃r) + 〈∇F (w∗), w∗ − w̃r〉 − ρ−ŝ 2 ‖w̃r − w∗‖2. (35)
Let φ = ( 12 ) rE[F (w̃0)− F (w∗)] + η1−ω‖∇ĨF (w ∗)‖2. Combining (34) and (35), we have
E[F (w̃r)− φ] ≤ F (w∗) ≤ F (w̃r) + 〈∇F (w∗), w∗ − w̃r〉 − ρ−ŝ 2 ‖w̃r − w∗‖2. (36)
Using the duality of norms, we have
E〈∇F (w∗), w∗ − w̃r〉 ≤ ‖∇F (w∗)‖∞E‖w̃r − w∗‖1 ≤ √ ŝ‖∇F (w∗)‖∞E‖w̃r − w∗‖. (37)
Combining (36), (37) and ‖Ew‖2 ≤ E‖w‖2, we have
ρ−ŝ 2 (E‖w̃r − w∗‖)2 ≤ √ ŝ‖∇F (w∗)‖∞E‖w̃r − w∗‖+ φ. (38)
Let a = E‖w̃r − w∗‖. From the above inequality, we solve the following quadratic function with respect to a,
ρ−ŝ 2 a2 − √ ŝ‖∇F (w∗)‖∞a− φ ≤ 0, (39)
which yields the bound
E‖w̃r − w∗‖ ≤
√ 2( 12 )
r[F (w̃0)− F (w∗)] ρ−ŝ +
√ 2η
1− ω ‖∇ĨF (w
∗)‖+ 2 √ ŝ‖∇F (w∗)‖∞
ρ−ŝ
≤
√ 2( 12 )
r[F (w̃0)− F (w∗)] ρ−ŝ +
( 2
ρ−ŝ +
2η
1− ω
)√ ŝ‖∇F(w∗)‖∞.
(40)
This completes the proof.
C.3 BOUND β < 1
Now we show that we can guarantee that β ≤ 1 by providing the appropriate choices of constants with η, k and m used in the theorem. More specifically, let η ≤ C3
ρ+ ŝ
≤ 1 48ρ+
ŝ , and |B| = 1, then we have
8ρ+ŝ η(n− |B|) |B|(1− ω)(n− 1) ≤ 8C3(n− |B|) |B|(n− 1)(1− 8C3(1 + n−|B||B|(n−1) )) ≤ 1 4 . (41)
If η ≥ C2 ρ+ ŝ with C2 ≤ C3, then we have
k ηδmρ−ŝ (1− ω) ≤ 1 2δmC2 3kκŝ . (42)
If we guarantee 14δmC2 5kκŝ ≤ 14 , then we have
m ≥ 6kκŝ δC2 = C4κŝ, (43)
where δ > 1. If we choose C2 = 160 , C3 = 1 48 , C4 = 1800, k = 10 and δ = 2, then we have β ≤ 1 2 . Note that the convergence rate β does not include α, therefore, the required value of k comes from the restricted constraint of the coefficient of ‖w̃r − w∗‖2, which implys that s = Ω(κŝs∗).

heading: C.4 PROOF OF COROLLARY 2
text: Proof.
E[F (w̃r)− F (w∗)] ≤ + η 1− ω ‖∇ĨF (w ∗)‖2. (44)
Let %1, %2, · · · be a non-negative sequence of random variables which is defined as
%r = max
{ F (w̃r)− F (w∗)− η
1− ω ‖∇ĨF (w
∗)‖2, 0 } . (45)
For a fixed > 0, it follows from (34) and the Markov inequality,
P(%r ≥ ) ≤ E%r ≤ ( 12 ) r[F (w̃r)− F (w∗)]
(46)
For a given ζ ∈ (0, 1), let the right side of (46) be not greater than ζ, which requires
r ≥ log2 F (w̃r)− F (w∗)
ζ . (47)
Therefore, the result in (44) holds with probability at least 1 − ζ. Thus, we need O(log( 1 )) outer iterations to get w̃r satisfying (44). Since within each outer-iteration, we need to calculate a full gradient and m stochastic variance reduced gradients with mini-batch size |B|, the overall oracle complexity is O((n + κŝ|B|k ) log 1 ). Note that in high-dimensional region, the effect of support of snapshot point can be ignored. On the other hand, we only perform a hard thresholding operation at each outer iteration. Then we have
log2
( 1 ) = log( 1 )
log(2) = O(log(1 )). (48)
Thus, the hard thresholding oracle complexity is O(log( 1 )). This completes the proof.

heading: D PROOF OF THEOREM OF SPARSE VARIANT
text: In this section, we theoretically analyze the convergence properties of our Serial Sparse SBCD-HTP (S2BCD-HTP).
Algorithm 3 Serial Sparse SBCD-HTP (S2BCD-HTP) Input: number of outer-loops R, number of inner-loops m, step size η, sparsity level s; Initialize: w̃0;
1: for r = 0, 1, . . . , R− 1 do 2: w0 = w̃ = w̃r; 3: ∇F (w̃) = 1n ∑n i=1∇fi(w̃); 4: G̃ = supp(w̃); 5: for t = 0, 1, . . . ,m− 1 do 6: Randomly sample it from{1, 2, . . . , n} uniformly; 7: Randomly sample jt from{1, 2, . . . , k} uniformly; 8: Tit :=support of sample it; 9: S = G̃ ∪ Gjt ;
10: ∇Sg([wt]Tit ) = ∇Sfit([w t]Tit )−∇Sfit([w̃]Tit ) +Dit∇SF (w̃); 11: wt+1 = wt − η∇Sg([wt]Tit ); 12: end for 13: w̃r+1 = HT (wm, s); 14: end for Output: w̃R.

heading: D.1 PROOF
text: Proof. We start with the iterate difference between wt+1 and w∗. By expanding iterate difference and taking expectation with respect to all randomness, we get
E‖wt+1 − w∗‖2 = E‖wt − η∇SgI(wt)− w∗‖2
= ‖wt − w∗‖2 + η2‖∇SgI(wt)‖2 + 2ηE〈∇SgI(wt), w∗ − wt〉
≤ ‖wt − w∗‖2 + η2‖∇SgI(wt)‖2 + 2η
k 〈∇IF (wt), w∗ − wt〉,
(49)
where the last inequality holds by the unbiasedness of the sparse gradient estimator ∇gI(wt) and the expectation with respect to randomized block. By using the ρ−ŝ -strongly convex condition, we get the bound for 〈∇F (wt), w∗ − wt〉 as follows:
〈∇F (wt), w∗ − wt〉 ≤ F (w∗)− F (wt)− ρ−ŝ 2 ‖w∗ − wt‖2. (50)
Due to the support set I, we need the extended RSC condition with δ > 1, which has been introduced in the proof of Lemma 5. Based on this underlying consensus, we can obtain the bound for 〈∇FI(wt), w∗ − wt〉.
〈∇IF (wt), w∗ − wt〉 ≤ δ[F (w∗)− F (wt)]. (51)
Using Lemma 6, we have the following sparse variance gradient bound,
E[‖∇SgI(wt)‖2] ≤ 1
k (12δρ+ŝ [F (w t)− F (w∗)] + (3 + 12(D2m − 1))‖∇IF (w∗)‖2
+ (24δρ+ŝ + 48δρ + ŝ (D 2 m − 1))[F (w̃)− F (w∗)]),
(52)
By combining the above inequalities, we have
E‖wt+1 − w∗‖2 ≤ ‖wt − w∗‖2 + η2‖∇SgI(wt)‖2 − 2ηδ
k [F (wt)− F (w∗)]
≤ ‖wt − w∗‖2 + 12η2δρ+ŝ − 2ηδ
k [F (wt)− F (w∗)]
+ 2δη2
k (12ρ+ŝ + 24ρ + ŝ (D 2 m − 1))[F (w̃)− F (w∗)]
+ η2
k (3 + 12(D2m − 1))‖∇IF (w∗)‖2,
(53)
where the second inequality holds by Lemma 6. Summing the above inequality over t = 0, · · · ,m− 1 and taking expectation with respect to all randomness in this epoch, we have
E‖wm − w∗‖2 ≤ E‖w̃r−1 − w∗‖2 − 2ηδm(1− 6ηρ+ŝ )
k [F (wm)− F (w∗)]
+ 2δη2m
k (12ρ+ŝ + 24ρ + ŝ (D 2 m − 1))[F (w̃)− F (w∗)]
+ η2m
k (3 + 12(D2m − 1))‖∇IF (w∗)‖2.
(54)
Using (28), we obtain
E‖wm − w∗‖2 ≤ E‖w̃r−1 − w∗‖2 − 2ηδm k (1− 6ηρ+ŝ )[F (w̃ r)− F (w∗)]
+ 2ηδm(1− 6ηρ+ŝ ) ρ+ŝ (d− s)
2kσ(d− s∗) ‖wm − w∗‖2
+ δη2m
k (24ρ+ŝ + 48ρ + ŝ (D 2 m − 1))[F (w̃r−1)− F (w∗)]
+ η2m
k (3 + 12(D2m − 1))‖∇IF (w∗)‖2.
(55)
Since w̃r = HT (wm, s), i.e., w̃r is the best s-sparse approximation of wm, and due to Lemma 2, then we have following result,
2ηδαm
k (1− 6ηρ+ŝ )[F (w̃
r)− F (w∗)] + (1− 2ηδm(1− 6ηρ+ŝ ) ρ+ŝ (d− s)
2kσ(d− s∗) )‖w̃r − w∗‖2
≤ αE‖w̃r−1 − w∗‖2 + δη 2mα
k (24ρ+ŝ + 48ρ + ŝ (D 2 m − 1))[F (w̃r−1)− F (w∗)]
+ η2mα
k (3 + 12(D2m − 1))‖∇IF (w∗)‖2
≤ 2α ρ−ŝ [F (w̃r−1)− F (w∗)] + δη 2mα k (24ρ+ŝ + 48ρ + ŝ (D 2 m − 1))[F (w̃r−1)− F (w∗)]
+ η2mα
k (3 + 12(D2m − 1))‖∇IF (w∗)‖2,
(56)
where α = 1 + 2 √ s∗√
s−s∗ , and we use ρ − ŝ -strongly convex of F (w) to bound ‖w̃ r−1 −w∗‖2 in the last inequality. Through proper simplification, we can obtain
2ηαδm
k (1− 6ηρ+ŝ )[F (w̃ r)− F (w∗)]
≤ ( 2α ρ−ŝ
+ (24ρ+ŝ + 48ρ + ŝ (D 2 m − 1))δη2mα
k )[F (w̃r−1)− F (w∗)]
+ η2mα
k (3 + 12(D2m − 1))‖∇IF (w∗)‖2,
(57)
where the second term of LHS of (56) is omitted due to the fact that the coefficient of E‖w̃r −w∗‖2 is larger than zero, when 6ρ+ŝ mη 2−mη+ kσ(d−s ∗)
δρ+ ŝ
(d−s) > 0, which can be naturally satisfied under our
parameter setting with s ≥ κŝs∗. Therefore, we can obtain the following result,
E[F (w̃r)− F (w∗)] ≤ 2kα ρ− ŝ + (24 + 48(D2m − 1))δmαρ+ŝ η 2
2ηαδm(1− 6ηρ+ŝ ) · [F (w̃r−1)− F (w∗)]
+ (3 + 12(D2m − 1))mαη2
2ηαδm(1− 6ηρ+ŝ ) ‖∇ĨF (w
∗)‖2, (58)
where Ĩ = supp(w∗) ∪ supp(HT (∇F (w∗), 2s)). We define ξ := [3+12(D 2 m−1)]η
2δ(1−6ρ+ ŝ η) ‖∇ĨF (w∗)‖2 .
In order to obtain a linear convergence rate, we have to bound
β :=
2kα ρ− ŝ + (24 + 48(D2m − 1))δmαρ+ŝ η 2
2ηαδm(1− 6ηρ+ŝ ) =
k
(η − 6ρ+ŝ η2)δmρ − ŝ
+ (12 + 24(D2m − 1))ρ+ŝ η
1− 6ρ+ŝ η < 1.
(59) If we choose η ≤ 1
48D2mρ + ŝ
, then the second term of β will be less than 12 , i.e., (12+24(D2m−1))ρ + ŝ η
1−6ρ+ ŝ η ≤ 1 2 . As for the first term, when m ≥ 500D4mκŝ D2m−1 with η ≥ 1 50D2mρ + ŝ
, δ = 4 and k = 10, we have 1
(η−6ρ+ ŝ η2)δmρ− ŝ
≤ 14 . Thus, we get the linear convergence for our serial sparse SBCD-HTP (S 2BCD-
HTP).
E[F (w̃r)− F (w∗)] ≤ 3 4 · [F (w̃r−1)− F (w∗)] + ξ, (60)
which means that the total gradient evaluation oracle complexity is O((n+ κŝk log 1 ), and the hard thresholding oracle complexity is O(log( 1 )). This completes the proof.

heading: E PROOF OF THEOREM 2
text: In this section, we analyze the convergence properties of the proposed Asynchronous Sparse SBCDHTP (ASBCD-HTP). We first have to specify the iterates labeling order used in our asynchronous analysis.
We use “After Read” labeling order (Leblond et al., 2017) in our proof, which enjoys a simpler analysis but requires the order of randomly sampling step to be an uniform distributed sample. In order to give a clear proof, we adopt the “After Read” labeling order and make the following assumptions: “After Read” labeling order: 1. Inconsistently read the iterate ŵt; 2. Increase iterates counter t+ 1 and sample a random sample it+1; 3. Compute the update −η∇g(ŵt); 4. Atomic write the update to shared variable coordinately.
Assumption 4 The labeling order increases after Step 8 in Algorithm 2 finished, and thus the future perturbation is not considered in the effect of asynchrony in the current step.
Assumption 5 We assume that uniform distributed samples and the independence of the sample it+1 with ŵt.
Following (Leblond et al., 2017), we can explicit the effect of asynchrony as follows:
ŵt − wt = η t−1∑
u=(t−τ)+
Gu∇g(ŵu), (61)
where wt is the t-th inner iteration solution. In the proof of Theorem 2, we use wt to represent the temporary solution of the t-th iteration for a clean representation. And Gu is a diagonal matrix with entries in {0,+1}. This explicitly defines the coordinate perturbation from the past updates. Here, τ denotes the maximum number of overlaps between concurrent threads. We also denote ∆ = maxj=1,··· ,d pj , which provides a measure of sparsity following (Leblond et al., 2017).

heading: E.1 PROOF
text: Proof. We analyze our asynchronous algorithm based on the “perturbed iterate analysis” framework (Mania et al., 2017) with “After Read” labeling order. By expanding the iterate difference and taking expectation with respect to the sample it+1, we get
Eit+1‖wt+1 − w∗‖2 = Eit+1‖wt − η∇SgI(ŵt)− w∗‖2
= ‖wt − w∗‖2 + η2‖∇SgI(ŵt)‖2 + 2ηEit+1〈w∗ − wt,∇SgI(ŵt)〉
≤ ‖wt − w∗‖2 + η2‖∇SgI(ŵt)‖2 + 2η
k Eit+1〈w∗ − ŵt,∇gI(ŵt)〉
+ 2ηEit+1〈ŵt − wt,∇SgI(ŵt)〉,
(62)
where ŵt is the “perturbed” iterate with perturbation ψ and the last inequality follows from the expectation with respect to randomized block. We need to bound the term, 2η〈w∗ − ŵt,∇gI(ŵt)〉. Since it+1 is independent to ŵt, we can write:
Eit+1〈w∗ − ŵt,∇gI(ŵt)〉 = 〈w∗ − ŵt,Eit+1∇gI(ŵt)〉 = 〈w∗ − ŵt,∇IF (ŵt)〉. (63)
We can use ρ−ŝ -strong convexity bound of F (ŵ) with constant factor δ > 1 as well as a squared triangle inequality to get:
−〈ŵt − w∗,∇F (ŵt)〉 ≤ −δ(F (ŵt)− F (w∗))− ρ+ŝ 2 ‖ŵt − w∗‖2, (64)
− ‖ŵt − w∗‖2 ≤ ‖ŵt − wt‖2 − 1
2 ‖wt − w∗‖2, (‖a+ b‖2 ≤ 2‖a‖2 + 2‖b‖2), (65)
〈w∗ − ŵt,∇IF (ŵt)〉 ≤ −δ(F (ŵt)− F (w∗)) + ρ−ŝ 2 ‖ŵt − wt‖2 − ρ−ŝ 4 ‖wt − w∗‖2. (66)
Then we make an identical deformation, and we have
2η
k Eit+1〈w∗ − ŵt,∇gI(ŵt)〉 ≤ − ηρ−ŝ 2k Eit+1‖wt − w∗‖2 + ηρ−ŝ k Eit+1‖ŵt − wt‖2
− 2ηδ k
(F (ŵt)− F (w∗)). (67)
Putting it all together, we get the initial recursive inequity, written here explicitly:
at+1 ≤ (1− ηρ−ŝ 2k )at+η2Eit+1‖∇SgI(ŵt)‖2 + ηρ−ŝ k Eit+1‖ŵt − wt‖2
+ 2ηEit+1〈ŵt − wt,∇SgI(ŵt)〉 − 2ηδ
k êt,
(68)
where at := Eit+1‖wt − w∗‖2 and êt := Eit+1 [F (ŵt)− F (w∗)].
From Lemma 1 in (Leblond et al., 2017), we can bound the asynchronous terms ‖ŵt−wt‖2, 〈ŵt−wt and ∇gI(ŵt)〉 as follows.
E〈ŵt − wt,∇SgI(ŵt)〉 = η t−1∑
u=(t−τ)+
E〈Gu∇SgI(ŵu),∇SgI(ŵu)〉
≤ η t−1∑
u=(t−τ)+
E |〈∇SgI(ŵu),∇SgI(ŵu)〉|
≤ η t−1∑
u=(t−τ)+
√ ∆
2 (E‖∇SgI(ŵu)‖2 + E‖∇SgI(ŵt)‖2)
≤ η √ ∆
2 t−1∑ u=(t−τ)+ E‖∇SgI(ŵu)‖2 + η √ ∆ 2 E‖∇SgI(ŵt)‖2.
(69)
E‖ŵt − wt‖2 ≤ t−1∑
v,u=(t−τ)+
|〈Gu∇SgI(ŵu), Gv∇SgI(ŵv)〉|
≤ η2 t−1∑
u=(t−τ)+
‖∇SgI(ŵu)‖2 + η2 t−1∑
u,v=(t−τ)+,u6=v
|〈Gu∇SgI(ŵu), Gv∇SgI(ŵv)〉|
≤ η2 t−1∑
u=(t−τ)+
E‖∇SgI(ŵu)‖2 + η2 √ ∆(τ − 1)+ t−1∑
u=(t−τ)+
E‖∇SgI(ŵu)‖2
= η2(1 + √ ∆(τ − 1)+) t−1∑
u=(t−τ)+
E‖∇SgI(ŵu)‖2
≤ η2(1 + √ ∆τ) t−1∑ u=(t−τ)+ E‖∇SgI(ŵu)‖2,
(70) where Gu are d × d diagonal matrices whose entries in {0,+1}. Each update in ŵt is already in wt– this is the case of 0. On the contrary, some updates might be late: this is the case of +1. Then by combining (68), (69) and (70), we get
at+1 ≤ (1− ηρ−ŝ 2k )at + η2‖∇SgI(ŵt)‖2 + η3ρ−ŝ (1 + √ ∆τ) t−1∑ u=(t−τ)+ Eit+1‖∇SgI(ŵu)‖2
+ η2 √ ∆ t−1∑ u=(t−τ)+ Eit+1‖∇SgI(ŵu)‖2 + η2 √ ∆τEit+1‖∇SgI(ŵt)‖2 − 2ηδ k êt
≤ (1− ηρ−ŝ 2k )at + η2C1Eit+1‖∇SgI(ŵt)‖2 + η2C2 t−1∑
u=(t−τ)+
Eit+1‖∇SgI(ŵu)‖2 − 2ηδ
k êt,
(71)
where C1 := 1 + √ ∆τ and C2 := √
∆ + ηρ−ŝ C1. τ represents the mawimum number of overlaps between concurrent threads. And ∆ = maxj=1,··· ,d pj , which is a measure of sparsity with 1/n ≤ ∆ ≤ 1. Using Lemma 6 to bound the sparse variance term, the above inequality becomes at+1 ≤ (1− ηρ−ŝ 2k )at + η2C1 k (12δρ+ŝ [F (ŵt)− F (w ∗)] + (24 + 48(D2m − 1))δρ+ŝ [F (w̃)− F (w ∗)]
+ (3 + 12(D2m − 1))‖∇IF (w∗)‖2) + η2C2 k t−1∑ u=(t−τ)+ Eit+1‖∇gI(ŵu)‖2 − 2ηδ k êt
≤ at + η 2C2 k t−1∑ u=(t−τ)+ Eit+1‖∇gI(ŵu)‖2 − 2ηδ(1− 6ρ+ŝ ηC1) k êt
+ (24ρ+ŝ + 48ρ + ŝ (D 2 m − 1))
δC1η 2
k ẽ+ (3 + 12(D2m − 1)) η2C1 k ‖∇IF (w∗)‖2,
(72) where ẽ := E[F (w̃) − F (w∗)], êt := E[F (ŵt) − F (w∗)]. Note that we have w̃ = w0 = w̃r−1. Summing (72) over t = 0, · · · ,m− 1 and taking expectation with all randomness in this epoch, we get
am ≤ a0 − 2ηδm k (1− 6ρ+ŝ ηC1)ê m + η2C2 k m∑ t=1 t−1∑ u=(t−τ)+ E[‖∇gI(ŵu)‖2]
+ (24ρ+ŝ + 48ρ + ŝ (D 2 m − 1))
δmC1η 2
k ẽr−1 + (3 + 12(D2m − 1)) mη2C1 k ‖∇IF (w∗)‖2,
(73)
where am = E‖wm − w∗‖2, and a0 = E‖w0 − w∗‖2. Then we focus on upper bounding the third term on the RHS of (73),
m∑ t=1 t−1∑ u=(t−τ)+ E[‖∇gI(ŵu)‖2] ≤ τ m−1∑ t=1 E[‖∇gI(ŵt)‖2] ≤ τ m∑ t=1 E[‖∇gI(ŵt)‖2]
≤ τ(12ρ+ŝ δmê m + δm(24ρ+ŝ + 48ρ + ŝ (D 2 m − 1))ẽr−1 +m(3 + 12(D2m − 1))‖∇IF (w∗)‖2).
(74)
Substituting the above inequality into (73), we obtain
am ≤ a0 − 2ηδm k (1− 6ρ+ŝ η(τC2 + C1))ê m +
δρ+ŝ mη 2
k (24 + 48(D2m − 1))(τC2 + C1)ẽr−1
+ η2m
k (3 + 12(D2m − 1))(C1 + τC2)‖∇IF (w∗)‖2.
(75)
Then we bring (28) into the above inequality, we have
am ≤ a0 − 2ηδm k (1− 6ρ+ŝ η(τC2 + C1))E[F (ŵm)− F (w̃ r) + F (w̃r)− F (w∗)]
+ δρ+ŝ mη 2
k (24 + 48(D2m − 1))(τC2 + C1)ẽr−1
+ η2m
k (3 + 12(D2m − 1))(C1 + τC2)‖∇IF (w∗)‖2
≤ a0 − 2ηδm k (1− 6ρ+ŝ η(τC2 + C1))ẽ r + 2ηδm(1− 6ρ+ŝ η(τC2 + C1)) ρ+ŝ (d− s) 2kσ(d− s∗) am
+ δρ+ŝ mη 2
k (24 + 48(D2m − 1))(τC2 + C1)ẽr−1
+ η2m
k (3 + 12(D2m − 1))(C1 + τC2)‖∇IF (w∗)‖2,
(76)
where the second inequality holds since we use (28) and make a proper simplification. Using Lemma 2, we have
(1− 2ηδm(1− 6ρ+ŝ η(τC2 + C1)) ρ+ŝ (d− s)
2kσ(d− s∗) )ãr
≤ 2α ρ−ŝ ẽr−1 − 2αηδm k (1− 6ρ+ŝ η(τC2 + C1))ẽ r
+ αδρ+ŝ mη 2
k (24 + 48(D2m − 1))(τC2 + C1)ẽr−1
+ αη2m
k (3 + 12(D2m − 1))(C1 + τC2)‖∇ĨF (w ∗)‖2,
(77)
where α = 1 + 2 √ s∗√
s−s∗ , ã r = ‖w̃r − w∗‖2, and Ĩ = supp(w∗) ∪ supp(HT (∇F (w∗), 2s)).
The first term of LHS of (77) is omitted due to the coefficient of ãr is larger than zero, when 6ρ+ŝ m(τC2 + C1)η 2 −mη + kσ(d−s ∗)
δρ+ ŝ
(d−s) > 0, which can be naturally satisfied under our parameter
settings with s ≥ κŝs∗. Therefore, we obtain the following recursive inequality,
ẽr ≤ 2kα ρ− ŝ + αδρ+ŝ mη 2(24 + 48(D2m − 1))(τC2 + C1)
2αηδm(1− 6ρ+ŝ η(τC2 + C1)) · ẽr−1
+ αη2m(3 + 12(D2m − 1))(τC2 + C1)
2αηδm(1− 6ρ+ŝ η(τC2 + C1)) · ‖∇ĨF (w
∗)‖2. (78)
We assume that the following settings has satisfied with the above constraint and then by choosing m = 120κŝ, η =
1 60ρ+
ŝ
, k = 10, we get
ẽr ≤ 20 + δ 75(ρ+ ŝ )2 (1 + 2(D2m − 1))(τC2 + C1)
4δ(1− 110 (τC2 + C1)) · ẽr−1 + ξ · ‖∇ĨF (w ∗)‖2, (79)
where ξ = η(3+12(D 2 m−1))(τC2+C1)
2δ(1−6ρ+ ŝ η(τC2+C1)) . In order to ensure linear speed up, τ needs to satisfy the following constraint,
θ , 20 + δ 75(ρ+ ŝ )2 (1 + 2(D2m − 1))(τC2 + C1)
4δ(1− 110 (τC2 + C1)) ≤ 1. (80)
We assume that 37(ρ+ŝ ) 2 ≥ D2m and δ = 10, then by simply setting τ ≤ min{ 35√∆ , 2κŝ,
√ 2κŝ√
∆ }, the
above constraint is satisfied with θ ≤ 0.7538, which implies that the oracle gradient evaluation and hard thresholding oracle complexities are O((n + κŝk ) log 1 ) and O(log( 1 )), respectively. We can simply infer the statistical estimation result as before. Here, we will omit this part. This completes the proof.

heading: E.2 DISCUSSION ABOUT SPARSE VARIANCE BOUND
text: In the dense update case (i.e., Dm = 1), Lemma 6 is very similar to Lemma 5 when |B| = 1. In the sparse update case, Lemma 6 highly correlates with the sparsity of datasets (∝ D2m), which could be much loose in some extreme cases. If Dm is very large (imaging a dataset with some dimensions contain only one entry among the n samples, so Dm = n), the step size of sparse (asynchronous) SBCD-HTP will be much smaller than that of SBCD-HTP. Our algorithms can still converge to the optimal solution up to the statistical error. Actually, it is still an open problem whether we can have a tighter variance bound in the sparse update setting that is uncorrelated with Dm.

heading: F MORE EXPERIMENTAL RESULTS

heading: F.1 REAL-WORLD DATASETS AND EQUIPMENTS
text: We summarize the detailed information of all the real-world datasets in Table 2, including rcv1train, rcv1-test, real-sim, news20 and E2006-TFIDF. All these datasets are provided in the LibSVM website1. The E2006-TFIDF dataset includes 16,087 training data points and 150,360 features.
1https://www.csie.ntu.edu.tw/˜cjlin/libsvm/
Instead of using all features, we randomly select 80,000 features for training. As for the news20 dataset, we randomly choose 30,000 features for training.

heading: F.2 BASELINE METHODS
text: We compare our serial dense SBCD-HTP method with the state-of-the-art methods:
• Fast Gradient with Hard Thresholding (FG-HT) (Yuan et al., 2014): This method is based on a standard gradient descent step and combined with hard thresholding to solve sparsity-constrained problems.
• Stochastic Gradient with Hard Thresholding (SG-HT) (Nguyen et al., 2017): This method combines stochastic gradient descent with hard thresholding for solving large-scale sparsity-constrained problems.
• Stochastic Variance Reduced Gradient with Hard Thresholding (SVRG-HT) (Li et al., 2016b): This method incorporates the variance reduction technique in (Johnson & Zhang, 2013) into the hard thresholding algorithm to improve the efficiency of stochastic optimization.
• Accelerated Stochastic Block Coordinate Gradient Descent with Hard Thresholding (ASBCDHT) (Chen & Gu, 2016): Block coordinate descent (BCD) is used for solving sparsity-constrained problems.
• Fast Newton Hard Thresholding Pursuit (FNHTP) (Chen & Gu, 2017): This method tries to iteratively approximate the inverse Hessian matrix and combines a Newton algorithm with hard thresholding.
For the proposed Serial Sparse SBCD-HTP (S2BCD-HTP) and Asynchronous sparse parallel SBCD-HTP (ASBCD-HTP) methods, we compare them with the following asynchronous parallel sparsity-constrained methods:
• Asynchronous Accelerated Stochastic Block Coordinate Descent Hard Thresholding (A2SBCD-HT): This algorithm is the asynchronous variant of ASBCDHT by using our parallel framework.
• Asynchronous Stochastic Gradient with Hard Thresholding (ASG-HT): This method is an extension of SG-HT from serial dense setting to asynchronous sparse setting.
• Asynchronous Stochastic Variance Reduced Gradient Hard Thresholding (ASVRGHT) (Li et al., 2016a): This method incorporates multicore structure to compute stochastic variance reduced gradient in a parallel way.

heading: F.3 FACE RECOGNITION DATASETS
text: Although there are many datasets available for face recognition, we choose a common dataset (i.e. the Extended Yale B database (Georghiades et al., 2001)). The Extended Yale B database contains 2,414 frontal-face images of 38 people under different controlled lighting conditions (Georghiades et al., 2001). For each individual, we randomly choose 26 images for training and 15 images for testing.

heading: F.4 EXPERIMENTAL SETUP
text: For each sparsity-constrained algorithm, the sparse level makes great difference to the solution of sparse representation, especially at different noise levels. Therefore, in order to approach the best performance of all these algorithms, we change the sparsity parameter within a certain range for each algorithm. Thus, we can make sure that all these algorithms achieve the best recognition rates in the parameter setting. In all settings of the experiments, the images are down-sampled to 32×32 pixels.
ϑ = 0.3 (left) and ϑ = 0.6 (right).
rcv1-train data.
Based on the sparse representation-based classification algorithm (SRC) (Wright et al., 2008), a series of processing operations are made to the above dataset. We first rescale the training matrix into [0,1] for the convenience of adding noise, and then add Gaussian noise with zero mean and standard deviation ϑ. Finally, we normalize the columns of the training matrix to have unit `2-norm. The parameter settings of comparison algorithms for face recognition task are the same as those described in Section 6.1.

heading: F.5 RESULTS
text: We first measure the effect of mini-batch size |B| in the convergence of SBCD-HTP on rcv1-train, as shown in Figure 6. It is obvious that the large mini-batch size can accelerate the convergence rate and make SBCD-HTP to obtain a better solution. This result is exactly in line with our theoretical analysis. However, when the mini-batch size is relatively large, it will deteriorate the performance of running time. One can simply choose the mini-batch size to 3 or 5 as suggested in Figure 6.

heading: F.6 FACE RECOGNITION TASKS
text: Then, we compare the performance of some hard thresholding algorithms, including SG-HT, SVRGHT and ASBCDHT, with our SBCD-HTP for solving face recognition tasks on the Extended Yale B dataset. As shown in Figure 5, SBCD-HTP can achieve the best performance among these stochastic hard thresholding methods in a short time. It is clear that SBCD-HTP can not only optimize the objective function efficiently but also obtain a better representative solution on real applications, which is consistent with our theoretical results as described in Section 4.2. This result also shows that our SBCD-HTP has a better generalization accuracy than other hard thresholding algorithms.
vertical axis shows the objective value minus the minimum, and the horizontal axis is the number of effective
passes over data or running time (seconds).

heading: F.7 LINEAR/LOGISTIC REGISTRATION
text: We show more practical evaluation results for sparse linear regression and sparse logistic regression problems. All hardware environment and parameter settings are the same as in Section 6. For SBCD-HTP, from Figures 7 and 8, we can come to the same conclusion as we did before, that is in both sparse linear regression and sparse logistic regression, our algorithm converges faster than the baseline algorithms in terms of effective passes and running time. However, we find an interesting
result, which is that for the E2006-TFIDF dataset, all the algorithms converge with high precision solutions, so that we have to use some methods to magnify their differences. Of course, our SBCDHTP has the best performance, but the gap is not that big. We think that this is due to ease with which E2006-TFIDF dataset can be optimized, since this phenomenon does not occur on other datasets.
For S2BCD-HTP and ASBCD-HTP, from Figure 9, our S2BCD-HTP and ASBCD-HTP far exceed the baseline algorithms in terms of running time. Therefore, we can confirmedly obtain the conclusion that our Sparse and Asynchronous algorithms outperform the state-of-the-art hard thresholding algorithms for sparse datasets. This is a significant improvement for hard thresholding algorithms to solve sparse learning problems in the high-dimensional asynchronous setting.

heading: F.7.1 UPPER BOUND
text: In order to conduct an experiment on the upper bound of Theorem 1. We first generate some n×d synthetic matrices X , each row of which is drawn independently from a d-dimensional Gaussian distribution with mean 0 and covariance matrix Σ ∈Rd×d. The response vector is generated from the model y=Xw∗+ , where w∗∈Rd is the s∗-sparse coefficient vector, and we need to generate the noise drawn from a multivariate normal distribution N(0, σ2I) with σ2 = 0.01. The nonzero entries in w∗ are sampled independently from a uniform distribution over the interval [−1, 1]. For the experiments, we construct the following synthetic data set: n = 5000, d = 10000, s∗= 500,
and the diagonal entries of the covariance matrix Σ are set to 1, and the other entries are set to 0.1.
sparse linear regression problems. The four asynchronous parallel algorithms (i.e., ASG-HT, ASVRH-HT,
A2SBCD-HT and ASBCD-HTP) run on 20 threads.
We set the step-size η = 148ρŝ and batch size |B| = 1 during our experiments. The upper bound is then computed according to Theorem 1. Under this condition, we obtain the result, as is shown in Figure 10. As we can see that the upper bound converges faster at first and then it becomes slow since the statistical error becomes the main factor. The convergence performance of SBCD-HTP is always bounded below the black line, which demonstrates the correctness of our theorem.

heading: F.8 LINEAR SVM
text: In this subsection, we consider the sparsity-constrained L2-SVM problem (i.e., minw 1 2n ∑n i=1 ( max{0, 1− y(i)w>x(i)} )2 + λ2 ‖w‖
2, subject to ‖w‖0 ≤ s), where y(i) is the known label and x(i) is the given sample. We fix the regularization parameter λ = 10−4. The other parameters are also same to the above experiments. The objective value gap and running time of the baselines are provided in Figure 11. We can see that we also make a similar observation as from the previous results. SBCD-HT performs the best among these considered baselines.
problems.

Review summaries:
This paper introduces a new optimization algorithm for convex functions with sparsity constraints, based on hard thresholding to enforce sparsity. It is argued that applying hard thresholding in every step is the computational bottleneck of previous works based on hard thresholding, and the paper proposes a variant that alleviates previous works’ convergence problems. This is achieved by having an inner and outer loop in which the hard thresholding is applied only in the outer loop, and in the inner loop block coordinate descent is used  (ie. applying updates to only a subset of the solution’s variables). The paper shows that the convergence bounds of the algorithms compares favourably to previous hard thresholding algorithms convergence bounds, and empirical results show the effectiveness of the algorithm compared to previous ones.