Some results on the Gaussian Markov Random Field construction problem based on the use of invariant subgraphs

The study of Gaussian Markov Random Fields has attracted the attention of a large number of scientific areas due to its increasing usage in several fields of application. Here, we consider the construction of Gaussian Markov Random Fields from a graph and a positive-definite matrix, which is closely related to the problem of finding the Maximum Likelihood Estimator of the covariance matrix of the underlying distribution. In particular, it is simultaneously required that the variances and the covariances between variables associated with adjacent nodes in the graph are fixed by the positive-definite matrix and that pairs of variables associated with non-adjacent nodes in the graph are conditionally independent given all other variables. The solution to this construction problem exists and is unique up to the choice of a vector of means. In this paper, some results focusing on a certain type of subgraphs (invariant subgraphs) and a representation of the Gaussian Markov Random Field as a Multivariate Gaussian Markov Random Field are presented. These results ease the computation of the solution to the aforementioned construction problem.


Introduction
A MRF is a random vector whose components are associated with the nodes of an undirected graph and the conditional distributions satisfy some properties (called Markov properties) over the structure of the graph. In this model, the dependence of two non-adjacent variables in the graph is explained in terms of the dependence between adjacent variables of the chains that connect these two non-adjacent variables (Rue and Held 2005). MRFs are a simple way to model spatial dependencies between the components of the random vector. For example, for modeling the evolution of a disease over some regions of a country, it is quite typical to consider an MRF over a graph in which two nodes are adjacent if the regions have a common border (see Section 4.4.2 in Rue and Held 2005).
The case in which the random vector has multivariate Gaussian distribution has attracted a lot of attention since; in such case, the MRF can be characterized by the null elements of the inverse of the covariance matrix. Furthermore, the modeling of a problem by means of a GMRF is, in general, computationally attractive, mainly because the inverse of the covariance matrix is positive-definite and typically sparse (Rue and Held 2005). The multivariate generalizations of GMRFs, called Multivariate Gaussian Markov Random Fields (MGMRFs), have also been widely studied in recent years (MacNab 2018).
In this paper, we follow the approach by Speed and Kiiveri (1986) in which they propose a GMRF construction problem over a graph. Formally, given a multivariate Gaussian distribution and an undirected graph, we search for another multivariate Gaussian distribution that keeps the same variances and the covariances between adjacent variables while satisfying the Markov properties. This problem is analogous to finding the Maximum Likelihood Estimator (MLE) of the covariance matrix of the GMRF given the sample covariance matrix of the distribution. If the GMRF model is reasonable, then this MLE will be similar to the sample covariance matrix but will also benefit from the computational advantages of having a sparse inverse. In this direction, several methods that allow to learn the underlying graph structure have been developed (see Furtlehner et al. 2021;Ferrer-Cid et al. 2021;Loftus et al. 2021 for some recent examples). Here, we solve this GMRF construction problem restricted to a certain type of subgraphs (which we will call invariant subgraphs) and then solve the whole GMRF construction problem by considering MGMRFs over forests. This will result in a significant reduction in computation time in comparison to directly solving the GMRF construction problem over the whole graph.
The remainder of the paper is organized as follows. In Sect. 2, we introduce the preliminary concepts necessary for presenting our results. The notions of invariant subgraph and complete separator are presented in Sect. 3. In Sect. 4, we define an MGMRF construction problem over a forest as a reduction in several GMRF construction problems over invariant subgraphs and we provide an easy formula for solving such construction problem. In Sect. 5, we study the complexity of the presented method and compare it with the original algorithms. Finally, in Sect. 6, we illustrate the presented results with three different examples and show the benefits of the considered approach. Some conclusions are presented in Sect. 7.

Simple undirected finite graphs
A simple undirected finite graph G = (V , E) consists of a couple of sets: a finite set of vertices (or nodes) V and a finite set of edges E, whose elements (i, j) are such that i, j ∈ V and satisfy that (i, j) ∈ E if and only if ( j, i) ∈ E and (i, i) / ∈ E for any i ∈ V . The cardinal of V is called the order of the graph and is usually denoted by n and the cardinal of E is called the size of the graph and is denoted by m. The subgraph (of If (i, j) ∈ E, then i and j are said to be adjacent. The number of adjacent nodes to a node is called the degree of incidence of such node. A sequence of nodes (a 1 , . . . , a k ) is called a chain from a 1 to a k if (a i , a i+1 ) ∈ E for all i ∈ {1, . . . , k − 1}. The number of edges involved in a chain is called the length of the chain. If a 1 = a k , then the chain is also called a cycle. Two chains that only share the first and the last node are called internally disjoint.
If there exists a chain between two nodes, both nodes are said to be connected. If there exist two internally disjoint chains between them, the nodes are called 2connected. A graph is called connected if all pairs of nodes are connected and is called 2-connected if all pairs of nodes are 2-connected. A subgraph that is connected (respectively, 2-connected) and maximal with respect to this property is called a connected (respectively, 2-connected) component.
A graph is called complete if all pairs of nodes are adjacent. A graph that does not contain any cycle is called a forest and, if it is also connected, then it is called a tree.
A subset V ⊂ V is called a separator of G = (V , E) if the subgraph induced by V \V is not connected. Given three pairwisely disjoint sets V , V 1 , V 2 ⊂ V , such that V 1 and V 2 are non-empty, we say that V separates V 1 and V 2 if any chain between a node in V 1 and a node in V 2 contains a node in V .
Given two graphs G 1 = (V 1 , E 1 ) and G 2 = (V 2 , E 2 ), the Cartesian product of G 1 and G 2 is the graph are adjacent if and only if it either holds that u 1 = u 2 and (v 1 , v 2 ) ∈ E 2 or v 1 = v 2 and (u 1 , u 2 ) ∈ E 1 .

Multivariate Gaussian distribution
A continuous random vector X has a multivariate Gaussian distribution if any linear combination of its components has a univariate Gaussian distribution (Mardia et al. 1979). Its joint density distribution has the following expression: where µ is the mean vector and Σ is the covariance matrix. We denote the set of positive-definite matrices by P and we require Σ ∈ P. Given two subsets of indices A and B, we denote the corresponding submatrix of Σ by Σ AB or simply by Σ A in case A = B. We also denote X\X A by X −A .
Given three continuous random vectors X, Y and Z of dimensions n X , n Y and n Z , respectively, with joint density function f (x, y, z), X and Y are said to be conditionally independent given Z (see, e.g., Rohatgi 1976) iff there exist h : R n X +n Z → [0, ∞] and g : R n Y +n Z → [0, ∞] such that f (x, y, z) = h(x, z)g(y, z) , ∀x ∈ R n X , ∀y ∈ R n Y , ∀z ∈ R n Z . We denote the fact that X A and X B are conditionally independent given X C by X A ⊥ X B | X C .
In particular, the conditional independence structure of a Multivariate Gaussian distribution can be characterized by using Σ −1 , often referred to as the precision matrix.
Theorem 1 (Rue and Held 2005) Let X be a multivariate Gaussian random vector with mean vector µ and covariance matrix Σ. For any i = j, it holds that

The GMRF model
Let G = (V , E) be a simple undirected finite graph where V = {1, ..., n} denotes the set of nodes and E ⊂ V × V denotes the set of edges. The neighborhood of a node i (i.e., the nodes that are adjacent to i) is denoted by N (i). Given a random vector X = (X 1 , ..., X n ), the Markov properties are defined as follows: -The pairwise Markov property: -The global Markov property: X A ⊥ X B | X C , for any pairwisely disjoint A, B, C ⊂ V with A, B = ∅ and where C separates A and B.
If X is a multivariate Gaussian random vector, then the three properties above are equivalent (Rue and Held 2005) and if any of the properties above is satisfied, X is called a Gaussian Markov Random Field (GMRF) over G. As a result of Theorem 1, given a GMRF, the Markov properties are characterized by the null elements of Σ −1 .
A multivariate version of the GMRFs can be also defined. Let X = (X 1 , . . . , X n ) be a multivariate Gaussian random vector and G = (V , E) with V = {1, ..., n}. Similarly to a GMRF, any node can be associated with a random vector instead of a random variable. This results in the notion of Multivariate Gaussian Markov Random Field (MGMRF). In this case, the Markov properties are characterized by null submatrices of Σ −1 .

The GMRF construction problem
We focus on the construction of a GMRF over a graph when the covariances between adjacent variables are fixed. We start with a positive-definite matrix P of dimension n × n and we allow to change the non-fixed values in order to satisfy the Markov properties over the graph. Given the matrix P as initial data, the search for a matrix F that coincides with P for adjacent variables and whose inverse F −1 has zeros at the positions associated with non-adjacent variables is referred to as the GMRF construction problem.
Theorem 2 (Speed and Kiiveri 1986) Let P,R ∈ P and G = (V , E) be a graph. There exists a unique F ∈ P that satisfies: In particular, setting the matrix R in the theorem above as the identity matrix, the existence of a unique solution to the GMRF construction problem is assured. The generalization of this construction problem to MGMRFs is immediate, with also a unique solution. Some algorithms to compute an approximate solution to this problem are provided in Speed and Kiiveri (1986) and Wermuth and Scheidt (1977). We highlight the importance of finding the solution to this problem, since it is related to several widely studied problems. Firstly, the problem above is equivalent to finding the MLE (Maximum Likelihood Estimation) of the covariance matrix of the GMRF given the sample covariance matrix (Dempster 1972;Xu et al. 2011;She and Tang 2019). More precisely, suppose that we have a sample of a GMRF and compute the sample covariance matrix. If we identify this sample covariance matrix with the matrix P in Theorem 2 and set R as the identity matrix, the matrix F corresponds to the MLE of the covariance matrix of the GMRF.
Proposition 1 Let X be a GMRF over G and consider a random sample of X. If the matrix P in Theorem 2 is the resulting sample covariance matrix, then the matrix F in Theorem 2 is the Maximum Likelihood Estimation of the covariance matrix of X.
This result is relevant to applications such as those in Furtlehner et al. (2021), Ferrer-Cid et al. (2021 and Loftus et al. (2021). If the assumption of the data being sampled from a GMRF is reasonable, then the MLE is similar to the sample covariance matrix while also having a sparse inverse matrix. In order to decide whether a GMRF model is reasonable or not, we can use a likelihood ratio test or apply some penalty procedures (Banerjee et al. 2008).
Secondly, the problem above is also equivalent to finding the distribution that maximizes the differential entropy among all the random vectors with the variances and some of the covariances specified (since it maximizes the determinant of the associated covariance matrix Grone et al. 1984).

The matrix F in Theorem 2 maximizes the determinant among all matrices in M.
This problem was also introduced by Dempster (1972) as a covariance selection model, which has been shown to be very useful for reducing the number of parameters in the estimation of the covariance matrix of a Multivariate Gaussian distribution (and, actually, of the exponential family) (Speed and Kiiveri 1986).
Note that there is a direct link between this problem and the construction of a GMRF from the marginal distribution of the cliques (complete subgraphs) of the graph by means of the positive definite completion of partial Hermitian matrices (Grone et al. 1984). In particular, the following result is a direct consequence of Theorem 2.

Proposition 3 Let G = (V , E) be a graph and C be the set of cliques of G. Consider
{X C | C ∈ C} to be a set of marginal multivariate normal distributions over the cliques of G satisfying that, for any C 1 , C 2 ∈ C, X C 1 restricted to C 2 and X C 2 restricted to C 1 are equally distributed. Let us denote by A the partial matrix such that A i, j = Cov(X i , X j ) if there exists C ∈ C with i, j ∈ C and all other elements unknown. It holds that a GMRF Y over G such that Y C = X C for all C ∈ C exists if and only if there exists a positive definite completion of A.
Another well-known approach to the construction of a GMRF starts from the full conditionals (see Section 2.2.4 in Rue and Held 2005), which is equivalent to determining the mean vector and the inverse of the covariance matrix (also known as the precision matrix). The main difference with our approach is that it is necessary to make additional considerations to assure the positive definiteness of the precision matrix. In our case, given a positive definite matrix, which in real-life applications typically is the sample covariance matrix, the solution to the GMRF construction problem always exists and is unique. However, if we are working with a degenerate GMRF, in which at least one variable may be expressed as a linear combination of the other ones, working with the precision matrix is a better option (see Sect. 3 in Rue and Held 2005).

Construction of invariant subgraphs from complete separators
Subvectors of a GMRF over a graph are not, in general, GMRFs over the subgraph induced by some of their components. For example, consider a GMRF over a tree with 3 nodes. The subgraph induced by the two nodes with degree of incidence 1 consists of two non-adjacent nodes. If the associated subvector were a GMRF over this subgraph, the associated variables would then be independent, a statement that is not true in general. For this very reason, we search for a type of subgraphs for which this property holds.
The term invariant subgraph refers to the fact that, given the submatrix of the initial data matrix P associated with the invariant subgraph, the submatrix of the solution matrix is invariant to the rest of values of P. This is a very important property when we are interested in finding the MLE estimator of the covariance matrix associated with an invariant subgraph, since it states that we can restrict our attention to the associated variables rather than to all the components of the random vector.
The simplest invariant subgraphs are the complete ones, but we are interested in more general ones. In particular, we will make use of complete separators to find these invariant subgraphs.
it is a separator of G and the subgraph of G induced by V is complete. The next theorem proves that, given a complete separator, two invariant subgraphs arise.

Theorem 3 Let X be a GMRF over G = (V , E) and A, B, C ⊂ V be pairwisely disjoint subsets satisfying A ∪ B ∪ C = V and A, B = ∅. If C separates A and B and the subgraph induced by C is complete, then X A∪C is a GMRF over the subgraph induced by A ∪ C and X B∪C is a GMRF over the subgraph induced by B ∪ C.
Proof We consider the global Markov property. For proving that X A∪C is a GMRF over the subgraph induced by A ∪ C, it suffices to prove that, if C ⊂ A ∪ C separates A and B in A ∪ C, then C separates A and B in G. Suppose that there exist a ∈ A and b ∈ B that are connected in G\C but are not connected in (A ∪ C)\C . We distinguish three cases: (i) If a , b ∈ C\C , since C is complete, then it holds that a and b are connected in (A ∪ C)\C , a contradiction. (ii) If a ∈ A\C and b ∈ C\C , then there exists a chain (a , ..., b ) between a and b of which at least an element is not contained in (A ∪ C)\C . Therefore, there exists b ∈ B with b ∈ (a , ..., b ). Since C separates A and B, the fact that b and a are connected implies that there exists c ∈ C\C connected to a . Since C\C is complete, c is also connected to b . The contradiction then follows from the fact that a and b are connected in A ∪ C. (iii) If a , b ∈ A\C , the contradiction is reached similarly to the previous case.
The proof for B ∪ C is identical to the one for A ∪ C.
It is concluded that, if we find a complete separator C that separates A, B = ∅, then the subgraphs induced by A ∪ C and B ∪ C are invariant subgraphs.

Example 2
The graph in Fig. 2

Complete separators of low order
Finding complete separators and invariant subgraphs is not always an easy task. However, it is easy to characterize complete separators of order 0 and 1. A complete separator of order 0 arises whenever the graph is partitioned into two connected components. Complete separators of order 1 are linked to Menger's Theorem.
The order of a minimal uv-separator of G is equal to the maximum number of internally disjoint chains that connect u and v in G.
As a consequence of this theorem, it is concluded that we can find some invariant subgraphs just by studying disjoint chains.

and G is maximal with respect to this property, then G is an invariant subgraph.
Proof If the maximal subgraph in which every two nodes are connected by two disjoint chains is the whole graph G, then the result follows. Otherwise, there exist two nodes connected by a maximum number of internally disjoint chains lower than 2. From Menger's Theorem, these two nodes are separated by another node or by the empty set. In both cases we have a complete separator, thus two invariant subgraphs, Theorem 3. We can repeat the same process on each of the obtained invariant subgraphs several times, until all obtained invariant subgraphs are such that for every u, v ∈ V there exist two disjoint chains from u to v in G and G is maximal with respect to this property.
Note that, if the order of all these subgraphs is greater than two, we obtain a decomposition into 2-connected components of the graph. It must be remarked that Menger's Theorem cannot be used to study a similar case associated with 3 disjoint chains because the minimal separator may be formed by two nodes that are not adjacent, so it might not be a complete separator. As a conclusion, we expect to find complete separators in sparse graphs.
We end this section by stressing that the existence of complete separators of order 0 and 1 is crucial when studying the complete separators of the Cartesian product of two graphs A B, that can been seen as a two-dimensional graph. This is the case for the ladder graph, which is considered in Sect. 6.1. In particular, Theorem 2.2 in Anand et al. (2012) can be reformulated as follows:

Theorem 5 Let A B be the Cartesian product of two graphs A and B. It holds that A B has a complete separator if and only if one of A and B is complete and the other
one has a complete separator of order 0 or 1.

Definition of the MGMRFs
From Theorem 3, we can find the elements of the submatrices Σ A∪C and Σ B∪C just by solving the GMRF construction problem over the associated subgraphs. The second step is then to find the rest of the values of Σ. Since for any a ∈ A and b ∈ B it holds that (a, b) / ∈ E, this can be done in an easy way by reducing the main problem to an MGMRF model. If C = ∅ we associate A, C, B with a 3-tree, being C associated with the node with degree of incidence 2 because X A and X B are conditionally independent given X C (since C separates A and B). Note that, if C = ∅, then the graph associated with the MGMRF only consists of two non-adjacent nodes, one for A and another one for B.
It is interesting to apply Theorem 3 iteratively whenever we identify a complete separator in A ∪C or in B ∪C. In such case, we may solve several GMRF construction problems over the obtained invariant subgraphs.

Solving MGMRF over forests
The next step is to solve the MGMRF construction problem defined above, ultimately resulting in a solution to the initial GMRF construction problem. We will focus on the case of trees, keeping in mind that, since the covariance matrix of an MGMRF over a forest is a diagonal block matrix where the blocks are associated with the connected components of the forest, the provided results are also valid for solving MGMRF construction problems over forests. Firstly, we provide a lemma concerning the MGMRF construction problem over a tree with 3 nodes.
Lemma 1 Let X = (X A , X C , X B ) be an MGMRF over a 3-tree G = (V , E) in which X C is associated with the node with degree of incidence 2. It holds that Σ AB = Σ AC Σ −1 C Σ C B . Proof By using the matrix formula to calculate the conditional distribution of (X A , X B ) given X C for a multivariate Gaussian distribution, we have that The result then follows from applying that Σ AB|C is a null matrix.
We can repeat the 3-tree structure all over the tree to calculate the remaining elements of Σ by just operating with the matrices obtained from the solved GMRF construction problems over the invariant subgraphs.

Proposition 5 Let X be an MGMRF over a tree G = (V , E). The covariance matrix between the subvectors associated with two non-adjacent nodes A and B is given by:
where K = (K 1 , ..., K l ) is the unique chain of length l − 1 from K 1 = A to K l = B.
Proof Since G is a tree, there exists a unique chain K = (K 1 , ..., K l ) from K 1 = A to K l = B (Kocay and Kreher 2016). The subvectors (X A , X K 2 , X K 3 ), (X A , X K 3 , X K 4 ), ..., (X A , X K l−1 , X B ) are MGMRFs over 3-trees. The result follows from applying l −2 times Lemma 1.

Corollary 1 Let X be a GMRF over a tree G = (V , E). Pearson's correlation coefficient between two variables is the product of Pearson's correlation coefficients of the adjacent variables in the unique chain that connects them.
The case of tree graphs, although is one of the simplest cases, is interesting due to its applicability. For instance, Gaussian Markov chains are examples of GMRFs over a tree graph, which are closely related to the Kalman Filter (Kalman 1960), which has been used in many real-life applications (see, e.g., Auger et al. 2013). Moreover, any GMRF over a non-acyclic graph can be approximated by a GMRF over a tree graph, which is specially relevant to spatial pattern classification and image restoration problems (Wu and Doerschuk 1995). Regarding the construction of the matrix F in Theorem 2, from the corollary above, it is immediate to determine the covariance between two variables in linear time O(n).

The proposed method and analysis of the computational complexity
As a conclusion of the previous results, a method to solve the GMRF construction problem by using invariant subgraphs may be defined. It can be summarized in the three following steps: 1. Decompose the graph into invariant subgraphs by using complete separators. 2. Solve the GMRF construction problem over every invariant subgraph. 3. Compute the solution for the whole graph.
For the second step, one of the original algorithms in Speed and Kiiveri (1986) and Wermuth and Scheidt (1977) is used to solve the GMRF construction problem over every invariant subgraph. These algorithms converge toward the solution by repeating a procedure iteratively. In each iteration, a loop over the non-adjacent nodes, the maximal cliques or the maximal cliques of the complementary graph, depending on the chosen algorithm, is performed. Fig. 1 and the initial matrix P defined by P i j = 0.5 if i = j and P i j = 1 if i = j, for any i, j ∈ {1, . . . , 9}. We recall the decomposition into invariant subgraphs of the graph provided in Fig. 2

Example 3 Consider the graph in
All other elements can be computed by applying Proposition 5:

Complexity of the proposed method
In order to study the complexity of the proposed method, it is necessary to examine the complexity of each one of the three sequential steps. In order to analyze the complexity we will assume that the number of adjacent nodes and the number of non-adjacent nodes grow with n 2 (see Tarjan 1985;Xu et al. 2011). This is a common assumption since the maximum number of adjacent nodes (and consequently the maximum number of non-adjacent nodes) in a graph of order n is n(n−1) 2 . The same is assumed for the number of maximal cliques of the graph, whose maximum number is the number of edges. This will lead to a common complexity for the first and second cyclic algorithm in Speed and Kiiveri (1986) and the variant in Wermuth and Scheidt (1977). The complexity of the three steps involves the complexity of inverting an n × n matrix. Its value depends on the method applied but is of the order O(n 2+ ), with > 0 (Xu et al. 2011).
The first step of our method is based on decomposing the graph into invariant subgraphs by using complete separators, which is a purely graph-theoretical matter. This decomposition in terms of complete separators was considered in a classical paper (Tarjan 1985), where an algorithm for obtaining this decomposition is proposed. The authors prove that the computation time is O(mn), being m the number of graph edges. In particular, their method is based on computing a minimal ordering of the graph (see Rose et al. 1976), and subsequently obtain the decomposition. More refined results in Coudert and Ducoffe (2018) reduce the computation time to the complexity of multiplying an n × n matrix; thus, the complexity of this step is O(n 2+ ).
For the second step, it suffices to apply one of the two cyclic algorithms in Speed and Kiiveri (1986) or the variant in Wermuth and Scheidt (1977) to all invariant subgraphs obtained in the previous step. Let us start by studying the complexity of solving the problem over the whole graph. In the three variants, an iteration of these algorithms visits all the cliques (second cyclic algorithm in Speed and Kiiveri 1986), the cliques of the complementary graph (first cyclic algorithm in Speed and Kiiveri 1986) or the edges of the complementary graph (algorithm in Wermuth and Scheidt 1977), which are of the order of m. For every element visited in the iteration, an inversion of an n × n matrix must be computed. Since the algorithm does not give the exact solution but converges toward this exact solution, the number of iterations may be determined by a maximum iteration number or a tolerance criterion. We conclude that the complexity of the algorithms is O(mn 2+ ). A similar study of the complexity has been carried out in Xu et al. (2011), reaching the same conclusions. For the algorithm over the invariant subgraphs, it is convenient to consider the first cyclic algorithm in Speed and Kiiveri (1986) or the variant in Wermuth and Scheidt (1977), since they focus on the complementary of the graph. It is easy to see that two nodes v, w ∈ V can be both contained in two different invariant subgraphs only if they are adjacent. Thus, considering an iteration over all the invariant subgraphs, the algorithms visit less cliques or edges in the complementary graph than if we consider the algorithm over the whole graph; thus, an upper bound will be m. In this case, we have to compute the inverse of the matrix associated with the invariant subgraph, with a maximum order of n. We conclude that the complexity for the second step has an upper bound of O(mn 2+ ).
For the last step, we need to invert the matrices associated with the complete separators and compute a matrix multiplication. A loose upper bound for the complexity can be obtained if we suppose that there are n complete separators of maximum order (which is actually not possible). Since any complete separator is always contained in an invariant subgraph, an upper bound of their order is n. In this case, the complexity has an upper bound of O(n 3+ ).
Thus, the complexity of each step of the sequential algorithm is, respectively, O(n 2+ ), O(mn 2+ ) and O(n 3+ ). Since m grows with n 2 , then the total complexity is O(mn 2+ ). Note that if we apply directly the algorithms in Speed and Kiiveri (1986) and Wermuth and Scheidt (1977), the complexity is also O(mn 2+ ). We conclude that the addition of the first and third steps to the method does not increase the complexity.

Further comments on computational aspects
On the one hand, if there exists a complete separator of the graph, the dimensions of the matrices in the second step are smaller than those for the original algorithm. Since this step is the most computationally expensive, a slight decrease in the dimension of the matrix can result in a notable decrease in the computation time of the algorithm. Thus, we expect the computation time to be reduced, even for a big graph in which we can only separate a few nodes (possibly in the corners or the outskirt of the graph) from the rest. The numerical results in Sects. 6.2 and 6.3 illustrate the improvement in the computation time for this case. In addition, the number of iterations in the algorithm used in step 2 is determined by a maximum number of iterations and a tolerance criterion. It seems reasonable that the number of iterations needed to reach the fixed tolerance criterion when we consider invariant subgraphs is smaller than in the general case.
On the other hand, if there does not exist any complete separator of the graph, no decomposition of the graph is reached and the second step consists in solving the GMRF construction problem over the whole graph, with the drawback of having computed step 1, losing time by trying to find complete separators that actually do not exist. This leads to an increase in the computation time with respect to the original algorithms. However, since the complexity of the first step is O(n 2+ ), which is smaller than the complexity of solving the problem over the whole graph, O(mn 2+ ), trying to find invariant subgraphs seems reasonable for graphs of large order.
We end this section by noticing that the decomposition step may be performed faster over some types of graphs. For instance, if the graph is planar then the complexity is O(n) and if the graph has a bounded treewidth then the complexity is O(n log(n)). We refer to (Sect. 6.1 in Coudert and Ducoffe 2018) for more information in this regard.

Examples
In this section, we consider three different scenarios in which the presented method simplifies the GMRF construction problem. Firstly, we consider a GMRF construction problem over the ladder graph. Secondly, we consider a generalization of the autoregressive model AR(k) in which the distribution is not necessarily stationary. Finally, we consider a real-life graph representing the (peninsular) regions of Spain as the nodes and in which two nodes are adjacent if the regions have a common border. We use real data concerning the mean temperature over the years 2011-2015 to find the Maximum Likelihood Estimation of the covariance matrix.

The ladder graph
Consider the random vector (X 1 , ..., X n , Y 1 , ..., Y n ) following a multivariate Gaussian distribution. We search for a model where, for any i ∈ {1, . . . , n}, X i is conditionally independent from all the other random variables but the three following ones: X i+1 , X i−1 and Y i . The same holds for Y i , being conditionally independent from all the other random variables but the three following ones: Y i+1 , Y i−1 and X i . We can express the multivariate Gaussian distribution as a GMRF over a ladder graph, see Fig. 3. This structure can be interpreted as a pair of random variables studied over a period of time and measured at specific instants of time in a way such that the values of the random variables only depend on the value of the variable at the previous instant of time and the value of the other variable at the same instant of time. We assume that the variance of all variables is equal to 1, the correlation between X i and X i−1 and between Y i and we obtain the correlation between Y i and X i+1 , denoted by φ(α, x), which is equal to the correlation between X i and Y i+1 . The value φ(α, x) is the only real solution of the equation φ(α, β) 3 − (α 2 + β 2 + 1)φ(α, β) + 2αβ = 0, obtained by inverting the covariance matrix and imposing the Markov properties. The following value is obtained: The computation of φ(α, β) is straightforward from the values of α and β. We may consider the MGMRF construction problem over the graph in Fig. 4, with A i = {X i , Y i } for any i ∈ {1, ..., n}. From Proposition 5 and after diagonalizing the matrix on the right side, it is easy to verify that, we obtain that the covariance matrix between A i and A i+d is:

Fig. 4 Graph associated with the MGMRF
This expression is easier to compute for any value of d than solving the GMRF construction problem over the whole graph, especially bearing in mind that the order of the graph is at least 2d.

Non-stationary AR(k) processes
In this example, we consider a generalization of the autoregression model AR(k). The AR(k) model considers a time series measured at specific instants of time that is stationary (i.e, it is invariant with respect to time translations) and such that the value at a certain instant of time only depends on the values of the k previous instants of time (Lindsey 2004). We consider a more general model, in which the time series is not necessarily required to be stationary. A non-stationary AR(k) model may be understood as a GMRF in which each node is adjacent to the k preceding nodes and the k succeeding nodes (or by all preceding/succeeding nodes in case there are less than k preceding/succeeding nodes). Notice that any subset of k consecutive nodes is a complete separator of the graph. Thus,all subsets {1,. . . ,k + 1},{2,. . . ,k + 2},. . . ,. . . ,n} are invariant subgraphs. For illustrative purposes, we consider a non-stationary AR(k) with n variables for all combinations of n ∈ {70, 110, 150} and k ∈ {3, 5, 10}. For setting the initial matrix, we choose a random matrix M of dimension n × n with elements ranging uniformly from 0 to 0.1. Subsequently, we compute the matrix P = M T + M + I n , where I n is the identity matrix of dimension n × n. The matrix M + M T is assured to be symmetric, whereas the addition of the term I n was incorporated in order to make the matrix M + M T positive definite.
As we have discussed earlier on, the matrix F in Theorem 2 maximizes the determinant among all matrices verifying F i j = P i j for any (i, j) ∈ E or i = j (Grone et al. 1984). Therefore, the evolution through time of the determinant of the matrix may be used as a tool for the comparison of the efficiency of the method. In particular, the faster the determinant increases, the faster we are approaching the solution to the GMRF construction problem.
We compute the determinant and the computation time (in seconds) of the algorithm presented by Wermuth and Scheidt (1977), both in case the decomposition into invariant subgraphs is considered and in case it is not considered for the aforementioned values of n and k. The result is illustrated in Fig. 5. It can be seen that for the method applied on the whole graph the determinant of the matrix converges smoothly toward the solution, whereas the here-presented method takes some initial time to decompose the graph into invariant subgraphs but attains the optimal solution immediately after this decomposition has been obtained. The difference between both methods seems to become larger for greater values of n.
In case the considered model is stationary, for the solution to the GMRF construction problem to be a stationary covariance matrix it is necessary for the initial matrix to be stationary. For this purpose, we may first estimate a stationary covariance matrix (see for instance Eldar et al. 2020) and then use the presented algorithm in order to impose the graph structure. Notice that in this case we can provide an explicit expression of the covariance between non-adjacent variables similarly as in Sect. 6.1.

A real-life graph
Let X be a multivariate Gaussian random vector associated with the mean temperature over a year of the (peninsular) regions of Spain. A map of the (peninsular) regions of Spain may be found on the left-hand side of Fig. 6. A reasonable assumption may be to consider that the mean temperature of a region only depends directly on the mean temperatures of the regions that share a border with it, just as in the example in Section 4.4.2 of Rue and Held (2005). The associated graph is presented on the right-hand side of Fig. 6.
Next, we consider the mean temperature through the years 2011-2015, as provided by the Instituto Nacional de Estadística (INE 2016), taking as reference the meteorological observatories located at the capital of the regions (with the exception of Extremadura, for which we have considered Badajoz). From these data, we construct the sample covariance matrix and solve the GMRF construction problem, both considering and not considering the decomposition into invariant subgraphs. The comparison of the computation time (in s) and the value of the determinant of the matrix at the current iteration for the classical algorithm in Wermuth and Scheidt (1977) and our proposed method is illustrated in Fig. 7. As can be observed, the convergence to the solution is faster in the case in which the graph is decomposed into invariant subgraphs.
We also want to stress that considering smaller regions (such as the provinces of Spain) makes it more difficult for complete separators to exist, mainly because in a planar graph any complete separator must have order 4 or less, as a consequence of Kuratowski's Theorem (Kuratowski 1930).

Conclusions
We have proposed a method to facilitate the search for the solution to the Gaussian Markov Random Field construction problem over a graph by independently considering several Gaussian Markov Random Field construction problems over invariant subgraphs and a Multivariate Gaussian Markov Random Field construction problem over forests. We have provided a way to find these invariant subgraphs by identifying complete separators of the graph, and we have studied the easiest cases in which the complete separators are of order 0 and 1. In particular, we have illustrated the utility of this proposal by considering a Gaussian Markov Random Field construction problem over a ladder graph, a non-stationary AR(k) process and a real-life problem regarding the mean temperature of the (peninsular) regions of Spain.
We have compared the complexity of the proposed method and the original algorithms in Speed and Kiiveri (1986) and Wermuth and Scheidt (1977), concluding that the decomposition of the graph in invariant subgraphs and the reconstruction of the solution for the initial problem do not increase the complexity. The numerical results illustrate that the proposed method has a lower computation time than the original algorithms for the considered examples. In case there does not exist a complete separator of the graph, both methods are equivalent, with the inconvenience of the here-presented method requiring to check the existence of complete separators at first. However, we believe that the cost of checking for the existence of complete separators is affordable, considering the potential gain. Furthermore, it should be borne in mind that the decomposition of the graph into complete separators only needs to be done once if several GMRF construction problems over the same graph are considered.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.