Efficient Localized Inference for Large Graphical Models

1 Introduction

Probabilistic graphical models such as Bayesian networks, Markov random fields, and conditional random fields are powerful tools for modeling complex dependencies over a large number of random variables ([Koller and Friedman, 2009; Wainwright et al., 2008]). Graphs are used to represent joint probability distributions, where nodes denote random variables, and edges represent dependency relationships between different nodes. With the specification of a graphical model, a fundamental problem is to calculate the marginal distributions of variables of interest. This problem is closely related to computing the partition function, or the normalization constant of a graphical model, which is known to be intractable and #P-complete. As a result, developing efficient approximation inference algorithms becomes a pressing need. The most popular algorithms include deterministic variational inference and Markov Chain Monte Carlo sampling.

However, many challenging practical problems involve very large graphs on which it is computationally expensive to use existing variational inference or Monte Carlo sampling algorithms. This happens, for example, when we use Markov random fields to represent the social network of Facebook or use a Bayesian network to model the knowledge graph that is derived from the entire Wikipedia, where in both cases the sizes of the graphical models can be prohibitively large

(e.g., millions or billions of variables). It is thus infeasible to perform traditional approximate inference such as message passing or Mote Carlo on these models because such methods need to traverse the entire model to make an inference. Despite the daunting sizes of large graphical models, in most real-world applications, users only want to make an inference on a set of query variables of interest. The distribution of a query variable is often only dependent on a small number of nearby variables in the graph. As a result, complete inference over the entire graph is not necessary and practical methods should perform inference only with the most relevant variables in local graph regions that are close to the query variables, while ignoring the variables that are weakly correlated and/or distantly located on the graph.

In this work, we develop a new localized inference method for very large graphical models. Our approach leverages the Dobrushin's comparison theorem that casts explicit bounds based on the correlation decay property in the graphs, in order to restrict the inference to a smaller local region that is sufficient for the inference of marginal distribution of the query variable. The use of the Dobrushin's comparison theorem allows us to explicitly bound the truncation error which guides the selection of localized region from the original large graph. Extensive experiments demonstrate both the effectiveness of our theoretical bounds and the accuracy of our inference algorithm on a variety of datasets.

RelatedWork Approximate inference algorithms of graphical models have been extensively studied in the past decades (see [Koller and Friedman, 2009; Wainwright et al., 2008; Dechter, 2013]). Query-specific inference, including [Chechetka and Guestrin, 2010] which proposed a focused belief propagation for query specific inference, and [Wick and McCallum, 2011; Shi et al., 2015] which study queryaware sampling algorithms, have recently been introduced for large graphical models. Compared with these methods, our work is theoretically motivated by the Dobrushin's comparison theorem and enables us to efficiently construct the localized region in a principled and practically efficient manner.

2 Background on Graphical Models

Graphical models provide a flexible framework for representing relationships between random variables [Heinemann and Globerson, 2014]. In graph G, we use X = (X1,X2, · · · ,Xn) to denote a finite collection of n random variables and we use x = (x1, x2, · · · , xn) to refer to an assignment. Suppose E is a set of edges and θ is a set of functions with θij(xi, xj) for edge 〈i j〉 ∈ E and θi(xi) for node i ∈ {1, 2, · · · , n}. We use µ(x; θ) to represent the joint distribution of the graphical model (G, E, θ) as following,

where Zµ(θ) is the normalization constant (also called the partition function).

In this work, we will focus on the Ising model, an extensively studied graphical model. The Ising model is a pairwise model with binary variables xi ∈ χ = {-1,+1}. The pairwise and singleton parameters are defined as follows

So the distribution of an Ising model is defined as,

Given a graphical model, marginal inference involves calculating the normalization constant, or the marginal probabilities of small subsets of variables. These problems require summation over an exponential number of configurations and are typically #P-hard in the worst case for loopy graphical models. However, practical problems can be often easier than the theoretically worst cases, and it is still possible to obtain efficient approximations by leveraging the special structures of given models. In this work, we focus on the query-specific inference, where the goal is to calculate the marginal distribution µ(xi) of given individual variable xi. For this task, it is possible to make good approximations based on a local region around xi, thus significantly accelerates the inference in very large graphical models.

3 Localized Inference and Correlation Decay

Given a large graphical model, it is usually not feasible to compute the exact marginal of a specific variable due to the exponential time complexity. Furthermore, it is even not practical to perform the variational approximation algorithms, such as mean field and belief propagation, when the graph is very large. This is because these traditional methods need to traverse the entire graph multiple times before convergence, and thus are prohibitively slow for very large models such as these built on social networks or knowledge bases.

On the other hand, it is relatively cheap to calculate exact or approximate marginals in small or medium size graphical models. In many applications, users are only interested in certain queries of node marginals. Because users' queries of interest often have strong associations with only a small number of nearby variables in the graph, the complete inference over the full graph is not necessary. This can be formally captured by the phenomenon of correlation decay, that is, when the graph G is large and sparse, the influence of a random variable on the distribution of another random variable decreases quickly as the distance of the shortest path between the corresponding nodes in the graph G increases.

The correlation decay property has been widely studied in statistical mechanics and graphical models [Rebeschini et al., 2015].

If a graphical model satisfies the property of correlation decay, it is possible that we can use only the local information in the graph to performmarginal inference, as the distant variables have little correlation with the query variable. This intuition allows us to use the information from the most relevant variables in the local region close to the queried variable to efficiently approximate its marginal distribution. Assume that µ is a large graphical model, and we want to calculate a marginal distribution µ(xi) of variable i. Localized inference constructs a much smaller model ν(xα), defined on a small subgraph α that includes i, such that ν(xi) ≈ µ(xi). The challenge here, however, is how to construct a good localized model and bound its approximation error. We address this problem via the Dobrushin's comparison theorem [F¨ollmer, 1982], and propose an efficient algorithm to find the local graph region for a given query node and provide an error bound between its approximate and true marginals. To get started, we first introduce the Dobrushin's comparison theorem, which is used to compare two Gibbs measures.

Theorem 1 [F¨ollmer, 1982] Dobrushin's comparison theorem Let µ be a Gibbs measure on a finite product space E = SI , where I is an index set. For i, j ∈ I, we define

where µi(xi|x) is the conditional distribution of the ith coordinate with respect to the σ-field generated by the coordinates with index j ̸= i, and ‖ · ‖ is the total variance distance. We compute

and assume c < 1. Let C = (Cij)i,j∈I and D = ∑∞n=0 C n = (I - C)-1 , then for any probability measure ν on the same place and any function f , we have

where b = (bj) is the singleton perturbation coefficient of node xj :

and δi(f) is the oscillation of f in the ith coordinate, that is,

In Theorem 1, µi(xi|x) is the probability of variable i conditioned on its adjacent variables whose assignments are the same as corresponding entries in x. According to the Markov property, calculating µi(xi|x) only requires information from the local star-shaped graph (the figure can be found in the full version). It is worth noting that a tighter bound can be obtained by defining b to be bj = 12 ∫ ‖µj (xj |x)-νj (xj |x)‖ν(dx). Here we use the definition in (3) for lower computational complexity. The matrix C is known as the Dobrushin's interaction matrix, and the inequality c = maxi∈I ∑j∈I Cij < 1 is the Dobrushin condition. If this condition holds, the theorem can give us a bound between two measures, which is the result of correlation decay.

In the following, we will apply Theorem 1 to undirected graphical models to derive an approximation bound of marginal distributions. We first denote I = {1, 2, · · · , n} by the index set of the variables and assume that we want to query the marginal distribution of variable xi. In order to apply Theorem 1, we set f(x) to be the indicator function of the variable xi, that is, f(x) = f(x1, x2, · · · , xn) = I[xi = k]. Then | ∫ fdµ- ∫ fdν| becomes the absolute marginal difference |µ(xi = k) - ν(xi = k)| between the two measures µ and ν. In addition, the oscillation of function f is thus reduced to δi(f) = 1 and ∀ j ̸= i, δj(f) = 0. With these simplifications, we obtain a bound of the maximum difference between marginals of the queried node i for two measures:

Corollary 1 Following the assumptions in Theorem 1 and the above text, we have

Note that the roles of µ and ν in (4) are not symmetric because the Dobrushin coefficient D is solely defined based on µ (and independent with ν). As a result, there are two ways to use bound (4) for localized inference, depending on whether we treat µ or ν as the original model that we want to query or the localized model that we use for approximation, respectively. We will next exploit both possibilities in the next sections. In Section 4, we take µ as the global model (or measure) and ν as the localized model (or measure) and derive a simple upper bound relates the approximation error to the distance between the query node and the boundary of the local region on the graph. In Section 5, we take ν as the global model (or measure) and µ as the localized model (or measure), we derive another upper bound that only involves the localized region, and leverage it to propose a greedy expansion algorithm to construct the localized model with guaranteed approximations.

4 Distance-based Upper Bound

In this section, we assume that µ = µ(x1, x2, · · · , xn) in Theorem 1 is defined by the original graphical model that we want to query, and ν is a simpler and more tractable distribution that we use to approximate the marginal of xi in µ.

For notational simplicity, we partition the node-set to two disjoint sets α and β = {1, 2, · · · , n} \ α, where α is the local subgraph that contains query node i ∈ α and β is the rest of the graph. We use ∂α and ˚α to represent the set of subscripts of nodes on the boundary and in the interior of α. Obviously, ∂α ⊆ α, ˚α ⊆ α, and ∂α ∪ ˚α = α. Similarly, ∂β ⊆ β, ˚ β ⊆ β, and ∂β ∪ ˚ β = β. In addition, we use xα to denote the variables in α and xβ to denote the variables in β. We will first apply the following lemma to obtain our first result on the relationship between the approximation error of marginals and the radius of the local subgraph α.

Lemma 1 [Rebeschini and van Handel, 2014] Assume I is a finite set and let m be a pseudo-metric on set I. C = (Cij)i,j∈I is a non-negative matrix. Suppose that

Then matrix D = ∞∑ n=0 C n = (I -C)-1 satisfies

In particular, this implies that

for every set J ⊆ I, where m(i, J) = maxj∈J m(i, j).

This lemma indicates that if Cij decays exponentially with the distance between i and j, the Dij , which is used in Theorem 1 and Corollary 1, also decays exponentially with the distance between i and j. The condition of this correlation decay lemma is usually mild in practice. When we choose m ≡ 0, which is naturally a pseudo-metric, and use Dobrushin's interaction matrix as C, the conditions of the lemma hold once the Dobrushin condition is satisfied, because matrix C in Theorem 1 is by definition a non-negative matrix and hence D is also non-negative, and every entry in b is less than 1/2. Applying Lemma 1, we can obtain the following result.

Theorem 2 Suppose µ is the probability measure for a graphical model for which we want to query the marginal distribution of node i. Let ν be the another probability on the same space, whose parameters of edges on subgraph α and parameters of nodes in ˚α are the same as µ. Assume the Dobrushin condition holds for µ (c = maxi∈I ∑j∈I Cij < 1). Let d(i, ∂α) denote the distance between node i and node-set ∂α on the Markov graph G of µ. If we assume

then ∀ ε > 0, we have

This theorem characterizes the error bound when approximating the global model µ using another model ν that matches µ locally in region α. Our result shows that in order to ensure an ε bound on the query node i, the distance d(i, ∂α) from the query node i to the boundary ∂α should be at least linear to ln(1/ε). In other words, the error ε decreases exponentially with d(i, ∂α). The proof of Theorem 2 can be found in the Appendix C.

As a result, given c and ε, we can get the minimum value of the lower bound of d(i, ∂α) by optimizing t. Theorem 2 gives a simple but general way to bound the local subgraph of variable xi, as we only need to check the Dobrushin condition and compute c on the whole true graphical model.

5 Localized Bound and Greedy Expansion

The bound in Theorem 2 requires computing the value c as defined in (2) for a given graphical model. However, since c is the maximum Cij of the entire graph, it can be very expensive to compute when the graph is large. In this section, we explore another approach of using the bound in Corollary 1, by setting ν = ν(x1, x2, · · · , xn) to be the distribution of the original graphical model and µ to be the localized model. In this way, we will derive a novel approximation approach by greedily constructing a local graph from the query variable i, with guaranteed upper bounds of the approximation error between marginal distributions of µ and ν.

To start with, we note that ν can be decomposed to

where ψα is the exponential of a potential function of xα, ψβ is the exponential of potential function of xβ, and ψ∂α∂β is the exponential of potential function defined on x∂α and x∂β.

We want to approximate ν with a simpler model µ in which the nodes in α and β are disconnected, so that the inference over i ∈ α can be performed locally within α, irrelevant to the nodes in α. Formally, we want to approximate µ by

which replaces the factor ψ∂α∂β with a product ˜ ψ∂α ˜ ψ∂β with approximations ˜ ψ∂α and ˜ ψ∂β. Therefore, the marginal distributions of xα and xβ get decoupled in µ, that is,

This decomposition thus allows us to approximately calculate marginal µ(xi) efficiently within subgraph α. The challenges here are 1) how to construct the factors ˜ ψ∂α and ˜ ψ∂β in µ to closely approximate ν, 2) how to decide the subgraph region and 3) how to bound the approximation error. We consider two methods for constructing ˜ ψ in this work:

1. [Dropping out] Simply remove the ψ∂α∂β in ν. To do so, we set

This corresponds to directly remove all the edges between ∂α and ∂β, which is also referred as the "dropping out" method in our experiments.

2. [Mean field] Find ˜ ψ∂α ˜ ψ∂β to closely approximate ψ∂α∂β by performing a mean field approximation, that is, we solve the following optimization problem:

where the KL(·||·) refers to the KL divergence of the corresponding normalized distributions. To apply the mean field approximation and reduce complexity, we further assume that the nodes are independent in ˜ ψ∂α and ˜ ψ∂β. By using the optimized approximation ˜ ψ∂α, we will be able to compensate the error of marginal of xi, which is introduced by simply removing the edges between ∂α and ∂β, as mentioned above.

Note that the potentials ˜ ψ∂β and ψβ in µ do not influence the calculation of µ(xi), for i ∈ α. Therefore, we remove all the edges in β. The marginal of node i ∈ α will not change.

By applying Corollary 1, we can now obtain an error bound which, remarkably, only involves the local region α. The details can be found in Appendix B.

Corollary 2 Assume µ(x) = µ(xα)µ(xβ), and the conditions in Theorem 1 holds, we have

where bj is defined in Eq 3, and D is defined by Dαα = (I - Cαα)-1 ; here C is defined in Theorem 1.

Note that the upper bound in (8) only involves the local region α and hence can be computed efficiently using mean field or belief propagation within the subgraph on α. The proof of Corollary 2 and the details on how to calculate C and D for Ising models in practice can be found in the full version.

Using the bound in (8), we propose a greedy algorithm to expand the local graph starting from query node i incrementally. At iteration, we add a neighboring node that yields the tightest bound using the above bound and repeat this process until the bound is tight enough or a maximum of graph size is reached. This process is summarized in Algorithm 1. After we complete the expanding phase, we can apply exact inference or on local region α to calculate the marginal of the query xi if the size of α is small or perform approximate inference methods if the size of α is medium. The actual size of α can vary in different graphical models, which is mainly determined by the correlation decay property near the query variable xi or the tightness of the upper bound in Eq (8).

Computational Complexity We find that both the complexity of expanding the local subgraph and the complexity of localized inference can be bounded. The details will be included in the Appendix A.

6 Experiments

We test our algorithm on both simulated and real-world datasets. The results indicate that our method provides an efficient localized inference technique.

Algorithm 1 Greedy expansion algorithm for localized inference

• 1: given a graphical model ν and a node i, approximate marginal probability ν(xi)
• 2: input: K =maximum number of nodes in the local subgraph α and δ = the improvement threshold
• 3: initialize local subgraph α = {i} and boundbest = 1
• 4: while |α| < K (|α| is the number of nodes in α) do
• 5: set β = {1, 2, · · · , n} \ α and ∂β the nodes in β that connects with α in ν.
• 6: for node k ∈ ∂β do
• 7: add node k to α and get a candidate local subgraph αnew k = α∪ {k}.
• 8: construct local model µ by setting ˜ ψ∂α = 1 (dropping out, Eq (6)) or estimating it using mean field as in Eq (7).
• 9: calculate the bound boundk in (8) where the boundk refers to ∑j∈∂αnew k Dij bj
• 10: end for
• 11: if mink∈∂β boundk < boundbest - δ then
• 12: update boundbest = mink∈∂β boundk
• 13: update α = α ∪ {argmink∈∂β boundk}
• 14: end if
• 15: end while

6.1 2D Ising Grid

In this section, we perform experiments on 2D-grid Ising models and regard the localized probability as ν and regard the true probability as µ. The graph is a 10×10 lattice and the coordinate of query node is (5, 5). The parameters in the Ising model is generated by drawing hi uniformly from [-I1, I1] for all nodes i and Jij uniformly from [-I2, I2] for all edges 〈i j〉. Here I1 and I2 control the locality and hardness of this Ising model.

Checking Dobrushin's condition We start with numerically checking the Dobrushin condition c = maxi ∑j Cij < 1. We find that c is smaller than one in most regions, but is larger than one when I2 is very large and I1 is very small, in which case the nodes are strongly coupled together (no correlation decay) and there is no significant local information. The hope, however, is that real problems tend to be easier because a large amount of information is available. The figure is shown in the full version.

Comparing Different Expansion Algorithms In this part, we compare the true approximation error maxk |µ(xi = k) - ν(xi = k)| to the bound ∑j∈∂α Dijbj given by our algorithm when we expand the local subgraph. The true error maxk |µ(xi = k) - ν(xi = k)| is evaluated using the brute-force algorithm. When removing the bipartite graph, we try both simply dropping edges and the mean field approximation. In all the experiments, we use the UGMMatlab package 1 for the mean field approximation.

In order to better compare the error, we also add two baselines. One baseline is that we expand the local subgraph in each step by randomly selecting a node in the boundary ∂β. Another baseline is that we expand the local subgraph greedily by choosing the node in ∂β that has the maximum L2 norm over the edge-set between such node and the subgraph α. The intuition is that when the the edges' weights are large, the node may be more related to the nodes in the subgraph.

1 http://www.cs.ubc.ca/ schmidtm/Software/UGM.html

In Figure 2, we compare our greedy expansion method stated in Algorithm 1 to the baselines stated above. For this experiment, we fix I1 = 1 and I2 = 0.25 and average on 100 random trials. We stop expanding the graph when the local subgraph contains 16 nodes. We calculate the mean value of the true errors and bounds in the 100 trials for a different number of nodes in the subgraph.

From Figure 2, we can see that, when combined the dropping out method for constructing ˜ ψ, our greedy expansion method significantly outperforms the two baselines. We also find that the mean field method for constructing ˜ ψ gives about the same true error as the dropping out method, but provides a tighter upper bound. It is interesting to note that the true errors of the two baseline expansion methods are sometimes even worse than the upper bounds of our greedy expansion, indicating the strong advantage of our method.

We further investigate how the parameters of the Ising model may influence the results of the algorithms and the tightness of the bound. For this purpose, we fix I2 = 0.25 and vary I1 in the range of {0, 0.5, 1, · · · , 10} in Figure 3. For each setting, we simulate 100 times and then calculate the mean error and bound. From Figure 3, we can find that the bound is again relatively tight, especially when the value of I1 is large. both the bounds and the true errors decrease as I1 increases because the correlation decay is stronger and the inference task is easier with strong local evidence on the singleton potentials (large I1).

6.2 Cora Data Set

We perform experimental evaluations on the Cora data set 2 . Cora consists of a large collection of machine learning papers with citation relations between the papers, in which each paper is labeled as one of seven classes. For our experiment, we binarize the labels by taking "Neural Networks" as label 1 and the remaining classes as label -1. We process the data by removing the hubs in the graph and truncate the graph to have a maximum degree of 15; this is done by randomly deleting edges of the nodes whose degree is larger than 15 until the whole graph is degree bounded by 15. We then experiment on the maximum connected subgraph, which consists of 2389 nodes and 4325 edges.

2 https://people.cs.umass.edu/ mccallum/data.html

In order to construct an Ising model, we random draw edge potentials by Jij ∼ N(0.25, 0.05) for each edge of the citation graph, and draw the singleton potentials by hi ∼ N(0.1I1, 1) for nodes with true label 1, and hi ∼ N(-0.1I1, 1) for nodes with true label -1. Here I1 is a parameter that we choose from {0, 1, · · · , 10}. When I1 , the node potentials increases so that marginal is more dominated by the status of the query node and the querying is easier.

Comparing local inference with global inference In this part, we want to compare the performance of inference on the local graph to the inference on the global graph. Since the global graph is too large, we can only use approximate inference algorithm. Here, we use mean field to do the global inference and use it as a baseline. For the local graph, we expand the graph greedily as stated in Algorithm 1 and choose a threshold of δ = 0.005 and stop expanding when the subgraph already has 16 nodes.

For I1 ∈ {0, 1, · · · , 10}, we query the same 500 nodes randomly selected out of the 2389 nodes and evaluate their marginal distributions. In global inference and local inference, we have the marginal on the each query node. If the marginal is larger than 0.5, we consider our inference algorithm give it label 1, whereas if the marginal is less than 0.5, we give it label -1.

In Figure 4, we report the accuracy of the labels given by the global and local inference evaluated w.r.t. the true labels, as well as the accuracy of the local inference evaluated w.r.t. the labels provided by global inference. We find that as I1 increases, both the accuracies of global and local inference w.r.t. the true labels increase significantly. In addition, the local inference gives similar result as the global inference (the green curve is high) and the accuracy increases as I1 increases as well. We also report in Figure 5 the precision, recall, and F-measure when comparing local inference with global inference, by treating label 1 as positives. Both figures show that when I2 is fixed and I1 increases, which means that the correlation decay is stronger, our local inference method achieves better results.

7 Conclusion

In this paper, we address query-specific marginal inference in large-scale graphical models using a new localized inference algorithm. We leverage the Dobrushin's comparison theorem to derive two error bounds for localized inference, including a simple bound based on graph distance and a localized bound from which we derive an efficient greedy expansion algorithm for constructing local regions for localized inference. Our experiments have shown that our bounds are practically useful and the algorithm works efficiently on various graphical models. Future directions include theoretical investigation on tighter bounds and development of more efficient greedy expansion algorithms.

A Computational Complexity

Here we consider the computational complexity of expanding the local subgraph and the complexity of localized inference. We always suppose that the maximum degree of the graph G is d and we define the maximum distance between the query node and any node in the subgraph to be the radius of the subgraph. First, given a threshold ε, from Theorem 2, we just need a subgraph with radius

where we recall that c is the Dobrushin coefficient c = maxi∈I ∑j∈I Cij. In particular, taking t = 2 shows that we just need r = ⌈- ln(ϵ(1 - c))/ ln 1+c 2c ⌉. It is worth noting that r decreases when c becomes small and/or the accuracy threshold ε becomes large. Since the size of the subgraph with radius r is no more than 1 + d + d2 + · · · + dr = (dr+1 - 1)/(d - 1), it can be much smaller than the whole graph. As a result, the inference over the subgraph is much more efficient.

Then, we discuss the computation complexity in each expansion step (Algorithm 1, line 6-10). We need to loop over the nodes in ∂β. In the loop, we need to calculate the vector b and the matrix D. The calculation for each element in b requires the enumeration of different assignments in the neighborhood of such node, which is bounded because it is not related to the size of the whole graph. In the calculation of matrix C, we only need to update the elements related to the new node. The number of such elements is no more than d and the calculation of each element is not related to the size of the whole graph. D can be derived from C and use historical information to calculate incrementally. The complexity is no more than computing the inverse of the whole matrix I -Cαα. If we use mean field approximation in the greedy expansion, the computation is also cheap because the sizes of ∂α and ∂β are small.

B Proof of Corollary 2

Proof 1 Note that the Dobrushin's interaction matrix C of µ is a block diagonal matrix. Since there are no edges between α and β, the corresponding blocks equal to zero. If the Dobrushin condition holds, D would also be a block-diagonal matrix and can be calculated easily from C. To see this, we have

Applying the bound in Corollary (1) gives the result.

C Proof of Theorem 2

Proof 2 Let m(i, j) = d(i, j) ln 1+(t-1)c tc , where d(i, j) represents the distance between node i and node j in the graph and t > 1. Then

Applying Lemma 1, we have

Substituting the inequality of d(i, ∂α) in the condition into right-hand side yields to the result.