Fast Projection onto the Capped Simplex with Applications to Sparse Regression in Bioinformatics
Year: 2021

Abstract:
We consider the problem of projecting a vector onto the so-called k-capped simplex, which is a hyper-cube cut by a hyperplane. For an n-dimensional input vector with bounded elements, we found that a simple algorithm based on Newton’s method is able to solve the projection problem to high precision with a complexity roughly about O(n), which has a much lower computational cost compared with the existing sorting-based methods proposed in the literature. We provide a theory for partial explanation and justification of the method. We demonstrate that the proposed algorithm can produce a solution of the projection problem with high precision on large scale datasets, and the algorithm is able to significantly outperform the state-of-the-art methods in terms of runtime (about 6-8 times faster than a commercial software with respect to CPU time for input vector with 1 million variables or more). We further illustrate the effectiveness of the proposed algorithm on solving sparse regression in a bioinformatics problem. Empirical results on the GWAS dataset (with 1,500,000 single-nucleotide polymorphisms) show that, when using the proposed method to accelerate the Projected Quasi-Newton (PQN) method, the accelerated PQN algorithm is able to handle huge-scale regression problem and it is more efficient (about 3-6 times faster) than the current state-of-the-art methods.

heading: 1 Introduction
text: The k-capped simplex [17] is defined as k := x 2 Rn x>1  k, 0  x  1
, where k 2 R is an input parameter, 0 and 1 are vectors of zeros and vector of ones in Rn, respectively (resp.), and the sign  is taken element-wise. Geometrically, for some k > 0, the set k represents an n-dimensional hyper-cube [0,1]n with a corner chopped off by the half-space x>1  k. If the half-space constraint x>1  k in k is replaced by the hyper-plane constraint x>1 = k, for some k > 0, this set represents the slice of the hyper-cube cut by the hyper-plane, and we denote this set as =k . In this paper, we consider the problem of projecting a vector y onto the k-capped simplex:
x⇤ = proj k(y) = argmin x n 1 2 kx yk22 s.t. x 2 k (or = k ) o . (1)
⇤Both authors contributed equally. 35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Remark 1. (On the input k) In this paper, we assume 0 < k < 1 is properly selected such that the set k (and =k ) has a nonempty interior, and therefore the projection is feasible. I.e., we assume p 2 Rn | p>1  k \ q 2 Rn | 0  q  1
6= ? (and similarly for =k ). In this setting, a
unique global minimizer x⇤ of Problem(1) exists, since the function kx yk22 is strongly convex. The value of k is closely related to the root of a piecewise-linear equation that plays a critical role in the solution of Problem (1), see discussion in § 2.
Problem (1) generalizes the projection onto the probability simplex [4, 5] which has many applications. In this paper, we consider the application in sparse regression in bioinformatics, see § 4. To the best of our knowledge, for an input vector y 2 Rn, the state-of-the-art algorithm [17] for solving the Projection (1) has a time complexity of O(n2). We found that a simple heuristic algorithm based on Newton’s method is able to solve the same problem to high precision with a complexity roughly O(n) on average, based on exploiting the 2nd-order information of a scalar minimization problem on the Lagrangian multiplier of Problem (1); see § 2.3.
Contributions We proposed a simple algorithm based on Newton’s method, as an alternative to the sorting-based method to solve Problem (1). We first convert the problem into a scalar minimization problem, then we show that such scalar problem can be solved using Newton’s iteration, see Theorem 1. We then show that the proposed method is superior to the current state-of-the-art method in terms of runtime. For a vector with dimension 106, the proposed method is about 6-8 times faster than a commercial solver. Furthermore, we give numerical results to show the effectiveness of the proposed method on solving sparse regression in a bioinformatics problem on a real-world dataset with one million data points, see § 4.
Scope, limitation, and impact We emphasize that the proposed method is not aimed at replacing all the well-established sorting-based methods in the literature, but to provide an alternative method for solving the projection problem in a special case. To be specific, we found that, when the input vector y has bounded elements (i.e., |yi|  ↵, for ↵ up to 1000), the proposed algorithm in this paper is able to solve the projection problem much faster than the existing methods in the literature. It is important to note that, when the bound ↵ becomes very huge (say a million), the proposed method sometimes converges slower than existing methods. Although limited in the scope of application, however, we argue that the proposed method is still very useful for applications that the input vectors are bounded. Using bioinformatics as an example, on an ordinary PC, we are able to reduce the runtime on identifying the genetic factors that contribute to dental caries (i.e., tooth decay) from half-hour to just 3 minutes; see § 4.
Paper organization In § 2, we give the background material, show the main theorem, and discuss the proposed algorithm. In § 3, we show experimental results to verify the efficiency of the proposed algorithm on solving Problem (1). In § 4, we illustrate the effectiveness of the proposed method on solving sparse regression in a large bioinformatics dataset. We conclude the paper in § 6 and provide supplementary material in the appendix.

heading: 2 Projection onto the capped simplex
text: In this section, we discuss how to solve Problem (1). By introducing a Lagrangian multiplier , we convert (1) to a scalar problem in the form of min !( ), see (5). Then, we provide a theorem for characterizing the function !( ). More specifically, we show that !( ) is n-smooth, i.e., the Lipschitz constant of the derivative !0( ) grows linearly as the dimension n, which makes the gradient method ineffective to solve Problem (5) due to very small stepsize for input vector y that has many variables. This motivates us to switch to the 2nd-order method, i.e., Newton’s method, to solve Problem (5). We end this section by discussing the properties of proposed Newton’s method.

heading: 2.1 Reformulate the projection problem to a scalar minimization problem
text: Consider the case of projection onto =k
2. Problem (1) is equivalent to Problem (2): by introducing a Lagrangian multiplier that corresponds to the equality constraint, we form a partial Lagrangian L(x, ), and arrive at the min-max problem (2), where we keep the inequality constraint on x as
min 0x1 max 2R
: n L(x, ) = 1
2 kx yk2 + (x>1 k)
o . (2)
2For projection onto k, we can use the same approach by adding nonnegativity constraint 0 in (2).
x⇤( ⇤) = argmin 0x1
( L(x) = nX
i=1
1 2 x2i + ( ⇤ yi)xi
) (4) = min n 1, [v( ⇤)]+ o , (3)
where [ · ]+ = max{ · , 0 }, both max,min in (3) are taken element-wise, and v( ) is defined as
v( ) = y 1. (4)
As stated in the introduction, under a sufficiently large k > 0, such x⇤ exists in the relative interior of =k , which implies Slater’s condition is true and thereby strong duality holds, i.e., we have the first equality sign in (3). Furthermore, we can swap the order of the min-max in (2) since it has a saddle point. By swapping the order of min-max, and using (3), we convert Problem (2) to a scalar minimization problem on :
min 2R :
⇢ !( ) := L ⇣ min
1, [v( )]+
, ⌘ , (5)
where the negative sign comes from flipping the maximization on to minimization. Now, it is clear that if we solve (5), we can construct the solution x⇤ as in (3) to solve Problem (1). Many papers [4, 5, 17] in the literature follow such a direction on solving the projection problem, in which they solve Problem (5) using a sorting-based technique to find the root of the piecewise-linear function !0( ) = 0 (see Theorem 1 for the expression of !0). In this paper, we solve (5) by Newton’s method. It is important to note that we can solve the same problem using gradient descent, but we found that it has a slow convergence, because !( ) can have a huge smoothness constant (see Theorem 1) and thus the feasible stepsize for the gradient step such that the function value will decrease is very small, which contributes to a slow convergence. We also observed that empirically the proposed method in this paper (see Algorithm 1) yields a much faster convergence (about 6-8 times faster) than the sorting-based methods, especially when the input vector y has a huge dimension; see § 3 for the experimental results. Before we proceed, we illustrate the idea of solving the projection on a toy example, see Fig.1.

heading: 2.2 The derivatives and the huge smoothness constant
text: Now we talk about the derivatives of !( ) with respect to (w.r.t.) . First, we show that !(·) is continuously differentiable and we derive its 1st and 2nd order derivatives. Then we show that the Lipschitz constant L of !0(·) is n, which is the dimension of the input vector y 2 Rn. With a huge L, we observed that gradient descent has a slow convergence on minimizing !( ), hence we propose to use the 2nd-order method, i.e., Newton’s method, to solve (2); see § 2.3. We now present the main theorem in this paper, which is about the derivatives of !( ). Theorem 1. The function !( ) is convex and twice differentiable with
!0( ) = k min ⇣ 1, [v( )]+ ⌘> 1 is n-Lipschitz, and !00( ) = nX
i=1
I0<vi<1 0, (6)
where I is an indicator function that IA = 1 if A is true and IA = 0 otherwise.
Now we present two lemmas to prove the Theorem 1. For the first lemma, recall that for all a, b 2 R,
max{a, b} = 1
2
⇣ a+ b+ |a b| ⌘ and min{a, b} = 1
2
⇣ a+ b |a b| ⌘ . (?)
Lemma 1. 8x,y 2 Rn, we have min
1, [x]+
min 1, [y]+ 2  kx yk2.
Proof. By triangle inequality and (?), for all a, b 2 R, we have min(a, 1) min(b, 1) < a b
and [a]+ [b]+ < a b . Therefore,
min(1, [x]+) min(1, [y]+) 2  [x]+ [y]+ 2  kx yk2.
The second lemma is related to the notion of optimal value function [3, Section 4]. Lemma 2. [Theorem 4.1 [3]] Let X be a metric space and U be a normed space. Suppose that for all x 2 X , the function L(x, ·) is differentiable and that L(x, ) and D L(x, ) (the partial derivative of L(x, ) w.r.t. ) are continuous on X ⇥ U . Let be a compact subset of X . Then the optimal value function ( ) := inf
x2 L(x, ), is (Hadamard) directionally differentiable. In addition,
if 8 2 U , L(·, ) has a unique minimizer x( ) over , then ( ) is differentiable at and the slope of ( ) w.r.t. is given by 0( ) = D L(x( ), ).
Now we are at the position to prove Theorem 1.
Proof of Theorem 1. To prove the differentiability of !( ), we apply Lemma 2 with X = Rn, U = R and = {x 2 X : 0  x  1}. Immediately we have that: i. L(x, ·) is differentiable; ii. L(x, ) and D L(x, ) = x>1 k are continuous on X ⇥ U ; iii. is a compact subset of X ; and iv. 8 2 U , L(x, ) has a unique minimizer x( ) = min
1, [y 1]+ over . The last one follows
that rxL(x, ) = 0 has a unique solution and further indicates that x( ) = min 1, [y 1]+ =
argmin x2 L(x, ). Now applying Lemma 2 gives ( ) = inf x2 L(x, ) = L min(1, [y 1]+), is differentiable, and the slope is
0( ) = min 1, [y 1]+ > 1 k = !0( ). (7)
We now show the derivative on !0( ) w.r.t. is the number of element of the vector v = y 1 that is strictly between 0 and 1. Let S0, S1, S be the sets of indices i that vi  0, vi 1, and 0 < vi < 1, resp.. Now, expressing 0 in (7) using the sets S0, S1, S gives
0( ) =
X
i2S0
0 + X
i2S1
1 + X
i2S (yi )
! k = |S1|+ X
i2S yi
! |S| k, (8)
which gives 0( ) = |S| +constant, hence, !00( ) = 00( ) = ( |S|) which gives (6). It is clear that L(x, ) in (2) is convex in x and the constraint set is closed and convex. For (y, k) that the feasible set for x is nonempty, by von Neumann Lemma [15], there exists a saddle point for Problem (2), as a result, ( ) is concave and !( ) = ( ) is convex. Lastly, we show the Lipschitz constant of !0( ) is n: let v( ) = y 1, then
!0( 1) !0( 2) =
min ⇣ 1, ⇥ v( 1) ⇤ + ⌘> 1 min ⇣ 1, ⇥ v( 2) ⇤ + ⌘> 1 2
 k1k2 · min 1, [v( 1)]+
min 1, [v( 2)]+
2
Lemma 1  p nkv( 1) v( 2)k2 = p nk 11 21k2 = n| 1 2|.

heading: 2.3 Solving the projection by Newton’s method
text: Now it is clear that we solve Problem (1) by solving the scalar minimization problem (5) using Newton’s method, see Algorithm 1. Before we discuss the case of projection onto k, we first give some remarks about Algorithm 1.
Algorithm 1: Newton’s method for solving Problem (1) onto =k 1 Input: a nonzero vector y, and initialize 0 2 (max(yi) 1,max(yi)) ; 2 for t = 1, 2, ... do 3 v = y t 11; 4 t = t 1 !0( t 1)
!00( t 1) , for !0( t 1) = k min 1, [v]+ > 1, !00( t 1) =
nX
i=1
I0<vi<1;
5 end 6 return x⇤ = min(1, [v]+)
Quantify the performance of the algorithm: feasibility gap As the goal of solving Problem (1) is to find a vector inside the k-crapped simplex =k , hence we can define the feasibility gap as
dist(x, =k ) 2 := dist
⇣ x, u 2 Rn | u>1 = k ⌘2 + dist ⇣ x, u 2 Rn | 0  u  1 ⌘2 .
By the definition of the algorithm, x⇤ is always inside the interval [0,1], so we can simplify the feasibility gap as dist x, u 2 Rn | u>1 = k = |x>1 k|. For the toy example in Fig.1, we can see that such a gap is monotonically decreasing for the proposed algorithm.
Range for ⇤ and !00( ⇤) for general problem In general, the minimizer ⇤ can take any value in [ 1,+1] (with 1 included) for any (y, k). Also, by definition, !00( ⇤) can take any value in [0, n]. Hence, cases like !00( ⇤) = 0 and/or ⇤ 2 {±1} are possible. For example, suppose k = n, then for any bounded y 2 Rn, one can find that a root of !0( ) = !00( ) = 0 is = 1. Furthermore, if the parameter k is too large (or too small) such that l1 and l2 in Fig.1 do not intersect each other, it is possible for the algorithm to produce ⇤ with !00( ⇤) = 0 and/or ⇤ 2 {±1}.
On k, and the feasible range of In practice, k is a user-defined value that is related to the sparsity of the desired solution vector, so k is usually (much) smaller than n, therefore we focus on the case k ⌧ n (in § 4, we will illustrate the experimental result on solving bioinformatics problems on real-world datasets, where we used k ⌧ n). In the following, we discuss , assuming 0 < k ⌧ n is selected such that the Problem (1) is feasible to solve. Now we use the toy example in Fig.1 as an illustration. As Problem (1) is feasible, we have l1 and l2 in Fig.1 intersect at a ⇤ /2 {±1} (in the toy example, ⇤ is unique). Now for the update of in Algorithm 1 to work, we need !00( ) 6= 0. Notice that !00( ) = 0 () {vj  0 or vj 1} 8j, i.e., the “bad region” such that !00( ) = 0 is the intersection of all Ij := {  yj 1 or yj  }. Therefore, the feasible range for such that !00( ) 6= 0 is the complement of such bad region, i.e,
2 R\ \
j
Ij ⇢
i min
yi 1,max yi h , (9)
where ]a, b[ denotes the open interval a < x < b, and min(yi),max(yi) are the smallest and largest value in y, resp.. We can make use of this to design an initialization of 0. For example, we can simply set 0 2 [min(yi) 1 + ,min(yi) ] for a small constant > 0, which has the cost of O(1). Or, we can search 0 within the interval (9), with a worst-case complexity of O(n log n), by sorting all the n elements of yi and trying all the n feasible intervals. Refer to the toy example in Fig.1, we see that once 0 is initialized within the feasible range, the algorithm will produce a sequence { t}t2N that converges to the minimizer ⇤. The convergence is guaranteed based on the convergence theory of Newton’s method. Furthermore, if ⇤ 6= 1, the sequence { t}t2N will not diverges to infinity since all !00( ) 6= 0, based on the following corollary. Corollary 2. If ⇤ 6= 1 is the minimizer of Problem (5), then !00( ⇤) 6= 0.
Proof. As !0( ⇤) = 0, then ⇤ (7)= k |S1|+ P
i2S yi |S| 6= 0 and hence |S| 6= 0 and ! 00( ⇤) 6= 0.
Computational cost The cost of the Algorithm 1 mainly comes from line 4, which has a periteration cost of O(n). Assuming it takes the algorithm T iterations to converge, then the total cost inside the loop is TO(n). Theoretically speaking, it can take Newton’s method many iterations to converge, but empirically we observed that the algorithm usually converges in a small number of iterations, so T is small, see Table 2 in § 3. In comparison, it costs O(n2) for the method proposed in [17] for solving the root of the piecewise-linear equation !0( ) = 0.
Comparisons with existing methods To the best of our knowledge, all existing works solve the projection problem (or a similar type of projection problem) by finding the root of !0( ) = 0 [4, 5, 17]3, in which they solve such a piecewise-linear equation by sorting-based techniques. As these methods do solve !0( ) = 0 exactly, so when the projection problem is feasible, they will eventually produce an exact result. In contrast, the proposed method is an iterative numerical algorithm, so it is possible that, due to numerical issues, the proposed method does not produce an exact solution but only an approximate solution. However, we argue that, for application with an input vector that has a huge dimension, the proposed method is good enough. If the equation !0( ) = 0 turns out to have no root for the input pair (y, k), in this case, all the sorting-based methods will fail, but the proposed algorithm is still able to produce an approximate solution, which can be useful in practice. Lastly, as the proposed method has a lower computational cost, the method is very suitable for the projection of high dimensional vectors, see § 3. The proposed method also finds application on large-scale bioinformatics problems, see § 4.
The case of projection onto k If we want to project onto k, the same algorithm can be used with the following twists: i. The feasibility gap becomes
dist ⇣ x, u 2 Rn | u>1  k ⌘ = ⇢ 0 if x>1  k x>1 k else , and
ii. has to be nonnegative, which can be done simply as = [ ]+. Note that now it is possible for !00( ) = 0, since the safe range (9) may become empty after = [ ]+, for some input vector y.

heading: 3 Experiment on comparing the projection algorithms
text: Now we show the experimental results to demonstrate the efficiency of the proposed projection algorithm. The algorithm is implemented in MATLAB4, and is compared with the MATLAB command quadprog (a MATLAB general-purpose solver for the quadratic program), the C++ implementation of [17], and a state-of-the-art commercial solver Gurobi [7] that uses the interior-point method [9] to solve the projection problem. All experiments were conducted on a PC with an Intel Xeon CPU (3.7GHz) and 32GB memory.
The result of Table 1: fast convergence of the proposed method across data size We compared the methods on the simulation data used in [17]. We generate vector y 2 Rn with yi uniformly sampled from [ 0.5, 0.5]. We pick k as an integer between 1 and n by random, and we pick the dimension n 2 {50, 102, 103, 104, 105, 106, 107, 108}. We repeat the experiments 100 times, and the comparison results are shown in Table 1. We see that as n increases, the efficiency of the proposed algorithm becomes more and more significant. For large n, many methods take more than 60 seconds (denoted as ’-’ in Table 1). The results confirmed that the proposed projection method is the most efficient one among these methods. In the following experiments, we only compare the proposed algorithm and the commercial solver Gurobi, as they are the two most efficient methods.
The result of Table 2: it only takes a small number of iterations for the algorithm to converge We currently do not have a solid theory on characterizing the exact convergent rate for Algorithm 1, hence we performed experiments to showcase the relationship between T (the number of iterations for Algorithm 1 to converge) and n (the number of dimension of the input vector y). From Table 2, we see an expected result that runtime and T increase with n, and the table empirically supports the claim that it only takes a small number of iterations for the proposed method to converge. For example, for n = 106, the empirical complexity of the proposed approach is about 10O(n), which is much lower than O(n log n) and O(n2) with such a large n.
The result of Fig. 2: fast convergence of the proposed method across range We compare the performance of the methods on solving the projection problem w.r.t. different range of the elements of y. Here, we generate y 2 Rn uniformly from [ ↵,+↵] with n = 105 and ↵ = {1, 10, 100, 1000}. As shown in Fig. 2, the proposed algorithm converges much faster than the commercial solver. The figure also shows that the performance of the proposed method is not sensitive to ↵.
3Note that [4, 5] do not tackle Problem (1) as they only solve the projection onto the unit simplex. 4The code is in the appendix.

heading: 4 Application of the projection in sparse regression
text: An application of the proposed projection method is solving the cardinality-constrained `22-regularized regression [8]:
min kwk0k
( F (w) := nX
i=1
(yi w >xi) 2 + 1
2 ⇢kwk22
) , (P ⇤)
where (xi, yi) 2 Rn ⇥ R m i=1
is a collection of m samples of data point with label. The sparsityinducing constraint kwk0  k enforces w 2 Rn to be at most k-sparse (with k non-zero elements). A recent work [10] studied the Boolean relaxation of P ⇤, which gives the following problem
min u2 k
( G(u) := y> ✓ 1
⇢ XD(u)X> + I
◆ 1 y ) , (PBR)
where u 2 k means u has to be inside the k-capped simplex. Both the original work of [10] and several follow-up papers [1, 2] showed that the solution of the Boolean relaxation PBR empirically outperforms the solution of other sparse estimation methods on recovering the sparse features, such as Lasso [14] and elastic net [18], especially when the sample size is small and the feature dimension is huge. These empirical results inspired us to develop an efficient and scalable algorithm for solving PBR, which in turn motivates the study of solving Problem (1). Now we show extensive experimental results to demonstrate the efficiency of the proposed method (Algorithm 1) on accelerating the Projected Quasi-Newton (PQN) method [11] on solving PBR on big
data. To be specific, we use both a simulation dataset and a real-world dataset with various sizes. We compare the method with a state-of-the-art method, which uses a subgradient method [2] (denoted as SS) implemented in Julia to solve a min-max problem that is equivalent to PBR. In addition, we compare with the classical elastic net (ENet) [18] as a reference. For a fair comparison between methods that are implemented in different programming languages, we report the computational time for each algorithm relative to the time needed to compute a Lasso estimator with Glmnet [6, 13] on the same programming language, and on the same dataset. Glmnet is a package that fits the data to a generalized linear model, based on a penalized maximum likelihood estimation method. Glmnet has been implemented in various programming languages, such as Julia, R, and MATLAB.
Sparse regression on simulation data First, we conduct experiments on simulation data generated as follows: let [m] = {1, 2, . . . ,m}, we draw xi ⇠ N (0n,⌃), i 2 [m] independently from a n-dimensional normal distribution with zero mean and co-variance matrix ⌃. We randomly sample a weight vector wtrue 2 { 1, 0,+1} with exactly ktrue non-zero coefficients. Then we generate noise vector ✏ 2 Rm, where ✏i, i 2 [m] are drawn independently from a normal distribution scaled according to a chosen signal-to-noise ratio (SNR) defined as p SNR = kXwtruek/k✏k. Finally, we form Y = Xwtrue + ✏, where X = [x1,x2, . . . ,xn]>. We evaluate the performance on simulation data using accuracy defined as Acc(w) :=
{ i : wi 6= 0, wtruei 6= 0 } / { i : wtruei 6= 0 }
. For the first simulation data, we generate 100 datasets by setting ⌃ij = p|i j| to a Toeplitz covariance matrix with p = 0.2, n = 1, 000, ktrue = 20, and SNR= 6. As shown in Fig. 3 (a), PQN + the proposed projection and PQN + Gurobi have exactly the same performance in in terms of Acc value, which is expected since the only difference between these methods is the projection algorithm. In addition, the two methods converge quickly to Acc= 1, outperforming all the other methods. Fig. 3 (b) shows that PQN + the proposed projection is the most efficient algorithm for solving PBR among all the tested methods. The second simulation data is generated similarly to the first one. The only difference here is that p = 0.7, indicting the higher correlation between the features, meaning that now it is harder to recover the true support of wtrue. As shown in Fig. 3 (c), PQN + the proposed projection and PQN + Gurobi again have the same performance, and they still outperform all the other methods. Again, PQN + the proposed projection is the most efficient algorithm for solving the problem here. For the results on comparing the convergence between PQN and SS, as well as the details of all the parameter settings for all the methods, see the appendix.
Sparse regression on GWAS data We conducted experiments to evaluate the efficiency of the proposed projection algorithm on a Genome-Wide Association Study (GWAS) dataset. The data was obtained from the dbGaP web site, under the dbGaP Study Accession phs000095.v3.p15. This study is part of the Gene Environment Association Studies initiative with the goal to identify novel genetic factors that contribute to dental caries through large-scale genome-wide association studies of well-characterized families and individuals at multiple sites in the United States. The data we used in this experiment contains 227 Caucasians with 11,626,696 Single-nucleotide polymorphisms
5https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000095. v3.p1
(SNPs) for all chromosomes. Here, the goal is to identify the important SNPs by performing a sparse regression. We solve the Boolean-relaxed sparse regression problem PBR on chromosome 20, 5, and 2, which have 252,257 SNPs, 750,860 SNPs, and 1,562,480 SNPs, resp.. Fig. 4 shows the comparison result between PQN + the projection and PQN + Gurobi. We find that PQN + the proposed projection outperforms in all three problems (with different numbers of SNPs) in terms of computational time. In addition, PQN + the proposed projection and PQN + Gurobi converge to the same objective function value, indicating the correctness of the proposed projection algorithm.
Results on comparing different methods on the simulated GWAS data To further demonstrate the practical impact of solving PBR using the proposed algorithm, here we further compare different methods using the simulated GWAS data described in [16]. We compared PBR with Lasso [14] and Elastic Net (ENet). We used the out-of-sample Mean Square Error (MSE) to select the parameter k for PBR. For Lasso and ENet, we used cross-validation with the minimum MES+1SE (Standard Error) to tune the hyper-parameters, which was suggested by [16]. There are three different simulated GWAS datasets as stated in [16], in which they have different levels of correlation (LD: High, Mixed, and Low) between SNPs. All the datasets have 50,000 SNPs and 25 significant SNPs. In Table 3, we showed the comparison when the sample size is 1000 and 100. Here we report the median and the standard deviation (in parentheses) over 100 replicates. We can see that PBR can find more correct SNPs with low false positive rate, especially when the sample size is only 100. For the comparisons when the sample size is 150 and 50, see the appendix.

heading: 5 Conclusion
text: We propose to use Newton’s method to solve the problem of projecting a vector onto the k-capped simplex. On the theory side, we transform the problem to a scalar minimization problem and provide a theorem to characterize the derivatives of the cost function of such a problem. On the practical side, we show that the proposed method outperforms existing methods on solving the projection problem, and show that it can be used to accelerate the method on solving sparse regressions in bioinformatics on a huge dataset.

heading: 6 Acknowledgments
text: The presented materials are based upon the research supported by the Indiana University’s Precision Health Initiative. Andersen Ang acknowledges the support by a grant from NSERC (Natural Sciences and Engineering Research Council) of Canada.

Review summaries:
This paper poses a new method for projecting vectors onto the k-capped simplex. This method reformulates the quadratic program into a form that can harness Newton method, a second-order iterative optimization algorithm, to solve the projection goal. The key novelty is using a second-order method to iteratively solve for a minimum, as opposed to sorting-based methods with a cubic complexity. Furthermore, first-order methods are shown to be too slow based on a proof that the Lipschitz constant scales with dimensionality. Compared to existing methods for k-capped simplex projection, the proposed method is several times faster, and in synthetic experiments achieves state-of-the-art accuracy. On experiments in high-dimensional bioinformatics, the proposed projection converges faster and equally well as other state-of-the-art methods.

In this paper the authors propose a method for projecting onto the capped simplex. The algorithm lies on a well-known framework, which transforms the optimization problem via Langrange and duality theory, in other equivalent problem. The contribution is on a different approach, using the Newton method, to solve one of the sub-problems. The algorithm is tested on a problem with genomic data. 

This paper proposes an algorithm based on Newton’s method for projecting a vector with bounded elements onto a hyper-cube cut by a hyperplane (also known as the k-capped simplex). The time complexity of the algorithm is roughly linear (in terms of the vector dimension) and theoretical insights for partial justification of the method are provided. Empirical evaluation demonstrates that the algorithm (1) finds a solution with high precision on large-scale datasets and (2) considerably outperforms state-of-the-art sorting-based methods in terms of run time. Finally, the proposed algorithm is also applied to solving a sparse regression problem from the bioinformatics domain on which it is able to accelerate the Projected Quasi-Newton (PQN) method from 3 to 6 times compared to its alternatives.

This paper addresses the problem of accelerating the projection of a vector onto the k-capped simplex. From the theoretical side, the authors showed that the given problem can be casted as a scalar minimization problem, which is tackled by a Newton’s method, unlike existing methods based on sorting algorithm. In computational experiments, the authors evaluated the proposed method in a sparse regression problem. The baselines are Matlab’s quadprog, reference[15], and Gurobi, a commercial interior point solver. The proposed method worked faster especially when the response vector y is high-dimensional. It is also shown to have a higher support recovery nature than that of Elastic net. In experiments using real dataset consisting up to 1.5 million SNPs, the proposed method worked much faster than Gurobi.