Faster Algorithms for Max-Product Message-Passing

1 Introduction

It is well-known that exact inference in tree-structured graphical models can be accomplished efficiently by message-passing operations following a simple protocol making use of the distributive law (Aji and McEliece, 2000; Kschischang et al., 2001). It is also well-known that exact inference in arbitrary graphical models can be solved by the junction-tree algorithm; its efficiency is determined by the size of the maximal cliques after triangulation, a quantity related to the treewidth of the graph.

Figure 1 illustrates an attempt to apply the junction-tree algorithm to some graphical models containing cycles. If the graphs are not chordal ((a) and (b)), they need to be triangulated, or made chordal (red edges in (c) and (d)). Their clique-graphs are then guaranteed to be junction-trees, and the distributive law can be applied with the same protocol used for trees; see Aji and McEliece (2000) for a beautiful tutorial on exact inference in arbitrary graphs. Although the models in this example contain only pairwise factors, triangulation has increased the size of their maximal cliques, making exact inference substantially more expensive. Hence approximate solutions in the original graph (such as loopy belief-propagation, or inference in a loopy factor-graph) are often preferred over an exact solution via the junction-tree Algorithm.

Even when the model’s factors are the same size as its maximal cliques, neither exact nor approximate inference algorithms take advantage of the fact that many factors consist only of latent variables. In many models, those factors that are conditioned upon the observation contain fewer latent variables than the purely latent cliques. Examples are shown in Figure 2. This encompasses a wide variety of models, including grid-structured models for optical flow and stereo disparity as well as chain and tree-structured models for text or speech.

Figure 1: The models at left ((a) and (b)) can be triangulated ((c) and (d)) so that the junction-tree algorithm can be applied. Despite the fact that the new models have larger maximal cliques, the corresponding potentials are still factored over pairs of nodes only. Our algorithms exploit this fact.

Figure 2: Some graphical models to which our results apply: *cliques containing observations have fewer latent variables than purely latent cliques*. White nodes correspond to the observation, gray nodes to the labeling. In other words, cliques containing a white node encode the *data likelihood*, whereas cliques containing only gray nodes encode *priors*.


(a)	(b)	(c)	(d)

In this paper, we exploit the fact that the maximal cliques (after triangulation) often have potentials that factor over subcliques, as illustrated in Figure 1. We will show that whenever this is the case, the expected computational complexity of exact inference can be improved (both the asymptotic upper-bound and the actual runtime).

Additionally, we will show that this result can be applied so long as those cliques that are conditioned upon an observation contain fewer latent variables than those cliques consisting of purely latent variables; the ‘purely latent’ cliques can be pre-processed offline, allowing us to achieve the same benefits as described in the previous paragraph.

We show that these properties reveal themselves in a wide variety of real applications. Both of our improvements shall increase the class of problems for which inference via max-product belief-propagation is tractable.

A core operation encountered in the junction-tree algorithm is that of finding the index that chooses the largest product amongst two lists of length $N$ :

^i=argmaxi∈{1…N}{va[i]×vb[i]}.

(1)

Our results stem from the realization that while (eq. 1) appears to be a linear time operation, it can be decreased to $O (\sqrt{N})$ (in the expected case) if we know the permutations that sort $v_{a}$ and $v_{b}$ .

A preliminary version of this work appeared in McAuley and Caetano (2010).

1.1 Summary of Results

A selection of the results to be presented in the remainder of this paper can be summarized as follows:

We are able to lower the asymptotic expected running time of the max-product belief-propagation for any graphical model whose cliques factorize into lower-order terms.
The results obtained are exactly those that would be obtained by the traditional version of the algorithm, i.e., no approximations are used.
Our algorithm also applies whenever cliques containing an observed variable contain fewer latent variables than purely latent cliques, as in Figure 2 (meaning that certain computations can be taken offline).
For any cliques composed of pairwise factors, we obtain an expected speed-up of at least $Ω (\sqrt{N})$ (assuming $N$ states per node; $Ω$ denotes an asymptotic lower-bound).
For example, in models with third-order cliques containing pairwise terms, message-passing is reduced from $Θ (N^{3})$ to $O (N^{2} \sqrt{N})$ , as in Figure 1(d). For models containing pairwise (but purely latent) cliques, message-passing is reduced from $Θ (N^{2})$ to $O (N \sqrt{N})$ , as in Figure 2.
For cliques composed of $K$ -ary factors, the expected speed-up generalizes to at least $Ω (\frac{1}{K} N^{\frac{1}{K}})$ , though it is never asymptotically slower than the original solution.
The expected-case improvement is derived under the assumption that the order-statistics of different factors are independent.
If the different factors have ‘similar’ order statistics, the performance will be better than the expected case.
If the different factors have ‘opposite’ order statistics, the performance will be worse than the expected case, but is never asymptotically more expensive than the traditional version of the algorithm.

Our results do not apply for every semiring $S (+, \cdot)$ , but only to those whose ‘addition’ operation defines an order (for example, $min$ or $max$ ); we also assume that under this ordering, our ‘multiplication’ operator satisfies

a < b \land c < d \Rightarrow a \cdot c < b \cdot d .

(2)

Thus our results certainly apply to the max-sum and min-sum semirings (as well as max-product and min-product, assuming non-negative potentials), but not for sum-product (for example). Consequently, our approach is useful for computing MAP-states, but cannot be used to compute marginal distributions. We also assume that the domain of each node is discrete.

We shall initially present our algorithm as it applies to models of the type shown in Figure 1. The more general (and arguably more useful) application of our algorithm to those models in Figure 2 shall be deferred until Section 4, where it can be seen as a straightforward generalization of our initial results.

1.2 Related Work

There has been previous work on speeding-up message-passing algorithms by exploiting some type of structure in certain graphical models. For example, Kersting et al. (2009) study the case where different cliques share the same potential function. In Felzenszwalb and Huttenlocher (2006), fast message-passing algorithms are provided for cases in which the potential of a 2-clique is only dependent on the difference

of the latent variables (which is common in some computer vision applications); they also show how the algorithm can be made faster if the graphical model is a bipartite graph. In

Kumar and Torr (2006), the authors provide faster algorithms for the case in which the potentials are truncated, whereas in Petersen et al. (2008) the authors offer speed-ups for models that are specifically grid-like.

The latter work is perhaps the most similar in spirit to ours, as it exploits the fact that certain factors can be sorted in order to reduce the search space of a certain maximization problem. In practice, this leads to linear speed-ups over a $Θ (N^{4})$ algorithm.

Another closely related paper is that of Park and Darwiche (2003). This work can be seen to compliment ours in the sense that it exploits essentially the same type of factorization that we study, though it applies to sum-product versions of the algorithm, rather than the max-product version that we shall study. Kjærulff (1998) also exploits factorization within cliques of junction-trees, albeit a different type of factorization than that studied here.

In Section 3, we shall see that our algorithm is closely related to a well-studied problem known as ‘funny matrix multiplication’ (Kerr, 1970). The worst-case complexity of this problem has been studied in relation to the all-pairs shortest path problem (Alon et al., 1997; Karger et al., 1993).

2 Background

The notation we shall use is briefly defined in Table 1. We shall assume throughout that the max-product semiring is being used, though our analysis is almost identical for any suitable choice.

Example	description
$A; B$	capital letters refer to sets of nodes (or similarly, cliques);
$A \cup B; A \cap B; A ∖ B$	standard set operators are used ( $A ∖ B$ denotes set difference);
$dom (A)$	the domain of a set; this is just the Cartesian product of the domains of each element in the set;
$P$	bold capital letters refer to arrays;
$x$	bold lower-case letters refer to vectors;
$x [a]$	vectors are indexed using square brackets;
$P [n]$	similarly, square brackets are used to index a row of a 2-d array,
$P [n]$	or a row of an $(\| n \| + 1)$ -dimensional array;
$P^{X}; v^{a}$	superscripts are just labels, i.e., $P^{X}$ is an array, $v^{a}$ is a vector;
$v_{a}$	constant subscripts are also labels, i.e., if $a$ is a constant, then $v_{a}$ is a constant vector;
$x_{i}; x_{A}$	variable subscripts define variables; the subscript defines the domain of the variable;
$n \|_{X}$	if $n$ is a constant vector, then $n \|_{X}$ is the restriction of that vector to those indices corresponding to variables in $X$ (assuming that $X$ is an ordered set);
$Φ_{A}; Φ_{A} (x_{A})$	a function over the variables in a set $A$ ; the argument $x_{A}$ will be suppressed if clear, given that ‘functions’ are essentially arrays for our purposes;
$Φ_{i, j} (x_{i}, x_{j})$	a function over a pair of variables $(x_{i}, x_{j})$ ;
$Φ_{A} (n \|_{B}; x_{A ∖ B})$	if one argument to a function is constant (here $n \|_{B}$ ), then it becomes a function over fewer variables (in this case, only $x_{A ∖ B}$ is free);

Table 1: Notation

MAP-inference in a graphical model $G$ consists of solving an optimization problem of the form

^x = a r g m a x x \prod C \in C Φ_{C} (x_{C}),

(3)

where $C$ is the set of maximal cliques in $G$ . This problem is often solved via message-passing algorithms such as the junction-tree algorithm, loopy belief-propagation, or inference in a factor graph (Aji and McEliece, 2000; Weiss, 2000; Kschischang et al., 2001).

Two of the fundamental steps encountered in message-passing algorithms are defined below. Firstly, the message from a clique $X$ to an intersecting clique $Y$ is defined by

(4)

(where $Γ (X)$ returns the neighbors of the clique $X$ ). If such messages are computed after $Y$ has received messages from all of its neighbors except $X$ (i.e., $Γ (X) ∖ Y$ ), then this defines precisely the update scheme used by the junction-tree algorithm. The same update scheme is used for loopy belief-propagation, though it is done iteratively in a randomized fashion.

Secondly, after all messages have been passed, the MAP-states for a subset of nodes $M$ (assumed to belong to a clique $X$ ) is computed using

m_{M} (x_{M}) = max x_{X ∖ M} ⎧ ⎨ ⎩ Φ_{X} (x_{X}) \prod Z \in Γ (X) m_{Z \to X} (x_{X \cap Z}) ⎫ ⎬ ⎭ .

(5)

Often, the clique-potential $Φ_{X} (x_{X})$ shall be decomposable into several smaller factors, i.e.,

Φ_{X} (x_{X}) = \prod F \subset X Φ_{F} (x_{F}) .

(6)

Some simple motivating examples are shown in Figure 3

: a model for pose estimation from

Sigal and Black (2006), a ‘skip-chain CRF’ from Galley (2006), and a model for shape matching from Coughlan and Ferreira (2002). In each case, the triangulated model has third-order cliques, but the potentials are only pairwise. Other examples have already been shown in Figure 1; analogous cases are ubiquitous in many real applications.

Figure 3: (a) A model for pose reconstruction from Sigal and Black (2006); (b) A ‘skip-chain CRF’ from Galley (2006); (c) A model for deformable matching from Coughlan and Ferreira (2002). Although the (triangulated) models have cliques of size three, their potentials factorize into pairwise terms.


(a)	(b)	(c)

The optimizations we suggest shall apply to general problems of the form

m_{M} (x_{M}) = max x_{X ∖ M} \prod F \subset X Φ_{F} (x_{F}),

(7)

which subsumes both (eq. 4) and (eq. 5), where we simply treat the messages as factors of the model. Algorithm 1 gives the traditional solution to this problem, which does not exploit the factorization of $Φ_{X} (x_{X})$ . This algorithm runs in $Θ (N^{| X |})$ , where $N$ is the number of states per node, and $| X |$ is the size of the clique $X$ (we assume that for a given $x_{X}$ , computing $\prod_{F \subset X} Φ_{F} (x_{F})$ takes constant time, as our optimizations shall not modify this cost).

0: a clique

X

whose max-marginal

m_{M} (x_{M})

(where

M \subset X

) we wish to compute; assume that each node in

X

has domain

{1 \dots N}

1: for

m \in dom (M)

{i.e.,

{1 \dots N}^{| M |}

} do

m a x := - \infty

3: for

y \in dom (X ∖ M)

4: if

\prod_{F \subset X} Φ_{F} (m |_{F}; y |_{F}) > m a x

then

m a x := \prod_{F \subset X} Φ_{F} (m |_{F}; y |_{F})

6: end if

7: end for {this loop takes

Θ (N^{| X ∖ M |})

}

m_{M} (m) := m a x

9: end for {this loop takes

Θ (N^{| X |})

}

10: Return:

m_{M}

Algorithm 1 Brute-force computation of max-marginals

3 Optimizing Algorithm 1

In order to specify a more efficient version of Algorithm 1, we begin by considering the simplest possible nontrivial factorization: a clique of size three containing pairwise factors. In such a case, our aim is to compute

m_{i, j} (x_{i}, x_{j}) = max x_{k} Φ_{i, j, k} (x_{i}, x_{j}, x_{k}),

(8)

which we have assumed takes the form

m_{i, j} (x_{i}, x_{j}) = max x_{k} Φ_{i, j} (x_{i}, x_{j}) \times Φ_{i, k} (x_{i}, x_{k}) \times Φ_{j, k} (x_{j}, x_{k}) .

(9)

For a particular value of $(x_{i}, x_{j}) = (a, b)$ , we must solve

m_{i, j} (a, b) = Φ_{i, j} (a, b) \times max x_{k} Φ_{i, k} (a, x_{k})      v_{a} \times Φ_{j, k} (b, x_{k})      v_{b},

(10)

which we note is in precisely the form shown in (eq. 1).

There is a close resemblance between (eq. 10) and the problem of multiplying two matrices: if the ‘ $max$ ’ in (eq. 10) is replaced by summation, we essentially recover traditional matrix multiplication. While traditional matrix multiplication is well known to have a sub-cubic worst-case solution (see Strassen, 1969), the version in (eq. 10) (often referred to as ‘funny matrix multiplication’, or simply ‘max-product matrix multiplication’) is known to be cubic in the worst case, assuming that only multiplication and comparison operations are used (Kerr, 1970). The complexity of solving (eq. 10) can also be shown to be equivalent to the all-pairs shortest path problem, which is studied in Alon et al. (1997). Not surprisingly, we shall not improve the worst-case complexity, but shall instead give far better expected-case performance than existing solutions. Just as Strassen’s algorithm can be used to solve (eq. 10) when maximization is replaced by summation, there has been work studying the problem of sum-product inference in graphical models, subject to the same type of factorization we discuss (Park and Darwiche, 2003).

As we have previously suggested, it will be possible to solve (eq. 10) efficiently if $v_{a}$ and $v_{b}$ are already sorted. We note that $v_{a}$ will be reused for every value of $x_{j}$ , and likewise $v_{b}$ will be reused for every value of $x_{i}$ . Sorting every row of $Φ_{i, k}$ and $Φ_{j, k}$ can be done in $Θ (N^{2} log N)$ (for $2 N$ rows of length $N$ ).

The following elementary lemma is the key observation required in order to solve (eq. 10) efficiently:

Lemma 1.

If the $p^{th}$ largest element of $v_{a}$ has the same index as the $q^{th}$ largest element of $v_{b}$ , then we only need to search through the $p$ largest values of $v_{a}$ , and the $q$ largest values of $v_{b}$ ; any corresponding pair of smaller values could not possibly be the largest solution.

This observation is used to construct Algorithm 2. Here we iterate through the indices starting from the largest values of $v_{a}$ and $v_{b}$ , stopping once both indices are ‘behind’ the maximum value found so far (which we then know is the maximum). This algorithm is demonstrated pictorially in Figure 4.

0: two vectors

v_{a}

and

v_{b}

, and permutation functions

p_{a}

and

p_{b}

that sort them in decreasing order (so that

v_{a} [p_{a} [1]]

is the largest element in

v_{a}

)

1: Initialize:

s t a r t := 1

{e n d}_{a} := p_{a}^{- 1} [p_{b} [1]]

{e n d}_{b} := p_{b}^{- 1} [p_{a} [1]]

{if

e n d_{b} = k

, then the largest element in

v_{a}

has the same index as the

k^{th}

largest element in

v_{b}

}

b e s t := p_{a} [1]

m a x := v_{a} [b e s t] \times v_{b} [b e s t]

3: if

v_{a} [p_{b} [1]] \times v_{b} [p_{b} [1]] > m a x

then

b e s t := p_{b} [1]

m a x := v_{a} [b e s t] \times v_{b} [b e s t]

5: end if

6: while

s t a r t < {e n d}_{a}

{in practice, we could also stop if

s t a r t < {e n d}_{b}

, but the version given here is the one used for analysis in Appendix A} do

s t a r t := s t a r t + 1

8: if

v_{a} [p_{a} [s t a r t]] \times v_{b} [p_{a} [s t a r t]] > m a x

then

b e s t := p_{a} [s t a r t]

10:

m a x := v_{a} [b e s t] \times v_{b} [b e s t]

11: end if

12: if

p_{b}^{- 1} [p_{a} [s t a r t]] < {e n d}_{b}

then

13:

{e n d}_{b} := p_{b}^{- 1} [p_{a} [s t a r t]]

14: end if

15: {repeat Lines 8–14, interchanging

a

and

b

}

16: end while {this takes expected time

O (\sqrt{N})

}

17: Return:

b e s t

Algorithm 2 Find

i

such that

v_{a} [i] \times v_{b} [i]

is maximized

$Step 1:{\parbox{333.8874pt}{\includegraphics[angle=-90, wid% th=260.172pt]{figures_fast/alg2_1}}$

$Step 2:{\parbox{333.8874pt}{\includegraphics[width=260.172% pt]{figures_fast/alg2_2}}$

$Step 3:{\parbox{333.8874pt}{\includegraphics[width=260.172% pt]{figures_fast/alg2_3}}$

$Step 4:{\parbox{333.8874pt}{\includegraphics[width=260.172% pt]{figures_fast/alg2_4}}$

$Step 5:{\parbox{333.8874pt}{\includegraphics[width=260.172% pt]{figures_fast/alg2_5}}$

Figure 4: Algorithm 2, explained pictorially. The arrows begin at

p_{a} [s t a r t]

and

p_{b} [s t a r t]

; the red line connects

{e n d}_{a}

and

{e n d}_{b}

, behind which we need not search; a dashed arrow is used when a new maximum is found. Note that in the event that

v_{a}

and

v_{b}

contain repeated elements, they can be sorted arbitrarily.

A prescription of how Algorithm 2 can be used to solve (eq. 8) is given in Algorithm 3. Determining precisely the running time of Algorithm 2 (and therefore Algorithm 3) is not trivial, and will be explored in depth in Appendix A. We note that if the expected-case running time of Algorithm 2 is $O (f (N))$ , then the time taken to solve Algorithm 3 shall be $O (N^{2} (log N + f (N)))$ . At this stage we shall state an upper-bound on the true complexity in the following theorem:

Theorem 1.

The expected running time of Algorithm 2 is $O (\sqrt{N})$ , yielding a speed-up of at least $Ω (\sqrt{N})$ in cliques containing pairwise factors.

0: a potential

Φ_{i, j, k} (a, b, c) = Φ_{i, j} (a, b) \times Φ_{i, k} (a, c) \times Φ_{j, k} (b, c)

whose max-marginal

m_{i, j} (x_{i}, x_{j})

we wish to compute

1: for

n \in {1 \dots N}

2: compute

P^{i} [n]

by sorting

Φ_{i, k} (n, x_{k})

{takes

Θ (N log N)

}

3: compute

P^{j} [n]

by sorting

Φ_{j, k} (n, x_{k})

{

P^{i}

and

P^{j}

are

N \times N

arrays, each row of which is a permutation;

Φ_{i, k} (n, x_{k})

and

Φ_{j, k} (n, x_{k})

are functions over

x_{k}

, since

n

is constant in this expression}

4: end for {this loop takes

Θ (N^{2} log N)

}

5: for

(a, b) \in {1 \dots N}^{2}

(v_{a}, v_{b}) := (Φ_{i, k} (a, x_{k}), Φ_{j, k} (b, x_{k}))

(p_{a}, p_{b}) := (P^{i} [a], P^{j} [b])

b e s t := A l g o r i t h m ??? (v_{a}, v_{b}, p_{a}, p_{b})

{takes

O (\sqrt{N})

}

m_{i, j} (a, b) := Φ_{i, j} (a, b) \times Φ_{i, k} (a, b e s t) \times Φ_{j, k} (b, b e s t)

10: end for {this loop takes

O (N^{2} \sqrt{N})

}{the total running time is

O (N^{2} log N + N^{2} \sqrt{N})

, which is dominated by

O (N^{2} \sqrt{N})

}

11: Return:

m_{i, j}

Algorithm 3 Use Algorithm 2 to compute the max-marginal of a 3-clique containing pairwise factors

3.1 An Extension to Higher-Order Cliques with Three Factors

The simplest extension that we can make to Algorithms 2 and 3 is to note that they can be applied even when there are several overlapping terms in the factors. For instance, Algorithm 3 can be adapted to solve

m_{i, j} (x_{i}, x_{j}) = max x_{k}, x_{m} Φ_{i, j} (x_{i}, x_{j}) \times Φ_{i, k, m} (x_{i}, x_{k}, x_{m}) \times Φ_{j, k, m} (x_{j}, x_{k}, x_{m}),

(11)

and similar variants containing three factors. Here both $x_{k}$ and $x_{m}$ are shared by $Φ_{i, k, m}$ and $Φ_{j, k, m}$ . We can follow precisely the reasoning of the previous section, except that when we sort $Φ_{i, k, m}$ (similarly $Φ_{j, k, m}$ ) for a fixed value of $x_{i}$ , we are now sorting an array rather than a vector (Algorithm 3, Lines 2 and 3); in this case, the permutation functions $p_{a}$ and $p_{b}$ in Algorithm 2 simply return pairs of indices. This is illustrated in Figure 5. Effectively, in this example we are sorting the variable $x_{k, m} := x_{k} \otimes x_{m}$ , which has state space of size $N^{2}$ .

$Step 1:{\parbox{333.8874pt}{\includegraphics[width=260.172% pt]{figures_fast/alg2_examp}}$

Figure 5: The reasoning applied in Algorithm 2 applies even when the elements of

p_{a}

and

p_{b}

are multidimensional indices.

As the number of shared terms increases, so does the improvement to the running time. While (eq. 11) would take $Θ (N^{4})$ to solve using Algorithm 1, it takes only $O (N^{3})$ to solve using Algorithm 3 (more precisely, if Algorithm 2 takes $O (f (N))$ , then (eq. 11) takes $O (N^{2} f (N^{2}))$ , which we have mentioned is $O (N^{2} \sqrt{N^{2}}) = O (N^{3})$ ). In general, if we have $S$ shared terms, then the running time is $O (N^{2} \sqrt{N^{S}})$ , yielding a speed-up of $Ω (\sqrt{N^{S}})$ over the naïve solution of Algorithm 1.

3.2 An Extension to Higher-Order Cliques with Decompositions Into Three Groups

By similar reasoning, we can apply our algorithm to cases where there are more than three factors, in which the factors can be separated into three groups. For example, consider the clique in Figure 6(a), which we shall call $G$ (the entire graph is a clique, but for clarity we only draw an edge when the corresponding nodes belong to a common factor). Each of the factors in this graph have been labeled using either differently colored edges (for factors of size larger than two) or dotted edges (for factors of size two), and the max-marginal we wish to compute has been labeled using colored nodes. We assume that it is possible to split this graph into three groups such that every factor is contained within a single group, along with the max-marginal we wish to compute (Figure 6, (b)). If such a decomposition is not possible, we will have to resort to further extensions to be described in Section 3.3.

Figure 6: Algorithm 4, explained pictorially. In this case, the most computationally intensive step is the marginalization of $Z$ (in step (c)), which takes $Θ (N^{5})$ . However, the algorithm can actually be applied *recursively* to the group $Z$ , resulting in an overall running time of $O (N^{4} \sqrt{N})$ , for a max-marginal that would have taken $Θ (N^{8})$ to compute using the naïve solution of Algorithm 1.


(a)	(b)


(c)	(d)	(e)


(f)	(g)

0: potentials

Φ_{G} (x) = Φ_{X} (x_{X}) \times Φ_{Y} (x_{Y}) \times Φ_{Z} (x_{Z})

; each of the factors should be contained in exactly one of these terms, and we assume that

M \subseteq X

(see Figure 6)

1: Define:

X^{'} := ((Y \cup Z) \cap X) \cup M

;

Y^{'} := (X \cup Z) \cap Y

;

Z^{'} := (X \cup Y) \cap Z

{

X^{'}

contains the variables in

X

that are shared by at least one other group; alternately, the variables in

X ∖ X^{'}

appear only in

X

(sim. for

Y^{'}

and

Z^{'}

)}

2: compute

Ψ^{X} (x_{X^{'}}) := {max}_{X ∖ X^{'}} Φ_{X} (x_{X})

{we are marginalizing over those variables in

X

that do not appear in any of the other groups (or in

M

); this takes

Θ (N^{S})

if done by brute force (Algorithm 1), but may also be done by a recursive call to Algorithm 4}

3: compute

Ψ^{Y} (x_{Y^{'}}) := {max}_{Y ∖ Y^{'}} Φ_{Y} (x_{Y})

4: compute

Ψ^{Z} (x_{Z^{'}}) := {max}_{Z ∖ Z^{'}} Φ_{Z} (x_{Z})

5: for

n \in dom (X \cap Y)

6: compute

P^{Y} [n]

by sorting

Ψ^{Y} (n; x_{Y^{'} ∖ X})

{takes

Θ (S_{∖} N^{S_{∖}} log N)

;

Ψ^{Y} (n; x_{Y^{'} ∖ X})

is free over

x_{Y^{'} ∖ X}

, and is treated as an array by ‘flattening’ it;

P^{Y} [n]

contains the

| Y^{'} ∖ X | = | (Y \cap Z) ∖ X |

-dimensional indices that sort it}

7: end for {this loop takes

Θ (S_{∖} N^{S} log N)

}

8: for

n \in dom (X \cap Z)

9: compute

P^{Z} [n]

by sorting

Ψ^{Z} (n; x_{Z^{'} ∖ X})

10: end for {this loop takes

Θ (S_{∖} N^{S} log N)

}

11: for

n \in dom (X^{'})

12:

(v_{a}, v_{b}) := (Ψ^{Y} (n |_{Y^{'}}; x_{Y^{'} ∖ X^{'}}), Ψ^{Z} (n |_{Z^{'}}; x_{Z^{'} ∖ X^{'}}))

{

n |_{Y^{'}}

is the ‘restriction’ of the vector

n

to those indices in

Y^{'}

(meaning that

n |_{Y^{'}} \in dom (X^{'} \cap Y^{'})

); hence

Ψ^{Y} (n |_{Y^{'}}; x_{Y^{'} ∖ X^{'}})

is free in

x_{Y^{'} ∖ X^{'}}

, while

n |_{Y^{'}}

is fixed}

13:

(p_{a}, p_{b}) := (P^{Y} [n |_{Y^{'}}], P^{Z} [n |_{Z^{'}}])

14:

b e s t := A l g o r i t h m ??? (v_{a}, v_{b}, p_{a}, p_{b})

{takes

O (\sqrt{S_{∖}})

}

15:

m_{X} (n) := Ψ^{X} (n) \times Ψ^{Y} (b e s t; n |_{Y^{'}}) \times Ψ^{Z} (b e s t; n |_{Z^{'}})

16: end for

17:

m_{M} (x_{M}) := A l g o r i t h m ??? (m_{X}, M)

{i.e., we are using Algorithm 1 to marginalize

m_{X} (x_{X})

with respect to

M

; this takes

Θ (N^{S})

}

Algorithm 4 Compute the max-marginal of

G

with respect to

M

, where

G

is split into three groups

Ideally, we would like these groups to have size $≃ | G | / 3$ , though in the worst case they will have size no larger than $| G | - 1$ . We call these groups $X$ , $Y$ , $Z$ , where $X$ is the group containing the max-marginal $M$ that we wish to compute. In order to simplify the analysis of this algorithm, we shall express the running time in terms of the size of the largest group, $S = max (| X |, | Y |, | Z |)$ , and the largest difference, $S_{∖} = max (| Y ∖ X |, | Z ∖ X |)$ . The max-marginal can be computed using Algorithm 4.

The running times shown in Algorithm 4 are loose upper-bounds, given for the sake of expressing the running time in simple terms. More precise running times are given in Table 2; any of the terms shown in Table 2 may be dominant. Some example graphs, and their resulting running times are shown in Figure 7.

Description	lines	time
Marginalization of $Φ_{X}$ , without recursion	2	$Θ (N^{\| X \|})$
Marginalization of $Φ_{Y}$	3	$Θ (N^{\| Y \|})$
Marginalization of $Φ_{Z}$	4	$Θ (N^{\| Z \|})$
Sorting $Φ_{Y}$	5–7	$Θ (\| Y^{'} ∖ X \| N^{\| Y^{'} \|} log N)$
Sorting $Φ_{Z}$	8–10	$Θ (\| Z^{'} ∖ X \| N^{\| Z^{'} \|} log N)$
Running Algorithm 2 on the sorted values	11–16	$O (N^{\| X^{'} \|} \sqrt{N^{\| (Y^{'} \cap Z^{'}) ∖ X^{'} \|}})$

Table 2: Detailed running time analysis of Algorithm 4; any of these terms may be asymptotically dominant

Graph:					${$ A complete graph $K_{M}$ , with pairwise terms $}$
	(a)	(b)	(c)	(d)	(e)
Algorithm 1:	$Θ (N^{5})$	$Θ (N^{3})$	$Θ (N^{11})$	$Θ (N^{6})$	$Θ (N^{M})$
Algorithm 4:	$O (N^{3} \sqrt{N})$	$O (N^{2} \sqrt{N})$	$O (N^{6} \sqrt{N})$	$O (N^{5})$	$O (N^{5 M / 6})$
Speed-up:	$Ω (N \sqrt{N})$	$Ω (\sqrt{N})$	$Ω (N^{4} \sqrt{N})$	$Ω (N)$	$Ω (N^{M / 6})$

Figure 7: Some example graphs whose max-marginals are to be computed with respect to the colored nodes, using the three regions shown. Factors are indicated using differently colored edges, while dotted edges always indicate pairwise factors. (a) is the region

Z

from Figure 6 (recursion is applied again to achieve this result); (b) is the graph used to motivate Algorithm 3; (c) shows a query in a graph with regular structure; (d) shows a complete graph with six nodes; (e) generalizes this to a clique with

M

nodes.

3.2.1 Applying Algorithm 4 Recursively

The marginalization steps of Algorithm 4 (Lines 2, 3, and 4) may further decompose into smaller groups, in which case Algorithm 4 can be applied recursively. For instance, the graph in Figure 7(a) represents the marginalization step that is to be performed in Figure 6(c) (Algorithm 4, Line 4). Since this marginalization step is the asymptotically dominant step in the algorithm, applying Algorithm 4 recursively lowers the asymptotic complexity.

Another straightforward example of applying recursion in Algorithm 4 is shown in Figure 8, in which a ring-structured model is marginalized with respect to two of its nodes. Doing so takes $O (M N^{2} \sqrt{N})$ ; in contrast, solving the same problem using the junction-tree algorithm (by triangulating the graph) would take $Θ (M N^{3})$ . Loopy belief-propagation takes $Θ (M N^{2})$ per iteration, meaning that our algorithm will be faster if the number of iterations is $Ω (\sqrt{N})$ . Naturally, Algorithm 3 could be applied directly to the triangulated graph, which would again take $O (M N^{2} \sqrt{N})$ .

Figure 8: In the above example, lines 2–4 of Algorithm 4 are applied recursively, achieving a total running time of $O (M N^{2} \sqrt{N})$ for a loop with $M$ nodes (our algorithm achieves the same running time in the triangulated graph).

	$O (N^{2})$
	$+$
	$O (2 N^{2} \sqrt{N})$
	$+$
	$O (4 N^{2} \sqrt{N})$ (by Algorithm 3)

3.3 A General Extension to Higher-Order Cliques

Naturally, there are cases for which a decomposition into three terms is not possible, such as

m_{i, j, k} (x_{i}, x_{j}, x_{k}) = max x_{m} Φ_{i, j, k} (x_{i}, x_{j}, x_{k}) \times Φ_{i, j, m} (x_{i}, x_{j}, x_{m}) \times Φ_{i, k, m} (x_{i}, x_{k}, x_{m}) \times Φ_{j, k, m} (x_{j}, x_{k}, x_{m})

(12)

(i.e., a clique of size four with third-order factors). However, if the model contains factors of size $K$ , it must always be possible to split it into $K + 1$ groups (e.g. four in the case of (eq. 12)).

Our optimizations can easily be applied in these cases simply by adapting Algorithm 2 to solve problems of the form

^i=argmaxi∈{1…N}{v1[i]×v2[i]×⋯×vK[i]}.

(13)

Pseudocode for this extension is presented in Algorithm 5. Note carefully the use of the variable $r e a d$ : we are storing which indices have been read to avoid re-reading them; this guarantees that our Algorithm is never asymptotically worse than the naïve solution. Figure 9 demonstrates how such an algorithm behaves in practice. Again, we shall discuss the running time of this extension in Appendix A

. For the moment, we state the following theorem:

Theorem 2.

Algorithm 5 generalizes Algorithm 2 to $K$ lists with an expected running time of $O (K N^{\frac{K - 1}{K}})$ , yielding a speed-up of at least $Ω (\frac{1}{K} N^{\frac{1}{K}})$ in cliques containing $K$ -ary factors. It is never worse than the naïve solution, meaning that it takes $O (min (N, K N^{\frac{K - 1}{K}}))$ .

K

vectors

v_{1} \dots v_{K}

; permutation functions

p_{1} \dots p_{K}

that sort them in decreasing order; a vector

r e a d

indicating which indices have been read, and a unique value

T \notin r e a d

{

r e a d

is essentially a boolean array indicating which indices have been read; since creating this array is an

O (N)

operation, we create it externally, and reuse it

O (N)

times; setting

r e a d [i] = T

indicates that a particular index has been read; we use a different value of

T

for each call to this function so that

r e a d

can be reused without having to be reinitialized}

1: Initialize:

s t a r t := 1

m a x := {max}_{p \in {p_{1} \dots p_{K}}} \prod_{k = 1}^{K} v_{k} [p [1]]

b e s t := {a r g m a x}_{p \in {p_{1} \dots p_{K}}} \prod_{k = 1}^{K} v_{k} [p [1]]

2: for

k \in 1 \dots K

{e n d}_{k} := {max}_{q \in {p_{1} \dots p_{K}}} p_{k}^{- 1} [q [1]]

4: end for

r e a d [p [1]] = T

6: while

s t a r t < max {{e n d}_{1} \dots {e n d}_{K}}

s t a r t := s t a r t + 1

8: if

r e a d [p [s t a r t]] := T

then

9: continue

10: end if

11:

r e a d [p [s t a r t]] := T

12:

m := {max}_{p \in {p_{1} \dots p_{K}}} \prod_{k = 1}^{K} v_{k} [p [s t a r t]]

13:

b := {a r g m a x}_{p \in {p_{1} \dots p_{K}}} \prod_{k = 1}^{K} v_{k} [p [s t a r t]]

14: if

m > m a x

then

15:

b e s t := b

16:

m a x := m

17: end if

18: for

k \in {1 \dots K}

19:

e_{k} := {max}_{q \in {p_{1} \dots p_{K}}} p_{k}^{- 1} [q [s t a r t]]

20: end for

21: for

k \in {1 \dots K}

22:

{e n d}_{k} := min (e_{k}, {e n d}_{k})

23: end for

24: end while {see Appendix A for running times}

25: Return:

b e s t

Algorithm 5 Find

i

such that

\prod_{k = 1}^{K} v_{k} [i]

is maximized

$Step 1:{\parbox{333.8874pt}{\includegraphics[angle=-90, wid% th=260.172pt]{figures_fast/alg3}}$

Figure 9: Algorithm 2 can easily be extended to cases including more than two sequences.

Using Algorithm 5, we can similarly extend Algorithm 4 to allow for any number of groups (pseudocode is not shown; all statements about the groups $Y$ and $Z$ simply become statements about $K$ groups ${G_{1} \dots G_{K}}$ , and calls to Algorithm 2 become calls to Algorithm 5). The one remaining case that has not been considered is when the sequences $v_{1} \dots v_{K}$ are functions of different (but overlapping) variables; naïvely, we can create a new variable whose domain is the product space of all of the overlapping terms, and still achieve the performance improvement guaranteed by Theorem 2; in some cases, better results can again be obtained by applying recursion, as in Figure 7.

As a final comment we note that we have not provided an algorithm for choosing how to split the variables of a model into $(K + 1)$ -groups. We note even if we split the groups in a naïve way, we are guaranteed to get at least the performance improvement guaranteed by Theorem 2, though more ‘intelligent’ splits may further improve the performance. Furthermore, in all of the applications we have studied, $K$ is sufficiently small that it is inexpensive to consider all possible splits by brute-force.

4 Exploiting ‘Data Independence’ in Latent Factors

While (eq. 3) gave the general form of MAP-inference in a graphical model, it will often be more convenient to express our objective function as being conditioned upon some observation, $y$ . Thus inference consists of solving an optimization problem of the form

^x (y) = a r g m a x x \prod C \in C Φ_{C} (x_{C} | y_{C}) .

(14)

When our objective function is written in this way, further factorization is often possible, yielding an expression of the form

^x (y) = a r g m a x x \prod F \in F Φ_{F} (x_{F} | y_{F})      data dependent \times \prod C \in C Φ_{C} (x_{C})      data independent,

(15)

where each $F \in F$ is a subset of some $C \in C$ . We shall say that those factors that do not depend on the observation are ‘data independent’.

By far the most common instance of this type of model has ‘data dependent’ factors consisting of a single latent variable, and conditioned upon a single observation, and ‘data independent’ factors consisting of a pair of latent variables. This was precisely the class of models depicted at the beginning of our paper in Figure 2, whose objective function takes the form

^x (y) = a r g m a x x \prod i \in N Φ_{i} (x_{i} | y_{i})      node % potential \times \prod (i, j) \in E Φ_{i, j} (x_{i}, x_{j})      edge potential

(16)

(where $N$ and $E$ are the set of nodes and edges in our graphical model). As in the Section 3, we shall concern ourselves with this version of the model, and explain only briefly how it can be applied with larger factors, as in Section 3.2.

Note that in (eq. 16) we are no longer concerned solely with exact inference via the junction-tree algorithm. In many models, such as grids and rings, (eq. 16) shall be solved approximately by means of either loopy belief-propagation, or inference in a factor graph.

Given the decomposition of (eq. 16), message-passing now takes the form

m_{A \to B} (q) = Φ_{i} (q) \times max y_{j} Φ_{j} (y_{j})      v_{a} \times Φ_{i, j} (q, y_{j})      v_{b}

(17)

(where $A = (i, j)$ and $B = (i, k)$ ). Just as we made the comparison between (eq. 10) and matrix multiplication, we can see (eq. 17) as being related to the multiplication of a matrix ( $Φ_{i, j}$ ) with a vector ( $Φ_{j}$ ), again with summation replaced by maximization. Given the results we have already shown, it is trivial to solve (eq. 17) in $O (N \sqrt{N})$ if we know the permutations that sort $Φ_{j}$ , and the rows of $Φ_{i, j}$ . The algorithm for doing so is shown in Algorithm 6. The difficultly we face in this instance is that sorting the rows of $Φ_{i, j}$ takes $Θ (N^{2} log N)$ , i.e., longer than Algorithm 6 itself.

This problem is circumvented due to the following simple observation: since $Φ_{i, j} (x_{i}, x_{j})$ consists only of latent variables (and not upon the observation), this ‘sorting’ step can take place offline, i.e., before the ‘data’ has been observed.

Two further observations mean that even this offline cost can often be avoided. Firstly, many models have a ‘homogeneous’ prior, i.e., the same prior is shared amongst every edge (or clique) of the model. In such cases, only a single ‘copy’ of the prior needs to be sorted, meaning that in any model containing $Ω (log N)$ edges, speed improvements can be gained over the naive implementation. Secondly, where an iterative algorithm (such as loopy belief-propagation) is to be used, the sorting step need only take place prior to the first iteration; if $Ω (log N)$ iterations of belief propagation are to be performed (or indeed, if the number of edges multiplied by the number of iterations is $Ω (log N)$ ), we shall again gain speed improvements even when the sorting step is done online.

In fact, the second of these conditions obviates the need for data independence altogether. In other words, in any pairwise model in which $Ω (log N)$ iterations of belief propagation are to be performed, the pairwise terms need to be sorted only during the first iteration. Thus these improvements apply to those models in Figure 1, so long as the number of iterations is $Ω (log N)$ .

0: a potential

Φ_{i, j} (a, b) \times Φ_{i} (a | y_{i}) \times Φ_{j} (b | y_{j})

whose max-marginal

m_{i} (x_{i})

we wish to compute, and a set of permutation functions

P

such that

P [i]

sorts the

i^{th}

row of

Φ_{i, j}

(in decreasing order).

1: compute the permutation function

p_{a}

by sorting

Ψ_{j}

{takes

Θ (N log N)

}

2: for

q \in {1 \dots N}

(v_{a}, v_{b}) := (Ψ_{j}, Φ_{i, j} (q, x_{j} | y_{i}, y_{j}))

b e s t := A l g o r i t h m ??? (v_{a}, v_{b}, p_{a}, P [q])

{

O (\sqrt{N})

}

m_{A \to B} (q) := Φ_{i} (q) \times Φ_{j} (b e s t) \times Φ_{i, j} (q, b e s t | y_{i}, y_{j})

6: end for {expected-case

O (N \sqrt{N})

}

7: Return:

m_{A \to B}

Algorithm 6 Solve (eq. 17) using Algorithm 2

4.1 Extension to Higher-Order Cliques

Just as in Section 3.2, we can extend Algorithm 6 to factors of any size, so long as the purely latent cliques contain more latent variables than those cliques that depend upon the observation. The analysis for this type of model is almost exactly the same as that presented in Section 3.2, except that any terms consisting of purely latent variables are processed offline.

As we mentioned in 3.2, if a model contains (non-maximal) factors of size $K$ , we will gain a speed-up of $Ω (\frac{1}{K} N^{\frac{1}{K}})$ . If in addition there is a factor (either maximal or non-maximal) consisting of purely latent variables, we can still obtain a speed-up of $Ω (\frac{1}{K + 1} N^{\frac{1}{K + 1}})$ , since this factor merely contributes an additional term to (eq. 13). Thus when our ‘data-dependent’ terms contain only a single latent variable (i.e., $K = 1$ ), we gain a speed-up of $Ω (\sqrt{N})$ , as in Algorithm 6.

5 Performance Improvements in Existing Applications

Our results are immediately compatible with several applications that rely on inference in graphical models. As we have mentioned, our results apply to any model whose cliques decompose into lower-order terms.

Often, potentials are defined only on nodes and edges of a model. A $D^{th}$

-order Markov model has a tree-width of

D

, despite often containing only pairwise relationships. Similarly ‘skip-chain CRFs’ (Sutton and McCallum, 2006; Galley, 2006), and junction-trees used in SLAM applications (Paskin, 2003) often contain only pairwise terms, and may have low tree width under reasonable conditions. In each case, if the tree-width is

D

, Algorithm 4 takes

O (M N^{D} \sqrt{N})

(for a model with

M

nodes and

N

states per node), yielding a speed-up of

Ω (\sqrt{N})

Models for shape matching and pose reconstruction often exhibit similar properties (Tresadern et al., 2009; Donner et al., 2007; Sigal and Black, 2006). In each case, third-order cliques factorize into second order terms; hence we can apply Algorithm 3 to achieve a speed-up of $Ω (\sqrt{N})$ .

Another similar model for shape matching is that of Felzenszwalb (2005); this model again contains third-order cliques, though it includes a ‘geometric’ term constraining all three variables. Here, the third-order term is independent of the input data, meaning that each of its rows can be sorted offline, as described in Section 4. In this case, those factors that depend upon the observation are pairwise, meaning that we achieve a speed-up of $Ω (N^{\frac{1}{3}})$ . Further applications of this type shall be explored in Section 6.2.

In Coughlan and Ferreira (2002), deformable shape-matching is solved approximately using loopy belief-propagation. Their model has only second-order cliques, meaning that inference takes $Θ (M N^{2})$ per iteration. Although we cannot improve upon this result, we note that we can typically do exact inference in a single iteration in $O (M N^{2} \sqrt{N})$ ; thus our model has the same running time as $O (\sqrt{N})$ iterations of the original version. This result applies to all second-order models containing a single loop (Weiss, 2000).

In McAuley et al. (2008), a model is presented for graph-matching using loopy belief-propagation; the maximal cliques for $D$ -dimensional matching have size $(D + 1)$ , meaning that inference takes $Θ (M N^{D + 1})$ per iteration (it is shown to converge to the correct solution); we improve this to $O (M N^{D} \sqrt{N})$ .

Interval graphs can be used to model resource allocation problems (Fulkerson and Gross, 1965); each node encodes a request, and overlapping requests form edges. Maximal cliques grow with the number of overlapping requests, though the constraints are only pairwise, meaning that we again achieve an $Ω (\sqrt{N})$ improvement.

Belief-propagation can be used to solve LP-relaxations in pairwise graphical models. In Sontag et al. (2008), LP-relaxations are computed for pairwise models by constructing several third-order ‘clusters’, which compute pairwise messages for each of their edges. Again, an $Ω (\sqrt{N})$ improvement is achieved.

Finally, in Section 6.2 we shall explore a variety of applications in which we have pairwise models of the form shown in (eq. 16). In all of these cases, we see an (expected) reduction of a $Θ (M N^{2})$ message-passing algorithm to $O (M N \sqrt{N})$ .

Table 3 summarizes these results. Reported running times reflect the expected case. Note that we are assuming that max-product belief-propagation is being used in a discrete model; some of the referenced articles may use different variants of the algorithm (e.g. Gaussian models, or approximate inference schemes). We believe that our improvements may revive the exact, discrete version as a tractable option in these cases.

Reference	description	running time	our method
McAuley et al. (2008)	$D$ -d graph-matching	$Θ (M N^{D + 1})$ (iter.)	$O (M N^{D} \sqrt{N})$ (iter.)
Sutton and McCallum (2006)	Width- $D$ skip-chain	$O (M N^{D})$	$O (M N^{D - 1} \sqrt{N})$
Galley (2006)	Width-3 skip-chain	$Θ (M N^{3})$	$O (M N^{2} \sqrt{N})$
Paskin (2003) (discrete case)	SLAM, width $D$	$O (M N^{D})$	$O (M N^{D - 1} \sqrt{N})$
Tresadern et al. (2009)	Deformable matching	$Θ (M N^{3})$	$O (M N^{2} \sqrt{N})$
Coughlan and Ferreira (2002)	Deformable matching	$Θ (M N^{2})$ (iter.)	$O (M N^{2} \sqrt{N})$
Sigal and Black (2006)	Pose reconstruction	$Θ (M N^{3})$	$O (M$

Faster Algorithms for Max-Product Message-Passing

Authors

Convergent message passing algorithms - a unifying view

A new look at reweighted message passing

Anytime Exact Belief Propagation

Feedback Message Passing for Inference in Gaussian Graphical Models

A Comparison of Lauritzen-Spiegelhalter, Hugin, and Shenoy-Shafer Architectures for Computing Marginals of Probability Distributions

Generalized sequential tree-reweighted message passing

Variational message passing (VMP) applied to LDA

1 Introduction

1.1 Summary of Results

1.2 Related Work

2 Background

3 Optimizing Algorithm 1

Lemma 1.

Theorem 1.

3.1 An Extension to Higher-Order Cliques with Three Factors

3.2 An Extension to Higher-Order Cliques with Decompositions Into Three Groups

3.2.1 Applying Algorithm 4 Recursively

3.3 A General Extension to Higher-Order Cliques

Theorem 2.

4 Exploiting ‘Data Independence’ in Latent Factors

4.1 Extension to Higher-Order Cliques

5 Performance Improvements in Existing Applications

Faster Algorithms for Max-Product Message-Passing

Authors

Related Research

Convergent message passing algorithms - a unifying view

A new look at reweighted message passing

Anytime Exact Belief Propagation

Feedback Message Passing for Inference in Gaussian Graphical Models

A Comparison of Lauritzen-Spiegelhalter, Hugin, and Shenoy-Shafer Architectures for Computing Marginals of Probability Distributions

Generalized sequential tree-reweighted message passing

Variational message passing (VMP) applied to LDA

This week in AI

1 Introduction

1.1 Summary of Results

1.2 Related Work

2 Background

3 Optimizing Algorithm 1

Lemma 1.

Theorem 1.

3.1 An Extension to Higher-Order Cliques with Three Factors

3.2 An Extension to Higher-Order Cliques with Decompositions Into Three Groups

3.2.1 Applying Algorithm 4 Recursively

3.3 A General Extension to Higher-Order Cliques

Theorem 2.

4 Exploiting ‘Data Independence’ in Latent Factors

4.1 Extension to Higher-Order Cliques

5 Performance Improvements in Existing Applications