An analytical solution to the multicommodity network flow problem with weighted random routing

We derive an analytical expression for the mean load at each node of an arbitrary undirected graph for the non-uniform multicommodity flow problem under weighted random routing. We show the mean load at each node, net of its demand and normalized by its (weighted) degree, is a constant equal to the trace of the product of two matrices: the Laplacian of the demand matrix and the generalized inverse of the graph Laplacian. For the case of uniform demand, this constant reduces to the sum of the inverses of the non-zero eigenvalues of the graph Laplacian. We note that such a closed-form expression for the network capacity for the general multicommodity network flow problem has not been reported before, and even though (weighted) random routing is not a practical procedure, it enables us to derive a (tight) upper bound for the capacity of the network under more standard routing policies. Using this new expression, we compute network capacity for a sample of demand matrices for some prototypical networks, including uniform demand (one unit between all node pairs) and broadcast demand (one unit between a source node and each other node as destination), and finally derive estimates of the mean load in some asymptotic cases.

computational schemes for the minimum cost routing of demands. There are no known closed-form solutions, or complete characterizations, of the solution to this problem under standard routing policies, such as shortest path or minimum cost routing. The reader is invited to consult multiple surveys of the state of the art in this literature (Kennington 1978;Assad 1978;Ouorou et al. 2000;Awerbuch and Leighton 2005;Wang 2018a, b). By contrast, in this paper, we report on a closed-form solution that provides a complete description of the mean load at each node and each link of a network under stochastic routing of demand governed by the weights of the edges of the network, as we describe below.
Single commodity or multicommodity network flow models in communication, transportation and numerous other settings typically assume shortest path routing. There are natural settings in which alternative routing not involving shortest paths may be required. For example, it may happen that longer routes are used for load balancing or, in the case of capacitated networks, to avoid network expansion (Minou 1989;Magnanti and Wong 1984). Or an inverse problem may be posed: to determine edge weights so that shortest path routes identified by these weights result in the smallest load across the network (Applegate and Cohen 2003). Given the universality of the network flow model, there is a vast number of applications of this model to numerous areas, in particular network design, bandwidth and server capacity computation (Davis et al. 2002;Bienstock et al. 2006), network coding (Bassoli et al. 2013;Jalali et al. 2016) and search algorithms (Chung 1997).
But a key application of the multicommodity flow problem remains the same: determining the network capacity that meets a specified demand matrix, given the connectivity of the network. This is a problem of general interest to which there are no closed-form solutions that compactly and efficiently leverage the intrinsic features of the network and the demand matrix to give a direct analytical answer, the focus of this paper. Determination of even (tight) upper or lower bounds today requires a use of numerical computation. In this report, we provide one such solution that gives (tight) upper bounds for nodal and edge capacities.
As stated above, there are few analytical results concerning the multicommodity flow problem with shortest path routing, in the sense of having a closed form solution as a function of a small number of parameters characterizing the network and the commodities or point to point demands. These include the characterization of the maximal load for hyperbolic graphs (Narayan and Saniee 2011;Jonckheere et al. 2011). In this setting, for a network of N nodes one assumes 1 unit of (directed) flow between all N (N − 1) node pairs, and then asks how the load scales due to shortest path routing as a function of N. This measure is sometimes referred to as the betweenness centrality, see Newman (2010) used for ranking of nodes for search algorithms.
In this paper, we study the near opposite of shortest path routing: when flows are routed in a uniformly or weighted random manner, each flow starting from its source and moving at each step randomly to a neighboring node and only stopping when the destination of the flow is reached. An analytical solution to this problem, if one could be found, would serve as an upper bound for the network capacity under various more realistic shortest path routing models discussed above. In a previous paper (Narayan et al. 2020), this problem was solved for the special case of the uniform multicommodity flow problem, that is, when all source-destination demands are present and are equal to one. Further, it was assumed that there were no biases in random routing. Here we generalize this result further by allowing (1) arbitrary and non-symmetric demand and (2) probabilistic routing that is driven by weights on the network edges. More specifically, we consider the case when an arbitrary amount of traffic d ij , in discrete units of messages, or packets from hereon, is injected into the network at every time step and at each node i for each possible destination node j = i. Thus there are i =j d ij units of traffic (packets/messages) injected into the N-node network at every time step, and this traffic is removed from the network when it reaches its destination. The network is assumed to be connected, i.e., have a single component.
We first demonstrate that a steady-state distribution is achieved, and then derive a closed-form expression for the expected flow, or the average number of packets or messages passing through each node, in terms of both the graph and demand (matrix) Laplacians. Next, we show how the result simplifies considerably when the demand matrix is symmetric. Following that, we derive a closed-form expression for when we each message only once for a node it passes through, even if it passes through multiple times, thus providing a tighter upper bound on the network capacity needed to process the required load. Finally, to illustrate the results more concretely, we use the derived analytical results to estimate the largest mean load for a few archetypal demand matrices over prototypical networks.

Weighted random routing with repetition
We shall use the following notation in the sequel.

Time evolution equations
As described in the previous section, we consider an undirected connected graph G(N, E) with N nodes, in which packets of traffic are injected at various nodes in a deterministic manner and move towards specified destinations. The dynamics are discrete time, i.e. packets of traffic move from node to node at time t = 0, 1, 2, 3 . . . . At each time step, therefore, the number of packets at a node depends on the number of packets at all its adjacent nodes at the previous time step, and any packets injected at their source or removed at their destination. We first obtain the equations for this time evolution.
At every time step, the number of packets injected at a node k with destination node l is given by T kl , which is the kl'th element of the demand matrix T. Furthermore, any packet at node i at time t that has not yet reached its destination (i.e. whose destination is not the node i) moves to one of the nodes, j, adjacent to node i at time t + 1. The move from node i to an adjacent node j is chosen randomly, with probability Here we use the notation i ∼ j if the nodes i and j have an edge connecting them, and i ∼ j if they do not. A is the weighted adjacency matrix defined in Table 1, with A ij = 0 if i ∼ j. If i ∼ j, the matrix element A ij is a positive real number which reflects the weight (1) of the edge (ij), i.e. whether the edge should be favored or disfavored when traffic flows through the network. We assume that the graph is undirected, i.e. that the matrix A is symmetric.
As stated before, any packet of traffic that is at its destination at time t is removed from the network and is no longer present at time t + 1. On the other hand, if a packet returns to its source as it moves around randomly, it continues as it would from any other node. The same is true if it returns to a node it has previously visited, i.e. if its path includes a closed loop. We return to this point in the next section of this paper.
We are interested in the expected value of the number of packets at each node. We expect that in steady state, if and when it exists, packets are removed from any node i at the rate of j T ji per time step. Recall the number of packets injected at the node i is j T ij per time step. We seek to find the steady state load, i.e. the average number of packets at each and every node of the network. This load may be interpreted as how fast on average each node needs to route packets it received net of those terminating at it. We first need to ensure there is a single and unique equilibrium distribution for the load in the network under the above conditions. The reader may refer to a simplified version of this model, in which T ij = 1 for all i = j and A ij = 1 for i ∼ j, which was discussed in earlier work (Narayan et al. 2020).
Theorem 1 For a connected graph G(N, E), with a deterministic injection rate of T ij packets per time step at node i destined for node j for each pair of nodes (ij), where each packet is routed randomly from its current node to its neighbors according to the weights in the adjacency matrix A until it reaches its destination, there exists a unique steady state number of packets at each node.
Proof Let us first consider the case when only one element in the demand matrix is non-zero, and this element is unity. That is, T ij = δ ik δ jl for some choice of k and l. With this demand matrix, let p kl i (t) be the ensemble average-i.e. averaged over many realizations of the same network-of the number of packets at node i at time t, except that p kl l (t) = 0, reflecting the fact that packets are immediately removed when they reach the destination node. The rate equation for the p i 's is then Table 1 Notation Here d j is the (weighted) degree of node j, as defined in Table 1. The boundary condition p kl l is an example of a Dirichlet boundary condition, where a function is defined in a region and is specified to be zero on the boundary of the region; in this case, the boundary is the node l and the region is all the other nodes in the graph.
We now show that, under the time evolution of Eq. (2), the function p kl i (t) reaches a t-independent unique steady state. Let p kl(1) i (0) and p kl(2) i (0) be two initial configurations at t = 0, that evolve according to Eq. (2). Define r kl Then r kl satisfies with the Dirichlet boundary condition at i = l. This is equivalent to where L is the Laplacian: Note that the matrix A is the weighted adjacency matrix. The real symmetric matrix L has a complete set of eigenvectors, with eigenvalues ≥ 0. The eigenvalue = 0 only exists if one can construct a function f on the graph for which f i = f j for all the nodes, in which case it is nondegenerate. This is not possible with Dirichlet boundary conditions, so that > 0. Thus the operator I − L (with Dirichlet boundary conditions) is a contraction. Therefore r kl (t → ∞) → 0, and as t → ∞ all initial configurations tend to the same t-independent steady state configuration.

Steady state solution
Having obtained the time evolution equation Eq. (2), we find its fixed point in this section, with Dirichlet boundary conditions as introduced in the proof of Theorem (1). As before, { α , α < N } represent the eigenvalues of the graph Laplacian.
In steady state, the arguments t and t + 1 in Eq.
(2) can be dropped. Moreover, we know that, of the traffic flowing from the node k to the node l, the load that flows into l at any time step must be equal to the load injected into the node k, i.e. unity. Therefore A lj p kl j /d j = 1, and we can extend Eq.
(2) as for all i, with the additional condition p kl l = 0. It may seem that we have gained nothing by restricting our analysis to the steady state configuration, since we still have to impose Dirichlet boundary conditions at the l'th node. However, as we shall see immediately, the solution to Eq. (6) can easily be found in terms of the eigenvectors of the Laplacian without the Dirichlet boundary condition, i.e. independent of k and l. We rewrite Eq. (6) in vector form as where p kl is a column vector whose i'th entry is p kl i , the matrix D is a diagonal matrix with D ij = d i δ ij , and v k and v l are unit column vectors whose i'th elements are δ ik and δ il respectively. Defining p kl = Dr kl , this equation is equivalent to with r kl l = 0. If the matrix L were invertible, we could obtain r kl = L −1 (v k − v l ), but L has a zero eigenvalue. So we define the matrix M = L + P, where P is a constant matrix with entries P ij = 1/N . Noting that the eigenvector ξ 0 of L with = 0 is a constant vector, with entries equal to 1/ √ N to ensure normalization, the matrix P projects onto ξ 0 . The eigenvalues of the matrix M are then equal to α + δ α0 in terms of the eigenvalues of the Laplacian. In terms of this matrix, and its inverse M −1 , where we have used the fact that MP = LP + P 2 = P 2 = P, i.e. P = M −1 P. The second term on the right hand side, Pr kl , is a column vector with all entries identical. Therefore in component form, For a general demand matrix, we arrive at the result:

Theorem 2 The total steady state load at the j'th node is
The term T j = k T kj is added to account for the processing of the traffic at its destination node, based on the assumption that this is the same as the amount of processing at any intermediate node on its path; this term can be easily modified if the assumption is changed. If the demand matrix is symmetric, this simplifies to a more compact form with a natural interpretation.

Theorem 3 When the demand matrix is symmetric, the total steady state load at the j'th node is
It is also possible to write Eq. (11) in terms of the eigenvectors and eigenvalues of the Laplacian, L. Using the result When the demand matrix is symmetric, this reduces to Further, when the demand matrix is uniform, i.e., T kl = 1 − δ kl , since the first summed expression in Eq. (14) is clearly zero when k = l, we can drop the factor of T kl . Since l ξ α l = 0 for α = 0 and l (ξ α l ) 2 = 1 by normalization, this reduces still further to which was derived in earlier work (Narayan et al. 2020).
Remark 1 For a symmetric demand matrix, the result that j is linearly dependent on d j can be obtained directly. The traffic from node k to node l can be represented as a stream of random walkers that diffuse through the network at discrete time steps. At every time step, in addition to the diffusive dynamics, T kl walkers are introduced at node k, and all the walkers at node l are removed. Comparing with Eq.
(2), the expected number of random walkers at node j at time t is equal to p kl j (t). If the random walks corresponding to all source destination pairs take place simultaneously, with each walker labelled with an index corresponding to its destination, we have random walkers with N different labels moving through the network. In addition to the random walk dynamics, walkers are created and destroyed at their sources and destinations respectively. In steady state, the number of walkers created and destroyed at the j'th node at any time step are equal to i T ji and i T ij respectively, but they have different labels. If we ignore the labels on the random walkers, the creation and destruction of random walkers can be ignored. The steady state solution for k l p kl j (t) is proportional to the steady state solution for a diffusion process on the graph with no sources or sinks. It is easy to verify that, in this steady state, the number of random walkers at the j'th node is proportional to i A ji = d j . Although this tells us that [ j −T j ] is a constant, independent of j, it does not tell us that this constant is equal to Tr[(T − T )M −1 ].
Remark 2 So far we have dealt with connected undirected graphs. We point out that when the graph is directed, then assuming that steady state distribution is achievable, Remark 1 implies that the expected load � j =T j + Cπ j where C is some constant independent of the node and (π j ) is the right principal eigenvector (with eigenvalue 1) of the random walk matrix for the directed graph, which for undirected graphs is equal to (d j ).

Remark 3
We observe that the proofs of both theorems carry through essentially unchanged if we replace the deterministic arrival of one packet at each source node for each destination node at each time step with a Poisson arrival process with a mean of one packet arrival per node per unit time for each destination node.

Remark 4
Specializing to the case of a uniform demand matrix, where Eq. (15) applies, in the large-N limit, the spectral density of the Laplacian α δ( − α ) tends to N ρ( ), where ρ( ) is smooth, and which is ∼ N 2 if the integral has a spectral gap or is otherwise convergent. On the other hand, if ρ( → 0) is non-zero, the spectral gap for large N is proportional to 1/N, and N 2 (ρ( )/ )d ∼ N 2 ln N . The maximum load at any node in the network is, up to an additive constant, Nd max 1/ α , which scales differently than the average load if d max diverges as N → ∞.

Random routing without repetition
With different edges in the graph weighted differently, random routing preferentially follows paths that have edges with high weights. However, this generalized model suffers from an important shortcoming: if A ij is large for an edge connecting a node i to a node j, then a packet that reaches node i is likely to move to j at the next time step, and vice versa. The packet may spend several time steps going to and fro between the two nodes. If A ij is increased with the weights of other edges held constant, this behavior becomes more pronounced. The multiple visits to nodes i and j increase the load at these nodes. Clearly, this is not how any reasonable routing scheme, even with randomness, would function. One way in which to eliminate the enhancement of the load at such nodes is to count each node on the path of a packet only once, regardless of how many times it is visited. In effect, this is as if the back and forth travel between the two adjoining nodes were to be eliminated.This is the strategy that we use in this section of the paper.
Counting each node on a path only once is an imperfect solution, since in addition to eliminating such back and forth paths, it also affects paths that include closed loops. The node at which the path enters and exits the loop, through which it goes twice, is only counted once even though the loop is not being eliminated. (If the loop were to be eliminated, all the other nodes in the loop would be counted zero times, not once.) However, this solution eliminates the greatest shortcoming of random routing for a weighted graph in the previous section of this paper. The analysis with an unweighted network and uniform demand between all source-destination pairs has been reported in Narayan and Saniee (2018), which we extend here. As before, we start with packets that are sourced at node k and have their destination as node l, and perform a random walk through the network. At each time step, a packet is removed from the network if it has reached the destination l or a query node m. The steady state equation is then a modification of Eq. (7): where µ kl m is determined by the condition that p kl l = p kl m = 0. The fraction µ kl m of the traffic between k and l that is removed at node m is the fraction that passes through the node m at some point in its path from k to l. In reality, the packet continues on after it encounters m, perhaps passing through m multiple times, before it reaches its destination l, but removing it from the network ensures that it is counted only once at the node m.
Defining r kl as before, we obtain The solution to this equation is a modification of Eq. (10): Applying the condition r kl l = r kl m = 0, and subtracting one from the other, we obtain Note that µ kl k = 1, but µ kl l is indeterminate because the terms with v l and v m in Eq. (17) merge together. Since all paths from k to l pass through both k and l, we fix µ kl l = 1. As noted earlier, µ kl m is the fraction of the traffic from node k to node l that passes through the node m, possibly multiple times. Summing over all source-destination pairs, the load at node m with repeated traversals counted only once is given by the following.

Theorem 4 The load at node m with repeated traversals counted only once is given by
As before, this can be expressed in terms of the eigenvectors and eigenvalues of the Laplacian, yielding

Numerical results
In this section, we provide numerical results for three different types of graphs: square lattice graphs, hyperbolic grids, and Erdös-Rényi random graphs (Erdös and Rényi 1959) for two types of demand regimes: broadcast with constant value and broadcast exponential decay. Numerical results for the uniform demand matrix for hyperbolic grids, Erdös-Rényi graphs, and scale-free networks (Barabási and Albert 1999); Dorogovtsev et al. (2000) have been presented earlier (Narayan et al. 2020).

Square lattice
We consider a L 0 × L 0 square lattice with N = L 2 0 nodes. Periodic boundary conditions are used, so that every node has exactly four neighbors. It is convenient to choose the normalized eigenfunctions of the Laplacian as the complex functions where r = (x, y) is the location of the node, with x and y integers, and q = (q x , q y ) = 2π(m, n)/L 0 , with n, m = 0, 1, 2 . . . L 0 − 1, and m = n = 0 not allowed. For such complex eigenfunctions, Eq. (13) is generalized to For the case of unit demand between all node pairs, from Eq. (16) and the following discussion, the load at any node grows as ∼ N 2 ln N for large N. If we change the normalization so that the total demand between all node pairs adds up to 1, i.e. T kl = (1 − δ kl )/[N (N − 1)], the load at any node grows as ∼ ln N for large N.
Instead, we consider the case when all the traffic is broadcast from a single node k, which we choose to have coordinates x = y = 0. (Because of periodic boundary conditions, all nodes are equivalent.) Eq. (24) then becomes where we have used the fact that q = −q .
If T 0l is independent of the destination node l, normalized to T 0l = (1 − δ l0 )/(N − 1) to ensure unit total demand, we can replace this with T 0l = 1/(N − 1) inside the sum, because 1 − exp[i(r l · q] = 0 when l is at the origin. The sum over l can easily be performed, and the l-dependent term cancels out when summed, resulting in In the large-N limit, this approaches where the sum diverges as N → ∞, and the integral is finite. It is possible to argue, and to verify numerically in Fig. 1, that the large-N behavior of the sum is (2/π ) ln N + const. Numerically, the constant is 0.3901 . . . . It is similarly possible to obtain that the integral approaches −(2/π ) ln |r j | + const for large |r j |. Noting that the maximum distance between any pair of nodes is O( √ N ), the load decreases as one moves away from the source, with a singular part that ranges from (2 ln N )/π near the source to (ln N /)π far from the source. This singular term is not present with geodesic routing.
Alternatively, if the demand falls off exponentially with distance from the source node, it is possible to take the large-N limit without any singular terms. We assume that the demand between the source node at (0, 0) and a destination node at (x, y) = (0, 0) is proportional to p |x|+|y| . A normalization constant of (1 − p) 2 /(4p) ensures unit total demand, matching the previous paragraph. Then where, as in the previous paragraph, it is not necessary to specify (l x , l y ) = (0, 0) in the sum. The sums over l x and l y are then independent and easy to perform. In the large-N limit, we can then replace N −1 q with 2π 0 2π 0 dq x dq y /(4π) 2 , since the zero in q is canceled by a zero in 1 − exp[i(l x q x + l y q y )] for small-q, resulting in the absence of singular terms in the large-N limit. The resulting integrals can be evaluated numerically, but an inspection of Eq. (28) shows that the j-dependent term 1 + exp[−ir j · q] is a rapidly oscillating function of q, averaging to unity, when |r j | → ∞. Thus j far from the source is half its value near the source, even though the demand is concentrated near the source. Fig. 1 a Load as a function of the logarithm of the distance r from the center of a square lattice. Demand is only broadcast from the center, uniformly to all other nodes. For definiteness, the load is shown for nodes along the x axis, for a lattice with N = 1600 2 . The dashed line has a slope of −2/π . b Load at the center as a function of N, the number of nodes in the square lattice, when demand is broadcast from the center equally to all nodes. The (2/π) ln N dependence is in agreement with the analytical results The r-dependence of the load can be understood using a continumm approximation, uaing an analogy to electrostatics. The source and destination of traffic map to positive and negative charges respectively, and the load maps to the electrostatic potential. The details of the analysis are outside the scope of this paper, but for demand broadcast equally to all nodes, the approximation predicts a ∼ ln r dependence of the load, while for demand decaying exponentially with distance, �(r) − �(∞) ∼ 1/r 4 for large r. Both of these are consistent with our numerical results.

Hyperbolic grids
A hyperbolic grid H p,q is an infinite regular planar graph with constant degree q, and each face a p-sided polygon, satisfying the condition (p − 2)(q − 2) > 4. If (p − 2)(q − 2) = 4, it is possible to make all the polygons as regular polygons of the same size. This is no longer possible if (p − 2)(q − 2) > 4, but one can interpret the graph as a planar projection of a graph on a surface of constant negative curvature, and on that surface, the polygons are regular and identical.
When the demand between all pairs of nodes (i, j) with i = j is equal to 1/(N (N − 1)), so that the total demand is 1, Eq. (15) is adjusted to The second term is negligible when N → ∞, while the first term was shown numerically Narayan et al. (2020) to scale linearly with ln N for large N when (p, q) = (3, 7).
Here, we show numerical results when demand is only broadcast from the node at the center, equally to all other nodes. To have the same total demand as in the previous case, the demand for each node is 1/(N − 1). As seen in Fig. 2, the load decays as �(r) = C exp[−cr] + �(∞), where r is the geodesic distance between the center and the node. Furthermore, the load at the center grows linearly with the radius of the grid, i.e. ∼ ln N , which is the same as for uniform demand on the hyperbolic grid, and also similar to the analytical results for the square lattice, when demand is broadcast from the center and is independent of the distance to the destination. We also consider the case when the demand is broadcast from the central node, but drops off exponentially as a function of the distance r from the center, i.e. proportional to p r , with the proportionality constant chosen so that the total demand is unity, irrespective of N. We choose p = 0.25; the total demand on an infinite grid is finite as long as p < (3 − √ 5)/2. The results for this case are not shown here, but within numerical accuracy, �(r) − �(∞) decreases exponentially as a function of r. The load is found to be almost independent of N, similar to the analytical results for exponentially decaying demand on the square lattice.
As with the square lattice, the r-dependence of the load can be understood in a continuum approximation with an electrostatic analogy, which predicts that �(r) − �(∞) will decay exponentially with r.

Random graphs
We performed numerical simulations on Erdös-Rényi random graphs, with the probability of an edge being connected being p = 2/N for a graph with N nodes. We only retain the giant component of the graph, which has O(N) nodes, numerically found to be ≈ 0.8N . For each value of N, several random graphs were generated, so that the total number of nodes in the giant components of all the graphs was approximately 10 4 . Demand was broadcast from the node at the center of each graph, either uniformly to all other nodes, or decaying as ∼ (1/4) r with the distance r from the center to the destination.
As shown in Fig. 3, rather surprisingly, the average load at the nodes at a distance r from the center has the same form when the demand is uniform or exponentially decaying, up to an overall proportionality constant. The average load �(r) for nodes at a distance r from the center does not have a simple functional form. There is also very little change in the plots as a function of N, consistent with our earlier result Narayan et al. Fig. 3 a Average load �(r) for nodes at a distance r from the center, from which demand is broadcast, for the Erdös-Rényi model with p = 2/N and N = 2500. Filled black circles are �(r) when the demand is the same, irrespective of the destination. Open red circles are 2.5�(r) when the demand decreases as (1/4) r ; the factor of 2.5 is chosen by eye to best match the two cases. The total demand for all destination nodes is normalized to 1. b Average load �(r) as a function of r for various values of N, for the broadcast model with the same demand for all destination nodes. The difference between the curves is small (2020) that the load at the center with uniform demand between all node pairs is independent of N.

Conclusions
We have shown that, for the general non-uniform multicommodity flow problem on an arbitrary connected graph under weighted random routing, the mean load at each node of the graph exists and is unique. We have also derived an explicit expression for it in terms of the graph and demand (matrix) Laplacians. Further, we have derived a closed form expression for the mean load without counting repeated visits to a transit node. Using these explicit expressions, we have obtained numerical estimates for the mean load for the regular and hyperbolic lattices in the large-size regime, and have numerically computed the mean load for Erdös-Rényi random graphs. This has been done for uniform and for broadcast demand. For the latter (broadcast demand), we have observed that the average load decreases as a function of distance from the broadcast source, ranging from exponential decrease (for regular hyperbolic graphs), to logarithmic decrease (for regular grids), to staged decrease (for Erdös-Rényi random graphs).