{'fraction_non_alphanumeric': 0.03942739562360334, 'fraction_numerical': 0.013076089120959356, 'mean_word_length': 4.2053437751929055, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Safe Active Feature Selection for Sparse Learning

We present safe active incremental feature selection~(SAIF) to scale up the computation of LASSO solutions. SAIF does not require a solution from a heavier penalty parameter as in sequential screening or updating the full model for each iteration as in dynamic screening. Different from these existing screening methods, SAIF starts from a small number of features and incrementally recruits active features and updates the significantly reduced model. Hence, it is much more computationally efficient and scalable with the number of features. More critically, SAIF has the safe guarantee as it has the convergence guarantee to the optimal solution to the original full LASSO problem. Such an incremental procedure and theoretical convergence guarantee can be extended to fused LASSO problems. Compared with state-of-the-art screening methods as well as working set and homotopy methods, which may not always guarantee the optimal solution, SAIF can achieve superior or comparable efficiency and high scalability with the safe guarantee when facing extremely high dimensional data sets. Experiments with both synthetic and real-world data sets show that SAIF can be up to 50 times faster than dynamic screening, and hundreds of times faster than computing LASSO or fused LASSO solutions without screening.


Introduction
LASSO has been a powerful tool for sparse learning to analyze data sets with p ≫ n, where p is the number of covariates or features and n the number of samples. LASSO screening methods provide efficient approaches to scale up sparse learning without solving the full LASSO problems, based on either sequential or dynamic screening methods (Ghaoui et al., 2012;Wang et al., 2014a;Ndiaye et al., 2015;Ren et al., 2017). However, the existing sequential screening requires the LASSO solution with a heavier regularization penalty parameter so that the range of dual variables can be estimated tightly to help effectively screen redundant features. Different from such static sequential screening methods, dynamic screening does not require the solution with the heavier penalty parameter but relies on duality gaps during optimization iterations for feature screening. To achieve high screening power, a significant number of optimization iterations have to be operated on the full-scale problems with the original high-dimensional feature set to compute the effective duality gap. Both sequential and dynamic screening require to update the original full-scale LASSO model. To further scale up the existing screening methods, we introduce an active incremental screening method-safe active incremental f eature selection (SAIF)-that can overcome all these shortcomings. SAIF dynamically manipulates the feature set according to the duality gap of the sub-problems with only a small number of active features in hand. In this way we can maximally reduce the redundant computation on inactive features that have zero coefficients in the optimal LASSO solution.
We first review the relevant literature on screening and homotopy methods for sparse learning, and then present the theoretical and experimental results for SAIF in the following sections.

Sequential Screening
Most sequential screening methods derive screening rules by leveraging the solutions to the LASSO model with a heavier regularization parameter. There are two broad categories of sequential screening methods for LASSO problems: heuristic and safe screening methods. The heuristic screening methods (Tibshirani et al., 2012;Fan and Lv, 2008) rely on heuristics to remove features. For example, the strong rule (Tibshirani et al., 2012) derives the screening rule based on the assumption that the absolute values of the inner products between features and the residue are non-expansive with respect to the parameter values. It is obvious that this assumption does not always hold. Such heuristic screening rules are not safe, meaning that they cannot guarantee that the removed features will have corresponding zero coefficients in the optimal LASSO solution to the original fullscale problem. Sequential screening methods, such as the ones proposed by Ghaoui et al. (2012); Wang et al. (2014a,b); Ren et al. (2017), do not take the unsafe assumptions that the heuristic screening methods use, but try to develop safe feature screening rules based on the structure of the problem. Most of these screening methods are inspired by the seminal work by Ghaoui et al. (2012) and derive screening rules with the help of the LASSO solution with a heavier regularization parameter.
Assume n data samples with p features form an n × p design matrix X, with an n × 1 label vector, y. Let the general loss function f in (1) be α-smooth, and γ-convex with respect to (w.r.t.) the L 2 norm (α > 0, γ > 0), and f * is the conjugate of f . According to Theorem 6 in Kakade et al. (2009), f * is 1 α -convex 1 γ -smooth w.r.t. the L 2 norm. The primal and dual forms of the general LASSO problem can be written as follows: where x i is the ith column of X, i.e, the ith feature vector; the row vector x j• denotes the jth sample; F denotes the index space of the features in the original LASSO problem; θ is the dual variable. The derivation of the dual problem (2) can be referred to Ndiaye et al. (2015); Wang et al. (2014a). The optimal primal and dual variable relationship is f ′ (x j• β * ) = −λθ * j , where f ′ is the first-order derivative of f . With KKT (Karush-Kuhn-Tucker) conditions (Tibshirani and Taylor, 2011;Wang et al., 2014a), we have According to (4), if |x T i θ * | < 1, we can have [β * ] i = 0, and i is an inactive feature. Most existing screening methods aim to estimate a convex or ball region B(θ, r) = {θ * | ||θ * − θ|| 2 ≤ r} as the range of θ * . If θ * ∈ B(θ, r), we can write θ * = θ + ρ with ||ρ|| 2 ≤ r. It follows that x T i θ * = x T i θ + x T i ρ, and we have x T i θ − ||x i || 2 r ≤ x T i θ * ≤ x T i θ + ||x i || 2 r. Thus |x T i θ| + ||x i || 2 r < 1 =⇒ i is an inactive feature.
Safe sequential screening methods estimate the ball region of the dual variables θ * at λ as B(θ * (λ ′ ), r). Here λ ′ > λ, and θ * (λ ′ ) can be computed based on the primal-dual relation for the solution to the LASSO problem with λ ′ as the regularization penalty parameter. For example, DPP (Dual Polytope Projection) (Wang et al., 2014a) takes the dual problem (2) as a projection problem and relies on the estimation of B(θ * (λ ′ ), r) for the range of θ * based on the non-expansiveness properties. Typically, sequential screening requires to solve a sequence of LASSO problems corresponding to a sequence of descending λ's to tighten the range estimates of θ * to achieve the high screening power. Such a sequential procedure is suitable and efficient when solving a sequence of sparse learning problems with different regularization parameters, for example, for parameter selection by cross validation. However, sequential screening methods are not absolutely safe since there is always computational error to the solution of LASSO with a heavier penalty parameter, which is pointed out by the authors of dynamic screening methods (Ndiaye et al., 2015) reviewed in the next subsection.

Dynamic Screening
Different from sequential screening methods that require the solution with a heavier penalty parameter, dynamic screening (Ndiaye et al., 2015;Bonnefoy et al., 2015) can scale up LASSO solutions by dual variable range estimation based on the duality gap during the algorithm iterations. The ball region for θ * is estimated based on the duality gap as a function of the primal and dual objective function values at iterative updates (Ndiaye et al., 2015;: Here Ω F = {θ | |x T i θ| ≤ 1, ∀i ∈ F} is the dual feasible space corresponding to the feature set F; β is the current estimation of primal variables; and θ is the projected feasible dual variables of β. P (β) and D(θ) are the primal and dual objective function values at β and θ, respectively. The tightness of the results depends on the duality gap [P (β) − D(θ)], determined by the quality of iterative updates for β and θ. Dynamic screening algorithms in Ndiaye et al. (2015);  iteratively update β and θ for the original LASSO problem with the whole feature set X to check the duality gap and apply screening rules to remove inactive features. Without the solution information from a heavier parameter, dynamic screening has to iterate the operations in optimization, such as sub-gradient computation, on the original whole feature set many times to gain a small duality gap. The computation cost of these operations dilute the screening benefits as the iterations have to be repeated many times to arrive at a duality gap that is small enough to achieve desired screening power.

Homotopy Methods
Homotopy methods have been applied for LASSO to compute the solution path when λ varies (Efron et al., 2004;Osborne et al., 2000;Malioutov et al., 2005;Garrigues and Ghaoui, 2008;Friedman et al., 2010;. This type of methods rely on a sequence of decreasing λ values and "warm start" (starting the active set with the solution from the previous λ) to achieve computational efficiency. Usually these methods have multiple iteration loops to incorporate the strong rule screening, active set, and path-wise coordinate descent. The inner loop performs coordinate descent and active set management. The outer-loop goes through a sequence of decreasing λ values and initializes the active set at each λ with the strong rule and warm start. Since they do not utilize safe convergence stopping criteria for the active set, they may miss some of the active features in the optimal solutions to the original LASSO formulation with the corresponding λ values. Furthermore, this type of methods do not employ any screening rule for the inner-loop sub-problems, and it may limit their scalability.
Besides screening and homotopy methods, working set methods (Johnson and Guestrin, 2015) maintain a working set according to some violation rules and solve a sub-problem regarding the working set at each step. By estimating an extreme feasible point based on the current solution, this method constructs the working set for the next step by the

At Rt
Sub-Prob. Iterations ADD DEL Figure 1: SAIF feature selection. A t stands for the Active set, while R t stands for the Remaining set at step t.
constraints that are closest to the feasible point. This kind of methods also start from solving the original full-scale problem as the existing LASSO screening methods. However, when p ≫ n, the basic assumption of sparse learning is that most of the given features are irrelevant and should be inactive for the optimal solutions. These existing algorithms may not be efficient due to redundant time-consuming operations on inactive features.

Our Contributions
In this paper, we propose a novel safe LASSO feature selection method to further scale up LASSO solutions by overcoming the issues in the existing methods. Instead of taking the whole feature set as the initial input, our method SAIF starts from a small set of features, which is taken as the active set. The features that are not in the active set are put in the remaining set ( Figure 1). Time-consuming iterations such as coordinate minimization with soft-thresholding are only performed on the features in the active set. Features are actively recruited or removed from the active set according to the estimated ranges of optimal dual variables. Based on duality properties, efficient feature operation rules and safe stopping criteria have been developed to keep most inactive and redundant features out of the active set. With a small active set, CPU time and memory operations can be tremendously reduced. Complexity analysis is provided for both dynamic screening and SAIF for comparison. Theoretical results show that SAIF converges to the optimal solution to the original LASSO problem and the running time of SAIF is only proportional to the active feature size, the number of features with non-zero coefficients in the optimal LASSO solution, rather than the original input feature size. Experiments on simulated and realworld data sets verified the advantages of the proposed method. Moreover we show that the proposed method can be extended to fused LASSO with tree dependency structures among features.
Data: Data matrix X, label y, penalty λ, stopping duality gap ǫ Result: Coefficient Vector β Choose ⌈c log( md+mx λ ) log(p)⌉ features from F in the descending order of |X T f ′ (0)|; δ = λ λmax , IsAdd = True; while True do Update β t with K iterations of soft-thresholding operations on A t ; Compute a ball region B(θ t , r t ) based on (11) or (12); IsAdd = False; Continue; end end ADD operation; end end Put β t in to β, and inflate the other entries with 0.

Safe Active Incremental Feature (SAIF) Selection for LASSO
We derive an innovative incremental feature screening algorithm, SAIF, in which we can iteratively solve much smaller sub-problems than the original LASSO problem, i.e., iteratively update the duality gap while adding or removing features by leveraging the active ball region estimates for the optimal dual variables of these sub-problems. The schematic illustration of SAIF is given in Figure 1. Let A t and R t denote the active feature index set and remaining feature index set at iteration step t, respectively. Instead of solving either the original full-scale LASSO primal problem or the corresponding dual problem, SAIF screening is different from the existing sequential and dynamic screening as it only needs to solve significantly reduced sub-problems and updates the screening rules based on the duality gap without solving these sub-problems exactly. More importantly, SAIF has the safe guarantee that only irrelevant or redundant features in the original LASSO problem will be removed. Algorithm 1 summarizes our SAIF screening procedure, which starts with A 0 and dynamically moves active features between R t and A t .

ADD and DEL Operations
Two operations in SAIF are ADD and DEL. Starting from an initial active set A 0 , whose features can be selected by some simple heuristics, for example, based on their correlation with the output, SAIF iteratively adds features (ADD) into or removes features (DEL) from the active set. At the tth iteration, we derive both ADD and DEL operations to dynamically update A t based on the primal sub-problem with only the current active features: Let Ω At be the dual feasible region and D(θ t ) denote the dual objective function value of the sub-problem at the dual variable θ t considering only the active features in A t with θ * t being the corresponding optimal dual solution. We use β t ∈ R |At| to represent the updated β values after t out-loop iterations in SAIF. β * t denotes the optimal active feature solution to the problem P t . P t (β) is the objective value of P t with inputβ, andβ can have a different set of features from A t ; we inflate the missing entries inβ with zeros and ignore the entries or features not in A t in the calculation of P t (β). Let S A represent the set of the optimal primal solutions for any feature set A, θ * the optimal dual solution with the full feature set F, andĀ for the optimal active feature set that {i : |x T i θ * | = 1}. Let B(θ t , r t ) = {θ * t ||θ * t − θ t || 2 ≤ r t } be an estimated ball region for θ * t at step t. SAIF carries out ADD and DEL operations as follows: We have the following theorem regarding ADD and DEL operations: Proof: a) From the dual form (8), if we add i to A t , there will be one more constraint for the dual problem at step t + 1, thus Ω A t+1 ⊆ Ω At . As we have a smaller feasible space at As Ω At is convex and closed, and f * is convex and smooth, the optimal dual solution for the active set A t is unique, which means D( Since the dual problem is convex and smooth, and the feasible set is closed and convex, θ * t = θ * . Hence, β * t ∈ S F as the primal solution may not be unique. Remark 1 Theorem 1-c) provides us with the stopping criterion for ADD operations in our SAIF algorithm. We can apply ADD and DEL operations in iterations to minimize This concludes that our stopping criterion for ADD operations ensures the safe feature screening.
The DEL operation is similar to the screening steps in dynamic screening. As we can see, at step t with the DEL operation, D(θ * t ) = D(θ * t−1 ). Theorem 1-a) implies D(θ * t ) ≤ D(θ * t−1 ). Thus the optimal dual objective value always goes down. Theorem 1-c) and Remark 1 suggest that after the stopping of ADD operations, A t already has all the active features for the original problem. The algorithm then stops once it reaches the pre-specified accuracy value of the duality gap. Such monotonicity leads to the convergence of SAIF detailed in Section 3.

Implementation
We first discuss how we derive a tighter ball region B(θ t , r t ) for the range estimate of θ * t , taking the advantages of the existing screening methods.
Dual variable range estimation: Accurately estimating the range of θ * t , B(θ t , r t ), for the sub-problem is critical for efficient SAIF screening with ADD and DEL operations at each iteration. With f as the vector form of the loss function with all of the samples, we provide the following theorem to estimate the ball region for θ * t with the similar idea from sequential screening.
Theorem 2 For the LASSO problem (1) with the loss function f , if f * is 1 α -strongly convex, and θ * 0 and θ * are the optimal solutions to the dual problem (2) at λ 0 and λ with λ < λ 0 , then If we have θ ∈ Ω F , the bound can be further improved by ). The proof for Theorem 2 can be found in Appendix A. At step t with the active set , we can use Theorem 2 to estimate θ * t . For linear regression, the estimation can be further improved based on the projection properties as in DPP (Wang et al., 2014a).
Theorem 2 provides a tight estimation when λ 0 is close to λ. When λ is far away from λ 0 , we can adopt the tighter dual variable range estimation with the following ball region by dynamic screening (Ndiaye et al., 2015;. At step t, we have For β t , with the primal-dual relation, we can project it to the dual feasible region Ω At to get a feasible dual variable θ t . With two ball regions from Theorem 2 and the duality gap, we can derive a tighter constrained region by computing the center and radius of a ball region B(θ t , r t ) that covers the intersection of two ball regions, B 1 (θ 1 , r 1 ) and B 2 (θ 2 , r 2 ): where B 1 can be derived from Theorem 2, and B 2 from (11). The resulting B(θ t , r t ) gives us a tighter bounded region at step t when r t < min{r1, r2}. When we do not have the solutions with other λ values, we simply set the bounded region for θ * t based on (11). Improve SAIF with an estimation factor: The estimation of dual variables may be inaccurate to have high enough screening power during the optimization iterations, especially at the beginning of the algorithm. We add a factor to the radius of the ball region to reduce redundant computation resulted from inaccurately recruited features. At the beginning of Algorithm 1, δ is a value smaller than 1. δ will be increased to 1 during the SAIF iterations to ensure the safe guarantee of SAIF algorithm.
ADD operation implementation details: The number of added features in each ADD operation can vary to reduce redundant iterations. Generally, the relationship between the screening power and this number depends on the regularization parameter λ and how well feature vectors x i , i ∈ F, correlate with the outcome label y. In this paper, we empirically set the number to be h = ⌈c log( md+mx λ ) log(p)⌉. Here mx and md are the maximum and median of |X T f ′ (0)| (|X T y| with linear regression). Many iterations may need to be operated to reach the dual space point that can distinguish h features, and this may reduce the efficiency of the algorithm. We can decrease the redundancy by relaxing the strict condition in Theorem 1-d). Let V i represent the set of features that violate the condition in Theorem 1-d) regarding feature i, i.e., Hereh = ⌈ζh⌉, and ζ > 0. Algorithm 2 summarizes the implementation of the ADD operation.
Algorithm 2: Algorithm for ADD operation

Convergence Analysis
In this section, we first discuss the convergence properties of SAIF and then provide the detailed complexity analysis of our SAIF algorithm. In this manuscript we present the complexity of SAIF with coordinate minimization (CM) as the inner base algorithm. We can derive the corresponding complexity in a similar way if an alternative base algorithm such as FISTA (Beck and Teboulle, 2009) is employed.

Algorithm Properties
Similar to dynamic screening, SAIF employs coordinate minimization (CM) in the primal variable space. Besides feature screening (by DEL operations), SAIF has feature recruiting operation (ADD). In this subsection, we first discuss the convergence of the base algorithm, then we show that the numbers of DEL and ADD operations are finite in SAIF.

Coordinate Minimization (CM)
The base algorithm we employ in SAIF is shooting algorithm (Fu, 1998), which is a cyclic block coordinate minimization method. Coordinate descent (CD) and coordinate minimization (CM) methods have been studied extensively by many researchers (Nesterov, 2012;Beck and Tetruashvili, 2013;Wright, 2015). Recently Li et al. (2017) have achieved faster convergence rates for CD and CM methods on convex problems. Based on the analysis from Li et al. (2017), we can prove the following lemma regarding CM for LASSO. 1 Lemma 1 (Adaptation of Li et al. (2017)) For the LASSO problem (1) with a γ-convex loss function, at most log ψ ε P (β 0 )−P (β * ) base operations are performed with cyclic coordinate minimization to arrive at β a such that P σ max is the largest eigenvalue of X T X, L is the Lipschitz constant of f ′ , and β 0 denotes the starting point of the primal variables.
The base operation (soft-thresholding) in CM is operated in the primal variable space. Feature screening or feature selection operations such as ADD and DEL operations rely on the dual variable estimation. We provide the following lemma to show that the accuracy of dual variable estimates is almost linearly bounded by the accuracy of primal variable estimates when the iteration number is large.
Lemma 2 For the primal problem (1) and dual problem (2), , and θ k = τ kθk , with a large enough k in coordinate minimization, we have The proofs of Lemmas 1 and 2 can be found in Appendix B. With Lemma 2, we can see that the estimation of dual variables relies on the accuracy of primal variable estimates. In SAIF, the starting point for each β t is already with relatively high accuracy as empirically there are only one or a few features different between steps t and t − 1.

Finite number of ADD and DEL Operations
With CM as the inner base algorithm, we prove that the outer loop can stop in a finite number of steps. The ADD operation recruits more features into the active set, and thus results in decreasing optimal objective value as shown in Theorem 1. Since the DEL operation does not change the optimal objective function value, the corresponding optimal dual objective function value of the sub-problem decreases monotonically and finally converges to the value of the original full-scale problem. Experimentally, for a given λ, the running time of SAIF is proportional to the size of the optimal active setĀ. The following theorem provides the guarantee for the convergence of SAIF.
Proof : a) IfĀ A t , from Remark 2 in Section 2.1, we can see that, ∃i ∈ R t , |x T i θ * t | ≥ 1. If max i∈Rt , |x T i θ * t | > 1, we can apply the ADD operation at step t to add the most active feature to A t+1 . We will have D( If max i∈Rt , |x T i θ * t | = 1, the optimal dual variable is already on the hyperplanes |x T i θ * t | = 1. From the SAIF algorithm, we can see that, with an ADD operation to move all x i : |x T i θ * t | = 1 to A t , the optimal dual solution will remain the same, i.e., θ * t+1 = θ * t . The ADD operation will stop at step t, as max i∈R t+1 , |x T i θ * t+1 | < 1. DEL does not remove as the optimal dual variable will remain the same, and the algorithm will stop.

Complexity Analysis
For complexity analysis, we split the SAIF algorithm into three phases: feature recruiting, inactive feature deletion, and accuracy pursuing. The inactive feature deletion phase is the same as the feature screening phase in dynamic screening. We first present the complexity analysis for dynamic screening, which is our additional contribution, and then based on that and previous results, we give the detailed complexity analysis for SAIF.

Complexity Analysis for Dynamic Screening
Dynamic screening (Ndiaye et al., 2015; starts its active set with the whole feature set. Let r i be the radius of the ball region for the screening of feature i, according to the same screening rule as the DEL operation, Here m is the feature with the value of max i |x T iθ t |, x m is the corresponding feature vector, . If feature i does not belong to the final active set A, then |x T i θ * | < 1. With large t, m belongs toĀ according to Theorem 1, and |x T m θ * | = 1. We have Thus the screening radius for feature i is determined by how closeθ t and θ * are, linearly determined by the primal variable accuracy according to Lemma 2. With ε as the prespecified objective function value accuracy, the following theorem gives the time complexity of the dynamic screening procedure. Theorem 4 Assuming that the time complexity for one coordinate minimization operation is O(u), the time complexity for dynamic screening is O uL Here G 0 = P (β 0 ) − P (β * ), and ε D is the accuracy of the objective function value for the last feature screening operation.
Remark 3 With coordinate minimization, the number of iterations to reach the accuracy of the objective function value ǫ is O L2 γ 2 p log 1 ε D +|Ā| log ε D ε . As p >> |Ā|, the computation cost in dynamic screening is mainly from the iterations to reach ǫ D .
The proof of Theorem 4 can be found in Appendix C. Experiments will confirm the conclusions from Theorem 4 and Remark 3 with the results presented in Section 5.

Complexity Analysis for SAIF
With the complexity analysis for dynamic screening, we now derive the complexity of SAIF and show its advantages over dynamic screening theoretically. SAIF starts the algorithm from the feature recruiting phase. The ADD operation recruit a feature with max i∈Rt |x T i θ t |.
Here we use θ * t rather than θ * as the algorithm has not reached the stopping point of ADD operations andĀ A t . In (16), the radius for adding feature i into the active set is determined by how much it can outperform the other features. We use T a to represent the running time consumed in the feature recruiting phase. The inactive feature deletion phase starts from setting IsADD = False in SAIF in Algorithm 1. Let p A be the total number of features involved in the ADD operation; after the dth feature (d in the sequence of {1, 2, ..., d, ..., p A }) has been added into the active set, we use P d , p d , and β * d to denote the primal objective function, the size of the active set, and the optimal primal solution of the sub-problem, respectively. Let Q d (β) = P d (β) − P d (β * d ). The time complexity for SAIF with CM is given by the following lemma and theorem.
Lemma 3 With O(u) as the complexity for the base operation of coordinate minimization of the LASSO problem with a γ-convex loss function, p A the total number of features involved in ADD operations, and p p A the size of the active set when IsADD is set to false, the complexity for the feature recruiting phase in SAIF is Theorem 5 With O(u) as the complexity for the base operation of coordinate minimization of the LASSO problem with a γ-convex loss function, the time complexity for SAIF is Herep is the maximum size of the active set during the algorithm iterations,Q is the geometric mean of the accuracies of the sub-problem objective function values corresponding to each ADD operation, and ε D is the accuracy of the objective function value for the last feature DEL operation.

Remark 4 With coordinate minimization, the number of iterations to reach the accuracy of the objective function value
.Q is a value much smaller than G 0 in dynamic screening (as the value of Q d for adding feature d usually is very small).
The proofs for Lemma 3 and Theorem 5 are given in Appendix C. According to our experiments,p is often close to the number of the actual active features in the optimal LASSO solution, |Ā|. The dominating factor for the computational complexity of SAIF is the second termpp A . The fewer features being added in the active set, the less time SAIF will consume. Experimentally, p A is often a value several times larger than |Ā|, and p A << p. We can conclude that SAIF takes much less time than dynamic screening based on the analysis of Theorems 4 and 5. With the theoretical safe and convergence guarantees, SAIF can work with extremely high-dimensional data to obtain optimal LASSO solutions.

SAIF for Fused LASSO
The active incremental philosophy can be applied to many convex sparse models. In this section, we show how to scale up fused LASSO solutions based on SAIF. The formulation for fused LASSO is where ||Dβ|| 1 = (a,b)∈E |β a − β b |, and each pair in E denotes an edge in a tree with F as the vertex set. The tree G(F, E) captures the dependency structures among features.
Here D is a matrix representation of the tree, and in each row of D, we have zero entries except two with values equal to 1 and −1, respectively. With M −p denotes the reduced matrix of a given matrix (vector) M without its corresponding pth column (entry), the fused LASSO problem can be further transformed into the equivalent LASSO formulation with the following theorem.
HereX =X −p , and H = I h , h = −x 1,p xn,p , ..., −x n−1,p xn,p .θ = Hθ −p , and θ = − , andx i is the ith column ofX. With the primal form (18) and dual form (19) in Theorem 6, we just need a transformation on the feature set to apply SAIF to fused LASSO problems. From the proof of Theorem 2 in Ndiaye et al. (2015), we can easily get ∀θ ∈ Ω,β = β b ∈ R p×1 , ||θ * −θ|| 2 2 ≤ 2 λ 2 P (β) −D(θ) . With the duality gap, we can derive the ADD and DEL operation for fused LASSO. The following Theorem shows how to project the current dual estimationθ in (19) to the feasible space Ω for regression with the least square loss function.
Theorem 7 For linear regression problems with fused LASSO regularization,θ = τθ is the scaled feasible dual variable vector ofθ closest toθ * in (19) with τ = min max{ y Tθ ||X Tθ ||∞ . The proofs for Theorems 6 and 7 are presented in Appendix D. The algorithm for fused LASSO is the same as what we presented for SAIF for LASSO with the transformation steps. As the transformation matrix is highly sparse and only has column operations on the feature matrix X, we can replace matrix multiplication with column operations to further improve computation efficiency.

Experiments
In this section, we present the experiments comparing SAIF with other existing LASSO and fused LASSO methods. We first evaluate the selected methods for the LASSO formulation based on a simulation study and then apply them to one real-world study. In the second subsection, we evaluate SAIF for logistic regression with two real-world data sets. We present the comparison between SAIF and sequential screening as well as homotopy methods in the third subsection. In the last subsection, two studies are performed to evaluate the selected methods for fused LASSO. The base algorithm (coordinate minimization) is implemented with C, and the main algorithms of SAIF, dynamic screening (Ndiaye et al., 2015), DPP (Wang et al., 2014a), and the homotopy method  are coded with Matlab. We use the provided package (Johnson and Guestrin, 2015) for the BLITZ method. The experiments are run on an iMac with OS Sierra version 10.12.1 and Intel Core i5. The implementation environment is the same for all the experiments unless specified.

Results for Linear Regression
Similar to sequential and dynamic screening algorithms, SAIF can be assembled with different kinds of LASSO solution methods. Shooting algorithm (CM) is chosen as the base algorithm in our experiments. Both dynamic screening (Ndiaye et al., 2015) and SAIF can do feature screening or selection without the help from a heavier parameter solution. We specifically focus on the performance comparison among (1) shooting algorithm without screening (No Scr.), (2) shooting algorithm with dynamic screening (Ndiaye et al., 2015) (Dyn. Scr), (3) working set method (Johnson and Guestrin, 2015) (BLITZ), and (4) shooting algorithm with SAIF screening (SAIF). All of these are safe methods for LASSO problems.

Simulation Study
First, we simulate the data sets with n = 100 samples and p = 5, 000 features according to a linear model y = Xβ + ǫ, where each column of X is a vector with random values uniformly sampled from the interval [−10, 10], and the white noise ǫ ∼ N (0, 1). For the linear coefficients β, 20% entries (0.2p) are randomly set to the values in [1, −1], and the rest (0.8p) to zero. For this data set, we can derive λ max = 2.183E4. The left plot in Figure 2 illustrates the running time for different methods in the logarithmic time scale with λ = 20, 100, and 1, 000 at different stopping accuracy 1.0E − 9 and 1.0E − 6 for the desired duality gap (DGap). We can see that, SAIF takes much less time than the other methods to reach the optimal solutions with the specified accuracy. The results also show that SAIF is more efficient compared to the existing safe methods when λ is small.

Breast Cancer Data
Breast cancer data set consists of gene expression data of 8,141 genes for 78 metastatic and 217 non-metastatic breast cancer patients from the study introduced in Chuang et al. (2007). In this set of experiments, the metastatic samples are labeled as 1 and nonmetastatic as -1 as the output of the LASSO linear regression problem. The right plot in Figure 2 compares the running time for the same four different methods at different λ's. Again, SAIF takes the least computation time at different duality gaps.
We further investigate the size of the active set along with the optimization iterations for dynamic screening and SAIF in Figures 3-a,c), with λ = 0.1 and 5. We can see that SAIF starts from a small active feature set and gradually increases its size with time, while λ λmax (x-axis) and log(100×t(sec.)) (y-axis) for dynamic screening (a), and SAIF (b) on breast cancer data. dynamic screening starts from the whole feature set and takes longer time to reach the point with screening power. Figures 3-c,d) illustrate the change of the dual objective function values D(θ t ) for SAIF during the optimization procedure. With the active feature set size increasing, D(θ t ) decreases and finally converges to a steady value D(θ * ), indicating the algorithm obtains the optimal solutions to the original LASSO problems.
Let p t be the feature number at step t for SAIF or dynamic screening. The left column in Figure 4 shows the change of pt p with respect to the regularization penalty (log 10 λ λmax on x-axis) and the optimization time (log(100 × t(sec.)) on y-axis). Similarly, we plot the change of log( pt p ′ ), where p ′ is the corresponding optimal active feature size in the right column of Figure 4. From Figure 4, it is clear that dynamic screening always takes longer time to reach the optimal active feature set size, especially when λ is small. Before reaching the point with screening power, the active feature set size is almost p. While the active feature set size for SAIF grows gradually from a small set. Due to the small active set size for the starting iterations, SAIF can more efficiently reach the optimal active set size with much shorter running time. All of these results confirms the theoretical complexity analysis for dynamic screening and SAIF. Furthermore, both Figures 3 and 4 illustrate that SAIF is more scalable than the existing methods as it always starts from a very small active set and iteratively focuses on a small subset of the features.

Results for Logistic Regression
We evaluate the proposed algorithms for sparse logistic regression with two data sets, Gisette and USPS, from LibSVM (Chang and Lin, 2011) Website. The Gisette data set has 5,000 features and 6,000 samples. There are 256 features, 7,291 samples, and 10 labels in the USPS data set; and we categorize the label values large than 4 as positive, and negative otherwise. For these two data sets, λ max is 932,575 and 992, respectively. We use the logistic regression code from L1General (Schmidt, 2006) as the inner iteration algorithm in this set of experiments. Figure 5 plots the running time at different λ values for dynamic screening, BLITZ, and SAIF. Although BLITZ may achieve comparable performance in a couple of cases when the active set is very small, SAIF consistently takes less computation at different λ values for both data sets. From all these results, SAIF can achieve more computational efficiency for both linear and logistic regression compared to the existing safe screening methods.

Comparison with Sequential Screening and Homotopy Methods
With a sequence of decreasing λ values, SAIF can be further improved with the warm start strategy. Given the simulation and the breast cancer data sets in Section 5.1, a decreasing sequence of λ values are evenly sampled from the logarithmic scale of the range [0.001λ max , λ max ]. The plots in Figure 6 present the running time for DPP (Wang et al., 2014a), the homotopy method , and SAIF with a different number of λ values on both data sets. In this set of experiments, we set the stopping criterion with the  . According to the recall values, the homotopy method always misses some of the active features at different numbers of λ values. Furthermore, the homotopy method leads to the inclusion of inactive features into the final solution as evidenced in Table 1 that the precision cannot reach 1 at different numbers of λ values. On the contrary, our SAIF has the safe guarantee: the recall and precision values regarding active features recovered by SAIF are always one. Clearly, the unsafe strategies employed by the homotopy method do not always reduce computation, and the employed inactive features may lead to larger CPU time consumption as shown in the left plot in Figure 6.

Results for Fused LASSO
We further present the experiments for fused LASSO with the formulation (17). There are a few solvers that are suitable for tree fused LASSO problems, such as CVX and the path solution method (Arnold and Tibshirani, 2016). Due to the scalability and solution accuracy issues with the path solution package, we only take CVX as the baseline for comparison in our experiments. We first compare the running time between SAIF and CVX on breast cancer data regarding fused LASSO linear regression; then we compare them on the FDG-PET data set (ADNI) by logistic regression.

Breast Cancer Data
For the same breast cancer data set as in the previous subsections, we would like to incorporate the interaction relationships among genes to formulate the fused LASSO problems for regression analysis. The largest connected component in the human protein-protein interaction (PPI) network was identified in Chuang et al. (2007) to capture the gene-gene relationships by a connected graph with 7,782 nodes. The left plot in Figure 7 gives the running time for both CVX and SAIF at different λ's for the specified duality gap 1.0E − 6.
The results show that SAIF can significantly reduce computation cost compared to the CVX solver without screening for fused LASSO.

FDG-PET Data Set
The FDG-PET data set has 74 Alzheimer's disease (AD) patients, 172 mild cognitive impairment (MCI) subjects, and 81 normal control (NC) subjects, downloaded from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database. 116 features (each feature corresponds to a brain region) can be derived for each subject after preprocessing. We further use the method described in Yang et al. (2012) to construct a correlation tree on these features. We take AD as positive (1) and NC as negative (0), and all of MCI samples are not used in this set of experiments in fitting to a fused LASSO logistic regression model. The right plot in Figure 7 gives the running time for three λ values at duality gap 1.0E − 6. Again SAIF takes much less time than CVX does on this data set.

Conclusions
In this paper, we have developed a novel feature selection method for LASSO-SAIF. From the experimental results, SAIF can achieve improved efficiency compared with existing methods. SAIF has the potential to scale up for data sets with high dimensional features due to its incremental property. Theoretical analysis reveals the safety guarantee and low algorithm complexity of the proposed method. Furthermore, SAIF can be potentially extended to group LASSO (Yuan and Lin, 2006) and other sparse models. SAIF provides us with a new direction for scaling up sparse learning. Given a data set with extremely high feature dimension, SAIF can be further improved with the multi-level active set and remaining set schema.

Appendix A. Proof of Theorem 2
Theorem 2 For the LASSO problem (1) with the loss function f , if f * is 1 α -strongly convex, and θ * 0 and θ * are the optimal solutions to the dual problem (2) at λ 0 and λ with λ < λ 0 , then If we have θ ∈ Ω, the bound can be further improved by .

Appendix B. Proof of Lemma 1 and Lemma 2
Lemma 1 (Adaptation of Li et al. (2017)) For the LASSO problem (1) with a γ-convex loss function, at most log ψ ε P (β 0 )−P (β * ) base operations are performed with cyclic coordinate minimization to arrive at β a such that P (β a )− P (β * ) ≤ ε, where ψ = pL 2 pL 2 +γ 2 ,L = √ σ max L, σ max is the largest eigenvalue of X T X, L is the Lipschitz constant of f ′ , and β 0 denotes the starting point of the primal variables.
Proof: With L as the Lipschitz constant of f ′ ,L = √ σ max L is the Lipschitz constant of X T f ′ . Following the proof of Theorem 8 by Li et al. (2017), we have We have Hence, Recursively applying (25) k times, we have The algorithm reaches the desired accuracy smaller than ε, and it means P (β k )−P (β * ) ≤ ε. We can set Hence, for any iteration number k ≥ log ψ ε P (β 0 )−P (β * ) , we always have the primal gap Lemma 2 For the primal problem (1) and dual problem (2) , and θ k = τ kθk , with a large enough k in coordinate minimization, we have ||θ k − θ * || 2 2 ≤ LM λ 2 ||β k − β * || 2 Σ , where Σ = X T X, and M is a constant value.
L is the Lipschitz continuous constant of f ′ . As f is convex and smooth, there always is a point, e.g.β, with f ′ (Xβ) = 0. Hence, With β k andβ being finite values, ||f ′ (Xβ k )|| 2 is also a finite value. We can make the following assumption: Hence, (Q.E.D.) Here η = 1 + 1 C + u Cn . With ǫ D = G Z , ignoring the last two items, the complexity of dynamic screening can be simplified as O uL 2 γ 2 p log G 0 ε D + |Ā| log ε D ε . (Q.E.D.) Lemma 3 With O(u) as the complexity for the base operation of coordinate minimization of the LASSO problem with a γ-convex loss function, p A the total number of features involved in ADD operations, and p p A the size of the active set when IsADD is set to false, the complexity for the feature recruiting phase in SAIF is Proof : Let p A be the total number of features involved in the ADD operations. Q d (β) = P d (β) − P d (β * d ). We add up the time complexity of the outer loops regarding each added feature. O(Ku) is the time complexity for K base CM operations; O(np) is the computation complexity for duality gap and the ADD operation in one iteration of outer loop. We have With ψ d+1 = γ 2 p d+1L 2 d+1 +γ 2 , by following the proof of Theorem 3 in Li et al. (2017), As ∀d,L d+1 ≤L, the time complexity for the feature recruiting phase in SAIF Theorem 5 With O(u) as the complexity for the base operation of coordinate minimization of the LASSO problem with a γ-convex loss function, the time complexity for SAIF is O uL 2 γ 2 p logQ ε D +pp A + |Ā| log ε D ε . Here p A is the total number of features involved in ADD operations,p is the maximum size of the active set during the algorithm iterations, Q is the geometric mean of the accuracies of the sub-problem objective function values corresponding to each ADD operation, and ε D is the accuracy of the objective function value for the last feature DEL operation. and Let ǫ D = G Z D , the time complexity for SAIF can be simplified as O uL Here the primal and dual relation for the corresponding optima is θ * = − f (Xβ * ) λ .
From the objective function in (34), we can derive τ = y Tθ λ||θ|| 2 2 to reach the minimum distance toθ * if no constraint is imposed on τ . Therefore we need to estimate the range of τ to determine the minimum. From the constrained region |x T i τθ| ≤ 1, ∀i ∈ {1, ..., p − 1} , the range for τ is −