# A Fast and Scalable Pathwise-Solver for Group Lasso and Elastic Net
Penalized Regression via Block-Coordinate Descent

James Yang  Department of Statistics, Stanford University  Trevor Hastie
Department of Statistics, Stanford University

###### Abstract

We develop fast and scalable algorithms based on block-coordinate descent to
solve the group lasso and the group elastic net for generalized linear models
along a regularization path. Special attention is given when the loss is the
usual least squares loss (Gaussian loss). We show that each block-coordinate
update can be solved efficiently using Newton’s method and further improved
using an adaptive bisection method, solving these updates with a quadratic
convergence rate. Our benchmarks show that our package adelie performs 3 to 10
times faster than the next fastest package on a wide array of both simulated
and real datasets. Moreover, we demonstrate that our package is a competitive
lasso solver as well, matching the performance of the popular lasso package
glmnet.

##  1 Introduction

The group lasso [Yuan and Lin, 2006] is an extension of the celebrated lasso
[Tibshirani, 1996] that aims to additionally capture group structure within
the features to induce group-wise sparsity in the solutions. Specifically, the
group lasso tends to zero-out _groups of coefficients_ rather than individual
coefficients like the lasso. The group lasso is a bona fide extension in that
if all groups are of size one, then we recover the lasso. Despite the wide
popularity of the lasso for its interpretability of its sparse solutions,
there is a drawback when multiple columns of a feature matrix represent one
entity. For example, a categorical feature, or a factor, is often represented
as indicators for each category, yielding many columns. In such cases, the
group lasso will select all or none of these columns whereas the lasso may
select any subset of them, leaving the user in an ambiguous situation to
decide whether the factor is informative or not [Meier et al., 2008].
Moreover, the result of lasso changes based on the encoding scheme of factors,
which is an undesirable property [Yuan and Lin, 2006]. For these reasons, the
group lasso has gained popularity over the years.

While the group lasso seems very close in spirit to the lasso (as both models
induce sparsity in their solutions), their computational aspects are vastly
different. A popular method for solving the lasso is coordinate-descent
[Friedman et al., 2010]. The authors demonstrate that coordinate-descent can
be made extremely efficient with a closed-form solution to each coordinate
update. As we shall see in Section 3, the group lasso does not enjoy such
properties as soon as a group size is greater than one. Still, inspired by the
success of coordinate-descent with the lasso, many authors pursue (block)
coordinate-descent as the primary workhorse of their algorithm [Yuan and Lin,
2006; Meier et al., 2008; Tseng, 2001; Liang et al., 2022]. To solve each
(block) coordinate update, people have gravitated towards proximal gradient
methods such as ISTA, FISTA, and FISTA with adaptive-restarts [Liang et al.,
2022; Beck and Teboulle, 2009; Klosa et al., 2020; Wright et al., 2009; Loris,
2009; Hastie et al., 2016; O’Donoghue and Candés, 2015]. These methods are
attractive for problems like the group lasso because the block-coordinate
updates can be solved iteratively where each iteration has a simple closed-
form solution, now reminiscient of the lasso.

In this paper, we develop fast and scalable algorithms for solving the group
lasso under the usual least squares loss (Gaussian loss) as well as any twice-
continuously differentiable convex, proper, and closed losses. In particular,
we can solve generalized linear models such as logistic regression, poisson
regression, and multinomial regression with the group lasso penalty. Following
the work of Zou and Hastie [2005], which extends the lasso to the elastic net,
we work with an analogous elastic net penalty for the group lasso. Our work
largely follows the work of Friedman et al. [2010] and our main result lies in
the block-coordinate update. Specifically, we propose a novel algorithm to
solve each block-coordinate update via Newton’s method, which can be shown to
have a convergence guarantee and a quadratic convergence rate. We provide a
publicly available Python package adelie with brief installation directions in
Section 7.

In Section 3, we present our main result for solving the group elastic net
under the Gaussian loss using block-coordinate descent with a special emphasis
on the block-coordinate update. Section 4 extends this result to general
convex losses using proximal quasi-Newton method. Section 5 demonstrates that
we may then use our results to easily fit multi-response data as well with no
special changes to the solver. Finally, Section 6 shows our benchmark results
comparing our package to existing R packages for solving the group lasso as
well as the lasso.

##  2 Preliminaries and Notations

Let $\mathbb{R}$, $\mathbb{R}^{n}$, and $\mathbb{R}^{m\times n}$ denote the
space of real numbers, real vectors of size $n$, and real matrices of shape
$m\times n$, respectively. Similarly, we define the analogous spaces for the
non-negative real numbers $\mathbb{R}_{+}$ and the positive real numbers
$\mathbb{R}_{++}$. We denote $\ell_{p}$-norms on $\mathbb{R}^{n}$ as
$\left\|x\right\|_{p}:=(\sum_{i=1}^{n}\left|x_{i}\right|^{p})^{1/p}$. We write
$\left\|x\right\|_{W}:=\sqrt{x^{\top}Wx}$ to be the $\ell_{2}$-norm scaled by
weights $W\in\mathbb{R}_{+}^{n\times n}$. We denote $\mathbf{0}$ and
$\mathbf{1}$ to be the vector of zeros and ones, respectively. The identity
matrix is denoted as $I_{n}\in\mathbb{R}^{n\times n}$ and $I$ if the
dimensions are implicitly defined based on context. We define
$[n]:=\left\\{1,\ldots,n\right\\}$ to be the set of integers ranging from $1$
to $n$. For any vector $v\in\mathbb{R}^{n}$ and a set of indices
$S\subseteq[n]$, let $v_{S}$ be the subset of $v$ along indices in $S$. If $S$
is a subset of $G$ groups, that is, $S\subseteq[G]$ and each group $g\in[G]$
has size $p_{g}$ such that $\sum_{g=1}^{G}p_{g}=p$, then $v_{S}$ is the subset
of $v$ along each group so that $v_{S}=\left\\{v_{g}\in\mathbb{R}^{p_{g}}:g\in
S\right\\}$. In either context, we denote $v_{-S}\equiv v_{S^{c}}$ as the
subset of $v$ along indices (groups) not in $S$. For a real-valued function
$f:\mathbb{R}^{n}\to\mathbb{R}$, we denote its gradient as $\nabla f$ and the
hessian as $\nabla^{2}f$. For two square matrices $A,B\in\mathbb{R}^{n\times
n}$, we denote $A\preceq B$ if $B-A$ is positive semi-definite and $A\prec B$
if $B-A$ is positive definite. We denote $A\otimes B$ to be the Kronecker
product of any two arbitrary matrices $A$ and $B$.

##  3 Algorithms for Group Lasso/Elastic Net under Gaussian Loss

Consider the usual regression setup where $X\in\mathbb{R}^{n\times p}$ and
$y\in\mathbb{R}^{n}$ are the feature matrix and response vector, respectively.
Let $G$ be the number of groups and $p_{g}$ be the size of group $g\in[G]$ so
that $X$ may be decomposed as $\left[X_{1}\,X_{2}\,\ldots\,X_{G}\right]$ where
$X_{g}\in\mathbb{R}^{n\times p_{g}}$ is the feature matrix corresponding to
group $g$ features only. In this section, we discuss our main algorithm to
solve the group elastic net problem under the Gaussian loss given by

$\displaystyle\operatorname*{minimize\;}_{(\beta_{0},\beta)\in\mathbb{R}^{1+p}}\quad$
$\displaystyle\frac{1}{2}\left\|y-X\beta-\beta_{0}\mathbf{1}\right\|_{W}^{2}+\lambda
P_{\alpha,\omega}(\beta)$ (1)

and

$\displaystyle
P_{\alpha,\omega}(\beta):=\sum\limits_{g=1}^{G}\omega_{g}\left(\alpha\left\|\beta_{g}\right\|_{2}+\frac{1-\alpha}{2}\left\|\beta_{g}\right\|_{2}^{2}\right)$
(2)

where $\beta_{g}\in\mathbb{R}^{p_{g}}$ are the coefficients corresponding to
group $g$. Here, $W\in\mathbb{R}_{+}^{n\times n}$ is assumed to be a diagonal
matrix where the diagonal specifies the observation weights. Moreover, we
assume that the sum of the diagonal $\mathbf{1}^{\top}W\mathbf{1}$ is exactly
$1$ so that the weights are normalized. Typically, we use the uniform weights,
which results in $W\equiv n^{-1}I$. The penalty factor
$\omega\in\mathbb{R}_{+}^{G}$ allows us to control the _relative_ penalty
across the groups and $\lambda\geq 0$ to control the _overall_ penalty. A
common strategy is to set $\omega_{g}\equiv\sqrt{p_{g}}$ so that the groups
are penalized relative to their size [Yuan and Lin, 2006]. The elastic net
parameter $\alpha\in[0,1]$ dictates the proportion of the group lasso penalty
($\left\|\cdot\right\|_{2}$) and the ridge penalty
($\left\|\cdot\right\|_{2}^{2}$) so that $\alpha$ closer to $1$ induces more
sparsity in the solutions. Note that the intercept $\beta_{0}$ is unpenalized.
Hence, without loss of generality, we assume for simplicity that $X$ is
column-centered and $y$ is centered (where the center is a weighted average
using the diagonal of $W$ as weights) so that we may replace $\beta_{0}\equiv
0$. This simplification is specific to the Gaussian loss and we will soon
remove this assumption when we move to general convex losses in Section 4.

Although the groups typically have less correlation _between_ groups than
_within_ groups, as most applications tend to group correlated features
together, it is nonetheless still possible to have strong correlation between
groups. In this case, the group lasso can still suffer from the same erratic
behavior observed with the lasso where the group lasso will tend to pick one
of the correlated groups and ignore the rest. For this reason, the group
elastic net penalty can help ameliorate this extreme behavior by applying some
ridge penalty to shrink every group towards each other.

Following [Friedman et al., 2010], we show in Figure 1 a comparison of the
solution paths for the group lasso, group elastic net, and ridge on the
Leukemia dataset from Golub et al.
[1999]111https://www.openintro.org/data/index.php?data=golub. This dataset
contains 72 samples and 7129 gene expression probes that are measured with DNA
microarrays. To introduce group structure, we take each column of gene probes
$X_{\cdot j}$ and create two more features $X_{\cdot,j}^{2}$,
$X_{\cdot,j}^{3}$, that is, the square and cube of the original column
(element-wise). We then group every three such features together. The response
$y$ is a binary variable indicating cancer or not for each observation.
However, for the sake of illustration, we treat it as a continuous variable.
We run the group lasso ($\alpha=1$), group elastic net ($\alpha=0.2$), and
ridge ($\alpha=0$) for $100$ values of $\lambda$, and plot the coefficient
profiles for the first $10$ values. The total time to run all three models for
all $100$ values of $\lambda$ was $0.5$ seconds on a standard M1 Macbook Pro.
It is clear that group lasso tends to pick sparse models with few large
coefficients whereas the ridge includes all the features but shrinks them
towards zero. The group elastic net is somewhere in between where we observe
more non-zero coefficients than the group lasso but fewer than the ridge.

Figure 1: A comparison of the coefficient profile for the group lasso, group
elastic net, and the ridge on the Leukemia dataset by [Golub et al., 1999].
From left to right, we gradually see more features entering the model with an
overall shrinkage on the coefficients towards zero until the ridge includes
every feature.

Our method to solve 1 is based on the block-coordinate descent (BCD) algorithm
[Tseng, 2001]. Since the objective is composed of a smooth, convex, closed,
and proper component and a separable non-differentiable component, BCD
provably attains the minimizer. The BCD algorithm cycles through each group
$g=1,\ldots,G$ and solves 1 only for $\beta_{g}$ while fixing all other
coefficients. In other words, we solve the following optimization problem:

$\displaystyle\operatorname*{minimize\;}_{\beta_{g}\in\mathbb{R}^{p_{g}}}\quad$
$\displaystyle\frac{1}{2}\beta_{g}^{\top}X_{g}^{\top}WX_{g}\beta_{g}-(y-\sum_{g^{\prime}\neq
g}X_{g^{\prime}}\beta_{g^{\prime}})^{\top}WX_{g}\beta_{g}+\lambda\omega_{g}\left(\alpha\left\|\beta_{g}\right\|_{2}+\frac{1-\alpha}{2}\left\|\beta_{g}\right\|_{2}^{2}\right)$
(3)

As we show in Section 3.1, it is advantageous to first perform an eigen-
decomposition of $X_{g}^{\top}WX_{g}=Q_{g}\Lambda_{g}Q_{g}^{\top}$ and
transform $\beta_{g}\mapsto Q_{g}^{\top}\beta_{g}$. Here,
$Q_{g}\in\mathbb{R}^{p_{g}\times p_{g}}$ is a square orthogonal matrix and
$\Lambda_{g}\in\mathbb{R}_{+}^{p_{g}\times p_{g}}$ is a non-negative diagonal
matrix, which may contain some zeros. Using the fact that the penalty is
invariant under rotation and letting
$\tilde{\beta}_{g}:=Q_{g}^{\top}\beta_{g}$, 3 changes to

$\displaystyle\operatorname*{minimize\;}_{\tilde{\beta}_{g}\in\mathbb{R}^{p_{g}}}\quad$
$\displaystyle\frac{1}{2}\tilde{\beta}_{g}^{\top}\Lambda_{g}\tilde{\beta}_{g}-(y-\sum_{g^{\prime}\neq
g}X_{g^{\prime}}\beta_{g^{\prime}})^{\top}WX_{g}Q_{g}\tilde{\beta}_{g}+\lambda\omega_{g}\left(\alpha\left\|\tilde{\beta}_{g}\right\|_{2}+\frac{1-\alpha}{2}\left\|\tilde{\beta}_{g}\right\|_{2}^{2}\right)$
(4)

We abstract this problem into the following general problem

$\displaystyle\operatorname*{minimize\;}_{x\in\mathbb{R}^{d}}\quad$
$\displaystyle\frac{1}{2}x^{\top}\Sigma
x-v^{\top}x+\lambda\left\|x\right\|_{2}$ (5)

where $\Sigma\in\mathbb{R}_{+}^{d\times d}$ is a non-negative diagonal matrix
and $\lambda\geq 0$. Relating back to 4, we have
$\Lambda_{g}+\lambda\omega_{g}(1-\alpha)I\mapsto\Sigma$,
$Q_{g}^{\top}X_{g}^{\top}W(y-\sum_{g^{\prime}\neq
g}X_{g^{\prime}}\beta_{g^{\prime}})\mapsto v$, and
$\lambda\omega_{g}\alpha\mapsto\lambda$ (with a slight abuse of notation). The
rest of the section focuses on solving 5. Since $\lambda>0$ is the only
interesting case, we assume without loss of generality that $\lambda>0$.

We make a few remarks on 5. Suppose for the moment that the observation
weights $W=n^{-1}I$ are equal weights to be consistent with the other
references. In the literature, we have found that most, if not all, methods
only consider isotropic $\Sigma$, that is, $\Sigma=tI$ for some scalar $t>0$.
Yuan and Lin [2006] directly assume that $X_{g}$ in 3 has orthonormal columns
so that $\Sigma\equiv I$. In the general case, they at least assume $X_{g}$ to
be full-rank so that it admits a QR decomposition with an invertible $R$.
Then, they solve 3 with $X_{g}$ replaced with the $Q$ matrix and naively
rescale the result by $R^{-1}$. This strategy has drawbacks because it assumes
full-rankness and, more fundamentally, the resulting coefficient is _not_ the
solution to the original problem 3. This was also observed and noted in Simon
and Tibshirani [2012] and the authors proposed to fix the issue by amending
the group lasso penalty from $\left\|\beta_{g}\right\|_{2}$ to
$\left\|X_{g}\beta_{g}\right\|_{2}$. While this allows the derivation to work
out nicely, we otherwise view this change merely as a mathematical convenience
and not a satisfactory solution to the original problem 3. Other authors such
as Liang et al. [2022]; Beck and Teboulle [2009]; Meier et al. [2008] majorize
$X_{g}^{\top}X_{g}$ using an isotropic matrix. They are able to show
convergence and a convergence rate of $O(1/k^{2})$ where $k$ is the number of
iterations. However, we noticed in practice that their methods converge very
slowly either when $X_{g}^{\top}X_{g}$ is near singular or the dimension of
the block increases. As the block-coordinate update is the most critical
component of the BCD algorithm, we desire a more reliable and a faster
algorithm. Our method directly solves 5 with no workarounds or approximations
and works in full generality for any non-negative diagonal $\Sigma$, or
equivalently, any _arbitary_ $X_{g}$. Moreover, we attain a quadratic
convergence rate and a stable algorithm in practice.

It was previously noted that we must perform an eigen-decomposition for
$X_{g}^{\top}WX_{g}$. This decomposition only needs to be done _at most once
per group_. Moreover, we will see later in Section 3.6 that the decomposition
only needs to be done on a small subset of the groups.

We conclude this section with the statement of Theorem 3.1 (see Appendix A for
the proof). We state this theorem for two reasons. First, the problem in 5 is
indeed more general than 4. If $\Sigma$ is not positive definite, then for a
specific choice of $v$, 5 may not have a solution due to an unbounded
objective. Since we will now restrict our attention to 5, we provide a
sufficient condition for the existence of a solution for completion. It is
easy to check via the SVD of $X_{g}$ that the conditions of Theorem 3.1 hold
for our original problem 3. Secondly, condition 6 will become crucial in
Section 3.2 where we develop our adaptive bisection method.

######  Theorem 3.1 (Sufficient Condition for the Existence of the Minimizer
of 5).

Consider the optimization problem of 5 with $\lambda>0$. If the following
condition holds

$\displaystyle v_{S}=\mathbf{0},\quad
S:=\left\\{i\in[d]:\Sigma_{ii}=0\right\\}$ (6)

then a finite minimizer exists. Moreover, the minimizer is given by

$\displaystyle x^{\star}$
$\displaystyle=\begin{cases}\left(\Sigma+\frac{\lambda}{\left\|x^{\star}\right\|_{2}}I\right)^{-1}v,&\left\|v\right\|_{2}>\lambda\\\
0,&\left\|v\right\|_{2}\leq\lambda\end{cases}$ (7)

###  3.1 Vanilla Newton’s Method-Based Algorithm

We now discuss how to solve 8 via Newton’s method. Assume 6 holds so that a
solution exists by Theorem 3.1. Without loss of generality, we only consider
the first case of 7 when $\left\|v\right\|_{2}>\lambda$. Without knowing
$\left\|x^{\star}\right\|_{2}$, there is no closed-form solution for
$x^{\star}$. In an endeavor to find an expression for
$\left\|x^{\star}\right\|_{2}$, we take the $\ell_{2}$-norm on both sides to
get that $\left\|x^{\star}\right\|_{2}$ must satisfy

$\displaystyle\sum\limits_{i=1}^{d}\frac{v_{i}^{2}}{(\Sigma_{ii}\left\|x^{\star}\right\|_{2}+\lambda)^{2}}=1$
(8)

Unfortunately, there is no closed-form solution for
$\left\|x^{\star}\right\|_{2}$ from 8 either. We see this as the main
difference from performing coordinate-descent for lasso. Namely, the
coordinate update in lasso has a simple closed-form solution (soft-
thresholding function) [Friedman et al., 2010]. However, as soon as a group
has size greater than 1, there is no closed-form solution to its block-
coordinate update in the group lasso. In spite of this difficulty, we shall
see shortly that a simple Newton’s method can numerically solve 8 efficiently
with guaranteed convergence. In turn, once we find
$\left\|x^{\star}\right\|_{2}$, we may substitute it into 7 to get the full
solution $x^{\star}$ in $O(d)$ time. Note that we can achieve this complexity
because $\Sigma$ is diagonal. Indeed, the subsequent arguments will also
heavily rely on this fact to achieve $O(dk)$ time complexity, where $k$ is the
number of Newton iterations, which would otherwise be $\Omega(d^{2}k)$ if
$\Sigma$ were not diagonal. This computation reduction was the primary reason
for transforming our problem from 3 to 4. We also note in passing that Hastie
et al. [2016, Ex. 4.7] also discovered the relationship between 7 and 8 from
which they proposed to solve for $\left\|x^{\star}\right\|_{2}$ via the
golden-section search.

Now, define $\varphi:[0,\infty)\to\mathbb{R}$ by

$\displaystyle\varphi(h)$
$\displaystyle:=\sum\limits_{i=1}^{d}\frac{v_{i}^{2}}{(\Sigma_{ii}h+\lambda)^{2}}-1$
(9)

so that 8 is now a one-dimensional root-finding problem for $\varphi$. In
Appendix C, we show that $\varphi$ is strictly decreasing, strictly convex,
and $\varphi(0)>0>\lim_{h\to\infty}\varphi(h)$, so that a root (uniquely)
exists. Our strategy is then to apply Newton’s method to search for the root
of $\varphi$. Newton’s method starts with an initial point $h^{(0)}$ and
iteratively updates

$\displaystyle
h^{(k+1)}=h^{(k)}-\varphi^{\prime}(h^{(k)})^{-1}\varphi(h^{(k)})$ (10)

Algorithm 1 summarizes Newton’s method. By Theorem 3.2, Newton’s method has a
guaranteed convergence to the root of $\varphi$ for any initial point
$h^{(0)}\geq 0$ such that $\varphi(h^{(0)})\geq 0$ (e.g. $h^{(0)}\equiv 0$).
For a proof of Theorem 3.2, see Appendix B. Moreover, since $\varphi$ is
strongly convex with Lipschitz second derivative, we also have a quadratic
convergence rate [Boyd and Vandenberghe, 2004].

######  Theorem 3.2 (Convergence Result of Newton’s Method).

Suppose $f:\mathbb{R}\mapsto\mathbb{R}$ is convex, continuously differentiable
with a non-vanishing derivative, and decreasing. Suppose further that there
exists a root $r\in\mathbb{R}$ such that $f(r)=0$. Let $x^{(0)}$ be any
initial point such that $f(x^{(0)})\geq 0$. Then, the Newton iterates
$\left\\{x^{(k)}\right\\}_{k=0}^{\infty}$ defined by

$\displaystyle x^{(k+1)}=x^{(k)}-f^{\prime}(x^{(k)})^{-1}f(x^{(k)})$ (11)

is an increasing sequence that converges to a root of $f$.

Data: $\Sigma$, $v$, $\lambda$, $\varepsilon$, $h^{(0)}$

1 $h\leftarrow h^{(0)}$;

2 while  _$\left|\varphi(h)\right| >\varepsilon$_ do

3 $h\leftarrow h-\frac{\varphi(h)}{\varphi^{\prime}(h)}$;

4

5return $h$;

Algorithm 1 Newton’s Method

It is clear from the discussion so far why our method, in principle, is faster
and more stable than the aforementioned proximal gradient methods such as
ISTA, FISTA, and FISTA with adaptive restarts. First, none of these methods
are able to give any convergence rate comparable to our quadratic convergence
rate largely because they are first-order methods. Secondly, while these
proximal gradient methods search the entire $\mathbb{R}^{d}$ space, we have
reduced our problem to a _one-dimensional_ root finding problem. Lastly, we
have established a well-behaving sequence of iterates (e.g. increasing and
bounded).

Figure 2 shows a comparison of the proximal gradient methods and our Newton
method, clearly demonstrating the superiority of our method. We implemented
all these methods as optimized as possible in C++. For FISTA and FISTA with
adaptive restarts, we follow the details laid out by the authors of Beck and
Teboulle [2009]; O’Donoghue and Candés [2015]. We generate the diagonal of
$\Sigma\in\mathbb{R}_{++}^{d\times d}$ from the uniform distribution on the
interval $(0,1)$, the linear term $v=\Sigma^{1/2}v^{\prime}$ where
$v^{\prime}\sim\mathcal{N}(0,I_{d})$, and $\lambda=0.1$. We measure accuracy
as

$\displaystyle\max\limits_{i=1,\ldots,d}\left|\hat{x}_{i}-\tilde{x}_{i}\right|$
(12)

where $\hat{x}$ is the solution returned by a solver and $\tilde{x}$ is
defined by 7 where the $\left\|x\right\|_{2}$ is replaced with
$\left\|\hat{x}\right\|_{2}$. Note that the proximal methods fail to converge
within 10000 iterations as $d$ increases, and even when they converge,
Newton’s method performs 10 to 1000 times faster. Furthermore, the accuracy of
the Newton’s method remains steady regardless of the dimension $d$, while the
proximal methods have trouble searching the space efficiently, suffering from
the lack of guidance from the curvature information.

Figure 2: Time and accuracy comparison between the proximal gradient methods
and our Newton’s method. Overall, the proximal gradient methods suffer from
large computational overhead due to large number of iterations, and fails to
converge as the dimension increases. Contrastingly, Newton’s method is quite
stable, properly converges, and runs 10 to 1000 times faster.

Despite the success of our Newton’s method approach, there is still room for
improvement. It is widely known that descent methods (including Newton’s
method) is highly sensitive to the initialization. A poor initialization may
lead to significantly slower convergence in practice while an initialization
close to the root will be rewarded with incredible speed. It then remains to
show how to properly choose the initial starting point $h^{(0)}$. First, it is
always advantageous to select $h^{(0)}$ such that $\varphi(h^{(0)})\geq 0$ to
preserve the increasing property of the iterates and the convergence
guarantees. We have chosen $h^{(0)}\equiv 0$ thus far, but it may not be the
best initial value as we may observe relatively larger number of iterations,
especially if $\Sigma$ is singular and $\lambda$ is small. The left plot in
Figure 3 shows an example of such a case. We see that $\varphi$ decays rapidly
near $h\approx 0$ and has an extremely flat tail. The iterates slowly move
from the origin until $h\approx 100$ and the speed picks up from there. One
possible remedy is to simply start Newton’s method from a large value of $h$
rather than $0$. However, as Figure 3 suggests, because the function usually
has a very flat tail, Newton’s method will likely shoot the iterate to the
left with a large magnitude. As a result, we may end up with an iterate that
is negative. Upon projecting back to $[0,\infty)$, we return to the situation
where we start Newton’s method at the origin. For these reasons, we develop in
Section 3.2 a special modification of our Newton-based approach, which
describes our most performant and robust method that reduces iterations in
nearly all configurations of the inputs.

Figure 3: Plot of the iterates for the vanilla Newton’s method (left) and
Newton-ABS (right). In both plots, we display the iterates by their value of
$\varphi$ and describe their path by the color gradient scheme. The iterates
are guaranteed to move from left to right. The solution is marked with a star.
For visualization purposes, the axes were scaled using the symmetrical log
transformation, which deceptively makes $\varphi$ look non-convex. It is clear
that Newton’s method struggles where there is a sharp decay near the origin as
the Newton iterates slowly exit the kink. Newton-ABS gets around this problem
by finding a good initial point sufficiently away from the origin using our
adaptive bisection strategy.

###  3.2 Newton’s Method with Adaptive Bisection Starts Algorithm

In this section, we improve our vanilla Newton’s method-based algorithm
described in Section 3.1. The key idea is to first apply a bisection method to
search for a good initial point. In particular, based on Figure 3, it is
important to avoid the region of fast decay when $h$ is close to $0$, as this
is the region with many Newton iterations. Once we find a point sufficiently
close to the root, we apply Newton’s method as a refinement procedure,
recovering all the convergence guarantees. Although one could use any
bisection method in theory, we propose an _adaptive bisection method_ that has
been most effective in our experience. We note in passing that we have
experimented with other bisection methods such as the naive one that bisects
evenly and other sophisticated methods such as Dekker’s method and Brent’s
method [Brent, 2013; Dekker, 1969]. Brent’s method, in particular, was the
only competitive method to our approach, however the runtime was quite
unpredictable as it was too dependent on the input values.

We first discuss how to find lower and upper bounds $h_{\star}$,
$h^{\star}\in[0,\infty)$, respectively, such that the root lies in
$[h_{\star},h^{\star}]$. To motivate the derivation, note that the root-
finding problem for $\varphi$ is equivalent to finding the largest value $h>0$
such that $\varphi(h)\geq 0$, or

$\displaystyle\sup\left\\{h>0:\varphi(h)\geq 0\right\\}$ (13)

The idea is to find a different function $\varphi_{\star}(\cdot)$ such that

$\displaystyle h_{\star}:=\sup\left\\{h>0:\varphi_{\star}(h)\geq
0\right\\}\implies\varphi(h_{\star})\geq 0$

and $h_{\star}$ is easily computable. Similarly, we find a function
$\varphi^{\star}(\cdot)$ such that

$\displaystyle h^{\star}:=\inf\left\\{h>0:\varphi^{\star}(h)\leq
0\right\\}\implies\varphi(h^{\star})\leq 0$

This shows that the root must lie in $[h_{\star},h^{\star}]$. In Appendix D,
we find suitable functions $\varphi_{\star}$ and $\varphi^{\star}$, and derive
the following expressions for $h_{\star}$ and $h^{\star}$:

$\displaystyle h_{\star}$
$\displaystyle=\left(\frac{-\lambda\mathbf{1}^{\top}\Sigma\mathbf{1}+\sqrt{(\mathbf{1}^{\top}\Sigma\mathbf{1})^{2}-(\mathbf{1}^{\top}\Sigma^{2}\mathbf{1})(\lambda^{2}d-\left\|v\right\|_{1}^{2})}}{\mathbf{1}^{\top}\Sigma\mathbf{1}}\right)_{+}$
(14) $\displaystyle h^{\star}$
$\displaystyle=\sqrt{\sum\limits_{i:\Sigma_{ii}>0}v_{i}^{2}\Sigma_{ii}^{-2}}$
(15)

We now describe our proposed adaptive bisection method. Ideally, we would like
to know if the root is closer to $h_{\star}$ or $h^{\star}$. If we believe
that the root is much closer to $h_{\star}$, we do not have to bisect
$[h_{\star},h^{\star}]$ at the mid-point, but perhaps at a point closer to
$h_{\star}$. Likewise, if the root were much closer to $h^{\star}$, we would
like to bisect closer to $h^{\star}$. Since we do not know the root a priori,
it is not possible to determine which side it is closer to. However, we can
make an informed guess based on a _prior_ of how likely the root is closer to
$h_{\star}$, for example. With this in mind, we construct our prior based on
the derivation of $h^{\star}$ in Appendix D. For completeness, we provide the
relevant part of the proof:

$\displaystyle\sum\limits_{i=1}^{d}\frac{v_{i}^{2}}{(\Sigma_{ii}h+\lambda)^{2}}$
$\displaystyle=\sum\limits_{i:\Sigma_{ii}>0}\frac{v_{i}^{2}}{(\Sigma_{ii}h+\lambda)^{2}}\leq
h^{-2}\sum\limits_{i:\Sigma_{ii}>0}\frac{v_{i}^{2}}{\Sigma_{ii}^{2}}$ (16)

It is easy to see that our $h^{\star}$ makes the right-side of 16 equal to 1.
We note that the inequality above becomes tighter as $h^{\star}$ approaches
the root. This motivates us to consider the worst approximation error (rate)

$\displaystyle
w:=\max\limits_{i:\Sigma_{ii}>0}\frac{\lambda}{\Sigma_{ii}h^{\star}+\lambda}=\frac{\lambda}{\Sigma_{\star}h^{\star}+\lambda}\in(0,1),\quad\Sigma_{\star}:=\min\limits_{i:\Sigma_{ii}>0}\Sigma_{ii}$

If $w\approx 0$, the inequality in 16 is tight, which suggests that the root
is close to $h^{\star}$. Hence, $w$ represents the prior that the root is
close to $h_{\star}$. We bisect at the new point
$h:=wh_{\star}+(1-w)h^{\star}$. If $\varphi(h)>\epsilon$ (sufficiently
positive), we use $h^{(0)}:=h$ as the initial point for our Newton’s method.
Otherwise if $\varphi(h)<-\epsilon$ (sufficiently negative), we set
$h^{\star}:=h$ and repeat the bisection argument. Finally, if
$\left|\varphi(h)\right|\leq\epsilon$ (sufficiently close to 0), we declare
$h$ to be the root. Algorithms 2 and 3 summarize this procedure, which we call
_Newton’s Method with Adaptive Bisection Starts_ (Newton-ABS). We emphasize
that we still have a convergence guarantee since we eventually apply Newton’s
method.

Data: $\Sigma$, $v$, $\lambda$, $\varepsilon$

1 Compute $h_{\star}$ using 14 and $h^{\star}$ using 15;

2 $\Sigma_{\star}\leftarrow\min\limits_{i:\Sigma_{ii}>0}\Sigma_{ii}$;

3 $h\leftarrow h^{\star}$;

4 while  _$\varphi(h) <-\varepsilon$_ do

5 $h^{\star}\leftarrow h$;

6 $w\leftarrow\frac{\lambda}{\Sigma_{\star}h^{\star}+\lambda}$;

7 $h\leftarrow wh_{\star}+(1-w)h^{\star}$;

8

return $h$;

Algorithm 2 Adaptive Bisection

Data: $\Sigma$, $v$, $\lambda$, $\varepsilon$

1 $h\leftarrow$ result of Algorithm 2;

2 if  _$\left|\varphi(h)\right| >\varepsilon$_ then

3 $h\leftarrow$ result of Algorithm 1 starting at $h^{(0)}=h$;

4

5$x^{\star}\leftarrow(\Sigma+\lambda h^{-1}I)^{-1}v$;

6 return $x^{\star}$;

Algorithm 3 Newton’s Method with Adaptive Bisection Starts (Newton-ABS)

Algorithm 3 can be further optimized. If $h^{\star}-h_{\star}$ is below some
threshold (e.g. $0.1$), then we may skip bisection entirely and start Newton’s
method at $h^{(0)}=h_{\star}$. This is because the adaptive bisection may move
too slowly if the range is too small. One may also enforce enough movement
towards $h_{\star}$ by taking the max of $w$ with a minimal probability (e.g.
$0.05$). This will ensure that at least some proportion of $h_{\star}$ is
taken if the prior suggests too strongly that the root is close to
$h^{\star}$. The idea is that it is better to overshoot towards $h_{\star}$
such that $\varphi$ is non-negative so that we can apply Newton’s method
sooner.

Figure 3 shows a plot of the Newton iterations for the vanilla Newton’s method
from Section 3.1 and our proposed Newton-ABS method. It is clear that Newton’s
method struggles where there is a sharp decay near the origin as the Newton
iterates slowly exit the kink. However, Newton-ABS gets around this problem
due to the fact that the adaptive bisection first moves in the opposite
direction, starting at the upper bound $h^{\star}$ and moving towards
$h_{\star}$, until $\varphi$ is non-negative. In this example, only 2 adaptive
bisections were required and 5 subsequent Newton iterations, which is a
noticeable improvement from the 34 vanilla Newton iterations. Empirically, we
have never observed a case where the initial value remains at the kink after
performing the adaptive bisection.

###  3.3 Naive and Covariance Updates

As in the coordinate descent algorithm for lasso, we have a “naive” and
“covariance” method to update the residual or gradient vector, respectively
[Friedman et al., 2010]. Denote $\beta$ as the current coefficient vector and
suppose in the BCD algorithm, we are at group $g$.

The “naive” method keeps track of the residual $r:=y-X\beta$. To solve the BCD
update 3, we must prepare the inputs to 4. Recall that we can cache the eigen-
decomposition of $X_{g}^{\top}WX_{g}$ so that we have access to its
eigenvalues $\Lambda_{g}$ and eigenvectors $Q_{g}$. Then, we may compute the
linear term in the following way

$\displaystyle Q_{g}^{\top}X_{g}^{\top}W(y-\sum_{g^{\prime}\neq
g}X_{g^{\prime}}\beta_{g^{\prime}})=Q_{g}^{\top}X_{g}^{\top}Wr+\Lambda_{g}Q_{g}^{\top}\beta_{g}$

In other words, we simply need to compute $\gamma_{g}:=X_{g}^{\top}Wr$,
express $\gamma_{g}$ and $\beta_{g}$ in the eigen-basis $Q_{g}$, and perform
simple elementwise operations. Altogether, the cost is $O((n+p_{g})p_{g})$.
Once the BCD update, $\tilde{\beta}_{g}$, is found by solving 3, we must
update the residual to reflect the change in $\beta_{g}$. Given the current
residual $r$, the change in residual, $\Delta r$, can then be expressed as

$\displaystyle\Delta r=-X_{g}\left(\tilde{\beta}_{g}-\beta_{g}\right)$

So, the residual update requires another $O(np_{g})$ operations. Sometimes, we
can avoid performing the residual update entirely if we observe that
$\tilde{\beta}_{g}\equiv\beta_{g}$, i.e. the BCD update did nothing.
Nevertheless, the “naive” method has a worst-case cost of
$O((n+p_{g}+k)p_{g})$ per BCD update where $k$ is the number of iterations for
the Newton-ABS method.

The “covariance” method keeps track of the gradient
$\gamma:=X^{\top}W(y-X\beta)$. This simplifies the input preparation for the
BCD update 3, since the linear term does not require the $O(np_{g})$
operations to compute $X_{g}^{\top}Wr$ as in the “naive” method so that the
complexity is $O(p_{g}^{2})$. However, it could potentially make the gradient
update more costly. The change in the gradient $\Delta\gamma$ can be expressed
as

$\displaystyle\Delta\gamma=-X^{\top}WX_{g}(\tilde{\beta}_{g}-\beta_{g})$

where we only have access to $X^{\top}WX_{g}$ as a whole. Hence, the gradient
update requires $O(p\cdot p_{g})$ operations. Therefore, the “covariance”
method has a worst-case cost of $O((p+p_{g}+k)p_{g})$, where $k$ is the number
of iterations for the Newton-ABS method.

Depending on the use-case, the user may find that one method is faster than
the other. Specifically, when $n\ll p$, the “naive” method is more efficient,
and when $p\ll n$, the “covariance” method is more efficient. However, in
terms of memory, the “naive” method is almost always preferred (unless $n\gg
p\cdot n_{\mathcal{A}}$ where $n_{\mathcal{A}}$ is the number of active
coefficients). This is because the “covariance” method requires saving
$X^{\top}WX_{g}\in\mathbb{R}^{p\times p_{g}}$ at the very least for the active
groups $g$ (with non-zero coefficients). With large $p$, this may be
unmanageable. While we could technically recompute $X^{\top}WX_{g}$ on-the-fly
if we have insufficient memory, this would significantly hurt the performance
and will fail to be competitive against the “naive” method. On the other hand,
the “naive” method only requires $O(n)$ extra memory for the residual.

###  3.4 Matrix Abstraction

The reader may notice that the BCD algorithm only ever interacts with the
matrix $X$ through a few operations, namely, inner-products with columns of
$X$. Our software adelie takes advantage of this fact and abstracts out the
matrix class. Any matrix-like class supporting a few methods to compute such
inner-products can be passed into our solver.

This abstraction has been proven to be quite useful for many applications. As
shown in Section 3.3, the bulk of the computation in each BCD update lies in
computing inner-products with $X$. Hence, any effort into speeding up these
inner-products will yield large performance gains for the solver. So, if $X$
is structured in some sense, it is extremely helpful to exploit its structure
to simplify the inner-product computations. A simple example is when $X$ is a
sparse matrix. A more motivating example is in genome-wide association studies
(GWAS) where the feature matrix of interest is the genotype matrix of SNP
information taking on values in $\left\\{0,1,2,\text{NA}\right\\}$. Both the
sparsity of this matrix as well as the restricted range of values allow us to
efficiently store this matrix with far more structure than a dense or a sparse
matrix. There is a plethora of references in this area, so we only mention a
few [Uffelmann et al., 2021; Witte, 2010; Mills and Rahal, 2019; He et al.,
2022, 2021, 2024; Qian et al., 2020; Li et al., 2020]. These genotype matrices
can be of daunting size. The UK biobank dataset [Sudlow et al., 2015], for
example, contains 500,000 samples and about 1.7 million SNPs (features).
Moreover, along with this genotype matrix, we also have a small set of dense
covariates such as age, sex, and 10 principal components of the genotype
matrix. Without the matrix abstraction, one must then represent the full
feature matrix $X$ as a dense matrix, which would require nearly $6$ terabytes
of floating point numbers, a size too large to load in the RAM of any modern
computers. Indeed, there has been some work to solve lasso on such large-scale
data [Qian et al., 2020; Li et al., 2020, 2021; He et al., 2024; Zeng and
Breheny, 2018], all implementing their own solvers for the purpose of
supporting their own representation of the large matrices. Using our method,
we can represent, for example, a dense matrix for the covariates, a structured
matrix for the genotype matrix, and a column-wise “concatenated” matrix of the
two. The sole purpose of the concatenated matrix is simply to delegate the
inner-products to the correct sub-matrix, which is an inexpensive operation.
If one wishes to represent the matrix differently, one only needs to implement
a small matrix class without implementing a separate solver.

We provide a last example specifically for the “covariance” method. The
covariance matrix $X^{\top}WX$ is usually a dense matrix, however, there are
cases when it can be approximated as a block-diagonal matrix [He et al.,
2024]. This structure was exploited to greatly speed up the lasso solver and
save large amounts of memory, however, they implemented a separate lasso
solver specifically assuming this structure. The same optimization can be done
with our software with few lines of code just to represent the block-diagonal
matrix. The user may otherwise reuse our general solver.

###  3.5 Pathwise Block-Coordinate Descent

We compute solutions along a decreasing sequence of the regularization
parameter $\lambda$. We start with the smallest $\lambda$, say
$\lambda_{\max}$, for which all groups with non-zero group lasso penalty have
coefficients identically $\mathbf{0}$. This value can be directly computed
once we fit all unpenalized variables. That is, if
$U:=\\{g\in[G]:\omega_{g}\alpha=0\\}$ denotes the set of groups that have
vanishing group lasso penalty, we first solve 1 only for groups in $U$ to get
their corresponding coefficients $\hat{\beta}_{U}$. Then, we compute
$\lambda_{\max}$ given by

$\displaystyle\lambda_{\max}:=\max_{g:g\notin
U}\frac{\left\|X_{g}^{\top}W(y-X_{U}\beta_{U})\right\|_{2}}{\alpha\omega_{g}}$
(17)

By the KKT conditions, we have our desired property that the zero solution on
groups not in $U$ is optimal at $\lambda_{\max}$. Note that if $U$ is empty
(which occurs if $\alpha=0$ or every $\omega_{g}=0$) so that every group is
included in the model, then we may arbitrarily take any value for
$\lambda_{\max}$. From $\lambda_{\max}$, we generate a decreasing sequence of
$K$ number of $\lambda$ values down to $\epsilon\lambda_{\max}$ evenly-spaced
on the log scale, where $\epsilon\in[0,1]$ is a user-specified value that
determines the depth of the path. Following Friedman et al. [2010], we
typically set $\epsilon\in[0.0001,0.01]$ and $K=100$.

The advantage of a pathwise solver is that the solver is warm-started at each
$\lambda$ (including $\lambda_{\max}$). This helps tremendously for smaller
values of $\lambda$ where a cold-start generally makes the solver wander
longer before reaching the optimum. This was also observed by Friedman et al.
[2010]. Similar ideas exist in the convex literature, in particular, with
interior-point methods and barrier methods [Boyd and Vandenberghe, 2004].

###  3.6 Screen and Active Sets

In the BCD algorithm, we can layer two types of optimizations that both reduce
computation by cycling over only a subset of the groups. The first
optimization is based on a _screening rule_ , which tries to guess the set of
groups that will likely contain the true active set (groups with non-zero
coefficients). We call such a set the _screen set_ and denote it as
$\mathcal{S}$. It is possible to modify the BCD algorithm to only ever iterate
over the screen set. Therefore, as mentioned in Section 3, we may only cache
the eigen-decomposition of $X_{g}^{\top}WX_{g}$ for $g\in\mathcal{S}$.
However, this optimization comes at a (small) price of checking the KKT
conditions. The second optimization further reduces most of the computation to
a subset of active groups in the screen set.

In the first optimization, the screen rule can be a heuristic or a guarantee.
The literature contains many types of such rules including SAFE rule, strong
rule, and EDPP rule [Ghaoui et al., 2011; Tibshirani et al., 2012; Wang et
al., 2015]. In particular, the SAFE rule and the EDPP rule are examples of
guarantees in the sense that at a fixed $\lambda$ if these rules discard
groups, then they guarantee that these groups will truly have zero
coefficients. On the other hand, the strong rule is an example of a heuristic,
which does not always guarantee that the groups that it discards will be
inactive. Despite the lack of a guarantee, we find that among these rules the
strong rule was the most effective one in practice as it includes much fewer
groups and is almost always correct. For this reason, we implement the strong
rule as part of our software.

We briefly describe the strong rule for the group elastic net problem. Suppose
we have a solution $\tilde{\beta}$ at $\tilde{\lambda}$ and the current screen
set is $\tilde{\mathcal{S}}$. We must now solve for $\lambda<\tilde{\lambda}$
to obtain the next solution $\beta$. For the group elastic net, the strong
rule searches for all remaining groups $g\notin\tilde{\mathcal{S}}$ such that

$\displaystyle\left\|X_{g}^{\top}W(y-X\tilde{\beta})\right\|_{2}<2\lambda-\tilde{\lambda}$
(18)

If this condition is satisfied, group $g$ is discarded from the screen set.
Note that $2\lambda-\tilde{\lambda}=\lambda-(\tilde{\lambda}-\lambda)$ so that
18 looks at the target $\lambda$ and a band below equal to the path gap at
$\lambda$. Once all such groups are discarded to obtain the next screen set
$\mathcal{S}$, we proceed to solve the group elastic net only over the groups
in $\mathcal{S}$, effectively solving the (constrained) group elastic net
problem with the constraint that $\beta_{-\mathcal{S}}=\mathbf{0}$; at the
start of the path, the screen set only contains unpenalized groups so that the
screen set is exactly the active set at $\lambda_{\max}$. It is easy to
imagine that the more accurate the screening rule is, the smaller the screen
set $\mathcal{S}$, and therefore, the less computation needed in the BCD
algorithm. Since we guess that groups outside the screen set are inactive, we
must verify the KKT conditions to confirm our guess. This amounts to checking
that

$\displaystyle\left\|X_{g}^{\top}W(y-X\beta)\right\|_{2}\leq\lambda,\quad\forall
g\notin\mathcal{S}$

If the condition holds, then we may declare $\beta$ as the solution at
$\lambda$.

The second optimization reduces the computation further by restricting most of
the BCD iterations to the current active set within the screen set. Namely, we
first iterate over the screen set once to determine the active groups.
Afterwards, if we have not converged, we iterate only over the active groups
until convergence. We then repeat the entire process until we reach
convergence on the screen set. This active-set strategy is also mentioned in
Friedman et al. [2010]; Meier et al. [2008]; Krishnapuram et al. [2005].

###  3.7 Convergence Criterion

In the BCD algorithm, we measure convergence by measuring the average amount
in which each block-update results in a change in the linear prediction and
take the largest such value. Concretely, suppose $\tilde{\beta}$ is the
current coefficient vector and we have just updated group $g$ so that the new
coefficient $\beta$ only differs from $\tilde{\beta}$ in the entries defined
by group $g$. Then, we compute
$p_{g}^{-1}\|X_{g}(\beta_{g}-\tilde{\beta}_{g})\|_{W}^{2}$. Note that since we
have already cached $X_{g}^{\top}WX_{g}=Q_{g}\Lambda_{g}Q_{g}^{\top}$, it is
equivalent to computing
$p_{g}^{-1}\|Q_{g}^{\top}(\beta_{g}-\tilde{\beta}_{g})\|_{\Lambda_{g}}^{2}$.
Since we must apply the rotation $Q_{g}^{\top}$ on $\tilde{\beta}_{g}$ and
$\beta_{g}$ anyways to solve 4, computing the convergence measure is an easy
task. After one full cycle of the groups, we check if

$\displaystyle\max\limits_{g}\quad
p_{g}^{-1}\left\|X_{g}(\beta_{g}-\tilde{\beta}_{g})\right\|_{W}^{2}\leq\epsilon$
(19)

where $\epsilon\geq 0$ is a user-specified tolerance level. If the condition
holds, we declare convergence.

###  3.8 Lasso Optimization

As mentioned before, groups of size $1$ have the same coordinate update as in
the lasso, which admits a closed-form solution. Hence, at every coordinate
update, our software checks for the group sizes and falls back to the lasso
update in this special case. Otherwise, we use our Newton-ABS method to solve
the update.

##  4 Algorithms for Regularized Generalized Linear Models

In Section 3, we discussed our main algorithm for solving the group elastic
net under the Gaussian loss. In this section, we extend our algorithm for
general convex losses.

Suppose $\ell(\cdot):\mathbb{R}^{n}\to\mathbb{R}$ is a convex, proper, and
closed function that is twice-continuously differentiable. We typically
consider $\ell(\cdot)$ as the negative log-likelihood of some GLM so that is
of the form

$\displaystyle\ell(\eta):=\sum\limits_{i=1}^{n}w_{i}\left(-y_{i}\eta_{i}+A(\eta_{i})\right)$

where $A(\cdot)$ is the log-partition function of an exponential family.
However, for the current discussion, we need not restrict ourselves to such
form. Our goal is to solve the more general optimization problem

$\displaystyle\operatorname*{minimize\;}_{(\beta_{0},\beta)\in\mathbb{R}^{1+p}}\quad$
$\displaystyle\ell(\eta)+\lambda P_{\alpha,\omega}(\beta)$ (20)
$\displaystyle\operatorname{subject\;to}\quad$
$\displaystyle\eta=X\beta+\beta_{0}\mathbf{1}+\eta^{0}$

where the group elastic net penalty $P_{\alpha,\omega}(\cdot)$ is defined by
2. The only new term is $\eta^{0}\in\mathbb{R}^{n}$, which we call the
_offset_ vector. The offset is a fixed vector that defines possibly different
baseline linear predictions for each observation so that
$X\beta+\beta_{0}\mathbf{1}$ is interpreted as the change in the linear
prediction from the given baseline. For the Gaussian loss

$\displaystyle\ell(\eta):=\sum\limits_{i=1}^{n}w_{i}\left(-y_{i}\eta_{i}+\frac{\eta_{i}^{2}}{2}\right)\equiv\frac{1}{2}\left(\left\|y-\eta\right\|_{W}^{2}-\left\|y\right\|_{W}^{2}\right)$

the offset can be managed in a trivial manner by replacing $y$ with
$y-\eta^{0}$ and running the algorithm described in Section 3 with no other
changes. However, with general loss functions, the offset can interact with
the objective in a non-trivial manner.

We use proximal quasi-Newton method to solve 20. This method is also mentioned
in Friedman et al. [2010]; Meier et al. [2008]. Moreover, Lee et al. [2014]
proves a superlinear convergence rate for proximal quasi-Newton methods (with
line search). However, as in Friedman et al. [2010], we do not perform the
line search since it will significantly impact the overall runtime and it does
not seem necessary in practice thanks to the warm-starts in our pathwise
solver. We observe fast convergence, usually only requiring 1 to 5 iterations.

We are now in position to describe the proximal quasi-Newton method for our
problem. We begin with an initial value
$(\beta_{0}^{(0)},\beta^{(0)})\in\mathbb{R}^{1+p}$. For each iteration $k\geq
0$ until convergence, we compute our current linear prediction
$\eta^{(k)}=X\beta^{(k)}+\beta_{0}^{(k)}\mathbf{1}+\eta^{0}$. We then perform
a second-order expansion of $\ell(\cdot)$ around $\eta^{(k)}$ and solve for
the proximal operator to get the next iterate
$(\beta_{0}^{(k+1)},\beta^{(k+1)})$. That is,

$\displaystyle(\beta_{0}^{(k+1)},\beta^{(k+1)})=\operatorname*{argmin}_{(\beta_{0},\beta)\in\mathbb{R}^{1+p}}\quad$
$\displaystyle\nabla\ell(\eta^{(k)})^{\top}\eta+\frac{1}{2}(\eta-\eta^{(k)})^{\top}W^{(k)}(\eta-\eta^{(k)})+\lambda
P_{\alpha,\omega}(\beta)$ (21)
$\displaystyle\quad\operatorname{subject\;to}\quad$
$\displaystyle\eta=X\beta+\beta_{0}\mathbf{1}+\eta^{0}$

Note that the linear term is indeed the gradient of $\ell(\cdot)$ at
$\eta^{(k)}$; however, $W^{(k)}$ need not be $\nabla^{2}\ell(\eta^{(k)})$. In
general, our strategy is to find a positive diagonal matrix $W^{(k)}$ that
majorizes $\nabla^{2}\ell(\eta^{(k)})$, i.e. $0\prec W^{(k)}$ diagonal and
$\nabla^{2}\ell(\eta^{(k)})\preceq W^{(k)}$. Of course, if
$\nabla^{2}\ell(\eta^{(k)})$ is already positive and diagonal, we may simply
take $W^{(k)}\equiv\nabla^{2}\ell(\eta^{(k)})$. The proximal operator 21 is
simply an instance of the group elastic net problem under the Gaussian loss
given by 1. Hence, we may use our algorithm in Section 3 to solve 21. We
summarize our algorithm in Algorithm 4. This algorithm is referred to as
_iteratively reweighted least squares_ (IRLS) in Friedman et al. [2010];
McCullagh and Nelder [1989].

Data: $X$, $y$, $\beta_{0}^{(0)}$, $\beta^{(0)}$, $\eta^{0}$

1 $k\leftarrow 0$;

2 while  _not converged based on 23_ do

3 $\eta^{(k)}=X\beta^{(k)}+\beta_{0}^{(k)}\mathbf{1}+\eta^{0}$;

4 Compute $\nabla\ell(\eta^{(k)})$;

5 Compute $W^{(k)}$ appropriate for the given loss $\ell(\cdot)$;

6 Compute $(\beta_{0}^{(k+1)},\beta^{(k+1)})$ by solving 21 using our method
in Section 3;

7 $k\leftarrow k+1$;

8

return $h$;

Algorithm 4 Proximal Quasi-Newton Method for General Smooth, Convex Losses

###  4.1 Application on Generalized Linear Models

Our primary application of Algorithm 4 is when $\ell(\cdot)$ represents the
negative log-likelihood of a GLM. Some examples include the Binomial family
with the logit or probit link function, Poisson family with the log link
function, and the Cox family. In terms of the solver, we do not make any
special changes based on the loss function. We note that this is one of the
more notable deviations from the approach of Friedman et al. [2010], as the
authors made amends to the solver specialized for each GLM, which we found
unnecessary.

The only manner in which the loss function plays a role in Algorithm 4 is in
computing the gradient $\nabla\ell(\eta^{(k)})$ and the hessian-like quantity
$W^{(k)}$. As a result, we have structured our software to abstract out the
loss function to a class that provides the gradient and the hessian, similar
to how we abstracted out the matrix classes (see Section 3.4). As an example,
for the Binomial family with the logit link function [McCullagh and Nelder,
1989] where the loss is given by

$\displaystyle\ell(\eta):=\sum\limits_{i=1}^{n}w_{i}\left(-y_{i}\eta_{i}+\log(1+\exp(\eta_{i}))\right)$

we implement a class to compute the gradient and hessian, respectively, given
by

$\displaystyle\nabla\ell(\eta)_{i}=-w_{i}(y_{i}-p(\eta_{i})),\quad\nabla^{2}\ell(\eta)_{ij}=\begin{cases}w_{i}p(\eta_{i})(1-p(\eta_{i})),&i=j\\\
0,&i\neq j\end{cases}$ (22)

where $p(t):=1/(1+\exp(-t))$ is the sigmoid function.

###  4.2 Screen Sets and the KKT Check

We may reuse the strong rule mentioned in Section 3.6 when solving 20. We
follow the notation as in Section 3.6. The only change we make is in computing
the scores so that we discard group $g\notin\tilde{\mathcal{S}}$ if

$\displaystyle\left\|X_{g}^{\top}\nabla\ell(X\tilde{\beta}+\tilde{\beta_{0}}\mathbf{1}+\eta^{0})\right\|_{2}<2\lambda-\tilde{\lambda},$

Similarly, we check the KKT condition by verifying that

$\displaystyle\left\|X_{g}^{\top}\nabla\ell(X\beta+\beta_{0}\mathbf{1}+\eta^{0})\right\|_{2}\leq\lambda,\quad\forall
g\notin\mathcal{S}$

###  4.3 Numerical Stability Issues

In 21, it is possible that $W^{(k)}$ is nearly semi-definite, which may lead
to divergence issues. For example, in the case of a Binomial family (logistic
regression), if the estimated probabilities
$p_{i}^{(k)}:=1/(1+\exp(-\eta_{i}^{(k)}))$ are near $0$ or $1$, then
$W_{ii}^{(k)}$ will be near $0$ (see 22). Hence, we majorize $W^{(k)}$ one
more time by taking the maximum of each diagonal entry with a small positive
number (e.g. $10^{-12}$) to guarantee positiveness of $W^{(k)}$. The authors
of Friedman et al. [2010] only make this change for the Binomial family, but
we see no reason why this should not be done for any other loss functions.

###  4.4 Convergence Criterion

We discuss our convergence criterion. At iteration $k\geq 0$, upon solving the
proximal operator 21, we check that

$\displaystyle\left|\left(\eta^{(k+1)}-\eta^{(k)}\right)^{\top}\left(\nabla\ell(\eta^{(k+1)})-\nabla\ell(\eta^{(k)})\right)\right|\leq\epsilon\cdot
n_{\mathcal{A}}$ (23)

and declare convergence if the condition holds. Here, $n_{\mathcal{A}}$
denotes the current number of coefficients (not groups) that were updated and
$\epsilon\geq 0$ is a user-specified tolerance level. We explain our intuition
behind this criterion below.

First, suppose we perform a first-order Taylor expansion of
$\nabla\ell(\cdot)$ around $\eta^{(k)}$,

$\displaystyle\nabla\ell(\eta)\approx\nabla\ell(\eta^{(k)})+\nabla^{2}\ell(\eta^{(k)})(\eta-\eta^{(k)})$

Then, we see that

$\displaystyle\left(\eta^{(k+1)}-\eta^{(k)}\right)^{\top}\left(\nabla\ell(\eta^{(k+1)})-\nabla\ell(\eta^{(k)})\right)\approx\left(\eta^{(k+1)}-\eta^{(k)}\right)^{\top}\nabla^{2}\ell(\eta^{(k)})\left(\eta^{(k+1)}-\eta^{(k)}\right)$
(24)

which measures the change in linear predictions weighted accordingly based on
the hessian at $\eta^{(k)}$. This way, a change in the prediction along a
direction for which the hessian is small is downweighted, as this suggests the
loss $\ell(\cdot)$ is only mildly affected by such a change. Conversely, a
change along a direction for which the hessian is large is upweighted, as this
suggests the loss $\ell(\cdot)$ may be largely affected.

Next, we argue that 23 is a generalization of the convergence metric in
Section 3.7. To see this, first replace $\eta=X\beta$, where we ignore the
intercept and the offset momentarily for simplicity. Then, from 24 we have
that

$\displaystyle\left(\eta^{(k+1)}-\eta^{(k)}\right)^{\top}\nabla^{2}\ell(\eta^{(k)})\left(\eta^{(k+1)}-\eta^{(k)}\right)=\Delta\beta^{\top}X^{\top}\nabla^{2}\ell(\eta^{(k)})X\Delta\beta$
(25)

where $\Delta\beta:=\beta^{(k+1)}-\beta^{(k)}$. Under the Gaussian loss where
$\nabla^{2}\ell(\eta^{(k)})\equiv W$, the right side of 25 becomes
$\left\|X\Delta\beta\right\|_{W}^{2}$. Now, suppose (as in the BCD algorithm)
that $\beta^{(k+1)}$ only differs in one group from $\beta^{(k)}$. That is,
for some $g\in[G]$, $\beta^{(k+1)}_{g}\neq\beta^{(k)}_{g}$ and
$\beta^{(k+1)}_{g^{\prime}}\equiv\beta^{(k)}_{g^{\prime}}$ for all
$g^{\prime}\neq g$. Then,
$\left\|X\Delta\beta\right\|_{W}^{2}=\left\|X_{g}\Delta\beta_{g}\right\|_{W}^{2}$.
In 19, we scale this quantity by $p_{g}$, which is the total number of
coefficients that were updated (since we have fixed all other coordinates).
Hence, we scale by the analogous quantity in 23.

##  5 Algorithms for Multi-Response Data

In this section, we discuss methods to fit multi-response data. In multi-
response data, we have a response _matrix_ $y\in\mathbb{R}^{n\times c}$ rather
than a response vector. As a result, we change the group elastic net objective
20 to

$\displaystyle\operatorname*{minimize\;}_{\beta_{0}\in\mathbb{R}^{c},B\in\mathbb{R}^{p\times
c}}\quad$ $\displaystyle\ell(\operatorname{vec}(\eta^{\top}))+\lambda
P_{\alpha,\omega}(\operatorname{vec}(B^{\top}))$ (26)
$\displaystyle\operatorname{subject\;to}\quad$
$\displaystyle\eta=XB+\mathbf{1}\beta_{0}^{\top}+\eta^{0}$

where $\operatorname{vec}(\cdot)$ is the operator that flattens a matrix into
a one-dimensional array by stacking the columns. Note that we must now fit a
matrix of coefficients along with an intercept vector
$\beta_{0}\in\mathbb{R}^{c}$ (one for each response).

While it is possible to specify any arbitrary group structure on $B$, we focus
only on two types: _grouped_ and _ungrouped_. The grouped method defines each
_row_ of $B$ as a group whereas the ungrouped method defines each
_coefficient_ of $B$ as a group. Hence, the grouped method yields a non-
trivial group lasso problem and the ungrouped method reduces to a lasso
problem. The penalty $P_{\alpha,\omega}$ simplifies to the following
respective forms

$\displaystyle
P_{\alpha,\omega}^{\text{grouped}}(\operatorname{vec}(B^{\top}))$
$\displaystyle=\sum\limits_{g=1}^{p}\omega_{g}\left(\alpha\left\|B_{g\cdot}\right\|_{2}+\frac{1-\alpha}{2}\left\|B_{g\cdot}\right\|_{2}^{2}\right)$
$\displaystyle
P_{\alpha,\omega}^{\text{ungrouped}}(\operatorname{vec}(B^{\top}))$
$\displaystyle=\sum\limits_{g=1}^{p}\omega_{g}\left(\alpha\left\|B_{g\cdot}\right\|_{1}+\frac{1-\alpha}{2}\left\|B_{g\cdot}\right\|_{2}^{2}\right)$

Although it may seem like 26 requires a different solver, we now demonstrate
that it can be easily solved with the tools built in Sections 3 and 4. It is
not difficult to show that the constraint in 26 can be rewritten as

$\displaystyle\operatorname{vec}(\eta^{\top})=(X\otimes
I_{c})\operatorname{vec}(B^{\top})+(\mathbf{1}\otimes
I_{c})\beta_{0}+\operatorname{vec}(\eta^{0\top})$ (27)

Written this way, with the data matrix $[(\mathbf{1}\otimes
I_{c})\quad(X\otimes I_{c})]\in\mathbb{R}^{nc\times(1+p)c}$ and the
coefficient vector
$(\beta_{0},\operatorname{vec}(B^{\top}))\in\mathbb{R}^{(1+p)c}$, we return to
the setting of 20 (with no intercept). As discussed in Section 3.4, as soon as
we implement a matrix to represent the Kronecker product with an identity
matrix, we are able to solve 26 with our existing tools. Such a matrix can be
implemented to compute inner-products efficiently due to its structure. For
the Multigaussian family where the loss function is given by

$\displaystyle\ell(\eta)=\sum\limits_{i=1}^{n}w_{i}\left(-\sum\limits_{j=1}^{c}y_{ij}\eta_{ij}+\frac{\left\|\eta_{i\cdot}\right\|_{2}^{2}}{2}\right)$

we may simply rely on our algorithm that solves the group elastic net under
the Gaussian loss (see 1). Otherwise, we must rely on our algorithm that
handles general convex losses (see 20).

We make a special remark on the Multinomial family where the loss is given by

$\displaystyle\ell(\eta):=\sum\limits_{i=1}^{n}w_{i}\left(-\sum\limits_{j=1}^{c}y_{ij}\eta_{ij}+\log\left(\sum\limits_{j=1}^{c}\exp(\eta_{ij})\right)\right)$

Unlike the GLMs discussed previously, this marks the first example where the
hessian $\nabla^{2}\ell(\eta)$ is not exactly positive and diagonal. Observe
that the model is not identifiable in the sense that for any fixed
$\left\\{\beta_{0,j},\beta_{\cdot,j}\right\\}_{j=1}^{c}$, we have that
$\left\\{\beta_{0,j}-d_{0},\beta_{\cdot,j}-d\right\\}_{j=1}^{c}$ yields the
same loss for any $d_{0}\in\mathbb{R}$ and $d\in\mathbb{R}^{p}$. As a result,
the hessian is never positive definite. Moreover, the hessian is only block-
diagonal. Consequently, Friedman et al. [2010] approach the Multinomial family
quite differently from the other GLMs. First, they tweak the solver to cycle
over the classes $j=1,\ldots,c$ while fixing all other class coefficients, so
that the gradient and hessian in 21 is only with respect to the linear
predictors for class $j$, namely, $\eta_{\cdot,j}$. Because of this
modification, given a candidate solution $(\beta_{0},\beta)$, they found it
necessary to then further minimize only the lasso penalty with respect to
$d_{0}$ and $d$ to improve the convergence speed. Despite these efforts, we
have found a real dataset on which their software glmnet fails to converge. We
pursue a different and much simpler approach. We simply need to construct a
positive, diagonal matrix that majorizes $\nabla^{2}\ell(\eta)$. Simon et al.
[2013b] show that $2\operatorname{diag}(\nabla^{2}\ell(\eta))$ is a valid
choice for dominating $\nabla^{2}\ell(\eta)$. Hence, we set
$W^{(k)}:=2\operatorname{diag}(\nabla^{2}\ell(\eta^{(k)}))$ in 21 and reuse
our work from Section 4. With this approach, our software is able to converge.
This method implicitly settles the identifiability issue precisely because
$W^{(k)}$ is now positive definite (or can be made positive definite as in
Section 4.3). Out of curiosity, we also experimented with minimizing the
penalty with respect to $d_{0}$ and $d$, but to our surprise, this resulted in
a less stable algorithm often diverging.

##  6 Benchmark

In this section, we show an assortment of benchmarks to demonstrate the power
of our method. We developed a highly efficient implementation of our
algorithms discussed in Sections 3, 4 and 5 in our package adelie where most
of the code is written in C++. Although we intend adelie to be a Python-first
package, we are currently developing a R package as
well222https://github.com/JamesYang007/adelie-r. We exported enough
functionalities to perform the benchmark comparisons in this section. We note
that both the Python and R packages share the same C++ code, which performs
most of the computation, and the only difference comes from the input
preparation for the C++ solver. Hence, there is only a negligible difference
in their runtimes. All benchmarks were run on an M1 Macbook Pro.

The packages we benchmark against are sparsegl [Liang et al., 2022; McDonald
et al., 2023], gglasso [Yang and Zou, 2015; Yang et al., 2024], and grpnet
[Helwig, 2024]. We note that these packages are pathwise-solvers as well. We
attempted to include the following packages, however, they either failed to
converge or had substantially longer runtimes: grplasso [Meier et al., 2008;
Meier, 2020], SGL [Simon et al., 2013a, 2019], and grpreg [Breheny and Huang,
2015].

For our real data analyses, we consider the following datasets:

  * •

Leukemia dataset [Golub et al., 1999]: the Leukemia dataset has been
previously described in Section 3. However, in this section, we will use it to
classify whether an individual has Leukemia.

  * •

UJIIndoorLoc dataset [Joaqun Torres-Sospedra, Raul Montoliu, Adolfo Martnez-
Us, Tomar Arnau, Joan Avariento, 2014]: the UJIIndoorLoc dataset contains WiFi
fingerprints of Received Signal Strength Intensity (RSSI) from location
trackers. We remove any redundant features from the dataset. Our task is to
predict the longitude coordinate of the sample.

  * •

Prostate dataset [Singh et al., 2002]: the Prostate dataset gene expressions
from a microarray study of prostate cancer. Our task is to classify whether an
individual has prostate cancer.

We summarize the datasets in Table 1.

Name | n | G | Task
---|---|---|---
Leukemia | 72 | 7129 | classification
UJIIndoorLoc | 19937 | 465 | regression
Prostate | 102 | 6033 | classification
Table 1: The number of samples ($n$) and groups ($G$) as well as the task for
each of the real datasets.

###  6.1 Group Lasso Comparisons

(a) Gaussian Loss

(b) Binomial Loss

Figure 4: Timing comparisons for solving the group lasso under the Gaussian
loss (Figure 4(a)) and the Binomial loss (Figure 4(b)) against existing R
packages. We study a small ($n=100$) and large ($n=1000$) sample size case,
and for each case we vary the (equi-)correlation of the features ($\rho$).
Figure 5: Timing comparisons for solving the group lasso for the real datasets
against existing R packages.

We show a comparison of solving the group lasso under the Gaussian and
Binomial losses using both simulated and real data.

We first describe our simulation setup. Our setup is similar to the one used
by Yuan and Lin [2006]. Let $n$ be the number of samples, $G$ be the number of
groups, and $\rho\in[0,1]$ be the equi-correlation between features. We
generate $Z_{i}\sim\mathcal{N}(0,I_{G})$ for $i\in[n]$ and
$W\sim\mathcal{N}(0,I_{n})$ all independent. We then define a sample
$Y_{i}=\sqrt{\rho}W_{i}\mathbf{1}+\sqrt{1-\rho}Z_{i}$ so that the features
have variance $1$ and equi-correlation $\rho$. Finally, we construct our
design matrix $X$ with a third-degree polynomial expansion of each $Y_{i,g}$
so that $X_{i,3g-2}=Y_{i,g}$, $X_{i,3g-1}=Y_{i,g}^{2}$, and
$X_{i,3g}=Y_{i,g}^{3}$ for all $i\in[n]$ and $g\in[G]$, so that
$X\in\mathbb{R}^{n\times 3G}$. Naturally, we define every 3 columns of $X$ as
a group. The natural parameter $\eta$ is generated from a linear model
$\eta=X\beta+\sigma\varepsilon$. The first 6 components of $\beta$ are
independently generated from a standard normal and the rest are set to 0. The
noise $\varepsilon\sim\mathcal{N}(0,I_{n})$ and $\sigma$ is such that the
signal-to-noise ratio is 3. For the Gaussian loss, our response $y$ is
identical to $\eta$, and for the Binomial loss, we generate $y$ as Bernoulli
random variables where the logit parameters are given by $\eta$. Finally,
prior to fitting, we center and standardize the columns of $X$ and $y$.

In Figure 4, we show the timings for solving the group lasso under the
Gaussian loss (Figure 4(a)) and the Binomial loss (Figure 4(b)). We study a
small ($n=100$) and large ($n=1000$) sample size case, and for each case we
vary the (equi-)correlation of the features ($\rho$). Each plot shows the
average times across $3$ trials to solve the same path of 100 regularization
values evenly-spaced on the log-scale starting at $\lambda_{\max}$, the
smallest $\lambda$ for which $\mathbf{0}$ is an optimal solution, until
$0.01\lambda_{\max}$. In all cases except one ($n=1000$, $\rho=0.95$,
$G=5000$, Binomial loss), adelie runs uniformly faster than the other
packages. adelie is up to 10 times faster than the next fastest package and,
more importantly, this difference factor remains to be the case as the data
scales both in $n$ and $G$. In the one case where adelie was not the fastest,
we observed that the screen set via the strong rule pulled in a large number
of groups from the very beginning. This resulted in a large number of
coordinate descent iterations and therefore a relatively long runtime. We are
currently working to improve the strong rule to ameliorate such cases.

We briefly mention that we tried the following setting with a larger group
size of $s=100$ rather than 3. In this setting, we generate
$Y\in\mathbb{R}^{n\times sG}$ as before such that every feature is Gaussian
with mean 0, variance 1, and equi-correlated with $\rho$. We then group every
$s$ features so that we have $G$ number of groups. The response vector $y$ is
generated in the same manner as before. Every package except adelie struggled
to converge as their objective values were unacceptably far from that of
adelie; we verified from the KKT conditions that adelie converged properly.
For this reason, we omit the timing comparison for this setting.

For the real datasets, we first expand each feature (column) $x$ into three
polynomial terms $x$, $x^{2}$, and $x^{3}$ just as with the Leukemia dataset
in Section 3 so that our feature matrix $X\in\mathbb{R}^{n\times 3G}$. Prior
to fitting, we center and standardize the columns of $X$ and $y$. In Figure 5,
we show the average timings to solve group lasso over 3 trials. Similar to the
simulation analysis, we solve for the same path of 100 regularization values
evenly-spaced on the log-scale for all packages. Once again, we observe a
dramatic speed improvement with adelie.

###  6.2 Lasso Comparisons

(a) Gaussian Loss

(b) Binomial Loss

Figure 6: Timing comparisons for solving the lasso under the Gaussian loss
(Figure 6(a)) and the Binomial loss (Figure 6(b)) against existing R packages.
We study a small ($n=100$) and large ($n=1000$) sample size case, and for each
case we vary the (equi-)correlation of the features ($\rho$).  Figure 7:
Timing comparisons for solving the lasso for the real datasets against
existing R packages.

In this section, we demonstrate that adelie is a competitive lasso solver. We
also include glmnet [Friedman et al., 2010] in the list of packages for
comparison. As in Section 6.1, we show a timings comparison of solving the
lasso under the Gaussian and Binomial losses using simulated and real data.

Our simulation setup is taken from Friedman et al. [2010]. For completion, we
describe it next. We begin by generating $n$ independent samples with $p$
features from a Gaussian distribution where every feature is equi-correlated
with $\rho$. The natural parameter $\eta=X\beta+\sigma\epsilon$ where
$\beta_{j}=(-1)^{j}\exp(-2(j-1)/20)$, $\epsilon\sim\mathcal{N}(0,1)$, and
$\sigma$ is chosen such that the signal-to-noise ratio is $3$. Note that the
coefficients have alternating signs and decrease exponentially in magnitude.
For the Gaussian loss, we let our response $y\equiv\eta$ and for the Binomial
loss, we generate $y$ from the Bernoulli distribution where the logit
parameter is given by $\eta$. Prior to fitting, we center and standardize the
columns of $X$ and $y$.

For the real datasets, we do not perform any special modifications to the data
as in the group lasso setting except for centering and standardizing the
columns of $X$ and $y$.

In Figures 6 and 7, we show the timings for solving the lasso under the
Gaussian loss (Figure 6(a)) and the Binomial loss (Figure 6(b)) for the
simulated and real data setting, respectively. We interpret the figures in the
same way as in Figures 4 and 5 for the group lasso comparisons in Section 6.1.
We note that adelie and glmnet have nearly identical timings and they clearly
outperform the other packages. This shows that adelie is not only a fast group
lasso solver, but also a fast lasso solver, suffering no loss of performance
compared to the highly optimized and specialized lasso solver glmnet.

##  7 Discussion

We demonstrate that block-coordinate descent is a powerful and effective
method for solving the group elastic net problem. We derive a novel algorithm
to solve the block-coordinate update for the Gaussian loss, outperforming
existing popular proximal gradient methods in both performance and accuracy.
Moreover, our method naturally generalizes to other convex losses via the
proximal (quasi)-Newton method. Through simulated and real data examples, we
demonstrate that our method is remarkably fast for both group lasso and lasso
problems.

Our Python package adelie is available under the MIT license from PyPI
(https://pypi.org/project/adelie) and GitHub
(https://github.com/JamesYang007/adelie). Our documentation is hosted on
GitHub Pages (https://jamesyang007.github.io/adelie). The user can install the
latest version by running pip install adelie in the terminal.

##  8 Acknowledgments

Trevor Hastie was partially supported by grants DMS-2013736 from the National
Science Foundation, and grant R01GM134483 from the National Institutes of
Health. The authors thank Jonathan Taylor for sharing his notes on the Cox
model and other helpful comments, and Balasubramanian Narasimhan and Rob
Tibshirani for helpful comments. We also thank Kevin Fry, Kevin Guo, and Anav
Sood for giving feedback on earlier drafts.

## References

  * Beck and Teboulle [2009] Beck, Amir and Teboulle, Marc.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  _SIAM Journal on Imaging Sciences_ , 2(1):183–202, 2009.  doi: 10.1137/080716542.  URL https://doi.org/10.1137/080716542.
  * Boyd and Vandenberghe [2004] Boyd, Stephen and Vandenberghe, Lieven.  _Convex optimization_.  Cambridge university press, 2004.
  * Breheny and Huang [2015] Breheny, Patrick and Huang, Jian.  Group descent algorithms for nonconvex penalized linear and logistic regression models with grouped predictors.  _Stat. Comput._ , 25(2):173–187, March 2015.
  * Brent [2013] Brent, Richard P.  _Algorithms for minimization without derivatives_.  Courier Corporation, 2013.
  * Dekker [1969] Dekker, Theodorus Jozef.  Finding a zero by means of successive linear interpolation.  _Constructive aspects of the fundamental theorem of algebra_ , pages 37–51, 1969.
  * Friedman et al. [2010] Friedman, Jerome, Hastie, Trevor, and Tibshirani, Robert.  Regularization paths for generalized linear models via coordinate descent.  _Journal of Statistical Software_ , 33(1), 2010.  ISSN 1548-7660.  doi: 10.18637/jss.v033.i01.  URL http://dx.doi.org/10.18637/jss.v033.i01.
  * Ghaoui et al. [2011] Ghaoui, Laurent El, Viallon, Vivian, and Rabbani, Tarek.  Safe feature elimination for the lasso and sparse supervised learning problems, 2011.
  * Golub et al. [1999] Golub, T. R., Slonim, D. K., Tamayo, P., Huard, C., Gaasenbeek, M., Mesirov, J. P., Coller, H., Loh, M. L., Downing, J. R., Caligiuri, M. A., Bloomfield, C. D., and Lander, E. S.  Molecular classification of cancer: Class discovery and class prediction by gene expression monitoring.  _Science_ , 286(5439):531–537, October 1999.  ISSN 1095-9203.  doi: 10.1126/science.286.5439.531.  URL http://dx.doi.org/10.1126/science.286.5439.531.
  * Hastie et al. [2016] Hastie, Trevor, Tibshirani, Rob, and Wainwright, Martin.  _Statistical Learning with Sparsity_.  Chapman & Hall, 2016.
  * He et al. [2021] He, Zihuai, Le Guen, Yann, Liu, Linxi, Lee, Justin, Ma, Shiyang, Yang, Andrew C., Liu, Xiaoxia, Rutledge, Jarod, Losada, Patricia Moran, Song, Bowen, Belloy, Michael E., Butler, Robert R., Longo, Frank M., Tang, Hua, Mormino, Elizabeth C., Wyss-Coray, Tony, Greicius, Michael D., and Ionita-Laza, Iuliana.  Genome-wide analysis of common and rare variants via multiple knockoffs at biobank scale, with an application to alzheimer disease genetics.  _The American Journal of Human Genetics_ , 108(12):2336–2353, December 2021.  ISSN 0002-9297.  doi: 10.1016/j.ajhg.2021.10.009.  URL http://dx.doi.org/10.1016/j.ajhg.2021.10.009.
  * He et al. [2022] He, Zihuai, Liu, Linxi, Belloy, Michael E., Le Guen, Yann, Sossin, Aaron, Liu, Xiaoxia, Qi, Xinran, Ma, Shiyang, Gyawali, Prashnna K., Wyss-Coray, Tony, Tang, Hua, Sabatti, Chiara, Candès, Emmanuel, Greicius, Michael D., and Ionita-Laza, Iuliana.  Ghostknockoff inference empowers identification of putative causal variants in genome-wide association studies.  _Nature Communications_ , 13(1), November 2022.  ISSN 2041-1723.  doi: 10.1038/s41467-022-34932-z.  URL http://dx.doi.org/10.1038/s41467-022-34932-z.
  * He et al. [2024] He, Zihuai, Chu, Benjamin, Yang, James, Gu, Jiaqi, Chen, Zhaomeng, Liu, Linxi, Morrison, Tim, Belloy, Michael E., Qi, Xinran, Hejazi, Nima, Mathur, Maya, Le Guen, Yann, Tang, Hua, Hastie, Trevor, Ionita-laza, Iuliana, Sabatti, Chiara, and Candès, Emmanuel.  In silico identification of putative causal genetic variants.  March 2024.  doi: 10.1101/2024.02.28.582621.  URL http://dx.doi.org/10.1101/2024.02.28.582621.
  * Helwig [2024] Helwig, Nathaniel E.  _grpnet: Group Elastic Net Regularized GLMs and GAMs_ , 2024.  URL https://CRAN.R-project.org/package=grpnet.  R package version 0.3.
  * Joaqun Torres-Sospedra, Raul Montoliu, Adolfo Martnez-Us, Tomar Arnau, Joan Avariento [2014] Joaqun Torres-Sospedra, Raul Montoliu, Adolfo Martnez-Us, Tomar Arnau, Joan Avariento.  UJIIndoorLoc, 2014.
  * Klosa et al. [2020] Klosa, J., Simon, N., Westermark, P.O., Liebscher, V., and Wittenburg, D.  Seagull: lasso, group lasso and sparse-group lasso regularization for linear regression models via proximal gradient descent.  _BMC Bioinformatics_ , 21(407), 2020.  doi: 10.1186/s12859-020-03725-w.
  * Krishnapuram et al. [2005] Krishnapuram, B., Carin, L., Figueiredo, M.A.T., and Hartemink, A.J.  Sparse multinomial logistic regression: fast algorithms and generalization bounds.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 27(6):957–968, June 2005.  ISSN 2160-9292.  doi: 10.1109/tpami.2005.127.  URL http://dx.doi.org/10.1109/TPAMI.2005.127.
  * Lee et al. [2014] Lee, Jason D., Sun, Yuekai, and Saunders, Michael A.  Proximal newton-type methods for minimizing composite functions.  _SIAM Journal on Optimization_ , 24(3):1420–1443, January 2014.  ISSN 1095-7189.  doi: 10.1137/130921428.  URL http://dx.doi.org/10.1137/130921428.
  * Li et al. [2020] Li, Ruilin, Chang, Christopher, Justesen, Johanne M, Tanigawa, Yosuke, Qian, Junyang, Hastie, Trevor, Rivas, Manuel A, and Tibshirani, Robert.  Fast lasso method for large-scale and ultrahigh-dimensional cox model with applications to uk biobank.  _Biostatistics_ , 23(2):522–540, September 2020\.  ISSN 1468-4357.  doi: 10.1093/biostatistics/kxaa038.  URL http://dx.doi.org/10.1093/biostatistics/kxaa038.
  * Li et al. [2021] Li, Ruilin, Chang, Christopher, Tanigawa, Yosuke, Narasimhan, Balasubramanian, Hastie, Trevor, Tibshirani, Robert, and Rivas, Manuel A.  Fast numerical optimization for genome sequencing data in population biobanks.  _Bioinformatics_ , 37(22):4148–4155, June 2021\.  ISSN 1367-4811.  doi: 10.1093/bioinformatics/btab452.  URL http://dx.doi.org/10.1093/bioinformatics/btab452.
  * Liang et al. [2022] Liang, Xiaoxuan, Cohen, Aaron, Heinsfeld, Anibal Solón, Pestilli, Franco, and McDonald, Daniel J.  sparsegl: An r package for estimating sparse group lasso, 2022.  URL https://arxiv.org/abs/2208.02942.
  * Loris [2009] Loris, Ignace.  On the performance of algorithms for the minimization of $\ell_{1}$-penalized functionals.  _Inverse Problems_ , 25(3):035008, jan 2009.  doi: 10.1088/0266-5611/25/3/035008.  URL https://dx.doi.org/10.1088/0266-5611/25/3/035008.
  * McCullagh and Nelder [1989] McCullagh, P. and Nelder, J. A.  _Generalized Linear Models_.  Springer US, 1989.  ISBN 9781489932426.  doi: 10.1007/978-1-4899-3242-6.  URL http://dx.doi.org/10.1007/978-1-4899-3242-6.
  * McDonald et al. [2023] McDonald, Daniel J., Liang, Xiaoxuan, Solón Heinsfeld, Anibal, and Cohen, Aaron.  _sparsegl: Sparse Group Lasso_ , 2023.  URL https://CRAN.R-project.org/package=sparsegl.  R package version 1.0.2.
  * Meier [2020] Meier, Lukas.  _grplasso: Fitting User-Specified Models with Group Lasso Penalty_ , 2020.  URL https://CRAN.R-project.org/package=grplasso.  R package version 0.4-7.
  * Meier et al. [2008] Meier, Lukas, Van De Geer, Sara, and Böhlmann, Peter.  The group lasso for logistic regression.  _Journal of the Royal Statistical Society B_ , 70(1):53–71, 2008.  doi: 10.1111/j.1467-9868.2007.00627.x.
  * Mills and Rahal [2019] Mills, Melinda C. and Rahal, Charles.  A scientometric review of genome-wide association studies.  _Communications Biology_ , 2(1), January 2019.  ISSN 2399-3642.  doi: 10.1038/s42003-018-0261-x.  URL http://dx.doi.org/10.1038/s42003-018-0261-x.
  * O’Donoghue and Candés [2015] O’Donoghue, Brendan and Candés, Emmanuel.  Adaptive restart for accelerated gradient schemes.  _Foundations of Computational Mathematics_ , 15:715–732, 2015.  doi: https://doi-org.stanford.idm.oclc.org/10.1007/s10208-013-9150-3.
  * Qian et al. [2020] Qian, Junyang, Tanigawa, Yosuke, Du, Wenfei, Aguirre, Matthew, Chang, Chris, Tibshirani, Robert, Rivas, Manuel A., and Hastie, Trevor.  A fast and scalable framework for large-scale and ultrahigh-dimensional sparse regression with application to the uk biobank.  _PLOS Genetics_ , 16(10):e1009141, October 2020\.  ISSN 1553-7404.  doi: 10.1371/journal.pgen.1009141.  URL http://dx.doi.org/10.1371/journal.pgen.1009141.
  * Simon and Tibshirani [2012] Simon, Noah and Tibshirani, Robert.  Standardization and the group lasso penalty.  _Statistica Sinica_ , 22(3):983–1001, 2012.  ISSN 10170405, 19968507.  URL http://www.jstor.org/stable/24309971.
  * Simon et al. [2013a] Simon, Noah, Friedman, Jerome, Hastie, Trevor, and Tibshirani, Robert.  A Sparse-Group lasso.  _J. Comput. Graph. Stat._ , 22(2):231–245, April 2013a.
  * Simon et al. [2019] Simon, Noah, Friedman, Jerome, Hastie, Trevor, and Tibshirani, Rob.  _SGL: Fit a GLM (or Cox Model) with a Combination of Lasso and Group Lasso Regularization_ , 2019.  URL https://CRAN.R-project.org/package=SGL.  R package version 1.3.
  * Simon et al. [2013b] Simon, Noah R., Friedman, Jerome H., and Hastie, Trevor J.  A blockwise descent algorithm for group-penalized multiresponse and multinomial regression.  _arXiv: Computation_ , 2013b.  URL https://api.semanticscholar.org/CorpusID:18386302.
  * Singh et al. [2002] Singh, Dinesh, Febbo, Phillip G, Ross, Kenneth, Jackson, Donald G, Manola, Judith, Ladd, Christine, Tamayo, Pablo, Renshaw, Andrew A, D’Amico, Anthony V, Richie, Jerome P, Lander, Eric S, Loda, Massimo, Kantoff, Philip W, Golub, Todd R, and Sellers, William R.  Gene expression correlates of clinical prostate cancer behavior.  _Cancer Cell_ , 1(2):203–209, March 2002.
  * Sudlow et al. [2015] Sudlow, Cathie, Gallacher, John, Allen, Naomi, Beral, Valerie, Burton, Paul, Danesh, John, Downey, Paul, Elliott, Paul, Green, Jane, Landray, Martin, Liu, Bette, Matthews, Paul, Ong, Giok, Pell, Jill, Silman, Alan, Young, Alan, Sprosen, Tim, Peakman, Tim, and Collins, Rory.  Uk biobank: An open access resource for identifying the causes of a wide range of complex diseases of middle and old age.  _PLOS Medicine_ , 12(3):e1001779, March 2015.  ISSN 1549-1676.  doi: 10.1371/journal.pmed.1001779.  URL http://dx.doi.org/10.1371/journal.pmed.1001779.
  * Tibshirani [1996] Tibshirani, Robert.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 58(1):267–288, 1996.  ISSN 00359246.  URL http://www.jstor.org/stable/2346178.
  * Tibshirani et al. [2012] Tibshirani, Robert, Bien, Jacob, Friedman, Jerome, Hastie, Trevor, Simon, Noah, Taylor, Jonathan, and Tibshirani, Ryan J.  Strong rules for discarding predictors in lasso-type problems.  _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ , 74(2):245–266, 2012.  ISSN 13697412, 14679868.  URL http://www.jstor.org/stable/41430939.
  * Tseng [2001] Tseng, P.  Convergence of a block coordinate descent method for nondifferentiable minimization.  _Journal of Optimization Theory and Applications_ , 109:475–494, 2001.  doi: https://doi-org.stanford.idm.oclc.org/10.1023/A:1017501703105.
  * Uffelmann et al. [2021] Uffelmann, Emil, Huang, Qin Qin, Munung, Nchangwi Syntia, de Vries, Jantina, Okada, Yukinori, Martin, Alicia R., Martin, Hilary C., Lappalainen, Tuuli, and Posthuma, Danielle.  Genome-wide association studies.  _Nature Reviews Methods Primers_ , 1(1), August 2021.  ISSN 2662-8449.  doi: 10.1038/s43586-021-00056-9.  URL http://dx.doi.org/10.1038/s43586-021-00056-9.
  * Wang et al. [2015] Wang, Jie, Wonka, Peter, and Ye, Jieping.  Lasso screening rules via dual polytope projection.  _Journal of Machine Learning Research_ , 16(32):1063–1101, 2015.  URL http://jmlr.org/papers/v16/wang15a.html.
  * Witte [2010] Witte, John S.  Genome-wide association studies and beyond.  _Annual Review of Public Health_ , 31(1):9–20, March 2010.  ISSN 1545-2093.  doi: 10.1146/annurev.publhealth.012809.103723.  URL http://dx.doi.org/10.1146/annurev.publhealth.012809.103723.
  * Wright et al. [2009] Wright, Stephen J., Nowak, Robert D., and Figueiredo, MÁrio A. T.  Sparse reconstruction by separable approximation.  _IEEE Transactions on Signal Processing_ , 57(7):2479–2493, 2009.  doi: 10.1109/TSP.2009.2016892.
  * Yang and Zou [2015] Yang, Yi and Zou, Hui.  A fast unified algorithm for solving group-lasso penalize learning problems.  _Stat. Comput._ , 25(6):1129–1141, November 2015\.
  * Yang et al. [2024] Yang, Yi, Zou, Hui, and Bhatnagar, Sahir.  _gglasso: Group Lasso Penalized Learning Using a Unified BMD Algorithm_ , 2024.  URL https://CRAN.R-project.org/package=gglasso.  R package version 1.5.1.
  * Yuan and Lin [2006] Yuan, Ming and Lin, Yi.  Model selection and estimation in regression with grouped variables.  _Journal of the Royal Statistical Society B_ , 68(1):49–67, 2006.  doi: 10.1111/j.1467-9868.2005.00532.x.
  * Zeng and Breheny [2018] Zeng, Yaohui and Breheny, Patrick.  The biglasso package: A memory- and computation-efficient solver for lasso model fitting with big data in r, 2018.
  * Zou and Hastie [2005] Zou, Hui and Hastie, Trevor.  Regularization and variable selection via the elastic net.  _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ , 67(2):301–320, 2005.  ISSN 13697412, 14679868.  URL http://www.jstor.org/stable/3647580.

##  Appendix A Proof of Theorem 3.1

Consider the optimization problem 5 with $\lambda>0$. Suppose that
$v_{S}=\mathbf{0}$ where $S$ is given in 6. Without loss of generality, we may
assume $S^{c}$ is non-empty so that there exists at least one $i$ such that
$\Sigma_{ii}>0$. Otherwise, by hypothesis, we must have $v=\mathbf{0}$ so that
5 reduces to minimizing $\left\|x\right\|_{2}$, which is uniquely minimized at
$x^{\star}=\mathbf{0}$.

By the KKT conditions, the (finite) minimizer $x^{\star}$ must satisfy

$\displaystyle\Sigma x-v+\lambda u=\mathbf{0}$

where $u$ is a sub-gradient of $\left\|\cdot\right\|_{2}$ at $x$. We first
prove the identity 7. If $x^{\star}=\mathbf{0}$, then we must have that
$\left\|v\right\|_{2}\leq\lambda$. Otherwise, $u\equiv x/\left\|x\right\|_{2}$
so that

$\displaystyle
x=\left(\Sigma+\frac{\lambda}{\left\|x\right\|_{2}}I\right)^{-1}v$ (28)

It remains to show that if $\left\|v\right\|_{2}\leq\lambda$, then
$x^{\star}=\mathbf{0}$, which will establish the identity 7.

In that endeavor, suppose $x^{\star}\neq\mathbf{0}$. Then, taking norms on
both sides of 28, we have that $x^{\star}$ satisfies

$\displaystyle\sum\limits_{i=1}^{d}\frac{v_{i}^{2}}{(\Sigma_{ii}\left\|x\right\|_{2}+\lambda)^{2}}=1$

Note that we then have that $\left\|v\right\|_{2}^{2}\lambda^{-2}\geq 1$, or
equivalently, $\left\|v\right\|_{2}\geq\lambda$. There must exist some $i$
such that $v_{i}\neq 0$ and $\Sigma_{ii}>0$. Otherwise, $v\equiv\mathbf{0}$ by
hypothesis that $v_{S}=\mathbf{0}$, which contradicts
$\left\|v\right\|_{2}\geq\lambda>0$. Consequently, for the element $i$ such
that $v_{i}\neq 0$ and $\Sigma_{ii}>0$, we have the strict inequality

$\displaystyle\frac{v_{i}^{2}}{(\Sigma_{ii}\left\|x\right\|_{2}+\lambda)^{2}}<v_{i}^{2}\lambda^{-2}$

Hence, we have that $\left\|v\right\|_{2}>\lambda$. This establishes the
equivalence $x^{\star}=\mathbf{0}\iff\left\|v\right\|_{2}\leq\lambda$.

##  Appendix B Proof of Theorem 3.2

Consider any $f$ satisfying the hypothesis of Theorem 3.2. Let $x^{(0)}$ be
any initial point such that $f(x^{(0)})\geq 0$. For a given $k\geq 0$, suppose
$f(x^{(k)})\geq 0$. Then, by convexity of $f$ and the definition of
$x^{(k+1)}$,

$\displaystyle f(x^{(k+1)})$ $\displaystyle\geq
f(x^{(k)})+f^{\prime}(x^{(k)})(x^{(k+1)}-x^{(k)})=0$

By induction, we have that $f(x^{(k)})\geq 0$ for every $k\geq 0$. Since $f$
is decreasing and has a non-vanishing derivative, we must have that
$f^{\prime}(x)<0$ for all $x$. In particular, $x^{(k+1)}\geq x^{(k)}$ so that
it is an increasing sequence. Note that if $f(x^{(k)})=0$ for some $k$, then
$x^{(k+j)}=x^{(k)}$ for all $j\geq 0$ so that the sequence converges in a
trivial fashion to a root of $f$. Hence, without loss of generality, suppose
$f(x^{(k)})>0$ for all $k\geq 0$. By the decreasing property of $f$, we must
then have that $x^{(k)}\leq r$ for every $k\geq 0$. Since the sequence is
increasing and bounded above, it converges to some value $x^{(\infty)}$. By
continuity, taking $k\to\infty$ in 11,

$\displaystyle
x^{(\infty)}=x^{(\infty)}-\frac{f(x^{(\infty)})}{f^{\prime}(x^{(\infty)})}\implies
f(x^{(\infty)})=0$

which shows that $x^{(\infty)}$ is a root of $f$.

##  Appendix C Properties of $\varphi$ in 9

Consider $\varphi$ as in 9. Recall by Theorem 3.1, we only consider $\varphi$
when $\left\|v\right\|_{2}>\lambda$. Similar to the proof of Theorem 3.1 in
Appendix A, we have that there exists some $i$ such that $v_{i}\neq 0$ and
$\Sigma_{ii}>0$. Then, differentiating $\varphi$, we get that

$\displaystyle\varphi^{\prime}(h)=-2\sum_{i=1}^{d}\frac{v_{i}^{2}\Sigma_{ii}}{(\Sigma_{ii}h+\lambda)^{3}}<0,\quad\varphi^{\prime\prime}(h)=6\sum_{i=1}^{d}\frac{v_{i}^{2}\Sigma_{ii}^{2}}{(\Sigma_{ii}h+\lambda)^{4}}>0$

This shows that $\varphi$ is strictly decreasing and strictly convex.

Next, rewriting $\varphi$ as

$\displaystyle\varphi(h)=\sum_{i\notin
S}\frac{v_{i}^{2}}{(\Sigma_{ii}h+\lambda)^{2}}-1$ (29)

we see that $\lim_{h\to\infty}\varphi(h)<0$. Note that
$\varphi(0)=\left\|v\right\|_{2}^{2}/\lambda^{2}-1>0$. Hence, a finite
positive root uniquely exists for $\varphi$.

##  Appendix D Derivation of $h_{\star}$ and $h^{\star}$

We begin with the lower bound $h_{\star}$. Let $a(h),b(h)\in\mathbb{R}^{d}$ be
defined by $a_{k}(h):=\Sigma_{kk}h+\lambda$ and
$b_{k}(h):=\left|v_{k}\right|/a_{k}(h)$ for each $k=1,\ldots,d$. Then, by
Cauchy-Schwarz,

$\displaystyle\left\|v\right\|_{1}$
$\displaystyle:=\sum\limits_{k=1}^{d}\left|v_{k}\right|=\sum\limits_{k=1}^{d}a_{k}(h)b_{k}(h)\leq\left\|a(h)\right\|_{2}\left\|b(h)\right\|_{2}$

Hence, if $\left\|a(h)\right\|_{2}\leq\left\|v\right\|_{1}$, then
$\left\|b(h)\right\|_{2}\geq 1$, or equivalently, $\varphi(h)\geq 0$. We see
that $\left\|a(h)\right\|_{2}\leq\left\|v\right\|_{1}$ if and only if

$\displaystyle\sum\limits_{i=1}^{d}(\Sigma_{ii}h+\lambda)^{2}\leq\left\|v\right\|_{1}^{2}$

By the quadratic formula, we may choose the largest $h$ that satisfies the
above condition, which is given by 14.

The upper bound $h^{\star}$ is even easier. Since

$\displaystyle\sum\limits_{i=1}^{d}\frac{v_{i}^{2}}{(\Sigma_{ii}h+\lambda)^{2}}$
$\displaystyle=\sum\limits_{i:\Sigma_{ii}>0}\frac{v_{i}^{2}}{(\Sigma_{ii}h+\lambda)^{2}}\leq
h^{-2}\sum\limits_{i:\Sigma_{ii}>0}\frac{v_{i}^{2}}{\Sigma_{ii}^{2}}$

we may set $h^{\star}$ to 15. Then, we immediately have that
$\varphi(h^{\star})\leq 0$.