Consider the problem of sparse linear regression:
where is the design matrix (also called feature matrix) for samples and features and
is the corresponding vector ofresponses or observations. We assume without loss of generality that the columns of the matrix are normalized to , and the response vector is also normalized.
Given a subset of features (denoted by ), it is easy to find the best regression coefficients by projecting on the span of . The statistic (coefficient of determination
) measures the proportion of the variance explained by this subset. Sparse regression can be therefore seen as maximizing aset function : for any given set of columns , measures how well these columns explain the observations .
This set function is monotone increasing: including more features can only increase the explained variance111However, adjusted is not monotonic.. Solving sparsity constrained maximization of would involve searching over all subsets of size up to and selecting the one that maximizes
. This combinatorial optimization problem is unfortunately NP-Hard[Kempe:2011ue] . A widely used approach is to use greedy forward selection: select one feature at a time, greedily choosing the one that maximizes in the next step. Orthogonal Matching Pursuit and variations [Tropp:2004gc] [Needell08cosamp:iterative] can also be seen as approximate accelerated greedy forward selection algorithms that avoid re-fitting the model for each candidate vector at each step.
When maximizing set functions using the greedy algorithm it is natural to consider the framework of submodularity. In a classical result, Nemhauser:1978ck show that for submodular monotone functions, the greedy -sparse solution is within of the optimal -sparse solution, i.e., greedy gives a constant multiplicative factor approximation guarantee. The concept of submodularity has led to several effective greedy solutions to problems such as sparse prediction [koyejo2014], sparse factor analysis [khanna2015], model interpretation [kim:2016ee], etc.
Unfortunately, it is easy to construct counterexamples to show that is not submodular [Kempe:2011ue, Elenberg:2016tq]. In their breakthrough paper Kempe:2011ue show that if the design matrix satisfies the Restricted Isometry Property (RIP), then the set function satisfies a weakened form of submodularity. This weak form of submodularity is obtained by bounding a quantity called a submodularity ratio defined subsequently. The authors further showed that a bounded submodularity ratio is sufficient for Nemhauser’s proof to go through with a relaxed approximation constant (that depends on ). Therefore, even weak submodularity implies a constant factor approximation for greedy and RIP implies weak submodularity.
The weak submodularity framework was recently extended beyond linear regression by Elenberg:2016tq to concave functions (for example, likelihood functions of generalized linear models). This generalization of the RIP condition is called Restricted Strong Convexity (RSC) and the general result is that RSC implies weak submodularity for the set function obtained from the likelihood of the model restricted to subsets of features. This shows that greedy feature selection can obtain constant factor approximation guarantees in a very general setting, similar to results obtained by Lasso but without further statistical assumptions.
Running the greedy algorithm can be computationally expensive for large datasets. This is because for each greedy step we need to refit the model using the previously selected choices and the new candidate. This has led to the development of faster variants. For example, DistributedGreedy [Mirzasoleiman:2013vi] distributes the computational effort across available machines, and StochasticGreedy [Mirzasoleiman:2015th] exploits a stochastic greedy step. Both algorithms have provable performance guarantees for submodular functions.
In this paper, we show that submodularity is not required to obtain approximation guarantees for these two algorithms and that a bounded submodularity ratio suffices. This extends the scope of these algorithms to significantly larger class of functions like sparse linear regression with RIP design matrices. Furthermore, as we shall discuss, submodularity ratio can be used to provide data dependent bounds. This implies that one can sometimes obtain tighter guarantees also for submodular functions.
Our contributions are as follows: (1) We analyze and obtain approximation guarantees for DistributedGreedy and StochasticGreedy using the submodularity ratio extending the scope of their application to non-submodular functions,
(2) We show that the submodularity ratio can give tighter data dependent bounds even for submodular functions,
(3) We derive explicit bounds for both DistributedGreedy and StochasticGreedy for the special case of sparse linear regression.
(4) We derive bounds for sparse support selection for general concave functions.
(5) We also present empirical evaluations of these algorithms vs several established baselines on simulated and real world datasets.
Kempe:2011ue defined submodularity ratio, and showed that for any function that has its submodularity ratio bounded away from , one can provide appropriate greedy guarantees. Elenberg:2016tq used the concept to derive a new relationship between submodularity and convexity, specifically stating that the Restricted Strong Concavity can be used to bound the submodularity ratio. This results in providing bounds for greedy support selection for general concave functions.
Another notion of approximate additive submodularity was explored by Krause:2010tj for the problem of dictionary selection. This was, however, superceded by [Kempe:2011ue] who showed that submodularity ratio provides tighter approximation bounds. Horel:2016mas consider another generalization from submodular functions – -approximate submodular functions which are functions within of some submodular function, and provide approximation guarantees for greedy maximization.
The DistributedGreedy algorithm was introduced by Mirzasoleiman:2013vi. They provide deterministic bounds for any arbitrary distribution of data onto the individual machines. Barbosa:2015vq showed that for sparsity constraints, and under the assumption that the data is split uniformly at random to all the machines, one can obtain a guarantee in expectation. kumar2010 also provide distributed algorithms for maximizing a monotone submodular function subject to a sparsity constraint. They extend the Threshold Greedy algorithm of gupta2010 by augmenting it with a sample and prune strategy. It runs the Threshold Greedy algorithm on a subset of data to obtain a candidate solution. The latter is then used to prune the remaining data to reduce its size. This process is repeated a constant number of times, and the algorithm provides a constant factor guarantee
Altschuler2016 provide approximation guarantees for greedy selection variants for column subset selection. They do not use the submodularity framework, and their results are not directly applicable or useful for other problem settings. Similarly, farahatEGK13 also use greedy column subset selection. However, their focus is not towards obtaining approximation guarantees, but rather on more efficient algorithmic implementation.
Notation: We represent vectors as small letter bolds e.g. . Matrices are represented by capital bolds e.g. . Matrix or vector transposes are represented by superscript . Identity matrices of size are represented by , or simply when the dimensions are obvious. is a column vector of all ones (zeroes). Sets are represented by sans serif fonts e.g. , complement of a set is . For a vector , and a set of support dimensions with , denotes subvector of supported on . Similarly, for a matrix , denotes the submatrix supported on . We denote as .
Throughout this manuscript, we assume the set function is monotone. Our goal is to maximize under a cardinality constraint :
We begin by defining the submodularity ratio of a set function .
Definition 1 (Submodularity Ratio [Kempe:2011ue] ).
Let be two disjoint sets, and . The submodularity ratio for with respect to is given by
The submodularity ratio of a set with respect to an integer is given by
It is straightforward to show that is submodular if and only if for all sets and . Generalizing to the functions with provides a notion of weak submodularity [Elenberg:2016tq]. For weakly submodular functions, even though the function may not be submodular, it still provides approximation guarantees for the greedy algorithm. We use the submodularity ratio to provide new bounds for DistributedGreedy and StochasticGreedy, thereby generalizing these algorithms to non-submodular functions.
2.1 Greedy Selection
We briefly go over the classic greedy algorithm for subset selection. A greedy approach to optimizing a set function is myopic – the algorithm chooses the element from the available choices that gives the largest incremental gain for the set of choices previously made. The algorithm is illustrated in Algorithm 1. The algorithm makes outer iterations, where is the desired sparsity. Each iteration is a full pass over the remaining candidate choices, wherein the marginal gain is calculated for each remaining candidate choice. Thus the greedy algorithm has the computational complexity of calls to the function evaluation oracle.
For a large , the linear scaling of the greedy algorithm for a fixed may be prohibitive. As such, algorithms that scale sublinearly are useful for truly large scale selections. The DistributedGreedy algorithm, for example, achieves this sublinear scaling by making use of multiple machines. The data is split uniformly at random to machines. Each machine then performs its own independent greedy selection (Algorithm 1), and outputs a sized solution. All of the greedy solutions are collated by a central machine, which performs another round of the greedy selections to output the final solution. The algorithm is illustrated in Algorithm 2, and is analyzed in Section 3. It has a computational complexity of The algorithm is easy to implement in parallel or within a distributed computing framework e.g. MapReduce.
An alternative to distributing the data, say when several machines are not available, is to perform the greedy selection stochastically. The StochasticGreedy algorithm for submodular functions was introduced by Mirzasoleiman:2015th. At any given iteration , instead of performing a function evaluation for each of the remaining candidates, a subset of a fixed size (where
is a pre-specified hyperparameter) is uniformly sampled from the availablechoices using the subroutine Subsample, and the function evaluation is made on those subsampled choices as if they were the only available candidates. This speeds up the greedy algorithm to function evaluations. The algorithm is presented in Algorithm 3, and its approximation bounds are discussed in Section 4.
Finally, we discuss a property of the greedy algorithm, which is fundamental to the analysis of the DistributedGreedy algorithm. The greedy algorithm belongs to a larger class of algorithms called -nice algorithms [Mirrokni:2015kk]. The following result allows us to remove or add unselected items from the choice set that is accessible to the algorithm.
[Mirrokni:2015kk] Say , and let be the -sized set returned by Algorithm 1. For any , .
Note that Lemma 1 is a property of the algorithm, and is independent of the function itself. Prior works [Mirrokni:2015kk], [Barbosa:2015vq] have exploited this property in conjunction with properties of submodular functions to obtain approximation bounds for the distributed algorithms. Our work extends these results to weakly submodular functions. As such, it is easy to see that our results are easily extensible to other nice algorithms – including distributed OMP and distributed stochastic greedy – that have closed form bounds for the respective single machine algorithm. For ease of exposition, we focus our discussion on the distributed greedy algorithm.
3 Distributed Greedy
In this section, we obtain approximation bounds for DistributedGreedy (Algorithm 2). The algorithm returns the best out of solutions : the local solutions (steps 2,4), and the final aggregated one (step 3). Our strategy to obtain the approximation bound for the algorithm is as follows. To obtain an overall approximation bound, we obtain individual bounds on each of the solutions in terms of the submodularity ratio (Definition 1) and use the subadditivity ratio (Definition 2) to show that one of the two shall always hold. For approximation bounds on the local solutions, we make use of the niceness of the Greedy (Lemma 1). The bound on the aggregated solution is more involved, since it involves tracking the split of the true solution across machines. The assumption of partitioning uniformly at random is vital here. This helps us lower bound the greedy gain in expectation by a probabilistic overlap with the true solution. The trick of tracking the split of the true solution across machines is similar to the one that has been used for analysis of submodular functions [Mirrokni:2015kk], [Barbosa:2015vq], but without the explicit connection to submodularity and subadditivity ratios. As we shall see in Sections 5, 6 elucidating these connections leads to novel bounds for support selection for linear regression and general concave functions.
We next define the subadditivity ratio, which helps us generalize subadditive functions in the way similar to how submodularity ratio generalizes submodular functions.
Definition 2 (Subaddivity ratio).
We define the subadditivity ratio for a set function w.r.t a set as:
We further define the subadditivity ratio of a function for an integer , , which takes a uniform bound over all sets of size :
By definition, the function is subadditive iff . Since submodularaity implies subadditivity (the converse is not always true), if the function is submodular, .
We next present some notation and few lemmas that lead up to the main result of this section (Theorem 1). Let be the entire set of available choices. Partition the set uniformly at random into . Let be the -sized solution returned by running the greedy algorithm on i.e. . Note that each induces a partition onto the optimal -sized solution as follows:
Having defined the notation, we start by lower bounding the local solutions in terms of value of the subset of that is not selected as part of the respective local solution.
The next lemma is used to lower the bound the value of the aggregated solution (step 4 in Algorithm 2) in terms of the value of the subset that is selected as part of the respective local solution.
We are now ready to present our main result about the approximation guarantee for Algorithm 2.
Let be the set returned by the distributed greedy (Algorithm 2). Let . Then,
There are machines, each with its local greedy solution . In addition, there is the aggregated solution set . The key idea is to show that atleast one of the solutions is good enough.
Say for some , then by Lemma 2, .
On the other hand, say for all , , then for all , . By Lemma 3, the result then follows. ∎
Theorem 2 generalizes the approximation guarantee of obtained by Barbosa:2015vq for submodular functions. Their analysis uses convexity of the Lovasz extension of submodular functions, and hence can not be trivially extended to weakly submodular functions. In addition to being applicable for a larger class of functions, our result can also provide tighter bounds for specific applications or datasets even for submodular functions, since they are also applicable for submodular functions, and bounding and away from 1 from domain knowledge will give tighter approximations than the generic bound of .
4 Stochastic Greedy
For analysis of Algorithm 3, we show that the subsampling parameter governs the tradeoff between the speedup and the loss in approximation quality vis-a-vis the classic Greedy. Before formally providing the approximation bound, we present an auxillary lemma that is key to proving the new approximation bound. The following result is a generalization of Lemma 2 from [Mirzasoleiman:2015th] for weakly submodular functions.
Let , with . Consider another set drawn randomly from with . Then,
We are now ready to present our result that shows that stochastic greedy selections (Algorithm 3) can be applied to weakly submodular functions with provable approximation guarantees. All the proofs missing from the main text are presented in the supplement.
Let be the optimum set of size , and be the set built by StochasticGreedy at step . Then, .
Define . Using Lemma 4 with and , we get at the -th step,
Define , . Note that . Taking expectation on both sides over , (6) becomes
Using and above, along with the fact that for ,
Note that is the tradeoff hyperparameter between the speedup achieved by subsampling and the corresponding approximation guarantee. A larger value of means the algorithm is faster with weaker guarantees and vice versa. As , we tend towards the bound which recovers the bound for weakly submodular functions obtained by Kempe:2011ue for the classic greedy selections (Algorithm 1), and also recovers the well known bound of for submodular functions.
5 Large Scale Sparse Linear Regression
In this section, we derive novel bounds for greedy support selections for linear regression using both Algorithms 2, 3. Recall that is the feature matrix, with samples and features, and is the vector of responses. We assume, without loss of generality, that the columns of the matrix are normalized to , and the response vector is also normalized to have norm . Let be the covariance matrix.
We know from standard linear algebra, that for a fixed set of columns , where is the index into columns of , . Minimizing the error in (1) is thus equivalent to maximizing the following set function (modulo a constant ):
where is the projection matrix for orthogonal projection onto the span of columns of . as defined above is also the statistic for the linear regression problem (1). The respective combinatorial maximization is:
The function defined in (8) is not submodular. However, submodularity is not required for giving guarantees for greedy forward selection. A bounded submodularity ratio is a weaker condition that is sufficient to provide such approximation guarantees. Kempe:2011ue analyzed the greedy algorithm for (1), and showed that this function was weakly submodular. Our goal is to maximize the statistic for linear regression using Algorithms 2, 3.
For a positive semidefinite matrix , and be the largest and smallest
-sparse eigenvalue of. We make use of the following result from Kempe:2011ue:
For the statistic (7), .
Recall that we need to only bound the submodularity ratio over the selected greedy set to obtain the approximation bounds (See Theorems 1 and 2). Lemma 5 provides a data dependent union bound for the submodularity ratio in terms of sparse eigenvalue of the covariance matrix. We now provide the corresponding approximation bounds for Algorithm 3 next.
Let be the optimal support set for sparsity constrained maximization of the statistic (8). Let be the solution returned by StochasticGreedy (). Then,
For the maximization of the (7), , where is the submatrix of with rows and columns indexed by .
Let be the solution returned by the DistributedGreedy algorithm, and let be the optimal solution for the sparsity constrained maximization of (8). Then,
6 Support Selection for general functions
The sparsity constraint problem with a given for a concave function is:
Similar to the developments in Section 5, we can define an associated set function as:
We recall that the submodular guarantee of is for normalized submodular functions. To extend the notion of normalization to general support selection, we subtract .
To bound the submodularity ratio for in (10), the concept of strong concavity and smoothness is required. For a function , define . We say is Restricted Strongly Concave (RSC) and Restricted Strongly Smooth (RSM) over a subdomain if
We make use of the following result that lower bounds the submodularity ratio for :
Lemma 7 (Elenberg:2016tq).
If the given function is -strongly concave on all sparse supports and -smooth over all sparse supports,
6.1 RSC implies Weak Subadditivity
In this section, we establish a lower bound on the subadditivity ratio in terms of only the restricted strong concavity (RSC) and smoothness (RSM) constants. This is analogous to lower bounding the submodularity ratio by Elenberg:2016tq.
Say the function is strongly concave and -smooth over all supported on . Then, the subadditivity ratio can be lower bounded as:
To prove Theorem 4, we make use of the following two results. Recall that for a set and vector , denotes the subvector of supported on .
For a support set , .
For any support set , .
Let be a partition of a given support set i.e , .
We can use Lemma 8 to lower bound the numerator of the subadditivity ratio as follows:
Since we have a bound for the subadditivity ratio for general strongly concave and smooth functions, we can now provide a novel approximation guarantee for support selection by DistributedGreedy.
Say the function is -strongly concave over all sparse support sets and -smooth over all sparse support sets. Let be the optimal support set that maximizes the sparsity constrained (9). Let be the solution set returned by DistributedGreedy. Then,
To the best of our knowledge, the bounds obtained in Corollary 3 are the first for a distributed algorithm for support selection for general functions. Note that we have taken a uniform bound for restricted strong concavity and smoothness to be over all sized support sets, though it is only required to be over the optimal support set.
7.1 Distributed Linear Regression
We consider sparse linear regression in a distributed setting. We generate a -sparse regression vector is generated by selecting random nonzero entries of ,
where is a standard i.i.d. Gaussian. Measurements are taken according to , where is i.i.d. Gaussian with variance set to be . Each row of the design matrix is generated by an autoregressive process,
where is i.i.d. Gaussian with variance . We take for the number of both training and test measurements. Results are averaged over iterations, each with a different partition of the features.
We evaluate four variants of DistributedGreedy. The two greedy algorithms are Greedy Forward Selection (FS) and Orthogonal Matching Pursuit (OMP). Lasso sweeps an regularization parameter using LARS [EfronLARS]. This produces nested subsets of features corresponding to a sequence of thresholds for which the support size increases by . Lasso uses this threshold, while Lasso+ fits an unregularized linear regression on the support set selected by Lasso.
Figure 1 shows the performance of all algorithms on the following metrics: log likelihood (normalized with respect to a null model), generalization to new test measurements from the same true support parameter, area under ROC, and percentage of the true support recovered for .
Next, we run a similar experiment on a large, real-world dataset. We sample time series measurements across customers from the ElectricityLoadDiagrams time series dataset [UCI]
. We consider the supervised learning experiment of predicting the electrical load at thenext time . We use half of the customers for training and the rest for testing. Figure 2 shows performance of the same algorithms with data distributed across partitions to select the top features. We see that distributed Forward Selection produces both largest likelihood and highest generalization score. OMP has second largest likelihood, but its generalization varies widely for different values of . This is likely due to the random placement of predictive features across a large number of partitions.
7.2 Stochastic Sparse Logistic Regression
In this section we demonstrate the applicability of Algorithm 3 for greedy support selection for sparse logistic regression. Note that the respective set function (10) when is the log likelihood for logistic regression is not submodular. As such one would not typically apply the StochasticGreedy algorithm for sparse logistic regression. However, the guarantees obtained in Section 6 suggest good practical performance which is indeed demonstrated in Figure 3. For the experiment we use the gisette dataset obtained from the UCI website [UCI]. The dataset is of a handwritten digit recognition problem to separate out digits ‘’ and ‘’. It has 13500 instances, and 5000 features. Figure 3 illustrates depicts the tradeoff between the time taken to learn the model and the respective training log likelihood for different values of as used in Algorithm 3. As shown, we obtain tremendous speed up with relatively little loss in the log likelihood value even for reasonably large values.
We provided novel bounds for two greedy algorithm variants for maximizing weakly submodular functions with applications to linear regression and general concave functions. Our research opens questions on what other known algorithms with provable guarantees for submodular functions can be extended to weakly submodular functions.
Research supported by NSF Grants CCF 1344179, 1344364, 1407278, 1422549, IIS 1421729, and ARO YIP W911NF-14-1-0258.
Appendix A Proofs
In this section, we provide detailed proofs of all the results used in this manuscript. For lemma and theorem statements repeated from the main text, we add an apostrophe to indicate that it is not a new lemma/theorem being introduced.
We make use of the following result that provides an approximation bound for greedy selections for weakly submodular functions.
a.1 Distributed Greedy
From Lemma 1, we know that running greedy on instead of will still return the set ,
For proving Lemma 3, we require another auxillary result.
For any .
We now prove Lemma 3.
a.2 Stochastic Greedy
Let , with . Consider another set drawn randomly from with . Then
a.3 Linear regression
For the maximization of the (7) , where is the submatrix of with rows and columns indexed by .
Say , is an arbitrary partition of . Consider,
where (A.3) results from the fact that all orthogonal projection matrices are symmetric and idempotent. Repeating a similar analysis for instead of , we get
A similar analysis also gives , which gives the desired result. ∎
a.4 RSC implies weak subadditivity
For a support set , .
For any with support in ,
Using , we get the desired result.
For any support set , .
By strong concavity,
where is an arbitrary vector that has support only on . Optimizing the RHS over gives the desired result. ∎
Appendix B Additional experiments
Figure 4 shows the performance of all algorithms on the following metrics: log likelihood (normalized with respect to a null model), generalization to new test measurements from the same true support parameter, area under ROC, and percentage of the true support recovered for . Recall that Figure 1 presents the results from the same experiment with