A greedy algorithm for sparse precision matrix approximation

TL;DR

This paper presents GISS$^{{ ho}}$, a fast greedy algorithm for sparse precision matrix estimation that combines $l_1$ minimization with computational efficiency, outperforming existing methods in accuracy and speed.

Contribution

The paper introduces GISS$^{{ ho}}$, a novel greedy algorithm for sparse precision matrix estimation, with proven convergence and superior numerical performance.

Findings

GISS$^{{ ho}}$ outperforms ADMM and HTP in accuracy.

The algorithm has favorable convergence properties.

Numerical experiments demonstrate efficiency and advantages over existing methods.

Abstract

Precision matrix estimation is an important problem in statistical data analysis. This paper introduces a fast sparse precision matrix estimation algorithm, namely GISS$^{{\rho}}$, which is originally introduced for compressive sensing. The algorithm GISS$^{{\rho}}$ is derived based on $l_1$ minimization while with the computation advantage of greedy algorithms. We analyze the asymptotic convergence rate of the proposed GISS$^{{\rho}}$ for sparse precision matrix estimation and sparsity recovery properties with respect to the stopping criteria. Finally, we numerically compare GISS$^{\rho}$ to other sparse recovery algorithms, such as ADMM and HTP in three settings of precision matrix estimation. The numerical results show the advantages of the proposed algorithm.

Tables

Table 1. Table 1: Case 1

Dimension	Algorithms	$L_0$ Norm (1e-8)	$L_0$ Norm (1e-4)	Time consuming (s)	Relative Error (Frobenius Norm)
200	GISS	598	598	0.1084	9.10e-10
	HTP	598		0.0826	9.10e-10
	L1-magic	1780		0.4679	1.00e-09
	${\text{ADMM}}_{{\lambda }_n}$	1388		0.2649	8.14e-10
	ADMM	1390		0.2615	8.99e-10
400	GISS	1198	1198	0.1002	9.09e-10
	HTP	1198		0.0975	9.09e-10
	L1-magic	3580		2.7838	1.00e-09
	${\text{ADMM}}_{{\lambda }_n}$	2788		1.5010	8.05e-10
	ADMM	2792		1.5014	8.98e-10
600	GISS	1798	1798	0.1479	9.09e-10
	HTP	1798		0.1567	9.09e-10
	L1-magic	5380		16.4614	1.04e-09
	${\text{ADMM}}_{{\lambda }_n}$	4188		9.7355	8.02e-10
	ADMM	4192		10.4155	8.97e-10
800	GISS	2398	2398	0.2518	9.09e-10
	HTP	2398		0.3942	9.09e-10
	L1-magic	7180		113.4222	1.04e-09
	${\text{ADMM}}_{{\lambda }_n}$	5588		36.3055	8.01e-10
	ADMM	5592		37.9462	8.97e-10
1000	GISS	2998	2998	0.4105	9.09e-10
	HTP	2998		0.7366	9.09e-10
	L1-magic	8980		392.0562	1.04e-09
	${\text{ADMM}}_{{\lambda }_n}$	6988		74.3292	8.00e-10
	ADMM	6992		77.3578	8.97e-10
2000	GISS	5998	5998	2.4252	9.09e-10
	HTP	5998		5.4236	9.09e-10
	L1-magic	17980		6.1685e+03	1.04e-09
	${\text{ADMM}}_{{\lambda }_n}$	13988		709.4917	7.98e-10
	ADMM	13992		733.7321	8.97e-10

Table 2. Table 2: Case 2, p = 200 𝑝 200 p=200 , λ 200 subscript 𝜆 200 \lambda_{200} = 0.7 ​ log ​ p n 0.7 log 𝑝 𝑛 0.7\sqrt{\frac{\mathrm{log}p}{n}} with sample number n 𝑛 n = 3000.

$p$=200	Frobenius norm	Matrix $L_1$ norm	Operator norm	$l_{\infty }$ Norm	Time consuming (s)
GISS	0.5457(0.008)	0.1084(0.016)	0.0125(0.003)	0.1314(0.003)	0.1219
HTP	0.2485(0.007)	0.0884(0.016)	0.0571(0.003)	0.1088(0.003)	0.3885
ADMM	0.3343(0.014)	0.2107(0.035)	0.0793(0.008)	0.1412(0.003)	0.1308
ADMM${}_{{\lambda }_n}$	1.2255(0.009)	0.5865(0.038)	0.2515(0.007)	0.1530(0.003)	0.4589

Table 3. Table 3: Case 2, p = 400 𝑝 400 p=400 , λ 400 subscript 𝜆 400 \lambda_{400} = 0.7 ​ log ​ p n 0.7 log 𝑝 𝑛 0.7\sqrt{\frac{\mathrm{log}p}{n}} with sample number n 𝑛 n = 5000.

$p$=400	Frobenius norm	Matrix $L_1$ norm	Operator norm	$l_{\infty }$ Norm	Time consuming (s)
GISS	1.4850(0.011)	0.2796(0.028)	0.1253(0.004)	0.1478(0.002)	0.8004
HTP	0.4452(0.009)	0.1192(0.018)	0.1187(0.003)	0.1056(0.003)	6.3441
ADMM	0.9107(0.023)	0.1878(0.025)	0.1896(0.007)	0.1940(0.004)	0.8222
ADMM${}_{{\lambda }_n}$	2.0134(0.012)	0.5611(0.030)	0.3289(0.005)	0.1606(0.002)	8.4315

Table 4. Table 4: Case 2, p = 600 𝑝 600 p=600 , λ 600 subscript 𝜆 600 \lambda_{600} = 0.7 ​ log ​ p n 0.7 log 𝑝 𝑛 0.7\sqrt{\frac{\mathrm{log}p}{n}} with sample number n 𝑛 n = 8000.

$p$=600	Frobenius norm	Matrix $L_1$ norm	Operator norm	$l_{\infty }$ Norm	Time consuming (s)
GISS	1.6309(0.010)	0.1536(0.022)	0.0171(0.003)	0.1290(0.001)	3.4359
HTP	0.4823(0.007)	0.2450(0.027)	0.2918(0.003)	0.0845(0.001)	43.1578
ADMM	2.1578(0.023)	0.3680(0.031)	0.0323(0.007)	0.1759(0.001)	12.7318
ADMM${}_{{\lambda }_n}$	2.2408(0.009)	0.5123(0.020)	0.2902(0.004)	0.1341(0.001)	13.6197

Table 5. Table 5: Case 2, p = 200 , 400 , 600 𝑝 200 400 600 p=200,~{}400,~{}600 .

	TP(%)			TN(%)
	200	400	600	200	400	600
GISS	99.93	98.43	94.29	99.71	99.93	100
HTP	99.95	97.95	92.35	99.81	99.99	100
ADMM	99.75	96.01	91.89	99.1	99.94	99.97
ADMM${}_{{\lambda }_n}$	99.97	99.09	96.52	98.69	99.84	99.99

Table 6. Table 6: Case 2, p = 60 𝑝 60 p=60 , λ 60 subscript 𝜆 60 \lambda_{60} = 0.001 ​ log ​ p n 0.001 log 𝑝 𝑛 0.001\sqrt{\frac{\mathrm{log}p}{n}} with sample number n 𝑛 n = 30.

$p$=60	Frobenius norm	Matrix $L_1$ norm	Operator norm	$l_{\infty }$ Norm
GISS	1.3570(0.14)	3.8780(0.21)	1.1073(0.14)	1.4990(0.15)
HTP	8.4734(0.01)	2.8927(0.01)	2.2451(0.01)	0.8689(0.01)
ADMM	6.8352(0.00)	5.8557(0.02)	0.4718(0.00)	1.0770(0.01)
ADMM${}_{{\lambda }_n}$	6.8379(0.00)	5.8608(0.02)	0.4719(0.00)	1.0770(0.01)

Table 7. Table 7: Case 2, p = 100 𝑝 100 p=100 , λ 100 subscript 𝜆 100 \lambda_{100} = 0.001 ​ log ​ p n 0.001 log 𝑝 𝑛 0.001\sqrt{\frac{\mathrm{log}p}{n}} with sample number n 𝑛 n = 50.

$p$=100	Frobenius norm	Matrix $L_1$ norm	Operator norm	$l_{\infty }$ Norm
GISS	1.6411(0.14)	5.1000(0.22)	0.9443(0.10)	1.4988(0.12)
HTP	11.2055(0.01)	4.1144(0.01)	2.7254(0.01)	0.9448(0.01)
ADMM	12.7463(0.00)	9.2365(0.02)	0.7048(0.00)	1.4200(0.01)
ADMM${}_{{\lambda }_n}$	12.7505(0.00)	9.2454(0.02)	0.7048(0.00)	1.4200(0.01)

Equations

$${\Sigma }_n=\frac{1}{n-1}\sum _{k=1}^n({\mathbf{X}}_k-\bar{\mathbf{X}})({\mathbf{X}}_k-\bar{\mathbf{X}}{)}^T$$

$$\Vert \Omega {\Vert }_1$$

$$|{\Sigma }_n\Omega -I{|}_{\infty }\le {\lambda }_n,\Omega \in {\mathbb{R}}^{p\times p}$$

$$\hat{\Omega }=({\hat{\omega }}_{ij}),\text{where}$$

$${\hat{\omega }}_{ij}={\hat{\omega }}_{ji}={\omega }_{ij}\mathbb{I}(|{\omega }_{ij}|\le |{\omega }_{ji}|)+{\omega }_{ji}\mathbb{I}(|{\omega }_{ij}|>|{\omega }_{ji}|)\text{for}i,j=1,2,\cdots ,p.$$

$$|{\mathbf{\beta }}_{\mathbf{i}}{|}_1$$

$$|{\Sigma }_n{\mathbf{\beta }}_{\mathbf{i}}-{\mathbf{e}}_{\mathbf{i}}{|}_{\infty }\le {\lambda }_n,{\beta }_i\in {\mathbb{R}}^p,$$

$${\Sigma }_{n,\gamma }$$

$$|{\mathbf{\beta }}_{\mathbf{i}}{|}_1$$

$${\Sigma }_{n,\gamma }{\mathbf{\beta }}_{\mathbf{i}}={\mathbf{e}}_{\mathbf{i}},{\mathbf{\beta }}_{\mathbf{i}}\in {\mathbb{R}}^p$$

$$|\mathbf{\beta }{|}_1$$

$$\Sigma \mathbf{\beta }=\mathbf{e}$$

$${\mathbf{\beta }}^{k+1}=\arg \min |\mathbf{\beta }{|}_1+\frac{\lambda }{2}|\Sigma \mathbf{\beta }-\mathbf{e}{|}_2^2-⟨{\mathbf{p}}^k,\mathbf{\beta }⟩,$$

$${\mathbf{p}}^{k+1}={\mathbf{p}}^k+\lambda {\Sigma }^T(\mathbf{e}-\Sigma {\mathbf{\beta }}^{k+1}),$$

$$\frac{{\mathbf{p}}^{k+1}-{\mathbf{p}}^k}{\lambda }={\Sigma }^T(\mathbf{e}-\Sigma {\mathbf{\beta }}^{k+1}).$$

$${\partial }_t\mathbf{p}(t)={\Sigma }^T(\mathbf{e}-\Sigma \mathbf{\beta }(t))$$

$$\mathbf{p}(t)\in \partial |\beta (t){|}_1.$$

$$\displaystyle\bm{p}\in\partial|\bm{\beta}|_{1}\Leftrightarrow\left\{\begin{array}[]{lll}\bm{p}_{i}=\mathrm{sign}(\bm{\beta}_{i}),&\mathrm{if}\ \bm{\beta}_{i}\neq 0,\\ |\bm{p}_{i}|\leq 1,&\mathrm{otherwise}.\end{array}\right.$$

$$|{\mathbf{\beta }}_{\mathbf{i}}{|}_1$$

$${\Sigma }_n{\mathbf{\beta }}_{\mathbf{i}}={\mathbf{e}}_{\mathbf{i}},{\mathbf{\beta }}_{\mathbf{i}}\in {\mathbb{R}}^p$$

$$\mathrm{E}{e}^{t(X_i-{\mu }_i{)}^2}\le K\text{textless}\infty \mathrm{for}\mathrm{all}|t|\le \eta ,\mathrm{for}\mathrm{all}i,$$

$$\mathrm{E}|X_i-{\mu }_i{|}^{4\tilde{\gamma }+4+\delta }\le K\mathrm{for}\mathrm{all}i.$$

$$\kappa ={\max }_{ij}\mathrm{E}[(X_i-{\mu }_i)(X_j-{\mu }_j)-{\sigma }_{ij}^0{]}^2=\text{mathchar}58{\max }_{ij}{\kappa }_{ij}.$$

$$U\text{mathchar}58=U(q,s_0(p))=\{\Omega \text{mathchar}58\Omega \succ 0,\Vert \Omega {\Vert }_{L_1}\le M,{\max }_{1\le i\le p}\sum _{j=1}^p|{\omega }_{ij}{|}^q\le s_0(p)\}$$

$$|\hat{\Omega }-{\Omega }_0{|}_{\infty }\le C_0\left(3+\frac{1+3p{\lambda }_n}{1-p{\lambda }_n}\right)M^2\sqrt{\frac{\log p}{n}},$$

$$|\hat{\Omega }-{\Omega }_0{|}_{\infty }\le D_0\left(3+\frac{1+3p{\lambda }_n}{1-p{\lambda }_n}\right)M^2\sqrt{\frac{\log p}{n}},$$

$${\sup }_{{\Omega }_0\in U}\mathrm{E}|{\hat{\Omega }}_{\gamma }-{\Omega }_0{|}_{\infty }^2=O\left(M^4\frac{\log p}{n}\right),$$

$$\mathbf{r}=\hat{\Sigma }{\mathbf{\beta }}^{\ast }-\mathbf{e}\sim N(0,{ϵ}^2{I}_p).$$

$$\mu \text{mathchar}58={\max }_{i,j}\left|\frac{1}{p}⟨{\hat{\Sigma }}_i,{\hat{\Sigma }}_j⟩\right|\text{textless}\frac{1}{2s-1},s=|S|.$$

$$\theta =1-\mu (s-1)$$

$$\vartheta =\frac{1-\mu (2s-1)}{1-\mu (s-1)}.$$

$$P(\mathbf{r}\in B_2)\ge 1-\frac{1}{p}$$

$$P(\mathbf{r}\in B_{\infty }(\varsigma ))\ge 1-\frac{1}{2p^{\varsigma }\sqrt{\pi \log p}}$$

$${\mathbf{\beta }}_{\min }^{\ast }\ge \frac{2ϵ}{\sqrt{\theta }}\left(\sqrt{1+2\sqrt{\frac{\log p}{p}}}+\sqrt{\frac{\log s}{p}}\right)$$

$${\mathbf{\beta }}_{\min }^{\ast }\ge \frac{2ϵ\sqrt{{\max }_i\Vert {\hat{\Sigma }}_i{\Vert }_p(1+\varsigma )s\log p}}{p\theta \Vert \hat{\Sigma }{\Vert }_{L_1}}+2ϵ\sqrt{\frac{\log s}{p\theta }}.$$

Full text

K****eywords Precision matrix estimation, CLIME estimator, sparse recovery, inverse scale space method, greedy methods.

1 Introduction

Covariance matrix and precision matrix estimation are two important problems in statistics analysis and data science. The problems become more challenging in high-dimensional setting where the number of variable dimension $p$ is larger than the sample size $n$, hence the need for a fast, accurate and stable precision/covariance matrix estimation is necessary. In the high-dimensional setting, classical methods and theoretical results with fixed $p$ and large $n$ are no longer applicable. Another huge challenge due to high dimension is high computational cost. Therefore, effective model and method are urgent facing high-dimensional data challenge.

Denote $\bm{X}=(X_{1},X_{2},\cdots,X_{p})^{T}$ by a $p$ variate random vector. The covariance matrix and precision matrix can be traditionally denoted by $\Sigma_{0}$ and $\Omega_{0}=\Sigma^{-1}_{0}$ respectively. Assume an independent and identically distributed $n$ random samples $\{\bm{X}_{1},\bm{X}_{2},\cdots,\bm{X}_{n}\}$ are from the distribution of $\bm{X}$. The sample covariance matrix is the common method among the estimators of covariance matrix, which is defined as follows,

[TABLE]

where $\bar{\bm{X}}=\frac{1}{n}\sum_{k=1}^{n}\bm{X}_{k}$ denotes the sample mean. When $p$ is larger than $n$, it is obvious that $\Sigma_{n}$ is singular and the estimation for $\Omega_{0}$ naturally becomes unstable and imprecise.

Estimation of precision matrix in high-dimensional setting has been studied for a long time. For example, when the random variable $\bm{X}$ follows a certain ordering structure, methods based on banding the Cholesky factor of the inverse of sample covariance matrix were studied in [37, 2]. Penalized likelihood methods such as $l_{1}$-MLE type estimators were studied in [21, 14, 40] and the convergence rate in Frobenius norm was given by [34]. In [39], the authors established the convergence rate for sub-Gaussian distribution cases. For more restrictive conditions, such as mutual incoherence or irrepresentable conditions, [33] showed the convergence rates in elementwise $l_{\infty}$ norm and spectral norm. To overcome the drawbacks that $l_{1}$ penalty inevitably leads to biased estimation, nonconvex penalty such as SCAD penalty [24, 17] was proposed [18, 43], although it often requires high computational cost.

Recently, [11] proposed a new constrained $l_{1}$ minimization approach called CLIME for sparse precision estimation. Convergence rates in spectral norm, elementwise $l_{\infty}$ norm and Frobenius norm were established under weaker assumptions and shown to be faster than those $l_{1}$-MLE estimators when the population distributions have polynomial-type tails. In addition, CLIME provides an effective computational performance as the columns of precision matrix estimation can be independently obtained and accelerated by parallelized computing. However, in [11], each column is obtained by a solving a linear programming, which could be still time-consuming for high dimension. More efficient approach are still needed to be developed for practical high dimensional applications.

In compressive sensing and sparse optimization community, many algorithms and related theoretical results are developed for $l_{1}$ minimization optimization problems [1, 6, 7, 23, 30, 38, 41]. Greedy inverse scale space flows (GISS) [27], originally stems from the adaptive inverse scale (aISS) method [5], is a new sparse recovery approach combining the idea of greedy approach and $\ell_{1}$ minimization. The advantage of GISS method was efficiency compared to the aISS method. GISSρ with $\rho$ being an acceleration factor, as a variant of GISS, can further accelerate sparse solution recovery by increasing the support of the current iterate by many indices at once.

In this article, we take the advantages of CLIME estimator framework and GISSρ algorithm and propose a new efficient approach for sparse precision matrix estimation. More specifically, we transfer the the constraint bound in CLIME estimator to a turning parameter in the stop criteria in GISSρ method. Convergence result in elementwise $l_{\infty}$ norm is established under weaker assumptions same as [11]. Moreover, under the Gaussian noise setting, the stopping time based on the data term is also obtained. The numerical experiments shows the competitive advantage of computation time, on obtaining the same level of sparsity and accuracy compared to other existing methods.

The rest of the paper is organized as follows. In Section 2, we firstly introduce basic notations and simply revisit the CLIME estimator. In Section 3, we present the basic idea of our new method derived from CLIME estimator and GISSρ algorithm. In Section 4, we establish the theoretical analysis with assumptions. Section 5 presents the numerical results including simulated experiments and application on real data. We provide the discussion and conclusion in Section 6 and the proof of the main results can be found in Appendix.

2 Preliminary

2.1 Notations and definitions

Before presenting our proposed precision matrix estimator, we first introduce some essential notations and definitions.

For a vector $\bm{a}=(a_{1},a_{2},\ldots,a_{n})^{T}\in\mathbb{R}^{n}$, we define $|\bm{a}|_{1}=\sum_{j=1}^{n}|a_{j}|$, $|\bm{a}|_{2}=\sqrt{\sum_{j=1}^{n}a_{j}^{2}}$ and $|\bm{a}|_{\infty}=\mathrm{max}_{1\leq j\leq n}~{}|a_{j}|$ to denote the vector norms. Define the two types of inner products between two vectors by $\langle\bm{a},\bm{b}\rangle=\bm{a}^{T}\bm{b}$ and $\langle\bm{a},\bm{b}\rangle_{n}=\frac{1}{n}\bm{a}^{T}\bm{b}$, and denote $\|\bm{\bm{a}}\|_{n}=\frac{1}{\sqrt{n}}|\bm{a}|_{2}$.

For a matrix $\bm{A}=(a_{ij})\in\mathbb{R}^{p\times p}$, we define the elementwise $l_{\infty}$ norm $|\bm{A}|_{\infty}=\mathrm{max}_{1\leq i\leq p,1\leq j\leq p}|a_{ij}|$, the spectral norm $\|\bm{A}\|_{2}=\mathrm{sup}_{|\bm{x}|_{2}\leq 1}|\bm{Ax}|_{2}$, the $L_{1}$ norm $\|\bm{A}\|_{L_{1}}=\mathrm{max}_{1\leq j\leq p}\sum_{i=1}^{p}|a_{ij}|$, the infinity norm $\|\bm{A}\|_{\infty}=\mathrm{max}_{1\leq i\leq p}\sum_{j=1}^{p}|a_{ij}|$, the Frobenius norm $\|\bm{A}\|_{F}=\sqrt{\sum_{i=1}^{p}\sum_{j=1}^{p}a_{ij}^{2}}$, and the elementwise $L_{1}$ norm $\|\bm{A}\|_{1}=\sum_{i=1}^{p}\sum_{j=1}^{p}|a_{ij}|$.

The $i$th row and $j$th column of the matrix $\bm{A}$ are denoted by $\bm{A}_{i,\cdot}$ and $\bm{A}_{\cdot,j}$ respectively and the transpose of $\bm{A}$ is denoted by $\bm{A}^{T}$. Let $\bm{A}^{*}=\frac{1}{p}A^{T}$ denote the adjoint operator of $\bm{A}$ with respect to the inner product $\langle\cdot,\cdot\rangle_{p}$. The inner product between two matrix $\bm{A}$ and $\bm{B}$ is denoted by $\langle\bm{A},\bm{B}\rangle=\sum_{i}(\bm{A}^{T}\bm{B})_{ii}$ with proper size. For two index sets $T$ and $T^{\prime}$, we use $\bm{A}_{TT^{\prime}}$ to denote the $|T|\times|T^{\prime}|$ matrix with rows and columns of $\bm{A}$ indexed by $T$ and $T^{\prime}$ respectively and $\bm{A}\succ 0$ means that $\bm{A}$ is positive definite. It also should be noted the $p\times p$ identity matrix is defined by $I$ while $\mathbb{I}$ denoted the indicator function.

For two real sequences $\{\epsilon_{n}\}$ and $\{\eta_{n}\}$, we write $\epsilon_{n}=O(\eta_{n})$ if there exists a constant $C$ such that $|\epsilon_{n}|\leq C|\eta_{n}|$ for large $n$, $\epsilon_{n}=o(\eta_{n})$ if $\mathrm{lim}_{n\rightarrow\infty}{\epsilon_{n}}/{\eta_{n}}=0$.

2.2 CLIME estimator

The CLIME estimator in [11] proposed to obtain a precision matrix estimation via solving the following $l_{1}$ minimization:

[TABLE]

where $\Sigma_{n}\in\mathbb{R}^{p\times p}$ is the sample covariance matrix generated by $n$ data samples, $\Omega$ is the precision matrix to be estimated and the turning parameter $\lambda_{n}$ is set for controlling the approximation error under the elementwise $l_{\infty}$ norm.

In general, the solution $\hat{\Omega}_{1}=(\omega_{ij})$ of (1) is asymmetric, hence a symmetry strategy as followed was adopted in [11],

[TABLE]

The above defined $\hat{\Omega}$ is the final estimated precision matrix through CLIME estimator.

It is easy to see that the convex program (1) can be decomposed into the following $p$ vector convex minimization problems:

[TABLE]

for all $1\leq i\leq p$, where $\bm{e_{i}}$ is a standard unit vector in $\mathbb{R}^{p}$. Denote $\hat{\bm{\beta_{i}}}$ as the solution of (3), then the assembling of the $p$ vectors in the form of $[\hat{\bm{\beta}}_{1},\hat{\bm{\beta}}_{2},\cdots,\hat{\bm{\beta}}_{p}]$ is a solution of (1), and vice visa.

For the case that the covariance matrix $\Sigma_{n}$ is not invertible, we can also consider a regularized

[TABLE]

for $\gamma>0$. Our GISSρ estimator can be obtained via solving the following minimization problem:

[TABLE]

for all $1\leq i\leq p$ with the $l_{\infty}$ norm stopping criterion $i.e.$ $|\Sigma_{n,\gamma}\bm{\beta_{i}}-\bm{e_{i}}|_{\infty}\leq\lambda_{n}$.

3 Sparse precision matrix estimation by GISSρ

Before presenting the GISSρ estimator for sparse precision matrix, we first introduce the basic ideas of GISS and GISSρ. It stems from the adaptive inverse scale (aISS) method [5], which was firstly proposed for solving compressive sensing problem. In the following, we present the algorithms in the context of precision matrix estimation for the consistency of notations.

3.1 Inverse scale space (ISS) based algorithms

The aISS algorithm was proposed to solve the equality constrained $l_{1}$ minimization (basic pursuit) problem. In the context of precision matrix estimation, we consider

[TABLE]

where $\Sigma\in\mathbb{R}^{p\times p}$ denotes the covariance matrix, $\bm{\beta}\in\mathbb{R}^{p}$ represents one column of the sparse precision matrix to be reconstructed and $\bm{e}\in\mathbb{R}^{p}$ denotes one column of identity matrix $I$.

For the $l_{1}$ minimization optimization problem for a general matrix $\Sigma$ and vector $\bm{e}$, Bregman iteration (BI) and augmented Lagrangian (AL) method are two efficient equivalent iterative methods [29, 16, 38]. By applying BI to solve (6), we obtain the following iterative sequence:

[TABLE]

where $k$ is the iteration step, $\bm{p}$ is the dual variable introduced in the optimization problem and $\lambda>0$ is step size for dual variable update.

In the above iterative scheme, the optimal condition of (7a) yields $\bm{p}^{k+1}\in\partial|\bm{\beta}^{k+1}|_{1}$ on combining with Equation (7b). Equation (7b) can be reformulated as

[TABLE]

If we consider both primal variable $\bm{\beta}$ and dual variable $\bm{p}$ as functions of time $t$ and interpret $\lambda$ as the time step, BI can be transformed into a dynamic equation:

[TABLE]

The above differential inclusion is called inverse scale space flow (ISS) [4].

The main idea of aISS [5] is that the support of the solution $\bm{\beta}(t^{k})$ can be determined by the subgradient $\bm{p}(t^{k})$ at some time $t^{k}$ based on the following relation

[TABLE]

for all $1\leq i\leq p$, i.e. the support of $\bm{\beta}(t^{k})$ is restricted onto the set $I=\{i\mid|\bm{p}_{i}(t^{k})|=1\}$. Furthermore, the step size $t^{k}$ can be determined as the change of the signal is piecewise linear, and the signal was obtained via solving a low dimensional optimization problem under the constraint of the sign consistence between $\bm{p}_{i}(t^{k})$ and $\bm{\beta}_{i}(t^{k})$.

3.2 GISSρ estimator

Motivated by the similarities between greedy algorithms (e.g. OMP [32, 26, 36], WOMP [35, 20], CoSaMP [28] and HTP [19]) and $l_{1}$ minimization algorithms (e.g. aISS), the authors in [27] proposed GISS method which approximates aISS flow closely. By removing the sign consistence constraint in aISS, GISS method achieves a more efficient scheme and provide a posterior method to estimate the distance to the $l_{1}$ minimizer. GISS method not only inherits remarkable properties from aISS, such as finite steps convergence but also shows much faster computation speed [27]. The modified version, namely GISSρ, can further accelerate the computation speed by including larger set of index with a factor $\rho\geq 1$.

Because of the efficiency of GISSρ algorithm for sparse recovery, we propose to apply it to estimate the columns of precision matrix, based on the formulation (1). More specifically, we propose GISSρ estimator by applying GISSρ algorithm via solving the following minimization problem:

[TABLE]

for all $1\leq i\leq p$ with the $l_{\infty}$ norm stopping criterion $i.e.$ $|\Sigma_{n}\bm{\beta_{i}}-\bm{e_{i}}|_{\infty}\leq\lambda_{n}$. The detail for solving each subproblem is shown in Algorithm 1.

We note that as CLIME estimator, the estimation of each column can be performed in parallel. Also the difference between GISSρ estimator and CLIME estimator is the role of $\lambda_{n}$. More precisely, $\lambda_{n}$ in CLIME estimator is used as the model tolerance. By controlling the tuning parameter, CLIME estimator can explore broader region of the solution space and the optimization problem is solved by a primal dual interior point algorithm applied on the equivalent linear programming problem. For GISSρ estimator, $\lambda_{n}$ plays the role of the stopping criterion and the algorithm produces a sequence of estimators. As expected, from our observation, both in the theoretical analysis and simulation experiments, $\lambda_{n}$ shows similar results in GISSρ estimator comparing to CLIME estimator, while GISSρ estimator is more efficient, especially for large scale problem. The detailed numerical performance will be shown in Section 5.

4 Theoretical Analysis

4.1 Asymptotic analysis

In this section, we will present the asymptotic results of GISSρ estimator under standard moments conditions of the random variables. We denote the sample covariance matrix as \textbf{\Sigma_{n}}=(\hat{\sigma}_{ij})=(\bm{\hat{\sigma}}_{1},\bm{\hat{\sigma}}_{2},\ldots,\bm{\hat{\sigma}}_{p}), where $\bm{\hat{}}{\sigma}_{i},~{}i=1,\cdots,p$ are column vectors. The true covariance matrix as \textbf{\Sigma_{0}}=(\sigma_{ij}^{0}) and the expectation of a random variable (r.v.) $\bm{X}$ as E$\bm{X}=(\mu_{1},\ldots,\mu_{p})^{T}$. The following two cases are usually considered according to the moment conditions on $\bm{X}$.

(C1) Exponential-type tails: Suppose that there exist some $0\textless\eta\textless 1/4$ such that $\mathrm{log}p/n\leq\eta$ and

[TABLE]

where $K$ is a upper bounded constant.

(C2) Polynomial-type tails: Suppose that for some $\tilde{\gamma},\ c_{1}\textgreater 0,\ p\leq c_{1}n^{\tilde{\gamma}},$ and for some $\delta\textgreater 0$,

[TABLE]

The following quantity is to measure the relation with the variance $\hat{\sigma}_{ij}$ which is useful in the theoretical property analysis,

[TABLE]

The maximum value $\kappa$ captures the overall variability of $\Sigma_{n}$. With the condition (C1) and (C2), $\kappa$ is a bounded constant depending on $\tilde{\gamma},\delta,K$ [11].

Through the following convergence analysis, we consider the uniformity class of matrices to estimate the precision matrix $\Omega_{0}$,

[TABLE]

for $0\leq q<1$, where $\Omega=\mathrel{\mathop{\mathchar 58\relax}}(\omega_{ij})=(\bm{\omega}_{1},\bm{\omega}_{2},\cdots,\bm{\omega}_{p})$. For the special case $q=0$, $U(0,s_{0}(p))$ is a class of $s_{0}(p)$-sparse matrices.

Theorem 4.1**.**

Suppose that $\Omega_{0}\in U(q,s_{0}(p))$ and $\hat{\Omega}$ is the solution solved through GISSρ process with $\rho\geq 1$.

• •

Assume that (C1) holds. Let $\lambda_{n}=C_{0}M\sqrt{\frac{\mathrm{log}p}{n}}\textless\frac{1}{p}-\varrho$ for some $\varrho\textgreater 0$, where $C_{0}=2\eta^{-2}(2+\tau+\eta^{-1}e^{2}K^{2})^{2}$ and $\tau>0$. Then

[TABLE]

with probability at least $1-4p^{-\tau}$.

• •

Assume that (C2) holds. Let $\lambda_{n}=D_{0}M\sqrt{\frac{\mathrm{log}p}{n}}\textless\frac{1}{p}-\varrho$ for some $\varrho\textgreater 0$, where $D_{0}=\sqrt{(5+\tau)(\kappa+1)}$. Then

[TABLE]

with probability at least $1-O\left(n^{-\delta/8}+p^{-\tau/2}\right)$.

From the Theorem 4.1, we find that the convergence rate of our estimator is same to the one of CLIME estimator [11] in the sense of elementwise $l_{\infty}$ norm under two types of moment conditions, which consequently outperform $l_{1}$-MLE type estimators in the case of polynomial-type tails [33].

Similarly, we can obtain the following convergence results:

Theorem 4.2**.**

Suppose $\Omega_{0}\in U(q,s_{0}(p))$ and ($C1$) holds. Let $\lambda_{n}=C_{0}M\sqrt{\frac{\mathrm{log}p}{n}}\textless\frac{1}{p}-\varrho$ for some $\varrho\textgreater 0$ with $C_{0}$ defined in Theorem 4.1 and $\tau$ is sufficiently large. Let $\gamma=\sqrt{\frac{\mathrm{log}p}{n}}$. If $p\geq n^{\Xi}$ for some $\Xi\geq 0$. Then

[TABLE]

4.2 Sparsity recovery properties

In this subsection, we provide some analysis from the point of view of sparsity recovery in compressive sensing [35, 42, 15] as the original idea of GISSρ algorithm from the signal processing for the noisy free case. Since our method can be computed in the parallelism form, we consider only one column recovery of the precision for simplicity. Denote $\bm{r}=\hat{\Sigma}\bm{\beta}^{*}-\bm{e}$ where $\bm{e}\in\mathbb{R}^{p}$ stands for one column of the identity matrix ${I}_{p}$, $\hat{\Sigma}$ is the $p\times p$ sample or true covariance matrix and $\bm{\beta}^{*}$ is the related column of the precision matrix to be recovered. In the algorithm, the algorithm is stopped with the criteria $|\bm{r}|_{\infty}\leq\lambda_{n}$. In the above section, we have analyzed the convergence of the approximation in terms of $\lambda_{n}$. In this section, we will show two theoretical results on the sparsity recovery guarantee with the assumption that the residual is a gaussian random variable. More precisely, we analyze the setting that the linear operator as the covariance matrix or the approximation of covariance matrix as a linear operator, the columns of the precision matrix is the sparse vector to be recovered. The problem is reformulated as a sparse recovery problem as followed,

[TABLE]

The assumption of the residual or observation error is to be Gaussian with standard variation $\epsilon$. We assume that each $\bm{\beta}^{*}$ has less than $s\leq p$ nonzero components. For convenience, let $S$ = supp($\bm{\beta}^{*}$) be the support of the nonzero index set of $\bm{\beta^{*}}$ and $T$ be its complement set. $\hat{\Sigma}_{S}$ denotes the submatrix of $\hat{\Sigma}$ formed by the columns of $\hat{\Sigma}$ in $S$, which are assumed to be linearly independent. One can similarity define $\hat{\Sigma}_{T}$.

Before providing the sparsity recovery guarantee, we present necessary assumptions on the design matrix $i.e.$ covariance matrix $\hat{\Sigma}$.

(A1) Mutual Incoherence Condition:

[TABLE]

It can be shown [35, 5] that once A1 holds, then

[TABLE]

and

[TABLE]

where $\theta$ and $\vartheta$ satisfies $\hat{\Sigma}_{S}^{*}\hat{\Sigma}_{S}\geq\theta I$ and $\|\hat{\Sigma}_{T}^{*}\hat{\Sigma}_{S}^{\dagger}\|_{\infty}\leq 1-\vartheta$ for $\hat{\Sigma}_{S}^{\dagger}\mathrel{\mathop{\mathchar 58\relax}}=\hat{\Sigma}_{S}(\frac{1}{p}\hat{\Sigma}_{S}^{T}\hat{\Sigma}_{S})^{-1}$ defined in Restricted Strong Convexity and Irrepresentable Conditions respectively. Moreover, the above two popular conditions can be verified with condition A1, while they are more difficult to be checked in practice [35, 5].

In the following, we will show that with large probability, the solution with the stopping rule $|\bm{r}|_{\infty}\leq\lambda_{n}$ recovers the true subset of the support index, under the assumption of Gaussian residual.

We further define the residual of the iterate $\bm{r}(t)=\bm{e}-\hat{\Sigma}\bm{\beta}(t)$ and then define the largest and the smallest nonzero magnitudes of $\bm{\beta}^{*}$ by $\bm{\beta}^{*}_{\mathrm{max}}\mathrel{\mathop{\mathchar 58\relax}}=\mathrm{max}(|\bm{\beta}_{i}^{*}|\mathrel{\mathop{\mathchar 58\relax}}i\in S)$ and $\bm{\beta}^{*}_{\mathrm{min}}\mathrel{\mathop{\mathchar 58\relax}}=\mathrm{min}(|\bm{\beta}_{i}^{*}|\mathrel{\mathop{\mathchar 58\relax}}i\in S)$ respectively.

The following Lemma show that the Gaussian noise is essentially bounded [10, 9].

Lemma 4.3**.**

For $\bm{r}\sim N(\bm{0},\epsilon^{2}I_{p})$ satisfies

[TABLE]

where $B_{2}=\left\{\bm{r}\mathrel{\mathop{\mathchar 58\relax}}|\bm{r}|_{2}\leq\epsilon\sqrt{p+2\sqrt{p\mathrm{log}p}}\right\}$ and $B_{\infty}(\varsigma)=\left\{\bm{r}\mathrel{\mathop{\mathchar 58\relax}}|\hat{\Sigma}^{T}\bm{r}|_{\infty}\leq\epsilon\sqrt{2(1+\varsigma)\mathrm{log}p}\right\}$ with $\varsigma>0$.

For the two essentially bounded cases, we have

Theorem 4.4**.**

Assume that (A1) holds and $\bm{r}\sim N(\bm{0},\epsilon^{2}I_{p})$, then

• •

For $l_{2}$ bounded Gaussian noise case (refer to (19a)),

[TABLE]

The GISSρ algorithm with the stopping rule $\|\bm{r}(t)\|_{\infty}\leq\epsilon\sqrt{1+2\sqrt{\frac{\mathrm{log}p}{p}}}$ selects the true subset of $S$ with probability at least $1-\frac{\mathrm{exp}{\left(-\frac{p}{2}+\frac{1}{2}\left(1+2\sqrt{\frac{\mathrm{log}p}{p}}\right)+\frac{p}{2}\mathrm{log}\left(2-\frac{1}{p}\left(1+2\sqrt{\frac{\mathrm{log}p}{p}}\right)\right)\right)}}{1-\frac{1}{p}\left(1+2\sqrt{\frac{\mathrm{log}p}{p}}\right)}-O(p^{-1})$.

• •

For $l_{\infty}$ bounded Gaussian noise case (refer to (19b)),

[TABLE]

The GISSρ with the stopping rule $\|\bm{r}(t)\|_{\infty}\leq 2\epsilon\|\hat{\Sigma}\|_{L_{1}}^{-1}{\sqrt{\mathrm{max}_{i}\|\hat{\Sigma}_{i}\|_{p}{(1+\varsigma)}\mathrm{log}p}}$ selects the true subset $S$ with probability at least $1-O(p^{-1})$.

Theorem 4.4 introduces a stopping rule with $l_{\infty}$ norm which only depends the residual level assumption and the conditions on $\beta^{*}_{\mathrm{min}}$ to ensure that all significant variables are selected [9].

5 Numerical Experiments

5.1 Experiments Setup

In the data simulation, we consider the following three cases.

• •

Case 1: Solve the sparse matrix with a designed covariance matrix: $\Sigma_{ij}=0.5^{|i-j|}$. The exact precision matrix $\Omega=\Sigma^{-1}$ is given as

[TABLE]

• •

Case 2: The second simulation refers to [34]. Assume that the true precision matrix $\Omega=B+\delta I$, where each off-diagonal entry in $B$ is generated independently and equals 0.5 with probability 0.1 or 0 with probability 0.9. $\delta$ is chosen such that conditional number of $\Omega$ is equal to $p$ and the matrix is normalized to have the unit diagonals. Then we sample a set of data from the $p-$dimensional normal distribution $N(0,\Sigma)$ where $\Sigma=\Omega^{-1}$. The estimated covariance matrix from this set of generated data is used as $\Sigma_{n}$ for the recovery.

• •

Case 3: Application to the fMRI data for ADHD. Attention Deficit Hyperactivity Disorder (ADHD) is a mental disorder causing 5%-10% of school-age children in behavior controlling through the United States. The ADHD-200 project released a resting-state fMRI dataset of healthy controls and ADHD children for scientific research. The data we processed is from Kennedy Krieger Institute (KKI), one of participating centers. The preprocessing pipeline can be found in https://www.nitrc.org/plugins/mwiki/ index.php/neurobureau:AthenaPipeline. After preprocessing, we get 148 time points from each of 116 brain regions for each subject. We use the data of each subject to generate the sample covariance matrix for estimating the precision matrix which is similar to the methodology of numerical examples in [25].

The above three simulation setting are tested with different algorithms for a comparison to ours. The first one is Hard Thresholding Pursuit (HTP) proposed in  [19] for solving the basis pursuit problem similar to (13). The algorithm is a combination of the Iterative Hard Thresholding algorithm and the Compressive Sampling Matching Pursuit algorithm (Cosamp). The other algorithms that we draw into comparisons are based on $l_{1}$ regularizer, such as L1-magic  [12] and Alternating direction method of multipliers (ADMM) [3, 22, 13]. We denote $\text{ADMM}_{\lambda_{n}}$ for solving inequality constrained model (3) and ADMM for the equality constrained model (13) for comparison.

The computation time will be reported with a computation environment : MATLAB R2016b, Intel(R) Xeon(R) CPU E5-2689 v4 @3.10GHz.

5.2 Results

The result for Case 1 with different dimension $p$ is shown in the Table 1. The third column of the table shows the non-zero entries number of the different methods with the thresholding value $10^{-8}$ applied on the recovered precision matrix for those $l_{1}$ based methods. The third column shows the result with the thresholding value $10^{-4}$, as we find all the methods successfully recover the true support. The last two columns show the computation time and the relative error to the true solution in terms of Frobenius norm. We find that both our method GISSρ (here $\rho=1$) have good computation time performance, especially in the high dimensional settings.

For Case 2, we vary the dimension $p$ from $200$ to $600$ and replicate 20 times for each algorithms. We adopt 20 cores for a parallelization implementation for this case. In this numerical setting, we take $\rho=1.2,1.5,2.0>1$ in GISSρ for different $p$ and set 0.05 for the thresholding value. Table 2, 3, 4 and 5 show the averages errors (Mean) and standard errors (SE) of losses in terms of different matrix norms: Frobenius norm, $L_{1}$ norm, $L_{2}$ norm and elementwise $l_{\infty}$ norm. True positive (TP) number and True negative (TN) number are also computed. The last column is the average time consuming of each algorithm. For $p\geq n$ situation, we show the numerical results in Table 6 and 7 for $p=60,~{}100$ respectively. We set 0.01 for the thresholding value and we found the Frobenius norm results outperform other methods, although the SE are larger than others.

In Theorem 4.1, we showed the convergence rate of the estimator with respect to the elementwise $l_{\infty}$ norm. From Figure 1, we can also numerically observe the linear convergence rate of $L_{2}$ norm and Frobenius norm, although this can not be shown theoretically in the current framework.

For Case 3 where the real fMRI dataset for ADHD is used, we compare the performance of support recovery using the data from all subjects, and the results suggest that GISSρ has competitive performance with CLIME. Figure 2 compares the running time of GISSρ with different $\rho$, HTP algorithm, and $\text{ADMM}_{\lambda_{n}}$ algorithm. Similar to the procedure described before, for each subject, we use the each algorithm to recover a designated edge percentages: 10%, 20%, 30%, 40%, 50%. This plot shows that the running time of GISSρ grows quickly for $\rho=1$ when we need to recover more connections. However, when we change the value by taking $\rho$=2, we find both HTP and GISSρ ($\rho=2$) shows the similar performance. $\text{ADMM}_{\lambda_{n}}$ algorithm is the slowest among above three methods. We note that we omit the application of ADMM method to this dataset as the computation time performance is much worse than any other algorithms in Figure 2.

6 Discussion and Conclusion

This paper proposed a fast sparse algorithm GISSρ for estimating precision matrix, especially for the high dimensional problem setting. The GISSρ is a algorithm based on the $l_{1}$ minimization. However, it combines the advantage of the greedy algorithms. The asymptotic convergence and the sparse reovery properties of this algorithm are analyzed. Numerical simulations show that the proposed GISSρ estimator shows a good performance compared to HTP, ADMM and other $l_{1}$ based estimators, especially for high dimensional problems.

Appendix

Lemma 6.1**.**

Suppose that $\Omega_{0}\in U(q,s_{0}(p))$ and $\hat{\Omega}_{n,\gamma}$ is the solution solved by GISSρ for $\rho\geq 1$ and $\hat{\Sigma}_{n,\gamma}$ defined as (4) and symmetry operation defined in (2). If $\|\Omega_{0}\|_{L_{1}}(\mathrm{max}_{ij}|\hat{\sigma}_{ij}-\sigma_{ij}^{0}|+\gamma)\leq\lambda_{n}\leq\frac{1}{p}$, then we have

[TABLE]

where $C_{1}\leq 2(1+2^{1-q}+3^{1-q})\left(3M+{\frac{1+3p\lambda_{n}}{1-p\lambda_{n}}M}\right)^{1-q}$.

Proof.

According to the proof of Theorem 6 in [11], we denote $\hat{\Omega}_{n,\gamma}\mathrel{\mathop{\mathchar 58\relax}}=\hat{\Omega}_{\gamma}$. Let $\hat{\bm{\beta}}$ be the one column of $\Omega_{0}$ and the $\bm{\beta}$ be the final result of GISSρ estimator satisfying $|\Sigma_{n,\gamma}\bm{\beta}-\bm{e}|_{\infty}\leq\lambda_{n}$, then we have

[TABLE]

Transform the result (28h) into matrix form and let $\hat{\Omega}^{1}_{\gamma}$ be the matrix form of $\beta$, it follows that

[TABLE]

since

[TABLE]

when $A$ is symmetric matrix.

Let $\tilde{\bm{p}}$ be the subgradient of $\|\cdot\|_{1}$ at $\bm{\beta}$, then

[TABLE]

and since

[TABLE]

we have

[TABLE]

In the matrix form, we have

[TABLE]

Therefore,

[TABLE]

Combining (35) and (29), we have

[TABLE]

Following the symmetric operation of $\|\hat{\Omega}_{\gamma}\|_{L_{1}}$ introduced in (2), we have

[TABLE]

This establishes (27) with (36). ∎

Combining the proof of Theorem 1(a) and 1(b) and Theorem 4(a), 4(b) in [11], Theorem 4.1 can be derived with Lemma 6.1. Similarly, Theorem 4.2 can be derived based on Theorem 2 and Theorem 5 in [11] and Lemma 6.1.

In the following, we prove the sparse recovery property Theorem 4.4.

Proof.

We first prove the first statement in Theorem 4.4. Following the proof of Theorem 6.1 and 6.2 in [31], the proof is divided into two parts. The first part shows that GISSρ will stop at the chosen stopping criteria with high probability and the second part is to assure the process will not stop if there still have supports to choose.

Let $b_{\infty}$ be the stopping point in the form of $l_{\infty}$ norm i.e. $b_{\infty}=\epsilon\sqrt{1+2\sqrt{\frac{\mathrm{log}p}{p}}}$. ${P}_{t}=\hat{\Sigma}_{S(t)}\left(\hat{\Sigma}^{*}_{S(t)}\hat{\Sigma}_{S(t)}\right)^{-1}\hat{\Sigma}_{S(t)}^{*}$. Obviously, $P_{t}$ is a symmetric projection matrix. It follows from Lemma 4 in [8], we have

[TABLE]

which means the algorithm stops when the residual is less than $b_{\infty}$.

The following part is to show that the algorithm will not stop whenever there exits the $i\in S$ such that $\bm{\beta}_{i}^{*}(t)=0$. Similar to the proof of the Theorem 6.1 in [31], we have

[TABLE]

provided that ${\bm{\beta}}^{*}_{\mathrm{min}}\geq\frac{2\epsilon}{\sqrt{\theta}}\left(\sqrt{1+2\sqrt{\frac{\mathrm{log}p}{p}}}\right)$. Since

[TABLE]

so it suffices to have

[TABLE]

∎

The second statement of Theorem 4.4 can be similarly proved.