The Function-on-Scalar LASSO
with Applications to Longitudinal GWAS

arXiv:1610.07403v1 [math.ST] 24 Oct 2016

Rina Foygel Barber1 , Matthew Reimherr

∗2

, and Thomas Schill

2

1

2

University of Chicago
Pennsylvania State University

Abstract
We present a new methodology for simultaneous variable selection and parameter
estimation in function-on-scalar regression with an ultra-high dimensional predictor
vector. We extend the LASSO to functional data in both the dense functional setting
and the sparse functional setting. We provide theoretical guarantees which allow for an
exponential number of predictor variables. Simulations are carried out which illustrate
the methodology and compare the sparse/functional methods. Using the Framingham
Heart Study, we demonstrate how our tools can be used in genome-wide association
studies, finding a number of genetic mutations which affect blood pressure and are
therefore important for cardiovascular health.

K eywords: Functional Data Analysis, High-Dimensional Regression, Variable Selection,
Functional Regression

1

Introduction

Over the last several decades, technological advances have supported and necessitated the
rapid growth of high-dimensional statistical methods. One the most important applications
for these methods are genome-wide association studies (GWAS). In these studies, researchers
search through hundreds of thousands or millions of genetic mutations known as SNPs, single
nucleotide polymorphisms, finding those which significantly impact an outcome or phenotype
of interest, e.g. blood pressure, diabetes, asthma, etc. GWAS have been hugely successful at
finding gene/disease associations as evidenced by the massive repositories of genetic studies
∗

Corresponding Author: Department of Statistics, Pennsylvania State University, University Park, PA,
16802

1

and findings, such as dbGaP (http://www.ncbi.nlm.nih.gov/gap). Next-generation or highthroughput sequencing technologies are capable of sequencing entire genomes, producing
ever-larger genetic datasets to explore. Thus, to further our understanding of the genetic
architecture of complex human diseases, there is a substantial and continuing need for powerful, high-dimensional association techniques.
The vast majority of GWAS are cross-sectional, examining only one point in time. The
aim of this paper is to present a new framework which combines functional data analysis,
FDA, and machine leaning for finding and estimating significant effects on longitudinally
measured outcomes. We refer to this methodology as Function-on-Scalar LASSO, or FSLASSO. As in Reiss et al. (2010), we use the term Function-on-Scalar to help distinguish
our setting from the wide array of regression problems now in existence. In GWAS, the
“function” is the time-varying phenotype of the individual, and it is being regressed on the
“scalar”, which is the individual’s constant (non-time-varying) genotype.
Our aim is to simultaneously exploit the sparse effect of the SNPs and the smooth nature
of the longitudinal outcomes. While we are strongly motivated by genetic studies, our
methods are general and allow for any setting with a longitudinal/functional outcome and
high-dimensional scalar predictors. Our primary goal is to select and estimate the effect of
predictors in the following functional linear model
(1)

Yn (t) = µ(t) +

I
X

Xni βi (t) + εn (t).

i=1

Here Yn (t) is the value of a quantitative outcome for subject n ∈ {1, . . . , N }, at time t ∈ T ⊂
R. The scalars Xni are real-valued, though in our application they take values in {0, 1, 2}
indicating the minor allele count for the SNP. The number of predictors, I, is allowed to
be much greater than the sample size N . This is known as a “scalar predictor/functional
response” model. Functional data analysis (FDA) now consists of two main branches: (1)
sparse FDA, where the outcomes, Yn , are observed at a relatively small number of time
points and are contaminated with noise, and (2) dense FDA, where outcomes are observed
at a large number of time points (possibly with a small level of noise). In practice, some data
clearly fall into one of these two categories, say with hundreds or thousands of observations
per subject, or with only a handful. There are also many settings where it is not clear which
scenario one is in. While establishing a clear cutoff still a relatively open problem, we refer
the interested reader to Li and Hsing (2010) and Zhang and Wang (2015+) who suggest the
line is when the number of points per curve is greater/less than the quad-root of the sample
size. We emphasize that this is just a basic rule of thumb, and care must be taken with each
application. We present our methodology for each scenario and we later explore how their
performances differ via simulations. In the functional case, as we will see, our methods and
theory will actually include a wide range of settings beyond just (1).
At the heart of our methodology is the now classic idea of combining basis expansions for
estimating the functions βi , with proper penalties which induce sparse estimates. While the
2

methods and theory are new, we can phrase the problem as a type of group LASSO, which
allows us to utilize existing computational tools.
The contributions of this paper include the following. First, we provide new asymptotic
results for the sparse setting. Our results can be viewed as an extension of the varying
coefficient methods discussed in Wei et al. (2011). In particular, we present a new restricted
eigenvalue result which is one of the cornerstones for developing convergence rates in high
dimensions, as well as the accompanying asymptotic convergence rates for our estimates.
Second, we provide a new methodology and accompanying asymptotic theory for a broad
class of dense settings, which, to the best of our knowledge, has not been explored before.
Interestingly, we show that the FS-LASSO applied to outcomes in any separable Hilbert
space achieves the exact same rates of convergence as in the scalar setting. This could have
broad implications for not only traditional functional outcomes, but also spatial processes,
functional panels, and imaging data such as fMRI. Finally, we demonstrate via simulations
some surprising results concerning the choice between the sparse and dense tool sets. In
particular, we find that a dense approach is comparable to the sparse in terms of variable selection even in traditionally sparse settings. Furthermore, the dense methods can be carried
out at a fraction of the computational cost (both in terms of power and memory). However,
the sparse methods produce estimates which are more accurate in traditionally sparse settings. This opens the doors to interesting two-stage procedures where variable screening is
done using dense tools, and final estimation is done via sparse ones.
Related literature For foundations on functional data analysis we refer to Ramsay and
Silverman (2005) and Horváth and Kokoszka (2012), while an overview of machine learning
can be found in Hastie et al. (2001) and James et al. (2013). The literature on functional
regression is now quite large, but we attempt to provide several key methodological papers,
which help outline the field. Functional regression methods for sparse data include the following. Hoover et al. (1998) examine spline based methods for estimating a simpler form
of the functional predictor/functional response model. Fan and Zhang (2000) provide a two
step method based on local polynomial smoothing to estimate a functional predictor/scalar
response model. Yao et al. (2005) present what has likely become the most common method
based on scatter plot smoothing (i.e. local polynomial smoothing) and functional principal
components for estimating a full functional predictor/functional response model. Zhu et al.
(2012) extend local linear smoothing methods for estimating scalar predictor/multivariate
functional response models. Methods for high dense FDA include the following. Cardot
et al. (2003) explore spline based methods for a functional predictor/scalar response model.
Kokoszka et al. (2008) use a PCA based approach for estimating a full functional predictor/functional response model. James et al. (2009) incorporate shrinkage methods in estimating a functional predictor/scalar response model (for the purposes of estimating the
”zero” parts of the regression function). Reiss et al. (2010) examine a B-spline approach
for scalar predictor/functional response models, while Reimherr and Nicolae (2014) explore
more direct least squares (of a functional norm) estimators for scalar predictor/functional re3

sponse models. However, all of these methods assume a fixed number of covariates. Related
work which incorporates an increasing number of predictors mainly comes from the literature on varying coefficient models, for example Wang et al. (2008), Zhao and Xue (2010) and
Wei et al. (2011). The only related work in the FDA literature concerns scalar-on-function
regression, the reverse of our focus. In particular, Matsui and Konishi (2011) combine basis
expansions with a SCAD style penalty, but give no theoretical guarantees. Lian (2013) also
uses a SCAD penalty, but combined with FCPA instead of a general basis function. Theory
is provided, but only in the case where the number of predictors is fixed. Gertheiss et al.
(2013) developed a procedure for a generalized scalar-on-function linear model again using
basis expansions, but combined with a LASSO style penalty which simultaneously shrinks
and smooths parameter estimates; no theoretical guarantees are given. Fan et al. (2014,
2015+) consider functional additive regression models with a large number of predictors.
Their methods allow for nonlinear relationships and theoretical properties are established.
Lastly, Kong et al. (2016) consider the problem of variable selection with a scalar outcome
and scalar predictors, but with a small set of functional predictors that must be accounted
for as well. We emphasize that all of these methods are for the scalar-on-function case; for
the function-on-scalar case the only work we are aware of is by Chen et al. (2016) who consider a basis expansion approach with a MCP style penalty and fixed number of covariates,
but their method cannot be applied to settings where I  N .
Outline The remainder of the paper is organized as follows. In Section 2 we provide
a framework for scalar predictor/functional response regression with sparsely observed outcomes. We combine basis expansions with a group LASSO style penalty to carry out variable
selection and parameter estimation simultaneously. In Section 3 we present the high dense
setting. We provide a very general framework where the response functions are allowed to
take values from a separable Hilbert space. This allows the functions to be from say L2 (T ),
which is commonly used, as well as more general spaces such as Sobelev spaces or product
spaces of function spaces (i.e. vectors of functions). In Section 4 we provide computational
details, including a pre-screening rule which allows one to reduce I substantially before fitting
the FS-LASSO, and details for making use of established computational machinery for the
group LASSO. We present a simulation study in Section 5 where we explore our procedure
in terms of variable selection, parameter estimation, and computational time. In Section 6
we apply our methods to the Framingham Heart Study (FHS), to identify genetic markers
that influence blood pressure, in hopes of gaining a better understanding of cardiovascular
disease. Genome–wide SNP data and phenotype information were downloaded from dbGaP
(http://www.ncbi.nlm.nih.gov/gap) study accession phs000007.v25.p9. Concluding remarks
are given in Section 7, while all theoretical results are proven in the Appendix.

4

2

The Sparse Setting

Sparse FDA occurs quite often in longitudinal studies where one has a relatively small number
of observed time points per subject or unit. In this case, it is common to work with the
raw measurements as opposed to the functional embeddings discussed in Section 3. A brief
introduction to sparse FDA can be found in Müller (2008).
We begin with the underlying model.
Assumption 1. Assume, for 1 ≤ n ≤ N and 1 ≤ m ≤ Mn ≤ M < ∞, that
Ynm =

I
X

Xni βi? (tnm ) + εn (tnm ),

tnm ∈ T .

i=1

The error processes {εn (t) : t ∈ T } are iid Gaussian and we define the Mn × Mn covariance

matrix Σn := Cov (εn (tn1 ), . . . , εn (tnMn )) . The design matrix, X, is deterministic.
Throughout, we will use a ? to denote true parameter values. For the moment, we do not
include assumptions on {tnm } or coefficient functions βi? . Note that the Gaussian assumption
is not crucial—what is really needed is that the processes have subgaussian tails, though we
make the Gaussian assumption to simplify the arguments.
Let {ej (·)} be a basis in L2 (T ). Approximating the functions βi? using this basis, we have
that
I X
J
X
Ynm =
Xni Bij? ej (tnm ) + Tnm + εn (tnm ),
i=1 j=1

where Tnm is the truncation error obtained after cutting off the basis expansion at J. The
following notation will be used repeatedly throughout this section and in the Appendix:
P

T = (T11 , T12 , . . . , TN MN )> ∈ R

n

Mn

,

Enm = (e1 (tnm ), . . . , eJ (tnm ))> ∈ RJ ,
En = (En1 , . . . EnMn )> ∈ RMn ×J ,
F=

N
1 X >
E En ∈ RJ×J .
N n=1 n

For a matrix B ∈ RI×J , define the sparse target function as
!2
N Mn
I X
J
X
1XX
L(B) =
Ynm −
Xni Bij ej (tnm ) + λ kBk`1 /`2 ,
2 n=1 m=1
i=1 j=1
where kBk`1 /`2 =

P

kBi∗ k2 promotes row-wise sparsity (here Bi∗ denotes the ith row of
b is the minimizer of the above expression. The estimated
the matrix B). The estimate B
i

5

P
bij ej (·). The target function can be
coefficient functions are then given by βbi (·) = Jj=1 B
rephrased so that traditional group LASSO machinery can be invoked. Notice that
I X
J
X

>
>
Xni Bij ej (tnm ) = Xn> BEnm = Enm
B> Xn = (Xn> ⊗ Enm
) vec(B> ).

i=1 j=1
>
One can then stack the (Xn> ⊗ Enm
) vectors into a matrix, A,
  >
 >

>
X1 ⊗ E11
X1 ⊗ E 1
 X > ⊗ E >   X > ⊗ E2 
P
12 
 1
 2

( n Mn )×IJ
A=
=
∈
R
.



..
..
 


.
.
>
XN> ⊗ EN
MN

XN> ⊗ EN

The target function can now be expressed as
1
L(B) = kY − A vec(B> )k22 + λkBk`1 /`2 ,
2

(2)
where

Y = (Y11 , . . . , Y1M1 , . . . , YN,MN )> .
b
Group LASSO computational tools can then be used to find B.
Since the matrix A will generally have more columns than rows, the linear system Y ≈
A vec(B> ) is underdetermined; in particular the (right) null space of A is large. However,
the grouped sparsity structure in B allows us to resolve this difficulty. We first recall the
notion of restricted eigenvalues Bickel et al. (2009), used for sparse regression (without group
structure):
Definition 2.1 (Restricted eigenvalue condition). A matrix A ∈ RN ×I satisfies the RE(I0 , α)
condition if, for all subsets S ⊂ {1, . . . , I} with |S| ≤ I0 ,
kAwk22 ≥ αN kwk22 for all w ∈ RI with kwS c k`1 ≤ 3 kwS k`1 .
The constant 3 here is somewhat arbitrary, but it is standard in the literature and we use it
here for convenience. Lounici et al. (2011) extend this definition to the group-sparse setting:
Definition 2.2 (Grouped restricted eigenvalue condition). A matrix A ∈ RN ×IJ satisfies
the REgroup (I0 , α) condition if, for all subsets S ⊂ {1, . . . , I} with |S| ≤ I0 ,
A vec(W> )

2
2

≥ αN kWk2F for all W ∈ RI×J with kWS c k`1 /`2 ≤ 3 kWS k`1 /`2 .

Here WS denotes the submatrix of W obtained by extracting the rows indexed by S, while
WS c contains only the rows in S c , the complement of S.

6

2.1

Bounding the error

With these definitions in place, our first result is a modified version of Lounici et al. (2011)’s
b − B? , where B? is
analysis of the group LASSO. It proves resulting bounds on the error B
the (truncated) true parameter matrix, as long as λ is sufficiently large and A satisfies the
grouped restricted eigenvalue condition.
Theorem 1. Suppose that Assumption 1 holds and that for some δ ∈ (0, 1)


q
q
(3)
λ ≥ 2 N kFkop kT k2 + max kΣn kop (2J + 3 log(2I/δ)) .
n

Assume also that B? has at most I0 nonzero rows, and that A satisfies the REgroup (I0 , α)
b of (2)
condition for some α > 0. Then with probability at least 1 − δ, any minimizer B
satisfies
√
3λ
I0
12λI0
?
b −B
b − B?
B
≤
≤
and B
.
αN
αN
F
`1 /`2
Remark 2.1. Results for sparse regression in a non-grouped setting commonly give error
bounds in both an `2 norm and an `1 norm; the former typically gives a more favorable scaling
of sample size with respect to sparsity, while the latter ensures that for a larger sample size
the error itself is nearly sparse, i.e. does not have a “long tail” of small errors on many
coordinates. Similarly, here we give results in both the Frobenius norm (with more favorable
sample size scaling) and the `1 /`2 norm.
Remark 2.2. As is standard in the LASSO literature, the consistency results given in this
theorem can yield a selection consistency result as well: if we assume that mini∈I0 kB?i k2 >
√
6λ I0
b − B? implies that the following estimated support,
, then the bound on B
αN

F


b i∗
b
S = i ∈ {1, . . . , I} : B

2

√ 
3λ I0
>
,
αN

is correct (i.e. Sb is equal to the row support of B? ) with high probability.
As is widely appreciated by researchers working on such results, Theorem 1 essentially
hinges on concentration inequalities and restricted eigenvalue type conditions. Neither of
these tools can be readily applied to our setting. Thus, much of our theoretical work is
focussed on showing how to extend these ideas to the functional setting. Below we discuss the
restricted eigenvalue condition at length, while discussion of the concentration inequalities
involved can be found at the end of Appendix B.1. Concentration inequalities allow us to
control the stochastic error in the model and largely dictate the rates of convergence, while
b − B? from above by
the restricted eigenvalue conditions allow us to bound the error B
F

b − B? )> k2 , and therefore to obtain convergence rates for the
a term involving kA vec (B
parameters of interest.
7

Note that the phrasing of Theorem 1 allows for the possibility of multiple minimizers of (2),
guaranteeing only that any minimizer must achieve the specified converge rates. However,
in practice it rare is for LASSO style estimates to not be unique. We direct the interested
reader to Tibshirani (2013) for a thorough discussion of this issue.
Example: We next give an example to illustrate the type of scaling that is obtained in
Theorem 1. Let W τ,2 [0, 1] denote the Sobelev space consisting of functions over [0, 1], which
have square integrable derivatives of up to order τ . Suppose that the coefficient functions,
βi , take values in a finite ball in W τ,2 . Let the basis e1 (·), e2 (·), . . . denote the Fourier basis
and assume that the time points tnm are iid U [0, 1]. Examining λ, there are now a number
to terms which can be made more explicit. First, since the basis is orthogonal we have
F ≈ IJ×J , so that kFkop is behaves like a constant. Next, examining the truncation term
P
T , there are n Mn ∼ N coordinates, and each coordinate consists of a sum of I0 functions
√
each of which lies in a Sobelev ball, thus we have that kT k ≈ N I0 J −τ . Examining the
second term in λ, we have that
q
p
p
max kΣn kop (2J + 3 log(2I/δ)) ∼ J + log(I) ≤ J log(I).
n

We therefore want to choose J to balance the two errors, which means taking J such that
1



p
p
N I0 1+2τ
−τ
N I0 J = J log(I) =⇒ J =
.
log(I)
Plugging this into our expression for the convergence rates, and assuming that the restricted
eigenvalue condition REgroup (I0 , α) holds for the matrix A with some α > 0 that we treat as
a constant, we have that
p

 τ
1/2
N 1/2 J log(I)I0
log(I) 1+2τ
= OP (1) · I0
,
k∆kF = OP (1) ·
N
N I0
b We see that the standard nonparametric rate of convergence appears
where ∆ := B − B.
which relates the convergence rates to the smoothness of the underlying parameter functions.
This rate also applies to the difference kβb − βkL2 . Notice that by Parceval’s identity and the
Triangle inequality
kβb − βkL2 =

I Z
X

!1/2
≤ k∆kF + kΠ⊥
J βkL2 ,

(βbi (t) − βi (t)2 ) dt

i=1

where Π⊥
J is the projection onto the remaining basis functions eJ+1 (·), . . . . By the assumed
smoothness of the β, we have that
−τ

kΠ⊥
J βkL2

τ



p
p  N I0  1+2τ
log(I) 1+2τ
−τ
∼ I0 J = I0
≤ I0
,
log(I)
N I0
8

Thus the rate for kβb − βkL2 is the same as for k∆kF .
For the `1 /`2 norm we then have
k∆k`1 /`2

λI0
3/2
= OP (1) ·
= OP (1) · I0
N



log(I)
N I0

τ
 1+2τ

.

By the same arguments as before, the above rate also applies to the sum of the normed
differences between the functions, i.e. kβb − βk`1 . As we will see, when τ → ∞, the rate
above converges to the rate found for the dense setting.

2.2

Restricted eigenvalue condition

Our next result concerns the restricted eigenvalue condition. It proves that, if the covariates
Xn are drawn from a Gaussian or subgaussian distribution with well-conditioned covariance
structure, then the grouped restricted eigenvalue condition will hold for A.
iid

Theorem 2. Suppose that Xn ∼ N (0, Σ) with kΣk ≤ ν 2 , or more generally, the Xn ’s are
>
iid where each Xn has mean 0, covariance Σ, and is ν-subgaussian meaning that E[ev Xn ] ≤
2
2
eν kvk2 /2 for any fixed vector v. Suppose that λmin (Σ) > 0. Consider any fixed sequence of
matrices E1 ∈ RM1 ×J , . . . , EN ∈ RMN ×J satisfying
sP
P
2
kE
wk
n
n
2
n kEn wk2
(4)
min
≥ γ0 > 0,
max
≤ γ1 ,
N kwk2
w∈RJ \{0}
w∈RJ \{0}
N kwk22
P

and define A ∈ R n Mn ×IJ as before. Then there exist c0 , c1 , c2 > 0 depending only on
ν, λmin (Σ), γ0 , γ1 , such that, if
(5)

N ≥ c0 · k · log(IJ) · (J + log(I)) ,

then with probability at least 1 − e−c1 N , A satisfies the REgroup (k, c2 ) condition.
Remark 2.3. This theorem is a primary distinction between our work discussed thus far
and that of Wei et al. (2011). Their work cites existing results which imply that the matrix A behaves appropriately for each approximately-row-sparse W in expectation, that is,

2
E A vec(W> ) 2 is lower-bounded for each W. In contrast, our theorem above proves the
2
much stronger statement that A vec(W> ) 2 can be lower-bounded simultaneously with high
probability for all approximately-row-sparse W, which is necessary for estimation in high
dimensions.
Remark 2.4. The conditions (4) on the basis matrices En are of course dependent on the
choice of basis; this type of basis-dependent assumption is standard in the nonparametric
literature. As a simple example of a choice of basis where (4) is satisfied, suppose that there
are J evenly spaced time points and our basis is given by J orthonormal functions (e.g. a
9

Fourier basis). Suppose measurements for each individual are taken at a random subset of
J 0 time points from the original J. Then each En is a J 0 × J submatrix chosen at random
from some larger J × Jporthogonal matrix (corresponding to e.g. the Fourier basis), in which
case γ1 = 1 and γ0 ∼ J 0 /J.
The assumption given in (4) is a condition number type constraint. The upper bound is
P
simply a bound on the largest eigenvalue of F = N1 n E>
n En , while the lower bound is a bit
stronger than a corresponding bound on the smallest eigenvalue. It is, in essence, saying that
the norms kEn wk22 should be on the order of kwk22 for many of the indices n. (In contrast,
upper and lower eigenvalue conditions would only ensure that kEn wk22 is on the order of
kwk22 on average, but would not prevent degenerate scenarios such as kE1 wk22 = N kwk22 and
kE2 wk22 = · · · = kEN wk22 = 0.)

3

The Dense Setting

When the underlying functions are observed at a relatively large number of time points, a
different approach than the one described in Section 2 is commonly employed. In particular,
for each subject n ∈ {1, . . . , N }, the observations Ynm = Yn (tnm ) are embedded into a function space, which are then treated as though they were fully observed functions. Implicitly,
one is assuming that the error from the embedding is negligible compared to other sources
of variability. More details and background can be found in Horváth and Kokoszka (2012).
We begin by defining the underlying model.
Assumption 2. Let Y1 , . . . , YN be independent random elements of a real separable Hilbert
space, H, satisfying the functional linear model
Yn =

I
X

Xni βi + εn .

i=1

Assume the N × I design matrix X = {Xni } is deterministic and has standardized columns,
the εn are iid square-integrable Gaussian random elements of H with mean 0 and covariance
operator C, and βi are deterministic elements of H. Let Λ = (λ1 , λ2 , . . . ) denote the vector
of eigenvalues of C.
Remark 3.1. We emphasize that, throughout this section, the responses Yn , coefficients βi ,
and errors εn are all elements of the Hilbert space (e.g. functions), rather than scalars as in
the usual high-dimensional regression setting. Our work in this section can be viewed as a
functional analogue of existing high-dimensional sparse signal recovery results for the scalar
setting (i.e. the ordinary LASSO).

10

The most common choice for the function space H is L2 (T ) where T is a closed and
bounded interval. However, the advantage of our phrasing is that we allow for a number of
other spaces, including Sobolev spaces (if one wants to better utilize the smoothness of the
data), product spaces (if the response is actually a vector of functions), or multidimensional
domains (which could arise in areas such as spatial statistics). Throughout this section,
when we write k · kH we will always mean the inner product norm on H. Norms which
deviate from this inner product norms will include alternative subscripts.
As before, when applying the FS-LASSO, we are assuming that the true underlying model
is sparse, with I0 denoting the number of true predictors and S0 denoting the true set
of predictors. The FS-LASSO estimate is then the solution to the following minimization
problem
βb = argminβ∈HI LF (β),
where

N

1X
LF (β) =
kYn − Xn> βk2H + λkβk`1 /H
2 n=1

for kβk`1 /H =

I
X

kβi kH .

i=1

The norm k · k`1 /H is a type of `1 norm on the product space HI , which encourages sparsity
among the list of functions β1 , . . . , βI (that is, the function βi will be uniformly zero for many
indices i). The target function LF (·) is a direct Hilbert space generalization of the LASSO
target function.
Next we introduce a functional restricted eigenvalue assumption, which is a direct analogue
of similar assumptions in the LASSO literature (e.g. Bickel et al. (2009)).
Definition 3.1. We say a matrix A ∈ RN ×I satisfies a functional restricted eigenvalue
condition, REF (I0 , α), if for all subsets S ⊂ {1, . . . , I} with |S| ≤ I0 , we have
kAxk2HN ≥ αN kxk2HI for all x ∈ HI that satisfy kxS c k`1 /H ≤ 3kxS k`1 /H .
In fact, this functional restricted eigenvalue assumption is no stronger than the usual
(scalar) restricted eigenvalue assumption—in the following theorem we show that any matrix
A satisfying the usual (scalar) restricted eigenvalue assumption, will also satisfy the functional version given in Definition 3.1. We do this by showing that the common inequality
used for proving that a matrix satisfies a scalar restricted eigenvalue condition immediately
implies the same inequality in the functional setting, and leads to the restricted eigenvalue
condition.
Theorem 3. For a fixed matrix X ∈ RN ×I , suppose that for some c1 , c2 > 0, X satisfies
p
√
kXzk2 ≥ c1 N · kzk2 − c2 log(I) · kzk`1 ,
for all z ∈ RI . Then the same inequality holds for all x ∈ HI , that is,
p
√
(6)
kXxkH ≥ c1 N · kxkH − c2 log(I) · kxk`1 /H .
11

Furthermore, if c1 > 4c2

q

I0 log(I)
,
N

then X satisfies the REF (I0 , α) property with
r

α=

c1 − 4c2

I0 log(I)
N

!2
.

We are now ready to present the main result of the section.
Theorem 4. If Assumption 2 holds, X satifies REF (I0 , α), and
p
√ q
λ ≥ 2 N kΛk1 + 2kΛk2 log(I/δ) + 2kΛk∞ (log(I/δ)),
then with probability at least 1 − δ, any minimizer βb of LF (β) satifies
√
16λI0
4λ
I0
?
and
kβb − β ? k`1 /H ≤
.
kX(βb − β )kH ≤ √
αN
αN
We mention that since Λ is vector of the eigenvalues of C in decreasing order, then one
P∞
P
2
λ
,
kΛ
k
=
has that kΛk22 = ∞
1
1
i
i=1 λi , and kΛk∞ = λ1 . The target function LF (β) is
i=1
convex and coercive, and thus has a least one solution. As in the scalar case, the function may
have multiple solutions, but this still depends heavily on the design X and does not happen
often in practice. In the next section we will discuss how this estimate can be computed.
Interestingly, the rates of convergence for H-valued response variables are exactly the
same as those for a scalar response. This is due to our ability to extend p
scalar concentration
inequalities to general Hilbert spaces, see Lemma
9.
If
we
take
λ
∼
N log(I) then the
p
?
prediction error becomes
kX(βb − β )kH ∼ I0 log(I) while the estimation error becomes
q
kβb − β ? k` /H ∼ I0 log(I) . Since this result applies to any separable Hilbert space, our
1

N

methodology is applicable to a wide range of applications beyond just the GWAS we consider
here. From random fields in spatial statistics to brain imaging in fMRI studies, all fall under
this umbrella if they are embedded into a Hilbert space. As a final note, we mention that the
λ given in Theorem 4 gives a very tight control of the probability that the two inequalities
hold, however, for the convergence rates all that matters is that λ has the right order, as is
typical in LASSO style results.

4

Computational Details

Here we present several computational tools which we utilize and can be found in Matlab code
available through the corresponding authors’s website. We begin by providing computational
details for the methods in Section 3. Traditionally, when handling dense functional data,
one constructs the functional objects by utilizing basis expansions. This is the cornerstone

12

of the FDA package in R. If e1 (·), e2 (·), . . . is a basis of H (e.g. B-splines or Fourier), then
we can approximate LF as
v
N X
J
I uX
X
X
u J
1
2
>
t hβi , ej i2 ,
LF (β) ≈
hYn − Xn β, ej i +
2 n=1 j=1
i=1
j=1
for J large (often over one hundred). Letting Bi,j = hβi , ej i and Ynj = hYn , ej i, we have that
the above can be expressed as
N

J

I

X
1 XX
(Ynj − Xn> B∗j )2 +
kBi∗ k
2 n=1 j=1
i=1
I

X
1
kBi∗ k,
= kY − AF vec(B> )k2 +
2
i=1
where AF = X ⊗ IJ×J . For large values of J, this can quickly become a substantial computational burden. However, one can use a data driven basis, such as FPCA, so that J can
be taken relatively small. However, we stress that dimension reduction is not our intent,
and thus we can choose the number of FPCs to explain nearly all of the variability of the
processes. Using such an approach, it is common to move from 200 B-spline basis functions
down to 5-10 functional principal components, with nearly no information loss. This phrasing now allows us to use the same group LASSO computational tools as in the sparse setting.
However, it is now possible (even in sparse data settings) that J is less than the number
observed time points, resulting in AF being smaller than A and with a far simpler form.
Since I, the number of groups (SNPs), will generally be extremely large, we implement a
screening rule that allows us to substantially reduce the potential number of groups (SNPs)
we consider when we minimize L(·) or LF (·) for some fixed λ. Using Wang et al. (2013,
Theorem 4), we do the following:
P

• (In parallel.) For each SNP i, define Ai ∈ R n Mn ×J when using the sparse algorithm
with entries
(Ai )nm,j = Xni ej (tnm ) .
For the dense algorithm define Ai ∈ RN J×J as Ai = X (i) ⊗ IJ×J , and compute
A>
i Y
• Find any λ0 ≥ maxi A>
i Y
bλ0 = 0.
will be B

2

2

and kAi kF .

. At penalty parameter λ0 , the group LASSO solution

• For any λ < λ0 , according to Wang et al. (2013),


1
1
1
>
Ai Y 2 +
−
kAi kF kY k2 < 1
λ0
λ λ0
13

=⇒

bλ )i∗ = 0 .
(B

Therefore, if our aim is to apply a convex optimizer to solve the group LASSO problem
with at most s SNPs, then after finding λ0 = maxi A>
i Y 2 we can choose any λ sufficiently
large so that no more than s SNPs violate the inequality above. SNPs which do violate
the above are then dropped since they will not enter the solution path. This allows us to
make a substantial reduction in the number of predictors. In our simulations and application
we could handle on the order of 10000 SNPs jointly when fitting the FS-LASSO using the
ADMM procedure, with the sparse tools requiring around 30gb of RAM, and the smooth
requiring around 4gb. By using this screening rule, we can reduce millions of predictors to
tens of thousands, and then do a final fit using any number of convex optimization routines.
Finally, we mention the choice of the smoothing parameter λ. In our simulations we utilize
the BIC, though we include several other options in the application section. We calculate
the BIC as
X
log(σ̂ 2 )
Mn + JIactive log(N )
or
log(σ̂ 2 )N Mpc + Mpc Iactive log(N ),
for the sparse and dense methods respectively. Here Iactive is the number of predictors in the
current model. The error σ̂ 2 is calculated as
1 XX
1 XX
P
(Yn (tnm ) − Ŷn (tnm ))2
or
(Ynj − Ŷnj )2 ,
Mn n m
N Mpc n j
for the sparse and dense methods respectively. The predicted values are computed by recomputing the corresponding least squares estimates, to eliminate the effect of the bias. Given
the dependence in the data, the BIC is an adequate, though not optimal choice. Further
work is needed on tuning parameter selection for functional models, but we leave this for
future research. Fitting the model for a grid of values for λ can be done relatively quickly
by utilizing a warm start, i.e. using the previous solution for a particular λ value as the
starting point for finding the next solution.

5

Simulations

In this section we present a simulation study to compare the performance of the discussed
methods. The predictors, Xni , are generated from a normally distribution with Cov(Xni , Xn0 j ) =
1n=n0 ρ|i−j| , which is an autoregressive covariance with the rate of decay controlled by ρ (with
independence across subjects, i.e. Xn∗ and Xn0 ∗ are independent for n 6= n0 ). We mention
that in our genetic application the predictors take values 0/1/2 and are thus not normal.
Additional simulations using Binomial predictors can be found in Appendix A, and while the
performance of all of the methods presented here decrease, the relative conclusions from comparing the methods stay the same, we thus focus on the normal case here. We take I = 1000
and I0 = 10. The nonzero functions {βi (t) : t ∈ [0, 1]} are randomly generated from the
Matérn process with parameters (0, I0−1 , 1/4, 0, 5/2) and the errors {εn (t) : t ∈ [0, 1]} are generated in the same way but with parameters (0, 1, 1/4, 0, 3/2), which results in errors that
14

are less smooth than the parameter functions. The Matérn process is very flexible family of
stationary processes which produce more realistic structures for biological applications, as
compared to Brownian motion or simple low dimensional structures. Each subject is observed
at 10 uniformly distributed locations which differ by curve. We consider N = 50, 100, 200
and ρ = 0.5 and 0.75, and use 1000 repetitions of each scenario. For the functional method,
FPCA is carried out using the PACE package in Matlab. For the sparse method we use
cubic Bsplines with J = 30.
To compare the methods without having to worry about tuning parameter selection, we
examine smoothed ROC curves which give the proportion of true positives found as a function
of the false positives. The curves for all scenarios are given in Figure 1. Surprisingly, the
ROC curves for the two method are about the same, meaning that in terms of variable
selection the sparse and functional methods have nearly equivalent performance on average,
though, as we will see in the application section, for a single iteration they can still disagree
on the selected subset. Next we examine the average prediction error in Figure 2, with the
penalty parameters chosen by BIC, AIC, and 2-fold cross validation. There we see that the
sparse method has an advantage, though this decreases for large sample sizes. The BIC and
AIC perform about the same, while the CV criteria seems to perform the worst. We note
that, since the prediction error compares linear combinations of the parameter estimates, we
use it as a single number proxy summarizing the estimation error as well, though clearly if
there were particular patterns for β(t) one was interested in examining (linear, sinusoidal,
etc.), it would be interesting to include them in the simulations and examine their estimation
error, however we don’t purse this further here.
ROC curves, ρ = 0.5

1.25

ROC curves, ρ = 0.75

1.25

0.75

Algorithm
Smoothing
Sparse

0.50

True positive rate

1.00

True positive rate

1.00

0.75

Smoothing
Sparse

0.50

N
50

0.25

Algorithm

N
50

0.25

100

100

200

200

0.00

0.00
0.000

0.025

0.050

0.075

0.100

0.000

False positive rate

0.025

0.050

0.075

0.100

False positive rate

Figure 1: ROC curves comparing the smooth and sparse methods.
Next we compare to two forms of scalar LASSO, one which “time-corrects” by incorporating a nonparametric estimate of the mean (TC-LASSO) and one which does not. Using the
nonparametric smoother is an attempt to account for the non-stationarity of the processes.
The resulting ROC curves are given in Figure 3 for ρ = 0.5 and ρ = 0.75, with the solid
lines indicating scalar LASSO and the dashed lines the time-corrected version. Examining
the plots we see a potentially interesting pattern. For smaller false positive rates, LASSO
is actually doing better than FS-LASSO (Figure 1). However, as one moves along the x15

Mean Squared Prediction Error

Smooth, ρ=.5
Smooth, ρ=.75
Sparse, ρ=.5
Sparse, ρ=.75

0.15

0.4

Smooth, ρ=.5
Smooth, ρ=.75
Sparse, ρ=.5
Sparse, ρ=.75

Computation Time [hours]

0.20

0.3

0.10

0.2

0.05

0.1

0.00
50

100

150

200

50

N

100

150

200

N

Mean Squared Prediction Error

Smooth, ρ=.5
Smooth, ρ=.75
Sparse, ρ=.5
Sparse, ρ=.75

0.15

0.10

0.20

Smooth, ρ=.5
Smooth, ρ=.75
Sparse, ρ=.5
Sparse, ρ=.75

Mean Squared Prediction Error

0.20

0.15

0.10

0.05

0.05

0.00

0.00
50

100

150

200

50

N

100

150

200

N

Figure 2: Prediction error and computation times for smooth and sparse methods. The top
row is the MSE and computation time for the BIC, while the first panel on the second row is
the MSE for AIC, and the second panel on the second row is the MSE for cross-validation.
axis (the false positive rate), the neither version of LASSO changes much, while FS-LASSO
quickly tops out, finding all significant predictors. This suggests that LASSO is able to find
a subset of predictors rather quickly, likely those which do not vary much over time, but is
unable to capture the more complex signals. Surprisingly, incorporating a nonparametric
estimate of the mean (TC LASSO) actually results in a large decrease in power. This further
supports the need for FS-LASSO, if one believes there is nonlinear structure in the data,
then patching LASSO is substantially worse than moving to a functional framework.
As a final comparison of statistical performance, we examine the sensitivity of the sparse
FS-LASSO on the choice of J. We consider J = 27 and J = 33, which constitutes a 10%
change in the number of basis functions. The resulting MSE curves are given in Figure 4.
Unfortunately, there does appear to be some sensitivity in the performance of FS-LASSO
relative to the number of basis functions. While J = 33 and J = 27 perform similarly, both
are slightly worse than J = 30, and in fact, worse than the smoothing algorithm. This can be
remedied by also choosing the number of basis functions by BIC, AIC, or CV, or by using the
smoothing method to select the predictors. Namely, one would ideally like to use the penalty
to also control the level of smoothing (not just the selection of the predictors). FS-LASSO
behaves very similarly to a ridge regression on the selected predictors, but a penalty which
acts more like a smoother would likely perform better (at least in terms of MSE). However,
16

ROC curves, ρ = 0.5

1.25

ROC curves, ρ = 0.75

1.25

0.75

True positive rate

1.00

True positive rate

1.00

0.75

Algorithm

0.50

LASSO
TC LASSO

N

0.25

Algorithm

0.50

LASSO
TC LASSO

N

0.25

50
100
200

0.00

50
100
200

0.00
0.000

0.025

0.050

0.075

0.100

0.000

0.025

0.050

False positive rate

0.075

0.100

False positive rate

Figure 3: ROC curves for scalar LASSO with (TC LASSO) and without (LASSO) a nonparametric estimate of the mean.
we do not explore this further here as this is an important topic of further research.
0.20

Mean Squared Prediction Error

Basis size 27
Basis size 30
Basis size 33
Sparse, ρ=.5
Sparse, ρ=.75

0.15

0.10

0.05

0.00
50

100

150

200

N

Figure 4: Mean squared prediction error for sparse FS-LASSO with differing numbers of
basis functions.
Lastly, we mention computation time and memory. Methods which scale well are critical
for applications such as genetic association studies. There, one needs to work with millions
of predictors and both computation time and memory become major issues. On a desktop
computer, the sparse method presented here can be applied with predictors of the order of
10,000 or so, but not much higher. The smoothing method, which allows for the inclusion
of an FPCA, can be applied with predictors of the order of 20,000. On clusters these
numbers can be increased depending on the available memory. In Figure 2 we plot average
computation times for the two methods when taking I = 10000, but keeping everything
else the same. We see that the smooth method has a substantial edge in computation time,
resulting in shorter times which also scale better with sample size.

17

rs10439990
0.010

Smooth
Sparse

200

Systolic Blood Pressure [mm Hg]

●

●

●

0.005
●

●
●

160

●

●

●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●●
●
●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●●
●●
●
●
●
●●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●● ● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●

0.000

120
−0.005

80

●

●

−0.010
20

40

60

80

Age [years]

20

40

60

80

Figure 5: Left Panel: Systolic blood pressure as a function of age from the FHS. The red
line is a local linear smoother estimating the mean. Right Panel: Plot of estimated effect
for SNP rs10439990 using the ‘smooth’ algorithm (blue) and the ‘sparse’ algorithm (red).

6

Framingham Heart Study

The Framingham Heart Study, FHS, is a long term longitudinal study with the goal of better understanding the risk factors for heart disease. Genome–wide SNP data and phenotype
information were downloaded from dbGaP (http://www.ncbi.nlm.nih.gov/gap) study accession phs000007.v25.p9. The study consists of three cohorts of subjects with the first cohort
recruited in 1958. FHS has had a tremendous impact on our understanding of health and
risk factors for heard disease. Since a full account would take far too much space, we refer the
interested reader to O’Donnell and Elosua (2008); Mahmood et al. (2014); Chen and Levy
(2016) for details on different findings and impacts due to FHS. Here we examine the second
cohort which consists of 1924 subjects. Each subject contributed up to 7 clinical exams over
29 years, though some subjects passed away during the study and contributed less. The subjects were genotyped using the Illumina Omni 5M platform resulting in approximately 4.3
million SNPs genotyped. Our goal here is to find a subset of SNPs impacting systolic blood
pressure, a primary risk factor for heart disease. As an illustration, in Figure 5 we plot all
of the blood pressure measurements versus age and include a local linear smoother. We see
that (as is well known) blood pressure increase monotonically, though nonlinearly, with age.
When we combine this fact with the high within-subject correlation of blood pressure, we see
that scalar procedures such as LASSO are inappropriate. When applying our procedures,
we first remove the effects of gender, height, and HDL cholesterol levels nonparametrically
using a local linear smoother.
Due to the size of the data, we cannot apply the methods to the entire set jointly. We
thus first use our method to screen the number of predictors down to a computationally
manageable size. One could take the top performers from marginal regression, but instead
we use our procedure to rank the SNPs according to when they enter the model (as one
18

varies the tuning parameter), and take the top ones. First, we divide up the genome into
segments of 10K SNPs for the sparse methods and 20K for the smooth. Using our algorithm,
we rank the SNPs in each segment according to the order in which they enter the model
(as one decreases the tuning parameter). We take the top ranking SNPS so that, when
pooled, we have screened down to approximately 10K SNPs for the sparse and 20K for the
smooth. A final application of our methods are applied to these SNPs with the smoothing
parameter chosen in four different ways: AIC, BIC, extended BIC with parameter 0.2, and
2-fold cross-validation. The top SNPs for each method are presented in Table 6. As we can
see, the smoothing method has found substantially more significant SNPs than the sparse.
However, there were 4 candidates which appeared in both lists. Given how similar the sparse
functional methods compared in simulations, it might at first be surprising that they only
overlap on four SNPs, however, given the huge number of predictors, it seems reasonable
that there would be substantial differences between the two. As an illustration, in Figure
5, we provide a plot of one of the estimated SNPs, rs10439990, which was selected using
both algorithms. There we see that the effect is negative, meaning that a subject with one
or two copies of the mutation has a lower chance of having high blood pressure. However,
the effect is not constant over time as it starts lower a younger ages, 20-40, and rises as
subjects age, peaking around 70 or so. While the two methods produce similar plots, they
disagree on the magnitude of the effect later in life, with the sparse algorithm indicating
that the protection is essentially gone, while the smooth algorithm indicates the protection
is still present, though diminished. However, we stress that the two submodels selected
were different, and thus the difference in estimates could be due in part to having different
predictors. This plot also illustrates a cautionary note about interpreting endpoints. With
so few observations at the ends, it is safer to focus interpretations on sections with the most
observations.
Quantifying the uncertainty of the output from these types of high-dimensional procedures
is an ongoing area of research even in the scalar setting. However, to test if either method is
dramatically over fitting, we also applied a 10-fold cross-validation to compare the predictive
performance between the sparse and smooth method. The resulting values were 0.00961 for
the sparse approach and 0.00972 for the smooth approach, thus showing that the sparse
approach is doing a slightly better job in terms of prediction. Given the inherently sparse
nature of the data, it is maybe surprising that the smooth approach was still so close in terms
of prediction. Regardless, as is now common in genetic studies, a follow up on a separate
data set would be required to validate the results. Note that for serious follow up studies,
it would be useful to extend techniques such as Shi (2015) to the functional setting, which
attempt to recover some of the relatively small effects that LASSO style procedures shrink
to zero. This is especially useful in genetic studies as the effect sizes are relatively small.
Lastly, an extension incorporating dependence between curves would useful as the FHS data
included some related individuals, though we do not pursue this further here.
As a final form of validation, we examined the different association results found in the
literature for our selected SNPs. This was accomplished using GWAS Central, http://www.
19

gwascentral.org/index, which provides a searchable database of GWAS results. There
were at least a few results that validate some of our findings. In particular, SNP rs10439990
is located the gene ZBTB20, which has been associated with Triglyceride levels (Kathiresan
et al., 2007), a common risk factor for heart disease. SNP kgp29965466 is located on the
gene ATRNL1, which is associated with several negative health outcomes and indicators
including nicotine dependence (Bierut et al., 2007). SNP rs10497371 is on gene MYO3B,
which has been associated with Diabetes (Scott et al., 2007). Lastly SNP rs7692492 is on
gene FAM190A, which has been associated with coronary heart disease (Samani et al., 2007;
Larson et al., 2007).
Top SNPs selected by sparse algorithm
Chromosome, Name AIC BIC EBIC CV
2, rs3845756
2, rs6414023
2, rs10497371
3, rs10439990
4, rs7692492
5, kgp30202888
7, kgp3808198
7, kgp4510449
8, kgp8137960
9, rs1702645
12, rs978561
13, rs1924783
15, kgp6228266

Top SNPs selected by smooth algorithm
Chromosome, Name AIC BIC EBIC CV
1, kgp5933154
1, rs17107710
2, kgp5982336
3, kgp26868438
3, rs10439990
3, kgp11928513
5, kgp30202888
7, kgp3808198
7, kgp4510449
10, kgp29965466
11, kgp10123049
12, kgp6953877
12, rs10859106

Figure 6: Top SNPs selected by the Sparse Algorithm and by the Smooth Algorithm. Gray
boxes indicated the SNP is selected when using the corresponding variable selection criteria.
Red indicates SNPs which were chosen by both the smooth and sparse methods.

7

Conclusion

We have provided powerful new tools for analyzing functional or longitudinal data with a
large number of predictors. In the sparse case we provided new theory in the form of a
restricted eigenvalue condition and accompanying asymptotic theory. While phrased as a
functional data method, the sparse case is closely related to varying coefficient models, and
thus this work can be viewed as an extension of the work in Wei et al. (2011). In the
dense case, we provide a completely new methodology and accompanying asymptotic theory
which allows the response functions to take values from any separable Hilbert space. Such
generality means the methods can be applied to a variety of settings including traditional
20

functional outcomes, functional panels, spatial processes, and image data such as fMRI. We
also provide accompanying computer code which takes advantage of the structure of group
LASSO so that our methods can be applied efficiently.
Our simulations suggest that the choice between sparse and smooth tool sets is not straightforward. The simulations were done in a traditionally sparse setting, but the two methods
were nearly equivalent in terms of variable selection and the functional method required
far fewer computational resources. However, estimation accuracy was better for the sparse
methods. This opens the door to future work which could involve multistage methods where
selection is done via smooth methods and estimation via sparse methods.
While our methods perform well, they also highlight the need for better tuning parameter
selection methods. Methods such as BIC and cross-validation can be applied in this setting
with good results, but such methods are not tailored to the dependence inherent in functional
data, and we hope to investigate possible adaptations to these selection criteria for the
functional setting in future work.

References
R. Baraniuk, M. Davenport, R. DeVore, and M. Wakin. A simple proof of the restricted isometry
property for random matrices. Constructive Approximation, 28(3):253–263, 2008.
P. Bickel, Y. Ritov, and A. Tsybakov. Simultaneous analysis of lasso and dantzig selector. The
Annals of Statistics, 37(4):1705–1732, 2009.
L. J. Bierut, P. A. Madden, N. Breslau, E. O. Johnson, D. Hatsukami, O. F. Pomerleau, G. E.
Swan, J. Rutter, S. Bertelsen, and L. Fox. Novel genes identified in a high-density genome wide
association study for nicotine dependence. Human molecular genetics, 16(1):24–35, 2007.
P. Bühlmann and S. Van de Geer. Statistics for high-dimensional data: methods, theory and
applications. Springer, 2011.
H. Cardot, F. Ferraty, and P. Sarda. Spline estimators for the functional linear model. Statistica
Sinica, 13:571–591, 2003.
G. Chen and D. Levy. Contributions of the framingham heart study to the epidemiology of coronary
heart disease. JAMA cardiology, 2016.
Y. Chen, J. Goldsmith, and R. T. Ogden. Variable selection in function-on-scalar regression. Stat,
2016.
J. Fan and J. Zhang. Two-step estimation of functional linear models with applications to longitudinal data. Journal of the Royal Statistical Society (B), 62:303–322, 2000.
Y. Fan, N. Foutz, G. M. James, W. Jank, et al. Functional response additive model estimation
with online virtual stock markets. The Annals of Applied Statistics, 8(4):2435–2460, 2014.

21

Y. Fan, G. James, and P. Radchenko. Functional Additive Regression. Annals of Statistics, 00
(000):0000–0000, 2015+.
J. Gertheiss, A. Maity, and A. Staicu. Variable selection in generalized functional linear models.
Stat, 2(1):86–101, 2013.
T. Hastie, R. Tibshirani, and J. Friedman. Elements of Statistical Learning. Springer, 2001.
D. Hoover, J. Rice, C. Wu, and L.-P. Yang. Nonparametric smoothing estimates of time–varying
coefficient models with longitudinal data. Biometrika, 85(4):809–822, 1998.
L. Horváth and P. Kokoszka. Inference for Functional Data with Applications. Springer, New York,
2012.
D. Hsu, S. Kakade, and T. Zhang. A tail inequality for quadratic forms of subgaussian random
vectors. Technical report, Technical report, Preprint, 2011.
G. James, T. Hastie, D. Witten, and R. Tibshirani. An Introduction to Statistical Learning:
With Applications in R. Springer Texts in Statistics. Springer London, Limited, 2013. ISBN
9781461471370. URL http://books.google.com/books?id=at1bmAEACAAJ.
G. M. James, J. Wang, and J. Zhu. Functional linear regression that’s interpretable. The Annals
of Statistics, 37:2083–2595, 2009.
S. Kathiresan, A. K. Manning, S. Demissie, R. B. D’agostino, A. Surti, C. Guiducci, L. Gianniny,
N. P. Burtt, O. Melander, and M. Orho-Melander. A genome-wide association study for blood
lipid phenotypes in the framingham heart study. BMC medical genetics, 8(1):1, 2007.
P. Kokoszka, I. Maslova, J. Sojka, and L. Zhu. Testing for lack of dependence in the functional
linear model. Canadian Journal of Statistics, 36:207–222, 2008.
D. Kong, K. Xue, F. Yao, and H. H. Zhang. Partially functional linear regression in high dimensions.
Biometrika, page asv062, 2016.
M. G. Larson, L. D. Atwood, E. J. Benjamin, L. A. Cupples, R. B. D’Agostino, C. S. Fox, D. R.
Govindaraju, C.-Y. Guo, N. L. Heard-Costa, and S.-J. Hwang. Framingham heart study 100k
project: genome-wide associations for cardiovascular disease outcomes. BMC medical genetics,
8(1):1, 2007.
B. Laurent and P. Massart. Adaptive estimation of a quadratic functional by model selection.
Annals of Statistics, 28:1302–1338, 2000.
Y. Li and T. Hsing. Uniform convergence rates for nonparametric regression and principal component analysis in functional/longitudinal data. The Annals of Statistics, pages 3321–3351, 2010.
H. Lian. Shrinkage estimation and selection for multiple functional regression. Statistica Sinica,
23:51–74, 2013.

22

K. Lounici, M. Pontil, S. Van de Geer, and A. Tsybakov. Oracle inequalities and optimal inference
under group sparsity. The Annals of Statistics, 39(4):2164–2204, 2011.
S. S. Mahmood, D. Levy, R. S. Vasan, and T. J. Wang. The framingham heart study and the
epidemiology of cardiovascular disease: a historical perspective. The Lancet, 383(9921):999–
1008, 2014.
H. Matsui and S. Konishi. Variable selection for functional regression models via the L1 regularization. Computational Statistics & Data Analysis, 55(12):3304–3310, 2011.
H. Müller. Functional modeling of longitudinal data. In Longitudinal Data Analysis, pages 223–252.
CRC Press, 2008.
C. J. O’Donnell and R. Elosua. Cardiovascular risk factors. insights from framingham heart study.
Revista Espanola de Cardiologia (English Edition), 61(3):299–310, 2008.
R. Oliveira. The lower tail of random quadratic forms, with applications to ordinary least squares
and restricted eigenvalue properties. arXiv preprint arXiv:1312.2903, 2013.
J. Ramsay and B. Silverman. Functional Data Analysis. Springer, 2005.
M. Reimherr and D. Nicolae. A functional data analysis approach for genetic association studies.
The Annals of Applied Statistics, 8:406–429, 2014.
P. Reiss, M. Mennes, and L. Huang. Fast function–on–scalar regression with penalized basis expansions. International Journal of Biostatistics, 6:Article 28, 2010.
N. J. Samani, J. Erdmann, A. S. Hall, C. Hengstenberg, M. Mangino, B. Mayer, R. J. Dixon,
T. Meitinger, P. Braund, and H.-E. Wichmann. Genomewide association analysis of coronary
artery disease. New England Journal of Medicine, 357(5):443–453, 2007.
L. J. Scott, K. L. Mohlke, L. L. Bonnycastle, C. J. Willer, Y. Li, W. L. Duren, M. R. Erdos, H. M.
Stringham, P. S. Chines, and A. U. Jackson. A genome-wide association study of type 2 diabetes
in finns detects multiple susceptibility variants. science, 316(5829):1341–1345, 2007.
P. Shi. Weak signal identification and inference in penalized model selection. PhD thesis, University
of Illinois at Urbana-Champaign, 2015.
R. Tibshirani. The lasso problem and uniqueness. Electronic Journal of Statistics, 7:1456–1490,
2013.
J. Wang, B. Lin, P. Gong, P. Wonka, and J. Ye. Lasso screening rules via dual polytope projection.
In Advances in Neural Information Processing Systems, pages 1070–1078, 2013.
L. Wang, H. Li, and J. Huang. Variable selection in nonparametric varying-coefficient models for
analysis of repeated measurements. Journal of the American Statistical Association, 103(484):
1556–1569, 2008.

23

F. Wei, J. Huang, and H. Li. Variable selection and estimation in high-dimensional varyingcoefficient models. Statistica Sinica, 21(4):1515–1540, 2011.
F. Yao, H.-G. Müller, and J.-L. Wang. Functional linear regression analysis for longitudinal data.
The Annals of Statistics, 33:2873–2903, 2005.
X. Zhang and J. Wang. From sparse to dense functional data and beyond. Technical report,
Preprint, 2015+.
P. Zhao and L. Xue. Variable selection for semiparametric varying coefficient partially linear errorsin-variables models. Journal of Multivariate Analysis, 101(8):1872–1883, 2010.
H. Zhu, R. Li, and L. Kong. Multivariate varying coefficient model for functional responses. Annals
of Statistics, 40(5):2634–2666, 2012.

24

A

Binomial Simulations

Here we carry out a nearly identical simulation study as in Section 5 except the normally
distributed predictors are replaced with Bin ∼ Binom(2, pi ) random variables to more closely
emulate the structure of the application. To generate Bin which are correlated, we use a
probit model based on the the same simulation structure as in Section 5; we take those
normal random variables as the input in probit model where pi = logit−1 (Xin ). The MSE
curves are summarized in Figure 7. If we compare these results with those from the bottom
right panel from Figure 2, here we also use CV, we see that the pattern is very similar. The
error does not change dramatically between the two correlations, but the sparse method has
a noticeably lower MSE.
Mean Squared Prediction Error

Smooth, ρ=.5
Smooth, ρ=.75
Sparse, ρ=.5
Sparse, ρ=.75

0.6

0.4

0.2

0.0
50

100

150

200

N

Figure 7: MSE curves comparing the smooth and sparse methods with Binom(2, 1/2) predictors. Tuning parameter is selected with cross-validation.

B
B.1

Proofs for Section 2
Proof of Theorem 1

Our analysis for this Theorem closely follows standard techniques for sparse and group sparse
regression under restriced eigenvalue conditions, such as the work by Lounici et al. (2011)
and Bickel et al. (2009).
b is a minimizer of
Recall that B
1
L(B) = kY − A vec(B> )k22 + λkBk`1 /`2 .
2
b ≤ L(B? ), that is,
Then in particular, we must have L(B)
1
b > )k2 + λkBk
b ` /` ≤ 1 kY − A vec(B?> )k2 + λkB? k` /` .
kY − A vec(B
2
2
1 2
1 2
2
2
After rearranging some terms,
D
E


1
b > − B?> )k2 ≤ Y − A vec(B?> ), A vec(B
b > − B? ) + λ kB? k` /` − kBk
b ` /` .
kA vec(B
2
1 2
1 2
2
25

b − B? . Then
Let S0 ⊂ [I] indicate the row support of B? , with |S| ≤ I0 . Write ∆ = B
standard arguments show that
b ` /` ≤ k∆S0 k` /` − k∆S c k` /` .
kB? k`1 /`2 − kBk
1 2
1 2
1 2
0
Furthermore,
D
E

b > − B?> ) = A> Y − A vec(B?> ) , vec(∆> )
Y − A vec(B?> ), A vec(B
X

?>
=
A>
Y
−
A
vec(B
)
, ∆i,∗
i
i

≤

X


?>
A>
Y
−
A
vec(B
)
i

· k∆i,∗ k2
2



i


?>
≤ max A>
i Y − A vec(B )
i

2

· k∆k`1 /`2 ,

where Ai is the matrix with entries (Ai )nm,j = Anm,ij = Xni · ej (tnm ).

?>
Y
−
A
vec(B
)
. On this event, we would then
Consider the event λ ≥ 2 maxi A>
i
2
have


1
?>
c k` /`
kA vec(∆> )k2 ≤ k∆k`1 /`2 · max A>
+
λ
k∆
k
−
k∆
S
S
`
/`
0
i Y − A vec(B )
1
2
1
2
0
2
i
2
≤ k∆S0 k`1 /`2 · 1.5λ − k∆S0c k`1 /`2 · 0.5λ .
Since the left-hand side is nonnegative, therefore, k∆S0c k`1 /`2 ≤ 3k∆S0 k`1 /`2 . Applying the
REgroup (I0 , α) condition, then,
kA vec(∆> )k2 ≥ N α k∆k2F .
Combining everything,
1
1
N α k∆k2F ≤ kA vec(∆> )k2 ≤ k∆S0 k`1 /`2 · 1.5λ − k∆S0c k`1 /`2 · 0.5λ ≤ k∆S0 k`1 /`2 · 1.5λ .
2
2
For the Frobenius norm result, we further write
p
p
k∆S0 k`1 /`2 · 1.5λ ≤ |S0 | · k∆S0 kF · 1.5λ ≤ |S0 | · k∆kF · 1.5λ
and so
k∆kF ≤

p

|S0 | · 1.5λ ·

2
Nα

while for the `1 /`2 -norm result we have
p
p
k∆k`1 /`2 = k∆S0 k`1 /`2 + k∆S0c k`1 /`2 ≤ 4k∆S0 k`1 /`2 ≤ 4 |S0 |k∆S0 kF ≤ 4 |S0 |k∆kF .
This proves the theorem, as long as we can show that with high probability,

?>
λ ≥ 2 max A>
.
i Y − A vec(B )
2
i

26

To prove this, take any i ∈ [I]. Then

?>
>
>
>
kA>
i Y − A vec(B ) k2 = kAi (T + ) k2 ≤ kAi T k2 + kAi k2 ,
where T is the truncation error defined earlier while  has entries nm = n (tnm ). Now we
bound the two pieces separately. First, we bound the operator norm of Ai . By definition,
(Ai )nm,j = Xni · ej (tnm ) and so we have
X
X
2
2
2
E>
E>
Xni
A>
n En = kXk∞ · N · F ,
n En  kXk∞
i Ai =
n

n

proving that kAi kop ≤ kXk∞ ·

p

N kFkop . Then

kA>
i T k2 ≤ kAi kop · kT k2 ≤ kXk∞ ·

q
N kFkop · kT k2 .

>
>
Next, to bound kA>
i k2 , observe that Ai  ∼ N (0, Ai ΣAi ), where where Σ is a blockdiagonal matrix with blocks given by Σ1 , . . . , ΣN . Then, applying (Hsu et al., 2011, Proposition 1.1),


q
>
> 2
>
>
2
P kAi k2 > trace(Ai ΣAi ) + 2 trace((Ai ΣAi ) ) · t + 2kAi ΣAi kop · t ≤ e−t

for
√ any t ≥ 0. Since trace(M ) ≤ JkM kop for any J × J matrix, and using the simple identity
2 Jt ≤ J + t, we simplify this to the weaker statement

2
>
−t
P kA>
i k2 > kAi ΣAi kop · (2J + 3t) ≤ e
We calculate
2
2
kA>
i ΣAi kop ≤ kAi kop · kΣkop ≤ kXk∞ · N kFkop · max kΣn kop .
n

Combining everything, and setting t = log(I/δ), for each I we have
n
o δ
2
> 2
P kAi k2 > kXk∞ · N kFkop · max kΣn kop · (2J + 3 log(I/δ)) ≤ .
n
I
Therefore, with probability at least 1 − δ, for all i = 1, . . . , I,
2
2
kA>
i k2 ≤ kXk∞ · N kFkop · max kΣn kop · (2J + 3 log(I/δ))
n

and so for all i,

?>
>
>
kA>
i Y − A vec(B ) k2 ≤ kAi T k2 + kAi k2


q
q
≤ kXk∞ N kFkop · kT k2 + max kΣn kop · (2J + 3 log(I/δ))
n

≤ λ/2 ,
as desired.
27

B.2

Proof of Theorem 2

We use the Transfer Principle of Oliveira (2013, Lemma 5.1), which connects sparse eigenvalues of A to a restricted eigenvalue property for A. This result is stated for sparsity but
can be straightforwardly extended to group-sparsity. We restate the group-sparse form here
as a Lemma without proof.
Lemma 1 (Adapted from Oliveira (2013, Lemma 5.1)). Suppose A satisfies that, for any
k1 -row-sparse W ∈ RI×J ,
2

A vec(W> )

(7)

2

≥ λ0 kWk2F .

Then for all W ∈ RI×J ,
>

A vec(W )

2
2

≥

λ0 kWk2F

−

maxi kAi k2op · kWk2`1 /`2
k1 − 1

.

To apply this result to our work, we need to first find a λ0 such that (7) is satisfied (with
high probability), and then we need to compute a bound (holding with high probability) on
each kAi kop .
Step 1: finding λ0 for condition (7) We first give another lemma, proved below:
Lemma 2. Suppose that the assumptions of Theorem 2 hold, and take any fixed sequence
{En } satisfying (4). Choose any δ > 0 and any k1 ≥ 1. Then there are constants a1 , a2 > 0
depending only on ν, λmin (Σ), γ0 , γ1 , such that if
(8)

N ≥ a1 (k1 J + k1 log(I) + log(1/δ)),

then with probability at least 1 − δ/2, for all k1 -row-sparse W,
P
>
E>
n W Xn 2
pn
≥ a2 .
N trace(W> ΣW)
To apply these results to our work, first note that by Lemma 2, if N satisfies (8), then
with probability at least 1 − δ/2, for all k1 -sparse W (where we specify k1 later),
kA · vec(W)k22 =

X

>
E>
n W Xn

n

2
2

1
≥
N

!2
X

>
E>
n W Xn

2

n

≥ N trace(W> ΣW) · (a2 )2 .
In particular, the assumption (7) in Lemma 1 holds with λ0 = N (a2 )2 λmin (Σ).

28

Step 2: bounding kAi kop Next we bound maxi kAi k2op . From the proof of Theorem 1, we
P
know that kAi k2op ≤ kXk2∞ · N kFkop , where F = N1 n E>
n En . By (4), we see that kFkop ≤
2
γ1 . Furthermore, Xij is mean zero and ν-subgaussian, and so, by standard subgaussian tail
bounds, with probability at least 1 − δ/2,
p
kXk∞ = max |Xni | ≤ ν 2 log(4IJ/δ) .
ni

Step 3: applying Lemma 1 Combining the results from Step 1 and Step 2, and applying
Lemma 1, we see that with probability at least 1 − δ, for all W ∈ RI×J ,
#
"
2
2
2
2ν
log(4IJ/δ)
·
γ
·
kWk
2
1
`
/`
2
1
2
.
A vec(W> ) 2 ≥ N (a1 )2 λmin (Σ) · kWkF −
k1 − 1
For any W with kWS c k`1 /`2 ≤ 3 kWS k`1 /`2 (where |S| ≤ k), we then have
√
√
kWk`1 /`2 ≤ 4 kWS k`1 /`2 ≤ 4 k kWS kF ≤ 4 k kWkF ,
and so the result above gives
>

A vec(W )

2
2

≥N

kWk2F



2ν 2 log(4IJ/δ) · γ12 · 16k
(a2 ) λmin (Σ) −
k1 − 1
2


.

Taking
k1 ≥ 1 +

2ν 2 log(4IJ/δ) · γ12 · 16k
,
0.5(a2 )2 λmin (Σ)

we obtain
A vec(W> )

2
2

≥ N kWk2F · 0.5(a2 )2 λmin (Σ)

for all W ∈ RI×J with kWS c k`1 /`2 ≤ 3 kWS k`1 /`2 for any |S| ≤ k, as desired. The sample
size requirement (5) in Theorem 2 ensures that the sample size assumption in Lemma 2 will
hold for the specified choice of k1 .
B.2.1

Proof of Lemma 2

For the proof of this lemma, we’ll need several supporting lemmas, proved below. The first
P
>
lemma shows that n E>
n W Xn 2 satisfies upper and lower bounds with high probability
for any fixed single matrix W.
Lemma 3. For any fixed W,
(
P
P c1 ≤

N

n

>
E>
n W Xn

2
p
≤ c2
trace(W> ΣW)

where c1 , c2 , c3 > 0 depend only on ν, λmin (Σ), γ0 , γ1 .
29

)
≥ 1 − 2e−c3 N ,

P
>
The next lemma shows that proving a lower bound on n E>
n W Xn 2 over all rowsparse W can be reduced to proving a lower bound on a finite covering set. This result
is a simple extension of Baraniuk et al. (2008, Lemma 5.1), from the sparse setting to the
group-sparse setting, and so we do not give the proof here.
Lemma 4 (Adapted from Baraniuk et al. (2008, Lemma 5.1)). Let

W = W ∈ RI×J : W is k1 -row-sparse, trace(W> ΣW) = 1

.

Let f (W) be any seminorm (that is, f (c · W) = |c| · f (W) and f (W1 + W2 ) ≤ f (W1 ) +
f (W2 )). Choose any  ∈ (0, 1/2). Then there exists W 0 ⊂ W with |W 0 | ≤ I k1 (3/)k1 J such
that

inf f (W) ≥ inf 0 f (W) −
sup f (W) .
W∈W
W∈W
1 −  W∈W 0
Proof of Lemma 2. By Lemma 3, for any fixed W, with probability at least 1 − 2e−c3 N ,
P
>
E>
n W Xn 2
≤ c2
c1 ≤ pn
N trace(W> ΣW)
where c1 , c2 , c3 > 0 depend only on ν, λmin (Σ), γ0 , γ1 . Now choose  = 4cc12 . In the notation
P
>
0
0
of Lemma 4, let f (W) = n E>
n W Xn 2 , and take W as in the lemma. Since |W | ≤
I k1 (3/)k1 J , the above bound is true for all W ∈ W 0 with probability at least
1 − 2e−c3 N +log(I

k1 (3/)k1 J )

≥ 1 − δ/2,

by our lower bound N ≥ a1 (k1 J + k1 log(I) + log(1/δ)) as long as the constant a1 is chosen
appropriately. Next, applying Lemma 4, for all W ∈ W,
P
>
E>

c1
n W Xn 2
pn
≥ c1 −
· c2 ≥ .
>
1−
2
N trace(W ΣW)
After defining a2 = c1 /2, this proves the lemma.
B.2.2

Proof of Lemma 3

We begin by stating a supporting lemma.
Lemma 5. Let B1 ∈ RM1 ×d , . . . , BN ∈ RMN ×d be fixed matrices with
N
X

kBn kF ≥ c1 N

n=1

and

v
u N
uX
√
t
kBn k2F ≤ c2 N
n=1

30

for some c1 , c2 > 0. Let X1 , . . . , XN ∈ Rd be iid random
i with E2[Xn ] = 0 and
h vectors


2
v > Xn
>
≤ eσ kvk2 /2 for all fixed
E Xn Xn = Id , which are σ-subgaussian, meaning that E e
v ∈ Rd . Then

(
P c3 N ≤

N
X

)
kBn Xn k2 ≤ c4 N

≥ 1 − 2e−c5 N ,

n=1

where c3 , c4 , c5 > 0 depend only on c1 , c2 , σ.
Proof of Lemma 3. First we calculate
X

> 1/2
E>
nW Σ

2
F

=

X

n

>
trace(E>
n W ΣWEn )

n

= trace(Σ1/2 W(

X

> 1/2
En E>
) ≤ trace(Σ1/2 WW> Σ1/2 ) · N
n )W Σ

n

1 X
En E>
n
N n

1 X
>
En E>
n ≤ γ1 N trace(W ΣW) .
N n

= trace(W> ΣW) · N
And,
X
n

> 1/2
E>
nW Σ
F

=

X

s
X

n

2

k(Σ1/2 W)i,∗ En k2 ≥

p
n

i

≥ γ0 N

X
i

(Σ1/2 W)i,∗

XX

(Σ1/2 W)i,∗

i

2
>
trace(W ΣW)

2
2

p
= γ0 N
trace(W> ΣW)

(Σ1/2 W)i,∗ En

p
trace(W> ΣW) .

en = Σ−1/2 Xn and defining
Now, writing X
>
E>
nW
,
Bn = p
trace(W> ΣW)

p
en ’s are iid with mean zero, identity covariance, and are ν/ λmin (Σ)we see that the X
p
>
en . Applying Lemma 5,
subgaussian. And, E>
trace(W> ΣW) · Bn Σ1/2 X
n W Xn 2 =
2

then, for some constants c1 , c2 , c3 depending only on ν, λmin (Σ), γ0 , γ1 ,
(
)
X
p
p
>
P c1 N trace(W> ΣW) ≤
E>
trace(W> ΣW) ≥ 1 − 2e−c3 N .
n W X n 2 ≤ c2 N
n

B.2.3

Proof of Lemma 5

Before we prove Lemma 5, we first state two additional supporting lemmas, proved below:

31

2

√
Lemma 6. Let x ∈ RN be a vector with kxk1 ≥ c1 N , kxk2 ≤ c2 N . Then
n
c1 o
(c1 )2
n : |xn | ≥
≥
N.
2
4(c2 )2
2
Lemma 7.
i Z1 , . . . , ZN are independent non-negative random variables with E [(Zn ) ] ≥
h Suppose
2
a1 and E e(Zn ) ≤ a2 for all n = 1, . . . , N , where a1 , a2 > 0. Then

P

(
X

)
Zn < a3 N

≤ e−a4 N

and

(
)
X
P
(Zn )2 > a5 N ≤ e−a6 N

n

n

where a3 , a4 , a5 , a6 > 0 are constants depending only on a1 , a2 .
Now we prove Lemma 5.
Proof of Lemma 5. For each n, define Zn =

kBn Xn k2
.
2σkBn kF

We have







E kBn Xn k22 = E Xn> Bn> Bn Xn = E trace(Bn> Bn Xn Xn> ) = trace(Bn> Bn ) = kBn k2F ,
and so E [(Zn )2 ] =

1
4σ 2

=: a1 . Next, it is known from Hsu et al. (2011, Remark 2.3) that

h
i
σ 4 trace((Bn> Bn )2 )η 2
ηkBn Xn k22
2
>
E e
≤ σ trace(Bn Bn )η +
1 − 2σ 2 kBn> Bn k η
for all 0 ≤ η <

1
.
>B
2σ 2 kBn
nk

Using basic facts about matrix norms we have

kBn k2 ≤ kBn k2F and trace(Bn> Bn ) = kBn k2F and trace((Bn> Bn )2 ) =
1
Setting η = 4σ2 kB
we obtain
k2

2
(Bn> Bn )2 F

Bn> Bn

≤ Bn> Bn

=
4
.
F

n F

h
i
h
i
2
2
2
2
E e(Zn ) = E ekBn Xn k2 /(4σ kBn kF ) ≤ e3/8 =: a2 .
Applying Lemma 7 to {Zn : n = 1, . . . , n}, we see that for constants a5 , a6 > 0 depending
only on a1 , a2 ,
)
(
X kBn Xn k2
2
≤ a5 N ≥ 1 − e−a6 N .
P
2 kB k2
4σ
n F
n
Furthermore, setting I = {n : kBn kF ≥ c1 /2}, by applying Lemma 6 to the vector
(c1 )2
(kBn kF )n=1,...,N , we see that |I| ≥ 4(c
2 N . Applying Lemma 7 to {Zn : n ∈ I} we see that
2)
for constants a3 , a4 > 0 depending only on a1 , a2 ,
(
)
(c )2
X kBn Xn k
−a4 1 2 N
−a4 |I|
2
4(c2 )
≥ a3 |I| ≥ 1 − e
≥1−e
.
P
2σ kBn kF
n∈I
Next, assume these events hold. We have
v
u
s
X
X kBn Xn k
X
p
√
uX kBn Xn k22
2
2
t
· kBn kF ≤
·
kB
k
≤
2σ
a
N
·
c
N.
kBn Xn k2 =
n
5
2
F
kBn kF
kBn k2F
n
n
n
n
32

Furthermore,
X

kBn Xn k2 ≥

n

X

kBn Xn k2 =

n∈I

X

2σ kBn kF ·

n∈I

≥ c1 ·

kBn Xn k2
2σ kBn kF

X kBn Xn k

2

n∈I

2σ kBn kF

≥ c1 · a3 |I| ≥ c1 · a3 ·

(c1 )2
N.
4(c2 )2

Setting


√
(c1 )2
(c1 )2
c 3 = c 1 · a3 ·
,
, c4 = 2σc2 a5 , c5 = min a6 , a4
4(c2 )2
4(c2 )2
we have proved the lemma.
Finally, we prove the two supporting results, Lemmas 6 and 7.
Proof of Lemma 6. We have
c1 N ≤
≤

N
X
n=1
N
X

|xn | =

N
X

|xn | · 1 {|xn | ≥ c1 /2} +

N
X

|xn | · 1 {|xn | < c1 /2}

n=1

n=1

|xn | · 1 {|xn | ≥ c1 /2} + N · c1 /2

n=1

v
v
u N
u N
uX
uX
t
≤
|xn |2 · t
1 {|xn | ≥ c1 /2}2 + N · c1 /2
n=1

≤ c2

√

n=1

p
N · |{n : |xn | ≥ c1 /2}| + N · c1 /2,

and therefore, |{n : |xn | ≥ c1 /2}| ≥

(c1 )2
N.
4(c2 )2

Proof of Lemma 7. For the upper bound, we have
(
)
h P
i
X
2
P
(Zn )2 > 2a2 N ≤ E e n (Zn ) −2a2 N =
n

!
i
Y h
2
E e(Zn )
· e−2a2 N = e−a2 N .
n

For the lower bound, we have, for each n,
p


 p

a1 ≤ E (Zn )2 = E (Zn )0.5 (Zn )1.5 ≤ E [Zn ] · E [(Zn )3 ] since Zn ≥ 0
q
p
p
p
p
√
≤ E [Zn ] · 4 E [(Zn )6 ] ≤ E [Zn ] · 4 E [6e(Zn )2 ] ≤ E [Zn ] · 4 6a2

33

and so E [Zn ] ≥


−tZn

E e



a2
√1 .
6a2

Next, taking any t ∈ [0, 1],



X (−1)k tk E (Zn )k
= 1 − tE [Zn ] +
k!
k≥2
!

X t2k E (Zn )2k
≤ −tE [Zn ] + 1 +
by removing negative terms
(2k)!
k≥1
h 2 2i
= −tE [Zn ] + E et (Zn )
since k! ≤ (2k)!
h
i
2 (Zn )2
2
≤ −tE [Zn ] + E t e
+ (1 − t ) · 1
by convexity of z 7→ ez

a2
≤ 1 − t √ 1 + t2 (a2 − 1).
6a2
o
n
a2
Setting t = min 1, 2(a2 −1)1 √6a2 we get
2



−tZn

E e



(a )
−t √1
(a1 )2
≤1−t √
≤ e 2 6a2 .
2 6a2

Then
(
)


P
(a )2
X
tN √1 −t n Zn
(a1 )2
4 6a2
Zn < √ N ≤ E e
P
4 6a2
n
!

 Y
(a )2
(a )2
(a )2
(a )2


tN √1
−tN √1
tN √1 −tN √1
−tZ
n
4 6a2
4 6a2
2 6a2
4 6a2
=E e
E e
=e
.
=e
n

Setting


(a1 )2
a21
(a1 )2
√
a3 = √ , a4 = min 1,
· √ , a5 = 2a2 , a6 = a2 ,
4 6a2
2(a2 − 1) 6a2
4 6a2
we have proved the lemma.

C

Proofs for Section 3

We begin this section with the proof of Theorem 3. After, we will divide the proof of Theorem
4 into a sequence of lemmas.

C.1

Proof of Theorem 3

Let e1 , e2 , . . . be an orthonormal basis of H. Then each coordinate xi of x ∈ HI , can be
expressed as
X
xi =
xi,k ek .
k

34

We denote xi,k ∈ R as the coordinates with respect to the e1 , e2 , . . . basis, and we will also
let x(k) represent the vector {x1,k , . . . , xI,k }. We can then write
kXxk2H = hXx, XxiH =

N
X

hXn x, Xn xi

n=1

=

N X
I X
I
X

Xni Xnj hxi , xj iH

n=1 i=1 j=1

Expressing the xi with respect to the orthonormal basis e1 , e2 , . . . ,
=

XXXXX

=

XXXX

n

i

n

i

=

X

=

X

j

j

k

Xni Xnj xi,k xj,l hek , el i

l

Xni Xnj xi,k xj,k

k

x(k)> X> Xx(k)

k

kXx(k) k22 .

k

Note that the last norm above is simply the Euclidean norm on RN . By assumption in the
Theorem, we have that the above is bounded from below as follows:
X
kXxk2H =
kXx(k) k22
k

2
X √
p
≥
c1 N kx(k) k2 − c2 log(I)kx(k) k1
k

= (c1 )2 N

X

kx(k) k22 − 2c1 c2

X
X
p
kx(k) k21 .
kx(k) k2 kx(k) k1 + (c2 )2 log(I)
N log(I)

k

k

k

Examining the first sum in this last line, we have by Parceval’s identity
X
XX
X
kxk2H =
kxi k2H =
|xi,k |2 =
kx(k) k2 .
i

i

k

k

For the second sum, the Cauchy-Schwartz inequality gives
sX
X
(k)
(k)
kx k2 kx k1 ≤ kxkH
kx(k) k21 .
k

k

So we have that

kXxk2H ≥ c1

√

2
sX
p
N kxkH − c2 log(I)
kx(k) k21 
k

35

or equivalently
kXxkH ≥ c1

√

sX
p
N kxkH − c2 log(I)
kx(k) k21 .
k

Now for the last piece observe that (by Cauchy-Schwarz)
X
XXX
kx(k) k21 =
|xi,k ||xj,k |
k

i

k

j

!1/2
≤

XX X
i

j

|xi,k |2

k

|xj,k |2

k

!1/2 2


=

!1/2
X

X X
i

|xi,k |2



k

!2
=

X

= kxk2`1 /H .

kxi kH

i

We can therefore conclude that
p
√
kXxkH ≥ c1 N kxkH − c2 log(I)kxk`1 /H
as claimed.
q
Next suppose that c1 > 4c2 I0 log(I)
. We will show that X satisfies the REF (I0 , α) property

2 N
q
with α = c1 − 4c2 I0 log(I)
. To see this, take any set S with |S| ≤ I0 , and any x with
N
kxS c k`1 /H ≤ 3 kxS k`1 /H . By the Cauchy-Schwarz inequality, which relates the `1 and `2
norms of any vector (including a vector of functions xS ), we have
p
p
kxS k`1 /H ≤ kxS kH · |Support(xS )| ≤ kxS kH · I0 .
By the triangle inequality,
kxk`1 /H ≤ kxS k`1 /H + kxS c k`1 /H ≤ 4 kxS k`1 /H ≤ 4 kxS kH ·

p
p
I0 ≤ 4 kxkH · I0 .

Then returning to the above,

p 
p
p
√
√
kXxkH ≥ c1 N kxkH − c2 log(I)kxk`1 /H ≥ c1 N kxkH − c2 log(I) 4 kxkH · I0
!
r
√
√
I0 log(I)
= kxkH · αN ,
= kxkH · N · c1 − 4c2
N
proving that X satisfies the REF (I0 , α) property.

36

C.2

Proof of Theorem 4

We follow a very similar structure to the one found in Chapter 6 of Bühlmann and Van de
Geer (2011), adapting the arguments for function spaces as needed. We break the proof into a
sequence of Lemmas which can be thought of as a basic inequality, a functional concentration
inequality, applying that inequality, a last lemma to set up the final proof. We then combine
all the lemmas to prove Theorem 4.
Lemma 8 (Basic Inequality). With probability one, we always have the inequality


1
?
2
>
(i)
b ` /H ≤ max kε X k kβb − β ? k` /H + λkβ ? k` /H .
kX(βb − β )kH + λkβk
1
1
1
1≤i≤I
2
Proof. We can set up a basic inequality using
1
b 2 + λkβk
b ` /H ≤ 1 kY − Xβ ? k2 + λkβ ? k` /H .
kY − Xβk
1
1
2
2
We can rewrite
1
b 2 = 1 kY − Xβ ? − X(βb − β ? )k2
kY − Xβk
2
2
1
1
2
= kY − Xβ ? k2 + kX(βb − β ? )k2 − hY − Xβ ? , X(βb − β ? )i.
2
2
2
This gives
1
b ` /H ≤ hY − Xβ ? , X(βb − β ? )i + λkβ ? k` /H
kX(βb − β ? )k2 + λkβk
1
1
2
I
X
hε> X(i) , βbi − βi? i + λkβ ? k`1 /H
=
i=1

≤

I
X

kε> X(i) kkβbi − βi? k + λkβ ? k`1 /H

i=1

≤ max kε> X(i) kkβb − β ? k`1 /H + λkβ ? k`1 /H ,
1≤i≤I

which is the desired result.
Lemma 9. Let X be an H valued Gaussian process with mean zero and covariance operator
C. Let Λ> = (λ1 , λ2 , . . . ) be a vector of the eigenvalues of C (in decreasing order). Then we
have the bound
o
n
√
P kXk2 ≥ kΛk1 + 2kΛk2 t + 2kΛk∞ t ≤ exp(−t).
Proof. Using the Karhunen-Loeve expansion, we can express
D

kXk2 =

∞
X
j=1

37

λj Zj2 ,

where {Zj } are iid standard normal. Note that since X is square integrable, we have that
P
i λi = kΛk1 < ∞, and therefore kΛk2 < ∞ as well. Define the events, for J = 1, 2, . . .
( J
)
X
√
AJ =
λj Zj2 ≥ kΛk1 + 2kΛk2 t + 2λ1 t .
j=1

Since kΛk1 ≥
(2000) that

PJ

2
j=1 λj and kΛk2 ≥

PJ

j=1

λ2j , it follows from Lemma 1 in Laurent and Massart

P (AJ ) ≤ exp(−t),
for all J. Since

PJ

j=1

λj Zj2 is a strictly increasing sequence of J we have that
A1 ⊆ A2 ⊆ A3 . . .

Using the continuity from below of probability measures, we can pass the limit J → ∞ to
obtain the desired result.
Lemma 10. Assume that ε = (ε1 , . . . , εN ) are iid Gaussian elements of H and X is a design
matrix with standardized columns. Define the event


> (i)
F = max kε X k ≤ λ0
1≤i≤I

where, for δ > 0,
p
√ q
λ0 = N kΛk1 + 2kΛk2 log(I/δ) + 2λ1 log(I/δ).
Then we have that
P (F ) ≥ 1 − δ.
Proof. Notice
that if {εn } are iid Gaussian with mean zero and covariance operator C, then
√
> (i)
ε X / N is also mean zero Gaussian, and the covariance operator is given by
N
1
1 X
Cov(ε> X(i) ) =
Cov(εn Xni ) = C,
N
N n=1

√
since the columns of X are standardized. Therefore, Vi = ε> X(i) / N is also Gaussian with
mean zero and covariance operator C. Now apply Lemma 9 with t = − log(δ) to obtain


p
2
P max kVi k ≥ kΛk1 + 2kΛk2 t + log(I) + 2kΛk∞ (t + log(I))
1≤i≤I

≤ I exp(−(t + log(I)) = exp(−t).
Since exp(−t) = δ, the claim holds.

38

Lemma 11. On F with λ ≥ 2λ0 we have that
b 2 + λkβbS c k` /H ≤ 3λk(βb − β ? )S k` /H .
kX(β ? − β)k
0
1
1
0
Proof. Applying Lemma 8 and multiplying both sides by 2, we have that, on F ,
b 2 + 2λkβk
b ` /H ≤ λkβb − β ? k` /H + 2λkβ ? k` /H .
kX(β ? − β)k
1
1
1
Analyzing the penalty term on the LHS, we can apply the reverse triangle inequality to
obtain
b ` /H = kβbS k` /H + kβbS c k` /H ≥ kβ ? k` /H − kβbS − β ? k` /H + kβbS c k` /H
kβk
0
0
S0 1
S0 1
1
1
1
1
0
0
and on the RHS we we have
kβb − β ? k`1 /H = kβbS0 − βS?0 k`1 /H + kβbS0c k`1 /H .
Combining these two calculations we arrive at the bound
b 2 ≤λkβb − β ? k` /H + 2λkβ ? k` /H − 2λkβk
b ` /H
kX(β ? − β)k
1
1
1
?
b
b
≤λ(kβS − β k` /H + kβS c k` /H ) + 2λkβ ? k`
0

S0

1

S0

1

0

1 /H

− 2λ(kβS?0 k`1 /H − kβbS0 − βS?0 k`1 /H + kβbS0c k`1 /H )
=3λkβbS − β ? k` /H − λkβbS c k` /H ,
0

S0

1

0

1

which is the desired result.
Proof of Theorem 4. We examine both terms together to obtain
kX(βb − β ? )k2 + λkβb − β ? k`1 /H = kX(βb − β ? )k2 + λk(βb − β ? )S0 k`1 /H + λkβbS0c k`1 /H .
Applying Lemmas 10 and 11 we have that, with probability 1 − δ,
p
kX(βb − β ? )k2 + λk(βb − β ? )S0 k`1 /H + λkβbS0c k ≤ 4λk(βb − β ? )S0 k`1 /H ≤ 4λ I0 k(βb − β ? )S0 k.
Applying Definition 3.1, we have that the RHS above is bounded by
√
p
4λ
I0
1
8λ2 I0
?
?
?
2
b
b
b
4λ I0 k(β − β )S0 k ≤ √
kX(β − β)k ≤ kX(β − β)k +
.
2
αN
αN
Note the last inequality follows from the simple bound 4uv ≤ 12 u2 + 8v 2 for any real u and
v (the LHS just completes the square). We therefore have that
2
b 2 + λk(βb − β ? )S k` /H + λkβbS c k` /H ≤ 1 kX(β ? − β)k
b 2 + 8λ I0 .
kX(β ? − β)k
0
1
1
0
2
αN

This implies
2
1
b 2 + λkβb − β ? k` /H ≤ 8λ I0 ,
kX(β ? − β)k
1
2
αN

which proves the claim.
39

