--- title: "Stochastic Second-Order Method for Large-Scale Nonconvex Sparse Learning Models" authors: - "Hongchang Gao" - "Heng Huang" doi: "10.24963/ijcai.2018/294" url: "https://ijcai.org/proceedings/2018/294" pdf: "https://ijcai.org/proceedings/2018/0294.pdf" venue: "Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence" year: 2018 pages: "2128-2134" track: "Main track" keywords: - "Machine Learning Applications: Big data ; Scalability" - "Machine Learning: Feature Selection ; Learning Sparse Models" --- # Stochastic Second-Order Method for Large-Scale Nonconvex Sparse Learning Models Hongchang Gao, Heng Huang *Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence*, pages 2128-2134. DOI: [10.24963/ijcai.2018/294](https://doi.org/10.24963/ijcai.2018/294) ## Abstract Sparse learning models have shown promising performance in the high dimensional machine learning applications. The main challenge of sparse learning models is how to optimize it efficiently. Most existing methods solve this problem by relaxing it as a convex problem, incurring large estimation bias. Thus, the sparse learning model with nonconvex constraint has attracted much attention due to its better performance. But it is difficult to optimize due to the non-convexity. In this paper, we propose a linearly convergent stochastic second-order method to optimize this nonconvex problem for large-scale datasets. The proposed method incorporates second-order information to improve the convergence speed. Theoretical analysis shows that our proposed method enjoys linear convergence rate and guarantees to converge to the underlying true model parameter. Experimental results have verified the efficiency and correctness of our proposed method. ## Full text ## 1 Introduction Sparse learning models, which play an important role in high dimensional machine learning applications [Lee et al., 2006; Gao et al., 2015; 2017], have attracted much attention in the past decade. Specifically, it assumes that only a few number of model parameters are responsible for the response. Thus, a straightforward way is to enforce sparsity on the model parameter by the ℓ0-norm constraint, which restricts the number of non-zero entries in the model parameter. Due to the nonconvexity of ℓ0-norm, most existing works [Tibshirani, 1996; Van de Geer, 2008; Yuan and Lin, 2006; Friedman et al., 2008; Banerjee et al., 2008] employ the relaxed ℓ1-norm regularization to enforce the sparsity of model parameters, since it is easy to solve due to the convexity of ℓ1-norm. For instance, the well-known Lasso [Tibshirani, 1996] solves a regression term to fit the data and an ℓ1-norm regularization term to pursue a sparse model parameter. However, such a convex relaxation usually degenerates the performance of the model. Thus, it is necessary and challenging to solve the ℓ0norm constraint problem directly. Formally, in this paper, we focus on the following sparsityconstrained optimization problem: where F(w) is a smooth and convex function, which measures how well the model fits the input space. ∥w∥0 ≤ s denotes the number of nonzero entries in w is not more than s, controlling the sparsity level of the model parameter. This model is very common in machine learning area. A representative case is the sparse linear regression problem, which is shown as follows: where y = [y1, · · · , yn]T ∈ Rn is the response vector, X = {x1, · · · , xn} ∈ Rn×d is the design matrix, and w ∈ Rd is the model parameter. The challenge to solve Eq. (1) is the nonconvex sparse constraint, which makes Eq. (1) as an NP-hard problem. In the past few decades, a large family of algorithms [Mallat and Zhang, 1993; Needell and Tropp, 2009; Tropp and Gilbert, 2007; Zhang, 2011] have been proposed to solve Eq. (1). Among them, there has been much progress towards the gradient-based method, such as gradient hard thresholding pursuit (GraHTP) [Yuan et al., 2014], iterative hard thresholding (IHT) [Blumensath and Davies, 2009], and so on. In particular, the gradient-based method updates the model parameter with the gradient descent method followed by HardThresholding. However, with the development of large-scale data in recent years, these algorithms fail to handle large-scale datasets. The reason is that they need to compute the gradient with respect to all data points at each iteration, making it prohibitive for large-scale datasets. To address this problem, some researchers resort to stochastic algorithms to solve Eq. (1) in recent years, such as SGHT [Nguyen et al., 2014], SVR-GHT [Li et al., 2016], ASBCDHT [Chen and Gu, 2016], and so on. Unlike the full gradient descent method whose complexity of each iteration is O(nd), stochastic methods have only O(d) complexity in each iteration so that they are efficient for large-scale datasets. Although the gradient-based method has achieved good performance when solving Eq. (1), yet it only considers the first-order information of the objective function, ignoring the second-order curvature information. As a result, its performance is far from satisfactory in some cases. For instance, when the condition number of the objective function in Eq. (1) is extremely large, the first-order method will converge very slowly. If incorporating the second-order curvature information of the objective function into the first-order gradient method, we can obtain a better searching direction in each iteration, making it converge fast. Thus, it is important to employ the second-order method to solve Eq. (1). In optimization community, there has been much progress towards the second-order method, such as [Byrd et al., 2016; Gower et al., 2016; Moritz et al., 2016; Zhao et al., 2017]. But most of them just focus on the convex problem. In [Yuan and Liu, 2014], a Newton greedy pursuit method was proposed to solve the nonconvex Eq. (1). However, it is not suitable for large-scale problems since it uses all data points in each iteration. Recently, some stochastic second-order methods have been proposed, such as [Moritz et al., 2016; Gower et al., 2016]. Although these algorithms enjoy a good converging property, the convergence analysis is based on the strongly convex condition, which is not applicable for the nonconvex Eq. (1). Thus, employing the second-order method to solve Eq. (1) is necessary and challenging. To incorporate the second-order curvature information and address the scalability problem, we propose a stochastic LBFGS method to solve the large-scale problem in Eq. (1). In particular, at each iteration, we first randomly sample a minibatch of data points to evaluate the approximated inverse Hessian matrix and the gradient, updating the model parameter without considering the sparsity constraint and then performing Hard-Thresholding on the updated model parameter to get the sparse one. Additionally, due to the stochastic sampling, the introduced variance will slow down the convergence rate. To address this problem, we incorporate the variance reduction technique as [Moritz et al., 2016]. One of the most important contributions of this paper is that we prove the linear convergence rate of the second-order stochastic L-BFGS for solving the nonconvex Eq. (1). As far as we know, this is the first work which proves stochastic L-BFGS has a linear convergence rate for the nonconvex problem with the sparsity constraint. Furthermore, the output estimator from our proposed method is guaranteed to converge to the unknown true model parameter, which is also the first work showing such a result for stochastic L-BFGS. The experiments on both synthetic and real-world datasets have shown the correctness and effectiveness of our proposed method. Notation Throughout this paper, the matrix is represented by the uppercase letter, the vector is denoted by the bold lowercase letter, and the scalar is represented by by the unbold lowercase letter. In particular, X = {x1, · · · , xn} ∈ Rn×d denotes the design matrix, y = [y1, · · · , yn]T ∈ Rn denotes the response vector, and w ∈ Rd denotes the model parameter. For the vector w ∈ Rd , we define ∥w∥1 = ∑di=1 |wi|, ∥w∥2 = ∑ni=1 w2 i , ∥w∥∞ = maxi |wi|. Additionally, supp(w) denotes the index of nonzero elements in w, and supp(w, s) denotes the index of the top s elements of w in regard to magnitude. w(t) denotes the vector in the t-th iteration. ## 2 Stochastic L-BFGS for Large-Scale Nonconvex Sparse Learning Models In this section, we will present the detail of our proposed method for large-scale nonconvex sparse learning models. The core idea is employing stochastic L-BFGS to solve Eq. (1). However, the naive stochastic L-BFGS converges slowly due to the introduced variance by random sampling. Inspired by [Moritz et al., 2016], we employ the variance reduction technique [Johnson and Zhang, 2013] to accelerate it. Meanwhile, unlike the traditional stochastic L-BFGS [Moritz et al., 2016] which is only applicable for strongly convex problems, our method can successfully handle the nonconvex problem with the sparsity constraint. The details of our proposed method are summarized in Algorithm 1. In Algorithm 1, there are two nested loops. In the outer loop, we calculate the full gradient ˜µ in Line 8 such that we can use it to reduce the variance of the stochastic gradient. In the inner loop, our algorithm combines the variance reduced gradient v(t) and the approximated inverse Hessian matrix H(r) to update the model parameter. More specifically, the gradient is an unbiased estimation to ∇F(w(t) ) = 1n ∑ni=1 ∇fi(w(t) ), where ∇fB(w(t) ) = 1 |B| ∑i∈B ∇fi(w(t) ). After that, we update the model parameter without considering the sparsity constraint as follows: where η is the step size, and ¯w(t+1) is the temporary model parameter. With this updating rule, the second-order curvature is incorporated by H(r) . Thus, it will converge faster than the first order approach. In the following, the HardThresholding is performed on ¯w(t+1) to obtain the solution satisfying the sparsity constraint as follows: where s is the sparsity level in Eq. (1), and the HardThresholding operator H(·, ·) is defined as follows: where supp(w, s) denotes the s largest non-zero values of the model parameter w. In Line 11-19 of Algorithm 1, after every L iterations, we update the approximated inverse Hessian matrix H(r) by the L-BFGS schema as follows: where r - M + 1 ≤ i ≤ r, M is the memory size, ρ(i) = 1 s(i)T y(i) , and H(r) r-M = s(r)T y( r) y(r)T y(r) I. Then, we set H(r) = H(r) r . Note that unlike traditional L-BFGS method, we update y(r) = ∇2 fB′ (θ(r) )s(r) since it works better in the stochastic setting [Moritz et al., 2016], where ∇2 fB′ (θ(r) ) = 1 |B′ | ∑i∈B′ ∇2 fi(θ(r) ). Practical Acceleration In Algorithm 1, we need to compute the Hessian matrix ∇2 fB′ (θ(r) ) and its inverse approximation H(r) . Both of them require O(d2 ) storage, which is prohibitive for high dimensionality problems. Instead of constructing H(r) explicitly, we employ the two-loop recursion method [Nocedal and Wright, 2006] to directly compute H(r)v(t) based on the correction pairs {s(i) , y(i) }ri=r-M+1 . For the Hessian matrix, we assume it can be represented as Actually, it is very common in many machine learning problems. For example, Ai(θ(r) ) in Eq. (2) is xi so that ∇2 fB′(θ(r) ) = 1 |B′ | ∑i∈B′ xixTi . Based on this representation, instead of computing ∇2 fB′(θ(r) ) explicitly, we can directly compute y(r) as follows: which will save much storage and computation since no explicit Hessian matrix needs to store. Algorithm 1 Stochastic L-BFGS Algorithm for Solving Eq. (1). Input: X ∈ Rn×d , Y ∈ Rn , s > 0. Output: w ∈ Rd ``` 1: Initialize r = 0, H0 = I 2: for k = 0, 1, 2, · · · do 3: ˜w = ˜w(k-1) , w(0) = ˜w(k-1) 4: ˜µ = 1n ∑ni=1 ∇fi( ˜w) 5: for t = 0, 1, 2, · · · ,m- 1 do 6: Randomly sample a subset B from {1, 2, · · · , n} 7: Compute gradient v(t) =∇fB(w(t) )-∇fB( ˜w)+˜µ 8: Update ¯w(t+1) =w(t) - ηH(r)v(t) 9: Hard-Thresholding w(t+1) =H( ¯w(t+1) , s) 10: if mod(t, L)=0 then 11: r = r +1 12: θ(r) = 1L ∑t-1 j=t-L w(j) 13: Randomly sample a subset B′ from {1, 2, · · · , n} 14: Compute ∇2 fB′ (θ(r) ) 15: Compute s(r) = θ(r) - θ(r-1) 16: Compute y(r) = ∇2 fB′ (θ(r) )s(r) 17: Compute H(r) with Eq. (7) 18: end if 19: end for 20: Set ˜w(k) as the randomly selected w(i) from {w(0) , · · · ,w(m-1) } 21: end for ``` - 22: return w ## 3 Convergence Analysis In this section, we will present the convergence analysis about Algorithm 1. At first, we will introduce the following assumptions that our analysis depends on. Assumption 1. (Restricted Strong Convexity) Function F satisfies restricted λ˜ s-strong convexity condition at sparse level ˜ s. Formally, we have Assumption 2. (Restricted Strong Smoothness) Function fi satisfies restricted L˜ s-strong smoothness condition at sparse level ˜ s. Formally, we have These two assumptions indicate that F(w) is strongly convex and fi(w) is smooth in the sparse subspace. Additionally, based on these two assumptions, we can define the restricted condition number as κ˜ s = L˜ s λ˜ s . Assumption 3. The gradient is bounded as follows: In the following, we present three lemmas for proving the main theorem. Lemma 1. Suppose Assumption 1 and 2 satisfy with the sparsity level ˜ s = 2s + s∗ . Then, for the sparse vector w∗ ∈ Rd such that ||w∗ ||0 ≤ s∗ , and the sparse vector w(t) ∈ Rd such that ||w(t) ||0 ≤ s, we have where S ⊇ (supp(w∗ ) ∪ supp(w(t) )). This lemma bounds the variance of the stochastic gradient v(t) S . The proof can be found in Lemma 3.5 [Li et al., 2016]. Thus, we do not include it due to the space limitation. Lemma 2. Given the sparse vector w∗ ∈ Rd such that ||w∗ ||0 ≤ s∗ , for s > s∗ and any w ∈ Rd , we have This lemma actually presents the projection error bound for the Hard-Thresholding operator. The proof can be referred to Lemma 3.3 of [Li et al., 2016]. Lemma 3. If Assumption 1 and 2 hold, the estimation of inverse Hessian matrix H(r) (for all r ≥ 1) is bounded by where 0 < γ ≤ Γ. This is a common condition for stochastic L-BFGS method, and the detailed proof can be found from Lemma 4 in [Moritz et al., 2016]. Based on these assumptions and lemmas, we turn to present the main result of our proposed method. Theorem 1. Suppose Assumption 1 and 2 satisfy with the sparsity level ˜ s = 2s + s∗ . Denote w∗ as the unknown true model parameter such that ||w∗ ||0 ≤ s∗ and ˜ S = supp(H(∇F(w∗ ), 2s)) ∪ supp(w∗ ). By choosing C2 Γ2L˜ s ≤ η ≤ C1 Γ2L˜ s and s ≥ (1 + 1/A2 )s∗ where A = C2-3κ˜ s(1+η)η 6κ˜ s(1+η)2 with valid constants C1, C2, and large m, such that where β = (1+ 2√s∗ √s-s∗ )(1+η). Then, for all k > 0, we have In the following, we present the detailed proof about Theorem 1. Proof. Due to ¯w(t+1) = w(t) - ηH (r)v(t) , conditioning on w(t) , we have where S = supp(w∗ ) ∪ supp(w(t) ) ∪ supp(w(t+1) ). The third step is due to E[H(r)v(t) S ] = H(r)∇SF(w(t) ). The fourth step follows from Assumption 1, Lemma 3, and the fact 2ab ≤ a2 + b2 . The fifth step is due to Lemma 1. According to Lemma 2, we have where β = α(1+ η) and α = 1 + 2√s∗ √s-s∗ . By summing over t = 0, · · · ,m- 1, we have where ˜w(k-1) = w(0) = ˜w and ˜ S = supp(H(∇F(w∗ ), 2s)) ∪ supp(w∗ ). By rearranging the above inequality, we have Now, we need to verify that At first, assume η ≤ C1 Γ2L˜ s < 1 18Γ2L˜ s , then Furthermore, assume s ≥ (1 + 1/A2 )s∗ where A = C2-3κ˜ s(1+η)η 6κ˜ s(1+η)2 , and η ≥ C2 Γ2L˜ s where C2 ≤ C1, to guarantee we get Hence, to guarantee θ < 1, we should have where B = C2 C2-3κ˜ s(1+η)(β-1) . In this way, we have θ < 1. By recursively applying Eq. (20), we can get the desired result as follows: which completes the proof. Remark 1. Theorem 1 indicates a linear convergence rate. In particular, to get a pre-defined accuracy ϵ > 0 with respect to the function value gap, we need O(log(1/ϵ)) outer iterations. Additionally, to have linear convergence rate, m should be set sufficiently large as in Eq. (25). In the following, we present the approximation error bound for the estimator obtained from Algorithm 1. Corollary 1. With the same conditions as Theorem 1, for all k > 0, we have the error bound for the estimator ˜w(k) as follows: Proof. Due to Assumption 1, we have Furthermore, denote then based on Eq. (16), we have Furthermore, we have Put Eq. (31) into Eq. (30), we have By solving this inequality with respect to E||w∗ - ˜w(k) ||2, we can obtain the desired result. Remark 2. This error bound for the estimator consists of three terms. The first term corresponds to the optimization error, the second term is a constant, and the third term corresponds to the statistical error. After sufficient iterations, the first term will approach to zero. Therefore, our algorithm can always converge to the unknown true parameter w∗ , up to the statistical error and a constant value. ## 4 Experiments In this section, we will present the performance of our proposed method on both synthetic and real-world datasets. Throughout the experiments, we compare it with two stateof-the-art stochastic methods. They are SGHT [Nguyen et al., 2014] and SVR-GHT [Li et al., 2016]. Specifically, - SGHT [Nguyen et al., 2014]: It employs stochastic gradient to update the model parameter and then performs Hard-Thresholding on the obtained model parameter. - SVR-GHT [Li et al., 2016]: This method adopts the variance reduced gradient to update model parameter, accelerating the converging speed. All of these methods belong to the stochastic method so that they are suitable for large-scale problems. Additionally, we set L = 10, M = 10, |B| = 10, and |B′ | = 50. The step length of each method is chosen to achieve the best performance. ## 4.1 Synthetic Data In this experiment, we focus on the sparse linear regression problem, just as shown in Eq. (2). For the synthetic data, each row of the design matrix X ∈ Rn×d is independently generated from a multivariate Gaussian distribution N(0, Σ) where Σ ∈ Rd×d . For the sparse regression coefficient w∗ , the nonzero entries are independently generated from a uniform distribution in [-1, 1]. The response vector is constructed by y = Xw∗+ϵ, where the noise ϵ is generated from a Gaussian distribution N(0, σ2 I), and we set σ2 = 0.01. With these settings, we construct two synthetic datasets. Toy1 is with n = 20000, d = 2000, s∗ = 100, s = 200, Σ = I. Toy-2 is with n = 50000, d = 5000, s∗ = 500, s = 1000, and diagonal entries of the covariance matrix Σ are set as 1, the other entries are set as 0.1. Figure 1: The objective value and estimation error of the model parameter about two synthetic datasets. ![Image](https://ijcai.org/proceedings/2018/fulltext/294/image_000000_4f241ab065ad96cef3dca50bba86e60a6fb714c6a4fff1cd2ff62cf19cf1dbce.png) In Figure 1, we show the logarithm of the objective function value and the estimation error with respect to the number of gradient evaluations. From Figure 1, we can find that our method converges faster than the state-of-the-art methods in terms of objective function value and the estimation error. In particular, although all these methods are theoretically shown with linear convergence rate for the sparse linear regression problem [Li et al., 2016], our method and SVR-GHT have better performance than SGHT due to the incorporation of the variance reduction technique. Furthermore, our method incorporates the second-order information so that it converges faster than SVR-GHT, which is consistent with the motivation of our method. ## 4.2 Real-World Data In this experiment, we will evaluate the performance of the sparse linear regression defined in Eq. (2) and the sparse logistic regression on the real-world data. Specifically, the sparse logistic regression is defined as follows: where yi ∈ {0, 1} is the class label. For the sparse linear regression model, we evaluate its performance on E2006-TFIDF dataset , which includes 16,087 training data points, 3,308 testing data points, and 150,360 features. Instead of using all features, we randomly select 20,000 features for both training set and testing set. Additionally, the sparsity level s is set as 2000. For the sparse logistic regression model, we evaluate its classification performance on the RCV1-Binary dataset , which includes 20,242 training samples and 677,399 testing samples from two classes. Additionally, it has 47,236 features totally. Here, we choose 5000 samples from each class of the testing samples as our testing set. At last, we set the sparsity level s as 500. Note that both datasets are available at the LIBSVMwebsite 1 . In Figure 2(a) and 2(b), we show the logarithm of the sparse linear regression's objective function value on the training set and testing set. We can find the similar result with the synthetic dataset. In particular, our method converges faster than the other state-of-the-art methods and has better performance on the testing set, which further confirms the effectiveness of our proposed method. In Figure 2(c) and 2(d), we show the classification accuracy of the sparse logistic regression on both training set and testing set. We can find that our method achieves the best classification result consistently during iterations. Additionally, the improvement is significant. In conclusion, our method converges fast and the obtained solution generalizes well. Figure 2: (a-b) show the logarithm of the sparse linear regression's objective function value. (c-d) show the classification accuracy of the sparse logistic regression. ![Image](https://ijcai.org/proceedings/2018/fulltext/294/image_000001_51106e066e6cd12fd89c162c6827af80a14f84dfd043e06f5072c4ad133be305.png) ## 5 Conclusion In this paper, we propose a stochastic L-BFGS method for solving large-scale nonconvex sparse learning problems. By theoretical analysis, the proposed method has shown a linear convergence rate for this kind of nonconvex problems. Meanwhile, it can guarantee to converge the underlying true model parameters. The extensive experiments have verified the efficiency of the proposed method. Thus, it can be applied to the real-world nonconvex large-scale problems. ## Acknowledgements This work was partially supported by the following grants: NSF-IIS 1302675, NSF-IIS 1344152, NSF-DBI 1356628, NSF-IIS 1619308, NSF-IIS 1633753, NIH R01 AG049371. 1 https://www.csie.ntu.edu.tw/∼cjlin/libsvmtools/datasets ## References 1. [Banerjee et al., 2008] Onureena Banerjee, Laurent El Ghaoui, and Alexandre dAspremont. Model selection through sparse maximum likelihood estimation for mul- tivariate gaussian or binary data. Journal of Machine learning research, 9(Mar):485-516, 2008. 2. [Blumensath and Davies, 2009] Thomas Blumensath and Mike E Davies. Iterative hard thresholding for compressed sensing. Applied and computational harmonic analysis, 27(3):265-274, 2009. 3. [Byrd et al., 2016] Richard H Byrd, Samantha L Hansen, Jorge Nocedal, and Yoram Singer. A stochastic quasi- newton method for large-scale optimization. SIAM Jour- nal on Optimization, 26(2):1008-1031, 2016. 4. [Chen and Gu, 2016] Jinghui Chen and Quanquan Gu. Ac- celerated stochastic block coordinate gradient descent for sparsity constrained nonconvex optimization. In Confer- ence on Uncertainty in Artificial Intelligence, 2016. 5. [Friedman et al., 2008] Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estima- tion with the graphical lasso. Biostatistics, 9(3):432-441, 2008. 6. [Gao et al., 2015] Hongchang Gao, Lin Yan, Weidong Cai, and Heng Huang. Anatomical annotations for drosophila gene expression patterns via multi-dimensional visual de- scriptors integration: Multi-dimensional feature learning. In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 339-348. ACM, 2015. 7. [Gao et al., 2017] Hongchang Gao, Feiping Nie, and Heng Huang. Local centroids structured non-negative matrix factorization. In Thirty-First AAAI Conference on Artifi- cial Intelligence, 2017. 8. [Gower et al., 2016] Robert Gower, Donald Goldfarb, and Peter Richtárik. Stochastic block bfgs: squeezing more curvature out of data. In International Conference on Ma- chine Learning, pages 1869-1878, 2016. 9. [Johnson and Zhang, 2013] Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. In Advances in Neural Information Processing Systems, pages 315-323, 2013. 10. [Lee et al., 2006] Su-In Lee, Honglak Lee, Pieter Abbeel, and Andrew Y Ng. Efficient l˜1 regularized logistic re- gression. In AAAI, 2006. 11. [Li et al., 2016] Xingguo Li, Raman Arora, Han Liu, Jarvis Haupt, and Tuo Zhao. Nonconvex sparse learning via stochastic optimization with progressive variance reduc- tion. arXiv preprint arXiv:1605.02711, 2016. 12. [Mallat and Zhang, 1993] Stéphane G Mallat and Zhifeng Zhang. Matching pursuits with time-frequency dic- tionaries. IEEE Transactions on signal processing, 41(12):3397-3415, 1993. 13. [Moritz et al., 2016] Philipp Moritz, Robert Nishihara, and Michael Jordan. A linearly-convergent stochastic l-bfgs algorithm. In Artificial Intelligence and Statistics, pages 249-258, 2016. 14. [Needell and Tropp, 2009] Deanna Needell and Joel A Tropp. Cosamp: Iterative signal recovery from incomplete and inaccurate samples. Applied and Computational Har- monic Analysis, 26(3):301-321, 2009. 15. [Nguyen et al., 2014] Nam Nguyen, Deanna Needell, and Tina Woolf. Linear convergence of stochastic iterative greedy algorithms with sparse constraints. arXiv preprint arXiv:1407.0088, 2014. 16. [Nocedal and Wright, 2006] Jorge Nocedal and Stephen J Wright. Numerical optimization (2nd edition), 2006. 17. [Tibshirani, 1996] Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B (Methodological), pages 267-288, 1996. 18. [Tropp and Gilbert, 2007] Joel A Tropp and Anna C Gilbert. Signal recovery from random measurements via orthogo- nal matching pursuit. IEEE Transactions on information theory, 53(12):4655-4666, 2007. 19. [Van de Geer, 2008] Sara A Van de Geer. High-dimensional generalized linear models and the lasso. The Annals of Statistics, pages 614-645, 2008. 20. [Yuan and Lin, 2006] Ming Yuan and Yi Lin. Model selec- tion and estimation in regression with grouped variables. Journal of the Royal Statistical Society: Series B (Statisti- cal Methodology), 68(1):49-67, 2006. 21. [Yuan and Liu, 2014] Xiao-Tong Yuan and Qingshan Liu. Newton greedy pursuit: A quadratic approximation method for sparsity-constrained optimization. In Proceed- ings of the IEEE Conference on Computer Vision and Pat- tern Recognition, pages 4122-4129, 2014. 22. [Yuan et al., 2014] Xiaotong Yuan, Ping Li, and Tong Zhang. Gradient hard thresholding pursuit for sparsity- constrained optimization. In Proceedings of the 31st In- ternational Conference on Machine Learning (ICML-14), pages 127-135, 2014. 23. [Zhang, 2011] Tong Zhang. Adaptive forward-backward greedy algorithm for learning sparse representations. IEEE transactions on information theory, 57(7):4689-4708, 2011. 24. [Zhao et al., 2017] Renbo Zhao, William B Haskell, and Vincent YF Tan. Stochastic l-bfgs revisited: Improved convergence rates and practical acceleration strategies. arXiv preprint arXiv:1704.00116, 2017.