arXiv:1712.09483v2 [] 20 Feb 2018

Minimax Estimation of Large Precision Matrices with
Bandable Cholesky Factor
Yu Liu and Zhao Ren
University of Pittsburgh1

Abstract
Last decade witnesses significant methodological and theoretical advances in estimating large precision matrices. In particular, there are scientific applications such
as longitudinal data, meteorology and spectroscopy in which the ordering of the variables can be interpreted through a bandable structure on the Cholesky factor of the
precision matrix. However, the minimax theory has still been largely unknown, as
opposed to the well established minimax results over the corresponding bandable covariance matrices. In this paper, we focus on two commonly used types of parameter
spaces, and develop the optimal rates of convergence under both the operator norm
and the Frobenius norm. A striking phenomenon is found: two types of parameter
spaces are fundamentally different under the operator norm but enjoy the same rate
optimality under the Frobenius norm, which is in sharp contrast to the equivalence
of corresponding two types of bandable covariance matrices under both norms. This
fundamental difference is established by carefully constructing the corresponding minimax lower bounds. Two new estimation procedures are developed: for the operator
norm, our optimal procedure is based on a novel local cropping estimator targeting
on all principle submatrices of the precision matrix while for the Frobenius norm, our
optimal procedure relies on a delicate regression-based block-thresholding rule. We
further establish rate optimality in the nonparanormal model. Numerical studies are
carried out to confirm our theoretical findings.

Keywords: optimal rate of convergence, precision matrix, local cropping, Cholesky factor, minimax lower bound, block-thresholding, operator norm, Frobenius norm
AMS 2000 Subject Classification: Primary 62H12; secondary 62F12, 62C20, 62G09.

1

Department of Statistics, University of Pittsburgh, Pittsburgh, PA 15260. e-mail: yul125@pitt.edu,
zren@pitt.edu

1

1

Introduction

Covariance matrix plays a fundamental role in many important multivariate statistical
problems. They include the principal component analysis, linear and quadratic discriminant analysis, clustering analysis, regression analysis and conditional dependence relationship studies in graphical models. During the last two decades, with the advances of
technology, it is very common that the datasets are high-dimensional (the dimension p can
be much larger than the sample size n) in many applications such as genomics, fMRI data,
astrophysics, spectroscopic imaging, risk management, portfolio allocation and numerical
weather forecasting [26, 20, 25, 33, 45, 39]. It has been well-known that the sample covariance matrix performs poorly and can yield to invalid conclusions in the high-dimensional
settings. For example, see [49, 50, 29, 21, 40, 30] for details on the limiting behaviors of
the spectra of sample covariance matrices when both n and p increase.
To avoid the curse of dimensionality, certain structural assumptions are almost necessary in order to estimate the covariance matrix or its inverse, the precision matrix,
consistently. In this paper, we consider large precision matrix estimation with bandable
Cholesky factor. Its connection with other structures are discussed at the end of introduction. Both the operator norm loss (kSkop = supkxk2 =1 kSxk2 ) and the Frobenius norm
P
1
loss (kSkF = ( i,j s2ij ) 2 ) are investigated.
We begin with introducing the bandable Cholesky factor of the precision matrix. Assume that X = (X1 , . . . Xp )T is a centered p-variate random vector with covariance matrix
Σ. Let ai = (ai1 , . . . , ai(i−1) )T be the coefficients of the population regression of Xi on
Pi−1
ait Xt =
its previous variables X1,i−1 = (X1 , X2 . . . Xi−1 )T . In other words, X̂i = t=1
T
X1,i−1 ai is the linear projection of Xi on X1,i−1 in population (Define X̂1 = 0). Set A
as the lower triangular matrix with zeros on the diagonal and zero-padded coefficients
(aTi , 0) arranged in the rows. Denote the residual  = X − X̂ = (I − A)X and D = Var().
The regression theory implies the residuals are uncorrelated, and thus the matrix D is
diagonal. The modified Cholesky decomposition of Ω is
Ω = Σ−1 = (I − A)T D−1 (I − A),

(1)

where I − A is the Cholesky factor of Ω. There is a natural order on the variables
based on the above Cholesky decomposition. Indeed, the well-known AR(k) model can
be characterized by the k-banded Cholesky factor A ≡ [aij ]p×p of the precision matrix
in which aij = 0 if i − j > k. Inspired by the auto-regression model, we consider the
bandable structures imposed on the Cholesky factor. More specifically, for M > 0, η > 1
we define the parameter space Pα (η, M ) of precision matrices by
n
Pα (η, M ) = Ω : η −1 ≤ λmin (Ω) ≤ λmax (Ω) < η,
o
X
(2)
max
|aij | < M k −α , k ∈ [p] .
i

j<i−k

2

Here, λmax (Ω), λmin (Ω) are the maximum and minimum eigenvalues of Ω and the index
set [p] = {1, 2, . . . , p}. We follow the convention that the sum over an empty set of
indices is equal to zero when i − k ≤ 1. This parameter space was first proposed in
[7]. The parameter α specifies how fast the sequence aij decays to zero as j goes away
from i. The covariance matrix estimation problem has been extensively studied when
a similar bandable structure is imposed on the covariance matrix (e.g., [7, 13]). Unlike
the order in these bandable covariance matrices, in which large distance |i − j| implies
nearly independence, the order in bandable Cholesky factor encodes nearly conditional
independence in the sense that the coefficients aij is close to zero when i − j > 0 is large.
Although several approaches have been developed to estimate the precision matrix
with bandable Cholesky factor, the optimality question remains mostly open, partially
due to the following two reasons. (i) Intuitively, one would expect the minimax rate of
convergence over Pα (η, M ) under the operator norm to be the same as that over the class
of bandable covariance matrices with the same decay parameter α. Under sub-Gaussian
−2α+1
assumptions, [13] established the optimal rate of convergence EkΩ̃ − Ωk2op  n 2α + logn p
uniformly for all bandable covariance matrices Σ = Ω−1 = [σij ]p×p with bounded spectra
P
such that maxi |j−i|>k |σij | < M k −α , k ∈ [p]. To establish such a rate of convergence
for Pα (η, M ), [34] provided a lower bound with the matching rate. However, we show a
surprising result in this paper that estimation over Pα (η, M ) is a much harder task than
that over bandable covariance matrices. Therefore, the lower bound in [34] is sub-optimal,
−2α+1
and all attempts on showing the same rate of convergence n 2α + logn p intrinsically
cannot succeed. (ii) From the methodological aspect, due to the regression interpretation
of the Cholesky decomposition (1), almost all existing methods rely on an intermediate
estimator of A obtained by running regularized regression of each variable against its
P
previous variables Xi ∼ i−1
j=1 aij Xj . For instance, [7] estimated each row of A by fitting
Pi−1
aij Xj with some bandwidth k. [51] used
the banded regression model Xi ∼ j=max{1,i−k}
an AIC or BIC penalty to pick the best bandwidth k. In addition, [28] proposed adding a
Lasso or Ridge penalty while [35] proposed using a nested Lasso penalty to the regression.
See, for instance, [3, 34] for Bayesian approaches following the similar idea. The typical
analysis for those estimation procedures in a row-wise fashion is to bound the operator
norm by its matrix `1 /`∞ norm. Although this analysis may provide optimal rates of
convergence under the operator norm loss for some sparsity structure (see, i.e., [14, 10]
for sparse covariance and precision matrices estimation), it might be sub-optimal for the
bandable structure as seen in bandable covariance matrix estimation [13, 7]. Therefore, in
order to obtain rate-optimality over Pα (η, M ), a novel analysis or even a new estimation
approach is expected.
Main Results. With regard to the above two issues, we provide satisfactory solutions
in this paper. We at the first time show that the rate of convergence under the operator
norm over Pα (η, M ) is intrinsically slower than that over the counterpart class of bandable
covariance matrices. This is achieved via a novel minimax lower bound construction.
Moreover, in order to obtain a rate-optimal estimator, we propose a novel local cropping
3

estimator which does not rely on any estimator of A, and thus requires a new analysis.
Our local cropping approach targets on accurate estimation of principal submatrices of the
precision matrix under the operator norm, which results in a tradeoff between one variance
term and two bias terms. The name comes after the idea of estimating each principal
submatrix of the precision matrix, which is to crop the center k by k submatrix of the
inverse of 3k by 3k sample covariance matrix using their neighbors in two directions of the
same size. (During the finalizing process of this paper, we realized that a similar estimator
is independently proposed to estimate precision matrices with a different structure [27].)
Since our procedure does not directly explore the structure on each row of A, the analysis
of bias terms is much more involved, requiring a block-wise partition strategy. More
details are discussed in Sections 2.1 and 3.1. In the end, besides Pα (η, M ), a similar type
of classes of parameter spaces with bandable Cholesky factor is considered as well,
n
Qα (η, M ) = Ω : η −1 ≤ λmin (Ω) ≤ λmax (Ω) < η,
o
(3)
|aij | < M (i − j)−α−1 , j ∈ [i − 1] .
We further establish another surprising result: the optimal rates of convergence of two
spaces, namely Pα (η, M ) and Qα (η, M ), are different under the operator norm. This
remarkable distinction is different from the comparison of two similar types of parameter
spaces for bandable covariance matrices in [13] and bandable Toeplitz covariance matrices
in [11]. The contrast of minimax results on Pα (η, M ) and Qα (η, M ) is summarized in
Theorem 1 below. We mainly focus on the high-dimensional setting, assuming that log p =
O(n) and n = O(p). Otherwise, one can always obtain trivial constant minimax rate (i.e.,
inconsistency) or the minimax rate as the smaller of p/n and the one in Theorem 1.
Theorem 1. Under normality assumption, the minimax risk of estimating the precision
matrix Ω over the parameter space Pα (η, M ) with α > 12 given in (2) under the operator
norm satisfies
2α−1
log p
inf sup EkΩ̃ − Ωk2op  n− 2α +
.
(4)
n
Ω̃ Pα (η,M )
The minimax risk of estimating the precision matrix Ω over the parameter space Qα (η, M )
with α > 0 given in (3) under the operator norm satisfies
inf

sup

Ω̃ Qα (η,M )

2α

EkΩ̃ − Ωk2op  n− 2α+1 +

log p
.
n

(5)

Moreover, we also consider the minimax rates of convergence of precision matrix estimation under the Frobenius norm loss over Pα (η, M ) and Qα (η, M ). This time, we prove
that two types of spaces enjoy the same optimal rate of convergence. Together with the
different rates of convergence under the operator norm loss, we demonstrate the intrinsic
difference between operator norm and Frobenius norm. The Frobenius norm of a p by p
matrix is defined as the `2 vector norm of all entries. Driven by this fact, our estimation
4

approach is naturally obtained by optimally estimating A and D in (1) separately. Due
to the decay structure in Pα (η, M ), which is defined in terms of nested `1 norm of each
row of A, our estimator is based on regression with a delicate block-thresholding rule.
The minimax procedure is motivated by wavelet nonparametric function estimation, although the space Pα (η, M ) cannot be exactly described by any Besov ball ([9, 19]). We
summarize the optimality result under the Frobenius norm in Theorem 2 below.
Theorem 2. Under normality assumption, the minimax risk of estimating the precision
matrix Ω over Pα (η, M ) and Qα (η, M ) given in (2) and (3) satisfies
inf
Ω̃

2α+1
1
1
EkΩ̃ − Ωk2F  inf sup
EkΩ̃ − Ωk2F  n− 2α+2 .
Ω̃ Qα (η,M ) p
Pα (η,M ) p

sup

(6)

Related Literature. During the last decade, various structural assumptions are
imposed in the literature of high-dimensional statistics in order to estimate the covariance/precision matrix consistently under various loss functions. While mostly driven by
the specific scientific applications, popular structures include ordered sparsity (bandable
covariance matrices, precision matrices with bandable Cholesky factor), unordered sparsity (sparse covariance matrices, sparse precision matrices) and other more complicated
ones such as certain combination of sparsity and low-rankness (spike covariance matrices,
covariance with tensor product, latent graphical models). Many estimation procedures
have been proposed accordingly to estimate high-dimensional covariance/precision matrices via taking advantages of these specific structures. For example, banding ([52, 7, 53, 8])
and tapering methods ([24, 13]) were developed to estimate bandable covariance matrices or precision matrices with bandable Cholesky factor; thresholding procedures were
used in [6, 22, 15] to estimate sparse covariance matrices; penalized likelihood estimation
([28, 56, 18, 2, 43, 32, 41]) and penalized regression methods ([37, 55, 16, 47, 42]) are
designed for sparse precision matrix estimation.
The fundamental difficulty of various covariance/precision matrices estimation problems have been carefully investigated in terms of the minimax risks under the operator
norm loss among other losses, especially for those ordered and unordered sparsity structures. Specifically, for unordered structures, [14] considered the problems of optimal
estimation of sparse covariance while [10] (see [42] as well) established the optimality
results for estimating sparse precision matrices. For ordered structures, [13] established
the optimal rates of convergence over two types of bandable covariance matrices. In addition, with an extra Toeplitz structure, [11] studied optimal estimation of two types of
bandable Toeplitz covariance matrices. However, it was still largely unknown about the
optimality results on estimating precision matrices with bandable Cholesky factor. See
an exposure paper with discussion [12] and references therein on minimax results of covariance/precision matrix estimation under some other losses. In this paper, we provide
a solution to this open problem by establishing the optimal rates of convergence over two
types of precision matrices with bandable Cholesky factor. Thus, this paper completes the
minimaxity results of all four sparsity structures commonly considered in literature.
5

Organization of the Paper. The rest of the paper is organized as follows. First,
we propose our estimation procedures for precision matrix estimation in Section 2. In
particular, local cropping estimators and regression-based block-thresholding estimators
are designed for estimating precision matrices under the operator norm and the Frobenius
norm respectively. Section 3 establishes the optimal rates of convergence under the operator norm for two commonly used types of parameter spaces Pα (η, M ) and Qα (η, M ). A
striking difference between two spaces are revealed when considering operator norm loss.
Section 4 considers rate-optimal estimation under the Frobenius norm The results reveal
that the fundamental difficulty of estimation for two parameter spaces are the same when
considering Frobenius norm loss. In Section 5, we extend the results to nonparanormal
models for inverse correlation matrix estimation by applying local cropping procedure to
rank-based estimators. Section 6 presents the numerical performance of our local cropping
procedure to illustrate the difference between two parameter spaces by simulation studies.
We also demonstrate the sub-optimality of banding estimators. All technical lemmas used
in proofs of main results are relegated to the supplement.
Notation. we introduce some basic notations that will be used in the rest of the
paper. 1 (·) indicates the indicator function while 1 indicates the all-ones vector. sgn(·)
indicates the sign function. bsc represents the largest integer which is no more than s. dse
represents the smallest integer which is no less than s. Define an  bn if there is a constant
C > 0 independent of n such that C −1 ≤ an /bn ≤ C. For any vector x, kxkp indicates its
`p norm. For any p by q matrix S = [sij ]p×q ∈ Rp×q , we use S T to denote its transpose.
The `p matrix norm is define as kSkp = supkxkp =1 kSxkp . The `2 matrix norm is also
called the the operator norm or the spectral norm, and denoted as kSkop . The Frobenius
P
1
norm is defined as kSkF = ( i,j s2ij ) 2 . λmax (S) and λmin (S) are the largest and smallest
singular values of S when S is not symmetric. When S is a real symmetric matrix, λmax (S)
and λmin (S) denote its largest and smallest eigenvalues. rowi (S) and coli (S) indicate the
i-th row and column of matrix S. a : b denotes the index set {a, a + 1, . . . , b}. [p] is short
for the set 1 : p. For the random vector X ∈ Rp×1 and the data matrix Z ∈ Rn×p , Xa:b
and Za:b indicates the (a : b)-th columns of XT and Z. For any square matrix S, diag(S)
denotes the diagonal matrix with diagonal entries being those on the main diagonal of S
while for any vector v, diag(v) denotes the diagonal matrix with diagonal entries being v.
In the estimation procedure under the operator norm, we use the matrix notation in the
(k)
form of Sm to facilitate the proof, where S is always a square matrix, m indicates the
(k)
location information, and (k) indicates that the size of Sm is k. Throughout the paper
we denote by C a generic positive constant which may vary from place to place but only
depends on α, η, M and possibly some sub-Gaussian distribution constant ρ in (17).

2

Methodologies

In this section, we introduce our methodologies for estimating precision matrices over
Pα (η, M ) and Qα (η, M ) under both the operator norm and the Frobenius norm. Assume
6

that X = (X1 , . . . Xp )T , a p-variate random vector with mean zero and precision matrix
Ωp . Our estimation procedures are based on its n i.i.d. copies Z ∈ Rn×p . We write Z =
(Z1 , . . . , Zp ), where each Zi consists of n i.i.d. copies of Xi . Our estimation procedures
are different under the operator norm and the Frobenius norm.

2.1

Estimation Procedure under the Operator Norm

We focus on the estimation problem under the operator norm first. As we discussed
in the introduction, almost all existing methodologies ([51, 28, 7]) directly appeal the
Cholesky decomposition of the precision matrix. They first estimate the Cholesky factor
A and D by auto-regression and then estimate the precision matrix according to Ω =
(I − A)T D−1 (I − A). The corresponding analysis in the row-wise fashion may not suitable
for the operator norm loss. In this paper, we propose a novel local cropping estimator,
which focuses on the estimation of Ω directly.
To facilitate the illustration of the estimation procedure, we define two matrix operators. The cropping operator is designed to crop the center block out of the matrix.
For a p by p matrix E ≡ [eij ]p×p , we define the k × k matrix Ckm (E) ≡ [cij ]k×k , where
1 ≤ m ≤ p − k + 1, with
cij = ei+m−1,j+m−1 , when 1 ≤ i, j ≤ k.

(7)

The parameter m indicates the location and k indicates the dimension. It is clear that
Ckm (E) is a principal submatrix of E. The expanding operator is designed to put a
small matrix onto a large zero matrix. For a k by k matrix, C ≡ [cij ]k×k , define the p × p
matrix Epm (C) ≡ [eij ]p×p , where 1 ≤ m ≤ p − k + 1, with
eij = ci−m+1,j−m+1 , when m ≤ i, j ≤ m + k − 1, otherwise eij = 0.

(8)

The parameter m indicates the location and p indicates the dimension. Note that for a k
by k matrix C, we have Ckm (Epm (C)) = C. An illustration of two operators is provided in
Figure 1.
In addition, for technical reasons (of obtaining rates of convergence in expectation
rather than in probability), we introduce a projection operator. For a real square
matrix S, let the singular value decomposition of S be S = U ΛV T with U U T = I,
V V T = I and Λ = diag(λi ). Let Λ∗ = diag(λ∗i ), where λ∗i = min{max{λi , η −1 }, η}, then
define
Pη (S) = U Λ∗ V T .
(9)
For a symmetric matrix S, we modify Pη (·) a little bit and define Pη (S) = U Λ∗ U T , where
S = U ΛU T is its eigen-decomposition. Since all eigenvalues of Pη (·) are in the interval
[η −1 , η], Pη (S) is always invertible and positive definite.
We are ready to construct the local cropping estimator Ω̃op
k with bandwidth k < p. At
a high level, we first propose an estimator of each principal submatrix of size k and 2k in
7

m

k

k
Ckm
Epm

p
Figure 1: An illustration of the cropping operator and the expanding operator.
Ω using cropping and extending operators. Then we arrange over those local estimators
to estimate Ω. Since the core idea of estimating those local estimators in our procedure is
to crop the inverse of sample covariance matrix with a relatively larger size, we call Ω̃op
k
in (12) the local cropping estimator.
(k)
Specifically, we first define an estimator Ω̃m of the principal submatrix Ckm (Ω) at each
location m. To this end, we select the sample covariance matrix with a relative larger
size, in this case, 3k. Let the modified local sample covariance matrix be
1 T 
(3k)
Σ̃m−k = Pη C3k
m−k ( Z Z) .
n

(10)

(3k)

Note that the operator Pη (·) guarantees Σ̃m−k to be invertible. Then we use the center
part of its inverse to estimate Ckm (Ω), i.e.,
(3k) −1 
k
.
Ω̃(k)
m = Ck+1 (Σ̃m−k )

(11)

(2k)

Similarly, we can define local estimators of Ω̃m via replacing k by 2k. Arranging over
these estimators in the form of weighted sum, we obtain the estimator of Ω, that is,
Ω̃op
k

= Pη

1
k

p
X

Epm (Ω̃(2k)
m )

m=2−2k

−

p
X

Epm (Ω̃(k)
m )



.

(12)

m=2−k

The operator Epm (·) makes these local estimators in the correct places. The final step (12)
is motivated by the analysis of optimal bandable covariance matrix estimation procedure
proposed in [13]. Indeed, the optimal tapering estimator in [13] can be rewritten as a
sum of many principal submatrices of the sample covariance matrix in a similar way as
(12). In contrast, our estimator is not in a form of tapering the sample covariance matrix.
However, in the analysis of our local cropping estimator in Section 3, the direct target of
Ω̃op
k is a certain tapered population precision matrix with bandwidth k. There are natural
bias and variance terms involved in the distance of Ω̃op
k and its direct target. Together
8

with the bias of the tapered population precision matrix, our analysis involves two bias
terms and one variance term, which critically determine the optimal choice of bandwidth.
In Section 3, we show that the local cropping estimator with an optimal choice of
bandwidth would achieve the minimax risk under the operator norm over parameter spaces
Pα (η, M ) in (2) and Qα (η, M ) in (3). However, the optimal choices of bandwidth are
fundamentally distinct between Pα (η, M ) and Qα (η, M ). Specifically, we show that the
1
1
optimal bandwidth over Pα (η, M ) is k  n 2α while that one over Qα (η, M ) is k  n 2α+1 .
Remark 1. Of note, the estimator Ω̃op
k depends on Z2−4k , . . . , Zp+4k−1 . The index of variable is clear most of the time, while we need to be careful when it is close to the boundary.
When the index is beyond the index set [p], we shrink the size of the corresponding block
by discarding the data with meaningless indexes.

2.2

Estimation Procedure under the Frobenius Norm

Under the Frobenius norm, our estimation procedure is based on the Cholesky decomposition of the precision matrix (1). More specifically, we estimate the matrix A and
D respectively by auto-regression, and then combine them together to construct the estimator of Ω. The following estimation procedure applies to both the parameter space
Pα (η, M ) and Qα (η, M ) as we will show that they enjoy the same optimal rate of convergence in Section 4.
Our estimator of the i-th row of A is based on the regression of Xi against its previous
variables. Unlike those existing methods ([51, 28, 7]) which rely on certain banding
or penalized approaches for such a regression problem, we apply a block-thresholding
procedure due to the decay structure in Pα (η, M ) which is defined in terms of nested
`1 norm. To this end, we first regress Xi against Xi−k1 :i−1 = (Xi−k1 , . . . , Xi−1 )T with
bandwidth k1 = d nc e with some sufficiently large c > 0. Recall that the n × 1 matrix Zi
consists of n observations of Xi , and the n × k1 matrix Zi−k1 :i−1 represents n observations
of Xi−k1 :i−1 . The empirical regression coefficients are
(âi(i−k1 ) , . . . âi(i−1) )T = (ZTi−k1 :i−1 Zi−k1 :i−1 )−1 ZTi−k1 :i−1 Zi .

(13)

We then further threshold the coefficients by taking advantages of the bandable structure
of the Cholesky factor A. Specifically, we define â∗i ∈ Ri−1 with coordinate â∗ij as follows,

â∗ij

=




âij ,

if i − k0 < j ≤ i − 1,

âij 1 (|âij | > λj ) , if i − k1 < j ≤ i − k0 ,


0,
if 1 ≤ j ≤ i − k1 ,

1

1

(14)

where k0 = dn 2α+2 e, λj = (dlog2i−j − logk20 eR) 2 , and R = 8ηk(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop .
Note that we keep the last k0 coefficients and apply a block-thresholding rule in which the
size of each block doubles backwards sequentially for the remaining coefficients in (14).
9

Our procedure is inspired by the optimal estimation procedure over Besov balls for many
nonparametric function estimation problems, or equivalently, the corresponding Gaussian
sequence models (See [9] the reference therein). We emphasize that any linear estimator
of the coefficients (âi(i−k1 ) , . . . âi(i−1) )T cannot yield to the optimal rates of convergence
in our setting under the Frobenius norm.
Our estimator of I − A can be constructed by arranging zero-padded â∗T
i , i ∈ [p]
∗
accordingly with an identity matrix. Specifically, set the ij-th entry of Â as â∗ij when
i ∈ [p], j ∈ [i−1], otherwise as zero. We also need to bound the singular values of (I − Â∗ ).
To this end, we define
I^
− A = Pη (I − Â∗ ),
as our estimator of (I − A), where Pη (·) is defined in (9).
The estimation of D is based on the sample variances of those empirical residuals
in the regression of Xi against Xi−k1 :i−1 = (Xi−k1 , . . . , Xi−1 )T . For each i, the sample
variance of the empirical residual is
1
ZT (I − Mi )T (I − Mi )Zi ,
(15)
dˆi =
n − k1 i
ˆ where dˆ =
where Mi = Zi−k1 :i−1 (ZTi−k1 :i−1 Zi−k1 :i−1 )−1 ZTi−k1 :i−1 . Let D̂ = diag(d),
T
(dˆ1 , . . . , dˆp ) . We define D̃ = Pη (D̂) as our estimator of D.
Finally, define our estimator of Ω as
^ T −1 ^
Ω̃F
k = (I − A) D̃ (I − A).

(16)

Remark 2. For the parameter space Qα (η, M ), a much simpler banding estimation
scheme on the empirical regression coefficients is able to achieve the minimax risk. Set
1
k = dn 2α+2 e. We use the empirical residuals and coefficients obtained by regressing each
Xi against Xi−k:i−1 to directly construct the estimators of A and D. It can be proved that
this estimator achieves the minimax risk over the parameter space Qα (η, M ).

3

Rate Optimality under the Operator Norm

In this section, we establish the optimal rates of convergence for estimating the precision
matrix over the parameter spaces Pα (η, M ) and Qα (η, M ) given in (2) and (3) under the
operator norm. We first derive the risk upper bound of the local cropping estimator in
Section 3.1 over parameter space Pα (η, M ). We provide a matching risk lower bound by
applying the Assouad’s lemma and the Le Cam’s method in Section 3.2 over Pα (η, M ).
The establishment of the rate optimality over the parameter space Qα (η, M ) is similar to
the one over Pα (η, M ), which is provided in Section 3.3.
Throughout this section, we assume that X = (X1 , . . . Xp )T follows certain subGaussian distribution with constant ρ > 0, that is,
P{|v T (X − EX)| > t} ≤ 2 exp(−t2 /(2ρ)),
for all t > 0 and all unit vectors kvk2 = 1.
10

(17)

3.1

Minimax Upper Bound under the Operator Norm over Pα (η, M )

In this section, we develop the following upper bound of our estimation procedure proposed
in Section 2.1.
1

Theorem 3. When dn 2α e ≤ p, the local cropping estimator defined in (12) of the precision
matrix Ω over Pα (η, M ) with α > 12 given in (2) satisfies
2
−2α+1
+C
sup EkΩ̃op
k − Ωkop ≤ Ck

Pα (η,M )

log p + k
.
n

1

When k = dn 2α e, we have
2
−
sup EkΩ̃op
k − Ωkop ≤ Cn

2α−1
2α

+C

Pα (η,M )

log p
.
n

1

The optimal choice of k  n 2α is due to the bias-variance tradeoff. Combining Theorem 3 with the minimax lower bound derived in Section 3.2, we immediately obtain that
the local cropping estimator is rate optimal.
Proof. As we discussed in Section 2.1, the direct target of our local cropping estimator
is certain tapered population precision matrix with bandwidth k, which can be written
as a weighted sum of many principal submatrices of the population precision matrix. We
construct this corresponding tapered population precision matrix Ω∗k as follows. Denote
∗]
the precision matrix Ω ≡ [ωij ]p×p . We define Ω∗k ≡ [ωij
p×p such that for i, j ∈ [p],
1
1
∗
ωij
= mij ωij , where mij = max{0, 2 − |i − j|} − max{0, 1 − |i − j|}.
(18)
k
k
The following lemma elucidates the decomposition of this tapered precision matrix Ω∗k .
Lemma 1. The Ω∗k defined in (18) can be written as
Ω∗k =

2k+1
1 X
k
m=2

bp/2kc

Epm+2kj C2k
m+2kj (Ω)

X



j=−1

−

k+1
X

bp/kc

m=2

j=−1

X

Epm+kj Ckm+2kj (Ω)



.

The proof of Lemma 1 can be found in [13] (refer to the proof of Lemma 1 with
covariance matrix therein replaced by the precision matrix), and thus omitted. Define
Ω̃∗k =

2k+1
1 X
k
m=2

bp/2kc

X

k+1
 X
(2k)
Epm+2kj (Ω̃m+2kj ) −

bp/kc

m=2

j=−1

j=−1

X

(k)

Epm+kj (Ω̃m+kj )



.

∗
−1
It is easy to check Ω̃op
k = Pη (Ω̃k ). Since the eigenvalues of Ω are in the interval [η , η], the
operator Pη (·) would not increase the risk much. Indeed, according to (53) in Lemma 21,
we have
2
∗
2
EkΩ̃op
k − Ωkop ≤ 4EkΩ̃k − Ωkop

≤ 8EkΩ̃∗k − Ω∗k k2op + 8kΩ∗k − Ωk2op .
11

(19)

The following lemma bounds the bias between our direct target Ω∗k and the population
precision matrix.
Lemma 2. For Ω in the parameter space Pα (η, M ) defined in (2) with α > 12 , Ω∗k defined
in (18), we have
kΩ∗k − Ωk2op ≤ Ck −2α+1 .
Remark 3. Unlike existing methods, our procedure does not directly utilize the decay
structure of Cholesky factor. Consequently, the proof of Lemma 2 is involved and requires
a block-wise partition strategy.
Then we turn to the analysis of EkΩ̃∗k − Ω∗k k2op .
EkΩ̃∗k − Ω∗k k2op
bp/2kc
2k+1 bp/2kc
X p
 2
1 X X p
(2k)
k
Em+2kj (Ω̃m+2kj ) −
≤2E
Em+2kj C2k
m+2kj (Ω) kop
k
m=2

j=−1

(20)

j=−1

bp/kc
k+1 bp/kc
X p
 2
1 X X p
(k)
Em+kj Ckm+kj (Ω) kop .
k
Em+kj (Ω̃m+kj ) −
+ 2E
k
m=2

j=−1

j=−1

These two terms can be bounded in the same way, we only focus on the second term.
bp/kc
k+1 bp/kc
X p
 2
1 X X p
(k)
E
k
Em+kj (Ω̃m+kj ) −
Em+kj Ckm+kj (Ω) kop
k
m=2

j=−1

j=−1

bp/kc

≤E max k
m

X

 
(k)
(Epm+kj (Ω̃m+kj ) − Epm+kj Ckm+kj (Ω) k2op

j=−1
(k)

≤E max kΩ̃m+kj − Ckm+kj (Ω)k2op



m,j

 
k
3k
−1 −1 2
≤2E max kΩ̃(k)
kop
m − Ck+1 (Cm−k (Ω ))
m∈[p]

(21)



−1 −1
+ 2 max kCkk+1 (C3k
− Ckm (Ω)k2op ,
m−k (Ω ))
m∈[p]

where we further have variance term and bias term of local estimators. For the variance
term in (21), we further have

k
3k
−1 −1
kΩ̃(k)
kop
m − Ck+1 (Cm−k (Ω ))
(3k)

−1 −1
≤k(Σ̃m−k )−1 − (C3k
m−k (Ω )) kop
(3k)

−1
≤η 2 kΣ̃m−k − C3k
m−k (Ω )kop
1
T
3k
−1
≤2η 2 kC3k
m−k ( ZZ ) − Cm−k (Ω )kop .
n

(22)
(3k)

The last two inequalities hold because of the fact that the eigenvalues of Σ̃m−k and
−1
−1
C3k
m−k (Ω ) are in the interval [η , η], and Lemma 21. The following concentration
inequality of sample covariance matrix facilitates our proof.
12

Lemma 3. For the observations Z following certain sub-Gaussian distribution with constant ρ and precision matrix Ω, we have

1
log p + k
T
3k
−1 2
E max kC3k
.
m−k ( ZZ ) − Cm−k (Ω )kop ≤ C
n
n
m∈[p]
Lemma 3 is an extension of the result in Chapter 2 of [44]. Its proof can be found in
Lemma 3 of [13].
Combining Lemma 3, (21) and (22), we have
 
log p + k
k
3k
−1 −1 2
.
E max kΩ̃(k)
kop ≤ C
m − Ck+1 (Cm−k (Ω ))
n
m∈[p]

(23)

We turn to bounding the bias term of local estimator in (21).
Lemma 4. Assume that Ω ∈ Pα (η, M ) defined in (2) with α > 12 . Then we have

−1 −1
− Ckm (Ω)k2op ≤ Ck −2α+1 .
kCkk+1 (C3k
m−k (Ω ))
Lemma 4, together with (23), (21) and (20), implies that
EkΩ̃∗k − Ω∗k k2op ≤ C

log p + k
+ Ck −2α+1 .
n

(24)

Plugging Lemma 2 and (24) into (19), we finish the proof of Theorem 3.

3.2

Minimax Lower Bound under the Operator Norm over Pα (η, M )

Theorem 3 in Section 3.1 proves that the local cropping estimator defined in (12) attains
−2α+1
the convergence rate of n 2α + logn p . In this section, we establish the following matching
lower bound, which proves the rate optimality of the local cropping estimator.
Theorem 4. The minimax risk of estimating the precision matrix Ω over Pα (η, M ) defined in (2) under the operator norm with α ≥ 21 satisfies
inf

sup EkΩ̃ − Ωk2op ≥

Ω̃ Pα (η,M )

τ 2  − 2α−1 log p 
n 2α +
,
32
n

(25)

1

where 0 < τ < min{M, 41 η −1 , η 2 − 1}.
Remark 4. Theorems 3 and 4 together show that the minimax risk for estimating the
precision matrices over Pα (η, M ) stated in (4) of Theorem 1. It is worthwhile to notice
that there is no consistent estimator over Pα (η, M ) under the operator norm, when α ≤ 12 .
Proof. The lower bound of parameter space Pα (η, M ) can be established by the lower
bounds over its subsets. We construct two subsets P1 and P2 and calculate the lower
bound over those two subsets separately. Let τ be a positive constant which is less than
1
min{M, 41 η −1 , η 2 − 1}.
13

1

First, we construct P1 . Set k = min{dn 2α e, p2 }. Set the index set Θ = {0, 1}k , i.e.,
for any θ ≡ {θi }1≤i≤k ∈ Θ, each θi is either 0 or 1. Then we define the k × k matrix
1
A∗k (θ) ≡ [aij ]k×k with aij = τ n− 2 θi 1 (j = k) and


0k×k
0k×k
0k×(p−2k)


A(θ) =  A∗k (θ)
0k×k
0k×(p−2k)  .
0(p−2k)×k 0(p−2k)×k 0(p−2k)×(p−2k)
We then define P1 as the collection of 2k matrices indexed by Θ,

P1 = Ω(θ) : Ω(θ) = (Ip − A(θ))T (Ip − A(θ)), θ ∈ Θ .

(26)

Next, we construct P2 as the collection of the diagonal matrices in the following equation,
n
P2 = Ω(m) ≡ [wij (m)]p×p :
o
(27)
−1
1
,m ∈ 0 : p ,
wij (m) = 1 (i = j) + τ a 2 1 (i = j = m)
where a = min{ logn p , 1}.
Lemma 5. P1 and P2 are subsets of Pα (η, M ).
Note that we assume log p = O(n) and n = O(p). Without loss of generality, we
further assume log p < n < p. For any estimator Ω̃ based on n i.i.d. observations, we
establish the lower bounds over those two subsets in Sections 3.2.1 and 3.2.2 respectively,
sup EkΩ̃ − Ωk2op ≥
P1

sup EkΩ̃ − Ωk2op ≥
P2

1
τ 2 −1
p
τ 2 2α−1
n min{n 2α , } ≥ n− 2α ,
16
2
16

(28)

τ 2 −1
τ 2 log p
n min{log p, n} ≥
.
16
16 n

(29)

According to Lemma 5, (P1 ∪ P2 ) ⊂ Pα (η, M ). Therefore, we obtain
sup EkΩ̃ − Ωk2op ≥ max{sup EkΩ̃ − Ωk2op ,
P1

Pα (η,M )

≥

τ2 
32

n−

2α−1
2α

+

log p 
n

sup EkΩ̃ − Ωk2op }
P2

,

which completes the proof of Theorem 4.
We introduce some further notation before establishing (28) using Assouad’s lemma in
P
Section 3.2.1 and (29) using Le Cam’s method in Section 3.2.2. Let H(θ, θ0 ) = ki=1 |θi −θi0 |
be the Hamming distance on {0, 1}k , which is the number of different elements between θ
R
and θ0 . The total variation affinity kP ∧ Qk = p ∧ q dµ, where p and q are the density
functions of two probability measure P and Q with respect to any common dominating
measure µ.
14

3.2.1

Assouad’s lemma in proof of (28)

Assouad’s lemma [1] is a powerful tool to provide the lower bound over distributions
indexed by the hypercube Θ = {0, 1}k . Let Pθ be the distribution generated from observations indexed by Ω(θ). The proof of Lemma 6 can be found in [54], and thus omitted.
Lemma 6 (Assouad). Let Ω̃ be an estimator based on observations from a distribution
in the collection {Pθ , θ ∈ Θ}, where Θ = {0, 1}k . Then
sup 22 Eθ kΩ̃ − Ω(θ)k22 ≥
θ∈Θ

kΩ(θ) − Ω(θ0 )k22 k
H(θ, θ0 )
2
H(θ,θ0 )≥1
min

min

H(θ,θ0 )=1

kPθ ∧ Pθ0 k.

Applying the Assouad’s lemma to the subset P1 , we have the following results.
Lemma 7. Let Pθ be the joint distribution of n i.i.d. observations from N (0, Ω(θ)−1 ),
where Ω(θ) ∈ P1 defined in (26). Then
min

H(θ,θ 0 )=1

kPθ ∧ Pθ0 k ≥ 0.5.

Lemma 8. Consider all Ω(θ) ∈ P1 defined in (26). Then
1
kΩ(θ) − Ω(θ0 )k22
≥ (τ n− 2 )2 .
0
0
H(θ, θ )
H(θ,θ )≥1

min

1

Lemmas 6, 7 and 8 together imply the desired (28), with the choice k = dn 2α e. The
proofs of the above lemmas can be found in the supplement.
3.2.2

Le Cam’s method in proof of (29)

Le Cam’s method can be used to establish the lower bound via testing a single distribution
against a convex hull of distributions. Set r = inf m∈[p] kΩ(0) − Ω(m)k2op . Let Pi be the
distribution generated from observations indexed by Ω(i), where 0 ≤ i ≤ p. Define
P
P̄ = pm=1 Pm . The proof of the following lemma can be found in [54], and thus omitted.
Lemma 9 (Le Cam). Let Ω̃ be an estimator based on observations from a distribution in
the collection {Pi , 0 ≤ i ≤ p}. Then
1
sup EkΩ̃ − Ω(m)k2op ≥ rkP0 ∧ P̄ k.
2
0≤m≤p

1 2
τa 2
Applying Le Cam’s method to P2 , we obtain that r =
≥ 14 τ 2 a and the
1
1+τ a 2

following results.
Lemma 10. Let Pm be the joint distribution of n i.i.d. observations from N (0, Ω(m)−1 ),
where Ω(m) ∈ P2 defined in (27). Then
7
kP0 ∧ P̄ k > .
8
Combining the above results in Lemmas 9 and 10, we obtain the desired (29), i.e.,
sup EkΩ̃ − Ωk2op ≥

0≤m≤p

7 2
τ2
log p
τ a≥
min{
, 1}.
64
16
n
15

3.3
3.3.1

Rate Optimality under the Operator Norm over Qα (η, M )
Minimax upper bound

In this section, we establish the optimal rate of convergence over Qα (η, M ) under the
operator norm. Our estimation procedure remains the local cropping estimator, except
1
that the bandwidth is k  n 2α+1 , which is due to smaller bias terms.
1

Theorem 5. When dn 2α+1 e ≤ p, the local cropping estimator defined in (12) of the
precision matrix Ω over Qα (η, M ) given in (3) satisfies
sup
Qα (η,M )

2
−2α
EkΩ̃op
+C
k − Ωkop ≤ Ck

log p + k
,
n

1

When k = dn 2α+1 e, we have
sup
Qα (η,M )

2α

− 2α+1
2
EkΩ̃op
+C
k − Ωkop ≤ Cn

log p
.
n

Proof. We employ the same proof strategy as that of Theorem 3. Only two lemmas
bounding bias terms need to be replaced. We only emphasize the differences here.
We replace Lemma 2 in the proof by Lemma 11, which bounds the distance of the
population precision matrix and its tapered one.
Lemma 11. For Ω in the parameter space Qα (η, M ) defined in (3), Ω∗k is defined in (18),
we have
kΩ∗k − Ωk2op ≤ Ck −2α .
In addition, we replace Lemma 4 by Lemma 12, which bounds the bias term of each
local estimator.
Lemma 12. For Ω ∈ Qα (η, M ) defined in (2) with α > 0, we have

−1 −1
− Ckm (Ω)k2op ≤ Ck −2α .
kCkk+1 (C3k
m−k (Ω ))
The remaining part of the proof remains the same, including a similar upper bound
for the variance term stated in Lemma 3. Therefore, we complete our proof.
3.3.2

Minimax lower bound

Theorem 6. The minimax risk for estimating the precision matrix Ω over Qα (η, M )
defined in (3) under the operator norm with α > 0 satisfies
inf

sup

Ω̃ Qα (η,M )

EkΩ̃ − Ωk2op ≥

16

2α
τ 2 − 2α+1
log p
(n
+
).
32
n

(30)

Remark 5. Theorems 5 and 6 together show that the minimax risk for estimating the
precision matrices over Qα (η, M ) stated in (5) of Theorem 1. In contrast to Pα (η, M ),
the optimal rate of convergence over Qα (η, M ) is faster. In particular, rate-optimal local
cropping estimators are always consistent as long as α > 0.
Proof. To establish the lower bound for Qα (η, M ) in which the decay of aij is in the
entry-wise fashion, we repeat the proof scheme in Section 3.2 with a few changes. Let τ
1
be a positive constant which is less than min{M, 41 η −1 , η 2 − 1}.
1

Set k = min{dn 2α+1 e, p2 } and the index set Θ = {0, 1}k , i.e., for any θ ∈ Θ, θ ≡
{θi }1≤i≤k , each θi is either 0 or 1. Define the k × k matrix Bk∗ (θ) ≡ [bij ]k×k with bij =
1
τ (nk)− 2 θi . Define


0k×k
0k×k
0k×(p−2k)


B(θ) =  Bk∗ (θ)
0k×k
0k×(p−2k)  .
0(p−2k)×k 0(p−2k)×k 0(p−2k)×(p−2k)
We construct the collection of 2k matrices as

P3 = Ω(θ) : Ω(θ) = (Ip − B(θ))T (Ip − B(θ)), θ ∈ Θ .

(31)

Lemma 13. P3 is a subset of Qα (η, M ).
Let Pθ be the joint distribution of n i.i.d. observations from N (0, Ω(θ)−1 ), where
Ω(θ) ∈ P3 defined in (31). Parallel to Lemmas 7 and 8 and the lower bound (28) for P1 ,
we establish the following lower bound for P3 .
Lemma 14. Consider all Ω(θ) ∈ P3 defined in (31). Then
min

H(θ,θ 0 )=1

kPθ ∧ Pθ0 k ≥ 0.5,

1
kΩ(θ) − Ω(θ0 )k22
≥ (τ n− 2 )2 .
0
0
H(θ, θ )
H(θ,θ )≥1

min

(32)

(33)

According to Assouad’s lemma, for any estimator Ω̃ based on n i.i.d. observations, we
have
1
p
τ2
(34)
sup EkΩ̃ − Ωk2op ≥ n−1 min{n 2α+1 , }.
16
2
P3
It is easy to show (P3 ∪ P2 ) ⊂ Qα (η, M ), where P2 is defined in (27). Therefore,
combining (34) and (29), we complete the proof of Theorem 6.
Remark 6. The estimation of the covariance matrix Σ is of significant importance as
well. We propose the estimator of Σ by inverting our estimator Ω̃op
k given in (12). The
results and the analysis given in Section 3 can be used to establish the minimax optimality
−1 − Σk
of our estimator under the operator norm. According to the inequality k(Ω̃op
op ≤
k )
op −1
op
op
−1
−1
−1
k(Ω̃k ) kop kΩ̃k − Ωkop kΩ kop and the fact that both kΩ̃k ) kop and kΩ kop are
17

−1
bounded by η, we establish the upper bound of our estimator (Ω̃op
k ) . Furthermore, considering the analog between the covariance matrix and the precision matrix in the subset
P1 and P2 defined in (26) and (27), the matching lower bound can be proved by a similar
argument in Section 3.2. Therefore, we have the following rate optimality of estimating the
−1
covariance matrix under the operator norm, which can be achieved by estimator (Ω̃op
k ) ,

inf

+

log p
,
n

EkΩ̃−1 − Ω−1 k2op  n− 2α+1 +

log p
.
n

sup EkΩ̃−1 − Ω−1 k2op  n−

2α−1
2α

Ω̃ Pα (η,M )

inf

sup

Ω̃ Qα (η,M )

4

2α

Rate Optimality under the Frobenius Norm

In this section, we establish that the optimal rates of convergence for estimating the
precision matrix over the parameter spaces Pα (η, M ) and Qα (η, M ) are identical under
the Frobenius norm. Since Qα (η, αM ) ⊂ Pα (η, M ), one immediately obtain that
inf
Ω̃

1
1
EkΩ̃ − Ωk2F ≥ inf sup
EkΩ̃ − Ωk2F .
Ω̃ Qα (η,αM ) p
Pα (η,M ) p
sup

(35)

In order to show (6) in Theorem 2, it suffices to establish the upper bound over the parameter space Pα (η, M ) and the matching lower bound over the parameter space Qα (η, M ).
We assume that X follows the p-variate Gaussian distribution, with mean zero and precision matrix Ω in this section.

4.1

Minimax Upper Bound under the Frobenius Norm

In this section, we establish the following risk upper bound of the regression-based blockthresholding estimation procedure we proposed in Section 2.2 under the Frobenius norm
over Pα (η, M ).
1

Theorem 7. Assume dn 2α+2 e ≤ p. The estimator defined in (16) of the precision matrix
1
Ω over Pα (η, M ) and Qα (η, αM ) given in (2) and (3) with k = dn 2α+2 e satisfies
2α+1
1
1
− 2α+2
2
2
EkΩ̃F
EkΩ̃F
.
k − ΩkF ≤ sup
k − ΩkF ≤ Cn
Qα (η,αM ) p
Pα (η,M ) p

sup

(36)

Proof. We focus on the second inequality since the first one is trivial. Note that Ω̃F
k =
T
−1
T
−1
^
^
(I − A) D̃ (I − A) according to (16) while Ω = (I − A) D (I − A). The risk upper
bound can be controlled by bounding I^
− A − (I − A) and D̃ − D. To this end, we first
provide some properties of our estimator.
Lemma 15. Assume that X follows the p-variate Gaussian distribution with mean zero
and precision matrix Ω = (I −A)T D−1 (I −A), which belongs to parameter space Pα (η, M )
18

defined in (2). For any fixed i, di is the i-th diagonal of D, ai ∈ Ri−1 corresponds the i-th
row of the lower triangle in A. dˆi is defined in (15), and â∗i ∈ Ri−1 corresponds the i-th
row of the lower triangle in Â∗ defined in (14). Then we have
2α+1

E|dˆi − di |2 ≤ Cn− 2α+2 ,
2α+1

Ekâ∗i − ai k22 ≤ Cn− 2α+2 .
We are ready to establish the upper bounds of I^
− A − (I − A) and D̃ − D sepa−1
−1
rately. Note that kD̃ kop ≤ η and kD kop ≤ η, which is due to Lemma 22. Therefore,
Lemma 21 yields EkD̃ − Dk2F ≤ 4EkD̂ − Dk2F , which further implies that
1
1
EkD−1 − D̃−1 k2F ≤ EkD̃−1 k2op kD̃ − Dk2F kD−1 k2op
p
p
1
≤ 4η 4 EkD̂ − Dk2F
p
1X ˆ
E|di − di |2 .
≤ 4η 4
p
i

Together with Lemma 15, it follows that
2α+1
1
EkD−1 − D̃−1 k2F ≤ Cn− 2α+2 .
p

(37)

2α+1
− A) − (I − A)k2F ≤ Cn− 2α+2 . Lemma 21 implies
Next, we turn to prove that p1 Ek(I^

1
4
4X
Ek(I^
− A) − (I − A)k2F ≤ EkÂ∗ − Ak2F ≤
Ekâ∗i − ai k22 .
p
p
p
i

Combining above equation with Lemma 15, we have
2α+1
1
Ek(I^
− A) − (I − A)k2F ≤ Cn− 2α+2 .
p

(38)

At last, we derive the risk upper bound of our estimator. It is clear that kI^
− Akop ≤ η,
≤ η. According to Lemma 22, kI − Akop ≤ η, kD−1 kop ≤ η. Combining these
facts with (37) and (38), we have
kD̃−1 kop

1
3
− A) − (I − A)k2F
EkΩ̃ − Ωk2F ≤ E kI − Ak2op kD−1 k2op k(I^
p
p
+ kI − Ak2op kD−1 − D̃−1 k2F kI^
− Ak2op

+ k(I^
− A) − (I − A)k2F kD̃−1 k2op kI^
− Ak2op
1
1
≤6η 4 Ek(I^
− A) − (I − A)k2F + 3η 4 EkD−1 − D̃−1 k2F
p
p
2α+1

≤Cn− 2α+2 .
Therefore, we finish the proof of Theorem 7.
19

4.2

Minimax Lower Bound under the Frobenius Norm
2α+1

In this section, we establish the matching lower bound n− 2α+2 over parameter spaces
Pα (η, M ) and Qα (η, M ).
Theorem 8. The minimax risk for estimating the precision matrix Ω over Pα (η, M ) and
Qα (η, αM ) under the Frobenius norm satisfies
inf
Ω̃

1
τ 2 2α+1
1
EkΩ̃ − Ωk2F ≥ inf sup
EkΩ̃ − Ωk2F ≥ n− 2α+2 .
32
Ω̃ Qα (η,αM ) p
Pα (η,M ) p
sup

Remark 7. The minimax risk for estimating the precision matrices over Pα (η, M ) and
Qα (η, M ) under the Frobenius norm in Theorem 2 immediately follows from Theorems 7
and 8.
Proof. It is sufficient to establish the lower bound over Qα (η, M ) since the first inequality immediately follows from (35). We construct a least favorable subset in Qα (η, M ).
1
p
Without loss of generality, we assume 2k
is an integer where k = min{dn 2α+2 e, p2 }. Dekp

p
fine the index set Θ0 = {0, 1} 2 . For each θ ∈ Θ0 , we further denote it as 2k
many k 2
dimensional vectors, i.e., θ = {θ(s)}1≤s≤d p e , where θ(s)ij is equal to 0 or 1. For such an
2k
index θ, there is a corresponding p × p block diagonal matrix C(θ) such that each k × k
1
p
e. We set τ as a positive
block Cs (θ(s)) ≡ [c(s)ij ]k×k , where c(s)ij = τ n− 2 θ(s)ij , s ∈ d 2k
1
1 −1
constant which is less than min{M, 4 η , η 2 − 1}.


0k
0k


02k
...
02k
 C1 (θ(1)) 0k







0k
0k


0
.
.
.
0


2k
2k

.
C
(θ(2))
0
2
k
C(θ) = 



..
..
..
..


.
.
.
.





0k
0k 


02k
02k
...
p
Cd p e (θ(d 2k
e)) 0k
2k

Finally, we define the subset of Qα (η, M ) indexed by Θ0 as follows

P4 = Ω(θ) : Ω(θ) = (Ip − C(θ))T (Ip − C(θ)), θ ∈ Θ0 .

(39)

Lemma 16. P4 is a subset of Qα (η, M ).
Applying Lemma 6 to P4 , we obtain that,
inf max 22 Eθ kΩ̃ − Ω(θ)k2F ≥
Ω̃ θ∈Ω(Θ0 )

kΩ(θ) − Ω(θ0 )k2F kp
min kPθ ∧ Pθ0 k
H(θ, θ0 )
4 H(θ,θ0 )=1
H(θ,θ0 )≥1
min

20

(40)

Lemma 17. Let Pθ be the joint distribution of n i.i.d. observations from N (0, Ω(θ)−1 ),
where Ω(θ) ∈ P4 defined in (39). Then
min

H(θ,θ 0 )=1

and

kPθ ∧ Pθ0 k ≥ 0.5,

(41)

kΩ(θ) − Ω(θ0 )k2F
≥ τ 2 n−1 .
H(θ, θ0 )
H(θ,θ0 )≥1
min

(42)

Applying Lemma 17 into (40), we obtain
inf
Ω̃

n 1 po
1
1
τ2
EkΩ̃ − Ωk2F ≥ inf sup EkΩ̃ − Ωk2F ≥ n−1 min n 2α+2 ,
,
32
2
Ω̃ P4 p
Qα (η,M ) p
sup

which completes the proof of Theorem 8, noting that n < p.

5

An Extension to Nonparanormal Distributions

In this section, we extend the minimax framework to the nonparanormal model. Assume that X = (X1 , X2 , . . . , Xp )T follows the p-variate Gaussian distribution with covariance matrix Σ. Instead of n i.i.d. copies X1 , X1 , . . . , Xn of X, we only observe
their transformations. Specifically, we denote the transformed variables of X by Y =
(f1 (X1 ), f2 (X2 ), . . . , fp (Xp ))T , where each fi is some unknown strictly increasing function. Then our observation is Z = (Y1 , Y2 , . . . , Yn )T ∈ Rn×p , where each Yi is the
transformed Xi . This is a form of the Gaussian copula model [5], or the nonparanormal
model [36]. To avoid the identifiability issue, we set diag(Σ) = I, which makes Σ the
correlation matrix. Here we consider the same structural assumption as in previous sections on the inverse of the correlation matrix, which is denoted by Ω. Based on Pα (η, M )
and Qα (η, M ) defined in (2) and (3), the following two types of parameter spaces are of
interest,
(
)
diag(Ω−1 ) = I, Ω ∈ Pα (η, M );
0
Pα (η, M ) = {Ω, {fi }} :
,
(43)
fi is strictly increasing, i ∈ [p].
and

(
Q0α (η, M ) =

{Ω, {fi }} :

diag(Ω−1 ) = I,

)
Ω ∈ Qα (η, M );

fi is strictly increasing, i ∈ [p].

.

(44)

Our goal is to estimate the latent correlation structure, the inverse of the correlation
matrix Ω, using the observation Z. We establish the minimax risk of estimating Ω over
the parameter spaces Pα0 (η, M ) and Q0α (η, M ) under the operator norm in the following
theorem.
Theorem 9. Assume log p = O(n), n = O(p). Then for the nonparanormal model, the
minimax risk of estimating Ω under the operator norm over Pα0 (η, M ) with α > 21 satisfies
inf

sup

0 (η,M )
Ω̃ {Ω,{fi }}∈Pα

EkΩ̃ − Ωk2op  n−
21

2α−1
2α

+

log p
.
n

(45)

The minimax risk of estimating Ω under the operator norm over Q0α (η, M ) satisfies
inf

2α

sup

Ω̃ {Ω,{fi }}∈Q0α (η,M )

EkΩ̃ − Ωk2op  n− 1+2α +

log p
.
n

(46)

Rank-based estimator are widely applied in the nonparanormal model. Progress has
been made in this field during the last decade especially for high-dimensional statistics.
For instance, see [38] for bandable correlation matrix estimation, [4] for Gaussian graphical
models, and [23] for multi-task regression via Cholesky decomposition.

5.1

Estimation Procedure

We introduce our rate-optimal estimation procedure over the parameter spaces Pα0 (η, M )
and Q0α (η, M ) under the operator norm. The approach to estimate the inverse of the
correlation matrix in nonparanormal model is almost the same as the estimation scheme
of the precision matrix under the operator norm in Section 2.1, except that the sample covariance matrix needs to be replaced by its rank-based nonparametric variant via
Kendall’s tau (τ ) [31] or Spearman’s correlation coefficient rho (ρ) [46].
Kendall’s tau is defined as
X
2
sgn(Zk1 i − Zk2 i )sgn(Zk1 j − Zk2 j ).
τ̂ij =
n(n − 1)
1≤k1 <k2 ≤n

Then define

π
Σ̂τ = [sin( τ̂ij )]p×p .
2

(47)

Spearman’s rho is defined as
Pn
(rki − (n + 1)/2)(rkj − (n + 1)/2)
,
ρ̂ij = pPn k=1
Pn
2
2
k=1 (rki − (n + 1)/2)
k=1 (rkj − (n + 1)/2)
where rij is the rank of Zij among Z1j , Z2j , . . . , Znj . Define
π
Σ̂ρ = [2 sin( ρ̂ij )]p×p .
6

(48)

It is well-known that both Σ̂τ and Σ̂ρ are unbiased estimators of the population correlation
matrix Σ. We adopt almost the same estimation procedure proposed in Section 2.1, but
replacing n1 ZT Z in (10) with either Σ̂τ or Σ̂ρ . In this way, we construct the nonparametric
local cropping estimators Ω̃τk and Ω̃ρk in replace of Ω̃op
k in (12). Note that the optimal
choices of the bandwidth k are picked differently over two types of parameter spaces
Pα0 (η, M ) and Q0α (η, M ) as we did over Pα (η, M ) and Qα (η, M ) in Section 2.1.

5.2

Minimax Optimality

In this section, we prove that both Ω̃τk and Ω̃ρk achieve the minimax optimality under
the operator norm over the parameter space Pα0 (η, M ). The minimax optimality over
Q0α (η, M ) can be established in the same way.
22

First, we derive the risk lower bound over Pα0 (η, M ). Following the same strategy in
Section 3.2, we construct the two subset P10 and P20 based on P1 and P2 given in (26) and
(27). Define the subset


1
1

Ω = diag(Ω0−1 ) 2 Ω0 diag(Ω0−1 ) 2 , Ω0 ∈ P1 ;
,
(49)
P10 = {Ω, {fi }} :
1


fi (x) = diag(Ω0−1 )i2 x, i ∈ [p].
Lemma 18. P10 is a subset of Pα0 (η 2 , M η).
One can easily check that the probability measure of the Gaussian distribution with
precision matrix Ω and transformation {fi }, where {Ω, {fi }} ∈ P10 is equivalent to the the
probability measure of the Gaussian distribution with precision matrix Ω0 , where Ω0 ∈ P1 .
Therefore, by this one-to-one correspondence of probability measure between index sets
P10 and P1 , we immediately have
sup
0 (η,M )
{Ω,{fi }}∈Pα

EkΩ̃ − Ωk2op ≥

sup
{Ω,{fi }}∈P10

EkΩ̃ − Ωk2op = sup EkΩ̃ − Ωk2op .

(50)

Ω0 ∈P1

Applying the same proof strategy in Section 3.2, we have the following result.
1

1

Lemma 19. We set Ω = diag(Ω0−1 ) 2 Ω0 diag(Ω0−1 ) 2 for each Ω0 (θ) ∈ P1 in (26). Then it
holds that
kΩ(θ) − Ω(θ0 )k22
− 21 2
min
≥
(τ
n
) ,
H(θ, θ0 )
H(θ,θ0 )≥1
where 0 < τ < min{M η, 41 η −2 , η − 1}.
Then by Assoaud’s lemma, we obtain that
inf sup EkΩ̃ − Ωk2op ≥
Ω̃ Ω0 ∈P1

τ 2 − 2α−1
n 2α .
32
1

We can use a similar strategy to construct P20 . Note that in this case we need to put τ a 2
on the first sub-diagonal in I − A, instead of the diagonal of Σ. Combined with the result
from the Le Cam’s lemma on the subset P20 , we have
sup
0 (η,M )
{Ω,{fi }}∈Pα

EkΩ̃ − Ωk2op ≥ C(n−

2α−1
2α

+

log p
).
n

(51)

Next, we turn to the risk upper bound of our rank-based local cropping estimators. The
risk can be decomposed into the bias terms and the variance term in the same fashion
in Section 3.1. Since the bias terms are deterministic and only due to the bandable
structure of the Cholesky factor of the inverse correlation matrix, the upper bounds of
two bias terms we derived in Lemma 2 and Lemma 4 still hold. For the variance term, a
simple extension of Theorem 1 in [38] provides the following result.

23

Lemma 20. For any Ω such that diag(Ω−1 ) = I and Ω ∈ Pα (η, M ) ∪ Qα (η, M ), we have

log p + k
τ
3k
−1 2
E max k(C3k
,
m−k (Σ̂ ) − Cm−k (Ω )kop ≤ C
n
m∈[p]

log p + k
ρ
3k
−1 2
E max k(C3k
.
m−k (Σ̂ ) − Cm−k (Ω )kop ≤ C
n
m∈[p]
Replacing Lemma 3 by Lemma 20 and following the rest of the proof in Theorem 3,
we finally obtain that
sup

(EkΩ̃τk − Ωk2op + EkΩ̃ρk − Ωk2op ) ≤ Cn−

2α−1
2α

+C

0 (η,M )
{Ω,{fi }}∈Pα

log p
.
n

(52)

The lower bound (51) and upper bound (52) together give the optimal rate of convergence
(45) in Theorem 9. The optimal rate of convergence (46) in Theorem 9 can be proved similarly. Therefore, we establish the minimax framework for nonparanormal distributions.

6

Numerical Studies

In this section, we turn to the numerical performance of the proposed rate-optimal estimators under the operator norm for Pα (η, M ) and Qα (η, M ) defined in (2) and (3) to further
illustrate the fundamental difference of Pα (η, M ) and Qα (η, M ). In addition, we compare
them with the banding estimator proposed in [6], which is based on the auto-regression
between variables. Specifically, for a given bandwidth k < n, the banding estimator is
defined as Ω̃BL = (I − ÃBL )T (D̃BL )−1 (I − ÃBL ). Here the i-th row of the lower triangular
matrix ÃBL is the vector âi in (13), i.e., the least square estimates of the coefficients for
the regression of Xi against Xi−k:i−1 . The i-th entry of the diagonal matrix D̃BL is the
estimate of the residual variance for the regression of Xi against Xi−k:i−1 .

6.1

Simulation in Qα (η, M ) under the operator norm

We first focus on the parameter space Qα (η, M ) and compare the performance of local
cropping estimator and the banding estimator. Specifically, we generate the precision
matrix in the following form:
Ω = (I − A)T D−1 (I − A),

A ≡ [aij ]p×p ,

D = Ip ,

where aij = −(i − j)−α−1 when i > j; otherwise aij = 0. It is easy to check that
Ω ∈ Qα (η, 1) with some large η > 0. The simulation is done with a range of parameter
values for p, n, α. Specifically, the decay rate α ranges from 0.5 to 2 with a step of 0.5,
the sample size n ranges from 500 to 4000, the dimension p ranges from 500 to 2000.
In this setting, we compare our local cropping estimator (denoted as cropping.Q.)
with the banding estimator (denoted as BL) proposed in [6]. According to [6], the bandwidth of banding estimator is chosen as k  (n/ log p)1/(2α+2) . The optimal bandwidth
24

over Qα (η, M ) is k  n1/(2α+1) . In the simulation, the bandwidth of BL estimator is
b(n/ log p)1/(2α+2) c and the bandwidth of crop.Q is bn1/(2α+1) c.
Table 1 reports the average errors of the banding estimator (BL) and local cropping
estimator (crop.Q) under the operator norm over 100 replications. The smaller errors in
each experiment are highlighted in boldface. Figure 2 displays the boxplots of the errors
of BL and crop.Q.
It can be seen from Table 1 that crop.Q outperforms BL in most cases with a few exceptions when n is small. As the sample size increases, the average errors of both methods
decrease, which matches our intuition. In addition, the dimension p has minor effect on
the errors of both estimators, which is partially reflected by the optimal rates (dominating
2α
term n− 2α+1 ) obtained in Theorem 1. For each fixed dimension p, the superiority crop.Q
over BL becomes more significant as the sample size n increases, which implies that BL
estimator is indeed sub-optimal.

6.2

Simulation in Pα (η, M ) under the operator norm

We demonstrate the fundamental difference between two types of parameter space Pα (η, M )
and Qα (η, M ) by numerical studies in this section. Of note, although local cropping estimators proposed in (12) are rate-optimal over both Pα (η, M ) and Qα (η, M ), the corresponding optimal choices of bandwidth are distinct. We generate precision matrices in
the following way to guarantee that Ω is always in Pα (η, M ) but not in Qα (η, M ) with
some fixed η and M . Considering
Ω = (I − A)T D−1 (I − A),

A ≡ [aij ]p×p ,

D = Ip ,

where the first column of A is ai1 = −2(i − 1)−α , 2 ≤ i ≤ p. The remaining entries are
all zeros. It is easy to check that Ω ∈ Pα (η, 2) with some large η > 0. The simulation
is carried out with a similar range of values for p, n, α as in Section 6.1. Note that the
consistent estimator exists only if α > 0.5. Therefore, in this setting, the decay rate α
varies among 1, 1.5 and 2.
The optimal choice of bandwidth of local cropping estimator over Pα (η, M ) is k 
1
2α
n , which is different from the one of crop.Q. We denote this rate-optimal estimator in
1
Pα (η, M ) by crop.P. In the simulation, the bandwidth of crop.P is bn 2α c. We also include
BL estimator as a reference.
Table 2 reports the average errors of the three procedures, crop.P, crop.Q and BL,
under the operator norm over 100 replications. The smallest errors in each experiment are
highlighted in boldface. Figure 3 plots the boxplots of their errors for p = 500, 1000, 2000.
Since Ω always belongs to Pα (η, M ) but not Qα (η, M ), the estimator crop.Q is suboptimal and thus expected to have an inferior performance. Table 2 shows this point,
i.e., for fixed p and α, the advantage of crop.P is more obvious as n increases. Especially,
crop.P outperforms the other two estimators when n = 4000. We also see a similar pattern
as in Table 1 that p has minor effect on the errors of all the estimators.
25

Table 1: The average errors under the operator norm of the banding estimator (BL) and
the local cropping estimator (crop.Q) over 100 replications.
α = 0.5
α=1
α = 1.5
α=2
p
n
crop.Q BL crop.Q BL crop.Q BL crop.Q BL
500
4.68
5.44
1.64
2.38
1.18
1.16
0.93
0.81
1000
3.29
4.89
1.17
1.72
0.82
1.08
0.66
0.69
500
2000
2.47
4.45
0.89
1.33
0.59
0.69
0.48
0.59
4000
1.84
3.80
0.62
1.07
0.41
0.64
0.34
0.53
500
4.96
5.74
1.75
2.40
1.30
1.19
0.99
0.84
1000
3.43
5.19
1.24
1.74
0.86
1.10
0.68
0.70
1000
2000
2.58
4.75
0.93
1.35
0.62
0.71
0.51
0.60
4000
1.93
4.10
0.66
1.33
0.44
0.65
0.36
0.55
500
5.14
5.97
1.85
2.41
1.33
1.21
1.06
0.89
1000
3.58
5.41
1.30
1.76
0.90
1.12
0.72
0.71
2000
2000
2.69
4.97
0.98
1.37
0.65
0.73
0.54
0.62
4000
2.01
4.32
0.69
1.34
0.45
0.66
0.38
0.55

Table 2: The average errors under the operator norm of the banding estimator (BL) and
the local cropping estimator (crop.P & crop.Q) with two different bandwidths over 100
replications.
α=1
α = 1.5
α=2
p
n
crop.P crop.Q BL crop.P crop.Q BL crop.P crop.Q BL
500
1.50
1.18
2.32 0.66
0.73
0.86 0.52
0.65
0.53
1000
1.09
0.96
1.80 0.47
0.56
0.83 0.38
0.56
0.45
500
2000
0.83
0.80
1.53 0.35
0.43
0.55 0.27
0.32
0.41
4000 0.64
0.68
1.33 0.26
0.35
0.54 0.19
0.24
0.38
500
1.50
1.20
2.36 0.68
0.74
0.91 0.57
0.68
0.59
1000
1.12
0.98
1.82 0.49
0.58
0.81 0.39
0.55
0.46
1000
2000
0.84
0.81
1.54 0.37
0.44
0.55 0.27
0.32
0.41
4000 0.65
0.68
1.52 0.26
0.35
0.53 0.19
0.24
0.38
500
1.51
1.21
2.39 0.69
0.75
0.96 0.62
0.71
0.63
1000
1.16
1.00
1.81 0.51
0.60
0.84 0.39
0.56
0.46
2000
2000
0.85
0.81
1.55 0.39
0.44
0.56 0.27
0.33
0.41
4000 0.65
0.69
1.70 0.26
0.35
0.53 0.19
0.24
0.38

26

Figure 2: The boxplot of the errors from the local cropping estimator with the optimal
bandwidth in Qα (η, M ) (cropping.Q) and the banding estimator (BL) over 100 replications.

References
[1] Patrice Assouad. Deux remarques sur l’estimation. Comptes rendus des séances de
l’Académie des sciences. Série 1, Mathématique, 296(23):1021–1024, 1983.
[2] Onureena Banerjee, Laurent El Ghaoui, and Alexandre dAspremont. Model selection
through sparse maximum likelihood estimation for multivariate Gaussian or binary
data. Journal of Machine learning research, 9(Mar):485–516, 2008.
[3] Sayantan Banerjee and Subhashis Ghosal. Posterior convergence rates for estimating

27

Figure 3: The boxplot of the errors from the local cropping estimator with the optimal
bandwidth in Pα (η, M ) (cropping.P), the local cropping estimator with the optimal bandwidth in Qα (η, M ) (cropping.Q) and the banding estimator (BL) over 100 replications.
large precision matrices using graphical models. Electronic Journal of Statistics,
8(2):2111–2137, 2014.
[4] Rina Foygel Barber and Mladen Kolar. Rocket: Robust confidence intervals via
Kendall’s tau for transelliptical graphical models. arXiv preprint arXiv:1502.07641,
2015.
[5] Peter J Bickel, Chris AJ Klaassen, Peter J Bickel, Y Ritov, J Klaassen, Jon A Wellner,
and YA’Acov Ritov. Efficient and adaptive estimation for semiparametric models.
Johns Hopkins University Press Baltimore, 1993.
28

[6] Peter J Bickel and Elizaveta Levina. Covariance regularization by thresholding. The
Annals of Statistics, 36(6):2577–2604, 2008.
[7] Peter J Bickel and Elizaveta Levina. Regularized estimation of large covariance
matrices. The Annals of Statistics, 36(1):199–227, 2008.
[8] Jacob Bien, Florentina Bunea, and Luo Xiao. Convex banding of the covariance
matrix. Journal of the American Statistical Association, 111(514):834–845, 2016.
[9] T Tony Cai. Minimax and adaptive inference in nonparametric function estimation.
Statistical Science, 27(1):31–50, 2012.
[10] T Tony Cai, Weidong Liu, and Harrison H Zhou. Estimating sparse precision matrix:
Optimal rates of convergence and adaptive estimation. The Annals of Statistics,
44(2):455–488, 2016.
[11] T Tony Cai, Zhao Ren, and Harrison H Zhou. Optimal rates of convergence for estimating Toeplitz covariance matrices. Probability Theory and Related Fields, 156(12):101–143, 2013.
[12] T Tony Cai, Zhao Ren, and Harrison H Zhou. Estimating structured high-dimensional
covariance and precision matrices: Optimal rates and adaptive estimation. Electronic
Journal of Statistics, 10(1):1–59, 2016.
[13] T Tony Cai, Cun-Hui Zhang, and Harrison H Zhou. Optimal rates of convergence
for covariance matrix estimation. The Annals of Statistics, 38(4):2118–2144, 2010.
[14] T Tony Cai and Harrison H Zhou. Optimal rates of convergence for sparse covariance
matrix estimation. The Annals of Statistics, 40(5):2389–2420, 2012.
[15] Tony Cai and Weidong Liu. Adaptive thresholding for sparse covariance matrix
estimation. Journal of the American Statistical Association, 106(494):672–684, 2011.
[16] Tony Cai, Weidong Liu, and Xi Luo. A constrained `1 minimization approach to
sparse precision matrix estimation. Journal of the American Statistical Association,
106(494):594–607, 2011.
[17] I Csiszár. Information-type indices of the divergence of distributions. i, ii. Magyar
Tud. Akad. Mat. Fiz. Oszt. Közl. 17 (1967), 123-149; ibid, 17:267–291, 1967.
[18] Alexandre d’Aspremont, Onureena Banerjee, and Laurent El Ghaoui. First-order
methods for sparse covariance selection. SIAM Journal on Matrix Analysis and Applications, 30(1):56–66, 2008.
[19] Bernard Delyon and Anatoli Juditsky. On minimax wavelet estimators. Applied and
Computational Harmonic Analysis, 3(3):215–228, 1996.
29

[20] Michael B Eisen, Paul T Spellman, Patrick O Brown, and David Botstein. Cluster
analysis and display of genome-wide expression patterns. Proceedings of the National
Academy of Sciences, 95(25):14863–14868, 1998.
[21] Noureddine El Karoui. On the largest eigenvalue of Wishart matrices with identity
covariance when n, p and p/n tend to infinity. arXiv preprint math/0309355, 2003.
[22] Noureddine El Karoui. Operator norm consistent estimation of large-dimensional
sparse covariance matrices. The Annals of Statistics, 36(6):2717–2756, 2008.
[23] Jianqing Fan, Lingzhou Xue, and Hui Zou. Multitask quantile regression under the
transnormal model. Journal of the American Statistical Association, 111(516):1726–
1735, 2016.
[24] Reinhard Furrer and Thomas Bengtsson. Estimation of high-dimensional prior and
posterior covariance matrices in Kalman filter variants. Journal of Multivariate Analysis, 98(2):227–255, 2007.
[25] Thomas M Hamill, Jeffrey S Whitaker, and Chris Snyder. Distance-dependent filtering of background error covariance estimates in an ensemble kalman filter. Monthly
Weather Review, 129(11):2776–2790, 2001.
[26] Mark H Heyer and F Peter Schloerb. Application of principal component analysis to
large-scale spectral line imaging studies of the interstellar medium. The Astrophysical
Journal, 475(1):173, 1997.
[27] Addison Hu and Sahand Negahban. Minimax estimation of bandable precision matrices. In Advances in Neural Information Processing Systems, pages 4893–4901,
2017.
[28] Jianhua Z Huang, Naiping Liu, Mohsen Pourahmadi, and Linxu Liu. Covariance matrix selection and estimation via penalised normal likelihood. Biometrika, 93(1):85–
98, 2006.
[29] Iain M Johnstone. On the distribution of the largest eigenvalue in principal components analysis. The Annals of Statistics, pages 295–327, 2001.
[30] Iain M Johnstone and Arthur Yu Lu. On consistency and sparsity for principal components analysis in high dimensions. Journal of the American Statistical Association,
104(486):682–693, 2009.
[31] Maurice G Kendall. A new measure of rank correlation. Biometrika, 30(1/2):81–93,
1938.
[32] Clifford Lam and Jianqing Fan. Sparsistency and rates of convergence in large covariance matrix estimation. The Annals of Statistics, 37(6B):4254, 2009.
30

[33] Olivier Ledoit and Michael Wolf. Improved estimation of the covariance matrix of
stock returns with an application to portfolio selection. Journal of Empirical Finance,
10(5):603–621, 2003.
[34] Kyoungjae Lee and Jaeyong Lee. Estimating large precision matrices via modified
Cholesky decomposition. arXiv preprint arXiv:1707.01143, 2017.
[35] Elizaveta Levina, Adam Rothman, and Ji Zhu. Sparse estimation of large covariance
matrices via a nested lasso penalty. The Annals of Applied Statistics, 2(1):245–263,
2008.
[36] Han Liu, John Lafferty, and Larry Wasserman. The nonparanormal: Semiparametric
estimation of high dimensional undirected graphs. Journal of Machine Learning
Research, 10(Oct):2295–2328, 2009.
[37] Nicolai Meinshausen and Peter Bühlmann. High-dimensional graphs and variable
selection with the lasso. The Annals of Statistics, pages 1436–1462, 2006.
[38] Ritwik Mitra and Cun-Hui Zhang. Multivariate analysis of nonparametric estimates
of large correlation matrices. arXiv preprint arXiv:1403.6195, 2014.
[39] Nikhil Padmanabhan, Martin White, Harrison H Zhou, and Ross O’Connell. Estimating sparse precision matrices. Monthly Notices of the Royal Astronomical Society,
460(2):1567–1576, 2016.
[40] Debashis Paul. Asymptotics of sample eigenstructure for a large dimensional spiked
covariance model. Statistica Sinica, 17(4):1617, 2007.
[41] Pradeep Ravikumar, Martin J Wainwright, Garvesh Raskutti, and Bin Yu. Highdimensional covariance estimation by minimizing `1 -penalized log-determinant divergence. Electronic Journal of Statistics, 5:935–980, 2011.
[42] Zhao Ren, Tingni Sun, Cun-Hui Zhang, and Harrison H Zhou. Asymptotic normality
and optimalities in estimation of large Gaussian graphical models. The Annals of
Statistics, 43(3):991–1026, 2015.
[43] Adam J Rothman, Peter J Bickel, Elizaveta Levina, and Ji Zhu. Sparse permutation
invariant covariance estimation. Electronic Journal of Statistics, 2:494–515, 2008.
[44] Leonas Saulis and VA Statulevicius. Limit theorems for large deviations, volume 73.
Springer Science & Business Media, 2012.
[45] Juliane Schäfer and Korbinian Strimmer. A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics. Statistical applications in genetics and molecular biology, 4(1), 2005.

31

[46] Charles Spearman. Spearmans rank correlation coefficient. Amer J Psychol, 15:72–
101, 1904.
[47] Tingni Sun and Cun-Hui Zhang. Sparse matrix inversion with scaled lasso. The
Journal of Machine Learning Research, 14(1):3385–3418, 2013.
[48] Roman Vershynin. Introduction to the non-asymptotic analysis of random matrices.
arXiv preprint arXiv:1011.3027, 2010.
[49] Kenneth W Wachter. Probability plotting points for principal components. In Ninth
Interface Symposium Computer Science and Statistics, pages 299–308. Prindle, Weber
and Schmidt, Boston, 1976.
[50] Kenneth W Wachter. The strong limits of random matrix spectra for sample matrices
of independent elements. The Annals of Probability, 6(1):1–18, 1978.
[51] Wei Biao Wu and Mohsen Pourahmadi. Nonparametric estimation of large covariance
matrices of longitudinal data. Biometrika, 90(4):831–844, 2003.
[52] Wei Biao Wu and Mohsen Pourahmadi. Banding sample autocovariance matrices of
stationary processes. Statistica Sinica, pages 1755–1768, 2009.
[53] Luo Xiao and Florentina Bunea. On the theoretic and practical merits of the banding
estimator for large covariance matrices. arXiv preprint arXiv:1402.0844, 2014.
[54] Bin Yu. Assouad, Fano, and Le Cam. Festschrift for Lucien Le Cam, 423:435, 1997.
[55] Ming Yuan. High dimensional inverse covariance matrix estimation via linear programming. Journal of Machine Learning Research, 11(Aug):2261–2286, 2010.
[56] Ming Yuan and Yi Lin. Model selection and estimation in the Gaussian graphical
model. Biometrika, 94(1):19–35, 2007.

32

7

Appendix

Theorem 1 immediately follows from Theorems 3, 4, 5 and 6 while Theorem 2 immediately
follows from Theorems 7 and 8. We provide key lemmas in this section in the proofs of
those main theorems. Recall that for 1 ≤ q ≤ ∞, the matrix `q norm of a matrix
S is defined by kSkq = maxkxkq =1 kSxkq . We use the following two facts repeatedly
p
√ F ≤ kSkop ≤
in this section: (i) kSk
kSk1 kSk∞ ; (ii) for a real symmetric matrix S,
p
kSkop ≤ kSk1 .

7.1

Preliminary Lemmas: Lemma 21, Lemma 22, and Lemma 23

In this section, we provide three preliminary lemmas and their proofs, which are important
in the proofs of the other lemmas.
7.1.1
Lemma 21. Operator Pη (·) is defined in (9). For any square matrix A ≡ [aij ]p×p such
that η −1 ≤ λmin (A) ≤ λmax (A) ≤ η, we have,
kA − Pη (S)kop ≤ 2kA − Skop ,

(53)

kA − Pη (S)kF ≤ 2kA − SkF .

(54)

Proof. Since
kA − Pη (S)k∗ ≤ ||S − Pη (S)||∗ + ||A − S||∗ ,
we only need to prove that
||S − Pη (S)||∗ ≤ ||A − S||∗ ,
where || · ||∗ is either the operator norm or Frobenius norm.
For an asymmetric matrix S, assume the singular value decomposition of S is U T SV =
W . Since operator norm and Frobenius norm are invariant to the orthogonal transformation, it is sufficient to prove
||U T SV − U T Pη (S)V ||∗ ≤ ||U T AV − U T SV ||∗ .
Here, U T SV and U T Pη (S)V are the diagonal matrices. η −1 ≤ λmin (U T AV ) ≤ λmax (U T AV ) ≤
η. Without loss of generality, we only need to prove that for any diagonal matrix
W ≡ diag(w), where w = (w1 , . . . , wp )T ,
||W − Wη ||∗ ≤ ||A − W ||∗ ,

33

where Wη = diag(max{min{wi , η}, η −1 }) and η −1 ≤ λmin (A) ≤ λmax (A) ≤ η. For the
operator norm, note the fact that,
kA − W kop =

kAv − W vk2

sup
v∈Rp ,kvk2 =1

≥ max kAei − W ei k2
i

≥ max max kAei k2 − kW ei k2 , kW ei k2 − kAei k2
i

≥ max max η −1 − wi , wi − η, 0
i

= max |wi − max{min{wi , η}, η −1 }|
i

= kWη − W kop .
where ei denotes the vector all of whose components are zero, except the i-th component
being one. For the Frobenius norm, we have
X
X
kA − W k2F =
a2ij +
(aii − wi )2
i

i6=j

XX
=
(
a2ij + wi2 − 2aii wi )
i

j

sX
XX
2
2
a2ij wi )
≥
(
aij + wi − 2
i

=

X

j

j

(wi −

sX

i

a2ij )2 .

j

qP
2
Since η −1 ≤ λmin (A) ≤ λmax (A) ≤ η, we obtain η −1 ≤
j aij ≤ η. It is easy to check
qP
2
−1
2
2
that for any i, (
j aij − wi ) ≥ (max{min{wi , η}, η } − wi ) . Then we have
kA − W k2F ≥ kWη − W k2F ,
that yields kA − W kF ≥ kWη − W kF .
The same result can be derived easily in the case of symmetric matrix S using similar
argument on the eigen-decomposition. Thus, we finish the proof.
7.1.2
Lemma 22. Assume precision matrix Ω satisfies that η −1 ≤ λmin (Ω) ≤ λmax (Ω) ≤ η.
Let Σ = Ω−1 and Ω = (I − A)T D−1 (I − A) (see (1)), then,
η −1 ≤ λmin (Σ) ≤ λmax (Σ) ≤ η,
η

−1

≤ λmin (D) ≤ λmax (D) ≤ η,

η −1 ≤ λmin (I − A) ≤ λmax (I − A) ≤ η.
34

(55)
(56)
(57)

Proof. Since Σ = Ω−1 , it is trivial that η −1 ≤ λmin (Σ) ≤ λmax (Σ) ≤ η. For the symmetric positive-definite matrices Ω ≡ [ωij ]p×p and Σ ≡ [σij ]p×p , ωii ≤ λmax (Ω) ≤ η; σii ≤
λmax (Σ) ≤ η for any i ∈ [p]. ωii−1 can be interpreted as the variance of the residual of
the regression of Xi on X−i = (X1 . . . Xi−1 , Xi+1 . . . Xp )T . σii can be interpreted as the
variance of the residual of the regression of Xi on 0. Therefore, the i-th entry of D, the
variance of the residual of the regression of Xi on X1:i−1 = (X1 . . . Xi−1 )T , has the bound
ωii−1 ≤ di ≤ σii , then η −1 ≤ λmin (D) ≤ λmax (D) ≤ η.
For any unit vector u, v such that (I − A)u = λv, λ > 0, we have
uT Ωu = uT (I − A)T D−1 (I − A)u = λ2 v T D−1 v,
we obtain that λ2 = (uT Ωu)/(v T D−1 v) ∈ [η −2 , η 2 ], that yields
η −1 ≤ λmin (I − A) ≤ λmax (I − A) ≤ η.
7.1.3
Lemma 23. Assume X is an i-variate random vector with covariance matrix Σ such that
η −1 ≤ λmin (Σ) ≤ λmax (Σ) ≤ η. Let the linear projection of Xi onto Xi−k:i−1 in population
hki
be X̂i . With the corresponding coefficients padded with i − k − 1 zeros in the front, we
hki
hki
hki
can rewrite it as X̂i = XT1:i−1 β i , where β i ∈ Ri−1 with its first i − k − 1 coordinates
hki

hki

being zeros. In addition, set i = Var(Xi − X̂i ).
hi−1i
hi−1i
hi−1i T
(i) Whenever β i
= βi1 , . . . , βi(i−1) satisfies
X

hi−1i

|βij

| < M k −α ,

k ∈ [i − 1],

(58)

j<i−k

we have,
hi−1i

hki

− β i k2 ≤ 2η 2 M k −α ,

kβ i

hi−1i

|i
hi−1i

(ii) Whenever β i

hi−1i

= βi1

hki

− i | ≤ 4η 4 M k −α .

hi−1i T

, . . . , βi(i−1)

hi−1i

|βij

(59)
(60)

satisfies

| < M (i − j)−α−1 ,

k ∈ [i − 1],

(61)

it holds that,
hi−1i

kβ i

1

hki

− β i k2 ≤ 2η 2 M (k − 1)−α− 2 ,

hi−1i
|i

−

hki
i |

4

−α− 12

≤ 4η M (k − 1)

35

.

(62)
(63)

hi−1i

hi−1i

hi−1i

hki

hki

hki

Proof. For any fixed i and k, let β i
= ((β A )T , (β B )T )T and β i = ((β A )T , (β B )T )T ,
hi−1i
hki
hi−1i
hki
where the sizes of β A
and β A are i − k − 1, the sizes of β B
and β B are k. Note
hki
that β A = 0. To facilitate the proof, we divide Σ into several block matrices,


Σ11 Σ12 Σ13


Σ = Σ21 Σ22 Σ23  ,
Σ31 Σ32 Σ33
where Σ11 , Σ22 , Σ33 are the covariance matrices of X1:i−k−1 , Xi−k:i−1 and Xi respectively.
By the definition of the linear projection, we have
hi−1i

Σ23 = Σ21 β A

hi−1i

+ Σ22 β B

hi−1i T

hi−1i = (−(β i

,

hi−1i T

) , 1)Σ(−(β i

) , 1)T ,

hki

Σ23 = Σ22 β B ,
hki

hki

hki = (−(β i )T , 1)Σ(−(β i )T , 1)T .
Condition in (58) implies
hi−1i

hi−1i

kβ A

k2 ≤ kβ A

hki

hi−1i

k1 ≤ M k −α .

Then we have,
hi−1i

kβ i

− β i k2 ≤ kβ A

hi−1i

≤ kβ A

hi−1i

− 0k2 + kβ B

hki

− β B k2

hi−1i

k2 + kΣ−1
22 Σ21 β A

k2

≤ M k −α + η 2 M k −α
≤ 2η 2 M k −α .
that yields the desired (59).
hi−1i T
Assume the decomposition of Σ−1 is (I − A)T D−1 (I − A). Note that (−(β i
) , 1)
corresponds to the ith row of the lower triangle of I − A. According to Lemma 21,
hi−1i T
kI − Akop ≤ η, then k(−(β i
) , 1)T k2 ≤ kI − Akop ≤ η. Applying the same argument
hki
on the covariance matrix of Xi−k:i , we have k(−(β i )T , 1)T k2 ≤ η. Moreover,
hi−1i

|i

hki

− i |

hi−1i T

hi−1i T

hki

hki

) , 1)T − (−(β i )T , 1)Σ(−(β i )T , 1)T |

hi−1i T
hki
hi−1i
hki
≤ k(−(β i
) , 1)T k2 + k(−(β i )T , 1)T k2 kΣkop kβ i
− β i k2
≤|(−(β i

) , 1)Σ(−(β i

≤4η 4 M k −α ,
that yields the desired (60).
Similarly, condition in (61) implies
hi−1i

kβ A

hi−1i

k2 ≤ kβ A

1

k1 ≤ M (k − 1)−α− 2 .

36

Then we have
hi−1i

kβ i

1

hki

− β i k2 ≤ 2η 2 M (k − 1)−α− 2 ,

hi−1i

|i

1

hki

− i | ≤ 4η 4 M (k − 1)−α− 2 ,

which completes the proofs of (62) and (63).

7.2

Proofs of Lemma 2 and Lemma 4 in analysis of Theorem 3

In this section, we prove Lemma 2 and Lemma 4 to establish Theorem 3.
7.2.1

Proof of Lemma 2

For Ω ≡ [ωij ]p×p , define
bdk (Ω) ≡ [ωij 1 (|i − j| ≤ k)]p×p .
P
It is easy to check Ω∗k = k1 2k−1
i=k bdi (Ω). Then we have
kΩ −

Ω∗k kop

(64)

2k−1
1 X
≤
(kΩ − bdi (Ω)kop ).
k
i=k

1

We turn to the analysis of kΩ − bdk (Ω)kop . Define D− 2 (I − A) as B ≡ [bij ]. We know
Ω = B T B, and
X
1
max
|bij | ≤ M η 2 k −α .
i

j<i−k

Set bdk (B) as Bk , and B − bdk (B) as B−k . Then we can rewrite Ω as
Ω = BT B
= (Bk + B−k )T (Bk + B−k )
T
T
B−k .
Bk + BkT B−k + B−k
= BkT Bk + B−k

Checking the entries of BkT Bk , we find that bdk BkT Bk = BkT Bk . Thus,

kΩ − bdk (Ω)kop

(65)

=kB T B − bdk (B T B)kop
T
T
T
T
=kB−k
Bk + BkT B−k + B−k
B−k − bdk (B−k
Bk + BkT B−k + B−k
B−k )kop
T
T
T
T
≤kB−k
Bk + BkT B−k + B−k
B−k kop + kbdk (B−k
Bk + BkT B−k + B−k
B−k )kop
T
T
T
T
=kB−k
B + B T B−k − B−k
B−k kop + kbdk (B−k
Bk + BkT B−k + B−k
B−k )kop
T
T
≤kB−k
Bkop + kB T B−k kop + kB−k
B−k kop
T
T
+ kbdk (B−k
Bk )kop + kbdk (BkT B−k )kop + kbdk (B−k
B−k )kop .

The key is to control kB−k kop . In addition, in order to get rid off the operator
T B )k , we adopt a technique which requires
bdk (·) when handling the term kbdk (B−k
−k op
controlling the operator norm of the matrix in which each entry is the absolute value of
the entry in B−k . To this end, we state the key lemma below.
37

Lemma 24. Assume that matrix B−k is defined above. We have
1

+
kB−k
kop ≤ Ck 2 −α ,

where X + is the matrix in which each entry is the magnitude of the corresponding entry
in X.
Proof. Assume Fi is the matrix composed by the (2i−1 k + 1)-th to (2i k)-th sub-diagonals
Pdlog2 (p/k)e
+
in matrix B + , i ∈ [dlog2 (p/k)e]. Note that B−k
=
(Fi ). For any i ≥ 1,
i=1
1

kFi k∞ ≤ M η 2 (2i−1 k)−α . Since there are at most 2i−1 k entries in each column of Fi ,
1
kFi k1 ≤ M η 2 (2i−1 k)1−α . The operator norm of Fi can be bounded by the two terms
above,
1

kFi kop ≤ (kFi k∞ kFi k1 ) 2
1

1

≤ M η 2 (2i−1 k) 2 −α
1

1

1

= M η 2 k 2 −α × (2 2 −α )i−1 .
Consequently,
dlog2 (p/k)e
+
kop ≤
kB−k

X

kFi kop

i=1
dlog2 (p/k)e

≤

X

1

1

1

(M η 2 k 2 −α × (2 2 −α )i−1 )

i=1
1

dlog2 (p/k)e

1

= M η 2 k 2 −α ×

X

1

(2 2 −α )i−1

i=1
1
2

≤ CM η k

1
−α
2

.

The bounds of the operator norms of many other matrices can be derived from the
above result combining the following lemma.
Lemma 25. Let X be a square matrix, then
kXkop ≤ kX + kop ,
and
kbdk (X)kop ≤ kX + kop ,
where X + is the matrix in which each entry is the magnitude of the corresponding entry
in X.
Proof. Assume X ≡ [xij ]p×p , let kXkop = uT Xv, then
X
X
kXkop = uT Xv =
ui xij vj ≤
|ui ||xij ||vj | = (u+ )T X + v + ≤ kX + kop .
i,j

i,j

38

Assume bdk (X) = [x∗ij ]p×p , kbdk (X)kop = uT (bdk (X))v, then
kbdk (X)kop = uT (bdk (X))v =

X

ui x∗ij vj ≤

i,j

X

|ui ||xij ||vj | ≤ kX + kop .

i,j
1

+
It follows from the previous two lemmas that kB−k kop ≤ kB−k
kop ≤ Ck 2 −α . Then we
have
1

kB T B−k kop ≤ kBkop kB−k kop ≤ Ck 2 −α ,
1

T
kB−k
Bkop ≤ kBkop kB−k kop ≤ Ck 2 −α ,

(66)

T
kB−k
B−k kop ≤ kB−k kop kB−k kop ≤ Ck 1−2α .

In addition, Lemma 25 implies that,
+ T +
T
T
) B−k kop ≤ Ck 1−2α .
B−k )+ kop ≤ k(B−k
B−k )kop ≤ k(B−k
kbdk (B−k

(67)

T B )k +kbd (B T B )k . We control kbd (B T B )k
Then we turn to bound kbdk (B−k
k
k op
k
−k k ∞
k −k op
T B )k first. For h ∈ [p], we have
and kbdk (B−k
1
k
T
T
B))k1
Bk )k∞ ≤ max krowh (bdk (B−k
kbdk (B−k

≤

h
h+k
X

h+2k
X

(

|bjh ||bij |)

i=h+1 j=h+k+1

=

h+2k
X

(

h+k
X

|bjh ||bij |)

j=h+k+1 i=h+1

≤

h+2k
X

M 2 η((j − h)−α (j − h − k)−α )

j=h+k+1

=

k
X

M 2 η((j + k)−α j −α )

j=1
2

≤M ηk

−α

k
X
j=1

≤Ck

1−2α

,

39

j −α

and
T
T
kbdk (B−k
Bk )k1 ≤ max kcolh (bdk (B−k
B))k1

≤

h
h−1
X

(

h+k
X

|bji ||bjh |)

i=h−k j=h+1

≤

h+k
X

(

h−1
X

|bji ||bjh |)

j=h+1 i=h−k

≤

h+k
X

M 2 η(k −α (j − h)−α )

j=h+1

≤

k
X

M 2 η(k −α j −α )

j=1

≤ Ck 1−2α .
Therefore, we obtain that
T
T
T
Bk )k1 )1/2 ≤ Ck 1−2α ,
Bk )k∞ kbdk (B−k
Bk )kop ≤ (kbdk (B−k
kbdk (B−k

kbdk (BkT B−k )kop ≤ (kbdk (BkT B−k )k∞ kbdk (BkT B−k )k1 )1/2 ≤ Ck 1−2α .

(68)

Combining (65), (66), (67) and (68), we prove that
1

kΩ − bdk (Ω)kop ≤ Ck 2 −α ,
which follows kΩ − Ω∗k k2op ≤ Ck 1−2α .
7.2.2

Proof of Lemma 4

The proof strategy in this lemma is not complicated, although the notation is quite in
−1 −1
volved. Our target is to bound the distance between Ckm (Ω) and Ckk+1 (C3k
m−k (Ω ))
under the operator norm. To this end, we introduce an intermediate term Ckm (Ω)∗ to facilitate our proof. Specifically, we break the target into two terms kCkm (Ω) − Ckm (Ω)∗ k2op and

−1 −1 k2 and derive their bounds respectively. The construckCkm (Ω)∗ −Ckk+1 (C3k
op
m−k (Ω ))
tion of Ckm (Ω)∗ with corresponding Cholesky decomposition are illustrated in Figure 4,

−1 −1 .
in contrast with those of Ckm (Ω) and Ckk+1 (C3k
m−k (Ω ))
To express the decomposition of Ckm (Ω) in the matrix form, we define the p × k matrix
Rpm,k ≡ [rij ]p×k ,

rij = 1 (i − m = j − 1) .

Assume Ω = (I − A)T D−1 (I − A). Set (I − A)Rpm,k as G. One can check
Ckm (Ω) = (Rpm,k )T ΩRpm,k
= (Rpm,k )T (I − A)T D−1 (I − A)Rpm,k
= GT D−1 G.
40

Ckm (Ω)

=

×

D−1

×

Ckm (Ω)∗

=
GT Rpm,2k (Rpm,2k )T×

D−1

×Rpm,2k (Rpm,2k )T G

Ckm (Ω)∗

=

−1 −1
=
Ckk+1 ((C3k
m−k (Ω )) )

GT

G

HT

×

−1
C2k
m (D )

×

H

KT

×

−1
C2k
k+1 (E )

×

K

Figure 4: An illustration of the proof strategy in Lemma 4.
Define
Ckm (Ω)∗ = (Rpm,k )T (I − A)T Rpm,2k (Rpm,2k )T D−1 Rpm,2k (Rpm,2k )T (I − A)Rpm,k
= GT Rpm,2k (Rpm,2k )T D−1 Rpm,2k (Rpm,2k )T G.
We first bound kCkm (Ω) − Ckm (Ω)∗ k2op . Since I − A is a lower triangular matrix,
consists of the first 2k columns of G. Then we have

Rpm,2k (Rpm,2k )T G

G − Rpm,2k (Rpm,2k )T G = (0, 0, . . . , 0, gm+2k , . . . , gp )T ,
where gi = rowi (G), and for i ∈ (m + 2k) : p,
kgi k2 ≤ kgi k1 ≤ M (i − m − k + 1)−α .

(69)

Consequently,
kG − Rpm,2k (Rpm,2k )T Gk2op
≤k(kgm+2k k2 , . . . , kgp k2 )T k22
≤k(M (k + 1)−α , . . . , M (p − m − k)−α )T k22
≤M 2 k −2α+1 .
Then we have
kCkm (Ω) − Ckm (Ω)∗ k2op
−1

(70)

≤k(G D G − G Rpm,2k (Rpm,2k )T D−1 Rpm,2k (Rpm,2k )T Gk2op
≤kG − Rpm,2k (Rpm,2k )T Gk2op kD−1 k2op (kGk2op + kRpm,2k (Rpm,2k )T Gk2op )
≤2kD−1 k2op kGk2op kG − Rpm,2k (Rpm,2k )T Gk2op
2
2 −2α+1
T

≤2η M k

T

.

41


−1 −1 k2 . Assume
Next, we turn to derive the bound of kCkm (Ω)∗ − Ckk+1 (C3k
op
m−k (Ω ))
−1 ))−1 = (I − B)T E −1 (I − B). One can also check
that (C3k
(Ω
m−k

−1 −1
Ckk+1 (C3k
m−k (Ω ))
k+1,k T
−1 −1 k+1,k
=(R3k
) (C3k
m−k (Ω )) R3k
k+1,k
k+1,k T
) (I − B)T E −1 (I − B)R3k
=(R3k
k+1,k T
k+1,2k
k+1,2k T −1 k+1,2k
k+1,2k T
k+1,k
=(R3k
) (I − B)T R3k
(R3k
) E R3k
(R3k
) (I − B)R3k
.
k+1,2k T
k+1,2k
To ease our notation, one can check (Rpm,2k )T DRpm,2k = C2k
) ER3k
=
m (D) and (R3k
m,2k T
m,k
k+1,2k T
k+1,k
2k
Ck+1 (E). In addition, we set (Rp ) (I − A)Rp
as H, (R3k
) (I − B)R3k
k
∗
T
2k
−1
k
as K. Then we can rewrite Cm (Ω) = H Cm (D )H and kCm (Ω) − Ckm (Ω)∗ k2op =
−1
2
2k
−1
2k
−1 2
K T C2k
k+1 (E )K. We bound kH − Kkop and kCm (D ) − Ck+1 (E )kop separately
below.
Referring to the instruction in Figure 4, one can check rowi (H) is part of the coefficients
of the regression Xm+i ∼ X1:m+i−1 and rowi (K) is part of the coefficients of the regression
Xm+i ∼ Xm−k:m+i−1 . According to Lemma 23, rowi (H −K) can be bounded by 2M η 2 (k+
i)−α . The dimension of the non-zero part of H − K is 2k × k. So we have

kH − Kk2op ≤ 8η 4 M 2 k −2α+1 .
2k
Similarly, the i-th elements of C2k
m (D) and Ck+1 (E) are the residuals of the above two
2k
2
8
2 −2α . Therefore,
regressions. According to Lemma 23, kC2k
m (D) − Ck+1 (E)kop ≤ 16η M k
we have
−1 2
−1
2k
kC2k
m (D ) − Ck+1 (E )kop
−1 2
2
2k
−1 2
2k
2k
≤kC2k
m (D )kop kCm (D) − Ck+1 (E)kop kCk+1 (E )kop

≤16M 2 η 12 k −2α .
Combing the above two results, we have

−1 −1 2
kop
kCkm (Ω)∗ − Ckk+1 (C3k
m−k (Ω ))
−1
2
−1
T 2k
≤kH T C2k
m (D )H − K Ck+1 (E )Kkop
−1
2k
−1 2
2
≤kH T k2op kC2k
m (D ) − Ck+1 (E )kop kHkop
−1 2
2
+ (kHk2op + kKk2op )kC2k
k+1 (E )kop kH − Kkop

≤η 4 × 16M 2 η 12 k −2α + 6η 4 × 8η 4 M 2 k −2α+1
≤96M 2 η 16 k −2α+1 .
In the end, based on the Equations (70) and (71), we have

−1 −1 2
kCkm (Ω) − Ckk+1 (C3k
kop
m−k (Ω ))

−1 −1 2
≤2kCkm (Ω) − Ckm (Ω)∗ k2op + 2kCkm (Ω)∗ − Ckk+1 (C3k
kop
m−k (Ω ))
≤2 × (2η 2 M 2 k −2α+1 + 96M 2 η 16 k −2α+1 )
≤200M 2 η 16 k −2α+1 .
We finish the proof of Lemma 4.
42

(71)

7.3

Proofs of Lemma 5, Lemma 7, Lemma 8, and Lemma 10 in the
analysis of Theorem 4

In this section, we prove Lemma 5, Lemma 7, Lemma 8, and Lemma 10 to establish
Theorem 4.
7.3.1

Proof of Lemma 5

First we prove that P1 ∈ Pα (η, M ). Let A∗k (θ) ≡ [aij ]k×k .
P
each entry. It is easy to check i−j>k |aij | ≤ M k −α . One
specific form of

Ik + (A∗k (θ))T A∗k (θ) −(A∗k (θ))T

Ω(θ) = 
−A∗k (θ)
Ik
0(p−2k)×k
0(p−2k)×k

We know the exact value of
can check Ω(θ) ∈ P1 has the

0k×(p−2k)

0k×(p−2k)  .
Ip−2k

Let 1 denote the vector with all 1’s in Θ. One can check that
2
2
λmax (Ω(θ)) = λmax (I − A(θ)) ≤ λmax (I + A(θ)) ≤(1 + λmax (A∗k (1)))2
1

1

=(1 + k 2 n− 2 τ )2 ≤ η.
The second inequality above is due to that the entries of A(θ) are all non-negative and
Lemma 25. Recall Σ(θ) = (I + A(θ))(I + A(θ))T . Thus, we can check
−1
−2
λmin (Ω(θ)) = λmax (Σ(θ))
= λmax (I + A(θ)) ≥(1 + λmax (A∗k (1)))−2
1

1

=(1 + k 2 n− 2 τ )−2 ≥ η −1 .
The eigenvalues of Ω(θ) are in the interval [η −1 , η]. So P1 ∈ Pα (η, M ).
Then we turn to prove that P2 ∈ Pα (η, M ). The Cholesky factor A of Ω(m) is the
1
zero matrix. The minimum eigenvalue of Ω(m) is (1 + τ a 2 )−1 , which is greater than η −1
and maximum one is 1, which is less than η. So P2 ∈ Pα (η, M ).
7.3.2

Proof of Lemma 7

Since kPθ ∧ Pθ0 k = 1 − 21 kPθ − Pθ0 k1 , we turn to control maxH(θ,θ0 )=1 kPθ − Pθ0 k1 . The
following First Pinsker’s inequality will facilitate our analysis.
Lemma 26 (First Pinsker’s Inequality [17]).
1
kPθ − Pθ0 k21 ≤ K(Pθ0 |Pθ )
2
n 1
1
p
= [ tr(Σ(θ0 )Σ(θ)−1 ) − log det(Σ(θ0 )Σ(θ)−1 ) − ]
2 2
2
2
where K(·|·) is the Kullback-Leibler divergence.

43

One can check Σ(θ) ∈ P1 has the form of
Σ(θ) = (I + A(θ))(I + A(θ)T )

Ik
(A∗k (θ))T
 ∗
=  (Ak (θ)) Ik + (A∗k (θ))(A∗k (θ))T
0(p−2k)×k
0(p−2k)×k


0k×(p−2k)

0k×(p−2k)  .
Ip−2k .

Set Σ(θ0 ) = D + Σ(θ) and d = θ0 − θ. Note that whenever max H(θ, θ0 ) = 1, D has the
form of
D = (A(θ0 ) − A(θ)) + (A(θ0 ) − A(θ))T + (A(θ0 )A(θ0 )T − A(θ)A(θ)T )
= A(d) + A(d)T + A(θ0 )A(d)T + A(d)A(θ)T ,
where A(d) is similarly defined as A(θ) except that θi is replaced by di . Let 1 denote the
all-ones vector, one can check
1

kDkop ≤ 2(1 + kA∗k (1)kop )kA∗k (d)kop ≤ 4τ n− 2 ,
kDkF ≤ 2(1 + kA∗k (1)kF )kA∗k (d)kF ≤ 4τ n
1

− 21

.

(72)
(73)

1

Furthermore, we have kDΣ(θ)−1 kop ≤ η4τ n− 2 ≤ n− 2 , kDΣ(θ)−1 kF ≤ kΣ(θ)−1 kop kDkF ≤
1
1
η × 4τ n− 2 ≤ n− 2 . One can easily check that Σ(θ0 )Σ(θ)−1 = I + DΣ(θ)−1 , 12 tr(I +
P
DΣ(θ)−1 ) = p2 + 21 tr(DΣ(θ)−1 ), log det(I +DΣ(θ)−1 ) = tr(DΣ(θ)−1 )+ i (log(1+λi )−λi ),
where λi ’s are the eigenvalues of DΣ(θ)−1 . Applying the First Pinsker’s inequality in this
case, we have
n 1
1
p
[ tr(Σ(θ0 )Σ(θ)−1 ) − log det(Σ(θ0 )Σ(θ)−1 ) − ]
2 2
2
2
n 1
1
p
−1
−1
= [ tr(I + DΣ(θ) ) − log det(I + DΣ(θ) ) − ]
2 2
2
2
nX
=
(λi − log(1 + λi )).
2

kPθ − Pθ0 k21 ≤

i

1

Since |λi | ≤ kDΣ(θ)−1 kop , all the λi are bounded by ±n− 2 . By Taylor expansion,
P 2
−1 2
i (λi − log(1 + λi )) ≤ 2
i λi = 2kDΣ(θ) kF , we have

P

kPθ − Pθ0 k21 ≤ nkDΣ(θ)−1 k2F ≤ 1,
which implies kPθ ∧ Pθ0 k ≥ 0.5.

44

7.3.3

Proof of Lemma 8

The proof is as follows,
kΩ(θ0 ) − Ω(θ)kop
≥

sup

(uT , 0T , 0T )(Ω(θ0 ) − Ω(θ))(0, v, 0)T

kuk2 =1,kvk2 =1

≥

sup
kuk2 =1,kvk2 =1

uT (A∗k (θ) − A∗k (θ0 ))T v

=k(A∗k (θ0 ) − A∗k (θ))kop
1

=(H(θ0 , θ))1/2 τ n− 2 ,
which immediately implies,
1
kΩ(θ) − Ω(θ0 )k22
≥ (τ n− 2 )2 .
0
0
H(θ, θ )
H(θ,θ )≥1

min

7.3.4

Proof of Lemma 10

Denote the density functions of Pi and P̄ by fi and f¯ respectively, where 0 ≤ i ≤ p. It is
R ¯2
R ¯2
1
sufficient to bound ff0 du − 1 because kP0 ∧ P̄ k ≥ 1 − 21 ( ff0 du − 1) 2 . To this end, note
that
Z ¯2
Z 2
X Z fi fj
fk
f
1 X
du − 1 = 2 (
du +
du) − 1
f0
p
f0
f0
1≤i6=j≤p

k∈[p]

n
1
(p(1 − τ 2 kn−1 )− 2 + p2 − p) − 1
2
p
1 2
1
≤
τ k≤ ,
2p
16

=

which further implies kP0 ∧ P̄ k ≥ 1 −

7.4

1
2

×

1
4

= 78 .

Proofs of Lemma 11 and Lemma 12 in the analysis of Theorem 5

In this section, we prove Lemma 11 and Lemma 12 to establish Theorem 5.
7.4.1

Proof of Lemma 11
P
Recall Ω∗k = k1 2k−1
h=k bdh (Ω), where bdh (·) is defined in Equation (64). Then,
kΩ∗k − Ωkop ≤

2k−1
1 X
kbdh (Ω) − Ωkop .
k
h=k

45

For any fix k ≤ h < 2k, define Uh ≡ [ωij 1 (i − j > h)]p×p , then bdh (Ω) − Ω = Uh + UhT .
Note that Ω = (I − A)T D−1 (I − A), where I − A ≡ [a0ij ]p×p . We have
kUh k∞ ≤ max krowi (Uh )k1 ≤ max η
i

i

≤ max ηM
i

≤ CηM k

p
X

p
X

|a0si |

s=i

i−k
X

|a0sj |

j=1

(s − i)−α−1 (s − i + k)−α

s=i
−α

.

Similarly, we have
kUh k1 ≤ max kcolj (Uh )k1 ≤ max ηM
j

j

≤ max ηM
j

p
X

p
X

|a0sj |

s=j+k

(s − j)−α−1

s=j+k

≤ CηM k −α .
With the bounds for kUh k1 and kUh k∞ , we have
kΩ∗k − Ωkop ≤

2k−1
1 X
kbdi (Ω) − Ωkop
k
i=k

≤ max kbdh (Ω) − Ωkop
k≤h<2k

≤ max (kUh kop + kUhT kop )
k≤h<2k
p
≤ max 2 kUh k∞ kUh k1
k≤h<2k

≤ CηM k −α .
7.4.2

Proof of Lemma 12

The proof of Lemma 12 is basically the same as the one of Lemma 4, except for a few
steps, which are highlighted below.
Since Ω ∈ Qα (η, M ), (69) should be updated by
kgi k2 ≤ M (i − m − k + 1)−α−1/2 ,
and consequently (70) should be replaced by
kCkm (Ω) − Ckm (Ω)∗ k2op ≤ 2M 2 k −2α .
According to Lemma 23, (71) is replaced by
−1 −1 2
2 16 −2α
kCkm (Ω)∗ − Ckk+1 ((C3k
.
m−k (Ω )) )kop ≤ 96M η k

Therefore, we have
−1 −1 2
2 16 −2α
kCkm (Ω) − Ckk+1 ((C3k
.
m−k (Ω )) )kop ≤ 200M η k

46

7.5

Proofs of Lemma 13 and Lemma 14 in the analysis of Theorem 6

In this section, we prove Lemma 13 and Lemma 14 to establish Theorem 6.
7.5.1

Proof of Lemma 13

Let Bk∗ (θ) ≡ [bij ]k×k , then we have the specific value of each entry. It is easy to check
bij ≤ M (i − j)−α−1 . One can check Ω(θ) ∈ P3 has the specific form of


Ik + Bk∗ (θ)T Bk∗ (θ) −Bk∗ (θ)T 0k×(p−2k)


Ω(θ) = 
(74)
−Bk∗ (θ)
Ik
0k×(p−2k)  .
0(p−2k)×k
0(p−2k)×k
Ip−2k
Let 1 denote the vector with all 1’s in Θ. Then one can check that
2
2
λmax (Ω(θ)) = λmax (I − B(θ)) ≤ λmax (I + B(θ)) ≤(1 + λmax (Bk∗ (1)))2
1

1

=(1 + k 2 n− 2 τ )2 ≤ η.
The second inequality above follows from that the entries of A(θ) are all non-negative and
Lemma 25. Recall Σ(θ) = (I + B(θ))(I + B(θ))T . Therefore, we have
−1
−2
λmin (Ω(θ)) = λmax (Σ(θ))
= λmax (I + B(θ)) ≥(1 + λmax (Bk∗ (1)))−2
1

1

=(1 + k 2 n− 2 τ )−2 ≥ η −1 .
The eigenvalues of Ω(θ) are in the interval [η −1 , η]. Therefore, P3 ∈ Pα (η, M ).
7.5.2

Proof of Lemma 14

The proof follows most part of the one of Lemma 7. Some inequalities in the proof
of Lemma 7 need to be rechecked. Assume that Σ(θ) = {Ω(θ)−1 : Ω(θ) ∈ P3 } and
D = Σ(θ0 ) − Σ(θ), one can verify that D has the same decomposition as that in the proof
of Lemma 7 except that A∗k is replaced by Bk∗ . We can show
1
1
1
kDΣ(θ)−1 kop ≤ ηkDk∞ ≤ η(τ (nk)− 2 k + (τ (nk)− 2 )2 k 2 ) ≤ ,
2
1

1

1

kDΣ(θ)−1 k2F ≤ η 2 (2k(τ (nk)− 2 )2 + (2k − 1)((τ (nk)− 2 )2 k)2 ) ≤ 4η 2 (τ n− 2 )2 .
Thus, (32) still holds in this case. As for (33), similarly,
1

1

kΩ(θ0 ) − Ω(θ)k2op ≥ kBk∗ (θ0 ) − Bk∗ (θ)k2op = H(θ, θ0 )(τ (nk)− 2 )2 k = H(θ, θ0 )(τ n− 2 )2 ,
1
kΩ0 (θ) − Ω0 (θ0 )k22
≥ (τ n− 2 )2 .
0
0
H(θ, θ )
H(θ,θ )≥1

min

Plugging the results of (32) and (33) into Lemma 6, we have
inf sup EkΩ̃ − Ωk2op ≥
Ω̃

P3

1
τ 2 −1
τ2
p
n k = n−1 min{n 2α+1 , }.
16
16
2

47

7.6

Proof of Lemma 15 in the analysis of Theorem 7

In this section, we prove Lemma 15 to establish Theorem 7.
7.6.1

Proof of Lemma 15
hki

Let the linear projection of Xi onto the linear span of Xi−k:i−1 in population be X̂i .
With the corresponding coefficients padded with i−k −1 zeros in the front, we can rewrite
hki
hki
hki
it as X̂i = XT1,i−1 ai , where ai ∈ Ri−1 with its first i − k − 1 coordinates being zeros.
hki

hki

hi−1i

hi−1i

hk i

In addition, we set di = Var(Xi − X̂i ). Note that ai = ai
, di = di
. Let âi 1
hk
i
and dˆi 1 be the empirical coefficient and residual of the regression Xi ∼ Xi−k1 :i−1 , where
hk i∗
k1 = d nc e with some sufficiently large c > 1. The coefficients with threshold âi 1 is
hk i∗
hk i
defined as (14). According to (14) and (15), we know â∗ = â 1 , dˆi = dˆ 1 .
i

2α+1

i

i

First, we prove that E|dˆi − di |2 ≤ Cn− 2α+2 . To this end, we decompose it as follow,
hk i
hi−1i 2
E|dˆi − di |2 = E|dˆi 1 − di
|
hk i
hk i
hk i
hi−1i 2
≤ 2E|dˆi 1 − di 1 |2 + 2|di 1 − di
| .
hk i

2α+1

hi−1i

|2 ≤ Cn−2α ≤ Cn− 2α+2 , noting that
According to Lemma 23 we know that |di 1 − di
hk i hk i
k1 = d nc e and α > 21 . Besides, the regression theory implies that, (n − k1 )dˆi 1 /di 1 ∼
χ2 (n − k1 ). So we have
hk1 i

P(|dˆi

hk1 i

/di

− 1| > t) ≤ 2 exp(−(n − k1 )t2 /8),

t ∈ (0, 1).

Then one can check that
hk1 i

E|dˆi

2α+1

hk1 i 2

| ≤ Cn− 2α+2 .

− di

2α+1

hk i
Therefore, we have shown that E|dˆi 1 − di |2 ≤ Cn− 2α+2 by combining the above two
equations together.
2α+1
Then we turn to prove Ekâ∗i − ai k22 ≤ Cn− 2α+2 .
hk1 i∗

Ekâ∗i − ai k22 =Ekâi

hi−1i 2
k2
hk1 i 2
ai k2

− ai

hk1 i∗

≤2Ekâi

−

hi−1i

According to Lemma 23, we know that kai
hk1 i∗

α > 21 . It is sufficient to prove that Ekâi

−

hk i
âi 1

hi−1i

+ 2kai

hk1 i 2
k2 .

− ai

hk1 i 2
k2 ≤ Cn−2α
2α+1
hk i
ai 1 k22 ≤ Cn− 2α+2 .

− ai

2α+1

≤ Cn− 2α+2 , since

We focus on the regression coefficients
first. The following analysis is conditioned
on Zi−k1 :i−1 . Note that with probability one that (ZTi−k1 :i−1 Zi−k1 :i−1 )−1 exists since k1 =
dn/ce. For any fixed i, we have
hk1 i

âi

hk1 i

|Zi−k1 :i−1 ∼ N (ai

, (ZTi−k1 :i−1 Zi−k1 :i−1 )−1 Var(Xi |Zi−k1 :i−1 )).

48

hk1 i

For each coordinate in âi

hk i

hk i

hk i

1
1
= (0, âi(i−k
. . . , âi(i−1)
)T , âij 1 |Zi−k1 :i−1 with j ∈ i−k1 : i−1
1)

hk i

can be represented as aij 1 + σj , where σj follows the normal distribution with variance
which can be bounded as follows since Ω ∈ Pα (η, M ),
Var(σj ) ≤ k(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop Var(Xi |Zi−k1 :i−1 )

(75)

≤ ηk(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop .
Next, we can show,
hk i∗

|âij 1



1
hk i
hk i 3
− aij 1 | ≤ min{|aij 1 |, λj } + |3σj |1 |σj | > λj .
2
2

(76)

To see this, one can check that
hk i∗

|âij 1





hk i
hk i
hk i
hk i
− aij 1 | = max{|σj |1 |aij 1 + σj | > λj , |aij 1 |1 |aij 1 + σj | ≤ λj }.

For the first term, we have


hk i
|σj |1 |aij 1 + σj | > λj


hk i
≤|σj |1 (|σj | > λj /2) + |σj |1 |σj | ≤ λj /2 ∩ |aij 1 | > λj /2
hk i

≤|σj |1 (|σj | > λj /2) + min{|aij 1 |, λj /2}.
Similarly, for the second term, we have


hk i
hk i
|aij 1 |1 |aij 1 + σj | ≤ λj


hk i
hk i
≤|aij 1 |1 |aij 1 + σj | ≤ λj ∩ |σj | > λj /2


hk i
hk i
+ |aij 1 |1 |aij 1 + σj | ≤ λj ∩ |σj | ≤ λj /2
hk i

≤|3σj |1 (|σj | > λj /2) + min{|aij 1 |, 3λj /2}.
We finish the proof of Equation (76).
Equation (76) further implies
hk1 i∗

E(kâi
≤2

i−1
X

hk1 i 2
k2 |Zi−k1 :i−1 )

− ai

hk i

min{|aij 1 |, 3/2λj }2 + 2

j=i−k1

≤4

i−1
X

hi−1i

|, 3/2λj }2 + 18

j=i−k1

i−1
X
j=i−k1

i−1
X

E(σj2 1 (|σj | > λj /2) |Zi−k1 :i−1 )

j=i−k1

hk1 i

≤6

E (3σj )2 1 (|σj | > λj /2) |Zi−k1 :i−1



j=i−k1

min{|aij

+ 4kai

i−1
X

hi−1i 2
k2

− ai
hi−1i

|aij

|λj + 18

i−1
X

2α+1

E(σj2 1 (|σj | > λj /2) |Zi−k1 :i−1 ) + Cn− 2α+2 .

j=i−k1

49

(77)

Set J0 as logk20 , J1 as logk21 . Due to that Ω ∈ Pα (η, M ) and (75), we can show
i−1
X

hi−1i

|aij

|λj

j=i−k1

≤

JX
1 +1 p

(k − J0 )R × M (2−(k−J0 )α − 2−(k+1−J0 )α )2−J0 α

k=J0

≤n

−1/2 −J0 α

2

q
M 8ηkn(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop

J1 −J
X0 +1

√ 
2−kα k

k=0
− 2α+1
2α+2

≤Cn

q

kn(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop
2α+1

≤C(kn(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop + 1)n− 2α+2 ,
and
i−1
X

E(σj2 1 (|σj | > λj /2) |Zi−k1 :i−1 )

j=i−k1

≤

i−1
X

E(ηk(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop 1 (|σj | > λj /2) |Zi−k1 :i−1 )

j=i−k1

≤ηk(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop

i−1
X

P(|σj | > λj /2|Zi−k1 :i−1 )

j=i−k1

≤n−1 Ckn(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop

i−1
X

exp(−

j=i−k1

≤n−1 Ckn(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop

JX
1 +1

λ2j
8ηk(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop

)

exp(−(k − J0 ))(2k )

k=J0
2α+1

≤Ckn(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop n− 2α+2 .
Plugging the above two equations in (77), we have
hk1 i∗

hk1 i 2
k2 )
hk i∗
hk i
=E(E(kâi 1 − ai 1 k22 |Zi−k1 :i−1 ))

E(kâi

− ai

2α+1

≤C(Ekn(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop + 1)n− 2α+2
2α+1

2α+1

+ CEkn(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop n− 2α+2 + Cn− 2α+2
2α+1

≤Cn− 2α+2 .
The last inequality holds since Ekn(ZTi−k1 :i−1 Zi−k1 :i−1 )−1 kop can be bounded by constant,
noting that k1 = dn/ce with some sufficiently large c > 1. Indeed, one can follow the
strategy in the proof of Theorem 5 in [13], together with the concentration inequality of
50

the sample covariance matrix under the operator norm (e.g., from Theorem 5.39 in [48]).
It follows that
2α+1
E(kâ∗i − ai k22 ) ≤ Cn− 2α+2 .

7.7

Proofs of Lemma 16 and Lemma 17 in the analysis of Theorem 8

In this section, we prove Lemma 16 and Lemma 17 to establish Theorem 8.
7.7.1

Proof of Lemma 16

Let Cs (θ(s)) ≡ [c(s)ij ]k×k , it is easy to check
Ω(θ) in P4 has the specific form of

E1 (θ(1))
02k

E2 (θ(2))
 02k
Ω(θ) = 
.
..

..
.

02k
02k

P

i−j>k

|cij | ≤ M k −α−1 . One can check

...
...
..
.



02k
02k
..
.



,


p
p
. . . Ed e (θ(d 2k e))
2k

where,
"
#
Ik + Cs (θ(s))T Cs (θ(s)) −Cs (θ(s))T
Es =
.
−Cs (θ(s))
Ik
Let 1 denote the vector with all 1’s in Θ. Then one can check
2
2
λmax (Ω(θ)) = λmax (I − C(θ)) ≤ λmax (I + C(θ)) ≤(1 + λmax (Cs (1)))2
1

=(1 + kn− 2 τ )2 ≤ η.
The second inequality above is due to fact that the entries of C(θ) are all non-negative
and Lemma 25. Recall Σ(θ) = (I + C(θ))(I + C(θ))T . Therefore, we have
−1
−2
λmin (Ω(θ)) = λmax (Σ(θ))
= λmax (I + C(θ)) ≥(1 + λmax (Cs (1)))−2
1

=(1 + kn− 2 τ )−2 ≥ η −1 .
The eigenvalues of Ω(θ) are in the interval [η −1 , η]. So P4 ∈ Pα (η, M ).
7.7.2

Proof of Lemma 17

For (41), the proof is almost the same as that of Lemma 7. Since Lemma 17 is about the
distributions in the subset P4 , we need to recheck the inequalities (72) and (73). Assume
that Σ(θ) = {Ω(θ)−1 : Ω(θ) ∈ P4 } and D = Σ(θ0 ) − Σ(θ), note that H(θ, θ0 ) = 1, in this
case, one can verify that we only need to replace A∗k by Cs . Specifically, Cs corresponds
to the one where θ and θ0 are different. Then we still have
1

kDkop ≤ 2(1 + kCs (1)kop )kCs (d)kop ≤ 4τ n− 2 ,
1

kDkF ≤ 2(1 + kCs (1)kF )kCs (d)kF ≤ 4τ n− 2 .
51

So (41) holds in this case. As for (42), similarly,
X
1
kΩ(θ0 ) − Ω(θ)k2F ≥
kCs (θ0 (s) − θ(s))k2F = H(θ, θ0 )(τ n− 2 )2 ,
s

kΩ0 (θ) − Ω0 (θ0 )k22
≥ τ 2 n−1 .
H(θ, θ0 )
H(θ,θ 0 )≥1
min

7.8

Proofs of Lemma 18 and Lemma 19 in the analysis of Theorem 9

In this section, we prove Lemma 18 and Lemma 19 to establish Theorem 9.
7.8.1

Proof of Lemma 18

It is easy to show diag(Ω−1 ) = I. Now we need to verify Ω ∈ Pα (η 2 , M η). Let S denote
1
(diag(Ω0−1 )) 2 . The Cholesky decomposition of Ω is:
Ω = SΩ0 S
= S(I − A0 )T D0−1 (I − A0 )S
= (S − A0 S)T D0−1 (S − A0 S)
= (I − S −1 A0 S)T SD0−1 S(I − S −1 A0 S).
1

1

According to Lemma 22, η − 2 ≤ λmin (S) ≤ λmax (S) ≤ η 2 and η −1 ≤ λmin (Ω0 ) ≤
λmax (Ω0 ) ≤ η. Therefore, we obtain η −2 ≤ λmin (Ω) ≤ λmax (Ω) ≤ η 2 . Let A ≡
[aij ]p×p = S −1 A0 S. The desired result Ω ∈ Pα (η 2 , M η) then immediately follows from
P
that maxi j<i−k |aij | < M ηk −α for k ∈ [i − 1].
1

Since diag(Ω0−1 )i2 > 0, fi is a strictly increasing function. Therefore, {Ω, {fi }} ∈
Pα0 (η 2 , M η).
7.8.2

Proof of Lemma 19

Set Ik + A∗k (θ)T A∗k (θ) as W (θ). One can check

1
W (θ)
−A∗k (θ)T diag(W (θ)) 2
1

Ω(θ) = −diag(W (θ)) 2 A∗k (θ)
diag(W (θ))
0(p−2k)×k
0(p−2k)×k


0k×(p−2k)

0k×(p−2k)  .
Ip−2k

Then,
1

1

kΩ(θ0 ) − Ω(θ)k2op ≥ kdiag(W (θ)) 2 A∗k (θ) − diag(W (θ0 )) 2 A∗k (θ0 )k2op
1

1

= H(θ, θ0 )(1 + (τ n− 2 )2 )(τ n− 2 )2 ,
which further implies
1
kΩ(θ) − Ω(θ0 )k22
≥ (τ n− 2 )2 .
H(θ, θ0 )
H(θ,θ0 )≥1

min

52

