VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS WITH
HIGH-DIMENSIONAL COVARIANCE MATRIX ESTIMATION

arXiv:1707.04145v1 [] 13 Jul 2017

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET
Abstract. In this paper, we propose a novel variable selection approach in the framework
of multivariate linear models taking into account the dependence that may exist between
the responses. It consists in estimating beforehand the covariance matrix Σ of the responses
and to plug this estimator in a Lasso criterion, in order to obtain a sparse estimator of the
coefficient matrix. The properties of our approach are investigated both from a theoretical
and a numerical point of view. More precisely, we give general conditions that the estimators
of the covariance matrix and its inverse have to satisfy in order to recover the positions of
the null and non null entries of the coefficient matrix when the size of Σ is not fixed and
can tend to infinity. We prove that these conditions are satisfied in the particular case of
some Toeplitz matrices. Our approach is implemented in the R package MultiVarSel available
from the Comprehensive R Archive Network (CRAN) and is very attractive since it benefits
from a low computational load. We also assess the performance of our methodology using
synthetic data and compare it with alternative approaches. Our numerical experiments show
that including the estimation of the covariance matrix in the Lasso criterion dramatically
improves the variable selection performance in many cases.

1. Introduction
The multivariate linear model consists in generalizing the classical linear model, in which
a single response is explained by p variables, to the case where the number q of responses
is larger than 1. Such a general modeling can be used in a wide variety of applications
ranging from econometrics (Lütkepohl (2005)) to bioinformatics (Meng et al. (2014)). In the
latter field, for instance, multivariate models have been used to gain insight into complex
biological mechanisms like metabolism or gene regulation. This has been made possible
thanks to recently developed sequencing technologies. For further details, we refer the reader
to Mehmood et al. (2012). However, the downside of such a technological expansion is to
include irrelevant variables in the statistical models. To circumvent this, devising efficient
variable selection approaches in the multivariate setting has become a growing concern.
A first naive approach to deal with the variable selection issue in the multivariate setting
consists in applying classical univariate variable selection strategies to each response separately. Some well-known variable selection methods include the least absolute shrinkage and
selection operator (LASSO) proposed by Tibshirani (1996) and the smoothly clipped absolute
deviation (SCAD) approach devised by Fan and Li (2001). However, such a strategy does
not take into account the dependence that may exist between the different responses.
In this paper, we shall consider the following multivariate linear model:
(1)

Y = XB + E,

Date: July 14, 2017.
1

2

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

where Y = (Yi,j )1≤i≤n,1≤j≤q denotes the n × q random response matrix, X denotes the n × p
design matrix, B denotes a p×q coefficient matrix and E = (Ei,j )1≤i≤n,1≤j≤q denotes the n×q
random error matrix, where n is the sample size. In order to model the potential dependence
that may exist between the columns of E, we shall assume that for each i in {1, . . . , n},
(2)

(Ei,1 , . . . , Ei,q ) ∼ N (0, Σ),

where Σ denotes the covariance matrix of the ith row of the error matrix E. We shall moreover
assume that the different rows of E are independent. With such assumptions, there is some
dependence between the columns of E but not between the rows. Our goal is here to design
a variable selection approach which is able to identify the positions of the null and non null
entries in the sparse matrix B by taking into account the dependence between the columns
of E.
This issue has recently been considered by Lee and Liu (2012) who extended the approach
of Rothman et al. (2010). More precisely, Lee and Liu (2012) proposed three approaches for
dealing with this issue based on penalized maximum likelihood with a weighted `1 regularization. In their first approach B is estimated by using a plug-in estimator of Σ−1 , in the second
one, Σ−1 is estimated by using a plug-in estimator of B and in the third one, Σ−1 and B are
estimated simultaneously. Lee and Liu (2012) also investigate the asymptotic properties of
their methods when the sample size n tends to infinity and the number of rows and columns
q of Σ is fixed.
In this paper, we propose to estimate Σ beforehand and to plug this estimator in a Lasso
criterion, in order to obtain a sparse estimator of B. Hence, our methodology is close to
the first approach of Lee and Liu (2012). However, there are two main differences: The
first one is the asymptotic framework in which our theoretical results are established and the
second one is the strategy that we use for estimating Σ. More precisely, in our asymptotic
framework, q is allowed to depend on n and thus to tend to infinity as n tends to infinity at a
polynomial rate. Moreover, in Lee and Liu (2012), Σ−1 is estimated by using an adaptation of
the Graphical Lasso (GLASSO) proposed by Friedman et al. (2008). This technique has also
been considered by Yuan and Lin (2007), Banerjee et al. (2008) and Rothman et al. (2008).
In this paper, we give general conditions that the estimators of Σ and Σ−1 have to satisfy in
order to be able to recover the support of B that is to find the positions of the null and non
null entries of the matrix B. We prove that when Σ is a particular Toeplitz matrix, namely
the covariance matrix of an AR(1) process, the assumptions of the theorem are satisfied.
Let us now describe more precisely our methodology. We start by “whitening” the observations Y by applying the following transformation to Model (1):
(3)

Y Σ−1/2 = XB Σ−1/2 + E Σ−1/2 .

The goal of such a transformation is to remove the dependence between the columns of Y .
Then, for estimating B, we proceed as follows. Let us observe that (3) can be rewritten as:
(4)

Y = X B + E,

with
(5)

Y = vec(Y Σ−1/2 ), X = (Σ−1/2 )0 ⊗ X, B = vec(B) and E = vec(E Σ−1/2 ),

where vec denotes the vectorization operator and ⊗ the Kronecker product.

VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

3

With Model (4), estimating B is equivalent to estimate B since B = vec(B). Then, for
estimating B, we use the classical LASSO criterion defined as follows for a nonnegative λ:

b
(6)
B(λ)
= ArgminB kY − X Bk22 + λkBk1 ,
where k.k1 and k.k2 denote the classical `1 -norm and `2 -norm, respectively. Inspired by Zhao
and Yu (2006), Theorem 1 established some conditions under which the positions of the null
b
and non null entries of B can be recovered by using B.
In practical situations, the covariance matrix Σ is generally unknown and has thus to be
b denote an estimator of Σ. Then, the estimator Σ
b −1/2 of Σ−1/2 is such that
estimated. Let Σ
b −1 = Σ
b −1/2 (Σ
b −1/2 )0 .
Σ
b −1/2 , (3) becomes
When Σ−1/2 is replaced by Σ
b −1/2 = XB Σ
b −1/2 + E Σ
b −1/2 ,
YΣ

(7)

which can be rewritten as follows:
e
Ye = XeB + E,

(8)
where
(9)

b −1/2 ), Xe = (Σ
b −1/2 )0 ⊗ X, B = vec(B) and Ee = vec(E Σ
b −1/2 ).
Ye = vec(Y Σ

In Model (8), B is estimated by
(10)

n
o
e
B(λ)
= ArgminB kYe − XeBk22 + λkBk1 .

By extending Theorem 1, Theorem 5 gives some conditions on the eigenvalues of Σ−1 and
b and its inverse to Σ and Σ−1 , respectively, under which the
on the convergence rate of Σ
e
positions of the null and non null entries of B can be recovered by using B.
We prove in Section 2.3 that when Σ is a particular Toeplitz matrix, namely the covariance
matrix of an AR(1) process, the assumptions of Theorem 5 are satisfied. This strategy has
been implemented in the R package MultiVarSel, which is available on the Comprehensive
R Archive Network (CRAN), for more general Toeplitz matrices Σ such as the covariance
matrix of ARMA processes or general stationary processes. For a successful application of
this methodology to particular “-omic” data, namely metabolomic data, we refer the reader
to Perrot-Dockès et al. (2017). For a review of the most recent methods for estimating highdimensional covariance matrices, we refer the reader to Pourahmadi (2013).
The paper is organized as follows. Section 2 is devoted to the theoretical results of the
paper. The assumptions under which the positions of the non null and null entries of B can
be recovered are established in Theorem 1 when Σ is known and in Theorem 5 when Σ is
unknown. Section 2.3 studies the specific case of the AR(1) model. We present in Section
3 some numerical experiments in order to support our theoretical results. The proofs of our
main theoretical results are given in Section 4.
2. Theoretical results
2.1. Case where Σ is known. Let us first introduce some notations. Let
(11)

C=

1 0
X X and J = {1 ≤ j ≤ pq, Bj 6= 0},
nq

4

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

where X is defined in (5) and where Bj denotes the jth component of the vector B defined in
(5).
Let also define
1
1
(12)
CJ,J =
(X•,J )0 X•,J and CJ c ,J =
(X•,J c )0 X•,J ,
nq
nq
where X•,J and X•,J c denote the columns of X belonging to the set J defined in (11) and to
its complement J c , respectively.
More generally, for any matrix A, AI,J denotes the partitioned matrix extracted from A
by considering the rows of A belonging to the set I and the columns of A belonging to the
set J, with • indicating all the rows or all the columns.
The following theorem gives some conditions under which the estimator Bb defined in (6) is
sign-consistent as defined by Zhao and Yu (2006), namely,


b = sign(B) → 1, as n → ∞,
P sign(B)
where the sign function maps positive entries to 1, negative entries to -1 and zero to 0.
Theorem 1. Assume that Y = (Y1 , Y2 , . . . , Ynq )0 satisfies Model (4). Assume also that there
exist some positive constants M1 , M2 , M3 and positive numbers c1 , c2 such that 0 < c1 + c2 <
1/2 satisfying:
(A1) for all n ≥ 1, for all j ∈ {1, . . . , pq}, n1 (X•,j )0 X•,j ≤ M1 , where X•,j is the jth column
of X defined in (5),
(A2) for all n ≥ 1, n1 λmin ((X 0 X )J,J ) ≥ M2 , where λmin (A) denotes the smallest eigenvalue
of A,
(A3) |J| = O(q c1 ), where J is defined in (11) and |J| is the cardinality of the set J,
(A4) q c2 minj∈J |Bj | ≥ M3 .
Assume also that the following strong Irrepresentable Condition holds:
(IC) There exists a positive constant vector η such that
(X 0 X )J c ,J ((X 0 X )J,J )−1 sign(BJ ) ≤ 1 − η,
where 1 is a (pq − |J|) vector of 1 and the inequality holds element-wise.
Then, for all λ that satisfies




1
λ
λ
= o q −(c1 +c2 ) , as n → ∞,
(L) q = qn = o n 2(c1 +c2 ) , √ → ∞ and
n
n
we have


b
P sign(B(λ))
= sign(B) → 1, as n → ∞,
b
where B(λ)
is defined by (6).
Remark 1. Observe that if c1 + c2 < (2k)−1 , for some positive k, then the first condition of
(L) becomes q = o(nk ). Hence for large values of k, the size q of Σ is much larger than n.
The proof of Theorem 1 is given in Section 4. It is based on Proposition 2 which is an
adaptation to the multivariate case of Proposition 1 in Zhao and Yu (2006).
b
Proposition 2. Let B(λ)
be defined by (6). Then


b
P sign(B(λ))
= sign(B) ≥ P(An ∩ Bn ),

VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

5

where

(13)

An =

−1

(CJ,J )

WJ <

√



λ
−1
nq |BJ | −
|(CJ,J ) sign(BJ )|
2nq

and

(14)

Bn =

−1

CJ c ,J (CJ,J )

WJ − WJ c



λ
−1
,
1 − CJ c ,J (CJ,J ) sign(BJ )
≤ √
2 nq

√
with W = X 0 E/ nq. In (13) and (14), CJ,J and CJ c ,J are defined in (12) and WJ and
WJ c denote the components of W being in J and J c , respectively. Note that the previous
inequalities hold element-wise.
The proof of Proposition 2 is given in Section 4.
We give in the following proposition which is proved in Section 4 some conditions on X
and Σ under which Assumptions (A1) and (A2) of Theorem 1 hold.
Proposition 3. If there exist some positive constants M10 , M20 , m1 , m2 such that, for all
n ≥ 1,
(C1) for all j ∈ {1, . . . , p}, n1 (X•,j )0 X•,j ≤ M10 ,
(C2) n1 λmin (X 0 X) ≥ M20 ,
(C3) λmax (Σ−1 ) ≤ m1 ,
(C4) λmin (Σ−1 ) ≥ m2 ,
then Assumptions (A1) and (A2) of Theorem 1 are satisfied.
Remark 2. Observe that (C1) and (C2) hold in the case where the columns of the matrix X
are orthogonal.
We give in Proposition 6 in Section 2.3 some conditions under which Condition (IC) holds
in the specific case where Σ is the covariance matrix of an AR(1) process.
2.2. Case where Σ is unknown. Similarly as in (11) and (12), we introduce the following
notations:
e = 1 Xe0 Xe
(15)
C
nq
and
eJ,J = 1 (Xe•,J )0 Xe•,J and C
eJ c ,J = 1 (Xe•,J c )0 Xe•,J ,
C
nq
nq

(16)

where Xe•,J and Xe•,J c denote the columns of Xe belonging to the set J defined in (11) and to
its complement J c , respectively.
A straightforward extension of Proposition 2 leads to the following proposition for Model (8).
e
Proposition 4. Let B(λ)
be defined by (10). Then


e
fn ∩ B
fn ),
P sign(B(λ))
= sign(B) ≥ P(A
where

(17)

fn =
A



λ e
√
−1 f
−1
e
(CJ,J ) WJ < nq |BJ | −
|(CJ,J ) sign(BJ )|
2nq

6

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

and

(18)

fn =
B

eJ c ,J (C
eJ,J )−1 W
fJ
C

fJ c
−W


λ 
−1
eJ c ,J (C
eJ,J ) sign(BJ )
,
≤ √
1− C
2 nq

e √nq. In (17) and (18), C
eJ,J and C
eJ c ,J are defined in (16) and W
fJ and
with W = Xe0 E/
fJ c denote the components of W
f being in J and J c , respectively. Note that the previous
W
inequalities hold element-wise.
The following theorem extends Theorem 1 to the case where Σ is unknown and gives
some conditions under which the estimator Be defined in (10) is sign-consistent. The proof of
Theorem 5 is given in Section 4 and is based on Proposition 4.
Theorem 5. Assume that Assumptions (A1), (A2), (A3), (A4), (IC) and (L) of Theorem 1
hold. Assume also that, there exist some positive constants M4 , M5 , M6 and M7 , such that
for all n ≥ 1,
(A5) k(X 0 X)/nk∞ ≤ M4 ,
(A6) λmin ((X 0 X)/n) ≥ M5 ,
(A7) λmax (Σ−1 ) ≤ M6 ,
(A8) λmin (Σ−1 ) ≥ M7 .
Suppose also that
b −1 k∞ = OP ((nq)−1/2 ), as n tends to infinity,
(A9) kΣ−1 − Σ
b = OP ((nq)−1/2 ), as n tends to infinity.
(A10) ρ(Σ − Σ)
e
Let B(λ)
be defined by (10), then


e
P sign(B(λ))
= sign(B) → 1, as n → ∞.
In the previous assumptions, λmax (A), λmin (A), ρ(A) and kAk∞ denote the largest eigenvalue,
the smallest eigenvalue, the spectral radius and the infinite norm (induced by the associated
vector norm) of the matrix A.
Remark 3. Observe that Assumptions (A5) and (A6) hold in the case where the columns of
the matrix X are orthogonal. Note also that (A7) and (A8) are the same as (C3) and (C4)
in Proposition 3.
In order to estimate Σ, we propose the following strategy:
• Fitting a classical linear model to each column of the matrix Y in order to have access
b of the random error matrix E. It is possible since p is assumed to
to an estimation E
be fixed and smaller than n.
b by assuming that Σ has a particular structure, Toeplitz for
• Estimating Σ from E
instance.
b defined in the first step is such that:
More precisely, E

b = IdRn − X(X 0 X)−1 X 0 E =: ΠE,
(19)
E
which implies that
(20)
where E is defined in (5).

b = [IdRq ⊗ Π] E,
Eb = vec(E)

VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

7

We prove in Proposition 7 below that our strategy for estimating Σ provides an estimator
satisfying the assumptions of Theorem 5 in the case where (E1,t )t , (E2,t )t , ..., (En,t )t are
assumed to be independent AR(1) processes.
2.3. The AR(1) case.
2.3.1. Sufficient conditions for Assumption (IC) of Theorem 1. The following proposition
gives some conditions under which the strong Irrepresentable Condition (IC) of Theorem 1
holds.
Proposition 6. Assume that (E1,t )t , (E2,t )t , ..., (En,t )t in Model (1) are independent AR(1)
processes satisfying:
Ei,t − φ1 Ei,t−1 = Zi,t , ∀i ∈ {1, . . . , n},
where the Zi,t ’s are zero-mean i.i.d. Gaussian random variables with variance σ 2 and |φ1 | < 1.
Assume also that X defined in (1) is such that X 0 X = νIdRp , where ν is a positive constant.
Moreover, suppose that if j ∈ J, then j > p and j < pq − p. Suppose also that for all j, j − p
or j + p is not in J. Then, the strong Irrepresentable Condition (IC) of Theorem 1 holds.
The proof of Proposition 6 is given in Section 4.
2.3.2. Sufficient conditions for Assumptions (A7), (A8), (A9) and (A10) of Theorem 5. The
following proposition establishes that in the particular case where the (E1,t )t , (E2,t )t , ...,
(En,t )t are independent AR(1) processes, our strategy for estimating Σ provides an estimator
satisfying the assumptions of Theorem 5.
Proposition 7. Assume that (E1,t )t , (E2,t )t , ..., (En,t )t in Model (1) are independent AR(1)
processes satisfying:
Ei,t − φ1 Ei,t−1 = Zi,t , ∀i ∈ {1, . . . , n},
where the Zi,t ’s are zero-mean i.i.d. Gaussian random variables with variance σ 2 and |φ1 | < 1.
Let


1
φb1 φb21 . . . φbq−1
1

 φb
1 φb1 . . . φbq−2
1


1
1 
.. 
..
..
..
..
b

.
.
.
Σ=
. 
,
 .
1 − φb21  .
.
.
.
.
..
..
.. 
..

 ..
q−1
... ... ...
1
φb1
where
Pn Pq b b
`=2 Ei,` Ei,`−1
b
(21)
φ1 = i=1
,
Pn Pq−1 b 2
E
i=1

`=1

i,`

b = (E
bi,` )1≤i≤n,1≤`≤q is defined in (19). Then, Assumptions (A7), (A8), (A9) and
where E
(A10) of Theorem 5 are valid.
The proof of Proposition 7 is given in Section 4. It is based on the following lemma.
Lemma 8. Assume that (E1,t )t , (E2,t )t , ..., (En,t )t in Model (1) are independent AR(1)
processes satisfying:
Ei,t − φ1 Ei,t−1 = Zi,t , ∀i ∈ {1, . . . , n},

8

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

where the Zi,t ’s are zero-mean i.i.d. Gaussian random variables with variance σ 2 and |φ1 | < 1.
Let
Pn Pq b b
`=2 Ei,` Ei,`−1
b
,
φ1 = i=1
Pn Pq−1 b 2
i=1
`=1 Ei,`
b = (E
bi,` )1≤i≤n,1≤`≤q is defined in (19). Then,
where E
√
nqn (φb1 − φ1 ) = Op (1), as n → ∞.
Lemma 8 is proved in Section 4. Its proof is based on Lemma 10 in Section 5.
3. Numerical experiments
The goal of this section is twofold: i) to provide sanity checks for our theoretical results in
a well-controlled framework; and ii) to investigate the robustness of our estimator to some
violations of the assumptions of our theoretical results. The latter may reveal a broader scope
of applicability for our method than the one guaranteed by the theoretical results.
We investigate i) in the AR(1) framework presented in Section 2.3. Indeed, all assumptions
made in Theorems 1 and 5 can be specified with well-controllable simulation parameters in
the AR(1) case with balanced design matrix X.
Point ii) aims to explore the limitations of our theoretical framework and assess its robustness. To this end, we propose two numerical studies relaxing some of the assumptions of our
theorems: first, we study the effect of an unbalanced design – which violates the sufficient condition of the irrepresentability condition (IC) given in Proposition 6 – on the sign-consistency;
and second, we study the effect of other types of dependence than an AR(1).
In all experiments, the performance are assessed in terms of sign-consistency. In other
words, we evaluate the probability for the sign of various estimators to be equal to sign(B).
We compare the performance of three different estimators:
• Bb defined in (6), which corresponds to the LASSO criterion applied to the data
whitened with the true covariance matrix Σ; we call this estimator oracle. Its theoretical properties are established in Theorem 1.
• Be defined in (10), which corresponds to the LASSO criterion applied to the data
b we refer to this estimator as
whitened with an estimator of the covariance matrix Σ;
whitened-lasso. Its theoretical properties are established in Theorem 5.
• the LASSO criterion applied to the raw data, which we call raw-lasso hereafter.
Its theoretical properties are established only in the univariate case in Alquier and
Doukhan (2011).
3.1. AR(1) dependence structure with balanced one-way ANOVA. In this section,
we consider Model (1) where X is the design matrix of a one-way ANOVA with two balanced
groups. Each row of the random error matrix E is distributed as a centered Gaussian random
vector as in Equation (2) where the matrix Σ is the covariance matrix of an AR(1) process
defined in Section 2.3.
In this setting, Assumptions (A1), (A2) and Condition (IC) of Theorem 1 are satisfied, see
Propositions 3 and 6. The three remaining assumptions (A3), (A4) and (L) are related to
more practical quantities: (A3) controls the sparsity level of the problem, involving c1 ; (A4)
basically controls the signal-to-noise ratio, involving c2 and (L) links the sample size n, q and
the two constants c1 , c2 , so that an appropriate range of penalty λ exists for having a large
probability of support recovery. This latter assumption is used in our experiments to tune

VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

9

the difficulty of the support recovery as follows: we consider different values of n, q, c1 , c2
and we choose a sparsity level |J| and a minimal magnitude in B such that Assumptions (A3)
and (A4) are fulfilled. Hence, the problem difficulty is essentially driven by the validity of
Assumption (L) where q = o(nk ) with c1 + c2 = 1/2k, and so by the relationship between n,
q and k.
We consider a large range of sample sizes n varying from 10 to 1000 and three different
values for q in {10, 50, 1000}. The constants c1 , c2 are chosen such that c1 + c2 = 1/2k with
c1 = c2 and k in {1, 2, 4}. Additional values of c1 and c2 have also been considered and
the corresponding results are available upon request. Finally, we consider two values for the
parameter φ1 appearing in the definition of the AR(1) process: φ1 ∈ {0.5, 0.95}.
Note that in this AR(1) setting with the estimator φb1 of φ1 defined in (21) , all the
assumptions of Theorem 5 are fulfilled, see Proposition 7.
The frequencies of support recovery for the three estimators averaged over 1000 replications
is displayed in Figure 1.

k= 1
1.00
●
●

●
●

●
● ●
● ●●

k= 2
●
●

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

● ●●
● ●●

k= 4
●
●

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

● ●
● ●

● ●●

●

● ●●

●

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

●

0.75

●

q= 10

●
●

●

0.50

● ●
● ●

●
●

●
●

0.25

φ1

0.00
1.00

●

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

●
● ●

●●

0.75

●

●

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

● ●
●

0.5

● ● ● ● ●●●●

0.95

●

●

●

●
●

●
●

●

●

●

●
●●
● ●
●

●

●

●
●

●

0.25

●
●

q= 50

0.50

●
● ●●
●

●

●

●

oracle

●

0.00

●
●

●

raw−lasso

1.00

●
●

●
●
●
● ●
●

●

●

●
●

●
●

●

0.00

10

● ●●

100

●
●

●

●

●

●

●
● ●

●
●

●

●

whithened−lasso

● ● ● ● ●●●●

●

●
● ●

●

●

●

●

●

●

●

0.25

●●
● ●
●
●

q= 1000

●

0.50

●

●

●

0.75

● ● ● ● ●●●●

●

1000 10

●
●

100

1000 10

100

1000

n
Figure 1. Frequencies of support recovery in a multivariate one-way ANOVA
model with two balanced groups and an AR(1) dependence.

10

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

We observe from Figure 1 that whitened-lasso and oracle have similar performance
since φ1 is well estimated. These two approaches always exhibit better performance than
raw-lasso, especially when φ1 = 0.95. In this case, the sample size n required to reach the
same performance is indeed ten time larger for raw-lasso than for oracle and whitened-lasso.
Finally, the performance of all estimators are altered when n is too small, especially in
situations where the signal to noise ratio (SNR) is small and the signal is not sparse enough,
these two characteristics corresponding to small values of k.
3.2. Robustness to unbalanced designs and correlated features. The goal of this
section is to study some particular design matrices X in Model (1) that may lead to violation
of the Irrepresentability Condition (IC).
To this end, we consider the multivariate linear model (1) with the same AR(1) dependence
as the one considered in Section 3.1. Then, two different matrices X are considered: First, an
one-way ANOVA model with two unbalanced groups with respective sizes n1 and n2 such that
n1 + n2 = n; and second, a multiple regression model with p correlated Gaussian predictors
such that the rows of X are i.i.d. N (0, ΣX ).
For the one-way ANOVA, violation of (IC) may occur when r = n1 /n is too different from
1/2, as stated in Proposition 6. For the regression model, we choose for ΣX a 9 × 9 matrix
X
(p = 9) such that ΣX
i,i = 1, Σi,j = ρ, when i 6= j. The other simulation parameters are fixed
as in Section 3.1.
We report in Figure 2 the results for the case where q = 1000 and k = 2 both for unbalanced
one-way ANOVA (top panels) and regression with correlated predictors (bottom panels).
For the one-way ANOVA, r varies in {0.4, 0.2, 0.1}. For the regression case, ρ varies in
{0.2, 0.6, 0.9}. In both cases, the gray lines correspond to the ideal situation (that is, either
unbalanced or uncorrelated) denoted Ideal in the legend of Figure 2. The probability of
support recovery is estimated over 1000 runs.
From this figure, we note that correlated features or unbalanced designs deteriorate the
support recovery of all estimators. This was expected for these LASSO-based methods which
all suffer from the violation of the irrepresentability condition (IC). However, we also note
that whitened-lasso and oracle have similar performance, which means that the estimation
of Σ is not altered, and that whitening always improves the support recovery.
3.3. Robustness to more general autoregressive processes. In this section, we consider
the case where X is the design matrix of a one-way ANOVA with two balanced groups
and where Σ is the covariance matrix of an AR(m) process with m in {5, 10}. Figure 3
displays the performance of the different estimators when q = 500. Here, for computing
b in whitened-lasso, the parameters φ1 , . . . , φm of the AR(m) process are estimated as
Σ
b defined in (19) the estimations
follows. They are obtained by averaging over the n rows of E
(i)
(i)
b by using standard estimation approaches for AR
φb1 , . . . , φbm obtained for the ith row of E
processes described in Brockwell and Davis (1990). As previously, we observe from this figure
that whitened-lasso and oracle have better performance than raw-lasso.
4. Proofs
Proof of Proposition 2. For a fixed nonnegative λ, by (6),

b
Bb = B(λ)
= ArgminB kY − X Bk22 + λkBk1 .

VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

r = 0.4
1.00

●

●
●

●

●

●

r = 0.2
●

●

●

●

●

●

●

●

r = 0.1

●

●

●

●

●
●
●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

0.50

●

●
●

0.25
0.00

●
●
●
●

●
●

●

●

●
●
●

●
●●●
●

●

●

●

0

250

500

●

●

●
●
●
●

●
●
●

●
●

●

●
●●
●
●
●

q = 1000

0.75

11

750 1000 0

250

500

750 1000 0

250

500

750 1000

n
ρ = 0.2
1.00

●

●
●
●
●●

●

●

●

●

ρ = 0.6
●

●

●

●

●
●

●●

0.50
0.25

●

●

●

●

●

●
●

●

●
●

●
●

●

●

●
●

●

●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●●

●

●

●

●●●
●
●●

0.00
0

250

500

750 1000 0

250

500

●

oracle

●

raw−lasso

●

whithened−lasso

●

oracle
raw−lasso

●

●
●●

●

Ideal

ρ = 0.8

q = 1000

0.75

●

●

750 1000 0

●

●

250

●

●

●

500

●

●

●

whithened−lasso

●

750 1000

n
Figure 2. Frequencies of support recovery in general linear models with unbalanced designs: one-way ANOVA and regression.

m= 5
1.00

●

m= 10
●

●

●

●
●

0.75

●

●

●

●

q= 500

●

0.50

●

●

●

●
●

●
●

0.25

●

0.00

0

●

100

200

300

400

500 0

raw−lasso
whithened−lasso

●
●
●

oracle

100

200

300

400

500

n
Figure 3. Frequencies of support recovery in one-way ANOVA with AR(m)
covariance matrix.
Denoting u
b = Bb − B, we get
b 2 + λkBk
b 1 = kX B + E − X Bk
b 2 + λkb
kY − X Bk
u + Bk1 = kE − X u
bk22 + λkb
u + Bk1
2
2
= kEk22 − 2b
u0 X 0 E + u
b0 X 0 X u
b + λkb
u + Bk1 .

12

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

Thus,
u
b = Argminu V (u),
where

√
√
√
V (u) = −2( nqu)0 W + ( nqu)0 C( nqu) + λku + Bk1 .

Since the first derivative of V with respect to u is equal to

√
√
2 nq C( nqu) − W + λ sign(u + B),
u
b satisfies
λ
λ
√
CJ,J ( nqb
uJ ) − WJ = − √ sign(b
uJ + BJ ) = − √ sign(BbJ ), if b
uJ + BJ = BbJ 6= 0
2 nq
2 nq
and

λ
√
CJ c ,J ( nqb
uJ ) − WJ c ≤ √ .
2 nq
Note that, if |b
uJ | < |BJ |, then BbJ 6= 0 and sign(BbJ ) = sign(BJ ).
Let us now prove that when An and Bn , defined in (13) and (14), are satisfied then there
exists u
b satisfying:
(22)
(23)
(24)

λ
√
CJ,J ( nqb
uJ ) − WJ = − √ sign(BJ ),
2 nq
|b
uJ | < |BJ |,
λ
√
CJ c ,J ( nqb
uJ ) − WJ c ≤ √ .
2 nq

Note that An implies:
(25) 



λ
λ
√
√
−1
−1
−1
(CJ,J ) sign(BJ ) < (CJ,J ) WJ < nq |BJ | +
(CJ,J ) sign(BJ ) .
nq −|BJ | +
2nq
2nq
By denoting
(26)

1
λ
u
bJ = √ (CJ,J )−1 WJ −
(CJ,J )−1 sign(BJ ),
nq
2nq

we obtain from (25) that (22) and (23) hold. Note that Bn implies:

λ
− √
1 − CJ c ,J (CJ,J )−1 sign(BJ )
2 nq

λ
≤ CJ c ,J (CJ,J )−1 WJ − WJ c ≤ √
1 + CJ c ,J (CJ,J )−1 sign(BJ ) .
2 nq
Hence,


λ
λ
−1
−1
CJ c ,J (CJ,J ) WJ − √ (CJ,J ) sign(BJ ) − WJ c ≤ √ ,
2 nq
2 nq
which is (24) by (26). This concludes the proof.
Proof of Theorem 1. By Proposition 2,


b
P sign(B(λ))
= sign(B) ≥ P(An ∩ Bn ) = 1 − P(Acn ∪ Bnc ) ≥ 1 − P(Acn ) − P(Bnc ),



VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

13

where An and Bn are defined in (13) and (14). It is thus enough to prove that P(Acn ) and
P(Bnc ) tend to zero as n tends to infinity.
By definition of An ,



λ
√
c
−1
−1
P(An ) = P (CJ,J ) WJ ≥ nq |BJ | −
|(CJ,J ) sign(BJ )|
2nq



λ
√
(27)
|bj |
,
≤ sup P |ξj | ≥ nq |Bj | −
2nq
j∈J
where
1
ξ = (ξj )j∈J = (CJ,J )−1 WJ = √ (CJ,J )−1 (X•,J )0 E =: HA E,
nq
and
b = (bj )j∈J = (CJ,J )−1 sign(BJ ).
By definition of Bn and (IC),



λ
−1
c
−1
1 − CJ c ,J (CJ,J ) sign(BJ )
P(Bn ) = P CJ c ,J (CJ,J ) WJ − WJ c > √
2 nq


λ
−1
c
c
≤ P CJ ,J (CJ,J ) WJ − WJ > √ η
2 nq


λ
(28)
≤ sup P |ζj | > √ η ,
c
2
nq
j∈J
where

1
ζ = (ζj )j∈J c = CJ c ,J (CJ,J )−1 WJ − WJ c = √
CJ c ,J (CJ,J )−1 (X•,J )0 − (X•,J c )0 E =: HB E.
nq
Note that, for all j in J,
|bj | ≤

X


1/2
X
p
p
|bj | ≤ |J| 
b2j  = |J|kbk2 .

j∈J

j∈J

Moreover,
kbk2 = k(CJ,J )−1 sign(BJ )k2 ≤ k(CJ,J )−1 k2

p
p
|J| := λmax ((CJ,J )−1 ) |J|,

where λmax (A) denotes the largest eigenvalue of the matrix A. Observe that
(29)

λmax ((CJ,J )−1 ) =

1
q
q
=
≤
,
0
λmin (CJ,J )
λmin ((X X )J,J )/n
M2

by Assumption (A2) of Theorem 1. Thus, for all j in J,
(30)

|bj | ≤

q|J|
.
M2

By Assumption (A4) of Theorem 1, we get thus that for all j in J,
(31) 






λ
λ
λq|J|
√
√
√
−1
−c2
nq |Bj | −
(CJ,J ) sign(BJ ) j = nq |Bj | −
|bj | ≥ nq M3 q
−
.
2nq
2nq
2nqM2

14

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

Thus,
P(Acn )

(32)


≤ sup P |ξj | ≥

√

j∈J



λq|J|
−c2
.
nq M3 q
−
2nqM2

Since E is a centered Gaussian random vector having a covariance matrix equal to identity,
ξ = HA E is a centered Gaussian random vector with a covariance matrix equal to:
0
HA HA
=

1
(CJ,J )−1 (X•,J )0 X•,J (CJ,J )−1 = (CJ,J )−1 .
nq

Hence, by (29), we get that for all j in J,
Var(ξj ) = (CJ,J )−1


jj

−1
≤ λmax (CJ,J
)≤

q
.
M2

Thus,
√






λq|J|
M2
λq|J|
√
−c2
−c2 √
P |ξj | ≥ nq M3 q
−
≤ P |Z| ≥ √
M3 q
nq − √
,
2nqM2
q
2 nqM2
where Z is a standard Gaussian random variable. By Chernoff inequality, we thus obtain that
for all j in J,




2 !
λq|J|
M
λq|J|
√
√
2
P |ξj | ≥ nq M3 q −c2 −
M3 q −c2 nq − √
≤ 2 exp −
.
2nqM2
2q
2 nqM2
By Assumption (A3) of Theorem 1, we get that under the last condition of (L),
λq|J|
√ 
= o q −c2 nq , as n → ∞.
√
nq

(33)
Thus,

P(Acn ) → 0, as n → ∞.

(34)

Let us now bound P(Bnc ). Observe that ζ = HB E is a centered Gaussian random vector
with a covariance matrix equal to:
1
0
−1
(CJ c ,J (CJ,J )−1 (X•,J )0 − X•,J
CJ,J c − X•,J c )
c )(X•,J (CJ,J )
nq

1
= CJ c ,J c − CJ c ,J (CJ,J )−1 CJ,J c =
(X•,J c )0 IdRnq − X•,J ((X•,J )0 X•,J )−1 (X•,J )0 X•,J c
nq


1
(X•,J c )0 IdRnq − ΠIm(X•,J ) X•,J c ,
=
nq

0
HB HB
=

where ΠIm(X•,J ) denotes the orthogonal projection onto the column space of X•,J . Note that,
for all j in J c ,



1 
Var(ζj ) =
(X•,J c )0 IdRnq − ΠIm(X•,J ) X•,J c
nq
jj



1
1
(X•,J c )0 X•,J c jj −
=
(X•,J c )0 ΠIm(X•,J ) X•,J c
nq
nq
jj

M
1
1
(X•,J c )0 X•,J c jj ≤
,
≤
nq
q

VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

15

where the inequalities come from Lemma 9 and Assumption (A1) of Theorem 1. Thus, for
all j in J c ,




√
λ q
λ
P |ζj | > √ η ≤ P |Z| > √ √ η ,
2 nq
2 M1 nq
where Z is a standard Gaussian random variable. By Chernoff inequality, for all j in J c ,
(
2 )



λ
λ
1
√ √ η
.
P |ζj | > √ η ≤ 2 exp −
2 nq
2 2 M1 n
Hence, under the following assumption
λ
√ → ∞,
n
which is the second condition of (L),
(35)

P(Bnc ) → 0, as n → ∞.


Proof of Proposition 3. Let us first prove that (C1) and (C3) imply (A1). For j ∈ {1, . . . , pq},
by considering the Euclidian division of j − 1 by p given by (j − 1) = pkj + rj , we observe
that
(X•,j )0 X•,j = (((Σ−1/2 )0 ⊗ X)•,j )0 ((Σ−1/2 )0 ⊗ X)•,j
= ((Σ−1/2 ) ⊗ X 0 )j,• )((Σ−1/2 )0 ⊗ X)•,j
= ((Σ−1/2 )kj +1,• ⊗ (X•,rj +1 )0 )(((Σ−1/2 )•,kj +1 )0 ⊗ X•,rj +1 )
= (Σ−1/2 )kj +1,• ((Σ−1/2 )•,kj +1 )0 ⊗ (X•,rj +1 )0 X•,rj +1
= (Σ−1 )kj +1,kj +1 ⊗ (X•,rj +1 )0 X•,rj +1
= (Σ−1 )kj +1,kj +1 (X•,rj +1 )0 X•,rj +1 .
Hence, using (C1), we get that for all j in {1, . . . , pq},
1
(X•,j )0 X•,j ≤ M10 (Σ−1 )kj +1,kj +1 ≤ M10
sup
((Σ−1 )k+1,k+1 )
n
k∈{0,...,q−1}
≤ M10 λmax (Σ−1 ) ≤ M10 m1 ,
where the last inequality comes from (C3), which gives (A1).
Let us now prove that (C2) and (C4) imply (A2). Note that
(X 0 X )J,J = (((Σ−1/2 )0 ⊗ X)0 ((Σ−1/2 )0 ⊗ X))J,J
= (Σ−1/2 (Σ−1/2 )0 ⊗ X 0 X)J,J
= (Σ−1 ⊗ X 0 X)J,J .
Then, by Theorem 4.3.15 of Horn and Johnson (1986),
λmin ((X 0 X )J,J ) = λmin ((Σ−1 ⊗ X 0 X)J,J )
≥ λmin (Σ−1 ⊗ X 0 X)
(36)

= λmin (X 0 X)λmin (Σ−1 ).

16

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

Finally, by using Conditions (C2) and (C4), we obtain
1
1
λmin (X 0 X )J,J ≥ λmin (X 0 X)λmin (Σ−1 ) ≥ M20 m2 ,
n
n
which gives (A2).



Proof of Theorem 5. By Proposition 4,


e
en ∩ B
en ) = 1 − P(A
ec ∪ B
e c ) ≥ 1 − P(A
ec ) − P(B
e c ),
P sign(B(λ))
= sign(B) ≥ P(A
n
n
n
n
en and B
en are defined in (17) and (18). By definition of A
en , we get
where A



λ e
√
c
−1 f
−1
e
e
P(An ) = P
(CJ,J ) WJ ≥ nq |BJ | −
|(CJ,J ) sign(BJ )|
.
2nq
Observing that


eJ,J )−1 W
fJ = (CJ,J )−1 WJ + (CJ,J )−1 W
fJ − WJ
(C


eJ,J )−1 − (CJ,J )−1 WJ
+ (C



fJ − WJ ,
eJ,J )−1 − (CJ,J )−1 W
+ (C


eJ,J )−1 sign(BJ ) = (CJ,J )−1 sign(BJ ) + (C
eJ,J )−1 − (CJ,J )−1 sign(BJ ),
(C
and using the triangle inequality, we obtain that


√ 
nq
λ
−1
−1
c
e
(CJ,J ) WJ ≥
(CJ,J ) sign(BJ )
P(An ) ≤ P
|BJ | −
5
2nq


√ 


nq
λ
−1 f
−1
WJ − W J ≥
+P
(CJ,J )
(CJ,J ) sign(BJ )
|BJ | −
5
2nq
 


√

nq
λ
−1
−1
−1
eJ,J ) − (CJ,J )
WJ ≥
(C
+P
|BJ | −
(CJ,J ) sign(BJ )
5
2nq
 

√ 


nq
λ
−1
−1
−1
e
f
WJ − WJ ≥
+P
(CJ,J ) − (CJ,J )
|BJ | −
(CJ,J ) sign(BJ )
5
2nq


√ 


nq
λ
λ
−1
−1
−1
e
(
C
)
−
(C
)
sign(B
)
≥
|B
|
−
(C
)
sign(B
)
.
+P
√
J,J
J,J
J
J
J,J
J
2 nq
5
2nq
(37)
The first term in the r.h.s of (37) tends to 0 by the definition of Acn and (34). By (31), the
last term of (37) satisfies, for all j ∈ J:




2nq
λ
−1
−1
−1
eJ,J ) − (CJ,J )
P
(C
sign(BJ ) ≥
|BJ | −
(CJ,J ) sign(BJ )
5λ
2nq
 




2nq
λq|J|
−1
−1
−c2
e
(CJ,J ) − (CJ,J )
≥
≤P
sign(BJ )
M3 q
−
.
5λ
2nqM2
j

VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

17

eJ,J )−1 − (CJ,J )−1 and s = sign(BJ ) then for all j in J:
Let U = (C
|(U s)j | =

(38)

X

Ujk sk ≤

p
|J|kU k2 .

k∈J

We focus on
eJ,J )−1 − (CJ,J )−1 k2 = k(C
eJ,J )−1 (CJ,J − C
eJ,J )(CJ,J )−1 k2 ≤ k(C
eJ,J )−1 k2 kCJ,J − C
eJ,J k2 k(CJ,J )−1 k2
k(C
≤

eJ,J )
eJ,J )
ρ(CJ,J − C
ρ(CJ,J − C
≤
,
eJ,J )λmin (CJJ )
eJ,J )(M2 /q)
λmin (C
λmin (C

where the last inequality comes from Assumption (A2) of Theorem 1, which gives that
q
(39)
k(CJ,J )−1 k2 ≤
.
M2
Using Theorem 4.3.15 of Horn and Johnson (1986), we get
eJ,J )−1 − (CJ,J )−1 k2 ≤
k(C

e
qρ(C − C)
.
e 2
λmin (C)M

e given in (11) and (15), respectively, we get
By definition of C and C
(40)

C=

0
b −1
Σ−1 ⊗ (X 0 X)
e = Σ ⊗ (X X) .
and C
nq
nq

By using that the eigenvalues of the Kronecker product of two matrices is equal to the product
of the eigenvalues of the two matrices, we obtain
b −1 )λmax ((X 0 X)/n)q
b −1 )λmax (Σ)λ
b max ((X 0 X)/n)q
ρ(Σ−1 − Σ
ρ(Σ−1 − Σ
≤
b −1 )λmin ((X 0 X)/n)M2
λmin ((X 0 X)/n)M2
λmin (Σ


b −1 ) ρ(Σ
b − Σ) + λmax (Σ) λmax ((X 0 X)/n)q
ρ(Σ−1 − Σ
≤
,
λmin ((X 0 X)/n)M2

eJ,J )−1 − (CJ,J )−1 k2 ≤
k(C

(41)

where the last inequality follows from Theorem 4.3.1 of Horn and Johnson (1986). Thus, by
Assumptions (A5), (A6), (A8), (A9) and (A10), we get that
(42)

eJ,J )−1 − (CJ,J )−1 k2 = OP (q(nq)−1/2 ), as n → ∞.
k(C

Hence, by (38), we get for all j in J that
 


eJ,J )−1 − (CJ,J )−1 sign(BJ )
P
(C



λq|J|
2nq
M3 q −c2 −
≥
5λ
2nqM2
j



√
p
2 nq
λq|J|
−1
−1
−c2 √
e
≤P
|J| k(CJ,J ) − (CJ,J ) k2 ≥
M3 q
nq − √
.
5λ
2 nqM2

By (33), (42) and (A3), it is enough to prove that

nq −c2 
P q c1 /2 q(nq)−1/2 ≥
q
→ 0, as n → ∞.
λ
By the last condition of (L),
nq −c2
λ q
→ ∞, as n → ∞,
q 1+c1

18

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

and the result follows since n tends to infinity. Hence, the last term of (37) tends to zero as
n tends to infinity.
Let us now study the second term in the r.h.s of (37).



 
1  e0 e
fJ − WJ = √1
W
Xe0 Ee − X 0 E J = √
X E − X 0E
nq
nq
J
J





i
1 h b −1/2
b −1/2 )0 ⊗ IdRn Vec(E) − Σ−1/2 ⊗ X 0 (Σ−1/2 )0 ⊗ IdRn Vec(E)
=√
Σ
⊗ X 0 (Σ
nq
J
(43)

o
i
1 hn b −1
d
=√
Σ − Σ−1 ⊗ X 0 Vec(E) = AZ,
nq
J
where Z is a centered Gaussian random vector having a covariance matrix equal to identity
and

on
oi
1 hn b −1
(44)
A= √
.
Σ − Σ−1 ⊗ X 0
(Σ1/2 )0 ⊗ IdRn
nq
J,•
By Cauchy-Schwarz inequality, we get for all K × nq matrix B, and all nq × 1 vector U that
for all k in {1, . . . , K},
|(BU )k | =

(45)

nq
X

Bk,` U` ≤ kBk2 kU k2 .

`=1

Thus, for all j in J, for all γ in R and all |J| × |J| matrix D,
 




fJ − WJ
(46)
P
D W
≥ γ = P (DAZ)j ≥ γ ≤ P (kDk2 kAk2 kZk2 ≥ γ) ,
j

where A is defined in (44) and Z is a centered Gaussian random vector having a covariance
matrix equal to identity. Hence, for all j in J,


√ 



nq
λ
−1 f
−1
WJ − WJ
P
(CJ,J )
≥
|Bj | −
(CJ,J ) sign(BJ ) j
5
2nq
j



√
nq
λq|J|
−c2
−1
M3 q
−
≤ P k(CJ,J ) k2 kAk2 kZk2 ≥
.
5
2nqM2
Let us bound kAk2 . Observe that

on
oi
hn
b −1 − Σ−1 ⊗ X 0
Σ
(Σ1/2 )0 ⊗ Id
=ρ

h


J,• 2

 

i
b −1 − Σ−1 Σ Σ
b −1 − Σ−1 ⊗ (X 0 X)
Σ




1/2

1/2
J,J

b −1 − Σ−1 Σ Σ
b −1 − Σ−1
Σ
λmax X 0 X



b −1 − Σ−1 λmax (Σ)1/2 λmax X 0 X 1/2 ,
≤ρ Σ

≤ρ

(47)

1/2

VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

19

where the first inequality comes from Theorem 4.3.15 of Horn and Johnson (1986). Hence,
by (A5), (A8) and (A9)

on
oi
hn
1
b −1 − Σ−1 ⊗ X 0
= OP (q −1/2 (nq)−1/2 ).
Σ
(Σ1/2 )0 ⊗ Id
(48)
kAk2 = √
nq
J,• 2
By (33), (39) and (48), it is enough to prove that
!
nq
X
P
Zk2 ≥ nq n q −2c2 → 0, as n → ∞.
k=1

The result follows from the Markov inequality and the first condition of (L).
Let us now study the third term in the r.h.s of (37). Observe that


i
1 h −1/2
WJ = √
Σ
⊗ X 0 (Σ−1/2 )0 ⊗ IdRn Vec(E)
nq
J
h

i

1
d
(49)
Z =: A1 Z,
=√
Σ−1 ⊗ X 0 (Σ1/2 )0 ⊗ IdRn
nq
J,•
where Z is a centered Gaussian random vector having a covariance matrix equal to identity
and
i

1 h −1
(50)
A1 = √
.
Σ ⊗ X 0 (Σ1/2 )0 ⊗ IdRn
nq
J,•
Using (45), we get for all j in J, for all γ in R and all |J| × |J| matrix D,




(51)
P (D WJ )j ≥ γ = P (DA1 Z)j ≥ γ ≤ P (kDk2 kA1 k2 kZk2 ≥ γ) ,
where A1 is defined in (50) and Z is a centered Gaussian random vector having a covariance
matrix equal to identity. Hence, for all j in J,
 

√ 



nq
λ
−1
−1
−1
e
WJ
P
(CJ,J ) − (CJ,J )
≥
|Bj | −
(CJ,J ) sign(BJ ) j
5
2nq
j


√ 
nq
λq|J|
−1
−1
−c2
e
≤ P (CJ,J ) − (CJ,J )
kA1 k2 kZk2 ≥
M3 q
−
.
5
2nqM2
2
Let us now bound kA1 k2 . Note that
h
i

Σ−1 ⊗ X 0 (Σ1/2 )0 ⊗ IdRn

J,• 2

=

h
i
Σ−1/2 ⊗ X 0

J,• 2


 1/2

1/2
= ρ Σ−1 ⊗ (X 0 X) J,J
≤ ρ Σ−1 ⊗ (X 0 X)
≤ λmax (Σ−1 )1/2 λmax (X 0 X)1/2 ,
where the first inequality comes from Theorem 4.3.15 of Horn and Johnson (1986). Hence,
by (A5) and (A7),
1
(52)
kA1 k2 ≤
λmax (Σ−1 )1/2 λmax (X 0 X)1/2 = OP (q −1/2 ).
nq
By (33), (42) and (52) it is thus enough to prove that
!
nq
X
2
−2c2
P
Zk ≥ nq n q
→ 0, as n → ∞.
k=1

The result follows from the Markov inequality and the first condition of (L).

20

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

Let us now study the fourth term in the r.h.s of (37). By (46), for all j in J,
 
P

≤P

eJ,J )−1 − CJ,J )−1
(C


√ 

nq
λ
≥
|Bj | −
(CJ,J )−1 sign(BJ ) j
5
2nq
j

√ 
nq
λq|J|
−c2
M3 q
−
,
kAk2 kZk2 ≥
5
2nqM2



fJ − WJ
W

eJ,J )−1 − CJ,J )−1
(C

2

where A is defined in (44).
By (33), (42) and (48), it is thus enough to prove that

P

nq
X

!
Zk2 ≥ (nq) n2 q 1−2c2

→ 0, as n → ∞.

k=1

The result follows from the Markov inequality and the fact that c2 < 1/2.
en ). By definition of B
en , we get that
Let us now study P(B
enc ) = P
P(B




λ 
−1 f
−1
e
f
e
e
e
c
c
c
CJ ,J (CJ,J ) WJ − WJ ≥ √
1 − |CJ ,J (CJ,J ) sign(BJ )|
.
2 nq

Observe that
eJ c ,J (C
eJ,J )−1 W
fJ − W
fJ c = CJ c ,J (CJ,J )−1 WJ − WJ c
C


fJ − WJ
+ CJ c ,J (CJ,J )−1 W


eJ,J )−1 − (CJ,J )−1 WJ
+ CJ c ,J (C



eJ,J )−1 − (CJ,J )−1 W
fJ − WJ
+ CJ c ,J (C


eJ c ,J − CJ c ,J (CJ,J )−1 WJ
+ C




fJ − WJ
eJ c ,J − CJ c ,J (CJ,J )−1 W
+ C



eJ c ,J − CJ c ,J (C
eJ,J )−1 − (CJ,J )−1 WJ
+ C




fJ − WJ
eJ c ,J − CJ c ,J (C
eJ,J )−1 − (CJ,J )−1 W
+ C
fJ c .
+ WJ c − W
Moreover,
eJ c ,J (C
eJ,J )−1 sign(BJ ) = CJ c ,J (CJ,J )−1 sign(BJ )
C


eJ,J )−1 − (CJ,J )−1 sign(BJ )
+ CJ c ,J (C


eJ c ,J − CJ c ,J (CJ,J )−1 sign(BJ )
+ C



eJ c ,J − CJ c ,J (C
eJ,J )−1 − (CJ,J )−1 sign(BJ ).
+ C

VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

By (IC) and the triangle inequality, we obtain that


λ
c
−1
e ) ≤ P CJ c ,J (CJ,J ) WJ − WJ c ≥ √ η
P(B
n
24 nq




λ
−1
f
WJ − WJ ≥ √ η
+ P CJ c ,J (CJ,J )
24 nq




λ
−1
−1
e
+P
CJ c ,J (CJ,J ) − (CJ,J )
WJ ≥ √ η
24 nq





λ
−1
−1
e
f
+P
CJ c ,J (CJ,J ) − (CJ,J )
WJ − WJ ≥ √ η
24 nq
 


λ
eJ c ,J − CJ c ,J (CJ,J )−1 WJ ≥ √
+P
C
η
24 nq

 



λ
−1
fJ − WJ ≥ √ η
eJ c ,J − CJ c ,J (CJ,J )
W
+P
C
24 nq
 



λ
−1
−1
e
e
+P
CJ c ,J − CJ c ,J (CJ,J ) − (CJ,J )
WJ ≥ √ η
24 nq

 



λ
−1
−1
f
e
e
WJ − WJ ≥ √ η
+P
CJ c ,J − CJ c ,J (CJ,J ) − (CJ,J )
24 nq


λ
fJ c ≥ √
η
+P
WJ c − W
24 nq
n


o
eJ,J )−1 − (CJ,J )−1 sign(BJ ) ≥ η
+P
CJ c ,J (C
12
n 

η o
−1
e
+P
CJ c ,J − CJ c ,J (CJ,J ) sign(BJ ) ≥
12
n 


η o
−1
−1
e
e
+P
sign(BJ ) ≥
.
CJ c ,J − CJ c ,J (CJ,J ) − (CJ,J )
12
(53)
The first term in the r.h.s of (53) tends to 0 by (35).
Let us now study the second term of (53). By (46), we get that for all j in J c ,




λ
−1
f
P
CJ c ,J (CJ,J )
WJ − WJ
≥ √ η
24 nq
j


λ
−1
≤ P kCJ c ,J k2 k(CJ,J ) k2 kAk2 kZk2 ≥ √ η .
24 nq
Observe that

k(X•,J c )0 k2 kX•,J k2
(X•,J c )0 X•,J (X•,J )0 X•,J c 1/2
1
kCJ c ,J k2 = ρ
=
(X•,J c )0 X•,J 2 ≤
√
√
nq
nq
nq
nq
nq

1/2 
1/2
0
0
(X•,J c ) X•,J c
(X•,J ) X•,J
ρ
≤ρ
nq
nq
λmax (Σ−1 )
(54)
= ρ(CJ c J c )1/2 ρ(CJ,J )1/2 ≤ ρ(C) =
λmax (X 0 X/n) = OP (q −1 ).
q


21

22

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

In (54) the last inequality and the fourth equality come from Theorem 4.3.15 of Horn and
Johnson (1986) and (40), respectively. The last equality comes from (A5) and (A7).
By (39), (48) and (54), it is thus enough to prove that
2 !

2 !

nq
nq
X
X
λ
λ
√
P
Zk2 ≥ (nq)1/2 q √
→ 0, as n → ∞,
=P
Zk2 ≥ (nq) √
nq
n
k=1

k=1

which holds true by the second condition of (L) and Markov inequality. Hence, the second
term of (53) tends to zero as n tends to infinity.
Let us now study the third term of (53). By (51), we get that for all j in J c ,





λ
−1
−1
e
P
CJ c ,J (CJ,J ) − (CJ,J )
WJ
≥ √ η
24 nq
j


λ
−1
−1
e
≤ P kCJ c ,J k2 k(CJ,J ) − (CJ,J ) k2 kA1 k2 kZk2 ≥ √ η .
24 nq
By (42), (52) and (54), it is thus enough to prove that
2 !

 !

nq
nq
X
X
λ
λ 2
2
2
1/2 √
q√
=P
Zk ≥ (nq) √
→ 0, as n → ∞,
P
Zk ≥ (nq)
nq
n
k=1

k=1

which holds true by the second condition of (L) and Markov inequality. Hence, the third term
of (53) tends to zero as n tends to infinity.
Let us now study the fourth term of (53). By (46), it amounts to prove that


λ
−1
−1
e
P kCJ c ,J k2 k(CJ,J ) − (CJ,J ) k2 kAk2 kZk2 ≥ √ η → 0, as n → ∞.
24 nq
By (54), (42) and (48) it is enough tho prove that

 !
nq
X
λ 2
2
P
Zk ≥ (nq) (nq) √
→ 0, as n → ∞,
n
k=1

which holds true by the second condition of (L). Hence, the fourth term of (53) tends to zero
as n tends to infinity.
Let us now study the fifth term of (53). By (51), proving that the fifth term of (53) tends
to 0 amounts to proving that


λ
−1
e
c
c
P CJ ,J − CJ ,J k(CJ,J ) k2 kA1 k2 kZk2 ≥ √ η → 0, as n → ∞.
24 nq
2
eJ c ,J k2 .
Let us now bound kCJ c ,J − C
CJ c ,J
≤
(55)

=



eJ c ,J
−C

e
C −C

2



=


J c ,J

1
b −1
Σ−1 − Σ
q



e
C −C

e
C −C




J c ,J 2
1/2

J c ,J ∞

X 0X
∞

n

≤

=ρ




e
C −C

e
C −C




J c ,J



e
C −C

e
C −C
1/2
∞

1/2


J c ,J

e
≤ C −C

= OP (q −1 (nq)−1/2 ),
∞

as n tends to infinity, where the last equality comes from (A5) and (A9).

∞

VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

23

By (39), (52) and (55), to prove that the fifth term of (53) tends to zero as n tends to
infinity, it is enough to prove that
 !

nq
X
λ 2
2
→ 0, as n → ∞,
P
Zk ≥ nq √
n
k=1

which holds using Markov’s inequality and the second condition of (L).
Using similar arguments as those used for proving that the second, third and fourth terms
of (53) tend to zero, we get that the sixth, seventh and eighth terms of (53) tend to zero, as
n tends to infinity, by replacing (54) by (55).
Let us now study the ninth term of (53). Replacing J by J c in (43), (44), (46, (47) and
(48) in order to prove that the ninth term of (53) tends to 0 it is enough to prove that

2 !
nq
X
λ
P
Zk2 ≥ nq √
→ 0, as n → ∞,
n
k=1

which holds using Markov’s inequality and the second condition of (L).
Let us now study the tenth term of (53). Using the same idea as the one used for proving
(38), we get that

n

o
eJ,J )−1 − (CJ,J )−1 sign(BJ ) ≥ η
P
CJ c ,J (C
12 
p
η
−1
−1
e
≤P
|J| kCJ c ,J k2 (CJ,J ) − (CJ,J )
≥
,
12
2
which tends to zero as n tends to infinity by (A3), (42), (54) and the fact that c1 < 1/2.
Let us now study the eleventh term of (53). Using the same idea as the one used for proving
(38), we get that

o
n 
eJ c ,J − CJ c ,J (CJ,J )−1 sign(BJ ) ≥ η
P
C
12 
p
−1
eJ c ,J − CJ c ,J k(CJ,J ) k2 ≥ η ,
≤P
|J| C
12
2
which tends to zero as n tends to infinity by (A3), (39) and (55) and the fact that c1 < 1/2.
Finally, the twelfth term of (53) can be bounded as follows:
n 


o
eJ c ,J − CJ c ,J (C
eJ,J )−1 − (CJ,J )−1 sign(BJ ) ≥ η
C
P
12

p
η
−1
−1
eJ c ,J − CJ c ,J
eJ,J ) − (CJ,J )
≤P
|J| C
(C
≥
,
12
2
2
which tends to zero as n tends to infinity by (A3), (42) and (55) and the fact that c1 < 1/2. 
Proof of Proposition 6. Observe that

1
−φ1
0
···
0
2
···
0
−φ1 1 + φ1 −φ1

..
..
..

−1
.
.
−φ1
.
(56)
Σ = 0
 .
.
.
 ..
..
. . 1 + φ2 −φ
1
1
0
0
···
−φ1
1





.



24

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

Let S = X 0 X = Σ−1 ⊗ X 0 X. Then,

nri +1



(1 + φ21 )nri +1
Si,j =
−φ1 nri +1



0

if j = i and ki ∈ {0, q − 1}
if j = i and ki ∈
/ {0, q − 1}
,
if j = i + p or if j = i − p
otherwise

where i − 1 = (p − 1)ki + ri corresponds to the Euclidean division of (i − 1) by (p − 1).
In order to prove (IC), it is enough to prove that
kSJ c ,J (SJ,J )−1 k∞ ≤ 1 − η,
where η ∈ (0, 1).
Since for all j, (j − p) ∈ J c or (j + p) ∈ J c ,
kSJ c ,J k∞ = ν|φ1 |.
Let A = SJ,J . Since A = (ai,j ) is a diagonally dominant matrix, then, by Theorem 1 of Varah
(1975),
1
X
kA−1 k∞ ≤
.
mink (ak,k −
ak,j )
1≤j≤|J|

j6=k

Using that for all j, (j − p) ∈

Jc

J c,

or (j + p) ∈
X
ak,j ≤ ν|φ1 |.
1≤j≤|J|

j6=k

If k ∈ J then k > p and k < pq − p. Thus,
ak,k ≥ ν(1 + φ21 ).
Hence,
kA−1 k∞ ≤

1
ν(1 + φ21 − |φ1 |)

and

|φ1 |
.
1 + φ21 − |φ1 |
Since |φ1 | < 1, the strong Irrepresentability Condition holds when
kSJ c ,J (SJ,J )−1 k∞ ≤ kSJ c ,J k∞ k(SJ,J )−1 k∞ ≤

|φ1 | ≤ (1 − η)(1 + |φ1 |2 − |φ1 |),
which is true for a small enough η.
Proof of Proposition 7. Since |φ1 | < 1,
kΣ−1 k∞ ≤ |φ1 | + |1 + φ21 | ≤ 3,
which gives (A7) by Theorem 5.6.9 of Horn and Johnson (1986).
Observe that
!


q−1
X
3 − |φ1 |
1
1
2
3
h
kΣk∞ ≤
1+2
|φ1 |
≤
=
≤
,
1+
2
2
2
1 − |φ1 |
1 − φ1
1 − φ1
1 − φ1
1 − φ21
h=1

which gives (A8) by Theorem 5.6.9 of Horn and Johnson (1986).



VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

25

c1
b −1 has the same expression as Σ−1 defined in (56) except that φ1 is replaced by φ
Since Σ
defined in (21), we get that
b −1
Σ−1 − Σ

∞


2
≤ 2 φ1 − φb1 + φ1 − φb1 ,

which implies Assumption (A9) of Theorem 5 by Lemma 8.
Let us now check Assumption (A10) of Theorem 5. Since, by Theorem 5.6.9 of Horn and
b ≤ kΣ − Σk
b ∞ , it is enough to prove that
Johnson (1986), ρ(Σ − Σ)
b
Σ−Σ

∞

= OP ((nq)−1/2 ), as n → ∞.

Observe that
b
Σ−Σ

∞

≤

q−1
X
1
φh1
φbh1
1
−
+
2
−
1 − φ21 1 − φb21
1 − φ21 1 − φb21
h=1

q−1
q−1
X
X
φh1 − φbh1
≤
+2
φbh1
+2
1 − φ21
(1 − φ21 )(1 − φb21 )
h=1
h=1

φ21 − φb21

1
1
−
2
1 − φ1 1 − φb21

!

!
q−1
q−1 

X
X
φh1 − φbh1
1
(φ1 − φb1 )(φ1 + φb1 )
1
≤
+2
+2
φbh1 − φh1
−
1 − φ21
1 − φ21 1 − φb21
(1 − φ21 )(1 − φb21 )
h=1
h=1
!
q−1
X
1
1
h
+2
φ1
−
1 − φ21 1 − φb21
h=1
! q−1


(φ1 − φb1 )(φ1 + φb1 )
(φ1 − φb1 )(φ1 + φb1 ) X bh
2
1
≤
1+
+2
φ1 − φh1 .
2| +
2
2
2
2
b
b
1
−
|φ
|
|1
−
φ
1
(1 − φ1 )(1 − φ1 )
(1 − φ1 )(1 − φ1 )
1
h=1
Moreover,
q−1
X

φbh1 − φh1 ≤ φb1 − φ1

h=1

≤ φb1 − φ1

q−1 X
h−1
X

k

|φ1 | |φb1 |h−k−1 ≤ φb1 − φ1

h=1 k=0

1
1 − |φb1 |

!

1
1 − |φ1 |

1 − |φb1 |q−1
1 − |φb1 |

!

1 − |φ1 |q−1
1 − |φ1 |




.

Thus, by Lemma 8,
b
Σ−Σ

∞

= OP ((nq)−1/2 ),

which implies Assumption (A10) of Theorem 5.
Proof of Lemma 8. In the following, for notational simplicity, q = qn . Observe that
Pn Pq b b
√1
i=1
√ b
`=2 Ei,` Ei,`−1
nq
nq φ1 =
.
P
P
q−1 b 2
n
1
i=1
`=1 Ei,`
nq



26

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

By (19),
q
n X
X

bi,` E
bi,`−1 =
E

i=1 `=2

q
X

b•,`−1 =
b•,` )0 E
(E

`=2

`=2

(57)

=

q
X

q
X
(ΠE•,` )0 (ΠE•,`−1 )

(φ1 ΠE•,`−1 + ΠZ•,` )0 (ΠE•,`−1 )

`=2
q−1
X

= φ1

0

(ΠE•,` ) (ΠE•,` ) +

`=1

q
X

(ΠZ•,` )0 (ΠE•,`−1 ),

`=2

where (57) comes from the definition of (Ei,t ).
Hence,
Pq
0
√1
√
`=2 (ΠZ•,` ) (ΠE•,`−1 )
nq
b
nq(φ1 − φ1 ) =
.
1 Pn Pq−1 b 2
i=1
`=1 Ei,`
nq
√
In order to prove that nq(φb1 − φ1 ) = OP (1), it is enough to prove that
n q−1
n q−1
1 XX 2
1 X X b2
Ei,` −
Ei,` = oP (1), as n → ∞,
nq
nq

(58)

i=1 `=1

i=1 `=1

by Lemma 10 and
q
1 X
(ΠZ•,` )0 (ΠE•,`−1 ) = OP (1), as n → ∞.
√
nq

(59)

`=2

Let us first prove (58). By (20),
Eb = [IdRq ⊗ Π] E := AE.
Note that
b = A(Σ ⊗ IdRn )A0 = Σ ⊗ Π.
Cov(E)
Hence, for all i
Var(Ebi ) ≤ λmax (Σ).
Since the covariance matrix of E is equal to Σ ⊗ IdRn , for all i
Var(Ei ) ≤ λmax (Σ).
By Markov inequality,
q
q
n q−1
n q−1
n
n
1 X X b2
1 XX 2
1 X X b2
1 XX 2
Ei,` −
Ei,` =
Ei,` −
Ei,` + oP (1)
nq
nq
nq
nq
i=1 `=1
i=1 `=1
i=1 `=1
i=1 `=1

1 
b 2 + oP (1).
=
kEk22 − kEk
2
nq
Observe that
b 2 = kEk2 − kAEk2 = E 0 E − E 0 A0 AE = E 0 (IdRnq − IdRq ⊗ Π) E
kEk2 − kEk
2

0

2

2

2

= E (IdRq ⊗ (IdRn − Π)) E =

pq
X
i=1

Eei2 ,

VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

27

where Ee = OE, where O is an orthogonal matrix. Using that
e i,i ≤ λmax (Σ),
E(Eei2 ) = Cov(E)
and Markov inequality, we get (58).
Let us now prove (59). By definition of (Ei,t ) and since |φ1 | < 1, E[(ΠZ•,` )0 (ΠE•,`−1 )] = 0.
Moreover,


2 
!2 
! n
q X
q
n
n
X
X
X
X
(ΠZ•,` )0 (ΠE•,`−1 )  = E 
Πi,k Zk,` 
E
Πi,j Ej,`−1  
`=2 i=1

`=2

=

X

X

k=1

j=1

Πi,k Πi0 ,k0 Πi,j Πi0 ,j 0 E Zk,` Zk0 ,`0 Ej,`−1 Ej 0 ,`0 −1



2≤`,`0 ≤q 1≤i,j,k,i0 ,j 0 ,k0 ≤n

=

X

X

Πi,k Πi0 ,k0 Πi,j Πi0 ,j 0

2≤`,`0 ≤q 1≤i,j,k,i0 ,j 0 ,k0 ≤n

X


φr1 φs1 E Zk,` Zk0 ,`0 Zj,`−1−r Zj 0 ,`0 −1−s ,

r,s≥0

since the (Ei,t ) are AR(1) processes with |φ1 | < 1. Note that E(Zk,` Zk0 ,`0 Zj,`−1−r Zj 0 ,`0 −1−s ) =
0 except when ` = `0 , k = k 0 , j = j 0 and r = s.
Thus,



!2 
q
q
X
X
X
X
0
4
2r 


(ΠZ•,` ) (ΠE•,`−1 )
=σ
E
φ1
Πi,k Πi0 ,k Πi,j Πi0 ,j
r≥0

`=2

=

`=2 1≤i,j,k,i0 ≤n

qσ 4
nqσ 4
Tr(Π) ≤
,
2
1 − φ1
1 − φ21

where Tr(Π) denotes the trace of Π, which concludes the proof of (59) by Markov inequality.

5. Technical lemmas
Lemma 9. Let A ∈ Mn (R) and Π an orthogonal projection matrix. For any j in {1, . . . , n}
(A0 ΠA)jj ≥ 0.
Proof of Lemma 9. Observe that
(A0 ΠA) = A0 Π0 ΠA = (ΠA)0 (ΠA),
since Π is an orthogonal projection matrix. Moreover,
(A0 ΠA)jj = e0j (ΠA)0 (ΠA)ej ≥ 0,
since (ΠA)0 (ΠA) is a positive semidefinite symmetric matrix, where ej is a vector containing
null entries except the jth entry which is equal to 1.

Lemma 10. Assume that (E1,t )t , (E2,t )t , ..., (En,t )t are independent AR(1) processes satisfying:
Ei,t − φ1 Ei,t−1 = Zi,t , ∀i ∈ {1, . . . , n},
where the Zi,t ’s are zero-mean i.i.d. Gaussian random variables with variance σ 2 and |φ1 | < 1.
Then,
n qn −1
σ2
1 XX 2 P
Ei,` −→
, as n → ∞.
nqn
1 − φ21
i=1
`=1

28

M. PERROT-DOCKÈS, C. LÉVY-LEDUC, L. SANSONNET, AND J. CHIQUET

2 ) =
Proof of Lemma 10. In the following, for notational simplicity, q = qn . Since E(Ei,`
σ 2 /(1 − φ21 ), it is enough to prove that
n q−1
 P
1 XX 2
2
Ei,` − E(Ei,`
) −→ 0, as n → ∞.
nq
i=1 `=1

Since

2


2
Ei,`
=

X

φj1 Zi,`−j  =

(60)

=

n q−1

1 XX 2
2
Ei,` − E(Ei,`
)
nq

1
(nq)2

i=1 `=1
n
X
X

X

0

φj1 φj1 Zi,`−j Zi,`−j 0 ,

j,j 0 ≥0

j≥0

Var

X

!

n
1 X
=
(nq)2

X

i=1

0

X

2
2
Cov(Ei,`
; Ei,`
0)

1≤`,`0 ≤q−1

0

φj1 φj1 φk1 φk1 Cov(Zi,`−j Zi,`−j 0 ; Zi,`0 −k Zi,`0 −k0 ).

i=1 1≤`,`0 ≤q−1 j,j 0 ≥0 k,k0 ≥0

By Cauchy-Schwarz P
inequality | Cov(Zi,`−j Zi,`−j 0 ; Zi,`0 −k Zi,`0 −k0 )| is bounded by a positive
constant. Moreover j≥0 |φ1 |j < ∞, hence (60) tends to zero as n tend to infinity, which
concludes the proof of the lemma.

References
Alquier, P. and P. Doukhan (2011). Sparsity considerations for dependent variables. Electron.
J. Statist. 5, 750–774.
Banerjee, O., L. E. Ghaoui, and A. D’aspremont (2008). Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data. Journal of Machine
Learning Research 9, 485–516.
Brockwell, P. J. and R. A. Davis (1990). Time Series: Theory and Methods. New York, NY,
USA: Springer-Verlag New York, Inc.
Fan, J. and R. Li (2001). Variable selection via nonconcave penalized likelihood and its oracle
properties. Journal of the American Statistical Association 96 (456), 1348–1360.
Friedman, J., T. Hastie, and R. Tibshirani (2008). Sparse inverse covariance estimation with
the graphical lasso. Biostatistics 9 (3), 432.
Horn, R. A. and C. R. Johnson (1986). Matrix Analysis. New York, NY, USA: Cambridge
University Press.
Lee, W. and Y. Liu (2012).
Simultaneous Multiple Response Regression and Inverse Covariance Matrix Estimation via Penalized Gaussian Maximum Likelihood.
J. Multivar. Anal. 111, 241–255.
Lütkepohl, H. (2005). New introduction to multiple time series analysis. Berlin: Springer.
Mehmood, T., K. H. Liland, L. Snipen, and S. Saebo (2012). A review of variable selection methods in partial least squares regression. Chemometrics and Intelligent Laboratory
Systems 118, 62 – 69.
Meng, C., B. Kuster, A. C. Culhane, and A. M. Gholami (2014). A multivariate approach to
the integration of multi-omics datasets. BMC Bioinformatics 15 (1), 162.
Perrot-Dockès, M., C. Lévy-Leduc, J. Chiquet, L. Sansonnet, M. Brégère, M. P. Étienne,
S. Robin, and G. Genta-Jouve (2017). A multivariate variable selection approach for analyzing LC-MS metabolomics data. arXiv:1704.00076.

VARIABLE SELECTION IN MULTIVARIATE LINEAR MODELS

29

Pourahmadi, M. (2013). High-Dimensional Covariance Estimation. Wiley Series in Probability and Statistics.
Rothman, A. J., P. J. Bickel, E. Levina, and J. Zhu (2008). Sparse permutation invariant
covariance estimation. Electron. J. Statist. 2, 494–515.
Rothman, A. J., E. Levina, and J. Zhu (2010). Sparse multivariate regression with covariance
estimation. Journal of Computational and Graphical Statistics 19 (4), 947–962.
Tibshirani, R. (1996). Regression shrinkage and selection via the Lasso. J. Royal. Statist. Soc
B. 58 (1), 267–288.
Varah, J. (1975). A lower bound for the smallest singular value of a matrix. Linear Algebra
and its Applications 11 (1), 3 – 5.
Yuan, M. and Y. Lin (2007). Model selection and estimation in the gaussian graphical model.
Biometrika 94 (1), 19–35.
Zhao, P. and B. Yu (2006). On model selection consistency of lasso. Journal of Machine
Learning Research 7, 2541–2563.
UMR MIA-Paris, AgroParisTech, INRA, Université Paris-Saclay, 75005, Paris, France
E-mail address: marie.perrot-dockes@agroparistech.fr
UMR MIA-Paris, AgroParisTech, INRA, Université Paris-Saclay, 75005, Paris, France
E-mail address: celine.levy-leduc@agroparistech.fr
UMR MIA-Paris, AgroParisTech, INRA, Université Paris-Saclay, 75005, Paris, France
E-mail address: laure.sansonnet@agroparistech.fr
UMR MIA-Paris, AgroParisTech, INRA, Université Paris-Saclay, 75005, Paris, France
E-mail address: julien.chiquet@agroparistech.fr

