Fast Projection onto the Capped Simplex with Applications to Sparse Regression in Bioinformatics
Authors: Andersen Ang
Jianzhu Ma
Nianjun Liu
Kun Huang
Yijie Wang
Venue: NeurIPS 2021 Poster (2021)
Keywords: Capped Simplex
Projection
Newton's Method
Sparse Regression via Boolean Relaxation
Bioinformatics
GWAS

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.

Fast Projection onto the Capped Simplex with Applications to Sparse Regression in Bioinformatics
Andersen Ang"
Jianzhu Ma Institute for Artificial Intelligence Peking University maj anzhu@pku. edu. cn
Dept. of Combinatorics and Optimization University of Waterloo
maSangduwaterloo.cr
Nianjun Liu
Kun Huang
Dept. of Epidemiology and Biostatistics
Dept. of Biostatistics and Health Data Science Indiana University
Indiana University Bloomington
liun an@indians edu
kunhuang@iu.co edu
Vijie Wang"
Dept. of Computer Science Indiana University Bloomington yi vijwangiiu. edu
Abstract
We consider the problem of projecting vector onto the so-called k-capped simplex, which is hyper-cube by a hyperplane. For n-dimensional input vector with bounded elements, we found that simple algorithm based on Newton's 's method is able solve the projection problem to high precision with complexity roughly about O(n), which has a much lower computational cost compared with the existing sorting-basee based methods proposed in the literature. We provide 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 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 I million variables more).
We further illustrate the effectiveness of the proposed algorithm solving sparse regression 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-an methods.
Introduction
The k-capped simplex is defined as A (xE x 0 <x an input parameter, 0 and are vectors of zeros and vector of ones in R", respectively (resp.), and the sign < is taken element-wise. Geometrically, for some > 0. the set A represents an n n-dimensional hyper-cube [0, 1]" with corner chopped off by the half-space x < k. If the in Aş is replaced by the hyper-plane constraint x T = k, for some A >> 0, this set represents the slice of the hyper-cube cut by the hyper-plane, and denote this set as A In this paper, we consider the problem of projecting a vector y onto the k-capped simplex:
1}. where & E R
half-space constraint x
= proj (y) = argmin { ||x - yl1 s.t. x E A (or A
(1)
Both authors contributed equally.
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

Remark 1. (On the input k) In this paper, we assume 0 < l < < 3C is properly selected such that the set D (and A has nonempty interior, and therefore the projection is feasible. Le., we assume E TR p k) nig E R. 10-0001 + Ø (and similarly for ST). In this setting, unique global minimizer x" of Problem (T) exists, since the Sunction ||x
xx--lio y strongly convex. The value of k is closely related to the root of piecewise-lineal equation that plays a critical role in the solution of Problem
① see discussion in $ $2
Problem I 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 $
To the best of our knowledge, for an input vector E 12". the state-of-the-art algorithm solving the Projection (T) has a time complexity of (n²)). We found that a simple heuristic algorithm based on Newton's 's method is able to solve the same problem to high precision with a complexity roughly (n) on average, based on exploiting the 2nd-order information of a scalar minimization
for
problem on the Lagrangian multiplier of Problem 1 see S 2.3
Contributions We proposed a simple algorithm based on Newton's 's method, an alternative to the sorting-based method to solve Problem ① We first convert the problem into a scalar minimization problem, then we show that such scalar problem can be solved using Newton's 's iteration, see Theo- rem | We then show that the proposed method is superior to the current state-of-the. the-the-art method in terms of runtime. For vector with dimension 10,, 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 solving sparse regression in bioinformatics problem on real-world dataset with one million data points, see $
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 special case. To be specific, we found that, when the input vector y has bounded elements (i.e., lyil S a, for a 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 a 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 example,
bounded. Using bioinformatics runtime on identifying the genetic factors that contribute to dental caries (i.e., tooth decay) from half-hour to just 3 minutes; see $
ordinary PC, we are able to reduce the
Paper organization In $ we give the background material, show the main theorem, and discuss 3 we show experimental results to verify the efficiency of the proposed
the proposed algorithm. In S algorithm on solving Problem [ In $ $ 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.
2 Projection onto the capped simplex
In this section, we discuss how solve Problem D.. By introducing Lagrangian multiplier y, we convert (1) to a scalar problem in the form of min, (7)), see ( Then, we provide theorem for characterizing the function w(y). More specifically, we show that ((7) is n-smooth, i.e., the Lipschitz constant of the derivative w' (7) grows linearly as the dimension
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 's method, to solve Problem ( We end this section by discussing the properties of proposed Newton's method.
2.1 Reformulate the projection problem to scalar minimization problem
Consider the case of projection onto A= Problem ① is equivalent to Problem 2 : by introducing
a Lagrangian multiplier 7 that corresponds to the equality constraint, we form a partial Lagrangian
L(x, 7), and arrive at the min-max problem 2 where we keep the inequality constraint on x as
min max {{ L(x, 7) = マ ||x y122+ kb) OSx51 TER
(2)
For projection onto Ac. we use the same approach by adding nonnegativity constraint 0
2

---

w'()] = xx(x))
accoss
w76)
(((n))
(R x(x))
4 k
x(7)
70
h ::
x(n)
1
*
7'
-0.5
w'(n)
-2
-2
-2
-1.5
2
A x
Iteration
Inenation
S on a toy example with y [0.1, 1.5, -1] [1] TRa, 1.5 and too 1.1. The line l2 intersects the piecewise-liner function 4 at y", which is the root of l2- The idea of the corting-based methods is to find such root by sorting the 3 elements of y and try each of the linear segment of I find the root. The proposed method (Algorithm [ solves the problem in 2 iterations, hence it gains speed-up. In the middle we plot and the iteration of 7/- The right-most plots show the iteration of w' and w" (see Theorem [ for their explicit expression). Note that 11, l2 may not intersect each other for some (y,k). in this case l2 has no root and the sorting-based methods do not work. However, the proposed algorithm can still produce an approximate solution.
Figure 1: Illustration on solving Problem
Closed-form solution on x
Let be the optimal yy in Q then the minimizer x* of ( is
x ()') = argmin [(x) = Σ 2 + ( -y:0 1=1
min
(3)
0<x$1
where [+]+ max{ 0, J. both max, min in 3 are taken element- and v(n) is defined as
v()) y
As stated in the introduction, under sufficiently large ls > O, such x" exists in the relative interior of Ag. which implies Slater's condition is true and thereby strong duality holds, i.e., we have the first equality sign in Bp. Furthermore, we can swap the order of the min-max in 2 since it has
saddle point. By swapping the order of min-max, and using we convert Problem 2 to a sealar minimization problem on y:
(u())=-((mi(()))))
min yER
(5)
where the negative sign comes from flipping the maximization on to minimization. Now, clear that if we solve ( we can construct the solution x" as in ( to solve Problem ( Many papers
4. [7] in the literature follow such a direction on solving the projection problem, in which they solve Problem 5 using sorting- based technique to find the root of the piecewise-ling linear function w'()) = 0 (see Theorem for the expression of w').
In this paper, we solve S by Newton's 's method. it important to note that we can solve the same problem using gradient descent, but we found that it has slow convergence, because w(>) can have a huge smoothness constant (see Theorem D and thus the feasible stepsize for the gradient step such that the function value will decrease is very small, which contributes to slow convergence. We also observed that empirically the proposed method in this paper (see Algorithm| [ 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 S for the experimental results.
Before we proceed, we illustrate the idea of solving the projection on a toy example, see Fig
2.2 The derivatives and the huge smoothness constant
Now we talk about the derivatives of w()) with respect (w.r.t) 7. First, we show that w(-) is continuously differentiable and we derive its 1st and 2nd order derivatives. Then we show that the Lipschitz constant of w'(-) is n, which is the dimension of the input vector y E R", With huge L. we observed that gradient descent has a slow convergence on minimizing w(7), hence we propose to use the 2nd-order method, i.c.. Newton's 's method, to solve @ see $ 22.3
We now present the main theorem in this paper, which is about the derivatives of w()).
Theorem 1. The function w(T) is convex and twice differentiable with
w'(7) = k min (1.Vv)))))' is n-piischir,, and w²())= = Σ Toenet 20, (6)
where 1 is an indicator function that LA = 1 if A is true and IA 0 otherwise.
3

---

Now we present two lemmas to prove the Theorem ! For the first lemma, recall that for all a, b E R.
max(a,b) b) (a bb+b+aa a -- bl)
and min b} 2 a a b- bb-l- blla-abb)
1
min (1, [x]+) min (1, (y]+))| S ||x yl12.
Lemma 1. Vx, E R", have
Proof. By triangle inequality and (+), for all a, 6 E 1. we have min((, 1) min(b, 1) < and |[a]+ [6]+] < aa b) Therefore,
min [x]+) min 1x)+)112 - ||[x]+ - 1x1+112 S |xx yl12.
The second lemma related to the notion of optimal value function
K Section 4].
Lemma 2. [Theorem 4.1 [3] Let X be a metric space and U be a normed space. Suppose that for all E X. the function [(x, is differentiable and that L(x, 7) and D,L(xx 7) (the partial derivative of [(x, 7) w.r.t. 7) are contimuous on X x xU. Let & be a compact subset of X. Then the optimal value function (6T)) inf
7(x, is (Hadamard) directionally differentiable. In addition, x
v/ E U, ((:-7) has a unique minimizer x(y) over $ then ((T)) is differentiable at 7 and the slope of ((7)) w.r.t fy is given by ('()) (x,)( 7).
Now are at the position to prove Theorem L
Proof of Theorem 7 To prove the differentiabilitty of (()), we apply Lemma 2 with X'
R T
and 4- {x X' 0 < x < 1}. Immediately we have that: i. [(x. -) is differentiable; ii [(x, 7) and D, 7) = x
& are continuous on X' x 24; iii. Tb is a compact subset of X:; and iv. by E U4 [(x, 7) has a unique minimizer x(7) = min (1. ly - -11]+ over t. The last one follows
that Vx[(x,7) 7) = O has a unique solution and further indicates that x(7) = min (1, - [--->1]+) = argmin L(x,7), Now applying Lemma 2 gives (²(1) = inf ((x,7) L( minn(, [y -11)))+) is XE# differentiable, and the slope is
x
('(r) = min (1, yy - - - w7))
We now show the derivative on w'(7) w.r.t. is the number of element of the vector v = =y y-- that strictly between 0 and Let So. S1, be the sets of indices that v; S 0, v 2 1, and 8/2 < 1, resp.. Now, expressing ø' in using the sets So, S1, S gives
('(7) = Σ
0+ Σ1
1
Σ (n -6--7)
7)
- k = =1811 +
Σ iES
y1
18 k
(8)
SE
SES,
SES
which gives d'(7) = - hence, (7) = (-) = (-1S)) which gives 6
It is clear that L(x, 7) in 2 is convex in x and the constraint set closed and convex. For (y k) that the feasible set for x is nonempty, by von Neumann Lemma [5). there exists a saddle point for
Problem 2 as result, o()) is concave and w(7) -6(y) is convex.
Lastly, we show the Lipschitz constant of w/ (1) is n: let v(r) = y 71, then
((m))-((xx)) = min 1. [v())))
min (1. [V((T>)) 1//22 ||||| min (1,[v(7))++ - min (1, [v(22))++
<
Lemma IF
vnllv()) v(12)l22 = -all - 722 1ellle l|2 = n//1 72l-
2.3 Solving the projection by Newton's method
Now it is clear that we solve Problem | by solving the scalar minimization problem 5 using Newton's 's method, see Algorithm Before we discuss the case of projection onto As, first give some remarks about Algorithm I

---

Algorithm 1: Newton's method for solving Problem onto A
Input: a nonzero vector y, and initialize "Yo E (max(ys) 1, max (yi)): for t 1. do
1 w' (7t-1) for w' (7-1) min (1,1+]+) 1, w w((n-1 = Σ
m
7t 7
wn (m)-))
Jo<e
end
return x" = min(1, [v]+)
Quantify the performance of the algorithm: feasibility gap As the goal of solving Problem is to find a vector inside the k-crapped simplex Ag hence we can define the feasibility gap as
dist(x, ) = dist fu R u = k))² + dist xx, fu EuER! R 0 (0Su<i)]
|
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 R. u = k)) |x k||. For the toy example in Fig we can see that such gap monotonically decreasing for the proposed algorithm.
Range for and w"("") for general problem In general, the minimizer >" can take any value in (-00, +oo] (with oo: included) for any (y, k). Also, by definition, wn (7") can take any value in [0, n]]- Hence, cases like (>*) = and/or " E { are possible. For example, suppose k = n, then for any bounded E IR" one can find that root of wn)) = w (n) 0 is y = - co. Furthermore,
the parameter k is too large (or too small) such that 4 and l2 in Fig D do not intersect each other, is possible for the algorithm to produce 7 with (n"()) = 0 and/or r {=oo).
On kt and the feasible range of In practice, A is 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 solving bioinformatics problems real-world datasets, where we used k < n). In the following, we discuss y, assuming 0 < A < n is
selected such that the Problem ①1 is feasible to solve. Now we use the toy example in Fig 1 as illustration. As Problem
① is feasible, we have I and 12 in I lintersect at a 7" & {=:00) } (in the
toy example, 7 is unique). Now for the update of y in Algorithm I to work, we need w" ()) + 0.
Notice that (n)) = > {uz < 0 or Uj 2 1} Vj. i.e., the "bad region" such that (>) ) = 0 is
the intersection of all I; {o y/ 1 or yj S 7. Therefore, the feasible range for 7 such that w" (7) * 0 is the complement of such bad region, i.e,
7 E TERR c min (v:) - 11, max (vi)
where Ja, bf denotes the open interval a < x < b, and min(y:). max(y:) are the smallest and largest value in y, resp.. We make use of this design initialization of to- For example, we can simply set to E [min(y;) - + ó. min(y: 5) for small constant 5 >> 0, which has the cost of (1). Or, can search tt within the interval ( with worst-case complexity of O(n log n), by sorting all the n elements of y, and trying all the n feasible intervals. Refer to the toy example in Fig T see that once To initialized within the feasible range, the algorithm will produce sequence {Tr hen that converges to the minimizer y* The convergence is guaranteed based on the convergence theory of Newton's method. Furthermore, 7 j xo, the sequence {7t (ee will not diverges to infinity since all w" (x) + O, based on the following corollary.
Corollary 2. 7f + do is the minimizer of Problem then w"(r") #0.
Proof. As w (r*) O. then 7 A= 511+ M6 [8]
# and hence |S| o and w w"(?") +
Computational cost The cost of the Algorithm I mainly comes from line 4, which has a per- iteration cost of (n). Assuming it takes the algorithm T iterations to converge, then the total cost inside the loop is (n)).. Theoretically speaking. 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 small, see Table 2 in 3 In comparison, it costs (n?)) for the method proposed in He for solving the root of the piecewise- linear equation w/ (7) 0.
5

---

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 w/ u/(?) in which they solve such a piecewise- methods do solve w w'(r) = 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 approximate solution. However, we argue that, for application with an input vector that has huge dimension, the proposed method is good enough.
equation by sorting-based techniques. As these
If the equation w' (7) 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, the proposed method has lower computational cost, the method is very suitable for the The proposed method also finds application on
projection of high dimensional vectors, see $ 3
large-scale bioinformatics problems, see $
The case of projection onto De If we want to project onto 4g. the same algorithm can be used with the following twists: i. The feasibility gap becomes
0
i' S k
dist x, (
and
- k
else
ii y has to be nonnegative, which can be done simply as 7 = hn] Note that now is possible for
w" (7) 0, since the safe range ( may become empty after y +. for some input vector y.
3 Experiment on comparing the projection algorithms
Now we show the experimental results to demonstrate the efficiency of the proposed projection algo- rithm. The algorithm is implemented in MATLAB and is compared with the MATLAB command quadprog (a MATLAB general-purposs solver for the quadratic program), the C++ implementation of [17) and offtt-o--hheart commercial solver Gurobi (7] that uses the interior-poin method [9] to solve the projection problem. All experiments were conducted on a PC with an Intel Xeon CPU (3. (3.7GHz) and 32GB memory.
The result of Table [ fast convergence of the proposed method across data size We compared
the methods on the simulation data used in [17. We generate vector y I²" with y; uniformly sampled from [-0.5, 0.5]. We pick l as an integer between 1 and n by random, and we pick the dimension {50, 102, 103, 10a 105 106, 10", 100). We repeat the experiments 100 times, and the comparison results are shown in Table J 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 D.. 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 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 [] showcase the relationship between T' (the number of iterations
hence we performed experiments for Algorithm to converge) and n (the number of dimension of the input vector y). From Table
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 10 the empirical complexity of the proposed approach is about 10O(n). which is much lower than (n log n.) and O (n²) with such large n.
The result of Fig. D 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 I uniformly from I-a, +a] with n 10 and {1, 10, 100, 1000). As shown in Fig. g 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 a.
3 Note that 5 not tackle Problem [ as they only solve the projection onto the unit simplex.
'The code is in the appendix.

---

Table 1: Runtime in seconds (in mean=ist) of different methods on various sizes n.
Method
5 x 10'
10a
10a
104
quadprog
0.0128 + 0.01058
0.0034 = 0.0095
0.14191 0.0183
C++
0.00003:1 0.00007
0.00004:: 0.00004
0.0011 + 0.0014
0.101 0.129
Gurobi
0.002 = 0.0048
0.0022: 0.0021
0.0041 0.0022
0.02+ 0.0032
Proposed
0.00013 + 0.00014
0.00016 + 0.00065
0.00039=: 0.00038
0.0039: 0.0037
Method
10"
10
10'
10
quadprog
C++ Gurobi Proposed
10.3985 13.703
0.2413++ 0.00117
2.8229 0.0955 0.4723 + 0.3629
32.26 + 0,4953 4.2067 + 2.5795
0.0406+ 0.0377
14.925 + 2.9893
Table 2: Runtime and iterations (in mean) on varying n, k (over 100 experiments).
Varying (k 10a)
Varying k (n 10%)
Time
# iter. T
Time
# iter. T
n 10
0.0011 0.00013
6.2 0.4
10
0.15 0.0096
11.4 + 0.49
10a
0.0108 + 0.00081
8.0 0.2
10a
0.13 = 0.0054
10 - 0
n 106
0.13 + 0.0054
10
A: 103
0.12 - 0.0084
0.5
n 107
2.1 + 0.03
12
104
0.11 0.0039
7
1
21.95 = 0.046
13
10S
0.10 =
0.2
Le. 100
100
130 100
123 321
20
20
20
10
10
13
30
300
-30
-20
-10
-20
-20
-20
-20
-20
-50
-30
-30
-10
0.1
02
Time (s)
Time (3)
Time 151
Time 60
Figure 2: Comparison between the proposed projection method and the Gurobi projection for 10, 100, 10003. The thick curves are the median of the results of over 100 experiments on 100 datasets. The figure shows the superior performance of the proposed method in terms of runtime.
Application of the projection in sparse regression
An application of the proposed projection method is solving the cardinality E Ge-reuularzzed regression {
1 2 pll 1
min [wlloSk
F(w)
Σ (w w -x,)2 +
(P*)
where {(x;, 9: E IR" x IR J inducing constraint ||w lo < k enforces w E R² to be at most k-sparse (with k non-zero elements). recent work [T0] studied the Boolean relaxation of pa, which gives the following problem
m is a collection of m samples of data point with label. The sparsity-
(((() (c()==y) y (<x)=>>(XXDDX +1) y
min uE
(PBR)
where E A means has to be inside the k-capped simplex. Both the original work of [TC and several follow-up papers [ showed that the solution of the Boolean relaxation PER empirically outperforms the solution of other sparse estimation methods recovering the sparse features, such as Lasso [I]] and elastic net [1]] especially when the sample size is small and the feature dimension huge. These empirical results inspired us develop efficient and scalable algorithm for solving
PAR. which in turn motivates the study of solving Problem |
Now we show extensive experimental results to demonstrate the efficiency of the proposed method (Algorithm I on accelerating the Projected Quasi-Newton (PQN) method [|] on solving PBR on big
7

---

p-02
B
p.02
(6)
(d)
Relative time wunt ginggeeeee de log)
1o"
0.3
0.6
163
0.2
0.2
40 200
100 200
100 120 140 102 180 220
80 1.0 120 142 100 150 200
306 120
100 200
Sample dite
Sample size
Sample $100
Sampl $220
-7 PON++our projection
PGN+-Gurbli
-1 witnem00
- 550teratinnn20000 --- ENet
Figure 3: Comparison between different methods. The results are the average over 100 datasets, and all the error bars are in med From left to right, the sub-figures are: (a) The Acc values between different methods when n 10 k 10, p 0.2, and SNR= 6; (b) The computation time of the algorithms in (a): (c) The Acc values between different methods when changed to 0.7; and (d) The computation time of the algorithms in (c).
both simulation dataset and real-world dataset with various sizes. We compare the method with a state-of-the-an method, which uses a subgradient method [2] (denoted as SS) implemented in Julia to solve min-max problem that is equivalent PBR. In addition, we compare with the classical elastic net (ENet) [8] as 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] [3] on the same programming language, and on the same dataset. Glmnet is a package that fits the data to generalized linear model. based on penalized maximum likelihood estimation method. Gimnet has been implemented in various programming languages, such as Julia, R, and MATLABB.
data. To be specific,
Sparse regression on simulation data First, we conduct experiments on simulation data generated as follows: let [m] {1.2.... m). we draw x, (O,, Σ), [m]] independently from n-dimensional normal distribution with zero mean and co- co-variance matrix Σ. We randomly sample weight vector w-nuee E 0,1,0,+1, with exactly kiue non-zero coefficients. Then we generate noise vector E RTm where 6, E [m] are drawn independently from a normal distribution scaled according chosen signal-to-noise ratio (SNR) defined as VSNR ||
wree ||/|||]| Finally, we form Y = Xw +,, where X = X21 Xm] We evaluate the performance on simulation data using accuracy defined as Acc( () == {i- wi F O, whte * 0 / [ wite f0]l-
For the first simulation data, generate 100 datasets by setting Σg p 5ー5l to Toeplitz covariance matrix with = 0.2, n = 1, 000, kinue 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 Acc= outperforming all the other methods. Fig. (b) shows that PQN + the proposed projection is the most efficient algorithm for solving PAR 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 winee, 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 Genome-Wide Association Study (GWAS) dataset. The data was obtained from the dbGaP web site, under the dbGaP Study Accession phs000095.v3.p This study is part of the Gene Environment Association Studies initiative with the goal identify novel genetic factors that contribute to dental caries through large-scale genome-wide association studies of well-characterried families and individuals at multiple sites in the United States. The data we used in this experiment contains 227 Caucasians with 1 11,26,,699 Singlc-nucleotide polymorphisms
htttps://ww..neee nebi. nlm, nih, www.://////oooosss//////g/gggge/ http:///c...oo.e//////////////oa wgwp::otyddddddddpppooooo
v3. p1

---

Chromosoms
Chromosome (750.860
Chromosome (1,362,480 SNP))
50
50
45
8.4 44
35
15
(30
10
28
00
20
10
10
10
50
100
200
230
200
400
600
800
1000
200
400
600
500
1000
1200
1400
1500 1800
Figure 4: Convergence comparison between PQN with the proposed projection algorithm and PQN with Gurobi projection. From left to right, the sub-figures are: (a) Comparison on chromosome 20; (b)Comparison on chromosome 5; and (c) Comparison chromosome 2.
(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 which have 252,257 SNPs, 750,860 SNPs, and 1,562,480 SNPs, resp..
Fig. f 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 PER using the proposed algorithm, here we further compare different methods using the simulated GWAS data described in [[6]
We compared PBR with Lasso [E] and Elastic Net (ENet). We used the sample Mean Square Error (MSE) to select the parameter & for PIR- For Lasso and ENet, used cross-validation with the minimum MES+ ISE (Standard Error) to tune the hyper-parameters, which was suggested by [16] There are three different simulated GWAS datasets 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 that PER 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.
Table 3: SNP=50,000, significant SNPs = 25.
Sample size
1000
Sample size = 100
Lasso
ENet
PaR
Lasso
ENet
PBR
High LD
Correct
3(0.91)
25(0)
25(0)
2(1.01)
20(3.47)
25(1.52)
False positive
0(0)
0(0.1)
0(0)
1(0.93)
4(1.4))
1(1.03)
Mixed LD
Correct
3(0.75)
18(1.47)
25(0)
2(1.35)
16(1.47)
24(2.03)
False positive
0(0)
0(1.19)
0(0)
1(0.68)
2(0.35)
2(1.24)
Low LD
Correct
17(2.12)
25(0.51)
25(0)
5(2.12)
17(3.51)
23(3.07)
False positive
0(0.39)
0(2.06)
0(0)
2(0.59)
5(2.06)
3(1.47)
5 Conclusion
We propose to use Newton's 's method to solve the problem of projecting vector onto the k-capped simplex. On the theory side, we transform the problem to a scalar minimization problem and provide 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 solving sparse regressions in bioinformatics on a huge dataset.
6 Acknowledgments
The presented materials are based upon the research supported by the Indiana University's Precision Health Initiative. Andersen Ang acknowledges the support by grant from NSERC (Natural Sciences and Engineering Research Council) of Canada.

---

References
[1] D. Bertsimas and B. V. Parys. Sparse high-dimensional regression: Exact scalable algorithms and phase transitions. The Annals of Statistics, 48(1):300- - 323,, 2020.
[2] D. Bertsimas, J. Pauphilet, and B. V. Parys. Sparse Regression: Scalable Algorithms and Empirical Performance. Statistical Science, 35(4):555 - 578, 2020.
[3] J. F. Bonnans and A. Shapiro. Optimization problems with perturbations: A guided tour. SIAM Review, 40(2):228 264, 1998.
[4] L. Condat. Fast projection onto the simplex and the II ball. Mathematical Programming, 158(1):575-585, 2016.
[5] J. Duchi, S. Shalev-Shwartz, Y. Singer, and T. Chandra. Efficient projections onto the 1-ball for learning in high dimensions. In Proceedings of the 25th international conference on Machine learning, pages 272-279, 2008
[6] Friedman, T. Hastie, and R. Tibshirani. Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software, 13(1):1-22, 2010.
[7] L. Gurobi Optimization. Gurobi optimizer reference manual, 2021.
[8] T. Hastie, R. Tibshirani, and J. Friedman. The Elements of Statistical Learning. Springer Series
Statistics. Springer New York Inc., New York, NY, USA, 2001.
[9] Y. Nesterov. Lectures on convex optimization, volume 137. Springer, 2018.
[10] M. Pilanci, M. J. Wainwright, and L. El Ghaoui. Sparse learning via boolean relaxations. Mathematical Programming. 151(1):63-87 2015.
[11] M. Schmidt, E. Berg, M. Friedlander, and K. Murphy, Optimizing costly functions with simple constraints: A limited- memory projected quasi-newton algorithm. In AISTATS 2009, pages 456-463, PMLR. 16-18 18 Apr 2009.
[12] M. Schmidt, D. Kim, and Sra. Projected Newton-type methods machine learning, pages 305-330, MIT Press, Cambridge, MA, USA, Dec. 2011.
[13] N. Simon, J. Friedman, T. Hastie, and R. Tibshirani. Regularization paths for cox's proportional hazards model via coordinate descent. Journal of Statistical Software, 39(5):1-13, 2011.
[14] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society (Series B), 58(1):867-288, 1996.
[15] Neumann. Zur theorie der gesellschaftsspiele. Mathematische Annalen, 100(1):295-32) 1928.
[16] P. Waldmann, G. Mészáros, B. Gredler, C. Fürst, and J. Sölkner. Evaluation of the lasso and the elastic net genome-wide association studies. Frontiers Genetics, 4:270, 2013.
[17] W. Wang and C. Lu. Projection onto the capped simplex. arXiv preprint arXiv:1503.0100 2015.
[18] H. Zou and T. Hastie. Regularization and variable selection via the elastic net. Journal the Royal Statistical Society, Series B, 67(2):301-320 2005.
10


Reviews:
[{"review_id": "cEdTt9ECmjL", "reviewer": "Reviewer_3a3C", "summary": "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.\n", "questions": "", "limitations": "", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The contribution of the paper is the use of Newton's method for solving the subproblem (5). The main part of the paper is focused on the computation of derivatives and results related to them, many of them very basic. I think this part is not relevant, and in any case might be included in the supplementary material. \n\nThe modification of the method seems to be of importance, given the complexity discussed (even though it might not work fast for some vectors), but I think the presentation should be dramatically improved.\n\nThe texts is repetitive, the toy example in Fig.1 is not explained or motivated at all, and there are some technical parts with details or omissions. For instance, in line 26 it should say that the vector to be projected is $y$; problem (1) should be actually argmin; on the line 30, the convexity of $\\Delta_k$ has to be used, otherwise the statement is false, Newton's algorithm doesn't have a stopping criteria, it runs for a fixed number $T$ of iterations, but $T$ is not an input, and is not discussed as it should.\n\nAlso, please proofread and check for typos, missing words, etc, and please be consistent with the notation, for instance in some parts $[]_+$ is used, and in other parts just $\\max(0,\\cdot)$, for the same argument.\n\nI tend to think that the result is important, but unfortunately I think the paper needs more work.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors claim that the bound of the input vector plays a role in the time-performance of their method, but more discussion is needed in that sense.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "Zl7Mr2fx5J0", "reviewer": "Reviewer_CRDH", "summary": "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.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Originality: The formulation of the optimization problem is well known; however, to the best of my knowledge, the idea of converting this problem into a scalar minimization problem and solving it by Newton’s iteration appears to be novel. Also, there are some theoretical insights supporting the proposed algorithm which I believe have not been presented in the literature. The related work appears to be comprehensive enough.\n\nQuality: The design and justifications for the proposed algorithm are technically sound. The observations made regarding the high precision and efficiency of the algorithm are empirically well supported. Overall, the work seems to be mature in terms of quality.\n\nClarity: The paper is well written, organized, and technically detailed to a satisfactory extent. The notation is clear and consistent throughout the paper.\n\nSignificance: The problem addressed in this paper is of considerable importance as high-impact bioinformatics applications such as the sparse regression problem of identifying single-nucleotide polymorphisms for chromosomes from the Genome-Wide Association Study (GWAS) dataset. Another significant aspect of this paper is the proposed algorithm’s superior efficiency over state-of-the-art methods, while being even 6-8 times more efficient than commercial software on highly dimensional data with million variables or more.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "One major concern is that it was not thoroughly examined if (and possibly how) $T$ is dependent on $n$, such investigation would be beneficial for the overall study.\n\nIn the experiments on simulated data, the Acc metric is defined so as to measure only the ‘overall’ misestimation of the $\\omega_i$ values. However, $\\omega_i$ can either be -1, 0, or 1, and the zero values in the entire $\\mathbf{w}$ vector are predominant I presume. If this is the case, it would be insightful to measure metrics such precision or recall w.r.t. to each of the possible values in $\\mathbf{w}$. Micro/macro or weighted variants of these metrics might also be applicable. \n\nThe input parameter $k$ seems to be fixed in the simulated data experiments. I am wondering if such a parameter has an impact on run time of the proposed algorithm, which should have been analyzed under different feature dimensions. If such analysis was indeed conducted, I would encourage the authors to include and discuss the findings in the paper.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "Mo18OuvytBK", "reviewer": "Reviewer_NzEp", "summary": "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.", "questions": "", "limitations": "", "rating": 7, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The method in this paper nicely extends the idea of capped simplex projection using Newton’s method. While the formulation in (5) appears elsewhere, the insight here is that approximate numerical methods can find a similarly high-quality solution to existing methods. Excluding typos, the writing throughout is clear, and the paper is well-organized and fairly easy to follow. The limitations are discussed in a fairly upfront manner, such as the bounded vector magnitude and lack of exact solution. The primary concern comes from a lack of applications that may be of interest to the community. The GWAS results are not analyzed for any biological significance, and while faster than previous methods, does not unlock any analysis that was not achievable prior. Furthermore, it would be helpful to have scaling results as a function of the number of samples, since the considered example is a very small GWAS. Modern GWAS are on the order of 100k’s of individuals. One such example of further impact would be to pick an appropriate sparsity level (k) based on cross-validation or other measures, since robustness and reproducibility are key. This might highlight the advantages of the proposed method better, and have further applications to hypothesis testing in GWAS.\n\nUpdate to author response: please see discussion with authors regarding expanded experiments on synthetic GWAS data highlighting the effectiveness of the proposed method. Review score has been updated.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have sufficiently addressed limitations and there are no clear societal impacts.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "GCvP7RTuFxC", "reviewer": "Reviewer_Rv7X", "summary": "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.\nIn 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.\nThe 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.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is easy to read, and the theoretical contribution is fair, and the computational experiments are carried out well.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors claimed that the computational complexity of the proposed method is T O(n), where T is the number of Newton iterations, and only a few iterations were required empirically.\nHowever, this issue should be surveyed more, since linear time computational complexity is supposed to be one of the main contributions of this paper. An analysis about the worst case or average case running times should be studied.\n\nFigure 1 is not illustrative enough due to the lack of explanation both in the caption and the main text.\n\n#Typos\nIn abstract, CUP -> CPU\n\nIn section 2.1 rootof -> “root of”\n\nIn Theorem 1, convex -> concave ? since the second derivative is nonnegative.\n\nIn Theorem 1, “a indicator” -> “an indicator”", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}]

Author rebuttal:

Meta-review: 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.

Out of the 4 reviewers, 3 support accepting the manuscript. Detailed technical comments are provided in the reviews and the authors responded to these reviews in detail and appropriately. One reviewer suggests to reject the paper based on the fact that the manuscript has typos and an opinion that some elements of the manuscript would be better moved to the supplement. In essence, I think the low score by that reviewer is not well justified and in my own assessment of this manuscript I have excluded it. What's remaining are three reviewers that support the paper quite strongly overall (scores 7,7,6). I therefore recommend to accept this manuscript. 
Decision: Accept (Poster)