Title: EFFICIENT HIGH-DIMENSIONAL DATA REPRESENTA-TION LEARNING VIA SEMI-STOCHASTIC BLOCK CO-ORDINATE DESCENT METHODS

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.

1. INTRODUCTION:

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 (Bü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:

EQUATION

where F (w) is a finite-sum convex and smooth function, each function f i (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

Table 1 : Comparison of some properties of the hard thresholding algorithms for solving sparsity-constrained problems. |B| is the mini-batch size of block coordinate descent methods, and k is the number of coordinate blocks. Since the sparsity level s = Ω(κ s s * ) required in our algorithms (i.e., SBCD-HTP and its parallel variant, ASBCD-HTP) is much smaller than those of other algorithms (e.g., Ω(κ 2 s s * )), the condition number κ s in our algorithms is also smaller than κ s in other algorithms, where s = 2Ω(κ s s * )+s * and s= 2Ω(κ 2 s s * )+s * . Î = supp(HT (∇F (w * ), 2s)) ∪ supp(w * ) with s = Ω(κ 2 s s * ) is a support set. 

Methods:



Required:

O(κ s log( 1 )) O( √ s ∇F (w * ) ∞ + ∇ Î F (x * ) ) ASBCDHT Ω(κ 2 s s * ) O((n+ κ s |B| k )log( 1 )) O(κ s log( 1 )) O( √ s ∇F (w * ) ∞ ) SBCD-HTP Ω(κ s s * ) O((n+ κ s |B| k )log( 1 )) O(log( 1 )) O( √ s ∇F (w * ) ∞ ) ASBCD-HTP Ω(κ s s * ) O((n+ κ s k )log( 1 )) O(log( 1 )) O( √ s ∇F (w * ) ∞ ) S 2 BCD-HTP Ω(κ s s * ) O((n+ κ s k )log( 1 )) O(log( 1 )) O( √ s ∇F (w * ) ∞ )

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 (ASBCD-HTP). 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 = Ω(κ 2 s 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., 1 n n i=1 ∇f i (w * ) ) due to the large magnitude of ∇f i (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Ω(κ 2 s 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 S 2 BCD-HTP) are much lower than those of the state-of-the-art hard thresholding algorithms, and that of S 2 BCD-HTP is similar to ASBCD-HTP under conditions on the delay. Actually, the gradient oracle complexity of S 2 BCD-HTP is better than SBCD-HTP when we deal with sparse data sets, since S 2 BCD-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) .

2. RELATED WORK:

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; Konečnỳ 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; Konečnỳ 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.

3. PRELIMINARIES:

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 ∈ R d 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 I C as the complement set of I, and v I ∈ R d , where 

[v I ] j = v j if j ∈ I and [v I ] j = 0 if j / ∈ I. Given an integer

w 0 = w = w r , ∇F ( w) = 1 n n i=1 ∇f i ( 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 j t from [k] uniformly, and S = G ∪ G jt ;

6:

∇ S g(w t ) = 1 |B| it∈B ∇ S f it (w t ) -∇ S f it ( w) + ∇ S F ( w); 7:

w t+1 = w t -η∇ S g(w t );

8:

end for 9:

w r+1 = HT (w m , 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 ρ - sstrongly convex at sparsity level s, if there exists a uniform constant ρ - s > 0 such that for any w, w ∈ R d with w -w 0 ≤ s, we have

EQUATION

)

Assumption 2 (Restricted Strong Smoothness) For any i ∈ [n], the differentiable function f i (•) is restricted ρ + s -strongly smooth at sparsity level s, if there exists a uniform constant ρ + s > 0 such that for any w, w ∈ R d with w -w 0 ≤ s, we have

EQUATION

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.

4. SEMI-STOCHASTIC BLOCK COORDINATE DESCENT HARD THRESHOLDING PURSUIT:

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.

4.1. OUR SBCD-HTP ALGORITHM:

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 G jt , where j t ∈ {1, 2, . . . , k}, {G 1 . . . G k } 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 G jt . 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 w t+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 G jt 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.

4.2. CONVERGENCE ANALYSIS:

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 f i (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

EQUATION

)

where Ĩ = 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 ( 1 2 ) r E[F ( w 0 ) -F (w * )] ≤ , we can obtain the following corollary.

Corollary 1 Suppose the conditions in Theorem 1 hold. To achieve ( 1 2 ) r E[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 G jt 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)) and O((n+κ s ) log( 1)), respectively), SBCD-HTP has a much lower oracle complexity. For ASBCDHT (Chen & Gu, 2016) , whose gradient oracle complexity is O((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 AS-BCDHT, 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 = Ω(κ 2 s 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

EQUATION

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 = Ω(κ 2 s 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 * ) ∞ + ∇ Ĩ F (w * ) ) with s = 2Ω(κ 2

s 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( 1 n n i=1 ∇f i (w * ) ))

, since the magnitude of the individual gradient norm ∇f i (w * ) can still be relatively large (Zhou et al., 2018b) .

5. ASYNCHRONOUS VARIANT OF SBCD-HTP:

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 T it of sample i t . Let P it be the projection matrix for the support T it , we can get D it since D it = P it D. Then we can find that E it [D it ∇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] Ti t ) to the intersection set of S and T it and obtain ∇ S g([ w] Ti t ). 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 SBCD-HTP (S 2 BCD-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 S 2 BCD-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: T it :=support of sample i t , and t = t + 1; //atomic increase counter t 8: 

EQUATION

∇F ( w) = 1 n n i=1 ∇f i ( w); //computed

S = G ∪ G jt , [ w] Ti t := inconsistent read of shared variable [ w r ] Ti t ; 9: ∇ S g([ w] Ti t ) = ∇ S f it ([ w] Ti t ) -∇ S f it ([ w] Ti t ) + D it ∇ S F ( w); 10: [ w r ] S Ti t = [ w r ] S Ti t -η∇ S g([ w]

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 f i (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

EQUATION

where τ denotes the maximum number of concurrent threads (Mania et al., 2017) and ∆ = max j=1•••d p j , which is a indicator to measure the sparsity of datasets (Leblond et al., 2017) .

Theorem 2 implies that ASBCD-HTP (or S 2 BCD-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.

6. EXPERIMENTS:

In this section, we evaluate the performance of SBCD-HTP, S 2 BCD-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 w r -w r+1 / w r ≤ 10 -6 , and then we use the minimum F (w r ) as the approximate optimal value F * for sub-optimality estimation in our experimental results Zhou et al. (2018b) .

6.1. DENSE ALGORITHMS:

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.

6.2. ASYNCHRONOUS ALGORITHMS:

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 A 2 SBCD-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 S 2 BCD-HTP significantly outperform than other algorithms. Moreover, ASBCD-HTP has a faster convergence performance than S 2 BCD-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 G jt . 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.

7. CONCLUSION:

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.

APPENDIX A KEY LEMMAS:

Lemma 1 ( (Jain et al., 2014 )) For any index set I, any w ∈ R I and HT (•, s) : R d → R d 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 * ∈ R I 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 * ∈ R d be the optimal sparse vector such that w * 0 ≤ s * and HT (•, s) : R d → R d 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 ∈ R d , 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 ∈ R d , 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 f i (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

EQUATION

)

Proof. See Section B.1 for the proof of Lemma 4.

Lemma 5 Suppose F (w) is ρ - s -strongly convex and each function f i (w) is ρ + s -strongly smooth with parameter s = 2s + s * . Let w * ∈ R d 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 w t , w ∈ R d and the sample i t , denote ∇ S g(w

t ) = 1 |B| it∈B ∇ S f it (w t )-∇ S f it ( w)+ ∇ S F ( w)

, then we can bound E ∇ S g I (w t ) 2 as follows:

E ∇ S g I (w t ) 2 ≤ 1 k (16δρ + s (1 + n -|B| |B|(n -1) )E[F (w t ) -F (w * )] + 4 ∇ I F (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 f i (w) is ρ + s -strongly smooth with parameter s = 2s + s * . Let w * ∈ R d 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

D m = max j=1,••• ,d 1 pj . For any w t , w ∈ R d and the sample i t , denote ∇ S g(w t ) = ∇ S f it (w t ) -∇ S f it ( w) + D it ∇ S F ( w),

we can bound E ∇ S g I (w t ) 2 as follows:

E ∇ S g I (w t ) 2 ≤ 1 k (12δρ + s [F (w t ) -F (w * )] + (3 + 12(D 2 m -1)) ∇ I F (w * ) 2 + (24δρ + s + 48δρ + s (D 2 m -1))[F ( w) -F (w * )])

, where p j is the probability that dimension j belonging to the support set of randomly sampling sample i t .

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.

B PROOFS OF KEY LEMMAS:

In this section, we prove three key lemmas.

B.1 PROOF OF LEMMA 4:

Proof. Let

EQUATION

)

This proof is similar to that of (Zheng & Kwok, 2016) .

B.2 PROOF OF LEMMA 5:

Proof. It is obvious that the stochastic variance reduced gradient satisfies

EQUATION

Thus, ∇g(w t ) is an unbiased estimator of the full gradient ∇F (w t ). Now, we bound E ∇g I (w t ) 2 . For any i ∈ {1, 2, • • • , n}, we define the following function

EQUATION

)

It is easy to find that ∇h i (w * ) = 0, which implies that h i (w * ) = min w h i (w). For any vector w, we have

EQUATION

)

where the second inequality follows from the RSS condition, and the second equality holds due to the fact that ∇h i (w), ∇ I h i (w) = ∇ I h i (w) 2 , and the last equality holds due to the fact that η = 1 ρ + s minimizes the function. Then, we have

EQUATION

Since the sampling i is chosen uniformly from {1, 2, • • • , n}, we have

EQUATION

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 * )) ≥ | ∇ I F (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

EQUATION

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 ∇ S g I (w t ) 2 ≤ 1/kE ∇g I (w t ) 2 where I = supp(w * ) ∪ supp( w r+1 ) ∪ supp( w r ) and S = G jt ∪ supp( w r ), we have

EQUATION

)

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.

B.3 PROOF OF LEMMA 6:

Proof. By the definition ∇g(w t ) = ∇f it (w t )-∇f it ( w)+D it ∇F ( w), we have E[D it ∇F ( w)] = ∇F ( w). Therefore, ∇g(w t ) is an unbiased estimator of ∇F (w t ).

Based on ( 16), we have the following result:

EQUATION

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[∇ I f it ( w) -∇ I f it (w * )] = ∇ I F ( w) -∇ I F (w * ), and E ∇ I f i (w t ) -∇ I f i (w * ) 2 ≤ 4δρ + s [F (w t ) -F (w * )] with δ > 1, which is also satisfied for w with δ > 1. Since D m = max v=1,••• ,d

1 pv , we bound the ∇ I F ( w) -D∇ I F ( w) 2 term as follows:

EQUATION

)

Due to the fact that E ∇ S g I (w t ) 2 ≤ 1 k E ∇g I (w t ) 2 , we have

EQUATION

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.

C PROOFS OF THEOREM 1 AND COROLLARIES 1 AND 2:

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.

C.1 PROOF OF THEOREM 1:

Proof. In this subsection, we provide the proof of Theorem 1. Let w t+1 = w t -η∇ S g I (w t ) and I = I * ∪I r ∪I r+1 , where I * = supp(w * ), I r = supp( w r ) and I r+1 = supp( w r+1 ). Conditioning on w t , we have the following expectation

E w t+1 -w * 2 = E w t -η∇ S g I (w t ) -w * 2 = E w t -w * 2 + η 2 E ∇ S g I (w t ) 2 -2η w t -w * , E∇ S g I (w t ) ≤ E w t -w * 2 + η 2 E ∇ S g I (w t ) 2 + 2η k w * -w t , ∇ I F (w t ) ≤ E w t -w * 2 + η 2 E ∇ S g I (w t ) 2 + 2ηδ k [F (w * ) -F (w t )] ≤ E w t -w * 2 + 1 k (16δρ + s η 2 (1 + n -|B| |B|(n -1) )E[F (w t ) -F (w * )] + 16δρ + s η 2 n -|B| |B|(n -1) E[F ( w) -F (w * )] -2ηδ[F (w t ) -F (w * )]) + 4η 2 k ∇ I F (w * ) 2 = E w t -w * 2 + 1 k ((16δρ + s η 2 (1 + n -|B| |B|(n -1) ) -2ηδ)E[F (w t ) -F (w * )] + 16δρ + s η 2 n -|B| |B|(n -1) E[F ( w) -F (w * )]) + 4η 2 k ∇ I F (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 = w 0 = w r-1 . Summing ( 22) over t = 0, • • • , m -1 and taking expectation with respect to all randomness, we have

E w m -w * 2 ≤ E w r-1 -w * 2 + 1 k ((16δρ + s η 2 m(1 + n -|B| |B|(n -1) ) -2mηδ)E[F (w m ) -F (w * )] + 16δρ + s η 2 m n -|B| |B|(n -1) E[F ( w r-1 ) -F (w * )]) + 4η 2 m k ∇ I F (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:

EQUATION

where the second inequality holds due to the fact that we have the following results: When F ( w r ) -F (w m ) > 0, by using Lemma 3 on w m and w r , we have

EQUATION

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

EQUATION

)

After simplification, we have

EQUATION

and

EQUATION

) When F ( w r )-F (w m ) < 0, we can omit -(16δρ + s η 2 m(1+ n-|B| |B|(n-1) )-2mηδ)E[F ( w r )-F (w m )] directly since 16δρ + s η 2 m(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:

EQUATION

Since w r = H k (w m ), i.e., w r is the best s-sparse approximation of w m , then we have the following result due to Lemma 2,

EQUATION

) Let α = 1 + 2 √ s * √

s-s * , and by combining ( 29) and ( 30), we have

(1 + 1 k (16δρ + s η 2 m(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 η 2 m(1 + n -|B| |B|(n -1) ) -2mηδ)E[F ( w r ) -F (w * )] + 16δρ + s η 2 m n -|B| |B|(n -1) E[F ( w r-1 ) -F (w * )]) + 4η 2 mα k ∇ I F (w * ) 2 ≤ 2α ρ - s E[F ( w r-1 ) -F (w * )] + α k ((16δρ + s η 2 m(1 + n -|B| |B|(n -1) ) -2mηδ)E[F ( w r ) -F (w * )] + 16δρ + s η 2 m n -|B| |B|(n -1) E[F ( w r-1 ) -F (w * )]) + 4η 2 mα k ∇ Ĩ F (w * ) 2 (31)

where the last inequality follows from the RSC condition and the definition of Ĩ, i.e., Ĩ = supp(HT (∇F (w * ), 2s)) ∪ supp(w * ).

Through proper simplification, we have the following result:

α k (2mηδ -16δρ + s η 2 m(1 + n -|B| |B|(n -1) ))E[F ( w r ) -F (w * )] ≤ ( 2α ρ - s + 16δρ + s η 2 mα n -|B| |B|k(n -1) )E[F ( w r-1 ) -F (w * )] + 4η 2 mα k ∇ Ĩ 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

EQUATION

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 β ≤ 1 2 < 1:

EQUATION

)

This completes the proof.

C.2 PROOF OF COROLLARY 1:

Proof.

We then give the proof for the statistical estimation analysis. Following from the RSC condition of F (•), we have 34) and ( 35), we have

EQUATION

) Let φ = ( 1 2 ) r E[F ( w 0 ) -F (w * )] + η 1-ω ∇ Ĩ F (w * ) 2 . Combining (

EQUATION

)

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

EQUATION

)

Let a = E w r -w * . From the above inequality, we solve the following quadratic function with respect to a, ρ -

EQUATION

which yields the bound

EQUATION

)

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

EQUATION

EQUATION

)

If we guarantee

EQUATION

where δ > 1. If we choose C 2 = 1 60 , C 3 = 1 48 , C 4 = 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 * ).

C.4 PROOF OF COROLLARY 2:

Proof.

EQUATION

)

Let 1 , 2 , • • • be a non-negative sequence of random variables which is defined as

EQUATION

)

For a fixed > 0, it follows from (34) and the Markov inequality,

EQUATION

For a given ζ ∈ (0, 1), let the right side of (46) be not greater than ζ, which requires

EQUATION

)

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

EQUATION

)

Thus, the hard thresholding oracle complexity is O(log( 1)). This completes the proof.

D PROOF OF THEOREM OF SPARSE VARIANT:

In this section, we theoretically analyze the convergence properties of our Serial Sparse SBCD-HTP (S 2 BCD-HTP).

Algorithm 3 Serial Sparse SBCD-HTP (S 2 BCD-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:

w 0 = w = w r ; 3: ∇F ( w) = 1 n n i=1 ∇f i ( w); 4: G = supp( w); 5:

for t = 0, 1, . . . , m -1 do 6:

Randomly sample i t from{1, 2, . . . , n} uniformly;

7:

Randomly sample j t from{1, 2, . . . , k} uniformly; 8:

T it :=support of sample i t ;

9:

S = G ∪ G jt ; 10: ∇ S g([w t ] Ti t ) = ∇ S f it ([w t ] Ti t ) -∇ S f it ([ w] Ti t ) + D it ∇ S F ( w);

11:

w t+1 = w t -η∇ S g([w t ] Ti t ); 12:

end for 13:

w r+1 = HT (w m , s); 14: end for Output: w R .

D.1 PROOF:

Proof. We start with the iterate difference between w t+1 and w * . By expanding iterate difference and taking expectation with respect to all randomness, we get

EQUATION

where the last inequality holds by the unbiasedness of the sparse gradient estimator ∇g I (w t ) and the expectation with respect to randomized block. By using the ρ - s -strongly convex condition, we get the bound for ∇F (w t ), w * -w t as follows:

EQUATION

)

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 ∇F I (w t ), w * -w t .

EQUATION

Using Lemma 6, we have the following sparse variance gradient bound,

EQUATION

By combining the above inequalities, we have

EQUATION

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

EQUATION

Using ( 28), we obtain

EQUATION

)

Since w r = HT (w m , s), i.e., w r is the best s-sparse approximation of w m , and due to Lemma 2, then we have following result,

EQUATION

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

EQUATION

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,

EQUATION

where Ĩ = supp(w * ) ∪ supp(HT (∇F (w * ), 2s)). We define ξ :=

[3+12(D 2 m -1)]η 2δ(1-6ρ + s η) ∇ Ĩ F (w * ) 2 .

In order to obtain a linear convergence rate, we have to bound

EQUATION

If we choose η ≤

1 48D 2 m ρ + s

, then the second term of β will be less than 1 2 , i.e.,

(12+24(D 2 m -1))ρ + s η 1-6ρ + s η ≤ 1 2 . As for the first term, when m ≥ 500D 4 m κ s D 2 m -1 with η ≥ 1 50D 2 m ρ + s

, δ = 4 and k = 10, we have

1 (η-6ρ + s η 2 )δmρ - s ≤ 1 4

. Thus, we get the linear convergence for our serial sparse SBCD-HTP (S 2 BCD-HTP).

EQUATION

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.

E PROOF OF THEOREM 2:

In this section, we analyze the convergence properties of the proposed Asynchronous Sparse SBCD-HTP (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 ŵt ; 2. Increase iterates counter t + 1 and sample a random sample i t+1 ; 3. Compute the update -η∇g( ŵ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 i t+1 with ŵt .

Following (Leblond et al., 2017) , we can explicit the effect of asynchrony as follows:

EQUATION

where w t is the t-th inner iteration solution. In the proof of Theorem 2, we use w t to represent the temporary solution of the t-th iteration for a clean representation. And G u 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 ∆ = max j=1,••• ,d p j , which provides a measure of sparsity following (Leblond et al., 2017) .

E.1 PROOF:

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 i t+1 , we get

EQUATION

where ŵ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 * -ŵt , ∇g I ( ŵt ) . Since i t+1 is independent to ŵt , we can write:

EQUATION

)

We can use ρ - s -strong convexity bound of F ( ŵ) with constant factor δ > 1 as well as a squared triangle inequality to get:

EQUATION

EQUATION

EQUATION

)

Then we make an identical deformation, and we have

EQUATION

Putting it all together, we get the initial recursive inequity, written here explicitly:

EQUATION

where a t := E it+1 w t -w * 2 and êt :

= E it+1 [F ( ŵt ) -F (w * )].

From Lemma 1 in (Leblond et al., 2017) , we can bound the asynchronous terms ŵt -w t 2 , ŵt -w t and ∇g I ( ŵt ) as follows.

EQUATION

) E ŵt -w t 2 ≤ t-1 v,u=(t-τ )+ | G u ∇ S g I ( ŵu ), G v ∇ S g I ( ŵv ) | ≤ η 2 t-1 u=(t-τ )+ ∇ S g I ( ŵu ) 2 + η 2 t-1 u,v=(t-τ )+,u =v | G u ∇ S g I ( ŵu ), G v ∇ S g I ( ŵv ) | ≤ η 2 t-1 u=(t-τ )+ E ∇ S g I ( ŵu ) 2 + η 2 √ ∆(τ -1) + t-1 u=(t-τ )+ E ∇ S g I ( ŵu ) 2 = η 2 (1 + √ ∆(τ -1) + ) t-1 u=(t-τ )+ E ∇ S g I ( ŵu ) 2 ≤ η 2 (1 + √ ∆τ ) t-1 u=(t-τ )+ E ∇ S g I ( ŵu ) 2 ,

(70) where G u are d × d diagonal matrices whose entries in {0, +1}. Each update in ŵt is already in w t -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

EQUATION

where C 1 := 1 + √ ∆τ and C 2 := √ ∆ + ηρ - s C 1 . τ represents the mawimum number of overlaps between concurrent threads. And ∆ = max j=1,••• ,d p j , which is a measure of sparsity with 1/n ≤ ∆ ≤ 1.

Using Lemma 6 to bound the sparse variance term, the above inequality becomes 

EQUATION

EQUATION

)

where a m = E w m -w * 2 , and a 0 = E w 0 -w * 2 . Then we focus on upper bounding the third term on the RHS of ( 73),

EQUATION

Substituting the above inequality into (73), we obtain

EQUATION

)

Then we bring (28) into the above inequality, we have

EQUATION

)

where the second inequality holds since we use (28) and make a proper simplification. Using Lemma 2, we have

EQUATION

where α = 1 + 2 √ s * √ s-s * , ãr = w r -w * 2 , and Ĩ = supp(w * ) ∪ supp(HT (∇F (w * ), 2s)). The first term of LHS of ( 77) is omitted due to the coefficient of ãr is larger than zero, when

6ρ + s m(τ C 2 + C 1 )η 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,

EQUATION

)

We assume that the following settings has satisfied with the above constraint and then by choosing

EQUATION

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,

EQUATION

)

We assume that 37(ρ + s ) 2 ≥ D 2 m and δ = 10, then by simply setting τ ≤ min{ 3 5 √ ∆ , 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 ) log1 ) 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.

E.2 DISCUSSION ABOUT SPARSE VARIANCE BOUND:

In the dense update case (i.e., D m = 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 (∝ D 2 m ), which could be much loose in some extreme cases. If D m is very large (imaging a dataset with some dimensions contain only one entry among the n samples, so D m = 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 D m . 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.

F MORE EXPERIMENTAL RESULTS:



F.2 BASELINE METHODS:

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. (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 (S 2 BCD-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 (A 2 SBCD-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 (ASVRG-HT) (Li et al., 2016a) : This method incorporates multicore structure to compute stochastic variance reduced gradient in a parallel way.

F.3 FACE RECOGNITION DATASETS:

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.

F.4 EXPERIMENTAL SETUP:

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. 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.

F.5 RESULTS:

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 .

F.6 FACE RECOGNITION TASKS:

Then, we compare the performance of some hard thresholding algorithms, including SG-HT, SVRG-HT 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. 

F.7 LINEAR/LOGISTIC REGISTRATION:

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 SBCD-HTP 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 S 2 BCD-HTP and ASBCD-HTP, from Figure 9 , our S 2 BCD-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. We set the step-size η = 1 48ρ 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.

F.8 LINEAR SVM:

In this subsection, we consider the sparsity-constrained L 2 -SVM problem (i.e., min w 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. 

:

https://www.csie.ntu.edu.tw/ ˜cjlin/libsvm/

Figures and Tables:

rcv1-train Figure1: Comparison of all the hard thresholding algorithms for solving sparse logistic regression problems. In each plot, the vertical axis shows the objective value minus the minimum, and the horizontal axis is the number of effective passes over data (left) or running time (seconds, right).

rcv1-train Figure2: Comparison of all the hard thresholding algorithms for solving sparse linear regression problems.

Figure 3: Comparison of ASG-HT, ASVRG-HT, A 2 SBCD-HT, S 2 BCD-HTP and ASBCD-HTP for solving sparse logistic regression problems. The four asynchronous parallel algorithms (i.e., ASG-HT, ASVRG-HT, A 2 SBCD-HT and ASBCD-HTP) run on 20 threads.

) where ẽ := E[F ( w) -F (w * )], êt := E[F ( ŵt ) -F (w * )]. Note that we have w = w 0 = w r-1 . Summing (72) over t = 0, • • • , m -1 and taking expectation with all randomness in this epoch, we get

Figure 5: Recognition rates of all the algorithms on the Extended Yale B with different levels of Gaussian noise:

Figure 7: Comparison of all the algorithms for solving sparse logistic regression problems. In each plot, the

Figure 8: Comparison of all the algorithms for solving sparse linear regression problems.

Figure 10: The upper bound of Theorem 1 vs the convergence performance.

Algorithm 1 Semi-stochastic Block Coordinate Descent Hard Thresholding Pursuit (SBCD-HTP)

Table content: <table/>

Ti t ); //atomic write

Table content: <table><tr><td>11:</td><td>end while</td></tr><tr><td>12:</td><td/></tr></table>

Summary of large-scale and high-dimensional datasets

Table content: <table><tr><td>Datasets</td><td># Data</td><td># Feature</td><td>Density</td></tr><tr><td>rcv1-train</td><td>20,242</td><td>47,236</td><td>0.16%</td></tr><tr><td>rcv1-test</td><td>677,399</td><td>47,236</td><td>0.15%</td></tr><tr><td>real-sim</td><td>72,309</td><td>20,958</td><td>0.024%</td></tr><tr><td>news20</td><td>19,996</td><td>30,000</td><td>0.49%</td></tr><tr><td>E2006-TFIDF</td><td>16,087</td><td>80,000</td><td>1.54%</td></tr><tr><td>F.1 REAL-WORLD DATASETS</td><td/><td/><td/></tr></table>

• 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)

Table content: <table/>

Peer review (iclr2020):
===== Update after author response Thanks for the clarifications and edits in the paper. I recommend acceptance of the paper. Other comments: Definition 1 in the updated version is still too vague (*difference of what?* -- function values? distance in norm between iterates?) -- this should be clarified. ======== This paper considers the problem of sparsity-constrained ERM and asks whether one can design a variant of the stochastic hard thresholding approaches where the hard-thresholding complexity does not depend on a (sparsity dependent) condition number, unlike all previous approaches (Table 1). It proposes a method which combines SVRG-type variance reduction, with block-coordinate updates, leaving the hard thresholding operation outside the inner loop, to accomplish this goal. It provides a convergence analysis which significantly improves the previous best rates (by having both the sparsity level shat which is significantly lower (kappa_shat vs. kappa_stilde^2) as well as a condition number independent hard thresholding complexity (Table 1). An asynchronous and sparse (in the features) variant is also proposed, with even better complexity. Some standard experiments on sparse linear regression and sparse logistic regression is presented showing an improvement in both number of iterations as well as CPU time. I think the clarity of the paper should be quite improved (see detailed comment), hence why I think the paper is borderline, but I am leaning towards an accept given the significant theoretical improvements over the past literature (and positive empirical results), even though the algorithmic suggestion is somewhat incremental. The proposed Algorithm 1 seems very close to the one of Chen & Gu (2016), the paper should be more clear about this. There seems to be mainly two changes: a) extending the support projection of the gradient to the union of the sampled block with the one of the support of the reference parameter wtilde (vs. just the sampled block in Chen & Gu (2016) and b) moving the hard-thresholding iteration outside of the SVRG-inner loop. These small tweaks to the algorithm yield a significant theoretical improvement, though. == Detailed comments == Clarity: the number of long abbreviations with only one letter change make it hard to follow the different algorithms; perhaps a better more differentiating naming scheme could be used. Moreover, I think more background on the sparse optimization setup should be provided in the introduction or at least in the preliminaries, as I do not think the wider ICLR community is very familiar with it (in particular, no cited paper was at ICLR). For example, define early the separation in optimization error and statistical error; and point out that F(w_t) might even be lower than F(w*) as the sparsity threshold s might be much higher than s*. This will make Table 1 more concrete and less abstract for people who not are not yet experts on this particular analysis framework. - Table 1: I would suggest to put the rate for S2BCD-HTP instead on the last row and mention instead that the rate for ASBCD is similar under conditions on the delay; as it is interesting to already have a better gradient complexity for S2BCD vs. SBCD. ** Questions: ** 1) In Corollary 1, how is the gradient oracle complexity defined or computed? And more specifically, how does one compare fairly the cost of doing a gradient update in Algorithm 1 on the *bigger set* S = Gtilde U G_jt vs. just G_jt for the Chen & Gu ASBCD algorithm? Is this accounted in the computation? 2) In Figure 1, which *minimum* is referred to and how is it found? I suspect it is not F(w*) (as it could be higher than F(w_t)), i.e. it is *not* the minimum of (1) with s*. One natural guess is that it might be min_w F(w) s.t. ||w||_0 <= s, though I do not see any guarantee in the main paper that running the algorithm would make F(w_t) converge to such a value (i.e. all we know from Thm 1 is that F(w_t) might be within O(||nabla_Itilde F(w*)||^2) of F(w*) ultimately. Please explain and clarify! == Potential improvement == The current result in Theorem 1, which is building on a similar proof technique as the original SVRG paper, has the annoying property of requiring the knowledge of the condition number in setting the size of the inner loop iteration. I suspect that this is an artifact of using an outdated version of the SVRG algorithm. This has been solved since then by considering a *loopless* version of SVRG which implicitly defines the size of the inner loop in a random manner using a quantity *which does not depend on the condition number*. This was proposed first by Hofmann et al. [2015], and then re-used by Lei & Jordan [2016] and more recently by Kovalev et al. [2019] e.g. Note that Leblond et al. (2017) that you cited profusely also used this variant of SVRG. I suspect that this technique could be re-used in your case to obtain a similar result with a loopless variant (which also gives cleaner complexity results). (Though I only skimmed through your proof.) Caveat: the sensibility of the theory in the main paper seems reasonable, but I did not check the proofs in the appendix. = References: - Hofmann et al. [2015]: Variance Reduced Stochastic Gradient Descent with Neighbors, Thomas Hofmann, Aurelien Lucchi, Simon Lacoste-Julien and Brian McWilliams, NeurIPS 2015 - Lei & Jordan [2016]: Less than a Single Pass: Stochastically Controlled Stochastic Gradient Method, Lihua Lei and Michael I. Jordan, AISTATS 2016 - Kovalev et al. [2019]: Don*t Jump Through Hoops and Remove Those Loops: SVRG and Katyusha are Better Without the Outer Loop, Dmitry Kovalev, Samuel Horvath and Peter Richtarik, arXiv 2019