#  sparsegl: An R Package for Estimating Sparse Group Lasso

Xiaoxuan Liang  Department of Statistics, The University of British Columbia,
Vancouver, BC Canada  Aaron Cohen  Department of Statistics, Indiana
University, Bloomington, IN USA  Anibal Solón Heinsfeld  Department of
Computer Science, The University of Texas, Austin TX USA  Franco Pestilli
Department of Psychology, The University of Texas, Austin TX, USA  Center for
Perceptual Systems, The University of Texas, Austin TX, USA  Daniel J.
McDonald  Department of Statistics, The University of British Columbia,
Vancouver, BC Canada

###### Abstract

The sparse group lasso is a high-dimensional regression technique that is
useful for problems whose predictors have a naturally grouped structure and
where sparsity is encouraged at both the group and individual predictor level.
In this paper we discuss a new R package for computing such regularized
models. The intention is to provide highly optimized solution routines
enabling analysis of very large datasets, especially in the context of sparse
design matrices.

11footnotetext: To whom correspondence should be addressed. E-mail:
<EMAIL_ADDRESS>
##  1 Introduction

Regularized linear models are now ubiquitous tools for prediction and,
increasingly, inference. When solving such high-dimensional learning problems,
adding regularization helps to reduce the chances of overfitting and improve
the model performance on unseen data. Sparsity inducing $\ell_{1}$-type
penalties such as the lasso (Tibshirani, 1996) or the Dantzig selector (Candes
and Tao, 2007) perform both variable selection and shrinkage, resulting in
near-optimal statistical properties. The group lasso (Yuan and Lin, 2006)
modifies the regularizer, replacing the $\ell_{1}$ penalty with a groupwise
sum of $\ell_{2}$ norms. When covariates have natural groupings, such as with
genomics data or one-hot encoded factors, this group penalty is preferable,
because the resulting estimate will include or exclude entire groups of
covariates. To simultaneously attain sparsity at both group and individual
feature levels, Simon et al. (2013) proposed the sparse group-lasso, a convex
combination of the $\ell_{1}$ lasso penalty and the group lasso penalty.

While a number of packages exist for solving the sparse group lasso, our R (R
Core Team, 2022) implementation in sparsegl (McDonald et al., 2022) is
designed to be fast, especially in the case of large, sparse covariate
matrices. This package focuses on finding the optimal solutions to sparse
group-lasso penalized learning problems at a sequence of regularization
parameters, implements risk estimators in an effort to avoid cross validation
if necessary, leverages a fast, compiled Fortran implementation, avoids
extraneous data copies, and undertakes a number of additional computational
efficiency improvements. In R, there are already excellent implementations of
sparse group lasso and group lasso, namely SGL (Simon et al., 2019), gglasso
(Yang et al., 2020; Yang and Zou, 2015), and biglasso (Zeng et al., 2022; Zeng
and Breheny, 2020). Of these, only SGL employs the additional $\ell_{1}$
sparsity-inducing penalty. However, it has a number of drawbacks that result
in much slower performance, even on small data. One major reason is the
omission of so-called “strong rules” (Tibshirani et al., 2012) that help
coordinate descent algorithms to avoid many of the groups which will turn out
to have zero coefficient estimates. The gglasso and biglasso packages are both
computationally fast. The former incorporates the strong rule, and the latter
involves a hybrid safe-strong rule along with scalable storage and a parallel
implementation in C++ and R that allows for data that exceeds the size of
installed random access memory. Unfortunately, neither allows within-group
sparsity (i.e., they perform group lasso, not sparse group lasso). Thus, the
estimated coefficients produced by these packages will have some active groups
and some inactive groups, but within an active group, generally all the
coefficients will be nonzero.

In python (Van Rossum et al., 2011), asgl (Civieta et al., 2021) implements
adaptive sparse group-lasso, which flexibly adjusts the weights in the
penalization on the groups of features. Additionally it incorporates quantile
loss. As with the other packages mentioned above, it can also solve the
special cases (lasso and group lasso). However, for all optimization problems,
it directly uses CVXPY (Diamond and Boyd, 2016; Agrawal et al., 2018), a
general purpose optimizer, without strong rules or other tricks to relate
solutions to each other across values of the tuning parameter.111A similar
implementation could be achieved in R using CVXR (Fu et al., 2020, 2022).

Our contribution, then, is to provide a package that performs sparse group
lasso and is faster than existing implementations. In particular, sparsegl has
the following benefits:

  * •

Performs Gaussian and logistic regression using fast, compiled code;

  * •

Allows arbitrary generalized linear models using R’s family object, though
with slightly less efficiency;

  * •

Allows for interval constraints and differential weights on the coefficients;

  * •

Accommodates a sparse design matrix and returns the coefficient estimates in a
sparse matrix;

  * •

Uses strong rules (and active set iteration) for fast computation along a
sequence of tuning parameters;

  * •

Uses dotCall64 to interface with low-level Fortran functions and avoid
unnecessary copying as well as allow for 64-bit integers (see Gerber et al.,
2017, 2018);

  * •

Provides information criteria as risk estimators (AIC/BIC/GCV) in addition to
cross validation.

A comparison of features of this and related R packages is shown in Table 1.
Figure 1 compares the speed of the dense and sparse implementations in
sparsegl with SGL across a number of different problem sizes, finding speedups
of 1.5 to 2.5 orders of magnitude.

In Section 2, we describe the algorithmic implementation in detail, paying
particular attention to the strong rule. In Section 3, we show how to use the
package, running through an example with simulated data. Section 4
demonstrates many of the unique features of sparsegl in two applications. We
summarize our contributions in Section 5.

| Regression & | Within group | Sparse | Strong | Avoids | Interval
---|---|---|---|---|---|---
| classification | sparsity | matrices | rules | copies | constraints
sparsegl | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
gglasso | ✓ |  |  | ✓ |  |
SGL | ✓ | ✓ |  |  |  |
Table 1: This table summarizes the features available in sparsegl and related
R packages. Figure 1: This figure shows the time required to compute sparse
group lasso solutions across a number of different problem sizes. In all
cases, we use $n=500$ observations and 100 values of the tuning parameter
$\lambda$. The median is taken across 5 replications for each method and
problem size. Note that both axes are on the log scale.

##  2 Methodology, estimation and prediction

Given a sample of $n$ observations of a univariate response $y_{i}$ and a
corresponding vector of features $\mathbf{x}_{i}\in\mathbb{R}^{p}$, the
standard linear regression setup has

$y_{i}=\mathbf{x}_{i}^{\mathsf{T}}\boldsymbol{\beta}+\sigma\epsilon_{i},\ \
i=1,\ldots,n,$

where $\epsilon_{i}$ is independent standard Gaussian noise and $\sigma>0$.
While ordinary least squares estimates the coefficient vector
$\boldsymbol{\beta}$ by solving
$\min_{\boldsymbol{\beta}}\frac{1}{2n}\sum_{i=1}^{n}(y_{i}-\mathbf{x}_{i}^{\mathsf{T}}\boldsymbol{\beta})^{2}$,
this method tends to behave poorly if $p\gg n$. In what follows, we will write
$\mathbf{y}=(y_{1},\ldots,y_{n})$ and let $\mathbf{X}$ be the rowwise
concatenation of
$\mathbf{x}^{\mathsf{T}}_{1},\ldots,\mathbf{x}_{n}^{\mathsf{T}}$.

The lasso adds an $\ell_{1}$ penalty to the optimization problem:

$\min_{\boldsymbol{\beta}}\left\\{\frac{1}{2n}\sum_{i=1}^{n}(y_{i}-\mathbf{x}_{i}^{\mathsf{T}}\boldsymbol{\beta})^{2}+\lambda\sum_{j=1}^{p}|\beta_{j}|\right\\}=\min_{\boldsymbol{\beta}}\left\\{\frac{1}{2n}\left\lVert\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\right\rVert_{2}^{2}+\lambda\left\lVert\boldsymbol{\beta}\right\rVert_{1}\right\\},$
(1)

where $\left\lVert\cdot\right\rVert_{2}$ is the Euclidean ($\ell_{2}$) norm
and $\left\lVert\cdot\right\rVert_{1}$ is the $\ell_{1}$ norm. The benefit of
this penalty is that it tends to allow only a subset of coefficient estimates
to be nonzero, hence performing variable selection. Here, $\lambda$ is a
hyperparameter that trades fidelity to the data—small $\lambda$ emphasizes
minimization of the squared-error—with desirable regularization that selects a
subset of variables and improves prediction accuracy.

A variant of this, the group lasso (Yuan and Lin, 2006) is appropriate when
there is a natural grouping structure for the features. That is, we assume
that both the design matrix $\mathbf{X}$ and the corresponding vector of
coefficients can be partitioned into interpretable non-overlapping groups,
and, by analogy with lasso regression, only a few of the groups are active,
i.e., have nonzero coefficients. The group lasso thus performs regularization
that has the effect of discarding groups of predictors rather than the
predictors themselves:

$\min_{\boldsymbol{\beta}}\left\\{\frac{1}{2n}\big{\lVert}\mathbf{y}-\sum_{g=1}^{G}\mathbf{X}^{(g)}\boldsymbol{\beta}^{(g)}\big{\rVert}_{2}^{2}+\lambda\sum_{g=1}^{G}\sqrt{w_{g}}\lVert\boldsymbol{\beta}^{(g)}\rVert_{2}\right\\}.$
(2)

Grouping may occur naturally—say with the inclusion of many categorical
predictors, groups of genes, or brain regions—or may be a design choice using
additive models and basis expansions. Note that in Equation 2, the grouping
structure is explicitly stated: the vector of coefficients,
$\boldsymbol{\beta}$, is thought of as a concatenation of the coefficient
subvectors of the various groups $\boldsymbol{\beta}^{(g)}$, and similarly the
data matrix $\mathbf{X}$ is the concatenation of submatrices, each submatrix
$\mathbf{X}^{(g)}$ being composed of the columns that correspond to that
particular group. Thus the first part of the equation,
$\mathbf{y}-\sum_{g=1}^{G}\mathbf{X}^{(g)}\boldsymbol{\beta}^{(g)}$, is
identical to the more simply-written equation
$\mathbf{y}-\mathbf{X}\boldsymbol{\beta}$, but the notation serves to
emphasise the partitioning.

However, the penalty,
$\sum_{g=1}^{G}\sqrt{w_{g}}\lVert\boldsymbol{\beta}^{(g)}\rVert_{2}$, is
different from the corresponding part in Equation 1, using instead the sum of
the (non-squared) $\ell_{2}$-norms of the coefficient vectors of the various
groups. It is the non-differentiability of this expression at
$\mathbf{0}\in\mathbb{R}^{|g|}$ (with $|g|$ meaning the size of group $g$)
that accounts for the group-discarding property of the solution, similar to
the way that the non-differentiability the absolute value at $0$ is
responsible for discarding individual predictors in the lasso.

As with Equation 1, there is only a single tuning parameter $\lambda$, whose
value determines the strength of regularization. Within the second summation
are the relative weights of the groups, $w_{g}$. These are often taken to be
the size of the corresponding group. For simplicity, this notation is
suppressed below where the meaning is clear.

Finally, in a group-structured problem as above, it may be desirable to
enforce sparsity not only among the groups but also within the groups. The
sparse group lasso (Simon et al., 2013) does this by combining the penalties
in Equations 1 and 2:

$\min_{\boldsymbol{\beta}}\left\\{\frac{1}{2n}\lVert\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\rVert_{2}^{2}+(1-\alpha)\lambda\sum_{g=1}^{G}\lVert\boldsymbol{\beta}^{(g)}\rVert_{2}+\alpha\lambda\sum_{g=1}^{G}\lVert\boldsymbol{\beta}^{(g)}\rVert_{1}\right\\}.$
(3)

There is now a second tuning parameter $\alpha$, which controls the relative
emphasis of intra- versus inter-group sparsity in the coefficient estimates.
In Equation 3, we have chosen to write the group dependence explicitly in the
$\ell_{1}$ component of the penalty, but note that
$\sum_{g=1}^{G}\lVert\boldsymbol{\beta}^{(g)}\rVert_{1}=\lVert\boldsymbol{\beta}\rVert_{1}$.
Similar to the weights $w_{g}$ in the group component, sparsegl also allows
individual predictor weights in the $\ell_{1}$ component,
$\sum_{j=1}^{p}\omega_{j}|\beta_{j}|$, but we suppress this generality for
clarity, setting $\omega_{j}=1$ for all $j$ below.

###  2.1 The group-wise majorization-minimization algorithm

There is no closed-form solution to the optimization problem in Equation 3, so
we require a numerical procedure. Because the problem is convex, a variety of
methods may be used. The general framework for our algorithm is the same as
the majorized block-wise coordinate descent algorithm developed in (Yang and
Zou, 2015; Simon et al., 2013). What this means is that, for a fixed value of
$\lambda$, we loop over the groups and update only those coefficients while
holding all other groups constant. In particular, instead of using the exact
Hessian to determine the step size and direction in every update step, we
update according to a simpler expression that majorizes the objective.

For the rest of this section, we describe this majorization algorithm,
focusing on a particular group $g$ and holding the coefficients for all other
groups fixed. We note here that, because the loss function in Equation 3 is
differentiable and the penalty terms are convex and separable (i.e., they can
be decomposed into a sum of functions each only involving a single group),
this block coordinate descent algorithm is guaranteed to converge to a global
optimum (Tseng, 2001).

To begin with, we introduce some notation. Let

$\mathbf{r}_{(-g)}=\mathbf{y}-\sum_{k\neq
g}\mathbf{X}^{(k)}\boldsymbol{\beta}^{(k)}$

be the partial residual without group $g$ where all the group fits besides
that of group $g$ are subtracted from $\mathbf{y}$. With all other groups held
fixed, we aim to solve:

$\min_{\boldsymbol{\beta}^{(g)}}\frac{1}{2n}\lVert\mathbf{r}_{(-g)}-\mathbf{X}^{(g)}\boldsymbol{\beta}^{(g)}\rVert_{2}^{2}+(1-\alpha)\lambda\lVert\boldsymbol{\beta}^{(g)}\rVert_{2}+\alpha\lambda\lVert\boldsymbol{\beta}^{(g)}\rVert_{1}.$
(4)

In what follows, we will suppress the $(g)$ notation, with the understanding
that we are really referring to only the $g^{\textrm{th}}$ group of the
coefficient vector and the partial residual $\mathbf{r}_{(-g)}$. We will also
define, the unpenalized loss function

$\ell(\boldsymbol{\beta})=\frac{1}{2n}\lVert\mathbf{r}-\mathbf{X}\boldsymbol{\beta}\rVert_{2}^{2},$

so that our objective function for the $g^{\textrm{th}}$ group becomes
$\ell(\boldsymbol{\beta})+(1-\alpha)\lambda\left\lVert\boldsymbol{\beta}\right\rVert_{2}+\alpha\lambda\left\lVert\boldsymbol{\beta}\right\rVert_{1}$,
and we are interested in finding an optimal value, $\hat{\boldsymbol{\beta}}$.
This enables the procedure to generalize easily to logistic loss or, in
principle, other exponential families.

Any global minimum must satisfy a subgradient equation, similar to a first-
derivative test for an optimum, except that $\left\lVert\cdot\right\rVert_{2}$
and $\left\lVert\cdot\right\rVert_{1}$ are non-differentiable at $\mathbf{0}$.
For Equation 4 above, taking the subdifferential and setting equal to zero
gives us the following first-order condition:

$\nabla\ell(\boldsymbol{\beta})=(1-\alpha)\lambda\mathbf{u}+\alpha\lambda\mathbf{v},$
(5)

where $\mathbf{u}$ is the subgradient of
$\left\lVert\boldsymbol{\beta}\right\rVert_{2}$ and $\mathbf{v}$ is the
subgradient of $\left\lVert\boldsymbol{\beta}\right\rVert_{1}$. The first is
defined to be
$\boldsymbol{\beta}/\left\lVert\boldsymbol{\beta}\right\rVert_{2}$ if
$\boldsymbol{\beta}$ is a nonzero vector, and is any vector in the set
$\\{\mathbf{u}:\left\lVert\mathbf{u}\right\rVert_{2}\leq 1\\}$ otherwise; the
second, $\mathbf{v}$, is defined coordinate-wise as
$v_{j}=\text{sign}(\beta_{j})$ if $v_{j}\neq 0$, and is any value
$v_{j}\in\\{v_{j}:|v_{j}|\leq 1\\}$ otherwise.

For the Gaussian case, the unpenalized loss $\ell(\boldsymbol{\beta})$ is a
quadratic function in $\boldsymbol{\beta}$, so it is equal to its second order
Taylor expansion about any point $\boldsymbol{\beta}_{0}$ in the parameter
space. We thus start with the following equality for any given
$\boldsymbol{\beta}_{0}$ (recalling that $\boldsymbol{\beta}_{0}$ here is only
for group $g$):

$\forall\boldsymbol{\beta},\
\boldsymbol{\beta}_{0},\quad\ell(\boldsymbol{\beta})=\ell(\boldsymbol{\beta}_{0})+(\boldsymbol{\beta}-\boldsymbol{\beta}_{0})^{\mathsf{T}}\nabla\ell(\boldsymbol{\beta}_{0})+\frac{1}{2}(\boldsymbol{\beta}-\boldsymbol{\beta}_{0})^{\mathsf{T}}\mathbf{H}(\boldsymbol{\beta}-\boldsymbol{\beta}_{0}),$

where the gradient $\nabla\ell$ is the first total derivative of $\ell$
(evaluated at $\boldsymbol{\beta}_{0}$) and $\mathbf{H}$, the Hessian, is the
second total derivative. For
$\ell(\boldsymbol{\beta})=\frac{1}{2n}\lVert\mathbf{r}-\mathbf{X}\boldsymbol{\beta}\rVert_{2}^{2}$,
a short computation shows that the Hessian is
$\mathbf{H}=\frac{1}{n}\mathbf{X}^{\mathsf{T}}\mathbf{X}$.

For the large-scale problems motivating this work, the matrix $\mathbf{X}$ is
large, so computing $\mathbf{X}^{\mathsf{T}}\mathbf{X}$, storing it in memory,
or inverting it, is computationally prohibitive. Instead, we replace this
matrix with a simpler one, $t^{-1}\mathbf{I}$, a diagonal matrix with the
value of $t$ selected to be such that this dominates the Hessian (in the sense
that $t^{-1}\mathbf{I}-\mathbf{H}$ is positive definite). For our algorithm we
choose the largest eigenvalue of the Hessian and use that for $t^{-1}$. Note
that this eigenvalue must be computed for each group $g\in G$, but this
computation is relatively simple using the power method or other techniques as
implemented with RSpectra (Qiu and Mei, 2022). This upper bound leads to the
following inequality:

$\forall\boldsymbol{\beta},\
\boldsymbol{\beta}_{0},\quad\ell(\boldsymbol{\beta})\leq\ell(\boldsymbol{\beta}_{0})+(\boldsymbol{\beta}-\boldsymbol{\beta}_{0})^{\mathsf{T}}\nabla\ell(\boldsymbol{\beta}_{0})+\frac{1}{2t}(\boldsymbol{\beta}-\boldsymbol{\beta}_{0})^{\mathsf{T}}(\boldsymbol{\beta}-\boldsymbol{\beta}_{0}).$
(6)

Replacing the loss function in the original minimization problem in Equation 4
with the right-hand side of Equation 6 leads to a majorized version of the
original problem

$\ell(\boldsymbol{\beta}_{0})+(\boldsymbol{\beta}-\boldsymbol{\beta}_{0})^{\mathsf{T}}\nabla\ell(\boldsymbol{\beta}_{0})+\frac{1}{2t}\left\lVert\boldsymbol{\beta}_{0}-\boldsymbol{\beta}\right\rVert_{2}^{2}+(1-\alpha)\lambda\left\lVert\boldsymbol{\beta}\right\rVert_{2}+\alpha\lambda\left\lVert\boldsymbol{\beta}\right\rVert_{1},$
(7)

which no longer involves operations with the Hessian matrix.

As before, the optimal value for Equation 7 is determined by its subgradient
equation, similar to that of Equation 5:

$\frac{1}{t}\Big{(}\boldsymbol{\beta}-\big{(}\boldsymbol{\beta}_{0}-t\nabla\ell(\boldsymbol{\beta}_{0})\big{)}\Big{)}+(1-\alpha)\lambda\mathbf{u}+\alpha\lambda\mathbf{v}=\mathbf{0},$

with $\mathbf{u}$ and $\mathbf{v}$ as defined above. Solving this for
$\boldsymbol{\beta}$ in terms of $\boldsymbol{\beta}_{0}$ results in the
following expression:

$\hat{\boldsymbol{\beta}}=U(\boldsymbol{\beta}_{0})=\left(1-\frac{t(1-\alpha)\lambda}{\left\lVert
S(\boldsymbol{\beta}_{0}-t\nabla\ell(\boldsymbol{\beta}_{0}),\
t\alpha\lambda)\right\rVert_{2}}\right)_{+}S\big{(}\boldsymbol{\beta}_{0}-t\nabla\ell(\boldsymbol{\beta}_{0}),\
t\alpha\lambda\big{)},$ (8)

where $(z)_{+}=\max\\{z,\ 0\\}$ and $S$ is the coordinate-wise soft threshold
operator, on a vector $\boldsymbol{\gamma}$ and scalar $b$,

$(S(\boldsymbol{\gamma},b))_{j}=\text{sign}(\gamma_{j})(|\gamma_{j}|-b)_{+},$

i.e., for each coordinate in the vector, it shrinks that coordinate in
magnitude by the amount $b$, and sets it to zero if the magnitude of that
coordinate was smaller than $b$ to begin with. It is this soft-thresholding
operation that encourages within-group sparsity.

An examination of Equation 8 shows that it is possible for the entire group to
be set to zero (made inactive) due to the (hard) threshold operator
$(\cdot)_{+}$ in the first part of the expression. It is also possible for
individual components of $\boldsymbol{\beta}^{(b)}$ to be zeroed out by the
coordinate-wise (soft) threshold operator $S$. Therefore, performing this
update step tends to enforce coefficient sparsity at both the group- and
individual-level.

Algorithm 1 Sparse group lasso solution for fixed $\lambda$, regression
version

1: Input: $\lambda\geq 0$, $\alpha\in[0,\ 1]$, set of groups $\mathcal{G}$,
initial coefficients $\boldsymbol{\beta}$,
$\mathbf{r}=\mathbf{y}-\mathbf{X}\boldsymbol{\beta}$

2: while Not converged do

3: for $g\in\mathcal{G}$ do

4: Update $\boldsymbol{\beta}^{(g)}=\left(1-\frac{t(1-\alpha)\lambda}{\lVert
S\big{(}\boldsymbol{\beta}^{(g)}-t\nabla\ell(\boldsymbol{\beta}^{(g)}),\
t\alpha\lambda\big{)}\rVert_{2}}\right)_{+}S\big{(}\boldsymbol{\beta}^{(g)}-t\nabla\ell(\boldsymbol{\beta}^{(g)}),\
t\alpha\lambda\big{)}$.

5: Update $\mathbf{r}=\mathbf{r}-\mathbf{X}^{(g)}\boldsymbol{\beta}^{(g)}$.

6: end for

7: end while

8: return $\boldsymbol{\beta}$

Above, we have focused on the Gaussian linear model with
$\ell(\boldsymbol{\beta})=\frac{1}{2n}\lVert\mathbf{r}-\mathbf{X}\boldsymbol{\beta}\rVert_{2}^{2}$,
$\nabla\ell(\boldsymbol{\beta})=-\frac{1}{n}\mathbf{X}^{\mathsf{T}}(\mathbf{r}-\mathbf{X}\boldsymbol{\beta})$,
and $\mathbf{H}\preceq t^{-1}\mathbf{I}$. In the case of logistic regression,
we use exactly the same procedure but with
$\ell(\boldsymbol{\beta})=\frac{1}{n}\sum_{i}\log(1+\exp\\{-r_{i}\mathbf{x}_{i}^{\mathsf{T}}\boldsymbol{\beta}\\})$,
$\nabla\ell(\boldsymbol{\beta})=-\frac{1}{n}\sum_{i}y_{i}\mathbf{x}_{i}^{\top}(1+\exp\\{-r_{i}\mathbf{x}_{i}^{\mathsf{T}}\boldsymbol{\beta}\\})^{-1}$,
and $\mathbf{H}(\boldsymbol{\beta})\preceq 4t^{-1}\mathbf{I}$. This procedure
is explicitly stated in Algorithm 1. For other exponential families (for
example, Poisson, Gamma, or Probit regression), we provide functionality to
pass an R family object. These will generally be much slower than the built-in
families described above because they require iteratively reweighted least
squares as an outer loop combined with inner majorization-minimization
iterations as described here.

While the procedure described so far solves Equation 3 for fixed choices of
$\lambda$ and $\alpha$, the data analyst does not typically know these ahead
of time. Rather, we would like to solve the problem for a collection of values
of $\lambda$ (and perhaps $\alpha$ as well). It turns out that the structure
of this optimization problem allows for some heuristics that can perform this
sequential optimization with a minimum of additional computational resources,
in some cases, solving Equation 3 faster for a sequence of values
$\lambda_{m}\in\\{\lambda_{1},\ldots,\lambda_{M}\\}$ than for a single choice
(Tibshirani et al., 2012). We describe our implementation of this procedure
next.

###  2.2 Sequential strong rule, KKT conditions, and active set iteration

For any fixed value of $\lambda$, many groups of coefficient estimates will
end up being equal to zero. If, somehow, we knew _which_ groups, we could
completely avoid visiting them in the blockwise coordinate descent updates,
and therefore avoid calculating Equation 8 for those groups. This would
significantly speed up computations.

Re-examining Equation 5, we can see that the first order condition implies
that, for each group $g$, any solution must satisfy

$\lVert S(\nabla\ell(\boldsymbol{\beta}_{g}),\
\lambda\alpha)\rVert_{2}\leq(1-\alpha)\lambda.$ (9)

This is because, as $\mathbf{u}$ is the subgradient of
$\lVert\boldsymbol{\beta}_{g}\rVert_{2}$, $\lVert\mathbf{u}\rVert_{2}\leq 1$.
Furthermore, if $\lVert\mathbf{u}\rVert_{2}<1$, then
$\hat{\boldsymbol{\beta}}_{j}=\mathbf{0}$. In the previous section, we used
the sufficiency of the Karush-Kuhn-Tucker (KKT) stationarity condition to
derive a solution, while Equation 9 is the necessary version. So given a
potential solution, it is easy to check its validity. Unfortunately, this is
not constructive.

The sequential strong rule (Tibshirani et al., 2012) begins from Equation 9
and makes use of the fact that we are solving for a sequence of parameters
$\\{\lambda_{1}>\lambda_{2}>\dots>\lambda_{M}\\}$ rather than a single value.
At each $\lambda_{m}$, we rely on the fact that we have already solved the
problem at $\lambda_{m-1}$ and use this information to quickly discard many
groups of predictors. Without loss of generality, for the rest of this
section, assume that the problem has been solved for $\lambda_{m-1}$.

Define $c_{g}(\lambda)=S(\nabla\ell(\boldsymbol{\beta}_{g}),\ \lambda\alpha)$.
Now, we make the assumption that $c_{g}(\lambda)$ is $(1-\alpha)$-Lipschitz,
i.e., that

$\forall\lambda,\ \lambda^{\prime}\ \ \ \lVert
c_{g}(\lambda)-c_{g}(\lambda^{\prime})\rVert_{2}\leq(1-\alpha)|\lambda-\lambda^{\prime}|.$

This Lipschitz assumption appears unintuitive, and in fact, is not always
true, but it turns out to be useful.

By Equation 9, if we knew that $\lVert
c_{g}(\lambda_{m})\rVert_{2}<(1-\alpha)\lambda_{m}$ then we could safely
ignore it. But we already have the solution at $\lambda_{m-1}$. By the
triangle inequality (first) and the Lipschitz assumption (second),

$\lVert c_{g}(\lambda_{m})\rVert_{2}\leq\lVert
c_{g}(\lambda_{m})-c_{g}(\lambda_{m-1})\rVert_{2}+\lVert
c_{g}(\lambda_{m-1})\rVert_{2}\leq(1-\alpha)(\lambda_{m-1}-\lambda_{m})+\lVert
c_{g}(\lambda_{m-1})\rVert_{2}.$

We want to be able to assert that
$(1-\alpha)(\lambda_{m-1}-\lambda_{m})+\lVert
c_{g}(\lambda_{m-1})\rVert_{2}\leq(1-\alpha)\lambda_{m}$, allowing us to
ignore group $g$, and this assertion holds precisely when

$\lVert
c_{g}(\lambda_{m-1})\rVert_{2}\leq(1-\alpha)(2\lambda_{m}-\lambda_{m-1}).$

Applying this logic to Equation 7 gives the sequential strong rule for the
sparse group lasso:

$\lVert S\big{(}\nabla\ell(\boldsymbol{\beta}_{g}),\
t\alpha\lambda_{m-1}\big{)}\rVert_{2}\leq
t(1-\alpha)(2\lambda_{m}-\lambda_{m-1}).$ (10)

For more details in related settings, see Tibshirani et al. (2012). If
Equation 10 holds, then we ignore group $g$ when solving the problem at
$\lambda_{m}$. That is to say, when we move from $\lambda_{m-1}$ to
$\lambda_{m}$, we first check this condition using the previously computed
solution for $\hat{\boldsymbol{\beta}}(\lambda_{m-1})$, and then perform
blockwise coordinate descent, using only those groups that failed this
inequality.

This discarding rule is fast, because it uses the previously computed solution
combined with a simple inequality, and, in practice, it tends to accurately
discard large numbers of groups. However, we should reiterate that it is
possible for the strong rule to fail. The Lipschitz assumption is not a
guarantee. Because of this, it is critical that, after discarding some of the
groups and running the algorithm on the others, the KKT condition is checked
on all discarded groups. If there are no violations, then we have the
solution.

To minimize gradient computations for groups that will eventually be
determined to be inactive, we actually keep track of two sets: the strong set
$\mathcal{S}$ and the active set $\mathcal{A}$. The active set collects all
groups that have ever had non-zero coefficients at previous values of
$\lambda$. We first iterate over previously active groups, then check the
strong set to see if we missed any, and finally check all the remaining
groups. When the number of groups is very large, this avoids onerous
computations for as many groups as possible. The complete algorithm including
this active set iteration is shown in Algorithm 2.

1: Input: $\mathbf{X}$, $\mathbf{y}$, $\mathcal{G}$, and
$\\{\lambda_{1},\ldots,\lambda_{M}\\}$. Output: $\hat{\boldsymbol{\beta}}$.

2: Initialize: $\mathcal{A}=\mathcal{S}=\varnothing$,
$\hat{\boldsymbol{\beta}}=0$.

3: for $m=1$ to $M$ do

4: Update
$\mathcal{S}\leftarrow\mathcal{S}\bigcup\left\\{g\in\mathcal{S}^{c}:\lVert
S(\nabla\ell(\hat{\boldsymbol{\beta}}_{g}),\
t\alpha\lambda_{m})\rVert_{2}>t(1-\alpha)\lambda_{m}\right\\}$.

5: Apply Algorithm 1 with $\mathcal{G}=\mathcal{A}$ (MM gradient update).

6: Update
$\mathcal{A}\leftarrow\mathcal{A}\bigcup\left\\{g\in\mathcal{S}\bigcap\mathcal{A}^{c}:\lVert
S(\nabla\ell(\hat{\boldsymbol{\beta}}_{g}),\
t\alpha\lambda_{m})\rVert_{2}>t(1-\alpha)\lambda_{m}\right\\}.$

7: If there were any violations, go to to Line 5.

8: Update
$\mathcal{A}\leftarrow\mathcal{A}\bigcup\left\\{g\in\mathcal{S}^{c}\bigcap\mathcal{A}^{c}:\lVert
S(\nabla\ell(\hat{\boldsymbol{\beta}}_{g}),\
t\alpha\lambda_{m})\rVert_{2}>t(1-\alpha)\lambda_{m}\right\\}.$

9: If there were any violations, go to to Line 5.

10: Set $\mathcal{S}=\mathcal{S}\bigcup\mathcal{A}$.

11: end for

Algorithm 2 Sequential strong rule and Majorization Minimization in sparsegl

###  2.3 Risk estimation

For many regularized prediction methods, tuning parameter selection is largely
performed with cross validation. However, cross validation can be
computationally expensive when the data set is large enough that the initial
fit is slow. For this reason, sparsegl provides information criteria as well
as cross validation routines.

In the Gaussian linear regression model given by Equation 2, if $\sigma$ is
unknown then a general form for a family of information criteria is given by

$\textrm{info}(C_{n},g)=\log\left(\frac{1}{n}\lVert\mathbf{y}-\mathbf{X}\hat{\boldsymbol{\beta}}\rVert_{2}^{2}\right)+C_{n}\;g(\mathrm{df}),$
(11)

where $C_{n}$ depends only on $n$, $g:[0,\infty)\rightarrow\mathbb{R}$ is a
fixed function, and the degrees of freedom ($\mathrm{df}$) measures the
complexity of the estimation procedure. The choices $C_{n}=2/n$ or
$C_{n}=\log(n)/n$ with $g(x)=x$ are commonly referred to as AIC (Akaike, 1973)
and BIC (Schwarz, 1978), respectively. Additionally, generalized cross
validation (Golub et al., 1979, GCV) is defined as

$\textrm{GCV}=\frac{\frac{1}{n}\lVert\mathbf{y}-\mathbf{X}\hat{\boldsymbol{\beta}}\rVert_{2}^{2}}{(1-\mathrm{df}/n)^{2}}.$

Written on the log scale, GCV takes the form of Equation 11 with
$g(x)=\log(1-x/n)$ and $C_{n}=-2$.

The key components for all three information criteria are the negative log
likelihood and the degrees of freedom. The first is simply a function of the
in-sample (training) error. On the other hand, the degrees of freedom, while
simple in the unregularized linear model (it is the number of parameters), is
less obvious for the sparse group lasso. In general, the degress of freedom
for any predictor $\hat{\mathbf{y}}$ of $\mathbf{y}$ is defined as (Efron,
1986)

$\mathrm{df}(\hat{\mathbf{y}})=\frac{1}{\sigma^{2}}\sum_{i=1}^{n}\mathrm{Cov}(\mathbf{y}_{i},\
\hat{\mathbf{y}}_{i}).$

In the case of any linear predictor, $\hat{\mathbf{y}}=\mathbf{A}\mathbf{y}$
for some matrix $\mathbf{A}$, it is easy to see that
$\mathrm{df}=\mathrm{tr}(\mathbf{A})$. Vaiter et al. (2012) gives an explicit
formula for the group lasso (without intragroup sparsity), but only minor
modifications are required for the sparse group lasso. We give a simplified
version of their result here.

######  Proposition 1 (Vaiter et al. 2012).

Suppose that for a fixed $\lambda>0$, the active set of $\hat{\beta}$ is
$\mathcal{A}$ and that $\mathbf{X}_{\mathcal{A}}$ is the set of columns
associated to $\mathcal{A}$. Assume that $\mathbf{X}_{\mathcal{A}}$ has full
column rank. Then,

$\mathrm{df}=\mathrm{tr}\left(\mathbf{X}_{\mathcal{A}}\left(\mathbf{X}_{\mathcal{A}}^{\mathsf{T}}\mathbf{X}_{\mathcal{A}}+\lambda\mathbf{K}\right)^{-1}\mathbf{X}_{\mathcal{A}}^{\mathsf{T}}\right).$

Here, $\mathbf{K}\in\mathbb{R}^{\mathcal{A}\times\mathcal{A}}$ is a block
diagonal matrix with each block corresponding to a group $g$ having at least 1
nonzero $\hat{\beta}$. For such a group $g$, denote
$\hat{\beta}_{g|\mathcal{A}}$ the subvector of nonzero coefficient estimates.
Then

$\mathbf{K}_{g}=\frac{1}{\lVert\hat{\boldsymbol{\beta}}_{g|\mathcal{A}}\rVert_{2}}\left(\mathbf{I}-\frac{\hat{\boldsymbol{\beta}}_{g|\mathcal{A}}\hat{\boldsymbol{\beta}}_{g|\mathcal{A}}^{\mathsf{T}}}{\lVert\hat{\boldsymbol{\beta}}_{g|\mathcal{A}}\rVert^{2}_{2}}\right).$

As long as the number of nonzero coefficients $|\mathcal{A}|$ is reasonably
small, the degrees of freedom can be efficiently calculated for each value of
$\lambda$. However, this calculation is generally cubic in $|\mathcal{A}|$. In
these cases, an approximation may be desired. We have found, in practice, that
$\lambda\mathbf{K}\approx\mathbf{0}$ is reasonably accurate, suggesting that
$\mathrm{df}\approx|\mathcal{A}|$ is also reasonable. This approximation is
exact for the lasso with $\alpha=1$ (Zou et al., 2007).

##  3 Example usage

This section provides a simple illustration of using the sparsegl package
(McDonald et al., 2022) to fit the regularization path for sparse group-lasso
penalized learning problems. We first examine the linear regression model when
the response variable is continuous and then briefly go over the logistic
regression case. The package is available from the Comprehensive R Archive
Network (CRAN) at https://CRAN.R-project.org/package=sparsegl and can be
installed and loaded in the usual manner:222The development version of the
package is hosted at https://github.com/dajmcdon/sparsegl with accompanying
documentation.

install.packages("sparsegl")

library("sparsegl")

We first create a small simulated dataset along with a vector indicating the
grouping structure.

set.seed(1010)

n <- 100

p <- 200

X <- matrix(rnorm(n*p), nrow = n, ncol = p)

beta <- c(

rep(5, 5), c(5, -5, 2, 0, 0), rep(-5, 5), c(2, -3, 8, 0, 0),

rep(0, (p - 20))

)

groups <- rep(1:(p / 5), each = 5)

eps <- rnorm(n, mean = 0, sd = 1)

y <- X %*% beta + eps

pr <- 1 / (1 + exp(-X %*% beta))

y0 <- rbinom(n, 1, pr)

The sparsegl package is mainly used with calls to 2 functions:

  * •

sparsegl(): fits sparse group regularized regression and classification
models;

  * •

cv.sparsegl(): repeatedly calls sparsegl() for the purposes of tuning
parameter selection via cross validation.

The interface is intended to closely mimic that available in other R packages
for regularized linear models, most notably glmnet (Friedman et al., 2022). To
perform the regularization path fitting at a sequence of regularization
parameters, sparsegl() takes as required inputs, only x, the design matrix,
and y, the response vector. Other optional arguments are the grouping vector
group, the family (either "gaussian" or "binomial"), a penalty vector for
group weights other than the size, the relative weight of lasso penalty,
desired lower or upper bounds for coefficient estimates, and other optional
configurations.

fit <- sparsegl(X, y, group = groups)

We include a number of S3 methods for sparsegl typical for linear models:
plot(), coef(), predict() and print(). The plot() function displays either the
coefficients or the group norms on the $y$-axis against either
$\\{\lambda_{m}\\}_{m=1}^{M}$ or the scaled penalty on the $x$-axis.333We have
chosen to implement plotting throughout the package using ggplot2 (Wickham et
al., 2022b). The resulting figures are shown in Figure 2.

plot(fit, y_axis = "coef", x_axis = "penalty", add_legend = FALSE)

plot(fit, y_axis = "group", x_axis = "lambda", add_legend = FALSE)

Figure 2: The left panel plots the estimated coefficients aganist the sparse
group penalty; while the right plots the $\ell_{2}$-norm of each group against
$\lambda$.

The coef() and predict() methods give the coefficients or predicted values for
a new design matrix $\widetilde{\mathbf{X}}$ at the requested $\lambda$’s,
potentially allowing for $\lambda$ values different from those used at the
fitting stage.

coef(fit, s = c(0.02, 0.03))[c(1, 3, 25, 29), ]

    ## 4 x 2 sparse Matrix of class "dgCMatrix"
    ##                      s1         s2
    ## (Intercept) -0.05189536 -0.1718156
    ## V2           4.71082485  4.6339817
    ## V24          .           .
    ## V28          .           .

predict(fit, newx = tail(X), s = fit$lambda[2:3])

    ##               s1        s2
    ##  [95,] -3.894658 -2.966973
    ##  [96,] -3.906349 -2.945468
    ##  [97,] -4.119689 -4.241786
    ##  [98,] -4.184564 -4.555082
    ##  [99,] -4.175593 -4.382721
    ## [100,] -4.071804 -4.091689

print(fit)

    ##
    ## Call:  sparsegl(x = X, y = y, group = groups)
    ##
    ## Summary of Lambda sequence:
    ##          lambda index nnzero active_grps
    ## Max.    0.62948     1      0           0
    ## 3rd Qu. 0.19676    26     20           4
    ## Median  0.06443    50     19           4
    ## 1st Qu. 0.02014    75     25           5
    ## Min.    0.00629   100    111          23

The cv.sparsegl() function implements $K$-fold cross validation and has a
similar signature to sparsegl(). It allows the user to choose the number of
splits and the loss function for measuring prediction or classification
accuracy on the held-out sets. Here, the S3 plot() method displays the cross-
validation curve with upper and lower confidence bounds calculated as $\pm 1$
standard error across the folds for each $\lambda$ in the regularization path
(Figure 3).

cv_fit <- cv.sparsegl(X, y, groups, nfolds = 15)

plot(cv_fit)

Figure 3: The cross validation estimate of out-of-sample prediction mean-
squared error is displayed against the sequence of $\lambda$ values.

The coef() and predict() methods work similarly to those above. The only
difference being that they can additionally accept the strings lambda.min or
lambda.1se, respectively the $\lambda$ that minimizes the average cross
validation error and the largest $\lambda$ such that the cross-validated
prediction error is within 1 standard error of the minimum.

coef(cv_fit, s = "lambda.1se")[c(1, 3, 25, 29), ]

    ## (Intercept)          V2         V24         V28
    ## 0.004435981 4.740139458 0.000000000 0.000000000

predict(cv_fit, newx = tail(X), s = "lambda.min") |> c()

    ## [1]  11.364725  39.985246   4.635314 -34.832413  -6.602096 -16.138344

For logistic regression, only a different family is required. Cross validation
can be implemented with misclassification or deviance loss (Figure 4).

fit_logit <- sparsegl(X, y0, groups, family = "binomial")

cv_fit_logit <- cv.sparsegl(

X, y0, groups, family = "binomial", pred.loss = "misclass"

)

plot(cv_fit_logit, log_axis = "none")

Figure 4: The cross validation estimate of out-of-sample error is displayed
against the sequence of $\lambda$ values. For logistic regression,
misclassification error may be used.

In some cases, when computations are at a premium, cross validation my be too
demanding for the purposes of risk estimation. For this reason, sparsegl
provides an estimate_risk() function. It can be used to compute any of AIC
(Akaike, 1973), BIC (Schwarz, 1978), and GCV (Golub et al., 1979). All three
are computed by combining the log-likelihood with a penalty term for model
flexibility. In addition to a fitted sparsegl model, estimate_risk() also
needs the original design matrix. Because the exact degrees-of-freedom can be
computationally expensive, setting approx_df = TRUE uses the number of non-
zero coefficient estimates, which can be reasonably accurate.

er <- estimate_risk(fit, X)

Figure 5: For Gaussian loss, AIC, BIC, and GCV (solid lines) along with their
minima (vertical dashed lines) can be estimated.

In this simulation, the $\lambda$ that minimizes AIC is 0.013 while the CV
minimizer is 0.01. The estimated risk curves are plotted against $\lambda$ in
Figure 5.

Additional documentation and examples are provided on the package website.

##  4 Applications

We examine two applications for which sparse group lasso is a natural
estimator. The first uses data regarding COVID-19 and trustworthiness of
information sources, which is included in the package. The second uses a very
large though sparse data set from neuroimaging. Finally, we briefly
investigate the accuracy of sparsegl relative to gglasso and CVXR.

###  4.1 Geographic distribution of trust in experts

Two typical uses for sparse group lasso are (1) additive models where
continuous predictors are expanded in a basis and (2) discrete factors as
predictors. Here we demonstrate an example using both at the same time. We
examine data from The Delphi Group at Carnegie Mellon University U.S. COVID-19
Trends and Impact Survey (CTIS), in partnership with Facebook. In particular,
we examine the publicly available contingency table reports, which break down
survey responses by age, race/ethnicity, gender, and other demographic
variables of interest. The necessary data to reproduce this analysis is
included in sparsegl as trust_experts.

In particular, we will focus on the “estimated percentage of respondents who
trust $\ldots$ to provide accurate news and information about COVID-19.” This
survey item is reported for a variety of different potential sources of
information—personal doctors/nurses, scientists, the Centers for Disease
Control (CDC), government health officials, politicians, journalists, friends,
and religious leaders. In this analysis, we average the first 4, characterize
the combination as “experts”, and use this as the response variable in a
linear model.

We regress “trust in experts” on 5 factor predictors representing month of
report, state of residence, age group, race/ethnicity, and self-reported
gender identity. We also include two continuous predictors: “estimated
percentage of people with Covid-like illness” and “estimated percentage of
people reporting Covid-like illness in their local community, including their
household” to control for the amount of exposure that respondents may have
been having to the pandemic. Both continuous predictors are incorporated with
a B-spline basis expansion and 10 degrees-of-freedom. The result is largely
similar to a generalized additive model as implemented with mgcv (Wood, 2017,
2022). The design matrix can be created using the following code:

library("dplyr")

library("splines")

data("trust_experts", package = "sparsegl")

trust_experts <- mutate(trust_experts, across(

where(is.factor),

~ `attr<-`(.x, "contrasts", contr.sum(nlevels(.x), contrasts = FALSE))

))

x <- Matrix::sparse.model.matrix(

trust_experts ~ 0 + region + age + gender + raceethnicity + period +

bs(cli, df = 10) + bs(hh_cmnty_cli, df = 10),

data = trust_experts, drop.unused.levels = TRUE)

After omitting both structural and otherwise missing data, the final model is
estimated with 9759 observations on 101 predictors. As shown in the code, we
did not use contrasts, fully expanding each factor in a one-hot encoding. This
allows all estimated coefficients to be interpreted as deviations from the
grand mean conditional on continuous predictors, which is natural. Such a
formulation (along with the group penalty) is closely related to Bayesian
linear models with separate Gaussian priors centered at 0 for each level of
the factor. Other contrasts could be used by modifying the above, but the
interpretation is more complicated. Encoded as a sparse matrix, this requires
about 2.1 MB of RAM to store, as opposed to 8.5 MB if it were dense. We
estimated the model using cv.sparsegl() and default arguments. Finally, we
chose $\lambda$ to be the largest lambda within 1 standard error of the CV
minimum (lambda.1se), resulting in a sparser model. Figure 6 displays the
estimated coefficients for the state-of-residence predictors. Even controlling
for age, race, gender, and the amount of circulating Covid-like illness, the
United States displays strong geographic disparities when it comes to
citizens’ trust in scientists and other health authorities.

Figure 6: This figure displays estimates for each state’s level of trust in
experts’ advice about Covid-19. The value displayed represents the change
relative the U.S.-wide average.

###  4.2 Estimating white matter connectivity

Figure 7: This graphic illustrates how the streamlines and voxels are
converted from a diffusion-weighted image to a linear model. Each voxel is
measured on 90 angles, so it occupies 90 rows in the data. When a streamline
(column of $\mathbf{X}$) passes through a voxel, the values within that voxel
are given by a physical model based on the direction of passage. Otherwise, if
the streamline does not cross the voxel, the respective rows are zero.

Pestilli et al. (2014) formulated an optimization model that takes as input a
set of brain connections generated using tractography algorithms and predicts
the MRI diffusion signal via a linear model (Pestilli et al., 2014; Daducci et
al., 2015). The Pestilli et al. (2014) model had no regularization, but
Aminmansour et al. (2019) extended the problem to include group-regularization
(this is an approach recently followed up by Schiavi et al., 2020). In this
study, we re-implemented the Aminmansour et al. (2019) formulation using
sparsegl to illustrate the feasibility and efficiency of the DWI modeling.

The neurological model predicts the DWI signal using the tractogram,
apportioning the image signal at each voxel to the streamlines according to
the measured gradient field. This preprocessing is shown pictorially Figure 7.
We estimate streamline weights using sparse group lasso, allowing the amount
of regularization applied to each group to be proportional to their
cardinality. For our study, we used one subject from the Human Connectome
Project (Van Essen et al., 2012). The full brain tractogram has 3M
streamlines, though we used only the streamlines identified as being part of
the Arcuate Fasciculus for illustration.444The processed data used to estimate
the sparse group lasso is available at
https://doi.org/10.6084/m9.figshare.20314917.

Aminmansour et al. (2019) used an algorithm based on Orthogonal Matching
Pursuit to estimate a related model. The data used in that study measures
diffusion in 11,823 voxels using 96 magnetic field angles and attempts to
reconstruct the image using the ENCODE method (Caiafa et al., 2017), resulting
in 1057 orientations and 868 fascicles. This results in a linear regression
problem with $n\approx 1M$ and $p=868$. In comparison, our data contains
77,630 voxels measured on 90 angles, with a tractography of 10,244
streamlines. The resulting linear model has $n\approx 6.9M$ and $p\approx
88K$, around $700\times$ the size of the previous analysis. The design matrix
would occupy over 500 GB if it were dense, but since it is only about 0.02%
non-sparse, it requires 1.8 GB of memory when stored in CSC format.

Figure 8: The group norm of the 12 groups based on neuroanatomical structure
is plotted against the magnitude of the penalty.

Estimating the group lasso using sparsegl with 12 groups and 100 values of
$\lambda$ required a little over 1 minute and about 6 GB of peak memory usage
on an Intel i7-11700K PC with 64 GB of RAM. The previous method required
nearly a day for a single value of $\lambda$. Figure 8 displays the group
norms of the 12 groups against the magnitude of the penalty.

###  4.3 Accuracy on synthetic problems

Figure 9: This figure shows the % change in the objective function for gglasso
and CVXR relative to sparsegl. The x-axis is scaled to make the $\lambda$
range comparable across conditions.

A small-scale simulated comparison illustrates that sparsegl is highly
accurate in many regimes of interest. We generate synthetic data from a linear
model and examine the objective function for the estimated model across a
range of $\lambda$. Specifically, we first generate the predictors
$\mathbf{X}$ by simulating each element $x_{ij}$ i.i.d. standard Gaussian. We
generate $n=100$ observations and $p=\\{50,\ 100,\ 150\\}$ predictors. We use
10 groups in all cases, with

$\beta=(\underbrace{1,\ldots,1,}_{\textrm{group }1}\
\underbrace{0,\dots,0}_{\textrm{group }2},\
\underbrace{1,\ldots,1}_{\textrm{group }3},\ \ldots,\
\underbrace{0,\ldots,0}_{\textrm{group }10}).$

The groups are all of size $p/10$. The expected signal is then equal to $p/2$.
We simulate the response from Equation 2 with $\sigma$ chosen to produce
expected signal-to-noise ratios (SNR) of $\\{0.1,\ 1,\ 10\\}$. We note that
this design produces group sparsity, even-numbered groups have no effect on
the response, but it does not produce within-group sparsity. Figure 9 shows
the results for the 9 combinations of these conditions across 20 values of
$\lambda$ determined automatically by sparsegl and reused for the other
packages. Note that gglasso is optimizing a different objective function than
CVXR or sparsegl which both use $\alpha=0.2$. Despite this discrepancy, the
objective values for gglasso and sparsegl are generally quite close, with
gglasso tending to be slightly higher, as expected, but by no more than 0.75%.
On the other hand CVXR is often much less accurate, especially for large
values of $\lambda$. This divergence is due to its inability to recover exact
0 solutions, tending to instead produce estimates nearly, but not exactly,
equal to 0. When $p=n=100$, CVXR is actually slightly more accurate,
especially for large or small $\lambda$. Inspecting the estimated
$\widehat{\beta}$ in this setting, it seems to have slightly less bias on the
non-zero groups, producing estimates with slightly larger magnitude.

##  5 Discussion

We developed this package for solving sparse group lasso optimization problems
using group $\ell_{2}$ and $\ell_{1}$ penalties with an eye toward
computational efficiency for very large, potentially sparse design matrices.
This efficiency is achieved through a customized Fortran implementation,
avoidance of deep copy behaviour, and the use of sequential strong rules for
the regularization parameter. We also provide heuristics for tuning parameter
selection without the need for refitting inherent in cross-validation and
enable some simple extensions such as differential weights in the $\ell_{1}$
penalty and boundary constraints on the coefficients.

## Acknowledgements

In addition to the R packages described above, this manuscript and sparsegl
benefitted from devtools (Wickham et al., 2022d), usethis (Wickham et al.,
2022a), testthat (Wickham, 2022), knitr (Xie, 2022), and rticles (Allaire et
al., 2022).

Franco Pestilli was supported with grants from the National Science Foundation
(OAC-1916518, IIS-1912270, IIS-1636893, BCS-1734853), National Institutes of
Health (NIMH R01MH126699, NIBIB R01EB030896, NIBIB R01EB029272), a Microsoft
Investigator Fellowship, and a gift from the Kavli Foundation. Daniel J.
McDonald was partially supported by the National Science Foundation
(DMS–1753171) and the National Sciences and Engineering Research Council of
Canada (RGPIN-2021-02618). The analysis in Section 4.1 is based on survey
results from Carnegie Mellon University’s Delphi Group.

## References

  * Agrawal et al. (2018) Agrawal, A., Verschueren, R., Diamond, S. and Boyd, S. (2018) A rewriting system for convex optimization problems.  Journal of Control and Decision, 5, 42–60.
  * Akaike (1973) Akaike, H. (1973) Information theory and an extension of the maximum likelihood principle.  In Proceedings of the 2nd International Symposium of Information Theory (eds. B. N. Petrov and F. Csaki), 267–281.
  * Allaire et al. (2022) Allaire, J., Xie, Y., Dervieux, C., R Foundation, Wickham, H., Journal of Statistical Software, Vaidyanathan, R., Association for Computing Machinery, Boettiger, C., Elsevier, Broman, K., Mueller, K., Quast, B., Pruim, R., Marwick, B., Wickham, C., Keyes, O., Yu, M., Emaasit, D., Onkelinx, T., Gasparini, A., Desautels, M.-A., Leutnant, D., MDPI, Taylor and Francis, Öğreden, O., Hance, D., Nüst, D., Uvesten, P., Campitelli, E., Muschelli, J., Hayes, A., Kamvar, Z. N., Ross, N., Cannoodt, R., Luguern, D., Kaplan, D. M., Kreutzer, S., Wang, S., Hesselberth, J. and Hyndman, R. (2022) rticles: Article Formats for R Markdown.  URL: https://CRAN.R-project.org/package=rticles.  R package version 0.24.
  * Aminmansour et al. (2019) Aminmansour, F., Patterson, A., Le, L., Peng, Y., Mitchell, D., Pestilli, F., Caiafa, C. F., Greiner, R. and White, M. (2019) Learning macroscopic brain connectomes via group-sparse factorization.  In Advances in Neural Information Processing Systems (eds. H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox and R. Garnett), vol. 32.  URL: https://proceedings.neurips.cc/paper/2019/hash/0bfce127947574733b19da0f30739fcd-Abstract.html.
  * Bache and Wickham (2022) Bache, S. M. and Wickham, H. (2022) magrittr: A Forward-Pipe Operator for R.  URL: https://CRAN.R-project.org/package=magrittr.  R package version 2.0.3.
  * Caiafa et al. (2017) Caiafa, C. F., Sporns, O., Saykin, A. and Pestilli, F. (2017) Unified representation of tractography and diffusion-weighted mri data using sparse multidimensional arrays.  In Advances in Neural Information Processing Systems (eds. I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan and R. Garnett), vol. 30.  URL: https://proceedings.neurips.cc/paper/2017/hash/ccbd8ca962b80445df1f7f38c57759f0-Abstract.html.
  * Candes and Tao (2007) Candes, E. J. and Tao, T. (2007) The Dantzig selector: Statistical estimation when $p$ is much larger than $n$.  The Annals of Statistics, 35, 2313–2351.
  * Civieta et al. (2021) Civieta, A. M., Aguilera-Morillo, M. C. and Lillo, R. E. (2021) Adaptive sparse group lasso in quantile regression.  Advances in Data Analysis and Classification, 15, 547–573.
  * Daducci et al. (2015) Daducci, A., Dal Palù, A., Lemkaddem, A. and Thiran, J. P. (2015) COMMIT: Convex optimization modeling for microstructure informed tractography.  IEEE Transactions on Medical Imaging, 34, 246–257.
  * Diamond and Boyd (2016) Diamond, S. and Boyd, S. (2016) CVXPY: A python-embedded modeling language for convex optimization.  Journal of Machine Learning Research, 17, 1–5.  URL: http://jmlr.org/papers/v17/15-408.html.
  * Efron (1986) Efron, B. (1986) How biased is the apparent error rate of a prediction rule?  Journal of the American Statistical Association, 81, 461–470.
  * Friedman et al. (2022) Friedman, J., Hastie, T., Tibshirani, R., Narasimhan, B., Tay, K., Simon, N. and Yang, J. (2022) glmnet: Lasso and Elastic-Net Regularized Generalized Linear Models.  URL: https://CRAN.R-project.org/package=glmnet.  R package version 4.1-4.
  * Fu et al. (2020) Fu, A., Narasimhan, B. and Boyd, S. (2020) CVXR: An R package for disciplined convex optimization.  Journal of Statistical Software, 94, 1–34.
  * Fu et al. (2022) Fu, A., Narasimhan, B., Kang, D. W., Diamond, S. and Miller, J. (2022) CVXR: Disciplined Convex Optimization.  URL: https://CRAN.R-project.org/package=CVXR.  R package version 1.0-11.
  * Gerber et al. (2017) Gerber, F., Moesinger, K. and Furrer, R. (2017) Extending R packages to support 64-bit compiled code: An illustration with spam64 and GIMMS NDVI3g data.  Computer & Geoscience, 104, 109–119.
  * Gerber et al. (2018) — (2018) dotCall64: An R package providing an efficient interface to compiled C, C++, and Fortran code supporting long vectors.  SoftwareX, 7, 217–221.
  * Golub et al. (1979) Golub, G. H., Heath, M. and Wahba, G. (1979) Generalized cross-validation as a method for choosing a good ridge parameter.  Technometrics, 21, 215–223.
  * McDonald et al. (2022) McDonald, D. J., Liang, X., Solón Heinsfeld, A. and Cohen, A. (2022) sparsegl: Sparse Group Lasso.  URL: https://CRAN.R-project.org/package=sparsegl.  R package version 1.0.1.
  * Pestilli et al. (2014) Pestilli, F., Yeatman, J. D., Rokem, A., Kay, K. N. and Wandell, B. A. (2014) Evaluation and statistical inference for human connectomes.  Nature Methods, 11, 1058–1063.
  * Qiu and Mei (2022) Qiu, Y. and Mei, J. (2022) RSpectra: Solvers for Large-Scale Eigenvalue and SVD Problems.  URL: https://CRAN.R-project.org/package=RSpectra.  R package version 0.16-1.
  * R Core Team (2022) R Core Team (2022) R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.  URL: https://www.R-project.org/.
  * Schiavi et al. (2020) Schiavi, S., Ocampo-Pineda, M., Barakovic, M., Petit, L., Descoteaux, M., Thiran, J.-P. and Daducci, A. (2020) A new method for accurate in vivo mapping of human brain connections using microstructural and anatomical information.  Science Advances, 6, eaba8245.
  * Schwarz (1978) Schwarz, G. (1978) Estimating the dimension of a model.  The Annals of Statistics, 6, 461–464.
  * Simon et al. (2013) Simon, N., Friedman, J., Hastie, T. and Tibshirani, R. (2013) A sparse-group lasso.  Journal of Computational and Graphical Statistics, 22, 231–245.
  * Simon et al. (2019) — (2019) SGL: Fit a GLM (or Cox Model) with a Combination of Lasso and Group Lasso Regularization.  URL: https://CRAN.R-project.org/package=SGL.  R package version 1.3.
  * Tibshirani (1996) Tibshirani, R. (1996) Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society B, 58, 267–288.
  * Tibshirani et al. (2012) Tibshirani, R., Bien, J., Friedman, J., Hastie, T., Simon, N., Taylor, J. and Tibshirani, R. J. (2012) Strong rules for discarding predictors in lasso-type problems.  Journal of the Royal Statistical Society B, 74, 245–266.
  * Tseng (2001) Tseng, P. (2001) Convergence of a block coordinate descent method for nondifferentiable minimization.  Journal of Optimization Theory and Applications, 109, 475–494.
  * Vaiter et al. (2012) Vaiter, S., Deledalle, C., Peyré, G., Fadili, J. and Dossal, C. (2012) The degrees of freedom of the group lasso for a general design.  Tech. rep., arXiv.
  * Van Essen et al. (2012) Van Essen, D., Ugurbil, K., Auerbach, E., Barch, D., Behrens, T., Bucholz, R., Chang, A., Chen, L., Corbetta, M., Curtiss, S., Della Penna, S., Feinberg, D., Glasser, M., Harel, N., Heath, A., Larson-Prior, L., Marcus, D., Michalareas, G., Moeller, S., Oostenveld, R., Petersen, S., Prior, F., Schlaggar, B., Smith, S., Snyder, A., Xu, J. and Yacoub, E. (2012) The Human Connectome Project: A data acquisition perspective.  NeuroImage, 62, 2222–2231.
  * Van Rossum et al. (2011) Van Rossum, G. et al. (2011) python Programming Language.  URL: https://www.python.org.
  * Wickham (2022) Wickham, H. (2022) testthat: Unit Testing for R.  URL: https://CRAN.R-project.org/package=testthat.  R package version 3.1.5.
  * Wickham et al. (2022a) Wickham, H., Bryan, J. and Barrett, M. (2022a) usethis: Automate Package and Project Setup.  URL: https://CRAN.R-project.org/package=usethis.  R package version 2.1.6.
  * Wickham et al. (2022b) Wickham, H., Chang, W., Henry, L., Pedersen, T. L., Takahashi, K., Wilke, C., Woo, K., Yutani, H. and Dunnington, D. (2022b) ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics.  URL: https://CRAN.R-project.org/package=ggplot2.  R package version 3.4.0.
  * Wickham et al. (2022c) Wickham, H., François, R., Henry, L. and Müller, K. (2022c) dplyr: A Grammar of Data Manipulation.  URL: https://CRAN.R-project.org/package=dplyr.  R package version 1.0.10.
  * Wickham et al. (2022d) Wickham, H., Hester, J., Chang, W. and Bryan, J. (2022d) devtools: Tools to Make Developing R Packages Easier.  URL: https://CRAN.R-project.org/package=devtools.  R package version 2.4.5.
  * Wood (2022) Wood, S. (2022) mgcv: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation.  URL: https://CRAN.R-project.org/package=mgcv.  R package version 1.8-41.
  * Wood (2017) Wood, S. N. (2017) Generalized Additive Models: An Introduction with R.  New York: Chapman and Hall/CRC, 2nd edn.
  * Xie (2022) Xie, Y. (2022) knitr: A General-Purpose Package for Dynamic Report Generation in R.  URL: https://CRAN.R-project.org/package=knitr.  R package version 1.40.
  * Yang and Zou (2015) Yang, Y. and Zou, H. (2015) A fast unified algorithm for solving group-lasso penalized learning problems.  Statistics and Computing, 25, 1129–1141.
  * Yang et al. (2020) Yang, Y., Zou, H. and Bhatnagar, S. (2020) gglasso: Group Lasso Penalized Learning Using a Unified BMD Algorithm.  URL: https://CRAN.R-project.org/package=gglasso.  R package version 1.5.
  * Yuan and Lin (2006) Yuan, M. and Lin, Y. (2006) Model selection and estimation in regression with grouped variables.  Journal of the Royal Statistical Society B, 68, 49–67.
  * Zeng and Breheny (2020) Zeng, Y. and Breheny, P. (2020) The biglasso package: A memory- and computation-effcient solver for lasso model fitting with big data in R.  The R Journal, 12, 6–19.
  * Zeng et al. (2022) Zeng, Y., Wang, C. and Breheny, P. (2022) biglasso: Extending Lasso Model Fitting to Big Data.  URL: https://CRAN.R-project.org/package=biglasso.  R package version 1.5.2.
  * Zou et al. (2007) Zou, H., Hastie, T. and Tibshirani, R. (2007) On the ‘Degrees of Freedom’ of the Lasso.  The Annals of Statistics, 35, 2173–2192.