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.

The study of network capacity, or load, for arbitrary demand matrices is over half a century old, and goes back to the pioneering work of Ford and Fulkerson (1956) and Elias et al. (1956) for the single commodity and to early attempts (Hu 1963; Lomonosov 1985; Shahrokhi and Matula 1990; Papernov 1976; Okamura and Seymour 1981) for the multicommodity flow solutions of the problem. This rather large literature aims to provide characterizations of the load, or the minimal capacity required to route the specified demand, in terms of nodal or link capacities needed, based on cut values. For the case of the single commodity model cut values provide both necessary and sufficient conditions, while for the multicommodity case, at best they provide necessary conditions.

The existing literature on multicommodity flow problem typically deals with either necessary or sufficient conditions for the existence of solutions or, more frequently, 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\ne i.\) Thus there are \(\sum _{i \ne 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.

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\ldots .\) 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\sim j\) if the nodes i and j have an edge connecting them, and \(i\not \sim j\) if they do not. A is the weighted adjacency matrix defined in Table 1, with \(A_{ij} = 0\) if \(i\not \sim j.\) If \(i\sim j,\) the matrix element \(A_{ij}\) is a positive real number which reflects the weight 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 \(\sum _j T_{ji}\) per time step. Recall the number of packets injected at the node i is \(\sum _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\ne j\) and \(A_{ij} = 1\) for \(i\sim 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} = \delta _{ik} \delta _{jl}\) for some choice of k and l. With this demand matrix, let \(p_i^{kl}(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_l^{kl}(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

for \(i\ne l,\) with the boundary condition \(p^{kl}_l(t+1) = 0.\) Here \(d_j\) is the (weighted) degree of node j,  as defined in Table  1. The boundary condition \(p_l^{kl}\) 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}_i(t)\) to be equal to \([p^{kl(1)}_i(t) - p^{kl(2)}_i(t)]/d_i.\) 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 \(\lambda \ge 0.\) The eigenvalue \(\lambda = 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 \(\lambda > 0.\) Thus the operator \(I - L\) (with Dirichlet boundary conditions) is a contraction. Therefore \(r^{kl}(t\rightarrow \infty )\rightarrow 0,\) and as \(t\rightarrow \infty\) all initial configurations tend to the same t-independent steady state configuration. \(\square\)

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, \(\{\lambda _{\alpha },\alpha <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 \(\sum 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 \(\mathbf{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 \delta _{ij},\) and \({\mathbf {v}}^k\) and \({\mathbf {v}}^l\) are unit column vectors whose i’th elements are \(\delta _{ik}\) and \(\delta _{il}\) respectively. Defining \(\mathbf{p }^{kl} = D {\mathbf {r}}^{kl},\) this equation is equivalent to

with \(r^{kl}_l = 0.\)

If the matrix L were invertible, we could obtain \({\mathbf {r}}^{kl} = L^{-1} ({\mathbf {v}}^k - {\mathbf {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 \(\varvec{\xi }^0\) of L with \(\lambda = 0\) is a constant vector, with entries equal to \(1/\sqrt{N}\) to ensure normalization, the matrix P projects onto \(\varvec{\xi }^0.\) The eigenvalues of the matrix M are then equal to \(\lambda _\alpha + \delta _{\alpha 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, \(P{\mathbf {r}}^{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

\(\square\)

The term \({\tilde{T}}_j = \sum _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

\(\square\)

It is also possible to write Eq. (11) in terms of the eigenvectors and eigenvalues of the Laplacian, L. Using the result \([M^{-1}]_{ij} = \sum _\alpha \xi ^\alpha _i (\lambda _\alpha + \delta _{\alpha 0})^{-1} \xi ^\alpha _j,\) which is in turn equal to \(\sum _{\alpha \ne 0} \xi ^\alpha _i \lambda _\alpha ^{-1} \xi ^\alpha _j + 1/N,\) we obtain

When the demand matrix is symmetric, this reduces to

Further, when the demand matrix is uniform, i.e., \(T_{kl} = 1 - \delta _{kl},\) since the first summed expression in Eq. (14) is clearly zero when \(k=l,\) we can drop the factor of \(T_{kl}.\) Since \(\sum _l \xi _l^\alpha = 0\) for \(\alpha \ne 0\) and \(\sum _l (\xi _l^\alpha )^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 \(\Lambda _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 \(\sum _i T_{ji}\) and \(\sum _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 \(\sum _k \sum _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 \(\sum _i A_{ji} = d_j.\) Although this tells us that \([\Lambda _j - {\tilde{T}}_j]\) is a constant, independent of j,  it does not tell us that this constant is equal to \(Tr[({\tilde{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 \(\Lambda _j = {\tilde{T}}_j + C \pi _j\) where C is some constant independent of the node and \((\pi _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 \(\sum _\alpha \delta (\lambda - \lambda _\alpha )\) tends to \(N\rho (\lambda ),\) where \(\rho (\lambda )\) is smooth, and

which is \(\sim N^2\) if the integral has a spectral gap or is otherwise convergent. On the other hand, if \(\rho (\lambda \rightarrow 0)\) is non-zero, the spectral gap for large N is proportional to 1/N,  and \(N^2 \int (\rho (\lambda )/\lambda ) d\lambda \sim N^2 \ln N.\) The maximum load at any node in the network is, up to an additive constant, \(N d_{max} \sum 1/\lambda _\alpha ,\) which scales differently than the average load if \(d_{max}\) diverges as \(N\rightarrow \infty .\)

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 \(\mu _m^{kl}\) is determined by the condition that \(p_l^{kl} = p_m^{kl} = 0.\) The fraction \(\mu _m^{kl}\) 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 \({\mathbf {r}}^{kl}\) as before, we obtain

The solution to this equation is a modification of Eq. (10):

Applying the condition \(r_l^{kl} = r_m^{kl} = 0,\) and subtracting one from the other, we obtain

Note that \(\mu _k^{kl} = 1,\) but \(\mu _l^{kl}\) is indeterminate because the terms with \({\mathbf {v}}^l\) and \({\mathbf {v}}^m\) in Eq. (17) merge together. Since all paths from k to l pass through both k and l,  we fix \(\mu _l^{kl} = 1.\)

As noted earlier, \(\mu _m^{kl}\) 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

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\times L_0\) square lattice with \(N = L_0^2\) 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 \(\mathbf{r } = (x, y)\) is the location of the node, with x and y integers, and \(\mathbf{q } = (q_x, q_y) = 2 \pi (m, n)/L_0,\) with \(n, m = 0, 1, 2\ldots 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 \(\sim 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 - \delta _{kl})/[N (N - 1)],\) the load at any node grows as \(\sim \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 \(\lambda _{\mathbf{q }} = \lambda _{-\mathbf{q }}.\)

If \(T_{0l}\) is independent of the destination node l,  normalized to \(T_{0l} = (1 - \delta _{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 ({\mathbf{r }}_l\cdot {\mathbf{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\rightarrow \infty ,\) 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/\pi ) \ln N + const.\) Numerically, the constant is \(0.3901\ldots .\) It is similarly possible to obtain that the integral approaches \(- (2/\pi ) \ln |{\mathbf{r }}_j| + const\) for large \(|{\mathbf{r }}_j|.\) Noting that the maximum distance between any pair of nodes is \(O(\sqrt{N}),\) the load decreases as one moves away from the source, with a singular part that ranges from \((2 \ln N)/\pi\) near the source to \((\ln N/)\pi\) far from the source. This singular term is not present with geodesic routing.

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/\pi .\) 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/\pi ) \ln N\) dependence is in agreement with the analytical results

Full size image

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)\ne (0, 0)\) is proportional to \(p^{|x| + |y|}.\) A normalization constant of \((1 - p)^2/(4 p)\) 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)\ne (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} \sum _{{\mathbf{q }}}\) with \(\int _0^{2\pi }\int _0^{2\pi } dq_x dq_y/(4\pi )^2,\) since the zero in \(\lambda _{\mathbf{q }}\) is canceled by a zero in \(\left( 1 - \exp [i (l_x q_x + l_y q_y)]\right)\) 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 \(\left( 1 + \exp [-i {\mathbf{r }}_j\cdot {\mathbf{q }}]\right)\) is a rapidly oscillating function of \({\mathbf{q }},\) averaging to unity, when \(|{\mathbf{r }}_j|\rightarrow \infty .\) Thus \(\Lambda _j\) far from the source is half its value near the source, even though the demand is concentrated near the source.

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 \(\sim \ln r\) dependence of the load, while for demand decaying exponentially with distance, \(\Lambda (r) - \Lambda (\infty ) \sim 1/r^4\) for large r. Both of these are consistent with our numerical results.

Hyperbolic grids

A hyperbolic grid \({\mathbb {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\ne 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\rightarrow \infty ,\) 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 \(\Lambda (r) = C \exp [-c r] + \Lambda (\infty ),\) 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. \(\sim \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.

a (Log scale) Load as a function of r,  the distance from the center, for a \({\mathbb {H}}_{3,7}\) grid with radius 8, i.e. with \(N=11{,}173.\) Demand is only broadcast from the center, uniformly to all other nodes. The plot of \(\ln [\Lambda (r) - \Lambda (\infty )]\) is a linear function of r,  implying exponential decay of the load as a function of r,  towards its asymptotic value at \(r\rightarrow \infty .\) b Load at the center as a function of the radius of the \({\mathbb {H}}_{3,7}\) grid. The load is a linear function of the radius, i.e. a linear function of \(\ln N\) for large N.

Full size image

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 - \sqrt{5})/2.\) The results for this case are not shown here, but within numerical accuracy, \(\Lambda (r) - \Lambda (\infty )\) 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 \(\Lambda (r) - \Lambda (\infty )\) 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 \(\approx 0.8 N.\) 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 \(\sim (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 \(\Lambda (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. (2020) that the load at the center with uniform demand between all node pairs is independent of N.

a Average load \(\Lambda (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 \(\Lambda (r)\) when the demand is the same, irrespective of the destination. Open red circles are \(2.5 \Lambda (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 \(\Lambda (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

Full size image

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).

Not applicable.

O. Narayan acknowledges useful conversations with Prof. Josh M. Deutsch.

Not applicable.

Authors and Affiliations

1. Department of Physics, University of California, Santa Cruz, CA, 95064, USA

   Onuttom Narayan

2. Mathematics & Algorithms Group, Bell Labs, Nokia, Murray Hill, NJ, 07974, USA

   Iraj Saniee

Contributions

The authors contributed equally to all aspects of this paper. Both authors read and approved the final manuscript.

Corresponding author

Correspondence to Iraj Saniee.

Competing interests

The authors declare that they have no competing interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions