A. Attractive Gaussian Markov Random Fields

We denote an undirected graph by G=(V,E), where V = {1,…,p} is the vertex set, and E is the edge set. Let y=(y1,…,yp) be a zero-mean p-dimensional random vector following y∼N(0,Σ). We associate the random variables y₁,...,y_p with the vertex set V. Then, the random vector yforms a Gaussian Markov random field with respect to a graph G=(V,E), where

View Source\begin{equation*}\begin{array}{l} {\Theta _{ij}} \ne 0 \Leftrightarrow (i,j) \in E\forall i \ne j, \\ {\Theta _{ij}} = 0 \Leftrightarrow {y_i} \perp {y_j}\left| {{{\mathbf{y}}_{[p]\backslash \{ i,j\} }},} \right. \end{array} \tag{1}\end{equation*}

In this paper, we focus on the attractive Gaussian Markov random fields, where the precision matrix Θ is a symmetric Mmatrix [24], i.e., [Θij]≤0, for any i ≠ j. In other words, all the partial correlations are nonnegative in attractive Gaussian Markov random fields.

We aim to estimate a sparse precision matrix under attractive Gaussian Markov random fields given independent and identically distributed observations y⁽¹⁾,…,y⁽ⁿ⁾ ∈ ℝ^p. Let S=1n∑ni=1y(i)(y(i))⊤ be the sample covariance matrix. The sparse precision matrix estimation under the attractive Gaussian Markov random fields can be formulated as the ℓ₁-norm regularized maximum likelihood estimation under the MTP₂ constraints,

View Source\begin{equation*}{X^ \star } = \arg \mathop {\min }\limits_{{\mathbf{X}} \in {\mathcal{M}^p}} - \log \det (X) + \operatorname{tr} (XS) + \lambda \sum\limits_{i \ne j} {\left| {{X_{ij}}} \right|} ,\tag{2}\end{equation*}

Mp:={X∈Sp++∣Xij≤0,∀i≠j}.(3)

View Source\begin{equation*}{\mathcal{M}^p}: = \left\{ {X \in \mathbb{S}_{ + + }^p\mid {X_{ij}} \leq 0,\forall i \ne j} \right\}.\tag{3}\end{equation*}

We propose a second-order algorithm in Section III to solve this problem.

B. Related Work

Sparse precision matrix estimation under Gaussian Markov random fields has been extensively studied in the literature. One popular approach is the ℓ₁-norm regularized Gaussian maximum likelihood estimation [11], [12], [13], and numerous algorithms have been proposed for solving this problem. A representative, yet not exhaustive, list of works available in the literature include: block coordinate ascent method [11], [25], Nesterov’s smooth gradient method [12], projected gradient method [26], projected quasi-Newton [27], augmented Lagrangian method [28], inexact interior point method [29], primal proximal-point with Newton-CG method [30], and Newton’s method with quadratic approximation [31], [32], [33]. However, all the methods mentioned above cannot be directly extended to estimate precision matrices in our problem because of MTP₂ constraints.

To estimate precision matrices under MTP₂ constraints, one option is using the Gaussian maximum likelihood estimator [5], [34], and the sign constraint can promote the sparsity of the solution implicitly. To solve this problem, a primal algorithm [2] and a dual algorithm [23] were proposed based on block coordinate descent. The authors in [6] proposed the ℓ₁-norm regularized Gaussian maximum likelihood estimation, and updated each column/row of the variable by solving a nonnegative quadratic program. In addition, there is growing interest in learning graphs under the generalized graph Laplacian models [6], [23], [35], where the precision matrices satisfy MTP₂ constraints. In this paper, we aim to propose an efficient algorithm for estimating sparse precision matrices under MTP₂ constraints, and establish the successful edge recovery guarantees.

A. Proposed Algorithm

To solve Problem (2), we propose a projected Newton-like algorithm. In each iteration, we partition the algorithm variables into two sets, i.e., free and restricted sets, and only update the variables in the free set while fixing the variables in the restricted set. The restricted set of variables is defined by

View Source\begin{equation*} {\mathcal{I}_k}: = \left\{ {(i,j) \in {{[p]}^2}\mid {{\left[ {{X_k}} \right]}_{ij}} = 0,{{\left[ {\nabla f\left( {{X_k}} \right)} \right]}_{ij}} < 0} \right\},\tag{4}\end{equation*}

Now we construct the projected Newton-like step as follows,

View Source\begin{equation*}{X_{k + 1}} = {\mathcal{P}_\Omega }\left( {{X_k} - {\gamma _k}{P_k}} \right),\tag{5}\end{equation*}

pveck(Pk)=Qk pvec k(∇f(Xk)),(6)

View Source\begin{equation*}{\operatorname{pvec} _k}\left( {{P_k}} \right) = {Q_k}{\text{ pvec}}{{\text{ }}_k}\left( {\nabla f\left( {{X_k}} \right)} \right),\tag{6}\end{equation*}

pveck(X):=[[X]Ick[X]Ik],(7)

View Source\begin{equation*}{\operatorname{pvec} _k}(X): = \left[ {\begin{array}{l} {{{[X]}_{\mathcal{I}_k^c}}} \\ {{{[X]}_{{\mathcal{I}_k}}}} \end{array}} \right],\tag{7}\end{equation*}

The approximate Newton direction in (6) may involve computing the inverse of the Hessian matrix with the dimension p² × p² or solving a system of linear equations of the same dimension, which is computationally challenging. In what follows, we aim to simplify the computation in (6). We construct the gradient scaling matrix Q_kas follows,

View Source\begin{equation*}{Q_k} = \left[ {\begin{array}{cc} {{{\left[ {H_k^{ - 1}} \right]}_{\mathcal{I}_k^c\mathcal{I}_k^c}}}&0 \\ 0&{{D_k}} \end{array}} \right],\end{equation*}

H−1k=Xk⊗Xk.(8)

View Source\begin{equation*}H_k^{ - 1} = {X_k} \otimes {X_k}.\tag{8}\end{equation*}

Following from the property of Kronecker product that vec(ABC) = (C^⊤⊗A)vec(B), we obtain

View Source\begin{equation*}H_k^{ - 1}\operatorname{vec} \left( {\nabla f\left( {{X_k}} \right)} \right) = \operatorname{vec} \left( {{X_k}\nabla f\left( {{X_k}} \right){X_k}} \right).\tag{9}\end{equation*}

As a result, we can get

View Source\begin{equation*}\begin{array}{l} {Q_k}{\operatorname{pvec} _k}\left( {\nabla f\left( {{X_k}} \right)} \right) = \left[ {\begin{array}{cc} {{{\left[ {H_k^{ - 1}} \right]}_{\mathcal{I}_k^c\mathcal{I}_k^c}}}&0 \\ 0&{{D_k}} \end{array}} \right]\left[ {\begin{array}{l} {{{\left[ {\nabla f\left( {{X_k}} \right)} \right]}_{\mathcal{I}_k^c}}} \\ {{{\left[ {\nabla f\left( {{X_k}} \right)} \right]}_{{\mathcal{I}_k}}}} \end{array}} \right] \\ = \left[ {\begin{array}{c} {{{\left[ {{X_k}{\mathcal{P}_{\mathcal{I}_k^c}}\left( {\nabla f\left( {{X_k}} \right)} \right){X_k}} \right]}_{\mathcal{I}_k^c}}} \\ {{D_k}{{\left[ {\nabla f\left( {{X_k}} \right)} \right]}_{{\mathcal{I}_k}}}} \end{array}} \right], \end{array} \end{equation*}

[PIck(A)]ij={Aij0if(i,j)∈Ick,otherwise.

View Source\begin{equation*}{\left[ {{\mathcal{P}_{\mathcal{I}_k^c}}(A)} \right]_{ij}} = \begin{cases} {{A_{ij}}}&{{\text{if}}\,(i,j) \in \mathcal{I}_k^c,} \\ 0&{{\text{otherwise}}{\text{.}}} \end{cases} \end{equation*}

Therefore, we have

View Source\begin{equation*}{\left[ {{P_k}} \right]_{\mathcal{I}_k^c}} = {\left[ {{X_k}{\mathcal{P}_{\mathcal{I}_k^c}}\left( {\nabla f\left( {{X_k}} \right)} \right){X_k}} \right]_{\mathcal{I}_k^c}}.\end{equation*}

In each iteration, we only update the variables in the free set, i.e., [X]Ick, and set the remaining variables in the restricted set to be zero. Finally, we can rewrite the projected Newton-like iteration in (5) as

View Source\begin{equation*}{\left[ {{X_{k + 1}}} \right]_{\mathcal{I}_k^c}} = {\mathcal{P}_\Omega }\left( {{{\left[ {{X_k}} \right]}_{\mathcal{I}_k^c}} - {\gamma _k}{{\left[ {{X_k}{\mathcal{P}_{\mathcal{I}_k^c}}\left( {\nabla f\left( {{X_k}} \right)} \right){X_k}} \right]}_{\mathcal{I}_k^c}}} \right),\end{equation*}

[Xk+1]Ik=0.

View Source\begin{equation*}{\left[ {{X_{k + 1}}} \right]_{{\mathcal{I}_k}}} = 0.\end{equation*}

Note that the Newton-type methods for solving our problem (2) are usually computationally expensive, because of the computations of the (approximate) inverse of the Hessian matrices with the dimension p² × p². However, it is observed that our constructed iterates as shown above only involves the matrix multiplication with the dimension p × p and gradient computation. The step size can be computed by the Armijo step-size rule.

Algorithm 1 Projected Newton-like method

B. Theoretical Results

In this subsection, we provide theoretical results to show that the estimator X^⋆ defined in (2) can recover the underlying graph edges correctly under irrepresentability condition.

Let S:={(i,j)∣Θij≠0} be the support set of the underlying precision matrix Θ, and d be the maximum number of nonzero elements in any row of Θ. Define

View Source\begin{equation*}{K_\Sigma }: = \mathop {\max }\limits_{i \in \{ 1, \ldots ,p\} } \sum\limits_{j = 1}^p {{\Sigma _{ij}}} ,{\text{ and }}{K_H}: = {\left| {\left\| {\left( {{{\mathbf{H}}_{\mathcal{S}\mathcal{S}}}} \right)} \right\|} \right|^{ - 1}}_\infty ,\tag{10}\end{equation*}

Assumption 1. There exists some α ∈ (0,1] such that

View Source\begin{equation*}{\left| {\left\| {{H_{{\mathcal{S}^c}\mathcal{S}}}{{\left( {{H_{\mathcal{S}\mathcal{S}}}} \right)}^{ - 1}}} \right\|} \right|_\infty } \leq 1 - \alpha .\tag{11}\end{equation*}

Assumption 2. The nonzero elements of the underlying precision matrix Θ satisfy

View Source\begin{equation*}\mathop {\min }\limits_{(i,j) \in \mathcal{S}} \left| {{\Theta _{ij}}} \right| \geq \left( {1 + \frac{\alpha }{2}} \right){K_H}\lambda .\tag{12}\end{equation*}

Assumption 1 presents the irrepresentability condition that is almost necessary for the ℓ₁-norm approaches to recover the supports correctly [14]. Assumption 2 imposes a lower bound on the minimum absolute value of nonzero elements of Θ.

Theorem 1. Suppose the sample covariance matrix satisfies ∥S−Σ∥max≤α4λ. Under Assumptions 1 and 2, if the sample size is lower bounded by n ≥ cd² logp, where c=(3(1+α2)cλKHKΣmax((1+2α)KHK2Σ,1))2, then X^⋆ obtained by solving Problem (2) with λ=cλlogpn−−−−√ can recover the support of Θ correctly, i.e.,

View Source\begin{equation*}\operatorname{supp} \left( {{X^ \star }} \right) = \operatorname{supp} ({\mathbf{\Theta }}),\tag{13}\end{equation*}

where c_λ is a positive constant.

Theorem 1 shows that the support of the minimizer X^⋆ of Problem (2) is the same with that of the underlying precision matrix Θ under the irrepresentability condition, implying that all the underlying graph edges can be identified correctly through X^⋆. In addition, the condition ∥S−Σ∥max≤α4λ can hold with overwhelming probability 1 − c₁ exp(−c₂ logp) for Gaussian observations, where c₁ and c₂ are two positive constants.

A. Synthetic Data

We consider to learn Barabasi-Albert graphs that are useful in many applications. The graph structure and its weights associated with the edges are generated randomly. We obtain a weighted adjacency matrix W, where W_ij denotes the graph weight between node i and node j. We construct the underlying precision matrix by Θ = D− W, where Dis a diagonal matrix, which is generated to ensure that Θ is an M-matrix and the irrepresentability condition in Assumption 1 holds.

Given the underlying precision matrix Θ ∈ ℝ^p×p, we generate n independent observations y(1),…,y(n)∼N(0,Θ−1), and compute the sample covariance matrix by S=1n∑ni=1y(i)(y(i))⊤.

We compare the computational time with the state-of-the-art BCD [2], and GGL [6] in solving the ℓ₁-norm regularized Gaussian maximum likelihood estimation. It is observed in Table I that the three methods lead to the same objective function value, because all of them can obtain the global minimizer. However, our proposed algorithm takes significantly less computational time than BCD and GGL. The results are averaged over 10 Monte Carlo realizations. In addition, GGL can incorporate connectivity constraints in graph learning.

TABLE I: Comparisons of computational time in learning Barabasi-Albert graphs. The "Objective" denotes the objective function value, and the unit of time is seconds.

Next, we show that the estimator defined in (2) can recover the graph edges correctly under irrepresentability condition. We use F-score to measure the performance of edge recovery, which is defined by

View Source\begin{equation*}{\text{ F - score }}: = \frac{{2{\text{ TP }}}}{{2{\text{ TP }} + {\text{ FP }} + {\text{ FN }}}},\tag{14}\end{equation*}

It is observed in Figure 1 that our algorithm can achieve the F-score to be 1, implying that the graph structure is identified perfectly under the irrepresentability condition, which is consistent with our theoretical results in Theorem 1. We compare our method with the well-known Glasso method [25], which does not impose the MTP₂ constraints. We can observe that our algorithm can obtain a higher F-score than Glasso in the region of small sample size ratios. This is expected because imposing more prior knowledge always helps to improve the estimation performance especially when the samples are limited. The results are averaged over 30 realizations.

Fig. 1:

F-score as a function of the sample size ratio of n/p in learning Barabasi-Albert random graphs. The regularization parameter for both methods is set as λ = 0.05.

Show All

B. Financial Time-series Data

The MTP₂ models are justified well on the stock data since the market factor leads to the positive dependencies among stocks [3]. The data is collected from 189 stocks composing the S&P 500 index during a period from Oct. 1st 2017 to Jan. 1st 2018, resulting in 62 observations per stock, i.e., p = 189 and n = 62. We construct the log-returns data matrix by

View Source\begin{equation*}{X_{i,j}} = \log {P_{i,j}} - \log {P_{i - 1,j}},\tag{15}\end{equation*}

It is observed in Figure 2 that the performance of our algorithm is better than Glasso, because the latter has more gray-colored connections, which are between stocks from distinct sectors and often spurious from a practical perspective. The modularity values for Glasso and the proposed method are 0.41 and 0.45, respectively. A higher modularity value indicates a better representation of the actual network of stocks.

Fig. 2:

Stock graphs learned via (a) Glasso, and (b) the proposed method.

Show All