{'fraction_non_alphanumeric': 0.044811979661395766, 'fraction_numerical': 0.011789685422137788, 'mean_word_length': 3.6974628171478567, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 65, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Message Passing Algorithm for the Problem of Path Packing in Graphs

We consider the problem of packing node-disjoint directed paths in a directed graph. We consider a variant of this problem where each path starts within a fixed subset of root nodes, subject to a given bound on the length of paths. This problem is motivated by the so-called kidney exchange problem, but has potential other applications and is interesting in its own right. We propose a new algorithm for this problem based on the message passing/belief propagation technique. A priori this problem does not have an associated graphical model, so in order to apply a belief propagation algorithm we provide a novel representation of the problem as a graphical model. Standard belief propagation on this model has poor scaling behavior, so we provide an efficient implementation that significantly decreases the complexity. We provide numerical results comparing the performance of our algorithm on both artificially created graphs and real world networks to several alternative algorithms, including algorithms based on integer programming (IP) techniques. These comparisons show that our algorithm scales better to large instances than IP-based algorithms and often finds better solutions than a simple algorithm that greedily selects the longest path from each root node. In some cases it also finds better solutions than the ones found by IP-based algorithms even when the latter are allowed to run significantly longer than our algorithm.

In this paper we consider the problem of packing node-disjoint directed paths into a directed graph, with each path starting within a designated subset of "root" nodes. For a given maximum path length, our objective is to include as many nodes as possible in paths. This problem is motivated by the kidney exchange problem (KEP) studied in [2], in which donors and recipients must be matched together to maximize the number of donations. In some situations these donations are performed in chains starting from a lone donor, and logistical concerns encourage limiting the length of chains. This problem also has other potential applications such as coordinating chained transactions, and is an interesting and hard combinatorial optimization problem in general.
We present a new algorithm for this problem based on the message passing/belief propagation technique (BP for short). This algorithm depends on a new representation of the problem that allows us to create a graphical model which embeds the optimization problem. After describing this representation we will show that a direct application of standard BP results in poor scaling, so we provide an efficient implementation that provides significant scaling advantages. Letting n be the number of nodes in the graph, K be the maximum path length, and ∆ be the maximum degree, the efficient implementation can perform a full iteration of the algorithm in O (n∆K) computations, while a direct implementation of belief propagation requires O n 2 ∆ 6 K 2 .
We compare our algorithm to several alternatives and show that it consistently scales well and often finds better solutions than other fast algorithms. In particular, we compare our efficiently implemented BP to a simple greedy algorithm, called Greedy from now on for simplicity, and an existing Integer Programming (IP) algorithm used in [2] which we call KEP. The Greedy algorithm which we consider makes greedy selections only between constructed paths and not during the construction of individual paths, as at every root node it searches exhaustively for the longest path which can be added to the solution starting from that node. We also use our new representation of the problem to derive a secondary novel algorithm based on IP, which we call Parent-Child-Depth (PCD). Comparisons are made on a variety of randomly created graphs as well as on several real world networks. Our random graphs are generated with each edge existing independently with probability c/n for various values of c. We find that for n = 1000 and n = 10000 with a path length bound of 5 nodes, with 10% to 25% of nodes designated as root nodes, and c between 2 and 5, our message passing algorithm consistently finds more nodes than Greedy while running in comparable time. In many cases BP and Greedy run significantly faster than the IP algorithms, which are cut off after a time limit if they have not yet found the optimal solution, and in several of these cases BP also finds solutions with more nodes than any other algorithm. We find that BP also performs well relative to alternative algorithms when the path length bound is increased to 10 or 15 nodes. In the regime where paths as long as 15 nodes are allowed, PCD often finds the optimal solution in reasonable time while KEP almost never does. We also provide numerical results for real world networks, with sizes ranging from 5518 nodes to 260982 nodes. In all but the smallest network the IP algorithms produce no useful results, and in most graphs BP finds better solutions on average with respect to a random choice of root nodes than Greedy, while running in comparable time.
Our approach to this problem is motivated by previous work applying message passing techniques to the problem of prize-collecting Steiner trees (PCST) done by Bayati, et. al. [4] (also used by Bailly, et. al. [3]). In particular, the PCST problem does not have a natural representation as a graphical model so the authors in [4] design a new representation of the problem that leads to a graphical model and then show how to implement BP on that model. This paper uses a similar approach for the path packing problem.
In another relevant paper, Altarelli, et. al. [1] apply message passing techniques to edge-disjoint paths. They directly derive a message passing algorithm for general graphs from an iterative cost calculation scheme valid for trees and apply it to locally tree-like graphs and find that their algorithm performs better in terms of paths found than various alternatives. Unlike this work, in our model paths may end at any non-root node (as we do not consider a lone root node to constitute a path in itself, every path must include at least one edge and thus one non-root node).
Another example of applying a message passing approach while making major modifications to standard methods is [5]. The authors introduce a new type of graphical model, called a memory factor network, along with a message passing style algorithm to perform inference on such a model. Designing the message passing algorithm to match the structure of the problem provides efficiency advantages over variants of belief propagation and also provides convergence guarantees.
The remainder of this paper is organized as follows: In section 2 we describe the problem formulation, in section 3 we detail the core of our message passing algorithm, in section 4 we detail our method for efficiently implementing this algorithm, in section 5 we discuss how to translate the result of the message passing into a solution to the original problem, in section 6 we discuss alternative algorithms, and in section 7 we provide numerical results comparing our algorithm to these alternatives.

Problem setup.
The problem description consists of a directed graph G with node set [n] = {1, 2, . . . , n}, an edge set E, and a positive integer K. Let E ′ be the undirected version of the edge set of G. We assume that G contains no isolated nodes, so for every i ∈ [n] there exists some j ∈ [n] such that (i, j) ∈ E ′ . The nodes [n] are partitioned into subsets U and V , and the edges are elements of U × V and V × V only. In particular, no edges start in V and end in U . We define a directed simple path in a directed graph as a sequence of distinct nodes {i 1 , i 2 , . . . , i k } such that (i l , i l+1 ) ∈ E for all l = 1, . . . , k − 1. All paths considered in this paper will be directed simple paths. For a path P = {i 1 , . . . , i k } we define the length of P , denoted |P |, to be the number of nodes in the path, namely k. Our goal is to find a collection P = {P 1 , . . . , P k } of node-disjoint paths with each path P j starting with a node in U and each having length at most K. Specifically, we define v(P) = P ∈P |P | to be the count of the nodes in a collection of paths. Let Π be the set of feasible collections P. Then our problem is We now introduce a formalization for the problem (1) that is conducive to message passing. We will use special symbols * and • to represent certain relationships to be described below. We will also introduce the notation ∂i to represent the set of nodes adjacent to the node i in the undirected version of the graph G. That is, ∂i is the set of nodes j ∈ [n] such that (i, j) ∈ E or (j, i) ∈ E. We also let ∂ X i = ∂ i ∩ X and let ∂i \ j stand for the more formally correct ∂i \ {j}.
For each node i ∈ [n], we introduce variables These variables intend to represent the depth of the node (d i ), its parent (p i ), and its child (c i ) in a given solution. The special value * indicates that the node does not participate in a path, in which case we will have d i = p i = c i = * . The special value • indicates that the node is at one end of a path. Thus if p i = • this means node i does not have a parent and thus starts a path, whereas if c i = • node i does not have a child and thus ends a path. Note that if c i is a node it must be a node in V because nodes in U must either act as the root of a path or not participate in a path.
To formulate the optimization problem, we need to define a feasible set and an objective function. We begin with the feasible set. We first define a function that acts as an indicator for the consistency of the variables (d i , p i , c i ) at a single node i.
Term (2) is the case in which node i does not participate, term (3) is the case in which node i is the root of a path, and term (4) is all cases where node i participates in a path but is not the root. Note that these three cases are mutually exclusive, so f i is an indicator function. It is equal to zero if the variables at node i are not consistent with any global configuration (d, p, c) = (d j , p j , c j ) j∈[n] representing a set of valid paths for the graph G, and equal to one otherwise.
For each pair of nodes i and j that are adjacent in the undirected version of G, we define a function g ij (d i , p i , c i , d j , p j , c j ) of the variables (d i , p i , c i , d j , p j , c j ) for the nodes i and j which is an indicator that these variables are consistent with each other and the existence or nonexistence of the edges (i, j) and (j, i). Namely, they correspond to a locally valid configuration, so if, e.g., node i reports that its parent is node j, then the edge (j, i) exists, node j reports that its child is node i, and we have Lines (5) and (6) capture the cases where j is the parent of i and i is the parent of j, respectively, and enforce the relationships that must occur in those cases. Line (7) captures all cases where there is no direct connection between nodes i and j. All three cases are mutually disjoint, so g ij is an indicator function. In summary, g ij is equal to one if parent/child relationships agree, the edges (i, j) and (j, i) are used only in the appropriate directions, and depth relationships are consistent.
Finally we define a indicator function for variable consistency at both the node and edge level, which we define for each i ∈ [n], j ∈ ∂i as (8) With this definition of h ij , we can define the set which will serve as the feasible set for our optimization. Note that each collection P ∈ Π of node disjoint paths is associated with a unique assignment of variables {(d i , p i , c i )} i∈ [n] , which we will denote by (d(P), p(P), c(P)). In fact the converse is also true: Proposition 1. For every (d, p, c) ∈ M there exists a unique P ∈ Π such that (d, p, c) = (d(P), p(P), c(P)).
Proof. Given a set of variables (d, p, c) ∈ M , recall that for every i ∈ [n] we have f i (d i , p i , c i ) = 1 and for every j ∈ ∂i we have g ij (d i , p i , c i , d j , p j , c j ) = 1.
We will construct a feasible set of node-disjoint paths by the following procedure. Find the set of nodes i such that d i = 1, which we denote R.
If R is empty then for every i ∈ U we must be in case (2) because case (3) is excluded. (2) and (4) we either have d i = * or d i > 1. But by (5) d i = 2 only if d j = 1 for some j ∈ ∂i. This implies that d i = 2 for any i ∈ V . Repeating this argument implies d i = l for any l ∈ [K], so we conclude d i = * . Then (2) implies p i = c i = * . Thus for every node i ∈ [n] we have Otherwise, if R is nonempty, we will construct one path for each i ∈ R. Consider some i 1 ∈ [n] such that d i1 = 1. By (3) we must have i 1 ∈ U , so it will be valid to start a path at i 1 . Again by (3) we know c i1 = * , •, so c i1 = i 2 for some i 2 ∈ ∂ V i. By (6) with i = i i and j = i 2 we have d i2 = d i1 + 1 = 2, (i 1 , i 2 ) ∈ E, p i2 = i 1 and c i2 = i 1 . Thus it is valid to start our path with the node sequence {i 1 , i 2 }. If c i2 = • then we end this path and add P = {i 1 , i 2 } to P. If c i2 = • then by d i2 = 2 and (4) we have c i2 = i 3 for some i 3 ∈ ∂ V i 2 . By (6) with i = i 2 and j = i 3 we have d i3 = 3, (i 2 , i 3 ) ∈ E, p i3 = i 2 , and c i3 = i 2 . Thus it is valid to start our path with the node sequence {i 1 , i 2 , i 3 }. If c i3 = • we end this path and add it to P. Otherwise we iterate, because (4) implies c i3 = i 4 for some i 4 ∈ ∂ V i 3 . In the generic step, we are considering c i l−1 = i l for some i l ∈ ∂ V i l−1 and we have d i l−1 = l − 1. Then (6) for i = i l−1 and j = i l implies d i l = l, (i l−1 , i l ) ∈ E, p i l = i l−1 and c i l = i l−1 . If c i l = • we terminate the path and add P = {i 1 , i 2 , . . . , i l } to P. Otherwise there exists some i l+1 ∈ ∂ V i l such that c i l = i l+1 so we iterate. This process is guaranteed to terminate because d i l = l and d i l ≤ K. Once we have added the path starting at node i 1 we repeat the process for each node i with d i = 1.
To see that (d(P), p(P), c(P)) = (d, p, c) after this process, first consider a node j which does not participate in any path in P. We have d j (P) = p j (P) = c j (P) = * so we want to show d j = p j = c j = * . Because every node i with d i = 1 starts a path in the above procedure, we must have d j = 1. Suppose d j = l for some l ∈ {2, . . . , K}. Then (4) implies p j = { * , •} so there exists some i ∈ ∂j such that p j = i and (6) implies d i = l − 1. If i participates in a path in P then our above procedure would also include j. If i does not participate then we have found a non-participating node with depth l − 1. Thus if there are were any non-participating j with d j ∈ [K] there would be a nonparticipating j with d j = 1, but we have already shown there is no such j. We conclude that d j = * for all non-participating d j . Then (2) implies p j = c j = * , as desired.
Next consider a node i l which participates in a path in P at depth l. If l = 1 then i l starts a path and has no parent and (by the above procedure) has child c i l . Thus we have d i l (P) = 1, p i l (P) = •, and c i l (P) = c i l . The only nodes i l that start paths are those with d i l = 1 and (3) implies p i l = •, so we have (d i l (P), p i l (P), c i l (P)) = (d i l , p i l , c i l ), as desired. If l > 1 then i l has d l = l, parent i l−1 = p i l (as noted during the procedure above) and either ends the path or has a child i l+1 = c i l . By the procedure the path ends only if c i l = •, so we have c i l (P) = • = c i l . Thus we have (d i l (P), p i l (P), c i l (P)) = (d i l , p i l , c i l ), as desired. If the path continues it continues to i l+1 = c i l so c i l (P) = c i l , and again we have (d i l (P), p i l (P), c i l (P)) = (d i l , p i l , c i l ), as desired.
Because we have a bijection between feasible collections of nodedisjoint paths in G and the set M , we can treat M as the feasible region for our optimization problem.
To define our objective, we define the function which assigns value 0 to the case where the node i does not participate in a path and value 1 to all other cases. We further define the function 3 Algorithm.
Belief propagation algorithms are used to estimate either the marginal distributions of variables or the solution to the maximum a posteriori (MAP) problem in graphical models. We will use the framework corresponding to the latter option. The MAP solution for a probability distribution is the most likely joint assignment of all the variables in the distribution. To utilize this approach, we transform the optimization problem (12) into a problem of finding a MAP assignment for a probability distribution. To achieve this transformation we introduce a corresponding undirected graphical model. A maximum likelihood assignment of all variables in this distribution will be an optimal solution to (12).
To create our undirected graphical model, we introduce node potentials where β > 0 is a system parameter to be chosen. We use h ij as defined in (8) as our edge potentials. Together these define the probability distribution P This probability distribution assigns zero probability to any global configuration (d, p, c) which does not correspond to an element of M , and among positive probability configurations it assigns higher probability to those which include more nodes.
We can now solve our original optimization problem by solving the MAP problem on this undirected graphical model. Note that we may have multiple solutions.
The core component of the algorithm is the set of messages, which we define for each ordered pair (j, i) with j ∈ [n] and i ∈ ∂j, representing the message from node j to node i. Thus each undirected edge (i, j) ∈ E ′ is associated with two sets of messages, representing the information sent from i to j and from j to i, respectively. The Our approach to this problem will be to implement parallel belief propagation in the standard Min-Sum form (which can be derived from the Max-Product form by taking the negative log of all messages; see, e.g., [6] or [8]), with messages sent according to the update rule We use the convention − log(0) = ∞. In other words, if h ij = 0 for some configuration then that configuration will not be included in the minimization. Otherwise the log(h ij ) term contributes log(1) = 0. Because our graphical model potentially has cycles, we must compute these messages iteratively from some starting values, which we set to be b j→i (d i , p i , c i ) = 1 for all j ∈ [n], i ∈ ∂j and all values of d i , p i , c i , with the goal of finding a fixed point solution. From the messages we can compute the max-marginals at each node as These max-marginals are, according to the Belief Propagation approach, intended to approximate the probability of a maximum-probability configuration of the whole system given the arguments d i , p i , c i are fixed. As such, they can be used in reconstructing the maximum global configuration. The non-uniqueness of solutions to (12) complicates this, however, which we discuss below in section 5. The update rule (14) and max-marginal equation (15) characterize the entirety of our algorithm. For some a priori maximum number of iterations T and initial messages b 0 If the messages reach a fixed point, i.e. b t = b t+1 for some t < T , we halt iterations at that point. Furthermore, for any 0 < t ≤ T we can use b t and (15) to compute estimates of the max-marginals for any node in any configuration at that step of the algorithm.
In this form, the number of messages is O K i∈[n] ∆ 3 i where ∆ i = |∂i| is the degree of node i in the undirected graph. Indeed, for each node i ∈ [n] there are ∆ i incoming messages, and for each message there are K choices of d i and O ∆ 2 i choices of (p i , c i ) that define messages. For ∆ = max i∈[n] ∆ i the number of messages is At each iteration of the algorithm a message from j to i must be updated by computing a minimization over all O K∆ 2 j choices of the variables (d j , p j , c j ). The argument of the minimization includes a sum over the neighbors of j, which can be performed once for each j ∈ [n] for a total computational cost O (n) per iteration. Thus a full iteration of all messages requires O nK 2 Rather than implementing this scheme directly we will instead introduce a more efficient representation of the messages, which is described in the next section.

Efficient representation.
In this section we introduce a more efficient representation of the message passing algorithm described above. This representation will require only O (nK∆) messages compared to O nK∆ 3 for the original description (14), and will require O (nK∆) computations to be performed per iteration compared to O n 2 K 2 ∆ 6 computations for the original representation. Before defining the messages for the new representation, we now introduce an intermediate set of messages and rewrite (14) as two steps, with If the set over which we are minimizing in (18) is empty we define the minimum to be ∞. The messages b j→i (d i , p i , c i ) evolve exactly as they do in (14). We will define a collection of messages in terms of ψ and b, intended as an efficient representation of the message passing algorithm, which we call the A-H form. After defining the messages we will provide the rules by which these messages are updated. We will also obtain algorithmic complexity bounds for our representation, which are summarized in Theorem 1. Then we will demonstrate in Theorem 2 the equivalence of the A-H form algorithm to the original description (14) by showing that iterating the A-H form allows us to compute the maxmarginals (15) as if we had done the iterations in the original form (14). This theorem will be proved in two parts, with Lemma 1 establishing a way to write all original messages b j→i (d i , p i , c i ) in terms of A-H form messages and Lemma 2 establishing how to update the A-H form messages.

Message definitions.
The messages to be defined fall into three categories based on the type of the node which sends the message, and the type of the node which receives it. We first define messages sent from a non-root node j ∈ V to another non-root node i ∈ V . Next we will describe messages sent from a non-root node j ∈ V to a root node u ∈ U , and finally we describe messages sent from a root node u ∈ U to a non-root node Fix (j, u) ∈ E ′ such that j ∈ V and u ∈ U .
• Let Fix (u, i) ∈ E ′ such that u ∈ U and i ∈ V .
• Let Note that when i ∈ U or j ∈ U there are different numbers of messages in each direction. To refer to the set of messages from one node to another we define for Note that in all cases the number of messages from i to j is O (K).

Update rules and the complexity of the A-H form.
The iteration update rules for these new messages are listed in Tables 1-3. The derivation of these update rules from the definitions of the messages is given in the proof of Theorem 2. These update rules allow us to provide complexity bounds for the system: In performing an iteration of the A-H form we proceed node-bynode. As we will show below, at a node j we will perform O (K∆ j ) total computations. This implies that we require O K j∈[n] ∆ j = O (n∆K) total computations for a full iteration of the A-H form, as claimed.
We now provide more detail on the computations required per node j. We consider the case j ∈ V ; if j ∈ U the analysis is similar. We will also focus our attention on the most complex message to compute, Table 1: Table 2: Messages from j ∈ V to u ∈ U , for (u, j) ∈ E ′ and 3 ≤ d ≤ K − 1 (26) Rather than compute these directly for each i, we first compute select the three neighbors l ∈ ∂ V j with the smallest values of and select the three neighbors l ∈ ∂ V j with the smallest values of We also record the values of (28) and (29) for the nodes we record. Once we have these quantities, for each i ∈ ∂j we can compute

Initialization and termination.
The equations in Tables 1-3 collectively define update rules for our algorithm. To initialize the system we set X 0 j→i = 1 for all (i, j) ∈ E ′ , by which we mean we set each element of X 0 j→i to 1. Given some a priori bound T on the number of iterations, we then compute X t from X t−1 for t ∈ [T ] according to the update rules in Tables 1-3. If the messages reach a fixed point with X t = X t+1 for some t < T we halt the algorithm.

Validity of the A-H Form.
Identities relating original form messages b j→i (d i , p i , c i ) to A-H form messages for various choices of d i , p i , and c i are given in Tables 4-6. It can be checked that for any choice of (d i , p i , c i ) that does not appear in a table we have b j→i (d i , p i , c i ) = ∞ by (18) because no choice of (d j , p j , c j ) can satisfy h ij (d i , p i , c i , d j , p j , c j ) = 1. The validity of these identities given the definitions of the A-H form variables is proved in Lemma 1 below.
To use the original form (14) to find solutions to (12) we need to compute the max-marginals (15) for different values of (d i , p i , c i ) from the messages b k→j (d j , p j , c j ). Specifically, we need to find the choices of (d i , p i , c i ) that maximizep i (d i , p i , c i ). Thus to show that we can iterate the A-H form to find the same approximate solutions we would find by iterating the original form (14) we prove that the max-marginal can be computed from the A-H form. In particular, we define We prove the following: when both sides are computed with appropriate messages b t j→i and X t j→i updated t times according to (14) and Tables 1-3, respectively. We will use two lemmas to prove this theorem: Lemma 1. The identities described in Tables 4-6 hold.
Proof. Fix i, j ∈ V . We focus on two representative messages, namely A d j→i for 3 ≤ d ≤ K − 1 and H j→i . In the first example we will prove the fifth identity of Table 4 and in the second will prove the sixth. The other identities described in Table 4 and all identities described in Tables 5-6 are proved similarly so we omit the details.
We will first show that the fifth identity of Table 4 holds. Namely,  Table 5: Classes of b j→u messages for j ∈ V , u ∈ U .
Suppose (i, j) ∈ E. Then h ij (d, p, j, d j , p j , c j ) = 0 for any choice of d j , p j , c j and thus b j→i (d, p, j) = ∞. By (20) we have A d+1 j→i = ∞, so b j→i (d, p, c) = A d+1 j→i , as desired. Now suppose (i, j) ∈ E. We now use h ij (d, p, j, d j , p j , c j ) = 1 to infer information about d j , p j , c j . From (6) we must have d j = d + 1, p j = i, and c j = i. Then (4) further implies c j = * . Substituting these facts into (31) shows Observe that because δ i→j = 0 the right hand is the definition of A d+1 j→i as given in (20), so we conclude b j→i (d, p, j) = A d+1 j→i , as desired. Now we will show the sixth identity of Table 4 holds. Namely, for We now want to break the minimization into cases based on how we satisfy h ij (d, p, c, d j , p j , c j ) = 1. First observe that p = j and c = j implies we must be in case (7), which implies p j = i and c j = i. Then there are two possible cases for how to satisfy f j (d j , p j , c j ) = 1.
First we can choose (2) and let d j = p j = c j = * . Otherwise we are in (4) so we have d j = d ′ for some 2 ≤ d ≤ K, p j ∈ { * , •} and c j ∈ { * , p j }. Combining these cases implies Observe that the right hand side contains the definitions of G j→i from (23) and F d ′ j→i from (22), so we have b j→i (d, p, c) = min G j→i , min as desired, where the last step is applying the definition of H j→i from (24). Tables 1-3

for updating A-H form messages are valid.
Proof. We will show how to derive one representative update equation, namely for 3 ≤ d ≤ K −1, which is the first equation of Table 1. The remainder of Table 1 is derived via similar methods, as are the relations in Tables  2 and 3.
Recall (21): We can split this minimization into two cases, namely c j = • and c j ∈ ∂j \ i, which allows us to write We now replace each instance of b k→j with the appropriate value according to the identities in Tables 4 and 6. First note that for k ∈ ∂ U j \ i we have b k→j (d, i, k) = ∞ because d ≥ 3 and k ∈ U implies h jk (d, i, k, d k , p k , c k ) = 0 since g jk (d, i, k, d k , p k , c k ) = 1 if and only if d k = d + 1 and if d k = d + 1 > 1 then f k (d k , p k , c k ) = 0. Thus minimizing over k ∈ ∂j \ i is equivalent to minimizing over k ∈ ∂ V j \ i so we can write We now replace these particular messages via b k→j (d, i, k) = A d+1 k→j from (30) (which is the fifth line of Table 4) and b l→j (d, i, k) = H l→j and b k→j (d, i, •) = H k→j from (32) (which is the sixth line of Table  4), so we have where the second line simply extracts the sum from the minimization. This equation now matches (33).
We are now prepared to prove Theorem 2.
Proof of Theorem 2. Tables 4-6 directly provide the functions ρ i required to compute max-marginals. This substitution is valid after any number of iterations because it is valid from the original definitions of the A-H form variables by Lemma 1 and Lemma 2 proves the A-H variables are iterated by rules that follow directly from the iteration rules for the original form variables.

Constructing solutions.
Because we have no guarantee that messages will converge, or that if they converge they will converge to an optimal solution, we elect to compute a feasible solution at every iteration of the algorithm and return the best solution found over all iterations. In particular we run our algorithm for a fixed number of iterations unless we converge to a fixed point before this limit. Finding all these solutions comes at an additional computational cost, but it gives consistently better solutions than simply using the result at the end of the iterations. To initialize our solution construction algorithm, let A = V be the set of non-root nodes which we have not yet assigned to a path. Nodes will be removed from this set as we determine their role in our constructed solution. For each node u ∈ U , we will construct a (possibly empty) path P u , and our full solution will be the collection of these paths. The first step is to compute the max-marginals for the possible choices at node u, which are either "do not participate in a path" (in which case P u = ∅) or "start a path by choosing a child from among available neighbors of u". We want to maximize the maxmarginal among these choices. Since our variables are in Min-Sum form we actually compute the negative log of the max-marginals and choose the minimizing configuration. To choose a configuration at the node u, we use the functions ρ i from Definition 1 to compute and find the minimum value. If the minimum is (34), we let P u = ∅ and proceed to the next root node. If the minimum is (35) for some neighbor j, we initialize our new path as P u = {u, j} and remove j from the set of available nodes by setting A = A \ {j}. Note that (35) considers only neighbors which are also elements of A insuring that P u does not include any nodes used in a previous path. Further note that we must have (u, j) ∈ E so this is a valid start to a path because j ∈ ∂ V u and E contains no edges in V × U . We next proceed to the node j to determine whether P u should be extended or added to the solution as is. This process is essentially the same as that for u, involving choosing the optimal value of the max-marginal at j, but we are now constrained to configurations with parent u. This means that we are choosing whether to end the path at j or continue to one of j's neighbors in V , so we compute and choose the minimum. If the minimum is (36) the path is terminated so we add P u to our solution as is and proceed to the next root node. Otherwise the minimum is achieved by a particular choice i, so we set P u = {u, j, i}, remove i from the set of available nodes by setting A = A \ {i}, and proceed to determine whether P u should be extended or added to the solution as is. Note again that i ∈ A guarantees we have not previously selected i to participate in any other path.
To confirm that (j, i) ∈ E, note that the first line of Table 1 implies (37) is infinite for that choice of i. Such an i cannot be the chosen minimum value because (36) is always finite as can be seen from the last two lines of Table 1 and the  first line of Table 6. Thus if we choose i to continue the path we must have (j, i) ∈ E which guarantees this is a valid choice.
After the first two nodes of a path the process becomes more generic, so let us suppose we have some P u such that |P u | = d ≥ 3 and the last two nodes of the path P u are are p and j. We still need to determine whether to extend the path P u or add it to the solution as is. If d = K then we know that we cannot extend the path, so we add it to the solution and move on to the next root node. Otherwise, we compute and choose the minimum. If the minimum is (38) we leave P u as is, add it to the solution and move on to the next root node. If the minimum is (39) for some i, then we append i to the end of P u and remove i from A. We now have a path of length d + 1 and can repeat the above procedure until we terminate the path. As above, we restrict our attention to nodes in A to insure no node is assigned to multiple paths. Furthermore, for i ∈ ∂ V j ∩ A if (j, i) ∈ E then A d+1 i→j = ∞ but (38) is always finite because H k→j is finite for any pair of nodes k and j and B d−1 p→j is finite because we are guaranteed to have (p, j) ∈ E by the fact that these are the last two nodes of the partially constructed path P u . Thus if we choose to continue the path we do so along an edge with the correct direction to a node that is not participating in another already constructed path, and so the end result of this process is a feasible collection of node-disjoint directed paths.
Notice, however, that the result depends on the order in which we process roots, so in practice we often repeat the procedure several times with different orders and choose the best solution among those we find.

Alternative algorithms.
Our numerical results in section 7 below will compare our efficient BP algorithm to three alternatives: a Greedy algorithm, an IP based algorithm for KEP from [2], and a novel IP based algorithm. We now describe these algorithms.
We propose the Greedy algorithm as a simple to implement approximation algorithm with good scaling behavior. It is greedy in the sense of adding paths to the solution one at a time and adding the longest possible path at each step, but when adding a path the search for the longest path is exhaustive. Given an arbitrary ordering of root notes, the Greedy algorithm proceeds as follows. Let A = V be the set of non-root nodes that have not yet been added to a path. From a root node u ∈ U , all possible directed paths with length at most K consisting of u and a collection of nodes in A are explored and the longest such path is chosen, with ties being broken arbitrarily. The longest possible path is found by recursively looking for the longest path from each neighbor, keeping track of the depth bound to insure a path longer than K is not chosen. The chosen path is then added to the solution and its nodes are removed from A. If the path consists only of the node u it is discarded. The algorithm now proceeds to the next root node in the given order.
This algorithm is guaranteed to generate a feasible collection of directed paths, but the size of that collection is highly dependent on the order in which we consider root nodes, so we repeat the process for a large number of different orderings of the root nodes, choosing the best solution. For the numerical results to follow we often used about 200 randomly generated orders. In most cases, the number of orderings is chosen so the running time of the Greedy algorithm is roughly equal to the message passing algorithm (which also includes a similar randomization step, as noted in section 5).
The KEP algorithm from [2] is based on a Traveling Salesman IP formulation approach to the kidney exchange problem and can find bounded length cycles as well as paths, which may be bounded or unbounded in length. For our comparisons we enforce our path length bound K and disallow all cycles.
Finally we have implemented a simple IP formulation of the problem based on the same description we used in the development of the message passing algorithm. We call this the Parent-Child-Depth (PCD) formulation. We let The formulation is as follows: There is a bijection between feasible solutions to the IP problem (40)-(53) and the set M defined in (9) and an optimal solution to the IP problem maximizes the objective function H(d, p, c) defined in (11).
Remark. This proposition states that the IP problem (40)-(53) is in fact the same problem we deal with elsewhere in this paper. Before providing a formal proof, we provide some intuition for the formulation.
The variable x i acts as an indicator for whether the node i participates in a path, with the various p ji and p •i indicators representing each possible choice of parent for node i (including the possibility of starting a path). Similarly the c ij variables represent the possible choices of child for the node i. The integer variable d i represents the depth of i in the path. This depth variable is decomposed into variables d ji and d •i which will be equal to d i if the parent of i is j or i starts a path, respectively, and zero otherwise. Thus at most one of d •i and the d ji will be nonzero, representing the actual choice of parent implied by the values of p •i and the p ji . The various constraints serve to guarantee that this interpretation of these variables is accurate and all feasibility constraints for the feasible set M .
More specifically, constraint (41) insures that a node participates if and only if it has a parent or starts a path. The remaining constraints insure that each node has at most one parent and one child, and that all variables agree locally with each other. Note (43), which requires a root to have a child (so each path contains at least one edge). The depth bound is enforced by requiring d i ≤ K and insuring that d i accurately represents the depth of i. Note in particular (49), which insures that depth increases from parent to child. This constraint is quadratic.
Proof of Proposition 2. Given (d ′ , p ′ , c ′ ) ∈ M , we can construct IP variables (x, d, p, c) node-by-node. For each i ∈ [n] there are three possible cases, represented by (2) then the node i does not participate in a path, so we let x i = p •i = d i• = 0, for j ∈ ∂ − i we let p ji = d ji = 0, and for j ∈ ∂ + i we let c ij = 0.
• If p ′ i = •, then the node i starts a path, so we let •} then the node i has a parent p ′ i ∈ ∂ − i, so we let x i = p p ′ i i = 1 and let p •i = 0 and for j ∈ ∂ − i \ p ′ i let p ji = 0. By (4) we have 2 ≤ d ′ i ≤ K, so we let d i = d p ′ i i = d ′ i and let d •i = 0 and for j ∈ ∂ − i \ p ′ i let d ji = 0. If c ′ i = • we let c ij = 0 for all j ∈ ∂ + i because i has no child. Otherwise we let c ic ′ i = 1 and for all j ∈ ∂ + i \ c ′ i let c ij = 0. With variables constructed in this way it is straightforward to check that all of the constraints are satisfied, so we omit the details. Now suppose we have a feasible set of IP variables (x, d, p, c). We will construct (d ′ , p ′ , c ′ ) as described in Section 2 and show that it is an element of the set M . We proceed node-by-node. For i ∈ [n]: For i ∈ [n] the three cases above correspond exactly to the cases (2)-(4) to guarantee f i (d ′ i , p ′ i , c ′ i ) = 1. Note in particular that (47) For i ∈ [n], j ∈ ∂i, we have three cases. If p ′ i = j then we know p ji = 1 and j ∈ ∂ − i, so (48) implies c ji = 1 and thus c ′ j = 1. (50) guarantees p ′ j = i and c ′ i = j. j ∈ ∂ − i guarantees (j, i) ∈ E and (49) Thus we are always in one of the cases (5)-(7), so To see that the optimization is the same, we need only note that x i = 1 if and only if p ′ i = * . Thus by (10) we have . In other words, both are simply measuring the number of nodes that participate in a path.
Note that this formulation is polynomial in n, both in the number of variables and in the number of constraints, and is therefore simple to implement. As we will see below, for a number of instances, with n = 1000 and K = 15, PCD produces optimal solutions in more than 50% of tested graphs (see Table 9), while KEP does not find optimal solutions in the same allowed time.

Numerical results.
We compare the performance of our algorithm, which we call BP, to the three alternative algorithms described in section 6. Note that the two IP-based algorithms (KEP and PCD) allow us to solve the problem to optimality given sufficient time, whereas the Greedy and BP algorithms have no optimality guarantees.
When running the BP algorithm we use the value β = 0.01 for the parameter defined in (13). Informal testing of different values of β indicated that this value had, in some cases, a small advantage of about 1% in the size of solutions over larger β values while smaller β values resulted in significantly worse solutions. When this advantage was not present the value of β had no noticeable effect on the solution size.
For all cases except the three largest real networks we use 200 orders of root nodes for the Greedy algorithm and 5 orders of root nodes for each iteration of BP. For the largest graphs these are decreased to 50 and 2, respectively. In all instances BP is run for a maximum of T = 50 iterations, terminating earlier only if it has converged to a fixed point or exceeded a preset time limit. In most cases BP did not  Table 7: Random graphs on n = 1000 nodes with maximum path length K = 5 and edge probability c/n. 100 samples were run for each row, with β = 0.01 for BP. Dark shaded KEP results are cases where KEP was optimal in at least 99% of samples. BP results are labeled by finding more nodes than Greedy and finding the most nodes of any method . consistently converge, though it did converge somewhat more often for smaller graphs with n = 1000 and for longer allowed paths with K = 15. We have not seen a significant correlation between instances where BP converges and those where it finds significantly better solutions, but convergence does improve running time simply by allowing the algorithm to terminate earlier.

Random graph comparisons.
We report our results on random graphs constructed as follows. We let each ordered pair of nodes in V × V and U × V have an edge independently with probability p. We set p = c/n and vary c. Thus for every pair (i, j) ∈ U ∪ V × V the directed edge (i, j) is present with probability c/n and absent with probability 1 − c/n. Table 1 shows results for some varieties of random graphs. The column "Root%" indicates the percentage of nodes which are in the set U , namely |U |/n. The column "c" indicates the value of the parameter c for the graphs tested. The remaining columns summarize the results of testing the algorithms on 100 random graphs, first showing the average number of nodes in the returned solution and then showing the average running time in seconds. For the IP based algorithms a time limit of 60 seconds was implemented by the IP solver and if the optimal solution was not found within that time the best feasible solution identified by the solver was used. While in many cases KEP returned a nontrivial feasible solution, if PCD did not find the optimal solution it almost always returned the trivial feasible solution zero. For the cases represented in Table 7 PCD never returned a nontrivial feasible solution so we omit PCD from the table. It does, however, appear in other cases below and performs well in some ranges outside those that appear in tables, as will be discussed later.
Testing over a wide range of parameter values we find that in many instances one or both of the IP algorithms runs to termination in time comparable to BP or Greedy. In these cases Greedy and BP both generally return suboptimal solutions, with both algorithms having ranges of parameters where they enjoy a small advantage over the other. This advantage ranges from Greedy finding about 5% more nodes than BP to BP finding around 10% more nodes than Greedy. Table 7 shows results for random graphs with 1000 nodes and a variety of percentages of root nodes. The algorithms were run with path length bound K = 5 and edge probability c/n for c = 2, 3, 4. For c = 2 KEP is optimal in every sample and has a running time comparable to BP. BP finds feasible solutions with 93%-98% of the nodes of the optimal solution, whereas the Greedy solutions have 85%-95%. As we increase c, KEP finds fewer optimal solutions, though it is still optimal in all but 1 of our 100 samples when 10% of the nodes of the graph are roots for c = 3. In most cases BP finds more nodes than Greedy does while falling short of even the suboptimal KEP solutions, but for c = 4 with 20% and 25% root nodes, BP performs the best among the compared algorithms, finding the most nodes in time comparable to Greedy and significantly shorter than KEP.
While we do not show full results here, we can provide some sense of the behavior of the algorithms for graphs just outside the regimes in Table 7. If the root percentage is decreased below 10% KEP is generally optimal and both BP and Greedy often find this optimal solution, likely because the scarcity of roots leads to less interaction between paths in feasible solutions. If the root percentage is increased above 25% BP still finds more nodes than Greedy but significantly fewer than KEP, which generally finds optimal solutions. Though PCD did not find nontrivial solutions for the cases in Table 7, it is often optimal for these cases with many roots. If the root percentage is kept in the 10%-25% range but c is decreased to 1, both KEP and PCD are always optimal, with Greedy missing a small percentage of nodes and BP slightly behind Greedy. If c is increased above 4, both Greedy and BP find significantly (10%-20%) more nodes than KEP, which is never optimal, and Greedy usually has a small (1%-3%) advantage over BP. In this regime PCD finds no nontrivial solutions. Table 8 shows tests over some of the same parameters as  Table 8: Random graphs on n = 10000 nodes with 20% root nodes, maximum path length K = 5, and edge probability c/n. 100 samples were run for each row, with β = 0.01 for BP. Dark shaded KEP results were optimal in all samples. BP results are labeled by finding more nodes than Greedy and finding the most nodes of any method .
limit for KEP and PCD is increased to 20 minutes, and again PCD is omitted because it never finds nontrivial solutions within the time limit. Most general trends remain the same, with KEP always finding optimal solutions within the time limit for c = 2 but BP performing better for c = 3 and c = 4. Specifically, BP now has an advantage over KEP for c = 3. Another factor to note is that BP has a significant advantage over Greedy in running time in the c = 4 case, which we can credit to the fact that BP converged in all samples for that case, running an average of 19 iterations as compared to 50 for the other cases.
In Table 9 we provide summary results for problems with longer allowed paths; that is with K = 10 and K = 15. Here KEP is almost never optimal, though it does often provide a feasible solution with more nodes than BP or Greedy, albeit at the cost of running to the cutoff time rather than terminating in a few seconds. In all cases for this regime BP finds more solutions than Greedy while running in comparable time, and in two of these cases it also on average finds more nodes than any other tested algorithm. In some cases the solution size advantage over Greedy is significant, with BP solutions more than 10% larger than Greedy solutions. Note that PCD performs rather well for K = 15, often finding the optimal solution, though the average size of the solution is reduced by the fact that PCD rarely returns a nontrivial feasible solution, generally either returning the optimal solution or an empty solution. PCD found the optimal solution in at least 50% of samples for all the K = 15 cases, while KEP only found any for c = 2, for which it was optimal in only 6% of cases. Thus there is an advantage to PCD over KEP for all K = 15 cases as it much more regularly finds the optimal solution. When it does find the optimal solution it is also relatively fast, with an average runtime of 23 seconds for cases across all 4 values of c in which PCD found the optimal solution. BP and  Table 9: Random graphs on n = 1000 nodes with 15% root nodes, maximum path length K, and edge probability c/n. 100 samples were run for each row, with β = 0.01 for BP. Shaded PCD results were optimal in 96% of samples, and in the other K = 15 cases PCD was optimal in at least 50% of samples. BP results are labeled by finding more nodes than Greedy and finding the most nodes of any method .
Greedy both remain significantly faster, and BP does on average find better solutions than either PCD or Greedy. For similar graphs just outside the regime presented in Table 9 we have a couple situations. As in the case of a shorter path bound, decreasing c to 1 makes KEP and PCD optimal in all cases, with Greedy missing a small percentage of nodes and BP slightly behind Greedy. Increasing c above 5 leads to Greedy and BP having essentially identical performance, both generally better than KEP and PCD. PCD still finds optimal solutions in a reasonable percentage of K = 15 cases, however. If the path length bound is kept as long as 15 or increased up to 25, PCD performs very well over a wide variety of c values and root percentages from 15% to 30%, generally finding the optimal solution. In many of these cases KEP also finds an optimal solution but does not terminate and confirm that the solution is optimal so PCD has a significant advantage in running time in those cases. BP and Greedy are suboptimal in these high-K cases, both finding 80%-90% of the nodes in the optimal solution.

Computational results for real world networks.
We compared the performance of KEP, PCD, Greedy, and BP algorithms on several real world networks taken from the Stanford Large Network Dataset collection [7]. We chose directed graphs and randomly selected a subset of nodes to act as roots (discarding directed edges pointing in to roots and discarding isolated nodes). We present performance comparisons for five networks, which are as follows: • Epinions: this network is drawn from the website Epinions.com and represents "trust" relationships between users, in which each user can indicate whether they trust the opinions of each other user. Thus an edge from i to j indicates that user i trusts user j.
• Gnutella: a snapshot of the Gnutella peer-to-peer file sharing network, with nodes representing hosts and edges representing the connections between those hosts.
• Slashdot: users of the technology news site Slashdot were allowed to tag other users as friends or foes; our network has users as nodes and tags as directed edges, drawn from a snapshot in 2008.
• Wiki-vote: A record of votes on administration positions for Wikipedia, in which administrators vote on the promotion of other potential administrators. Nodes are wikipedia users and a directed edge from node i to node j represents that user i voted on user j.
• Amazon: product co-purchasing on Amazon.com, based on the "customers who bought this item also bought feature". A directed edge from product i to product j indicates that j appears on the page for i as a frequently co-purchased item. The data is a snapshot from March 2003.
The results of our tests on the full graphs are summarized in Table  10. The "Nodes" column indicates the size of the given graph, which is given as an average over the five tests because the selection of root nodes results in a different number of nodes becoming isolated and being dropped in each test. The column "∆" is the average degree of the graph, which is also averaged across the choices of root nodes for each trial. Remaining columns show the number of nodes in the average solution and the running time in minutes for each algorithm. With the exception of the graph Gnutella, KEP in all cases crashed rather than running to the stopping time and so provided no solution, which we indicate by "-". We believe these crashes resulted from KEP exceeding the available memory. They also caused the running times for KEP to be erratic, with the algorithm crashing before the time limit in most cases. The extreme running time for Amazon, the largest graph, can be explained by the external software we used failing to adhere to the requested time limits before crashing. In all cases PCD failed to return a nontrivial solution, so it is omitted from the table. For all graphs except Slashdot, BP enjoys an advantage in solution size over Greedy, finding 4%-15% more nodes while running only slightly longer on average. On Slashdot Greedy has an advantage in both solution size and running time. We also note that even on the graph Gnutella,  Table 10: Real world networks with 20% root nodes (randomly selected) and maximum path length K = 5. We use β = 0.01 for BP. ∆ is the average degree of the sampled graph. 5 samples were run for each graph, with − representing graphs for which KEP crashed without returning any solution. Shaded cells represent the best solutions, with lighter shading for the case when BP provides a better solution than Greedy.
the one case where KEP returned feasible solutions, BP has a huge advantage in running time for a solution with 90% of the nodes of that returned by KEP.