Positive Definite $\ell_1$ Penalized Estimation of Large Covariance Matrices

Introduction

Estimating covariance matrices is of fundamental importance for an abundance of statistical methodologies. Nowadays, the advance of new technologies has brought massive high-dimensional data into various research fields, such as fMRI imaging, web mining, bioinformatics, climate studies and risk management, and so on. The usual sample covariance matrix is optimal in the classical setting with large samples and fixed low dimensions (Anderson, 1984), but it performs very poorly in the high-dimensional setting (Johnstone, 2001). In the recent literature, regularization techniques have been used to improve the sample covariance matrix estimator, including banding (Wu and Pourahmadi, 2003; Bickel and Levina, 2008a), tapering (Furrer and Bengtsson, 2007; Cai, Zhang, and Zhou, 2010) and thresholding (Bickel and Levina, 2008b; El Karoui, 2008; Rothman, Levina, and Zhu, 2009). Banding or tapering is very useful when the variables have a natural ordering and off-diagonal entries of the target covariance matrix decays to zero as they move away from the diagonal. On the other hand, thresholding is proposed for estimating permutation-invariant covariance matrices. Thresholding can be used to produce consistent covariance matrix estimators when the true covariance matrix is bandable (Bickel and Levina, 2008b; Cai and Zhou, 2011a). In this sense, thresholding is more robust than banding/tapering for real applications.

However, there is no guarantee that the thresholding estimator is always positive definite. Although the positive definite property is guaranteed in the asymptotic setting with high probability, the actual estimator can be an indefinite matrix, especially in real data analysis. To illustrate this issue, we consider the Michigan lung cancer gene-expression data (Beer et al., 2002) which have $86$ tumor samples from patients with lung adenocarcinomas and $5217$ gene expression values for each sample. More details about this dataset are referred to Beer et al. (2002) and Subramaniana et al. (2005). We randomly choose $p$ genes ($p=200,500$), and obtain the soft-thresholding sample correlation matrix for these genes. We repeat the process ten times for $p=200$ and $500$ respectively, and each time the thresholding parameter $\lambda$ is selected via the 5-fold cross validation. We found that none of the soft-thresholding estimators would become positive definite for both $p=200$ and $500$. On average, there exist $22$ and $124$ negative eigenvalues for the soft-thresholding estimator for $p=200$ and $p=500$, respectively. Figure 1 displays the $30$ smallest eigenvalues for $p=200$ and the $130$ smallest eigenvalues for $p=500$.

Note that the solution to (2) could be positive semidefinite. To obtain a positive definite covariance estimator, we can consider the positive definite constraint $\{\boldsymbol{\Sigma}\succeq\epsilon\boldsymbol{I}\}$ for some arbitrarily small $\epsilon>0$. Then the modified $\hat{\boldsymbol{\Sigma}}^{+}$ is always positive definite. In this work, we focus on solving the positive definite $\hat{\boldsymbol{\Sigma}}^{+}$ as follows

Despite its natural motivation, (3) is actually a very challenging optimization problem due to the positive semidefinite constraint. To our best knowledge, the first attempt for solving (3) was recently proposed by Rothman (2011) who added the log-determinant barrier function to (3):

Alternating Direction Algorithm

We use an alternating direction method to solve (3) directly. The alternating direction method is closely related to the operator-splitting method that has a long history back to 1950s for solving numerical partial differential equations, see e.g., Douglas and Rachford (1956); Peaceman and Rachford (1955). Recently, the alternating direction method has been revisited and successfully applied to solving large scale problems arising from different applications. For example, Scheinberg, Ma, and Goldfarb (2010) introduced the alternating linearization methods to efficiently solve the graphical lasso optimization problem. We refer to Fortin and Glowinski (1983); Glowinski and Le Tallec (1989) for more details on operator-splitting and alternating direction methods.

The solution to (5) gives the solution to (3). To deal with the equality constraint in (5), we shall minimize its augmented Lagrangian function for some given penalty parameter $\mu$, i.e.

where $\boldsymbol{\Lambda}$ is the Lagrange multiplier. We iteratively solve

and then update the Lagrangian multiplier $\boldsymbol{\Lambda}^{i+1}$ by

For (7) we do it by alternatingly minimizing $L(\boldsymbol{\Theta},\boldsymbol{\Sigma};\boldsymbol{\Lambda}^{i})$ with respect to $\boldsymbol{\Theta}$ and $\boldsymbol{\Sigma}$.

To sum up, the entire algorithm proceeds as follows:

For $i=0,1,2,\ldots$, solve the following three sub-problems sequentially till convergence

Next, define an entry-wise soft-thresholding rule for all the non-diagonal elements of a matrix $\boldsymbol{Z}$ as $\boldsymbol{S}(\boldsymbol{Z},\tau)=\{s(z_{ij},\tau)\}_{1\leq i,j\leq p}$ with

Then the $\boldsymbol{\Sigma}$ step has a closed-form solution given below

Algorithm 1 shows the complete details of our alternating direction method for (3). In Section 4 we provide the convergence analysis of Algorithm 1 and prove that Algorithm 1 always converges to the optimal solution of (5) from any starting point.

In our implementation we use the soft-thresholding estimator as the initial value for both $\boldsymbol{\Theta}^{0}$ and $\boldsymbol{\Sigma}^{0}$, and we set $\boldsymbol{\Lambda}^{0}$ as a zero matrix. The value for $\mu$ is 2. Before invoking Algorithm 1, we always check whether the soft-thresholding estimator is positive definite. If yes, then the soft-threhsolding estimator is the final solution to (3).

Numerical Examples

Before delving into theoretical analysis of the algorithm and the resulting estimator, we first use simulation to show the competitive performance of our proposal. In all examples we standardize the variables to have zero mean and unit variance. In each simulation model, we generated $100$ independent datasets, each with $n=50$ independent $p$-variate random vectors from the multivariate normal distribution with mean and covariance matrix $\boldsymbol{\Sigma}_{0}=(\sigma^{0}_{ij})_{1\leq i,j\leq p}$ for $p=100,200~{}\&~{}500$. We considered two covariance models with different sparsity patterns:

partition the indices $\{1,2,\ldots,p\}$ into $K=p/20$ non-overlapping subsets of equal size, and let $i_{k}$ denote the maximum index in $I_{k}$.

Model 1 has been used in Bickel and Levina (2008a) and Cai and Liu (2011), and Model 2 is similar to the overlapping block diagonal design used in Rothman (2011).

First, we compare the run times of our estimator $\hat{\boldsymbol{\Sigma}}^{+}$ with the log-barrier estimator $\breve{\boldsymbol{\Sigma}}^{+}$ by Rothman (2011). As shown in Table 1, our method is much faster than the log-barrier method.

In what follows, we compare the performance of $\hat{\boldsymbol{\Sigma}}^{+}$, $\breve{\boldsymbol{\Sigma}}^{+}$ and the soft-thresholding estimator $\hat{\boldsymbol{\Sigma}}$. For all three regularized estimators, the thresholding parameter was chosen by $5$-fold cross-validation (Bickel and Levina, 2008b; Rothman et al., 2009; Cai and Liu, 2011). The estimation performance is measured by the average losses under both the Frobenius norm and the spectral norm. The selection performance is examined by the false positive rate

Moreover, we compare the average number of negative eigenvalues and the percentage of positive-definiteness to check the positive-definiteness.

Table 2 and Table 3 show the average metrics over 100 replications. The soft-thresholding estimator $\hat{\boldsymbol{\Sigma}}$ is positive definite in 19 or fewer out of 100 simulation runs, while $\hat{\boldsymbol{\Sigma}}^{+}$ and $\breve{\boldsymbol{\Sigma}}^{+}$ can always guarantee a positive-definite estimator. The larger the dimension, the less likely for the soft-thresholding estimator to be positive definite. In terms of estimation, both $\hat{\boldsymbol{\Sigma}}^{+}$ and $\breve{\boldsymbol{\Sigma}}^{+}$ are more accurate than $\hat{\boldsymbol{\Sigma}}$. As for the selection performance, $\hat{\boldsymbol{\Sigma}}^{+}$ and $\breve{\boldsymbol{\Sigma}}^{+}$ achieve a slightly better true positive rate than $\hat{\boldsymbol{\Sigma}}$. Overall, $\hat{\boldsymbol{\Sigma}}^{+}$ is the best among all three regularized estimators.

2 Real data

To demonstrate our proposal we further consider two gene expression datasets: one from a small round blue-cell tumors microarray experiment (Khan et al., 2001) and the other one from a cardiovascular microarray study (Efron, 2009, 2010). The first dataset has 64 training tissue samples with four types of tumors (23 EWS, 8 BL-NHL, 12 NB, and 21 RMS), and 6567 gene expression values for each sample. We applied the pre-filtering step used in Khan et al. (2001) and then picked the top 40 and bottom 160 genes based on the F-statistic as done in Rothman et al. (2009). The second dataset has 63 subjects with 44 healthy controls and 19 cardiovascular patients, and 20426 genes measured for each subject. We used the F-statistic to pick the top 50 and bottom 150 genes. By doing so, it is expected that there is weak dependence between the top and the bottom genes. We considered the soft-thresholding estimator (Bickel and Levina, 2008b), the log-barrier estimator (Rothman, 2011) and our estimator. For all three estimators, the thresholding parameter was chosen by 5-fold cross validation.

As evidenced in Plot 2, the soft-thresholding estimator yields an indefinite matrix for both real examples whereas the other two regularized estimators guarantee the positive-definiteness. The soft-thresholding estimator contains $37$ negative eigenvalues in the small round blue-cell data, and $46$ negative eigenvalues in the cardiovascular data. Regularized correlation matrix estimation has a natural application in clustering when the dissimilarity measure is constructed using the correlation among features. For both datasets we did hierarchical clustering using the three regularized estimators. The heat maps are shown in Figure 3 in which the estimated sparsity pattern well matches the expected sparsity pattern.

Finally, we compared the average run times over $5$ cross validations for both $\hat{\boldsymbol{\Sigma}}^{+}$ and $\breve{\boldsymbol{\Sigma}}^{+}$, as shown in Table 4. It is obvious that our proposal is much more efficient.

Theoretical properties

In this section, we prove that the sequence $(\boldsymbol{\Theta}^{i},\boldsymbol{\Sigma}^{i},\boldsymbol{\Lambda}^{i})$ produced by the alternating direction method (Algorithm 1) converges to $(\hat{\boldsymbol{\Theta}}^{+},\hat{\boldsymbol{\Sigma}}^{+},\hat{\boldsymbol{\Lambda}}^{+})$, where $(\hat{\boldsymbol{\Theta}}^{+},\hat{\boldsymbol{\Sigma}}^{+})$ is an optimal solution of (5) and $\hat{\boldsymbol{\Lambda}}^{+}$ is the optimal dual variable. This automatically implies that Algorithm 1 gives an optimal solution of (3).

We define some necessary notation for ease of presentation. Let $G$ be a $2p$ by $2p$ matrix defined as

Define the norm $\|\cdot\|_{G}^{2}$ as $\|U\|_{G}^{2}=\langle U,GU\rangle$ and the corresponding inner product $\langle\cdot,\cdot\rangle_{G}$ as $\langle U,V\rangle_{G}=\langle U,GV\rangle$. Before we give the main theorem about the global convergence of Algorithm 1, we need the following lemma.

Assume that $(\hat{\boldsymbol{\Theta}}^{+},\hat{\boldsymbol{\Sigma}}^{+})$ is an optimal solution of (5) and $\hat{\boldsymbol{\Lambda}}^{+}$ is the corresponding optimal dual variable associated with the equality constraint $\boldsymbol{\Sigma}=\boldsymbol{\Theta}$. Then the sequence $\{(\boldsymbol{\Theta}^{i},\boldsymbol{\Sigma}^{i},\boldsymbol{\Lambda}^{i})\}$ produced by Algorithm 1 satisfies

Now we are ready to give the main convergence result of Algorithm 1.

The sequence $\{(\boldsymbol{\Theta}^{i},\boldsymbol{\Sigma}^{i},\boldsymbol{\Lambda}^{i})\}$ produced by Algorithm 1 from any starting point converges to an optimal solution of (5).

2 Statistical analysis of the estimator

Define $\boldsymbol{\Sigma}^{0}$ as the true covariance matrix for the observations $\boldsymbol{X}=(X_{ij})_{n\times p}$, and define the active set of $\boldsymbol{\Sigma}^{0}=(\sigma^{0}_{jk})_{1\leq j,k\leq p}$ as $A_{0}=\{(j,k):\sigma^{0}_{jk}\neq 0{,j\neq k}\}$ with the cardinality $s=|A_{0}|$. Denote by $\boldsymbol{B}_{A_{0}}$ the Hadamard product $\boldsymbol{B}_{p\times p}\circ(I_{\{(j,k)\in A_{0}\}})_{1\leq j,k\leq p}=(b_{jk}\cdot I_{\{(j,k)\in A_{0}\}})_{1\leq j,k\leq p}$. Define $\sigma_{\max}=\max_{j}\sigma^{o}_{jj}$ as the maximal true variance in $\boldsymbol{\Sigma}^{0}$.

Assume that the true covariance matrix $\boldsymbol{\Sigma}^{0}$ is positive definite.

Under the exponential-tail condition that for all $|t|\leq\eta$ and $1\leq i\leq n,1\leq j\leq p$

we also assume that $\log p\leq n$. For any $M>0$, we pick the thresholding parameter as

With probability at least $1-3p^{-M}$, we have

Under the polynomial-tail condition that for all $\gamma>0$, $\varepsilon>0$and $1\leq i\leq n,1\leq j\leq p$

we also assume that $p\leq cn^{\gamma}$ for some $c>0$. For any $M>0$, we pick the thresholding parameter as

With probability at least $1-O(p^{-M})-3K_{2}p(\log n)^{2(1+\gamma+\varepsilon)}n^{-\gamma-\varepsilon}$, we have

Define $d=\max_{j}\sum_{k}I_{\{\sigma_{jk}\neq 0\}}$ and assume that $\sigma_{\max}$ is bounded by a fixed constant, then we can pick $\lambda=O((\log p/n)^{1/2})$ to achieve the minimax optimal rate of convergence under the Frobenius norm as in Theorem 4 of Cai and Zhou (2011b) that

However, to attain the same rate in the presence of the log-determinant barrier term, Rothman (2011) instead would require that $\sigma_{\textrm{min}}$, the minimal eigenvalue of the true covariance matrix, should be bounded away from zero by some positive constant, and also that the barrier parameter should be bounded by some positive quantity. We would like to point out that if $\sigma_{\textrm{min}}$ is bounded away from zero, then the soft-thresholding estimator $\hat{\boldsymbol{\Sigma}}_{st}$ will be positive-definite with an overwhelming probability tending to $1$, (Bickel and Levina, 2008b; Cai and Zhou, 2011a, b). Therefore the theory requiring a lower bound on $\sigma_{\textrm{min}}$ is not very appealing.

Conclusions

The soft-thresholding estimator has been shown to enjoy good asymptotic properties for estimating large sparse covariance matrices. But its positive definiteness property can be easily violated, which means the soft-thresholding estimator could be in principle an inadmissible estimator for covariance matrices. In this paper we have put the soft-thresholding estimator in a convex optimization framework and considered a natural modification by imposing the positive definiteness constraint. We have developed a fast alternating direction method to solve the constrained optimization problem and the resulting estimator retains the sparsity and positive definiteness properties simultaneously. The algorithm and the new estimator are supported by numerical and theoretical results.

Acknowledgement

We thank Adam Rothman for sharing his code. Shiqian Ma’s research is supported by the National Science Foundation postdoctoral fellowship through Institute for Mathematics and Its Applications at University of Minnesota. Hui Zou’s research is supported in part by grants from the National Science Foundation and the Office of Naval Research.

Appendix: Technical Proofs

Since $(\hat{\boldsymbol{\Theta}}^{+},\hat{\boldsymbol{\Sigma}}^{+},\hat{\boldsymbol{\Lambda}}^{+})$ is optimal to (5), it follows from the KKT conditions that the followings hold.

Note that the optimality conditions for the first subproblem in Algorithm 1, i.e. the subproblem with respect to $\boldsymbol{\Theta}$ in (8), are given by

Using the updating formula for $\boldsymbol{\Lambda}^{i}$ in Algorithm 1, i.e.,

Now by letting $\boldsymbol{\Theta}=\boldsymbol{\Theta}^{i+1}$ in (16) and $\boldsymbol{\Theta}=\hat{\boldsymbol{\Theta}}^{+}$ in (19), we can get that

The optimality conditions for the second subproblem in Algorithm 1, i.e., the subproblem with respect to $\boldsymbol{\Sigma}$ in (8) are given by

Note that by using (18), (23) and (24) can be respectively rewritten as:

Using the fact that $\partial|\cdot|$ is a monotone function, (12), (13), (25) and (26) imply

Combining (28) with $\boldsymbol{\Theta}^{i+1}=\mu(\boldsymbol{\Lambda}^{i}-\boldsymbol{\Lambda}^{i+1})+\boldsymbol{\Sigma}^{i+1}$ and $\hat{\boldsymbol{\Theta}}^{+}=\hat{\boldsymbol{\Sigma}}^{+}$ leads to

Simple algebraic derivation from (29) yields the following inequality:

Rearranging the terms on the left hand side of (30) using $\hat{\boldsymbol{\Theta}}^{+}-\boldsymbol{\Theta}^{i+1}=(\hat{\boldsymbol{\Theta}}^{+}-\boldsymbol{\Theta}^{i})+(\boldsymbol{\Theta}^{i}-\boldsymbol{\Theta}^{i+1})$ and $\hat{\boldsymbol{\Sigma}}^{+}-\boldsymbol{\Sigma}^{i+1}=(\hat{\boldsymbol{\Sigma}}^{+}-\boldsymbol{\Sigma}^{i})+(\boldsymbol{\Sigma}^{i}-\boldsymbol{\Sigma}^{i+1})$, then (28) can be reduced to

Using the notation of $U^{i}$ and $U^{*}$, (31) can be rewritten as

Combining (32) with the following identity

Now, using (25) and (26) for $i$ instead of $i+1$, we get,

Combining (25), (26), (34), (35) and using the fact that $\partial|\cdot|$ is a monotone function, we obtain,

By substituting (36) into (33), we get the desired result (11). ∎

$\|U^{i}-U^{*}\|_{G}^{2}$ is monotonically non-increasing and thus converges.

It follows from (i) that $\boldsymbol{\Lambda}^{i}-\boldsymbol{\Lambda}^{i+1}\rightarrow 0$ and $\boldsymbol{\Sigma}^{i}-\boldsymbol{\Sigma}^{i+1}\rightarrow 0$. Then (18) implies that $\boldsymbol{\Theta}^{i}-\boldsymbol{\Theta}^{i+1}\rightarrow 0$ and $\boldsymbol{\Theta}^{i}-\boldsymbol{\Sigma}^{i}\rightarrow 0$. From (ii) we obtain that, $U^{i}$ has a subsequence $\{U^{i_{j}}\}$ that converges to $\bar{U}=(\bar{\boldsymbol{\Lambda}},\bar{\boldsymbol{\Sigma}})$, i.e., $\boldsymbol{\Lambda}^{i_{j}}\rightarrow\bar{\boldsymbol{\Lambda}}$ and $\boldsymbol{\Sigma}^{i_{j}}\rightarrow\bar{\boldsymbol{\Sigma}}$. From $\boldsymbol{\Theta}^{i}-\boldsymbol{\Sigma}^{i}\rightarrow 0$ we also get that $\boldsymbol{\Theta}^{i_{j}}\rightarrow\bar{\boldsymbol{\Theta}}:=\bar{\boldsymbol{\Sigma}}$. Therefore, $(\bar{\boldsymbol{\Theta}},\bar{\boldsymbol{\Sigma}},\bar{\boldsymbol{\Lambda}})$ is a limit point of $\{(\boldsymbol{\Theta}^{i},\boldsymbol{\Sigma}^{i},\boldsymbol{\Lambda}^{i})\}$.

Note that (25) and (24) respectively imply that

(37), (38) and (39) together with $\bar{\boldsymbol{\Theta}}=\bar{\boldsymbol{\Sigma}}$ mean that $(\bar{\boldsymbol{\Theta}},\bar{\boldsymbol{\Sigma}},\bar{\boldsymbol{\Lambda}})$ is an optimal solution to (5). Therefore, we showed that any limit point of $\{(\boldsymbol{\Theta}^{i},\boldsymbol{\Sigma}^{i},\boldsymbol{\Lambda}^{i})\}$ is an optimal solution to (5). ∎

Without loss of generality, we may always assume that $E(X_{ij})=0$ for all $1\leq i\leq n,1\leq j\leq p$. By the condition that $\boldsymbol{\Sigma}^{0}$ is positive definite, we can always choose some very small $\epsilon>0$ such that $\epsilon$ is smaller than the minimal eigenvalue of $\boldsymbol{\Sigma}^{0}$. We introduce $\boldsymbol{\Delta}=\boldsymbol{\Sigma}-\boldsymbol{\Sigma}^{0}$, and then we can write (3) in terms of $\boldsymbol{\Delta}$ as follows,

Note that it is easy to see that $\hat{\boldsymbol{\Delta}}=\hat{\boldsymbol{\Sigma}}^{+}-\boldsymbol{\Sigma}^{0}$.

Note that $\hat{\boldsymbol{\Delta}}$ is also the optimal solution to the following convex optimization problem

Under the same probability event, $\|\hat{\boldsymbol{\Delta}}\|_{F}\leq 5\lambda{(s+p)}^{1/2}$ would always hold. Otherwise, the fact that $G(\boldsymbol{\Delta})>0$ for $\|\boldsymbol{\Delta}\|_{F}=5\lambda{(s+p)}^{1/2}$ should contradict with the convexity of $G(\cdot)$ and $G(\hat{\boldsymbol{\Delta}})\leq G(\boldsymbol{0})=0$. Therefore, we can obtain the following probability bound

Now we shall prove the probability bound under the exponential-tail condition. First it is easy to verify two simple inequalities that $1+u\leq\exp(u)\leq 1+u+\frac{1}{2}u^{2}\exp(|u|)$ and $v^{2}\exp(|v|)\leq\exp(v^{2}+1)$. The first inequality can be proved by using the Taylor expansion, and the second one can be easily derived using the obvious facts that $\exp(v^{2}+1)\geq\exp(2|v|)$ and $\exp(|v|)\geq v^{2}$.

Let $t_{0}=(\eta\frac{\log p}{n})^{1/2}$, $c_{0}=\frac{1}{2}eK_{1}\eta^{1/2}+\eta^{-1/2}(M+1)$ and $\varepsilon_{0}=c_{0}(\frac{\log p}{n})^{1/2}$. For any $M>0$, we can apply the Markov inequality to obtain that

where we apply $\exp(u)\leq 1+u+\frac{1}{2}u^{2}\exp(|u|)$ in the second inequality and $1+u\leq\exp(u)$ in the third inequality, and then use $v^{2}\exp(|v|)\leq\exp(v^{2}+1)$ in the fourth inequality. Moreover, the simple facts that $E[X_{ij}]=0$ ($1\leq i\leq n$) and $t_{0}^{2}=\eta\frac{\log p}{n}\leq\eta$ are also used.

Let $t_{1}=\frac{1}{2}\eta(\frac{\log p}{n})^{1/2}$ and $c_{1}=2K_{1}(\eta^{-1}+\frac{1}{4}\eta\sigma_{\max}^{2})\exp(\frac{1}{2}\eta\sigma_{\max})+2\eta^{-1}(M+2)$. Define $\varepsilon_{1}=c_{1}(\frac{\log p}{n})^{1/2}$. For any $M>0$, we first apply the Cauchy inequality to obtain that

where we use the simple inequality $\exp(|v|)\geq v^{2}$ in the third inequality. Then, combining this result with the Cauchy inequality again yields that

where we use the fact that $t_{1}=\frac{1}{2}\eta(\frac{\log p}{n})^{1/2}\leq\frac{1}{2}\eta<\eta$ in the first inequality, and then use $|\sigma^{0}_{jk}|\leq(\sigma^{0}_{jj}\sigma^{0}_{kk})^{1/2}\leq\sigma_{\max}$ in the third inequality. Now, we can apply the Markov inequality to obtain the following probability bound

where we apply $\exp(u)\leq 1+u+\frac{1}{2}u^{2}\exp(|u|)$ and $E[X_{ij}X_{ik}]=\sigma^{0}_{jk}$ for $i=1,2,\cdots,n$ in the second inequality, and we use $1+u\leq\exp(u)$ in the third inequality.

Recall that $\lambda=c_{0}\frac{\log p}{n}+c_{1}(\frac{\log p}{n})^{1/2}=\varepsilon^{2}_{0}+\varepsilon_{1}$ and

Therefore, we can complete the probability bound under the exponential-tail condition as follows

In the sequel we shall prove the probability bound under the polymonial-tail condition. First, we define $c_{2}=8(K_{2}+1)(M+1)$ and $\varepsilon_{2}=c_{2}(\frac{\log p}{n})^{1/2}$. Define $\delta_{n}=n^{1/4}(\log n)^{-1/2}$, $Y_{ij}=X_{ij}I_{\{|X_{ij}|\leq\delta_{n}\}}$ and $Z_{ij}=X_{ij}I_{\{|X_{ij}|>\delta_{n}\}}$. Then we have $X_{ij}=Y_{ij}+Z_{ij}$ and $E[X_{ij}]=E[Y_{ij}]+E[Z_{ij}]$. By construction, $|Y_{ij}|\leq\delta_{n}$ are bounded random variables, and $E[Z_{ij}]$ are bounded by $o(\varepsilon_{2})$ due to the fact that $|E[Z_{ij}]|\leq\delta_{n}^{-3}E[{|X_{ij}|^{4}}I_{\{|X_{ij}|>\delta_{n}\}}]\leq K_{2}\delta_{n}^{-3}=o(\varepsilon_{2}).$ Now we can apply the Bernstein’s inequality (Bernstein, 1946; Bennett, 1962) to obtain that

where the fact that $var(Y_{ij})\leq E[X^{2}_{ij}]\leq E[X^{2}_{ij}I_{\{|X_{ij}|\geq 1\}}]+E[X^{2}_{ij}I_{\{|X_{ij}|\leq 1\}}]\leq K_{2}+1$ is used in the second inequality. Besides, we can apply the Markov inequality to obtain that

Then, we can derive the following probability bound

Let $c_{3}=8(K_{2}+1)(M+2)$ and $\varepsilon_{3}=c_{3}(\frac{\log p}{n})^{1/2}$. Recall that $\delta_{n}=(\frac{n}{\log(n)})^{1/4}$, and define $R_{ijk}=X_{ij}X_{ik}I_{\{|X_{ij}|>\delta_{n}\textrm{~{}or~{}}|X_{ik}|>\delta_{n}\}}$. Then we have $X_{ij}X_{ik}=Y_{ij}Y_{ik}+R_{ijk}$ and $\sigma^{0}_{jk}=E[X_{ij}X_{ik}]=E[Y_{ij}Y_{ik}]+E[R_{ijk}]$. By construction, $|Y_{ij}Y_{ik}|\leq\delta_{n}^{2}$ are bounded random variables, and $E[R_{ijk}]$ is bounded by $o(\varepsilon_{3})$ due to the fact that

Again, we can apply the Bernstein’s inequality to obtain that

where the fact that $var(Y_{ij}Y_{ik})\leq E[X^{2}_{ij}X^{2}_{ik}]\leq(E[X^{4}_{ij}]E[X^{4}_{ik}])^{1/2}\leq K_{2}+1$ is used.

Recall that $\lambda=c_{2}\frac{\log p}{n}+c_{3}(\frac{\log p}{n})^{1/2}=\varepsilon^{2}_{2}+\varepsilon_{3}$. Therefore, we can prove the desired probability bound under the polynomial-tail condition as follows

References