Ensemble Estimation of Large Sparse Covariance
Matrix Based on Modified Cholesky Decomposition
Xiaoning Kang† and Xinwei Deng‡∗
†

International Business College, Dongbei University of Finance and Economics, China

arXiv:1801.00380v1 [] 1 Jan 2018

‡

Department of Statistics, Virginia Tech, USA
(January 3, 2018)
Abstract

Estimation of large sparse covariance matrices is of great importance for statistical
analysis, especially in the high dimensional setting. The traditional approach such
as sample covariance matrix could perform poorly due to the high dimensionality. In
this work, we propose a positive-definite estimator for the covariance matrix based
on the modified Cholesky decomposition. The modified Cholesky decomposition relies
on the order of variables, which provides the flexibility to obtain a set of covariance
matrix estimates under different orders of variables. The proposed method considers
an ensemble estimator as the “center” of such a set of covariance matrix estimates with
respect to the Frobenius norm. The proposed estimator is not only guaranteed to be
positive definite, but also can capture the underlying sparse structure of the covariance
matrix. Under some weak regularity conditions, both algorithmic convergence and
asymptotical convergence are established. The merits of the proposed method are
illustrated through simulation studies and one real data example.

1

Introduction

High-dimensional data are widely occurred in modern scientific applications, such as bioinformatics, imaging recognition, weather forecasting, and financial service. Estimation of
large covariance matrix from the high-dimensional data is thus an important and challenge
problem in the multivariate data analysis. For example, dimension reduction using principal
∗

Address for correspondence: Xinwei Deng, Associate Professor, Department of Statistics, Virginia Tech,

Blacksburg, VA 24060 (E-mail: xdeng@vt.edu).

1

component analysis usually relies on accurate estimation of covariance matrix. Under the
context of graphical models, the estimation of covariance matrix or its inverse is often used
to infer the network structure of the graph. However, conventional estimation of covariance
matrix is known to perform poorly due to the high dimension problems when the number of
variables is close to or larger than the sample size (Johnstone, 2001). To overcome the curse
of dimensionality, a variety of methods proposed in literature often assumes certain patterns
of sparsity for the covariance matrices.
In this work, our focus is on the estimation of sparse covariance matrix for high-dimensional
data. Early work on covariance matrix estimation includes shrinking eigenvalues of the sample covariance matrix (Dey and Srinivasan, 1985; Haff, 1991), a linear combination of the
sample covariance and a proper diagonal matrix (Ledoit and Wolf, 2004), improving the
estimation based on matrix condition number (Aubry et al., 2012; Won et al., 2013), and
regularizing the eigenvectors of the matrix logarithm of the covariance matrix (Deng and
Tsui, 2013). However, the above mentioned methods do not explore the sparse structure
of the covariance matrix. A sparse covariance matrix estimate can be useful for subsequent
inference, such as inferring the correlation pattern among the variables. Bickel and Levina
(2008) proposed to threshold the small entries of the sample covariance matrix to zeroes and
studied its theoretical behavior when the number of variables is large. Rothman, Levina and
Zhu (2009) considered to threshold the sample covariance matrix with more general thresholding functions. Cai and Yuan (2012) proposed a covariance matrix estimation through
block thresholding. Their estimator is constructed by dividing the sample covariance matrix into blocks and then simultaneously estimating the entries in a block by thresholding.
However, the threshold-based estimator is not guaranteed to be positive definite. To make
the estimate being sparse and positive-definite, Bien and Tibshirani (2011) considered a
penalized likelihood method with a Lasso penalty (Tibshirani, 1996) on the entries of the
covariance matrix. Their idea is similar to the Graphical Lasso for inverse covariance matrix estimation in the literature (Yuan and Lin, 2007; Friedman, Hastie, and Tibshirani,
2008; Rocha, Zhao, and Yu, 2008; Rothman et al., 2008; Yuan, 2008; Deng and Yuan, 2009;
and Yuan, 2010), but the computation is much more complicated due to the non-convexity

2

of the objective function. Xue, Ma and Zou (2012) developed a sparse covariance matrix
estimator for high-dimensional data based on a convex objective function with positive definite constraint and L1 penalty. They also derived a fast algorithm to solve the constraint
optimization problem. Some other work on the estimation of high-dimensional covariance
matrix can be found in Fan, Liao and Mincheva (2013), Liu, Wang and Zhao (2014), Xiao
et al. (2016), Cai, Ren and Zhou (2016). A comprehensive review of recent development
of covariance matrix estimation can be found in Pourahamdi (2013) and Fan, Liao and Liu
(2016).
Another direction of sparse covariance matrix estimation is to take advantage of matrix
decomposition. One popular and effective decomposition is the modified Cholesky decomposition (MCD) (Pourahmadi, 1999; Wu and Pourahmadi, 2003; Pourahmadi, Daniels and
Park, 2007; Rothman, Levina and Zhu, 2009; Dellaportas and Pourahmadi, 2012; Xue, Ma
and Zou, 2012; Rajaratnam and Salzman, 2013). It assumes that the variables have a natural order based on which the variables can be sequentially orthogonalized to re-parameterize
the covariance matrix. By imposing certain sparse structures on the Cholesky factor, it
can result in certain sparse structure on the estimated covariance matrix. For example,
Wu and Pourahmadi (2003) proposed a banded estimator of Cholesky factor, which can be
obtained by regressing each variable only on its closest k predecessors. Bickel and Levina
(2008) showed that banding the Cholesky factor can produce a consistent estimator in the
operator norm under weak conditions on the covariance matrix. Huang et al. (2006) considered to impose an L1 (Lasso) penalty on the entries of the Cholesky factor for estimating
the sparse covariance matrix. However, the MCD-based approach for estimating covariance
matrix depends on the order of variables X1 , . . . , Xp . Such an pre-specification on the order
of variables may not hold in practice. A natural order of variables is often not available or
cannot be pre-determined before in many applications such as the gene expression data and
stock marketing data.
In this paper, we adopt the MCD approach for estimating the large covariance matrix,
but alleviate the drawback of order dependency of the MCD method. Different from most
existing MCD-based approaches, the proposed method takes advantage of the fact that one

3

can always obtain a MCD-based covariance matrix estimate for a given order of variables.
By considering a set of covariance matrix estimates under different orders of variables, we
can obtain an ensemble estimator for the covariance matrix with positive-definite and sparse
properties. Specifically, using the permutation idea, we first obtain a number of estimates of
covariance matrix based on different orders of variables used in the MCD approach. With
such a set of estimates, the proposed ensemble estimator is obtained as the “center” of them
under the Frobenius norm through an L1 penalized objective function. The L1 regularization
is imposed to achieve the sparsity of the estimate. Such an estimator takes advantage
of multiple orders of variables due to the ensemble effort. An efficient algorithm is also
developed to make the computation attractive for solving the estimator. The numerical
convergence of this algorithm is guaranteed with theoretical justification. Furthermore, the
consistent properties of the proposed estimator are also established under Frobenius norm
with some regularity conditions. The simulation and real-data analysis elaborate that the
proposed method performs much better than several existing approaches in terms of both
estimation accuracy and sparsity accuracy.
The remainder of this work is organized as follows. Section 2 briefly reviews the MCD
approach to estimate the covariance matrix. Section 3 introduces the proposed method by
addressing the order issue. An efficient algorithm is also developed to solve the objective
function. In Section 4, the theoretical properties are presented. The simulation study and
one real data example are reported in Section 5 and 6, respectively. We conclude the paper
in Section 7.

2

The Modified Cholesky Decomposition

Without loss of generality, suppose that X = (X1 , . . . , Xp )0 is a p-dimensional random
vector with mean 0 and covariance matrix Σ. Let x1 , . . . , xn be n independent and identically distributed observations following N (0, Σ). Pourahmadi (1999) proposed the modified
Cholesky decomposition (MCD) for the estimation of a covariance matrix, which is statistically meaningful and grantees the positive definiteness of the estimate. This decomposition

4

arises from regressing each variable Xj on its predecessors X1 , . . . , Xj−1 for 2 ≤ j ≤ p.
Specifically, consider to fit a series of regressions
j−1
X
Xj =
(−tjk )Xk + j = X̂j + j ,
k=1

where j is the error term for the jth regression with Ej = 0 and V ar(j ) = d2j . Let 1 = X1
and D = diag(d21 , . . . , d2p ) be the diagonal covariance matrix of  = (1 , . . . , p )0 . Construct
the unit lower triangular matrix T = (tjk )p×p with ones on its diagonal and regression
coefficients (tj1 , . . . , tj,j−1 )0 as its jth row. Then one can have
D = V ar() = V ar(X − X̂) = V ar(T X) = T ΣT 0 ,
and thus
−1

Σ = T −1 DT 0 .

(2.1)

The MCD approach reduces the challenge of modeling a covariance matrix into the task of
modeling (p−1) regression problems, and is applicable in high dimensions. However, directly
imposing the sparse structure on Cholesky factor matrix T in (2.1) does not imply the sparse
pattern of covariance matrix Σ since it requires an inverse of T . Thus the formulation (2.1)
is not convenient to impose a sparse structure on the estimation of Σ. Alternatively, one
can consider a latent variable regression model based on the MCD. Writing X = L would
lead to
V ar(X) = V ar(L)
Σ = LDL0 .

(2.2)

This decomposition can be interpreted as resulting from a new sequence of regressions,
where each variable Xj is regressed on all the previous latent variable 1 , . . . , j−1 rather
than themselves. It gives a sequence of regressions
Xj = lTj  =

X

ljk k + j ,

j = 2, . . . , p,

k<j

where lj = (ljk ) is the jth row of L. Here ljj = 1 and ljk = 0 for k > j.
5

(2.3)

With the data matrix X = (x1 , . . . , xn )0 , define its jth column to be x(j) . Let e(j) denote
the residuals of x(j) for j ≥ 2, and e(1) = x(1) . Let Z(j) = (e(1) , . . . , e(j−1) ) be the matrix
containing the first (j − 1) residuals. Now we consider a sparse covariance matrix estimate
Σ̂ by encouraging the sparsity on L̂. One approach to achieving such sparsity is to use the
Lasso for the regression (Tibshirani, 1996)
l̂j = arg min kx(j) − Z(j) lj k22 + ηj klj k1 , j = 2, . . . , p,
lj

(2.4)

where ηj ≥ 0 is a tuning parameter and selected by cross validation. k · k1 stands for the
vector L1 norm. e(j) = x(j) − Z(j) lj is used to construct the residuals for the last column of
Z(j+1) . Then d2j is estimated as the sample variance of e(j)
dˆ2j = Vd
ar(ê(j) ) = Vd
ar(x(j) − Z(j) l̂j )

(2.5)

ˆ0 will be a sparse covariwhen constructing matrix D̂ = diag(dˆ21 , . . . , dˆ2p ). Hence, Σ̂ = L̂D̂ L
ance matrix estimate.

3

The Proposed Method

ˆ0 depends on the order of random variables X , . . . , X .
Clearly, the estimate Σ̂ = L̂D̂ L
1
p
It means that different orders would lead to different sparse estimates of Σ. To address
this order-dependent issue, we consider an ensemble estimation of Σ by using the idea of
permutation. Specifically, we generate M different permutations of {1, . . . , p} as the orders
of the random variables, denoted by πk , s, k = 1, . . . , M . Let P πk be the corresponding
permutation matrix. Under a variable order πk , the estimate is obtained as
0

Σ̂πk = L̂πk D̂ πk L̂πk ,

(3.1)

where L̂πk and D̂ πk are calculated based on (2.4) and (2.5). Then transforming back to the
original order, we have
Σ̂k = P πk Σ̂πk P 0πk .

6

(3.2)

To obtain an ensemble estimator for Σ, a naive solution would be Σ̄ =

1
M

PM

k=1

Σ̂k . However,

such an estimate may not be sparse since the sparse structure in Σ̂k is destroyed by the
average.
In order to simultaneously achieve the positive definiteness and sparsity for the estimator,
we propose to consider the estimate
M
1 X
Σ̂ = arg min
kΣ − Σ̂k k2F + λ|Σ|1 ,
Σν I 2M k=1

(3.3)

where k · kF stands for the Frobenius norm, λ ≥ 0 is a tuning parameter, and | · |1 is L1
norm for all the off-diagonal elements. Here ν is some positive arbitrarily small number.
The constraint Σ  νI is to guarantee the positive-definiteness of the estimate. The penalty
term is to encourage the sparse pattern in Σ̂. It is worth pointing out that, if λ = 0 in (3.3),
PM
the solution of Σ̂ would be Σ̄ = M1
k=1 Σ̂k ; if without the constraint Σ  νI in (3.3), the
solution of Σ̂ would be the soft-threshold estimate of Σ̄. Therefore, the proposed estimate is
to pursue the “center” of the multiple estimates Σ̂k , while maintain the properties of being
positive-definite and sparse.
To efficiently solve the optimization in (3.3), we employ the alternating direction method
of multipliers (ADMM) (Boyd et al., 2011). The ADMM technique has been widely used in
solving the convex optimization under the content of L1 penalized covariance matrix estimation (Xue, Ma and Zou, 2012). The ADMM technique does not require the differentiability
assumption of the objective function and it is easy to implement. Specifically, let us first
introduce a new variable Φ and an equality constraint as follows
M
1 X
kΣ − Σ̂k k2F + λ|Σ|1 : Σ = Φ, Φ  νI}.
(Σ̂, Φ̂) = arg min {
2M
Σ ,Φ
k=1

(3.4)

Note that the solution of (3.4) gives solution to (3.3). To solve (3.4), we minimize its
augmented Lagrangian function for some given penalty parameter τ as
L(Σ, Φ; Λ) =

M
1 X
1
kΣ − Σ̂k k2F + λ|Σ|1 − hΛ, Φ − Σi + kΦ − Σk2F ,
2M k=1
2τ

(3.5)

where Λ is the Lagrangian multiplier. The notation h·, ·i stands for the matrix inner product
P
as hA, Bi = i,j aij bij , where aij and bij are the elements of matrices A and B. The ADMM
7

iteratively solves the following steps sequentially for i = 0, 1, 2, . . . till convergence
Φ step : Φi+1 = arg min L(Σi , Φ; Λi )
Φν I

(3.6)

Σ step : Σi+1 = arg min L(Σ, Φi+1 ; Λi )
Σ

(3.7)

1
(3.8)
Λ step : Λi+1 = Λi − (Φi+1 − Σi+1 ).
τ
Pp
0
Assume the eigenvalue decomposition of a matrix Z is
i=1 λi ξ ξ, and define (Z)+ =
Pp
0
i=1 max(λi , ν)ξ i ξ i . Then we develop the closed form for (3.6) as
∂L(Σi , Φ; Λi )
1
= −Λi + (Φ − Σi ) , 0
∂Φ
τ
Φ = Σi + τ Λi
Φi+1 = (Σi + τ Λi )+ .
Next, define an element-wise soft threshold for each entry zij in matrix Z as s(Z, δ) =
{s(zij , δ)}1≤i,j≤p with
s(zij , δ) = sign(zij )max(|zij | − δ, 0)I{i6=j} + zij I{i=j} .
Then the solution of (3.7) is derived as
M
∂L(Σ, Φi+1 ; Λi )
1 X
1
=
(Σ − Σ̂k ) + Λi + (Σ − Φi+1 ) + λsign∗ (Σ) , 0
∂Σ
M k=1
τ
M
1 X
(τ + 1)Σ = τ (
Σ̂k − Λi ) + Φi+1 − λτ sign∗ (Σ)
M k=1

Σ

i+1

M
1 X
= {s(τ (
Σ̂k − Λi ) + Φi+1 , λτ )}/(τ + 1),
M k=1

where sign∗ (Σ) means sign(Σ) with the diagonal elements replaced by 0 vector. Algorithm
1 summarizes the developed procedure for solving (3.3) by using the ADMM technique.
Algorithm 1.
Step 1: Input initial values Σinit , Λinit and τ .
8

Step 2: Φi+1 = (Σi + τ Λi )+ .
PM
i
i+1
, λτ )}/(τ + 1).
Step 3: Σi+1 = {s(τ ( M1
k=1 Σ̂k − Λ ) + Φ
Step 4: Λi+1 = Λi − τ1 (Φi+1 − Σi+1 ).
Step 5: Repeat Step 2 - 4 till convergence.
This algorithm converges fast and produces the optimal solution of arg min L(Σ, Φ; Λ) in
PM
(3.5). In practice, the initial estimate Σinit is set to be the naive estimate Σ̄ = M1
k=1 Σ̂k .
The Λinit is set to be zero matrix, and τ = 2 and ν = 10−4 (Xue, Ma and Zou, 2012).
The optimal value of tuning parameter λ in (3.5) is chosen based on Bayesian information
criterion (BIC) (Yuan and Lin, 2007)
−1

−1

BIC(λ) = − log |Σ̂λ | + tr[Σ̂λ S] +

log n X
êij (λ),
n i≤j

(λ)

where S is the sample covariance matrix, Σ̂λ = (σ̂ij )p×p indicates the estimate of Σ obtained
(λ)

by applying our algorithm with tuning parameter λ. êij (λ) = 0 if σ̂ij = 0, and êij (λ) = 1
otherwise.
Note that the implementation of the proposed method also requires the choice of M , the
number of permutation orders. Obviously, the number of all possible permeation orders is p!,
which increases rapidly as the number of variables p increases. To get an appropriate value
for M for efficient computation, we have tried M = 10, 30, 50, 100 and 150 as the potential
number of random orders. The performances were quite comparable for M larger than 30.
Hence, we choose M = 100 as the number of permutation orders for the proposed method
in this paper. The more detailed discussion and justification on the choice of M
can be found in Kang and Deng (2017).

4

Theoretical Properties

In this section, Theorem 1 states that the sequence (Σi , Φi , Λi ) generated by Algorithm
+

+

+

1 from any starting point numerically converges to an optimal minimizer (Σ̂ , Φ̂ , Λ̂ ) of
+

(3.5), where Λ̂ is the optimal dual variable. Theorem 2 provides the asymptotic convergence
rate of the estimator Σ̂k in (3.2) under the Frobenius norm. Theorem 3 and 4 demonstrate
9

the asymptotic properties of the proposed estimator under some weak regularity conditions.
The details of the proofs of Theorems 1- 4 are in the Appendix. To facilitate the presentation
of the proofs, we first introduce some notations. Define a 2p by 2p matrix J as


τ I p×p
0
.
J =
−1
0
τ I p×p
Let the notation k · k2J be kU k2J = hU , J U i and hU , V iJ = hU , J V i. Let Σ0 = (σij0 )p×p =
L0 D 0 L00 be the true covariance matrix for the observations X = (xij )n×p , and define the
number of nonzero off-diagonal elements of Σ0 as s0 . Denote the maximal true variance in
(π )

Σ0 by σmax . Let Zπk = {(j, k) : k < j, l0jkk 6= 0} be the collection of nonzero elements in the
lower triangular part of matrix L0πk . Denote by s1 the maximum of the cardinality of Zπk
for k = 1, 2, . . . , M . Now we present the following lemma and theory.
+

+

+

Lemma 1. Assume that (Σ̂ , Φ̂ ) is an optimal solution of (3.4) and Λ̂ is the corresponding optimal dual variable with the equality constraint Σ = Φ, then the sequence (Σi , Φi , Λi )
generated by Algorithm 1 satisfies
kW + − W i k2J − kW + − W i+1 k2J ≥ kW i − W i+1 k2J ,
+

+

where W + = (Λ̂ , Σ̂ )0 and W i = (Λi , Σi )0 .
Theorem 1. Suppose x1 , . . . , xn are n independent and identically distributed observations
from Np (0, Σ). Then the sequence (Σi , Φi , Λi ) generated by Algorithm 1 from any starting
point converges to an optimal minimizer of the objective function in (3.5).
Theorem 1 demonstrates the convergence of Algorithm 1. It automatically indicates that
the sequence Σi , i = 1, 2, . . ., produced by Algorithm 1 converges to an optimal solution of
the objective (3.3).
In order to achieve the consistent property of Σ̂k in (3.2) obtained under each order using
the MCD approach, one needs a basic assumption that there exists a constant θ > 1 such
that the singular values of the true covariance matrix are bounded as
1/θ < svp (Σ0 ) ≤ sv1 (Σ0 ) < θ,
10

(4.1)

where we use sv1 (A), sv2 (A), . . . , svp (A) to indicate the singular values of matrix A in a
decreasing order. They are the square root of the eigenvalues of matrix AA0 . This assumption is also made in Rothman (2008), Lam and Fan (2009) and Guo et al. (2011). The
assumption guarantees the positive definiteness property and makes inverting the covariance
matrix meaningful. The following lemma and theorem give the convergence rate of Σ̂k in
(3.2) under the Frobenius norm.
0

Lemma 2. Let Σ0 = L0 D 0 L0 be the MCD of the true covariance matrix. Since the singular
values of Σ0 are bounded, there exist constants θ1 and θ2 such that 0 < θ1 < svp (Σ0 ) ≤
sv1 (Σ0 ) < θ2 < ∞, then there exist constants h1 and h2 such that
h1 < svp (L0 ) ≤ sv1 (L0 ) < h2 ,
and
h1 < svp (D 0 ) ≤ sv1 (D 0 ) < h2 .
0

Lemma 3. Let Σ0πk = L0πk D 0πk L0πk be the MCD of the true covariance matrix regarding a
P
variable order πk . Under (4.1), assume the tuning parameters ηj in (2.4) satisfy pj=1 ηj =
O(log(p)/n); s1 and p satisfy (s1 + p) log(p)/n = o(1), then L̂πk and D̂ πk in (3.1) have the
following consistent properties
p
kL̂πk − L0πk kF = Op ( s1 log(p)/n),
p
kD̂ πk − D 0πk kF = Op ( p log(p)/n).
Theorem 2. Assume that all the conditions in Lemma 3 hold. Then the estimator Σ̂k in
(3.2) has the following consistent property
kΣ̂k − Σ0 kF = Op (

p
(s1 + p) log(p)/n).

Theorem 2 establishes the consistent property of the estimate Σ̂k under variable order
πk , k = 1, 2, . . . , M . From this result, the convergence rate of the proposed estimator for
estimating a sparse covariance matrix can be obtained based on Frobenius norm by Theorem
3 and 4 as follows.
11

Theorem 3. Assume all the conditions in Lemma 3 hold and log p ≤ n. Under the
exponential-tail condition that for all |t| ≤ ρ and 1 ≤ i ≤ n, 1 ≤ j ≤ p
E{exp(tx2ij )} ≤ K1 .
For any m > 0, set
λ=

log
c20

p

n


+ c1

log p
n

1/2
,

where
1
c0 = eK1 ρ1/2 + ρ−1/2 (m + 1)
2
and
1 2
1
c1 = 2K1 (ρ−1 + ρσmax
) exp( ρσmax ) + 2ρ−1 (m + 2).
4
2
With probability at least 1 − 3p−m , we have
+

||Σ̂ − Σ0 ||F ≤ 5λ(s0 + p)1/2 .
Theorem 4. Assume all the conditions in Lemma 3 hold and p ≤ cnγ for some c > 0.
Under the polynomial-tail condition that for all γ > 0, ε > 0 and 1 ≤ i ≤ n, 1 ≤ j ≤ p
E{|xij |4(1+γ+ε) } ≤ K2 .
For any m > 0, set
log p
λ = 8(K2 + 1)(m + 1)
+ 8(K2 + 1)(m + 2)
n



log p
n

1/2
.

With probability at least 1 − O(p−m ) − 3K2 p(log p)2(1+γ+ε) n−γ−ε , we have
+

||Σ̂ − Σ0 ||F ≤ 5λ(s0 + p)1/2 .

5

Simulation Study

In this section, we conduct a comprehensive simulation study to evaluate the performance
of the proposed method. Suppose that data x1 , . . . , xn are generated independently from

12

the normal distribution N (0, Σ). Here we consider the following four covariance matrix
structures.
Model 1. Σ1 = MA(0.5, 0.3). The diagonal elements are 1 with the first sub-diagonal
elements 0.5 and the seconde sub-diagonal elements 0.3.
Model 2. Σ2 = AR(0.5). The conditional covariance between any two random variables
Xi and Xj is fixed to be 0.5|i−j| , 1 ≤ i, j ≤ p.
Model 3. Σ3 is generated by randomly permuting rows and corresponding columns of
Σ1 .
Model 4. Σ4 is generated by randomly permuting rows and corresponding columns of
Σ2 .
Note that Models 1-2 consider the banded or nearly-banded structure for the covariance
matrix. While the sparse pattern of the covariance matrix in Models 3-4 is not structured
due to the random permutation. For each case, we generate the data set with three settings
of different sample sizes and variable sizes: (1) n = 50, p = 30; (2) n = 50, p = 50 and (3)
n = 50, p = 100.
The performance of the proposed estimator is examined in comparison with several other
approaches, which are divided into three classes. The first class is the sample covariance
matrix S that serves as the benchmark. The second class is composed of three methods that
deal with the variable order used in the MCD, including the MCD-based method with BIC
order selection (BIC) (Dellaportas and Pourahmadi, 2012), the best permutation algorithm
(BPA) (Rajaratnam and Salzman, 2013) and the proposed method (Proposed). The BPA
selects the order of variables used in the MCD approach (2.2) such that ||D||2F is minimized.
While the BIC method determines the order of variables in the MCD approach (2.3) in
a forward selection fashion. That is, in each step, it selects a new variable having the
smallest value of BIC when regressing it on its previous residuals. For example, suppose that
C = {Xi1 , . . . , Xik } is the candidate set of variables and there are (p − k) variables already
chosen and ordered. By regressing each Xj , j = i1 , . . . , ik on the residuals [e1 , . . . , ep−k ], we
can assign the variable corresponding to the minimum BIC value among the k regressions
to the (p − k + 1)th position of the order.

13

The third class of competing methods consists of four approaches that estimate the
sparse covariance matrix directly without considering the variable order, including Bien and
Tibshirani’s estimate (BT) (Bien and Tibshirani, 2011), Bickel and Levina’s estimate (BL)
(Bickel and Levina, 2008), Xue, Ma and Zou’s estimate (XMZ) (Xue, Ma and Zou, 2012)
and Rothman et al.’s estimate (RLZ) (Rothman et al., 2010). The BT estimate minimizes
the negative log-likelihood function with L1 penalty on the entries of the covariance matrix,
that is,
n

−1

−1

o

Σ̂BT = arg min − log |Σ | + tr(Σ S) + η|Σ|1 ,
Σ0
where η ≥ 0 is the tuning parameter. The optimization is solved by using the majorizationminimization algorithm (Lange, Hunter and Yang, 2000). The BL estimate is obtained by
imposing hard thresholding (Bickel and Levina, 2008) on the entries of the sample covariance
matrix, so the resultant estimate may not be positive definite. The XMZ estimate is obtained
from encouraging the sparsity on the sample covariance matrix S as well as maintaining the
property of positive definiteness. The resultant estimate is
1
Σ̂XM Z = arg min kΣ − Sk2F + η|Σ|1 .
Σν I 2
The RLZ estimate is to introduce sparsity in the Cholesky factor matrix L by estimating
the first k sub-diagonals of L and setting the rest to zeroes. It means that each variable is
only regressed on the k previous residuals in (2.3). The tuning parameter k can be chosen
by AIC or cross validation.
To measure the accuracy of covariance matrix estimates Σ̂ = (σ̂ij )p×p obtained from
each approach, we consider the entropy loss (EN), quadratic loss (QL), L1 norm and mean
absolute error (MAE), defined as follows:
EN = tr[Σ−1 Σ̂] − log |Σ−1 Σ̂| − p,
−1
1
tr[Σ̂ Σ − I]2 ,
p
X
L1 norm = max
|σ̂ij − σij |,

QL =

j

i
p

MAE =

p

1 XX
|σ̂ij − σij |.
p i=1 j=1
14

Here we have not included the Kullback-Leibler loss (Kullback and Leibler, 1951), since it
is more suitable in measuring inverse covariance matrix estimates rather than covariance
matrix estimates (Levina et al., 2008). In addition, to gauge the performance of capturing
sparse structure, we consider the false selection loss (FSL), which is the summation of false
positive (FP) and false negative (FN). Here we say a FP occurs if a nonzero element in the
true matrix is incorrectly estimated as a zero. Similarly, a FN occurs if a zero element in
the true matrix is incorrectly identified as a nonzero. The FSL is computed in percentage as
(FP + FN) / p2 . For each loss function above, Table 1 to Table 4 report the averages of the
performance measures and their corresponding standard errors in the parentheses over 100
replicates. For each model, the two methods with lowest averages regarding each measure are
shown in bold. Dashed lines in the tables represent the corresponding values not available
due to matrix singularity.
For a short summary of Tables 1-4, the numerical results show that the proposed method
generally provides a superior performance over other approaches in comparison. It is able
to accurately catch the underlying sparse structure of the covariance matrix. When the
underlying covariance matrix is banded or tapered, the proposed method is comparable to
the RLZ method, and performs better than other approaches. Note that the RLZ method
targets on the banded covariance matrix. When the underlying structure of covariance matrix
is more general without any specification, the proposed method still performs well. As in
the high-dimensional cases, the advantage of the proposed method is even more evident.
Specifically, Table 1 and 2 summarize the comparison results for Models 1 and 2, respectively. From the perspective of competing methods, the sample covariance matrix S,
serving as a benchmark approach, does not give the sparse structure and performs poorly
under all the loss measures. The BIC and BPA in the second class of approaches provide
sparse covariance matrix estimates compared with S, but their false selection loss (FSL)
are considerably larger than the proposed method. Moreover, the proposed method greatly
outperforms the BIC and BPA regarding QL, L1 and MAE for all settings of p = 30, 50 and
100. Although the proposed method is comparable to the BIC and BPA methods under EN
criterion when p = 30, it performs slightly better when p = 50 and much better in the case

15

Table 1: The averages and standard errors (in parenthesis) of estimates for Model 1.
EN

QL

L1

MAE

FSL (%)

S

12.42 (0.08)

88.11 (2.83)

5.21 (0.07)

3.51 (0.03)

83.96 (0.01)

BIC

7.22 (0.08)

8.59 (0.59)

2.96 (0.04)

1.76 (0.02)

52.32 (0.55)

BPA

6.02 (0.09)

4.38 (0.23)

2.74 (0.05)

1.55 (0.02)

46.53 (0.68)

BT

7.78 (0.03)

19.44 (0.29)

2.14 (0.01)

1.85 (0.00)

6.92 (0.10)

BL

-

-

2.33 (0.05)

1.17 (0.02)

6.81 (0.14)

XMZ

10.61 (0.17)

0.13 (0.02)

1.88 (0.01)

1.25 (0.01)

7.13 (0.18)

RLZ

9.36 (0.07)

0.35 (0.03)

1.55 (0.01)

1.04 (0.01)

6.22 (0.00)

Proposed

7.10 (0.11)

0.38 (0.03)

1.92 (0.02)

1.22 (0.01)

6.75 (0.15)

S

-

-

8.26 (0.07)

5.75 (0.03)

90.19 (0.01)

BIC

15.77 (0.21)

125.73 (28.51)

3.85 (0.06)

1.98 (0.01)

43.30 (0.47)

BPA

12.98 (0.16)

15.67 (1.46)

3.62 (0.08)

1.84 (0.02)

41.49 (0.63)

BT

15.07 (0.28)

106.19 (8.57)

2.42 (0.02)

1.93 (0.02)

10.68 (0.46)

BL

-

-

2.41 (0.05)

1.27 (0.01)

4.80 (0.06)

XMZ

20.45 (0.29)

0.14 (0.02)

1.98 (0.01)

1.36 (0.01)

5.13 (0.07)

RLZ

16.20 (0.07)

0.62 (0.05)

1.63 (0.01)

1.06 (0.00)

3.84 (0.00)

Proposed

13.68 (0.10)

1.10 (0.06)

2.02 (0.01)

1.35 (0.01)

4.36 (0.06)

S

-

-

16.15 (0.12)

11.43 (0.03)

95.01 (0.00)

BIC

42.56 (0.45)

177503 (149636)

5.26 (0.07)

2.24 (0.01)

32.75 (0.36)

BPA

35.68 (0.38)

3223.74 (1749.81)

5.29 (0.11)

2.20 (0.02)

33.60 (0.38)

BT

28.78 (0.33)

97.08 (3.21)

2.40 (0.03)

1.87 (0.02)

4.08 (0.33)

BL

-

-

2.67 (0.06)

1.42 (0.01)

2.83 (0.02)

XMZ

364.15 (0.18)

Inf (Inf)

16.14 (0.12)

11.42 (0.03)

94.96 (0.01)

RLZ

33.40 (0.12)

1.28 (0.06)

1.69 (0.01)

1.08 (0.00)

1.96 (0.00)

Proposed

31.28 (0.14)

4.04 (0.12)

2.12 (0.01)

1.49 (0.00)

2.38 (0.02)

p = 30

p = 50

p = 100

16

Table 2: The averages and standard errors (in parenthesis) of estimates for Model 2.
EN

QL

L1

MAE

FSL (%)

S

12.58 (0.08)

92.44 (2.69)

5.10 (0.07)

3.49 (0.03)

46.66 (0.01)

BIC

5.35 (0.09)

5.55 (0.40)

2.94 (0.04)

1.91 (0.01)

43.47 (0.30)

BPA

4.60 (0.07)

3.64 (0.22)

2.82 (0.04)

1.78 (0.01)

42.53 (0.31)

BT

5.40 (0.03)

19.36 (0.32)

2.52 (0.01)

2.19 (0.00)

43.82 (0.08)

BL

-

-

2.61 (0.04)

1.57 (0.01)

43.70 (0.20)

XMZ

5.82 (0.11)

0.12 (0.01)

2.27 (0.01)

1.64 (0.01)

42.24 (0.20)

RLZ

3.16 (0.03)

0.43 (0.04)

1.89 (0.01)

1.34 (0.01)

43.56 (0.00)

Proposed

4.09 (0.06)

0.48 (0.04)

2.30 (0.01)

1.66 (0.01)

41.82 (0.16)

S

-

-

8.54 (0.08)

5.84 (0.03)

65.59 (0.01)

BIC

11.08 (0.20)

42.50 (7.50)

3.93 (0.08)

2.18 (0.01)

42.28 (0.24)

BPA

9.17 (0.15)

11.93 (0.86)

3.70 (0.06)

2.08 (0.02)

41.73 (0.28)

BT

10.07 (0.18)

52.46 (3.38)

2.69 (0.02)

2.26 (0.01)

30.34 (0.13)

BL

-

-

2.91 (0.05)

1.68 (0.01)

29.45 (0.07)

XMZ

11.26 (0.20)

0.12 (0.01)

2.36 (0.01)

1.76 (0.01)

28.85 (0.07)

RLZ

5.48 (0.04)

0.64 (0.04)

1.96 (0.01)

1.38 (0.00)

28.48 (0.00)

Proposed

7.22 (0.06)

0.98 (0.05)

2.40 (0.01)

1.77 (0.01)

28.60 (0.07)

S

-

-

16.04 (0.12)

11.43 (0.04)

81.86 (0.00)

BIC

29.70 (0.55)

17694 (11352)

5.32 (0.08)

2.47 (0.01)

33.77 (0.25)

BPA

23.65 (0.36)

854.23 (232.83)

5.16 (0.11)

2.44 (0.02)

34.41 (0.31)

BT

21.09 (0.36)

95.23 (3.68)

2.80 (0.03)

2.24 (0.02)

17.01 (0.23)

BL

-

-

3.04 (0.05)

1.82 (0.01)

16.07 (0.02)

XMZ

369.27 (0.19)

Inf (Inf)

16.03 (0.12)

11.42 (0.04)

81.83 (0.00)

RLZ

11.03 (0.07)

1.33 (0.06)

2.05 (0.02)

1.41 (0.00)

15.12 (0.00)

Proposed

16.29 (0.09)

3.46 (0.12)

2.49 (0.01)

1.90 (0.00)

15.63 (0.02)

p = 30

p = 50

p = 100

17

Table 3: The averages and standard errors (in parenthesis) of estimates for Model 3.
EN

QL

L1

MAE

FSL (%)

S

12.46 (0.09)

91.29 (3.34)

5.08 (0.06)

3.48 (0.02)

83.96 (0.01)

BIC

7.28 (0.10)

9.37 (0.70)

2.93 (0.04)

1.76 (0.01)

52.71 (0.51)

BPA

5.87 (0.09)

4.54 (0.38)

2.60 (0.05)

1.52 (0.02)

46.59 (0.78)

BT

7.77 (0.04)

19.38 (0.34)

2.14 (0.01)

1.85 (0.00)

6.94 (0.10)

BL

-

-

2.24 (0.05)

1.14 (0.01)

6.62 (0.15)

XMZ

10.53 (0.14)

0.16 (0.02)

1.87 (0.01)

1.23 (0.01)

7.00 (0.15)

RLZ

16.75 (0.11)

0.37 (0.03)

2.28 (0.02)

1.68 (0.00)

15.55 (0.00)

Proposed

6.96 (0.10)

0.38 (0.03)

1.89 (0.02)

1.21 (0.01)

6.84 (0.13)

S

-

-

8.28 (0.08)

5.75 (0.03)

90.19 (0.01)

BIC

16.00 (0.21)

113.06 (22.87)

3.85 (0.08)

1.98 (0.02)

43.07 (0.50)

BPA

12.92 (0.16)

15.30 (1.73)

3.55 (0.08)

1.83 (0.02)

40.90 (0.65)

BT

15.58 (0.26)

114.94 (9.63)

2.45 (0.02)

1.96 (0.01)

11.21 (0.47)

BL

-

-

2.46 (0.05)

1.27 (0.01)

4.79 (0.06)

XMZ

20.40 (0.32)

0.17 (0.02)

1.97 (0.01)

1.36 (0.01)

5.18 (0.07)

RLZ

31.67 (0.17)

0.69 (0.05)

2.43 (0.01)

1.90 (0.00)

11.36 (0.00)

Proposed

13.74 (0.11)

1.13 (0.07)

2.00 (0.01)

1.36 (0.01)

4.39 (0.07)

S

-

-

16.18 (0.10)

11.43 (0.03)

95.01 (0.00)

BIC

42.88 (0.52)

29998 (13821)

5.36 (0.10)

2.26 (0.01)

32.65 (0.37)

BPA

35.37 (0.43)

2794.90 (1780.94)

5.12 (0.12)

2.19 (0.02)

33.30 (0.41)

BT

28.78 (0.37)

99.03 (3.40)

2.36 (0.02)

1.89 (0.02)

3.91 (0.30)

BL

-

-

2.58 (0.05)

1.41 (0.01)

2.83 (0.02)

XMZ

364.17 (0.16)

Inf (Inf)

16.17 (0.10)

11.42 (0.03)

94.96 (0.00)

RLZ

64.68 (0.22)

1.28 (0.06)

2.50 (0.01)

1.92 (0.00)

5.72 (0.00)

Proposed

31.33 (0.15)

3.81 (0.11)

2.10 (0.01)

1.49 (0.00)

2.39 (0.02)

p = 30

p = 50

p = 100

18

Table 4: The averages and standard errors (in parenthesis) of estimates for Model 4.
EN

QL

L1

MAE

FSL (%)

S

12.48 (0.08)

89.99 (2.78)

5.16 (0.07)

3.49 (0.03)

46.66 (0.01)

BIC

5.22 (0.08)

5.30 (0.39)

2.89 (0.04)

1.90 (0.01)

43.03 (0.33)

BPA

4.48 (0.08)

3.55 (0.23)

2.79 (0.04)

1.76 (0.01)

42.66 (0.35)

BT

5.36 (0.04)

19.32 (0.35)

2.50 (0.01)

2.19 (0.00)

43.79 (0.08)

BL

-

-

2.69 (0.06)

1.58 (0.02)

43.86 (0.20)

XMZ

5.86 (0.09)

0.10 (0.01)

2.24 (0.01)

1.64 (0.01)

42.31 (0.19)

RLZ

11.57 (0.10)

0.46 (0.04)

2.67 (0.02)

2.15 (0.00)

50.67 (0.00)

Proposed

4.06 (0.05)

0.48 (0.04)

2.27 (0.01)

1.66 (0.01)

42.06 (0.16)

S

-

-

8.24 (0.08)

5.74 (0.03)

65.58 (0.01)

BIC

11.16 (0.21)

352.49 (280.45)

3.95 (0.08)

2.17 (0.01)

42.10 (0.24)

BPA

9.22 (0.15)

21.30 (8.30)

3.72 (0.08)

2.07 (0.01)

41.28 (0.30)

BT

10.41 (0.21)

60.15 (5.09)

2.71 (0.01)

2.27 (0.01)

30.31 (0.14)

BL

-

-

2.79 (0.06)

1.68 (0.01)

29.50 (0.07)

XMZ

10.96 (0.19)

0.10 (0.01)

2.37 (0.01)

1.75 (0.01)

28.78 (0.08)

RLZ

19.21 (0.11)

0.72 (0.04)

2.77 (0.01)

2.24 (0.00)

33.60 (0.00)

Proposed

7.29 (0.06)

1.13 (0.06)

2.41 (0.01)

1.78 (0.01)

28.56 (0.07)

S

-

-

16.10 (0.13)

11.44 (0.04)

81.85 (0.00)

BIC

30.15 (0.51)

21511.89 (6496.05)

5.39 (0.10)

2.46 (0.01)

33.63 (0.25)

BPA

23.53 (0.35)

7354.89 (5272.40)

5.39 (0.13)

2.43 (0.01)

33.92 (0.29)

BT

20.62 (0.33)

90.39 (3.13)

2.73 (0.02)

2.24 (0.02)

16.43 (0.17)

BL

-

-

3.03 (0.06)

1.82 (0.01)

16.08 (0.02)

XMZ

369.42 (0.22)

Inf (Inf)

16.09 (0.13)

11.43 (0.04)

81.82 (0.01)

RLZ

40.19 (0.21)

1.30 (0.08)

2.89 (0.01)

2.31 (0.00)

18.44 (0.00)

Proposed

16.39 (0.08)

3.41 (0.13)

2.49 (0.01)

1.90 (0.00)

15.64 (0.02)

p = 30

p = 50

p = 100

19

of p = 100.
For the BT method in the third class of approaches, the proposed method significantly
outperforms it in capturing the sparse structure for the cases of p = 50 and p = 100.
Furthermore, the proposed method gives superior performance to the BT with respect to
all the other loss criteria, especially QL measure. In comparison with the BL method, the
proposed method performs similarly. It is known that the BL method is asymptotically
optimal for sparse covariance matrix (Bickel and Levina, 2008). However, the estimated
covariance matrix obtained from the BL method does not guarantee to be positive definite,
which would result in matrix singularity in computing EN and QL losses. Compared with
XMZ approach, the proposed method is superior or comparable with respect to all the loss
measures except QL criterion in the settings of p = 30 and 50. In the high dimensional
case when p = 100, the proposed method performs much better than the XMZ approach.
Finally, it is seen that the performance of the proposed method is comparable to that of
the RLZ approach for Models 1-2 regarding these loss criteria. This is not surprising, since
the covariance matrices of Models 1-2 are banded and tapered respectively, and the RLZ
approach is designated to estimate such covariance matrix structures.
Table 3 and 4 present the comparison results for Models 3 and 4, respectively. Different
from Models 1-2, the covariance matrices under Models 3-4 are unstructured. This implies
that the RLZ approach does not have the advantage. Hence, it is clearly seen that the
proposed method performs much better than the RLZ approach, especially at capturing
the sparse structure and with respect to EN loss. Generally, the proposed method provides
superior performance to other approaches, with similar comparison results as described under
Models 1-2.

6

Application

In this section, a real prostate cancer data set (Glaab et al., 2012) is used to evaluate
the performance of the proposed method in comparison with other approaches described in
Section 5. It contains two classes with 50 normal samples and 52 prostate cancer samples,

20

and 2135 gene expression values recorded for each sample. Data are available online at
http://ico2s.org/datasets/microarray.html. Since it includes a large number of variables,
the variable screening procedure is performed through two sample t-test. Specifically, for
each variable, t-test is conducted against the two classes of the prostate cancer data such
that the variables corresponding to large values of test statistics are ranked as significant
variables. Then the top 50 significant variables as group 1 and the top 50 nonsignificant
variables as group 2 are selected for data analysis. As mentioned in Xue, Ma and Zou
(2012), there is supposed to be some correlations within each group of variables, but no
correlations between group 1 variables and group 2 variables. Data are centered within each
class and then used for the analysis. In the section, to make each variable at the same scale,
we focus on the correlation matrix rather than the covariance matrix.
S

AVE

BIC

BPA

Proposed

BT

BL

XMZ

RLZ

Figure 1: Heatmaps of the absolute values of the correlation matrices obtained from the
proposed method and other approaches for prostate cancer data. Darker colour indicates
higher density; lighter colour indicates lower density.

21

35

35

30

Proposed
BIC
BT
RLZ

25
20
0

0

5

5

10

15

Eigenvalues

20
15
10

Eigenvalues

25

30

Sample
AVE
BPA
XMZ
BL

0

20

40

60

80

100

0

index

20

40

60

80

100

index

Figure 2: Scree plot for correlation matrices obtained from the proposed method and other
approaches for prostate cancer data.

Figure 1 shows the heatmaps of the absolute values of the estimated correlation matrices
obtained from each method. It is clear to see that, overall, the estimated sparse pattern
of the proposed method, BT, and XMZ well matches the expected sparse pattern. They
shrink the nonzero correlations whereas other methods (except RLZ) do not. The proposed
method and XMZ perform the best in capturing the sparse structure with two diagonal
blocks, followed by BT producing a little denser structure. All the rest approaches either
result in a much sparser matrix as diagonal matrix (i.e. RLZ) or fail to identify the sparsity
pattern (i.e. S, BIC, BPA and BL).
In addition, Figure 2 displays the eigenvalues in the decreasing order for each estimate.
We see that the sample covariance matrix has the most spread out eigenvalues, followed by
BL and BPA. The eigenvalues from the proposed method and RLZ have the least spread.
The largest and smallest eigenvalues obtained from each approach are summarized in Table

22

Table 5: The largest and smallest eigenvalues of correlation matrices obtained from each
approach for prostate data.
Method

S

BIC

BPA

BT

BL

XMZ

RLZ

Proposed

Largest

33.7

11.1

25.0

12.6

33.6

14.8

5.1

5.5

3.3e-03

6.8e-03 1.6e-02 -2.9e-01 4.9e-04 2.9e-02

Smallest 1.2e-05

3.0e-01

5. BT estimator yields a negative definite matrix with negative eigenvalues, while the other
estimators guarantee the positive definiteness. Although RLZ performs slightly better than
the proposed method in terms of eigenvalues spread, it produces an incorrect sparsity pattern
shown in Figure 1.

7

Discussion

In this paper, we have introduced a positive definite ensemble estimate of covariance matrix
for the high-dimensional data. The proposed method takes advantage of the multiple estimates obtained from the modified Cholesky decomposition (MCD) approach under different
variable orders. The positive definite constraint and L1 penalty are added to the objective
function to guarantee the positive definiteness and encourage the sparse structure of the
estimated covariance matrix. An efficient algorithm is developed to solve the constraint
optimization problem. The proposed estimator does not require the prior knowledge of the
variable order used in the MCD, and performs well in the high-dimensional cases. Simulation
studies and real application demonstrate the superiority of the proposed method to other
existing approaches.
The ensemble idea plays an important role in constructing the proposed estimator in this
work. This idea may also be applied into other estimates, such as the inverse covariance
matrix estimate. However, one potential issue in the ensemble idea in practice is that the
variables may have relations among themselves, i.e. causal relationship or spatial information. It means that some orders of variables are meaningful and reflect such relations, while

23

others not. This is can be clearly seen through BIC and BPA methods in the numerical
study. Hence, ruling out the meaningless orders and only using the meaningful orders of
variables would improve the performance of the proposed method. How to implement this
idea in the real data needs further study.

8

Appendix

Proof of Lemma 1
Since (Σ+ , Φ+ , Λ+ ) is the optimal minimizer of (3.5), based on the Karush-Kuhn-Tucker
conditions (Karush, 1939; Kuhn and Tucker, 1951) we have
M
+
+
1 X
Σ̂k − Λ̂ )jl ∈ λ∂|Σ̂jl |, j = 1, . . . , p, l = 1, . . . , p, and j 6= l
(−Σ̂ +
M k=1
+

M
+
1 X
(−Σ̂ +
Σ̂k )jj + Λ̂jj = 0, j = 1, . . . , p
M k=1
+

+

Φ̂ = Σ̂

+

+

Φ̂  νI,

(8.1)

(8.2)

(8.3)

(8.4)

and
+

+

hΛ̂ , Φ − Φ̂ i ≤ 0, ∀Φ  νI.

(8.5)

The expressions in (8.1) and (8.2) result from the stationarity, and (8.3) and (8.4) are valid
because of the primal feasibility. By the optimality conditions of the problem (3.6) with
respect to Φ, we obtain
1
hΛi − (Φi+1 − Σi ), Φ − Φi+1 i ≤ 0, ∀Φ  νI.
τ
This, together with Λ step (3.8), yields
1
hΛi+1 − (Σi+1 − Σi ), Φ − Φi+1 i ≤ 0, ∀Φ  νI.
τ
24

(8.6)

+

Now by setting Φ = Φi+1 in (8.5) and Φ = Φ̂ in (8.6) respectively, it leads to
+

+

hΛ̂ , Φi+1 − Φ̂ i ≤ 0,

(8.7)

+
1
hΛi+1 − (Σi+1 − Σi ), Φ̂ − Φi+1 i ≤ 0.
τ

(8.8)

and

Summing (8.7) and (8.8) gives
+
+
1
h(Λi+1 − Λ̂ ) − (Σi+1 − Σi ), Φi+1 − Φ̂ i ≥ 0.
τ

(8.9)

On the other hand, by the optimality conditions of the problem (3.7) with respect to Σ, we
have
M
1 X i+1
1
0∈[
(Σ − Σ̂k ) + Λi + (Σi+1 − Φi+1 )]jl + λ∂|Σi+1
jl |, j = 1, . . . , p, l = 1, . . . , p, and j 6= l,
M k=1
τ

(8.10)
and
M
1 X i+1
1
[
(Σ − Σ̂k ) + Λi + (Σi+1 − Φi+1 )]jj = 0, j = 1, . . . , p,
M k=1
τ

(8.11)

Plugging Λ step (3.8) into (8.10) and (8.11) respectively results in
(−Σi+1 +

M
1 X
Σ̂k − Λi+1 )jl ∈ λ∂|Σi+1
jl |, j = 1, . . . , p, l = 1, . . . , p, and j 6= l,
M k=1

(8.12)

and
(Σi+1 −

M
1 X
Σ̂k )jj + Λi+1
jj = 0, j = 1, . . . , p.
M k=1

(8.13)

Since ∂| · | is monotonically non-decreasing, (8.1) and (8.12) yield for j 6= l

M
 ≥ (−Σ̂+ + 1 PM Σ̂ − Λ̂+ ) , if Σi+1 ≥ Σ̂+
1 X
k
jl
jl
jl
k=1
M
i+1
i+1
(−Σ +
Σ̂k − Λ )jl
P
+
+
+ ,
M
i+1
 ≤ (−Σ̂ + 1
M k=1
k=1 Σ̂k − Λ̂ )jl , if Σjl < Σ̂jl
M
that is,

 ≥ 0, if Σi+1 ≥ Σ̂+
+
+
jl
jl
i+1
i+1
(Σ̂ − Σ + Λ̂ − Λ )jl
.
+
 ≤ 0, if Σi+1 < Σ̂
jl
jl
25

As a result, we obtain
+

+

+

(Σi+1 − Σ̂ )jl (Σ̂ − Σi+1 + Λ̂ − Λi+1 )jl ≥ 0, j = 1, . . . , p, l = 1, . . . , p, and j 6= l. (8.14)
In addition, subtracting (8.13) from (8.2) implies
+

+

(Σ̂ − Σi+1 + Λ̂ − Λi+1 )jj = 0, j = 1, . . . , p.

(8.15)

Then combining (8.14) and (8.15) leads to
+

+

+

hΣi+1 − Σ̂ , Σ̂ − Σi+1 + Λ̂ − Λi+1 i ≥ 0.

(8.16)

By summing (8.9) and (8.16), we have
i
+
+
+
+
+
+
1
hΣi+1 − Σ̂ , Λ̂ − Λi+1 i + hΛi+1 − Λ̂ , Φi+1 − Φ̂ i − hΣi+1 − Σ̂ , Φi+1 − Φ̂ i ≥ kΣi+1 − Σ̂ k2F .
τ

This, together with (8.3) and Φi+1 = τ (Λi − Λi+1 ) + Σi+1 from Λ step (3.8), gives
+
+
1
τ hΛi+1 − Λ̂ , Λi − Λi+1 i + hΣi+1 − Σ̂ , Σi − Σi+1 i
τ
+

≥ kΣi+1 − Σ̂ k2F − hΛi − Λi+1 , Σi − Σi+1 i
+

+

+

(8.17)

+

By Φ̂ − Φi+1 = (Φ̂ − Φi ) + (Φi − Φi+1 ) and Σ̂ − Σi+1 = (Σ̂ − Σi ) + (Σi − Σi+1 ), (8.17)
is reduced to
+
+
1
τ hΛi − Λ̂ , Λi − Λi+1 i + hΣi − Σ̂ , Σi − Σi+1 i ≥ τ kΛi − Λi+1 k2F
τ
+
1 i
i+1 2
+ kΣ − Σ kF + kΣi+1 − Σ̂ k2F − hΛi − Λi+1 , Σi − Σi+1 i
τ

Using the notations W + and W i , the left hand side of (8.18) becomes
+
+
1
h(Λi − Λ̂ , Σi − Σ̂ )0 , [τ (Λi − Λi+1 ), (Σi − Σi+1 )]0 i
τ
+

+

= h(Λi , Σi )0 − (Λ̂ , Σ̂ )0 , J [(Λi , Σi )0 − (Λ̂
= hW i − W + , J (W i − W i+1 )i
= hW i − W + , W i − W i+1 iJ .

26

i+1

i+1 0

, Σ̂

) ]i

(8.18)

The first two terms on the right side of (8.18) becomes
1
1
τ kΛi − Λi+1 k2F + kΣi − Σi+1 k2F = τ hΛi − Λi+1 , Λi − Λi+1 i + hΣi − Σi+1 , Σi − Σi+1 i
τ
τ
1
i
i+1
i
i+1 0
= h(Λ − Λ , Σ − Σ ) , [τ (Λi − Λi+1 ), (Σi − Σi+1 )]0 i
τ
= h(Λi , Σi )0 − (Λi+1 , Σi+1 )0 , J [(Λi , Σi )0 − (Λi+1 , Σi+1 )0 ]i
= hW i − W i+1 , J (W i − W i+1 )i
= kW i − W i+1 k2J .
As a result, (8.18) can be rewritten as
+

hW i − W + , W i − W i+1 iJ ≥ kW i − W i+1 k2J + kΣi+1 − Σ̂ k2F − hΛi − Λi+1 , Σi − Σi+1 i.
(8.19)
Note a fact that
kW + − W i+1 k2J = kW + − W i k2J − 2hW + − W i , W i+1 − W i iJ + kW i − W i+1 k2J .
Therefore,
kW + − W i k2J − kW + − W i+1 k2J
= 2hW + − W i , W i+1 − W i iJ − kW i − W i+1 k2J
+

≥ 2kW i − W i+1 k2J + 2kΣi+1 − Σ̂ k2F − 2hΛi − Λi+1 , Σi − Σi+1 i − kW i − W i+1 k2J
+

= kW i − W i+1 k2J + 2kΣi+1 − Σ̂ k2F + 2hΛi+1 − Λi , Σi − Σi+1 i.

(8.20)

Hence, next we only need to show hΛi+1 − Λi , Σi − Σi+1 i ≥ 0. Now replacing i instead of
i + 1 in (8.12) and (8.13) yields
M
1 X
Σ̂k − Λi )jl ∈ λ∂|Σijl |, j = 1, . . . , p, l = 1, . . . , p, and j 6= l,
(−Σ +
M k=1
i

(8.21)

and
M
1 X
(Σ −
Σ̂k )jj + Λijj = 0, j = 1, . . . , p.
M k=1
i

27

(8.22)

So (8.12), (8.13), (8.21) and (8.22), together with the monotonically non-decreasing property
of ∂| · |, imply
hΣi − Σi+1 , Λi+1 − Λi + Σi+1 − Σi i ≥ 0.

(8.23)

After a simple algebra of (8.23), we have
hΣi − Σi+1 , Λi+1 − Λi i ≥ kΣi+1 − Σi k2F ≥ 0.
Hence the last two terms on the right hand side of (8.20) are both non-negative, which proves
Lemma 1.



Proof of Theorem 1
According to Lemma 1, we have
(a) kW i − W i+1 k2J → 0, as i → +∞;
(b) kW + − W i k2J is non-increasing and thus bounded.
The result (a) indicates that Σi − Σi+1 → 0 and Λi − Λi+1 → 0. Based on (3.8) it is easy to
see that Φi − Σi → 0. On the other hand, (b) indicates that W i lies in a compact region.
Accordingly, there exists a subsequence W ij of W i such that W ij → W ∗ = (Λ∗ , Σ∗ ). In
addition, we also have Φij → Φ∗ , Σ∗ . Therefore, lim (Σi , Φi , Λi ) = (Σ∗ , Φ∗ , Λ∗ ).
i→∞

∗

∗

∗

Next we show that (Σ , Φ , Λ ) is an optimal solution of (3.3). By letting i → +∞ in
(8.12), (8.13) and (8.6), we have
M
1 X
Σ̂k − Λ∗ )jl ∈ λ∂|Σ∗jl |, j = 1, . . . , p, l = 1, . . . , p, and j 6= l,
(−Σ +
M k=1
∗

(8.24)

M
1 X
Σ̂k )jj + Λ∗jj = 0, j = 1, . . . , p,
(Σ −
M k=1

(8.25)

hΛ∗ , Φ − Φ∗ i ≤ 0, ∀Φ  νI.

(8.26)

∗

and

(8.24), (8.25) and (8.26), together with Φ∗ = Σ∗ , imply that (Σ∗ , Φ∗ , Λ∗ ) is an optimal
solution of arg min L(Σ, Φ; Λ) in (3.5). Hence, we prove that the sequence produced by
28

Algorithm 1 from any starting point converges to an optimal minimizer of (3.5).



Proof of Lemma 2
The proof is very similar to that of Lemma A.2 in Jiang (2012), so we omit them.



Proof of Lemma 3
Since the conclusions of this lemma are invariant to the variable orders, we prove it under the
original order without the symbol πk to simplify the notations. From  = L−1 X ∼ N (0, D),
P
the negative log likelihood can be written as ni=1 [log |D|+tr(x0i L0−1 D −1 L−1 xi )], up to some
constant. Consequently, adding the penalty terms to the negative log likelihood leads to the
following objective function
n
X

[log |D| +

tr(x0i L0−1 D −1 L−1 xi )]

i=1

+

p
X

ηj

j=1

X

|ljk |.

k<j

Denote
Q(D, L) = (log |D| + tr(L

0−1

−1

−1

D L S) +

p
X
j=1

ηj

X

|ljk |.

k<j

Define G(∆L , ∆D ) = Q(D 0 + ∆D , L0 + ∆L ) − Q(D 0 , L0 ). Let AU1 = {∆L : k∆L k2F ≤
U12 s1 log(p)/n} and BU2 = {∆D : k∆D k2F ≤ U22 p log(p)/n}, where U1 and U2 are constants.
We will show that for each ∆L ∈ ∂AU1 and ∆D ∈ ∂BU2 , probability P (G(∆L , ∆D ) > 0)
is tending to 1 as n → ∞ for sufficiently large U1 and U2 , where ∂AU1 and ∂BU2 are the
boundaries of AU1 and BU2 , respectively. Additionally, since G(∆L , ∆D ) = 0 when ∆L = 0
and ∆D = 0, the minimum point of G(∆L , ∆D ) is achieved when ∆L ∈ AU1 and ∆D ∈ BU2 .
That is k∆L k2F = Op (s1 log(p)/n) and k∆D k2F = Op (p log(p)/n).
Assume k∆L k2F = U12 s1 log(p)/n and k∆D k2F = U22 p log(p)/n. From assumption (4.1)
and by Lemma 2, without loss of generality, there exists a constant h such that 0 < 1/h <
svp (L0 ) ≤ sv1 (L0 ) < h < ∞ and 0 < 1/h < svp (D 0 ) ≤ sv1 (D 0 ) < h < ∞. Write
D = D 0 + ∆D and L = L0 + ∆L , then we decompose G(∆L , ∆D ) into three parts and then

29

consider them separately.
G(∆L , ∆D ) = Q(D, L) − Q(D 0 , L0 )
−1 −1
= log |D| − log |D 0 | + tr(L0−1 D −1 L−1 S) − tr(L0−1
0 D 0 L0 S)

+

p
X

ηj

j=1

X

|ljk | −

p
X

ηj

j=1

k<j

X

|l0jk |

k<j

−1
= log |D| − log |D 0 | + tr[(D −1 − D −1
− D −1
0 )D 0 ] − tr[(D
0 )D 0 ]

+ tr(L

0−1

−1

−1

D L S) −

−1 −1
tr(L0−1
0 D 0 L0 S)

+

p
X
j=1

ηj

X

|ljk | −

k<j

p
X

ηj

j=1

X

|l0jk |

k<j

= M1 + M2 + M3 ,
where
M1 = log |D| − log |D 0 | + tr[(D −1 − D −1
0 )D 0 ],
−1 −1
−1
M2 = tr(L0−1 D −1 L−1 S) − tr(L0−1
− D −1
0 D 0 L0 S) − tr[(D
0 )D 0 ],

M3 =

p
X
j=1

ηj

X
k<j

|ljk | −

p
X
j=1

ηj

X

|l0jk |.

k<j

Based on the proof of Theorem 3.1 in Jiang (2012), we can have M1 ≥ 1/8h4 k∆D k2F . For
the second term,
−1
0−1
−1
M2 = tr(L0−1 D −1 L−1 S) − tr(L0−1 D −1
D −1
0 L S) + tr(L
0 L S)
−1 −1
−1
− tr(L0−1
− D −1
0 D 0 L0 S) − tr[(D
0 )D 0 ]
−1
0−1
−1
0−1
0−1
= tr(D −1 − D −1
] + trD −1
− L−1
0 )[L (S − Σ0 )L
0 (L SL
0 SL0 )
−1
0−1
+ tr(D −1 − D −1
− D0 )
0 )(L Σ0 L
−1
0−1
−1
0−1
0−1
= tr(D −1 − D −1
] + tr[D −1
− L−1
0 )[L (S − Σ0 )L
0 (L (S − Σ0 )L
0 (S − Σ0 )L0 )]
−1
0−1
0−1
−1
−1
0−1
− L−1
− D −1
− D0 )
+ tr[D −1
0 (L Σ0 L
0 Σ0 L0 )] + tr(D
0 )(L Σ0 L
−1
0−1
−1
0−1
0−1
= tr(D −1 − D −1
] + tr[D −1
− L−1
0 )[L (S − Σ0 )L
0 (L (S − Σ0 )L
0 (S − Σ0 )L0 )]
0−1
+ tr[D −1 (L−1 Σ0 L0−1 − L−1
0 Σ0 L0 )]
(1)

(2)

(3)

= M2 + M2 + M2 ,

30

0

(1)

0−1
0−1
where the fourth equality uses L−1
= L−1
= D 0 . The notations M2 ,
0 Σ0 L0
0 (L0 D 0 L0 )L0
(2)

M2

(3)

and M2

are defined in the following
(1)

−1
0−1
M2 = tr(D −1 − D −1
],
0 )[L (S − Σ0 )L
(2)

−1
0−1
0−1
M2 = tr[D −1
− L−1
0 (L (S − Σ0 )L
0 (S − Σ0 )L0 )],
(3)

0−1
M2 = tr[D −1 (L−1 Σ0 L0−1 − L−1
0 Σ0 L0 )].

Based on the proof of Theorem 3.1 in Jiang (2012), for any  > 0, there exists V1 > 0 and
V2 > 0 such that
(1)

|M2 | ≤ V1

p
p log(p)/nk∆D kF

and
(3)

(2)

M2 − |M2 | > 1/2h4 k∆L k2F − V2

X
p
p
log(p)/n
|ljk | − V2 h s1 log(p)/nk∆L kF ,
(j,k)∈Z c

where Z = {(j, k) : k < j, l0jk 6= 0}, and l0jk represents the element (j, k) of the matrix L0 .
Next, for the penalty term,
M3 =

p
X

X

ηj

j=1

|ljk | +

(j,k)∈Z c

p
X

ηj

j=1

X

(1)

(2)

(|ljk | − |l0jk |) = M3 + M3 ,

(j,k)∈Z

where
(1)

M3 =

p
X

ηj

X

|ljk |,

(j,k)∈Z c

j=1

and
(2)
|M3 |

=|

p
X
j=1

ηj

X

(|ljk | − |l0jk |)| ≤

p
X

ηj

j=1

(j,k)∈Z

≤

p
X
j=1
p

≤

X
j=1

31

X

|(|ljk | − |l0jk |)|

(j,k)∈Z

ηj

X

|ljk − l0jk |

(j,k)∈Z

√
ηj s1 k∆L kF ,

where the last inequality uses the fact that (a1 + a2 + · · · + am )2 ≤ m(a21 + a22 + · · · + a2m ).
Combine all the terms above together, with probability greater than 1 − 2, we have
|G(∆L , ∆D )|
(1)

(3)

(2)

(1)

(2)

≥ M1 − |M2 | + M2 − |M2 | + M3 − |M3 |
X
p
p
≥ 1/8h4 k∆D k2F − V1 p log(p)/nk∆D kF + 1/2h4 k∆L k2F − V2 log(p)/n
|ljk |
(j,k)∈Z c
p

p

X
X
X √
p
ηj
ηj s1 k∆L kF
|ljk | −
− V2 h s1 log(p)/nk∆L kF +
j=1

=

(j,k)∈Z c

j=1

X
p
U22
U12
log(p)/n
|ljk |
p
log(p)/n
−
V
U
p
log(p)/n
+
s
log(p)/n
−
V
1
2
1
2
8h4
2h4
c
(j,k)∈Z

p

− V2 U1 hs1 log(p)/n +

X
j=1

p

ηj

X

|ljk | − s1 U1

p

log(p)/n

X

(j,k)∈Z c

ηj

j=1

Pp

=

U1 s1 p log(p) U1
U2 p log(p) U2
j=1 ηj
( 4 − V1 ) +
( 4−p
− V2 h)
n
8h
n
2h
log(p)/n
p
X
X
p
+
|ljk |(
ηj − V2 log(p)/n).
(j,k)∈Z c

j=1

Here V1 and V2 are only related to the sample size n and . Assume

Pp

j=1

ηj = K(log(p)/n)

where K > V2 and choose U1 > 2h4 (K + hV2 ), U2 > 8h4 V1 , then G(∆L , ∆D ) > 0. This
establishes the lemma.



Proof of Theorem 2
Based on the proof of Theorem 3.2 in Jiang (2012), we can have
kΣ̂πk − Σ0πk k2F = Op (kL̂πk − L0πk k2F ) + Op (kD̂ πk − D 0πk k2F )
= Op (s1 log(p)/n) + Op (p log(p)/n)
= Op ((s1 + p) log(p)/n),

32

where the second inequality is provided by Lemma 3. Then
kΣ̂k − Σ0 k2F = kP πk Σ̂πk P 0πk − P πk Σ0πk P 0πk k2F
= kP πk (Σ̂πk − Σ0πk )P 0πk k2F
= kΣ̂πk − Σ0πk k2F
= Op ((s1 + p) log(p)/n),
where the third equality uses the fact that the Frobenius norm of a matrix is invariant on
the permutation matrix.



Proof of Theorem 3 and 4
As Σ0 is positive definite, there exists  > 0 such that  < λmin (Σ0 ), where λmin (Σ0 ) is
the smallest eigenvalue of Σ0 . By introducing ∆ = Σ − Σ0 , the expression of (3.3) can be
rewritten in terms of ∆ as
ˆ = arg
∆

M
1 X
min
k∆ + Σ0 − Σ̂k k2F + λ|∆ + Σ0 |1 (, F(∆)).
0
2M
∆=∆ ,∆+Σ0 I
k=1

ˆ = Σ̂+ − Σ0 . Now consider ∆ ∈ {∆ : ∆ = ∆0 , ∆ + Σ0 
Note that it is easy to see that ∆
√
I, k∆kF = 5λ s0 + p}. Define the active set of Σ0 as A0 = {(i, j) : σij0 6= 0, i 6= j}, and
B A0 = (bij · I {(i,j)∈A0 } )1≤i,j≤p . Let Ac0 be the complement set of A0 . Denote the element (i, j)
of matrix ∆ by ∆ij . Under the probability event {|σ̂ijk − σij0 | ≤ λ} where Σ̂k = (σ̂ijk )p×p , we
have

33

M
M
1 X
1 X
2
F(∆) − F(0) =
k∆ + Σ0 − Σ̂k kF −
kΣ0 − Σ̂k k2F + λ|∆ + Σ0 |1 − λ|Σ0 |1
2M k=1
2M k=1
M
1
1 X
2
= k∆kF +
< ∆, Σ0 − Σ̂k > +λ|∆Ac0 |1 + λ(|∆A0 + (Σ0 )A0 |1 − |(Σ0 )A0 |1 )
2
M k=1
X
1
≥ k∆k2F − λ(|∆|1 +
∆ii ) + λ|∆Ac0 |1 − λ|∆A0 |1
2
i
X
1
∆ii )
≥ k∆k2F − 2λ(|∆A0 |1 +
2
i
√
1
≥ k∆k2F − 2λ s0 + pk∆kF
2
5 2
= λ (s0 + p)
2

> 0.
ˆ is also the optimal solution to the convex optimization problem
Note that ∆
ˆ = arg
∆

min
F(∆) − F(0).
∆=∆0 ,∆+Σ0 I

The rest of proof is the same as that of Theorem 2 in Xue, Ma and Zou (2012), so we omit
them.



References
Aubry, A., De Maio, A., Pallotta, L., and Farina, A. (2012). Maximum Likelihood Estimation of a Structured Covariance Matrix with a Condition Number Constraint. Signal
Processing, IEEE Transactions on, 60(6), 3004–3021.
Bickel, P. J., and Levina, E. (2004). Some Theory of Fisher’s Linear Discriminant Function,Naive Bayes, and Some Alternatives When There Are Many More Variables than
Observations. Bernoulli, 10(6), 989–1010.
Bickel, P. J., and Levina, E. (2009). Covariance Regularization by Thresholding. The
Annals of Statistics, 36(6), 2577–2604.
34

Bien, J., and Tibshirani, R. J. (2011). Sparse Estimation of a Covariance Matrix. Biometrika,
98(4), 807–820.
Boyd, S., Parikh, N., Chu, E., Peleato, B., and Eckstein, J. (2011). Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers.
Foundations and Trends in Machine Learning, 3(1), 1–122.
Cai, T. T., and Yuan, M. (2012). Adaptive Covariance Matrix Estimation through Block
Thresholding. The Annals of Statistics, 40(40), 2014–2042.
Cai, T. T., Zhang, C. H., and Zhou, H. H. (2010). Optimal Rates of Convergence for
Covariance Matrix Estimation. The Annals of Statistics, 38(4), 2118–2144.
Cai, T. T., Ren, Z., and Zhou, H. H. (2016). Estimating Structured High-Dimensional Covariance and Precision Matrices: Optimal Rates and Adaptive Estimation. Electronic
Journal of Statistics, 10(1), 1–59.
Dellaportas, P., and Pourahmadi M. (2012). Cholesky-GARCH Models with Applications
to Finance. Statistics and Computing, 22(4), 849–855.
Deng, X., and Tsui, K. W. (2013). Penalized Covariance Matrix Estimation Using a MatrixLogarithm Transformation. Journal of Computational and Graphical Statistics, 22(2),
494–512.
Deng, X., and Yuan, M. (2009). Large Gaussian Covariance Matrix Estimation with Markov
Structure, Journal of Computational and Graphical Statistics, 18(3), 640–657.
Dey, D. K., and Srinivasan, C. (1985). Estimation of a Covariance Matrix under Stein’s
Loss. The Annals of Statistics, 13(4), 1581–1591.
Fan, J., Liao, Y., and Mincheva, M. (2013). Large Covariance Estimation by Thresholding
Principal Orthogonal Complements. Journal of the Royal Statistical Society: Series
B, 75(4), 603–680.

35

Fan, J., Liao, Y., and Liu, H. (2016). An Overview of the Estimation of Large Covariance
and Precision Matrices. The Econometrics Journal, 19(1), C1–C32.
Friedman, J., Hastie, T., and Tibshirani, T. (2008). Sparse Inverse Covariance Estimation
with the Graphical Lasso. Biostatistics, 9(3), 432–441.
Glaab, E., Bacardit, J., Garibaldi, J. M., and Krasnogor, N. (2012). Using Rule-Based
Machine Learning for Candidate Disease Gene Prioritization and Sample Classification
of Cancer Gene Expression Data. PloS one, 7(7), e39932.
Guo, J., Levina, E., Michailidis, G., and Zhu, J. (2011). Joint Estimation of Multiple
Graphical Models. Biometrika, 98(1), 1–15.
Haff, L. R. (1991). The Variational Form of Certain Bayes Estimators. The Annals of
Statistics, 19(3), 1163–1190.
Huang, J. Z., Liu, N., Pourahmadi, M., and Liu, L. (2006). Covariance Matrix Selection
and Estimation via Penalised Normal Likelihood. Biometrika, 93(1), 85–98.
Jiang, X. (2012). Joint Estimation of Covariance Matrix via Cholesky Decomposition. Ph.D
Dissertation. Department of Statistics and Applied Probability, National University
of Singapore.
Johnstone, I. M. (2001). On the Distribution of the Largest Eigenvalue in Principal Components Analysis. The Annals of Statistics, 29(2), 295–327.
Kang, X., and Deng, X. (2017) An Improved Modified Cholesky Decomposition Method
for Inverse Covariance Matrix Estimation. https://arxiv.org/abs/1710.05163.
Karush, W. (1939). Minima of Functions of Several Variables with Inequalities as Side
Conditions. Master Dissertation. Department of Mathematics, University of Chicago,
Chicago, Illinois.

36

Kuhn, H., and Tucker, A. (1951). Nonlinear Programming. Proceedings of the 2nd Berkeley Symposium on Mathematica Statistics and Probabilistics, University of California
Press, 481–492.
Kullback, S., and Leibler, R. A. (1951). On Information and Sufficiency. The Annals of
Mathematical Statistics, 22(1), 79–86.
Lam, C., and Fan, J. (2009). Sparsistency and Rates of Convergence in Large Covariance
Matrix Estimation. The Annals of Statistics, 37(6B), 4254–4278.
Lange, K., Hunter, D. R., and Yang, I. (2000). Optimization Transfer Using Surrogate
Objective Functions. Journal of Computational and Graphical Statistics, 9(1), 1–20.
Ledoit, O., and Wolf, M. (2004). A Well-Conditioned Estimator for Large-Dimensional
Covariance Matrices. Journal of Multivariate Analysis, 88(2), 365–411.
Levina, E., Rothman, A., and Zhu, J. (2008). Sparse Estimation of Large Covariance
Matrices via a Nested Lasso Penalty. The Annals of Applied Statistics, 2(1), 245–263.
Liu, H., Wang, L., and Zhao, T. (2014). Sparse Covariance Matrix Estimation with Eigenvalue Constraints. Journal of Computational and Graphical Statistics, 23(2), 439–459.
Pourahmadi, M. (1999). Joint Mean-Covariance Models with Applications to Longitudinal
Data: Unconstrained Parameterisation. Biometrika, 86(3), 677–690.
Pourahmadi, M. (2013). High-Dimensional Covariance Estimation: with High-Dimensional
Data. John Wiley & Sons, Chichester, UK.
Pourahmadi, M., Daniels, M. J., and Park, T. (2007). Simultaneous Modelling of the
Cholesky Decomposition of Several Covariance Matrices. Journal of Multivariate Analysis, 98(3), 568–587.
Rajaratnam, B., and Salzman, J. (2013). Best Permutation Analysis. Journal of Multivariate Analysis, 121(10), 193–223.

37

Rocha, G. V., Zhao, P., and Yu, B. (2008). A Path Following Algorithm for Sparse PseudoLikelihood Inverse Covariance Estimation. Technical Report.
Rothman, A., Bickel, P., Levina, E., and Zhu, J. (2008). Sparse Permutation Invariant
Covariance Estimation. Electronic Journal of Statistics, 2(3), 494–515.
Rothman, A. J., Levina, E., and Zhu, J. (2009). Generalized Thresholding of Large Covariance Matrices. Journal of the American Statistical Association, 104(485), 177–186.
Rothman, A. J., Levina, E., and Zhu, J. (2010). A New Approach to Cholesky-Based
Covariance Regularization in High Dimensions. Biometrika, 97(3), 539–550.
Tibshirani, R. (1996). Regression Shrinkage and Selection via the Lasso. Journal of the
Royal Statistical Society, Series B, 58(1), 267–288.
Wu, W. B., and Pourahmadi, M. (2003). Nonparametric Estimation of Large Covariance
Matrices of Longitudinal Data. Biometrika, 90(4), 831–844.
Won, J. H., Lim, J., Kim, S. J., and Rajaratnam, B. (2013). Condition - Number Regularized Covariance Estimation. Journal of the Royal Statistical Society: Series B,
75(3), 427–450.
Xiao, L., Zipunnikov, V., Ruppert, D., and Crainiceanu, C. (2016). Fast Covariance Estimation for High-Dimensional Functional Data. Statistics and Computing, 26(1-2),
409–421.
Xue, L., Ma, S., and Zou, H. (2012). Positive-Definite L1 -Penalized Estimation of Large Covariance Matrices. Journal of the American Statistical Association, 107(500), 1480–
1491.
Yuan, M., and Lin, Y. (2007). Model Selection and Estimation in the Gaussian Graphical
Model. Biometrika, 94(1), 19–35.
Yuan, M. (2008). Efficient Computation of the `1 Regularized Solution Path in Gaussian
Graphical Models. Journal of Computational and Graphical Statistics, 17(4), 809–826.
38

Yuan, M. (2010). High Dimensional Inverse Covariance Matrix Estimation via Linear
Programming. The Journal of Machine Learning Research, 11(12), 2261–2286.

39

