# Ultrametricity of optimal transport substates for multiple interacting paths over a square lattice network

###### Abstract

We model a set of point-to-point transports on a network as a system of polydisperse interacting self-avoiding walks (SAWs) over a finite square lattice. The ends of each SAW may be located both at random, uniformly distributed, positions or with one end fixed at a lattice corner. The total energy of the system is computed as the sum over all SAWs, which may represent either the time needed to complete the transport over the network, or the resources needed to build the networking infrastructure. We focus especially on the second aspect by assigning a concave cost function to each site to encourage path overlap. A Simulated Annealing optimization, based on a modified BFACF algorithm developed for polymers, is used to probe the complex conformational substates structure at zero temperature. We characterize the average cost gains (and path-length variation) for increasing polymer density with respect to a Dijkstra routing and find a non-monotonic behavior as recently found for random networks. We observe the emergence of ergodicity breaking and of non-trivial overlap distributions among replicas when switching from a convex to a concave cost function (e.g., , where represents the node overlap). Finally we show that the space of ground states for is compatible with an ultrametric structure as seen in many complex systems such as some spin glasses.

###### pacs:

Valid PACS appear hereThe problem of optimal transport over various kinds of networks is important both for theoretical and practical reasons Jayanth R et al. (1999). Areas of application range from river networks Rinaldo et al. (2014); Ijjász-Vásquez et al. (1993) to vascular systems in animals and plants Katifori et al. (2010); Corson (2010), and from electric energy distribution systems Bohn and Magnasco (2007) to communication networks Banavar et al. (2000); Yeung et al. (2013). The adoption of a cost function minimization scheme has allowed a unified approach to very diverse research fields. Cost functions may be thought as energy dissipation for electricity grids, time delay and/or resources needed to build the networking infrastructure. In recent years, the relation between the properties of the cost function and the associated optimal solutions has been extensively studied Banavar et al. (2000). For instance, when multiple sources are connected to a single destination (as in drainage basins), it is known that a concave cost leads to multiple (nearly equivalent) spanning trees, whereas a convex behavior shows a unique redundant solution with many loops. When the character of the cost function is not well defined, no a priori conclusions may be drawn Yeung and Saad (2012).

It is well known Banavar et al. (2000) that, for concave cost functions, a multiplicity of local optimal solutions exists. Since a hierarchical organization of the states is observed, one may hypothesize an ultrametric relation among them. Ultrametricity (UM) is one of the key features of the mean-field Parisi picture for spin glassesMezard et al. (1987): the states of the system obey an UM distance which translates into a hierarchical organization. This behavior has been hypothesized or observed for different polymer systems with noise, such as directed polymers in random media (DPRM) Zhang (1987); Kardar and Zhang (1987); Pang and Halpin-Healy (1993) and for self-interacting SAWs in external fields Fernández (1991).

In this paper we propose a model having three main differences with respect to the DPRM formulation: 1) the (locally) minimum cost is achieved by collectively optimizing several interacting chains; 2) polymers are polydisperse and have at least one end which is randomly located on the lattice meaning that quenched disorder is achieved through random topology rather than noisy bonds; 3) our polymers may not be directed: they are free to wander backwards/sideways to achieve a global cost gain. This behavior is seldom observed for weakly interacting polymers or very dilute systems.

Our analysis reveals a non-monotonic behavior of the optimized cost gain, with respect to a Dijkstra routing, when increasing polymer density as previously found for random networks Yeung et al. (2013). Moreover, we observe the emergence, only for concave costs, of ergodicity breaking and of non-trivial overlap distributions among replicas. Finally, we present evidence that the space of ground states is compatible with an ultrametric structure.

## Model and Numerical Methods –

We consider a square lattice network of nodes and side , each node connected to its four nearest neighbors via uniform links with adjacency matrix , zero elsewhere. A set of communications, modeled as polymers with fixed ends, compete on the network for the available resources, each occupying a path described by an interacting SAW. The self-avoidance condition enforces that no path uses the same node more than once, whereas distinct polymers can use the same node. The occupation number of each path on the node is denoted by . The total occupation number on the -th node is . The interaction among polymers is regulated by a Hamiltonian of the form

(1) |

where the concave/convex character of the cost function makes the system behave in qualitatively different ways Banavar et al. (2000). The occupation number is normalized by in order to have a uniform temperature behavior with respect to polymer multiplicity. The are node-dependent weights, in general set to unity, which may be used to model spatial non-uniformity. In this paper we consider only and the simplest functional form for the cost : the power function . This functional dependence leads to polymer repulsion for , whereas encourages overlap. For the polymers do not interact, and it is known that the ground state of is attained by (highly degenerate) shortest-path routing Yeung et al. (2013); Banavar et al. (2000).

To explore the energy landscape of for , we adopt a Simulated Annealing (SA) scheme in which temperature is gradually decreased to zero within a canonical Monte Carlo (MC). The basic MC move follows the BFACF algorithm developed for lattice polymers Berg and Foerster (1981); Aragão de Carvalho and Caracciolo (1983). With respect to the original scheme, at each iteration one polymer is randomly selected and instead of applying the basic move at a single random site, we perform multiple moves. The number of moves on the SAW is randomly chosen within one and the average polymer length. This random choice of basic move multiplicity guarantees good MC acceptance rates Caracciolo et al. (1990). Since in this work we are dealing with several interacting polymers, we rely on the Gibbs acceptation factor of the MC to extract a chain of states from the canonical ensemble instead of direct generation as in Ref. Berg and Foerster (1981). Assigning the same probability both to path-enlarging and path-shrinking BFACF moves leads to very low acceptance rates. We solve this problem by tying the probability of path-enlarging moves to the MC temperature Binder (1995).

The square lattice with uniform factors induces on this problem some peculiarities not found in a continuous representation of space. In particular, the ground states for basic routing problems (i.e., involving non-interacting polymers) in two dimensions are intrinsically degenerate because solutions with the same energy exist between any pair of points with horizontal and vertical distance Mernik (2009). We restrict our study to two-dimensional lattices with no periodic boundary conditions (PBC) to better model realistic network topologies. Preliminary numerical results show that qualitatively similar conclusions may be drawn when PBC are imposed at the boundaries.

Detecting UM in finite-volume systems can be very difficult due to finite-size effects especially with no PBC Rammal et al. (1986); Hed et al. (2004); Katzgraber and Hartmann (2009); Katzgraber et al. (2012); Ciliberti and Marinari (2004).
To measure differences between replicas for the same quenched disorder, we define a path overlap , computed as the ratio of the common visited nodes (node-overlap ) or common visited links (link-overlap ) with respect to path length ^{1}^{1}1Whenever two strings possess different lengths, we concatenate a padding string to the short one..
For homologue instances of the SAW belonging to replicas and , we define the node overlap as

(2) |

and the link overlap

(3) |

We obtain compatible results for and , but the data shown in this paper is computed by using the link-overlap, so we define . From the overlap we get the normalized Hamming distance as .

In an UM space Rammal et al. (1986) the triangle inequality valid for metric spaces, is replaced by a stronger version . This inequality is equivalent to imposing that any triple of points should form an acute isosceles triangle or, at most, an equilateral one. In order to discern between trivial UM due to equilaterals or true UM due to acute isosceles, after performing a standard UM test Hed et al. (2004), we propose a new procedure that consists of analyzing the frequency distribution of each triple of ordered distances , by keeping only two transformed components as defined in Ref.Contucci et al. (2007): vs .

We compute the distance for all pairs of states within the same quenched disorder. This is achieved by selecting pairs of homologue SAWs belonging to two states and , and finally computing their normalized Hamming distance. By applying a clustering algorithm Hed et al. (2004); Katzgraber and Hartmann (2009); Ciliberti and Marinari (2004) to each disorder, we obtain a dendrogram such as the one in the left part of Fig. 1. The procedure starts with each state in a separate cluster, then iteratively the nearest clusters are merged. At each step the inter-cluster distances are recomputed by averaging among all pairs of their member states. The procedure ends when a root common cluster appears. We then sort all states in the same order as the dendrogram and finally plot the adjacency matrix: a well-visible block-diagonal complex structure emerges for (we found similar results for ), along with a deep hierarchy visible in the dendrogram. On the other hand, nearly flat hierarchies, somewhat uniform distance matrices and pseudo-normal overlap distributions are observed for (see Fig. S7). To probe for an UM space structure, we randomly select three configurations from the hierarchical cluster structure (see Ref. Hed et al. (2004)), resulting in three mutual distances that we sort to get and finally compute the correlator . If the phase space is UM, we expect for . Thus should converge to a delta function in for and the variance of the distribution . By this approach, UM would be detected even in the case of equilateral triangles (trivial UM). The alternative definition of in Ref. Katzgraber and Hartmann (2009) is difficult to apply for polydisperse SAWs (each SAW length has a distinct distance distribution), so we devised a supplementary test to rule out trivial UM.

## Results and discussion –

All ground states have been obtained by performing a SA energy minimization with an exponential cooling converging to at one third of the total simulation length. To characterize cost gain and path-length variation for the system with both random ends, we performed several minimizations in which the number of SA timesteps varied with system size from for and to for and . The maximum number of basic BFACF moves per timestep was proportional to . For each lattice size we generated quenched disorders with random uniform configurations and for each disorder we produced local ground states. We plan to expand this number for future works. Some examples of the ground states we obtained are shown in Figs. S1-S4.

The first goal has been to characterize the system as regards the attainable cost gain with respect to the shortest path routing that is widely used for many transport applications. In Fig. 2 we plot the energy difference ratio both for concave (, continuous blue online) and convex (, dashed red online) cost functions.

Since the average polymer length increases with , we observe in Fig. 2 the tendency of the curves to superimpose for . After a steep cost gain growth, maximum efficiency with respect to Dijkstra is reached for both values, then the value slowly decreases since most nodes are already busy and the advantage of longer-than-Dijkstra detours is weaker. The peaks are shifted for and : and , respectively. Cost gain ratios are quite different in the two scenarios: for the convex case, the gain ratio is relatively constant for any at nearly , while in the concave situation, it rarely goes beyond and decreases more markedly for higher values. By comparing the peak values of Fig. 2 with their associated path-length variations in Fig. 3, we show that large cost gains may be obtained by employing paths slightly longer than Dijkstra ( for both values). These results are qualitatively very similar to those obtained by an alternative optimization approach presented in Refs. (Yeung and Saad, 2012; Yeung et al., 2013) for random paths on a random graph with constant connectivity . It should be stressed that the regular lattice is not tractable with that method since too many loops exist, leading to severe ground state degeneration. The present method may be exploited to minimize network construction resources by obtaining a set of lean structures with (resource sharing encouraged) and then selecting the best performing ground states for (paths competing for resources).

Let us focus for the rest of the paper on the concave cost case. In the Supplementary Materials to Ref. Yeung et al. (2013), there is a brief discussion regarding the possibility of a Replica Symmetry Breaking (RSB) scenario for . This led us to question whether an UM structure among the ground states of our closely related system exists. We consider two polymer distributions: one in which both SAW ends are uniformly distributed and another in which one end is constrained to a lattice corner. The RSB scenario is apparent when looking at non-trivial overlaps Franz and Parisi (2000) among replicas at : we observe multimodal overlap distributions for roughly half of the SAWs. This fraction grows with . In the SM we compare individual dendrograms, distance matrices and distributions. For randomized paths and for the overlap distributions are always quasi-normal.

It has become customary to show a tendency towards UM by plotting the distribution for several system sizes along with its variance . In Fig. 4 the for the random-ends polymers (top) and for the fixed-end polymers (bottom) are shown. There is a visible trend of both to diverge at when . Both polymer topologies share the same overall behavior.
In the two insets of Fig. 4(top and bottom) we can compare the variance of the for both systems compared to their randomized counterparts ^{2}^{2}2The randomized system is obtained by reshuffling every SAW sequence with the same quenched disorder. The randomized paths do maintain legal connection between polymer ends, but self-avoidance may be violated.: they both tend to zero as lattice size grows, so we cannot still conclude in which case whether UM is attained in the thermodynamic limit.

To further investigate this issue, we study the distribution of distance triples for the fixed-end polymer system. In Fig. 5 we plot the (considering all quenched disorders) both for optimized and randomized polymers with increasing lattice size. Plots for the single quenched disorders are shown in Fig. S8. All optimized systems reveal higher concentrations (dark red) of triangles along the -axis that is the signature of true UM; for larger , the high region gets progressively depleted (green-yellow). The randomized plots show a that decreases with a constant gradient starting from a maximum in the origin (equilaterals). The mode is in all six cases in , but for the optimized systems its value is roughly half the sum of bins along the -axis. This is better shown by the white circles representing the center of mass (CM) for each distribution: with respect to the randomized states, in which the CM shifts toward the origin with , the optimized states tend, on average and for each quenched disorder, to stay just above the -axis without converging to . The fraction of trivial UM is due to “pathological”SAWs with very few ground states. Scalene triangles are produced by the shortest-path degeneracy due to local low-polymer density on the lattice (see the individual overlap distributions within the Supplemental Materials for details).

## Conclusions –

In this letter we presented a novel approach to explore the ground states of an important class of transport optimization problems on a regular square lattice. We showed that with this method one is able to obtain solutions for the optimal transport of a set of interacting communications spread over the lattice with different topological constraints. The interaction among paths is obtained within the unifying framework of concave/convex cost functions. The fact that this method works on a lattice, allows the possibility of discovering a hierarchy of the most inexpensive network infrastructures encouraging transport coalescence () and then to optimize the distilled graph for performance, fault tolerance and congestion resistance with a repulsive cost function () as in Ref. Yeung et al. (2013). We tested our optimization procedure by characterizing the global cost gain and path-length variation for a system of randomly spread point-to-point communications over a square lattice that has been recently studied on random graphs via an unconventional use of the replica and cavity methods Yeung et al. (2013) obtaining qualitatively similar results. The appearance, for , of families of hierarchically related solutions (tree-like as predicted by Banavar et al.Banavar et al. (2000)) led us to investigate whether RSB and an UM structure hold. Similarities and differences with respect to spin systems allowed us to borrow a standard approach to probe for UM, which we slightly extended by plotting the distribution of triangle types for growing lattice sizes. In conclusion we found evidence supporting the RSB scenario (non-trivial overlap distributions) and UM at the level of single interacting polymers as hypothesized for DPRM systems Zhang (1987). Here the equivalent of noise (highly correlated) is apparently played by all polymers minimizing the total energy thus forming a rough landscape Buldyrev et al. (2006); Braunstein et al. (2002). It should be further investigated whether and how this phenomenon depends on polymer density/dispersity and values. Finally, we plan to further explore the possibility of defining a global similarity measure, encompassing all SAWs, to assess whether a system-wise UM structure exists.

###### Acknowledgements.

We thank Francesco Versaci for useful discussions. We also acknowledge the contribution of Sardinian Regional Authorities under projects ABLE and Predict.## References

- Jayanth R et al. (1999) A. B. Jayanth R, M. Amos, and Rinaldo, Nature 399, 130 (1999).
- Rinaldo et al. (2014) A. Rinaldo, R. Rigon, J. R. Banavar, A. Maritan, and I. Rodriguez-Iturbe, Proceedings of the National Academy of Sciences of the United States of America 111, 2417 (2014).
- Ijjász-Vásquez et al. (1993) E. J. Ijjász-Vásquez, R. L. Bras, I. Rodríguez-Iturbe, R. Rigon, and A. Rinaldo, Advances in Water Resources 16, 69 (1993).
- Katifori et al. (2010) E. Katifori, G. J. Szöllosi, and M. O. Magnasco, Physical Review Letters 104, 048704 (2010), arXiv:0906.0006 .
- Corson (2010) F. Corson, Physical Review Letters 104, 048703 (2010), arXiv:0905.4947 .
- Bohn and Magnasco (2007) S. Bohn and M. O. Magnasco, Physical Review Letters 98, 088702 (2007), arXiv:0607819 [cond-mat] .
- Banavar et al. (2000) J. R. Banavar, F. Colaiori, A. Flammini, A. Maritan, and A. Rinaldo, Physical Review Letters 84, 4745 (2000).
- Yeung et al. (2013) C. H. Yeung, D. Saad, and K. Y. M. Wong, Proceedings of the National Academy of Sciences of the United States of America 110, 13717 (2013), arXiv:arXiv:1309.0745v1 .
- Yeung and Saad (2012) C. H. Yeung and D. Saad, Physical Review Letters 108, 208701 (2012), arXiv:1202.0213 .
- Mezard et al. (1987) M. Mezard, G. Parisi, and M. Virasoro, Spin Glass Theory and Beyond, Lecture Notes in Physics Series (World Scientific, 1987).
- Zhang (1987) Y.-C. Zhang, Phys. Rev. Lett. 59, 2125 (1987).
- Kardar and Zhang (1987) M. Kardar and Y. C. Zhang, Physical Review Letters 58, 2087 (1987).
- Pang and Halpin-Healy (1993) N.-N. Pang and T. Halpin-Healy, Phys. Rev. E 47, R784 (1993).
- Fernández (1991) A. Fernández, International Journal of Theoretical Physics 30, 83 (1991).
- Berg and Foerster (1981) B. Berg and D. Foerster, Physics Letters B 106, 323 (1981).
- Aragão de Carvalho and Caracciolo (1983) C. Aragão de Carvalho and S. Caracciolo, J. Phys. France 44, 323 (1983).
- Caracciolo et al. (1990) S. Caracciolo, A. Pelissetto, and A. Sokal, Journal of Statistical Physics 60, 1 (1990).
- Binder (1995) K. Binder, Monte Carlo and Molecular Dynamics Simulations in Polymer Science (Oxford University Press, 1995).
- Mernik (2009) L. Mernik, J. of Pure and Applied Math 56, 589 (2009).
- Rammal et al. (1986) R. Rammal, G. Toulouse, and M. A. Virasoro, Reviews of Modern Physics 58, 765 (1986).
- Hed et al. (2004) G. Hed, A. P. Young, and E. Domany, Physical Review Letters 92, 157201 (2004).
- Katzgraber and Hartmann (2009) H. G. Katzgraber and A. K. Hartmann, Physical Review Letters 102, 037207 (2009), arXiv:arXiv:0807.3513v1 .
- Katzgraber et al. (2012) H. G. Katzgraber, T. Jorg, F. Krzakala, and A. K. Hartmann, Physical Review B - Condensed Matter and Materials Physics 86, 184405 (2012), arXiv:arXiv:1205.4200v1 .
- Ciliberti and Marinari (2004) S. Ciliberti and E. Marinari, Journal of Statistical Physics 115, 557 (2004).
- (25) Whenever two strings possess different lengths, we concatenate a padding string to the short one.
- Contucci et al. (2007) P. Contucci, C. Giardinà, C. Giberti, G. Parisi, and C. Vernia, Physical Review Letters 99, 057206 (2007).
- Franz and Parisi (2000) S. Franz and G. Parisi, European Physical Journal B 18, 485 (2000).
- (28) The randomized system is obtained by reshuffling every SAW sequence with the same quenched disorder. The randomized paths do maintain legal connection between polymer ends, but self-avoidance may be violated.
- Buldyrev et al. (2006) S. V. Buldyrev, S. Havlin, and H. E. Stanley, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 73, 036128 (2006).
- Braunstein et al. (2002) L. A. Braunstein, S. V. Buldyrev, S. Havlin, and H. E. Stanley, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 65, 056128 (2002).