Stochastic Second-Order Method for Large-Scale Nonconvex Sparse Learning Models

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.

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.

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