0

0

Elastic Net; Entropy adjustment; Gaussian Graphical Models; Lasso; Maximum Likelihood Estimation

1 Introduction↩︎

Inverse of the covariance matrix (also known as precision matrix) has an important role in modern statistics and machine learning. Its accurate estimation is required in different research fields including finance [1], psychology [2], genetics [3], brain studies [4], etc. The estimated precision matrix is fundamental in discriminant analysis, forecasting and several other statistical methodologies [5].

Precision matrix is closely related to the partial correlations among variables. Suppose the true precision matrix is \(\Omega=[\omega_{ij}]_{1\le i,j\le p} \in {\mathbb{R}}^{p\times p}\), which is positive definite. The entries of partial correlation matrix \(P=[\rho_{ij}]_{1\le i,j\le p}\) can be written in terms of the precision matrix entries: \[\begin{align} \label{partialcorr} \rho_{ij} &=& - \dfrac{\omega_{ij}}{\sqrt{\omega_{ii} \omega_{jj}}} \end{align}\tag{1}\]

Therefore, the \({(i,j)}\) entry of the precision matrix is zero if and only if the partial correlation of the variables \(X^i\) and \(X^j\) is zero. Under the assumption of multivariate normality, the entry \(\omega_{ij}=0\) of the precision matrix indicates the conditional independence between the variables \(X^i\) and \(X^j\), given all the other variables [6]. In this way, the precision matrix is closely related to the Gaussian Graphical Models (GGM) which can represent the conditional independence of multivariate normally distributed variables in a convenient form as a graph (e.g., gene interaction networks). The GGM is an undirected graph defined as \({G=(N,E)}\). Here, the set of nodes \({N=\{1,...,p\}}\) represents the variables and the set of edges \(E\subseteq N \times N\) consists of the pair indexes \({(i,j)}\) corresponding to the ‘active’ entries \({\omega_{ij}\not = 0}\) for \(1\le i,j \le p\) [7]. The selection of GGM is an important problem in studying genetic interactions, human brain activity network, psychological networks, etc.

In recent decades, substantial research has focused on estimating the precision matrix under high dimensional settings. The most popular approach under high dimensional settings is applying Lasso (Least absolute shrinkage and selection operator) or \(\ell_1\) norm penalty on an objective loss function. This approach is proposed in the regression context by [8] and is widely used for estimating precision and covariance matrices. In particular, [9] and [10] independently proposed the \(\ell_1\) norm penalized log-likelihood maximization approach, also known as Graphical Lasso or GLasso estimator. Methods that employ \(\ell_1\) norm penalization include the Neighborhood Selection [11], Sparse Partial Correlation Estimation [12], constrained \(\ell_1\) norm minimization for inverse matrix estimation or CLIME [13], \(\ell_1\) norm penalized D-trace loss minimization [14], Sparse Column-Wise Inverse Operator or SCIO [15], among several others. Despite the popularity of \(\ell_1\) norm penalty, other penalties (both convex and non-convex) are also considered, including Adaptive \(\ell_1\) norm [16], [17], SCAD [16], Ridge or \(\ell_2\) norm [18]–[20], Elastic-Net [21], [22], among others. We provide more details on estimation techniques in the next section.

In this article, we study the estimation of precision matrix \(\Omega\) without assuming a specific structure or a sparsity pattern. In line with above-cited literature, our proposed method uses \(\ell_1\) norm penalized log-likelihood optimization problem. Although this penalty produces a sparse estimate, it can also lead to an increased uncertainty in capturing the true relationships between variables. This is particularly important when Lasso incorrectly shrinks some of the entries to zero. Moreover, Lasso penalization tends to shrink the eigenvalues of the estimated precision matrix. In general, smaller eigenvalues of a precision matrix correspond to higher uncertainty in the data structure. In order to diminish the induced uncertainty, we consider an additional log-determinant penalty on our optimization problem. This term is known as an appropriate measure of uncertainty in a context of multivariate normal distribution and is directly related to entropy [23], [24]. Our introduced method therefore induces sparsity through Lasso penalty term and adjusts uncertainty through the log-determinant term. The new augmented penalty is the convex combination of these two terms, motivated by the Elastic-Net penalization framework of [21] and [22].

We evaluate the performances of the proposed estimator through an extensive numerical study. In particular, for our simulation analyses we consider different patterns (i.e., models) of the true precision matrix \(\Omega\), including those used in the experiments of the previous studies. We measure the performance of the considered methods in terms of statistical losses and GGM prediction measures. Finally, for the proposed method, we establish the convergence rate in the Frobenius norm under standard asymptotic theory (i.e., we assume that the number of variables is fixed).

The manuscript is organized as follows. In Section 2, we describe the related work and our proposed methodology. In Section 3, we evaluate the statistical performance of the proposed methodology and compare it with that of other approaches in the literature. In Section 4, we provide real data applications: classification of prostate cancer patients using Linear Discriminant Analysis and computation of an optimal financial portfolio. We provide our conclusions in Section 5. We demonstrate technical details in Appendix A, simulation models in Appendix B and additional numerical results in Appendix C. Finally, Appendix D develops the analytical convergence rate of the proposed method.

2 Methodology↩︎

2.1 Background and existing research↩︎

We use the following mathematical notations throughout this paper. For any symmetric matrix \(\boldsymbol{A}=[a_{ij}]_{1\le i,j\le p} \in \mathbb{R}^{p\times p}\), we denote the Frobenius or \(\ell_2\) norm by \(||\boldsymbol{A}||_2=\sqrt{\sum\limits_{i=1}^{p}\sum\limits_{j=1}^{p}a_{ij}^2}\), the \(\ell_{\infty}\) norm by \(||\boldsymbol{A}||_{\infty}=\max\limits_{1\le i,j\le p}|a_{ij}|\), the matrix \(\ell_1\) norm by \(||\boldsymbol{A}||_{\ell_1}=\max\limits_{1\le j\le p}\sum\limits_{i=1}^p|a_{ij}|\), the elementwise \(\ell_1\) norm by \(||\boldsymbol{A}||_{1}=\sum\limits_{i=1}^{p}\sum\limits_{j=1}^{p}|a_{ij}|\), the spectral norm by \(||\boldsymbol{A}||_{\text{sp}}=\sup\limits_{||x||_2\le 1}||\boldsymbol{A}x||_2\). Here, we define the \(\ell_{2}\) norm by \(||\boldsymbol{a}||_2=\sqrt{\sum\limits_{j=1}^{p}a_j^2}\) for any \(p\)-dimensional vector \(\boldsymbol{a}=(a_1,...,a_p)^T\in \mathbb{R}^p\).

Without loss of generality, we assume that \(\boldsymbol{X}_{n\times p}\) is mean-centered, observed sample data matrix, where each row \(X_{i} =\left(X_{i1},...,X_{ip}\right)\) is a \(p\)-dimensional normal random vector, i.i.d. for \(i=1,...,n\) and has a covariance matrix \(\Sigma=\Omega^{-1} \in \mathbb{R}^{p\times p}\). We define the sample covariance matrix as \(S=\dfrac{1}{n}\sum\limits_{i=1}^{n}X_i^{T}X_i = \dfrac{1}{n}\boldsymbol{X}^T\boldsymbol{X}\).

The estimation of a precision matrix under high dimensional settings has been widely studied in prior scholarly work. Among the proposed methods, the well-explored estimator is based on minimizing penalized (or constrained) negative log-likelihood function of a multivariate normal distribution. In particular, Graphical Lasso or GLasso [9], [10], [25] is the most popular approach for estimating the precision matrix and the corresponding GGM. This estimator is defined as: \[\begin{align} \label{GLASSO}\widehat{\Omega}_{\text{GLasso}}&=&\arg\min_{\Omega}\;-\log\det(\Omega)+\text{trace}(\Omega S)+\gamma||\Omega||_{1}, \end{align}\tag{2}\] where the Lasso or \(\ell_1\)-norm penalty is used to encourage sparse solutions. Here, \(\gamma>0\) is the associated penalty (i.e., tuning or calibration) parameter that determines the degree of applied penalization on \(\hat{\Omega}\). In case of GLasso, this parameter adjusts the accuracy and the sparsity of precision matrix estimator. When \(\gamma = 0\), the problem (2 ) leads to the Maximum Likelihood Estimate of the covariance matrix \(\hat{\Sigma} = S\).

Extant literature has studied the convergence rates of the GLasso precision matrix estimate \(\widehat{\Omega}_{\text{GLasso}}\) under certain assumptions. In particular, [26] derived the convergence rate of the Glasso estimate in terms of Frobenius norm \(\ell_2\), whereas [27] established the convergence rate in terms of \(\ell_\infty\) and spectral norms. Moreover, [27] prove the model selection consistency, i.e., show the estimate \(\widehat{\Omega}_{\text{Lasso}}\) specifies correctly the sparsity pattern of \(\Omega\), with probability converging to one.

In the context of penalization-driven approaches, a natural alternative to the Lasso penalization is the Ridge (or \(\ell_2\)-norm) penalization. The most straightforward version of this estimator is defined as: \[\begin{align} \label{GRidge}\widehat{\Omega}_{\text{GRidge}}&=&\arg\min_{\Omega}\;-\log\det(\Omega)+\text{trace}(\Omega S)+\gamma||\Omega||_{2}^2. \end{align}\tag{3}\] This estimator was initially discussed by [28] and [26] and was further developed independently by [18] and [19]. [19] refer to it as ROPE (Ridge-type Operator for the Precision matrix Estimation). On the other hand, [18] suggested an alternative Ridge-type estimator, defined as: \[\begin{align} \label{GRidgeT}\widehat{\Omega}_{\text{T-GRidge}}&=&\arg\min_{\Omega}\;-\log\det(\Omega)+\text{trace}(\Omega S)+\dfrac{\gamma}{2}||\Omega - {T}||_{2}^2, \end{align}\tag{4}\] where \({T} \in \mathbb{R}^{p\times p}\) is a symmetric positive definite target matrix. Clearly, the penalty term of \((\ref{GRidge})\) is the special case of \((\ref{GRidgeT})\) with \({T}=\boldsymbol{0}\). Note that Ridge-type estimators are less complex than GLasso estimator (3 ) because they can be calculated by an explicit closed-form solutions. However, Ridge penalty does not induce exact sparsity on the estimated precision matrix, similar to Lasso penalty. In this way, these estimators are more attractive, if the sparsity of \(\Omega\) is not an important property for the precision matrix estimate. For the sake of simplicity, analogously to the Graphical Lasso, we use the term Graphical Ridge (or GRidge) for the Ridge-type estimator \((\ref{GRidge})\) with no target, and Targeted Graphical Ridge (or T-GRidge) for the Ridge-type estimator \((\ref{GRidgeT})\) with a specified target.

More recently, scholarly work has studied the convex combination of Lasso and Ridge penalties. In the regression context, this penalty is known as Elastic-Net [29]. [22] consider the following augmented estimation problem for estimating the precision matrix \(\Omega\): \[\begin{align} \label{ElNet}\widehat{\Omega}_{\text{GEN}}&=&\arg\min_{\Omega}\;-\log\det(\Omega)+\text{trace}(\Omega S)+\gamma\left(\alpha ||\Omega||_{1} + (1-\alpha) ||\Omega||_{2}^2\right), \end{align}\tag{5}\] where the parameter \(\alpha \in [0;1]\) controls the balance between Lasso and Ridge penalization. For \(\alpha = 1\), we have the Graphical Lasso approach (2 ), and for \(\alpha = 0\), we have the Rope approach (4 ). Consistent with the terminology in regression context, this approach is called Graphical Elastic-Net (or GEN). The Elastic-Net penalization often brings the benefits of both Lasso and Ridge penalties: it leads to sparse estimators and simultaneously encourages stability in a presence of highly correlated variables.

Similar to the targeted Ridge penalization (4 ), [21] suggested a targeted version for the Elastic-Net penalty: \[\begin{align} \label{ElNet-T}\widehat{\Omega}_{\text{T-GEN}}&=&\arg\min_{\Omega}\; -\log\det(\Omega)+\text{trace}(\Omega S) \nonumber\\ &+&\gamma\left(\alpha ||\Omega - T||_{1} + \dfrac{1-\alpha}{2} ||\Omega - T||_{2}^2\right). \end{align}\tag{6}\] Several target matrices are suggested in [21]. The penalty term of \((\ref{ElNet})\) is a special case of \((\ref{ElNet-T})\) with \({T}=\boldsymbol{0}\), and the penalty term of \((\ref{GRidgeT})\) is a special case of \((\ref{ElNet-T})\) with \(\alpha=0\).

Note that the popularity of \(\ell_1\) and \(\ell_2\) norm penalties is due to their convexity. The literature also considers applying nonconcave penalties on the log-likelihood function, such as Smoothly Clipped Absolute Deviation or SCAD [16], \(\ell_0\) norm penalty [30], [31], \(\ell_q\) norm penalty with \(0<q<1\) [32]. For a further review on approaches for estimating high dimensional precision matrices, we refer to [33].

2.2 Proposed methodology↩︎

Following the Elastic-Net-type penalization (5 ), we propose a new approach for estimating the precision matrix. Consider the following optimization problem: \[\begin{align} \label{EAGL}\widehat{\Omega}_{\text{EAGL}}&=&\arg\min_{\Omega}\; -\log\det(\Omega)+\text{trace}(\Omega S) \nonumber\\ &+&\gamma\left( \alpha ||\Omega||_{1} + (1-\alpha) \log\det(\Omega^{-1}) \right), \end{align}\tag{7}\] where \(\alpha\in (0;1)\) is a combination parameter analogous to its use in the Elastic-Net penalization (5 ). The proposed approach is an augmented version of the GLasso method (2 ) with an additional \(H_C(\boldsymbol{X})=\log\det(\Omega^{-1})=\log(\det(\Sigma))\) penalty term. Here, we still use the Lasso penalty, because our goal is to obtain a sparse estimate. The motivation behind the \(\log\det\) term is twofold.

First, we can show that this term represents the entropy of multivariate Gaussian distribution, up to an additive constant. We show this in Appendix A. In the context of machine learning and statistics, entropy is viewed as a measure of uncertainty or randomness within a dataset. High entropy indicates more randomness, while low entropy suggests a more organized data and thus, better measurability of the relationships among the variables. The objective of the proposed technique is therefore to impose sparsity on the precision matrix estimate while adjusting the overall uncertainty, which leads to a more efficient estimate.

Second, the term \(H_C(\boldsymbol{X})\) produces an additional calibration for the eigenvalues of estimated precision matrix. Note that \(H_C(\boldsymbol{X})=\log\det(\Omega^{-1})=-\sum_{i=1}^{p}log(\lambda_{i})\), where \(\lambda_{i}\) is the \(i\)-th eigenvalue of the precision matrix \(\Omega\). The additional \(\log\det\) penalty forces the precision matrix eigenvalues to shrink less which occurs due to the Lasso penalization (such shrinkage pattern is in line with the findings by [25] and [34]). As a result, the eigenvalues of the estimated precision matrix \(\widehat{\Omega}_{\text{EAGL}}\) become higher and less shrunk than those of \(\widehat{\Omega}_{\text{GLasso}}\). In Figure 1 (a), we provide a simple example, which shows the distribution of eigenvalues of \(\widehat{\Omega}_{\text{GLasso}}\) and \(\widehat{\Omega}_{\text{EAGL}}\) for different values of the penalization parameter \(\gamma\) from a sequence varying from 0.1 to 2. For the data simulation, we use the settings described in Section 3 and we assume that the true precision matrix has a structure given in Model 1 for \(p=200\) and \(n=100\). We can clearly see that the eigenvalues of \(\widehat{\Omega}_{\text{EAGL}}\) shrink slower than those of \(\widehat{\Omega}_{\text{GLasso}}\). Therefore, we can expect a more accurate eigenspectrum of the estimated precision matrix with the additional penalty term.

Note that large precision matrix eigenvalues indicate strong concentration in the principal axes and thus, lower uncertainty. This behavior aligns with the intuition that as the matrix becomes sparser due to the Lasso penalization, it reflects to fewer restrictions on the relationships between the variables (i.e., many variables become conditionally independent given the others). This will lead to a greater variability in the joint distribution (i.e., it will be more spread out) and, consequently, to a higher entropy.

The use of \(\log\det\) penalty is not new in the literature for other purposes. For example, a similar \(\log\det\) penalty is used by [26] for estimating sparse covariance matrix. However, this penalty was used as a logarithmic barrier to ensure that the obtained estimator is positive definite. To the best of our knowledge, this term has not been used in the context of entropy adjustment.

We call the proposed approach Entropy Adjusted Graphical Lasso or EAGL. Note that when \(\alpha=1\), the problem (7 ) reduces to (2 ). Interestingly, when \(\alpha = 0\), the problem (7 ) leads to the covariance matrix estimator \(\hat{\Sigma} = \dfrac{1}{1 + \gamma} S\), which is a special case of the James-Stein estimator [35] and has a better performance than the traditional estimator \(\hat{\Sigma} = S\) in terms of the MSE [36]. This estimator can also be seen as the Ledoit-Wolf shrinkage covariance estimator [37] towards a zero target.

Remark: Following the targeted version of the Elastic-Net estimator (5 ), we can also employ an additional target matrix \(T\) in the Lasso penalty. However, in this paper we focus on the analysis of EAGL estimator without using target matrices.

3 Simulation analysis↩︎

In this section, we evaluate the numerical performance of the proposed method using simulated data. The data generation process is based on different scenarios (i.e., models) for the true precision matrix \(\Omega\). Specifically, we consider several popular methods for estimating the precision matrix. We compare our proposed estimator EAGL with Graphical Lasso, Graphical Elastic-Net (with and without a target), Graphical Ridge (with and without a target) and Graphical SCAD. Our study is implemented in \(\mathbf{R}\) using the following packages: glassoFast for EAGL and GLasso, GLassoElnetFast for Graphical Elastic-Net estimators, rope for Graphical Ridge estimators, ggmncv for Graphical SCAD. All packages are available at http://cran.r-project.org/web/packages. Following [19], we use scalar matrix \(T = \nu \boldsymbol{I}\) for the target matrix, which is used in all targeted estimators (here, \(\nu = p / \text{trace}(S)\) and \(\boldsymbol{I}\in \mathbb{R}^{p\times p}\) is identity matrix). As indicated earlier, Graphical Ridge approaches do not provide sparse estimators. Therefore, we perform an additional sparsification on these estimates using sparsify() function of rags2ridges package, with ‘localFDR’ thresholding argument [38].

3.1 Performance evaluation↩︎

We evaluate the performance of a precision matrix estimator based on several popular statistical distances and errors. We consider the relative entropy losses, such as the Kullback-Leibler loss and the Reverse Kullback-Leibler loss [34], [39], [40], the Relative Trace Error (RTE) and matrix losses, such as the Frobenius \(\ell_2\) norm, the spectral \(\ell_{sp}\) norm and the matrix \(\ell_1\) norm [13], [14], [18]. These measures are defined as: \[\begin{align} \tag{8} \text{KLL}(\widehat{\Omega},\Omega)&=&\text{trace}(\Omega^{-1} \widehat{\Omega})-\log\det(\Omega^{-1} \widehat{\Omega})-p,\\ \tag{9} \text{RKLL}(\widehat{\Omega},\Omega)&=&\text{trace}(\Omega \widehat{\Omega}^{-1})-\log\det(\Omega \widehat{\Omega}^{-1})-p,\\ \tag{10} \text{RTE}(\widehat{\Omega},\Omega)&=&\left|1-\dfrac{\text{trace}(\widehat{\Omega})}{\text{trace}(\Omega)}\right|\\ \tag{11} \ell_2(\widehat{\Omega},\Omega)&=&||\widehat{\Omega}-\Omega||_2,\\ \tag{12} \ell_{\text{sp}}(\widehat{\Omega},\Omega)&=&||\widehat{\Omega}-\Omega||_{\text{sp}},\\ \tag{13} \ell_1(\widehat{\Omega},\Omega)&=&||\widehat{\Omega}-\Omega||_{\ell_1}. \end{align}\]

Furthermore, we evaluate the graphical modeling (i.e., sparsity pattern prediction of the GGM) based on overall summary metrics: Matthews Correlation Coefficient [41] and Balanced Accuracy [42]. These metrics are respectively defined as: \[\begin{align} \tag{14} \text{MCC}&=&\dfrac{\text{TP}\times \text{TN}-\text{FP}\times \text{FN}}{\sqrt{(\text{TP}+\text{FP})(\text{TP}+\text{FN})(\text{TN}+\text{FP})(\text{TN}+\text{FN})}},\\ \tag{15}\text{BA}&=&\dfrac{1}{2} \left( \dfrac{TP}{TP + FN } + \dfrac{TN}{TN + FP} \right). \end{align}\]

Here, TP is the number of correctly selected non-zero entries (i.e., true positives), TN is the number of correctly selected zero entries (i.e., true negatives), FP is the number of incorrectly selected non-zero entries (i.e., false positives), and FN is the number of incorrectly selected zero entries (i.e., false negatives). In case of MCC, the coefficient of 1 and \(-1\) indicate the perfect and the worst classification, respectively. In case of the Balanced Accuracy, the coefficient 1 indicates the perfect selection, whereas the coefficient \(0\) indicates the worst selection. For convenience, we further rescale the MCC to the range of the Balanced Accuracy. \[\begin{align} \label{uMCC} \text{uMCC}&=&\dfrac{\text{MCC}+1}{2}. \end{align}\tag{16}\] Both metrics are popular techniques in statistics and machine learning as standard accuracy metrics for evaluating binary classifications [43].

3.2 Simulation settings↩︎

We generate multivariate normal random samples with zero mean and covariance matrix \(\Sigma=\Omega^{-1}\). The following models are considered for the true precision matrix \(\Omega = [\omega_{ij}]_{1\le i, j \le p}\):

• Model 1. Band matrix, with \(\omega_{ii}=1\), \(\omega_{i,i-1}=\omega_{i-1,i}=0.45\) and other values are 0 [10], [25]

• Model 2. Band matrix, with \(\omega_{ii}=1\), \(\omega_{i,i-1}=\omega_{i-1,i}=0.5\), \(\omega_{i,i-2}=\omega_{i-2,i}=0.35\) and other values are 0 [19], [34].

• Model 3. A block-diagonal matrix, with four equally sized blocks along the diagonal. Each block is defined as \(\omega_{ij}=0.6^{|i-j|}\) [13], [16]

• Model 4. A random positive definite matrix, with approximately \(50\%\) of non-zero entries. This matrix is generated using Matlab command sprandsym with a parameter \(0.5\) [34].

• Model 5. Erdős-Rényi random graph. We define \(\Theta = [\theta_{ij}]_{1\le i, j \le p}\), with \(\theta_{ij} = u_{ij} \delta_{ij}\), where \(u_{ij} = u_{ji}\) is a uniform random variable on an interval \([0.4, 0.8]\) and \(\delta_{ij}=\delta_{ji}\) is a Bernoulli random variable with a success probability \(0.05\). Next, we take \(\Omega = \Theta + (|\lambda_{\text{min}}(\Theta) + 0.05|) \boldsymbol{I}\), where \(\lambda_{\text{min}}(\Theta)\) is the smallest eigenvalue of \(\Theta\) [21].

• Model 6. Scale-free graph, based on the Barabási–Albert algorithm. The resulting graph has \(p\) edges. This matrix is generated using R package huge [21], [22].

• Model 7. Hub graph with 10 hubs, each hub containing \(p/20\) nodes. The resulting graph has \(p-10\) edges. This matrix is generated using R package huge with parameters \(\nu=0.3\) and \(u=0.1\) [21], [22].

The generated matrices are further standardized to have unit diagonal. We provide the graphical models corresponding to each considered precision matrix structure in Appendix B. Finally, we set the sample size \(n=100\) and the number of variables \(p=200\). The number of replications is 100.

3.3 Penalty parameter selection↩︎

Selecting the penalty parameter is a crucial step for any penalized optimization problem. In our numerical analysis, we select the penalty parameter \(\gamma\) using 5-fold Cross-Validation or CV [44] approach. This is a commonly used technique in precision matrix estimation methods and is consistent with studies of [21], [19], etc. Following [22] and [21], we set \(\alpha = 0.5\) for the estimators GEN, T-GEN and the proposed EAGL.

Another popular technique for estimating the penalty parameter of the precision matrix estimation methods is Bayesian Information Criterion or BIC [10]. In this paper, we put CV above BIC, because CV leads to lower statistical losses. However, according to [45], CV may provide many false positives compared to BIC. Following [22], we conduct additional analyses for demonstrating the performance of our proposed methodology where the penalty parameter \(\gamma\) is selected using BIC approach. We provide these results in Appendix C. These findings provide further insights into the comparison of the two selection techniques.

3.4 Computational time↩︎

Table 1 shows the average computational time required for each method. The reported time includes the required time to select the penalty parameter \(\gamma\) using CV and the required time to run the method. We assume that the true precision matrix has a structure given in Model 1 for \(p=200\) and \(n=100\). The analyses are performed on Intel Core i7-1280P CPU, 1.80 GHz. The computational time of EAGL is relatively low, compared to the other methods. Note that, in general, Ridge-type approaches are not complex and require little computation (usually less than a second). However, the time provided in Table 1 is high due to the sparsification procedure.

3.5 Discussion of results↩︎

Tables 5-11 report the averages of the measures over 100 replications. Corresponding standard deviations are provided in the parentheses. First, we observe that EAGL outperforms GLasso in terms of RKLL, RTE, \(\ell_2\), \(\ell_{\text{sp}}\) and \(\ell_1\) norms for all models, in terms of KLL for models 1, 3, 6 and 7, and in terms of uMCC for models 1, 2, 3, and 7. On the other hand, GLasso method performs better than EAGL in terms of KLL for models 2, 4 and 5 and in terms of uMCC for models 4, 5 and 6. Table 2 provides a summary for the comparison between EAGL and GLasso estimators, indicating scenarios, where EAGL outperforms GLasso.

Comparing our proposed method with other methods, we see that in general EAGL provides better results, especially in terms of the statistical losses. However, SCAD provides the best overall results in terms of the KLL for models 1, 2, 3, 4, in terms of the \(\ell_1\) norm for model 7 and in terms of uMCC for model 5. We observe that T-GEN estimator provides the best overall results in terms of KLL for models 5 and 6, in terms of \(\ell_2\) norm in terms for model 6, in terms of the \(\ell_1\) norm for model 6 and in terms of BA for model 4 and 5. On the other hand, GEN estimator provides the best overall results in terms of BA for model 6. T-GRidge estimator provides the best overall results in terms of the \(\ell_1\) norm for model 3, in terms of BA for model 1 and in terms of uMCC for model 7. Finally, GRidge estimator provides the best overall results in terms of the \(\ell_1\) norm for models 1 and 2, in terms of uMCC for models 1, 2, 3 and 6.

.

In sum, the proposed EAGL estimator in general provides favorable performance than GLASSO, GEN, T-GEN, GRidge, T-GRidge and SCAD methods for most of the models in terms of most statistical losses and GGM prediction measures. Moreover, EAGL shows a suitable trade-off between the matrix prediction and sparsity pattern identification, i.e., the outperformance in terms of one criterion does not diminish the other one. This suggests that entropy adjustment enhances both the statistical accuracy of the estimated precision matrix and the performance of GGM selection.

In addition, we conduct the same study where the penalty parameter is selected using BIC (see Appendix C). The results support our discussion above and the comparison of the considered methods remains comparatively the same.

4 Real data application↩︎

In this section, we conduct an empirical analysis of the proposed method through two real-data applications. The first application aims at predicting prostate cancer patients and the second one aims at selecting a large financial portfolio.

4.1 Prostate cancer study↩︎

In this application, we focus on the problem of predicting prostate cancer patients using Linear Discriminant Analysis (LDA) with different precision matrix estimates. We use a dataset analysed by [46], which contains measurements of the gene expression levels of 6033 genes for 102 specimens. Out of these specimens, 52 are taken from prostate cancer patients, and 50 are taken from healthy patients. The dataset singh2002 is available in R package sda.

In order to evaluate the performance of the proposed EAGL method, we first divide the data into a training set and a testing set with sizes 60 and 42, respectively. Second, for the training set we apply two sample t-test between the two groups (cancer and healthy) in order to select the most significant 100 genes with the smallest p-values. Based on the selected genes, we obtain the estimated precision matrix \(\hat{\Omega}\) using the methods considered in section 3. The penalty parameters are estimated using the 5-fold CV technique. The estimated precision matrix is then used in the following LDA procedure.

Let’s assume \(\boldsymbol{\mu}_1\) and \(\boldsymbol{\mu}_2\) are the population means of the gene expression levels for cancer and healthy patients, respectively, and \(\Omega\) is the corresponding population precision matrix. We use Mahalanobis distance \(\boldsymbol{a}^T(\boldsymbol{x}-\boldsymbol{\mu})\), where \(\boldsymbol{a} = \Omega (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_2)\) and \(\boldsymbol{\mu} = \dfrac{\boldsymbol{\mu}_1+\boldsymbol{\mu}_2}{2}\). We estimate the means \(\boldsymbol{\mu}_1\) and \(\boldsymbol{\mu}_2\) using within group averages \(\bar{\boldsymbol{x}}_1\) and \(\bar{\boldsymbol{x}}_2\), respectively, calculated with the training date. Next, we assign the testing observation \(\boldsymbol{x}\) to the cancer group if \(\hat{\boldsymbol{a}}^T(\boldsymbol{x}-\hat{\boldsymbol{\mu}})> 0\), where \(\hat{\boldsymbol{a}} = \hat{\Omega} (\bar{\boldsymbol{x}}_1 - \bar{\boldsymbol{x}}_2)\) and \(\hat{\boldsymbol{\mu}} = \dfrac{\bar{\boldsymbol{x}}_1+\bar{\boldsymbol{x}}_2}{2}\). We repeat this process 100 times and calculate the missclassification rate for each precision matrix estimator. In addition, we calculate the uMCC and BA scores, where we assume that TP and TN are the number of correctly predicted cancer and healthy patients, respectively, and FP and FN are the number of incorrectly classified cancer and healthy patients, respectively.

Table 3 reports the average prediction measures for each method. We observe that the proposed EAGL shows lower missclassification rate and higher uMCC and BA than the GLasso, GEN, GRidge, T-Gridg and SCAD estimators. We observe that EAGL and the targeted T-GEN have roughly the same classification measures. However, we note that in contrast to the T-GEN approach, the proposed EAGL approach does not require selecting a target matrix and is computationally more efficient.

4.2 S&P 500 portfolio optimization↩︎

In our second application, we focus on developing an optimal stock portfolio with minimum risk [22], [47], [48]. Note that estimating the precision matrix has a crucial role in computing optimal portfolio weights [49]. In the mean-variance optimization, the risk of a \(p\)-dimensional (global) portfolio \(w=\left(w_1, ..., w_p\right)\) is measured by the standard deviation of its returns, given as \(\sqrt{w^T\Sigma w}\) [50], where \(\Sigma\) is the global covariance matrix. In the absence of short-sale constraints, the estimated weights \(\hat{w}_{MVP}\) of the single-period minimum variance portfolio are defined as: \[\begin{align} \hat{w}_{MVP} &=& \arg \min_w w^T\hat{\Sigma} w \\ & &\text{subject to } w^T\boldsymbol{1} = 1 , \nonumber \end{align}\] where 1 is a \(p\)-dimensional vector of ones. [51] demonstrate that the minimum variance optimization problem has an explicit solution, defined as: \[\begin{align} \hat{w}_{MVP} &=& \dfrac{\hat{\Sigma}^{-1}\boldsymbol{1}}{\boldsymbol{1}^T\hat{\Sigma}^{-1}\boldsymbol{1}}. \end{align}\] This means that the minimum variance portfolio depends on the precision matrix estimate \(\hat{\Omega} = \hat{\Sigma}^{-1}\). Therefore, we expect that an accurate estimate of the precision matrix will lead to a better portfolio.

We use historical monthly returns of \(p=267\) stock constituents of S&P 500 index for a total of \(n = 364\) months, covering the period from January 1991 to April 2021. The constituents come from 11 major market sectors. The dataset is available at the R package probstats4econ. Following [37], we assume that stock returns are independent and identically distributed. In order to evaluate the performance of selected portfolio, we perform the ‘rolling-horizon’ procedure of [51]. First, we select a rolling window with length \(m=150\) months, which is used to estimate the precision matrix \(\Omega\) and the corresponding portfolio weights. Next, we compute the out-of-sample portfolio returns of the next period. This procedure is repeated \(n-m\) times by including the returns of next period and dropping the earliest one (thus, keeping the rolling window with the same size). Finally, we evaluate the performance of the portfolio by computing the out-of-sample mean (i.e., reward), the standard deviation (i.e., risk) and the Sharpe Ratio, based on the obtained \(n-m\) out-of-sample returns \(\hat{R}_t\): \[\begin{align} \hat{\mu} &=& \frac{1}{n-m}\sum_{t = m}^{n-1}\hat{R}_{t+1}, \\ \hat{\sigma}^2 & =& \frac{1}{n-m-1}\sum_{t = m}^{n-1}(\hat{R}_{t+1}-\hat{\mu})^2, \\ \widehat{\text{SR}} &=& \frac{\hat{\mu}}{\hat{\sigma}} . \end{align}\] The Sharpe Ratio, also known as the reward-to-risk ratio, measures the performance of a portfolio relative to its associated risk. A higher Sharpe Ratio indicates better investment performance in terms of risk-adjusted returns.

In our analysis, in addition to the selected methods, we also evaluate the performance of two alternatives, commonly used in the literature. The first one is the inverse of Ledoit-Wolf covariance estimator \(\hat{\Omega}_{\text{LW}} = \hat{\Sigma}_{\text{LW}}^{-1}\) [37], which is a popular approach for selecting an optimal financial portfolio. The second one is the naive approach \(\hat{\Omega}_{\text{Naive}} = \boldsymbol{I}\), which leads to allocating \(1/p\) weight of wealth to each of the \(p\) assets available for investment [52].

Table 4 presents the out-of-sample performances for portfolios obtained using different precision matrix estimates. The penalty parameter \(\gamma\) is selected using 5-fold CV technique. First, our results show that the proposed EAGL method provides the lowest out-of-sample portfolio risk. EAGL, GLASSO and LW estimators lead to the highest out-of-sample mean (i.e., reward) compared to all the other methods. Finally, our proposed approach demonstrate the highest Sharpe Ratio. This means that the portfolios obtained using the proposed precision matrix estimate are the most efficient in terms of the risk-adjusted return. The results show that the naive approach provides the lowest out-of-sample mean and the highest risk. Finally, despite its popularity, the LW estimator leads to relatively high risk.

In addition to the out-of-sample performance above, we also examine the interaction structure between the companies. For demonstrative purposes, we apply our proposed EAGL approach on the whole data. We observe that the partial correlation matrix (not shown), corresponding to the precision matrix, is roughly block-diagonal, where each block represents a market sectors. The interactions (i.e., partial correlations) are stronger among the companies within each sector. This indicates that there is a strong community structure, with communities made up of companies from the same sector. Figure 2 demonstrates the Gaussian Graphical Model illustrating the interactions among the companies. The GGM reveals groups of nodes (i.e., companies) such that within group interactions are stronger and more frequent than between group interactions.

5 Conclusions↩︎

In this article, we introduce a novel method for estimating sparse precision matrices in high dimensional settings. The proposed approach is an augmented version of the popular Graphical Lasso method, which incorporates Lasso penalization framework with an additional log-determinant penalty. This combined penalty allows for entropy adjustment of the multivariate Gaussian distribution, thereby reducing the uncertainty.

Through extensive numerical analyses, using both simulated and real datasets, we demonstrate that the proposed entropy adjustment helps achieve better performance of the estimated precision matrix without increasing the computational costs. We evaluate our method using different loss functions and prediction performance measures. Our method provides lower statistical losses and better model selection metrics for Gaussian Graphical Models compared to other established approaches in the literature. We calibrate the penalty parameter using 5-fold Cross-Validation and BIC, ensuring practical usability across different datasets. Our proposed approach does not rely on the selection of a target matrix, though it can be extended to include a target matrix in the Lasso norm. Finally, we establish the convergence rate of the proposed EAGL estimator in the Frobenius norm under standard asymptotic conditions.

For researchers, we suggest incorporating the proposed EAGL in situations, where uncertainty reduction and sparsity are of high importance. This recommendation is particularly pertinent for scholars and practitioners who employ the precision matrix in finance, discriminant analysis and network modeling, because the improved accuracy of the estimated precision matrix can significantly improve their analyses.

While our methodology has notable advantages, it is important to acknowledge certain limitations. First, in our study, we choose \(\alpha=0.5\), i.e., halfway between Lasso and \(\log\det\) penalties. However, the optimal value for \(\alpha\) may differ from this choice. Second, \(\alpha\) should be sufficiently far from zero, because EAGL estimate does not exist for \(\alpha=0\), when \(p\) is close to or exceeds \(n\). Finally, we may encounter some convergence problems for EAGL for a relatively large norm of the true covariance matrix.

Our results pave the way for future research. First, future research should further investigate the role of parameter \(\alpha\) on numerical and theoretical performance of our methodology. Second, it might be interesting to explore the relationship between the parameters \(\gamma\) and \(\alpha\) in greater detail. Third, a promising research direction involves applying the proposed entropy adjustment to other precision matrix estimators, such as the Entropy Adjusted Graphical Ridge. Finally, studying the proposed estimator from a Bayesian perspective could provide valuable insights and deepen our understanding of its theoretical properties.

Despite the limitations, we hope our paper contributes valuable insights into precision matrix estimation and simulates future research in this area.

Supplementary Materials↩︎

6 Technical details↩︎

Entropy is a crucial scientific concept closely related to uncertainty and randomness. It originates from the field of information theory and is introduced by [53]. It is borrowed in statistics for understanding the randomness in the data. Suppose \(\boldsymbol{X}\) follows a multivariate Gaussian distribution with mean \(\mu\) and covariance matrix \(\Sigma = \Omega^{-1}\), i.e., \(\boldsymbol{X}\sim N({\mu}, \Omega^{-1})\). The entropy of a continuous random variable is defined as: \[H(\boldsymbol{X}) = - \int p(\boldsymbol{X}) \log p(\boldsymbol{X})\,dx,\] where \(p(\boldsymbol{X})\) probability density function and is defined for the Gaussian distribution as: \[p(\boldsymbol{X}) = (2\pi)^{-p/2} \det(\Omega)^{1/2} \exp\left( -\dfrac{1}{2}(\boldsymbol{X}-\mu)^T\Omega(\boldsymbol{X}-\mu)\right)\] Therefore, we have: \[\begin{align} H(\boldsymbol{X}) & = & - \boldsymbol{E}[log (p(\boldsymbol{X}))] \\ & = & -\boldsymbol{E}\left[\log\left[(2\pi)^{-p/2} \det(\Omega)^{1/2} \exp\left( -\dfrac{1}{2}(\boldsymbol{X}-\mu)^T\Omega(\boldsymbol{X}-\mu)\right)\right]\right] \\ & = & \dfrac{p}{2}\log(2\pi) - \dfrac{1}{2}\log\det(\Omega) + \dfrac{1}{2}\boldsymbol{E}\left[(\boldsymbol{X}-\mu)^T\Omega(\boldsymbol{X}-\mu)\right]\\ & = & \dfrac{p}{2}(1+\log(2\pi)) - \dfrac{1}{2}\log\det(\Omega)\\ & = & \dfrac{p}{2}(1+\log(2\pi)) + \dfrac{1}{2}\log\det(\Omega^{-1})\\ & = & \dfrac{p}{2}(1+\log(2\pi)) + \dfrac{1}{2}H_C(\boldsymbol{X})\\ & = & C + \dfrac{1}{2}H_C(\boldsymbol{X}), \end{align}\] where \(C\) is a constant. Here, we used the fact that \(\boldsymbol{E}\left[(\boldsymbol{X}-\mu)^T\Omega(\boldsymbol{X}-\mu)\right]=p\).

7 Simulation models↩︎

8 Additional simulation results↩︎

In this section, we present additional numerical results based on the same simulation settings described in Section 3, with the penalty parameter \(\gamma\) now selected using the BIC criterion. For a comprehensive comparison, we provide the results for both 5-fold CV and BIC. The results indicate that, in the majority of cases, the precision matrices estimated using CV tend to have lower statistical losses compared to those estimated with BIC. This suggests that CV may offer a more effective approach for obtaining lower statistical losses in these scenarios.

When examining Gaussian Graphical Model (GGM) selection, the results are more mixed. For certain models, such as models 4 and 5, CV leads to higher values for both uMCC and BA, compared to BIC, suggesting better sparsity pattern selection performance. However, in other cases, such as models 2 and 6, CV results in higher BA but lower uMCC than BIC. Therefore, we do not observe a consistent advantage of BIC approach in terms of sparsity prediction measures.

9 Asymptotic analysis↩︎

In this section, we analyze the convergence rate of the proposed estimator \(\widehat{\Omega}_{\text{EAGL}}\). We make the following assumptions to guarantee the existence of the true precision matrix \(\Omega\):

• : \(\lambda_{\min}(\Omega)\ge \underline{\lambda}>0\),

• : \(\lambda_{\max}(\Omega)\le \bar{\lambda}\),

where \(\lambda_{\min}(\Omega)\) and \(\lambda_{\max}(\Omega)\) are the smallest and largest eigenvalues of \(\Omega\), respectively, and \(\bar{\lambda}\) and \(\underline{\lambda}\) are some positive values. We define the support set \(S=\{(i,j):\Omega_{ij}\not=0\}\) and assume that \(\text{card}(S)\le s\).

Under the assumptions \(A1\) and \(A2\), if \(\gamma \asymp \sqrt{\dfrac{\log p}{n}}\), then \[\label{convergence} ||\widehat{\Omega}_{\text{EAGL}}-\Omega||_2=O_P\left(\sqrt{\dfrac{(p+s)\log p}{n}}\right),\tag{17}\] under the standard asymptotics (i.e., assuming that \(p\) remains fixed, while \(n\) converges to infinity).

The proof of (17 ) is motivated by [54]. First, we rewrite the optimization problem of our proposed methodology as: \[\begin{align} \label{EAGL1}\widehat{\Omega}_{\text{EAGL}}&=&\arg\min_{\Omega}\; -\left(1+(1-\alpha)\gamma\right) \log\det(\Omega)+\text{trace}(\Omega S) + \gamma\alpha ||\Omega||_{1}. \end{align}\tag{18}\]

Consider the following function: \[\begin{align}Q(\Theta)& = & -\left(1+(1-\alpha)\gamma\right) \log\det(\Theta) + \text{trace}(\Theta S) + \gamma\alpha||\Theta||_1- \\ & - & \left( - \left(1+(1-\alpha)\gamma\right)\log\det(\Omega) + \text{trace}(\Omega S)+ \gamma\alpha||\Omega||_1 \right) \\ & = & -\left(1+(1-\alpha)\gamma\right) \left(\log\det(\Theta)-\log\det(\Omega)\right) + \text{trace}(\Theta-\Omega)S+\gamma \alpha \left(||\Theta||_1-||\Omega||_1\right), \end{align}\] where \(\Omega\) is the true precision matrix. Note that the estimator \(\widehat{\Omega}_{\text{EAGL}}\) minimizes the function \(Q(\Theta)\). This means that the difference matrix \(\widehat{\Delta}=\widehat{\Omega}_{\text{EAGL}} -\Omega\) minimizes the function \(G(\Delta)=Q(\Omega+\Delta)\). Moreover, \(G(\widehat{\Delta})\le G(0) = Q(\Omega) = 0\) and \(G(\Delta)\) is a convex function.

Define the following set: \[\begin{align} \Phi_n(M)=\{\Delta: \Delta=\Delta^T, ||\Delta||_2=M r_n\}, \end{align}\] where \(r_n=\sqrt{\dfrac{(p+s)\log p}{n}}\to 0\). If we show that \(\inf\{G(\Delta): \Delta\in \Phi_n(M)\}>0\), then the minimizer \(\widehat{\Delta}\) must be inside the sphere defined by \(\Phi_n(M)\), i.e., \(||\widehat{\Delta}||_2\le Mr_n\). \[\begin{align} \label{functionG} G(\Delta) & = & Q(\Omega + \Delta) = -(1+(1-\alpha)\gamma) \left( \log\det(\Omega + \Delta)-\log\det(\Omega) \right) \\ \nonumber & + & \text{trace}(\Omega + \Delta -\Omega)S + \gamma \alpha \left(||\Omega + \Delta||_1-||\Omega||_1\right) \\\nonumber & = & -\left(1+(1-\alpha)\gamma\right) \left(\log\det(\Omega + \Delta)-\log\det(\Omega)\right) + \text{trace} (S - \Sigma) \Delta\\ & + & \text{trace} \Sigma \Delta + \gamma \alpha \left(||\Omega + \Delta||_1-||\Omega||_1\right) \nonumber \end{align}\tag{19}\] From the Taylor expansion of the function \(g(t)=\log\det(\Omega + t\Delta)\), we get that: \[\begin{gather} \log\det(\Omega+\Delta)-\log\det(\Omega)=\text{trace}(\Omega^{-1}\Delta)-\\ \text{vec}(\Delta)^T\left[\int\limits_{0}^{1}(1-\tau)(\Omega+\tau\Delta)^{-1}\otimes(\Omega+\tau\Delta)^{-1}d\tau\right]\text{vec}(\Delta), \end{gather}\] where \(\otimes\) is the Kronecker product, and \(\text{vec}(\Delta)\) is a vectorization of matrix \(\Delta\). For \(\Delta\in \Phi_n(M)\), we have: \[\begin{gather} \text{vec}(\Delta)^T\left[\int\limits_{0}^{1}(1-\tau)(\Omega+\tau\Delta)^{-1}\otimes(\Omega+\tau\Delta)^{-1}d\tau\right]\text{vec}(\Delta) \ge \\ \lambda_{\min}\left(\int\limits_{0}^{1}(1-\tau)(\Omega+\tau\Delta)^{-1}\otimes(\Omega+\tau\Delta)^{-1}d\tau\right)||\Delta||_2^2 \ge \\ \int\limits_{0}^{1}(1-\tau)\lambda_{\min}^2(\Omega+\tau\Delta)^{-1}d\tau \;||\Delta||_2^2 \;\;\ge \;\; \frac{1}{2}\min_{0\le \tau\le 1}\lambda_{\min}^2(\Omega+\tau\Delta)^{-1} ||\Delta||_2^2. \end{gather}\] Furthermore: \[\begin{align} \min_{0\le \tau\le 1}\lambda_{\min}^2(\Omega + \tau \Delta)^{-1} \ge \lambda_{\max}^{-2}(\Omega + \Delta)\ge(||\Omega||_{\text{sp}}+||\Delta||_{\text{sp}})^{-2}\ge \dfrac{1}{\bar{\lambda}+||\Delta||_{\text{sp}}^{2}} \ge \dfrac{1}{2} \bar{\lambda}^{-2}, \end{align}\] since \(||\Delta||_{\text{sp}}\le||\Delta||_2=o(1)\), with probability tending to 1.

The equation \((\ref{functionG})\) can be rewritten in the following form: \[\begin{gather} \label{functionG1} G(\Delta)=-\left(1+(1-\alpha)\gamma\right)\left(\text{trace} \Sigma \Delta - \dfrac{1}{4} \bar{\lambda}^{-2}||\Delta||_2^2 \right) + \\ \text{trace}\left(\Delta(S-\Sigma)\right) + \text{trace}\left(\Sigma \Delta\right) + \gamma \alpha \left(||\Omega+\Delta||_1-||\Omega||_1\right) = \\ \left(1+(1-\alpha)\gamma\right)\dfrac{1}{4} \bar{\lambda}^{-2}||\Delta||_2^2 -(1-\alpha)\gamma \text{trace} \Sigma \Delta + \\ \text{trace}\left(\Delta(S-\Sigma)\right) + \gamma \alpha \left(||\Omega+\Delta||_1-||\Omega||_1\right). \end{gather}\tag{20}\]

For an index set \(S\) and a matrix \(A=[a_{ij}]\), we denote \(A_S=[a_{ij}I((i,j)\in S)]\), where \(I(\cdot)\) is an indicator function. Recall \(S=\{(i,j):\Omega_{ij}\not=0\}\), and \(\bar{S}\) is its complement. It is straightforward that \(||\Omega_{\bar{S}}||_1 = 0\). From the triangular inequality, we have \[\begin{gather} ||\Omega+\Delta||_1-||\Omega||_1= ||\Omega_S+\Delta_S||_1+||\Delta_{\bar{S}}||_1 - ||\Omega_S||_1 \ge \\ ||\Omega_S||_1 - ||\Delta_S||_1+||\Delta_{\bar{S}}||_1 - ||\Omega_S||_1 = ||\Delta_{\bar{S}}||_1 - ||\Delta_S||_1 \end{gather}\] Next, we consider the term \(\text{trace}\left(\Delta(S-\Sigma)\right)\). Following the results of [55], we obtain: \[\begin{align} \text{trace}\left(\Delta(S-\Sigma)\right) \le ||S-\Sigma||_\infty ||\Delta||_1 = O_P\left(\sqrt{\dfrac{\log p}{n}}\right)||\Delta||_1 \le C_1\sqrt{\dfrac{\log p}{n}}||\Delta||_1, \end{align}\] with probability tending to 1, where \(C_1\) is some positive value. Note that \(||\Delta||_1 = ||\Delta_{\bar{S}}||_1 + ||\Delta_S||_1\).

The lower bound of the term \(\text{trace} \Sigma \Delta\) is obtained using Cauchy-Schwartz inequality: \[\begin{align} \text{trace} \Sigma \Delta \le ||\Sigma||_2 ||\Delta||_2, \end{align}\] also with probability tending to 1. Now, take \(\gamma = \dfrac{C_1}{\epsilon}\sqrt{\dfrac{\log p}{n}}\), where \(\epsilon\) is a small positive value. Combining the obtained lower bounds, we get from (20 ): \[\begin{gather} \label{functionG2} G(\Delta)\ge \left(1+(1-\alpha)\dfrac{C_1}{\epsilon}\sqrt{\dfrac{\log p}{n}} \right)\dfrac{1}{4} \bar{\lambda}^{-2}||\Delta||_2^2 -(1-\alpha)\dfrac{C_1}{\epsilon}\sqrt{\dfrac{\log p}{n}} ||\Sigma||_2 ||\Delta||_2 -\\ C_1\sqrt{\dfrac{\log p}{n}}\left(||\Delta_{\bar{S}}||_1 + ||\Delta_S||_1\right) + \dfrac{C_1}{\epsilon}\sqrt{\dfrac{\log p}{n}} \alpha \left(||\Delta_{\bar{S}}||_1 - ||\Delta_S||_1 \right) = \\ \left(1+(1-\alpha)\dfrac{C_1}{\epsilon}\sqrt{\dfrac{\log p}{n}} \right)\dfrac{1}{4} \bar{\lambda}^{-2}||\Delta||_2^2 -(1-\alpha)\dfrac{C_1}{\epsilon}\sqrt{\dfrac{\log p}{n}} ||\Sigma||_2 ||\Delta||_2 -\\ ||\Delta_{\bar{S}}||_1C_1\sqrt{\dfrac{\log p}{n}} \left(1-\dfrac{\alpha}{\epsilon}\right) - ||\Delta_S||_1 C_1\sqrt{\dfrac{\log p}{n}} \left(1 + \dfrac{\alpha}{\epsilon}\right) \ge \\ \left(1+(1-\alpha)C_1\sqrt{\dfrac{\log p}{n}} \right)\dfrac{1}{4} \bar{\lambda}^{-2}||\Delta||_2^2 -(1-\alpha)C_1\sqrt{\dfrac{\log p}{n}} ||\Sigma||_2 ||\Delta||_2 - ||\Delta_S||_1 C_1\sqrt{\dfrac{\log p}{n}} \left(1 + \dfrac{\alpha}{\epsilon}\right). \end{gather}\tag{21}\] Note that \(||\Delta_S||_1\le\sqrt{s}||\Delta_S||_2\le\sqrt{s}||\Delta||_2\), which gives us \[\begin{gather} G(\Delta)\ge \left(1+(1-\alpha)\dfrac{C_1}{\epsilon}\sqrt{\dfrac{\log p}{n}} \right)\dfrac{1}{4} \bar{\lambda}^{-2}||\Delta||_2^2 -(1-\alpha)\dfrac{C_1}{\epsilon}\sqrt{\dfrac{\log p}{n}} ||\Sigma||_2 ||\Delta||_2 - \\ C_1\sqrt{\dfrac{s\log p}{n}} \left(1 + \dfrac{\alpha}{\epsilon}\right) ||\Delta||_2 \ge \\ ||\Delta||_2^2 \biggl( \Bigl(1+(1-\alpha)\dfrac{C_1}{\epsilon}\sqrt{\dfrac{\log p}{n}} \Bigl)\dfrac{1}{4} \bar{\lambda}^{-2} -(1-\alpha)\dfrac{C_1}{\epsilon}\sqrt{\dfrac{(s+p)\log p}{n}} ||\Sigma||_2 ||\Delta||_2^{-1} - \\ C_1\sqrt{\dfrac{(s+p)\log p}{n}} \left(1 + \dfrac{\alpha}{\epsilon}\right) ||\Delta||_2^{-1} \biggl) = ||\Delta||_2^2 \biggl( \Bigl(1+(1-\alpha)\dfrac{C_1}{\epsilon}\sqrt{\dfrac{\log p}{n}} \Bigl)\dfrac{1}{4} \bar{\lambda}^{-2} - \\ (1-\alpha)\dfrac{C_1}{M\epsilon} ||\Sigma||_2 - \dfrac{C_1}{M}\left(1 + \dfrac{\alpha}{\epsilon}\right) \biggl), \end{gather}\] which is positive for a sufficiently large \(M\). This establishes the convergence rate \((\ref{convergence})\). An open question remains regarding the degree of stringency of the positivity of \(G(\Delta)\) depending on \(\alpha\) and the limitations on \(||\Sigma||_2\). We leave this problem to be studied in the future.

References↩︎