A Highly-Efficient Group Elastic Net Algorithm with an Application to Function-On-Scalar Regression
Authors: Tobia Boschi
Matthew Reimherr
Francesca Chiaromonte
Venue: NeurIPS 2021 Poster (2021)
Keywords: Feature selection
Group Elastic Net
Function-on-Scalar regression
Functional Principal Components
Augmented Lagrangian
GWAS

Abstract: Feature Selection and Functional Data Analysis are two dynamic areas of research, with important applications in the analysis of large and complex data sets. Straddling these two areas, we propose a new highly efficient algorithm to perform Group Elastic Net with application to function-on-scalar feature selection, where a functional response is modeled against a very large number of potential scalar predictors. First, we introduce a new algorithm to solve Group Elastic Net in ultra-high dimensional settings, which exploits the sparsity structure of the Augmented Lagrangian to greatly reduce the computational burden. Next, taking advantage of the properties of Functional Principal Components, we extend our algorithm to the function-on-scalar regression framework. We use simulations to demonstrate the CPU time gains afforded by our approach compared to its best existing competitors, and present an application to data from a Genome-Wide Association Study on childhood obesity.

A Highly-Efficient Group Elastic Net Algorithm with an Application to Function-On-Scalar Regression
Tobia Boschi Department of Statistics Penn State University tub37@psu. edu
Matthew Reimherr Department of Statistics Penn State University m1r360psu. edu
Francesca Chiaromonte Department of Statistics Penn State University,
EMbeDS
Sant' Anna School of
Advanced Studies fxc110psu. edu
Abstract
Feature Selection and Functional Data Analysis are two dynamic areas of research, with important applications in the analysis of large and complex data sets. Strad- dling these two areas, we propose new highly efficient algorithm to perform Group Elastic Net with application to function-on-see on-scalar feature selection, where functional response is modeled against a very large number of potential scalar predictors. First, we introduce a new algorithm to solve Group Elastic Net in ultra- high dimensional settings, which exploits the sparsity structure of the Augmented Lagrangian to greatly reduce computational burden. Next, taking advantage of the properties Functional Principal Components, we extend our algorithm to the function-on-scall regression framework. We use simulations to demonstrate the CPU time gains afforded by our approach compared its best existing competitors, and present application data from Genome Wide Association Study on childhood obesity.
Introduction
As problems involving very large and potentially structured data become ever more ubiquitous, attention is being devoted to the integration of approaches and techniques from the areas of Feature Selection and Functional Data Analysis (FDA). Indeed, more and more regression applications comprise large number number of variables some of which scalar and some of which suitable for a functional representation, such as longitudinal measurements or biomedical images (Sorensen et al. [2033] (llah and Finch 2013] Cremona et al. [2019) A great deal of recent work has been concerned with feature selection in these applications. Matsui and Konishi (201 Gertheiss et al. (2013); Fan al. (2015) study the case where the response scalar and the features are functional. Chen et al. (2016) 6 Fan and Reimherr (2016); [Babber et al. )); Parodi et al. (2018); Mirshani and Reimherr (2019) tackle the called function-on-scalar case, where the response is functional and the features are scalar focusing on settings in which the number of features is bigger than the number of observations. However, recent developments in optimization have demonstrated that substantial computational gains can still be made when the number of features is massive (e.g ~16)). In this article, we present Functional Group Elastic Net (fgen), novel and highly efficient method settings where
to solve the function-on-scan feature selection problem in Altra-hig-ddimension the number of features indeed massive and much larger than the number of observations. The ability to solve these problems with lower computational burden is increasingly critical. Given the complex, noisy nature much contemporary data, changing some aspects of their pre-processing or some of the tuning parameters involved in the analysis can lead to completely different results (Krawezyk and Cano 2018 Murdoch et al. 2019). For this reason, repeating an analysis multiple
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

times (e.g., with different choices of data preprocessing pipelines, or to tune certain meta- parameters) paramount to capture significant signals and ensure the stability of outcomes (Yu and Kumbier 2020). Substantial reductions in computational burden enable such repetition, allowing scientists and practitioners to conduct truly meaningful and reproducible analyses.
Group Elastic Net incorporates the group structure (Yuan and Lin 2006) and the Elastic Net penalty (Zou and Hastie, 2005) into a penalized regression framework. The former allows one represent each feature (or component) by group of variables. The latter induces sparsity and regularizes the estimates. We consider the case where all groups have the same size k. The minimization problem is formulated as follows:
P min (1/2) (XB + A1 E1BBBll (A22222(22/2) (22/2) Σ B.ll
(1)
B
Let p be the number of features, the number of statistical units, and ll the 12 norm for matrices, i.e the Frobenius norm, and vectors. Then, X E R"XP is the design matrix (that we assume to have standardized columns), Y E T the response matrix, and B E TR the coefficient matrix. In other words, (I) describes a sparse multi-task model where the response and each of the features are represented by a group of & coefficients (Zhang and Yang 2018). Throughout this article, we follow the notation in Johnson et al. (2014) and we use the subscripts and (i) to indicate the i-th row and the i-th column of a matrix, respectively. Thus, B; E IR are the coefficient values associated with the i-th group, and X()) E 12d are the observed values relative the i-th feature. Before proceeding. note that ① can be expressed as
min (h(XB) + +B)). B
where = = - Y12 is the least-squares loss function and TTB) - Ei (B,)
_ is the Group Elastic Net penalty function. The first term in π not differentiable and creates sparsity at the group level, i.e, if component selected, then all its coefficients are selected and vice-versa. The second term is Ridge-type penalty which reduces model complexity and tries to control variance inflation due to multicollinearity feature selection models are indeed known to be less effective and not reliable in scenarios characterized by very high collinearity among features (Katrutsa and Strijov 2015). a1 and X2 are penalty parameters > 0 and control the weight of the two penalties with respect to the least square loss.
A1 E B (>e/2 [
To solve M.. we develop a new Semi-smooth Newton Augmented Lagrangian (SsNAL) algorithm. We then extend it to the function-on -oniiiol----ccllrr regression framework by means of Functional Principal Components (FPC) (James et al.] 22000] [Chiou al.) [004]] Hall and Hosseini-Nasab) (006). S&NAL exploits the sparsity induced by the augmented Lagrangian second order information to guarantee a super-linear convergence and greatly reduce the computational cost. This methodology, first introduced by Tomioka and Sugiyama (2009) and Tomioka et al. (201T), has been recently used in several applications, e.g., to regular Lasso (Li et al.) 2018), constrained Lasso (Deng and So 2019), and Elastic-Net (Boschi et al. 2020). However, incorporating the group structure significantly increases the dimension of the problem. Indeed, (I) is not separable and the optimization must be carried out jointly across the coordinates of the outcome. Therefore, considering new group penalty while preserving the efficiency of the method requires that we carefully redefine a set of all-new mathematical operators and the theory behind them.
We implemented an efficient version of fgen in python and benchmarked against the two best Group Elastic Net solvers we found in the literature: the python package sklearn (Pedregosa) jet al. T) and the R package glmnet (Friedman et al.) 2010), which is written in fortran Both of these solvers implement a highly optimized coordinate descent algorithm (Friedman et al.
2010
Breheny and Huang 2015) and outperform competitors such us FISTA (Beck and Teboulle
2009 Bonnefoy et al.) 2015), ADMM (Deng et al. 2013 Zhu 2017), and proximal gradient (Chen ct al. 2010) by at least one order of magnitude in terms of CPU time. Our simulation results
demonstrate that in sparse scenarios fgen is at least 3 times faster than glmnet and more than 10 times faster than skl earn. We also applied fgen to the Intervention Nurses Start Infants Growing on Healthy Trajectories (INSIGHT) study (Paul et al.) 2014),, which investigates risk factors for childhood obesity. Specifically, we examined the association between hundreds of thousand of Single Nucleotide Polymorphisms (SNPs) and growth eurves, which represent functional outcome.
The remainder of the article organized as follows. In Section we describe the Group Elastic
Net problem and introduce some preliminary results. In Section 3 we present our new methodology
2

---

and illustrate how to extend it to a function-on-se -n-alr feature selection problem. In Section 4 we investigate the performance of our method on simulated data and apply data from INSIGHT. In Section 5 we provide final remarks and discuss future developments. Proofs of theoretical results and additional simulations are included in the Supplement. The fgen code available at https://gg gitggtttuu..... www//oobbbabschh/eoeeeee.e..
2
Preliminaries
In this section we define the Group Elastic Net problem and we introduce some results related to Fenchel conjugate functions and proximal operators, which are essential tools in our developments.
2.1 Fenchel conjugate function and proximal operator of TTB)
Fenchel conjugate functions (Fenceel 1949) allow one to more readily define the dual problem
(Boyd and Vandenberghe 2004) of (P) which is called the primal problem. Let X' S R² be convex set and f X R. Then, the conjugate function of f X - E defined as
(=) suP=e> ((2, the inner product. i.e. the dot product. If X' C P i.e. if is matrix, the definition is still valid but ( is the Frobenius inner product. In our first proposition we provide closed form solution for the Group Elastic Net penalty conjugate function (see Supplemental Section A.1 for a proof).
- where X = {2 suPzex ((2,z) - f(z)) < <0 (..)) indicates
Proposition 1. Given Z E RRxx,, the conjugate function of π has the form
x (2)
- =2x))
Σ (110112~1144)
where [ ]+ is the positive part operator; [s]+ = 8 if > 0 and 0 otherwise.
Note that x² (Z) is continuous differentiable function. This is more general results than the one presented in Li et al. (2018); Boschi et al. 2020). because we extend the definition of T to the case where 2 is matrix and not just vector. In the simple scenario where k 1, i.e. when every group consists of just one variable, we obtain again the conjugate function of the standard Elastic Net penalty. Notably, starting from non-separable objective function, derive x² which does separate. As we will see in Section 3.1 this allows one to induce a new level of sparsity in the Lagrangian problem which actually key for the massive computational advantage offered by fgen.
Proximal operators (Rockafellar 1976a b) are fundamental in many optimization algorithms. Given a lower semi-continune convex function f IP + R. the proximal operator of f at x with parameter >> 0 is denoted prox,, R - R and defined prox, (z) arg min, (f( (2a)-1 =112). If f Rpxk + IR, i.e. if i is a matrix, then prox,; IR + ggxk and ||-11a is the Frobenius norm. Parikh al. (2014) and [eck (2017) (Chapter 6) provide numerous examples and properties. Combining their results, one can easily find the form of the proximal operator of (B) provided in our second proposition (see Supplemental Section A.2 for a proof).
Proposition 2. The proximal operator of T (BB) is: prox, n (B) (prox,, (Bi).... (Bro),((Bi.... where
prox, (Bp))
prox (B,) = (1 (12 1 |B,11 ad1 B,
(3)
Note that prox, (B) : Rxpx + Rpxk, To implement fgen, one also needs the proximal operator of T²" which can be obtained through the Moreau decomposition:
propppoxxxxo (x) + (=0/=+/o(1/0) 0> 0.
(4)
2.2 Dual formulation and Augmented Lagrangian
Here we introduce the dual Group Elastic Net problem and its augmented Lagrangian. From Boyd
and Vandenberghe (2004), a possible dual formulation of (P is
min (v)) s.t. x =
(D)
V.2
3

---

Algorithm Semi-smooth Augmented Lagrangian (SsNAL) method
Augmented Lagrangian method
Semi-smooth Newton method for 7
Start from the initial values V - B 0
solve 7 and find (V 2
while not converged do
while not converged do
(I) Given B", find V*+1 and = which ap- proximately solve the inner subproblem
(1) Find the descent direction solving exactly
by conjugate grudient the linear system
(Vk+1, z
ss min Lo (V, Z B²) (7) V,2
8 vec")) = vee (Vé(V*))
(9)
(2) Line search (Li et all 2018); choose P
(2) Update the Lagrangian multiplier B and the parameter o:
(0, 1/2) and reduce the step size until
1 (V* + + D') S (VV) + us (VN(V>), D))
Bh+ B² -a -x
z2
(3) Update V: V:1+1 = V D'D
(8)
t 00
(4) Update z prox ( - x
end while
end while
where V E R and Z E Rp*k the dual variables matrices. In particular V,, Zj R are the dual variables associated with the i-th group. h" and π are the Fenchel conjugate functions of h and t, respectively. Specifically, A (VV (1/2) | + (Y, V) (Dünner et al. 2016) and x" xZ) is given in Proposition We can now define the augmented Lagrangian function and the Ktarush-Kuhin Kahn-Tucker (KKT) system associated with (). The augmented Lagrangian is given by
C. (V, +x)) Σ (B. VT X(0 + Z.) + (a/2) Σ C Zill (5)
where a > 0. both the primal variable and the Lagrangian multiplier which penalizes the constraints' violations. The KKT system is given by the following three equations:
V XB= 0, 0, 0= 0= (Z) B=0,
X 7
(6)
Note that (V) V + Y. A closed form of (Z) is not essential for SsNAL method. The KKT equations will be useful to determine the convergence of our algorithm, since the set (** B*) solves the KKT (6) if and only if (V*, 2²) and 8 are the optimal solutions of (D) and
respectively (Boyd and Vandenberghe 2004).
3
Methodology
In this section we present our new methodology. First, introduce SsNAL algorithm to solve the Group Elastic Net problem. Next, we illustrate how to extend it to the function cecction-onscaat on-scalar regression framework. Finally, we describe how implement solution path over different values of Å₁.
3.1 SsNAL method
The S&NAL method is summarized in Algorithm I consists an Augmented Lagrangian method characterized by inner subproblem. The subproblem is solved with Semi-smooth Newton method which exploits the sparsity of the augmented Lagrangian second order information and greatly reduces computational costs. We now provide the details of its implementation and some important theoretical results. From Rookafellar ( (1976)). one can find the optimal solution of (D) by solving the Augmented Lagrangian method described in Algorithm[] [ The essential part of the algorithm is the subproblem
7 As described in Li et al. (018), an approximate solution (V, 2) for a given B can be found as
= arg min L. (V 2, B), V
= arg min C. V 2
(10)
With slight abuse of notation, we indicate by L. (VIZ, B) the function L, (V, B) where the parameter Z and B fixed. Similarly for L, V, B). Our third proposition provides explicit forms for C. (V Z, B) and Z (see Supplemental Section A.3 for proof).

---

Proposition 3. Define (VV = Lo (V Z, B). Then, for the Group Elastic Net problem we have
(a) (v) =((() 1 σ2
Σ
p
prox,x
oV (Bi-oV"X, (1))
1
p
|B:112
20
20 i=1
(11)
(b) 2 - ppoxx+ /a (B/o - V)),
where po (B/o X = (prox,- (Bi/o - V*Xa)), po (Bp/a V
Z has closed form. To find V one has to minimize to or, equivalently, find the solution of V 0. Note that v is continuous and differentiable, and thus Vi is well defined.
To solve the subproblem
7 we propose the Semi-smooth Newton method in Algorithm | and
according to the rule in Proposition 3 and V by minimizing & through one Newton step. The main computational cost is solving the linear system ( This leads to our next crucial result (see Supplemental Section A.4 for a proof).
Z updated iteratively
Theorem 1. Let T = B - oXTV, X = X @ ly (the nk x pk Kronecker product between X and the k x k identity matrix). O2rb be the generalized Hessian of v, and 0 prox, be the Clarke sub-differential of prox,= (Clarke) 1990). Then we have
(V)(()) = =V VY X ox (T) (())&(V))
Ink prox (T)
(12)
Pi
Moreover, let Q E Rokkpp be the block- diagonal matrix Q =
where each P, is
Pp.
squared is x k matrix defined
(1+0 (1- ||T. aA1) - 1 OATTTT))
1T > σd1
(ii)
(13)
0
o.w.
Then Q E ô prox,- (T) and 8 (V) vec(D) = (Ink + =xQXT) vec(D) for every D Bx in the domain of V - where vec( D) E Rin obtained by stacking all the columns of D.
Note that, while in Li et al. 2018): Deng and So (2019); Boschi et al. (2020) V and 8² are vector and a matrix, respectively. here the dimensions of these operators increase due to the group nature of the problem. In particular, Vei becomes a matrix and 82 higher order tensor which we express as an nk X nık matrix by stacking its dimensions. Moreover, Q is not simple diagonal matrix as in the previous SsNAL algorithms, but is now characterized by blocks associated to the different groups of variables.
Theorem is critical for preserving the efficiency of fgen, while integrating groups into the problem. First, it states that solving (9) is equivalent solving (Ink + xxQX") vec( vec (V)( Second, the form of Q still allows one to induce sparsity in the linear system and drastically reduce the computational cost. Indeed, let J {i : [22 2 σ and let [J] be the cardinality of J. Then the linear system (9) is equivalent to
+ Tok+oxxQu0
Qx vec(D) vec (VV(V)).
(14)
Ls. with Xy gaxr being the sub-matrix of X restricted to the columns in J. In addition, Q3 E Brrxrkk is the block-diagonal matrix formed by all the P. such that E J. Using the Cholesky factorization (the generalized Hessian is positive
X, Here, X- lknk is defined as x==
semidefinite) the total cost of solving the linear system reduces from O (nk (n² + p² + ni)) to O (nk33 + + nr)). This includes computing X which is O (nrk3 (n + r)), and the Cholesky factorization, which ( Because the sparsity induced by the Group Elastic Net penalty, is usually much smaller than p implying substantial computational gain. Even when p is very large (~ 10%) one can still solve the linear system efficiently, as long as the dimension k of each group is relatively small (< 102). Furthermore, if r < n, which is often the case when
5

---

Algorithm 2 Functional Group Elastic Net method
(1) Perform FPC of B and find the first & basis components (76, 7w)) with their eigenvalues (p pk)
(2) Find the first & FPC scores for each response function Vy: Y; = (i)i,,)........ (Di, 7h)2)
(3) Using Y as response matrix, apply SsNAL to solve ① and find the coefficient scores estimates B.
(4) Project B into the FPC basis to find the coefficient curve estimates: B, Σラ= B(,,5)
the solution of the Group Elastic Net problem is sparse, one can factorize an rk x rk (instead of nk x nk) matrix using the Sherman-Morrison Woodddbryy formula (Van Loan and Golub 1983):
(Ii ( + +xx3xx3333
Inn -Xa ((095)'+x3383)
X3.
(15)
The total cost is further reduced from O (n1=3(n² 22 nr)) O (rk- + r2 + + 1)), in- cluding the computation of Q77 which can be done with a cost of ((rk2) by inverting each one of the P. blocks independently. Finally, if in the first iterations of the algorithm and r are both larger than 10², one can solve ( approximately using the conjugate grudient method (Polyak) 1969)
To determine the convergence of the Augmented Lagrangian and the Semi-smooth Newton methods, we check the residuals of the third and first KKT in (6) respectively, i.c.
E E+
VT (1) + 2:112
res(kk)) Di=1 1 + Y. - X,BI1, Σ |2
res(kkts)
(16)
1+ 11
VI 1
Taking the 12-norm of the KKT residuals, normalizing them and using them to assess convergence a common procedure the literature (Li et al., 2018 Deng and So. 2019). Both methods have super-linear convergence rate. Accordingly, the convergence rate of the entire algorithm, which the sum of the convergence rate of the two sub-problems (Tomioka and Sugiyama 2009), still super-linear. Thus, we show in Section 4 fgen typically converges in very few iterations. The convergence analysis here follows directly from that in Tomioka et al. (20TT) and Boschi et al. (2020), where SsNAL convergence is proved for the standard Elastic Net. The proof leverages results in Rockafellar 1976a b): Luque [984);; [L al. (2018) and the fact that π is a continuous differentiable function (wiich is also true for the Group Elastic Net).
3.2 Extension to function-on-scalar regression
We now extend fgen to the function-on-scalas features selection problem. In function-on scalar regression, a functional response is regressed on a set of scalar predictors. Assuming the response belongs the Hilbert Space [[a b)), the optimization problem ①
becomes
min
(1/2) VI + A1 Σ |B + (A2/2) Σ | B
(17)
1=1
where y and B are functional objects with and p rows, respectively. Each row )1 is response function and each row B; is coefficient function. The squared L of function f is II.
(f.,,,, where the inner product between two functions ff and g is (f.g) = So fg.
Applying SsNAL directly (17) is not straightforward and would substantially hinder its efficiency. First, the definition of conjugate functions and proximal operators in functional spaces would require theoretical background. Second, and perhaps most important from practical standpoint, computing integrals is much more expensive than computing euclidean norms. For these reasons, in Algorithm 2 we take advantage of Functional Principal Components (FPC) (Horváth and Kokoszka (0017) to solve an optimization problem of the same type as (I) which (7). In particular, we build response matrix Y. where each group is formed by the first l FPC scores of the the function V, (Fan and Reimherr (016)). The level of approximation of fgen thus depends on the number of FPC scores k, i.e the dimension of each group. Indeed, given function f and its FPC basis (%a 10c 1. we have
2012 Kokoszka and Reimherr
fact a very close approximation
11 Σ IIF rheell
(18)

---

This property- which is true for every orthonormal basis system -plays crucial role in the extension of our SsNAL approach to the function-on catioon-o--ccllat regression, since it allows one to approximate the 1.2 function norm with the standard I2 matrix norm, Consequently, one can use the FPC scores to construct the response matrix Y and the coefficients matrix B in (T starting from the response functions Y' and the coefficient functions B. The number of selected FPC scores determines the dimension &or of each group in the Group Elastic Net problem. In many applications just few FPC scores allow one to obtain a very close approximation of the original functions. Indeed, among the many orthonormal bases one could envision, FPC has the advantage of being the most parsimonious allowing one to reconstruct the response curves using fewer coefficients than any other orthonormal basis. In scenarios investigated by simulation in Section 4 k = 5 is sufficient to capture more than the 99% of the -oono This produces an almost perfect approximation of (17) while fully preserving f fgen efficiency.
3.3 Solution path implementation
To evaluate different values of the penalty parameter '1 we implement an efficient solution path search. We compute the solution for decreasing sequence sequence of λ starting from Xmar (xx; TYll which selects O active features. When we move to the next Y value, we use the solution obtained the previous value for initialization (warm start). The two consecutive solutions tend to be close. and fgen converges in very few iterations usually just one. We also allow the user to specify a maximum number of selected features; when this number is reached the path search is stopped, further reducing computation.
To guide the choice of (2) we propose two quantitative criteria: k-fold Cross Validation (cv) and an Extended Bayesian Information Criterion (E-bic) (Chen and Chen 2012). which modifies the standard BIC also include the number of features p. In symbols, we have
e-bic B) = k log (rss(B)/(mkk)) + (kv) (log(nk) + log /n
(19)
where rss rss(B) is the residual sum of squares associated with the solution B, and v are the Group Elastic Net degrees of freedom. From Tibshirani al. 2012))
tr(X, (X3X, >2I,)) -xx))))x3) Note that cv can be very computationally expensive because it requires to run fgen multiple times for each value of Åι and x2 under consideration. In contrast, e-bic can be computed directly from the original solution. Before evaluating both criteria, we de-bias the fgen estimates following the approach suggested by Belloni et al. (2014); Zhao al. (2077). First, we fgen, then, we fit a standard least squares on the selected features. In the next section, following standard practice in the literature e.g.. Friedman al. (2010); Pedregoss et al. (201T)-
rewrite \₁ and X2 as A = Cy X and x2- = (1-a)cxAthar Amaa with cx E (0, 1] and a E (0, ex determines the reduction with respect of Xmaz or controls the relative weight of the two penalties.
4 Simulation study and INSIGHT data
In this section we use synthetic data to illustrate the computational efficiency of fgen, and apply our new method to Genome Wide Association Study (GWAS) on childhood obesity. In the simulations, we benchmark fgen against the two best Group Elastic Net solvers we found in the literature: the python package sklearn and R package glmnet, which is written in fortran. Other functional-od -onsclll feature selection methods, such as the ones proposed by Barber et al. (0017); Parodi et al. (2018); MMrshani and Reimherr (2019), have computational burden more than two order of magnitudes larger than fgen and could not complete instances with p > 10
4.1 Simulation results
We generate synthetic data as follows. The entries of the design matrix X E Raxp are each drawn independently from a standard normal distribution. The response curves are created as y XB + e. B contains Po non-zero curves. These and the errors c are generated from a O mean Gaussian process with a Matérn covariance function (Cressie and Huang] 1999) of the form
C(t, s) = w2 (f(m)2--1)~ -1 (10- (20) s1) " K, (2)
(20)
where K2 is a modified Bessel function. In particular, we set the point-wise variance 2² and the range | = 0.25 (this determines how fast the curves dependency decays). The smooth parameter
7

---

Table I: b and report CPU time in seconds for fgen, sklearn and glmnet, respectively. For fgen we also report the number of iterations parenthesis. is the number of selected features, I is the range parameter of the Matern process used to generate the coefficients.
a = 0.8,
7
0.25
n=500
n= 1000
n=50000
pi Po
Cx
a
r
b
a
C
0.8
2
0.1(3)
0.5
0.3
(.1
0.9
0.6
0.5
9.1
3.2
5
0.4
10
0.2(4)
0.5
0.3
4
0.3(3)
1.1
0.6
1.1(3)
10.9
3.1
2(10"); 10
21
0.4(4)
0.3(4)
1.2(3)
10.4
3,4
0.8
2
0.2
0.9
0.5
0.2(2)
1.9
1.2
1.3(2)
18.8
5.4
10
0.4
10
(43
0.8
0.5
0.4(3)
2.1
1.0
6
2.0(3)
22.3
4.9
0.2
21
(.9(4)
0.9
0.5
8
0.5(4)
2.2
1.2
10
1.9(3)
23.2
5.3
0.8
(23(2)
4.5
1.5
0.7(2)
9.9
3.3
2.6(2)
97.9
17.2
5
0.6
29
0.6(2)
4.7
1.6
17
0.7(2)
9.6
3.3
12
3.0(2)
85.7
19.8
10"; 10
0.4
486
(3.7(3)
4.8
0.2
53
1.4(2)
10.2
3.1
42
4.1(2)
99.7
17.8
0.8
0.6(2)
10.2
2.6
1.4(2)
22.6
5.2
3.5
175.2
27.9
10
29
(29(2)
10
2.4
17
(.5
19.6
4.8
3.8(2)
183.7
27.6
0.4
486
79.2(3)
9.1
2.4
53
1.8(2)
20.2
5.3
42
6.7(2)
190.7
28.0
0.5,
= 0.25
n=500
n= 1000
n=500
p. po
k
Cx
a
b
c
r
a
b
c
a
b
C
0.8
2
0.1
1.0
0.6
0.2(2)
1.8
1.3
1
0.9
18.3
5.2
2(10");
(32(3)
(.2(2)
0.9
7
(.9
20.6
0.2
73
1.1(3)
0.9
0.5
0.4(3)
2.0
1.3
(39
22.7
5.2
0.8,
= 0.10
n=500
n= 1000
n= 5000
p: Po
CA
r
b
r
a
C
0.8
3
0.2(3)
0.9
0.6
0.2
1.8
1.3
0.9
19.7
5.1
2(10*);
10
0.4
(4)
0.6
0.4
2.6
10
1.9
19.8
5.2
0.2
1.5(4)
1.2
0.5(4)
2.4
1.2
2.0(3)
22.2
4.0
responses
errors
coefficients
Lup estimated
T
49
?
Figure 1: Curves related to a simulation scenario described in Table ( row, 500). The left and center panels display sample of 10 response curves and IO linear model errors terms, respectively. The right panel depicts the true non-zeros coefficients curves (dashed lines) and the fgen de-biased estimates (solid lines).
equal 3.5 for and to 1.5 for c, i.e the errors are rougher than the coefficients. Each curve is sampled at 1000 evenly spaced points between 0 and 1. Figure shows instances of the response and error curves, and e, and the true non-zero coefficient curves in along with their de-bia estimates produced by fgen (the underlying simulation parameters are those in Table D second
row, n 500). all scenarios, fgen is rum with both the tolerances in (T6) set to 10 (we set the same tolerance for sklearn and glmnet) and and increase by a factor of every iteration. we start from smaller values of o, the algorithm needs more iterations to converge, while if o is too large, fgen does not converge to the optimal solution. We set D cyX and X2 = (1 a)cxAmm where cx E (0. 1), a E (0, 1). and mar max, [Xx, (X,) TYI. Note that for glmnet and sklearn need divide A mar by since both solvers divide the least squares loss in ( by the number of observations.
(2) set to 0.2. We start from の0 po/p
Table I reports CPU times for fgen, sklearn and glmnet under different simulation settings. fgen the fastest solver in almost every instance. When both n and p are large and the solution is sparse, fgen is approximately 6 times faster than glmnet and more than 30 times faster than sklearn. Note

---

response
10-fold
e-bic
0.18
==
9
m0
0.1.4
0.00608
D =0
7
0.10
10
0.00602
0.06
10
20
30
0.99
0.9
0.8
0.7 0.99
0.9
0.8
0.7
weeks
C)
Cy
Figure 2: Plots related to the INSIGHT study. The left panel displays the growth curves. The center and right panels depict values of the 10-fold co Mean Squared Error and the e-bic, respectively, against CA These are obtained from fgen run with 3 different values of a; 0.8 (green line), 0.6 (red line), 0.4 (blue line). fgen estimates are de-biased prior computing both criteria.
that the super-linear convergence rate allows fgen to converge in very few iterations (no more than 4 in all cases). The CPU time increases with k for all solvers. However, A: = 5 already captures more than the 99% of the -orm in all the scenarios considered. we decrease a from 0.8 to 0.5, fgen need even fewer iterations to converge, increasing its computational gain with respect to the competitors. Considering rougher coefficients (created with a Matern process with range parameter I 1= 0.1) does not affect the relative performance of the algorithms. The instance with an active set of 486 features is the only one where fgen performs worse than its competitors. As expected, in the presence of non-sparse solutions fgen looses some of its efficiency. However, tune the penalty parameters in practice, one evaluates a sequence of сΛ values starting from very sparse solutions. In the first steps of the solution path, fgen exploits sparsity and is very efficient. In the following steps, still converges very quickly thanks to the warm-slarl approach described in Section 33. In Supplemental Table B.2 compare the solution path computing time. fgen outperforms the other solvers in every scenario, being approximately 2 times faster than glmnet and from 10 to more than 30 times faster than sklearn. Finally, gauge uncertainty in our CPU time evaluations, we replicated a subset of the instances explored in Table independent times. Mean CPU times and standard errors over such replicates are reported in Supplemental Table B.I Results agree with those obtained considering just one replication. Furthermore, one can notice that fgen has also smaller variability CPU time when n 5000.
Taking into account that glmnet (written in fortran) and sklearn are highly optimized packages, the results above provide strong evidence in support of method. We also tracked prediction performance for all methods and in all simulation settings considered, but we do not report them here since all three solvers solve the same convex minimization problem and therefore converged the same solution in all settings.
4.2 INSIGHT study
Here, apply fgen to data from the Intervention Nurses Start Infants Growing Healthy Trajec- tories (INSIGHT) study (Paul et al.) 2014).. In particular, we focus on data collected to investigate genetic variants that may affect the risk of childhood obesity. As the prevalence of obesity increases also among children, examining possible causes and risk factors has become an essential public health concern. INSIGHT provides genome-wide Single Nucleotide Polymorphisms (SNPs) information for a cohort of very young children, along with longitudinal information on their growth. Selecting SNPs that may affect growth is thus GWAS Genome-Win Association Study) where the outcome is a growth curve. In recent years, many GWASs have identified SNPs strongly associated with obesity phenotypes (Locke et al. 2015). Before proceeding with the analysis, we point out that due high feature collinearity, low -ggallto-noise ratio, and ultra-high dimensionality. GWAS data are very hard examine and users should be very careful in interpreting results; e.g., selected SNPs may just be proxies for other causal SNPs in their vicinity. While being well aware of all its complexities and potential pitfalls, our main aim presenting GWAS analysis is to show the efficiency and broad applicability of fgen. Our functional outcome captures the evolution of weight/heighl ratios (Daniels et al.] 2015) measured at birth and at 4, 16, 28, and weeks for total of n 210 children. The growth curves shown in the left panel of Figure 2 are fitted as in Craig et al. (2019) using

---

Principal Analysis by Conditional Estimation (Chen et al., 2017). Building smooth curve for each child allows one to capture information along the entire time domain and at the same time de-noise and mitigate the effect of outlyng//aoomalon are available upon request at affibGaaP using the access number (as000149.vvl The growth curves are based on privacy protected data and cannot be made publicly available.
raw measurements. The SNPs collected in INSIGHT
Craig al. 2019) used flame (Parodi et al. 2018) solve the function-on-seit feature selection problem. To do so, they had to reduce the analysis from p 342325 to 10000 SNPs with various preliminary screening steps. The computational efficiency of fgen allows us to inspect all the 342325 SNPs simultaneously. The center and left panels of Figure 2 display 10-fold cV and e-bic for different values of a and ex. Both criteria identify just one (cx- 0.99) dominant SNP, rs79187646. Without drawing any strong domain conclusion, we remark that the selected SNP appears to be relevant in the literature. Notably, this the same SNP selected by [Boschi et al. (2020), where the same SNP data were associated to BMI at age 3 a scalar response. Also notably, this dominant SNP was not among the SNPs identified by Craig et al. (2019). Based on the U.S. National Library of Medicine rs79187646 is located in NTM, a well known gene. According to the NHGRI-EBI GWAS Catalog many GWASs, including two recent studies (Kichaev et al. 2019 Pulit et al. 2019), have connected NTM to body mass, food addiction, intake of sweet substances and other obesity-related traits.
5 Conclusions
In this article we proposed new Function Group Elastic Net method (fgen) to solve the function-on- scalar feature selection problem. Our proposal starts with the development of a novel, highly-efficient SsNAL algorithm solve the Group Elastic Net which is then extended to the function-on-scall regression framework using Functional Principal Components representation. Though we could rely on critical prior results (Tomioka and Sugiyama 2009: Li al., 2018 Boschi et al. 2020), in order to integrate the group structure into SsNAL, we had to tackle more complex mathematical operators and redefine the theoretical foundation. Our simulations show a substantial reduction in CPU time with respect to the best existing Group Elastic Net solvers. Finally, we applied fgen to a GWAS study detecting SNP that may affect obesity risk in children.
The current version of fgen is limited to the case where each group has the same size k and to the funetional-on-seallr feature selection problem. In the future, we plan to further extend work investigating more complex optimization problems (e.g., allowing each group to have a different size) and functional regression frameworks. In particular, we aim to adapt our methodology the function-on-functee feature selection scenario, where both the response and the predictors can be represented functional curves.
Acknowledgments and Disclosure of Funding
We thank Kateryna Makova and her laboratory Penn State for access to the INSIGHT data, Ana Kenney for help with the data and useful discussions. The work of Tobia Boschi and Francesca Chiaromonte was partially supported by the Huck Institutes of the Life Sciences at Penn State, the work of Matthew Reimher was partially supported by the Grant NSF SES-1853209.
References
Barber, R. F., M. Reimherr, T. Schill, et al. (2017). The function-on-ssallrt lasso with applications to longitudinal gwas. Electronic Journal of Statistics 11(1), 1351-1389,
Beck, A. (2017). First-order methods in optimization, Volume 25. SIAM.
Beck, A. and M. Teboulle (2009). A fast iterative serinkage-thresooldin algorithm for linear inverse problems. SIAM journal on imaging sciences 2(1), 183-202.
Belloni, A., V. Chernozhukov, and C. Hansen (2014). Inference on treatment effects after selection among nigh-dimensional controls. The Review of Economic Studies 81(,), 608-650
Bonnefoy, A., V. Emiya, L. Ralaivola, and R. Gribonval (2015). Dynamic screening: Accelerating first-order algorithms for the lasso and gassp- IEEE Transactions on Signal Processing 63(19), 5121-5132.
10

---

Boschi, T., M. Reimherr, and F. Chiaromonte (2020). An efficient semi-smooth newton augmented lagrangian method for elastic net. arXiv preprint 20000.0000..3970.
Boyd, S. and L. Vandenberghe (2004). Convex optimization. Cambridge university press.
Breheny, P. and J. Huang (2015). Group descent algorithms for nonconvex penalized linear and logistic regression models with grouped predictors. Statisties and computing 25(2), 173-187.
Chen, and Z. Chen (2012). Extended bic for small-n-large-p sparse glm. Statistica Sinica, 555-574.
Chen, K., X. Zhang, A. Petersen, and H.-G. Müller (2017). Quantifying infiite-diimensione data: Functional data analysis in action. Statistics in Biosciences 9(2), 582-604.
Chen, X., Q. Lin, S. Kim, J. G. Carbonell, and E. P. Xing (2010). An efficient proximal gradient method for general structured sparse learning. stat 1050.
Chen, J. Goldsmith, and R. T. Ogden (2016). Variable selection function-on-scalar regression. Stat 5(1), 88-101.
Chiou, J.-M.. H.-G. Müller, and J.-L. Wang (2004). Functional response models. Statistica Sinica, 675-693.
Clarke, F. H. (1990). Optimization and nonsmooth analyvis, Volume 5. Siam.
Craig, S. A. M. Kenney, J. Lin, M. Paul, L. L. Birch, J. Savage, M. E. Marini, F. Chiaromonte, M. L. Reimherr, and K. D. Makova (2019). Polygenic risk score based on weight gain trajectories is strong predictor of childhood obesity. bioRxiv, 606277.
Cremona, M. A., H. Xu, K. D. Makova, M. Reimherr, Chiaromonte, and P. Madrigal (2019). Functional data analysis for computational biology. Bioinformatics (Oxford, England) 35(17), 3211.
Cressie, N. and H.-C. Huang (1999). Classes of nonseparable, spatio-temporal stationary covariance functions. Journal of the American Statistical Association 94(448), 1330-1339,
Daniels. S. R., S. G. Hassink, et al. (2015). The role of the pediatrician in primary prevention of obesity. Pediatrics 136(1), e275-e292.
Deng, W., W. Yin, and Y. Zhang (2013). Group sparse optimization by alternating direction method. In Wavelets and Sparsity XV, Volume 8858, pp. 88580R. International Society for Optics and Photonics.
Deng, Z. and A. M.-C. So (2019), An efficient augmented lagrangian based method for constrained lasso. urXiv preprint 1rrXiv1903..00
Dünner, C., S. Forte, M. Takáč, and M. Jaggi (2016). Primal-dual rates and certificates. arXiv preprint
Fan, Y., G. M. James, P. Radchenko, et al. (2015). Functional additive regression. The Annals of Statistics 43(5), 2296-2325.
Fan. Z. and M. Reimherr (2016). High- digh-dimmnsional function-on function-on-scale regression via the alternating direction method of multipliers. 2016 3rd International Conference on Information Science and Control Engineering (ICISCE), pp. 397-399. IEFF.
Fenchel, W. (1949). On conjugate convex functions. Canadian Journal of Mathematics I(1), 73-77.
Friedman, J., T. Hastie, and R. Tibshirani (2010). Regularization paths for generalized linear models via coordinate descent. Journal of statistical software 33(1), 1.
Gertheiss, J., A. Maity, and A.-M. Staicu (2013). Variable selection in generalized functional linear models. Stat 2(1), 86-101
Hall, P. and M. Hosseini-Nasab (2006). On properties of functional principal components analysis.
Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(1), I), 109-126.

---

Horváth, L. and P. Kokoszka (2012). Inference for functional data with applications, Volume 200. Springer Science Business Media.
James, G. M., T. J. Hastie, and C. A. Sugar (2000). Principal component models for sparse functional data. Biometrika 87(3), 587-602.
Johnson, R. A., D. W. Wichern, et al. (2014). Applied multivariate statistical analysis, Volume 6. Pearson London, UK:
Katrutsa, A. and V. Strijov (2015). Stress test procedure for feature selection algorithms. Chemomet- rics and Intelligent Laboratory Systems 142, 172-183.
Kichaev, G., G. Bhatia, P.-R. Loh, S. Gazal, K. Burch, M. K. Freund, A. Schoech, B. Pasaniuc, and A. L. Price (2019). Leveraging polygenic functional enrichment to improve gwas power. The American Journal of Human Genetics 104(1), 65-75.
Kokoszka, P. and M. Reimherr (2017). Introduction functional data analysis. CRC Press.
Krawezyk, B. and A. Cano (2018). Online ensemble learning with abstaining classifiers for drifting and noisy data streams. Applied Soft Computing 68, 677-692.
Li, X., D. Sun, and K.-C. Toh (2018). A highly efficient semismooth newton augmented lagrangian method for solving lasso problems. SIAM Journal Optimization 28(1), 433-458.
Locke, A. E., B. Kahali, S. Berndt, A. E. Justice, T. H. Pers, F. Day, C. Powell, S. Vedantam, M. L. Buchkovich, J. Yang, et al. (2015). Genetic studies of body mass index yield new insights for obesity biology. Nature 518(7538), 197.
Luque, F. J. (1984). Asymptotic convergence analysis of the proximal point algorithm. SIAM Journal on Control and Optimization 22(2), 277-293
Matsui, H. and S. Konishi (2011). Variable selection for functional regression models via the I1 regularization. Computational Statistics & Data Analysis 55(12), 3304-3310.
Mirshani, A. and M. Reimherr (2019). Adaptive function-on-scalar regression with a smoothing elastic net. arXiv preprint arXiv: 1905.09881
Murdoch, W. J., C. Singh, K. Kumbier, R. Abbasi-Asl, and B. Yu (2019). Definitions, methods, and applications in interpretable machine learning. Proceedings of the National Academy of Sciences 16(44), 22071-22080.
Parikh, N., S. Boyd, al. (2014). Proximal algorithms. Foundations and Trends in Optimiza- tion 1(3), 127-239.
Parodi, A., M. Reimher,, et al. (2018). Simultancous variable selection and smoothing for high- dimensional function-on-scalar regression. Electronic Journal of Statistics /2(2), 4602-4639.
Paul, I. M., S. Williams, S. Anzman-Frasca, S. Beiler, K. D. Makova, M. E. Marini, L. B. Hess, S. E. Rzucidlo, N. Verdiglione, J. A. Mindell, et al. (2014). The intervention nurses start infants growing healthy trajectories (insight) study. BMC pediatrics 14(1), 184.
Pedregosa, F., G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Pretten- hofer, R. Weiss, V. Dubourg, et al. (201 1). Scikit-learn: Machine learning in python. the Journal of machine Learning research 12, 2825-2823
Polyak, B. T. (1969). The conjugate gradient method in extremal problems. USSR Computational Mathematics and Mathematical Physics 9(4,, 94-
Pulit, S. L.. C. Stoneman, A. P. Morris, A. R. Wood, C. A. Glastonbury, J. Tyrrell, L. Yengo, T. Ferreira, F. Marouli, Ji, et al. (2019). Meta-analysis genome-wide association studies for body fat distribution in 694 649 individuals of european ancestry. Human molecular genetics 28(1), 166-174.
Rockafellar, R. T. (1976a). Augmented lagrangians and applications of the proximal point algorithm in convex programming. Mathematics of operations research I(2), 97-116.
12

---

Rockafellar, R. T. (1976b). Monotone operators and the proximal point algorithm. SIAM journal on control and optimization 14(5), 877-898.
Sørensen, H., Goldsmith, and L. M. Sangalli (2013). An introduction with medical applications to functional data analysis. Statistics in medicine 32(30), 1222-52400
Tibshirani, R. J., J. Taylor, et al. (2012). Degrees of freedom in lasso problems. The Annals of Statistics 40(2), 1198-1232.
Tomioka, and M. Sugiyama (2009). Dual-augmented lagrangian method for efficient sparse reconstruction. IEEE Signal Processing Letters 16(12), 1067-1070.
Tomioka, R., T. Suzuki, and M. Sugiyama (2011). Super-linear convergence of dual augmented lagrangian algorithm for sparsity regularized estimation. Journal of Machine Learning Re- search 12(5).
Ullah, S. and C. F. Finch (2013). Applications of functional data analysis: systematic review. BMC medical research methodology I3(1), 1-12.
Van Loan, C. F. and G. H. Golub (1983). Matrix computations. Johns Hopkins University Press.
Yu, B. and K. Kumbier (2020). Sciences 117(8). 3920-3929.
Veridical data science. Proceedings of the National Academy of
Yuan, M. and Y. Lin (2006). Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68 1), 49-67.
Zhang, Y. and Q. Yang (2018). An overview of multi-tast learning. National Science Review 5(1), 40-43.
Zhao, S., A. Shojaie, and D. Witten (2017). In defense of the indefensible: very naive approach to high-dimensional inference. arXiv preprint (00dvv:7755.05550
7hu, Y. (2017). An augmented admm algorithm with application to the generalized lasso problem. Journal of Computational and Graphical Statistics 26(1), 195-204.
Zou, and Hastie (2005). Regularization and variable selection via the elastic net. Journal of the royal statistical society: series B (statistical methodology) 67(2), 301-320.
NeurIPS paper 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 describe the limitations of your work? [Ye]] We briefly describe some
limitations our methodology at the end of Section 5
(c) Did you discuss any potential negative societal impacts of your work? [No] We do not the potential for adverse societal impacts
(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 of assumptions all theoretical results? [Yes]
(b) Did you include complete proofs of all theoretical results? [Yes] All proofs are included
in Supplemental Section A
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed reproduce the main ex- perimental results (either in the supplemental material or URL)? [Yes] Code to reproduce the simulation results included the Supplemental Material.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] See Section 4
13

---

(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [Yes] We rum subset of experiments multiple times computing mean and standard errors. These results are included in Supplemental Section
(d) Did you include the total amount of compute and the type of resources used (e.g.. type of GPUs, internal cluster, or cloud provider)? [Yes] See Supplemental Section
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] We cited the creators of data from INSIGHT in Section 4.2
(b) Did you mention the license of the assets? [N/A]
(c) Did you include any assets either in the supplemental material URL? [Yes] The code to run fgen is included in the Supplementary Material and will be made available GitHub
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [No] The INSIGHT study operated under an appropriate IRB. We did not reference it or provide details here. We do thank the individuals who provided access to the data in the Acknowledgments Section.
(e) Did you discuss whether the data you using/curating contains personally identifiable information or offensive content? [Yes] The INSIGHT data we used for our application is privacy potected. We discuss how parts of it are available upon request through dbGaP, and parts of are not available, in Section 4.2
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent participant compensation? [N/A]
14


Reviews:
[{"review_id": "outB8URmNYr", "reviewer": "Reviewer_AsoJ", "summary": "The paper presents an optimization algorithm for functional data.\nIt is based on another optimization algorithm for very high-dimensional data.\nThe objective function is least squares with the usual group-lasso penalization.\nThe novelty claimed in the paper is that the presented algorithm is for matrices not just vectors.\nThe paper then extends the algorithm to functional data. This is done by first reducing the functional data to vector data through Principal Components (PC).", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "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": "\nORIGINALITY\n\nThe objective function is least squares with the usual group-lasso penalization.\nThe novelty claimed in the paper is that the presented algorithm is for matrices not just vectors.\nHowever, the distance to minimize is the Frobenius distance:\nlet A = (A_1 | A_2 | ... | A_p) be a nxp matrix, where A_i is a column vector of dimension n.\nThe square of the Frobenius norm is sum_{ij} (a_{ij})^2 = sum_j sum_i (a_{ij})^2 = sum_j || A_j ||^2\nThe objective function is || Y - XB||^2 = sum_j || Y_j - XB_j ||^2\nHence, this is just the sum of square Euclidean norms of the columns. If we add to this the group-lasso\nterms we get a sum of separate terms to penalize, one for each column:\n\nf(B) = sum_j f_j(B_j); and min{ f(B) : B}  = sum_j min{ f_j(B_j) : B_j}.\n\nThere is no joint minimization. Therefore, as presented, the algorithm is not a new algorithm; and\nall what is known for the vector least-squares with group lasso applies directly to this problem.\n\nSIGNIFICANCE\n\nThe paper extends the algorithm to functional data. This is done by first reducing the functional data\nto vector data through Principal Components (PC). This is done is several works with functional data. One of\nthe main issues with PC is to decide what is the intrinsic dimension k of the data, that is, the number of vectors\nneeded to represent the functional data. This is a hard problem. There are several works on estimating the intrinsic\ndimension while doing PC for other tasks (e.g., clustering). In the present paper, k is estimated beforehand to be about\n5. This is done by looking at the mount of variance explained by the k vectors. Another issue is what is the \"right\"\namount of variation that will be used as threshold. This has not been seriously addressed in the paper.\nOnce k is fixed, the problem becomes the above group-lasso optimization.\n\nAs the title states it, this paper should have focused on this latter problem (functional data) in more\ndepth, addressing all issues arising from functional data. Instead, the paper emphasizes the matrix group-lasso\noptimization, and dedicates a handful of paragraphs to the functional data problem.\n\n\nQUALITY\n\nProposition 1: the authors say that\n(a) pi*(Z) is a differentiable function. This is not true. The obvious case is the case\nk=1, where the Euclidean norm becomes absolute value. But the function is not differentiable for any k.\nIt is not differentiable at Z_i = 0, but it is not differentiable even if Z_i=(0,0,..., z_j, 0, ...,0).\n(b) that they have extended pi* to the case of matrices; but as said above, with the squared of the Frobenius norm\nthe matrix is treated as a vector; hence, there is no generalization here.\n\nLine 183: it is stated that pi* is a differentiable function. As shown above, this is not true.\n\nlines 31, 32: it is stated \"Given the complex, noisy nature of much contemporary data\"\nThis statement should be supported with reasons of references.\n\nSection 3.1: the efficiency of the algorithm presented lies in the two minimization steps of equation (10),\nwhich are realized as one Newton step for each. These are rough approximations to the minimization steps.\nThey will certainly speed up the algorithm, but it will not guarantee convergence to the optimal solution.\n\nTable 1: here the computation times are shown. The residuals should also be shown. They are more important than the computation times.\nOne gain nothing having a fast algorithm that performs poorly.\n\nSection 4.2: the application to the INSIGHT study is pretty lean. It looks as just an example of application.\n\n\n\nCLARITY\n\nlines 33,34: it is stated \"...is paramount to capture significant signals and ensure the stability of outcomes\"\nPlease explain clearly what you want to say here and why.\n\nlines 38, 39: define l_1 and l_2 norms (notation).\n\nline 39: it is stated \"...and regularize the estimates\" What is the meaning of \"to regularize estimates\" ?\n\nline 40: define k. What groups are you taking about?\n\nEquation (1) none of the terms in this equation are defined. You should define most of the terms before.\nInstead , the definitions of the notation comes after a period '.'\nIt would be better to change the period in the equation with a ',' and start the following sentence  with\n\"where...\"\n\nYou also talk about groups here. What groups? What is a \"component\"?\nDefine lambda_1, lambda_2\n\nline 107: please provide references to Moreau's decomposition.\n\nEquation (D): please explain notation '|'. Do you mean 'subject to' ?\n\nSection 3.1: the functions L_sigma(V | Z, B) and L( Z | V, B) have not been defined. I think they are never defined in the paper.\n\nProposition 3: psi(V)= L_sigma(V | Z, B) is supposed to be a function of Z; Z is given in (b), but it would be better if\nZ enters explicitly in the definition of psi(V).\n\nTheorem 1: Please define objects clearly: \"generalized Hessian\", \"Clarke sub-differential\".\n\nLine 151: What is the meaning of the symbol \"belongs to\" in 'Q belongs-to partial_prox(T)'.\nIs 'partial_prox(T)' a subset? From equation (13) (iii) P_i is an expression and not a subset.\n\nEquation (16): What kind of residuals are these? Please explain more. What are they good at? Why should we look at them?\n\nlines 224 to 231: it is stated that  \"de-biased estimates\" (why are they de-biased? where is the bias?) are obtained by fitting standard least squares on the selected features, which, following the authors, is an \"effective approach when p >> n\".\nHowever, it is know that least-squares are not are not to trust when p >> n. The whole objective of doing lasso when p >> n is to reduce\nthe dimension of the coefficients to some q << n, so as to get reliable estimates through, say, least-squares. Please comment on this issue.\n\nFigure 2: the image in the middle panel has a strange line at the bottom left corner (I think it is for alpha= 0.8).\nPlease comment on that.\n\nlines 253, 254: Why is sigma estimated? I thought this was a fixed parameter of the kernel.\n\n\nSection 4: I wonder if PC is also applied in the competitor methods sklearn and glmnet.\nIf it is not, then the comparison in times is not that fair, since the competitors would be working with larger objects.\nAlso, glmnet works with a Gaussian model, not just least-squares, so it has to estimate a covariance matrix.\nPlease give more details on what models (and commands) of glmnet and Sklearn were used in the comparison.\n\nAlso even though glmnet is apparently mostly written in FORTRAN, it is still a package of R, and there must be\nR scripts behind. These probably slows the running time.\n\n\n\nOTHER\n\nline 99: change \"denetoded\" to \"denoted\"\n\nline 279: change 'skleanr' to 'sklearn'\n\n\n\nREBUTTAL\n\nI have read the authors response, and the other reviews. I thank the authors for their detailed answers to all my questions. They have helped me understand better what they are proposing, which, in effect is a Group Elastic net. As the authors correctly claim, the optimization problem cannot be separated. I also now understand pi*: the authors are right, this function is differentiable.\nI would like to note that there is no l1-norm used in Group elastic net; it is just the l2 norm; the other term is the l2-norm-squared. So it is a bit confusing when the authors talk about the l1-norm in their procedure. It is the sparse group lasso that uses the l1-norm and the l2-norm, and the elastic net that uses the l1-norm and the l2-norm-squared.\n\nThe Group elastic net optimization problem is not new. The authors contribution is an efficient implementation of it. Therefore, the main result is Table 1 with the timing of the methods compared. In any case, it was not clear that to me that the methods were producing the same results.\n\n\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "\nLimitations have been mentioned in the conclusions as well as some future extensions of the model.\n\nSocietal impact: The authors see no potential for adverse societal impact of their work.\nHowever, every algorithm that searches for relations between disease and genes may be have unintended consequences,\naffecting certain individuals (e.g., health insurance, discrimination, etc.)\n\n", "review_text": "", "extra_scores": {"ethical_concerns": "No ethical concerns have been raised.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "10"}, "extra_text": {}}, {"review_id": "h8K84Hzy1ng", "reviewer": "Reviewer_zxwr", "summary": "The authors propose a semi-smooth augmented Lagrangian method for solving the least-squares vector-valued regression with group elastic net regularization in the high-dimensional setting. The proposed algorithm is able to exploit the second-order sparsity as well as achieve a super-linear convergence rate.  The application of this algorithm to function-on-scalar regression is also presented.", "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": "The semi-smooth augmented Lagrangian (SsNAL) algorithm is an efficient method for solving large-scale high-dimensional least-squares regression including Lasso and elastic net problems, where the augmented Lagrangian algorithm can guarantee the super-linear convergence and the semi-smooth Newton method can exploit the second-order sparsity. An interesting question is to extend the algorithms to solve so-called function-on-scalar regression problems with a functional response and a set of scalar covariates. In this paper, the authors propose the use of SsNal for vector-valued regression by imposing the group elastic net regularization in the high-dimensional setting. By virtue of functional principal component analysis (FPCA),  the proposed algorithm is extended to the function-on-scalar features selection problems. Numerical experiments on both simulated data and real data sets are presented to verify the efficiency of the proposed algorithm in solving function-on-scalar regression problems. Overall, the paper is well-written and easy to read. The presented method for solving function-on-scalar regression is of interest, but I have several concerns as follows:\n\n1. The contribution of the paper may be not novel enough for NeurIPS as it just extends the results of scalar-valued regression in Boschi et al. (2020) to vector-valued regression, which is usually quite straightforward in linear regression.\n\n2. The paper adopts FPCA to reduce function-on-scalar regression into vector-valued regression and then applies the proposed algorithm. The benefit of employing functional data analysis is not considered and discussed. Moreover, numerical comparisons of accuracy with other existing methods for function-on-scalar regression are not provided. Comparing only CPU time may not be enough to illustrate the performance of the proposed algorithm.\n\n3. For the functional outcome which captures the evolution of weight/high ratios, how about directly using the five-dimension response vector instead of functional data? Any comparison or discussion about this?\n\n\nAfter Rebuttal\n\nThank the authors for their detailed feedback. I appreciate and agree with the clarification about the novelty of the paper.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "RNjMU6FJc0O", "reviewer": "Reviewer_kboM", "summary": "This paper proposes a new efficient algorithm to solve the group elastic net regularized multivariate regression problem in ultra-high-dimensional settings. The algorithm (SsNAL) combines the augmented Lagrangian method with semi-smooth Newton method to achieve super-linear convergence and reduce the computational cost. Experiments with simulated data and a real dataset show a substantial reduction in CPU time of SsNAL with respect to the best existing group elastic net solvers.", "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": "Comments:\n\n1. It is unclear to me what are the motivation and real advantages of using Algorithm 1 to solve problem (D), compared with using the alternating direction method of multipliers (ADMM). Algorithm 1 looks complicated since it has two loops (the inner loop \"Semi-smooth Newton method\" is nested within the outer loop of the augmented Lagrangian method).\n\nDifferent from the augmented Lagrangian method which jointly update $V$ and $Z$ as in Eq. (7), ADMM updates $V$ and $Z$ in an alternating or sequential fashion:\n$$ V^{k+1} = \\arg\\min_V \\mathcal{L}_\\sigma\\left(V, Z^k, B^k\\right)$$\n\n$$Z^{k+1} = \\arg\\min_Z \\mathcal{L}_\\sigma\\left(V^{k+1}, Z, B^k\\right)$$\n\n$$B^{k+1} = B^k - \\sigma \\left(X^T V^{k+1} + Z^{k+1}\\right)$$\n\nand moreover, it can be shown that both $V$-update and $Z$-update subproblems admit a closed-form solution:\n\n$$V^{k+1} = \\left(I_n + \\sigma X X^T\\right)^{-1} \\left(X B^k - Y - \\sigma X Z^k\\right)$$\n\n$$Z_i^{k+1} = {\\rm prox}_{\\frac{\\pi^*}{\\sigma}}(T_i) = \\left[1 - \\frac{1}{\\sigma\\lambda_2 + 1}\\left(1-\\frac{\\lambda_1}{\\|T_i\\|_2}\\right)_+\\right] T_i, \\; i=1,2,\\dots,p,$$\n\nwhere $I_n$ denotes the identity matrix of size $n$, and $T = B^k/\\sigma - X^T V^{k+1}$. As a result, the ADMM algorithm to solve problem (D) contains one loop only. Note that $n$ is relatively small in high-dimensional settings and thus $V^{k+1}$ can be efficiently computed using e.g., Gaussian elimination. Therefore, the ADMM algorithm is presumably computationally more efficient than Algorithm 1.\n\nIt would be interesting to empirically compare the running time of Algorithm 1 and ADMM for varying dimensions, to better understand the pros and cons of Algorithm 1 when used to solve problem (4).\n\n\n2. Mistakes, notation abuse, etc.:\n\nLine 41: For a matrix (e.g., XB-Y and V), it would be clearer to use the symbol \\\\|\\cdot\\\\|_F rather than \\\\|\\cdot\\\\|_2 to denote its Frobenius norm. The symbol \\\\|\\cdot\\\\|_2 when applied on a matrix usually means the spectral norm, which is a matrix norm induced by the L2-norm for vectors.\nEq. (D) in line 110: It should reverse the sign of the objective function, or equivalently, change the \\min operator to \\max.\nLine 113: \\sum_{i=1}^k (Y^T)_i V_i should be written as \\sum_{i=1}^n (Y^T)_i^T V_i or \\langle Y, V\\rangle.\nEq. (5): (X^T)_i V should be written as V^T (X^T)_i. Note that to be consistent throughout the paper (and also to follow the convention), a vector (e.g., the i-th row of a matrix) is assumed to be a column vector, not a row vector.\nThe \\nabla \\psi(V) and \\hat{\\partial}^2 \\psi(V) and their sizes are not consistently defined (as a vector/matrix/tensor):\n- in Eq. (12) \\nabla \\psi(V) is defined as a matrix of size n \\times k but in Eq. (14) it is used as a vector of dimension nk.\n- in Eq. (12) \\hat{\\partial}^2 \\psi(V) is defined as a matrix of size nk \\times nk but in line 154 it says \\partial^2 \\psi is \"a tensor of order 3\".\nEq. (16): The i-1 should be i=1.\nEq. (A11) in the supplement, right-hand side of the second equality: The first line is missing T_i and in the second line the 0 should be T_i/\\sigma.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"ethical_concerns": "None.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes an efficient algorithm to solve the group elastic net regularized multivariate regression problem. The proposed SsNAL algorithm combines the augmented Lagrangian method with semi-smooth Newton method to achieve super-linear convergence. Several reviewers had problems understanding part of the paper. The presentation of the paper should be improved in the final version of the paper.
Decision: Accept (Poster)