# The Generalized Elastic Net for least squares regression with network-
aligned signal and correlated design

Huy<EMAIL_ADDRESS>[ Sansen
<EMAIL_ADDRESS>[ Claire<EMAIL_ADDRESS>[ Department of Statistics
The University of Chicago
Chicago, Illinois 60637
USA
E-mail: e2 E-mail: e3

###### Abstract

We propose a novel $\ell_{1}+\ell_{2}$-penalty, which we refer to as the
Generalized Elastic Net, for regression problems where the feature vectors are
indexed by vertices of a given graph and the true signal is believed to be
smooth or piecewise constant with respect to this graph. Under the assumption
of correlated Gaussian design, we derive upper bounds for the prediction and
estimation errors, which are graph-dependent and consist of a parametric rate
for the unpenalized portion of the regression vector and another term that
depends on our network alignment assumption. We also provide a coordinate
descent procedure based on the Lagrange dual objective to compute this
estimator for large-scale problems. Finally, we compare our proposed estimator
to existing regularized estimators on a number of real and synthetic datasets
and discuss its potential limitations.

62J05,

62J07,

62F12,

elastic net,

Fused Lasso,

network alignment,

Smooth Lasso,

sparsity,

total variation,

###### keywords:

[class=MSC]

###### keywords:

, and

##  1 Introduction

###  1.1 Problem formulation and the proposed penalty

Consider the usual linear regression model

$Y=X\beta^{*}+\epsilon$ (1)

where the design matrix $X\in\mathbb{R}^{n\times p}$ is random with
independent and identically distributed (i.i.d.) rows,
$\beta^{*}\in\mathbb{R}^{p}$ is the unknown true parameter, and
$\epsilon=(\epsilon_{1},\dots,\epsilon_{n})^{T}\in\mathbb{R}^{n}$ are i.i.d.
zero-mean Gaussian variables with (unknown) variance $\sigma^{2}$ and are
independent of the design matrix $X$. In addition to observing the responses
$y=(y_{1},\dots,y_{n})^{T}\in\mathbb{R}^{n}$, we also observe an undirected
simple graph $G=(V,E)$ with $p$ vertices and $m$ edges. Here, the $p$ vertices
index the entries of $\beta^{*}$ as well as the columns of $X$ (which we can
think of as feature vectors). This situation typically entails significant
correlation between feature vectors, thus leading to an ill-conditioned design
matrix. For simplicity, we assume throughout this paper that the rows of $X$
are i.i.d. $N(0,\Sigma)$-distributed. In our setting, the minimum eigenvalue
of $\Sigma\in\mathbb{R}^{p\times p}$ may be very small and $\beta^{*}$ might
be nearly unidentifiable. Although the addition of an unpenalized intercept
should present no difficulty, we assume no intercept in our model to simplify
the theoretical analysis.

We further assume that $\beta^{*}$ is structured with respect to the graph $G$
so that prediction and estimation can be done with small error, even in the
high-dimensional setting where $p\gg n$. As the entries of $\beta^{*}$ are
indexed by the vertices of $G$, a natural assumption is that $\beta^{*}_{i}$
and $\beta^{*}_{j}$ should be similar if $i$ and $j$ are adjacent vertices on
the graph $G$. This assumption is related to the notion of network cohesion as
discussed in Chapter 4 of Kolaczyk (2009): vertices may display similar
characteristics because they are connected (contagion), or they may be
connected because they have similar characteristics (homophily). Note that,
however, many prior works such as Li et al. (2019) discuss network cohesion in
the context where observations (the responses $y_{1},\dots,y_{n}$ and the rows
$x_{1},\dots,x_{n}$ of $X$) are indexed by a graph’s vertices and thus may no
longer be i.i.d., whereas we focus on the case where the features (the columns
of $X$) are indexed by the graph’s vertices. Following Li et al. (2019), we
also use the term network cohesion to cover both homophily and contagion,
without distinguishing the difference in causal direction between them.

More specifically, the notion of network cohesion encourages us to assume
either that the number of edges $(i,j)\in E$ where $\beta_{i}\neq\beta_{j}$ is
small (sparse signal jumps), or that $\beta^{*}$ is smooth enough so that
$\Gamma\beta^{*}$ lies in an $\ell_{q}$-ball, where $\Gamma$ is the edge-
incidence matrix of the graph $G$ and $0<q\leq 1$ (note that when $q\in(0,1)$,
an $\ell_{q}$-ball is not convex - see Figure 7.1 of Wainwright (2019) for an
illustration of what this ”ball” looks like). Mathematically, in our
theoretical analysis we assume either

$\|\Gamma\beta^{*}\|_{0}\leq s$ (2)

or

$\sum_{j=1}^{m}|(\Gamma\beta^{*})_{j}|^{q}\leq R_{q}$ (3)

for some $s\ll m$ or some $R_{q}>0$, respectively (see Section 1.5 for the
precise definition of any mathematical symbol). Assumption (2) means that the
number of edges with nonzero signal jumps is small and the true signal has
several piecewise constant regions on the graph, whereas Assumption (3) means
the signal is smooth over the graph in the $\ell_{q}$-sense. We use the term
network alignment to refer to either Assumption (2) or Assumption (3). In our
experiments, we sometimes also consider the notion that $\beta^{*}$ is smooth
over the graph in the sense that $\|\Gamma\beta^{*}\|_{\infty}$ is small. We
emphasize that we allow for the possibility that the entries $\beta^{*}$ are
all nonzero, as long as $\beta^{*}$ satisfies (or can be well approximated by
an oracle $\beta$ that satisfies) either Assumption (2) or Assumption (3).

Under Model (1) and either Assumption (2) or (3), we study the prediction and
estimation errors of the following estimator

$\hat{\beta}:=\arg\min_{\beta\in\mathbb{R}^{p}}\frac{1}{n}\|Y-X\beta\|_{2}^{2}+\lambda_{1}\|\Gamma\beta\|_{1}+\lambda_{2}\|\Gamma\beta\|_{2}^{2}$
(4)

which can also be rewritten as

$\hat{\beta}:=\arg\min_{\beta\in\mathbb{R}^{p}}\frac{1}{n}\|Y-X\beta\|_{2}^{2}+\lambda_{1}\sum_{(i,j)\in
E}|\beta_{i}-\beta_{j}|+\lambda_{2}\sum_{(i,j)\in E}(\beta_{i}-\beta_{j})^{2}$
(5)

Note that our focus is mainly on the penalty
$\lambda_{1}\|\Gamma\beta\|_{1}+\lambda_{2}\|\Gamma\beta\|_{2}^{2}$ where
$\Gamma$ is the incidence matrix of a general graph. Following the naming
conventions in Zou and Hastie (2005) and Tibshirani and Taylor (2011), we
refer to this penalty as the Generalized Elastic Net (GEN) penalty. The
estimator (4) can be easily extended to the generalized linear model (GLM)
setting, by replacing the term $\frac{1}{n}\|Y-X\beta\|_{2}^{2}$ with another
negative log-likelihood function from an exponential family distribution (see
Chapter 9 of Wainwright (2019) for more examples). For instance, if we have
binary responses that can be modeled with the logistic GLM, then using the
logistic log-likelihood function gives

$\hat{\beta}_{\text{logistic}}:=\frac{1}{n}\sum_{i=1}^{n}\log(1+e^{\langle
x_{i},\beta\rangle})-\bigg{\langle}\frac{1}{n}\sum_{i=1}^{n}y_{i}x_{i},\beta\bigg{\rangle}+\lambda_{1}\|\Gamma\beta\|_{1}+\lambda_{2}\|\Gamma\beta\|_{2}^{2}$
(6)

where $x_{1},\dots,x_{n}$ are the rows of $X$ and $y_{1},\dots,y_{n}$ are the
entries of $Y$ which are binary. For simplicity, we only focus on analyzing
the estimator (4) under Model (1), but analogous theoretical results for the
GLM setting should follow by adapting the theoretical framework of Chapter 6
of Bühlmann and van de Geer (2011).

###  1.2 Motivating applications

As network-linked features are quite common and we do not restrict our
attention to any particular type of graph, our proposed penalty is potentially
applicable to a wide variety of settings. We provide below a non-exhaustive
list of concrete examples where our penalty may be relevant.

Example 1: Structural MRI analysis. We consider the use of structural magnetic
resonance images (sMRI) of the brain in diagnosing Alzheimer’s disease, as in
Xin et al. (2014). In this case, the rows $x_{1},\dots,x_{n}$ of $X$ might
represent sMRI features of $n$ human subjects and the responses
$y_{1},\dots,y_{n}$ are binary variables indicating each subject’s disease
status. The estimator (6) can thus be applied using a 3D grid graph
representing contiguous brain voxels. In Xin et al. (2014), the Generalized
Fused Lasso penalty $\lambda_{L}\|\beta\|_{1}+\lambda_{1}\|\Gamma\beta\|_{1}$
is used, and this penalty leads to a solution that is both sparse and smooth.
However, it may be more reasonable to assume only that the true signal aligns
with the graph, in which case the estimator (6) may fare better for the
purpose of predicting Alzheimer’s disease.

Example 2: Microarray analysis with prior information. Following Segal et al.
(2003), we can also consider a microarray dataset with $X=[x_{ij}]$ where
$x_{ij}$ is the expression level of the $j^{\text{th}}$ gene for the
$i^{\text{th}}$ test subject, and $y_{i}$ is an outcome measure for subject
$i$ which can be continuous or discrete. Often, we have prior knowledge from
previous biomedical research in the form of gene regulatory pathways which can
serve as our graph $G$ (see Li and Li (2008) for specific examples). We can
incorporate this prior information using our GEN penalty. In Li and Li (2008),
the penalty $\lambda_{L}\|\beta\|_{1}+\lambda_{2}\beta^{T}\tilde{L}\beta$ is
used instead, where $\tilde{L}$ is the normalized Laplacian matrix. Assuming
the vast majority of genes has no effect on the outcome may make it easier to
interpret the estimated parameters (in terms of which genes may be responsible
for the outcome). However, if many of these genes can be grouped into clusters
with small (but nonzero) baseline effects on the outcome, using our penalty
may lead to better predictions.

Example 3: Microarray analysis without prior information. In the previous
example, without any prior information about gene regulatory pathways, we can
take $G$ to be the complete graph in our GEN penalty. The penalty
$\lambda_{L}\|\beta\|_{1}+\lambda_{1}\|\Gamma\beta\|_{1}$, where $\Gamma$ is
the incidence matrix of a complete graph, has been studied in She (2008) under
the name Clustered LASSO.

Example 4: Temporal data. Given a time series $\\{X_{t}\\}_{t\in\mathbb{N}}$,
we consider fitting an autoregressive model of the form
$X_{t}=\sum_{j=1}^{p}\beta_{j}X_{t-j}+\epsilon_{t}$. If the time points $t$
are sampled sufficiently far apart such that our data points
$(X_{t},X_{t-1},\dots,X_{t-p})$ can be considered independent across $t$, it
may be reasonable to apply our method with $G$ being a $p$-vertex chain graph.

###  1.3 Comparison with related works

The standalone $\ell_{1}$ penalty $\lambda_{1}\|\Gamma\beta\|_{1}$, which is
often known as the total variation penalty on graphs, has been studied
extensively in the context of the graph trend filtering problem where the
design matrix is the identity. More precisely, given the model
$Y=\beta^{*}+\epsilon$, the trend filtering estimator for $\beta^{*}$ is

$\hat{\beta}_{\text{tf}}:=\arg\min_{\beta\in\mathbb{R}^{n}}\frac{1}{n}\|Y-\beta\|_{2}^{2}+\lambda_{1}\|\Gamma\beta\|_{1}$
(7)

This estimator is also known as the analysis estimator, in the terminology of
Elad et al. (2007); see Hütter and Rigollet (2016), Wang et al. (2015),
Ortelli and van de Geer (2021) and Guntuboyina et al. (2020) for results on
prediction error bounds for the estimator (7) and its constrained form when
$\Gamma\beta^{*}$ is sparse. The graph considered in the trend filtering
problem is usually a chain or grid graph due to applications such as image
denoising, but results for other types of graphs such as trees and star graphs
are also available in the literature. The analysis matrix $\Gamma$ in (7) can
be generalized to higher order total variation operators (as defined in Wang
et al. (2015)). In comparison, we focus solely on the case where $\Gamma$ is
the incidence matrix defined in (9), and our design matrix $X$ is random with
i.i.d. rows rather than a pre-specified matrix consisting of fixed vectors
from some dictionary.

When the design matrix is general, the estimator

$\hat{\beta}_{\text{GL}}:=\arg\min_{\beta\in\mathbb{R}^{p}}\frac{1}{n}\|Y-X\beta\|_{2}^{2}+\lambda_{1}\|\Gamma\beta\|_{1}$
(8)

has been proposed by Tibshirani and Taylor (2011) (under the name Generalized
LASSO estimator) as well as Land and Friedman (1997) (where the penalty is
called variable fusion). These works mainly address computational techniques
for the estimator (8), rather than theoretical guarantees when $\beta^{*}$
aligns with the graph. The idea of working with the dual objective to derive
our algorithms comes from Kim et al. (2009) and Tibshirani and Taylor (2011).
Our analysis of the prediction and estimation errors for the estimator (4) is
also applicable to (8), and to our knowledge no similar analysis with random
design is available in the literature. However, the error bounds for our
estimator (4) are better due to the improved minimum eigenvalue in the
denominator of the bounds in Theorem 2.1.

Two previously introduced penalties which involve the Lasso penalty to induce
sparsity but are closely related to GEN have also been studied in the context
where the design matrix can be non-identity; they serve as the main benchmarks
in both our theoretical results and our experiments. The Smooth Lasso penalty
$\lambda_{L}\|\beta\|_{1}+\lambda_{2}\|\Gamma\beta\|_{2}^{2}$ was first
proposed by Hebiri and van de Geer (2011), in which the theoretical analysis
assumes fixed design and thus relies on a restricted eigenvalue assumption
(Assumption $B(\Theta)$ in Hebiri and van de Geer (2011)) on the expanded Gram
matrix $n^{-1}\tilde{X}^{T}\tilde{X}$ (see Section 1.5 for definition of
$\tilde{X}$). The Fused Lasso penalty
$\lambda_{L}\|\beta\|_{1}+\lambda_{1}\|\Gamma\beta\|_{1}$ was first proposed
by Tibshirani et al. (2005) for the chain graph. These two methods implicitly
assume that the true signal is both sparse and aligned with the graph. Such an
assumption can be overly restrictive, and sparsity of $\beta^{*}$ may not
always be a natural assumption in the general graph setting. When
$\|\beta^{*}\|_{0}=p$, error bounds proven for these estimators usually
involve the term $\frac{p\log p}{n}$. In comparison, our penalty only assumes
network alignment and should also work well in the sparse-and-smooth case when
the zero entries of $\beta^{*}$ form large contiguous blocks on the graph. The
Fused Lasso and the Smooth Lasso should only perform better than ours when
sparsity holds but the network alignment assumption is significantly violated.
Empirically, when $\beta^{*}$ aligns with the graph but is not sparse,
choosing the tuning parameters by cross-validation often results in
$\lambda_{L}$ being set to almost zero for both the Fused Lasso and the Smooth
Lasso.

In Li et al. (2018), the penalty
$\lambda_{2}\|\hat{\Gamma}\beta\|_{2}^{2}+\lambda_{1}\|\hat{\Gamma}\beta\|_{1}+\lambda_{L}\|\beta\|_{1}$
is introduced and referred to as the Graph Total Variation (GTV) method, which
involves three hyperparameters that require tuning. Unlike our penalty, the
incidence matrix $\hat{\Gamma}$ is obtained by first estimating $\Sigma$ with
$\hat{\Sigma}$ (which can depend on the design $X$ or side information) and
then treating $\hat{\Sigma}$ as the adjacency matrix of a graph $\hat{G}$ with
weighted edges. Note that this is a two-step process, and the graph $\hat{G}$
here also differs from our setting in that we do not consider non-binary edge
weights, since in many applications only a graph structure is provided.
Computationally, since we need to use 3D grid search for hyperparameter tuning
and the matrix $\hat{\Gamma}$ is very dense, the estimator introduced in Li et
al. (2018) does not scale well. Furthermore, even when we use the true
covariance $\Sigma$ to form $\hat{\Gamma}$, the performance of GTV in most of
our synthetic experiments does not compare favorably with that of our method,
Fused Lasso or Smooth Lasso. The theoretical analysis in Li et al. (2018) does
not account for the error in estimating $\Sigma$ with $\hat{\Sigma}$, which we
believe cannot be overlooked.

###  1.4 Organization of the paper

In Section 2, we provide upper bounds on the prediction and estimation errors
of our estimator (4). We specialize these bounds for the cases when $\Gamma$
is the incidence matrix of the $r$-dimensional grid graph, the star graph or
the complete graph, which are also discussed in Hütter and Rigollet (2016) and
may be encountered in practice. In Section 3, we describe our coordinate
descent algorithm in order to compute the estimator (4) for a particular
choice of $(\lambda_{1},\lambda_{2})$. We also provide runtime comparisons
between various techniques to compute (4), such as the interior point method
and the Alternating Direction Method of Multipliers (ADMM). In Section 4, we
use experiments on synthetic data to provide empirical evidence of some
properties of GEN that are suggested by our theoretical results, and compare
its performance with estimators from prior works that we discuss in Section
1.3. We also illustrate the usefulness of our proposed penalty in three real
data examples, namely COVID-19 trend prediction, Alzheimer’s disease detection
and estimation of Chicago’s crime patterns.

###  1.5 Notations and definitions

For any positive integer $n$, we denote $[n]$ as the set $\\{1,\dots,n\\}$.
For any matrix $A$, we denote by $A^{\dagger}$ the Moore-Penrose inverse of
$A$. For any vector $v$, $\|v\|_{0}$ refers to the number of nonzero entries
of $v$, and $\|v\|_{p}$ for $1\leq p\leq\infty$ refers to the usual
$\ell_{p}$-norm of $v$. We write $\mathbf{1}(\cdot)$ for the indicator
function. For a vector $v\in\mathbb{R}^{k}$ and any set $S\subseteq[k]$, we
denote by $v_{S}\in\mathbb{R}^{k}$ to be the vector with the $j^{\text{th}}$
coordinate given by $(x_{S})_{j}=x_{j}\mathbf{1}(j\in S)$. For any vector
$\theta\in\mathbb{R}^{m}$, we write $S_{\theta}$ to refer to the support
$\\{j\in[m]:\theta_{j}\neq 0\\}$ of $\theta$. We use $s$ to denote
$\|\Gamma\beta^{*}\|_{0}$. For any positive semi-definite matrix $M$, let
$\gamma_{\text{max}}(M)$ and $\gamma_{\text{min}}(M)$ denote its maximum and
minimum eigenvalues respectively, and $\text{ker}(M)$ the null space of $M$.
$I_{k}$ denotes the identity matrix of size $k$-by-$k$.

The notation $\lesssim$ means that the left-hand side (LHS) is bounded by the
right-hand side (RHS) multiplied by an absolute constant (not dependent on any
parameter of interest) that is omitted. The notation $\gtrsim$ is similarly
defined. The notation $\asymp$ means that both $\lesssim$ and $\gtrsim$ hold.
The constants $C$, $c$, $c_{1},c_{2}$ are absolute constants which are allowed
to change line by line.

Throughout this paper, the graph $G=(V,E)$ we consider is undirected and has
no self-loops. We identify the set of vertices $V$ with $[p]$ and the set of
edges $E$ with $[m]$; note that $m\leq p^{2}$, and $p\lesssim m$ if the graph
has no isolated vertices. We also denote the maximum degree of the graph $G$
by $d$ and the number of connected components of $G$ by $n_{c}$ (which is also
the dimension of the null space of $\Gamma$). The edge-vertex incidence matrix
of the graph $G$ is denoted by $\Gamma\in\\{-1,0,1\\}^{m\times p}$, which is
defined as follows: each edge $e=(i,j)\in E$ is represented by a row
$\Gamma_{e,\cdot}\in\\{-1,0,1\\}^{p}$ of $\Gamma$ whose $k^{\text{th}}$ entry
is given by

$\Gamma_{e,k}=\left\\{\begin{array}[]{rl}1&\text{if}\ k=\min(i,j)\\\
-1&\text{if}\ k=\max(i,j)\\\ 0&\text{otherwise.}\end{array}\right.$ (9)

The unnormalized Laplacian matrix of the graph $G$ (see Chung and Graham
(1997)) is then defined by $L:=\Gamma^{T}\Gamma$. We denote by
$\Pi\in\mathbb{R}^{p\times p}$ the projection matrix onto the kernel of
$\Gamma$. Note that we will use the facts $\Pi=\Pi^{T}$, $\Pi^{2}=\Pi$ and
$\Pi+\Gamma^{\dagger}\Gamma=I_{p}$ throughout the proofs.

In our theoretical analysis, we frequently make use of some definitions and
conventions from Hütter and Rigollet (2016). We denote $s_{1},\dots,s_{m}$ to
be the columns of $\Gamma^{\dagger}\in\mathbb{R}^{p\times m}$. The inverse
scaling factor of $\Gamma$ is defined as

$\rho(\Gamma):=\max_{j\in[m]}\|s_{j}\|_{2}$ (10)

while the compatibility factor of $\Gamma$ for a nonempty set $S\subseteq[m]$
is defined as

$k_{S}:=\inf_{\beta\in\mathbb{R}^{p}}\frac{\sqrt{|S|}\|\beta\|_{2}}{\|(\Gamma\beta)_{S}\|_{1}}$
(11)

Following Hebiri and van de Geer (2011), we also employ the notations

$\tilde{Y}:=\begin{pmatrix}Y\\\
0\end{pmatrix},\quad\tilde{X}:=\begin{pmatrix}X\\\
\sqrt{\lambda_{2}n}\Gamma\end{pmatrix},\quad\tilde{\epsilon}:=\begin{pmatrix}\epsilon\\\
-\sqrt{\lambda_{2}n}\Gamma\beta^{*}\end{pmatrix}$ (12)

Note that $\tilde{Y}=\tilde{X}\beta^{*}+\tilde{\epsilon}$ and we can write our
estimator as

$\hat{\beta}=\arg\min_{\beta\in\mathbb{R}^{p}}\frac{1}{n}\|\tilde{Y}-\tilde{X}\beta\|_{2}^{2}+\lambda_{1}\|\Gamma\beta\|_{1}$
(13)

##  2 Theoretical results

In this section, we aim to provide non-asymptotic bounds showing that the
estimator (4) is consistent in prediction and estimation under a network
alignment assumption, even in the high-dimensional setting where $p\gg n$. We
also show that the $\ell_{2}$ component of the penalty helps alleviate the
effects of an ill-conditioned covariance matrix $\Sigma$. Note that the tuning
parameters $\lambda_{1}$ and $\lambda_{2}$ in our theoretical analysis are
dependent on unobserved quantities $\beta^{*}$, $\Sigma$ and $\sigma$;
therefore, we cannot use the theoretical values for $\lambda_{1}$ and
$\lambda_{2}$ in practice and must in general rely on cross-validation. We do
not attempt to optimize the constants in our bounds, as our focus is on
understanding how the performance of our estimator depends on the quantities
$n$, $p$, $s$ (or $R_{q}$), $\Sigma$ and the graph $G$. All proofs are
deferred to the Appendix.

###  2.1 Main theorem

We begin by introducing bounds for the prediction and estimation errors that
are applicable to all graphs. However, these bounds may not be optimal for
some graphs, especially the $p$-vertex chain graph as in that case
$\rho(\Gamma)=\sqrt{p}$. The proof of Theorem 2.1 relies on the projection
argument used in Hütter and Rigollet (2016) to derive error bounds for the
trend filtering estimator (7). For simplicity, in the discussion of our
theoretical results, we assume that $\gamma_{\text{max}}(\Sigma)$, $n_{c}$ and
$\sigma^{2}$ are of constant order as $n$ goes to infinity. Recall that
$n_{c}$ is the dimension of $\ker(\Gamma)$, $d$ is the maximum degree of all
vertices of $G$, $L:=\Gamma^{T}\Gamma$, and $k_{S}$ is defined in (11).

######  Theorem 2.1 (Main theorem).

Fix $\delta>0$ and choose
$\lambda_{1}=32\sigma\rho(\Gamma)\sqrt{\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}}$, $\lambda_{2}\leq\frac{\lambda_{1}}{8\|\Gamma\beta^{*}\|_{\infty}}$.
Given any set $S$ satisfying both

$\frac{144\gamma_{\text{max}}(\Sigma)(\sqrt{n_{c}}+\delta)^{2}}{n}+\frac{36\lambda_{1}^{2}|S|k_{S}^{-2}}{\sigma^{2}}\leq\frac{1}{2}\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)$
(14)

and

$\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}\leq\frac{\sigma^{2}}{18}$ (15)

with probability at least $1-c_{1}\exp(-nc_{2})-\frac{2}{m}-e^{-\delta^{2}/2}$
we have

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{n_{c}+\delta^{2}}{n}+\frac{\lambda_{1}^{2}|S|k_{S}^{-2}}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}+\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}$
(16)
$\|\hat{\beta}-\beta^{*}\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}^{2}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{n_{c}+\delta^{2}}{n}+\frac{\lambda_{1}^{2}|S|k_{S}^{-2}}{\gamma_{\text{min}}^{2}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}+\frac{\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}$
(17)

Note that Theorem 2.1 is actually valid for any matrix $\Gamma$. However, if
$\Gamma$ is the incidence matrix of the graph $G$, we can further bound
$k_{S}^{-2}$ by applying Lemma 3 of Hütter and Rigollet (2016), which states
that $k_{S}^{-2}\lesssim\min(d,|S|)$.

Denote $\beta^{*}=\beta^{*}_{1}+\beta^{*}_{2}$, where
$\beta^{*}_{1}\in\text{ker}(\Gamma)$ and
$\beta^{*}_{2}\in\text{ker}(\Gamma)^{\perp}$. Note that the first term in the
RHS of (17) represents the error from estimating $\beta^{*}_{1}$, which is the
unpenalized component of $\beta^{*}$. The latter two terms represent the error
from estimating the penalized component $\beta_{2}^{*}$, and given a
particular graph $G$ we need to further bound $\rho(\Gamma)$ and $k_{S}^{-2}$
for that graph.

The estimation error bound (17) is only different from the prediction error
bound (16) by a factor of
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)$ in the
denominator. This means we have to make a stronger assumption about how fast
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)$ may decay to
zero in order to ensure the estimation error, rather than just the prediction
error, is also small. For example, when we specialize our bounds for the 3D
grid with $p$ vertices, the prediction error bound (30) only requires
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\gg\frac{s\log
p}{n}$ but the estimation error bound for this graph requires
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\gg\sqrt{\frac{s\log
p}{n}}$.

The conditions (14) and (15) on $S$ are the result of using Lemma 2.1. They
are equivalent to requiring that the RHS of (16) is sufficiently small
(smaller than $C\sigma^{2}$ for some absolute constant $C>0$). Assuming
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)$ is not too
small, it is reasonable to expect the prediction error to converge to zero as
$n$ become sufficiently large.

Theorem 1 is applicable to the estimator $\hat{\beta}_{\text{GL}}$ in (8)
(which corresponds to setting $\lambda_{2}=0$). However, when
$\gamma_{\text{min}}(\Sigma)$ is small, we may not have a meaningful error
bound for $\hat{\beta}_{\text{GL}}$. Generally, we want $\lambda_{2}$ to be as
large as possible to improve the minimum eigenvalue term without introducing
additional errors, and thus the choice
$\lambda_{2}=\frac{\lambda_{1}}{8\|\Gamma\beta^{*}\|_{\infty}}$ is
appropriate. When $\|\Gamma\beta^{*}\|_{\infty}$, the maximum signal
difference between adjacent vertices, is small (which is reasonable under the
assumption of network cohesion on $\beta^{*}$) and
$\gamma_{\text{min}}(\Sigma)$ is very close to zero, the improvement of the
minimum eigenvalue term can be significant. In contrast, in Theorem 3 of
Hebiri and van de Geer (2011) and Theorem 1 of Li et al. (2018), similar proof
ideas are used but the dependence between the $\ell_{2}$ and $\ell_{1}$ tuning
parameters is such that
$\lambda_{2}\propto\frac{\lambda_{1}}{\|L\beta^{*}\|_{\infty}}$. Since $L$ is
the second-order graph difference operator (see Wang et al. (2015) for the
definitions of higher-order total variation operators), the quantity
$\|L\beta^{*}\|_{\infty}=\|\Gamma^{T}\Gamma\beta^{*}\|_{\infty}$ is not as
related to Assumption (2) or (3) and can be much larger than
$\|\Gamma\beta^{*}\|_{\infty}$ for graphs with some high-degree nodes. For
example, for the star graph with $p$ nodes where the entries of $\beta^{*}$
are 0 at the central node and 1 at the leaves,
$\|\Gamma\beta^{*}\|_{\infty}=1$ but $\|L\beta^{*}\|_{\infty}$ is of order
$p$. The choice of $\lambda_{2}$ in Theorem 2.1, however, suggests that the
regularization effects of the $\ell_{2}$ component of the penalty may be
diminished if $\|\Gamma\beta^{*}\|_{\infty}$ is large. This is consistent with
what we observe in our synthetic experiments: when
$\|\Gamma\beta^{*}\|_{\infty}$ is large, cross-validation often yields
$\lambda_{2}\approx 0$.

The proof of Theorem 2.1 relies on the following lemma to relate the empirical
quadratic form $\frac{1}{n}\|Xv\|_{2}^{2}$ to the corresponding theoretical
quantity $\|\Sigma^{1/2}v\|_{2}^{2}$, uniformly for all $v\in\mathbb{R}^{p}$.
This lemma is an extension of the main result in Raskutti et al. (2010) for
our setting and may be of independent interest.

######  Lemma 2.1 (Restricted eigenvalue property for random Gaussian design).

If $X\in\mathbb{R}^{n\times p}$ has i.i.d. $N(0,\Sigma)$ rows and $m\geq 2$,
$n\geq 10$, then the event

$\bigg{\\{}\forall
v\in\mathbb{R}^{p}:\frac{\|Xv\|_{2}}{\sqrt{n}}\geq\frac{1}{4}\|\Sigma^{1/2}v\|_{2}-3\sqrt{\frac{\gamma_{\text{max}}(\Sigma)n_{c}}{n}}\|v\|_{2}-6\sqrt{2}\rho(\Gamma)\sqrt{\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}}\|\Gamma v\|_{1}\bigg{\\}}$

holds with probability at least $1-c_{1}\exp(-nc_{2})$, for some universal
constants $c_{1},c_{2}>0$.

By setting $S=S_{\Gamma\beta^{*}}$ and applying
$k_{S}^{-2}\lesssim\min(d,|S|)$, we obtain the following bounds which are
applicable when $\beta^{*}$ is piecewise constant on the graph $G$. When
$\rho(\Gamma)\gtrsim 1$, the second term in (18) and (19) should dominate.

######  Corollary 2.1.

If $\|\Gamma\beta^{*}\|_{0}=s$, with probability at least
$1-c_{1}\exp(-nc_{2})-\frac{2}{m}-e^{-\delta^{2}/2}$ we have

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\left(\frac{n_{c}+\delta^{2}}{n}+\rho^{2}(\Gamma)\min(d,s)\frac{s\log
p}{n}\right)$ (18)
$\|\hat{\beta}-\beta^{*}\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}^{2}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\left(\frac{n_{c}+\delta^{2}}{n}+\rho^{2}(\Gamma)\min(d,s)\frac{s\log
p}{n}\right)$ (19)

provided that the RHS of (18) is smaller than $C\sigma^{2}$.

On the other hand, if we set $S=\emptyset$, we obtain the following bounds
that are applicable when $\|\Gamma\beta^{*}\|_{1}$ is small. When $\beta^{*}$
is smoothly varying over $G$ and $\|\Gamma\beta^{*}\|_{0}$ is large, these
bounds are more helpful in explaining our estimator’s good performance.

######  Corollary 2.2.

With probability at least
$1-c_{1}\exp(-nc_{2})-\frac{2}{m}-e^{-\delta^{2}/2}$,

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{n_{c}+\delta^{2}}{n}+\sigma\rho(\Gamma)\sqrt{\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}}\|\Gamma\beta^{*}\|_{1}$ (20)
$\|\hat{\beta}-\beta^{*}\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}^{2}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{n_{c}+\delta^{2}}{n}+\frac{\sigma\rho(\Gamma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\sqrt{\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}}\|\Gamma\beta^{*}\|_{1}$ (21)

provided that the RHS of (20) is smaller than $C\sigma^{2}$.

We can also consider the notion that $\Gamma\beta^{*}$ is $\ell_{q}$-sparse
($0<q<1$), in the sense that $\sum_{j=1}^{m}|(\Gamma\beta^{*})_{j}|^{q}\leq
R_{q}$ (Assumption (3)). This notion of weak sparsity has been considered in
Raskutti et al. (2011) (where $\beta^{*}$ is assumed to lie in an
$\ell_{q}$-ball) and Cai and Zhou (2012) (where, in the context of covariance
estimation, the columns of the covariance matrix are assumed to lie in an
$\ell_{q}$-ball). In contrast, Hebiri and van de Geer (2011) defines the
smoothness of the true signal using $\ell_{2}$-norm, in the sense that
$\sum_{j=1}^{m}|(\Gamma\beta^{*})_{j}|^{2}\leq R_{2}$ for some $R_{2}>0$. If
there exists an edge with a large signal difference, $R_{2}$ can be very
large. For smaller values of $q$, we can more easily accommodate the
occasional large signal jump with a reasonably small $R_{q}$, which appears in
the bound (22).

By choosing $S$ to trade off the last two terms in the RHS of (16), we obtain
the following bound for the prediction error. The proof is routine and is thus
omitted.

######  Corollary 2.3.

With probability at least
$1-c_{1}\exp(-nc_{2})-\frac{2}{m}-e^{-\delta^{2}/2}$, if Assumption (3) holds
for some $q\in(0,1)$, we have

$\begin{split}&\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{n_{c}+\delta^{2}}{n}\\\
&+\min\left(\frac{[\sigma\rho(\Gamma)]^{2-q}\left(\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}\right)^{1-q/2}R_{q}d^{1-q}}{[\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)]^{1-q}},\frac{[\sigma\rho(\Gamma)]^{\frac{2}{1+q}}\left(\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}\right)^{\frac{1}{1+q}}R_{q}^{\frac{2}{1+q}}}{\left[\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\right]^{\frac{1-q}{1+q}}}\right)\end{split}$
(22)

provided that the RHS of (22) is smaller than $C\sigma^{2}$.

###  2.2 Discussion of the quantity
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)$

When $\Gamma=I_{p}$, our penalty is just the original Elastic Net penalty. In
that case, since $\rho(\Gamma)=1$ and $k_{S}^{-2}\leq 1$, the corresponding
estimator $\hat{\beta}_{\text{EN}}$ satisfies with high probability

$\|\Sigma^{1/2}(\hat{\beta}_{\text{EN}}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}I_{p}\right)}\frac{\|\beta^{*}\|_{0}\log
p}{n}$ (23)

Here, it is clear the minimum eigenvalue term is bounded below by
$\lambda_{2}$. When $\Gamma$ is an incidence matrix of a graph, however,
$\Gamma$ has a nontrivial kernel and so the behavior of the quantity
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)$ is less
clear.

We conjecture that under reasonable assumptions about $(\Sigma,L)$,
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)$ is bounded
below by $c\lambda_{2}$ for some absolute constant $c$, at least when
$\lambda_{2}$ is in a neighborhood of zero. We emphasize that the proof of
Theorem 2.1 makes no assumption about how $\Sigma$ is related to the graph $G$
or its Laplacian $L$. When we only have
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\geq
c\lambda_{2}$, (18) yields

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\sigma\sqrt{\gamma_{\text{max}}(\Sigma)}\|\Gamma\beta^{*}\|_{\infty}\left(\frac{n_{c}+\delta^{2}}{\rho(\Gamma)\sqrt{n\log
p}}+\rho(\Gamma)\min(d,s)s\sqrt{\frac{\log p}{n}}\right)$ (24)

but we fail to obtain any theoretical guarantee of consistency in estimation
when $\Gamma\beta^{*}$ is sparse. If we can assume
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\geq
c\sqrt{\lambda_{2}}$, however, we obtain from (19) that

$\|\hat{\beta}-\beta^{*}\|_{2}^{2}\lesssim\sigma\sqrt{\gamma_{\text{max}}(\Sigma)}\|\Gamma\beta^{*}\|_{\infty}\left(\frac{n_{c}+\delta^{2}}{\rho(\Gamma)\sqrt{n\log
p}}+\rho(\Gamma)\min(d,s)s\sqrt{\frac{\log p}{n}}\right)$ (25)

The bounds (24) and (25) may be more applicable when $\Sigma$ is ill-
conditioned and $\gamma_{\text{min}}(\Sigma)$ cannot be assumed to be bounded
away from zero. Unfortunately, characterizing the spectrum of the sum of two
symmetric matrices in terms of the spectra of the summands is known to be a
difficult problem, and we leave as an open problem the question of identifying
a reasonable assumption on $(\Sigma,L)$ (which may both have nontrivial
kernels) under which
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\geq
c\lambda_{2}$ holds. In comparison to our work, Corollary 1 of Hebiri and van
de Geer (2011) (which assumes fixed design) assumes that its restricted
eigenvalue constant $\phi_{\mu_{n}}$, defined with respect to the matrix
$\tilde{X}^{T}\tilde{X}/n$, may be greater than $\mu_{n}$ or $\sqrt{\mu_{n}}$
without further justification; here $\mu_{n}$ plays a similar role as our
$\lambda_{2}$. In order to prove
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\geq
c\lambda_{2}$, Lemma 1 of Li et al. (2018) assumes that, for some absolute
constant $c_{l}>0$, $\min_{j\in[p]}\sum_{k=1}^{p}|\Sigma_{jk}|\geq c_{l}$ and
$\max_{j\in[p]}\sum_{k=1}^{p}|\hat{\Sigma}_{jk}-\Sigma_{jk}|\leq c_{l}/4$;
again, note that $\hat{\Sigma}$ acts as the adjacency matrix of the graph
considered in Li et al. (2018). Such assumptions may be too restrictive as the
same absolute constant $c_{l}$ is used in both assumptions.

In Section 4.2, we provide empirical evidence to show that, in many situations
where the true covariance matrix $\Sigma$ reflects the structure of the graph
$G$ (that is, features indexed by adjacent or nearby nodes are more
correlated) and $\Sigma$ is degenerate, the improvement of the minimum
eigenvalue term is significant and can be better than $c\lambda_{2}$ (or even
$c\sqrt{\lambda}_{2}$).

###  2.3 Error bounds for specific types of graphs

In this section, we apply our results to some specific graph structures that
are also explored in Hütter and Rigollet (2016). Throughout this section, $s$
denotes $\|\Gamma\beta^{*}\|_{0}$ and $R_{q}$ denotes the bound on
$\sum_{j=1}^{m}|(\Gamma\beta^{*})_{j}|^{q}$. We only present prediction error
bounds here as the estimation error bounds are different only by a factor of
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)$ in the
denominator. We mainly assume
$\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}(\frac{1}{64}\Sigma+\lambda_{2}L}$
is of constant order, but we also specialize the bound (24) assuming
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\gtrsim\lambda_{2}$
for the case when $\Gamma\beta^{*}$ is sparse to illustrate the effects of the
$\ell_{2}$ component in our penalty when $\gamma_{\text{min}}(\Sigma)$ is very
small. In that situation, the bounds for the standalone $\ell_{1}$ penalty
provide no control on the errors.

The 2D grid. From Proposition 4 of Hütter and Rigollet (2016) as well as our
lower bound result on $\rho(\Gamma)$ for the 2D grid (proven in the Appendix),
we have the following lemma.

######  Lemma 2.2.

If $\Gamma$ is the incidence matrix of the 2D grid with $p$ vertices, then

$1\lesssim\rho(\Gamma)\lesssim\sqrt{\log p}$

We therefore obtain the following corollary for the 2D grid.

######  Corollary 2.4.

Let $\Gamma$ be the incidence matrix of the 2D grid with $p$ vertices. With
the same choice of $\delta$, $\lambda_{1}$ and $\lambda_{2}$ as in Theorem
2.1, with high probability we have

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\left(\frac{1+\delta^{2}}{n}+\frac{s(\log
p)^{2}}{n}\right)$ (26)
$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{1+\delta^{2}}{n}+\sigma\sqrt{\frac{\gamma_{\text{max}}(\Sigma)(\log
p)^{2}}{n}}\|\Gamma\beta^{*}\|_{1}$ (27)
$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{1+\delta^{2}}{n}+\frac{\sigma^{2-q}R_{q}\left(\frac{\gamma_{\text{max}}(\Sigma)(\log
p)^{2}}{n}\right)^{1-q/2}}{[\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)]^{1-q}}$
(28)

provided that the RHS of the bounds above are smaller than $C\sigma^{2}$. If
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\gtrsim\lambda_{2}$,
we also have

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\sigma\sqrt{\gamma_{\text{max}}(\Sigma)}\|\Gamma\beta^{*}\|_{\infty}\left(\frac{1+\delta^{2}}{\sqrt{n\log
p}}+\frac{s\log p}{\sqrt{n}}\right)$ (29)

The rates obtained in (26) and (29) are good if $s$ is of small order relative
to $n$. For example, if there is a small island of size $k$-by-$k$ where
$\beta^{*}$ attains a value distinct from its background value outside that
island (this situation can correspond to finding abnormal spots on an MRI
scan), then (26) gives us the rate $\frac{k(\log p)^{2}}{n}$, provided that
$\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}$
is of constant order. This can be compared with the rate obtained by the Lasso
estimator, which is $\frac{k^{2}\log p}{n}$ if the background value outside
the island is zero (but it fails to achieve this rate if the background value
is nonzero). However, in the situation where the 2D grid can be divided in the
middle into a left island and a right island and $\beta^{*}$ is constant on
each of these islands, then $s\asymp\sqrt{p}$ and our rates are meaningful
only in the $p\ll n$ setting.

The $r$-dimensional grid ($r\geq 3$). From Proposition 6 of Hütter and
Rigollet (2016) as well as our lower bound result on $\rho(\Gamma)$ for the
$r$-dimensional grid, we can conclude that $\rho(\Gamma)$ in this case is of
constant order, assuming $r$ is fixed.

######  Lemma 2.3.

If $\Gamma$ is the incidence matrix of the $r$-dimensional grid with $p$
vertices and $r\geq 3$, then

$c(r)\leq\rho(\Gamma)\leq C(r)$

for some constants $c(r),C(r)$ that only depend on $r$.

We obtain the following corollary for the $r$-dimensional grid.

######  Corollary 2.5.

Let $\Gamma$ be the incidence matrix of the $r$-dimensional grid with $p$
vertices, where $r\geq 3$ is fixed. With the same choice of $\delta$,
$\lambda_{1}$ and $\lambda_{2}$ as in Theorem 2.1, with high probability we
have

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\left(\frac{1+\delta^{2}}{n}+\frac{s\log
p}{n}\right)$ (30)
$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{1+\delta^{2}}{n}+\sigma\sqrt{\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}}\|\Gamma\beta^{*}\|_{1}$ (31)
$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{1+\delta^{2}}{n}+\frac{\sigma^{2-q}R_{q}\left(\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}\right)^{1-q/2}}{[\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)]^{1-q}}$
(32)

provided the RHS of the bounds above are smaller than $C\sigma^{2}$. If
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\gtrsim\lambda_{2}$,
we also have

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\sigma\sqrt{\gamma_{\text{max}}(\Sigma)}\|\Gamma\beta^{*}\|_{\infty}\left(\frac{1+\delta^{2}}{\sqrt{n\log
p}}+s\sqrt{\frac{\log p}{n}}\right)$ (33)

If we consider $r=3$ and there is a small island of size $k$-by-$k$-by-$k$
where $\beta^{*}$ attains a value distinct from its background value outside
that island, then (30) gives us the rate $\frac{k^{2}\log p}{n}$, whereas the
Lasso gives us the rate $\frac{k^{3}\log p}{n}$ if we further assume the
background value is zero. This suggests that if the signal is both sparse and
smooth over the graph, in some situations using our estimator is preferable to
using the Lasso. More generally, if the island is not cubic but rather has an
arbitrary shape, $\|\Gamma\beta^{*}\|_{0}$ should be the island’s surface
area, whereas $\|\beta^{*}\|_{0}$ should be the island’s volume.

The complete graph. As previously mentioned, we can consider regularization
with the complete graph when there is no prior structural information
available.

######  Lemma 2.4.

If $\Gamma$ is the incidence matrix of the complete graph with $p$ vertices,
$\rho(\Gamma)\asymp\frac{1}{p}$.

Proof In Proposition 10 of Hütter and Rigollet (2016), replace any ‘$\leq$’
sign with ‘$=$’. ∎

If we replace the term $\min(d,s)$ by $p$ (since $d\asymp p$), we obtain the
following corollary:

######  Corollary 2.6.

Let $\Gamma$ be the incidence matrix of the complete graph with $p$ vertices.
With the same choice of $\delta$, $\lambda_{1}$ and $\lambda_{2}$ as in
Theorem 2.1, with high probability we have

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\left(\frac{1+\delta^{2}}{n}+\frac{s\log
p}{pn}\right)$ (34)
$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{1+\delta^{2}}{n}+\frac{\sigma}{p}\sqrt{\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}}\|\Gamma\beta^{*}\|_{1}$ (35)
$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{1+\delta^{2}}{n}+\frac{\sigma^{2-q}\frac{R_{q}}{p}\left(\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}\right)^{1-q/2}}{[\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)]^{1-q}}$
(36)

provided that the RHS of the above bounds are smaller than $C\sigma^{2}$. If
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\gtrsim\lambda_{2}$,
we also have

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\sigma\sqrt{\gamma_{\text{max}}(\Sigma)}\|\Gamma\beta^{*}\|_{\infty}\left(\frac{p(1+\delta^{2})}{\sqrt{n\log
p}}+s\sqrt{\frac{\log p}{n}}\right)$ (37)

In the case when the signal takes $k\ll p$ different values, with $k-1$ of
those attained on small islands of size $l\ll p$, $s$ is of order $klp$, and
(34) yields the rate $\frac{kl\log p}{n}$, provided that
$\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}$
is of constant order. This is the same as the rate we obtain for the Lasso if
the complement of the small islands has value zero. However, if there are two
large components with two different values, $s$ is of order $p^{2}$ and so
(34) only guarantees some control when $p\ll n$. If
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)$ is of order
$\lambda_{2}$, then (37) only gives us a meaningful bound when $p\ll\sqrt{n}$,
provided that $\|\Gamma\beta^{*}\|_{\infty}$ is of constant order.

The star graph. Here, we consider the graph with $p$ nodes and with one center
node connected to $p-1$ leaves. A similar penalty has been considered by
Ollier and Viallon (2017) to model stratified data, and this penalty is useful
particularly when most outer nodes share the same value as the central node.

######  Lemma 2.5.

If $\Gamma$ is the incidence matrix of the star graph with $p$ vertices, then
$\rho(\Gamma)\asymp 1$.

Proof From Proposition 12 in Hütter and Rigollet (2016), any column $s_{j}$ of
$\Gamma^{\dagger}$ has $\|s_{j}\|_{2}^{2}=1-\frac{1}{p}$. ∎

######  Corollary 2.7.

Let $\Gamma$ be the incidence matrix of the star graph with $p$ vertices. With
the same choice of $\delta$, $\lambda_{1}$ and $\lambda_{2}$ as in Theorem
2.1, with high probability we have

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\left(\frac{1+\delta^{2}}{n}+\frac{s^{2}\log
p}{n}\right)$ (38)
$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{1+\delta^{2}}{n}+\sigma\sqrt{\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}}\|\Gamma\beta^{*}\|_{1}$ (39)
$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{1+\delta^{2}}{n}+\frac{\left(\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)R_{q}^{2}\log
p}{n}\right)^{\frac{1}{1+q}}}{[\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)]^{\frac{1-q}{1+q}}}$
(40)

provided that the RHS of the above bounds are smaller than $C\sigma^{2}$. If
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\gtrsim\lambda_{2}$,
we also have

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\sigma\sqrt{\gamma_{\text{max}}(\Sigma)}\|\Gamma\beta^{*}\|_{\infty}\left(\frac{1+\delta^{2}}{\sqrt{n\log
p}}+s^{2}\sqrt{\frac{\log p}{n}}\right)$ (41)

For the star graph, we obtain meaningful bounds for the prediction error, even
in the high-dimensional setting where $p\gg n$.

The chain graph. When $\Gamma$ is the $p$-vertex chain graph (1D grid graph),
$\rho(\Gamma)=\sqrt{p}$ and Theorem 2.1 does not yield an error bound that is
meaningful in the $p\ll n$ setting. We modify the proof of Theorem 2.1 using
an idea in Theorem 6 of Wang et al. (2015) to obtain the following bound when
$\|\Gamma\beta^{*}\|_{1}$ is small.

######  Theorem 2.2.

Let $\Gamma$ be the incidence matrix of the $p$-vertex chain graph, and fix
$\delta>0$. With an appropriate choice of $\lambda_{1}$ and
$\lambda_{2}\leq\frac{\lambda_{1}}{8\|\Gamma\beta^{*}\|_{\infty}}$, with high
probability we have

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}(\frac{1}{64}\Sigma+\lambda_{2}L)}\frac{1+\delta^{2}}{n}+\frac{(\sigma^{2}\gamma_{\text{max}}(\Sigma)\|\Gamma\beta^{*}\|_{1})^{2/3}}{\gamma_{\text{min}}^{1/3}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\sqrt[3]{\frac{p\log
p}{n^{2}}}$ (42)

provided that the bound above is smaller than $C\sigma^{2}$.

The bound above is meaningful when $n\gg\sqrt{p\log p}$ and thus sufficient to
justify the use of our estimator when $\Gamma$ is the chain graph. Optimal
error bounds under the assumption of hard sparsity on $\Gamma\beta^{*}$ are
available in the literature if $X$ is identity (see for example Ortelli and
van de Geer (2021) and Guntuboyina et al. (2020)). However, such bounds are
often derived under a “minimum length” condition, which requires that the
distances between jumps for the true signal are roughly of the same order. The
bound (42), on the other hand, requires minimal assumptions. We leave open for
future work the analysis of our estimator (4) under the assumption of hard
sparsity on $\Gamma\beta^{*}$.

##  3 Computation

In this section, we describe our coordinate descent procedure to compute the
estimator (4). For convenience, we will work with the following definition of
$\hat{\beta}$ where we replace the loss $\frac{1}{n}\|Y-X\beta\|_{2}^{2}$ in
(4) by $\frac{1}{2}\|Y-X\beta\|_{2}^{2}$. Note that this simply corresponds to
a different scaling of $\lambda_{1}$ and $\lambda_{2}$.

$\hat{\beta}:=\arg\min_{\beta\in\mathbb{R}^{p}}\frac{1}{2}\|Y-X\beta\|_{2}^{2}+\lambda_{1}\|\Gamma\beta\|_{1}+\lambda_{2}\|\Gamma\beta\|_{2}^{2}$
(43)

Again, let $\tilde{Y}:=\begin{pmatrix}Y\\\ 0\end{pmatrix}\in\mathbb{R}^{n+m}$,
$\tilde{X}:=\begin{pmatrix}X\\\
\sqrt{2\lambda_{2}}\Gamma\end{pmatrix}\in\mathbb{R}^{(n+m)\times p}$ so that
we can write

$\hat{\beta}:=\arg\min_{\beta\in\mathbb{R}^{p}}\frac{1}{2}\|\tilde{Y}-\tilde{X}\beta\|_{2}^{2}+\lambda_{1}\|\Gamma\beta\|_{1}$
(44)

If we fix $\lambda_{2}$, the solution path in terms of $\lambda_{1}$ for (44)
as well as its dual objective is piecewise linear, and a path-finding
algorithm for the dual objective yielding the entire solution path in terms of
$\lambda_{1}$ has been proposed in Tibshirani and Taylor (2011). However, for
the purpose of selecting tuning parameters, this is of limited usefulness
since $\lambda_{2}$ needs to be fixed. The solution path in terms of
$\lambda_{2}$ is not piecewise linear, and so we cannot use a LARS-like
algorithm to get the entire path in terms of both $(\lambda_{1},\lambda_{2})$.
Also, as mentioned in Tibshirani and Taylor (2011), the set of knots in the
solution path becomes very large as the problem size increases, and at each
knot we must solve a large least squares problem (especially at the
regularized end of the path, which is typically the region of interest in this
paper) in order to compute the whole path. If we want to compute (44) for a
small set of candidate $(\lambda_{1},\lambda_{2})$ values, then the path
algorithm in Tibshirani and Taylor (2011) is unlikely to be the most
efficient.

###  3.1 Coordinate descent on the dual objective

Our coordinate descent algorithm builds upon the dual problem derived in
Tibshirani and Taylor (2011) for the Generalized Lasso. Tibshirani and Taylor
(2011) suggests, without explicit derivations, that we can use coordinate
descent on the dual problem to compute the solution of (44) for a fixed value
of $(\lambda_{1},\lambda_{2})$. Coordinate descent cannot be directly applied
to the primal objective (44) as the $\ell_{1}$-penalty here is not separable
in terms of $\beta$; in such a situation, coordinate descent does not
necessarily converge. However, the dual objective (47) has a non-smooth
component that is separable, and thus convergence is guaranteed (since
conditions (A1), (B1)-(B3) and (C2) from Tseng (2001) hold). For completeness,
we fully derive this coordinate descent algorithm on the dual and provide
experiments to convince the reader that our estimator can be efficiently
computed.

Define $\check{Y}:=\tilde{X}\tilde{X}^{\dagger}\tilde{Y}\in\mathbb{R}^{m+n}$,
$\check{\Gamma}:=\Gamma\tilde{X}^{\dagger}\in\mathbb{R}^{m\times(m+n)}$. From
Equation (36) of Tibshirani and Taylor (2011), the dual problem is:

$\hat{u}=\arg\min_{u\in\mathbb{R}^{m}}\frac{1}{2}\|\check{Y}-\check{\Gamma}^{T}u\|_{2}^{2}\quad\text{subject
to }\|u\|_{\infty}\leq\lambda_{1},\Gamma^{T}u\in\text{row}(\tilde{X})$ (45)

and the primal-dual relation, as in Equation (37) of Tibshirani and Taylor
(2011), is:

$\hat{\beta}=\tilde{X}^{\dagger}(\check{Y}-\check{\Gamma}^{T}\hat{u})+z$ (46)

where $z\in\text{ker}(\tilde{X})$. In most situations, the augmented matrix
$\tilde{X}:=\begin{pmatrix}X\\\ \sqrt{2\lambda_{2}}\Gamma\end{pmatrix}$ has a
trivial kernel, in which case $\text{row}(\tilde{X})=\mathbb{R}^{p}$ and we
can ignore $z$ as well as the constraint
$\Gamma^{T}u\in\text{row}(\tilde{X})$. Now if we let
$Q:=\check{\Gamma}\check{\Gamma}^{T}\in\mathbb{R}^{m\times m}$ and
$b:=\check{\Gamma}\check{Y}\in\mathbb{R}^{m}$, then we can write the dual
objective as:

$\hat{u}=\arg\min_{u\in\mathbb{R}^{m}}\frac{1}{2}u^{T}Qu-b^{T}u\quad\text{subject
to }\|u\|_{\infty}\leq\lambda_{1}$ (47)

We denote the projection map from $\mathbb{R}$ onto $[-\lambda,\lambda]$ by
$T_{\lambda}(\cdot)$:

$T_{\lambda}(x):=\left\\{\begin{array}[]{rl}\lambda&\text{if}\ x>\lambda\\\
x&\text{if}\ -\lambda\leq x\leq\lambda\\\ -\lambda&\text{if}\
x<-\lambda\end{array}\right.$ (48)

Our coordinate descent algorithm is presented below.

Input: $\lambda_{1},\lambda_{2},\Gamma,Y,X$, tolerance $\epsilon$

Output: $\hat{\beta}$ as defined in (43)

1 Compute
$Q=\check{\Gamma}\check{\Gamma}^{T}=(\Gamma\tilde{X}^{\dagger})(\Gamma\tilde{X}^{\dagger})^{T}$
and $b=\check{\Gamma}\check{Y}=\Gamma\tilde{X}^{\dagger}\tilde{Y}$

2 Initialize $\hat{u}_{i}^{(0)}\leftarrow 0$ for all $i\in[m]$
$\triangleright$ 0 is an arbitrary choice in $[-\lambda_{1},\lambda_{1}]$

3 while  _$\|\hat{u}^{(k)}-\hat{u}^{(k-1)}\|_{2} >\epsilon$ _ do

4 $\hat{u}_{i}^{(k+1)}\leftarrow
T_{\lambda_{1}}\left(\frac{b_{i}-\sum_{j<i}Q_{ij}\hat{u}_{j}^{(k+1)}-\sum_{j>i}Q_{ij}\hat{u}_{j}^{(k)}}{Q_{ii}}\right)$

5Compute
$\hat{\beta}\leftarrow\tilde{X}^{\dagger}(\check{Y}-\check{\Gamma}^{T}\hat{u})$

6 Return $\hat{\beta}$

Algorithm 1 Coordinate descent on the dual objective

For general GLM loss functions, we can also derive the dual problem with a
separable non-smooth constraint; however, we may not be able to write the
coordinate descent updates in closed form (we can only do so in Algorithm 1
because the dual objective (47) is quadratic). In this case, we can use
coordinate proximal gradient descent, in which we apply the projection
operator to the gradient descent update for each coordinate.

In the Appendix, we also provide an alternative algorithm to compute (43),
based on the interior point method applied to the dual objective (as in Kim et
al. (2009)). This algorithm will be denoted as IP in the following section.

###  3.2 Runtime comparisons

We compare the runtimes for computing the estimator (4) using Algorithm 1
(CD), IP, ADMM, and the Embedded Conic Solver (ECOS) from Domahidi et al.
(2013) applied to the primal objective. ECOS is a generic solver for second-
order cone programs (SOCP) that performs well for small or medium-sized
problems. We use the highly optimized ECOS implementation in the Python
package CVXPY to serve as a benchmark for comparing the runtimes of our
algorithms. Figure 1 shows the growth of empirical runtimes as $n$ or $p$
increases for signals over the chain graph (where $m=p-1$) with
$\|\Gamma\beta^{*}\|_{\infty}=0.3$ fixed; here, the hyperparameters
$\lambda_{1},\lambda_{2}$ are chosen according to our theory so as to satisfy
$\lambda_{2}=\frac{\lambda_{1}}{8\|\Gamma\beta^{*}\|_{\infty}}$.

As we can see from Figure 1, our coordinate descent algorithm scales well as
$n$ and $p$ increase, and its runtime does not exceed 10 seconds if $n$ and
$p$ are both smaller than 1,000. More generally, when $\lambda_{2}$ is not too
close to zero, the matrix $Q=\check{\Gamma}\check{\Gamma}^{T}$ is not ill-
conditioned and our coordinate descent algorithm performs quite well. We note
that this is the setting where our estimator (4) should be preferred over the
Generalized Lasso estimator $\hat{\beta}_{\text{GL}}$ in (8), whose accuracy
is impeded by the ill-conditioned nature of the matrix $Q$ when $\lambda_{2}$
is equal to zero. While our estimator requires a two-dimensional grid search
to choose $(\lambda_{1},\lambda_{2})$, Algorithm 1 can significantly reduce
the time it takes to perform hyperparameter tuning, even for large-scale
problems where $p$ and $n$ are both in the thousands. Note that when both $n$
and $p$ are not too large, the generic SOCP solver ECOS can also be
competitive.

As for our interior point method, the main computational bottleneck is the
cost of solving a linear equation involving the Hessian matrix; in other
words, we need to solve the problem $Ax=b$ for each iteration, where $A$ is an
$m$-by-$m$ matrix. Solving it requires $O(m^{3})$ operations, and thus IP can
do well only if the number of iterations required is small. Figure 1(b) shows
that in the case of the chain graph, when we fix $n$ and increase $p$, IP
still performs better than the generic solver ECOS and scales well with $p$.

(a) $p=2000$

(b) $n=2000$

Figure 1: Runtimes of different algorithms (reported on the log scale) when
(a) $p$ is fixed but $n$ increases, or (b) $n$ is fixed but $p$ increases. The
tolerance levels for IP, CD, and ECOS are set at $10^{-4}$. The tolerance
level for ADMM is $10^{-3}$. Signals are defined on a 1D chain graph with $p$
vertices. In both situations, CD has the best runtime scaling, and IP scales
better than ECOS.

We also examine the runtimes for the 2D grid as well as the star graph when
$n$ is fixed but $p$ increases. For these graphs, IP no longer scales well
with $p$ whereas CD still has the best scaling, and ECOS is also competitive
for small problem sizes.

(a) $n=2000$, 2D grid

(b) $n=2000$, star graph

Figure 2: Runtimes of different algorithms (reported on the log scale) when
$n$ is fixed but $p$ increases. (a) Signals are defined on a $p$-vertex 2D
grid graph ($m=2p-2\sqrt{p}$) with $\|\Gamma\beta^{*}\|_{\infty}=0.66$. (b)
Signals are defined on a $p$-vertex star graph ($m=p-1$) with
$\|\Gamma\beta^{*}\|_{\infty}=0.5$. The tolerance levels for IP, CD, and ECOS
are set at $10^{-4}$. The tolerance level for ADMM is $10^{-3}$. As before,
$(\lambda_{1},\lambda_{2})$ are chosen according to theory. In both
situations, CD has the best runtime scaling.

##  4 Experiments

In this section, we present the empirical performance of our penalty
$\lambda_{1}\|\Gamma\beta\|_{1}+\lambda_{2}\|\Gamma\beta\|_{2}^{2}$ and
compare with some existing penalized M-estimators in the literature, under
several synthetic settings where we vary the true signal structure and graph
topology. Particularly, we focus on the case where $\beta^{*}$ is not sparse
but aligns with the graph $G$. The design matrix is also allowed to be
correlated in a way such that two vertices have more correlated feature
vectors if they are adjacent or nearby on the graph $G$. Such a covariance
structure is natural for node-indexed feature vectors and is in line with the
notion of network cohesion discussed in Section 1. We list the methods to
which we compare our estimator below.

Graph-independent methods that do not take into account the graph provided.
These methods usually do not perform well in the setting we describe above,
and they mainly serve as benchmarks for comparison.

  1. 1.

The ordinary least squares (OLS) estimator, which is a standard method in the
setting when $p<n$ and the underlying signal is dense. It often does not
perform well when we are in the high-dimensional setting ($p>n$) or the design
is highly correlated and $\gamma_{\text{min}}(\Sigma)$ is close to zero.

  2. 2.

The Lasso (L) penalty $\lambda_{L}\|\beta\|_{1}$ from Tibshirani (1996), which
can perform well in the $p\gg n$ setting if the true signal is known to be
sparse. In the $p>n$ case, however, it has been shown to select at most $n$
variables before it saturates. As discussed in Zou and Hastie (2005), the
Lasso lacks the ability to select groups of correlated variables, and it is
empirically observed to suffer from unstable selections in the presence of
high correlation between features.

  3. 3.

The Elastic Net (EN) penalty
$\lambda_{L}\|\beta\|_{1}+\lambda_{E}\|\beta\|_{2}^{2}$, which was developed
in Zou and Hastie (2005) to deal with highly correlated predictors. The
Elastic Net tends to encourage strongly correlated predictors to be in or out
of the model together while also preserving sparsity of representation like
the Lasso. It is a suitable candidate in our setting due to our assumption of
highly correlated design.

Graph-based methods that utilizes information from the given graph $G$ (except
for possibly the GTV method). We have described these methods in Section 1.3.

  1. 4.

The Fused Lasso (FL) penalty
$\lambda_{1}\|\Gamma\beta\|_{1}+\lambda_{L}\|\beta\|_{1}$ proposed in
Tibshirani et al. (2005) encourages the resulting estimate to be both sparse
and piecewise constant with respect to $G$. This penalty may be suitable if we
believe the true signal is sparse and also forms clusters on $G$ (that is, in
each cluster the true signal attains a single value). When the true signal is
not sparse, the tuning parameter $\lambda_{L}$ is often set to zero if we use
cross-validation (CV) for hyperparameter selection, and FL degenerates into
our GEN penalty with $\lambda_{2}=0$.

  2. 5.

The Smooth Lasso (SL) penalty
$\lambda_{2}\|\Gamma\beta\|_{2}^{2}+\lambda_{L}\|\beta\|_{1}$ in Hebiri and
van de Geer (2011) results in an estimate that is smooth, in the sense that
$\|\Gamma\hat{\beta}_{\text{SL}}\|_{\infty}$ is small. It is useful when
$\beta^{*}$ is sparse and we also believe $\|\Gamma\beta^{*}\|_{\infty}$ is
small. When the true signal is not sparse and we use CV for hyperparameter
selection, $\lambda_{L}$ for SL is often set to zero, in which case SL also
degenerates into our GEN penalty with $\lambda_{1}=0$.

  3. 6.

The Graph Total Variation (GTV) penalty
$\lambda_{1}\|\hat{\Gamma}\beta\|_{1}+\lambda_{2}\|\hat{\Gamma}\beta\|_{2}^{2}+\lambda_{L}\|\beta\|_{1}$
in Li et al. (2018) estimates $\Sigma$ with some covariance estimator
$\hat{\Sigma}$ and then treats $\hat{\Sigma}$ as the weighted adjacency matrix
of some graph $\hat{G}$ with incidence matrix $\hat{\Gamma}$. In our
experiments, as suggested by Li et al. (2018), the estimator $\hat{\Sigma}$ is
obtained by hard-thresholding the sample covariance matrix (see Bickel and
Levina (2008) for details). This choice of $\hat{\Sigma}$ means that we also
need to tune the covariance threshold in addition to the 3 hyperparameters
that appear in the GTV penalty. In general, however, $\hat{\Sigma}$ can be any
covariance estimator and can also incorporate side information such as the
graph $G$ provided in our setting.

  4. 7.

We also denote by GTV-oracle the GTV penalty based on using the unobserved
covariance matrix $\Sigma$ (rather than $\hat{\Sigma}$) to construct the
corresponding incidence matrix $\hat{\Gamma}_{\text{oracle}}$. Using the true
covariance matrix should eliminate any error from covariance estimation.
However, if all entries of $\Sigma$ are nonzero, computation of the GTV-oracle
estimator can be especially challenging, since the graph used in the GTV
penalty here is a weighted complete graph.

###  4.1 Experiments on synthetic data

We repeatedly generate training and testing data from the model
$y=X\beta^{*}+\epsilon$, where the rows of $X$ are generated i.i.d. from
$N(0,\Sigma)$ and independent of $\epsilon$ which is generated from
$N(0,\sigma^{2}I_{n})$. Hyperparameter selection via CV is performed using a
separate validation data set. We report the estimation error
$\|\hat{\beta}-\beta^{*}\|_{2}$, as well as the prediction error
$\frac{1}{n}\|X_{\text{test}}(\hat{\beta}-\beta^{*})\|_{2}^{2}$ computed using
the testing data.

####  4.1.1 Choices of $\Sigma$ and the graph $G$

We consider the chain graph, the 2D grid graph and the barbell graph in our
experiments. The first two graphs allow for easier visualization of the true
and estimated signals defined on them. The barbell graph, which consists of
two non-overlapping cliques connected by a single path that has an endpoint in
each clique, allows us to test the performance of our method on a denser graph
with a less homogenous degree distribution.

As previously mentioned, $\Sigma$ is constructed so that nearby nodes have
more correlated feature vectors. For the chain graph, we use the Toeplitz
covariance structure with $\Sigma_{ij}=\rho^{|i-j|}$ where, if not stated
otherwise, we typically choose $\rho=0.5$ (moderate correlation). For the 2D
grid and barbell graph, we construct $\Sigma$ by inverting the matrix
$L+0.5I_{p}$ (recall that $L$ denotes the Laplacian of the graph $G$) and then
normalize $\Sigma$ so that all covariates have unit variance. The resulting
covariance matrices obtained from this process are illustrated in Figure 3.

Figure 3: Left: the covariance matrix obtained for a 2D grid graph with
$p=3\times 3$ vertices. Right: the covariance matrix obtained for a barbell
graph with two cliques $\\{1,2,3\\}$ and $\\{7,8,9\\}$ connected by the path
$\\{3,4,5,6,7\\}$. Note that correlation is higher for adjacent or nearby
vertices.

####  4.1.2 Hyperparameter selection and tuning time

We select hyperparameters based on 5-fold CV using a fine grid search, where
each hyperparameter is chosen from a list of at least 20 values. The scorer
for CV is the negative mean squared error (MSE)
$-\frac{1}{n}\|Y-X\hat{\beta}\|_{2}^{2}$, which tends to select for
hyperparameters with better prediction performance.

Hyperparameter tuning for GEN is computationally manageable. When $G$ is a
chain graph, the tuning time for GEN is comparable to that of other methods
with two hyperparameters, namely the Elastic Net, the Fused Lasso and the
Smooth Lasso. If we disregard the covariance thresholding parameter, the GTV
penalty still involves three hyperparameters, and the graph $\hat{\Gamma}$
used in the GTV penalty, computed using the covariance estimate
$\hat{\Sigma}$, has more nonzero weighted edges compared to the given graph
$\Gamma$. These factors contribute to longer tuning time for the GTV penalty.
The tuning time is much worse for the GTV-oracle penalty since the true
covariance matrix $\Sigma$ is denser than $\hat{\Sigma}$, and so
$\hat{\Gamma}_{\text{oracle}}$ has many more nonzero weighted edges than
$\hat{\Gamma}$ does. We present the tuning times for a toy example where each
hyperparameter is selected from a small grid search in Table 1. Here, since
$n$ and $p$ are both not too large, we use the SOCP solver ECOS (see Section
3.2) for all methods.

Table 1: Tuning times with ECOS when $G$ is the chain graph. $p=110,m=109,n=210,\sigma=1$, and $\Sigma$ is constructed as in Section 4.1.1. The GTV penalty is based on $\hat{\Gamma}$ which has around 200 nonzero weighted edges. The GTV-oracle penalty is based on $\hat{\Gamma}_{\text{oracle}}$ which has almost 6000 nonzero weighted edges. 5-fold CV is performed for each method on a small grid with 5 candidate values $[0,0.1,1,10,100]$ for each hyperparameter. |  L |  EN |  FL |  SL |  GTV |  GTV-oracle |  GEN
---|---|---|---|---|---|---|---
# hyperparameters |  1 |  2 |  2 |  2 |  3 |  3 |  2
time [seconds] |  0.45 |  3.22 |  2.85 |  2.30 |  14.31 |  134.08 |  2.46

When the graph $G$ contains more edges, we can expect the tuning time for GEN
to increase relative to other two-hypterparameter methods, as both the
$\ell_{1}$ and $\ell_{2}$ components of the GEN penalty depend on $\Gamma$.
Table 2 repeats the above experiment but with $G$ being the barbell graph and
$\Sigma$ reflecting the structure of this graph. The tuning time with ECOS for
GEN is roughly double that of FL or SL, whose penalties contain only one
component depending on $\Gamma$.

Table 2: Tuning times with ECOS when $G$ is the barbell graph. $p=110,m=2461,n=210,\sigma=1$, and $\Sigma$ is constructed as in Section 4.1.1. The GTV penalty is based on $\hat{\Gamma}$ which has around 2500 nonzero weighted edges. As $\Sigma$ for the barbell graph is denser than $\Sigma$ in Table 1, $\hat{\Sigma}$ here is also denser than $\hat{\Sigma}$ in Table 1. |  L |  EN |  FL |  SL |  GTV |  GTV-oracle |  GEN
---|---|---|---|---|---|---|---
# hyperparameters |  1 |  2 |  2 |  2 |  3 |  3 |  2
time [seconds] |  0.44 |  2.51 |  4.62 |  4.57 |  52.47 |  131.91 |  9.68

####  4.1.3 Comparisons between GEN, FL and SL when $\beta^{*}$ is dense but
aligns with $G$

In this section, we focus on the case when $\beta^{*}$ is not sparse but
$\Gamma\beta^{*}$ is sparse or $\beta^{*}$ is otherwise smooth with respect to
$G$. As all parts of the true signals constructed in this section are far from
zero, the component $\lambda_{L}\|\beta\|_{1}$ in the FL and SL penalties is
of little use, and setting $\lambda_{L}>0$ worsens both prediction and
estimation errors in this setting. Consequently, CV yields $\lambda_{L}$
values that are almost identically zero for both FL and SL. Essentially, in
this section, FL refers to the standalone $\lambda_{1}\|\Gamma\beta\|_{1}$
penalty, whereas SL refers to the standalone
$\lambda_{2}\|\Gamma\beta\|_{2}^{2}$ penalty.

We observe that FL performs well when $\beta^{*}$ has few signal jumps on $G$,
regardless of whether there exists large jumps (i.e.
$\|\Gamma\beta^{*}\|_{\infty}$ is large). SL, on the other hand, tends to
perform well when the signal is smooth with respect to $G$, in the sense that
$\|\Gamma\beta^{*}\|_{\infty}$ is small, even if the number of signal jumps
$\|\Gamma\beta^{*}\|_{0}$ might be large. To demonstrate these observations,
we construct signals with varying smoothness ($\|\Gamma\beta^{*}\|_{\infty}$)
and numbers of jumps ($\|\Gamma\beta^{*}\|_{0}$). Figure 4 illustrates the
true signals on the 1D chain graph, whereas Figure 5 illustrates the true
signals on the 2D grid graph (note that $p$ is fixed for these graphs). For
the barbell graph, we let the signal values be constant (at 5 and 20
respectively) on each clique. The lengths of the path connecting the two
cliques are chosen from $\\{1,4,7,10,13,16\\}$ (and so $p$ has to vary), and
we let the signal decrease from 20 to 5 gradually on the connecting path, so
that $\|\Gamma\beta^{*}\|_{\infty}$ decreases from 15 to 1.46 while
$\|\Gamma\beta^{*}\|_{0}$ increases from 1 to 16.

Figure 4: True signals defined on the chain graph with $p=110$. The top left
signal is piecewise constant and has the smallest $\|\Gamma\beta^{*}\|_{0}=3$
but the largest $\|\Gamma\beta^{*}\|_{\infty}=5$. The bottom right signal is
the smoothest with the largest $\|\Gamma\beta^{*}\|_{0}=99$ and the smallest
$\|\Gamma\beta^{*}\|_{\infty}=0.24$. The intermediate signals are constructed
such that $\|\Gamma\beta^{*}\|_{0}$ decreases but
$\|\Gamma\beta^{*}\|_{\infty}$ increases gradually. All 6 signals have
$\|\Gamma\beta^{*}\|_{1}=15$.

Figure 5: True signals defined on the 2D grid with $p=15\times 15$. The top
left signal is piecewise constant and has the smallest
$\|\Gamma\beta^{*}\|_{0}=28$ but the largest $\|\Gamma\beta^{*}\|_{\infty}=3$.
The bottom right signal is the smoothest with the largest
$\|\Gamma\beta^{*}\|_{0}=412$ and the smallest
$\|\Gamma\beta^{*}\|_{\infty}=0.24$. All 6 signals have
$\|\Gamma\beta^{*}\|_{1}$ between 84 and 120.

Figure 6 illustrates the performances of FL, SL and GEN in terms of estimation
and prediction errors. When $\|\Gamma\beta^{*}\|_{\infty}$ is small, CV yields
$\lambda_{2}$ values that are larger relative to $\lambda_{1}$, which is
consistent with our theory in that $\lambda_{2}$ can be chosen up to
$\frac{C\lambda_{1}}{\|\Gamma\beta^{*}\|_{\infty}}$ without incurring
additional errors. As can be seen from Figure 6, our GEN penalty adapts well
to true signals of various smoothness levels, thus demonstrating the
importance of having both the $\lambda_{1}\|\Gamma\beta\|_{1}$ and
$\lambda_{2}\|\Gamma\beta\|_{2}^{2}$ components in our penalty. From the
performances of FL and SL, we can see that the
$\lambda_{1}\|\Gamma\beta\|_{1}$ penalty ensures good performance when the
true signal is piecewise constant, whereas the
$\lambda_{2}\|\Gamma\beta\|_{2}^{2}$ penalty ensures good performance when the
true signal is smooth over $G$.

Note again that FL and SL in this setting correspond to the GEN penalty with
$\lambda_{2}$ or $\lambda_{1}$ set to zero, respectively. Therefore, GEN’s
superior performance over FL and SL in terms of prediction error is not
surprising, given that the scorer used for CV
$-\frac{1}{n}\|Y-X\hat{\beta}\|_{2}^{2}$ selects for hyperparameters with
stronger prediction performance. However, GEN is also consistently better than
FL or SL in terms of estimation error. This can be understood better by
examining the signal estimates obtained from the three procedures. Figure 7
compare the estimated signals with the true signals defined on the 1D chain
graph. FL recovers the constant regions well but struggles with the smoothly
increasing region, whereas the SL estimate is better in the smoothly
increasing region but cannot reproduce the constant regions of the true
signal. GEN, on the other hand, is able to recover the true signal in all
regions. We can also make the same observations when $G$ is the 2D grid graph
(but they are harder to visualize).

(a) 1D chain

(b) 2D grid

(c) Barbell

Figure 6: Prediction and estimation errors for three graphs as
$\|\Gamma\beta^{*}\|_{\infty}$ and $\|\Gamma\beta^{*}\|_{0}$ vary. Results are
based on 500 resamplings. Vertical bars for each true signal connect the
$25^{\text{th}}$ and $75^{\text{th}}$ percentiles. The lines labeled by FL, SL
and GEN connect the medians of errors. 1D chain: $n=70,p=110$. 2D grid:
$n=150,p=225$. Barbell graph: $n=\frac{2}{3}p$, and $p=60,63,66,69,72,75$.

Figure 7: Left: estimated signals obtained from FL, SL and GEN. Right: true
signal. GEN recovers the true signal well in both the constant and the
smoothly increasing regions.

We also examine the performances of FL, SL and GEN as features become more
correlated and thus $\Sigma$ becomes more ill-conditioned. Figure 8 shows the
estimation errors for the chain graph when $\Sigma$ is the identity matrix
(which is the limit of the Toeplitz covariance matrix as $\rho\to 0$) and when
$\Sigma$ is the Toeplitz covariance matrix with $\rho=0.95$. From our
theoretical results, we expect that when the features are highly correlated,
the performance of the standalone $\lambda_{1}\|\Gamma\beta\|_{1}$ penalty
(FL) should be negatively affected by the ill-conditioned nature of $\Sigma$.
However, the $\lambda_{2}\|\Gamma\beta\|_{2}^{2}$ penalty should improve the
minimum eigenvalue term in the denominators of our error bounds, especially
when $\|\Gamma\beta^{*}\|_{\infty}$ is small and $\lambda_{2}$ can be chosen
to be larger. Such an improvement is not as noticeable when $\Sigma$ is the
identity, which is already well-conditioned.

Figure 8: Side-by-side comparison of the estimation errors for the chain graph
when $\Sigma$ is the identity matrix (left) and when $\Sigma$ has the Toeplitz
structure with $\rho=0.95$ (right). $\|\Gamma\beta^{*}\|_{1}$ is fixed at
$15$. Note the greater divergence between the estimation errors of FL and GEN
when there is higher correlation.

####  4.1.4 Performance comparisons as $n$ and $p$ vary

This section examines the performance of GEN relative to all other methods as
$n$ is fixed and $p$ increases, or as $p$ is fixed and $n$ increases. The
covariance matrix $\Sigma$ is constructed as in Section 4.1.1, and the graphs
we use are again the chain graph, the 2D grid and the barbell graph. The true
signal $\beta^{*}$ is again not sparse, but contains a mix of piecewise
constant regions and smoothly varying regions on the graph $G$ (similar to the
true signals with intermediate values of $\|\Gamma\beta^{*}\|_{\infty}$ and
$\|\Gamma\beta^{*}\|_{0}$ in Figure 4 and Figure 5). Note that we do include
the high-dimensional setting when $n$ is smaller than $p$. As described in
Section 4.1.2, hyperparameters for all methods are chosen based on the best
prediction scores in cross-validation. We only report the estimation errors in
Figure 9, since the prediction errors for all methods show the same trends.

(a) 1D chain

(b) 2D grid

(c) Barbell

Figure 9:  Estimation errors (reported on the log scale) based on 500
resamplings for all estimators as $p$ is fixed ($p=110$ for chain graph,
$p=121$ for 2D grid, $p=66$ for barbell graph) but $n$ increases (left), and
as $n=90$ is fixed but $p$ increases (right). $\sigma=1$ is fixed, and in each
plot $\|\Gamma\beta^{*}\|_{\infty}$ is kept roughly constant. CV yields
$\lambda_{L}$ identically equal to zero for the Lasso estimator, and thus its
performance coincides with that of OLS.

As we can see from Figure 9, GEN consistently has the best performance in
terms of estimation errors (except for the barbell graph when GTV-oracle
performs slightly better for some values of $(n,p)$). With regard to graph-
independent methods, OLS and the Lasso clearly fail to perform well in our
setting, and EN only provides limited improvements in terms of estimation
errors. FL and SL, whose penalties do take into account the graph $G$, perform
significantly better than the previous three methods but never better than
GEN. The performance of GTV, whose penalty depends on the covariance estimate
$\hat{\Sigma}$, is not consistent; it can be reasonable for the barbell graph
but, for the other two graphs, is not very different from OLS, EN and the
Lasso for certain values of $(n,p)$. Interestingly, the performance of GTV-
oracle can surpass that of GEN for the barbell graph; this can be attributed
to the fact that $\Sigma$ is constructed to reflect the graph structure and
hence is a good estimate for the graph $G$ itself. The divergence between the
estimation errors of GTV and GTV-oracle therefore suggests that the covariance
estimation error is not negligible, especially when $p$ is small relative to
$n$. Note that GTV requires much more time than other methods for
hyperparameter selection and model training, as we have discussed in Section
4.1.2.

####  4.1.5 Performance comparisons when $\beta^{*}$ is both sparse and smooth
over $G$

In Section 4.1.3 and Section 4.1.4, we have compared the performances of
various estimators when $\beta^{*}$ is dense. We now consider the case when
$G$ is the chain graph, and $\beta^{*}$ is sparse and has small variations in
its successive entries, as illustrated in the left plot of Figure 10. Such a
signal structure should be more favorable to either FL or SL, and we expect at
least one of them to outperform GEN in this case.

However, as can be seen in Table 3, GEN still has the best performance
compared to all other estimators. FL and SL perform better than the Lasso
estimator, which in turn is better than OLS as expected. Certainly, such a
strong performance relative to the other estimators may depend on the
smoothness and sparsity levels of the true signal. Nonetheless, this example
clearly demonstrates that effectively leveraging the true signal’s smoothness
over $G$ can be more important in reducing prediction and estimation errors
than exploiting its sparsity structure.

Figure 10: Left: Sparse and smooth signal with $p=100$,
$\|\beta^{*}\|_{0}=40$, $\|\Gamma\beta^{*}\|_{\infty}=0.39$. Right: The left
signal is modified to include a spike, so that $\|\Gamma\beta^{*}\|_{\infty}$
increases to 5. We use $\sigma=1,n=80$ and the Toeplitz covariance matrix with
$\rho=0.5$ for $\Sigma$ in this section.

We also consider a slight modification to the previous example, so that we
have a sharp spike in the zero region of the signal. Adding a single spike
should not significantly change the radius $R_{q}$ of the $\ell_{q}$-ball to
which $\Gamma\beta^{*}$ belongs. However, since $\|\Gamma\beta^{*}\|_{\infty}$
is now much larger, CV yields $\lambda_{2}$ identically equal to zero, and GEN
degenerates into FL with $\lambda_{L}=0$. As a result, FL performs better than
GEN (and so does GTV-oracle), although the deterioration of GEN’s performance
is not drastic and GEN still performs better than EN, SL, GTV and the Lasso.
It is therefore a question of interest for future research whether we can
replace the $\ell_{2}$ component of GEN with another penalty that is more
robust to signal spikes, while retaining the benefits of having the $\ell_{2}$
component as discussed in Section 4.1.3.

Table 3: Prediction and estimation errors for the true signals in Figure 10; ‘L’ and ‘R’ denote errors for the left and right true signals respectively. The mean and standard deviation of the errors based on 500 resamplings are shown below. Errors better than GEN’s errors are shown in orange. |  OLS |  L |  EN |  FL |  SL |  GTV |  GTV-oracle |  GEN
---|---|---|---|---|---|---|---|---
L: Est. errors |  $6.92\pm 1.05$ |  $1.98\pm 0.36$ |  $1.98\pm 0.36$ |  $0.56\pm 0.08$ |  $0.40\pm 0.06$ |  $0.87\pm 0.13$ |  $0.38\pm 0.06$ |  $0.27\pm 0.05$
L: Pred. errors |  $38.12\pm 15.33$ |  $2.94\pm 1.35$ |  $2.94\pm 1.35$ |  $0.28\pm 0.12$ |  $0.33\pm 0.14$ |  $0.72\pm 0.26$ |  $0.21\pm 0.11$ |  $0.15\pm 0.08$
R: Est. errors |  $7.25\pm 1.10$ |  $2.95\pm 0.59$ |  $2.95\pm 0.59$ |  $0.88\pm 0.17$ |  $1.75\pm 0.30$ |  $1.31\pm 0.22$ |  $0.75\pm 0.12$ |  $0.97\pm 0.18$
R: Pred. errors |  $41.75\pm 15.95$ |  $6.19\pm 2.93$ |  $6.19\pm 2.93$ |  $0.67\pm 0.26$ |  $2.37\pm 0.93$ |  $1.44\pm 0.50$ |  $0.49\pm 0.19$ |  $0.87\pm 0.32$

###  4.2 Empirical study of the quantity
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)$

In Section 2, if $\Sigma$ is ill-conditioned and we cannot assume
$\gamma_{\text{min}}(\Sigma)$ is bounded away from zero, then we assume that
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)$ may be
greater than $c\lambda_{2}$ or $c\sqrt{\lambda_{2}}$. These assumptions lead
to the bounds (24) and (25), which may allow for consistency in prediction and
estimation respectively.

(a) Chain graph; $\rho=0.8$

(b) Chain graph; $\rho=0.99$

(c) Chain graph; $\rho=0.9985$

(d) Complete graph; $\rho=0.99$

(e) 2D grid; $\Sigma$ from $L+10^{-4}I_{p}$

(f) Barbell; $\Sigma$ from $L+3\times 10^{-4}I_{p}$

Figure 11: Growth of $\gamma_{\min}(\frac{1}{64}\Sigma+\lambda_{2}L)$ as a
function of $\lambda_{2}$ for various choices of $\Sigma$ and $G$ ($p=100$ for
all plots). In (a), (b) and (c), $G$ is the chain graph and $\Sigma$ has
Toeplitz structure with varying $\rho$. In (d), $G$ is the complete graph and
$\Sigma$ has Toeplitz structure. In (e) and (f), we use the 2D grid and
barbell graph, with corresponding (and highly correlated) covariance
structures as in Section 4.1.1.

We conjecture that
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\geq\frac{1}{64}\lambda_{2}$
holds for all $\lambda_{2}\in[0,1]$ under reasonable assumptions about
$(\Sigma,L)$; this implies
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\geq\frac{1}{64}\min(\lambda_{2},1)$.
Figure 11 shows the growth of the quantity
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)$ as a
function of $\lambda_{2}$, for the various types of graphs and covariance
matrices we have considered in Section 4. When $G$ is the chain graph and
$\Sigma$ has the Toeplitz structure, we generally have
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\geq\frac{1}{64}\sqrt{\lambda_{2}}$
for all $\lambda_{2}\in[0,1]$ unless $\rho>0.99$. When $G$ is the 2D grid or
barbell graph and $\Sigma$ is constructed accordingly as in Section 4.1.1, we
can also observe the same trends. Overall, when $\Sigma$ is ill-conditioned
and $\lambda_{2}$ can be chosen to be sufficiently large, the $\ell_{2}$
component of the GEN penalty can significantly improve our error upper bounds.

###  4.3 Real data analysis

In this section, we apply the GEN penalty to a number of real datasets.

####  4.3.1 COVID-19 trend prediction

We consider the problem of predicting the number of COVID-19 cases 14 days in
advance for a given county in California, using a New York Times-curated
COVID-19 dataset. This problem may be of importance for hospitals and local
authorities, as they may wish to anticipate potential spikes in COVID-19 cases
based on current, local data. It is reasonable to assume that the number of
cases $Y_{tc}$ on day $t$ in county $c$ is Poisson-distributed, as in Agosto
and Giudici (2020), Bu et al. (2021) and Cori et al. (2013). In this case, we
can apply the variance-stabilizing Anscombe transform $x\mapsto
2\sqrt{x+\frac{3}{8}}$ to form $\tilde{Y}_{tc}$. Following the modeling
approach in Cori et al. (2013), we may then consider the Gaussian model

$\tilde{Y}_{tc}=\sum_{s=14}^{21}\alpha_{s}Y_{t-s,c}+\epsilon_{tc}$ (49)

where $\epsilon_{tc}\sim N(0,\sigma^{2})$. In order to reduce temporal
correlation between observations, the days $t$ are sampled such that
consecutive time points are at least 7 days apart. We restrict our analysis
(a) to the period from June 2020 to July 2021 to avoid non-stationary effects
in the evolution of the pandemic due to the appearance of new virus strains,
and (b) to the 25 densest counties in California where linear models are
typically a better fit. As in Ngonghala et al. (2022), we use cross-validation
to evaluate the accuracy of our model; 6/7 of our data is used for fitting and
the remaining data is for performance evaluation (i.e. 2 months of data).
Fitting an OLS model based on (49) usually results in a satisfactory fit with
an $R^{2}$ score above 0.8 (see the Appendix).

We hypothesize that for densely populated counties, rising cases in
neighboring counties may further explain a significant fraction of the
remaining variance in the data due to population movements between counties.
To test this hypothesis, we consider a model that incorporates the number of
cases from nearby counties within a two-hop radius of the given county $c$:

$\tilde{Y}_{tc}=\sum_{k\in
N_{2}(c)}\sum_{s=14}^{21}\alpha_{sk}Y_{t-s,k}+\epsilon_{tc}$ (50)

Fitting model (50) is a high-dimensional problem, where the number of
parameters $p$ can be up to 3 times the number of observations $n$, depending
on the county. OLS therefore is not a suitable method for model (50).
Consequently, in this experiment we fit the penalty-based methods (except GTV)
based on (50) and compare with the performance of OLS computed based on (49).
The graph $G$ we consider here is such that two feature vectors are connected
if they are indexed by the same day $t$ and by two adjacent counties, or if
they are indexed by the same county $k$ and two consecutive time points.

We perform this prediction task for each of 25 most densely populated counties
in California. For each county, we report the median root mean square error
(RMSE) computed on the test set in Table 4. As expected, incorporating the
numbers of cases in neighboring counties allows us to outperform OLS based on
(49) in 21 out of 25 counties. The graph-dependent methods FL, SL and GEN
perform better than other methods in 19 out of 25 counties. Among these 19
counties, GEN has the best performance in 7 of them and is therefore a
competitive candidate for this prediction task. The improvement it yields can
be quite substantial; for the county Sutter in particular, GEN reduces the
RMSE by 50% compared to OLS and at least 25% compared to FL and SL.

Table 4: Median RMSE achieved by various methods for 25 counties. OLS is fitted based on model (49), and all other methods are based on (50). The best performances for each county are highlighted in bold. County | OLS | L | EN | FL | SL | GEN
---|---|---|---|---|---|---
Alameda | 1.08 | 1.18 | 1.14 | 0.96 | 0.85 | 0.99
Butte | 3.20 | 1.81 | 1.83 | 1.46 | 1.88 | 2.10
Contra Costa | 1.21 | 1.79 | 1.69 | 3.47 | 2.49 | 3.35
Fresno | 8.22 | 7.06 | 9.77 | 5.25 | 7.71 | 5.95
Los Angeles | 4.92 | 6.92 | 7.34 | 5.28 | 6.06 | 4.56
Marin | 6.11 | 3.92 | 5.49 | 5.47 | 3.38 | 4.18
Merced | 7.94 | 9.75 | 9.03 | 9.26 | 9.68 | 6.08
Napa | 3.93 | 3.97 | 5.53 | 3.92 | 3.90 | 2.96
Orange | 1.84 | 4.44 | 3.97 | 3.34 | 2.08 | 2.56
Placer | 2.20 | 1.35 | 1.89 | 1.27 | 1.55 | 1.53
Riverside | 3.46 | 3.32 | 3.62 | 3.21 | 2.79 | 3.73
Sacramento | 2.23 | 3.11 | 2.53 | 3.61 | 2.31 | 1.66
San Diego | 1.43 | 1.67 | 1.00 | 0.98 | 0.88 | 0.91
San Francisco | 1.41 | 2.23 | 1.05 | 1.23 | 1.33 | 1.39
San Joaquin | 3.64 | 3.43 | 3.43 | 3.93 | 5.24 | 5.41
San Mateo | 1.44 | 2.34 | 2.45 | 1.68 | 1.56 | 1.75
Santa Barbara | 2.84 | 2.02 | 2.02 | 2.02 | 2.01 | 3.71
Santa Clara | 1.14 | 2.04 | 1.85 | 1.05 | 1.10 | 1.03
Santa Cruz | 6.56 | 3.86 | 4.62 | 3.55 | 4.17 | 4.59
Solano | 2.10 | 3.86 | 3.72 | 2.03 | 2.93 | 2.62
Sonoma | 1.74 | 3.47 | 3.60 | 2.78 | 2.62 | 2.67
Stanislaus | 9.29 | 6.41 | 9.17 | 4.55 | 4.76 | 4.88
Sutter | 4.07 | 7.51 | 7.51 | 4.33 | 2.58 | 1.94
Ventura | 2.02 | 1.22 | 1.20 | 1.23 | 1.16 | 1.36
Yolo | 3.43 | 5.13 | 4.54 | 1.93 | 2.45 | 1.79

####  4.3.2 Detection of Alzheimer’s disease

We test the GEN penalty’s performance in detecting Alzheimer’s disease, using
an MRI dataset available on Kaggle. The task is to classify whether the MRI
images in the dataset show signs of dementia. Since the responses are binary,
we need to consider the logistic extension (6) of our method as well as that
of all other methods. The original dataset has images labeled with moderate,
mild, very mild and no dementia, but we exclude the moderate cases due to the
small number of training samples. We also exclude the very mild cases since
the images may be too similar to those with no dementia, thus leading to lower
prediction accuracy for all methods.

Since the features are 2D MRI images, it is natural to use the 2D grid graph
as our graph $G$, which is of size $p=32\times 32=1024$ (we compress the
original images to this size for computational convenience). We use the first
800 images with no dementia and 400 images with mild dementia in the original
dataset. Out of these 1200 images, $n=480$ images are used as training data
(note that $n<p$), 480 images are use for hyperparameter tuning, and the other
240 images constitute our testing data. For computation, we use ECOS for all
methods. Since GTV requires at least a 3D grid search for hyperparameter
tuning, it is too slow to be considered for this experiment. All other methods
take at most 5 seconds of training time.

The classification accuracies for all methods except GTV are reported in Table
5. As expected, GEN shows better prediction performance than all other methods
in consideration.

Table 5: Prediction accuracies for classification of Alzhemer’s disease status. Here, OLS is replaced by logistic regression (LR), and the logistic extensions of all penalty-based methods (except GTV) are used. |  LR |  L |  EN |  FL |  SL |  GEN
---|---|---|---|---|---|---
Accuracy |  82.08% |  90.0 % |  92.50% |  91.25% |  92.08% |  92.92%

####  4.3.3 Estimation of crime patterns in Chicago

Consider the task of uncovering crime trends over time across the 77
communities of Chicago (which we denote by the set $\mathcal{C}$). Statistics
on the number of crimes per community between 2004 and 2022 are available on
the city’s data portal. The monthly crime rates (which are defined here as the
number of crimes per 100,000 inhabitants) vary over the years and across the
communities, and they are also subject to significant seasonal effects.
Additional details on the nature of the data and preprocessing are provided in
the Appendix.

Let $Y_{my}^{(c)}$ denote the crime rate for community $c\in\mathcal{C}$,
month $m$ and year $y$. Since we are working with count data, it is reasonable
to pre-process the data by applying the Anscombe transform to $Y_{my}^{(c)}$
to form $\tilde{Y}_{my}^{(c)}$. We then consider the following additive
Gaussian model

$\tilde{Y}_{my}^{(c)}=\sum_{i=1}^{12}\alpha_{i}\mathbf{1}[m=i]+\sum_{j=2004}^{2022}\beta_{j}\mathbf{1}[y=j]+\sum_{c\in\mathcal{C}}\gamma_{k}\mathbf{1}[k=c]+\epsilon_{my}^{(c)}$
(51)

where $\epsilon_{my}^{(c)}\sim N(0,\sigma^{2})$. While our design matrix here
is not equal to identity as in the trend filtering case, note that it contains
“one-hot” encodings rather than i.i.d. rows from some distribution
$\mathbb{P}$. The parameters $(\alpha,\beta)$ naturally exhibit temporal
smoothness, since we expect them to vary smoothly over time. The community
offset parameter $\gamma$, on the other hand, should exhibit spatial
smoothness, as we expect neighboring communities to have similar offsets. To
define our GEN penalty, we encode these prior beliefs in a regularizing graph
$G$ with 3 disconnected components: one chain graph reinforcing the temporal
smoothness of the month coefficients, another chain graph for that of the
years, and a third component encoding neighborhood adjacency. We compare our
method’s performance with all other methods except GTV.

Figure 12 compares the prediction performance (reported using RMSE computed on
held-out data across 40 independent trials) for all methods. Here, performance
is assessed for different data regimes: while the original dataset contains
17,094 observations, we use a fraction $\alpha\in\\{0.5\%,1\%,2\%,5\%\\}$ of
data for estimation of the $p=108$ parameters in our model ($\alpha=0.5\%$ and
$\alpha=1\%$ correspond to $p>n$ and $p\approx n$ respectively). As shown in
Figure 12, GEN performs consistently better than all other methods, especially
in the data-sparse regime. While in this example we are more interested in the
estimation of crime patterns rather than prediction (note that the model (51)
cannot be used to predict crime rates beyond 2022), Figure 12 provides
evidence for GEN’s superior performance and can be of interest if we are given
a dataset with many missing values that require data imputation.

Figure 12: RMSE achieved by different estimators as the proportion $\alpha$ of
data used for training varies. Figure 13: Visualization of the community
offsets $\gamma$ produced by different estimators. Note that GEN produces
smoother estimates with greater magnitudes.

Figure 13 and Figure 14 visualize the estimates of the community offsets
$\gamma$ and the temporal parameters $(\alpha,\beta)$ respectively. Note that
the estimate of $\gamma$ obtained from GEN contains fewer zero entries and is
significantly smoother when compared with other methods. We can clearly see
that GEN divides the communities into clusters with similar community offsets.
Only Smooth Lasso provides an estimate of $\gamma$ that is close to GEN’s, but
interestingly the estimates obtained by GEN tend to be greater in magnitude.
From Figure 14, we can see that GEN, FL and SL produce smooth estimates to
show that crime rates tend to decrease in the colder months and that there is
a general reduction in crime rates between 2002 and 2022.

(a) Estimates of $\alpha$

(b) Estimates of $\beta$

Figure 14: Visualization of the estimates of the monthly parameters $\alpha$
and the yearly parameters $\beta$’s produced by different methods. GEN, FL and
SL produce smoother estimates relative to graph-independent methods.

##  5 Conclusion

In this paper, we have introduced an $\ell_{1}+\ell_{2}$ penalty in which both
the $\ell_{1}$ and $\ell_{2}$ components are defined with respect to a given
graph $G$. We have shown that the assumption that $\beta^{*}$ is piecewise
constant or smooth with respect to $G$ can be sufficient for our estimator to
be consistent in prediction and estimation, even in the high-dimensional
setting. The $\ell_{1}$ component of our penalty improves estimation of the
piecewise constant regions of the true signal, whereas the $\ell_{2}$
component improves estimation of the true signal’s smooth regions and
alleviates the negative effects of a highly correlated design. We also show
that our estimator can be computed efficiently, even for large-scale problems
where choosing hyperparameters via cross-validation may be challenging.

## References

  * Agosto and Giudici (2020) Agosto, A. and Giudici, P. (2020) A Poisson autoregressive model to understand COVID-19 contagion dynamics.  Risks, 8, 77.
  * Anderson (1955) Anderson, T. W. (1955) The integral of a symmetric unimodal function over a symmetric convex set and some probability inequalities.  Proceedings of the American Mathematical Society, 6, 170–176.
  * Bickel and Levina (2008) Bickel, P. J. and Levina, E. (2008) Covariance regularization by thresholding.  Annals of Statistics, 36, 2577–2604.
  * Boyd et al. (2004) Boyd, S., Boyd, S. P. and Vandenberghe, L. (2004) Convex optimization.  Cambridge University Press.
  * Bu et al. (2021) Bu, F., Aiello, A. E., Volfovsky, A. and Xu, J. (2021) Likelihood-based inference for partially observed stochastic epidemics with individual heterogeneity.  arXiv preprint arXiv:2112.07892.
  * Bühlmann and van de Geer (2011) Bühlmann, P. and van de Geer, S. (2011) Statistics for high-dimensional data: methods, theory and applications.  Springer Science & Business Media.
  * Cai and Zhou (2012) Cai, T. T. and Zhou, H. H. (2012) Optimal rates of convergence for sparse covariance matrix estimation.  Annals of Statistics, 2389–2420.
  * Chung and Graham (1997) Chung, F. R. and Graham, F. C. (1997) Spectral graph theory, vol. 92.  American Mathematical Society.
  * Cori et al. (2013) Cori, A., Ferguson, N. M., Fraser, C. and Cauchemez, S. (2013) A new framework and software to estimate time-varying reproduction numbers during epidemics.  American Journal of Epidemiology, 178, 1505–1512.
  * Domahidi et al. (2013) Domahidi, A., Chu, E. and Boyd, S. (2013) ECOS: An SOCP solver for embedded systems.  In 2013 European Control Conference (ECC), 3071–3076. IEEE.
  * Elad et al. (2007) Elad, M., Milanfar, P. and Rubinstein, R. (2007) Analysis versus synthesis in signal priors.  Inverse problems, 23, 947.
  * Guntuboyina et al. (2020) Guntuboyina, A., Lieu, D., Chatterjee, S. and Sen, B. (2020) Adaptive risk bounds in univariate total variation denoising and trend filtering.  Annals of Statistics, 48, 205–229.
  * Hebiri and van de Geer (2011) Hebiri, M. and van de Geer, S. (2011) The Smooth-Lasso and other $\ell_{1}+\ell_{2}$-penalized methods.  Electronic Journal of Statistics, 5, 1184–1226.
  * Hütter and Rigollet (2016) Hütter, J.-C. and Rigollet, P. (2016) Optimal rates for total variation denoising.  In Conference on Learning Theory, 1115–1146. PMLR.
  * Kim et al. (2009) Kim, S.-J., Koh, K., Boyd, S. and Gorinevsky, D. (2009) $\ell_{1}$ trend filtering.  SIAM review, 51, 339–360.
  * Kolaczyk (2009) Kolaczyk, E. (2009) Statistical analysis of network data: Methods and models.  Springer Series In Statistics, 386.
  * Land and Friedman (1997) Land, S. R. and Friedman, J. H. (1997) Variable fusion: A new adaptive signal regression method.
  * Li and Li (2008) Li, C. and Li, H. (2008) Network-constrained regularization and variable selection for analysis of genomic data.  Bioinformatics, 24, 1175–1182.
  * Li et al. (2019) Li, T., Levina, E. and Zhu, J. (2019) Prediction models for network-linked data.  Annals of Applied Statistics, 13, 132–164.
  * Li et al. (2018) Li, Y., Mark, B., Raskutti, G. and Willett, R. (2018) Graph-based regularization for regression problems with highly-correlated designs.  In 2018 IEEE Global Conference on Signal and Information Processing (GlobalSIP), 740–742. IEEE.
  * Ngonghala et al. (2022) Ngonghala, C. N., Taboe, H. B., Safdar, S. and Gumel, A. B. (2022) Unraveling the dynamics of the Omicron and Delta variants of the 2019 coronavirus in the presence of vaccination, mask usage, and antiviral treatment.  Applied Mathematical Modelling.
  * Ollier and Viallon (2017) Ollier, E. and Viallon, V. (2017) Regression modelling on stratified data with the lasso.  Biometrika, 104, 83–96.
  * Ortelli and van de Geer (2021) Ortelli, F. and van de Geer, S. (2021) Prediction bounds for higher order total variation regularized least squares.  Annals of Statistics, 49, 2755–2773.
  * Osgood (2000) Osgood, D. W. (2000) Poisson-based regression analysis of aggregate crime rates.  Journal of Quantitative Criminology, 16, 21–43.
  * Raskutti et al. (2010) Raskutti, G., Wainwright, M. J. and Yu, B. (2010) Restricted eigenvalue properties for correlated Gaussian designs.  Journal of Machine Learning Research, 11, 2241–2259.
  * Raskutti et al. (2011) — (2011) Minimax rates of estimation for high-dimensional linear regression over $\ell_{q}$-balls.  IEEE transactions on information theory, 57, 6976–6994.
  * Rudin (1974) Rudin, W. (1974) Real and Complex Analysis.  McGraw-Hill Inc.,.
  * Segal et al. (2003) Segal, M. R., Dahlquist, K. D. and Conklin, B. R. (2003) Regression approaches for microarray data analysis.  Journal of Computational Biology, 10, 961–980.
  * She (2008) She, Y. (2008) Sparse regression with exact clustering.  Stanford University.
  * Tibshirani (1996) Tibshirani, R. (1996) Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 58, 267–288.
  * Tibshirani et al. (2005) Tibshirani, R., Saunders, M., Rosset, S., Zhu, J. and Knight, K. (2005) Sparsity and smoothness via the fused lasso.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67, 91–108.
  * Tibshirani and Taylor (2011) Tibshirani, R. J. and Taylor, J. (2011) The solution path of the generalized lasso.  Annals of Statistics, 39, 1335–1371.
  * Toharudin et al. (2020) Toharudin, T., Caraka, R. E., Chen, R. C., Nugroho, N. T., Sueb, M., Jaya, I. and Pontoh, R. (2020) Bayesian Poisson model for COVID-19 in West Java, Indonesia.  Sylwan, 164, 279–290.
  * Tseng (2001) Tseng, P. (2001) Convergence of a block coordinate descent method for nondifferentiable minimization.  Journal of Optimization Theory and Applications, 109, 475–494.
  * Wainwright (2019) Wainwright, M. J. (2019) High-dimensional statistics: A non-asymptotic viewpoint, vol. 48.  Cambridge University Press.
  * Wang et al. (2015) Wang, Y.-X., Sharpnack, J., Smola, A. and Tibshirani, R. (2015) Trend filtering on graphs.  In Artificial Intelligence and Statistics, 1042–1050. PMLR.
  * Xin et al. (2014) Xin, B., Kawahara, Y., Wang, Y. and Gao, W. (2014) Efficient generalized fused lasso and its application to the diagnosis of Alzheimer’s disease.  In Proceedings of the AAAI Conference on Artificial Intelligence, vol. 28.
  * Zou and Hastie (2005) Zou, H. and Hastie, T. (2005) Regularization and variable selection via the elastic net.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67, 301–320.

##  Appendix A Proofs of theoretical results

We restate our theorems in this appendix for convenience.

######  Theorem A.1 (Theorem 2.1).

Fix $\delta>0$ and choose
$\lambda_{1}=32\sigma\rho(\Gamma)\sqrt{\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}}$, $\lambda_{2}\leq\frac{\lambda_{1}}{8\|\Gamma\beta^{*}\|_{\infty}}$.
Given any set $S$ satisfying both

$\frac{144\gamma_{\text{max}}(\Sigma)(\sqrt{n_{c}}+\delta)^{2}}{n}+\frac{36\lambda_{1}^{2}|S|k_{S}^{-2}}{\sigma^{2}}\leq\frac{1}{2}\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)$
(52)

and

$\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}\leq\frac{\sigma^{2}}{18}$ (53)

with probability at least $1-c_{1}\exp(-nc_{2})-\frac{2}{m}-e^{-\delta^{2}/2}$
we have

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)\frac{n_{c}+\delta^{2}}{n}+\lambda_{1}^{2}|S|k_{S}^{-2}}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}+\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}$

and

$\|\hat{\beta}-\beta^{*}\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)\frac{n_{c}+\delta^{2}}{n}+\lambda_{1}^{2}|S|k_{S}^{-2}}{\gamma_{\text{min}}^{2}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}+\frac{\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}$

Proof By definition,

$\hat{\beta}=\arg\min_{\beta}\frac{1}{n}\|Y-X\beta\|_{2}^{2}+\lambda_{1}\|\Gamma\beta\|_{1}+\lambda_{2}\|\Gamma\beta\|_{2}^{2}$

We can also rewrite our estimator as:

$\hat{\beta}=\arg\min_{\beta}\frac{1}{n}\|\tilde{Y}-\tilde{X}\beta\|_{2}^{2}+\lambda_{1}\|\Gamma\beta\|_{1}$

Using subdifferential calculus, we can see that $\hat{\beta}$ must satisfy

$\frac{2\tilde{X}^{T}(\tilde{Y}-\tilde{X}\hat{\beta})}{n}=\lambda_{1}\Gamma^{T}\text{sign}(\Gamma\hat{\beta})$

where

$[\mathop{\rm sign}(x)]_{i}=\left\\{\begin{array}[]{rl}1&\text{if}\
x_{i}>0,\\\ \text{any value in }\left[-1,1\right]&\text{if}\ x_{i}=0,\\\
-1&\text{if}\ x_{i}<0\,.\end{array}\right.$

Hence, we obtain

$\frac{2}{n}\hat{\beta}^{T}\tilde{X}^{T}(\tilde{Y}-\tilde{X}\hat{\beta})=\lambda_{1}\hat{\beta}^{T}\Gamma^{T}\mathop{\rm
sign}(\Gamma\hat{\beta})=\lambda_{1}\|\Gamma\hat{\beta}\|_{1}$

and for any $\beta\in\mathbb{R}^{p}$,

$\frac{2}{n}\beta^{T}\tilde{X}^{T}(\tilde{Y}-\tilde{X}\hat{\beta})=\lambda_{1}\beta^{T}\Gamma^{T}\mathop{\rm
sign}(\Gamma\hat{\beta})\leq\lambda_{1}\|\Gamma\beta\|_{1}$

By subtracting the previous equality from the inequality above, for any
$\beta\in\mathbb{R}^{p}$ we have

$\frac{2}{n}(\beta-\hat{\beta})^{T}\tilde{X}^{T}(\tilde{Y}-\tilde{X}\hat{\beta})\leq\lambda_{1}(\|\Gamma\beta\|_{1}-\|\Gamma\hat{\beta}\|_{1})$

Since $\tilde{Y}=\tilde{X}\beta^{*}+\tilde{\epsilon}$,

$\displaystyle\frac{2}{n}(\hat{\beta}-\beta)^{T}\tilde{X}^{T}\tilde{X}(\hat{\beta}-\beta^{*})$
$\displaystyle\leq\frac{2}{n}\tilde{\epsilon}^{T}\tilde{X}(\hat{\beta}-\beta)+\lambda_{1}(\|\Gamma\beta\|_{1}-\|\Gamma\hat{\beta}\|_{1})$
$\displaystyle=\frac{2}{n}\epsilon^{T}X(\hat{\beta}-\beta)-2\lambda_{2}(\beta^{*})^{T}\Gamma^{T}\Gamma(\hat{\beta}-\beta)+\lambda_{1}(\|\Gamma\beta\|_{1}-\|\Gamma\hat{\beta}\|_{1})$
$\displaystyle\leq\frac{2}{n}\epsilon^{T}X(\hat{\beta}-\beta)+2\lambda_{2}\|\Gamma\beta^{*}\|_{\infty}\|\Gamma(\hat{\beta}-\beta^{*})\|_{1}+\lambda_{1}(\|\Gamma\beta\|_{1}-\|\Gamma\hat{\beta}\|_{1})$
$\displaystyle\leq\frac{2}{n}\epsilon^{T}X(\hat{\beta}-\beta)+\frac{\lambda_{1}}{4}\|\Gamma(\hat{\beta}-\beta)\|_{1}+\lambda_{1}(\|\Gamma\beta\|_{1}-\|\Gamma\hat{\beta}\|_{1})$

where the last inequality follows if we choose
$\lambda_{2}\leq\frac{\lambda_{1}}{8\|\Gamma\beta^{*}\|_{\infty}}$.

We wish to bound $\frac{2}{n}\epsilon^{T}X(\hat{\beta}-\beta)$. As
$\Pi\in\mathbb{R}^{p\times p}$ denotes the projection matrix onto the kernel
of $\Gamma$, we have $I_{p}=\Pi+\Gamma^{\dagger}\Gamma$. Hence,

$\displaystyle\frac{2}{n}\epsilon^{T}X(\hat{\beta}-\beta)$
$\displaystyle=\frac{2}{n}\epsilon^{T}X\Pi(\hat{\beta}-\beta)+\frac{2}{n}\epsilon^{T}X\Gamma^{\dagger}\Gamma(\hat{\beta}-\beta)$
(54) $\displaystyle\leq\frac{2}{n}\|\Pi
X^{T}\epsilon\|_{2}\|\hat{\beta}-\beta\|_{2}+\frac{2}{n}\|(\Gamma^{\dagger})^{T}X^{T}\epsilon\|_{\infty}\|\Gamma(\hat{\beta}-\beta)\|_{1}$
$\displaystyle\leq\frac{2}{n}\|\Pi
X^{T}\epsilon\|_{2}\|\hat{\beta}-\beta\|_{2}+\frac{\lambda_{1}}{4}\|\Gamma(\hat{\beta}-\beta)\|_{1}$

where the last inequality follows if we choose
$\lambda_{1}\geq\frac{8}{n}\|(\Gamma^{\dagger})^{T}X^{T}\epsilon\|_{\infty}$
(with high probability).

We obtain the bound:

$\begin{split}\frac{2}{n}(\hat{\beta}-\beta)^{T}\tilde{X}^{T}\tilde{X}(\hat{\beta}-\beta^{*})\leq&\frac{2}{n}\|\Pi
X^{T}\epsilon\|_{2}\|\hat{\beta}-\beta\|_{2}+\frac{\lambda_{1}}{2}\|\Gamma(\hat{\beta}-\beta)\|_{1}\\\
&+\lambda_{1}\|\Gamma\beta\|_{1}-\lambda_{1}\|\Gamma\hat{\beta}\|_{1}\end{split}$
(55)
$\frac{2}{n}(\hat{\beta}-\beta)^{T}\tilde{X}^{T}\tilde{X}(\hat{\beta}-\beta^{*})\leq\frac{2}{n}\|\Pi
X^{T}\epsilon\|_{2}\|\hat{\beta}-\beta\|_{2}+\frac{\lambda_{1}}{2}\|\Gamma(\hat{\beta}-\beta)\|_{1}+\lambda_{1}\|\Gamma\beta\|_{1}-\lambda_{1}\|\Gamma\hat{\beta}\|_{1}$

For any $S\subseteq[m]$:

$\displaystyle\frac{\lambda_{1}}{2}\|\Gamma(\hat{\beta}-\beta)\|_{1}+\lambda_{1}\|\Gamma\beta\|_{1}-\lambda_{1}\|\Gamma\hat{\beta}\|_{1}$
$\displaystyle\leq\frac{\lambda_{1}}{2}\|(\Gamma\hat{\beta}-\Gamma\beta)_{S}\|_{1}+\frac{\lambda_{1}}{2}\|(\Gamma\hat{\beta})_{-S}\|_{1}+\frac{\lambda_{1}}{2}\|(\Gamma\beta)_{-S}\|_{1}+\lambda_{1}\|\Gamma\beta\|_{1}-\lambda_{1}\|\Gamma\hat{\beta}\|_{1}$
$\displaystyle\leq\frac{3\lambda_{1}}{2}\|(\Gamma\hat{\beta}-\Gamma\beta)_{S}\|_{1}+\frac{3\lambda_{1}}{2}\|(\Gamma\beta)_{-S}\|_{1}-\frac{\lambda_{1}}{2}\|(\Gamma\hat{\beta})_{-S}\|_{1}$
$\displaystyle\leq\frac{3\lambda_{1}}{2}\|(\Gamma\hat{\beta}-\Gamma\beta)_{S}\|_{1}+2\lambda_{1}\|(\Gamma\beta)_{-S}\|_{1}-\frac{\lambda_{1}}{2}\|(\Gamma\hat{\beta}-\Gamma\beta)_{-S}\|_{1}$
$\displaystyle\leq
2\lambda_{1}\|(\Gamma\hat{\beta}-\Gamma\beta)_{S}\|_{1}+2\lambda_{1}\|(\Gamma\beta)_{-S}\|_{1}-\frac{\lambda_{1}}{2}\|\Gamma\hat{\beta}-\Gamma\beta\|_{1}$

and so we have

$\displaystyle\frac{2}{n}(\hat{\beta}-\beta)^{T}\tilde{X}^{T}\tilde{X}(\hat{\beta}-\beta^{*})+\frac{\lambda_{1}}{2}\|\Gamma\hat{\beta}-\Gamma\beta\|_{1}$
$\displaystyle\leq\frac{2}{n}\|\Pi
X^{T}\epsilon\|_{2}\|\hat{\beta}-\beta\|_{2}+2\lambda_{1}\|(\Gamma\hat{\beta}-\Gamma\beta)_{S}\|_{1}+2\lambda_{1}\|(\Gamma\beta)_{-S}\|_{1}$
$\displaystyle\leq 2\left(\frac{1}{n}\|\Pi
X^{T}\epsilon\|_{2}+\frac{\lambda_{1}\sqrt{|S|}}{k_{S}}\right)\|\hat{\beta}-\beta\|_{2}+2\lambda_{1}\|(\Gamma\beta)_{-S}\|_{1}$
$\displaystyle\leq
2\left(\sqrt{2\sigma^{2}\gamma_{\text{max}}(\Sigma)}\frac{\sqrt{n_{c}}+\delta}{\sqrt{n}}+\frac{\lambda_{1}\sqrt{|S|}}{k_{S}}\right)\|\hat{\beta}-\beta\|_{2}+2\lambda_{1}\|(\Gamma\beta)_{-S}\|_{1}$

with high probability, where we used the definition of $k_{S}$ and Lemma A.3.
If we set $\beta=\beta^{*}$, we obtain

$\begin{split}\frac{1}{n}\|\tilde{X}(\hat{\beta}-\beta^{*})\|_{2}^{2}+\frac{\lambda_{1}}{4}\|\Gamma\hat{\beta}-\Gamma\beta^{*}\|_{1}\leq&\left(\sqrt{2\sigma^{2}\gamma_{\text{max}}(\Sigma)}\frac{\sqrt{n_{c}}+\delta}{\sqrt{n}}+\frac{\lambda_{1}\sqrt{|S|}}{k_{S}}\right)\|\hat{\beta}-\beta^{*}\|_{2}\\\
&+\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}\end{split}$ (56)

which implies

$\lambda_{1}\|\Gamma\hat{\beta}-\Gamma\beta^{*}\|_{1}\leq
4\left(\sqrt{2\sigma^{2}\gamma_{\text{max}}(\Sigma)}\frac{\sqrt{n_{c}}+\delta}{\sqrt{n}}+\frac{\lambda_{1}\sqrt{|S|}}{k_{S}}\right)\|\hat{\beta}-\beta^{*}\|_{2}+4\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}$

or that

$\displaystyle 576\rho^{2}(\Gamma)\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}\|\Gamma\hat{\beta}-\Gamma\beta^{*}\|_{1}^{2}=\frac{576}{1024}\frac{\lambda_{1}^{2}}{\sigma^{2}}\|\Gamma\hat{\beta}-\Gamma\beta^{*}\|_{1}^{2}$
$\displaystyle\leq
18\frac{\lambda_{1}^{2}}{\sigma^{2}}\left(\sqrt{2\sigma^{2}\gamma_{\text{max}}(\Sigma)}\frac{\sqrt{n_{c}}+\delta}{\lambda_{1}\sqrt{n}}+\frac{\sqrt{|S|}}{k_{S}}\right)^{2}\|\hat{\beta}-\beta^{*}\|_{2}^{2}+18\frac{\lambda_{1}^{2}}{\sigma^{2}}\|(\Gamma\beta^{*})_{-S}\|_{1}^{2}$
$\displaystyle\leq\left(72\gamma_{\text{max}}(\Sigma)\frac{(\sqrt{n_{c}}+\delta)^{2}}{n}+36\frac{\lambda_{1}^{2}|S|k_{S}^{-2}}{\sigma^{2}}\right)\|\hat{\beta}-\beta^{*}\|_{2}^{2}+18\frac{\lambda_{1}^{2}}{\sigma^{2}}\|(\Gamma\beta^{*})_{-S}\|_{1}^{2}$
$\displaystyle\leq\left(72\gamma_{\text{max}}(\Sigma)\frac{(\sqrt{n_{c}}+\delta)^{2}}{n}+36\frac{\lambda_{1}^{2}|S|k_{S}^{-2}}{\sigma^{2}}\right)\|\hat{\beta}-\beta^{*}\|_{2}^{2}+\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}$
(57)

where we used the condition (53). Now if we apply Corollary A.1 to (56), we
have

$\begin{split}&(\hat{\beta}-\beta^{*})^{T}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)(\hat{\beta}-\beta^{*})\leq\left(\sqrt{2\sigma^{2}\gamma_{\text{max}}(\Sigma)}\frac{\sqrt{n_{c}}+\delta}{\sqrt{n}}+\frac{\lambda_{1}\sqrt{|S|}}{k_{S}}\right)\|\hat{\beta}-\beta^{*}\|_{2}\\\
&+\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}+\frac{72\gamma_{\text{max}}(\Sigma)n_{c}}{n}\|\hat{\beta}-\beta^{*}\|_{2}^{2}+576\rho^{2}(\Gamma)\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}\|\Gamma\hat{\beta}-\Gamma\beta^{*}\|_{1}^{2}\end{split}$

which, by (57) and the inequality
$\frac{n_{c}}{n}\leq\frac{(\sqrt{n_{c}}+\delta)^{2}}{n}$, implies

$\begin{split}&(\hat{\beta}-\beta^{*})^{T}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)(\hat{\beta}-\beta^{*})\leq\left(\sqrt{2\sigma^{2}\gamma_{\text{max}}(\Sigma)}\frac{\sqrt{n_{c}}+\delta}{\sqrt{n}}+\frac{\lambda_{1}\sqrt{|S|}}{k_{S}}\right)\|\hat{\beta}-\beta^{*}\|_{2}\\\
&+2\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}+\left(\frac{144\gamma_{\text{max}}(\Sigma)(\sqrt{n_{c}}+\delta)^{2}}{n}+\frac{36\lambda_{1}^{2}|S|k_{S}^{-2}}{\sigma^{2}}\right)\|\hat{\beta}-\beta^{*}\|_{2}^{2}\end{split}$

If we now apply the condition (52), we obtain

$\begin{split}(\hat{\beta}-\beta^{*})^{T}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)&(\hat{\beta}-\beta^{*})\leq\left(\sqrt{2\sigma^{2}\gamma_{\text{max}}(\Sigma)}\frac{\sqrt{n_{c}}+\delta}{\sqrt{n}}+\frac{\lambda_{1}\sqrt{|S|}}{k_{S}}\right)\|\hat{\beta}-\beta^{*}\|_{2}\\\
&+2\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}+\frac{1}{2}\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\|\hat{\beta}-\beta^{*}\|_{2}^{2}\end{split}$

which, by using
$\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\|\hat{\beta}-\beta^{*}\|_{2}^{2}\leq(\hat{\beta}-\beta^{*})^{T}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)(\hat{\beta}-\beta^{*})$,
implies both

$\begin{split}\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)\|\hat{\beta}-\beta^{*}\|_{2}^{2}\leq&2\left(\sqrt{2\sigma^{2}\gamma_{\text{max}}(\Sigma)}\frac{\sqrt{n_{c}}+\delta}{\sqrt{n}}+\frac{\lambda_{1}\sqrt{|S|}}{k_{S}}\right)\|\hat{\beta}-\beta^{*}\|_{2}\\\
&+4\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}\end{split}$ (58)

and

$\begin{split}(\hat{\beta}-\beta^{*})^{T}&\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)(\hat{\beta}-\beta^{*})\leq
4\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}\\\
&+2\frac{\sqrt{2\sigma^{2}\gamma_{\text{max}}(\Sigma)}\frac{\sqrt{n_{c}}+\delta}{\sqrt{n}}+\frac{\lambda_{1}\sqrt{|S|}}{k_{S}}}{\sqrt{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}}\sqrt{(\hat{\beta}-\beta^{*})^{T}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)(\hat{\beta}-\beta^{*})}\end{split}$
(59)

The error bounds follow from (58) and (59) if we note that $x^{2}-bx-c\leq 0$
implies $x^{2}\leq 4\max(b^{2},c)\leq 4(b^{2}+c)$, for $b,c>0$. ∎

######  Theorem A.2 (Theorem 2.2).

Let $\Gamma$ be the incidence matrix of the $p$-vertex chain graph, and fix
$\delta>0$. With an appropriate choice of $\lambda_{1}$ and
$\lambda_{2}\leq\frac{\lambda_{1}}{8\|\Gamma\beta^{*}\|_{\infty}}$, with high
probability we have

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}(\frac{1}{64}\Sigma+\lambda_{2}L)}\frac{1+\delta^{2}}{n}+\frac{(\sigma^{2}\gamma_{\text{max}}(\Sigma)\|\Gamma\beta^{*}\|_{1})^{2/3}}{\gamma_{\text{min}}^{1/3}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{(p\log
p)^{1/3}}{n^{2/3}}$ (60)
$\|\hat{\beta}-\beta^{*}\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)}{\gamma_{\text{min}}^{2}(\frac{1}{64}\Sigma+\lambda_{2}L)}\frac{1+\delta^{2}}{n}+\frac{(\sigma^{2}\gamma_{\text{max}}(\Sigma)\|\Gamma\beta^{*}\|_{1})^{2/3}}{\gamma_{\text{min}}^{4/3}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{(p\log
p)^{1/3}}{n^{2/3}}$ (61)

provided that the RHS of (60) is smaller than $C\sigma^{2}$.

Proof The proof is identical to that of Theorem A.1 up to (54). However, we
need to bound
$\frac{2}{n}\epsilon^{T}X\Gamma^{\dagger}\Gamma(\hat{\beta}-\beta)$
differently.

Let $\Gamma=U\Xi V^{T}$ be the singular value decomposition of $\Gamma$, and
let $\xi_{1},\dots,\xi_{p-1}$ be the nonzero singular values of $\Gamma$. Let
$u_{1},\dots,u_{m}$ and $v_{1},\dots,v_{p}$ denote the columns of $U$ and $V$.
Denote by $V_{[k]}\in\mathbb{R}^{p\times k}$ the matrix containing the first
$k$ columns of $V$ ($k$ is to be specified later) and
$V_{-[k]}\in\mathbb{R}^{p\times(p-k)}$ the matrix containing the other $p-k$
columns of $V$. Define the projection matrix
$P_{[k]}:=V_{[k]}V_{[k]}^{T}\in\mathbb{R}^{p\times p}$.

Noting that $\Gamma^{\dagger}\Gamma$ is a projection matrix, we have:

$\displaystyle\frac{2}{n}\epsilon^{T}X\Gamma^{\dagger}\Gamma(\hat{\beta}-\beta)$
$\displaystyle=\frac{2}{n}\epsilon^{T}XP_{[k]}\Gamma^{\dagger}\Gamma(\hat{\beta}-\beta)+\frac{2}{n}\epsilon^{T}X(I_{p}-P_{[k]})\Gamma^{\dagger}\Gamma(\hat{\beta}-\beta)$
$\displaystyle\leq\frac{2}{n}\|P_{[k]}X^{T}\epsilon\|_{2}\|\Gamma^{\dagger}\Gamma(\hat{\beta}-\beta)\|_{2}+\frac{2}{n}\|(\Gamma^{\dagger})^{T}(I_{p}-P_{[k]})X^{T}\epsilon\|_{\infty}\|\Gamma(\hat{\beta}-\beta)\|_{1}$
$\displaystyle\leq\frac{2}{n}\|P_{[k]}X^{T}\epsilon\|_{2}\|\hat{\beta}-\beta\|_{2}+\frac{\lambda_{1}}{2}\|\Gamma(\hat{\beta}-\beta)\|_{1}$
(62)

if we choose
$\lambda_{1}\geq\frac{8}{n}\|(\Gamma^{\dagger})^{T}(I_{p}-P_{[k]})X^{T}\epsilon\|_{\infty}$
with high probability.

In order to choose $k$, we need to bound
$\frac{8}{n}\|(\Gamma^{\dagger})^{T}(I_{p}-P_{[k]})X^{T}\epsilon\|_{\infty}$.
Let $s_{1}^{\prime},\dots,s_{m}^{\prime}$ be the columns of
$(I_{p}-P_{[k]})\Gamma^{\dagger}$. Let $e_{j}$, $j\in[m]$, denote the
$j^{\text{th}}$ canonical basis element. As in the proof of Theorem 6 of Wang
et al. (2015), we have:

$\displaystyle\|s_{j}^{\prime}\|_{2}^{2}$
$\displaystyle=\|(I_{p}-V_{[k]}V_{[k]}^{T})V\Xi^{\dagger}U^{T}e_{j}\|_{2}^{2}$
$\displaystyle=\|\begin{bmatrix}0&V_{-[k]}\end{bmatrix}\Xi^{\dagger}U^{T}e_{j}\|_{2}^{2}=\left\|\sum_{i=k+1}^{p-1}\xi_{i}^{-1}\langle
u_{i},e_{j}\rangle v_{i}\right\|_{2}^{2}$
$\displaystyle=\sum_{i=k+1}^{p-1}\xi_{i}^{-2}\langle
u_{i},e_{j}\rangle^{2}\leq\frac{2}{p}\sum_{i=k+1}^{p-1}\xi_{i}^{-2}$

where we made use of the fact that the left singular vectors
$\\{u_{i}\\}_{i=1}^{m}$ of $\Gamma$, when $\Gamma$ is the incidence matrix of
the chain graph with $p$ vertices, satisfy $\forall
i\in[m]:\|u_{i}\|_{\infty}\leq\sqrt{\frac{2}{p}}$.

For the chain graph, the nonzero singular values $\xi_{i}$ are such that

$\xi_{i}^{2}=4\sin^{2}\left(\frac{\pi i}{2p}\right)=2-2\cos\left(\frac{\pi
i}{p}\right),\text{ for }i=1,\dots,p-1$

Hence, as in Wang et al. (2015),

$\displaystyle\max_{j\in[m]}\|s_{j}^{\prime}\|_{2}^{2}\leq\frac{2}{p}\sum_{i=k+1}^{p-1}\xi_{i}^{-2}$
$\displaystyle=\frac{1}{2p}\sum_{i=k+1}^{p-1}\sin^{-2}\left(\frac{\pi
i}{2p}\right)$
$\displaystyle\leq\frac{1}{2p}\int_{k}^{p}\sin^{-2}\left(\frac{\pi
x}{2p}\right)dx=\frac{\cos\left(\frac{\pi
k}{2p}\right)}{\pi\sin\left(\frac{\pi k}{2p}\right)}\leq\frac{4p}{\pi^{2}k}$

where we used $\sin(x)\geq x/2$ and $\cos(x)\leq 1$ for $x\in[0,\pi/2]$.

Similar to Lemma A.4, we can then select
$\lambda_{1}=\frac{64}{\pi}\sigma\sqrt{p/k}\sqrt{\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}}$. We also have $\frac{2}{n}\|P_{[k]}X^{T}\epsilon\|_{2}\leq
4\sqrt{\frac{2\sigma^{2}\gamma_{\text{max}}(\Sigma)k}{n}}$ with probability at
least $1-e^{-n/8}-e^{-k^{2}/2}$, as in Lemma A.3. The rest of the proof is
again identical to that of Theorem 2.1, and we obtain for any $S$ that

$\|\Sigma^{1/2}(\hat{\beta}-\beta^{*})\|_{2}^{2}\lesssim\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)\left[\frac{1+\delta^{2}}{n}+\frac{k}{n}\right]+\lambda_{1}^{2}|S|k_{S}^{-2}}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}+\lambda_{1}\|(\Gamma\beta^{*})_{-S}\|_{1}$

with high probability. By setting $S=\emptyset$ and choosing $k$ such that

$\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)k}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)n}\asymp\lambda_{1}\|\Gamma\beta^{*}\|_{1}\asymp\sigma\sqrt{\frac{p\gamma_{\text{max}}(\Sigma)\log
p}{kn}}\|\Gamma\beta^{*}\|_{1}$

we obtain

$k\asymp\left(\frac{p\log
p\|\Gamma\beta^{*}\|_{1}^{2}\gamma_{\text{min}}^{2}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)n}{\sigma^{2}\gamma_{\text{max}}(\Sigma)}\right)^{1/3}$

and with this choice of $k$

$\frac{\sigma^{2}\gamma_{\text{max}}(\Sigma)k}{\gamma_{\text{min}}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)n}\asymp\frac{(\sigma^{2}\gamma_{\text{max}}(\Sigma)\|\Gamma\beta^{*}\|_{1})^{2/3}}{\gamma_{\text{min}}^{1/3}\left(\frac{1}{64}\Sigma+\lambda_{2}L\right)}\frac{(p\log
p)^{1/3}}{n^{2/3}}$

∎

We will often use the following lemma to compare probabilities involving two
Gaussian vectors.

######  Lemma A.1 (Anderson’s Gaussian comparison inequality Anderson (1955)).

Let $X$ and $Y$ be two zero-mean Gaussian vectors with covariance $\Sigma_{X}$
and $\Sigma_{Y}$ respectively. If $\Sigma_{Y}-\Sigma_{X}$ is positive semi-
definite, then for any convex set $C$ satisfying $C=-C$,

$\mathbb{P}(X\in C)\geq\mathbb{P}(Y\in C)$

######  Lemma A.2 (Lemma 2.1).

If $X\in\mathbb{R}^{n\times p}$ has i.i.d. $N(0,\Sigma)$ rows and $m\geq 2$,
$n\geq 10$, then the event

$\bigg{\\{}\forall
v\in\mathbb{R}^{p}:\frac{\|Xv\|_{2}}{\sqrt{n}}\geq\frac{1}{4}\|\Sigma^{1/2}v\|_{2}-3\sqrt{\frac{\gamma_{\text{max}}(\Sigma)n_{c}}{n}}\|v\|_{2}-6\sqrt{2}\rho(\Gamma)\sqrt{\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}}\|\Gamma v\|_{1}\bigg{\\}}$

holds with probability at least $1-c_{1}\exp(-nc_{2})$, for some universal
constants $c_{1},c_{2}>0$.

Proof We follow the proof outline of Raskutti et al. (2010). First note that
we can restrict our attention to $v\in\mathbb{R}^{p}$ satisfying
$\|\Sigma^{1/2}v\|_{2}=1$, as the inequality that defines the event above is
invariant to scaling of $v$. Define

$V(r,s):=\\{v\in\mathbb{R}^{p}:\|\Sigma^{1/2}v\|_{2}=1,\|\Gamma v\|_{1}\leq
r,\|v\|_{2}\leq s\\}$

and

$M(r,s,X):=\sup_{v\in V(r,s)}\left(1-\frac{\|Xv\|_{2}}{\sqrt{n}}\right)$

Bounding the expectation $\mathbb{E}(M(r,s,X))$: By an application of Gordon’s
inequality (see Section 4.2 of Raskutti et al. (2010) for the details),

$\displaystyle\mathbb{E}\left(\sup_{v\in
V(r,s)}(-\|Xv\|_{2})\right)=\mathbb{E}\left(\sup_{v\in V(r,s)}\inf_{u\in
S^{n-1}}u^{T}Xv\right)$ $\displaystyle\leq\mathbb{E}\left(\sup_{v\in
V(r,s)}\inf_{u\in S^{n-1}}g^{T}u+h^{T}\Sigma^{1/2}v\right)$
$\displaystyle=-\mathbb{E}\|g\|_{2}+\mathbb{E}\left(\sup_{v\in
V(r,s)}h^{T}\Sigma^{1/2}v\right)$

where $g\sim N(0,I_{n})$ independent of $h\sim N(0,I_{p})$. We know that
$\mathbb{E}\|g\|_{2}\geq\frac{3}{4}\sqrt{n}$ when $n\geq 10$, so we just need
to upper bound $\mathbb{E}\left(\sup_{v\in V(r,s)}h^{T}\Sigma^{1/2}v\right)$.

Since $\Pi+\Gamma^{\dagger}\Gamma=I_{p}$,

$h^{T}\Sigma^{1/2}v=h^{T}\Sigma^{1/2}(\Pi+\Gamma^{\dagger}\Gamma)v\leq\|\Pi\Sigma^{1/2}h\|_{2}\|v\|_{2}+\|(\Gamma^{\dagger})^{T}\Sigma^{1/2}h\|_{\infty}\|\Gamma
v\|_{1}$

and by definition of $V(r,s)$ we have $\|v\|_{2}\leq s$ and $\|\Gamma
v\|_{1}\leq r$ for all $v\in V(r,s)$, so we obtain

$\mathbb{E}\left(\sup_{v\in V(r,s)}h^{T}\Sigma^{1/2}v\right)\leq
s\mathbb{E}\|\Pi\Sigma^{1/2}h\|_{2}+r\mathbb{E}\|(\Gamma^{\dagger})^{T}\Sigma^{1/2}h\|_{\infty}$

Note that the spectral decomposition of $\Pi=U\Lambda U^{T}$, where $U$ is an
orthogonal matrix, is such that $\Lambda$ is a diagonal matrix with $n_{c}$
ones and $p-n_{c}$ zeros on the diagonal. Since
$\gamma_{\text{max}}(\Sigma)\Pi-\Pi\Sigma\Pi$ is positive semi-definite, by
Lemma A.1 we know that $\|\sqrt{\gamma_{\text{max}}(\Sigma)}\Pi h\|_{2}$
stochastically dominates $\|\Pi\Sigma^{1/2}h\|_{2}$, and hence

$\displaystyle\mathbb{E}\|\Pi\Sigma^{1/2}h\|_{2}$
$\displaystyle\leq\sqrt{\gamma_{\text{max}}(\Sigma)}\mathbb{E}\|\Pi h\|_{2}$
$\displaystyle=\sqrt{\gamma_{\text{max}}(\Sigma)}\mathbb{E}\|U\Lambda
U^{T}h\|_{2}$
$\displaystyle=\sqrt{\gamma_{\text{max}}(\Sigma)}\mathbb{E}\|\Lambda h\|_{2}$
$\displaystyle=\sqrt{\gamma_{\text{max}}(\Sigma)}\mathbb{E}\sqrt{h_{1}^{2}+\cdots+h_{n_{c}}^{2}}$
$\displaystyle\leq\sqrt{\gamma_{\text{max}}(\Sigma)n_{c}}$

where we have used Jensen’s inequality and the rotational invariance of the
standard Gaussian distribution in the above derivations.

By Exercise 2.12 b) of Wainwright (2019), we also have for $m\geq 2$:

$\displaystyle\mathbb{E}\|(\Gamma^{\dagger})^{T}\Sigma^{1/2}h\|_{\infty}$
$\displaystyle=\mathbb{E}\max_{j\in[m]}|\langle s_{j},\Sigma^{1/2}h\rangle|$
$\displaystyle\leq 2\sqrt{\gamma_{\text{max}}(\Sigma)}\rho(\Gamma)\sqrt{\log
m}\leq 2\sqrt{2}\sqrt{\gamma_{\text{max}}(\Sigma)}\rho(\Gamma)\sqrt{\log p}$

since {$\langle s_{j},\Sigma^{1/2}h\rangle:j=1,\cdots,m\\}$ is a collection of
$m$ zero-mean Gaussian variables with variance at most
$\gamma_{\text{max}}(\Sigma)\max_{j\in[m]}\|s_{j}\|_{2}^{2}=\gamma_{\text{max}}(\Sigma)\rho(\Gamma)^{2}$
(and in the last inequality we used $m\leq p^{2}$).

We can therefore conclude

$\mathbb{E}\left(-\inf_{v\in
V(r,s)}\|Xv\|_{2}\right)\leq-\frac{3}{4}\sqrt{n}+s\sqrt{\gamma_{\text{max}}(\Sigma)n_{c}}+2\sqrt{2}r\sqrt{\gamma_{\text{max}}(\Sigma)}\rho(\Gamma)\sqrt{\log
p}$

Dividing by $\sqrt{n}$ and adding 1 on both sides, we obtain

$\mathbb{E}\left(M(r,s,X)\right)\leq\frac{1}{4}+s\sqrt{\frac{\gamma_{\text{max}}(\Sigma)n_{c}}{n}}+2\sqrt{2}r\sqrt{\gamma_{\text{max}}(\Sigma)}\rho(\Gamma)\sqrt{\frac{\log
p}{n}}$

Concentration around the mean for $M(r,s,X)$: As $M(r,s,X)$ is a Lipschitz
function of a Gaussian vector (see Section 4.3 of Raskutti et al. (2010) for
details), for all $t>0$ we have:

$\mathbb{P}(|M(r,s,X)-\mathbb{E}M(r,s,X)|\geq t/2)\leq 2\exp(-nt^{2}/8)$

Substituting
$t=t(r,s):=\frac{1}{4}+s\sqrt{\frac{\gamma_{\text{max}}(\Sigma)n_{c}}{n}}+2\sqrt{2}r\sqrt{\gamma_{\text{max}}(\Sigma)}\rho(\Gamma)\sqrt{\frac{\log
p}{n}}$, we obtain

$\mathbb{P}\left(M(r,s,X)\geq\frac{3t(r,s)}{2}\right)\leq
2\exp(-nt(r,s)^{2}/8)$

Peeling: This part is adapted from Section 4.4 of Raskutti et al. (2010). We
have shown that

$\displaystyle\mathbb{P}\left(\sup_{\begin{subarray}{c}\|v\|_{2}\leq
s,\|\Gamma v\|_{1}\leq r\\\
\|\Sigma^{1/2}v\|_{2}=1\end{subarray}}\left(1-\frac{\|Xv\|_{2}}{\sqrt{n}}\right)\geq\frac{3}{8}+\frac{3}{2}\sqrt{\frac{\gamma_{\text{max}}(\Sigma)n_{c}}{n}}s+3\sqrt{2}\rho(\Gamma)\sqrt{\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}}r\right)$ $\displaystyle\leq
2\exp\left(-\frac{n}{18}\left(\frac{3}{8}+\frac{3}{2}s\sqrt{\frac{\gamma_{\text{max}}(\Sigma)n_{c}}{n}}+3\sqrt{2}\rho(\Gamma)r\sqrt{\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}}\right)^{2}\right)$

Let
$g_{1}(s):=\frac{3}{16}+\frac{3}{2}\sqrt{\frac{\gamma_{\text{max}}(\Sigma)n_{c}}{n}}s$
and
$g_{2}(r):=\frac{3}{16}+3\sqrt{2}\rho(\Gamma)\sqrt{\frac{\gamma_{\text{max}}(\Sigma)\log
p}{n}}r$. We can rewrite the above as

$\mathbb{P}\left(\sup_{\begin{subarray}{c}\|v\|_{2}\leq s,\|\Gamma v\|_{1}\leq
r\\\
\|\Sigma^{1/2}v\|_{2}=1\end{subarray}}\left(1-\frac{\|Xv\|_{2}}{\sqrt{n}}\right)\geq
g_{1}(s)+g_{2}(r)\right)\leq
2\exp\left(-\frac{n}{18}[g_{1}(s)+g_{2}(r)]^{2}\right)$

Note that $g_{1}\geq\mu$ and $g_{2}\geq\mu$ where $\mu:=\frac{3}{16}$. For
$i=1,2,\cdots,$ and $j=1,2,\cdots$, we define the sets

$A_{ij}:=\\{v\in\mathbb{R}^{p}:\|\Sigma^{1/2}v\|_{2}=1,2^{i-1}\mu\leq
g_{1}(\|v\|_{2})<2^{i}\mu,2^{j-1}\mu\leq g_{2}(\|\Gamma v\|_{1})<2^{j}\mu\\}$

Also, we define the events