Efficient Projections onto the 1 -Ball for Learning in High Dimensions John Duchi Google, Mountain View, CA 94043 Shai Shalev-Shwartz Toyota Technological Institute, Chicago, IL, 60637 Yoram Singer Tushar Chandra Google, Mountain View, CA 94043 J D U C H I @ C S . S TA N F O R D . E D U SHAI@TTI-C.ORG SINGER@GOOGLE.COM TUSHAR@GOOGLE.COM Abstract We describe efficient algorithms for projecting a vector onto the 1 -ball. We present two methods for projection. The first performs exact projection in O(n) expected time, where n is the dimension of the space. The second works on vectors k of whose elements are perturbed outside the 1 -ball, projecting in O(k log(n)) time. This setting is especially useful for online learning in sparse feature spaces such as text categorization applications. We demonstrate the merits and effectiveness of our algorithms in numerous batch and online learning tasks. We show that variants of stochastic gradient projection methods augmented with our efficient projection procedures outperform interior point methods, which are considered state-of-the-art optimization techniques. We also show that in online settings gradient updates with 1 projections outperform the exponentiated gradient algorithm while obtaining models with high degrees of sparsity. 1. Introduction A prevalent machine learning approach for decision and prediction problems is to cast the learning task as penalized convex optimization. In penalized convex optimization we seek a set of parameters, gathered together in a vector w, which minimizes a convex objective function in w with an additional penalty term that assesses the complexity of w. Two commonly used penalties are the 1norm and the square of the 2-norm of w. An alternative Appearing in Proceedings of the 25 th International Conference on Machine Learning, Helsinki, Finland, 2008. Copyright 2008 by the author(s)/owner(s). but mathematically equivalent approach is to cast the problem as a constrained optimization problem. In this setting we seek a minimizer of the objective function while constraining the solution to have a bounded norm. Many recent advances in statistical machine learning and related fields can be explained as convex optimization subject to a 1-norm constraint on the vector of parameters w. Imposing an 1 constraint leads to notable benefits. First, it encourages sparse solutions, i.e a solution for which many components of w are zero. When the original dimension of w is very high, a sparse solution enables easier interpretation of the problem in a lower dimension space. For the usage of 1 -based approach in statistical machine learning see for example (Tibshirani, 1996) and the references therein. Donoho (2006b) provided sufficient conditions for obtaining an optimal 1 -norm solution which is sparse. Recent work on compressed sensing (Candes, 2006; Donoho, 2006a) further explores how 1 constraints can be used for recovering a sparse signal sampled below the Nyquist rate. The second motivation for using 1 constraints in machine learning problems is that in some cases it leads to improved generalization bounds. For example, Ng (2004) examined the task of PAC learning a sparse predictor and analyzed cases in which an 1 constraint results in better solutions than an 2 constraint. In this paper we re-examine the task of minimizing a convex function subject to an 1 constraint on the norm of the solution. We are particularly interested in cases where the convex function is the average loss over a training set of m examples where each example is represented as a vector of high dimension. Thus, the solution itself is a high-dimensional vector as well. Recent work on 2 constrained optimization for machine learning indicates that gradient-related projection algorithms are more efficient in approaching a solution of good generalization than second-order algorithms when the number of examples and Efficient Projections onto the 1 -Ball for Learning in High Dimensions the dimension are large. For instance, Shalev-Shwartz et al. (2007) give recent state-of-the-art methods for solving large scale support vector machines. Adapting these recent results to projection methods onto the 1 ball poses algorithmic challenges. While projections onto 2 balls are straightforward to implement in linear time with the appropriate data structures, projection onto an 1 ball is a more involved task. The main contribution of this paper is the derivation of gradient projections with 1 domain constraints that can be performed almost as fast as gradient projection with 2 constraints. Our starting point is an efficient method for projection onto the probabilistic simplex. The basic idea is to show that, after sorting the vector we need to project, it is possible to calculate the projection exactly in linear time. This idea was rediscovered multiple times. It was first described in an abstract and somewhat opaque form in the work of Gafni and Bertsekas (1984) and Bertsekas (1999). Crammer and Singer (2002) rediscovered a similar projection algorithm as a tool for solving the dual of multiclass SVM. Hazan (2006) essentially reuses the same algorithm in the context of online convex programming. Our starting point is another derivation of Euclidean projection onto the simplex that paves the way to a few generalizations. First we show that the same technique can also be used for projecting onto the 1 -ball. This algorithm is based on sorting the components of the vector to be projected and thus requires O(n log(n)) time. We next present an improvement of the algorithm that replaces sorting with a procedure resembling median-search whose expected time complexity is O(n). In many applications, however, the dimension of the feature space is very high yet the number of features which attain non-zero values for an example may be very small. For instance, in our experiments with text classification in Sec. 7, the dimension is two million (the bigram dictionary size) while each example has on average one-thousand non-zero features (the number of unique tokens in a document). Applications where the dimensionality is high yet the number of "on" features in each example is small render our second algorithm useless in some cases. We therefore shift gears and describe a more complex algorithm that employs redblack trees to obtain a linear dependence on the number of non-zero features in an example and only logarithmic dependence on the full dimension. The key to our construction lies in the fact that we project vectors that are the sum of a vector in the 1 -ball and a sparse vector--they are "almost" in the 1 -ball. In conclusion to the paper we present experimental results that demonstrate the merits of our algorithms. We compare our algorithms with several specialized interior point (IP) methods as well as general methods from the literature for solving 1 -penalized problems on both synthetic and real data (the MNIST handwritten digit dataset and the Reuters RCV1 corpus) for batch and online learning. Our projection based methods outperform competing algorithms in terms of sparsity, and they exhibit faster convergence and lower regret than previous methods. 2. Notation and Problem Setting We start by establishing the notation used throughout the paper. The set of integers 1 through n is denoted by [n]. Scalars are denoted by lower case letters and vectors by lower case bold face letters. We use the notation w b to designate that all of the components of w are greater than b. We use · as a shorthand for the Euclidean norm · 2. The other norm we use throughout the paper is the 1n norm of the vector, v 1 = i=1 |vi |. Lastly, we consider order statistics and sorting vectors frequently throughout this paper. To that end, we let v(i) denote the ith order statistic of v, that is, v(1) v(2) . . . v(n) for v Rn . In the setting considered in this paper we are provided with a convex function L : Rn R. Our goal is to find the minimum of L(w) subject to an 1 -norm constraint on w. Formally, the problem we need to solve is minimize L(w) s.t. w w 1 z . (1) Our focus is on variants of the projected subgradient method for convex optimization (Bertsekas, 1999). Projected subgradient methods minimize a function L(w) subject to the constraint that w X , for X convex, by generating the sequence {w(t) } via ( w (t) 2) - t (t) w(t+1) = X where (t) is (an unbiased estimate of) the (sub)gradient of L at w(t) and X (x) = argmin y { x - y | y X } is Euclidean projection of x onto X . In the rest of the paper, the main algorithmic focus is on the projection step (computing an unbiased estimate of the gradient of L(w) is straightforward in the applications considered in this paper, as is the modification of w(t) by (t) ). 3. Euclidean Projection onto the Simplex For clarity, we begin with the task of performing Euclidean projection onto the positive simplex; our derivation naturally builds to the more efficient algorithms. As such, the most basic projection task we consider can be formally described as the following optimization problem, in 1 wi = z , wi 0 . (3) w - v 2 s .t. 2 2 =1 minimize w Efficient Projections onto the 1 -Ball for Learning in High Dimensions where R is a Lagrange multiplier and Rn is a + vector of non-negative Lagrange multipliers. Differentiating with respect to wi and comparing to zero gives the dL optimality condition, dwi = wi - vi + - i = 0. The complementary slackness KKT condition implies that whenever wi > 0 we must have that i = 0. Thus, if wi > 0 we get that wi = vi - + i = vi - . (4) When z = 1 the above is projection onto the probabilistic simplex. The Lagrangian of the problem in Eq. (3) is n - i 1 L(w, ) = w-v 2+ wi - z ·w , 2 =1 I N P U T: A vector v Rn and a scalar z > 0 Sort v into µ : µ1 µ2 . . . µp j D j r 1 µr - z Find = max [n] : µj - j =1 O i 1 µi - z efine = =1 > 0 U T P U T: w s.t. wi = max {vi - , 0} Figure 1. Algorithm for projection onto the simplex. 4. Euclidean Projection onto the 1 -Ball We next modify the algorithm to handle the more general 1 -norm constraint, which gives the minimization problem minimn ze w - v 2 i 2 w R All the non-negative elements of the vector w are tied via a single variable, so knowing the indices of these elements gives a much simpler problem. Upon first inspection, finding these indices seems difficult, but the following lemma (Shalev-Shwartz & Singer, 2006) provides a key tool in deriving our procedure for identifying non-zero elements. Lemma 1. Let w be the optimal solution to the minimization problem in Eq. (3). Let s and j be two indices such that vs > vj . If ws = 0 then wj must be zero as well. Denoting by I the set of indices of the non-zero components of the sorted optimal solution, I = {i [n] : v(i) > 0}, we see that Lemma 1 implies that I = [] for some 1 n. Had we known we could have simply used Eq. (4) to obtain that in wi = in w(i) = i i w(i) = i v . (i) s .t. w 1 z . (7) We do so by presenting a reduction to the problem of projecting onto the simplex given in Eq. (3). First, we note that if v 1 z then the solution of Eq. (7) is w = v. Therefore, from now on we assume that v 1 > z . In this case, the optimal solution must be on the boundary of the constraint set and thus we can replace the inequality constraint w 1 z with an equality constraint w 1 = z . Having done so, the sole difference between the problem in Eq. (7) and the one in Eq. (3) is that in the latter we have an additional set of constraints, w 0. The following lemma indicates that each non-zero component of the optimal solution w shares the sign of its counterpart in v. Lemma 3. Let w be an optimal solution of Eq. (7). Then, for all i, wi vi 0. Proof. Assume by contradiction that the claim does not ^ hold. Thus, there exists i for which wi vi < 0. Let w be a vector such that wi = 0 and for all j = i we have ^ wj = wj . Therefore, w 1 = w 1 - |wi | z and hence ^ ^ ^ w is a feasible solution. In addition, w-v 2- w-v 2 ^ 2 2 = 2 2 = wi - 2wi vi > wi > 0 . =1 =1 =1 =1 - and therefore 1 = v(i) - z = z (5) =1 Given we can characterize the optimal solution for w as wi = max {vi - , 0} . (6) (wi - vi )2 - (0 - vi )2 We are left with the problem of finding the optimal , and ^ We thus constructed a feasible solution w which attains an the following lemma (Shalev-Shwartz & Singer, 2006) proobjective value smaller than that of w. This leads us to the vides a simple solution once we sort v in descending order. desired contradiction. Lemma 2. Let w be the optimal solution to the minimization problem given in Eq. (3). Let µ denote the vector obBased on the above lemma and the symmetry of the obtained by sorting v in a descending order. Then, the numjective, we are ready to present our reduction. Let u be a ber of strictly positive elements in w is vector obtained by taking the absolute value of each comj > . ponent of v, u = |v |. We now replace Eq. (7) with j i i 1r µr - z 0 (z , µ) = max [n] : µj - j minimize - u 2 s.t. 1 z and 0 . (8) 2 =1 n R The pseudo-code describing the O(n log n) procedure for solving Eq. (3) is given in Fig. 1. Once we obtain the solution for the problem above we construct the optimal of Eq. (7) by setting wi = sign(vi ) i . Efficient Projections onto the 1 -Ball for Learning in High Dimensions I N P U T A vector v Rn and a scalar z > 0 I N I T I A L I Z E U = [n] s = 0 = 0 WHILE U = P I C K k U at random PA RT I T I O N U : G = {j U | vj vk } L = {j U | vj < vk } j C A L C U L AT E = |G| ; s = vj IF G 2001). The algorithm computes partial sums just-in-time and has expected linear time complexity. The algorithm identifies and the pivot value v() without sorting the vector v by using a divide and conquer procedure. The procedure works in rounds and on each round either eliminates elements shown to be strictly smaller than v() or updates the partial sum leading to Eq. (9). To do so the algorithm maintains a set of unprocessed elements of v. This set contains the components of v whose relationship to v() we do not know. We thus initially set U = [n]. On each round of the algorithm we pick at random an index k from the set U . Next, we partition the set U into two subsets G and L. G contains all the indices j U whose components vj > vk ; L contains those j U such that vj is smaller. We now face two cases related to the current summation of entriejs in v greater than the hypothesized v() (i.e. vk ). If :vj vk (vj - vk ) < z then by Eq. (9), vk v() . In this case we know that all the elements in G participate in the sum defining as given by Eq. (9). We can discard G and set U to be L as we still need to further identify the remaining elements in L. If j :vj vk (vj - vk ) z then the same rationale implies that vk < v() . Thus, all the elements in L are smaller than v() and can be discarded. In this case we can remove the set L and vk and set U to be G \ {k }. The entire process ends when U is empty. Along the process we also keep track of the sum and the number of elements in v that we have found thus far to be no smaller than v() , which is required in order not to recalculate partial sums. The pseudo-code describing the efficient projection algorithm is provided in Fig. 2. We keep the set of elements found to be greater than v() only implicitly. Formally, at each iteration of the algorithm we maintain a variable s, which is the sum of the elements in the set {vj : j U, vj v() }, and overload to designate the cardinality of the this set throughout the algorithm. Thus, when the algorithms exits its main while loop, is the maximizer defined in Lemma 1. Once the while loop terminates, we are left with the task of calculating usij g Eq. (10) and performing the actual projection. Since n :vj µ vj is readily available to us as the variable s, we simply set to be (s - z )/ and perform the projection as prescribed by Eq. (6). Though omitted here for lack of space, we can also extend theaalgorithms to handle the more general constraint that i |wi | z for ai 0. ELSE U G \ {k } ENDIF S E T = (s - z )/ O U T P U T w s.t. vi = max {vi - , 0} (s + s) - ( + )vk < z s = s + s ; = + ; U L Figure 2. Linear time projection onto the simplex. 5. A Linear Time Projection Algorithm In this section we describe a more efficient algorithm for performing projections. To keep our presentation simple and easy to follow, we describe the projection algorithm onto the simplex. The generalization to the 1 ball can straightforwardly incorporated into the efficient algorithm by the results from the previous section (we simply work in the algorithm with a vector of the absolute values of v, replacing the solution's components wi with sign(vi ) · wi ). For correctness of the following discussion, we add another component to v (the vector to be projected), which we set to 0, thus vn+1 = 0 and v(n+1) = 0. Let us start by examining again Lemma 2. The lemma implies that the index is the largest integer that still satisfies > 1 0. After routine algebraic v() - r = 1 v( r ) - z manipulations the above can be rewritten in the following somewhat simpler form: i v - v() < z and +1 i =1 (i) =1 v (i) - v(+1) z . (9) Given and v() we slightly rewrite the value as follows, 1 j = vj - z . (10) : v j v ( ) The task of projection can thus be distilled to the task of finding , which in turn reduces to the task of finding and the pivot element v() . Our problem thus resembles the task of finding an order statistic with an additional complicating factor stemming from the need to compute summations (while searching) of the form given by Eq. (9). Our efficient projection algorithm is based on a modification of the randomized median finding algorithm (Cormen et al., 6. Efficient Projection for Sparse Gradients Before we dive into developing a new algorithm, we remind the reader of the iterations the minimization algorithm takes from Eq. (2): we generate a sequence {w(t) } Efficient Projections onto the 1 -Ball for Learning in High Dimensions I N P U T A balanced tree T and a scalar z > 0 I N I T I A L I Z E v = , = n + 1, s = z C A L L P I VOT S E A R C H(root(T ), 0, 0) P R O C E D U R E P I VOT S E A R C H(v , , s) C O M P U T E = + r(v ) ; s = s + (v ) ^ ^ I F s < v + z ^ ^ // v pivot I F v > v v = v ; = ; s = s ^ ^ ENDIF I F leafT (v ) R E T U R N = (s - z )/ ENDIF ^^ C A L L P I VOT S E A R C H(leftT (v ), , s) // v < pivot ELSE I F leafT (v ) R E T U R N = (s - z )/ EE N D I F r C A L L P I VOT S E A R C H ightT (v ), , s NDIF E N D P RO C E D U R E Figure 3. Efficient search of pivot value for sparse feature spaces. step in which we deduct t from all the elements of the vector implicitly, adhering to the goal of performing a sublinear number of operations. As before, we assume that the goal is to project onto the simplex. Equipped with these variables, the j th component of the projected vector after t projected gradient steps can be written as max{vj - t , 0}. The second substantial modification to the core algorithm is to keep only the non-zero components of the weight vector in a red-black tree (Cormen et al., 2001). The red-black tree facilitates an efficient search for the pivot element (v() ) in time which is logarithmic in the dimension, as we describe in the sequel. Once the pivot element is found we implicitly deduct t from all the non-zero elements in our weight vector by updating t . We then remove all the components that are less than v() (i.e. less than t ); this removal is efficient and requires only logarithmic time (Tarjan, 1983). The course of the algorithm is as follows. After t projected gradient iterations we have a vector v(t) whose non-zero elements are stored in a red-black tree T and a global deduction value t which is applied to each non-zero component just-in-time, i.e. when needed. Therefore, each non-zero weight is accessed as vj - t while T does not contain the zero elements of the vector. When updating v with a gradient, we modify the vector v(t) by adding to it the gradientbased vector g(t) with k non-zero components. This update is done using k deletions (removing vi from T such that (t) (t) gi = 0) followed by k re-insertions of vi = (vi + gi ) into T , which takes O(k log(n)) time. Next we find in O(log(n)) time the value of t . Fig. 3 contains the algorithm for this step; it is explained in the sequel. The last step removes all elements of the new raw vector v(t) + g(t) which become zero due to the projection. This step is discussed at the end of this section. In contrast to standard tree-based search procedure, to find t we need to find a pair of consecutive values in v that correspond to v() and v(+1) . We do so by keeping track of the smallest element that satisfies the left hand side of Eq. (9) while searching based on the condition given on the right hand side of the same equation. T is keyed on the values of the un-shifted vector vt . Thus, all the children in the left (right) sub-tree of a node v represent values in vt which are smaller (larger) than v . In order to efficiently find t we keep at each node the following information: (a) The value of the component, simply denoted as v . (b) The number of elements in the right sub-tree rooted at v , denoted r(v ), including the node v . (c) The sum of the elements in the right sub-tree rooted at v , denoted (v ), including the value v itself. Our goal is to identify the pivot element v() and its index . In the previous section we described a simple condition for checking whether an element in v is greater or smaller than the pivot value. We now rewrite this expression yet one more time. A component with value v is not by iterating w w(t+1) = W here g(t) = -t (t) , W = {w | w projection onto this set. w (t) + g(t) 1 z } and W is In many applications the dimension of the feature space is very high yet the number of features which attain a non-zero value for each example is very small (see for instance our experiments on text documents in Sec. 7). It is straightforward to implement the gradient-related updates in time which is proportional to the number of non-zero features, but the time complexity of the projection algorithm described in the previous section is linear in the dimension. Therefore, using the algorithm verbatim could be prohibitively expensive in applications where the dimension is high yet the number of features which are "on" in each example is small. In this section we describe a projection algorithm that updates the vector w(t) with g(t) and scales linearly in the number of non-zero entries of g(t) and only logarithmically in the total number of features (i.e. non-zeros in w(t) ). The first step in facilitating an efficient projection for sparse feature spaces is to represent the projected vector as a "raw" vector v by incorporating a global shift that is applied to each non-zero component. Specifically, each projection step amounts to deducting from each component of v and thresholding the result at zero. Let us denote by t the shift value used on the tth iteration of the algorithm and by s t the cumulative sum of the shift values, t = t s . The representation we employ enables us to perform the Efficient Projections onto the 1 -Ball for Learning in High Dimensions smaller than the pivot iff the following holds: j vj > |{j : vj v }| · v + z . (11) * 10 2 10 1 f-f f-f : v j v Coordinate L1 - Line L1 - Batch L1 - Stoch IP * 10 3 10 2 Coordinate L1 - Batch L1 - Stoch IP 10 1 The variables in the red-black tree form the infrastructure for performing efficient recursive computation of Eq. (11). Note also that the condition expressed in Eq. (11) still holds when we do not deduct t from all the elements in v. The search algorithm maintains recursively the number and the sum s of the elements that have been shown to be greater or equal to the pivot. We start the search with the root node of T , and thus initially = 0 and s = 0. Upon entering a new node v , the algorithm checks whether the condition given by Eq. (11) holds for v . Since and s were computed for the parent of v , we need to incorporate the number and the sum of the elements that are larger than v itself. By construction, these variables are r(v ) and (v ), which we store at the node v itself. We let = + r(v ) ^ and s = s + (v ), and with these variables handy, Eq. (11) ^ distills to the expression s < v + z . If the inequality holds, ^ ^ we know that v is either larger than the pivot or it may be the pivot itself. We thus update our current hypothesis for µ and (designated as v and in Fig. 3). We continue searching the left sub-tree (leftT (v )) which includes all elements smaller than v . If inequality s < v + z does not ^ ^ hold, we know that v < µ , and we thus search the right subtree (rightT (v )) and keep and s intact. The process naturally terminates once we reach a leaf, where we can also calculate the correct value of using Eq. (10). Once we find t (if t 0) we update the global shift, t+1 = t + t . We need to discard all the elements in T smaller than t+1 , which we do using Tarjan's (1983) algorithm for splitting a red-black tree. This step is logarithmic in the total number of non-zero elements of vt . Thus, as the additional variables in the tree can be updated in constant time as a function of a node's child nodes in T , each of the operations previously described can be performed in logarthmic time (Cormen et al., 2001), giving us a total update time of O(k log(n)). 10 0 10 0 10 -1 -1 10 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 x 10 8 1 2 3 4 5 6 x 10 9 Approximate Flops Approximate Flops Figure 4. Comparison of methods on 1 -regularized least squares. The left has dimension n = 800, the right n = 4000 methods for both least squares and logistic regression (Koh et al., 2007; Kim et al., 2007). The algorithms we use are batch projected gradient, stochastic projected subgradient, and batch projected gradient augmented with a backtracking line search (Koh et al., 2007). The IP and coordinatewise methods both solve regularized loss functions of the form f (w) = L(w) + w 1 rather than having an 1 domain constraint, so our objectives are not directly comparable. To surmount this difficulty, we first minimize L(w) + w 1 and use the 1-norm of the resulting solution w as the constraint for our methods. To generate the data for the least squares problem setting, we chose a w with entries distributed normally with 0 mean and unit variance and randomly zeroed 50% of the vector. The data matrix X Rm×n was random with entries also normally distributed. To generate target values for the least squares problem, we set y = X w + , where the components of were also distributed normally at random. In the case of logistic regression, we generated data X and the vector w identically, but the targets yi were set to be sign(w · xi ) with probability 90% and to -sign(w · xi ) otherwise. We ran two sets of experiments, one each for n = 800 and n = 4000. We also set the number of examples m to be equal to n. For the subgradient methods in these e eriments and throughout the remainder, we set xp t = 0 / t, choosing 0 to give reasonable performance. (0 too large will mean that the initial steps of the gradient method are not descent directions; the noise will qui ly ck disappear because the step sizes are proportional to 1/ t). Fig. 4 and Fig. 5 contain the results of these experiments and plot f (w) - f (w ) as a function of the number of floating point operations. From the figures, we see that the projected subgradient methods are generally very fast at the outset, getting us to an accuracy of f (w) - f (w ) 10-2 quickly, but their rate of convergence slows over time. The fast projection algorithms we have developed, however, allow projected-subgradient methods to be very competitive with specialized methods, even on these relatively small problem sizes. On higher-dimension data sets interior point methods are infeasible or very slow. The rightmost graphs in Fig. 4 and Fig. 5 plot f (w) - f (w ) as functions of floating point operations for least squares and logistic re- 7. Experiments We now present experimental results demonstrating the effectiveness of the projection algorithms. We first report results for experiments with synthetic data and then move to experiments with high dimensional natural datasets. In our experiment with synthetic data, we compared variants of the projected subgradient algorithm (Eq. (2)) for 1 -regularized least squares and 1 -regularized logistic regression. We compared our methods to a specialized coordinate-descent solver for the least squares problem due to Friedman et al. (2007) and to very fast interior point Efficient Projections onto the 1 -Ball for Learning in High Dimensions L1 - Line L1 - Batch L1 - Stoch IP * 10 -1 10 -1 10 10 -2 -1 L1 - Stoch L1 - Full EG - Full EG - Stoch f-f * 10 -1 10 -2 L1 - 1 EG - 1 L1 - 100 EG - 100 * f-f f-f 10 -3 f-f * 10 -2 10 -2 10 -3 10 -3 10 -4 -4 L1 - Batch L1 - Stoch IP 10 10 -3 -4 10 10 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 x 10 8 -5 1 2 3 4 5 6 7 8 x 10 9 0 20 40 60 80 100 120 0 50 100 150 200 250 300 350 400 450 Approximate Flops Approximate Flops Time (CPU seconds) Time (CPU seconds) Figure 5. Comparison of methods on 1 -regularized logistic regression. The left has dimension n = 800, the right n = 4000 Figure 6. EG and projected subgradient methods on RCV1. gression with dimension n = 4000. These results indicate that in high dimensional feature spaces, the asymptotically faster convergence of IP methods is counteracted by their quadratic dependence on the dimension of the space. We also ran a series of experiments on two real datasets with high dimensionality: the Reuters RCV1 Corpus (Lewis et al., 2004) and the MNIST handwritten digits database. The Reuters Corpus has 804,414 examples; with simple stemming and stop-wording, there are 112,919 unigram features and 1,946,684 bigram features. With our preprocessing, the unigrams have a sparsity of 1.2% and the bigrams have sparsity of .26%. We performed 1 -constrained binary logistic regression on the CCAT category from RCV1 (classifying a document as corporate/industrial) using unigrams in a batch setting and bigrams in an online setting. The MNIST dataset consists of 60,000 training examples and a 10,000 example test set and has 10-classes; each image is a gray-scale 28 × 28 image, which we represent as xi R784 . Rather than directly use the input xi , however, we learned weights wj using the following Kernel-based "similarity" function for each class j {1, . . . , 10}: 1 i if yi = j k (x, j ) = wj i j i K (xi , x), j i = -1 otherwise. S mensional spaces, and it very quickly identifies and shrinks weights for irrelevant features (Kivinen & Warmuth, 1997). At every step of EG we update -Z (t) wi exp t i f (w(t) ) (t+1) (13) wi = t where Zt normalizes so that wi = z and i f denotes the ith entry of the gradient of f , the function to be minimized. EG can actually be viewed as a projected subgradient mi thod using generalized relative ene i tropy (D(x y) = xi log xi - xi + yi ) as the distance y function for projections (Beck & Teboulle, 2003). We can ^ replace i f with i f in Eq. (13), an unbiased estimator of e gradient of f , to get stochastic EG. A step size t th l 1/ t guarantees a convergence rate of O( og n/T ). For each experiment with EG, however, wexperimented with e learning rates proportional to 1/t, 1/ t, and constant, as well as different initial step-sizes; to make EG as competitive as possible, we chose the step-size and rate for which EG performed best on each individual test.. Results for our batch experiments learning a logistic classifier for CCAT on the Reuters corpus can be seen in Fig. 6. The figure plots the binary logistic loss of the different algorithms minus the optimal log loss as a function of CPU time. On the left side Fig. 6, we used projected gradient descent and stochastic gradient descent using 25% of the training data to estimate the gradient, and we used the algorithm of Fig. 2 for the projection steps. We see that 1 projections outperform EG both in terms of convergence speed and empirical log-loss. On the right side of the figure, we performed stochastic descent using only 1 training example or 100 training examples to estimate the gradient, using Fig. 3 to project. When the gradient is sparse, updates for EG are O(k ) (where k is the number of non-zeros in the gradient), so EG has a run-time advantage over 1 projections when the gradient is very sparse. This advantage can be seen in the right side of Fig. 6. For MNIST, with dense features, we ran a similar series of tests to those we ran on the Reuters Corpus. We plot the multiclass logistic loss from Eq. (12) over time (as a function of the number gradient evaluations) in Fig. 7. The left side of Fig. 7 compares EG and gradient descent using i (t+1) In the above, K is a Gaussian kernel function, so that K (x, y) = exp(- x - y 2/25), and S is a 2766 element support set. We put an 1 constraint on each wj , giving us the following multiclass objective with dimension 27,660: s 1 m r 1 k(xi ,r )-k(xi ,yi ) minimizew m i=1 log + =yi e .t. wj 1 z , wj 0. (12) As a comparison to our projected subgradient methods on real data, we used a method known in the literature as either entropic descent, a special case of mirror descent (Beck & Teboulle, 2003), or exponentiated gradient (EG) (Kivinen & Warmuth, 1997). EG maintains a weight vector w subi ject to the constraint that wi = z and w 0; it can easily be extended to work with negative weights under a 1-norm constraint by maintaining two vectors w+ and w- . We compare against EG since it works well in very high di- Efficient Projections onto the 1 -Ball for Learning in High Dimensions EG L1 0 EG L1 ECAT (versus more than 15% and 20% respectively for EG). 10 f-f f-f 10 0 Acknowledgments -1 * * 10 10 -1 We thank the anonymous reviewers for their helpful and insightful comments. 50 100 150 200 250 300 350 400 2 4 6 8 10 12 14 16 18 20 Gradient Evaluations Stochastic Subgradient Evaluations Figure 7. MNIST multiclass logistic loss as a function of the number of gradient steps. The left uses true gradients, the right stochastic subgradients. x 10 3.5 5 References Beck, A., & Teboulle, M. (2003). Mirror descent and nonlinear projected subgradient methods for convex optimization. Operations Research Letters, 31, 167­175. Bertsekas, D. (1999). Nonlinear programming. Athena Scientific. Candes, E. J. (2006). Compressive sampling. Proc. of the Int. Congress of Math., Madrid, Spain. Cormen, T. H., Leiserson, C. E., Rivest, R. L., & Stein, C. (2001). Introduction to algorithms. MIT Press. Crammer, K., & Singer, Y. (2002). On the learnability and design of output codes for multiclass problems. Machine Learning, 47. Donoho, D. (2006a). Compressed sensing. Technical Report, Stanford University. Donoho, D. (2006b). For most large underdetermined systems of linear equations, the minimal 1 -norm solution is also the sparsest solution. Comm. Pure Appl. Math. 59. Friedman, J., Hastie, T., & Tibshirani, R. (2007). Pathwise coordinate optimization. Annals of Applied Statistics, 1, 302­332. Gafni, E., & Bertsekas, D. P. (1984). Two-metric projection methods for constrained optimization. SIAM Journal on Control and Optimization, 22, 936­964. Hazan, E. (2006). Approximate convex optimization by online game playing. Unpublished manuscript. Kim, S.-J., Koh, K., Lustig, M., Boyd, S., & Gorinevsky, D. (2007). An interior-point method for large-scale 1 -regularized least squares. IEEE Journal on Selected Topics in Signal Processing, 4, 606­617. Kivinen, J., & Warmuth, M. (1997). Exponentiated gradient versus gradient descent for linear predictors. Information and Computation, 132, 1­64. Koh, K., Kim, S.-J., & Boyd, S. (2007). An interior-point method for large-scale 1 -regularized logistic regression. Journal of Machine Learning Research, 8, 1519­1555. Lewis, D., Yang, Y., Rose, T., & Li, F. (2004). Rcv1: A new benchmark collection for text categorization research. Journal of Machine Learning Research, 5, 361­397. Ng, A. (2004). Feature selection, l1 vs. l2 regularization, and rotational invariance. Proceedings of the Twenty-First International Conference on Machine Learning. Shalev-Shwartz, S., & Singer, Y. (2006). Efficient learning of label ranking by soft projections onto polyhedra. Journal of Machine Learning Research, 7 (July), 1567­1599. Shalev-Shwartz, S., Singer, Y., & Srebro, N. (2007). Pegasos: Primal estimated sub-gradient solver for SVM. Proceedings of the 24th International Conference on Machine Learning. Tarjan, R. E. (1983). Data structures and network algorithms. Society for Industrial and Applied Mathematics. Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. Royal. Statist. Soc B., 58, 267­288. Zinkevich, M. (2003). Online convex programming and generalized infinitesimal gradient ascent. Proceedings of the Twentieth International Conference on Machine Learning. 7 3 Cumulative Loss 2.5 EG - CCAT EG - ECAT L1 - CCAT L1 - ECAT % Sparsity 6 5 2 4 1.5 3 1 2 0.5 1 % of Total Features % of Total Seen 0 1 2 3 4 5 6 7 x 10 8 5 0 1 2 3 4 5 6 7 8 x 10 5 0 Training Examples Training Examples Figure 8. Online learning of bigram classifier on RCV1. Left is the cumulative loss, right shows sparsity over time. the true gradient while the right figure compares stochastic EG and stochastic gradient descent using only 1% of the training set to estimate the gradient. On top of outperforming EG in terms of convergence rate and loss, the 1 projection methods also gave sparsity, zeroing out between 10 and 50% of the components of each class vector wj in the MNIST experiments, while EG gives no sparsity. As a last experiment, we ran an online learning test on the RCV1 dataset using bigram features, comparing 1 projections to using decreasing step sizes given by Zinkevich (2003) to exponential gradient updates. The 1 projections are computationally feasible because of algorithm 3, as the dimension of our feature space is nearly 2 million (using the expected linear-time algorithm of Fig. 2 takes 15 times as long to compute the projection for the sparse updates in online learning). We selected the bound on the 1-norm of the weights to give the best online regret of all our experiments (in our case, the bound was 100). The results of this experiment are in Fig. 8. The left figure plots the cumulative log-loss for the CCAT and ECAT binary prediction problems as a function of the number of training examples, while the right hand figure plots the sparsity of the 1 -constrained weight vector both as a function of the dimension and as a function of the number of features actually seen. The 1 -projecting learner maintained an active set with only about 5% non-zero components; the EG updates have no sparsity whatsoever. Our online 1 -projections outperform EG updates in terms of the online regret (cumulative log-loss), and the 1 -projection updates also achieve a classification error rate of 11.9% over all the examples on the CCAT task and 14.9% on