Scalable Inference of Sparsely-changing Gaussian Markov Random Fields
Year: 2021

Abstract:
We study the problem of inferring time-varying Gaussian Markov random fields, where the underlying graphical model is both sparse and changes sparsely over time. Most of the existing methods for the inference of time-varying Markov random fields (MRFs) rely on the regularized maximum likelihood estimation (MLE), that typically suffer from weak statistical guarantees and high computational time. Instead, we introduce a new class of constrained optimization problems for the inference of sparsely-changing Gaussian MRFs (GMRFs). The proposed optimization problem is formulated based on the exact `0 regularization, and can be solved in near-linear time and memory. Moreover, we show that the proposed estimator enjoys a provably small estimation error. We derive sharp statistical guarantees in the high-dimensional regime, showing that such problems can be learned with as few as one sample per time period. Our proposed method is extremely efficient in practice: it can accurately estimate sparsely-changing GMRFs with more than 500 million variables in less than one hour.

heading: 1 Introduction
text: Contemporary systems are comprised of massive numbers of interconnected components that interact according to a hierarchy of complex, unknown, and time-varying topologies. For example, with billions of neurons and hundreds of thousands of voxels, the human brain is considered as one of the most complex physiological networks [18, 22, 28, 30, 37]. The temporal behavior of today’s interconnected systems can be captured via time-varying Markov random fields (MRF). Time-varying MRFs are associated with a temporal sequence of undirected Markov graphs Gt(V,Et), where V and Et are the set of nodes and edges in the graph at time t. The node set V represents the random variables in the model, while the edge set Et captures the conditional dependency between these variables at time t. A popular approach for the inference of MRFs is based on the maximum-likelihood estimation (MLE): to obtain a model based on which the observed data is most probable to occur [42].
Despite being known as theoretically powerful tools [20, 39], MLE-based methods suffer from several fundamental drawbacks which render them impractical in realistic settings. First, they often suffer from notoriously high computational cost in massive problems, where the number of variables to be inferred is in the order of millions, or more. Second, they struggle to incorporate sparsity amongst
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
their components, which is pervasive in large-scale systems. In particular, while sparsely-changing MRFs can in theory be accurately estimated using sparsity-promoting regularizers (such as `0 penalty), most of the existing methods resort to relaxed or weaker variants of such regularization (such as `1 penalty), thereby suffering from inferior statistical guarantees.
To address the aforementioned challenges, we propose a class of constrained optimization problems that achieve superior statistical and computational guarantees, compared to the regularized MLE, for the inference of time-varying Gaussian MRFs (GMRFs). Our approach departs from the usual wisdom in statistics and machine learning that inference problems with nonconvex `0 terms are intractable, and convex proxies should be used instead. In particular, we show that the inference of sparsely-changing Gaussian MRFs can be solved efficiently via nonconvex `0 penalties.
Notations. The ith element of a time-series vector vt is denoted as vt;i; the (i, j)th element of a time-series matrix Vt is denoted as Vt;ij . For a vector v, the notation vi:j is used to denote the subvector of v from index i to j. For a vector v, the notations kvk1, kvk2, kvk0 denote the `1 norm, `2 norm, and the number of nonzero elements, respectively. Moreover, for a matrix M , the notations kMk2, kMk1, kMk1/1, kMk1/1 refer to the induced 2-norm, induced1-norm, `1/`1 norm, and `1/`1 norm, respectively. Moreover, we define kMko↵ = kMk1/1 Pd i=1 |Mii|. For a vector v and matrix M , the notations supp(v) and supp(M) are defined as the sets of their nonzero elements. Given two sequences f(n) and g(n), the notation f(n) . g(n) implies that there exists a constant C <1 that satisfies f(n)  Cg(n), and f(n) ⇣ g(n) implies that f(n) . g(n) and g(n) . f(n). All proofs are deferred to the supplementary file.

heading: 1.1 Warm-up: Regularized MLE for Sparsely-changing GMRFs
text: Consider a multivariate zero-mean Gaussian process {Xt}Tt=0 with distribution
P(Xt)=exp ⇢ 1
2 h⇥t, XtX
> t i+h⌘t, Xti A(µt,⇥t)
(1)
for t = 0, . . . , T where A : Rd⇥d ! R is the log-partition function used to normalize the distribution. Without loss of generality, we assume that the mean is zero. At any given time t, a sequence of data samples {X(i)t } Nt i=1 is collected from (1). The inference of time-varying GMRFs reduces to estimating the time-varying precision matrix ⇥t from the data samples {X (i) t } Nt i=1. Moreover, the edge set of the Markov graph Gt coincides with the off-diagonal nonzero elements of ⇥t [45].
We first illustrate the fundamental drawbacks of the `1-regularized MLE for time-varying GMRFs with sparsely-changing structures. The sparse precision matrices can be estimated via the following regularized MLE, also known as time-varying Graphical Lasso (GL) [11, 16]:
{b⇥t}Tt=0 = argmin ⇥t
TX
t=0
⇣ h⇥t, b⌃ti log det(⇥t) ⌘
+ 1
TX
t=0
k⇥tko↵ + 2
TX
t=1
k⇥t ⇥t 1k1/1 (2a)
s.t. ⇥t 0 t = 0, 1, . . . , T (2b)
where b⌃t 2 Rd⇥d is the sample covariance matrix at time t. Without loss of generality, we assume that the samples have zero mean. Example 1 below shows that (2) may lead to poor estimates.
Example 1. Consider a scenario where {⇥t}4t=0 2 R25⇥25 are randomly generated symmetric and sparse matrices. At each time t = 0, . . . , 4, the precision matrix ⇥t has exactly 30 off-diagonal elements with value one in its upper-triangular part, and the remaining off-diagonal entries are set to zero. Moreover, the diagonal entries ⇥t;ii are chosen as 1+ P j 6=i ⇥t;ij . At every time, 5 nonzero offdiagonal elements are changed to zero, and 5 zero elements are set to one. The sample covariance b⌃t is obtained by collecting 500 samples from the Gaussian distribution with the constructed precision
matrices. Figure 1a illustrates a heatmap of the mismatch error, i.e., the total number of mismatches
in the sparsity patterns of the true and estimated precision matrices and their differences, for different
values of the regularization coefficients. It can be seen that after an exhaustive search over the
regularization coefficient space, the best achievable mismatch error is in the order of 50. Thus, the
estimated parameters reveal little information about the true structure of the time-varying GMRF.
Moreover, Figure 1b depicts the concatena-
tion of the nonzero elements in the true pre-
cision matrices (dashed red line), and their
corresponding values in the estimated matri-
ces (blue curve). It can be seen that, even when
the sparsity pattern of the elements is correctly
recovered, the estimated nonzero entries are
“shrunk” toward zero, incurring a substantial
bias.
The above example shows the inferior statistical performance of the time-varying GL as an instance of a regularized MLE method for sparsely-changing GMRFs. In addition to its subpar statistical performance, time-varying
GL suffers from expensive computational complexity: a general-purpose interior-point algorithm for solving time-varying GL has a prohibitive per-iteration complexity of O(Td6) [10]. More recent algorithms for solving time-varying GL have lower per-iteration complexity of O(Td3) [16, 26, 32]. However, these methods suffer from a slow (sublinear) convergence rate of O(1/✏), increasing the overall complexity to O(Td3/✏) in order to obtain an ✏-accurate solution. Thus, there is an inherent tradeoff between the quality of the solution found and the computational time required, with the performance deteriorating sharply as the number of precision digits increase. Solvers with such computational complexity may fall short of practical use in the large-scale settings. We now discuss the proposed method, which finds optimal solutions to the relevant optimization problems in strongly polynomial time.

heading: 2 Proposed Approach
text: The proposed framework is based on exact solutions to a class of tractable discrete `0-problems, thus circumventing bias and other drawbacks of the standard `1-approximations, while guaranteeing the scalability of the proposed method. As a general framework, we study the optimization problem:
min (1 ) TX
t=0
k⇥tk0 + TX
t=1
k⇥t ⇥t 1k0 (3a)
s.t. k⇥t eF ⇤(b⌃t)k1/1  t t = 0, 1, . . . , T (3b) ⇥t 2 Rd⇥d t = 0, 1, . . . , T (3c)
where the optimal solutions {b⇥t}Tt=0 are the estimates of the precision matrices of the sparselychanging GMRF, b⌃t 2 Rd⇥d is the sample covariance matrix at time t, and eF ⇤(·) is an approximate backward mapping of the model. In particular, we use the approximate backward mapping proposed in [46], see §4.1 for a formal definition.
First, we establish a deterministic guarantee on the estimation error of the optimal solution to (3).
Theorem 1 (Estimation error and sparsistency). Let 0 < < 1. For every t = 0, . . . , T , define St as the set of indices corresponding to the nonzero elements of the true precision matrix ⇥⇤t , and define Dt as the set of indices corresponding to the nonzero elements of ⇥⇤t ⇥ ⇤ t 1. Assume that
• ⇥⇤t eF ⇤(b⌃t) 1/1 < t, 80  t  T , • 2 t  min(i,j)2St |⇥ ⇤ t;ij |, 80  t  T , • 2 t + 2 t 1  min(i,j)2Dt |⇥ ⇤ t;ij ⇥ ⇤ t 1;ij |, 80  t  T .
Then, the following statements hold for every 0  t  T : Sparsistency supp ⇣ b⇥t ⌘ = supp (⇥⇤t ) and supp ⇣ b⇥t b⇥t 1 ⌘ = supp ⇥⇤t ⇥ ⇤ t 1 . Estimation error kb⇥t ⇥⇤t k1/1  2 t and kb⇥t ⇥⇤t k2  2 p |St| t.
Theorem 1 presents a set of conditions under which the proposed estimation method achieves sparsistency and small estimation error. The first condition entails that the true precision matrix ⇥⇤t is a feasible solution to (3). The second and third conditions imply that there is a non-negligible gap between the zero and nonzero elements of the true parameters and their temporal changes. In §4 we present additional bounds specific to our choice of backward mapping.
Theorem 2 (Computational complexity). Given eF ⇤(b⌃t), the optimization problem (3) can be solved to optimality in at most O((dT )2) time and memory on a single thread.
Theorem 2 shows that the optimization problem (3) can be solved efficiently and in strongly polynomial time despite its non-convex nature. As will be explained later, our choice of approximate backwards mapping [46] requires inverting T + 1 matrices, each with size d⇥ d, thereby increasing the overall complexity to O((dT )2 + Td3). If T  d, the complexity of the our method is dominated by that of the matrix inversion, which is unavoidable, even if the sample covariance matrix coincides with its true analog (since we still need to invert them to obtain the true precision matrices).
Our solution method for (3) relies on the element-wise decomposability of (3): we decompose (3) into smaller subproblems over different coordinates of {⇥t}Tt=0. Then, we show that the optimal solution to each subproblem can be obtained by solving a shortest path problem on an auxiliary weighted directed acyclic graph (DAG). The details of the solution method are presented in §5. Moreover, the proposed algorithm is easily parallelizable, leading to better runtimes in practice.
Example 1 (continued). Figure 2 depicts the
performance of the proposed method, compared
to the regularized MLE with 1 = 0.14 and 2 = 0.16 (corresponding to the smallest mismatch error) for the instances generated in Ex-
ample 1. The regularization parameter in the objective function of (3) is set to 0.2. Moreover, for simplicity, we set 0 = · · · = 4 = . Figure 2a demonstrates that the proposed method
enjoys a significantly smaller mismatch error,
for a wide range of . On the other hand, Fig-
ure 2b shows that the synthetic bias caused by
the `1 penalty in the regularized MLE is allevi-
ated via the proposed method.

heading: 3 Related Works
text: Time-varying MRF. In addition to the time-varying GL introduced in §1.1, a recent line of works have studied the inference of smoothly-changing GMRFs [23, 15, 48], where a kernel averaging technique combined with Graphical Lasso is used to estimate the smoothly-changing precision matrices. However, these methods do not leverage the prior information about the sparsity of the parameter differences. With the goal of addressing this deficiency, several works have studied the inference of sparsely-changing MRF (also known as sparse differential networks) [43, 47, 27]. However, the main drawback of these methods is that they only estimate the parameter differences, and their theoretical guarantees are restricted to problems with two time steps (T = 0, 1). Similarly, regression-based approaches have been proposed for change point detection problems [24, 36] with MRFs and two time periods, assuming the sparsity pattern of all entries of the precision matrices change at the same time. In contrast, [44] studies the inference of sparse MRFs given an index variable under the assumption that the sparsity pattern is invariant, whereas [14] assumes that the precision matrix is a linear function of the index variable.
Sparsity-promoting optimization. Optimization problems with `0 terms are often deemed to be intractable, and approximations are solved instead. Perhaps the most popular approach is the fused lasso [34, 40, 38, 41], which calls for replacing terms k⇥tk0 and k⇥t ⇥t 1k0 with their `1- approximations. Nonetheless, such approximations result in subpar statistical performance when compared with exact `0 methods [19, 29].
Exact or near-optimal methods for optimization problems of the form
min ✓2[`,u]p
k✓k0 + pX
i=1
pX
j=i+1
gij(✓i ✓j), (4)
for given one-dimensional functions gij : R! R, have also been studied in the literature. If functions gij are convex, then problem (4) admits pseudo-polynomial time algorithm [3, 7]. Moreover, convex relaxations that deliver near-optimal solutions for (4) were proposed for the special case of convex quadratic g functions [6]; if, additionally, we have ` = 0 and u = 1, then problem (4) is in fact solvable in strongly polynomial time [5]. On the other hand, problem (4) is much more challenging for non-convex g: if g(x) = {x 6= 0}, as is the case in (3), then problem (4) is NP-hard even if the term k✓k0 is dropped from the objective [17]. Nonetheless, as we show in this paper, problem (4) can be solved efficiently in the context of time-varying MRFs, where gij(x) = 0 whenever j > i+ 1.

heading: 4 Statistical analysis
text: Theorem 1 presents a set of deterministic conditions under which the estimates from (3) enjoy zero mismatch error and small estimation error. However, the formulation of (3) is contingent upon the availability of an accurate backward mapping, and a choice of t that satisfies the conditions of Theorem 1. In this section, we show how to efficiently design sample-efficient approximate backward mappings, and select t accordingly for the class of sparsely-changing GMRFs. Moreover, we use the deterministic conditions of Theorem 1 to arrive at a non-asymptotic probabilistic guarantee for the inference of time-varying GMRFs under different prior knowledge on their temporal behavior, such as sparsity and smoothness.

heading: 4.1 Sparsely-changing GMRFs
text: Given the true covariance matrix ⌃t, the backward mapping of time-varying GMRFs as defined in (1) takes the form F ⇤(⌃t) = ⌃ 1t . In light of this closed-form expression for the backward mapping, a commonly-used approximation is eF ⇤(b⌃t) = b⌃ 1t , where b⌃t = 1Nt PNt i=1 X (i) t X (i) t > is the sample covariance matrix. However, in the high-dimensional settings where d Nt, this approximate backward mapping is not well-defined, since the sample covariance matrix is highly rank-deficient.
To address this issue, [46] propose a proxy backward mapping for high-dimensional settings: consider the soft-thresholding operator ST⌫(M) : Rd⇥d ! Rd⇥d, where ST⌫(M)ij = Mij sign(Mij)min{|Mij |, ⌫} if i 6= j, and ST⌫(M)ij = Mij if i = j. The approximate backward mapping is then given by eF ⇤(b⌃t) = [ST⌫(b⌃t)] 1, which is well-defined, even in the high-dimensional setting, with an appropriate choice of the threshold ⌫ [46]. We make the following assumption. Assumption 1 (Bounded norm). There exist constant numbers 1 <1, 2 > 0, and 3 <1 such that, for every t = 0, . . . , T ,
k⇥tk1  1, inf w:kwk1=1 k⌃twk1 2, k⌃tk1/1  3.
Assumption 1 implies that the true covariance matrices and their inverses have bounded norms. Another key notion which plays an important role in our statistical guarantees is the weak sparsity of the covariance matrices, as defined next.
Definition 1 (Weak sparsity). Given 0  q  1 and dimension d, define sd(q, d) def = maxi Pd j=1 |[⌃t]ij | q .
Assuming that sd(0, d)⌧ d, Assumption 1 reduces to the covariance matrix being sparse. Moreover, in many cases, a sparse inverse covariance matrix leads to weakly sparse covariance matrices. For instance, if ⇥t has a banded structure with small bandwidth, then it is known that the elements of ⌃t = ⇥ 1 t enjoy exponential decay away from the main diagonal elements [12, 21]. Under such circumstances, one can verify that sd(q, d)  C1 ⇢q for some constant C > 0 and ⇢ < 1. More generally, a similar statement holds for a class of inverse covariance matrices whose support graphs have large average path length [8, 9]; a large class of inverse covariance matrices with row- and column-sparse structures satisfy this condition. Theorem 3 states that the proposed method results in
high-quality solutions provided that the number of samples Nt is sufficiently large with respect to the weak sparsity of the covariance matrices. Theorem 3. Consider a sparsely-changing GMRF and let ⇣ = max{logd(T + 1), 1}. Given an arbitrary ⌧ > ⇣ + 2, let Nt & sd(q, d) 2 1 q ⌧ log d, for some 0  q < 1. Then, with the approximate
backward mapping eF ⇤(b⌃t) = [ST⌫t(b⌃t)] 1, and parameters ⌫t ⇣ q ⌧ log d Nt and t ⇣ q ⌧ log d Nt , the estimates {b⇥t}Tt=0 obtained from (3) satisfy the following statements for all t = 0, 1, . . . , T , with probability of at least 1 4d ⌧+⇣+2: Sparsistency We have supp ⇣ b⇥t ⌘ = supp (⇥⇤t ) and supp ⇣ b⇥t b⇥t 1 ⌘ = supp ⇥⇤t ⇥ ⇤ t 1 . Estimation error kb⇥t ⇥⇤t k1/1. q ⌧ log d Nt and kb⇥t ⇥⇤t kF . q ⌧ |St| log d Nt .
Theorem 3 is a direct consequence of Theorem 1 and provides, to the best of our knowledge, the first non-asymptotic guarantee on the inference of sparsely-changing GMRFs with an arbitrary length of time horizon T . In particular, it shows that the proposed optimization (3) guarantees small estimation error and zero mismatch error for sparsely-changing GMRFs, provided that Nt scales logarithmically with the dimension of the precision matrices. In the static setting (T = 0), the derived bound recovers the existing results on the sample complexity of learning static GMRFs [25, 33, 35].

heading: 4.2 Sparsely-and-smoothly-changing GMRFs
text: In many applications, such as financial markets and motion detection in video frames, the associated graphical model should be learned “on-the-go”, as the data arrives with a continuously changing graphical model. Under such circumstances, one may have access to few (or even one) samples at each time.
Suppose that the precision matrices change smoothly over time. Such smooth changes can be modeled via a continuous function ⇥(x) : [0, 1] ! Rd⇥d with uniformly bounded element-wise second derivatives [⇥(x)ij ]00 = d2⇥(t)ij dx2 , such that ⇥ ⇤ t = ⇥(t/T ) [15, 48]. If ⇥(t) ⌫ aI for every t 2 [0, 1] and some a > 0, then the covariance matrix ⌃(t) = ⇥(t) 1 is well-defined and smooth. Then, the problem of inferring the time-varying GMRF reduces to estimating a sequence of precision matrices {⇥(0),⇥(1/T ), . . . ,⇥(1)} given the samples Xt ⇠ N (0,⌃(t/T )). To alleviate the scarcity of samples, [15, 48] leverage the smoothness of the precision matrices, by taking the weighted average of the samples over time, where the weights are obtained from a nonparametric kernel. In particular, consider the weighted sample covariance matrix b⌃wt :
b⌃wt = tX
s=0
w(s, t)⌃s, where w(s, t) = 1
Th K
✓ s t
Th
◆ (5)
and K(·) is a symmetric nonnegative kernel that satisfies a set of mild conditions which hold for most standard kernels, including the (truncated) Gaussian kernel. These conditions are delineated in the appendix. Next, we present the counterparts of Assumption 1 and Definition 1 for sparsely-andsmoothly-changing GMRFs. Assumption 2 (Bounded norm). There exist constant numbers 1 <1, 2 > 0, and 3 <1 such that, for every t 2 [0, T ],
k⇥(t/T )k1  1, inf w:kwk1=1 k⌃(t/T )wk1 2, k⌃(t/T )k1/1  3.
Definition 2 (Weak sparsity). Given 0  q  1 and dimension d, define sc(q, d) def = maxi,t2[0,T ] Pd j=1 |[⌃(t/T )]ij | q .
We now present the analog of Theorem 3 for sparsely-and-smoothly-changing GMRFs. Theorem 4. Consider a sparsely-and-smoothly-changing GMRF with one sample per time, let
⇣ = max{logd(T + 1), 1}, and suppose that the sample covariance matrices are constructed according to (5) with h ⇣ T 1/3. Given an arbitrary ⌧ > ⇣ + 2, let T & sc(q, d) 3 1 q (⌧ log d)3/2. Then, with the approximate backward mapping eF ⇤(b⌃t) = ST⌫t(b⌃wt ) and parameters ⌫t ⇣ p ⌧ log d T 1/3
and t ⇣ p ⌧ log d T 1/3
, the estimates {b⇥t}Tt=0 obtained from (3) satisfy the following statements for all T = 0, 1, . . . , T with probability of at least 1 d ⌧+⇣+2:
Sparsistency supp ⇣ b⇥t ⌘ = supp (⇥⇤t ) and supp ⇣ b⇥t b⇥t 1 ⌘ = supp ⇥⇤t ⇥ ⇤ t 1 .
Estimation error kb⇥t ⇥⇤t k1/1. q ⌧ log d T 2/3 and kb⇥t ⇥⇤t kF . q ⌧ |St| log d T 2/3 .
Theorem 4 shows how the smoothness assumption on the true covariance matrix can be used to construct the backward mappings using the samples collected during the entire time horizon, thereby significantly reducing the sample complexity of learning time-varying GMRFs. In particular, leveraging the smoothness of the covariance matrix can reduce the minimum required number of samples from O(T log d) to O((log d)1.5). On the other hand, Theorem 3 does not impose any lower bound on T , and its estimation error decays faster in terms of the sample size.

heading: 5 Solution Method
text: In this section, we describe the proposed algorithm for solving (3). For the simplicity of notation, we define the lower bound and upper bound matrices lt and ut as lt;ij = [ eF ⇤(b⌃t)]ij t and ut;ij = [ eF ⇤(bµt)]ij + t, for every 1  i, j  d. The following fact plays a key role in our analysis.
Fact 1. An optimal solution n b⇥t oT t=0 of (3) satisfies for every 1  i < j  d,
n b⇥t;ij oT t=0 2 argmin {⇥t;ij}Tt=0 (1 ) TX t=0 {⇥t;ij 6= 0}+ TX t=1 {⇥t;ij ⇥t 1;ij 6= 0} (6a)
s.t. lt;ij  ⇥t;ij  ut;ij 80  t  T. (6b)
Fact 1 implies that (3) decomposes into the smaller subproblems (6). Therefore, our main focus is devoted to solving each subproblem independently. To further simplify the notation, we drop the subscript ij from (6), whenever it is chosen arbitrarily, and use ✓t instead of ⇥t;ij .
Let OPTi!j( ) denote the truncated problem from time i to j with the regularization coefficient :
f ⇤ i!j = min
{✓t}jt=i (1 )
jX
t=i
{✓t 6= 0}+ jX
t=i+1
{✓t ✓t 1 6= 0} (7a)
subject to lt  ✓t  ut 8i  t  j. (7b)
Let the objective function for a candidate solution ✓ be denoted as fi!j(✓); by convention, we let fi!j(✓) = 0 whenever j < i. Moreover, the optimal objective value and the set of optimal solutions to OPTi!j( ) are respectively denoted as f⇤i!j and X ⇤i!j . Similarly, b✓i!j 2 X ⇤i!j is used to denote an optimal solution to OPTi!j( ). We omit the subscript i ! j whenever i = 0 and j = T . The tth feasible interval is defined as t = [lt, ut]. Accordingly, the notation \t!s refers to \t!s def = t \ t+1 \ · · · \ s.

heading: 5.1 Special case: = 1
text: As the first step, we consider the special case = 1, and provide an efficient algorithm (Algorithm 1) for solving OPT0!T (1), where the sparsity is only promoted on the parameter differences (and not on the individual parameters). As will be shown later, Algorithm 1 will be used as a subroutine in our proposed algorithm for the general case 0 < < 1. At a high level, Algorithm 1 recursively performs the following operations: at any given time ⌧ , the algorithm looks into the future to find a nonempty interval that is feasible for the longest possible time . Then, it sets the subvector ✓⌧ : to an arbitrarily chosen element from this nonempty interval.
Proposition 1. Greedy(l, u, 0, T ) returns an optimal solution {✓Greedyt } T t=0 to OPT0!T (1). Moreover, the truncated solution {✓ Greedy t } j t=0 is optimal for OPT0!j(1).

heading: 5.2 General case: 0 < < 1
text: Now, we present our main algorithm for the general case 0 < < 1.
Algorithm 1 Greedy(l, u, ⌧, T )
1: Output: Solution ✓Greedy⌧!T , the objective value f Greedy ⌧!T to OPT⌧!T (1), and the index set of
maximal nonempty intervals
2: Find largest such that \⌧! 6= 0; 3: Set ✓Greedy⌧ : = ⌘ for some ⌘ 2 \ ⌧! ; 4: if  T 1 then 5: Execute Greedy(l, u, + 1, T ); 6: end if 7: Set fGreedy⌧!T = PT t=⌧+1 {✓ Greedy t ✓ Greedy t 1 6= 0}; 8: Set [ { }; 9: Return {✓Greedyt } T t=⌧ , f Greedy ⌧!T , and ;
Definition 3. The set Zi!j = {i, i + 1, . . . , j} is called a zero-feasible sequence if lk  0  uk for every k 2 Zi!j . Moreover, the zero-feasible sequence Zi!j is called maximal if it is not strictly contained within another zero-feasible sequence.
Let Zi1!j1 ,Zi2!j2 , . . . ,ZiZ!jZ be the set of all maximal zero-feasible sequences such that 0  i1  j1 < i2  j2 < · · · < iZ  jZ  T , where Z is the number of maximal zero-feasible sequences. If Z = 0, i.e., there is no zero-feasible sequence, then it is easy to see that PT t=0 {✓t 6= 0} = T + 1 for every feasible solution, and hence, Greedy(l, u, 0, T ) leads to an optimal solution to (7). Another special case is when i1 = 0 and j1 = T , in which case the optimal solution is b✓ = 0. Therefore, without loss of generality, suppose that Z > 0 and either i1 6= 0 or j1 6= T .
Our goal now is to obtain an optimal solution to (7) by solving a shortest path problem over a weighted directed acyclic graph (DAG) whose nodes correspond to the maximal zero-feasible sequences. In particular, consider a weighted DAG G with the vertex set V = {0, 1, . . . , Z, Z + 1}, where the vertices k and l are connected via a directed arc (k, l) if k < l. Moreover, for every arc (k, l), the weight W (k, l) = 0 if (k, l) = (0, 1), i1 = 0 or (k, l) = (Z,Z + 1), jZ = T , and
W (k, l) =(1 )(il jk 1) + f Greedy jk+1!il 1 + {k 6= 0}+ {l 6= Z + 1} (8)
otherwise, where we define j0 = 1 and iZ+1 = T + 1.
Algorithm 2 Algorithm for solving (7)
1: Output: Optimal solution b✓ and the objective value f⇤ to OPT0!T ( ) 2: Find the maximal zero-feasible sequences; 3: Construct the DAG G with weights defined as (8); 4: Find the shortest path p = (v1, v2, . . . , vr) between the vertices 0 and Z + 1 in G; 5: Set b✓jvi+1:ivi+1 1 = ✓ Greedy jvi+1!ivi+1 1
and b✓ivl :jvl+1 = 0 for every l = 1, 2, . . . , r; 6: Return {b✓}Tt=0 and f⇤
Theorem 5. The shortest path from 0 to Z + 1 on G has value f⇤.
The above theorem implies that the optimal solution to (3) can be obtained via Algorithm 2. Theorem 6. Problem (7) can be solved in O(ZT ) time and memory.
Since Z = O(T ) and a solution to (3) requires solving O(d2) instances of (7), we find the total complexity stated in Theorem 2. Note however that if Z = O(1), then the overall complexity reduces to O(Td2), which is linear in the total number of variables. In the next section, we will show that the practical runtime of the proposed algorithm is near-linear with respect to the number of variables.

heading: 6 Numerical Analysis
text: In this section, we evaluate the performance of the proposed estimator in synthetically generated massive datasets, and a case study on the correlation network inference in stock markets. We refer
the reader to the supplementary file for an extensive discussion of our simulations (e.g., the choice of parameters, robustness analysis, and a comparison with other state-of-the-art methods). In all of our simulations, the parameters ⌫t and t are chosen directly from the data samples, i.e., without prior knowledge of the true solution, via Bayesian Inference Criterion (BIC) [31, 13].
Case Study on Large Datasets We consider randomly generated instances of sparsely-changing GMRFs, where the true inverse covariance matrix is constructed as follows: at time t = 0, we set ⇥0 = Id + P (i,j)2S A (i,j), where A(i,j) is a sparse positive semidefinite matrix with exactly two nonzero off-diagonal elements. For every (i, j) 2 S, we set A(i,j)ij = A (i,j) ji = 0.4 and A (i,j) ii = A (i,j) jj = 0.4. Clearly, A (i,j) ⌫ 0, and hence, ⇥0 0. In the first experiment, we fix T = 10 and change the values of d. The number of nonzero elements in the individual precision matrices and their differences are set to 3d and 0.04d, respectively. We evaluate the performance of the proposed method in the high dimensional settings, where Nt = d/2 for every t = 0, . . . , T . Moreover, define TPR and FPR for the individual parameters and their differences as the true positive and false positive values, normalized by the total number of nonzero and zero elements in the true precision matrices and their differences, respectively. Clearly, both TPR and FPR are between 0 and 1, with TPR = 1 and FPR = 0 corresponding to the perfect recovery of the sparsity patterns. Figure 7 depicts TPR, FPR, and the max-norm error of the estimated parameters, as well as the runtime of our algorithm for different values of d with and without parallelization. It can be seen that both TPR and FPR values improve with the dimension for the estimated parameters and their differences. Moreover, with a single core, the runtime of our algorithm scales almost linearly with d2, which is in line with the result of Theorem 2. Using 5 cores, the runtime of our algorithm is improved by 40% on average. Using 10 cores deteriorates the performance due to the shared memory limitations. Using our algorithm, we reliably infer instances of sparsely-changing GMRFs with more almost 500
million variables in less than one hour.
Case Study on Stock Market Finally, we illustrate the performance of our algorithm for the inference of stock correlation network. In particular, we consider an investor seeking to understand relationships between securities over time. Sparsity of the precision matrices guarantees interpretability,
while the sparsely-changing structure is imposed to identify sharp changes in market conditions, and a need to rebalance the portfolio. We consider the daily changes for 214 securities from 1990/01/04 to 2017/08/10, with the total number of 6990 days (d = 214 and T = 6990). Due to the continuously changing nature of the stock correlation network, we will use the kernel averaging approach introduced in Subsection 4.1. Using the constructed sample covariance matrices, we estimate the sparsely-changing precision matrix ⇥(t/T ) at discrete times t 2 {30, 60, 90, . . . , 6990}.
A drastic change in the correlation network signals a spike in the stock market, which may reflect the market’s response to unexpected events. Figure 11 shows the number of changes in the estimated network, for the choices of ⌫0 = 3, 0 = 0.16, and = 0.9, together with the historical chart of National Association of Securities Dealers Automated Quotations (NASDAQ) [1]. It can be seen that the major spikes in the estimated network can be attributed to the historical stock market crashes. For instance, the spikes A, B, and C respectively correspond to the “early 1990s recession”, “dot-com bubble”, and “global financial crisis”; see [4] for more details. Interestingly, the estimated network can also detect other historical (but less severe) downturns in 2011 (point D) and 2016 (point E). In the supplementary materials, we provide a more detailed analysis with different choices of parameters.
Acknowledgments and Disclosure of Funding
S.F. is supported by MICDE Catalyst Grant and MIDAS PODS Grant. A.G. is supported, in part, by the National Science Foundation under grant CIF 2006762.

Review summaries:
In this paper, the authors consider the problem of learning Sparsely changing Gaussian Markov Random Fields. They pose the problem as a non-convex optimization problem that can be solved exactly. The authors provide statistical guarantees. Furthermore, they also provide a computationally efficient algorithm to recover the graphical model.

  This paper study the problem of inferring time-varying Gaussian Markov random fields, where the underlying graphical model is both sparse and changes {sparsely} over time. A new class of constrained optimization problems was introduced for the inference of sparsely-changing Gaussian MRFs (GMRFs). The proposed optimization problem is formulated based on the exact regularization and can be solved in near-linear time and memory.   

This submission considers the problem of estimating the time-varying precision matrices corresponding to a Gaussian graphical model. The authors propose an \ell_0-regularized optimization problem (that runs in polynomial time), and establish theoretical guarantees as well as its performance on simulated and real datasets. The paper is quite well-written and seems like an important contribution to this active research area.