The cavity approach for Steiner trees Packing Problems
Alfredo Braunstein∗
DISAT, Politecnico di Torino, Corso Duca Degli Abruzzi 24, Torino, Italy
Human Genetics Foundation, Via Nizza 52, Torino, Italy and
Collegio Carlo Alberto, Via Real Collegio 1, Moncalieri, Italy

Anna Paola Muntoni†

arXiv:1712.07041v2 [] 11 Jan 2018

DISAT, Politecnico di Torino, Corso Duca Degli Abruzzi 24, Torino, Italy
The Belief Propagation approximation, or cavity method, has been recently applied to several
combinatorial optimization problems in its zero-temperature implementation, the Max-Sum algorithm. In particular, recent developments to solve the Edge-Disjoint paths problem and the Prize
collecting Steiner tree Problem on graphs have shown remarkable results for several classes of graphs
and for benchmark instances. Here we propose a generalization of these techniques for two variants
of the Steiner trees packing problem where multiple “interacting” trees have to be sought within a
given graph. Depending on the interaction among trees we distinguish the Vertex-Disjoint Steiner
trees Problem, where trees cannot share nodes, from the Edge-Disjoint Steiner trees Problem, where
edges cannot be shared by trees but nodes can be members of multiple trees. Several practical problems of huge interest in network design can be mapped into these two variants, for instance, the
physical design of Very Large Scale Integration (VLSI) chips.
The formalism described here relies on two components edge-variables that allows us to formulate
a massage-passing algorithm for the V-DStP and two algorithms for the E-DStP differing in the
scaling of the computational time with respect to some relevant parameters. We will show that one of
the two formalisms used for the edge-disjoint variant allow us to map the Max-Sum update equations
into a weighted maximum matching problem over proper bipartite graphs. We developed a heuristic
procedure based on the Max-Sum equations that shows excellent performance in synthetic networks
(in particular outperforming standard multi-step greedy procedures by large margins) and on large
benchmark instances of VLSI for which the optimal solution is known, on which the algorithm found
the optimum in two cases and the gap to optimality was never larger than 4%.

CONTENTS
I. Introduction
II. Two Steiner packing problems

3
4

III. An arborescent representation
1. Branching model
2. The flat formalism
A. Constraints for the Vertex-Disjoint Steiner trees problem
B. Constraints for the Edge-Disjoint Steiner trees problem

4
4
5
5
5

IV. Boltzmann distribution and marginals

6

V. The cavity equations
A. Vertex-disjoint Steiner trees Problem
B. Edge-disjoint Steiner trees problem
1. Neighbors occupation formalism
2. Mapping into a matching problem
C. The parameter D

9
9
10
10
12
14

VI. Max-Sum for loopy graphs
A. Reinforcement
B. Max-Sum based heuristics
1. Shortest Path Tree
2. Minimum Spanning Tree

14
15
15
15
16

VII. Numerical results
A. Fully connected graphs
1. Uncorrelated edge weights
2. Correlated edge weights

16
17
17
17

2
3. Fixed fraction of terminals
B. Regular graphs
1. Energy as a function of the degree
2. Energy as a function of the number of communications
C. Grid graphs
1. Branching and flat models for the V-DStP e E-DStP (neighbors occupation formalism)
2. V-DStP for VLSI circuits
VIII. Summary of results

∗
†

18
18
19
19
19
20
20
22

Acknowledgments

23

References

23

A. Message-Passing equations for V-DStP

25

B. Recursive expression of Z q for the E-DStP

26

alfredo.braunstein@polito.it
anna.muntoni@polito.it

3
I.

INTRODUCTION

The minimum Steiner tree problem (MStP) is an important combinatorial problem that consists in finding a
connected sub-graph within a given weighted graph, able to span a subset of vertices (called terminal s) with minimum
cost. It is easy to see that if weights are strictly positive the sub-graph satisfying all these constraints must be a tree.
The decisional problem of determining whether a solution within a given cost bound exists is NP-complete (it was
one of Karp’s original 21 problems [13]). The large difficulty of the MStP can be seen to arise from the large space of
subsets of non-terminal vertices (Steiner nodes). There exist several variants and generalizations of the MStP. One of
the most studied is the Prize collecting Steiner problem (PCStP) that have many applications in network technologies,
such as optimal heating and optical fibers distribution [14], in biology, e.g. in finding signal pathways in cell [4]. In
the Prize-Collecting variant the notion of terminals is relaxed so that every vertex has an associated prize (or reward).
The prize of included nodes is counted negatively in the solution cost (so that profitable vertices with positive reward
lower the total cost). In this variant the cost of the optimal tree will be the best trade-off between prizes of included
nodes and the cost of their connections given by edge-weights.
In this work we will address the Packing of Steiner Trees problem; where we aim at finding, within the same graph,
multiple Steiner trees which span disjoints sets of terminals and its prize-collecting version. We consider two different
variants regarding the interaction between trees. In the Vertex-disjoint Steiner trees problem (V-DStP), different trees
cannot share vertices (and consequently they cannot share edges either); in the Edge-disjoint Steiner tree problem
(E-DStP) only edge sets are pairwise disjoint but nodes can be shared by different trees. Being generalizations of
PCStP, both problems are NP-hard; but from a practical point of view the packing problems are more difficult than
their single-tree counterpart as it can be seen from the fact that even finding feasible solutions, i.e. trees satisfying
the interaction constraints (regardless the cost), is NP-hard [12]. In addition to its mathematical interest, a lot of
attention is devoted to the practical solution of Packing of Steiner trees problems since several layout design issues
arising from Very Large Scale Integration (VLSI) circuits [9, 10, 15] can be mapped into these variants of the MStP
[11, 12]. Integrated systems are composed by a huge number of logical units, called cells, typically arranged in 2D or
3D grids. Some specific cells, forming the so-called net, must be connected to one another in order to satisfy some
working conditions. The physical design phase of these circuits addresses the problem of connecting each element of
a net minimizing some objective function, namely the power consumption induced by the wires of the connection. It
can be easily seen that connecting the cells of a net at minimum power consumption is equivalent to solve a MStP on
a 2D or 3D grid graph. Thus, the problem of concurrently connecting multiple and disjoint sets of nets can be easily
mapped into a V-DStP or an E-DStP. The most common approaches to these combinatorial optimization problems
rely on linear programming formulations, for instance, the multi-commodity flow model [12].
In this work we devise six different models to represent these two problems; one for V-DStP and two for E-DStP,
the first one more suitable for graphs where the density of terminals is low and the second for instances with low
graph connectivity. We attempt their solution through the Cavity Method of statistical physics; and its algorithmic
counterpart, the Belief Propagation (BP) iterative algorithm (or rather, its zero-temperature limit, the Max-Sum
(MS) algorithm [16, 17]). This technique is an approximation scheme first developed to study disordered systems, and
nowadays applied to a wide range of optimization problems. Once a proper set of variables is defined, the optimization
problem, i.e. the constrained minimization of a cost function, can be mapped into the problem of finding the groundstate of a generalized system with local interactions. Ground-states can be investigated through observables related
to the Boltzmann-Gibbs distribution at zero temperature but, in most of the interesting cases, its exact computation
involves impractical computations. MS consists in iterating closed massage-passing equations on a factor graph that is
closely related to the original graph that, at convergence, provide an estimate of the marginal probability distribution
of the variables of interest. The cavity method can be proven to be exact on tree graphs or in some models on random
networks in the asymptotic limit) but nevertheless in practice reaches notable performances on arbitrary graphs. It
should be noted that in a simplified version of the problem (the minimum spanning tree), fixed points of Max-Sum
can be proven to parametrize the optimal solution [6]. As usual, the iterative solution of the Max-Sum equations
involve the solution of a related problem in a local star-shaped subgraph; which for some of these models is not trivial
(i.e. its naive solution is exponentially slow in the degree). We devise a mapping of the problem into a minimum
matching problem that can be solved in polynomial time in the degree (leading e.g. to linear time per iteration on
Erdos-Renyi random graphs).
In combination with these three initial models, a variant called the flat formalism, borrowed from [8] can be
independently included, leading to six different model combinations for the two problems. The flat formalism is more
suitable for graphs with large diameter and/or few terminals as it allows to reduce considerably the solution space.
Interestingly, the resulting flat models can be seen as generalizations of both [8] and [1, 3]; as the edge-disjoint path
problems in the last two publications can be seen as a packing of steiner trees problem in which each tree has exactly
two terminals.
With these algorithmic tools on hand, we perform numerical simulations of complete, Erdos-Renyi and random

4
regular graphs and on benchmark instances of V-DStP arising from the VLSI design problem.

II.

TWO STEINER PACKING PROBLEMS

Given a graph G = (V, E) whose vertices have non-negative real prizes {cµi : i ∈ V, µ = 1, . . . , M } and whose edges
µ
have real positive weights wij
: (i, j) ∈ E, µ = 1, . . . , M , we consider the problem of finding M connected sub-graphs
Gµ = (Vµ , Eµ ) spanning disjoint sets of terminals {Tµ ⊆ Vµ , µ = 1, . . . , M } that minimize the following cost or energy
function


X
X
X µ
µ

(1)
H=
wij
ci +
µ

i∈V \Vµ

(i,j)∈Eµ

This definition of the cost is extremely general: node prizes and edge costs can depend on sub-graph µ. For directed
µ
µ
by considering oriented trees (the trees we will consider will be ultimately rooted and
6= wji
graphs, we can admit wij
thus oriented). In the following we refer to vertices with strictly positive prizes as (generalized) terminals in analogy
with the MStP. This particular case can be integrated in our formalism imposing cµi = +∞ if node i is a (true)
terminal of tree µ (it suffices to have a large enough value for cµi instead of +∞), and cµj = 0, ∀µ for any non-terminal
node j ∈ V . Since we interpret the solution-trees as networks that allow terminals to “communicate” we will refer to
each sub-graph Gµ as a “communication” µ flowing within the graph.
Subsets Gµ must satisfy some interaction constraints depending on the packing variant we are considering. In the
Vertex-disjoint Steiner trees Problem (V-DStP), vertex-sets Vµ must be pairwise disjoint, i.e. Vµ ∩ Vν = ∅ if µ 6= ν
and, consequently, also edge sets will be pairwise disjoint. In the Edge-disjoint Steiner trees Problem (E-DStP), only
edge sets must be pairwise disjoints, i.e. Eµ ∩ Eν = ∅ if µ 6= ν, but vertex sets can overlap.
III.

AN ARBORESCENT REPRESENTATION

To deal with these two combinatorial optimization problems we will define a proper set of interacting variables
defined on a factor graph which is closely related to the original graph G. The factor graph is the bipartite graph of
factors and variables, in which an edge between a factor and a variable exists if the factor depends on the variable.
More precisely, to each vertex i ∈ V we associate a factor node with factor ψi and to each edge (i, j) ∈ E we associate
a two components variable (dij , µij ) ∈ {−D, . . . , 0, . . . , D} × {0, . . . , M }. Our choice of the edge-variables is similar
to the one adopted in [8] but here, in addition to a “depth” component, we introduce a “communication” variable µij
by which we label edges forming different trees.
.
Compatibility functions ψi are defined in a way that allowed configurations of variables (d, µ) = {(dij , µij ) : (i, j) ∈ E}
are in one to one correspondence to feasible solutions of the Vertex-disjoint or Edge-disjoint variant of the Steiner
trees problem. In particular, in order to ensure Steiner sub-graphs to be trees, i.e. to be connected and acyclic,
we impose local constraints on variables di = {dij : j ∈ ∂i} and µi = {µij : j ∈ ∂i} through compatibility functions
ψi (di , µi ) that will be equal to one if the constraints are satisfied or zero otherwise.
Consider a solution to the V-DStP or the E-DStP. Each variable µij takes value from the set {0, 1, . . . , M } and
denotes to which sub graph, if any, does the edge (i, j) belongs; the state µij = 0 will conventionally mean that no
tree employs the edge (i, j). Components dij ∈ {−D, . . . , 0, . . . , D} have a meaning of “depth” or “distance” within
the sub-graph. Value dij = 0 conventionally means that such edge is not employed by any communication and thus
it is admitted if and only if the associated µij = 0.
Being the interactions among nodes different as we deal with the V-DStP or the E-DStP, we will define two different
compatibility functions, ψiV and ψiE , for the two problems. Both functions will be written with the help of a single-tree
compatibility function ψiµ for two different formulations of the constraints, the branching and the flat model.
1.

Branching model

Let us consider a sub-graph Gµ constituting part of the solution for the V-DStP or the E-DStP. For each node
i ∈ Vµ , the variable dij measures the length, in “steps”, of the unique path from node i to root rµ passing through
j ∈ ∂i. Variable dij will be strictly positive (negative) if j is one step closer (farther) then i to root rµ . Thus, every
edge will satisfy µij = µji and the anti-symmetric condition dij = −dji . A directed tree structure is guaranteed if
there exists only one neighbor j ∈ ∂i such that µij 6= 0 and dij > 0 and all remaining neighbors k ∈ ∂i\j can either

5
not enter in this solution or be member of tree Gµ at the distance dki = dij + 1 from root rµ . Mathematically, we can
define the following compatibility function

ψiµ, b (di , µi ) =

Y

δµji ,0 δdji ,0 +

j∈∂i

XX


Y

δµ,µji δdji ,−d

d>0 j∈∂i

(δµ,µki δdki ,d+1 + δµki ,0 δdki ,0 )

(2)

k∈∂i\j

where δx,y is the discrete Kronecker delta function equal to 1 if x = y and 0 otherwise.

2.

The flat formalism

The diameter of solutions representable using the branching model strongly depends on the value of the parameter D.
A small value of the depth parameter can certainly prevent the representation of more elongated and, possibly, more
energetically favored solutions but, at the same time, a big value of D will significantly slow down the computation
of the compatibility function in (2). The flat model relaxes the depth-increasing constraint in the sense that, under
certain conditions, it allows chains of nodes, within the solution, with equal depth. According to a flat representation,
the depth variable increases of one unity, if a node i is a terminal node or there exist two or more neighbors connected
to i within a sub-graph Gµ , i.e. the degree of node i within the sub-graph is more then two. It can be shown [8] that
for D = T , where T = |Tµ | is the number of terminals per communication, these constraints admits all feasible trees
plus some extra structures which contain disconnected cycles with no terminals and thus are energetically disfavored.
The indicator function of these constraints is
ψiµ, f (di , µi ) = δcµi ,0

XX

δµ,µki δ−d,dki

d>0 k∈∂i

X

Y

δµ,µli δdli ,d

l∈∂i\k

δµmi ,0 δdmi ,0

(3)

m∈∂i\{k,l}




Finally, the single-tree compatibility function ψiµ = 1 − 1 − ψiµ, b 1 − ψiµ, f of configuration satisfying exactly
one of the two constraints can be written as
ψiµ (di , µi ) = ψiµ, b (di , µi ) + ψiµ, f (di , µi )

(4)

since ψiµ, b ψiµ, f ≡ 0.

A.

Constraints for the Vertex-Disjoint Steiner trees problem

In the V-DStP a node can belong to none or at most one sub-graph Gµ ; thus, if a vertex i is member of a Steiner
tree, it neighbor edges can either participate to the same communication or be unused. Being nodes sets of the
solution non-overlapping, we can consider separately the ownership of a vertex to a particular tree. The compatibility
function ψiV can be then expressed as the sum over all possible trees of a single-tree compatibility function ψiµ in (4).
ψiV (di , µi ) =

M
X

ψiµ (di , µi )

(5)

µ=1

B.

Constraints for the Edge-Disjoint Steiner trees problem

Differently to the V-DStP, in the E-DStP a vertex can belong to an arbitrary number of communications (including
zero) with the constraint that the local tree structure must be concurrently satisfied for each communication. If
the node does not participate in the solution we must admit configurations in which di = 0 if µi = 0. For the
remaining cases, if some neighbors k ∈ ∂i is a members of a Steiner tree µ, its distances dki will be different from
zero if µki = µ, and, additionally, they will satisfy the topological constraints. We can mathematically express such
conditions through the compatibility functions

6

ψiE,b (di , µi ) =

M
Y




Y


µ=1

δdki δµki ,µ ,0 +

k∈∂i

XX

δdki δµki ,µ ,−d

d>0 k∈∂i

Y 



δdli δµli ,µ ,d+1 + δdli δµli ,µ ,0 

(6)

l∈∂i\k

for the branching model and

ψiE,f (di , µi ) =

Y
µ



δcµ ,0
i

XX

δ−d,dki δµki ,µ

d>0 k∈∂i

X

Y

δµ,µli δdli δµli ,µ ,d

l∈∂i\k

δdmi δµmi ,µ ,0 

(7)

m∈∂i\{k,l}

for the flat model. Notice that we can express
ψiE (di , µi ) = ψiE,b (di , µi ) + ψiE,f (di , µi )

(8)

or eventually, if we define d˜ki = dki δµki ,µ , we can rephrase it as a product over single-tree compatibility functions
(compare with 5)
ψiE (di , µi ) =

M
Y



ψiµ d̃i , µi

(9)

µ=1

Some examples of feasible assignments of variables for both branching and flat models are shown in Figure 1 on
page 7. On the top-left (Figure 1 on page 7 (a)) we see one instance of the V-DStP containing two sub-graphs, the
“red” having root “4” and the “green” rooted at node “3”, and on the right two “red” and “green” edge-disjoint Steiner
trees, rooted at “0” and “5” respectively. Roots are represented as square nodes in contrast to circle colored nodes
that are terminals. Edges employed in both solutions are figured as arrows whose labels denote the value of the
(positive) depth component while the color mirrors the communication component. In agreement with our branching
representation we see that depth components increase as we cover the solution from the root to the leaves for both
problems. In the top-left figure, nodes of the vertex-disjoint trees are members either of the “red” or the “greed” trees
but such constraint is relaxed on the top-right of Figure 1 on page 7 (b) for E-DStP. To underline it, we allow node
“9” to be a terminal of communication “green” and a Steiner node of the “red” tree as two incident edges, (0, 9) and
(9, 8), belong to the “red” solution.
To show an example of representation through the flat formalism, we picture in Figure 2 on page 8 (a) and (b) two
representations, one for each model, of the same solution to the V-DStP on a grid graph. According to the branching
formalism, we see in Figure 2 on page 8 (a) that we need a minimum depth of D = 6 to allow all terminals of the
”blue” communication to reach the root node “1” but, since the tree in actually a chain of nodes, we notice that in
Figure 2 on page 8 (b) the same solution can be represented in the flat formalism using D = 3. In fact, only each
time we reach a terminal node the depth variable increases of one unit. Depth variables must increase in another
condition, precisely when we reach a branching point: this is exactly what happens in the neighborhood of node “21”
of the orange solution, Figure 2 on page 8 (b).

IV.

BOLTZMANN DISTRIBUTION AND MARGINALS

The formalism introduced above allows us to map each solution of the packing of Steiner trees to a certain assignment
of variables d = {dij : (i, j) ∈ E} and µ = {µij : (i, j) ∈ E} of the associated factor graph. The cost function in (1)
can be then expressed in terms of the new variables as

H (d, µ) =



M
X
X
X µ


c
I
[µ
=
6
µ]
+
wij 
i
i



µ=1

i∈V

(10)

dij >0:
µij =µ

µ
where, for sake of simplicity, we consider the “homogeneous” case wij
= wij ∀µ and ψi can be either equal to ψiV or
E
ψi . The function I [· ] is the indicator function that is equal to one if its argument is true and zero otherwise. The
expression µi 6= µ means that none of the neighbors k ∈ ∂i satisfies µki = µ.
The Boltzmann-Gibbs distribution associated with the energy H (d, µ) is given by:

7
(a)

(b)

6

11
6
5

2
1

0

7

0

9

1

4

1

2
5

2

3

7

8

2

8

1

3

10

2

2

4

3

2

1

2

1
5

3

3

1

4

Figure 1. Figures (a) and (b) show a feasible assignment of the variables for the V-DStP (left) and the E-DStP (right) using
the branching model.

Q

i

P (d, µ) =

ψi (di , µi ) e−βH(d,µ)
Zβ

(11)

in which β is a positive parameter, called the “inverse temperature” as in the statistical mechanics framework, and
the normalization constant
Zβ =

XY

ψi (di , µi ) e−βH(d,µ)

i

d,µ

is the partition function. Configurations of the variables that do not satisfy the topological constraints will have
zero probability measure, whereas all other configuration will be weighted according to the sum of weights of used
edges and of the penalties of non-employed nodes. In the limit β → +∞ the distribution will be concentrated in
the configuration(s) that minimizes H (d, µ) (the ground state of the system) that are exactly the solutions of the
optimization problems. Thus we are interesting in determining, for each edge (i, j) ∈ E, the assignment of variables
that, in the β → +∞, limit maximizes the marginal probability distribution Pij defined as:

 X
Pij d˜ij , µ̃ij =
P (d, µ) δdij ,d˜ij δµij ,µ̃ij

(12)

d,µ

Unfortunately the computation of (12) is impractical as it would require the calculation of a sum of an exponential
number of terms. We seek to estimate these marginals via the cavity method approach. We report here a standard
formulation of the cavity equations and we refer the interested reader to [16] for the detailed derivation. At finite β
the BP equations on our factor graph are:


m (d , µ )

 ij ij ij

=





= e−βwki I[dki >0] mki (dki , µki )

nki (dki , µki )

1
Zij

P

{dki ,µki }:
k∈∂i\j

ψi (di , µi ) e−β

P

µ

cµ
i I[µi 6=µ]

k∈∂i\j

nki (dki , µki )
(13)

where
Zij =

Q

X
{dij ,µij }

mij (dij , µij )

8
1

1

2

3

4

5

6

11

12

16

17

18

22

23

24

4

5

6

11

12

16

17

18

22

23

24

2

7

8

3

4

9

10

5

6

13

14

15

1

19

3

20

2

2

21

(a)
1

1

2

3

2

7

8

2

3

9

10

3

3

13

14

15

1

(b)

19

2

20

2

21

2

Figure 2. Figures (a) and (b) picture the same solution to a V-DStP using the branching formalism (figure (a)) and the flat
representation (figure (b))

is the normalization constant or “partial” partition function. The functions mij are called cavity marginals or “messages”, suggesting that some information is flowing on edge (i, j) within the factor graph from node i to node j. In
fact, the values of the messages mij are in some sense proportional to the probability of a particular assignment
(dij , µij ) for edge (i, j) if the node j were temporarily erased from the graph.
The system of equations in (13) can be seen as fixed point equations that can be solved iteratively. Starting from
a set of initial cavity marginals at time t = 0, we iterate the right-hand-side of (13) until numerical convergence to a
fixed point is reached. At convergence we calculate an approximation to marginals in (12) via the cavity fields defined
as
Mij (dij , µij ) ∝ nij (dij , µij ) nji (−dij , µij )

(14)

where the proportional sign denotes that a normalization constant is missing.
Cavity equations for optimization problems can be easily obtained by substituting the mij and Mij with the
variables hij (dij , µij ) = limβ→∞ β1 log n (dij , µij ) and Hij (dij , µij ) = limβ→∞ β1 log Mij (dij , µij ) into (13) and (14)
that play the role of cavity marginals and fields in the zero-temperature limit; the resulting closed set of equations
is known as the Max-Sum algorithm. At convergence we can extract our optimal assignment of variables by the
computation of the decisional variables

d∗ij , µ∗ij ∈ arg max Hij (dij , µij )

(15)

(dij ,µij )

Hij (dij , µij ) = hij (dij , µij ) + hji (−dij , µij ) − C 0

(16)

where C 0 is an additive constant that guarantees that normalization condition in the zero-temperature limit, i.e.
max(dij ,µij ) Hij (dij , µij ) = 0, is satisfied. In practice, converge is reached when the decisional variables computed

9
as in (15) do not change after a predefined number of successive iterations (often 10 ÷ 30). Notice that taking the
β → +∞ limit of the message-passing equations at finite β is not equivalent to the zero-temperature limit of the
Boltzmann distribution in (11).
In the following section we will show how to derive equations for the cavity marginals and cavity fields, for finite β
and in the limit β → +∞, depending on we are dealing with the V-DStP or the E-DStP problem.

V.
A.

THE CAVITY EQUATIONS

Vertex-disjoint Steiner trees Problem

To derive the Belief Propagation equations for the V-DStP problem suffices to impose ψi (di , µi ) = ψiV (di , µi ) in
(13). By a change of variables, we will determine a Max-Sum algorithm for this variant.
Equations for messages can be easily obtained by using the properties of Kronecker delta functions in ψiV (di , µi );
the explicit derivation in reported in appendix A. We can differentiate three cases depending on we are updating
messages nij for positive, negative or null depth dij :


b
f

mij (d, µ) = m+ (d, µ) + m (d, µ)
b
mij (d, µ) = m− P
(d, µ) + mf (d, µ)

f
m (0, 0) = e−β µ cµi Q
b
ij
k∈∂i\j nki (0, 0) + m0 + m0

∀d > 0, µ 6= 0
∀d < 0, µ 6= 0

(17)

where mb+ (d, µ) , mb− (d, µ) , mf (d, µ) , and mb0 , mf0 are defined as

mb+ (d, µ) =

Y

[nki (d + 1, µ) + nki (0, 0)]

k∈∂i\j

X

mf (d, µ) = δcµi ,0

k∈∂i\j

mb− (d, µ) =

X

nki (d + 1, µ)

XX

nli (0, 0)

l∈∂i\{j,k}

k∈∂i\j

mb0 =

Y

nki (d, µ)

Y

[nli (d, µ) + nli (0, 0)]

l∈∂i\{j,k}

mb− (d, µ)

µ6=0 d<0

mf0

=

XX X
µ6=0 d<0 k∈∂i\j

nki (d, µ)

X
l∈∂i\{j,k}

Y

nli (−d, µ)

nmi (0, 0)

m∈∂i\{k,l,j}

Replacing hij (dij , µij ) = limβ→+∞ nij (dij , µij ) in (17) we obtain the Max-Sum equations:


n
o
f
b

h
(d,
µ)
=
max
h
(d,
µ)
,
h
(d,
µ)

ij
+
+


n
o
hij (d, µ) = max hb− (d, µ) , hf− (d, µ)

o
n


f
b
hij (0, 0) = max − P cµ + P
µ i
k∈∂i\j hki (0, 0) , h0 , h0
for

∀d > 0, µ 6= 0
∀d < 0, µ 6= 0

(18)

10

X

hb+ (d, µ) = −wij +

max {hki (d + 1, µ) , hki (0, 0)}

k∈∂i\j

hf+ (d, µ) = −wij + log δcµi ,0 + max




k∈∂i\j 

hki (d, µ) +

X
l∈∂i\j



hli (0, 0)





X

hb− (d, µ) = max hki (d + 1, µ) − wik +
k∈∂i\j

max {hli (d, µ) , hli (0, 0)}

l∈∂i\{j,k}

hf− (d, µ) = log δcµi ,0 + max




k∈∂i\j 

X

hki (d, µ) − wik +

hli (0, 0)





l∈∂i\j

hb0 = max max hb− (d, µ)
µ6=0 d<0


hf0 = max max  max
µ6=0 d<0

k∈∂i\j ,
l∈∂i\{j,k}

B.


X

hki (d, µ) − wik + hli (−d, µ) +


hmi (0, 0)

m∈∂i\{j,k,l}

Edge-disjoint Steiner trees problem

As for the V-DStP, the Belief Propagation equations for the E-DStP can be computed imposing ψi (di , µi ) =
ψiV (di , µi ) into (13):
X

mij (dij , µij ) =

ψiE (di , µi ) e−β

P

µ

cµ
i I[µi 6=µ]

{dki ,µki }:
k∈∂i\j

Y

nki (dki , µki )

(19)

k∈∂i\j

Instead of considering the cavity messages as in V A, to compute (19) we will first define a partial partition function

Zi =

X

ψiE (di , µi ) e−β

di ,µi

P

µ

cµ
i I[µi 6=µ]

Y

nki (dki , µki )

(20)

k∈∂i

and then calculate the set of messages mij (dij , µij ) through (20) by temporarily setting nji (dji , µji ) = δ−dij ,dji δµij ,µji .
Due to the explicit expression of ψiE message-passing equations become intractable and, therefore, the update step
of the algorithm cannot be efficiently implemented. In the following subsections we overcome this issue by proposing
two different approaches for the computation of (20) where we make use of two different sets of auxiliary variables.
The first formalism relies on “binary occupation” variables that denote, for each node of the factor graph, if edges
incident on it are used or not by any communication; as we will see the associated computation scales exponentially in
the degree of the nodes. The second one consists in a mapping between the E-DStP update equation and a weighted
matching problem over bipartite graphs, that, in the β → +∞, becomes a weighted maximum matching problem
which can be solved efficiently. This implementation scales exponentially with respect to M but it may be more
efficient for vertices with large degrees with respect to the first algorithm.

1.

Neighbors occupation formalism

Suppose of associating with each vertex i ∈ V a vector x = {0, 1} |∂i| . A feasible assignment of these auxiliary
variables is guaranteed if, for every link (i, k) ∈ E incident on i, we will impose xk = 1 if this edge belongs to a tree
(i.e. dki 6= 0 and consequently
Qµki 6= 0) or xk = 0 otherwise (for µki = 0, dki = 0). Variables (di , µi ) must locally
satisfy the following identity k∈∂i I [xk = 1 − δdki ,0 ] = 1 for every node i ∈ V . If we insert this expression in (20)

11
and we sum over all possible assignments of x variables we obtain
P µ
X
XY 

Zi =
ψiE (di , µi ) e−β µ ci I[µi 6=µ]
I xj = 1 − δdji ,0 nji (dji , µji )
=

X

(21)

x j∈∂i

di ,µi

ZxM

(22)

x

where ZxM is defined by computing the following expression for q = M
P µ
X
Y
Zxq ≡
ψiE (di , µi ) e−β µ ci I[µi 6=µ]
I [xk = 1 − δdki ,0 ] nki (dki , µki )
di ,µi
µki ≤q

(23)

k∈∂i

The computation of Zxq is then performed using the following recursion (the equivalence of (24) to (23) is proven in
appendix B)
Zxq =

X


gy0 + gyb + gyf Zyq−1

y≤x
P
−β µ cµ
i

Y

Zx0 = e

(24)

δxj ,0 nji (0, 0)

(25)

j∈∂i

where the auxiliary functions gy0 , gyb , gyf are defined as
q Y
gy0 = e−βci
nki (0, 0)
k∈∂i
yk =0
xk =1

gyb =

XX
d>0 j∈∂i
yj =0
xj =1

gyf = δcqi ,0

Y

nji (−d, q)

[nki (d + 1, q) + nki (0, 0)]

k∈∂i\j
yk =0
xk =1

XX

nji (−d, q)

d>0 j∈∂i
yj =0
xj =1

X

nki (d, q)

k∈∂i\j
yk =0
xk =1

Y

nli (0, 0)

l∈∂i\{j,k}
yl =0
xl =1

Notice that with y ≤ x we denote all possible vectors y = {0, 1} |∂i| satisfying
(
yk =

yk ≤ xk
0

if µki =
6 q
if µki = q

(26)

We can also write the expressions above in the Max-Sum formalism. Define Fi = limβ→+∞
function of Max-Sum messages hij (dij , µij ) = limβ→+∞ β1 log nij (dij , µij ) as

1
β

log Zi and express it as

"
Fi =

max

di ,µi
ψiE (di ,µi )=1

max
x

#
X

log I [xk = 1 − δdki ,0 ] + hki (dki , µki ) −

X

cµi I [µi

6= µ]

(27)

µ

k∈∂i

P
where the function k∈∂i log I [xk = 1 − δdki ,0 ] takes value zero if variables satisfy the constraints or minus infinity
otherwise. As in the BP formulation, we rewrite it as
Fi = max FxM

(28)

x

with
"
FxM

=

max

di ,µi
ψiE (di ,µi )=1

X
k∈∂i

#
log I [xk = 1 − δdki ,0 ] + hki (dki , µki ) −

X
µ

cµi I [µi

6= µ]

12
It is computed recursively from

Fxq = max Fyq−1 + max {h0 , hb , hf }
y≤x
X
X µ
Fx0 = −
hki (0, 0)
ci + log I [x = 0] +
µ

(29)
(30)

k∈∂i

where
h0 =

X

hki (0, 0) − cqi

(31)

k∈∂i
yk =0
xk =1







X



h
(−d,
q)
+
max
[h
(d
+
1,
q)
,
h
(0,
0)]
hb = max max 
ki
li
li


d>0 k∈∂i 

l∈∂i\k
yk =0
xk =1

(32)

yl =0
xl =1




hf = log δ

cqi ,0



X



max
hki (−d, q) + hli (d, q) +
hmi (0, 0)
+ max 

d>0 k∈∂i, l∈∂i, k6=l

m∈∂i\{k,l}
yk =0, yl =0
ym =0
xm =1

xk =1 xl =1

2.

(33)

Mapping into a matching problem

We will develop an alternative method for the computation of the messages of BP and MS update equations, that
M
can lead to an exponential speedup in some cases. Let us introduce an auxiliary vector s ∈ {0, 1, . . . , D} associated
with each vertex of the graph. Components sµ take value in the set of the possible positive depths {1, . . . , D} if
this node is member of communication µ or 0 otherwise. For a node i that is not a root but a member of the
communication µ, there exists exactly one neighbor k such that dik > 0, dki = −sµki µki = µ and for the remaining
ones, dli δµli ,µ = sµki + 1 or dli δµli ,µ = 0, l ∈ ∂i\k. The compatibility function for E-DStP can be expressed as a
function of the new variables as


M
 Y
Y
X X
Y 

δd˜ki ,−sµ
(34)
ψiE (di , µi ) =
δd˜li ,sµ +1 + δd˜li ,0 +
δd˜ki ,0 
µ=1

=

sµ >0 k∈∂i


M
X Y
s



k∈∂i

l∈∂i\k


1 − δsµ ,0

X

µ=1


Y 

δ ˜
dki ,−sµ

k∈∂i



δd˜li ,sµ +1 + δd˜li ,0  +

M
Y
µ=1

l∈∂i\k

δsµ ,0

Y
k∈∂i

δd˜ki ,0




(35)



If now make explicit the dependency on µi we obtain

ψiE (di , µi ) =


X Y X
s

+



δµki ,µ δdki ,−sµ

µ:sµ >0 k∈∂i

Y

Y

Y 

δµli ,µ δdli ,sµ +1 + (1 − δµli ,µ ) +

(36)

l∈∂i\k

(1 − δµki ,µ ) δdki ,sµki




(37)



µ:sµ =0 k∈∂i

Now, introduce:
Zi =

X

ψiE (di , µi ) e−β

d,µi

=

X
s

P

µ

cµ
i I[µi 6=µ]

Y

nki (dki , µki )

(38)

k∈∂i

Qs

(39)

13
where
Qs =

X
di

X

Y

nki (dki , µki )


 Y

e−β

ν
ν6=µ ci I[µi 6=ν]



µ:sµ >0

{µki :sµki >0∨µki =0} k∈∂i

X

P

fkµ + e−β

µ
µ ci

P

k∈∂i

Y
k∈∂i

δµki ,0




Y 

δµli ,µ δdli ,sµ +1 + (1 − δµli ,µ )

fkµ = δµki ,µ δdki ,−sµ

(40)



(41)

l∈∂i\k

Let us concentrate in the computation of Qs for a fixed s. For simplicity of notation, we will assume, unless explicitly
noted, that µ indices run over the set {µ : sµ > 0}. Now as fkµ fkν = 0 if µ 6= ν (because δµki ,µ δµki ,ν = 0), we have
that
δµki ,µ δdki ,−sµ [δµki ,ν δdki ,sν +1 + (1 − δµki ,ν )] = δµki ,µ δdki ,−sµ

(42)

and equivalently
Y

[δµki ,ν δdki ,sν +1 + (1 − δµki ,ν )] =

X

ν

δµki ,ν δdki ,sν +1 + δµki ,0 δdki ,0

(43)

ν

Thus
YX

fkµ =

XY

µ k∈∂i

π

fπµ µ

µ

P
where the sum π runs over all the possible coupling between communications and neighbors of node i. Mathematically we have defined the one-to-one functions π
π : {µ : sµ > 0} → ∂i
with π : µ 7→ πµ . In the following, we will switch to an alternative representation of functions π. If we denote by
tkµ = δk,πµ , for a fixed π we obtain
Y
Y
Y 

fπ µ µ =
δµπµ i ,µ δdπµ i ,−sµ
δµli ,µ δdli ,sµ +1 + (1 − δµli ,µ )
µ

µ

l∈∂i\πµ

!1−Pν tkν
=

Y

X

k∈∂i

ν

Y

δµki ,ν δdki ,sν +1 + δµki ,0 δdki ,0

δµki ,µ δdki ,−sµ

tkµ

µ

0

with P
the convention that 0 = 1. Note that
P the vector t and the function π contain the same information: we have
that k∈∂i tkµ = 1 − δsµ ,0 for each µ and µ tkµ ≤ 1 for each k ∈ ∂i. These two conditions are complete; for a vector
t that satisfies these two constraints, the corresponding function π can be defined naturally. We will have then
"
Zi =

XXY
s

t

k∈∂i

X

#
X

I

µ

"
Y

tkµ = 1 − δsµ ,0

k∈∂i

"
Y X

×
=

e

−βcµ
i I[µi 6=µ]

I

#
X

tkµ ≤ 1 ×

(44)

µ

k∈∂i

#1−Pν tkν
Y

nki (sν + 1, ν) + nki (0, 0)

ν

tkµ

nki (−sµ , µ)

(45)

µ

Rs Z s

(46)

s

where
Rs =

"
Y X

#
nki (sν + 1, ν) + nki (0, 0)
"

Zs =

XY
t

(47)

ν

k∈∂i

µ

e

−βcµ
i I[sµ =0]

I

#
X

k∈∂i

tkµ = 1 − δsµ ,0

"
Y
k∈∂i

I

#
X
µ

tkµ ≤ 1

Y
P
k∈∂i

nki (−sµ , µ)
n
(s
ν + 1, ν) + nki (0, 0)
ν ki

tkµ
(48)

14
The term Zs is the partition function of a matching problem on the complete bipartite graph G = (V = A ∪ B, E = A × B)
with A = ∂i and B = {µ : sµ > 0}, and the energy of a matching is
 (t) =

X

tkµ log P

kµ

nki (−sµ , µ)
− βcµi I [sµ = 0]
ν nki (sν + 1, µ) + nki (0, 0)

In general, the partition function Zs is hard to compute exactly, because it corresponds to the calculation of a
matrix permanent which is computationally intractable. Fortunately, the situation is much easier in the β → ∞ limit:
using hki (−sµ , µ) = β1 log nki (−sµ , µ) and taking the limit β → ∞, the computation of Fi = β1 log Zi reduces to the
evaluation of


1
Fi = max
(log Rs + log Zs )
(49)
s
β
(
)


X
= max
max max hki (sµ + 1, µ) , hki (0, 0) + Fs
(50)
s

µ

k∈∂i

1
β

log Zs we need to solve a weighted maximum matching problem on a bipartite


2
graph which can be done in polynomial time (precisely, it can be performed in O (M + |∂i|) M |∂i| ). Indeed, for

To evaluate the second term Fs =

each assignment of the s we can define the weights wkµ associated with each edge (k, µ) as:

wkµ

(
hki (−sµ , µ) − maxν max {hki (sν + 1, ν) , hki (0, 0)}
=
−cµi

if sµ > 0
if sµ = 0

(51)

and solve

P

s = max
F
(k,µ) wkµ tkµ
P
k∈∂i tkµ ≤ 1

P t ≤ 1
µ kµ

:
∀µ
∀k ∈ ∂i

(52)

That the system in (52) for tkµ ∈ {0, 1} corresponds to a maximum matching problem on a bipartite graph, and
can be solved efficiently through e.g. its linear program relaxation.

C.

The parameter D

The branching formalism introduced in III relies on a parameter D that denotes the maximum allowed distance
between the root and the leaves of any tree. This parameter limits the depth of solution-trees and therefore the
goodness of the results: a small value for D may prevent the connection of some terminals but a large value of D
will slow down the algorithm affecting the converge. Thus this parameter needs to be carefully designed to ensure
good performances. Although there is not a technique able to predict the best value of D, some heuristics have
been proposed in recent works to determine a minimum feasible value of D for the MStP and PCStP [7]. In this
work, we adopt methods described in [8] to find a minimum value of Dµ for each communication µ and we than set
D = maxµ Dµ .
It is clear that the computing cost of both V-DStP and E-DStP strongly depends on the value of D, more precisely
linearly for the V-DStP and the binary occupation formalism and polynomially for the matching problem formulation
for the E-DStP, and could be still prohibit for graph with large diameter. Fortunately, the use of the flat formalism
allows us to reduce the parameter D to D = maxµ |Tµ | being |Tµ | the number of terminals of communication µ. A
proof of this property is reported in [8] for the single tree problem.

VI.

MAX-SUM FOR LOOPY GRAPHS

The goodness of the approximation of marginals is strictly related to the properties of the factor graph over which
we run the Belief Propagation algorithm. BP is exact on tree graphs but nevertheless benefits from nice convergence
properties even on general, loopy, graphs that are locally tree-like [18]. In the framework of the PCStP and multiple

15
trees variants, there are several instances of practical interest, such as square or cubic lattices (2D or 3D graphs)
modelling VLSI circuits, where many very short loops exist and the assumption of negligible correlation among
variables is not satisfied. In many of these cases MS fails to converge in most of the trials or it requires a prohibitive
run-time [8].
We employ here a reinforcement scheme [5, 7] that is able to make the algorithm converge on a tunable amount
of time with the drawback that the solution may be sub-optimal in terms of cost. From the viewpoint of the factor
graph it adds an extra factor to edge-variables that acts as an external field oriented in the direction of the cavity
fields of past iterations. It slightly modifies the original problem into an easier one where a feasible assignment of
variables is more likely to occur. The strength of this perturbation increases linearly in time in a way that, after
few iterations, first inaccurate predictions will be neglected but, after many iterations of MS, it let the algorithm
converge to, hopefully, an assignments of variables satisfying all the constraints. We report in VI A how to modify
the Max-Sum equations for the V-DStP and E-DStP for including the reinforcement factor.
The reinforcement scheme described in the following is generally sufficient to guarantee convergence on random
networks. In practice, however, MS did not converge in some benchmark instances, even adopting the bootstrapping
procedure. In [8] we have shown how to complement the MS equations with heuristics to solve PCStP instances in an
efficient and competitive way. At each iteration we perform a re-weight of node prizes and edge weights according to
temporarily Max-Sum predictions and we then apply a heuristics to find a tree connecting all nodes of the modified
graph. After a pruning procedure, we obtain a pruned minimum spanning tree which is surely a feasible candidate
solution for the PCStP. The motivation is based on the fact that although Max-Sum often outputs inconsistent
configurations while trying to reach the optimal assignment of variables, it still contains some valuable information.
Heuristics have the responsibility of adjusting the assignments of the temporarily decisional variables guaranteeing a
tree-structured solution for any iteration of the main algorithm. Furthermore heuristics results do not depend on the
parameter D of the model and they can provide solution-trees of any diameter. We show in VI B how to generalize
the combination of Max-Sum and heuristics in the case of multiple trees for the V-DStP and the E-DStP.

A.

Reinforcement

At each iteration t of the algorithm we modify the update equations for cavity marginals incident on node i as
t−1
h̃tji (dji , µji ) ∝ htji (dij , µji ) + γt Hji
(dji , µji )

(53)

and we compute the cavity fields as
t−1
t
Hij
(dij , µij ) ∝ htij (dij , µij ) + htji (−dij , µij ) + γt Hij
(dij , µij )

(54)

where γt = tγ0 is linearly proportional to γ0 , the reinforcement factor that governs the strength of the bootstrap. It
is usually very small, of the order of 10−5 not to deviate the dynamics towards the minimum of the energy and thus
affect the goodness of the solution.

B.

Max-Sum based heuristics

At each iteration of the main algorithm we perform a re-weighting of the graph to favor MS temporarily predictions
and we then apply two fast heuristics to find as many spanning trees as the number of communications that we want
to pack. These trees will be carefully pruned in order to decrease the cost of the solution. In this work we design
two different re-weighting schemes for two different heuristics and we refer the interested reader to the single-tree
heuristics explained in [8] for additional details. For each sub-graph µ we apply one of the two schemes as follows.

1.

Shortest Path Tree

For any MS iteration t we compute the auxiliary weights as
t
t
wij
= max Hij
(d, µ)
d6=0

(55)

16
t
t
Notice that since marginals Hij
are normalized, there exists only one assignment of the variables such that Hij
(d∗ , µ) =
∗
0 in correspondence of the most probable state (d , µ). According to (55) edges that are likely to be exploited within
communication µ, i.e. such that d∗ 6= 0, will have zero weights; differently, we penalize edges that, according to MS,
must not be used (for which d∗ = 0) imposing strictly positive weights equal to the MS marginals computed at the
most probable non-zero depth. We then compute the Shortest Paths Tree (SPT) of the modified graph and we prune
the solution tree removing a leaf i if it is not a terminal (for the MStP), and edges (i, j) satisfying wij > cµi (for the
PCStP); we repeat this procedure until we do not find such leaves.

2.

Minimum Spanning Tree

In this scheme we assign auxiliary costs to nodes of the graph according to MS prediction. Let us consider the two
auxiliary functions

(

n
o
P
hi (d, µ) = maxk∈∂i htik (−d, µ) + l∈∂i\k max [htli (d + 1, µ) , htli (0, µ)]
P
hi (0, µ) = k∈∂i htki (0, µ) − cµi

for d > 0

(56)

A node satisfying maxd>0 hi (d, µ) < hi (0, µ) will be penalized assigning to edges incident on it a large cost C. We
then apply the Minimum Spanning Tree (MST) algorithm to the modified graph and we prune the solution as in the
case of the SPT.
Heuristics are applied to the graph for all the communications providing, as a feasible solution for E-DStP or
V-DStP, a superposition of single-tree solutions. Notice that heuristics are sequentially applied, i.e. we consider one
communication at the time, and depending on we are dealing with V-DStP or E-DStP, edges (and Steiner nodes for
the V-DStP) selected in the first spanning trees cannot be further used by the successive applications. To overcome
this problem, we add an erasing step before the application of each heuristics in which we delete edges (and eventually
Steiner nodes) used by other communications. For V-DStP we only need to cut edges incident on terminals of other
sub-graphs to satisfy nodes-disjoint constraints. Unfortunately such strong edge cutting procedure may lead to a
graph with disconnected components or a graph in which the terminals that we aim at connecting may be isolated.
In these scenarios we cannot find further trees able to span the modified graph and thus this heuristic approach fails.
One way of preventing this problem is to randomize the order of the trees over which we apply the heuristics.

VII.

NUMERICAL RESULTS

In this section we report results for several experiments on synthetic networks and on benchmark, real-world,
instances for the VLSI. In all the cases we will solve the V-DStP or the E-DStP where terminals have infinite prizes,
i.e. the MStP variant, and a predefined root is selected for each sub-graph. The synthetic networks we chose are
fully connected, regular or grid graphs, whose properties will allow us to underline the main features of the models
and formalisms introduced in this work. In particular, by means of the fully connected graphs we will illustrate
the improvements carried by Max-Sum against a “greedy” search of the solutions in which we will make use of the
single-tree MS algorithm for the MStP. Furthermore, regular graphs allow us to verify the different scaling of the
running time with respect to the degree of the graph of the two algorithms presented for the E-DStP. Motivated
by their importance on technological applications, namely in the the design of VLSI, we also show some results
on grid, both synthetic and real-word, graphs: here we will underline the improvements carried by the flat model.
Generally, energies are averaged over several instances, meaning different realizations of the weighting of the edges
and assignment of the terminals, of the same graph. To measure the energy gap of the solutions found by the two
E −E
different procedures, for instance “x” and “y” algorithms, we measure the quantity xEy y assuming that Ex and Ey
are the energies of the solutions found by algorithm “x” and “y” respectively. If the gap is positive (negative) the “x”
(“y”) algorithm outperforms the other one.
We underline that, due to the intrinsic difficulty of the problem, there are very few (exact or approximate) results in
literature and few available algorithms to use for the comparison. In the case of VLSI circuits, we report the solution
costs of a state-of-the art linear programming technique for the V-DStP.

17
A.

Fully connected graphs

Here we report results for the V-DStP on fully connected graphs where we aim at packing M = 3 trees. We
compare our performances against the following “greedy” procedure: we apply MS algorithm for the single-tree MStP
to each communications, one at the time, and we compute the “greedy” energy as the sum of energies of single-tree
solutions. As in the case of the heuristics described in VI B we must perform a pre-processing of the graph before
the application of the single-tree algorithm in order not to use edges in more than one sub-graphs. Notice that this
“greedy” procedure is actually as hard as the packing problem, since even solving the MStP belongs to NP-hard class
of problem; nevertheless this procedure will be useful to underline the benefits carried by the parallel (packing) search
against the “greedy” and sequential one.
We deal with fully connected graphs because here the existence of a trivial solution of the packing problem, consisting
in a chain of terminal nodes, is always guaranteed. We perform two different experiments: we first fix the size of
the graphs (500 nodes) and we study how energies and gaps change for an increasing number of terminals nodes.
T
Secondly, we fix the fraction of terminals per communication, more precisely for α = Nµ = 0.08, µ ∈ {1, 2, 3} and we
compare the performances as we increase the size of the graphs (from 100 to 700 nodes). We run both algorithms
with fixed parameters D = {3, 5, 10} and fixed reinforcement factor γ0 = 10−5 .
1.

Uncorrelated edge weights

These experiments are performed on fully connected graphs where weights associated with edges are independently
and uniformly distributed random variables in the interval (0, 1). In this scenario, energies obtained by the greedy
procedure are always larger then the ones achieved by the parallel search, for all values of the number of terminals
and for any value of the parameter D used, as it is suggested by the plot of the gaps (right plot) in Figure 3 on page
17. Notice that the gaps are slightly greater then zero suggesting that solutions found by the two methods are very
similar in terms of energy cost, as reported in the plot in Figure 3 on page 17, left panel.
(b) 8

(a) 7

Energy

5

D3
D5
D 10

7

6

Gap(Greedy, MS) %

6

Greedy D 3
Max-Sum D 3
Greedy D 5
Max-Sum D 5
Greedy D 10
Max-Sum D 10

4

3

2

1

5

4

3

2

1

0
20

40

60

80

100

120

140

0
20

Number of terminals

40

60

80

100

120

140

Number of terminals

Figure 3. Energy (a) and energy gap (b) of the solutions of Max Sum and Greedy algorithm as functions of the number of
terminals. The test instances are fully connected graphs of 500 nodes with uncorrelated edge weights. Gaps reported in panel
(b) are always positive suggesting that solutions found by the global search are cheaper in terms of cost then the greedy ones.

2.

Correlated edge weights

To underline the benefits carried by the optimized strategy, we run reinforced and greedy reinforced Max-Sum on
complete graphs with correlated edge weights. With each node i we assign a uniformly distributed random variable xi

18
in the interval (0, 1) and for each edge (i, j) we pick a variable yij ∈ (0, 1). Then an edge (i, j) will be characterized by
a weight wij = xi xj yij . In this scenario, we expect that the cheapest edges will be chosen by the “greedy” algorithm
for the solution of the first trees and, as we proceed with the sequential search, the algorithm will become the more
and more forced to use the remaining expensive edges. In fact, as shown in Figure 4 on page 18, the gaps notably
increase of one order of magnitude for most of the number of terminals considered in these experiments.
Notice that energies encountered for D = {5, 10} are very close to one another suggesting that a further increasing
of the parameter D, and thus of the solution space, will not lead to a significant improvement of the solutions.
(b) 80

0.25

0.2

Greedy D 3
Max-Sum D 3
Greedy D 5
Max-Sum D 5
Greedy D 10
Max-Sum D 10

60

0.15

Energy

D3
D5
D 10

70

Gap(Greedy, MS) %

(a)

0.1

50

40

30

20

0.05

10

0

0
40

60

80

100

120

140

40

Number of terminals

60

80

100

120

140

Number of terminals

Figure 4. Energy (a) and energy gap (b) for Max Sum results against Greedy results as functions of the number of terminals
for correlated edge weighting. The energy gaps of panel (b) are positive and notably large.

3.

Fixed fraction of terminals

To study the performances in the asymptotic limit, namely for N → +∞ ,Tµ → +∞ for each communication µ and
constant α, we attempted the solution of V-DStP on complete graphs having a fixed fraction of terminals α = 0.08
and for an increasing number of nodes N . Although non-rigorous, this procedure can suggest us the behavior of the
energies and the energy gaps in the large N limit. As reported in Figure 5 on page 19 panel (a), when the number of
nodes reaches N ∈ [500, 700], the energy of both Max Sum and greedy solutions, for all values of D, seems to stabilize
to a constant value. As a consequence, as plotted in Figure 5 on page 19 panel (b), also energy gaps fluctuates around
a fixed value that seems to be different if one considers D = 3 or D = {5, 10}.

B.

Regular graphs

In V B we have seen how to deal with the update equations of BP and MS algorithms for the E-DStP with the
help of two different auxiliary set of variables. Although the final expressions of the equations are very different,
the energies obtained by both algorithms must be identical; the only differences rely on the computational cost
that strongly depends on the properties of the graph, precisely on the degree of the nodes of the graph and on the
number of communications. To underline these two features of the neighbors occupation formalism and matching
problem mapping, that from now will be denoted as NeighOcc and Matching algorithms, we perform two different
experiments on regular, fixed degree, graphs for different values of the degree and of the number of sub-graphs. For
these simulations we have fixed the values of the parameter D = 10 and the reinforcement factor γ0 = 10−4 .

19
(a)

(b) 80

0.07
0.065

75

0.06
70

Gap(Greedy, MS) %

0.055

Energy

0.05
0.045
0.04
0.035
Greedy D 3
Max-Sum D 3
Greedy D 5
Max-Sum D 5
Greedy D 10
Max-Sum D 10

0.03
0.025

100

200

60

55

50

45

D3
D5
D 10

40

0.02
0

65

300

400

500

600

700

800

0

100

200

300

400

500

600

700

800

Number of nodes

Number of nodes

Figure 5. Energy (a) and energy gap (b) for Max Sum against Greedy results as functions of the number of nodes for a fixed
fractions of terminals (per communication)α = 0.08.

1.

Energy as a function of the degree

Similarly to the experiments in VII B 2, here we consider regular graphs of N = 50 nodes containing M = 3 subgraphs for four possible degrees d ∈ {3, 4, 5, 6}. The energies provided by NeighOcc and Matching and plotted in
Figure 6 on page 20, panel (a), can be statistically considered the same, as for the fixed degree experiment shown
before. Here the computational costs (panel (b) and (c) of Figure 6 on page 20) scales exponentially only for the
NeighOcc (as it is remarked by the linear trend in the semi-log plot) while it scales polynomially for the Matching
formalism as predicted by the analysis on the update equations in V B.

2.

Energy as a function of the number of communications

In this experiment we try to solve the E-DStP on two sets of regular graphs of N = 50 nodes having fixed degree
4, for an increasing number of trees. Each communication has the same number of terminals T = 3. As shown in
Figure 6 on page 20, panel (d), the energy costs of the solutions provided by NeighOcc and Matching algorithms are
almost identical as we expected. At the same time, the computing time plotted in Figure 6 on page 20, panels (e)
and (f), shows that the Matching procedure needs a time that scales exponentially, i.e. linearly in an log-scale plot,
on the number of sub-graphs while it becomes polynomial for the NeighOcc algorithm.

C.

Grid graphs

This section is devoted to the illustration of results of both V-DStP and E-DStP on 2D and 3D lattices. The first
experiments are performed on synthetic 3D lattices of dimension 5 × 5 × 5 containing N = 125 nodes. Here we fix
the number of communications M ∈ {2, 3, 4} and we study how energies behave when the number of terminals T per
communication changes in the range [3, 6]. For the V-DStP and the E-DStP (only in the NeighOcc formalism) we
compare the results provided by the branching and flat models. While the parameter D can be arbitrary large for
the branching model, we keep the value D = T for the flat one since, as discussed in III 2, it is sufficient to explore
all the solution space. In the second part of this section we comment the performances of the MS algorithm and of
the MS-based heuristics presented in VI B applied to several benchmark instances for the design problem of VLSI
circuits.

20
(a) 10

7
6
5
4
3
3

4

5

6

7

Time [s] - NeighOcc

1 10 4

8

2

7 10 3

5 10

3

3 10

3

2

3

4

5

6

7

2

8
7
6
5
4
3
4

3

5

6

10

6

10

5

10

4

10

3

10

2

4

5

6

7

Degree 4

Time [s] - NeighOcc

Time [s] - Matching

9

3

2

Degree 4

10

2

10

(f) 3 10 3

Degree 4 Matching
Degree 4 NeighOcc

1

10 3

Degree

(e) 10 7

(d) 12
11

10 4

Degree

Degree

Energy

M3

M3

Time [s] - Matching

9

Energy

(c) 10 5

(b)
M 3 Matching
M 3 NeighOcc

1

2

3

4

5

M

M

6

2 10 3

1 10 3

1

2

3

4

5

6

M

Figure 6. In panel (a), energy of the solutions for the E-DStP on regular graphs of 3 communications as a function of the
degree. In (b) and (c) running time of Matching and NeighOcc algorithms as a functions of degree. In panel (d): energy of the
solutions for the E-DStP on regular graphs as a function of the number of packed trees. Panels (e) and (f): running time of
the Matching and NeighOcc algorithms as a functions of the number of communications. In all plots energies of the solutions
are almost the same, but the computing time dramatically differ as we are using the Matching formalism or the NeighOcc
algorithm.

1.

Branching and flat models for the V-DStP e E-DStP (neighbors occupation formalism)

As shown in Figure 7 on page 21, left panel, the energies of the solutions found by the flat model for the V-DStP
are always smaller than the energies found by the branching one. We underline that, as plotted in the right panel
of Figure 7 on page 21, the flat version of MS equations has the advantage of converging in a running time that is
always smaller then the one needed by the branching model. This is reasonable as the parameter D, which linearly
influences the computation time of both algorithms, is often greater (on average D = 8) for the branching model than
the one fixed for the flat representation.
A different behavior is observed for the resolution of the E-DStP on grids using our two models. As remarked in
Figure 8 on page 21, left panel, energies found by the flat and branching representations are comparable; here the
depth used by the branching model, on average equal to D = 8 and D = 9 for T ∈ [3, 4] and T ∈ [5, 6] respectively,
probably suffices to explore the same solution space considered by the flat formalism for smaller D. Still, the flat
model is preferable as it requires a computing time that is smaller then the one needed by the branching model for
all the cases we have considered.

2.

V-DStP for VLSI circuits

In this section we report several results for standard benchmark instances of circuit layout where we solve the
V-DStP. Instances are 3D grid graphs modelling VLSI chips where we pack relatively many trees, usually 19 or 24,
each of which typically contains few terminal nodes (3 or 4). Such grid graphs can be seen as multi-layers graphs

21
(a) 25

(b) 10 2

M 2 MS branching
M 2 MS flat
M 3 MS branching
M 3 MS flat
M 4 MS branching
M 4 MS flat

20

Time [s]

Energy

15

10

10

1

M 2 MS branching
M 2 MS flat
M 3 MS branching
M 3 MS flat
M 4 MS branching
M 4 MS flat

5

0
2

3

4

5

6

10

7

0

2

3

4

T

5

6

7

T

Figure 7. Energy (a) and computational time (b) as a function of the number of terminals per communications for 3D grid
graphs, V-DStP variant.
(b) 10 4

(a) 22

M 2 MS branching
M 2 MS flat
M 3 MS branching
M 3 MS flat
M 4 MS branching
M 4 MS flat

20

18

Time [s]

Energy

16

14

12

10

3

10

M 2 MS branching
M 2 MS flat
M 3 MS branching
M 3 MS flat
M 4 MS branching
M 4 MS flat

8

6

4
2

3

4

5

T

6

7

10

2

2

3

4

5

6

7

T

Figure 8. Energy (a) and computational time (b) as a function of the number of terminals per communications for 3D grid
graphs, E-DStP variant.

where we allow two different kinds of connections. In the multi-crossed layers, each node is connected to all its
possible neighbors in all directions: the resulting graphs are cubic lattices. The multi-aligned layers are similar to the
multi-crossed ones but in each layer we allow only connections in one direction, either east-to-west or north-to-south
[12]. For sake of simplicity, consider a cubic lattice in a three dimensional Cartesian coordinate system: depending
on the value of the z−coordinate, the allowed connections will be present in directions parallel to the x or to the
y axes. In Table I on page 22 we first report some information (type of the layers, size, number of sub-graphs and

22
total number of terminals) concerning each instance and our results. We show the energies achieved by reinforced
Max Sum along with the ones of the two heuristics described in VI B; in analogy with [8], we label as “J” heuristics
that performs a modified SPT and as “N” if instead we use the MST. Energies obtained using the flat model are
labeled as “(f)” while if nothing is specified or “(b)” is used, we made use of the branching representation. Results are
compared with respect to the ones obtained through state-of-the-art linear programming (LP) techniques [12] which
is able, for these particular instances, to find the optimal solutions. The sign “-” denotes that no solution has been
found. As shown in Table I on page 22, the gaps are always smaller than 4% and in two cases, for the multi-aligned
augmenteddense-2 and terminalintensive-2 instances, we reach the same performances of LP, obtaining the optimal
solutions. We stress that these graphs are very loopy and far from being locally tree-like but nevertheless we achieve
good performances thanks to the reinforcement procedure along with the introduction of the modified heuristics.

augmenteddense-2
augmenteddense-2
dense-3
difficult-2
difficult-2x
difficult-2y
difficult-2z
modifieddense-3
moredifficult-2
pedabox-2
terminalintensive-2

Type
Multi-aligned
Multi-crossed
Multi-crossed
Multi-aligned
Multi-aligned
Multi-aligned
Multi-aligned
Multi-crossed
Multi-aligned
Multi-aligned
Multi-aligned

Size
16x18x2
16x18x2
15x17x3
23x15x2
23x15x2
23x15x2
23x15x2
16x17x3
22x15x2
15x16x2
23x16x2

M
19
19
19
24
24
24
24
19
24
22
24

Ttot
Heur. “J”
Heur. “N” Rein. Max Sum
59 504 (b) 506 (f)
507
508 (f) 504 (b)
59
503
59
487
488
485
66
535
538
538
66
560
66
4776
4829
4816
66
1060
1063
1061
59
492
496
495
65
542
542
546
56
405
405
405
77 596 (f) 599 (b)
617
620

LP (opt)
504
498
464
526
unknown
unknown
unknown
479
522
390
596

Gap %
0%
1.0 %
4.0 %
1.7 %

2.6 %
3.8 %
3.8 %
0%

Table I. Results for circuit layout instances

VIII.

SUMMARY OF RESULTS

Using Max-Sum algorithm, we have explored through simulations some interesting theoretical questions in random
graphs which we summarize here. Simulations (up to N = 700, or around 2 × 105 edges) suggest that for the Steiner
Tree Packing problem on complete graphs with uniform independent weights, the energy converges to a constant value
if the fraction of terminal vertices is kept constant, in agreement with known results for single Steiner trees [2].
We have observed a non-negligible gap (up to 7% in the solution energy and increasing with tree depth) between a
greedy solution (which is commonly used by practitioners and consists in sequentially optimizing each communication
and removing its used components from the graph) and the joint optimum computed by MS. Interestingly this gap
is greatly expanded (up to 80% in experiments) with weights that are positively correlated. For the edge-disjoint
problem, we have compared all model variants on random regular graphs with various parameters (degree, number of
terminals, number of trees), confirming the convenience of each of them in a different parameter region. Simulations
on regular lattice graphs give qualitatively similar results.
Finally, we have attempted to optimize a set of publicly available benchmark problems (including 2d and 3d tree
packing problems), some of which have known optimum. Results are encouraging, as the solutions provided by the
Max-Sum algorithm show a gap no larger than 4% in all cases (0% in some cases) when the optimum is known. We
expect this gap to be generally independent of the problem size, which suggests that this strategy could be extremely
useful for large-scale industrial problems.

23

(b) terminalintensive-2

(a) augmenteddense-2 multi-aligned

(c) augmenteddense-2 multi-crossed

(d) dense-3

Figure 9. Examples of solutions for the V-DStP on VLSI circuits for multi-aligned, (a) and (b) figures, and multi-crossed, (c)
and (d) figures, layouts

ACKNOWLEDGMENTS

We acknowledge Fondazione CRT for project SIBYL under the initiative “La ricerca dei Talenti”. We warmly thank
F. Altarelli for interesting discussions and R. Zecchina for suggesting this topic.

[1] Fabrizio Altarelli, Alfredo Braunstein, Luca Dall’Asta, Caterina De Bacco, and Silvio Franz. The edge-disjoint path
problem on random graphs by message-passing. CoRR, abs/1503.00540, 2015.
[2] Omer Angel, Abraham D. Flaxman, and David B. Wilson. A sharp threshold for minimum bounded-depth and boundeddiameter spanning trees and Steiner trees in random networks. Combinatorica, 32(1):1–33, January 2012.
[3] C. De Bacco, S. Franz, D. Saad, and C. H. Yeung. Shortest node-disjoint paths on random graphs. J. Stat. Mech.,
2014(7):P07009, July 2014.
[4] M. Bailly-Bechet, C. Borgs, A. Braunstein, J. Chayes, A. Dagkessamanskaia, J.-M. François, and R. Zecchina. Finding
undetected protein associations in cell signaling by belief propagation. Proceedings of the National Academy of Sciences,
108(2):882–887, January 2011.
[5] M. Bayati, C. Borgs, A. Braunstein, J. Chayes, A. Ramezanpour, and R. Zecchina. Statistical Mechanics of Steiner Trees.

24
Physical Review Letters, 101(3):037208, July 2008.
[6] Mohsen Bayati, A. Braunstein, and Riccardo Zecchina. A rigorous analysis of the cavity equations for the minimum
spanning tree. Journal of Mathematical Physics, 49(12):125206, 2008. Cited by 0012.
[7] Indaco Biazzo, Alfredo Braunstein, and Riccardo Zecchina. Performance of a cavity-method-based algorithm for the
prize-collecting Steiner tree problem on graphs. Phys. Rev. E, 86:026706, August 2012.
[8] Alfredo Braunstein and Anna Muntoni. Practical optimization of Steiner trees via the cavity method. Journal of Statistical
Mechanics: Theory and Experiment, 2016(7):073302, 2016.
[9] M. Burstein and R. Pelavin. Hierarchical Wire Routing. IEEE Transactions on Computer-Aided Design of Integrated
Circuits and Systems, 2(4):223–234, October 1983.
[10] J. P. Cohoon and P. L. Heck. BEAVER: a computational-geometry-based tool for switchbox routing. IEEE Transactions
on Computer-Aided Design of Integrated Circuits and Systems, 7(6):684–697, June 1988.
[11] M. Grötschel, A. Martin, and R. Weismantel. The steiner tree packing problem in VLSI design. Mathematical Programming,
78(2):265–281.
[12] Nam-Dūng Hoáng and Thorsten Koch. Steiner tree packing revisited. Mathematical Methods of Operations Research,
76(1):95–123, 2012.
[13] Richard M Karp. Reducibility among combinatorial problems. Springer, 1972.
[14] Ivana Ljubić, René Weiskircher, Ulrich Pferschy, Gunnar W. Klau, Petra Mutzel, and Matteo Fischetti. An Algorithmic
Framework for the Exact Solution of the Prize-Collecting Steiner Tree Problem. Math. Program., 105(2-3):427–449, October
2005.
[15] Wing Kwong Luk. A greedy switch-box router. INTEGRATION, the VLSI journal, 3(2):129–149, 1985.
[16] Marc Mézard and Andrea Montanari. Information, Physics, and Computation. Oxford University Press, January 2009.
[17] Marc Mézard and Giorgio Parisi. The Cavity Method at Zero Temperature. Journal of Statistical Physics, 111(1-2):1–34,
April 2003.
[18] Y. Weiss and W.T. Freeman. On the optimality of solutions of the max-product belief-propagation algorithm in arbitrary
graphs. Information Theory, IEEE Transactions on, 47(2):736–744, February 2001.

25
Appendix A: Message-Passing equations for V-DStP

Consider the compatibility function for node i
Y

ψiV (di , µi ) =

δµki ,0 δdki ,0 +

X

δcµi ,0

µ

XX

Y

δµki ,µ δdki ,−d

µ=1 d>0 k∈∂i

k∈∂i

+

M X X
X

(δµli ,µ δdli ,d+1 + δµli ,0 δdli ,0 ) +

l∈∂i\k

X

δµ,µki δ−d,dki

d>0 k∈∂i

(A1)
Y

δµ,µli δdli ,d

l∈∂i\k

δµmi ,0 δdmi ,0

m∈∂i\{k,l}

For sake of simplicity we split ψiV (d, µ) in:
(1)

(2)

(3)

ψiV (di , µi ) = ψi (di , µi ) + ψi (di , µi ) + ψi (di , µi )

(A2)

where
(1)

ψi (di , µi ) =

Y

δµki ,0 δdki ,0

(A3)

k∈∂i
(2)

ψi (di , µi ) =



M X X
X



δµki ,µ δdki ,−d

µ=1 d>0 k∈∂i
(3)
ψi

(di , µi ) =

X

δcµi ,0

µ

Y

(A4)

l∈∂i\k

XX

X

δµ,µki δ−d,dki

d>0 k∈∂i

(1)

(δµli ,µ δdli ,d+1 + δµli ,0 δdli ,0 )

l∈∂i\k

(2)

Y

δµ,µli δdli ,d

δµmi ,0 δdmi ,0

(A5)

m∈∂i\{k,l}

(3)

Consider the action of ψi (di , µi ), ψi (di , µi ) and ψi (di , µi ) inside (13) and compute
(1)
(2)
(3)
mij (dij , µij ) = mij (dij , µij ) + mij (dij , µij ) + δcµi ,0 mij (dij , µij ).
(1)

X

mij (dij , µij ) =

e−β

{dki ,µki }:
k∈∂i\j
P
−β µ cµ
i

=e

P

µ

cµ
i I[µi 6=µ]

Y

δµli ,0 δdli ,0

l∈∂i

Y

δµij ,0 δdij ,0

Y

nki (dki , µki )

(A6)

k∈∂i\j

mki (0, 0)

(A7)

k∈∂i\j
(2)

mij (dij, µij ) =

X

e−β

µ
µ ci I[µi 6=µ]

P

µ


X
d>0

Y

δdji ,−d δµji ,µ



[nki (d + 1, µ) + nki (0, 0)] +

(A8)

k∈∂i\j

+ δdji ,d+1 δµij ,µ + δdji ,0 δµji ,0

 X

Y

nki (−d, µ)

k∈∂i\j

l∈∂i\{j,k}



[nli (d + 1, µ) + nli (0, 0)]


(A9)


(3)

mij (dij , µij ) =

XX

δµ,µji δdji ,−d

µ d>0

+δµ,µji δdji ,d

X
k∈∂i\j

X

Y

nki (d, µ)

nki (−d, µ)

k∈∂i\j

nli (0, 0) +

(A10)

l∈∂i\{j,k}

Y

nli (0, 0) +

(A11)

l∈∂i\{j,k}


+δµ,µji δdji ,0

X
k∈∂i\j

nki (d, µ)

X
l∈∂i\{j,k}

nli (−d, µ)

Y

nmi (0, 0)

(A12)

m∈∂i\{k,l,j}

If now we use that dji = −dij and µij = µji we can write the following set of equations:


Q
P
Q
+ nki (0, 0)] + δcµi ,0 k∈∂i\j nki (d, µ) l∈∂i\{j,k} nli (0, 0)
∀d > 0

mij (d, µ) = Pk∈∂i\j [nki (d + 1, µ)Q
P
Q
mij (d, µ) = k∈∂i\j nki (d + 1, µ) l∈∂i\{j,k} [nli (d, µ) + nli (0, 0)] + δcµi ,0 k∈∂i\j nki (d, µ) l∈∂j\{j,k} nli (0, 0)
∀d < 0
P

Q
P
P
P
Q

−β µ cµ
i
mij (0, 0) = e
< −AN N
k∈∂i\j nki (0, 0) +
µ6=0
d<0
k∈∂i\j nki (d + 1, µ)
l∈∂i\{j,k} [nli (d, µ) + nli (0, 0)]
(A13)

26
For d = µ = 0

mij (0, 0) = e−β

P

cµ
i

µ

Y

nki (0, 0) +

XX X

X

nki (d, µ)

µ6=0 d<0 k∈∂i\j

Y

nki (d + 1, µ)

µ6=0 d<0 k∈∂i\j

k∈∂i\j

+

XX X

[nli (d, µ) + nli (0, 0)] +

l∈∂i\{j,k}

Y

nli (−d, µ)

l∈∂i\{j,k}

nmi (0, 0)

m∈∂i\{k,l,j}

Appendix B: Recursive expression of Z q for the E-DStP

From Eq. (23)

X

Zxq =

ψiE (di , µi ) e−β

P

µ

cµ
i I[µi 6=µ]

di ,µi
µki ≤q, ∀k∈∂i

Y

I [xk = 1 − δdki ,0 ] nki (dki , µki )

(B1)

k∈∂i

we underline the possible contribution to a communication q from at least one on the neighbors of i as

Zxq =

X

ψiE (di , µi ) e−β

cµ
i I[µi 6=µ]

P

µ

d,i µi
µki ≤q

Y

I [xk = 1] I [dki 6= 0] nki (dki , µki )

k∈∂i:
µki =q

Y

I [xk = 1 − δdki ,0 ] nki (dki , µki )

k∈∂i:
µki ≤q−1

(B2)
Consider a vector x such that there exists at least one component xk = 1 for dki 6= 0, µki = q and eventually other
components different from zero assigned to one of the possible sub-graph µ ≤ q − 1. This vector can be seen as the
superposition of all vectors y ≤ x, that is, all vectors having at most the same number of non-zeros of x and the
component yk = 0 each time µki = q; all remaining components must satisfy yk0 = 1 − δdk0 i ,0 for µk0 i ≤ q − 1. Thus:
X X

Zxq =

q

e−βci I[µi 6=q] ψiq (di , µi )

di ,µi y≤x
µki ≤q

×

Y

nki (dki , µki ) δµki ,q ×

(B3)

k∈∂i:
yk =0,
xk =1
p

Y

Y

p≤q−1

k∈∂i:
µki ≤q−1

e−βci I[µi 6=p] ψip (di , µi ) I [yk = 1 − δdki ,0 ] nki (dki , µki )

(B4)

where we have made use of the expression of ψiE in (9). If we now collect the sum over y ≤ x and we explicitly use
the constraints on depth and communication variables we find

Zxq =





X
X

q

e−βci I[µi 6=q] ψiq (d, µ)



X

Y

p

e−βci I[µi 6=p] ψip (d, µ)

di ,µi p≤q−1
µki ≤q−1

=

X
y≤x

nki (dki , µki ) δµki ,q ×

(B5)

k∈∂i
yk =0
xk =1

di ,µ
y≤x 


µki ≤qi

×

Y


gy0 + gyb + gyf Zyq−1

Y
k∈∂i

I [yk = 1 − δdki ,0 ] nki (dki , µki )






(B6)




(B7)

27
where
q

gy0 = e−βci

Y

nki (0, 0)

k∈∂i
yk =0
xk =1

gyb =

XX

Y

nji (−d, q)

d>0 j∈∂i
yj =0
xj =1

gyf = δcqi ,0

XX

[nki (d + 1, q) + nki (0, 0)]

k∈∂i\j
yk =0
xk =1

nji (−d, q)

d>0 j∈∂i
yj =0
xj =1

X

nki (d, q)

k∈∂i\j
yk =0
xk =1

Y

nli (0, 0)

l∈∂i\{j,k}
yl =0
xl =1

In the special case in which no communications is flowing within the graph, that is for q = 0, we must impose the
value of Zx0 through
Zx0 = e−β

P

µ

cµ
i

I [x = 0]

Y
j∈∂i

nji (0, 0)

(B8)

