Zhaorong Zhang and Minyue Fu,

Abstract—Gaussian belief propagation algorithm (GaBP) is one of the most important distributed algorithms in signal processing and statistical learning involving Markov networks. It is well known that the algorithm correctly computes marginal density functions from a high dimensional joint density function over a Markov network in a finite number of iterations when the underlying Gaussian graph is acyclic. It is also known more recently that the algorithm produces correct marginal means asymptotically for cyclic Gaussian graphs under the condition of walk summability (or generalised diagonal dominance). This paper extends this convergence result further by showing that the convergence is exponential under the generalised diagonal dominance condition,and provides a simple bound for the convergence rate. Our results are derived by combining the known walk summability approach for asymptotic convergence analysis with the control systems approach for stability analysis.

I. Introduction

BELIEF propagation (BP) algorithm is a well-celebrated distributed algorithm for Markov networks that has been widely utilized in many disciplines, ranging from statistical learning and artificial intelligence to distributed estimation,distributed optimisation, networked control and digital communications [1]–[13].

Initially introduced by Pearl [1] in 1988, the BP algorithm is also known as Pearl’s algorithm, message-passing algorithm and sum-product algorithm. It is designed to compute the marginal probability densities of random variables from the joint probability density function over a large Markov network with sparse connections among individual random variables.The significance of the algorithm stems from the facts that it is fully distributed (i.e., only local information is needed for iteration computation) and that a wide range of application problems can be formulated as a BP problem. It is well known that the BP algorithm produces correct marginal probability densities in a finite number of iterations when the underlying graph for the joint density function is acyclic (i.e., no cycles or loops). But the properties of the algorithm for cyclic (loopy)graphs have been a major research topic over several decades.

The Gaussian BP algorithm (GaBP), a special version of the BP algorithm for Markov networks with Gaussian distributions (also known as Gaussian graphical model), has received special attention for the study of its convergence properties. In [2], it was shown that GaBP produces asymptotically the correct marginal means under the assumption that the joint information matrix is diagonal dominance. It was relaxed in [3] that the same asymptotic convergence holds when the joint information matrix is walksummable, which is equivalent to the condition of generalised diagonal dominance. The convergence property in [3] was generalised by [5] to allow an alternative decomposition of the optimizing function and more flexible message initialization.In [7], [8], necessary and sufficient conditions for asymptotic convergence of GaBP are studied. In [5], convergence properties of the BP algorithm for convex optimisation(including quadratic optimisation) are studied. A pertinent result in [5] is a bound on the convergence rate of the BP algorithm under a scaled diagonal dominance assumption and a particular decomposition of optimizing function. Many other interesting properties of GaBP can be found in, e.g., [14]–[16]and the references therein.

The purpose of this paper is to study the convergence rate of GaBP. Under the generalised diagonal dominance condition,we provide a simple bound for the exponential convergence rate of the marginal means. This bound is simply the spectral radius of the matrix related to the information matrix. This bound also coincides with that in [5] but under weaker conditions (see details later). Our results are derived by combining the walk summability approach in [3] for asymptotic convergence analysis with the control systems approach for stability analysis.

In the rest of the paper, we introduce GaBP in Section II and discuss the walk summability condition in Section III,followed by convergence rate analysis in Section IV,illustrating examples in Section V and conclusions in Section VI.

II. Problem Formulation

G={V,E}, where V={1,2,...,n} represents the set of nodes

Fig. 1. An example of Markov network.

A Gaussian graphical model is a Markov network with Gaussian distributions, characterised by an undirected graph and E is the set of edges (or unordered pairs {i,j}?V), with each nodei∈V being associated with a random variablexi.Fig. 1 shows an example of Markov network. The joint probability density forx=col{x1,x2,...,xn} is given by the following Gaussian density function

whereA={aij} is a sparse information matrix withai j=0 for allwhich is a symmetric and positive definite matrix, andbis the potential vector. It is straightforward to verify that the mean vector μ=E{x} and covariance matrixP=E{(x?μ)(x?μ)T}are given respectively by

The problem of concern is for each nodei∈V to compute,in a distributed fashion, the marginal density functionpi(xi) ofxi, defined by

wherex?idenotes the vectorxwith the componentxiremoved. It is well known that this amounts to computing the marginal mean μi(theith term of μ) and marginal variancepii(theith diagonal term ofP).

Using the Gaussian graphical model,p(x) can be factorized into

with

The BP algorithm is an iterative message-passing algorithm for computingpi(xi). In each iterationk, each nodei∈V computes and transmits to each nodej∈Ni(the set of neighbouring nodes ofi) the message

For a Gaussian graphical model, the messagecan be expressed as

This results in GaBP below

with

The initialization is done by takingand

The marginal mean and marginal variance ofare then given by, respectively

It is well known that, when the graph G is acyclic, GaBP converges inditerations with μi(k)=μiandpii(k)=piifor alli, wheredis the diameter of G (i.e., the largest distance between any two nodes in G). Actually, for each nodei,diiterations are sufficient to yield the above convergence, wherediis the largest distance from any node in G to nodei[2].(The distance of two nodes is the minimum path length between the nodes.)

For cyclic (or loopy) graphs, GaBP produces the correct marginal means asymptotically under certain conditions. In particular, it has been established in [3] that μi(k) converges touifor alliasymptotically under the so-called walk summability condition. This condition is also known to be equivalent to requiring the matrixAto be generalised diagonally dominant [17].

The goal of this paper is to study the convergence rate of GaBP under the same walk summability condition.

III. Walk Summability

Walk-sum analysis is an elegant approach introduced in [3](and their earlier references thereof) for studying the convergence of GaBP. Here we provide a quick summary of this approach.

Given a matrixR={ri j}∈Rn×nand its induced graph G=(V,E), a walkwin the graph is a node sequence

and its length isl. The weight of the walk is defined to be

As a convention, a single nodei∈V is regarded as a special(zero-length) walk with its weight ?(i)=1. A walkwfrom nodeitojis also denoted byw:i→j, and such a walk with lengthlis denoted byThe set of all walks from nodeito nodejis denoted by {i→j}, and the set of all lengthlwalks from nodeito nodejis denoted by. The walksum of a set of weightsWis denoted by

The importance of walk sums is revealed in the relationship that (i,j)th element of matrixRlis equal to

which can be verified by matrix multiplication. Now we give the definition of walk summability [3].

Definition 1:A matrixA∈Rn×nwithaii=1 for alliis said to be walk-summable if all the walk-sums ?({i→j}) converge absolutely, i.e.,| c onverges for alli,j. This is the same as the unordered sumis well defined (i.e.,converges to the same value for every possible summation order) for alli,j. Further, a linear systemAx=bis said to be walk-summable ifAis walk-summable.

DefiningR={ri j}=I?Aandthe following properties are known for walk-summable systems [3].

Lemma 1 [3]:The following conditions are equivalent.

1)A∈Rn×nwithaii=1 for alliis walk-summable;

4)

Using the walk-sum interpretation, the Gaussian variancePand mean μ in (2) can be expressed by walk sums under the assumption of walk summability [3]. More specifically, using Lemma 1,which implies that

and (10), we get

The connection between walk summability and diagonal dominance is revealed in the result below [17]. Recall [18]that a matrixA={aij} is called diagonally dominant ifaii>0 andf or alli.

Lemma 2 [3], [17]:A matrixA∈Rn×nwithaii=1 for alliis walk-summable (i.e.,<1) if and only ifAis generalised diagonally dominant, i.e., there exists a diagonal matrixD>0 such thatD?1ADis diagonally dominant.

IV. Convergence Rate Analysis

This section presents the main result of this paper on the convergence rate of GaBP. The key to this analysis is the socalled unwrapped tree graph proposed in [2], which is a computation tree graph, associated with the GaBP. Using this tool, the asymptotic convergence of GaBP was proved in [2]under the assumption of diagonal dominance. This tool was further used in [3] to relax the diagonal dominance assumption to walk summability (or equivalently, generalised diagonal dominance). Here we use the same tool for convergence rate analysis.

A. Unwrapped Tree Graph

Following the work of [2], we construct an unwrapped tree with deptht>0 for a loopy graph G [2]. Take nodeito be the root and then iterate the following procedurettimes:

1) Find all leaves of the tree (start with the root);

2) For each leaf, find all the nodes in the loopy graph that neighbor this leaf node, except its parent node in the tree, and add all these nodes as the children to this leaf node.

The variables and weights for each node in the unwrapped tree are copied from the corresponding nodes in the loopy graph. It is clear that taking each node as root node will generate a different unwrapped tree. Fig. 2 shows the unwrapped tree around root node 1 for a loopy graph. Note,for example, that nodes 1′,1′′,1′′′,1‘,1“,1“‘ all carry the same valuesb1anda11. Similarly, if node 1' is the parent (or child)of nodej′in the unwrapped tree, and node 1 and nodejare a wrapped version of nodes 1 andj, thena1′j′=a1j(or. A similar comment applies to unwrappedbi.

Fig. 2. Left: a loopy graph; Right: the unwrapped tree for root node 1 with 4 layers (t =4).

We have the following key property.

Lemma 3 [3]:There is a one-to-one correspondence between finite-length walks in G that end ati, and walks inThat is, every finite-length walk in G has a counterpart in somewith somei∈V and some sufficiently largek, and every finite-length walk infor anyi∈V andk≥0 corresponds to a finite-length walk in G.

B. Main Result

We first establish a relationship between μi(k) in (6)(obtained by GaBP) and the walks in

Lemma 4:Under the assumption that the information matrixAin (1) is walk summable, we have, for anyi∈V andk≥0

Proof:Without loss of generality, we assumei=1. For the unwrapped graphconsider the corresponding matrixand vectorDefineThen,can be solved by applying GaBP onAs noted in Section II, sinceis a tree graph, it is well known [2] that applying GaBP toresults in a correct solution for(the first component ofinkiterations because every node inis no more thankhops away from node 1. On the other hand, due to the fact that the parameters inandare all copied fromAandb,applying GaBP to the original graph G forkiterations is identical to applying it to G1(k). That is,in (6), which is obtained by applying GaBP on G forkiterations, is equal toNow, it is also known that every tree graph is walksummable [3]. Thus, we can apply (13) toto obtain

Now we can state the main result, derived by combining the walk summability approach in [3] (i.e., Lemmas 1, 2 and 4)with the control systems approach. That is, we view the evolution of μi(k)?μias a dynamic system and examine its exponential stability property.

Theorem 1:Suppose the information matrixAin (1) is generalised diagonally dominant. Then, the convergence rate of GaBP is at leasti.e.,

for alli∈V andk≥0 , whereCis a constant (independent ofk).

Proof:Firstly, using Lemma 2, we know that the assumption of general diagonal dominance is equivalent to assuming walk summability. In particular,Ais invertible.

From (2), we have μ=A?1b. Using the walk summability condition and Lemma 1, we get (13), i.e.,

On the other hand, μi(k) is given by (14), according to Lemma 4. Combining the above, we get

Denote byWi(k) the set of all the walks that end at nodeiwith walk length greater thank, and bythe subset of all the walks containing nodes not in Gk. It is clear that every walk inhas length greater thank. Then, the above expression can be rewritten as

Remark 1:The bound in Theorem 1 coincides with that in[5]. But our assumptions are weaker. More specifically, [5]

V. Illustrating Examples

Fig. 3. 13-node graph.

Fig. 4. 1000-node graph.

wherethe term l og10ρ2represents the slope.

Fig. 5. GaBP iterations for the 13-node graph.

Fig. 6. GaBP iterations for the 1000-node graph.

The simulation results for both two examples have shown that the error decreases exponentially with the increase of the iteration number. The slope for the 13-node example is measured to be roughly ?1.0502, corresponding to the convergence rate of 10?1.0502/2≈0.2985. The slope for the 1000-node example is measured to be roughly ?1.0642,corresponding to the convergence rate of 1 0?1.0642/2≈0.2937.In comparison, for the 13-node graph, the spectral radius ofis 0.9671. We see that in both examples, theupper is 0.6100; and for the 1000-node graph, the spectral radius of bounds the actual convergence rate of GaBP.

VI. Conclusions

In this paper, we have analysed the convergence property of GaBP for Markov networks and provided a simple bound on the convergence rate. The bound is characterised byand is guaranteed to be less than 1 under the walk summability (or generalised diagonal dominance) assumption. This result gives a simple extension to the known asymptotic convergence property of GaBP under the same assumption [2], [3]. We see in the simulation results that the actual convergence rate is faster than predicted byIt would be interesting to see how this bound can be further improved. Other future directions include relaxing the walk summability assumption (possibly by following the work of [15]), and generalising GaBP to wider distributed estimation and distributed optimisation problems.