Generalized belief propagation for probabilistic systems

[0001] This invention relates generally to modeling probabilistic systems, and more particularly, to modeling probabilistic systems using belief propagation in a Markov network.

[0002] Computer models are frequently used to study the behavior of complex probabilistic systems. When the systems contain many inter-dependent random variables, Markov networks are often used. In a Markov network, nodes of the network represent the possible states of a part of the system, and links between the nodes represent statistical dependencies between the possible states of those nodes.

[0003] By the Hammersly-Clifford theorem, from the study of Markov networks, the probability of any set of states at the nodes of the network can be written as the product of compatibility functions between clusters of nodes.

[0004] FIG. 1 shows a simple network with four nodes labeled a, b, c, and d. The links between the nodes represent the statistical dependencies between the possible states of the nodes. For the case of pairwise probabilistic interactions between nodes of the network, the overall joint probability of the system can be expressed as the product of compatibility functions for each linked pair of nodes:
P(s_a, s_b, s_c, s_d)=φ_ab(s_a, s_b)φ_bc(s_b, s_c)φ_ca(s_c, s_a)φ_bd(s_b, s_d),  [1]
where φ_ab is the compatibility matrix between nodes a and b, s_a is a random variable describing the state at node a, and similarly for the other nodes and variables.

[0005] Often, Markov networks for practical applications are very large. For example, an image acquired from a scene by a camera may be represented by a Markov network between all small neighboring patches, or even pixels, of the acquired image. Similarly, the well known “travelling salesman problem” can map onto a Markov network where the maximum probability state corresponds to the shortest path of the salesman's route. This network has as many nodes as cities to be visited. In some Markov networks, the nodes can represent measured input signals, such as visual input data. Markov models are also extensively used in speech recognition systems.

[0006] To analyze the probabilistic system modeled by a Markov network, one typically wants to find the marginal probabilities of certain network variables of interest. (The “marginal probability” of a variable signifies the probability of that variable ignoring, the state of any other network variable.) For example, it may be useful to examine the probability of a variable that represents an underlying explanation for some measured data, such as the probability of particular words used to vocalize particular speech sounds. To find those probabilities, the Markov network is marginalized over all the other variables in the network. This gives the probability of the variable representing the explanation, given the measured input data values. This marginalization is thus a form of inference.

[0007] One may also want to find states of the nodes, which maximize the network probabilities. For example, for the Markov network corresponding to the travelling salesman problem, it is desired to find the state at each node which maximize the probability of the Markov network. These states, which minimize the length of the salesman's route, are known as the maximum a posteriori probability (MAP) states.

[0008] In the example of FIG. 1, it is possible to determine the marginal probability P(s_a) of the variable at node a by summing the random values at nodes b, c, and d: $\begin{array}{cc}P\left(s_a\right)=\sum _{s_b,s_c,s_d}{\varphi }_{\mathrm{cb}}\left(s_a,s_b\right){\varphi }_{\mathrm{bc}}\left(s_b,s_c\right){\varphi }_{\mathrm{ca}}\left(s_c,s_a\right){\varphi }_{\mathrm{bd}}\left(s_b,s_d\right).& \left[2\right]\end{array}$

[0009] In general, especially for large networks, these marginal probabilities are infeasible to determine directly. The joint sum over all possible states of all the nodes can be of too high a dimension to sum numerically, particularly when the network has closed loops.

[0010] FIGS. 2a-b show examples of Markov networks with many loops for which it is difficult to find either the marginal probability at a node, or the state of the node which maximizes the overall probability of the Markov network. Both networks are in the form of lattices, which are commonly used to describe the joint probabilities of variables spatially, distributed over two dimensions. FIG. 2a shows a rectangular lattice, and FIG. 2b shows a triangular lattice. These type of lattice networks are used to model many systems.

[0011] Techniques to approximate the marginal probabilities for such structures are known, but these techniques are typically very slow. Simulated annealing can be used, or Gibbs sampling, see Geman et al. “Stochastic relaxation, Gibbs distribution, and the Bayesian restoration of images,” IEEE Pattern Analysis and Machine Intelligence, 6:721-741, 1984. Another class of approximation techniques are variational methods, see Jordan, “Learning in graphical models,” MIT Press, 1998. However, these methods require an appropriate class of variational approximation functions for a particular problem. It is not obvious which functions, out of all possible ones, to use for the approximations.

[0012] For the special case of Markov networks that form chains or trees, there is a local message-passing method that calculates the marginal probabilities at each node, see Pearl, “Probabilistic reasoning in intelligent systems: networks of plausible inference,” Morgan Kaufmann, 1988. The later method is now in widespread use, and is equivalent to the “forward-backward” and Viterbi methods for solving one dimensional Hidden Markov Models (HMM), and to Kalman filters and their generalization to trees, see Luettgen et al. in “Efficient multiscale regularization with applications to the computation of optical flow,” IEEE Trans. Image Processing, 3(1):41-64, 1994. This message-passing method gives the exact marginal probabilities for any Markov network that does not have loops. This is referred to as the “standard” belief propagation, or message-passing method below.

[0013] Unfortunately, many Markov networks of practical interest do contain loops. For example, an image, modeled as a Markov network of local image patches connected to their nearest neighbors, gives a lattice structured Markov network as shown in FIGS. 2a-b, also called a Markov random field. This type of network contains many loops.

[0014] Another method for inference in Markov networks applies the local message-passing rules derived for trees and chains in a network, even though the network may contain loops, see Weiss, “Belief propagation and revision in networks with loops,” Technical Report 1616, MIT AI Lab, 1997. This is referred to as the “loopy” belief propagation method in the description below, although it should be clearly understood that the “loopy” method is nothing more than the “standard” belief propagation method applied to a network with loops. When such a procedure converges, it can yield an approximate determination of the marginal probabilities. However, the loopy method sometimes gives too poor an approximation to the marginal probabilities, and often does not even converge. In the latter case, the approximation gives no single answer for the desired marginal probabilities.

[0015] Therefore, it is desired to provide a method for determining marginal probabilities in Markov networks that is both relatively fast and more accurate than the loopy method. Furthermore, it is desired to provide a method for networks with loops that converges more reliably than the prior art loopy belief propagation method.

[0016] The present invention provides a method for determining the probabilities of nodes in a network model of a complex probabilistic system. More particularly, the method determines desired marginal or maximum a posteriori probabilities in networks with loops. The method uses a message-passing scheme, appropriate for networks with loops, which is more accurate and typically converges more reliably and in fewer iterations than prior art loopy belief propagation methods.

[0017] The invention describes a class of procedures in which computational cost and accuracy can be traded off against each other, allowing a user of the invention to select more computation for more accuracy in a particular application of the invention.

[0018] The invention has two major advantages over the prior art loopy belief propagation method. First, the invented method normally gives much more accurate answers for the desired marginal probabilities. Second, the invented method can converge to a single answer in cases where the loopy method does not.

[0019] Instead of finding the marginal probability at a node, one embodiment of the invention finds the states at each node which approximately maximize the probability of the entire network. Thus, the invention provides a novel way to approximate both the marginal probabilities and MAP states in Markov networks.

[0020] Many Markov network problems of interest are known to be NP-hard problems. A problem is NP-hard when it is intrinsically harder than those problems that can be solved by a Turing machine in nondeterministic polynomial time. When a decision version of a combinatorial optimization problem belongs to the class of NP-complete problems, which includes the traveling salesman problem described above, an optimization version is NP-hard. The invented method yields fast, approximate solutions for some of these very difficult optimization problems.

[0021] More particularly, the invention provides a method for determining marginal probabilities of states of a system represented by a model. The model includes nodes connected by links. Each node represents possible states of a corresponding part of the system, and each link represents statistical dependencies between possible states of related nodes.

[0022] The nodes are grouped into arbitrary sized clusters such that every node is included in at least one cluster and each link is completely contained in at least one cluster. Based on the arbitrary sized cluster, messages are defined. Each message has associated sets of source nodes and destination nodes, a value for each state of the destination nodes and a rule that depends on other messages and on selected links connecting the source nodes and destination nodes. The values of each message are updated using the associated rule until a termination condition is reached, at which point the probabilities of the states of the system are determined from the values of the messages.

Introduction

[0043] Our invention builds on the “cluster variational method” for analyzing cooperative phenomena in statistical physics, see for example, the special issue in honor of R. Kikuchi, Progress in Theoretical Physics Supplements, vol. 115, 1994.

[0044] A physical system is often modeled in statistical mechanics as a Markov network, where the states of the nodes in the network correspond to the physical states of the particles or spins in the physical system. The cluster variational method derives approximations to the “Gibbs free energy” of a physical system. We refer to these as “Kikuchi approximations” to the Gibbs free energy, or the “Kikuchi free energy.” From a Kikuchi approximation to the Gibbs free energy, one can determine properties of a physical system, such as magnetization, specific heat, or the critical temperature for a magnetic spin system.

[0045] The simplest Kikuchi approximation derived with the cluster variational method corresponds to an approximation first proposed by Bethe, see Bethe, Proc. Roy. Soc. (London) A 150 (1935) 552. The loopy belief propagation method, when it converges, gives marginal probabilities identical to those obtained from the Bethe approximation.

[0046] Solutions obtained from propagation rules according to our invention are equivalent to solving the appropriate physical system under the corresponding Kikuchi approximation. The method of minimizing the Kikuchi free energy directly will therefore give identical results to our method, but such a direct minimization will be extremely slow and be plagued by the problem of convergence to inappropriate local minimal in the free energy surface. Thus, physics calculations made with more complicated Kikuchi approximations have typically only been tractable for homogeneous models of physical systems. As an advantage, our belief propagation method is also practical for heterogeneous models of physical systems where the nodes of the model have arbitrary topology and interdependencies.

Calculation of Marginal Probabilities by Kikuchi Free Energy Minimization

[0047] As shown in FIG. 3, we begin with a network 300 representing a system or model. Typically, the network is constructed ahead of time by an “expert” in the art field in which the model operates. It is not our job to construct the model; instead, we provide a method for using the model, even if the model is heterogeneous and includes many loops. For example, in a medical diagnosis system, the nodes can represent clinical findings, physiological and psychological conditions, laboratory results, pathologies, etc. Here, the goal may be to use the model to make the best (most probable) diagnosis given observations made for a particular patient. In an error-correcting coding system, the nodes represent source bits or parity bits, and the links represent the required relationships between the bits. Here, the goal may be to use the model to decode a block of bits that has been corrupted by noise.

[0048] The model includes a plurality of nodes 301. A random variable s at each node takes on one of a discrete set of states, or it may take on a continuum of states. The nodes are connected to others by links 302 to form the network model 300. Typically, a link connects two nodes, however higher-order links can connect more than two nodes. The network model 300 can represent a magnetic spin system, a medical diagnosis system, an image-processing system, a speech recognition system, a genetic regulatory network, a decoding system, or other complex real-world systems.

[0049] By definition, for a Markov network with pair-wise statistical dependencies between nodes, the probability of a particular assignment of values s at the nodes is given by: $\begin{array}{cc}P\left(s_1,s_2,\dots s_N\right)=\frac{1}{Z}\prod _{i,j}{\varphi }_{\mathrm{ij}}\left(s_i,s_j\right)\prod _i\psi \left(s_i\right),& \left[3\right]\end{array}$
where the first product runs over all the linked neighboring nodes, i and j.

[0050] The φ compatibility matrices represent the statistical dependencies between the nodes, as represented by the links 302 in the network model 300. For example, in a medical diagnosis system, one node might represent a symptom, and the nodes that it is linked to are related symptoms or diseases. The numerical values of the compatibility matrices φ, represent the strength of the relation between specific nodes. For example, in a medical diagnosis model, some symptoms might be strongly correlated with certain diseases, i.e. a strong statistical dependency, while other symptoms might only have weak statistical correlation with other symptoms or diseases.

[0051] For some systems, it may be desired to represent the statistical dependencies between the states of the nodes by compatibility functions (“links”) that are higher-order tensors. For example, if the states of nodes i, j, and k are mutually statistically dependent, but those dependencies cannot be expressed in terms of pair-wise interactions, then we introduce the tensor φ_ijk(s_i,s_j,s_k) to represent those dependencies. The over-all probability of the system would then also include the product over all higher-order tensors as well. The method we describe here will also apply to systems with such higher-order links between nodes. All that is necessary is to include terms that arise from the higher-order links in the “energy” and “region energies” that we describe below.

[0052] The ψ functions represent the “evidence” that each node is in any particular one of its states before we consider any of the other nodes. While the φ functions will normally never change, the ψ functions will typically change from one use of the model to another. For example, in the case of a medical diagnosis model, the ψ_i(s_i) function represents, for a specific patient, the results of a given test. To use the same model for a different patient, we use different evidence, represented through the ψ functions, but we do not change the φ functions, which represent the relationships between different symptoms, diseases, tests, etc.

[0053] In equation [3], Z is a normalization constant that insures that the sum of the probabilities of all possible states of the overall system is equal to one.

[0054] In order to describe the “Kikuchi free energy” of the model, it is helpful to define some additional functions, which are motivated by common conventions from statistical physics. We define the “bond strength” J_ij(s_i,s_j) between two nodes i and j by:
φ_ij(s_i,s_j)=e^{Jij (si,sj)/T}  [4]
where T is a parameter called the “temperature.” We define the “field” h_i(s_i) at the node i by:
ψ_i(s_i)=e^−hi(si)/T  [5]

[0055] In terms of these variables, we have: $\begin{array}{cc}P\left(s_1,s_2,\dots ,s_N\right)=\frac{1}{Z}{e}^{-E\left(s_1,s_2,\dots ,s_N\right)/T},& \left[6\right]\end{array}$
where E is the “energy”: $\begin{array}{cc}E\left(s_1,s_2,\dots ,s_N\right)=\sum _{\mathrm{ij}}{J}_{\mathrm{ij}}\left(s_i,s_j\right)+\sum _i{h}_i\left(s_i\right).& \left[7\right]\end{array}$

[0056] In statistical physics, equation [6], which is equivalent to the original definition [3] of the Markov model, is known as “Boltzmann's Law,” and the normalization constant Z is called the “partition function.” The partition function is defined explicitly by $\begin{array}{cc}Z=\sum _{s_1,s_2,\dots ,s_N}{e}^{-E\left(S_1,S_2,\dots ,S_N\right)/T}.& \left[8\right]\end{array}$

[0057] Following the conventions of statistical physics, we define a quantity called the “Gibbs free energy.” The Gibbs free energy G is a function of the probability distribution P(s₁,s₂, . . . ,s_N). By design, G is defined so that if it is minimized as a function P, then Boltzmann's Law is automatically obeyed. G is defined by: $\begin{array}{cc}\begin{array}{c}G\left(P\left(s_1,s_2,\dots ,s_N\right),\gamma \right)=U-TS+\\ \gamma \left[1-\sum _{s_1,s_2,\dots ,s_N}P\left(s_1,s_2,\dots ,s_N\right)\right].\end{array}& \left[9\right]\end{array}$

[0058] Here, U is the expected value of the energy, also called the “internal energy”: $\begin{array}{cc}U=\sum _{s_1,s_2,\dots ,s_N}E\left(s_1,s_2,\dots ,s_N\right)P\left(s_1,s_2,\dots ,s_N\right)& \left[10\right]\end{array}$
and S is the entropy: $\begin{array}{cc}S=-\sum _{s_1,s_2,\dots s_N}P\left(s_1,s_2,\dots ,s_N\right)\ln \left(P\left(s_1,s_2,\dots ,s_N\right)\right).& \left[11\right]\end{array}$

[0059] The last term in equation [9] is a Lagrange multiplier term designed to enforce the constraint that the sum of the probabilities over all states is equal to one. In order to insure that this constraint is satisfied, we set the partial derivative of G with respect to the Lagrange multiplier γ equal to zero.

[0060] If we differentiate the Gibbs free energy G with respect to the probability distribution P and the Lagrange multiplier y, and set those partial derivatives equal to zero, then we recover Boltzmann's Law. Thus, to understand physical systems, physicists often try to calculate and differentiate the Gibbs free energy.

[0061] For a large system, it is intractable to calculate G exactly, so Bethe, Kikuchi, and many others have developed approximation schemes. In these schemes, one determines exactly the Gibbs free energy over some small clusters of nodes of the network, and then approximates the total Gibbs free energy as a sum over many small clusters of nodes, correcting for over-counting in regions where the different clusters overlap. In the Bethe approximation, the clusters are precisely those pairs of nodes that are linked in the Markov network, while Kikuchi's more accurate approximations, also called the “cluster variational method,” allow for clusters of a larger size.

[0062] The cluster variational method errs by missing contributions to the Gibbs free energy that arise from high-order correlations between nodes that cannot be described in terms of the chosen clusters or intersections of clusters. The larger the size of the small clusters, the smaller the error.

[0063] In principle, the cluster variational method allows for arbitrary clusters on arbitrary Markov networks. In practice, however, most previous computations involving minimization of a Kikuchi free energy were restricted to symmetrically placed clusters over a homogenous network. Computations on inhomogeneous networks have been performed, but only for the Bethe approximation, see Morita, Physica 98A (1979) 566.

[0064] By “homogenous” networks, we mean networks for which the nodes are arranged in some regular repeated pattern, and the interactions between the nodes depends only on their position in that pattern. Many magnetic spin systems are “homogenous” in this sense, but most of the other types of systems that are modeled with Markov networks, such as the example of a medical diagnosis system, are not homogenous.

[0065] Previous prior art Kikuchi free energy minimization calculations were restricted to symmetrically placed clusters over a homogeneous network because, before our invention, there was no practical way to quickly minimize a Kikuchi free energy defined for arbitrary clusters over an inhomogeneous Markov network of arbitrary topology.

[0066] Using the Kikuchi approximate form for the total Gibbs free energy, we can approximately calculate the marginal probabilities of the variables at the nodes, and also the marginal probabilities of collections of nodes. We extend this technique to develop a message-passing method that converges to solutions corresponding to a Kikuchi approximation of the Gibbs free energy. A message-passage method, as an advantage over the prior art Kikuchi free energy calculations, is very much faster and more practical for inhomogeneous networks.

Kikuchi Free Energy Calculation for Simple Example

[0067] We use FIGS. 4a-c to illustrate our generalized message-based belief propagation method. This example network shows how we develop a message-passing method that has a corresponding Kikuchi approximation result as a fixed point of the messages. That is to say, our message-passing method will converge to a solution that is equivalent to the solution of the Kikuchi approximation.

[0068] FIG. 4a is a graph of many nodes grouped into three clusters. The clusters intersect each other. The nodes and their links are not shown, and indeed, the number of nodes per cluster can be quite large. In real world physical systems, one usually chooses many more than three clusters for a Kikuchi approximation. FIGS. 4b-c depict two possible instantiations of the three clusters and their intersections.

[0069] As described above with respect to FIG. 3, we group 310 the nodes of the network 300 into arbitrary intersecting clusters 311. A cluster is said to be intersecting if a single node appears in more than one cluster. Because we do an arbitrary clustering, we can focus on regions of the model that might be more significant to an exact solution. We require that every node is included in at least cluster. We also require that every link in the Markov network is completely included in at least one cluster. The first constraint ensures that our clusters represent and consider all nodes of the model, and the second constraint ensures that the clusters will intersect.

[0070] FIG. 4a shows abstractly three clusters of nodes (regions 1, 2, and 3), and their intersections (regions 12, 31, and 23), and intersections of intersections (region 123).

[0071] We use the following terminology. The largest regions are referred to as “clusters.” The clusters, along with their intersections, and the intersections of those intersections, and so on, are referred to collectively as “regions.” A “sub-region” of region r is a region that consists of a proper subset of the nodes in region r. Region r is a “super-region” of region s if and only if region s is a sub-region of region r.

[0072] The number of nodes in each cluster can vary, and be different from each other. As stated before larger clusters will yield better results than small clusters albeit at a greater computational cost. Note that region 12 is interpreted to include region 123, just as region 1 also includes region 123.

[0073] In equation [12] below, we approximate the total Gibbs free energy for the system 300 by the sum of the free energies of each of the clusters. Then, we include correction terms by subtracting the free energy of intersection regions where we over-counted the free energy. We may include additional correction terms for the over-counted intersections of the intersection regions, etc. In addition, we impose, with Lagrange multipliers, the additional constraints that the probabilities of each region sum to one, and that the probabilities have the correct marginalization relationships with respect to each other.

[0074] We express the Gibbs free energy in terms of the probabilities of each of the regions. We let α_i be a particular assignment of node values within the cluster region i. This is what we called s₁, s₂, . . . , s_N in equation [3] above. Our notation for region intersections is to concatenate the indices of the intersecting regions, so α_ijk is a particular assignment of node values within the region of intersection of regions i, j, and k.

[0075] The “region energy” E(α_i) includes those terms in the overall energy that can be attributed to nodes or links contained entirely in the region i, when the nodes in the region are set to the particular values denoted by α_i. Note that if there are higher-order tensor “links” in the definition of the system, the region energy will still contain those links that are contained.

[0076] The “region probability” P(α_i) is the marginal probability that the nodes in cluster i take on the particular values denoted by α_i. In our notation, the normalization constraint for region i is Σ_αiP(α_i)=1. That is, the sum over the probabilities of all possible value assignments to the nodes of region i is one. An analogous equation holds for all regions. These are the Lagrange multiplier constraints γ in equation [12], below.

[0077] The second set of constraints for the Gibbs free energy in equation [12] are the marginalization constraints, associated with Lagrange multipliers λ. The Lagrange multipliers impose consistency constraints between the probabilities of the regions and their intersecting sub-regions. The probability of a region, marginalized over all the nodes in the region not in some intersection region, must equal the marginal probability of the nodes in that intersection region. In our notation, Σ_αi\jP(α_i)=P(α_ij), where Σ_αi\j means the sum over all nodes in region i that are not in region j.

[0078] For the three clusters of nodes shown in FIGS. 4b-c, we express the Gibbs free energy in the Kikuchi approximation as G_1,2,3. We have three sets of U−ST terms, first a sum over states α_i for the cluster region free energies, second the terms involving α_ij where we subtract over-counted intersection regions, and a third term involving α₁₂₃ where we correct for over-subtracting that intersection region. After all these corrections are made, each subset of each region is counted exactly once, as it should be. $\begin{array}{cc}\begin{array}{c}{G}_{1,2,3}=\left(U-ST\right)\mathrm{terms}+\gamma \mathrm{constraints}+\lambda \mathrm{constraints}\\ =\sum _{{\alpha }_1}P\left({\alpha }_1\right)E\left({\alpha }_1\right)+\sum _{{\alpha }_2}P\left({\alpha }_2\right)E\left({\alpha }_2\right)+\sum _{{\alpha }_3}P\left({\alpha }_3\right)E\left({\alpha }_3\right)+\\ T\sum _{{\alpha }_1}P\left({\alpha }_1\right)\ln \left(P\left({\alpha }_1\right)\right)+T\sum _{{\alpha }_2}P\left({\alpha }_2\right)\ln \left(P\left({\alpha }_2\right)\right)+\\ T\sum _{{\alpha }_3}P\left({\alpha }_3\right)\ln \left(P\left({\alpha }_3\right)\right)-\sum _{{\alpha }_{12}}P\left({\alpha }_{12}\right)E\left({\alpha }_{12}\right)-\\ \sum _{{\alpha }_{23}}P\left({\alpha }_{23}\right)E\left({\alpha }_{23}\right)-\sum _{{\alpha }_{31}}P\left({\alpha }_{31}\right)E\left({\alpha }_{31}\right)-\\ T\sum _{{\alpha }_{12}}P\left({\alpha }_{12}\right)\ln \left(P\left({\alpha }_{12}\right)\right)-T\sum _{{\alpha }_{23}}P\left({\alpha }_{23}\right)\ln \left(P\left({\alpha }_{23}\right)\right)-\\ T\sum _{{\alpha }_{31}}P\left({\alpha }_{31}\right)\ln \left(P\left({\alpha }_{31}\right)\right)+T\sum _{{\alpha }_{123}}P\left({\alpha }_{123}\right)E\left({\alpha }_{123}\right)+\\ T\sum _{{\alpha }_{123}}P\left({\alpha }_{123}\right)\ln \left(P\left({\alpha }_{123}\right)\right)+{\gamma }_1\left[1-\sum _{{\alpha }_1}P\left({\alpha }_1\right)\right]+\\ {\gamma }_2\left[1-\sum _{{\alpha }_2}P\left({\alpha }_2\right)\right]+{\gamma }_3\left[1-\sum _{{\alpha }_3}P\left({\alpha }_3\right)\right]+\\ {\gamma }_{12}\left[1-\sum _{{\alpha }_{12}}P\left({\alpha }_{12}\right)\right]+{\gamma }_{23}\left[1-\sum _{{\alpha }_{23}}P\left({\alpha }_{23}\right)\right]+\\ {\gamma }_{31}\left[1-\sum _{{\alpha }_{31}}P\left({\alpha }_{31}\right)\right]+{\gamma }_{123}\left[1-\sum _{{\alpha }_{123}}P\left({\alpha }_{123}\right)\right]+\\ \sum _{{\alpha }_{12}}{\lambda }_{1\mathrm{\backslash 12}}\left({\alpha }_{12}\right)\left[P\left({\alpha }_{12}\right)-\sum _{{\alpha }_{1\mathrm{\backslash 12}}}P\left({\alpha }_1\right)\right]+\\ \sum _{{\alpha }_{12}}{\lambda }_{2\mathrm{\backslash 12}}\left({\alpha }_{12}\right)\left[P\left({\alpha }_{12}\right)-\sum _{{\alpha }_{2\mathrm{\backslash 12}}}P\left({\alpha }_2\right)\right]+\\ \sum _{{\alpha }_{23}}{\lambda }_{2\mathrm{\backslash 23}}\left({\alpha }_{23}\right)\left[P\left({\alpha }_{23}\right)-\sum _{{\alpha }_{2\mathrm{\backslash 23}}}P\left({\alpha }_2\right)\right]+\\ \sum _{{\alpha }_{23}}{\lambda }_{3\mathrm{\backslash 23}}\left({\alpha }_{23}\right)\left[P\left({\alpha }_{23}\right)-\sum _{{\alpha }_{3\mathrm{\backslash 23}}}P\left({\alpha }_3\right)\right]+\\ \sum _{{\alpha }_{31}}{\lambda }_{3\mathrm{\backslash 31}}\left({\alpha }_{31}\right)\left[P\left({\alpha }_{31}\right)-\sum _{{\alpha }_{3\mathrm{\backslash 31}}}P\left({\alpha }_3\right)\right]+\\ \sum _{{\alpha }_{31}}{\lambda }_{1\mathrm{\backslash 31}}\left({\alpha }_{31}\right)\left[P\left({\alpha }_{31}\right)-\sum _{{\alpha }_{1\mathrm{\backslash 31}}}P\left({\alpha }_1\right)\right]+\\ \sum _{{\alpha }_{123}}{\lambda }_{\mathrm{\backslash 12}\backslash 123}\left({\alpha }_{123}\right)\left[P\left({\alpha }_{123}\right)-\sum _{{\alpha }_{12\backslash 123}}P\left({\alpha }_{12}\right)\right]+\\ \sum _{{\alpha }_{123}}{\lambda }_{23\backslash 123}\left({\alpha }_{123}\right)\left[P\left({\alpha }_{123}\right)-\sum _{{\alpha }_{23\backslash 123}}P\left({\alpha }_{23}\right)\right]+\\ \sum _{{\alpha }_{123}}{\lambda }_{31\backslash 123}\left({\alpha }_{123}\right)\left[P\left({\alpha }_{123}\right)-\sum _{{\alpha }_{31\backslash 123}}P\left({\alpha }_{31}\right)\right].\end{array}& \left[12\right]\end{array}$

[0079] One point about the marginalization constraints. One might wonder why we do not also have terms like $\sum _{{\alpha }_{123}}{\lambda }_{1\backslash 123}\left({\alpha }_{123}\right)\left[P\left({\alpha }_{123}\right)-\sum _{{\alpha }_{1\backslash 123}}P\left({\alpha }_1\right)\right]$
which enforce that region 1 properly marginalizes down to region 123. We omit these terms because they are redundant with the marginalization equations already included in G_1,2,3.

[0080] Marginalization is a linear operator. Of the three possible constraints relating regions i, ij, and ijk, that region i marginalize down to region ij, that ij marginalize down to ijk, and that region i marginalize down to ijk, we only need two, and the third is linearly dependent on the other two equations.

[0081] The Kikuchi approximation to the Gibbs free energy can be generalized to the case of an arbitrary collection of regions in a relatively straightforward way, but we need to establish some additional more notation. Let r denote a region chosen from the set of all relevant regions R. Let S(r) denote the set of sub-regions of r and T(r) denote the set of super-regions of r. Define a “direct sub-region” of r as a sub-region that does not have a super-region that is also a sub-region of r. Thus, in our three cluster example, regions 12 and 13 are direct sub-regions of region 1, but region 123 is not. Let S_d(r) denote the set of direct sub-regions of r and T_d(r) denote the set of direct super-regions of r. Associate with each region r a “over-counting number” d_r, defined recursively by $\begin{array}{cc}{d}_r=1-\sum _{t\in T\left(r\right)}{d}_t.& \left[13\right]\end{array}$

[0082] The largest regions with no super-regions are defined to have a double-counting number of one. In the three-cluster case, we have d₁=d₂=d₃=1, d₁₂=d₂₃=d₃₁=−1, and d₁₂₃=1. The over-counting number is needed to insure that each link is ultimately counted exactly once.

[0083] The generalization of equation [12] for the Kikuchi approximation to the Gibbs free energy for the case of an arbitrary collection of regions is: $\begin{array}{cc}\begin{array}{c}G=\sum _{r\in R}[d_r\left(U_r-T{S}_r\right)+{\gamma }_r\left[1-\sum _{{\alpha }_r}P\left({\alpha }_r\right)\right]+\\ \sum _{s\in S_d\left(r\right)}\sum _{{\alpha }_s}{\lambda }_{r\backslash s}\left({\alpha }_s\right)\left[P\left({\alpha }_s\right)-\sum _{{\alpha }_{r\backslash s}}P\left({\alpha }_r\right)\right]],\end{array}& \left[14\right]\end{array}$
where U_r and S_r are the internal energy and entropy of the region r.

[0084] We want to find the cluster region probabilities that minimize the Kikuchi free energy while obeying the marginalization and normalization constraints. These probabilities are the correct marginal probabilities under the Kikuchi approximation.

[0085] To find equations for those probabilities, we differentiate equation [14] with respect to each of the region probabilities for each possible assignment of values to a region and set the result to zero. The resulting equations, illustrated for the three-cluster example of FIG. 4a, are:
E(α₁)+T+T ln P(α₁)+γ₁=λ_1\12(α₁₂+λ_1\31(α₃₁)  [15]
E(α₂)+T+T ln P(α₂)+γ₂=λ_2\12(α₁₂)+λ_2\23(α₂₃)  [16]
E(α₃)+T+T ln P(α₃)+γ₃=λ_3\23(α₂₃)+λ_2\31(α₃₁)  [17]
−E(α₁₂)−T−T ln P(α₁₂)+γ₁₂=−λ_1\12(α₁₂)−λ_2\12(α₁₂)+λ_12\123(α₁₂₃)  [18]
−E(α₂₃)−T−T ln P(α₂₃)+γ₂₃=−λ_2\23(α₂₃)−λ_3\23(α₂₃)+λ_23\123(α₁₂₃)  [19]
−E(α₃₁)−T−T ln P(α₃₁)+γ₃₁=−λ_1\31(α₃₁)−λ_3\31(α₃₁)+λ_31\123(α₁₂₃)  [20]
E(α₁₂₃)+T+T ln P(α₁₂₃)+γ₁₂₃=−λ_12\123(α₁₂₃)−λ_12\123(α₁₂₃)−λ_31\123(α₁₂₃)  [21]
Exponentiating both sides and rearranging terms of each equation, we have the following set of equations for the marginal posterior probabilities: $\begin{array}{cc}P\left({\alpha }_1\right)=k_1\exp \left(\frac{-E\left({\alpha }_1\right)}{T}\right)\exp \left(\frac{{\lambda }_{1/12}\left({\alpha }_{12}\right)+{\lambda }_{1\mathrm{\backslash 31}}\left({\alpha }_{31}\right)}{T}\right)& \left[22\right]\\ P\left({\alpha }_2\right)=k_2\exp \left(\frac{-E\left({\alpha }_2\right)}{T}\right)\exp \left(\frac{{\lambda }_{2/12}\left({\alpha }_{12}\right)+{\lambda }_{2\mathrm{\backslash 23}}\left({\alpha }_{23}\right)}{T}\right)& \left[23\right]\\ P\left({\alpha }_3\right)=k_3\exp \left(\frac{-E\left({\alpha }_3\right)}{T}\right)\exp \left(\frac{{\lambda }_{3/23}\left({\alpha }_{23}\right)+{\lambda }_{3\mathrm{\backslash 31}}\left({\alpha }_{31}\right)}{T}\right)& \left[24\right]\\ P\left({\alpha }_{12}\right)=k_{12}\exp \left(-\frac{E\left({\alpha }_{12}\right)}{T}\right)\exp \left(\frac{\begin{array}{c}{\lambda }_{1/12}\left({\alpha }_{12}\right)+{\lambda }_{2\mathrm{\backslash 12}}\left({\alpha }_{12}\right)-\\ {\lambda }_{12\backslash 123}\left({\alpha }_{123}\right)\end{array}}{T}\right)& \left[25\right]\\ P\left({\alpha }_{23}\right)=k_{23}\exp \left(-\frac{E\left({\alpha }_{23}\right)}{T}\right)\exp \left(\frac{\begin{array}{c}{\lambda }_{2/23}\left({\alpha }_{23}\right)+{\lambda }_{3\mathrm{\backslash 23}}\left({\alpha }_{23}\right)-\\ {\lambda }_{23\backslash 123}\left({\alpha }_{123}\right)\end{array}}{T}\right)& \left[26\right]\\ P\left({\alpha }_{31}\right)=k_{31}\exp \left(-\frac{E\left({\alpha }_{31}\right)}{T}\right)\exp \left(\frac{\begin{array}{c}{\lambda }_{3/31}\left({\alpha }_{31}\right)+{\lambda }_{1\mathrm{\backslash 31}}\left({\alpha }_{31}\right)-\\ {\lambda }_{31\backslash 123}\left({\alpha }_{123}\right)\end{array}}{T}\right)& \left[27\right]\\ P\left({\alpha }_{123}\right)=k_{123}\exp \left(-\frac{E\left({\alpha }_{123}\right)}{T}\right)\exp \left(\frac{\begin{array}{c}-{\lambda }_{12/123}\left({\alpha }_{123}\right)-\\ {\lambda }_{23\backslash 123}\left({\alpha }_{123}\right)-\\ {\lambda }_{31\backslash 123}\left({\alpha }_{123}\right)\end{array}}{T}\right),& \left[28\right]\end{array}$
where the k's are normalization constants, different for each equation, that normalize each marginal probability. Each of these equations is actually a collection of equations. Thus, if we focus on equation [22], we get one equation for each instantiation of α₁. On the right hand side of this equation, α₁₂ and α₃₁ are completely determined by α₁, because regions 12 and 31 are sub-regions of region 1.

[0086] Of course, we also need to differentiate the Gibb's free energy with respect to the λ and γ Lagrange multipliers. Doing that gives us the desired marginalization and normalization conditions for the region probabilities.

From Lagrange Multipliers to Messages

[0087] In the previous section, we derived a set of self-consistent equations for region probabilities and Lagrange multipliers by differentiating the Kikuchi approximation of the Gibbs free energy. These equations are a solution to the problem of finding marginal posterior probabilities at the level of the Kikuchi approximation.

[0088] In our new method, we develop equivalent and more convenient formulations of the same solution, in terms of a convergent “message-passing” method. By defining 320 the Lagrange multipliers, see step 320, in terms of another set of quantities, called “messages,” we obtain a different but equivalent set of equations for the cluster probabilities.

[0089] Our method works by treating the self-consistent equations as message update rules that are iterated until a termination condition is reached, for example, a fixed number of iterations, or the solution converges 330, at which point the marginal posterior probabilities 332 can be determined from the converged messages 340. Before the first update, we assign initial values 331 to the messages 321. The initial values can be random or special values, e.g., all random positive numbers, or all ones.

[0090] In general, there are many equivalent sets of message definitions and message-update rules that correspond to the same Kikuchi approximation. Some particular message definitions and message-update rules, which we describe later, have particularly nice properties and reduce to the already known ordinary message-passing methods when applied to the Bethe approximation. In this section, we want to describe the general properties that should be satisfied by messages so that any particular message-passing method is equivalent to the corresponding Kikuchi approximation.

[0091] We have used the marginalization relations to relate region probabilities to probabilities of their intersection regions. For the three-cluster example, equations [22]-[28] describe how the probabilities are related to the region energies, and the normalization and marginalization Lagrange multipliers. We ignore the normalization Lagrange multipliers (the γ's) for now, knowing that these simply give a multiplicative scalar which normalizes the final probabilities.

[0092] Our goal is to express messages in terms of the λ Lagrange multipliers. It will turn out that the marginalization constraint relations will then give us a set of fixed point equations for the messages, which we will interpret as the message update rules 321.

[0093] We stack all of the λ Lagrange multipliers to create a vector {right arrow over (λ)}. In this vector, there is one λ Lagrange multiplier for each linearly independent marginalization constraint, times the number of different possible variable assignments in the cluster region being marginalized down. For our three cluster example, the vector {right arrow over (μ)} includes nine sets of Lagrange multipliers: λ_1\12(α₁₂), λ_1\31(α₃₁), λ_2/12(α₁₂), λ_2/23(α₂₃), λ_3/23(α₂₃), λ_3/31(α₃₁), λ_12/123(α₁₂₃), λ_23/23(α₁₂₃), and λ_31/123(α₁₂₃), which we take to be stacked in that order. The total dimensionality of the vector {right arrow over (λ)} is D_λ=2D₁₂+2D₂₃+2D₃₁+3 D₁₂₃, where, for example, D₁₂ is the dimensionality of α₁₂.

[0094] In our example, equations [22-28] give the cluster probabilities in terms of the Lagrange multipliers {right arrow over (λ)}. We take the logarithm of all these equations, and define the local quantities

 L(α)=E(α)+T ln P(α)  [29]

for each cluster. Then we summarize these equations by:
{right arrow over (L)}=Aλ,  [30]
where {right arrow over (L)} is a vector of the L's in the seven regions taken in the order 1, 2, 3, 12, 23, 31, and 123, and {right arrow over (L)} has dimensionality:
D_L≡D₁+D₂+D₃+D₁₂+D₂₃+D₃₁+D₁₂₃.
A is a matrix with D_L rows and D_λ columns, but it has a special block form, which we express as: $\begin{array}{cc}A=\left(\begin{array}{ccccccccc}\tilde{1}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}\\ \tilde{0}& \tilde{0}& \tilde{1}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}\\ \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}\\ \tilde{1}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}& -\tilde{1}& \tilde{0}& \tilde{0}\\ \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{1}& \tilde{0}& \tilde{0}& -\tilde{1}& \tilde{0}\\ \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& -\tilde{1}\\ \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& -\tilde{1}& -\tilde{1}& -\tilde{1}\end{array}\right).& \left[31\right]\end{array}$

[0095] Each entry in this matrix actually denotes a sub-matrix. For example, the {tilde over (1)} in the upper left corner denotes a sub-matrix of D₁ rows and D₁₂ columns which is equal to one when the nodes in α₁₂ agree with the corresponding nodes in α₁, and is equal to zero otherwise. This becomes clear in a more specific example.

[0096] If region 1 includes two nodes S_a and S_b, and region 2 includes the single node S_b, and each node can be in two states that we label “1” and “2,” then α₁ can be in any of the four states S_a=1, S_b=1; S_a=1, S_b=2; S_a=2, S_b=1; S_a=2, S_b=2 and α₁₂ can be any of the two states S_b=1; S_b=2, and the sub-matrix {tilde over (1)} in the top left corner of A is expressed as: $\left(\begin{array}{cc}1& 0\\ 0& 1\\ 1& 0\\ 0& 1\end{array}\right).$

[0097] In general, the self-consistent equations derived from a Kikuchi approximation can be summarized in the form of a matrix equation {right arrow over (L)}=A{right arrow over (λ)} and a set of normalization and marginalization conditions.

[0098] We desire an alternative formulation for our solution in terms of “messages” and self-consistent message equations. We still use the same normalization and marginalization conditions on the cluster probabilities, but we redefine the Lagrange multipliers A in terms of messages.

[0099] We store the logarithms of all the messages between the regions in an ordered vector {right arrow over (m)}. To be commensurate with the Lagrange multipliers that we seek to summarize, we normally let the dimensionality of {right arrow over (m)} be the same as the dimensionality of {right arrow over (λ)}. In general, the vector {right arrow over (m)} can have different dimensionality than {right arrow over (λ)}, but the vector needs to span the same sub-space of constraints on the lattice so that we can re-write the Lagrange multiplier {right arrow over (λ)} in terms of messages {right arrow over (m)}.

[0100] Thus, for our example, we generate sets of logarithms of messages commensurate with {right arrow over (λ)} as m_1\12(α₁₂), m_1/31(α₃₁), m_2\12(α₁₂), m_2\23(α₂₃), m_3\23(α₂3), m_3\31(α₃₁), M_12\123(α₁₂₃), M_23\123(α₁₂₃), and m_31\123(α₁₂₃). The vector m stores these messages in that order.

[0101] There is a set of self-consistent equations for the cluster probabilities in terms of the messages where the equations have the matrix form:
L=B{right arrow over (m)},  [32]
or equivalently:
A{right arrow over (λ)}=B{right arrow over (m)},  [33]
where B is a “selection” matrix that determines the choice of the message-passing method. This matrix is specified below.

[0102] The matrix B and the vector {right arrow over (m)} constitute a valid reformulation of the Kikuchi solution as long as the system described by equation [33] has at least one solution that can be written in the form of a matrix equation:
{right arrow over (m)}=C{right arrow over (λ)}.  [34]

[0103] Any valid reformulation of the Kikuchi solution can be used as the basis for our message-passing method.

[0104] To illustrate this idea, we return to our three-cluster example of FIG. 4a, for which we have already chosen the vector {right arrow over (m)}. We can now choose many different matrices B, but we use one that arises from a scheme that we describe in more detail below.

[0105] Recall that the m's defined above are the logarithms of messages, which we denote by M's. We now set the marginal posterior probabilities of the regions to obey the set of equations: $\begin{array}{cc}P\left({\alpha }_1\right)=k_1\exp \left(-\frac{E\left({\alpha }_1\right)}{T}\right)M_{3\mathrm{\backslash 31}}\left({\alpha }_{31}\right)M_{2\mathrm{\backslash 12}}\left({\alpha }_{12}\right)M_{23\backslash 123}\left({\alpha }_{123}\right)& \left[35\right]\\ P\left({\alpha }_2\right)=k_2\exp \left(-\frac{E\left({\alpha }_2\right)}{T}\right)M_{1\mathrm{\backslash 12}}\left({\alpha }_{12}\right)M_{3\mathrm{\backslash 23}}\left({\alpha }_{23}\right)M_{31\backslash 123}\left({\alpha }_{123}\right)& \left[36\right]\\ P\left({\alpha }_3\right)=k_3\exp \left(-\frac{E\left({\alpha }_3\right)}{T}\right)M_{1\mathrm{\backslash 31}}\left({\alpha }_{31}\right)M_{2\mathrm{\backslash 23}}\left({\alpha }_{23}\right)M_{23\backslash 123}\left({\alpha }_{123}\right)& \left[37\right]\\ P\left({\alpha }_{12}\right)=k_{12}\exp \left(-\frac{E\left({\alpha }_{12}\right)}{T}\right)M_{1\mathrm{\backslash 12}}\left({\alpha }_{12}\right)M_{2\mathrm{\backslash 12}}\left({\alpha }_{12}\right)M_{23\backslash 123}\left({\alpha }_{123}\right)M_{31\backslash 123}\left({\alpha }_{123}\right)& \left[38\right]\\ P\left({\alpha }_{23}\right)=k_{23}\exp \left(-\frac{E\left({\alpha }_{23}\right)}{T}\right)M_{3\mathrm{\backslash 23}}\left({\alpha }_{23}\right)M_{3\mathrm{\backslash 23}}\left({\alpha }_{23}\right)M_{12\backslash 123}\left({\alpha }_{123}\right)M_{31\backslash 123}\left({\alpha }_{123}\right)& \left[39\right]\\ P\left({\alpha }_{31}\right)=k_{31}\exp \left(-\frac{E\left({\alpha }_{31}\right)}{T}\right)M_{1\mathrm{\backslash 31}}\left({\alpha }_{31}\right)M_{3\mathrm{\backslash 31}}\left({\alpha }_{31}\right)M_{123\backslash 123}\left({\alpha }_{123}\right)M_{23\backslash 123}\left({\alpha }_{123}\right)& \left[40\right]\\ p\left({\alpha }_{123}\right)=k_{123}\exp \left(-\frac{E\left({\alpha }_{123}\right)}{T}\right)M_{12\backslash 123}\left({\alpha }_{23}\right)M_{23\backslash 123}\left({\alpha }_{123}\right)M_{31\backslash 123}\left({\alpha }_{123}\right),& \left[41\right]\end{array}$
where again the k's are normalization constants that vary for each equation, as for equations [22-28] above.

[0106] Taking the logarithm of these equations, and recalling the order of the m's given above, we see that equations [35-41] are equivalent to the choice $\begin{array}{cc}B=\left(\begin{array}{ccccccccc}\tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{1}& \tilde{0}\\ \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}\\ \tilde{0}& \tilde{1}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}\\ \tilde{1}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{1}\\ \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{1}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{1}\\ \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{1}& \tilde{1}& \tilde{0}\\ \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{1}& \tilde{1}\end{array}\right).& \left[42\right]\end{array}$

[0107] The choice of B and {right arrow over (m)} give an equivalent reformulation of the region probabilities when there is a solution to the equation:
{right arrow over (m)}=C{right arrow over (λ)}  [43]
where
C=B⁺A,  [44]
and B⁺ is the pseudo-inverse of B. Of course, we can always compute a C matrix as above, but for the matrix to represent a true solution, we must verify that matrix equation
BC=A  [45]
is satisfied.

[0108] If a solution to A{right arrow over (λ)}=B{right arrow over (m)} exists for a particular choice of matrix B and message vectors {right arrow over (m)}, we discard the original Lagrange multipliers λ and work entirely with the messages. In this case, we guarantee that the self-consistent equations for the messages will ultimately give the same cluster probabilities as those given by the Kikuchi approximation.

[0109] The next step is to determine the self-consistent equations for the messages. To do this, we need to combine the equations for the cluster probabilities in terms of messages with the marginalization conditions on the cluster probabilities. In terms of our example, we use, for instance, the equation P(α₁₂)=Σ_α1\12P(α₁) along with the equations for P(α₁) and P(α₁₂) in terms of messages to obtain $\begin{array}{cc}{M}_{1\mathrm{\backslash 12}}{M}_{3\backslash 123}=k\sum _{{\alpha }_{1\mathrm{\backslash 12}}}\exp \left(\frac{E\left({\alpha }_{12}\right)-E\left({\alpha }_1\right)}{T}\right)M_{3\mathrm{\backslash 31}},& \left[46\right]\end{array}$
where k is again a normalization constant. In all, there will be one such equation for each marginalization constraint.

[0110] We use the D_λ equations of the form of equation [46] in our message-passing method, by interpreting them as message-update rules. In this interpretation, any messages on the right-hand-side of an equation are treated as the messages at iteration t, while the messages on the left-hand-side of the equation are the desired messages at iteration t+1. We start, for example, with random messages 331, and at each iteration replace the “old” messages on the right hand side of the equation with “new” messages from the left-hand side. In practice, the messages usually converge to a fixed point. When they converge, the fixed point corresponds to the solution in which we are interested.

[0111] As for all methods of this type, we can update the messages synchronously or asynchronously, and we can update the messages with some linear combination of the “old” and “new” messages. Our method is an iterative method for the solution of nonlinear equations, and as such, its convergence properties may be improved by using standard prior-art techniques for such methods, see “Iterative Methods for Linear and Nonlinear Equations”, C. T. Kelley, 1995.

[0112] Regarding the normalization constants k, we could continue to set these normalization constants by requiring that the total cluster probabilities sum to one. However, another procedure that works equally well, and is easy to implement, temporarily scales the messages by requiring that the sum of all the components of a message to equal a given arbitrary constant while iterating the self-consistent equations. Reasonable choices for the constant are one, or the dimensionality of the message. Thus, for instance, we can require that: $\begin{array}{cc}\sum _{{\alpha }_{12}}{M}_{1\mathrm{\backslash 12}}\left({\alpha }_{12}\right)=1.& \left[47\right]\end{array}$

[0113] After iterating the message equations to convergence, we re-scale the messages to guarantee that the region probabilities are properly normalized.

[0114] We have described the messages as being associated with a region r and one of its direct sub-regions s. We can equivalently describe the message as going from a set of “source” nodes to a set of “destination” nodes. The “source” nodes are those nodes in region r that are not in direct sub-region s, while the “destination” nodes are those nodes in s.

[0115] To obtain probability functions defined on a subset of the nodes in a region, one simply marginalizes the region probability by summing it over those nodes that are in the region but not in the subset of interest.

[0116] We can also determine the maximum a posteriori probability values of the nodes in a region rather than the region probabilities. For the prior-art “standard” belief propagation method, a simple modification suffices to obtain the approximate MAP values of the nodes. The necessary modification is to replace all sums over node values in the message-update rules with the value of the maximum, i.e., the “argmax,” over the same node values. The same modification can also be used to obtain the MAP values of regions for our method. Alternatively, we can use the physics formulation of the Markov network, and set the temperature T to a very small positive constant. In the limit T→0, we recover region probabilities that are entirely dominated by the MAP values of the nodes.

[0117] The message-update equations for some arbitrary B matrix is not, in general, convenient to solve. With an arbitrary B matrix, one might find that a message between two nodes on one side of the network would depend on messages between other nodes that are very far away. This is an undesirable property, both because it seems unintuitive, and because it means that solving the message-update equations will involve matrix multiplications with very large matrices. Finding the free energy for arbitrary B matrices also involve large matrix multiplications or inversions.

[0118] In order that iterating the message update equations be an efficient operation, we desire that the matrix B is sparse and local, i.e., only a few non-zero elements for each row or column, and only nearby messages influence each other.

[0119] In the next section, we describe a “canonical” message-passing method that yields these desired properties.

A Canonical Method

[0120] In the previous section, we described the properties that must be satisfied by messages and their associated message-update rules, which we call the “message-passing method,” in order that the messages give results equivalent to the Kikuchi approximation. In general, there may be many different message-passing methods that are all equivalent to the same Kikuchi approximation. In this section, we describe how to select one particular message-passing method that we call the “canonical” method.

[0121] Our canonical method gives the same message-update rules as the standard loopy message-passing method when applied at the level of the Bethe approximation. More generally, and as an advantage, our canonical method leads to local message-update rules that are relatively easy and convenient to solve because they correspond to a sparse B matrix.

[0122] We start with a Markov network with local compatibility functions φ_ij and evidence functions ψ_i. The overall probability of a configuration of states s_a, s_b, . . . is expressed as. $\begin{array}{cc}P\left(s_a,s_b,\dots \right)=\frac{1}{Z}\prod _{i,j}{\varphi }_{\mathrm{ij}}\left(s_i,s_j\right)\prod _i{\psi }_i\left(s_i\right),& \left[48\right]\end{array}$
where the first product runs over all linked neighbors, i and j.

[0123] Using FIG. 5, we illustrate this procedure with a specific example. In this example, nine nodes labeled a through i, are arranged in a 3×3 lattice 500. The only statistical dependencies or links φ_ij(s_i,s_j) included are the ones connecting nearest neighbors in the lattice.

General Procedure

[0124] The detailed steps of a general method is described with respect to FIG. 6. In step 601, a network 600 of nodes and links is grouped 601 into clusters or “super-nodes” 650 to be used in the approximation. These clusters can have a simple geometric shape, repeated over the lattice. Alternatively, the clusters can be triangles or squares of an irregular network, or the clusters can be some other arbitrary choice of node clusters best suited for the approximation. The set of clusters must cover all the links in the Markov network, and the clusters must be small enough so that we are able to perform the generalized belief propagation operations described below.

[0125] In our example, we group the network 500 into clusters: (abde), (bcef), (degh), and (efhi) 501-504 of FIG. 5. This set of clusters “covers” the lattice 500. That is, every link between nodes which influence each other is completely included in at least one of the four clusters we have selected. Of course, we can select a different collection of clusters which completely covers the lattice. That would just be a different Kikuchi approximation.

[0126] As stated before, the size of the clusters is proportional to the accuracy of the approximation. If we chose a cluster (abcdefghi), i.e., all the nodes in the system, then our approximation is exact. While that grouping may be feasible with such a small example, in general, we will perform summations whose number grows exponentially with the size of the cluster. Therefore, in a practical application, the number of nodes in a cluster is limited.

[0127] The prior art Bethe approximation is obtained when we select as the clusters all the pairs of nodes that influence each other. In our example, that would mean partitioning the nodes pairwise into clusters (ab), (bc), (de), (ef), (gh), (hi), (ad), (be), (cf), (dg), (eh), and (fi). Because our four clusters are larger than the twelve clusters of the Bethe approximation, we obtain more accurate results.

[0128] In step 602, we identify a collection of relevant regions 651 by recursively generating intersections of clusters, for our example, the regions, (be), (de) (fe), (he), and (e). We discard duplicate regions.

[0129] In step 603, we arrange the regions into a top-to-bottom hierarchy 652 of their intersections. The hierarchy for our specific example is shown in more detail in FIG. 7. This hierarchy has the property that each region is connected only to its “direct sub-regions” and “direct super-regions.” A sub-region s is a “direct sub-region” of another region r if there are no other regions in our collection that are super-regions of s and sub-regions of r. In our example, regions (be) and (de) are both direct sub-regions of region (abde), but region (e) is not. Region r is a direct super-region of region s if and only if region s is a direct sub-region of region r.

[0130] In step 604, we identify messages M 653 which connect all regions with their direct sub-regions. In our example, we would identify messages M_abde\be, M_abde\de, M_bcef\be, M_bcef\fe, M_degh\de, M_{degh he}, M_efhi\fe, M_efhi\he, M_be\e, M_de\e, M_fe\e, and M_he\e.

[0131] We can identify these messages using a more compact notation that we illustrate with one example:
M_be^ad≡M_abde\be.  [47]

[0132] For this message, we send a message from the upper nodes a and d to the lower nodes b and e. We can consider the upper nodes to be the “source” nodes for the message, and the lower nodes to be the “destination” nodes. Each of the messages between a region and its direct sub-region has dimensionality corresponding to the dimensionality of the sub-region. Thus, M_be^ad has the same dimensionality as region (be). When we want to emphasize that fact, we express the message as M_be^ad(S_b, s_e).

[0133] In step 605, we construct a probability function P(α_R) and a local potential Φ(α_R) for each region R. We construct a vector {right arrow over (L)} out of the quantities L(α)≡ln P(α)−ln Φ(α). For our example, we need to construct probabilities P(s_a,s_b,s_d,s_e), P(s_b,s_c,s_e,s_f), P(s_d,s_e,s_g,s_h), P(s_e,s_f,s_h,s_i), P(s_b,s_e), P(s_d,s_e), P(s_f,s_e), P(s_h,s_e) and P(s_e). The local potentials Φ(α_R) are constructed by combining all the original potentials φ_ij and ψ_i that belong to the region R. Thus, for example,
Φ(s_a,s_b,s_d,s_e)=φ_ab(s_a,s_b)φ_ad(s_a,s_d) φ_be(s_b,s_e)φ_de(s_d,s_e)ψ_a(s_a)ψ_b(s_b)ψ_d(s_d)ψ_e(s_e),  [48]
while Φ(s_b,s_e)=φ_be(s_b,s_e)ψ_b(s_b)ψ_e(s_e) and Φ(s_e)=ψ(s_e)

[0134] Each probability function has the dimensionality of the region it describes. For example, if node b can be in any of three states and node e can be in any of four states, then P(α_be) is described by 3*4=12 numbers.

[0135] In step 606, we specify the matrix B 656 that determines the selection of the message-passing method. This matrix will operate on the vector {right arrow over (m)} created from the logarithms of the messages specified in step 604. The matrix B has a block form that we explained when describing equation [33]. This special block form “hides” the details of the different possible values of the nodes within each message or region. Using this block form, the matrix B 655 has N_R rows and N_M columns, where N_R is the number of regions, and N_M is the number of messages. Each entry in the matrix is either a {tilde over (0)} matrix or a {tilde over (1)} matrix. We determine the entries of the matrix B by the following sub-step 606 for each entry:

[0136] (a) Label the region corresponding to the row of the entry R_R, and the region and sub-region of the message corresponding to the column of the entry R_C and S_C respectively. By R_C\S_C, we denote the set of nodes that are in R_C but not in S_C.

[0137] (b) If the nodes in SC are the same as the nodes in R_R, or if the nodes are a sub-set of the nodes in R_R, and if all the nodes in R_C\S_C are not in R_R, then the entry is a {tilde over (1)} matrix. Otherwise, the entry is a {tilde over (0)} matrix.

[0138] Explained more intuitively, this rule says that every region gets messages that start at source nodes that are outside of the region and end at destination nodes that are inside of the region. In our example, assuming the regions are listed in the order (abde), (bcef), (degh), (efhi), (be), (de), (fe), (he), (e), and the messages are ordered as we listed them above, the matrix B is: $\begin{array}{cc}B=\left(\begin{array}{cccccccccccc}\tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{1}\\ \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}\\ \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{1}& \tilde{0}& \tilde{1}& \tilde{0}\\ \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{1}& \tilde{0}& \tilde{0}\\ \tilde{1}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{1}& \tilde{1}\\ \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{1}& \tilde{1}\\ \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{1}& \tilde{1}& \tilde{0}& \tilde{1}\\ \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{0}& \tilde{1}& \tilde{1}& \tilde{1}& \tilde{1}& \tilde{0}\\ \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{0}& \tilde{1}& \tilde{1}& \tilde{1}& \tilde{1}\end{array}\right).& \left[49\right]\end{array}$

[0139] We explain one of these entries: the {tilde over (1)} matrix in the first row, third column. The first row corresponds to the region (abde), while the third column corresponds to the message M_be^cf. Both the nodes c and f are outside of the region (abde) and both the nodes b and e are inside the region (abde), so the conditions for choosing a {tilde over (1)} matrix are met for this entry.

[0140] One important advantage of this rule is that when it is applied to the messages and regions obtained from the Bethe approximation, the rule will always result in a matrix B that will ultimately yield the standard “loopy” belief propagation. Another important advantage is that this rule will normally give back a matrix B that results in a system of equations A{right arrow over (λ)}==B{right arrow over (m)} that is solvable. Another advantage is that this rule will result in a B matrix that is sparse and local.

[0141] In step 607, we exponentiate the equations represented by the matrix equation {right arrow over (L)}=B{right arrow over (m)} to obtain the equations 656 relating the region probabilities to the messages. For our example, we obtain nine equations from the nine rows of the matrix equation. As an example, we show how to construct the equation for the first row. The first row of the B matrix corresponds to the region (abde). Reading across the first row, we see that there are {tilde over (1)} sub-matrices in the third, fifth, eleventh, and twelfth columns, which correspond to the messages M_be^cf, M_de^gh, M_e^f, and M_e^h Therefore, the equation for the first row is:
P(s_a, s_b,s_d,s_e)=Φ(s_a,s_b,s_d,s_e)M_be^cf(s_b,s_e)M_de^gh(s_d,s_e)M_e^f(s_e)M_e^h(s_e).  [50]
The other equations are constructed in the same way. For example, the equation for the fifth row is
P(i s_b,s_e)=Φ(i s_b,s_e)M_be^ad(i s_b,s_e)M_be^cf(s_b,s_e)M_e^d(s_e)M_e^f(s_e)M_e^h(s_e)  [51]
In step 608, we express the N_m constraint equations on the region probabilities, and substitute in equations relating the region probabilities to the messages in order to obtain a set of self-consistent equations M_c (update rules) for messages 657.

[0142] In our example, we have twelve constraint equations. In the canonical method, there are always exactly as many constraint equations as there are messages. One of the constraint equations in our example relates P(s_a,s_b,s_d,s_e) to P(s_b,s_e): $\begin{array}{cc}P\left(s_b,s_e\right)=\sum _{s_b,s_e}P\left(s_a,s_b,s_d,s_e\right).& \left[52\right]\end{array}$

[0143] Using the results from the previous step, we obtain the self-consistent equation for the messages from the marginalization (abde)→(be): $\begin{array}{cc}{M}_{\mathrm{be}}^{\mathrm{ad}}\left(s_b,s_e\right)M_e^d\left(s_e\right)=\sum _{s_a{s}_d}{\varphi }_{\mathrm{ab}}\left(s_a,s_b\right){\varphi }_{\mathrm{ad}}\left(s_a,s_d\right){\varphi }_{\mathrm{be}}\left(s_b,s_e\right){\psi }_a\left(s_a\right){\psi }_d\left(s_d\right)M_{\mathrm{de}}^{\mathrm{gh}}\left(s_d,s_e\right).& \left[53\right]\end{array}$

[0144] Eleven other similar equations can be obtained from the constraint equations. From the marginalization of region (abde)→(de) we obtain: $\begin{array}{cc}{M}_{\mathrm{de}}^{\mathrm{ad}}\left(s_d,s_e\right)M_e^b\left(s_e\right)=\sum _{s_a{s}_b}{\varphi }_{\mathrm{ad}}\left(s_a,s_d\right){\varphi }_{\mathrm{ab}}\left(s_a,s_b\right){\varphi }_{\mathrm{de}}\left(s_d,s_e\right){\psi }_a\left(s_a\right){\psi }_b\left(s_b\right)M_{\mathrm{be}}^{\mathrm{cf}}\left(s_d,s_e\right)& \left[54\right]\\ \mathrm{From}\left(\mathrm{bcef}\right)->\left(\mathrm{be}\right):& \\ M_{\mathrm{be}}^{\mathrm{cf}}\left(s_b,s_e\right)M_e^f\left(s_e\right)=\sum _{s_c{s}_f}{\varphi }_{\mathrm{bc}}\left(s_b,s_c\right){\varphi }_{\mathrm{cf}}\left(s_c,s_f\right){\varphi }_{\mathrm{be}}\left(s_b,s_e\right){\psi }_c\left(s_c\right){\psi }_f\left(s_f\right)M_{\mathrm{fe}}^{\mathrm{ik}}\left(s_f,s_e\right)& \left[55\right]\\ \mathrm{From}\left(\mathrm{bcef}\right)->\left(\mathrm{fe}\right):& \\ M_{\mathrm{fe}}^{\mathrm{cb}}\left(s_f,s_e\right)M_e^b\left(s_e\right)=\sum _{s_c{s}_b}{\varphi }_{\mathrm{fc}}\left(s_f,s_c\right){\varphi }_{\mathrm{cb}}\left(s_c,s_b\right){\varphi }_{\mathrm{fe}}\left(s_f,s_e\right){\psi }_b\left(s_b\right){\psi }_f\left(s_f\right)M_{\mathrm{be}}^{\mathrm{ad}}\left(s_b,s_e\right)& \left[56\right]\\ \mathrm{From}\left(\mathrm{degh}\right)->\left(\mathrm{de}\right):& \\ M_{\mathrm{de}}^{\mathrm{gh}}\left(s_d,s_e\right)M_e^h\left(s_e\right)=\sum _{s_g{s}_h}{\varphi }_{dg}\left(s_d,s_g\right){\varphi }_{\mathrm{gh}}\left(s_g,s_h\right){\varphi }_{\mathrm{hc}}\left(s_h,s_c\right){\psi }_g\left(s_g\right){\psi }_h\left(s_h\right)M_{\mathrm{he}}^{\mathrm{if}}\left(s_h,s_c\right)& \left[57\right]\\ \mathrm{From}\left(\mathrm{degh}\right)->\left(\mathrm{he}\right):& \\ M_{\mathrm{he}}^{\mathrm{gd}}\left(s_h,s_c\right)M_c^d\left(s_e\right)=\sum _{s_g{s}_d}{\varphi }_{\mathrm{gh}}\left(s_g,s_h\right){\varphi }_{\mathrm{gd}}\left(s_g,s_d\right){\varphi }_{\mathrm{de}}\left(s_d,s_e\right){\psi }_d\left(s_d\right){\psi }_g\left(s_g\right)M_{dc}^{\mathrm{ab}}\left(s_d,s_c\right)& \left[58\right]\\ \mathrm{From}\left(\mathrm{efhi}\right)->\left(\mathrm{fe}\right):& \\ M_{\mathrm{fe}}^{\mathrm{ih}}\left(s_f,s_e\right)M_e^h\left(s_c\right)=\sum _{s_i{s}_h}{\varphi }_{\mathrm{fi}}\left(s_f,s_i\right){\varphi }_{\mathrm{hi}}\left(s_h,s_i\right){\varphi }_{\mathrm{fe}}\left(s_f,s_e\right){\psi }_h\left(s_h\right){\psi }_i\left(s_i\right)M_{\mathrm{he}}^{\mathrm{gd}}\left(s_h,s_e\right)& \left[59\right]\\ \mathrm{From}\left(\mathrm{efhi}\right)->\left(\mathrm{he}\right):& \\ M_{\mathrm{he}}^{\mathrm{fi}}\left(s_h,s_e\right)M_e^f\left(s_e\right)=\sum _{s_i{s}_j}{\varphi }_{\mathrm{hi}}\left(s_h,s_i\right){\varphi }_{\mathrm{fi}}\left(s_f,s_i\right){\varphi }_{\mathrm{fe}}\left(s_f,s_e\right){\psi }_f\left(s_f\right){\psi }_i\left(s_i\right)M_{\mathrm{fe}}^{\mathrm{cb}}\left(s_f,s_e\right)& \left[60\right]\\ \mathrm{From}\left(\mathrm{be}\right)->\left(e\right):& \\ M_c^b\left(e\right)=\sum _{s_b}{\varphi }_{\mathrm{be}}\left(s_b,s_e\right){\psi }_b\left(s_b\right)M_{\mathrm{be}}^{\mathrm{ad}}\left(s_b,s_e\right)M_{\mathrm{be}}^{\mathrm{cf}}\left(s_b,s_e\right)M_b^a\left(s_b\right)M_b^c\left(s_b\right)& \left[61\right]\\ \mathrm{From}\left(\mathrm{de}\right)->\left(e\right):& \\ M_e^d\left(e\right)=\sum _{s_d}{\varphi }_{\mathrm{de}}\left(s_d,s_e\right){\psi }_d\left(s_d\right)M_{\mathrm{de}}^{\mathrm{ab}}\left(s_d,s_e\right)M_{\mathrm{de}}^{\mathrm{gf}}\left(s_d,s_e\right)M_d^a\left(s_d\right)M_d^g\left(s_d\right)& \left[62\right]\\ \mathrm{From}\left(\mathrm{fe}\right)->\left(e\right):& \\ M_e^f\left(e\right)=\sum _{s_f}{\varphi }_{\mathrm{fe}}\left(s_f,s_e\right){\psi }_f\left(s_f\right)M_{\mathrm{fe}}^{\mathrm{cb}}\left(s_f,s_e\right)M_{\mathrm{fe}}^{\mathrm{ih}}\left(s_f,s_e\right)M_f^c\left(s_f\right)M_f^i\left(s_f\right)& \left[63\right]\\ \mathrm{From}\left(\mathrm{he}\right)->\left(e\right):& \\ M_c^h\left(e\right)=\sum _{s_h}{\varphi }_{\mathrm{he}}\left(s_h,s_e\right){\psi }_h\left(s_h\right)M_{\mathrm{he}}^{\mathrm{gd}}\left(s_h,s_e\right)M_{\mathrm{he}}^{\mathrm{if}}\left(s_h,s_e\right)M_h^g\left(s_h\right)M_h^i\left(s_h\right).& \left[64\right]\end{array}$
In step 609, we solve the self-consistent message equations 657 by iteratively updating the values according to the rules, beginning with some initial message values 658 and some “evidence” ψ659. The initial values 658 can be random or special values, e.g., all random positive numbers, or all ones. The “evidence” is given by the values of the Variables.