Beyond L1: Faster and Better Sparse Models with skglm
Authors: Quentin Bertrand
Quentin Klopfenstein
Pierre-Antoine Bannier
Gauthier Gidel
Mathurin Massias
Venue: NeurIPS 2022 Accept (2022)
Keywords: nonsmooth optimization
cooridnate descent
Anderson acceleration

Abstract: We propose a new fast algorithm to estimate any sparse generalized linear model with convex or non-convex separable penalties. Our algorithm is able to solve problems with millions of samples and features in seconds, by relying on coordinate descent, working sets and Anderson acceleration.  It handles previously unaddressed models, and  is extensively shown to improve state-of-art algorithms. We provide a flexible, scikit-learn compatible package, which easily handles customized datafits and penalties.

Beyond L1: Faster and Better Sparse Models with skglm
Quentin Bertrand Mila & UdeM, Canada
Quentin Klopfenstein
Luxembourg Centre for Systems Biomedicine
quentin. bertranddmila. quebec
University of Luxembourg
Asch-sur-Alzette, Luxembourg
Pierre-Antoine Bannier
Gauthier Gidel
Mathurin Massias
Independent Researcher
Mila & UdeM, Canada Canada CIFAR Al Chair
Univ. Lyon, Inria, CNRS, ENS de Lyon, UCB Lyon 1, LIP UMR 5668, F-69342 Lyon, France
Abstract
We propose a new fast algorithm estimate any sparse generalized linear model with convex non-convex separable penalties. Our algorithm is able to solve problems with millions of samples and features in seconds, by relying on coor- dinate descent, working sets and Anderson acceleration. It handles previously unaddressed models, and is extensively shown to improve state-of-art algorithms. We release skglm, a flexible, scikit-learn compatible package, which easily handles customized datafits and penalties.
Introduction
Sparse generalized linear models play central role in modern machine learning and signal processing. The Lasso (Tibshirani, 1996) and its derivatives (Zou and Hastie, 2005; Ng, 2004; Candes et al., 2008; Simon et al., 2013) have found numerous successful applications to large scale tasks in genomics (Ghosh and Chinnaiyan, 2005), vision (Mairal, 2010), or neurosciences (Strohmeier et al., 2016). This impact was made possible by two key factors: efficient algorithms and software implementations.
State-of-the-arr algorithms for "smooth non-smooth separable" problems predominantly rely on coordinate descent (CD, Tseng and S. Yun 2009; Nesterov 201 2), which, when it can be applied, is more efficient than full gradient methods (Richtárik and Takáč, 2014, Sec. 6.1). Coordinate descent can even be improved with Nesterov-like accelcration, to obtain improved convergence rates (Lin et al., 2014; Fercoq and Richtárik, 2015). However, these better rates may fail to reflect in practical accelerations. On the contrary, Bertrand and Massias (2021) relied on Anderson acceleration (Anderson, 1965) to provide both better rates and practical acceleration for coordinate descent.
Even with efficient algorithms such coordinate descent, the practical use of sparsity hits putational barrier for problems with more than millions of features (Le Morvan and Vert, 201 8). Multiple techniques have been proposed make coordinate descent scale huge problems. Notably, algorithms can be accelerated by reducing the number of variables to optimize over, using screening rules or working sets. Screening rules discard features from the problem in advance (El Ghaoui et al. 2010; Bonnefoy et al. 2015) or dynamically (Fercoq et al., 2015; Ndiaye et al., 201 7). On the other side, working sets (Johnson and Guestrin, 2015; Massias et al., 2018) iteratively solve larger subproblems and progressively include variables identified as relevant.
For the Lasso and few convex models, coordinate descent has been broadly disseminated to practitioners in off-the-shelf packages such as glmnet (Friedman et al., 2007) or scikit-learn
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

Lasso
MCP, 7 =
1
22
F1-score
RMSE
10
10
1
10
A/ 4/max
Figure Regularization paths computed with our algorithm. Non-convex sparse penalties behave better than the L1 norm. Due to their lower bias, they achieve perfect support recovery, lower prediction error and their optimal regularization strength À in estimation (top) and prediction (bottom) correspond.
(Pedregosa et al., 2011). More recently, celer, a sfate-of-the-art convex working set algorithm (Massias et al., 2020) allowed for successful applications the Lasso large scale problems in medicine (Reidenbach et al., 2021; Kim et al., 2021) or seismology (Muir and Zhan, 2021).
Yet the Lasso limited: non-convex sparse models enjoy better theoretical and empirical properties (Breheny and Huang, 2011; Soubies et al., 2015). As illustrated in Figure 1, they yield sparser solutions than convex penalties and mitigate the intrinsic Lasso bias. Yet, they have not often been applied to huge scale applications. This is mostly an algorithmic barrier: while coordinate descent can be applied non-convex penalties (Breheny and Huang, 2011; Mazumder al., 201 Bolte et al., 2014), screening rules and working sets are heavily dependent on convexity or quadratic datafits (Rakotomamonjy et al., 2019, 2022).
In this work, we solve this issue by designing a state-of-the-art generic algorithm to solve a wide range of sparse generalized linear models. The contributions are the following:
We propose non-convex converging working algorithm relying on Anderson accelerated coordinate descent. For a specific class non-convex penalties, we show:
(a) Convergence of the proposed working set algorithm (Proposition 5).
(b) Support identification coordinate descent (Proposition 10).
(c) Local convergence rates for the Anderson extrapolation (Proposition 13).
We provide an extensive experimental comparison and we show state-of-the th-the-arr improve- ments a wide range of convex and non-convex problems. In addition we release an efficient and modular python implementation, with scikit-lears API, for practitioners to apply non-convex penalties to large scale problems.
2 Framework and proposed algorithm
2.1 Problem setting
In this paper, we consider problems of the form:
argmin (BB) A F(X3) + BER"
Σ gg (,(
(1)
((())
where F' is smooth, and the functions 9j proper and lower semicontinuous but not necessarily convex, whose proximal operator can be computed exactly. We write = -; 9j- Instances of Problem (1) include convex estimators:
Lasso, the elastic net, the sparse logistic regression, the dual of SVM with hinge loss. They also include non-convex penalties: lo.5 and 22/3 penalties (Foucart and Lai, 2009), the minimax concave penalty (MCP, Zhang 2010) or SCAD (Zhang, 2010), both with regression and classification losses. Formally, the assumptions are the following.

---

R' Assumption 1. ::Rr - IR is convex and differentiable and for all j E [p]. the restriction of v, the j-th coordinate is Lipsschhtt::: for all (z, E RP R. fV he,) v L, L
Assumption 2. For any j E pp.g : R + IR is proper, closed, and lower bounded.
Following Attouch and Bolte 2009; Bolte et al. 2014 we focus on finding a critical point of
Definition 3. Using the Fréchet subdifferential (Kruger, 2003), a critical point a E R is a point which satisfies f(x 0g(z).
Assumptions and 2 are usual, and, under boundedness of the iterates, ensure convergence of and coordinate descent algorithms to a critical point (Attouch et al. 2013, Thm 5.1, Bolte et al. 2014, Thm. 3.1). In addition, our work focuses the case where 9j 's present non-differentiabiiiy points, leading to the following extended notion of sparsity.
Definition 4 (Generalized support). The generalized support 3 E RP is the set of indices j E [p] such that 9j is disserentiable at By: gsupp(3) = (j (p] Ogy (3,) is singleton].
Penalties such as L,, bg 9 1), MCP SCAD only not differentiable at O, and this corresponds to the usual notion of sparsity. But Definition 4 goes beyond sparsity and extends to estimators such as SVM, where gj ljo,C] and the generalized support is the complement of the support vectors' set {j [p] BJ 0 or 33 C). The generalized support of critical point is usually of cardinality much smaller than p. and its knowledge makes the problem easier and faster to solve. Our working set algorithm exploits this structure in order converge faster.
2.2 Proposed algorithm
The proposed algorithm exploits two main ideas:
A working set strategy, able handle large class of convex and non-convex penalties (Algorithm 1).
An Anderson accelerated coordinate descent for non-convex problems (Algorithm 2). The building blocks of Algorithm 2, coordinate descent (CD, Algorithm 3) and Anderson extrap- olation (Anderson, Algorithm 4), can be found in Appendix A.
To avoid wasting computation on features outside the generalized support, working set algorithms iteratively select subset of coordinates deemed important (the working set), and solve Prob- lem (1) restricted to them. The key question is thus the notion of important features. Stem- ming from Definition 3, we rank features by their violation of the optimality condition: score dist - (3),, 0g, (8)) For example, the MCP Fréchet subdifferential at 0g, (0) [-X, A].
and the proposed score reads
score 3 max (0, (,f(()) A} if 3, 0
(2)
Vg, otherwise
To control the working set growth, we use score rank the features. Then, with n& max(n*-1,2 2 gsupp(3(())) we take the ny largest of them in the working set, while retaining features currently in the working set. This growth quickly rises the unknown size of the generalized support while avoiding overshooting, as backed up by recent theory in Ndiaye and Takeuchi (2021).
Proposition 5. Let Wt be the 1-th working set. Suppose that Algorithm 2 converges toward critical point, and for all 2 0, W, C WH4,, then the iterates of Algorithm I converge towards a critical point of Problem (1).
3

---

Algorithm 1 skglm (proposed)
Algorithm 2 inner inner_solver
input: X, IR? Nout N,
input: :X, (0)) E R", ws C [p], Rin 6, M = 5
C,
7hin E N, ws_ size E N 0
for k = ...., Rinn do
for / 1, ...., Mout do
g(k) CD( , 3(=-1), X3, ws) // Algo. (3)
score dist (- V (3)) (g,) jelp ws_size = max wss_iize, 2 x g3)p
if mod M = then
// Algo. (4), OMM + Ma)
// ws_size features with largest scores
fext < Andersonn( (()-M)))))))))
(M
1/ test objective O(n[ws))
ws arg 1ss_topK(scoree, k= ws size)
if D (B
(((ws) then
if' m dist ( (V)
e
3 gextr, X,3 + w gextr
then stop
maxjews dist (gf(3). S
else // accelerated on working set
3 4 ieniinnerssee Sonner_solver((e (X, B, ws, M)mn
then stop
return (()
8
return :3
Proof of Proposition can be found in Appendix B. The second key ingredient our algorithm is to use state-of-the-an Anderson accelerated coordinate descent for non-convex problems. In Section 2.3 we show that coordinate descent yields finite time support identification for large class of non- convex problems (Proposition 10), which leads to acceleration (Proposition 13). As experiments demonstrate in Section 3, this rate allows our algorithm to surpass state-of- sfate-of-the-art t-t-hhe-at solvers.
2.3 Anderson accelerated coordinate descent analysis for a-semi-convex penalties
We now turn to our main technical contributions: we show that Algorithm 2 achieves finite time support identification (Proposition 10) of the generalized support (Definition 4) for specific class of non-smooth non-convex penalties (Assumption 6), which includes the MCP (Proposition 7). Based on Proposition 10, we are able to derive convergence rates for Anderson acceleration (Proposition 13).
We study our inner solver (Algorithm 2); for convenience we still refer to 3 and X for their counter- parts restricted to the working set. The following assumptions are required.
Assumption fe--smi-cconeer For all j E [pl 9 Ly is a a-semi-convex, i.e., g3/ al1-113/2 convex, with < 1.
Note that in statistics, the admissible value range of hyperparameters for MCP and SCAD are datafit-
dependent, (see Breheny and Huang 2011, Sec. 2.1, normalized columns and y > 1 = 1/ [X,]| = 1,L, or Soubies et al. 2015, Eq. 4.2) and yields a-semi- conveiit for MCP and SCAD
A/E 27
i < 7X if
Proposition 7 (a-semi-convexil of MCP). Ler MCPA,
1
Ify > 1, then MCP Ly is a-semi-convex with a = [1 + Tt,)) (i.c., Assumption 6 holds).
Note that Assumption 6 does not hold for the 4,-penalties (0 < q < 1), for which we propose an alternative in Appendix C.
Assumption (Existence). Problem (1) admits at least critical point.
In Proposition 10, convergence of Algorithm 2 toward a critical point 3 is assumed, and the following assumption is made on this critical point.
Assumption 9 (Non degeneracy). The considered critical point 3 E IR is non-degenerated: for all j & g5)
-5) (3)e
Assumption is a generalization of qualification constraints (Hare and Lewis, 2007, Sec. 1), and is usual in the machine learning literature (Zhao and Yu, 2006; Bach, 2008; Vaiter et al., 2015). For
Howeeer MCP and SCAD not a-semiconvex for all hyperparameter values.

---

the nooro,m if the entries of the design matrix X are drawn from an i.idd normal distribution, then Assumption holds with high probability (Candes and Tao, 2005; Rudelson and Vershynin, 2008).
Equipped with the previous assumptions we show that coordinate descent achieves model identifica- tion for this class of non-convex problems.
Proposition 10 (Model identification of CD). Suppose
1. Assumptions 1. 2, 6 and hold.
2. The sequence (3th) )22 generated by coordinate descent (Algorithm 2 without extrapola- tion) converges toward a critical point 3.
3. Assumption 9 holds for B.
Then, Algorithm (without extrupolation) identifies the model in finitely many iterations: there exists K 0 such that for all k 2 K, 33 3sc.
In other words, for k large enough, ( shares the generalized support of 3. The identification property was proved for a proximal gradient descent algorithm in the non-convex case (Liang et al., 2016) under the assumption that the non-smooth function g is partly smooth (Lewis, 2002). For ourselves, Proposition not rely on the partly smooth assumption to ensure identification property. Authors not aware of previous identification results for coordinate descent in the non-convex case.
In addition, if f and g are locally regular on the generalized support at the considered critical point, algorithm enjoys local acceleration when combined with Anderson extrapolation (Proposition 13).
Assumption 11 (Locally (²)). For all E S _ gsupp(3), 9j is locally C around By, and f is locally C around B.
Assumption 11 on the function f is mild and holds for usual machine learning datafitting terms. Assumption 11 on the functions 9j. E S. is stronger: for instance, for the MCP, it implies 3, * 7 for all j E S. However this assumption is standard the literature, see Liang et 2016, Sec. 3.3
Assumption 12. (Local strong convexity) The Hessian of f at the considered critical point IR,, restricted to its generalized support S. is positive definite, i.e., ( V3,5g(3) > 0.
Assumption 12 requires local strong convexity restricted to the generalized support S, which is standard in the MCP SCAD literature (Breheny and Huang 2011, Section 4.1) and is usual derive local linear rates of convergence (Liang al., 2016, Section 3.3). For instance, for the Lasso, if the entries of the design matrix X are drawn from continuous distribution, then Assumption 12 holds with probability one (Tibshirani, 2013, Lemma 4).
Proposition 13. Consider critical point 9 and suppose
1. Assumptions 1, 2 and hold.
2. The functions f and gjz j E (p) are piecewise quadratic (which is the case for the MCP regression).
3. The sequence (g(4)) k22 generated by Anderson accelerated coordinate descent with updates from to and to 1 (Algorithm with extrapolation) converges to a critical point B.
4. Assumptions 9, 11 and 12 hold for 3.
Then there exists K N, and a C¹ function vb RISI
- RSSl such that, for all le E N,N K:
g = By for all S Let & Jo(B). A V2 ( + V (B). S
(4)
(1-((T) and
B B (T- Id). Then P(T) < 1 and the iterates of Anderson extrapolation enjoy local accelerated convergence rate:
1 (k-K) A < (VK(H) 201-1 (k-K (4-K)/
14(20M-13
(1 -
(5)
The proof can be found in Appendix B.5.
5

---

Table 1: Most popular packages for sparse generalized linear models.
Name
Acceleration
Huge scale
Nnevx
Modular
glmnet (Friedman et al., 2010)
X
X
X
(Fortran)
scikit-learn (Pedregosa et al., 2011)
XXX
X
X
X
X
(Cython)
lightning (Blondel and Pedregosa, 2016)
X
X
V (Cython)
celer (Massias al., 2018)
J
X
X (Cython)
picasso (Ge et al., 2019)
X
X
V
X (C++)
PyGLMnet (Jas et al., 2020)
XX
X
(Python)
fireworks (Rakotomamoniy et al., 2022)
X
/
N.A. (Python)
skglm (ours)
(Python)
Related work. Most Anderson acceleration convergence results are shown for quadratic objectives for specific algorithms: gradient descent (Golub and Varga, 1961; Anderson, 1965), ADMM (Poon and Liang, 2019), coordinate descent (Bertrand et al., 2020). Outside of the quadratic case, convergence results are usually significantly weaker (Scieur et al., 2016; Sidi, 2017; Brezinski et al., 2018; Mai and Johansson, 2019; Ouyang et al., 2020). Regarding the smooth non-convex case, Wei et al. (2021) proposed a stochastic Anderson acceleration and proved convergence towards critical point. Proposition 13 generalizes Scieur et al. (2020, Prop 2.1) and Bertrand and Massias (2021, Prop. 4) to the proximal convex and c-semi-convex cases. To our knowledge this is one of the first quantitative results for Anderson accelcration in a non-convex setting.
2.4 Comparison with existing work
In this section we compare our contribution to existing algorithms and implementations, which are summarized Table 1. Huge scale refers to the fact that the algorithm can run on problems with millions of variables. Non-convex tells if the algorithm handles non-convex penalties. Modular indicates that it is easy to add a new model, through a different datafitting term or penalty.
The packages glmnet (Friedman et al., 2010), scikit-leand (Pedregosa et al., 2011) and lightning (Blondel and Pedregosa, 2016) implement coordinate descent (cyclic or random). They rely on compiled code such as Fortran or Cython, making it very difficult to implement new models' or faster algorithms like working set They do not handle non-convex penalties.
More recent algorithms such as blitz (Johnson and Guestrin, 2015), celer (Massias et al., 2018), 1)), picasso (Ge et al., 2019) fireworks (Rakotomamonjy al., 2022) use working set strategies. celer and blitz are state-of-the-art algorithms for the Lasso, but their score to prioritize features relies on duality. fireworks extends blitz convex functions), with scorefirework: subdifferential of g at the current point, but which a coarse information. Finally, fireworks, building upon the seminal non convex working set solver of Boisbunon et al. (2014), does not provide accelerated convergence rates and does not come with a public implementation, pi pasao (Ge et al., 2019) lacks modularity (penalties are hardcoded), and the solver is not suited for huge scale (it does not support large sparse matrices). Deng and Lan (2019) proposed an algorithm based on inertially accelerated coordinate descent, which fails to provide practical speedups according to Bertrand and Massias (2021).
non-convex penalties (writing as difference of dist V (3). 0g, (0)). Yet this rule does not consider the
Contrary to these algorithms, ours is generic and relies only on the knowledge of V and prox, 9. For any new penalty, this information can be written few lines of Python code, compiled with numba (Lam et al., 2015) for speed efficiency. We therefore improve state-of-the-ar algorithms in the convex case, and generalize virtually any datafit and penalty, even nonconvex.
Experiments
Our package relying numpy and numba (Lam et al., 2015; Harris et al., 2020) is attached in the supplementary material. An open source, fully tested and documented version of the code can be
years old)
3 wwwpscikitttiiieee
(5 years old)

---

blitz
celer
skglm (ours)
sklearn
Amax / 10
A 1100
Xmnx 11000
10
news20
10
0.0
0.5
1.0
0.0
2.5
5.0
0
20
10
kdda
10
20
50
100
0
200
400
10
url
10
0
100 Time (s)
200
0
500
1000 0
2000
Time (s)
Time (s)
Figure 2: Lasso, duality gap. Normalized duality gap as a function of time for the Lasso multiple datasets, for multiple values of X.
found https://githu:
www./oociiiieee. contriio/ntibbllnnn. We use datasets from 1ibsvm (Fan et al. 2008, see table 2).
We compare multiple algorithms solve popular Machine Learning and inverse problems: Lasso, Elastic net, multitask sparse regression, MCP regression. The compared algorithms are the following:
scikit-learn (Pedregosa et al., 201 which implements coordinate descent in Cython, celer (Massias al., 2020), which combines working sets, screening rules, coordinate descent, and Anderson acceleration in the dual, in Cython,
blitz (Johnson and Guestrin, 2015), which combines working sets with prox-Newton iterations (Lee al., 2012) in C++,
coordinate descent (CD, Tseng and S. Yun 2009),
skglm (Algorithm 1, ours), using M iterates for the Anderson extrapolation.
Other solvers. Experiment per experiment, there exist niche solvers (such aggressive Gap Safe Rules, Ndiaye et al. 2020). Since our goal is a general purpose algorithm able to deal with many models, we do not include them in the comparison. In addition, we focus on solving single instance of Problem (1), rather than a regularization path (i.e., a sequence of problems for multiple regularization strengths). As glmnet is designed to compute regularization paths, we could not include in the comparison. The reader can refer Johnson and Guestrin (2015, Fig. 4) or Figure in Appendix E for comparisons on single optimization problems with glmnet; glmnet and additional algorithms are discussed in Appendix E.
How to do fair comparison between solvers? To plot the convergence curves, we use the benchopt benchmarking package (Moreau et al., 2022). In order automate and reproduce optimization benchmarks it treats solvers black boxes. It launches them several times with increasing maximum number iterations, and stores the resulting objective values and times to reach it. As each point on solver curve is obtained in different run, the curves are not monotonic, and there may be several points corresponding to the same time. This merely reflects the variability in solvers running time across runs; we refer to Figure 10 in Appendix E.6 for the inevitability of this phenomenon with black box solvers.
csie. wwwp://..ei...t..u..uu....... 'https://gitan htt.t///..ccc/bcccccp/////oooooe
7

---

skglm (ours)
coordinate descent
sklearn
X 10
Amax /100
mmax/ Amax/1000 /1000
10
finance
10 8
50
100
0
200
400
Time (s)
Time (s)
Time (s)
Figure 3: Elastic net, duality gap. Normalized duality gap as a function of time for the elastic net for multiple values of X, p 0.5.
(a) 62.1-
(b) l2.0..
(c) Block MCP.
Figure 4: Real data, brain source locations recovered by convex and non-convex penalties after right auditory stimulation. 4(a) shows that a convex penalty fails at identifying one source in each hemisphere, while 4(b) and 4(c) demonstrates the capability of non-convex penalties to recover the correct solution.
3.1 Convex problems
Lasso. In Figure 2 we compare solvers for the Lasso: (f = 2m y/ - 1x. = A).). We
parametrize À as fraction of Amax || //l.o/nn /n, smallest regularization strength for which 3 0. For large scale datasets (rcvl, news20), skglm yields performances better or similar to the otate-of- the-art algorithms blitz and celer. For huge scale datasets (kdda and url), skglm yields significant speedups over them. The improvement over the popular scikit-lears can be of two orders of magnitude. Thus, while dealing with many more models, our algorithm still yields ofate-of-the-art speed for basic ones.
Elastic net. Our approach easily generalizes to other problems, such as the clastic net (f In - X-11,, 9j A(pl-1 + - Figure 3 shows the duality gap as function of time for skglm (ours), sklearn, and our numba implementation of coordinate descent. The proposed algorithm is orders of magnitude faster than scikit-learn and vanilla coordinate descent, in particular for large datasets and low regularization parameter values (finance, Amaxx/ Amaxx/1000) 11000). Note that blitz does not implement a solver for the elastic net. Many Lasso solvers would casily handle the elastic net, but relying on Cython/C++ code makes the implementation cinecconu,,ng, By contrast, it takes 40 lines of code to define an & + Es-squreed penalty with our implementation. An additional experiment on the dual of SVM with hinge loss Appendix E.4.
3.2 Non-convex problems
In this subsection we propose a comparison on two non convex problems.
MCP regression. MCP regression is Problem (1) with f = 1 11 - X x-1 9j MCP x, for 2n.
> 1. As usual for this problem, we scale the columns of X to have norm vn. On Figure 5, we compare our algorithm to pi passo on a dense dataset (n 1000, p = 5000); as this package does not support large sparse design matrices, for the rev/ dataset we use an iterative reweighted I.I algorithm (Candes et al., 2008). Since the derivative of the MCP vanishes for values bigger than Xy, this approach requires solving weighted Lassos with some 0 weights. Up our knowledge, our
8

---

coordinate descent
picasso
skgim (ours)
Amax/ /10 10
Amax/ /100
Xmax 11000
objective
1000
1000
1000
500
0-
0
opt. cond.
10-6
10
6
10-6
2
Time (s)
Time (s)
Time (s)
coordinate descent
iterative reweighted L1
skglm (ours)
Amax/ 10
Amax/ /100
Amax/ / /1000
opt. cond. objective
0.50
0.50
0.5
0.25
0.25
0.0
10
10
8
10
0.0
0.2
0.4
0
20
100
Time (s)
Time (s)
Time (s)
Figure MCP, objective value and violation of first order condition. Objective value and violation of optimality condition of the iterates, dist(-V/(3th for the MCP for multiple values of À (y = 3) on simulated dense dataset (top) and the revl dataset (normalized columns).
ag(3()))) (,g(3)))))) as function of time
algorithm is the only efficient one with such a property. Our algorithm handles problems of large size, converges to a critical point, and, due to its progressive inclusion of features, is able to reach sparser critical point than it competitors.
Application to neuroscience To demonstrate the usefulness of our algorithm for practitioners, we apply it the aegencto/elcctoennophallogggg
(M/EEG) inverse problem. It consists in reconstructing the spatial cortical current density at the origin of M/EEG measurements made at the surface of the scalp. Non-convex penalties (Strohmeier et al., 2015) exhibit several advantages over convex ones (Gramfort et al., 2013): they yield sparser physiologically plauusbble solutions and mitigate the l amplitude bias. Here the setting is multitask: Y xx and thus we use block penalties (details in Appendix D). We use real data from the mne software (Gramfort et al., 2014); the experiment is right auditory stimulation, with two expected neural sources recover in each auditory cortex. In Figure 4, while the &,1 penalty fails at localizing one source in each hemisphere, the non-convex penalties recover the correct locations. This emphasizes on the critical need for fast solvers for non-convex sparse penalties as well as our algorithm's 's ability to handle the latter. In this work we focused on optimization- based estimators to solve the inverse problem, note that one could have resort to other techniques, such as Bayesian techniques (Ghosh and Doshi-Voss Velez, 2017; Fang et al., 2020).
Ablation study. To evaluate the influence of the two components of Algorithm 1, an ablation study (Figure 6) is performed. Four algorithms are compared: with/without working sets and with/without Anderson acceleration. Figure represents the duality gap of the Lasso as function of time for multiple datasets and values of the regularization parameters (parametrized as a fraction of Amax). First, Figure 6 shows that working sets always bring significant speedups. Then, when combined with working set, Anderson acceleration bring significant speed-ups, especially for hard problems with low regularization parameters. An interesting observation is that on large scale datasets (news20 and

---

CD Warking sets and Acceleration
CD
Acceleration
CD Working sets
CD
10"
Amas/ 110
Amas
1100
Amax
1000
rcv1
10 10-10 10
0.0
0.2
2
0
5
10
10
news20
10 10
10
0
20
40
0
20
40
100
finance
10 10
0
20
40
0
20
0
20
40
Time (s)
Time (s)
Time (s)
Figure 6: Lasso, duality gap. Normalized duality gap as a function of time for the Lasso.
finance) and for low regularization parameters (Amax/ 1000 and Amax/ 1000) Anderson acceleration without working set does not bring acceleration. This highlights the importance of combining Anderson acceleration with working sets.
Conclusion and broader impact. In this paper, we have proposed an accelerated versatile algo- rithm for specific class of non-smooth non-convex problems. Based on working sets, coordinate descent and Anderson acceleration, we have improved state of the art on convex problems, and han- dled previously out-of-reach problems. Thorough experiments demonstrated the speed and interest of our approach. A limitation of this work is the considered function class (a-semi-convex) which can seen as restrictive. One possible extension would be weakly convex functions (Davis and Drusvyatskiy, 2019, Sec. 1). We deeply believe that the high quality code provided will benefit to practitioners, and ease the use of non-convex penalties for real world problems, from neuroimaging to genomies. We proposed optimization algorithm and do not see potential negative societal impacts.
Acknowledgements
The experiments were ran on the CBP cluster of ENS de Lyon (Quemener and Corvellec, 2013). QB would like to thank Samsung Electronics Co., Ldt. for funding this research. GG is supported by an IVADO grant.
10

---

References
D. G. Anderson. Iterative procedures for nonlinear integral equations. Journal of the ACM, 12(4): 547-560, 1965.
H. Attouch and J. Bolte. On the convergence of the proximal algorithm for nonsmooth functions involving analytic features. Mathematical Programming, 116(1):5-16, 2009.
H. Attouch, J. Bolte, and B. F. Svaiter. Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods. Mathematical Programming, 137(1):91-129 2013.
Bach. Consistency of the group Lasso and multiple kernel learning. Mach. Learn. Res., 1179-1225, 2008.
Q. Bertrand and M. Massias. Anderson acceleration of coordinate descent. In AISTATS, 2021.
Q Bertrand, Q Klopfenstein, M. Blondel, S. Vaiter, A. Gramfort, and J. Salmon. Implicit differentia- tion of lasso-type models for hyperparameter optimization. ICML, 2020.
M. Blondel and Pedregosa. Lightning: large-scale linear classification, regression and ranking in python, 2016.
A. Boisbunon, R. Flamary, and A. Rakotomamoniy. Active set strategy for high-dimensional non- convex sparse optimization problems. In 2014 IEFE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 1517-1521. IEEE, 2014.
J. Bolte, S. Sabach, and M. Teboulle. Proximal alternating linearized minimization for nonconvex and nonsmooth problems. Mathematical Programming, 146(1):459-494, 2014.
A. Bonnefoy, Emiya, L. Ralaivola, and R. Gribonval. Dynamic screening: accelerating first-order algorithms for the Lasso and Group-Lasso. IEEE Trans. Signal Process., 63(19)200, 2015.
S. Boyd, N. Parikh, and F. Chu. Distributed optimization and stafistical learning via the alternating direction method of multipliers. Now Publishers Inc, 2011.
Breheny and Huang. Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection. The annals of applied statistics, 5(1):232, 2011.
C. Brezinski, M. Redivo-Zaglia, and Y. Saad. Shanks sequence transformations and anderson acceleration. SIAM Review, 6063):646-669, 2018.
E. J. Candes and Tao. Decoding by linear programming. IEEE transactions on information theory, 51(12):4203-422 2005.
E. J. Candes, M. B. Wakin, and S. P. Boyd. Enhancing sparsity by reweighted ε minimization. Journal of Fourier unalysis and upplications, 14(5):877-905, 2008
D. Davis and D. Drusvyatskiy, Stochastic model-based minimization weakly convex functions. SIAM Journal on Optimization, 29(1) 29(1):2072239, 2019.
Q Deng and C. Lan. Efficiency of coordinate descent methods for structured nonconvex optimization. arXiv preprint arXiv: 00900, 2019.
Ghaoui, Viallon, and T. Rabbani. Safe feature elimination for the lasso and sparse supervised learning problems. arXiv preprint arXiv: arXiv:1009 (09,,4219, 2010.
R. E. Fan, K. W. Chang, C. J. Hsieh, X. R. Wang, and C. Lin. Liblinear: A library for large linear classification. JMLR, 9: 18877,1 2008.
S. Fang, Zhe, K.-C. Lee, K. Zhang, and J. Neville. Online bayesian sparse learning with spike and slab priors. In 2020 IEEE International Conference Data Mining (ICDM), pages 142-151. IEEE, 2020.
O. Fercoq and Richtárik. Accelerated, parallel, and proximal coordinate descent. SIAM Journal on Optimization, 25(4):1997-2022) 2015.

---

O. Fercoq, A. Gramfort, and J. Salmon. Mind the duality gap: safer rules for the lasso. In ICML, pages 333-342. PMLR, 2015.
Foucart and M.-J. Lai. Sparsest solutions of underdetermined linear systems via l mi-iiiiinzntion for 0 < 9 I 1. Applied and Computational Harmonic Analysis, (6(3):395-400 2009.
J. Friedman, J. Hastie, H. Höfling, and R. Tibshirani. Pathwise coordinate optimization. Ann. Appl. Stat., 1(2):302-332 202, 2007.
Friedman, T. Hastie, and Tibshirani. Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33(1):1, 2010.
J. Ge, X. Li, H. Jiang, H. Liu, T. Zhang, M. Wang, and T. Zhao. Picasso: A sparse learning library for high dimensional data analysis in and python. The Journal of Machine Learning Research, 20 (1):1692-1668) 2019.
D. Ghosh and A. M. Chinnaiyan. Classification and selection of biomarkers in genomic data using lasso. Journal of Biomedicine and Biotechnology. 2005(2):147, 2005.
S. Ghosh and Doshi- Model selection in bayesian neural networks via horseshoe priors. arXiv preprint arXiv: 10388, 2017.
G. H. Golub and R. S. Varga. Chebyshev methods, and second order richardson iterative methods. Numerische Mathematik, 3(1::147-156 1961.
methods, successive overrelaxation iterative
A. Gramfort, D. Strohmeier, J. Haueisen, M. S. Hämäläinen, and M. Kowalski. Time-frequency mixed- norm estimates: Sparse M/EEG imaging with statioati source activations. Neurolmage, 70: 410-422, 2013.
A. Gramfort, M. Luessi, E. Larson, D. A. Engemann, D. Strohmeier, C. Brodbeck, L. Parkkonen, and M. S. Hämäläinen. MNE software for processing MEG and EEG data. Neurolmage, 86:446 4600, 2014.
W. I.. Hare and A. S. Lewis. Identifying active manifolds. Algorithmic Operations Research, 2(2):
75-75, 2007.
C. R. Harris, K. Millman, S. J. van der Walt, Gommers, P. Virtanen, D. Cournapea,,
E. Wieser, J. Taylor, S. Berg, N. J. Smith, et al. Array programming with numpy. arXiv preprint arXiv:2006. 10256, 2020.
M. Jas, Achakulvisut, A. Idrizović, D. Acuna, M. Antalek, V. Marques, T. Odland, R. Garg, M. Agrawal, Y. Umegaki, P. Foley, H. Fernandes, D. Harris, B. Li. O. Pieters, S. Otterson, G. De Toni, C. Rodgers, E. Dyer, M. Hamalainen, K. Kording, and Ramkumar. Pyglmnet: Python implementation of elastic-net regularized generalized linear models. Journal of Open Source Software, 547): 1959, 2020.
B. Johnson and C. Guestrin. Blitz: A principled meta-algorithm for scaling sparse optimization. In ICML, volume 37, pages 1171-1179, 2015.
Y. J. Kim, N. Brackbill, E. Batty, J. Lee, C. Mitelut, W. Tong, EJ Chichilnisky, and L.. Paninski. Nonlinear decoding of natural images from large-scale primate retinal ganglion recordings. Neural Computation, (3(7):1719-177 2021.
Q. Klopfenstein, Q. Bertrand, A. Gramfort, J. Salmon, and S. Vaiter. Model identification and local linear convergence of coordinate descent. arXiv preprint arXiv:2010.11825, 2020.
A. Y. Kruger. On Fréchet subdifferentials. Journal of Mathematical Sciences, 116(3):3325-3358 2003.
S. K. Lam, A. Pitrou, and S. Seibert. Numba: A llym-based python jit compiler. In Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, pages 1-6, 2015.
M. Morvan and J.-P. Vert. WHInter: working set algorithm for high-dimensional sparse second order interaction models. In ICML, pages 3635-3644, PMLR, 2018.
12

---

J. D. Lee, Y. Sun, and M. Saunders. Proximal Newton-type methods for convex optimization. Advances in Neural Information Processing Systems, 25:827-835, 2012.
A. S. Lewis. Active sets, nonsmoothness, and sensitivity. SIAM Journal on Optimization, 13(3): 702-725, 2002.
J. Liang, J. Fadili, and G. Peyré. A multi-step inertial forward-backward splitting method for non-convex optimization. Advances in Neural Information Processing Systems, 29:4035-4044 2016.
Q. Lin, Z. Lu, and L. Xiao. An accelerated proximal coordinate gradient method. NeurlPS, pages 3059-3067. 2014.
D. C. Liu and J. Nocedal. On the limited memory BFGS method for large scale optimization. Mathematical programming, 555))5803-528 1989.
V. Mai and M. Johansson. Anderson acceleration of proximal gradient methods. JCMI.. 2019.
Mairal. Sparse coding for machine learning. image processing and computer vision. PhD thesis, Acolle normale supérieure de Cachan, 2010.
M. Massias, A. Gramfort, and J. Salmon. Celer: fast solver for the lasso with dual extrapolation. 2018.
M. Massias, S. Vaiter, A. Gramfort, and Salmon. Dual extrapolation for sparse generalized linear models. J. Mach. Learn. Res., 2020.
R. Mazumder, J. H. Friedman, and T. Hastie. Sparsenet: Coordinate descent with nonconvex penalties. Journal of the American Statistical Association, 106(495):1 1138, 2011.
Moreau, M. Massias, A. Gramfort, P. Ablin, B. Charlier, P.-A. Bannier, M. Dagréou, T. Dupré la Tour, G. Durif, C. F. Dantas, Q Klopfenstein, et al. Benchopt: Reproducible, efficient and collaborative optimization benchmarks. arXiv preprint arXiv:2206. 13424, 2022.
J. B. Muir and 7. Zhan. Seismic wavefield reconstruction using pre-conditioned cavelet-curvelet compressive sensing approach. Geophysical Journal International, 227(1):303-315, 2021.
E. Ndiaye and Takeuchi. Continuation path with linear convergence rate. arXiv preprint
12.05:211..05504, 2021.
E. Ndiaye, O. Fercoq, A. Gramfort, and J. Salmon. Gap safe screening rules for sparsity enforcing penalties. J. Mach. Learn. Res., 18(128):1-33 2017.
E. Ndiaye, O. Fercoq, and Salmon. Screening rules and complexity for active set identification. Journal of Convex Analysis, 2020.
Y. Nesterov. Efficiency of coordinate descent methods on huge-scale optimization problems. SIAM Journal on Optimization, 22(2):341 2012.
A. Y. Ng. Feature selection, vs. 12 regularization, and rotational invariance. In ICML, page 78, 2004.
J. Nutini. Greed is good: greedy optimization methods for large-scale structured problems. PhD thesis, University of British Columbia, 2018.
J. Nutini, M. W. Schmidt, I. H. Laradji, M. P. Friedlander, and H. A. Koepke. Coordinate descent converges faster with the Gauss-Souther rule than random selection. In ICML, pages 1641, 2015.
W. Ouyang, Y. Peng. Y. Yao, J. Zhang, and B. Deng. Anderson acceleration for nonconvex ADMM based on Douglas- Rachford splitting. Computer Graphics Forum, volume 39, pages 221-239. Wiley Online Library, 2020.
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Pretten- hofer, R. Weiss, V. Dubourg, Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learm: Machine learning in Python. JMLR, 12:2825-2834 2011.
13

---

C. Poon and J. Liang. Trajectory of alternating direction method of multipliers and adaptive accelera- tion. In NeurlPS, pages 7357-7365, 2019.
E. Quemener and M. Corvellec. Sidus- the solution for extreme deduplication of an operating system. Limux Journal, 2013(235):3, 2013.
A. Rakotomaaonij, G. Gasso, and J. Salmon. Screening rules for lasso with non-convex sparse regularizers. In ICML, pages 5341-5350, 2019.
A. Rakotomamonjy, R. Flamary, G. Gasso, and J. Salmon. Provably convergent working set algorithm for non-convex regularized regression. In AISTATS, 2022.
D. A. Reidenbach, A. Lal. L. Slim, O. Mosafi, and J. Isracli. Gepsi: A python library to simulate gwas phenotype data. bioRxiv, 2021.
P. Richtárik and M. Takáč. Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function. Mathematical Programming, 144(1-2):1-38, 2014.
M. Rudelson and R. Vershynin. On sparse reconstruction from fourier and gaussian measurements. Communications Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 61(5):1025-1045 2008.
D. Scieur, d' Aspremont, and F. Bach. Regularized nonlinear acceleration. In Advances Neural Information Processing Systems, pages 712-720, 2016.
D. Scieur, A. d' A'Apemmmot,, and F. Bach. Regularized nonlinear acceleration. Mathematical Program- ming, 179(1):47-83, 2020.
A. Sidi. Vector extrapolation methods with applications. SIAM, 2017.
N. Simon, J. Friedman, T. Hastie, and R. Tibshirani. A sparse-group lasso. J. Comput. Graph. Statist., 22(2):231-245 2013.
E. Soubies, L. FéracdFéraud, and G. Aubert. A continuous exact &n penalty (cel0) for least squares regularized problem. SIAM Journal on Imaging Sciences, 83)):1077-669,, 2015.
D. Strohmeier, A. Gramfort, and J. Haueisen. MEG/EEG source imaging with non-convex penalty in the time-frequency domain. In Pattern Recognition in Neurvimaging, 2015 International Workshop 2015.
D. Strohmeier, Y. Bekhti, J. Haueisen, and A. Gramfort. The iterative reweighted mixed-nome estimate for spatio-temporal MEG/EEG source reconstruction. IEEE transactions on medical imaging, 35(10):2218-222) 2016.
R. Tibshirani. Regression shrinkage and selection via the lasso. J. R. Stat. Soc. Ser. B Stat. Methodol., 58(1):267-288, 1996.
R. Tibshirani, J. Bien, Friedman, T. Hastie, N. Simon, J. Taylor, and R. J. Tibshirani. Strong rules for discarding predictors in lasso-type problems. Journal of the Royal Statistical Society: Series B (Statistical Methoddoogy), 74(2):245-266. 2012.
R. J. Tibshirani. The lasso problem and uniqueness. Electronic Journal of statistics, 7:1456-1490, 2013.
Tseng and S. Yun. A coordinate gradient descent method for nonsmooth separable minimization. Mathematical Programming, 117(1):8774423, 2009.
S. Vaiter, G. Peyré, and J. Fadili. Low complexity regularization of linear inverse problems. In Sampling Theory, a Renaissance, pages 103-153. Springer, 2015.
Wei, C. Bao, and Y. Liu. Stochastic Anderson mixing for nonconvex stochastic optimization.
Advances in Neural Information Processing Systems, 34, 2021.
F. Wen, L. Chu, P. Liu, and R. Qiu. A survey nonconvex regularization-basee sparse and low-rank recovery in signal processing, statisties, and machine learning. IEEE Access, 6:69883 2018.
14

---

C.-H. Zhang. Nearly unbiased variable selection under minimax concave penalty, The Annals of statistics, 9822):894-94,, 2010.
Zhao and B. Yu. On model selection consistency of lasso. Mach. Learn. Res., 7:2541-2563, 2006.
H. Zou and T. J. Hastie. Regularization and variable selection via the elastic net. J. R. Stat. Soc. Ser. Stat. Methodol,, 67(2))301-324 -300,, 2005.
15

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you descrihe the limitations of your work? [Yes] See limitations paragraph
(c) Did you discuss any potential negative societal impacts of your work? [Yes]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Proposi- tions 10, 13 and 14.
(b) Did you include complete proofs of all theoretical results? [Yes] See Appendix B.
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or URL)? [Yes]
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] See Section 3.
(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [N/A]
(d) Did you include the total amount of compute and the type of resources used (c.g., type of GPUs, internal cluster, cloud provider)? [N/A]
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] In particular we acknowledge the python ecosystem.
(b) Did you mention the license of the assets? [N/A]
(c) Did you include any new assets either in the supplemental material or as a URL? [Yes]
(d) Did you discuss whether and how consent was obtained from people whose data you're 're using/curating? [N/A]
(e) Did you discuss whether the data you are using/curating contains personally identifiable information offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
16


Reviews:
[{"review_id": "QILO_2fYNmC", "reviewer": "Reviewer_naRr", "summary": "The paper proposed a new algorithm to estimate generalized linear model with separable penalties. The algorithm is based on coordinate descent, working sets and Anderson acceleration. It is implemented in a new flexible sklearn compatible package. The numerical experiments show time improvement compared to existing methods.", "questions": "See \"Strengths and Weaknesses\"", "limitations": "See \"Strengths and Weaknesses\"", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The authors proposed a generic algorithm that uses coordinate descent, working sets and Anderson acceleration for generalized linear model with sparsity-induced separable penalties (which could be nonconvex). For a certain class of nonconvex penalties, they provided the theoretical analysis for the convergence, support identification, and local convergence rate of the proposed algorithm under certain regularity conditions. They implemented the algorithm in a flexible way that is compatible with scikit-learn package, which can also accommodate other penalties. They also perform extensive numerical experiments for both convex ($\\ell_1$, $\\ell_1-\\ell_2$) and nonconvex penalties (MCP), on huge-scale real datasets synthetic. They also conduct ablation study on working sets and Anderson acceleration. The paper is well written. \nThe following are some questions/comments:\n1. There are many assumptions and propositions in the paper, and the important propositions 10 and 13 uses multiple assumptions that may not hold in reality. Some of these assumptions like convergence to a critical point may be standard for the convergence analysis. But some of them such as Assumptions 11 and 12 might be too strong and may not hold for a general formulation, or even MCP. I think it is better to illustrate how strong these assumptions are. For example, for Proposition 13, Assumptions 1,2, and 8 and convergence towards a critical point might be standard; and piecewise quadratic holds for MCP regression; what about assumptions 9, 11 and 12 for $\\hat\\beta$? Because I do not think Assumption 11 and 12 will hold for any $\\hat\\beta$. Especially when $|\\hat\\beta_j|=\\gamma\\lambda$, Assumption 11 will not hold?\n2. Nonmonotonicity of the curves: when I read the figures for the first time, I noticed this nonmonotonicity in terms of the duality gap vs time curve. After carefully reading through the paper, I found that the authors explained it in line 219 on page 7 as well as in the last paragraph of supplement. It is probably better to point out this to reader and refer to this paragraph in some of these plot captions, so as to improve readability. Furthermore, I am wondering if it is possible to fix different levels of duality gaps and plot the median time of achieve this duality gap with the package to ensure the curve’s monotonicity?\n3. For the MCP plot (Figure 5), how is optimality condition in the y-axis of the bottom panel defined?\n", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ["I don’t know"]}, "extra_text": {}}, {"review_id": "PcGy10Mq_Pc", "reviewer": "Reviewer_uXjd", "summary": "Paper proposes a new generic algorithm for solving sparse optimization problems with very large number of variables, with convex or non-convex penalties (regularization). Their algorithm 'skglm' uses working-set strategy on the outer loop and coordinate descent with Anderson acceleration in the inner loop. Working set strategy enables scaling to millions of variables.", "questions": "- in the ablation study you note that the Anderson acceleration does not help unless working set strategy is used. This seems counterintuitive, isn't working set strategy quite orthogonal technique? Do you have any insight on this.\n- ", "limitations": "Yes.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n- broad area of application for the proposed algorithm, surely relevant to the community\n- work is placed well into the literature\n- algorithm shows promising experimental results\n- code is available and should be easy to use\n\nWeaknesses\n- I find the paper overly technical and hard to read. For example, would be good to include definitions of MCP and SCAD in the paper, as especially MCP is referred to many times. Similarly, it would be good to include the definition or intuition of Anderson acceleration in the paper itself.\n- The intuition, story, about how how authors arrived to the solution could be clearer. As the proposed algorithm seems combination of existing methods, what is the great insight, research contribution, that leads to this solution?\n- all plots are way too small and hard to understand. Consider different ways to present the results in a clearer way.\n\n\nAs a summary, the results look promising, even though I am not an expert of this particular area. By making the paper more readable, self-contained and intuitive, it would be more impactful. Consider limiting the other packages you compare to (you can include in the appendix) and highlighting the experiments where the algorithm particularly beats the competition,\n", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "B4vdDshVC09", "reviewer": "Reviewer_783Y", "summary": "The paper proposes a fast co-ordinate descent algorithm for sparse linear models with \\alpha-semi-convex penalties. The algorithm includes two key step, one is to introduce a score to rank the variables and to get a working set; the other is to user Anderson acceleration in the inner-loop. The authors give a series of theoretical guarantees. The numerical experiments show the effectiveness of the proposed algorithm. ", "questions": "See above ", "limitations": "N/A", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "I am not an expert of optimization theory and I might not be authoritative enough to evaluate the theoretical contributions of the paper. \n\nPros:\n1. the theoretical results seem rich and solid.\n\nCons:\n1.  the algorithmic contributions are not clearly articulated. It seems working set and Anderson acceleration are known ideas, right? How novel is it to apply them in the co-ordinate descent framework? Is there any other CD algorithm that uses the two ideas? Also, is the proposed score function brand-new? The authors should clarify their contributions in the two critical steps. \n2. The paper does not introduce the details about Anderson acceleration, not even including a brief summary. \n3. Line 93, why n_k is bounded by 2 | gsupp(\\beta^t)|? The authors do not explain the motivation and intuition. Can we change the factor 2? Since the size of the working set is important for the efficiency, the authors should give a clarification about the size choice. \n4. The authors should be aware of Bayesian sparse learning works, such as based on the spike-and-slab prior and horse-shoe prior. They are proven to be much more effective than L1/L1 +L2, and can  scale up to a large number of features/variables. They can provide uncertainty quantification, and well support non-convex loss functions. See the reference algorithms in [1][2]\n\n[1] Ghosh, Soumya, Jiayu Yao, and Finale Doshi-Velez. \"Model Selection in Bayesian Neural Networks via Horseshoe Priors.\" J. Mach. Learn. Res. 20.182 (2019): 1-46.\n[2] Fang, S., Zhe, S., Lee, K. C., Zhang, K., & Neville, J. (2020, November). Online bayesian sparse learning with spike and slab priors. In 2020 IEEE International Conference on Data Mining (ICDM) (pp. 142-151). IEEE.", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ["I don’t know"]}, "extra_text": {}}]

Author rebuttal:

Decision: Accept