Precision Matrix Estimation with Noisy and Missing Data

TL;DR

This paper introduces an ADMM algorithm for estimating precision matrices from noisy or incomplete data, addressing optimization challenges and comparing its performance with existing methods, with applications to political voting networks.

Contribution

It develops a novel ADMM-based method for precision matrix estimation with indefinite inputs and nonconvex penalties, filling a gap in high-dimensional noisy data analysis.

Findings

The proposed algorithm effectively estimates precision matrices with noisy and missing data.

Empirical comparisons show tradeoffs between different estimation methods.

Application to US senators' voting data reveals meaningful network structures.

Abstract

Estimating conditional dependence graphs and precision matrices are some of the most common problems in modern statistics and machine learning. When data are fully observed, penalized maximum likelihood-type estimators have become standard tools for estimating graphical models under sparsity conditions. Extensions of these methods to more complex settings where data are contaminated with additive or multiplicative noise have been developed in recent years. In these settings, however, the relative performance of different methods is not well understood and algorithmic gaps still exist. In particular, in high-dimensional settings these methods require using non-positive semidefinite matrices as inputs, presenting novel optimization challenges. We develop an alternating direction method of multipliers (ADMM) algorithm for these problems, providing a feasible algorithm to estimate precision matrices with indefinite input and potentially nonconvex penalties. We compare this method with existing alternative solutions and empirically characterize the tradeoffs between them. Finally, we use this method to explore the networks among US senators estimated from voting records data.

Tables

Table 1. Table 1: The relative norm error and FPR + FNR FPR FNR \text{FPR}+\text{FNR} performance of the Kronecker sum estimator using different methods. Here we set A 𝐴 A to be from an AR ​ ( 0.5 ) AR 0.5 \text{AR}(0.5) model and choose B 𝐵 B from an Erdos-Renyi random graph. We set m = 400 𝑚 400 m=400 and let τ B = 0.5 subscript 𝜏 𝐵 0.5 \tau_{B}=0.5 . Metrics are reported as the minimum value over a range of penalty parameters λ 𝜆 \lambda . The MCP penalty is chosen with a = 2.5 𝑎 2.5 a=2.5 , and we set R = 1.5 ​ ∥ A ∥ 2 𝑅 1.5 subscript delimited-∥∥ 𝐴 2 R=1.5\lVert A\rVert_{2} .

n	${\tau }_B$	method	penalty	Frobenius	Spectral	Nuclear	FPRFNR
80	0.3	Nonproj	${\ell }_1$	0.422	0.598	0.406	0.107
			MCP	0.450	0.613	0.422	0.106
		Proj	${\ell }_1$	0.424	0.610	0.411	0.113
			MCP	0.444	0.616	0.429	0.111
		Nodewise	${\ell }_1$	0.391	0.517	0.383	0.130
160	0.3	Nonproj	${\ell }_1$	0.342	0.509	0.327	0.013
			MCP	0.363	0.518	0.345	0.013
		Proj	${\ell }_1$	0.356	0.525	0.343	0.016
			MCP	0.341	0.493	0.321	0.015
		Nodewise	${\ell }_1$	0.288	0.429	0.280	0.017
80	0.5	Nonproj	${\ell }_1$	0.469	0.642	0.452	0.174
			MCP	0.481	0.659	0.458	0.177
		Proj	${\ell }_1$	0.464	0.651	0.450	0.194
			MCP	0.483	0.658	0.467	0.197
		Nodewise	${\ell }_1$	0.466	0.600	0.455	0.250
160	0.5	Nonproj	${\ell }_1$	0.389	0.573	0.369	0.052
			MCP	0.422	0.596	0.393	0.054
		Proj	${\ell }_1$	0.407	0.593	0.384	0.056
			MCP	0.399	0.587	0.377	0.055
		Nodewise	${\ell }_1$	0.358	0.538	0.349	0.083

Table 2. Table 2: The relative norm error and FPR + FNR FPR FNR \text{FPR}+\text{FNR} performance of the missing data estimator using different methods. Here we set A 𝐴 A to be from an AR ​ ( 0.6 ) AR 0.6 \text{AR}(0.6) model and set m = 400 𝑚 400 m=400 . Recall that ζ 𝜁 \zeta is the sampling rate. Metrics are reported as the minimum value over a range of penalty parameters λ 𝜆 \lambda . The MCP penalty is chosen with a = 2.5 𝑎 2.5 a=2.5 , and we set R 𝑅 R to be 1.5 1.5 1.5 times the oracle value for each method.

$A$ Model	n	$\zeta$	method	penalty	Frobenius	Spectral	Nuclear	FPRFNR
$\text{AR}\left(0.6\right)$	80	0.9	Nonproj	${\ell }_1$	0.367	0.506	0.363	0.0089
				MCP	0.308	0.533	0.296	0.0088
			Proj	${\ell }_1$	0.377	0.520	0.375	0.0085
				MCP	0.308	0.527	0.284	0.0083
			Nodewise	${\ell }_1$	0.292	0.487	0.280	0.0097
	130	0.7	Nonproj	${\ell }_1$	0.397	0.597	0.388	0.017
				MCP	0.384	0.632	0.363	0.016
			Proj	${\ell }_1$	0.417	0.599	0.407	0.019
				MCP	0.348	0.626	0.326	0.018
			Nodewise	${\ell }_1$	0.356	0.592	0.347	0.029
	250	0.5	Nonproj	${\ell }_1$	0.420	0.619	0.403	0.028
				MCP	0.457	0.680	0.436	0.026
			Proj	${\ell }_1$	0.437	0.626	0.429	0.031
				MCP	0.391	0.600	0.369	0.032
			Nodewise	${\ell }_1$	0.412	0.632	0.400	0.078
	700	0.3	Nonproj	${\ell }_1$	0.431	0.633	0.411	0.043
				MCP	0.505	0.718	0.470	0.040
			Proj	${\ell }_1$	0.450	0.644	0.431	0.034
				MCP	0.422	0.664	0.391	0.031
			Nodewise	${\ell }_1$	0.555	0.704	0.517	0.131

Table 3. Table 3: The relative norm error and FPR + FNR FPR FNR \text{FPR}+\text{FNR} performance of the missing data estimator using different methods. Here we set A 𝐴 A to be from an Erdos-Renyi random graph and set m = 400 𝑚 400 m=400 . Recall that ζ 𝜁 \zeta is the sampling rate. Metrics are reported as the minimum value over a range of penalty parameters λ 𝜆 \lambda . The MCP penalty is chosen with a = 2.5 𝑎 2.5 a=2.5 , and we set R 𝑅 R to be 1.5 1.5 1.5 times the oracle value for each method.

$A$ Model	n	$\zeta$	method	penalty	Frobenius	Spectral	Nuclear	FPRFNR
ER	80	0.9	Nonproj	${\ell }_1$	0.398	0.426	0.369	0.133
				MCP	0.379	0.444	0.355	0.132
			Proj	${\ell }_1$	0.405	0.420	0.375	0.129
				MCP	0.367	0.383	0.346	0.126
			Nodewise	${\ell }_1$	0.349	0.357	0.334	0.160
	130	0.7	Nonproj	${\ell }_1$	0.409	0.495	0.372	0.137
				MCP	0.410	0.562	0.372	0.137
			Proj	${\ell }_1$	0.423	0.497	0.385	0.135
				MCP	0.388	0.465	0.354	0.131
			Nodewise	${\ell }_1$	0.372	0.463	0.346	0.194
	250	0.5	Nonproj	${\ell }_1$	0.421	0.556	0.379	0.163
				MCP	0.463	0.680	0.401	0.170
			Proj	${\ell }_1$	0.437	0.556	0.394	0.163
				MCP	0.406	0.535	0.364	0.171
			Nodewise	${\ell }_1$	0.431	0.654	0.376	0.241
	700	0.3	Nonproj	${\ell }_1$	0.427	0.604	0.383	0.193
				MCP	0.485	0.701	0.415	0.189
			Proj	${\ell }_1$	0.445	0.575	0.401	0.184
				MCP	0.423	0.638	0.380	0.191
			Nodewise	${\ell }_1$	0.500	0.719	0.413	0.276

Table 4. Table 4: Measures of the indefiniteness of the input matrix Γ ^ n subscript ^ Γ 𝑛 \hat{\Gamma}_{n} . σ i subscript 𝜎 𝑖 \sigma_{i} denote the eigenvalues of Γ ^ n subscript ^ Γ 𝑛 \hat{\Gamma}_{n} , while σ i + superscript subscript 𝜎 𝑖 \sigma_{i}^{+} denote the eigenvalues of Γ ^ n + superscript subscript ^ Γ 𝑛 \hat{\Gamma}_{n}^{+} as defined in Section  4.1 . We set m = 400 𝑚 400 m=400 . For data generated from each model, we report the most negative eigenvalue, the maimum eigenvalues of both the nonprojected and projected sample covariances, the sum of the negative eigenvalues, and the number of negative eigenvalues.

Model	$A$	n	$\min {\sigma }_i$	$\max {\sigma }_i$	$\max {\sigma }_i^{+}$
KS	$\text{AR}\left(0.5\right)$	$n=80,{\tau }_B=0.3$	-0.51	17.0	15.3	-100.5	320
		$n=160,{\tau }_B=0.3$	-0.42	10.3	9.6	-74.1	240
		$n=80,{\tau }_B=0.5$	-0.93	21.3	18.1	-170.1	320
		$n=160,{\tau }_B=0.5$	-0.78	12.0	10.7	-124.6	243
MD	$\text{AR}\left(0.6\right)$	$n=80,\zeta =0.9$	-0.26	14.2	13.6	-36.2	320
		$n=130,\zeta =0.7$	-0.63	12.3	11.0	-116.6	270
		$n=250,\zeta =0.5$	-1.19	11.4	9.7	-183.6	218
		$n=700,\zeta =0.3$	-2.17	9.2	7.5	-228.9	188
	ER	$n=80,\zeta =0.9$	-0.26	13.4	12.7	-36.6	320
		$n=130,\zeta =0.7$	-0.62	11.7	10.4	-116.7	270
		$n=250,\zeta =0.5$	-1.20	10.3	8.7	-180.7	214
		$n=700,\zeta =0.3$	-2.17	8.5	6.9	-223.0	184

Table 5. Table 5: Highest estimated partial correlations among Republicans (top) and Democrats (bottom).

Senators		States		Tea Party
Enzi	Barrasso	WY	WY	No-No
Crapo	Risch	ID	ID	No-No
Chambliss	Isakson	GA	GA	No-No
Snowe	Collins	ME	ME	No-No
DeMint	Lee	SC	UT	Yes-Yes
Collins	Brown	ME	MA	No-No
Lee	Paul	UT	KY	No-No
Cochran	Hoeven	MS	ND	No-No
Coburn	Johnson	OK	WI	Yes*-Yes
Cochran	Wicker	MS	MS	No-No
Reed	Whitehouse	RI	RI
Murray	Cantwell	WA	WA
Wyden	Merkley	OR	OR
Udall	Bennet	CO	CO
Leahy	Sanders	VT	VT
Carper	Coons	DE	DE
Mikulski	Akaka	MD	HI
Feinstein	Boxer	CA	CA
Baucus	Tester	MT	MT
Akaka	Gillibrand	HI	NY

Table 6. Table 6: Most negative estimated partial correlations among senators.

Senators		States		Party
Akaka	McCaskill	HI	MO	D
Hatch	Hoeven	UT	ND	R
Hoeven	Toomey	ND	PA	R
Alexander	Rubio	TN	FL	R
Levin	Tester	MI	MT	D
McCaskill	Udall	MO	NM	D
Corker	Wicker	TN	MS	R
Collins	Johnson	ME	WI	R
Cardin	Hagan	MD	NC	D
Bingaman	Casey	NM	PA	D

Equations

$$\hat{\Theta }\in \underset{\Theta ⪰0}{\mathrm{argmin}}\left\{\mathrm{tr}({\hat{\Gamma }}_n\Theta )-\log \det (\Theta )+g_{\lambda }(\Theta )\right\},$$

$$\hat{\Theta }\in \underset{\Theta ⪰0,\Vert \Theta {\Vert }_2\le R}{\mathrm{argmin}}\left\{\mathrm{tr}({\hat{\Gamma }}_n\Theta )-\log \det (\Theta )+g_{\lambda }(\Theta )\right\}.$$

$$f(I+t{v}_1{v}_1^T)$$

$$=\mathrm{tr}({\hat{\Gamma }}_n)+t\mathrm{tr}({\hat{\Gamma }}_n{v}_1{v}_1^T)-\log (1+t)+g_{\lambda }(t{v}_1{v}_1^T)$$

$$=\mathrm{tr}({\hat{\Gamma }}_n)+t{\sigma }_1-\log (1+t)+g_{\lambda }(t{v}_1{v}_1^T).$$

$$f(\Theta )=\mathrm{tr}({\hat{\Gamma }}_n\Theta )-\log \det (\Theta )+g_{\lambda }(\Theta )+\text{mathds}{1}_{{\mathcal{X}}_R}(\Theta )$$

$$T_{\rho }(A)$$

$$\text{ where }{\tilde{D}}_{ii}=\min \left\{\frac{M_{ii}+(M_{ii}^2+\frac{4}{\rho }{)}^{1/2}}{2},R\right\},$$

$${\hat{\beta }}_j$$

$${\hat{a}}_j$$

$${\hat{\Gamma }}_n=\frac{1}{n}{X}^{T}X\oslash M\text{ where }{M}_{k\ell }=\{\begin{array}{ll}{\zeta }_k& \text{if }k=\ell \\ {\zeta }_k{\zeta }_{\ell }& \text{if }k\ne \ell \end{array}$$

$$g_{\lambda }(w)=\{\begin{array}{ll}\lambda |w|& \text{if }|w|\le \lambda \\ -\frac{w^2-2a\lambda |w|+{\lambda }^2}{2(a-1)}& \text{if }\lambda <|w|\le a\lambda \\ \frac{(a+1){\lambda }^2}{2}& \text{if }a\lambda <|w|\end{array}$$

$$g_{\lambda }(w)=\text{sign}(w)\lambda {\int }_0^{|w|}{\left(1-\frac{z}{\lambda a}\right)}_{+}dz$$

$$f(\Theta )=\mathrm{tr}({\hat{\Gamma }}_n\Theta )-\log \det (\Theta )+g_{\lambda }(\Theta )+\text{mathds}{1}_{{\mathcal{X}}_R}(\Theta )$$

$$\hat{\Theta }$$

$$\displaystyle\qquad\qquad\qquad\qquad+\mathds{1}_{\mathcal{X_{R}}}(\Theta)+g_{\lambda}(V)\Big{\}}\;\text{ s.t. }\Theta=V$$

$${\mathcal{L}}_{\rho }(\Theta ,V,\Lambda )$$

$$+g_{\lambda }(V)+\frac{\rho }{2}\Vert \Theta -V{\Vert }_F^2+⟨\Lambda ,\Theta -V⟩$$

$$\displaystyle\begin{aligned} V^{k+1}=\operatorname*{arg\,min}_{V\in\mathbb{R}^{m\times m}}\Big{\{}&g_{\lambda}(V)+\frac{\rho}{2}\lVert\Theta^{k}-V\rVert_{F}^{2}\\ &\quad+\langle\Lambda^{k},\Theta^{k}-V\rangle\Big{\}}\end{aligned}$$

$$\displaystyle\begin{aligned} \Theta^{k+1}=\operatorname*{arg\,min}_{\Theta\in\mathbb{R}^{m\times m}}\Big{\{}&-\log\det\Theta+\operatorname*{tr}(\hat{\Gamma}_{n}\Theta)\\ &\quad+\mathds{1}_{\mathcal{X}_{R}}(\Theta)+\frac{\rho}{2}\lVert\Theta-V^{k+1}\rVert_{F}^{2}\\ &\quad+\langle\Lambda^{k},\Theta-V^{k+1}\rangle\Big{\}}\end{aligned}$$

$${\Lambda }^{k+1}={\Lambda }^k+\rho ({\Theta }^{k+1}-V^{k+1})$$

$$V^{k+1}=\underset{V\in {\mathbb{R}}^{m\times m}}{\mathrm{argmin}}\left\{\frac{1}{\rho }{g}_{\lambda }(V)+\frac{1}{2}{\Vert V-\frac{\rho {\Theta }^k+{\Lambda }^k}{\rho }\Vert }_F^2\right\}.$$

$${\text{Prox}}_{g_{\lambda }/\rho }(w)=\{\begin{array}{ll}0& \text{if }|w|\le \nu \lambda \\ w-\text{sign}(w)\nu \lambda & \text{if }\nu \lambda \le |w|\le (\nu +1)\lambda \\ \frac{w-\text{sign}(w)\frac{a\nu \lambda }{a-1}}{1-\frac{\nu }{a-1}}& \text{if }(\nu +1)\lambda \le |w|\le a\lambda \\ w& \text{if }a\lambda \le |w|\end{array}$$

$${\text{Prox}}_{g_{\lambda }/\rho }(w)=\{\begin{array}{ll}0& \text{if }|w|\le \nu \lambda \\ \frac{w-\text{sign}(w)\nu \lambda }{1-\nu /a}& \text{if }\nu \lambda \le |w|\le a\lambda \\ w& \text{if }a\lambda \le |w|\end{array}$$

$$\displaystyle\Theta^{k+1}=\operatorname*{arg\,min}_{\Theta\in\mathbb{R}^{m\times m}}\Big{\{}$$

$$\displaystyle+\frac{\rho}{2}\left\lVert\Theta-\frac{\rho V^{k+1}-\hat{\Gamma}_{n}-\Lambda^{k}}{\rho}\right\rVert_{F}^{2}\Big{\}}$$

$$T_{\rho }(A)=T_{\rho }(UM{U}^T)$$

$$=\underset{\Theta }{\mathrm{argmin}}\left\{-\log \det \Theta +\text{mathds}{1}_{{\mathcal{X}}_R}(\Theta )+\frac{\rho }{2}\Vert \Theta -A{\Vert }_F^2\right\}$$

$$=U\tilde{D}{U}^T$$

$$\text{ where }{\tilde{D}}_{ii}=\min \left\{\frac{M_{ii}+(M_{ii}^2+\frac{4}{\rho }{)}^{1/2}}{2},R\right\}$$

$$V^{k+1}$$

$${\Theta }^{k+1}$$

$${\Lambda }^{k+1}$$

$$T_{\rho }(A)$$

$$=\underset{\Theta }{\mathrm{argmin}}\left\{-\log \det \Theta +\text{mathds}{1}_{{\mathcal{X}}_R}(\Theta )+\frac{\rho }{2}\Vert \Theta -A{\Vert }_F^2\right\}$$

$$T_{\rho }(A)$$

$$-\rho \mathrm{tr}(\Theta A)+\text{mathds}{1}_{{\mathcal{X}}_R}(\Theta )$$

$$=\underset{\Theta =WD{W}^T}{\mathrm{argmin}}-\log \det D+\frac{\rho }{2}\mathrm{tr}(DD)$$

$$-\rho \mathrm{tr}(WD{W}^{T}UM{U}^T)+\text{mathds}{1}_{{\mathcal{X}}_R}(D)$$

Full text

1 Introduction

Undirected graphs are often used to describe high-dimensional distributions. Under sparsity conditions, these graphs can be estimated using penalized methods such as

[TABLE]

where $\hat{\Gamma}_{n}$ is the sample covariance or correlation matrix and $g_{\lambda}$ is a separable (entry-wise) sparsity-inducing penalty function. Although this approach has proven successful in a variety of application areas such as neuroscience and genomics, its soundness hinges on the positive semidefiniteness (PSD) of $\hat{\Gamma}_{n}$. If $\hat{\Gamma}_{n}$ is indefinite, the objective may be unbounded from below.

In order to ensure this penalized $M$-estimator is well-behaving, Loh and Wainwright (2015) impose a side constraint of the form $\rho(\Theta)<R$, where $\rho$ is a convex function. Here we focus on the estimator using the operator norm as a side constraint

[TABLE]

Loh and Wainwright (2017) adopt this method and show in theory the superior statistical properties of this constrained estimator. Their results suggest that the addition of a side constraint is not only sufficient but also almost necessary to effectively untangle the aforementioned complications.

Unfortunately, this additional constraint precludes using existing methods to solve the penalized objective with non-PSD input. To close this gap, we develop an alternating direction method of multipliers (ADMM) algorithm to implement (2) efficiently. We conduct empirical studies comparing this new method to several other precision matrix estimators. Our simulation study reveals several trends that are not present in the fully observed case. Finally, we illustrate the performance of our methods in analyzing the US senate voting data, uncovering both known and novel phenomena from the modern political landscape.

The remainder of this paper is organized as follows. In Section 2, we provide an overview of existing related work and describe in detail the optimization issues that arise from indefinite inputs and nonconvex penalties. In Section 3, we present the proposed ADMM algorithm and present some convergence results. Section 4 provides numerical examples and comparisons. Section 5 presents an exploratory analysis of US Senate voting records data using this method and details several interesting conclusions that can be drawn from the estimated graphs. Finally, we summarize the empirical results and their practical implications regarding choice of method in Section 6.

2 Problem formulation and existing work

There is a wide body of work proposing methods to perform precision matrix estimation in the fully observed case, including Meinshausen and Bühlmann (2006), Yuan and Lin (2007), Rothman et al. (2008), Friedman et al. (2008), Banerjee et al. (2008), and Zhou et al. (2010), most of which are essentially a $\ell_{1}$-penalized likelihood approach (1) which we will refer to as the graphical Lasso.

Recent work has focused on using nonconvex regularizers such as SCAD and MCP for model selection in the regression setting (Fan and Li, 2001; Zhang, 2010; Breheny and Huang, 2011; Zhang and Zhang, 2012). Loh and Wainwright (2015, 2017) extend this analysis to general $M$-estimators, including variants of the graphical Lasso objective, and show their statistical convergence and support recovery properties. Estimators with these penalties have been shown to attain model selection under weaker theoretical conditions, but require more sophisticated optimization algorithms to solve, such as the local linear approximation (LLA) method of Fan et al. (2014).

In a fully observed and noiseless setting, $\hat{\Gamma}_{n}$ is the sample covariance and guaranteed to be at least positive semidefinite. Then, if $g_{\lambda}$ is the $\ell_{1}$-penalty, the objective of (1) is convex and bounded from below. In this setting, one can show that for $\lambda>0$ a unique optimum $\hat{\Theta}$ exists with bounded eigenvalues and that the iterates for any descent algorithm will also have bounded eigenvalues (for example, see Lemma 2 in Hsieh et al., 2014).

When working with missing, corrupted, and dependent data, the likelihood is nonconvex, and the expectation-maximization (EM) algorithm has traditionally been used to perform statistical inference. However, in these noisy settings, the convergence of the EM algorithm is difficult to guarantee and is often slow in practice. For instance, Städler and Bühlmann (2012) implement a likelihood-based method for inverse covariance estimation with missing values, but their EM algorithm requires solving a full graphical Lasso optimization problem in each M-step.

An alternative approach is to develop $M$-estimators that account for missing and corrupted data. For graphical models, Loh and Wainwright (2015) establish that the graphical Lasso, including a version using nonconvex penalties, can be modified to accommodate noisy or missing data by adjusting the sample covariance estimate.

These modified estimators depend on the observation that statistical theory for the graphical Lasso generally requires that $\lVert\hat{\Gamma}_{n}-\Sigma\rVert_{\infty}$ converges to zero at a sufficiently fast rate (e.g. Rothman et al., 2008; Zhou et al., 2010; Loh and Wainwright, 2017). When considering missing or corrupted data, it is often possible to construct covariance estimates $\hat{\Gamma}_{n}$ that satisfy this convergence criteria but are not necessarily positive semidefinite. In fact, in high-dimensional settings $\hat{\Gamma}_{n}$ may even be guaranteed to be indefinite. Attempting to input these indefinite covariance estimates into the graphical Lasso, however, presents novel optimization issues.

Unbounded objective. When attempting to move beyond the $\ell_{1}$ penalized case with positive semidefinite input, the problem in (1) becomes unbounded from below, so an optimum may not necessarily exist. This issue comes from two potential sources: 1) negative eigenvalues in $\hat{\Gamma}_{n}$, or 2) zero eigenvalues combined with the boundedness of the nonconvex penalty $g_{\lambda}$. For example, consider the restriction of the objective in (1) to a ray defined by an eigenvalue-vector pair $\sigma_{1},v_{1}$ of $\hat{\Gamma}_{n}$:

[TABLE]

If $\sigma_{1}<0$, we see that $f$ is unbounded from below due to the $t\sigma_{1}$ and $-\log(1+t)$ terms. In fact, if $\sigma_{1}=0$ and $g_{\lambda}$ is bounded from above, as is the case when using standard nonconvex penalties, the objective is also unbounded from below.

So unboundedness can occur anytime there is a negative eigenvalue in the input matrix, or whenever there are zero eigenvalues combined with a nonconvex penalty function $g_{\lambda}$. Unboundedness creates optimization issues, as an optimum no longer necessarily exists.

Handling unboundedness. In order to guarantee that an optimum exists for (1), an additional constraint of the form $\rho(\Theta)\leq R$ can be imposed, where $\rho$ is some convex function. In this paper, we consider the estimator (2), which uses a side constraint of the form $\lVert\Theta\rVert_{2}\leq R$. Loh and Wainwright (2017) show the rates of convergence of this estimator (2) and show that it can attain model selection consistency and spectral norm convergence without the incoherence assumption when used with a nonconvex penalty (see Appendix E therein), but do not discuss implementation or optimization aspects of the problem.

To our knowledge, there is currently no feasible optimization algorithm for the estimator defined in (2), particularly when the input is indefinite. Loh and Wainwright (2015) present a composite gradient descent method for optimizing a subset of side-constrained versions of (1). However, their algorithm requires a side constraint of the form $\rho(\Theta)=\frac{1}{\lambda}(g_{\lambda}(\Theta)+\frac{\mu}{2}\lVert\Theta\rVert_{F}^{2})$, which does not include the spectral norm constraint and therefore cannot attain the better theoretical results it achieves (Section C.7 compares the performance of different side constraints). It may be possible to develop heuristic algorithms that alternate performing a proximal gradient update ignoring the side constraint and projecting to the constraint set, but as far as we know there has not been any analysis of algorithms of this type (we discuss this in more detail in Section C.4).

An alternative approach to solving this unbounded issue is to project the input matrix $\hat{\Gamma}_{n}$ to the positive semidefinite cone before inputting into (1). We discuss this further in Section 4.1, but this only solves the unbounded issue when using the $\ell_{1}$ penalty; nonconvex penalties still require a side constraint to have a bounded objective and therefore our algorithm is still useful even for the projected methods.

3 ADMM Algorithm

Our algorithm is similar to the algorithm in Guo and Zhang (2017), which applies ADMM to the closely related problem of condition number-constrained sparse precision matrix estimation using the same splitting scheme as below. We discuss their method in more detail in Section A.6. The following algorithm is specialized to the case where the spectral norm is used as the side constraint. In Section B we derive a similar ADMM algorithm that can be used for any side constraint with a computable projection operator.

Rewrite the objective from (2) as

[TABLE]

where $\mathcal{X}_{R}=\{\Theta:\Theta\succeq 0,\lVert\Theta\rVert_{2}\leq R\}$ and $\mathds{1}_{\mathcal{X}}(\Theta)=0\text{ if }\Theta\in\mathcal{X}\text{ and }\infty\text{ otherwise}$.

Let $\rho>0$ be a penalty parameter and let $\text{Prox}_{g_{\lambda}/\rho}$ be the prox operator of $g_{\lambda}/\rho$. We derive these updates for SCAD and MCP in Section A.2. Let $T_{\rho}(A)$ be the following prox operator for $-\log\det\Theta+\mathds{1}_{\mathcal{X}_{R}}(\Theta)$, which we derive in Section A.3,

[TABLE]

where $UMU^{T}$ is the eigendecomposition of $A$. Then the ADMM algorithm for solving (4), which we derive in Section A.2, is described in Algorithm 1. Computationally this algorithm is dominated by the eigendecomposition used to evaluate $T_{\rho}$, and therefore has a complexity of $O(m^{3})$, which matches the scaling of other graphical Lasso solvers (e.g. Meinshausen and Bühlmann, 2006; Friedman et al., 2008; Hsieh et al., 2014).

3.1 Convergence

The following proposition applies standard results on the convergence of ADMM for convex problems to show convergence when the $\ell_{1}$ penalty is used. Details are in Section A.4.

Proposition 1**.**

If the penalty is convex and satisfies the conditions in Section A.1, Algorithm 1 converges to a global minimum of (4).

Remark**.**

Regarding the nonconvex penalty, recent work has established ADMM convergence results in some nonconvex settings (see Hong et al., 2016; Wang et al., 2015), but to our knowledge there is no convergence result that encompasses this nonsmooth and nonconvex application. We can show convergence if a fairly strong assumption is made on the iterates, but we are currently working on extending existing results to this case.

Proposition 2 shows that any limiting point of Algorithm 1 is a stationary point of the original objective (4). This is proved in Section A.5. When using the $\ell_{1}$ penalty or a nonconvex penalty with $R\leq\sqrt{2/\mu}$, where $\mu$ is the weak convexity constant of $g_{\lambda}$, the objective $f$ is convex and therefore any stationary point is unique and also the global optimum. See Section C.5 for a more detailed discussion.

Proposition 2**.**

Assume that the penalty $g_{\lambda}$ satisfies the conditions in Section A.1. Then for any limit point $(\Theta^{*},V^{*},\Lambda^{*})$ of the ADMM algorithm defined in Algorithm 1, $\Theta^{*}$ is also a stationary point of the objective $f$ as defined in (4).

The assumptions on $g_{\lambda}$ in Section A.1 are the same as those assumed in Loh and Wainwright (2015, 2017), and are satisfied by the Lasso, SCAD, and MCP functions.

Note that if a limiting point is found to exist when using a nonconvex penalty the result in Proposition 2 will still hold. Empirically we find that the algorithm performs well and converges consistently when used with nonconvex penalties, but there is no existing theoretical guarantee that a limiting point of ADMM will exist in that setting.

4 Simulations

We evaluate the proposed estimators using the relative Frobenius norm and the sum of the false positive rate and false negative rate ($\text{FPR}+\text{FNR}$). We present results over a range of $\lambda$ values, noting that all the compared methods would use similar techniques to perform model tuning. Section C.1 presents an example of how to use BIC or cross-validation to tune these methods. We present results using covariance matrices from auto-regressive and Erdős-Rényi random graph models. See Section C for descriptions of these models as well as additional simulation results.

4.1 Alternative methods

When faced with indefinite input, there are two alternative graphical Lasso-style estimators that can be used besides (2), which involve either $\ell_{\infty}$ projection to the positive semidefinite cone or nodewise regression in the style of Meinshausen and Bühlmann (2006).

Projection. Given an indefinite input matrix $\hat{\Gamma}_{n}$, Park (2016) and Greenewald et al. (2017) propose performing the projection $\hat{\Gamma}_{n}^{+}=\operatorname*{arg\,min}_{\Gamma\succeq 0}\lVert\Gamma-\hat{\Gamma}_{n}\rVert_{\infty}$. They then input $\hat{\Gamma}_{n}^{+}$ into the optimization problem (1). This is similar to the projection done in Datta and Zou (2017). In terms of the upper bound on statistical convergence rates, this method pays a constant factor cost, though in practice projection may result in a loss of information and therefore a decrease in efficiency.

After projecting the input, existing algorithms can be used to optimize (1) with the $\ell_{1}$ penalty. However, as mentioned in Section 2, using a nonconvex penalty still leads to an unbounded objective and therefore still requires using our ADMM algorithm to solve (2).

Nodewise regression. Loh and Wainwright (2012) and Rudelson and Zhou (2017) both study the statistical and computational convergence properties of using errors-in-variables regression to handle indefinite input matrices in high-dimensional settings. Following the nodewise regression ideas of Meinshausen and Bühlmann (2006) and Yuan (2010), we can perform $m$ Lasso-type regressions to obtain estimates $\hat{\beta}_{j}$ and form estimates $\hat{a}_{j}$, where

[TABLE]

and combine to get $\tilde{\Theta}$ with $\tilde{\Theta}_{-j,j}=\hat{a}_{j}\hat{\beta}_{j}$ and $\tilde{\Theta}_{j,j}=-\hat{a_{j}}$. Finally, we symmetrize the result to obtain $\hat{\Theta}=\operatorname*{arg\,min}_{\Theta\in S^{m}}\lVert\Theta-\tilde{\Theta}\rVert_{1}$, where $S^{m}$ is the set of symmetric matrices.

These types of nodewise estimators have gained popularity as they require less restrictive incoherence conditions to attain model selection consistency and often perform better in practice in the fully observed case. They have not, however, been as well studied when used with indefinite input.

4.2 Data models

We test these methods on two models that result in indefinite covariance estimators, the non-separable Kronecker sum model from Rudelson and Zhou (2017) and the missing data graphical model described in Loh and Wainwright (2015). In the main paper we focus on the missing data model, but Section C contains a detailed description of the Kronecker sum model as well as simulation results using it.

Missing data (MD). As discussed above, Loh and Wainwright (2013); Loh and Wainwright (2015) propose an estimator for a graphical model with missing-completely-at-random observations.

Let $W\in\mathbb{R}^{n\times m}$ be a mean-zero subgaussian random matrix. Let $U\in\{0,1\}^{n\times m}$ where $U_{ij}\sim\text{Bernoulli}(\zeta_{j})$ are independent of $W$. This corresponds to entries of the $j$th column of the data matrix being observed with probability $\zeta_{j}$. Then we have an unobserved matrix $Z$ and observed matrix $X$ generated by $Z=WA^{1/2}$ and $X=U\circ X$, where $\circ$ denotes the Hadamard, or element-wise, product. Here the covariance estimate for $A$ is

[TABLE]

where $\oslash$ denotes element-wise division. As we divide off-diagonal entries by smaller values, $\hat{\Gamma}_{n}$ will not necessarily be positive semidefinite.

4.3 Simulation results

Optimization performance. Figure 1 shows the optimization performance of Algorithm 1 using nonprojected input matrices from the missing data model with both $\ell_{1}$ and nonconvex penalties (MCP). The top two panels present an “easy” scenario with a higher sampling rate, while the bottom two have a more challenging scenario with significant missing data. Blue lines report the optimization error while red lines are the statistical error.

All the plots in Figure 1 have their optimization error quickly converge to below the statistical error. These plots also suggest that our algorithm can attain linear convergence rates. We find that the algorithm consistently converges well over a range of tested scenarios.

Comparing the statistical error of the top two plots, we see that MCP achieves significantly lower error for the easier scenario. But in the bottom two plots, where there is more missing data, it struggles relative to the $\ell_{1}$ penalty. This is a common trend through our simulations, as the performance of estimators using MCP degrades as missingness increases while the $\ell_{1}$-penalized versions are more robust.

Method comparisons. Figure 2 demonstrates the statistical performance along the full regularization path. Across the panels from left to right, the sampling rate decreases and therefore the magnitude of the most negative eigenvalue increases (see Table 4).

In terms of Frobenius error, both projected methods and the nonprojected estimator with the $\ell_{1}$ penalty get slightly worse across panels, but the nodewise regression and the nonprojected MCP estimator react much more negatively to more indefinite input. The nodewise regression in particular goes from being among the best to among the worst estimators as the sampling rate decreases.

Comparing the projected and nonprojected curves in Figure 2, we see that the optimal value of $\lambda$, as well as the range of optimal values, shrinks for the projected method as the sampling rate decreases. This pattern is consistently repeated across models and scenarios, likely because the $\ell_{\infty}$ projection is shrinking the off-diagonal entries of the input matrix. We find that the nonprojected graphical Lasso performs slightly better than the projected version when used with the $\ell_{1}$ penalty, likely due to the information lost in this shrinkage.

Figure 2 also shows how these methods perform in terms of model selection. We can see that the nonconvex penalties perform essentially identically to their $\ell_{1}$ penalized counterparts. In particular, the degradation of the nonprojected MCP estimator in terms of norm error does not seem to affect its model selection performance. The nodewise regression, however, still demonstrates this pattern, as its model selection performance degrades across the panels. For scenarios with more missing data, the nonprojected estimators seem to be easier to tune, maintaining a wider range of $\lambda$ values where they perform near-optimally. In Section C of the supplement we perform similar experiments in a variety of different noise and model settings.

Sensitivity to $R$. Figure 3 demonstrates the sensitivity of the nonprojected estimators to the choice of $R$, the size of the side constraint. We can see that all these methods are sensitive to the choice of $R$ for small values of $\lambda$ in terms of norm error. None of the methods are sensitive in terms of model selection.

The nonprojected graphical Lasso with MCP is the most sensitive to $R$ and is also sensitive for larger choices of $\lambda$, which is important since it never reaches its oracle minimum norm errors when $R$ is chosen to be larger than the oracle. The nonprojected graphical Lasso with $\ell_{1}$ and the projected graphical Lasso with MCP both still achieve the same best-case performance when $R$ is misspecified, though tuning $\lambda$ becomes more difficult.

The nodewise regression results are also plotted here. Here $R$ is the $\ell_{1}$ side constraint level in (5). For smaller values of $\lambda$ the nodewise estimator levels off, corresponding to when the side constraint becomes active over the penalty. Different values of $R$ change when this occurs and, if $R$ is chosen large enough, do not significantly affect ideal performance. Note that these use a stronger oracle that knows each column-wise $\ell_{1}$ norm, but do show that this method can be improved with careful tuning.

5 Senate voting analysis

Based on the missing data model from Section 4.2, we estimate the conditional dependence graph among senators using the ADMM algorithm from Section 3. The dataset includes voting records from the United States Senate during the 112th Congress (2011-2013). We drop senators who serve partial terms and unanimous votes, resulting in a dataset of voting records for 99 senators over 426 votes. Appendix D contains further details regarding data processing and the methods used as well as additional analysis.

Missing values in this data correspond to votes that are missed by senators and consist of roughly 2.6% of total votes. Note that only 109 of the votes are fully observed, so some type of correction or imputation should be used instead of omitting rows.

A major story at this time was the rise of the tea party movement in the Republican party. Across the US government tea party challengers rose to prominence. Though it was not an official party, politicians associated with the tea party movement tended to be more conservative and less likely to compromise than establishment Republicans, leading to a particularly politically polarized period of government.

Figure 4 plots the estimated graph among senators. As expected the distinction between Republicans and Democrats is stark. Both independent senators caucus with the Democrats, so as expected they are part of the Democratic component of the graph.

We identify senators who were present at the inaugural meeting of the unofficial Senate Tea Party Caucus as well as those elected in 2010 with significant tea party support.111The marked tea party senators are Marco Rubio, Mike Lee, Jerry Moran, Jim DeMint, Rand Paul, Ron Johnson, and Pat Toomey. These senators are colored in black, and we can see that within the Republican party they are clustered together.

In Figure 4(a) we can see that the sole connection between parties runs through the tea party (Rand Paul) and Jeff Merkley, a Democratic senator. This may be surprising, as Rand Paul is one of the most conservative senators and Merkley one of the most liberal. Paul is, however, regarded as a relatively libertarian conservative. So though he is extremely conservative in some dimensions, he may share liberal views with Merkeley on others.

Figure 4(b) plots the same graph estimated at a lower penalization level. The Republicans who have cross-party connections include some of both the most conservative (Paul) and the most moderate (Thad Cochran, Lisa Murkowski).222Here we are measuring ideology by NOMINATE, a standard method in political science for assessing a representative’s position on the political spectrum (Poole, 2005). See Appendix D for more details. On the Democratic side the cross-connected senators also include both the most liberal (Sanders, Merkley, Tom Udall) and relatively moderate (Claire McCaskill). As expected, moderates are among those most connected opposing party, but this shows that the most extreme members of a party can also be linked to the opposing party. Appendix D discusses these cross-party links in more detail.

Figure 4(c) shows the Republican subgraph from Figure 4(a). Here we can identify other senators who are closely associated with the tea party. In particular, two nodes near the tea party cluster are marked ‘H’ and ‘C,’ corresponding to Senators Orrin Hatch and Tom Coburn. Both have been linked to the tea party in the media, either as candidates supported by it or as being supportive of the movement.

It is also of interest that one marked senator is not clustered with the others, Jerry Moran. This suggests that he is not as closely connected to the tea party movement as the others we have identified.

6 Summary and discussion

In this paper, we study the estimation of sparse precision matrices from noisy and missing data. To close an existing algorithmic gap, we propose an ADMM algorithm that allows for fast optimization of the side-constrained graphical Lasso, which is needed to implement the graphical Lasso with either indefinite input and/or nonconvex penalties. We investigate its convergence properties and compare its performance with other methods that handle the indefinite sample covariance matrices that arise with dirty data.

We find that methods with nonconvex penalties are quite sensitive to the indefiniteness of the input covariance estimate, and are particularly sensitive to the magnitude of its negative eigenvalues. They may have better existing theoretical guarantees, but in practice we find that with nontrivial missingness or noise they perform worst than or, at best, recover the performance of their $\ell_{1}$-normalized counterparts. The nonconvex methods can outperform the $\ell_{1}$-penalized ones when there is a small amount of missingness or noise, but in these cases we often find the nodewise estimator to perform best.

In difficult settings with significant noise or missingness, the most robust and efficient method seems to be using the graphical Lasso with nonprojected input and an $\ell_{1}$ penalty. As the application becomes easier – with more observations or less missing data – the nodewise estimator becomes more competitive, just as it is understood to be with fully observed data.

The projected graphical Lasso estimator with an $\ell_{1}$ penalty seems to be slightly worse than its nonprojected counterpart. Projection does, however, allow for the use of nonconvex penalties in more difficult settings without the large degradation in performance we have observed. This may be desired in some scenarios when the nonzero off-diagonal precision matrix entries are expected to be large.

Finally, we also use this new algorithm to estimate conditional dependence graphs among US senators using voting records data. We identify several interesting patterns in these graphs, especially regarding the rise of the tea party movement and cross-party connections between senators.

Acknowledgements

The research is supported in part by the NSF under grants DMS-1316731 and DMS-1830247.

Appendix A Auxiliary Results

A.1 Nonconvex penalties

The nonconvex penalties we will focus on are the SCAD and MCP functions, introduced in Fan and Li (2001) and Zhang (2010), respectively. Following Loh and Wainwright (2015), we make the following assumptions regarding the (univariate) penalty function $g_{\lambda}\colon\mathbb{R}\to\mathbb{R}$.

1. (i)

$g_{\lambda}(0)=0$ and $g_{\lambda}(t)=g_{\lambda}(-t)$. 2. (ii)

$g_{\lambda}(w)$ is nondecreasing for $w>=0$. 3. (iii)

$g_{\lambda}(w)/w$ is nonincreasing for $w>0$. 4. (iv)

$g_{\lambda}^{\prime}(w)$ exists for all $w\neq 0$ and $\lim_{w\to 0^{+}}g_{\lambda}^{\prime}(w)=\lambda$. 5. (v)

$g_{\lambda}$ is weakly convex, i.e. there exists $\mu>0$ such that $g_{\lambda}(w)+(\mu/2)w^{2}$ is convex.

Note that Loh and Wainwright (2017) show stronger model selection results under the following additional assumption.

1. (vi)

There exists a constant $\gamma<\infty$ such that $g_{\lambda}^{\prime}(w)=0$ for all $w>\gamma\lambda$.

This excludes the $\ell_{1}$ penalty, but is satisfied by the nonconvex penalties we consider.

The SCAD penalty takes the form

[TABLE]

for some parameter $a>2$. Note that this penalty is weakly convex with constant $\mu=1/(a-1)$.

The MCP penalty has the form

[TABLE]

for some parameter $a>0$. This penalty is weakly convex with $\mu=1/a$.

A.2 Derivation of Algorithm 1

Recall that we can rewrite the objective as

[TABLE]

where $\mathcal{X}_{R}=\{\Theta:\Theta\succeq 0,\lVert\Theta\rVert_{2}\leq R\}$ and $\mathds{1}_{\mathcal{X}}(\Theta)=0\text{ if }\Theta\in\mathcal{X}\text{ and }\infty\text{ otherwise}$.

We then introduce an auxiliary optimization variable $V\in\mathbb{R}^{m\times m}$ and reformulate the problem as

[TABLE]

For a penalty parameter $\rho>0$ and Lagrange multiplier $\Lambda\in\mathbb{R}^{m\times m}$, we consider the augmented Lagrangian

[TABLE]

The ADMM algorithm is then, given current iterates $\Theta^{k}$, $V^{k}$, and $\Lambda^{k}$,

[TABLE]

Considering the $V$-subproblem, we can show that the minimization problem in (10) is equivalent to

[TABLE]

Which is a prox operator of $g_{\lambda}/\rho$. Let $W=\frac{\rho\Theta^{k}+\Lambda^{k}}{\rho}$ and $\nu=1/\rho$. If $g_{\lambda}$ is the $\ell_{1}$ penalty then these updates simply soft-threshold the elements of $W$ at level $\lambda/\rho$. For SCAD, these updates have the element-wise form

[TABLE]

While for MCP the updates are

[TABLE]

See Loh and Wainwright (2015) for the derivations of these updates.

For the $\Theta$-subproblem, we can similarly show that (11) is equivalent to

[TABLE]

For any matrix $A$ with corresponding eigendecomposition $A=RMR^{T}$ let us define the operator

[TABLE]

whose solution is derived in Section A.3. Then the solution to (11) is $T_{\rho}((\rho V^{k+1}-\hat{\Gamma}_{n}-\Lambda^{k})/\rho)$.

Using these results, the algorithm in (10)-(12) becomes

[TABLE]

A.3 Solution of $T_{\rho}$

Recall that in (16) we define

[TABLE]

Let $\Theta=WDW^{T}$ and $A=UMU^{T}$ be the eigen-decompositions of the optimization variable and $A$. Then, similar to the derivation in Guo and Zhang (2017), we can rewrite this problem as

[TABLE]

The final line is since, if we denote $O(m)$ to be the set of $m\times m$ orthonormal matrices,

[TABLE]

Which holds with equality when $W=U$. Note that the last equality here is from Theorem 14.3.2 of Farrell (1985).

We therefore get that $T_{\rho}(A)=U\tilde{D}U^{T}$ where

[TABLE]

We can see that this is separable by element. Let

[TABLE]

So $\tilde{D}_{ii}=\operatorname*{arg\,min}_{d}q(d;M_{ii})+\mathds{1}(0\leq d\leq R)$. Ignoring the constraints in the indicator function for now, we can set the derivative of $q$ equal to zero to get that

[TABLE]

Which we can solve with the quadratic formula to show that $q(d;M_{ii})$ has a unique minimizer over $d>0$ at

[TABLE]

Adding $\mathds{1}(0\leq d\leq R)$ back and noting that $q(d;M_{ii})$ is strictly convex over $d>0$, we get that we simply need to truncate this value at $R$. Therefore we get that

[TABLE]

A.4 Proof of Proposition 1

Proof.

The optimization problem (4) is equivalently

[TABLE]

where $f_{1}(\Theta)=\operatorname*{tr}(\hat{\Gamma}_{n}\Theta)-\log\det(\Theta)+\mathds{1}_{\mathcal{X}_{R}}(V)$, $f_{2}(V)=g_{\lambda}(V)$, $A=-I_{m^{2}}$, and $B=I_{m^{2}}$.

Boyd et al. (2010) show that if $f_{1}$ and $f_{2}$ are proper convex functions and if (18) is solveable then ADMM converges in terms of the objective value $\phi(\Theta^{t},V^{t})\to\phi^{*}$ and dual variable $\Lambda^{t}\to\Lambda^{*}$. Bertsekas and Tsitsiklis (1989, Proposition 4.2) and Mota et al. (2011) show that if in addition $A$ and $B$ have full column rank then we get convergence of the primal iterates $\Theta^{t}\to\Theta^{*}$ and $V^{t}\to V^{*}$, where $(\Theta^{*},V^{*})$ is the solution to (18). ∎

A.5 Proof of Proposition 2

Before we prove Proposition 2, we first define directional derivatives and stationary points.

Definition**.**

The directional derivative of a lower semi-continuous function $h$ at $\Theta$ in the direction $\Delta$ is

[TABLE]

Note that we allow $h^{\prime}(\Theta;\Delta)=+\infty$. We say that $\Theta$ is a stationary point of $h$ if it satisfies the first-order necessary conditions to be a local extrema, i.e.

[TABLE]

Note that this coincides with the definition of stationary point used in Loh and Wainwright (2017), though they use slightly different notation. Also note that $h^{\prime}(\Theta;\Delta)=\langle\nabla h(\Theta),\Delta\rangle$ when $h$ is continuously differentiable.

Proof.

From the first-order necessary conditions of the subproblems (10)-(11), we get that, for all $\Delta\in\mathbb{R}^{m\times m}$,

[TABLE]

And recall that

[TABLE]

We can rewrite (19)-(20) as

[TABLE]

Now consider a fixed point $(\Theta^{*},V^{*},\Lambda^{*})$ and consider (21)-(23) evaluated at this limit point. From (23) we get that $\Theta^{*}=V^{*}$. This combined with (21) gives us that, for all $\Delta\in\mathbb{R}^{m\times m}$,

[TABLE]

Finally, (22) gives us that

[TABLE]

Using the above and recalling the objective $f$ as defined in (4), we get that, for all $\Delta\in\mathbb{R}^{m\times m}$,

[TABLE]

So $\Theta^{*}$ is a stationary point of $f$ by definition. ∎

A.6 Comparison to Guo and Zhang (2017)

Guo and Zhang (2017) study the problem of condition number-constrained precision matrix estimation, where they consider the estimator

[TABLE]

Note that this is quite similar to the estimators we consider in (2), as they simply replace the maximum eigenvalue constraint with a constraint on the ratio of the maximum to minimum eigenvalues.

However, they do not study the application of their estimator to cases with indefinite input or its performance in noisy and missing data situations. In particular, constraining the condition number does not necessarily guarantee that the graphical Lasso objective (1) will be lower bounded, especially when using nonconvex penalties.

As a simple example, consider the case with an input matrix and iterates

[TABLE]

In this case the objective is

[TABLE]

which is unbounded below as $t$ grows even though the condition numbers of the iterates are constant.

More generally, whenever $\hat{\Gamma}_{n}\in\mathbb{R}^{m\times m}$ has eigenvalues $\sigma_{1},\dots,\sigma_{m}$, where $\sigma_{1}\geq\dots\geq\sigma_{m_{1}}\geq 0$ and $0>\sigma_{m_{1}+1}\geq\dots\geq\sigma_{m}$. Denote $S_{1}=\sum_{i=1}^{m_{1}}\sigma_{i}$ and $S_{2}=\sum_{i=m_{1}+1}^{m}-\sigma_{i}$. Let $VDV^{T}=\hat{\Gamma}_{n}$ be the eigendecomposition of the covariance estimate. Then for some condition number bound $\kappa$, we can consider iterates of the form $\Theta^{t}=tVMV^{T}$, where $M$ is a diagonal matrix with entries

[TABLE]

Which we note has a condition number of $\kappa$. Then we can see that the objective becomes

[TABLE]

So if $\kappa>S_{1}/S_{2}$ then this objective is still unbounded below.

Using a spectral norm bound $\lVert\Theta\rVert_{2}\leq R$ as the side constraint with a indefinite input guarantees a lower bound on the graphical Lasso objective regardless of the choice of $R$ and is therefore a more natural side constraint to use.

Appendix B ADMM for general side constraints

In this section we develop an ADMM algorithm for general side constraints, i.e. the following variant of (2).333Note that we switch the notation of the side constraint function from $\rho$ to $h$ to avoid confusion with the ADMM penalty parameter $\rho$.

[TABLE]

This algorithm has the same convergence guarantees as Algorithm 1, but in practice we find that Algorithm 1 converges faster and more consistently when the spectral norm side constraint is used.

B.1 Derivation

We first rewrite the objective as

[TABLE]

where $\mathcal{X}_{h,R}=\{\Theta:\Theta\succeq 0,h(\Theta)\leq R\}$ and

[TABLE]

We can then introduce auxiliary optimization variables $V_{1},V_{2}\in\mathbb{R}^{m\times m}$ and reformulate the optimization problem as

[TABLE]

For a penalty parameter $\rho>0$ and Lagrange multiplier matrices $\Lambda_{1},\Lambda_{2}\in\mathbb{R}^{m\times m}$, we consider the augmented Lagrangian of this problem

[TABLE]

The ADMM algorithm is then, given current iterates $\Theta^{k}$, $V_{1}^{k}$, $V_{2}^{k}$, $\Lambda_{1}^{k}$, and $\Lambda_{2}^{k}$,

[TABLE]

Considering the $V_{1}$-subproblem, we can show that the minimization problem in (27) is equivalent to

[TABLE]

Which is a prox operator of $g_{\lambda}/\rho$. These have the same form as described in Section A.2.

For the $V_{2}$-subproblem, we similarly see that (28) is equivalent to

[TABLE]

which is equivalent to the projection operator

[TABLE]

Note that if directly projecting onto $\mathcal{X}_{h,R}$ does not have an closed-form solution, we can perform this step using Dykstra’s alternating projection algorithm.

Finally, for the $\Theta$-subproblem, we can again show that (29) is equivalent to

[TABLE]

Let us define the operator

[TABLE]

whose solution is derived in Section A.3 if we set $R=\infty$. Then the solution to (29) is $\tilde{T}_{\rho}((\rho V_{1}^{k+1}+\rho V_{2}^{k+1}-\hat{\Gamma}_{n}-\Lambda_{1}^{k}-\Lambda_{2}^{k})/(2\rho))$.

Using these results, the algorithm in (27)-(31) becomes

[TABLE]

B.2 Convergence

Analogues to Propositions 1 and 2 can also be shown for this algorithm using similar methods. To do this, we first note that we can rewrite the optimization problem (25) as

[TABLE]

where

[TABLE]

and

[TABLE]

This results in the following augmented Lagrangian that is equivalent to (26).

[TABLE]

Even though we present our algorithm as a three-block ADMM in Section B.1, this formulation makes it clear that we are using a two-block splitting scheme where (27) and (28) are the separable subproblems of the $V$-step.

Showing similar convergence results to Propositions 1 and 2 can then be done using the same techniques as in Sections A.4 and A.5.

Appendix C Additional simulation results

C.1 Tuning parameter selection

Note that in practice tuning parameters must be selected for all these methods. In particular, we must tune $\lambda$ and possibly the side-constraint $R$. Note that one often has a reasonable prior for the magnitude of the spectral norm of the true precision matrix, so if that is the case a multiple of that can often be used to choose $R$. Also, as noted in Section 4.3, when using the $\ell_{1}$ penalty the choice of $R$ primarily affects how difficult tuning $\lambda$ will be. Though it is important to tune correctly when using nonconvex penalties, we do not recommend those methods when there is significant missing data. Therefore we will focus on tuning $\lambda$ here, though the same methods can be used to choose $R$ as well.

Two possible methods are to use cross-validation or a modified BIC criterion. Though the particular implementation of both of these will depend on the data model that is being used, as these methods can be applied to any method that generates an indefinite initial estimate of the covariance.

For the missing data case we can follow Städler and Bühlmann (2012), which uses the same data model. Recall the notation in Section 4.2, where $X_{ij}$ denotes the $i$th value of variable $j$ and $U_{ij}$ tracks if that value is observed. Here, we define the observed log-likelihood of an observation $X_{i}$ given a precision matrix estimate $\hat{\Theta}$ as

[TABLE]

where $X_{i,U_{i}}$ is the vector of values that are observed for observation $i$, $\hat{\Sigma}=\hat{\Theta}^{-1}$, and $\phi$ is the multivariate normal density. The BIC criterion, which we minimize, is therefore

[TABLE]

To cross-validate, we can divide the data into $V$ folds, where the $v$th fold contains indices $N_{v}$. The cross-validation score, which we maximize, is therefore

[TABLE]

where $\hat{\Sigma}_{-v}=\hat{\Theta}_{-v}^{-1}$ and $\hat{\Theta}_{-v}$ is the estimate based on the sample omitting the observations in $N_{v}$.

Figure 5 presents an example of parameter tuning on a simulated scenario. We see that both BIC and CV select slightly higher-than-optimal levels of penalization in terms of model selection, but that selected model still achieves fairly good model selection.

C.2 Kronecker sum (KS) model

Park et al. (2017) present a graphical model with additive noise that is dependent across observations. This noise structure was first studied in the regression setting in Rudelson and Zhou (2017) with a Kronecker sum covariance.

Let $W_{1},W_{2}\in\mathbb{R}^{n\times m}$ be independent mean-zero subgaussian random matrices. The data matrix is generated as $X=W_{1}A^{1/2}+B^{1/2}W_{2}\sim\mathcal{M}_{n,m}(0,A\oplus B)$, where $\mathcal{M}_{n,m}$ is the matrix variate normal distribution and for covariance matrices $A\in\mathbb{R}^{m\times m}$ and $B\in\mathbb{R}^{n\times n}$. Note that $A\oplus B=A\otimes I_{n}+I_{m}\otimes B$, where $\otimes$ denotes the Kronecker product. Here $X_{0}=W_{1}A^{1/2}$ contains the signal and has independent rows, while $W=B^{1/2}W_{2}$ is the noise matrix with independent columns but dependent rows. We are interested in estimating the signal precision matrix $\Theta=A^{-1}$, which has sparse off-diagonal entries. For our simulations, we normalize $B$ so that $\operatorname*{tr}(B)=n\tau_{B}$, where $\tau_{B}$ is a measure of the noise level. Then the initial covariance estimate for $A$ is given by

[TABLE]

as shown in Rudelson and Zhou (2017). Note that, in this model, $\hat{\Gamma}_{n}$ is guaranteed to not be positive semidefinite when $m>n$, as $X^{T}X$ will have zero eigenvalues.

C.3 Covariance models

We look into three different models from which $A$ Let $\Omega=A^{-1}=(\omega_{ij})$. We consider simulation settings using the following covariance models, which are also used in Zhou (2014).

• •

AR1($r$): The covariance matrix is of the form $A=(r^{\lvert i-j\rvert})_{ij}$.

• •

Star-Block (SB): Here the covariance matrix is block-diagonal, where each block’s precision matrix corresponds to a star-structured graph with $A_{ii}=1$. For the corresponding edge set $E$, then $A_{ij}=r$ if $(i,j)\in E$ and $A_{ij}=r^{2}$ otherwise.

• •

Erdos-Renyi random graph (ER): We initialize $\Omega=0.25I$ then randomly select $d$ edges. For each selected edge $(i,j)$, we randomly choose $w\in[0.6,0.8]$ and update $\omega_{ij}=\omega_{ji}\to\omega_{ij}-w$ and $\omega_{ii}\to\omega_{ii}+w$, $\omega_{jj}\to\omega_{jj}+w$.

C.4 Optimization performance

Figure 6 shows the convergence behavior for several initializations in terms of objective value. Our algorithm seems to attain a linear convergence rate in terms of the objective values even with a nonconvex penalty regardless of the initialization. We find that the algorithm consistently converges well over a range of tested scenarios.

Comparison to gradient descent. Figure 7 compares the optimization performance of our ADMM algorithm to gradient descent. Note that since proximal gradient descent is difficult to do in this setting, requiring an interior optimization step, we use a heuristic version similar to that suggested by Agarwal et al. (2012) that does the proximal gradient step ignoring the side-constraint then projects back to the side-constraint space. Note that since $\rho$ in ADMM is roughly equivalent to the inverse step size in gradient descent, we compare for difference values of $\rho$. These methods also take roughly the same computational time per iteration, as they are both dominated by either an eigenvalue decomposition or matrix inversion.

We can see that for large enough values of $\rho$, these methods are nearly identical. Although there is no known theoretical guarantee of convergence, it seems that this heuristic gradient descent still convergence well for small enough step sizes.

But for smaller values, i.e. larger step sizes, ADMM still performs well and obtains faster convergence rates while gradient descent is unstable and inconsistent. This combined with the convergence guarantee of ADMM leads us to recommend this algorithm.

C.5 Penalty nonconvexity and $R$

Suppose $g_{\lambda}$ is $\mu$-weakly convex and $R\leq\sqrt{\frac{2}{\mu}}$. Then, as shown in Lemma 6 of Loh and Wainwright (2017), the overall objective function is strictly convex over the feasible set, and Proposition 2 therefore shows that any limiting point of ADMM algorithm corresponds to the unique global optimum of the objective. However, this choice of $R$ radius on the $\lVert\cdot\rVert_{2}$ side constraint is quite restrictive. In particular, since we also require $R\geq\lVert\Theta^{*}\rVert_{2}$ we therefore need to choose large values of $a$ in the SCAD or MCP penalties to make $\mu$ small enough, which means in practice we simply recover the performance of the $\ell_{1}$ penalized methods. Though Loh and Wainwright (2017) show statistical properties for when the parameters are chosen satisfying this condition, in practice we can often do better by allowing the objective to be nonconvex even though no global optimum will exist.

Once we relax this condition ($R>\sqrt{2/\mu}$), the objective becomes nonconvex, and Proposition 2 simply shows that any limiting point of our ADMM algorithm will be a stationary point of the objective. In our simulations, we generally set $\mu$ and $R$ such that this condition is violated, and yet we show that our algorithm still results in good estimators. In fact, Figure 8 demonstrates how, in practice, choosing $\mu$ such that this condition is met tends to eliminate the advantages that nonconvex penalties provide. Here the choice of $a=8$ is the only one that satisfies the condition, and this choice has identical performance as the $\ell_{1}$ penalty. Using a smaller value of $a$ violates this condition but allows the estimator to take advantage of the unbiasedness of the penalty, resulting in better performance in this setting.

Note that for both of these cases, our ADMM algorithm provides a new feasible method of implementing estimation of this type of side-constrained graphical Lasso objective. This consideration is related to tuning, where satisfying the $(R,\mu)$ condition allows the support recovery without incoherence statistical results of Loh and Wainwright (2017) but in practice results in suboptimal performance, as the nonconvex penalties have to be chosen such that they lose their unbiased advantage over the $\ell_{1}$ penalty.

C.6 Method comparisons

Tables 1-3 present more detailed comparison based on the models from the Kronecker sum (KS) and the missing data (MD) models. We compared performance in terms of relative Frobenius and nuclear norm to the true precision matrix, as well as false positive rate plus false negative rate (FPRFNR). The Kronecker sum results are reported for two sample sizes and two values of the noise parameter $\tau_{B}$, while the missing data results are reported for two covariance models and three settings of the sample size and sampling rate $\zeta$.444Note that in the initial covariance estimator for the missing data model the effective sample size for estimating an off-diagonal element of the covariance is $n\zeta^{2}$; four settings are designed to keep this effective sample size roughly constant while changing the sampling rate $\zeta$. The effective sample sizes for the $n=80$, $n=130$, and $n=250$ settings are $64.8$, $63.7$, and $62.5$, respectively.

Comparing the projected and nonprojected methods, we see that these two methods are fairly competitive. In terms of model selection, the nonprojected methods tend to perform similarly or better than the projected methods. This improvement is particularly evident in the $n=80$ settings in Table 1. If we focus on the methods using the $\ell_{1}$ penalty, the nonprojected method performs at least similarly and sometimes significantly better than the projected method in terms of norm error. The lower sampling rate regime in Tables 2 and 3 shows this trend as well. Overall these results suggest a small but sometimes significant advantage for the nonprojected methods, supporting the idea that the projected methods pay a cost in terms of efficiency due to the loss of information in the projection.

There is no significant difference in model selection between MCP and the $\ell_{1}$ penalty. In fact, the different penalties perform almost identically across scenarios regardless of the $\ell_{\infty}$-projection step. Intuitively, the primary benefit of nonconvex penalties is their ability to more accurately estimate large entries, which are easy for the estimators to select.

In terms of norm error, however, there are significant differences depending on the indefiniteness of the optimization problem. Table 4 reports some statistics on the eigenspectrum of the input matrix. Nonprojected methods with MCP tends to perform relatively better than its $\ell_{1}$ counterpart if the input matrix is close to the positive semidefinite space. Simulation results from the missing data model Tables 2 and 3 further support this relationship between the most negative eigenvalue and the relative performance. Here we see how the MCP nonprojected estimator goes from being significantly better than its $\ell_{1}$ counterpart in terms of Frobenius error in the $\zeta=0.9$ case to significantly worse when $\zeta=0.5$. In the projected case, which projects away this indefinite issue, the MCP estimator consistently outperforms its $\ell_{1}$ counterpart in terms of Frobenius error.

The nonconvexity of the penalty interacts poorly with indefiniteness of the input matrix. When the $\ell_{1}$ penalty is used, it is better able to “control” the indefiniteness of the input due to its linear scaling, resulting in better norm error performance. The nonconvex penalty’s inability to resolve the indefiniteness issue results in a degradation of its relative performance as the input matrix becomes more indefinite.

Turning to the nodewise estimator, we see similar patterns. Again referring to Table 4, it seems that the relative performance of the nodewise estimator varies significantly with the indefiniteness of the input matrix. When the input matrix is closer to positive semidefinite, such as the $n=160$ situations in Table 1 or the $\zeta=0.9$ cases in Tables 2 and 3, it performs comparably in terms of model selection and significantly better in terms of norm error. But when the input matrix is very indefinite, such as the $\zeta=0.5$ cases in Tables 2 and 3 its relative performance quickly degrades.

Figure 9 demonstrates the patterns that we observed in Figures 2. Again, we vary the sampling rate $\zeta$ and fix the effective sample size for estimating off-diagonal entries of the covariance matrix ($n\zeta^{2}$), so the $\ell_{\infty}$ rate of $\hat{\Gamma}_{n}$ is kept constant. As the sampling rate decreases, the magnitude of the most negative eigenvalue in the covariance estimate increases, which we can see negatively affects the relative performance of the nodewise and nonprojected MCP methods. These methods are the best for high sampling rates but the worst when there is more missing data. The other methods are not as sensitive.

We can also see how the size of the off-diagonal entries in the precision matrix affect the potential benefits of the nonconvex penalties. In the top panel, which has small off-diagonal entries, the MCP estimators are consistently worse. But in the bottom panel, which has larger off-diagonal entries, the nonconvex penalties have better Frobenius norm performance when the sampling rate is high, though this advantage goes away as the sampling rate drops.

C.7 Comparison of side constraints

Here we compare using the operator norm side-constraint in (2) to the $\ell_{1}$-side constrained version considered in Loh and Wainwright (2015). Note that theoretically Loh and Wainwright (2017) show that under certain conditions the former can attain model selection without incoherence and spectral norm convergence under the scaling $n>d^{2}\log p$ (where $d$ is the maximal node degree), which has not been shown with the latter.

Figure 10 shows the performance in terms of relative norm error for various missing data model scenarios. For large values of $\lambda$ we can see that the two estimators are identical, as the side-constraints are not active.

As the penalty $\lambda$ shrinks, when the $\ell_{1}$ side constraint is used, selecting $R$ is akin to performing selection on the minimum value of $\lambda$ to use, since only one of the penalty and side-constraint can be active at a time. Following the green lines, we can see that the regularization from the operator norm side-constraint can improve results. This means that $\ell_{1}$ side-constrained estimator misses this additional improvement. At worst, for larger values of $R$ the operator norm-constrained simply has identical performance to the $\ell_{1}$-constained version as long as $\lambda$ is appropriately selected.

Appendix D Additional data analysis

As discussed in Section 5, we collected voting records data from the Senate during the 112th US Congress, which was from January 3, 2011 to January 3, 2013. This data is part of the public record and open-source code to download and process the data can be found at https://github.com/unitedstates/congress.

Due to changes in membership, there are data on 102 senators, which we drop three of due to serving incomplete terms. The data contains 486 votes in total. We drop votes that are unanimous or unanimous within both parties, resulting in 426 votes. Roughly 2.6% of values are missing in this data.

We use the ADMM algorithm from Section 3 to estimate the nonprojected version of (2) with an $\ell_{1}$ penalty to estimate conditional dependence graphs among senators. Since this is an exploratory analysis, we ran estimators using various levels of the penalization parameter $\lambda$ and have chosen plots to display based on the number of estimated edges and maintaining visual clarity.

For our preliminary analysis, we use a modified version of the missing data estimator as described in Zhou (2019), where bills have varying missing probabilities while we estimate the edges among senators. We also demean each vote by political party, similar to the demeaning done in Hornstein et al. (2018). See our future work for a more detailed study of this estimator and its properties.

Figure 11 plots the subgraph of senators with cross-party connections or links to those with cross-party connections from Figure 4(b). We look at the NOMINATE scores of these senators (Figure 12) to determine their positions on the political spectrum. NOMINATE is a probabilistic geometric model that places each senator in a two-dimensional space representing their ideological beliefs (Poole, 2005).

Although most of the linked senators are are either on the extremes or are moderates (see Section 5), the main exception to this is Dean Heller, who is linked to Mark Warner. This outlier connection is perhaps worth further investigation as to why they are linked.

The McCaskill-Cochran connection is unsurprising, as both are among the most moderate senators from their respective parties. The Udall-Murkowski connection is more interesting since Tom Udall is viewed as a relatively liberal Democrat while Murkowski is a moderate Republican. Murkowski connecting across the aisle is expected, but why she would be linked to Udall is unknown.

Paul-Sanders and Paul-Merkley are both interesting connections since they are between left-wing and right-wing senators. In fact, these three senators form a triangle of positively correlated edges. This could be because they are all ideological extremists and therefore vote together against bills with moderate support. Or this could be explained by the fact that all three senators are liberal on NOMINATE’s second dimension, and therefore vote together on those issues.

Figure 13 shows the estimated Republican and Democratic subgraphs with edge weights and directions as determined by the partial correlation esimates, while Table 5 lists the top ten edges by strongest positive partial correlation.

We can see that for both parties most of the strongest links are between senators from the same state, which is expected since these senators must appeal to the same constituents and therefore are likely to vote similarly. Of the 33 states where both senators are in the same party, 19 of the pairs are linked and 13 of those appear in the top connections in Table 5.

Looking at the exceptions closely, we can often identify why these links are particularly strong. For instance, DeMint-Lee and Coburn-Johnson are both between senators we have identified as being tea party-linked (Coburn due to his close proximity and support in the media). Also, Collins-Brown is between two moderate Republicans who are both from Democratic New England states.

Table 6 also contains the most negatively estimated partial correlations. Note that the negative relationships we find are in general of much smaller magnitude than the positive correlations. For instance, the most negative correlation (Akaka-McCaskill at -0.13) has roughly same magnitude as the 34th-most positive correlation (DeMint-Johnson).

D.1 Nodewise regression

We also demonstrate the usage of nodewise regression in Figure 14. Since the sampling rate for this dataset is fairly high, we expect nodewise regression to perform well, and the results are overall similar to those in Figure 4(b). Of the four cross-party links, three (Cochran-McCaskill, Sanders-Paul, and Warner-Heller) are also present in the previously estimated graph. The link between Democrat Mark Begich and Lisa Murkowski is new, a natural one since they both represent Alaska in the Senate.

Figure 14(b) exhibits similar patterns to those identified in Figure 4(c). Hatch and Coburn, marked as ‘H’ and ‘C,’ still appear to be closely connected to the tea party cluster. Moran is also still disconnected from the rest of the tea party despite attending the inaugural meeting of the Senate Tea Party Caucus.