††thanks: Equal contribution††thanks: Equal contribution

# Integer Traffic Assignment Problem: Algorithms and Insights on Random Graphs

Rayan Harfouche  Giovanni Piccioli<EMAIL_ADDRESS>Lenka Zdeborová
École Polytechnique Fédérale de Lausanne (EPFL)
Statistical Physics of Computation Laboratory

###### Abstract

Path optimization is a fundamental concern across various real-world
scenarios, ranging from traffic congestion issues to efficient data routing
over the internet. The Traffic Assignment Problem (TAP) is a classic
continuous optimization problem in this field. This study considers the
Integer Traffic Assignment Problem (ITAP), a discrete variant of TAP. ITAP
involves determining optimal routes for commuters in a city represented by a
graph, aiming to minimize congestion while adhering to integer flow
constraints on paths. This restriction makes ITAP an NP-hard problem. While
conventional TAP prioritizes repulsive interactions to minimize congestion,
this work also explores the case of attractive interactions, related to
minimizing the number of occupied edges.

We present and evaluate multiple algorithms to address ITAP, including the
message passing algorithm of [54], a greedy approach, simulated annealing, and
relaxation of ITAP to TAP. Inspired by studies of random ensembles in the
large-size limit in statistical physics, comparisons between these algorithms
are conducted on large sparse random regular graphs with a random set of
origin-destination pairs. Our results indicate that while the simplest greedy
algorithm performs competitively in the repulsive scenario, in the attractive
case the message-passing-based algorithm and simulated annealing demonstrate
superiority. We then investigate the relationship between TAP and ITAP in the
repulsive case. We find that, as the number of paths increases, the solution
of TAP converges toward that of ITAP, and we investigate the speed of this
convergence.

Depending on the number of paths, our analysis leads us to identify two
scaling regimes: in one the average flow per edge is of order one, and in
another the number of paths scales quadratically with the size of the graph,
in which case the continuous relaxation solves the integer problem closely.

††preprint: APS/123-QED

##  I Introduction

###  I.1 Motivation of the study

Ranging from congestion problems in traffic networks to routing problems for
data packages over the internet, path optimization is ubiquitous in real-world
problems. The traffic assignment problem (TAP) is one of the most prominent
routing problems, with a long history dating back to the 50’ [50]. In this
work, we consider the integer traffic assignment problem (ITAP), a variant of
TAP where the flows on edges are constrained to be integers. We are given a
graph representing the roads in a city and a list of origin-destination (OD)
pairs, each one representing a commuter that must move from one point to the
other. For each OD pair we must find a path connecting the origin node to the
destination node, while minimizing some congestion measure on the edges of the
graph. Intuitively we want to distribute the paths (commuters) so that no edge
(road) becomes too congested. While in TAP, for each OD pair, we are allowed
to split the flow across several paths each carrying a fractional amount of
the flow, in ITAP, each path must carry an integer amount of flow, that is, an
integer number of commuters. This restriction makes ITAP an NP-hard problem.
From a physics standpoint one can see this system as a system of variable-
length interacting polymers, where each polymer is clamped at its endpoints
and is either attracted to or repelled from others [13]. In TAP, the
interaction between paths is commonly considered repulsive, the idea being
that we want to minimize the congestion. In our work, we also consider the
attractive case, where it is advantageous to merge paths as much as possible.
This setting is motivated by network-building games in which one must build a
graph that connects the OD pairs efficiently while minimizing the number of
links used [43]. Think, for example, of a city that has to build road
infrastructure that reaches every house efficiently while minimizing the
length of constructed roads. In such formulation one clearly incentivizes
paths to share edges, that is, avoiding building a road between every pair of
houses. In another example, one aims to route internet traffic through as few
nodes as possible so as to turn off the rest and save energy [23].

Our study aims to provide a systematic comparison of various algorithms for
both the repulsive and attractive versions of the problem. We perform this
study on random instances of the problem where the underlying graph is sparse
and random with a varying average degree, the origin-destination pairs are
selected uniformly at random among all the nodes, and we vary the number of
the paths considered. Further, we vary the size of the graph while keeping the
average degree constant. Considering such random benchmarks is common in the
field of statistical physics of disordered systems where physical systems [32]
are often modelled using a random realization of the disorder (here, the graph
and the origin-destination pairs) and increasing the size of the system.
Studying these random ensembles of instances then aims to identify properties
that are prone to be universal in the sense that classes of real-world
instances will exhibit analogous properties. The study of random instances of
the ITAP problem was initiated in [54], and here we extend it by examining
more algorithms.

From an algorithmic point of view, we are particularly interested in the
performance gap between simple greedy heuristics and more sophisticated
algorithms, such as the message passing algorithm from [54] or simulated
annealing, that aim to take into account the interactions between the paths in
a more detailed manner. Another very common manner to solve an integer
optimization problem is to relax it to a version with continuous variables and
investigate the performance gap between the relaxed and integer problems.

###  I.2 Problem definition

In this section, we set the notation and define ITAP. Let $G=(V,E)$ be an
undirected graph. This is the graph over which paths will be routed. Let
$N=|V|$ be the number of nodes, which we index with numbers in
$[N]\coloneqq\\{1,2,\dots,N\\}$. The origin-destination pairs
$\\{(s^{\mu},t^{\mu})\\}_{\mu\in[M]}$ are also given. A path $\mu$ goes from
its origin node $s^{\mu}$ to its destination node $t^{\mu}$. There are a total
of $M$ paths, $\mu=1\dots M$, routed over the network. We name $\pi^{\mu}$ the
path joining $s^{\mu}$ and $t^{\mu}$. Each path $\pi^{\mu}$ is represented as
a sequence of adjacent nodes in $G$:
$\pi^{\mu}=(\pi^{\mu}(1)=s^{\mu},\pi^{\mu}(2),\dots,\pi^{\mu}(L(\pi))=t^{\mu})$,
where $L(\pi^{\mu})$ is the number of nodes in the path. Further, we use the
expression $e\in\pi$ to indicate that $\pi$ traverses edge $e$. The flow
$I_{e}$ through an edge $e\in E$ is the number of paths that pass through $e$.
Mathematically we can write this as

$I_{e}\coloneqq\sum_{\mu\in[M]}\mathbb{I}[e\in\pi^{\mu}]\,.$

We aim to minimize the following energy function

$H(\pmb{I})=\sum_{e\in E}\phi(I_{e}),$ (1)

We refer to $\phi:\operatorname{\mathbb{R}}\mapsto[0,\infty)$ as the
nonlinearity. $\phi(I_{e})$ represents the total cost of having flow $I_{e}$
on edge $e$. We always assume $phi$ to be an increasing function. The cost is
equally distributed among paths that pass through the edge: each one incurs in
a cost $\phi(I_{e})/I_{e}$.

When $\phi$ is convex, minimizing $H$ will have a repulsive effect on the
paths (this is the usual traffic scenario in TAP). If instead $\phi$ is
concave, minimizing $H$ has an attractive effect on the paths, leading many
paths to share the same edge. To see this, suppose $\phi(0)=0$, and consider
the cost of two paths (each carrying one unit of flow) sharing an edge vs.
passing through different edges. In the first case the contribution to the
energy is $\phi(2)$, while in the second case it is $2\phi(1)$. If $\phi$ is
convex then $\phi(2)\geq 2\phi(1)$ therefore it is convenient to not share the
edge. Vice versa if $\phi$ is concave the configuration with the paths
overlapping is energetically convenient.

Several algorithms that we use to minimize the energy (1) rely on defining a
Gibbs measure $\propto e^{-\beta H}$ over all possible paths and then studying
its properties to find minima of $H$. More precisely, the measure is defined
as

$P_{\beta}\left(\\{\pi^{\mu}\\}_{\mu\in[M]}|\\{s^{\mu},t^{\mu}\\}_{\mu\in[M]}\right)=\frac{e^{-\beta
H(\pmb{I})}}{Z\left(\beta,\\{s^{\mu},t^{\mu}\\}_{\mu\in[M]}\right)}\left[\prod_{\mu\in[M]}\operatorname{\mathbb{I}}[\pi^{\mu}(1)=s^{\mu}]\operatorname{\mathbb{I}}[\pi^{\mu}(L(\pi^{\mu}))=t^{\mu}]\operatorname{\mathbb{I}}[\pi^{\mu}\text{
is a SAW}]\right].$ (2)

$\beta$ is an inverse temperature parameter: for $\beta=0$ the distribution is
uniform across all self-avoiding walks (SAWs) with constrained endpoints,
while for $\beta\to\infty$ it is supported on the global minima of $H$. The
three constraints ensure respectively that the path has the correct origin and
destination and that it is self-avoiding, i.e., it does not visit the same
node more than once. Finally
$Z\left(\beta,\\{s^{\mu},t^{\mu}\\}_{\mu\in[M]}\right)$ is a normalization
constant.

The joint distribution of all the paths is very complex, and hard to work with
in practice. Therefore many algorithms rely instead on considering one path
(say $\pi^{\nu}$) at a time and approximating in different ways the
conditional distribution of the path when all the others are fixed. We now
derive the form of the conditional probability
$P_{\beta}\left(\pi^{\nu}|\\{\pi^{\mu}\\}_{\mu\in[M]\backslash\nu},\\{s^{\mu},t^{\mu}\\}_{\mu\in[M]}\right)$.

We first isolate the contribution of path $\nu$ in $H$. Define
$I_{e}^{\backslash\nu}$ to be the flow on edge $e$ in absence of path $\nu$.
Mathematically we have
$I_{e}^{\backslash\nu}=\sum_{\mu=1,\mu\neq\nu}^{M}\operatorname{\mathbb{I}}[e\in\pi^{\mu}]$.
We can write $H$ in the following form.

$\displaystyle H(\pmb{I})$ $\displaystyle=\sum_{e\in E}\phi(I_{e})=\sum_{e\in
E}\phi(I^{\backslash\nu}+\mathbb{I}[e\in\pi^{\nu}])=\sum_{e\in
E}\phi(I_{e}^{\backslash\nu})+\sum_{e\in
E}\mathbb{I}[e\in\pi^{\nu}]\left(\phi(I_{e}^{\backslash\nu}+1)-\phi(I_{e}^{\backslash\nu})\right)=$
(3) $\displaystyle=H(\pmb{I}^{\backslash\nu})+\Delta
H(\pi^{\nu};\pmb{I}^{\backslash\nu}),$ (4)

with $\Delta H(\pi^{\nu};\pmb{I}^{\backslash\nu})=\sum_{e\in
E}\mathbb{I}[e\in\pi^{\nu}]\left(\phi(I_{e}^{\backslash\nu}+1)-\phi(I_{e}^{\backslash\nu})\right)$.
Notice that the only dependence of $\pi^{\nu}$ is contained in the term
$\Delta H(\pi^{\nu};\pmb{I}^{\backslash\nu})$. This implies that the
conditional distribution has the form

$\displaystyle
P_{\beta}\left(\pi^{\nu}|\\{\pi^{\mu}\\}_{\mu\in[M]\backslash\nu},\\{s^{\mu},t^{\mu}\\}_{\mu\in[M]}\right)=\frac{e^{-\beta\Delta
H(\pi^{\nu};\pmb{I}^{\backslash\nu})}}{Z^{\nu}\left(\beta,\pmb{I}^{\backslash\nu},s^{\nu},t^{\nu}\right)}\operatorname{\mathbb{I}}[\pi^{\nu}(1)=s^{\nu}]\operatorname{\mathbb{I}}[\pi^{\nu}(L(\pi^{\nu}))=t^{\nu}]\operatorname{\mathbb{I}}[\pi^{\nu}\text{
is a SAW}],$ (5)

where $Z^{\nu}\left(\beta,\pmb{I}^{\backslash\nu},s^{\nu},t^{\nu}\right)$ is
again a normalization constant. In words, $\Delta
H(\pi^{\nu};\pmb{I}^{\backslash\nu})$ is the energy landscape seen by path
$\pi^{\nu}$ when all other paths are fixed.

###  I.3 Relaxation of the ITAP

In this section, we first explain what relaxation procedure transforms ITAP
into TAP and then illustrate some properties of TAP that we will exploit in
the numerical analysis. We start by defining the demand matrix
$D\in\operatorname{\mathbb{R}}^{N\times N}$, with elements
$D_{xy}=\sum_{\mu\in[M]}\operatorname{\mathbb{I}}[s^{\mu}=x]\operatorname{\mathbb{I}}[t^{\mu}=y]$.
In words, $D_{xy}$ counts the number of origin-destination (OD) pairs that are
equal to $(x,y)$. For simplicity, we say that each path corresponds to one
unit of flow. $D_{xy}$ then tells us how many units of flow must be routed
from $x$ to $y$. An instance of TAP or ITAP is uniquely defined by giving the
graph $G$ and the demand matrix $D$. For each pair of nodes $(x,y)$ we define
$\Pi_{xy}$ to be the set of paths going from $x$ to $y$. We indicate elements
of $\Pi_{xy}$ with $\pi_{xy}^{a}$, $a\in\\{1,\dots,|\Pi_{xy}|\\}$. Similarly
the flow carried by $\pi_{xy}^{a}$ is $h_{xy}^{a}$. We will refer to
$\\{h_{xy}^{a}\\}$ as the path flows, and to $\\{I_{e}\\}_{e\in E}$ as the
edge flows. We now give the definition of TAP and ITAP.

$\displaystyle\underset{}{\text{minimize}}\qquad H(\pmb{I})$ (6) subject to
$\displaystyle I_{e}=\sum_{(x,y)\in(V\times
V)}\sum_{a\in[|\Pi_{xy}|]}h^{a}_{xy}\operatorname{\mathbb{I}}[e\in\pi^{a}_{xy}]\quad\forall
e\in E$
$\displaystyle\sum_{a\in[|\Pi_{xy}|]}h_{xy}^{a}=D_{xy}\quad\forall(x,y)\in(V\times
V)$ $\displaystyle h_{xy}^{a}\in\\{0,1,\dots\\}\quad\forall(x,y)\in(V\times
V),\;\forall a\in[|\Pi_{xy}|]\textbf{ ITAP constraint}$ $\displaystyle
h_{xy}^{a}\in[0,\infty)\quad\forall(x,y)\in(V\times V),\;\forall
a\in[|\Pi_{xy}|]\textbf{ TAP constraint}$

The minimization is over all path flows $\\{h_{xy}^{a}\\}_{(x,y)\in V\times
V,a\in[|\Pi_{xy}|]}$.

ITAP and TAP are identical except for the last constraint. In ITAP the flow on
each path is constrained to be integer, while in TAP it can be any positive
real number. Define $R_{xy}\coloneqq|\\{a\in[|\Pi_{xy}|],\text{ s.t.
}h_{xy}^{a}>0\\}|$. This is the number of paths from $x$ to $y$ that carry a
strictly positive flow. From the constraints, one sees that $R_{xy}\leq
D_{xy}$ in ITAP, while it is potentially unbounded111that is bounded only by
$|\Pi_{xy}|$ which generally scales super-exponentially in $N$ in TAP, as one
can split the flow across all the paths. In the following we will use the term
’at optimality’ to mean that the paths and flow configuration corresponds to a
global minimum of $H$.

####  I.3.1 TAP is easy in the repulsive case

Let us first consider the case in which $\phi$ is convex. Below we list a few
properties of TAP that make it a useful relaxation. See [35, 7] for proofs.

######  Theorem I.1.

Suppose $\phi:\operatorname{\mathbb{R}}\mapsto\operatorname{\mathbb{R}}$ is
convex. Then:

  1. 1.

TAP is a convex optimization problem

  2. 2.

There exist polynomial (in $N,M$) time algorithms that solve TAP (see for
example [38])

  3. 3.

At optimality, the path flows $\\{h_{xy}^{a}\\}_{(x,y)\in V\times
V,a\in[|\Pi_{xy}|]}$ are not in general unique

  4. 4.

At optimality, the edge flows $\\{I_{e}\\}_{e\in E}$ are unique

The first two properties guarantee that TAP can be solved efficiently. The
third property establishes that there can be several path flows that all give
the minimal $H$. Last, the edge flows are instead unique when $H$ is minimal.

####  I.3.2 TAP=ITAP in the attractive case

In the previous section, we saw how ITAP relaxes to a convex optimization
problem when $\phi$ is convex (repulsive case). This is no longer the case
when $\phi$ is concave, as $H$ is not a convex function anymore and can, in
fact, have multiple local minima. As a result, TAP is NP-hard for non-convex
$\phi$.

Additionally, we now show that if $\phi$ is concave, then the solutions of TAP
and ITAP coincide, meaning that even in the relaxed problem the optimal path
flows are always integer. First we define a local minimum in the problem TAP.
We do so at the level of edge flows, since there can be several path flows
corresponding to one edge flow.

We say $\pmb{I}$ is a TAP (ITAP) feasible flow, if it satisfies the three TAP
(ITAP) constraints in (6) for some choice of the path flows. Informally, a
feasible edge flow, is one that is induced by some choice of path flows that
respect the demand and positivity (and integrality) constraints.

######  Definition I.2 (TAP local minimum).

A TAP feasible edge flow $\pmb{I}$ is a local minimum of $H$ if for every
feasible edge flow $\pmb{\tilde{I}}$, there exists an $\epsilon>0$ such that
$H(\pmb{I})\leq H(\epsilon\pmb{\tilde{I}}+(1-\epsilon)\pmb{I})$.

Notice that $\epsilon\pmb{\tilde{I}}+(1-\epsilon)\pmb{I}$ is feasible if both
$\pmb{\tilde{I}},\pmb{I}$ are; so our definition says that at a minimum,
moving in any feasible direction from $\pmb{I}$ results in an energy increase.
We now give a characterization of local minima of $H$ when $\phi$ is concave.

######  Proposition I.3.

Consider an instance $G,D$ of TAP, with $\phi$ concave. Then for every pair of
nodes $x,y$ with $D_{xy}>0$, at any local minimum of $H$ all the flow going
from $x$ to $y$ must be routed through a single path. Put differently
$R_{xy}=1$.

###### Proof.

Let $\pmb{I}^{\backslash xy}$ be a feasible flow for the problem with
$D_{xy}=0$ and the rest of the demand matrix unchanged. Put differently,
$\pmb{I}^{\backslash xy}$ is a flow that results from routing all the flow
except for that going from $x$ to $y$. Let $\pi^{0},\pi^{1}$ be two distinct
paths joining $x$ and $y$. Consider the following two flows

$I^{0}_{e}=I^{\backslash xy}_{e}+D_{xy}\mathbb{I}[e\in\pi^{0}],\quad
I^{1}_{e}=I^{\backslash xy}_{e}+D_{xy}\mathbb{I}[e\in\pi^{1}],\qquad e\in E$
(7)

Both $\pmb{I}^{0}$ and $\pmb{I}^{1}$ are feasible flows for $(G,D)$. All the
flows $\pmb{I}^{\lambda}$ in the interpolating sequence
$\pmb{I}^{\lambda}=\lambda\pmb{I}^{1}+(1-\lambda)\pmb{I}^{0},\;\lambda\in[0,1]$
are also feasible. In $\pmb{I}^{\lambda}$, $\lambda D_{xy}$ units of flow are
sent through $\pi^{1}$ and the rest is sent through $\pi^{0}$, thus the flow
is split among two paths. Using concavity we have that $H(\pmb{I}^{\lambda})$
is concave as a function of $\lambda\in[0,1]$. This implies that it can have
local minima only for $\lambda\in\\{0,1\\}$ corresponding to all of the flow
being routed through a single path. ∎

Proposition I.3 immediately implies that also on global optima all the flow
between any two nodes gets routed through a single path. As a consequence,
since in our case the matrix $D$ has integer entries, a solution to TAP also
satisfies the ITAP constraint in (6) and is therefore integer.

###  I.4 Related Literature

When the objective function is convex, ITAP is related to TAP, a problem
consisting of routing paths (users) over a network in a way that minimizes the
average time taken to reach the destination. The study of TAP dates back to
the seminal work of Wardrop [50] who formulated a mathematical model of
traffic assignment and introduced the concepts of user equilibrium and system
optimum. These correspond respectively to the case where each path
egoistically chooses the route that minimizes _its own_ travel time, and the
case where the route taken by each path is mandated so to minimize the total
travel time of all paths. Wardrop’s principles state the mathematical
conditions for user equilibrium to be satisfied and establish that both
problems can be formulated as optimization problems but with different
objectives. Subsequently, [4, 12] studied TAP as a convex optimization
problem. The first algorithm employed to solve it is the Frank-Wolfe method
[21], followed by increasingly efficient algorithms [2, 14, 51, 52, 33]. See
[35, 7] for a comprehensive review of TAP, and [38] for a comparison of all of
the common optimization algorithms to solve it.

The integer version (ITAP) of the problem also received some interest, for
example in [16] Wardrop’s principles and the optimization formulation were
extended to ITAP. ITAP is also related to the multicommodity flow problem
(MCF) and the integer multicommodity flow problem (IMCF), in which,
analogously to ITAP, one cannot split the flow among several paths [3, 17].
The difference with TAP is that in MCF and IMCF the interaction between paths
is given exclusively by a capacity constraint, i.e., the number of paths
crossing an edge cannot exceed its capacity [44, 45]. Variants of ITAP, such
as the routing and wavelength assignment problem, arise in optics and
telecommunications, where one must establish connections between pairs of
nodes in a fibre optics network [34, 30, 10].

TAP traditionally only considers the case in which the paths repel each other.
From a traffic perspective this makes sense, since we do not want too many
users on the same road. The ITAP problem we consider also allows for
attractive interactions between paths. This line of work has been studied in
the setting of network building games. For a review of congestion and network
building games see [43]. The game theory analysis is more concerned with
quantifying the price of anarchy in games, that is the inefficiency of user
(Nash) equilibria in the game where each user acts egoistically, with respect
to the system optimum [42, 39, 28, 1]. While finding Nash equilibria is a
seemingly different problem, Rosenthal proved that one can build a potential
function with the property that every minimum of the potential corresponds to
a Nash equilibrium [41]. Therefore finding a Nash equilibrium consists of
finding a local minimum of the ITAP objective function, for an appropriate
choice on the nonlinearity. The game theory approach can also be useful in
designing algorithms for ITAP, for example the greedy algorithm is inspired to
the best response dynamics in congestion games [18, 19]. Also the logit
dynamics for congestion games, which has been studied in [6, 25] has some
similarity with the simulated annealing algorithm we propose.

The statistical physics approach to ITAP has been initiated in [54]. In this
work the authors study ITAP on sparse random regular graphs (RRGs). They
identify one of the relevant parameters of the problem (i.e. they find the
$N,M$ dependence of $\rho$ in (21)) and propose the conditional belief
propagation (CBP) algorithm (our name). CBP is a belief propagation [36, 37,
47] based algorithm derived using the replica method [32, 9]. See [31, 55, 5]
for a comprehensive treatment of message passing and replica techniques. In
our analysis, we will borrow the RRG setting and provide a shorter and more
transparent derivation of CBP. The simulated annealing [24] algorithm we use
updates the paths by proposing at each time step a randomly generated SAW with
fixed endpoints. Several Markov chain Monte Carlo (MCMC) methods for sampling
unconstrained SAWs have been developed [49, 29, 46]. Knuth [27] devised an
algorithm to sample SAWs with fixed endpoints on a grid, but it samples from a
distribution that is nonuniform and difficult to characterize. We devise a
variant of Knuth’s algorithm that allows us to sample SAWs from a range of
distributions on a graph. In particular, by tuning a temperature parameter, we
can go from sampling small fluctuations around the shortest path, to sampling
SAWs uniformly.

To conclude, the integrality of the solutions to the relaxed problem, has been
studied in a game theoretical setting in [26]. However the integrality results
hold under quite restrictive conditions, which are not satisfied in our case.

###  I.5 Our contribution

Below we list our main contributions.

  * •

We propose two novel algorithms to solve ITAP. The first, RITAP, is a simple
algorithm based on relaxing ITAP to TAP, solving TAP, and then projecting the
solution back on the integer constraints. Its strength lies in the
computational efficiency when $M$ is large and in the possibility to study it
theoretically. The second algorithm is based on simulated annealing, a
standard technique to solve non-convex optimization problems. The advantage of
this method originates from the design of the underlying MCMC, which is based
on a novel algorithm that samples SAWs with fixed endpoints. Numerical results
show that simulated annealing is often the most effective method, especially
in the attractive case. Finally, the SAW sampler is also of independent
interest to the community studying techniques to generate SAWs randomly.

  * •

We present a greatly simplified (with respect to [54]) derivation of CBP that
emphasizes the relation with the greedy algorithm.

  * •

We compare the performance of the considered algorithms and study the
properties of ITAP in a random instance setting. The random instances of ITAP
consist of sparse RRGs, uniformly random OD pairs and nonlinearities of the
form $\phi(x)=x^{\gamma}$. The comparison is carried out in the regime where
the flow per edge is of order one. We find that in the repulsive case, most
algorithms including the simple greedy one, exhibit comparable performance. In
the attractive case the CBP and simulated annealing algorithms have an edge
over the others.

  * •

We study the relation between TAP and ITAP, as the number of paths grows. In
the random instance setting, we observe that for convex $\phi$, when $M\gg
N^{2}$, TAP and ITAP reach identical energies, moreover TAP solutions almost
satisfy the integrality constraints. This implies that ITAP is easy to solve
in this limit. We conjecture that this phenomenology holds for convex
nonlinearities, provided that the distribution of OD pairs is sufficiently
uniform. A formal statement of this conjecture is presented in III.2.

  * •

In random instances, as we change the number of paths $M$ on the graph, we
identify and study two intrinsic regimes of the problem, defined by the
parameters $\rho$ and $\eta$ (see respectively (21),(22)). The first, partly
identified in [54], arises when $M$ is such that the average flow per edge is
of order one. In this regime, the effective interaction between paths arises
and is usually the strongest, in the sense that the optimized energy found by
the algorithms differs the most from the energy of the shortest paths, which
are the optimal paths in the absence of interactions. The second regime is
encountered when $M=\Theta(N^{2})$, and thus the flow per edge is very large.
When $M\approx N^{2}$ we observe the onset of an asymptotic behavior that
continues for higher $M$. It is precisely in this scaling that the TAP
solution (in the convex case) approaches the ITAP one. The presence of the two
regimes is verified also in the case of graphs topologies taken from real-
world traffic networks.

##  II Algorithms

Below we state the algorithms we studied to solve ITAP. We release the code of
each algorithm at github.com/SPOC-group/algorithms_integer_traffic_assignment
222To request the code of the CBP algorithm as implemented in [54], please
write to Bill Chi Ho Yeung. and the code and data to reproduce our results at
github.com/SPOC-group/numerics_integer_traffic_assignment.

###  II.1 Greedy Algorithm

The greedy algorithm is perhaps the simplest algorithm that can be devised to
minimize the energy (1). First the paths are arbitrarily initialized. At every
time step, a path $\nu\in[M]$ is chosen. Then the path is updated to the value
that minimizes $H$, leaving all other paths fixed. Since we want to minimize
$H$ with respect to $\pi^{\nu}$, by the decomposition (3) this is equivalent
to minimizing $\Delta H(\pi^{\nu};\pmb{I}^{\backslash\nu})$. In turn this
coordinate minimization is equivalent to finding the shortest path from
$s^{\nu}$ to $t^{\nu}$ in a graph where the weight of edge $e$ is
$\left(\phi(I_{e}^{\backslash\nu}+1)-\phi(I_{e}^{\backslash\nu})\right)$. This
is easily doable for example via Dijkstra algorithm [15, 11]. The greedy
algorithm proceeds by iteratively selecting an index $\nu$, then setting
$\pi^{\nu}$ to be the shortest path in the aforementioned graph, until
$H(\pmb{I})$ stops decreasing. In appendix A we provide the pseudocode of an
efficient implementation of this algorithm. One can see the greedy algorithm
as a method that samples from the conditional distribution (5) after setting
$\beta=\infty$. In the infinite $\beta$ limit the conditional distribution
becomes in fact concentrated on the shortest paths.

###  II.2 Conditional Belief Propagation (CBP)

Belief propagation (BP) is a message passing algorithm that allows to
approximate the marginals of a given probability measure. It does so by
exchanging messages between the nodes of the factor graph associated to the
measure. One can attempt to apply BP to the joint Gibbs measure of all paths
(2). In this case BP would allow to approximate the probability that a path
passes through a given node, from which the paths can then be reconstructed.
Already in [54], it was recognized that this approach is unfeasible: the fact
that potentially all the paths interact on each edge results in an exponential
time complexity in $M$. To circumvent this problem the CBP algorithm was
proposed in [54]. The idea is to apply BP to the conditional distribution (5)
of a single path, while fixing all the others, and then iterating through
paths. In this section we sketch CBP algorithm; Appendix B contains the
details of its derivation. We indicate with $i\to j$ the directed edge from
node $i$ to $j$. In the graph $G$, every undirected edge is treated as a pair
of opposite directed edges. Path $\pi^{\nu}$ is treated as a directed sequence
of edges starting at $s^{\mu}$ and terminating at $t^{\nu}$. We first
introduce the binary variables $\sigma_{i\to j}^{\nu}\coloneqq\mathbb{I}[i\to
j\in\pi^{\nu}]$, which indicate the directed edges through which $\pi^{\nu}$
passes. For an undirected edge $ij$ we define
$H_{ij}^{\nu}\coloneqq\phi(I_{ij}^{\backslash\nu}+1)-\phi(I_{ij}^{\backslash\nu})$333
Recall that we do not separate flow on an edge by direction. $I_{e}$ is the
total bidirectional flow on the undirected edge $e$. This way we can write the
conditional (5) as

$\displaystyle P\left(\\{\sigma^{\nu}_{i\to j}\\}_{(ij)\in
E}|\\{H^{\nu}_{ij}\\}_{(ij)\in
E},s^{\nu},t^{\nu}\right)\propto\left[\prod_{i\in[N]\backslash
s^{\nu},t^{\nu}}\mathbb{I}\left[\sum_{\begin{subarray}{c}k,l\in\partial i\\\
k\neq l\end{subarray}}\sigma^{\nu}_{k\to i}+\sigma^{\nu}_{i\to l}=(0\text{ OR
}2)\right]\exp{\bigg{(}-\beta\sum_{k\in\partial
i}H^{\nu}_{ik}\sigma^{\nu}_{i\rightarrow k}\bigg{)}}\right]\times$
$\displaystyle\times\mathbb{I}\left[\sum_{k\in\partial
s^{\nu}}\sigma^{\nu}_{s^{\nu}\to k}=1\right]\mathbb{I}\left[\sum_{k\in\partial
s^{\nu}}\sigma^{\nu}_{k\to s^{\nu}}=0\right]\mathbb{I}\left[\sum_{k\in\partial
t^{\nu}}\sigma^{\nu}_{k\to t^{\nu}}=1\right]\mathbb{I}\left[\sum_{k\in\partial
t^{\nu}}\sigma^{\nu}_{t^{\nu}\to k}=0\right]\prod_{k\in\partial
s^{\nu}}\bigg{(}\exp{(-\beta
H^{\nu}_{s^{\nu}k}\sigma^{\nu}_{s^{\nu}\rightarrow k})}\bigg{)},$ (8)

where $\partial i$ denotes the set of neighbours of node $i$. The constraint
in the first line ensures that the path passes at most once in each node. The
constraints in the second line take care of the origin and destination nodes,
which represent special cases. From this probability distribution, BP
equations are derived following the procedure described in [31, 55]. We take
the limit $\beta\rightarrow\infty$ and obtain the following update equations
for two sets of messages $a^{\nu}_{i\rightarrow j}$ and $b^{\nu}_{i\rightarrow
j}$.

$\displaystyle a^{\nu}_{i\to j}$
$\displaystyle=\left\\{\begin{array}[]{lr}\min_{k\in\partial i\backslash
j}[a^{\nu}_{k\rightarrow
i}]+H^{\nu}_{ij}-\min\bigg{[}0,\min_{\begin{subarray}{c}k,l\in\partial
i\backslash j\\\ k\neq l\end{subarray}}[a^{\nu}_{k\rightarrow
i}+b^{\nu}_{l\rightarrow i}]\bigg{]},&\text{if
}i\not\in\\{s^{\nu},e^{\nu}\\}\\\ \infty&\text{if }i=e^{\nu}\\\ \\\
H^{\nu}_{s^{\nu}j}-\min_{k\in\partial s^{\nu}\backslash
j}[b^{\nu}_{k\rightarrow s^{\nu}}],&\text{if }i=s^{\nu}\end{array}\right.$
(13) $\displaystyle b^{\nu}_{i\to j}$
$\displaystyle=\left\\{\begin{array}[]{lr}\min_{k\in\partial i\backslash
j}[b^{\nu}_{k\rightarrow
i}]+H^{\nu}_{ji}-\min\bigg{[}0,\min_{\begin{subarray}{c}k,l\in\partial
i\backslash j\\\ k\neq l\end{subarray}}[a^{\nu}_{k\rightarrow
i}+b^{\nu}_{l\rightarrow i}]\bigg{]},&\text{if
}i\not\in\\{s^{\nu},e^{\nu}\\}\\\ \infty,&\text{if }i=s^{\nu}\\\ \\\
H^{\nu}_{je^{\nu}}-\min_{k\in\partial e^{\nu}\backslash
j}[a^{\nu}_{k\rightarrow e^{\nu}}],&\text{if }i=e^{\nu}\end{array}\right.$
(18)

From the messages the flow is computed as
$I^{\backslash\nu}_{ij}=\sum_{\mu\in[M]\backslash\nu}(\sigma^{\mu}_{i\rightarrow
j}+\sigma^{\mu}_{j\rightarrow i})$, with

$\displaystyle\sigma^{\nu}_{i\rightarrow j}=$
$\displaystyle\operatorname{\mathbb{I}}[i=s^{\nu}]\Theta\bigg{(}\min_{l\in\partial
j\backslash i}(b^{\nu}_{l\rightarrow j})-b^{\nu}_{i\rightarrow j}\bigg{)}$
$\displaystyle+(1-\operatorname{\mathbb{I}}[i=s^{\nu}])\Theta\bigg{(}\min{\bigg{(}0,\min_{\begin{subarray}{c}l\in\partial
j\backslash i\end{subarray}}{(a^{\nu}_{l\rightarrow
j}+\min_{\begin{subarray}{c}k\in\partial j\backslash
l\end{subarray}}b^{\nu}_{k\rightarrow j})}\bigg{)}}-(a^{\nu}_{i\rightarrow
j}+\min_{\begin{subarray}{c}k\in\partial j\backslash
i\end{subarray}}b^{\nu}_{k\rightarrow j})\bigg{)}.$ (19)

In [54] the quantity
$H^{\nu}_{ij}=\phi(I_{ij}^{\backslash\nu}+1)-\phi(I_{ij}^{\backslash\nu})$ is
Taylor approximated with $\phi^{\prime}(I_{ij}^{\backslash\nu})$. The CBP
algorithm alternates between updating all messages once using (13), (18) and
then updating the flows using (19). This is repeated until convergence is
reached, at which point, the paths are determined by plugging the final values
of the messages into equation (19). For our experiments we use the original
implementation of CBP from [54]. Its pseudocode, presenting some minor
differences with what we derived, is provided in Appendix B.

By picking a different update schedule one can essentially retrieve the greedy
algorithm. Suppose we pick a path $\nu$ and update the messages $\\{a_{i\to
j}^{\nu},b_{i\to j}^{\nu}\\}_{(i\to j)\in E}$ several times until they reach
convergence. Then BP approximates the conditional
$P_{\infty}\left(\pi^{\nu}|\\{\pi^{\mu}\\}_{\mu\in[M]\backslash\nu},\\{s^{\mu},t^{\mu}\\}_{\mu\in[M]}\right)$
which is supported on shortest paths from $s^{\nu}$ to $t^{\nu}$. Therefore,
with this schedule, BP approximates the shortest path computation described in
II.1. If one then updates the variables $\\{\sigma_{i\to j}^{\nu}\\}_{(i\to
j)\in E}$ using (32) and then repeats the whole procedure for other values of
$\nu$ this reproduces the greedy algorithm.

###  II.3 Simulated annealing

Given an energy function $H$, to be minimized, the simulated annealing
technique consists of building a probability measure $\propto e^{-\beta H}$
and then drawing samples from it while increasing $\beta$. The sampling part
is usually accomplished using iterative methods such as MCMC. For not-too-high
values of $\beta$, the noise in the distribution will allow the dynamics to
jump possible barriers in $H$ and explore the landscape. When increasing
$\beta$, the measure will become more and more concentrated on the global
minimum of $H$. In practice $\beta$ is changed with the iteration number $t$,
according to a schedule $\beta(t)$. For ITAP, the simulated annealing aims to
draw samples from the ($t$ dependent) probability measure
$P_{\beta(t)}\left(\\{\pi^{\mu}\\}_{\mu\in[M]}|\\{s^{\mu},t^{\mu}\\}_{\mu\in[M]}\right)$,
where $\beta(t)$ is an increasing function. To accomplish this task we use a
MCMC based on Gibbs sampling. Gibbs sampling [22] is an iterative MCMC
algorithm that samples from the joint distribution (2) by sampling each
variable from its conditional distribution (5) in a sequence. The details of
annealing schedule and the pseudocode for the algorithm are provided in
Appendix D. Here we explain how it works using a mock example with three
paths. To lighten the notation define the joint measure
$P_{\beta}(\pi^{1},\pi^{2},\pi^{3})\coloneqq
P_{\beta}\left(\\{\pi^{\mu}\\}_{\mu\in[3]}|\\{s^{\mu},t^{\mu}\\}_{\mu\in[3]}\right)$,
and the conditional measure $P(\pi^{1}|\pi^{2},\pi^{3})\coloneqq
P_{\beta}\left(\pi^{1}|\\{\pi^{\mu}\\}_{\mu\in[3]\backslash
1},\\{s^{\mu},t^{\mu}\\}_{\mu\in[3]}\right)$. Initially, given the
configuration $\pi^{1}(t),\pi^{2}(t),\pi^{3}(t)$, the algorithm samples
$\pi^{1}(t+1)$ from $P_{\beta(t)}(\pi^{1}|\pi^{2}(t),\pi^{3}(t))$, followed by
sampling $\pi^{2}(t+1)$ from $P_{\beta(t)}(\pi^{2}|\pi^{1}(t+1),\pi^{3}(t))$,
and finally $\pi^{3}(t+1)$ from
$P_{\beta(t)}(\pi^{3}|\pi^{1}(t+1),\pi^{2}(t+1))$. If $\beta(t)=\beta$ were
constant in $t$, then repeating this process, the samples
$(\pi^{1}(t),\pi^{2}(t),\pi^{3}(t))$ would converge to the distribution
$P_{\beta}(\pi^{1},\pi^{2},\pi^{3})$ [8, 40]. In the simulated annealing,
since $\beta(t)$ increases towards infinity, there is no stationary
distribution and the above stochastic process will lead to a local minimum of
$H$.

How do we sample from the conditional distribution (5)? The conditional is a
probability measure over a single self-avoiding path with constrained
endpoints and subject to the Gibbs measure $\propto e^{-\beta\Delta H}$. The
self-avoiding constraint makes this measure hard to sample directly. Instead,
we resort to a second, Metropolis based, MCMC that allows us to sample self-
avoiding paths. The method, thoroughly described in appendix C, consists of
generating SAWs using a proposal distribution which is easy to sample and that
approximates the conditional. Then a Metropolis accept/reject step is carried
out. The proposal distribution is chosen to capture the leading order behavior
of the conditional when $\beta\to\infty$. This means that the sampler is
particularly efficient (low rejection rate) at high $\beta$.

###  II.4 Relaxed ITAP solver (RITAP)

The last algorithm we study is based on relaxing ITAP to TAP, solving TAP, and
finally projecting the solution of TAP to a solution of ITAP. When having to
solve TAP we will assume that we have an algorithm `TAPsolver`$(G,D)$, that
given a graph $G=(V,E)$ and a demand matrix $D$, returns a collection of sets
$\operatorname{\mathcal{P}}=\\{\operatorname{\mathcal{P}}_{xy},(x,y)\in(V\times
V)\\}$. Every set $\operatorname{\mathcal{P}}_{xy}$ contains the paths from
$x$ to $y$ that at a local minimum (in the space of edge flows, see I.2) of
TAP carry some positive flow, as well as the flows carried by each path.
Mathematically we have

$\operatorname{\mathcal{P}}_{xy}=\left\\{(\pi^{a}_{xy},h^{a}_{xy}),\;(x,y)\in(V\times
V),\,a\in\\{1,\dots,R_{xy}\\}\right\\}.$ (20)

We indicate with $\operatorname{\mathcal{P}}_{xy}.\pi$ the set of paths
$\\{\pi_{xy}^{a}\\}_{a=1}^{R_{xy}}$, and with
$\operatorname{\mathcal{P}}_{xy}.h$ the set of path flows
$\\{h_{xy}^{a}\\}_{a=1}^{R_{xy}}$. $R_{xy}$ is the number of paths between $x$
and $y$ carrying a positive flow in the solution. In the repulsive case,
thanks to convexity, there is only one (global) minimum in the space of edge
flows. Therefore `TAPsolver` outputs some global optimal path flows
corresponding to the unique optimal edge flows.444Recall that while the global
minimum is unique in the space of edge flows, it is not in the space of path
flows. The particular path flows outputted depend on the implementation of
TAPsolver In contrast TAP is NP-hard in the attractive case; the particular
local minimum reached by `TAPsolver` depends on the implementation. Our
implementation is based on the Frank-Wolfe algorithm. This is neither the sole
nor the fastest algorithm to solve TAP [38], however we use it due to its
simplicity.

To obtain a valid solution to the integer problem we devise a procedure to
project the relaxed solution to an integer one. In the attractive case,
proposition I.3 guarantees that `TAPsolver` already outputs integer paths,
thus the projection step is superfluous. In the repulsive case we use the
following procedure. Looping over $\mu\in[M]$ we consider
$\operatorname{\mathcal{P}}_{s^{\mu}t^{\mu}}$. We set $\pi^{\mu}$ (the $\mu$th
integer path) to be the path in $\operatorname{\mathcal{P}}_{s^{\mu}t^{\mu}}$
that carries the most flow: mathematically
$a^{\star}=\operatorname{\arg\max}_{a}h^{a}_{s^{\mu}t^{\mu}}$ and
$\pi^{\mu}=\pi_{xy}^{a^{\star}}$. After assigning $\pi^{\mu}$ we update
$h^{a^{\star}}_{s^{\mu}t^{\mu}}\leftarrow h^{a^{\star}}_{s^{\mu}t^{\mu}}-1$,
this way, if there are multiple paths with the same origin and destination,
one takes into account the flow that is already assigned. We report below the
pseudocode describing this procedure. Since it contains a loop over $M$ the
complexity of this algorithm is at least linear in $M$. In appendix E.2 we
present a smarter implementation of RITAP whose time complexity depends on $M$
as $\min(M,N^{2})$, and is therefore more suited to very high $M$. The same
appendix also contains additional details about the Frank-Wolfe algorithm.

Algorithm 1 Relaxed ITAP solver (RITAP)

1:Input:  $G$ graph over which paths are sampled, $D$ demand matrix.

2:Output: $\\{\pi^{\mu}\\}_{\mu\in[M]}$ the paths outputted by the algorithm

3:$\operatorname{\mathcal{P}}\leftarrow\verb|solveTAP|(G,D)$. $\triangleright$
Solving TAP

4:for $\mu=1,\dots,M$ do $\triangleright$ Loop over OD pairs

5:
$a^{\star}\leftarrow\operatorname{\arg\max}_{a}\operatorname{\mathcal{P}}_{s^{\mu}t^{\mu}}.h$
$\triangleright$ Find path with maximum flow

6:
$\pi^{\mu}\leftarrow\operatorname{\mathcal{P}}_{s^{\mu}t^{\mu}}.\pi[a^{\star}]$

7:
$\operatorname{\mathcal{P}}_{s^{\mu}t^{\mu}}.h[a^{\star}]\leftarrow\operatorname{\mathcal{P}}_{s^{\mu}t^{\mu}}.h[a^{\star}]-1$$\triangleright$
Update the flow on $\pi_{s^{\mu}t^{\mu}}^{a^{\star}}$

8:end for

##  III Results on random instances

To test the algorithms and study ITAP, we consider it on undirected sparse
random regular graphs. Origin-destination pairs are sampled i.i.d. from the
uniform distribution over pairs of distinct vertices. Mathematically we sample
$M$ times independently from $P(x,y)=\mathbb{I}[x\neq y]/(N(N-1))$, where
$x,y\in(V\times V)$. We consider nonlinearities of the form
$\phi(x)=x^{\gamma}$, with $\gamma>0$. In the repulsive case we set
$\gamma=2$, while in the attractive case $\gamma=1/2$. We run extensive
numerical experiments with different algorithms to study their ability to
minimize the ITAP energy $H(\pmb{I})$. Additional details about the
experimental settings are given in appendix G.

Given a graph $G=(V,E)$ and a set of OD pairs
$\\{s^{\nu},t^{\nu}\\}_{\nu\in[M]}$, each algorithm outputs a set of paths
$\\{\pi^{\nu}\\}_{\nu\in[M]}$ connecting the respective OD pairs. We name $H$
the energy reached by the algorithm. To establish a baseline, we consider a
very basic algorithm that, for each OD pair $(x,y)$, outputs the shortest path
(in the topological distance on $G$) between $x$ and $y$. We call the
resulting energy $H_{\text{SP}}$.555When the shortest path is not unique we
output one of them; this arbitrariness does not significantly affect our
results. We use the shortest paths as a benchmark since they represent the
optimal routing for the non-interacting case (every path is routed
independently of others), corresponding to $\gamma=1$. In fact in this case
the energy can be written as
$H(\pmb{I})=\sum_{\mu\in[M]}\sum_{e\in\pi^{\mu}}1=\sum_{\mu\in[M]}\text{Length}(\pi^{\mu})$.
Comparing $H$ to $H_{\text{SP}}$ also allows us to gauge the effective
strength of interactions, that is how much gain in energy one achieves by
taking the interaction into account.

###  III.1 Regime with the average flow per edge $O(1)$

For every considered algorithm we plot the fraction of energy saved when
compared to the shortest paths, mathematically expressed as
$1-H/H_{\text{SP}}$. To minimize the energy we employ several algorithms:
Greedy, CBP, simulated annealing, and RITAP. In each curve we fix the degree
$d$ and the number of nodes $N$, and vary the number of OD pairs $M$. We plot
quantities as a function of the effective parameter

$\rho\coloneqq\frac{2M\log N}{Nd\log d}.$ (21)

This will allow us to compare different sizes or degrees. The choice of this
scaling is justified both by the empirical collapse of the curves with
different $N,d$ and by the fact that $\rho$ is approximately the average the
flow per edge therefore it is a proxy for how congested the graph is. To see
this consider the shortest path algorithm. In an RRG on average the shortest
path between two nodes has an average length of approximately $\log N/\log d$
[48]. Then the total flow $\sum_{e\in E}I_{e}$ present on the graph is about
$M\log N/\log d$, while the number of edges is $Nd/2$. This gives a flow per
edge of $\rho$. While the average length of the outputted paths changes with
the algorithm, the difference is small enough not to compromise the validity
of the approximation. The agreement between $\rho$ and the flow per edge is
further verified with numerical experiments in Appendix G. To compare
algorithms, we focus initially on the regime where $\rho$ is of order one,
since for higher $\rho$, most algorithms have similar performances, and for
$\rho\to 0$, the system is non-interacting. Furthermore in most cases, the
quantity $1-H/H_{\text{SP}}$ peaks for $\rho\approx 1$, indicating that this
is where the effective coupling among paths is the greatest. Error bars in all
plots correspond to a $\pm 1$ standard deviation interval around the mean.

####  III.1.1 Attractive paths case ($\gamma=1/2$)

We first consider the case in which paths attract each other, corresponding to
a concave $\phi$. In this case, TAP local optima coincide with ITAP ones, and
both problems are NP-hard. We start by comparing the algorithms for $N=500$
and $d=3$, as portrayed in figure 1.

Figure 1: Relative energy difference between the paths obtained using each
algorithm and the shortest paths. The higher the better. Algorithms used are:
greedy, CBP, simulated annealing, and RITAP. The CBP line is semi-transparent
when CBP converges on less than 80% of the instances. Results are averages
over 200 instances of ITAP, with $d=3,\,N=500,\,\gamma=1/2$.

We immediately remark that RITAP is less effective than other algorithms. This
is not surprising as in the attractive case, the cost function to minimize is
highly non-convex, and RITAP gets easily stuck in a local minimum. Simulated
annealing is observed to be the best-performing algorithm in this case. CBP
instead is affected by convergence issues at high $\rho$. When it converges
for less than 80% of the instances we use a semi-transparent line in the
plots. For high $\rho$, averaging over converged instances, CBP achieves
marginally better performance compared to simulated annealing.

To consolidate our findings, we repeat the experiments for different values on
$N,d$. Figure 2 shows the curves of energy for a fixed $d=3$ and for several
graph sizes.

Figure 2: Relative energy difference between the paths obtained using each
algorithm and the shortest paths. The higher the better. The algorithms used
are: (a) greedy, (b) CBP, (c) simulated annealing, and (d) RITAP. The CBP line
is semi-transparent when CBP converges on less than 80% of the instances. We
fix $\gamma=1/2$, and consider graphs of degree $d=3$ with different sizes
$N\in\\{100,200,400,1000\\}$. Results are averages over 200 instances of ITAP.

Notice that the larger the size, the more energy is saved compared to the
shortest paths. The peaks of curves with different $N$ align, confirming the
validity of $\rho$ as a size-independent measure of congestion. We remark that
for higher $N$, CBP becomes less stable, but at high $\rho$, if it converges,
it is usually the best-performing algorithm. Apart from CBP, simulated
annealing maintains a small but consistent margin over other algorithms. In
the last set of experiments, whose results are depicted in figure 3, we set
$N=200,\gamma=1/2$ and vary $d$.

Figure 3: Relative energy difference between the paths obtained using each
algorithm and the shortest paths. The higher, the better. The algorithms used
are (a) greedy algorithm, (b) CBP, (c) simulated annealing, and (d) RITAP.
Transparent parts of a line correspond to values for which the algorithm
converges only for less than 80% of instances. Results are averages over 200
instances of ITAP, with $N=200,\gamma=1/2$ and $d\in\\{3,6,12\\}$.

We remark that the higher the degree of the RRG, the higher will be the energy
saved when comparing to the shortest paths. We also observe that the
performance gap between RITAP and other algorithms closes at high $d$ and $M$.
For higher $d$, the proportion of instances for which CBP converges changes
haphazardly with $M$. For instance, for $N=200$ and $d=12$, the proportion is
above 80% for the two first points, then decreases to around 60% (transparent
part) and goes back to higher values (opaque part).

####  III.1.2 Repulsive paths case ($\gamma=2$)

We now turn to the repulsive case, where $\phi$ is convex. The relaxed
problem, TAP, in this case is a convex, and hence easy, optimization problem.
Figure 4 shows the dependence of $1-H/H_{\text{SP}}$ on $\rho$ for different
algorithms. The degree and the size of the graphs are fixed to $d=3$ and
$N=500$. We also plot the curve corresponding to the TAP optimal energy, which
constitutes an upper bound to the energy saving that can be achieved by any
algorithm that minimizes the energy under the ITAP constraints.

Figure 4: Relative energy difference between the paths obtained using each
algorithm and the shortest paths. The higher the better. Algorithms used are:
greedy, CBP, simulated annealing, RITAP and, as an upper bound, the TAP
optimal energy. Results are averaged over 200 instances of ITAP, with
$d=3,\,N=500,\,\gamma=2$.

Greedy, CBP, and simulated annealing display approximately the same
performance in terms of saved energy, instead RITAP is slightly less
efficient. Increasing $\rho$ however all algorithms perform comparably and the
gap between saved energy in ITAP and the TAP bound progressively decreases.
Considering very small $M$, all the curves (except for the bound) go to zero.
This is because when only a few paths are present, they are unlikely to share
edges and therefore do not interact. In absence of interactions, the shortest
paths are the global optimal solution. Conversely, the TAP bound follows a
different behavior. For small values of $\rho$, the quantity
$1-H/H_{\text{SP}}$ is close to 0.8. Indeed, since the paths are not
constrained to be integer, the flow is split across several paths, allowing a
considerable advantage of energy compared to any integer configuration.

In figure 5 we vary $N$ and keep $d=3$ fixed. We observe that the size does
not modify the phenomenology observed in 4. Moreover the curves for different
$N$ overlap with each other.

Figure 5: Relative energy difference between the paths obtained using each
algorithm and the shortest paths. The algorithms used are: (a) greedy , (b)
CBP, (c) simulated annealing, and (d) RITAP. The higher the better. Graphs of
degree $d=3$ with different sizes $N\in\\{100,200,400,1000\\}$ are used.
Results shown are averages over 200 realizations with $\gamma=2$.

Next we fix $N=200,\gamma=2$ and vary the degree of the RRG. Figure 6
illustrates the results. Similarly to the $\gamma=1/2$ case, as the degree
increases, the quantity $1-H/H_{\text{SP}}$ reaches larger values, indicating
that the shortest path algorithm is increasingly suboptimal. Additionally, the
simulated annealing gains some slight advantage over other algorithms for
higher $d$. In the same figure we also plot the TAP upper bound; while for
$d=3,6$ the gap between TAP and ITAP efficiencies narrows when increasing
$\rho$, it appears to widen for $d=12$. However, for even higher $\rho$, the
gap decreases and ultimately approaches zero (see figure 10). Results for CBP
have not been reported for $d=6,12$ since respectively for $\rho\geq 0.43$ and
$\rho\geq 0.14$, more than half of the 200 instances fail to converge. Last,
we remark that $\rho$ yields again a good collapse of the curves with
different degrees, by aligning their peaks.

We end this section with a note on the speed of different algorithms. In
general the greedy algorithm stands out as the most computationally efficient
algorithm that achieves good performance over the considered range of $M$.
RITAP is also fast computationally, but not as effective unless $\rho$ is
large. All the other algorithms incur in higher computational costs, posing
challenges to scalability for large sizes.

Figure 6: Relative energy difference between the paths obtained with each
algorithm and the shortest paths. Lower curves are ITAP solutions outputted
by: (a) greedy algorithm, (b) simulated annealing, (c) RITAP. The higher the
better. The upper curves represent the TAP bound. Graphs of size $N=200$ with
$\gamma=2$ and different degrees $d\in\\{3,6,12\\}$ are used. The results
shown are averages over 200 realizations.

###  III.2 Regime of large average flow per edge $M=O(N^{2})$: Integrality of
relaxed solutions

In the previous section we saw that the parameter $\rho$ is the relevant
quantity to consider when describing systems where the average flow per edge
of order one. Here we argue that the problem exhibits a second intrinsic
scaling regime when $M=\Theta(N^{2})$. We then investigate whether TAP
solutions are also ITAP solutions when $M$ is large and $\phi$ is convex.
Defining an effective parameter

$\eta\coloneqq\frac{M}{N(N-1)},$ (22)

we observe that a stable asymptotic behavior arises for $\eta\approx 1$ and
continues for higher $M$.

We have $\eta=\operatorname{\mathbb{E}}_{D}[D_{xy}]$, hence this regime
corresponds to a dense demand matrix. We run RITAP for several different
values of $d,N$ and plot quantities as a function of $\eta$. In this section
we report the results for fixed $d=3,\gamma=2$ and variable $N$. Appendix F
illustrates the additional results when varying $d$ and when using more
realistic nonlinearities.

Figure 7: Properties of the TAP and ITAP solutions for large $M$. Experiments
with RRG with $d=3$, $\gamma=2$. From darker to lighter the curves correspond
to $N=32,64,128,256,512$. The shaded area indicates a $\pm 1$ standard
deviation interval around the mean. In insets, the same quantities as the main
plots are plotted as a function of $\rho$. (a) Efficiency with respect to
shortest paths for TAP (upper curves) and ITAP (lower curves) as a function of
$\eta$. (b) Dependence on $\rho$ of the ratio between the TAP optimal energy
and the energy obtained after projecting the TAP solution onto an integer one
using RITAP. (c) Fraction of non-integer paths in TAP as a function of $\eta$.
The dashed line has a slope of $1/\eta$ (equivalently $1/M$), to illustrate
asymptotic behavior. (d) Average degeneracy of support paths as a function of
$\eta$.

Panel (a) in figure 7 shows the quantity $1-H/H_{SP}$ (averaged over
instances) as a function of $\eta$, for several $N$s. In the upper curves $H$
is the optimal energy in TAP (let us call it $H_{\text{TAP}}$). In the lower
curves it is the energy of the ITAP solution obtained with RITAP; we shall
name this quantity $H_{\text{ITAP}}$. We notice that for low $M$ the two
curves take quite different values, however as $M$ increases the gap between
the two shrinks towards zero. To closely examine the convergence we also plot
the ratio $H_{\text{ITAP}}/H_{\text{TAP}}-1$ in the right panel. The fact that
this quantity goes to zero confirms that $H_{\text{ITAP}}\to H_{\text{TAP}}$
when $\eta\to\infty$. Since $H_{\text{TAP}}$ is a lower bound to
$H_{\text{ITAP}}$, this result shows that RITAP is basically reaching the
global minimum of ITAP, and thus in the high $\eta$ regime, ITAP is easy to
solve approximately.

Next we explore how the convergence in energy reflects on the paths through
which flow is routed in TAP and ITAP solutions. We say a path is integer if it
carries an integer amount of flow at optimality in TAP. Using the notation of
section II.4, consider a TAP solution $\operatorname{\mathcal{P}}$. We define
the fraction of integer paths in TAP as
$F_{\text{int}}\coloneqq\frac{1}{M}\sum_{x,y\in(V\times V)}\sum_{a\in
R_{xy}}\lfloor h_{xy}^{a}\rfloor$. 666$\lfloor h_{xy}^{a}\rfloor$ is indeed
the integer part of the flow being routed on $\pi_{xy}^{a}$. In other words,
$F_{\text{int}}$ is the fraction of TAP optimal paths that already satisfy
ITAP constraints. Figure 7(c) depicts the dependence of $1-F_{\text{int}}$ on
$\eta$. $1-F_{\text{int}}$ decreases as $\eta$ increases, moreover when
$\eta\approx 1$, we observe the onset of a power law behavior $1/\eta$ as
indicated by the dashed black line. Overall this shows that also on the path
level, for $\eta\gg 1$, a solution to TAP is almost a solution to ITAP, in the
sense that almost all TAP paths are integer ($F_{\text{int}}\to 1$). In
Appendix H we prove that $F_{\text{int}}\to 1$ implies that
$H_{\text{ITAP}}/H_{\text{TAP}}\to 1$, therefore rigorously establishing the
relation between path integrality and energy difference in TAP and ITAP.

The main drawback of studying TAP paths at optimality is that these are not
unique, thus we could be observing some peculiarity of the Frank-Wolfe
algorithm rather than a fundamental property of the problem. To erase this
doubt and to cast light onto the behavior of $F_{\text{int}}$ we consider a
quantity which is unique at optimality in TAP. We start from the observation
that, although not unique, path flows at optimality are not arbitrary either.
Based on the optimal edge flows one can determine a set of paths on which
positive flow is possible. We call these the support paths. The following
result by Wardrop [50] characterizes the support paths.

######  Proposition III.1 (Wardrop’s second principle).

Let $G,D$ be a TAP instance and $\\{h_{xy}^{\star a}\\}_{(x,y)\in V\times
V,a\in[|\Pi_{xy}|]}$ be any path flows at optimality. Also let
$\pmb{I^{\star}}$ be the edge flows at optimality. Then for every $(x,y)\in
V\times V$, for every $a\in[|\Pi_{xy}|]$, a necessary condition for
$h_{xy}^{\star a}$ to be strictly positive is that

$\sum_{e\in\pi_{xy}^{a}}\phi^{\prime}(I_{e}^{\star})=\min_{\pi\in\Pi_{xy}}\sum_{e\in\pi}\phi^{\prime}(I_{e}^{\star})$
(23)

A proof of this fact is presented in Appendix H. Wardrop’s second principle
says that the support paths are the shortest paths in the graph whose edges
$e\in E$ are weighed with $\phi^{\prime}(I_{e}^{\star})$. If the optimal flow
between nodes $x$ and $y$ is split across several paths, it means all of these
paths are the shortest paths in the aforementioned graph. We indicate with
$\mathcal{S}_{xy}(\pmb{I^{\star}})\coloneqq\operatorname{\arg\min}_{\pi\in\Pi_{xy}}\sum_{e\in\pi}\phi^{\prime}(I_{e}^{\star})$
the set of support paths starting and terminating respectively at nodes $x$
and $y$. Notice that for a given instance of TAP, the support paths are
unique, as they depend only on the optimal edge flows. The larger the
cardinality of $\mathcal{S}_{xy}(\pmb{I^{\star}})$, the more potential split
across paths, the less integrality in TAP paths. Let $Q\coloneqq\sum_{x,y\in
V\times V}\mathbb{I}[D_{xy}>0]$ be the number of nonzero entries in the demand
matrix. We define the average TAP path degeneracy to be
$K\coloneqq\frac{1}{Q}\sum_{x,y\in V\times
V}\mathbb{I}[D_{xy}>0]\left|\mathcal{S}_{xy}(\pmb{I^{\star}})\right|$. The
behavior of $K-1$ as a function on $\eta$ is plotted in figure 7 (d). For
small $\eta$, there is a huge degeneracy in the support paths, with the flow
being possibly split across hundreds of paths. When $\eta$ increases the
degeneracy decreases and plateaus to a fixed value when $\eta$ becomes larger
than one.

Let us explain why the degeneracy plateauing plays a crucial role in
guaranteeing that TAP converges to ITAP. Fix two nodes $x,y$ and suppose
$\left|\mathcal{S}_{xy}(\pmb{I^{\star}})\right|$ is constant in $M$: when $M$
increases so does $D_{xy}$, however the flow will always be routed through the
same number of paths joining $x$ and $y$. In turn this implies that an
increasing fraction of the flow is integer, resulting in TAP solutions being
close to an ITAP one.777a moment of thought reveals that the fractional part
of the flow is smaller than $\left|\mathcal{S}_{xy}(\pmb{I^{\star}})\right|$,
which is constant in $M$. Therefore the proportion of fractional flow is at
most $\left|\mathcal{S}_{xy}(\pmb{I^{\star}})\right|/D_{xy}$ which goes to
zero. Translating this reasoning in formulas we arrive at the following
identity which gives a lower bound to $F_{\text{int}}$ in terms of the paths
degeneracies.

$F_{\text{int}}\geq 1-\frac{1}{\eta}\left(K-1\right)$ (24)

We give a proof of it in Appendix E.2. For example if
$S_{xy}(\pmb{I^{\star}})=1$ (no degeneracy) for all $(x,y)$, then the bound
gives $F_{\text{int}}\geq 1$. Therefore if $K$ is constant in $M$,
$1-F_{\text{int}}$ goes to zero like $1/\eta$ if not faster, reproducing the
rate observed in figure 7 (c).

One question stands: why when $M$ increases the degeneracy decreases and
ultimately plateaus? We conjecture that this phenomenon is related to the
demand matrix becoming more dense and uniform (in the sense that all entries
are of the same order) for higher $M$. In particular, the distribution of OD
pairs, rather than their number, plays a crucial role in reducing the
degeneracy. We support this claim by showing that the number of OD pairs can
be increased while leaving the degeneracy unchanged. Consider two TAP
instances with the same $G$ and demand matrices respectively $D$ and $2D$. Let
$\phi(x)=x^{\gamma}$ with $\gamma\geq 1$, and suppose we know the path flows
at TAP optimality in the instance $G,D$. Then we show in Appendix H that by
keeping the paths the same and doubling the flow on each of them, we obtain a
TAP optimal solution to the instance $G,2D$. This implies that the support
paths do not change. We conclude that the degeneracy is not influenced by the
norm of $D$. The above results and discussion lead us to formulate the
following conjecture.

######  Conjecture III.2.

Let $G\sim\text{RRG}(N,d)$ be an RRG with $N$ nodes and degree $d$, and $D$ be
a demand matrix obtained by sampling $M$ OD pairs uniformly at random. Suppose
$\phi$ is convex and let $K(G,D)$ be the average degeneracy in the TAP
instance $G,D$. Define
$\kappa(\eta)\coloneqq\lim\limits_{\begin{subarray}{c}M,N\to\infty\\\
M/(N(N-1))=\eta\end{subarray}}\operatorname{\mathbb{E}}_{G,D}K(G,D)$. Then
$\lim\limits_{\eta\to\infty}\kappa(\eta)$ exist and is finite.

Plugging this result into $F_{\text{int}}$ one sees that
$\lim\limits_{\begin{subarray}{c}M,N\to\infty\\\
M/(N(N-1))=\eta\end{subarray}}\operatorname{\mathbb{E}}_{G,D}F_{\text{int}}$,
approaches $1$ as $\eta\to\infty$.

##  IV Two scales present also on real-world graphs

###  IV.1 A tale of two scales

In our above analysis, we identified two main regimes in which the problem’s
behavior is interesting. These correspond respectively to when the effective
flow per edge $\rho$ and $\eta$, defined in (21) and (22), are of order one.
Figure 7 illustrates the presence of the two regimes by showing the same
quantities plotted as a function of $\eta$ and $\rho$, respectively, in the
main panels and in the insets. One observes that when $\rho=\Theta(1)$ the
curves for different sizes collapse when plotted against $\rho$. The
importance of the regime $\rho=\Theta(1)$ is that this is where the effective
interactions among paths arise and are the strongest. Panel (a) shows that the
greatest relative difference between the ITAP energy and the shortest paths
(non-interacting case) one is indeed attained for $\rho\approx 1$. Instead for
$\rho\to 0$ the effective interaction goes to zero and the ITAP energy
approaches that of the shortest paths. 888the fact that the peak of
$1-H/H_{\text{SP}}$ happens for $\rho\approx 1$ is found also in all other
plots, except for figure 3.

Increasing $M$ we move into a setting there $\rho\gg 1$ and $\eta\approx 1$.
In panels (b,c) one observes that the ’knee’ of the curve, where the power law
arises, indeed corresponds to $\eta\approx 1$, but different values of $\rho$,
depending on the size. Similarly, in panel (d), the degeneracy plateaus at
fixed $\eta$ but not at fixed $\rho$.

###  IV.2 Validation on real-world graphs

To prove the robustness of our results we consider two real-world traffic
networks:

  * •

The Anaheim city road network, composed by $N=416$ nodes. The edges are
directed and average total degree (in degree plus out degree) is 4.39.

  * •

The Eastern Massachusetts highway network, with $N=74$ nodes. The edges are
directed and average total degree is 6.97.

Both networks are taken from [20]. We use $\phi(x)=x^{\gamma}$ as nonlinearity
and sample the OD pairs uniformly at random, always with the constraint that
the origin and destination must be distinct. Figure 8 depicts the results of
these experiments.

Figure 8: Experiments on Anaheim and Eastern Massachusetts networks with OD
pairs picked uniformly at random, and nonlinearity $\phi(x)=x^{2}$. The shaded
area indicates a $\pm 1$ standard deviation interval around the mean. The
randomness is over 20 realizations of the OD pairs. In insets, the same
quantities as the main plots are plotted as a function of $\rho$. (a) Relative
energy gain with respect to the shortest path routing. The upper curves
correspond to TAP optimality, the lower ones to the RITAP solution. (b)
Dependence on $\rho$ of the ratio between the TAP optimal energy and the
energy obtained after projecting the TAP solution onto an integer one using
RITAP. (c) Fraction of non-integer paths in TAP as a function of $\eta$. The
dashed line has a slope of $1/\eta$ (equivalently $1/M$), to illustrate
asymptotic behavior. (d) Average degeneracy of support paths as a function of
$\eta$. In the case of Ahaheim the computation of the support paths was too
lengthy computationally, therefore we instead computed the number of paths
over which TAP flow is routed in the FW solution. This is found to agree well
with the actual degeneracy in the cases (low $\eta$) in which it could be
computed.

The quantities that are plotted are the same as in figure 7. As in the RRG
case the ITAP energy approaches the TAP one (panels (a,b)) when $\eta$
increases. Moreover the degeneracy $K-1$ plateaus at high $M$ (panel (d)) thus
producing the $1/\eta$ rate in the quantity $1-F_{\text{int}}$ (panel (c)).
Looking at the insets in panels (a,b,c) we remark that when $\rho=\Theta(1)$
the curves overlap when plotted as a function of $\rho$. Instead, considering
the regime $\eta=\Theta(1)$ we see that curves especially in panels (b,c)
collapse when plotted against $\eta$, with power law behavior arising for
$\eta\approx 1$. This confirms that $\eta$ and $\rho$ and are the relevant
parameters when studying the asymptotic and fixed flow per edge regimes. To
conclude, we also ran additional numerical experiments with typical
nonlinearities borrowed from traffic assignment. We remark that the above
phenomenology is quite robust with respect to the functional form of the
nonlinearity. See Appendix F for the detailed results.

##  V Conclusions

Our work provides an examination of the integer traffic assignment problem on
sparse graphs, with origin-destination pairs picked uniformly at random. We
first devised several algorithms, most of them based on the idea of optimizing
the position of one path while keeping the others fixed. The only exception is
RITAP which is based on relaxing the integral flow constraint. We found a more
transparent derivation of CBP, and proposed a novel method to sample self-
avoiding paths on arbitrary graphs.

To compare the algorithms we focused on the regime where the flow per edge is
of order one, $\rho\coloneqq 2M\log N/(Nd\log d)=O(1)$. In the attractive case
($\phi$ concave), we see that simulated annealing and CBP are the most
effective. CBP, however, fails to converge in some settings.

In the repulsive case ($\phi$ convex) we find that all algorithms except for
RITAP attain approximately the same performance, indicating that elaborated
methods are not significantly more performant than simple ones like the greedy
algorithm. Only simulated annealing exhibits a marginally better performance
when the degree of the RRG increases.

We conclude that in the absence of computational cost constraints, simulated
annealing represents a valid choice to solve ITAP. When taking speed into
consideration, we find that the greedy algorithm and RITAP are orders of
magnitude faster while sacrificing only a few per cent in energy. RITAP in
particular benefits from an implementation that makes its running time
independent of $M$, when $M>N^{2}$. It is thus better suited to cases where
$M$ is very large. Finally, we remark that CBP is also effective, but its
slowness and unreliable convergence limit its applicability.

Next, using RITAP, we investigate the regime $M=\Theta(N^{2})$ where every
entry of the demand matrix is of order one. We do so in the case
$\phi(x)=x^{2}$. We observe that as $\eta=M/(N(N-1))$ grows, TAP approaches
ITAP. Therefore, if we are interested in solving the NP-hard problem ITAP in
this limit, we can instead solve TAP, which is a convex problem and the
relaxed solution will be already very close to an integer one. The convergence
of TAP to ITAP is well explained by the fact that the degeneracy of the
support paths plateaus when $\eta>1$. We conjecture that the degeneracy
becoming constant for higher $M$ originates from the demand matrix becoming
more uniform (in the sense that all entries have similar magnitude). In
particular the distribution of OD pairs, rather than their number plays a
crucial role in controlling the degeneracy.

By fixing a graph and gradually increasing $M$ our study characterizes the
phases encountered by the system. One first finds the regime $\rho=\Theta(1)$,
where paths start interacting. Then, as $\eta$ becomes of order one, the
system enters in an asymptotic regime where most relevant quantities either
plateau (e.g. $1-H/H_{\text{SP}}$, $K$) or exhibit a power law behavior (e.g.
$H_{\text{ITAP}}/H_{\text{TAP}}-1$, $1-F_{\text{int}}$). This phase is
characterized by the aforementioned convergence of TAP to ITAP. The relevance
and properties of the two scales, $\rho=O(1)$ and $\eta=O(1)$, were further
verified in the case of real-world graphs. In future work it would be
important to consider realistic and not random origin-destination pairs and
investigate how that changes the properties of the problem.

## Acknowledgements

We thank Bill Chi Ho Yeung for sharing with us the code of the CBP algorithm.

## References

  * [1] Elliot Anshelevich, Anirban Dasgupta, Jon Kleinberg, Éva Tardos, Tom Wexler, and Tim Roughgarden.  The price of stability for network design with fair cost allocation.  SIAM Journal on Computing, 38(4):1602–1623, 2008.
  * [2] Hillel Bar-Gera.  Origin-based algorithm for the traffic assignment problem.  Transportation Science, 36(4):398–417, 2002.
  * [3] Cynthia Barnhart, Christopher A Hane, and Pamela H Vance.  Using branch-and-price-and-cut to solve origin-destination integer multicommodity flow problems.  Operations Research, 48(2):318–326, 2000.
  * [4] Martin Beckmann, Charles B McGuire, and Christopher B Winsten.  Studies in the economics of transportation.  Technical report, 1956.
  * [5] C Bishop.  Pattern recognition and machine learning.  Springer google schola, 2:531–537, 2006.
  * [6] Lawrence E Blume.  The statistical mechanics of strategic interaction.  Games and economic behavior, 5(3):387–424, 1993.
  * [7] Stephen D Boyles, Nicholas E Lownes, and Avinash Unnikrishnan.  Transportation network analysis.  Volume I, Version 0.85, 2020.
  * [8] George Casella and Edward I George.  Explaining the gibbs sampler.  The American Statistician, 46(3):167–174, 1992.
  * [9] Tommaso Castellani and Andrea Cavagna.  Spin-glass theory for pedestrians.  Journal of Statistical Mechanics: Theory and Experiment, 2005(05):P05012, 2005.
  * [10] Imrich Chlamtac, Aura Ganz, and Gadi Karmi.  Lightpath communications: An approach to high bandwidth optical wan’s.  IEEE transactions on communications, 40(7):1171–1182, 1992.
  * [11] Thomas H Cormen, Charles E Leiserson, Ronald L Rivest, and Clifford Stein.  Introduction to algorithms.  MIT press, 2022.
  * [12] Stella C Dafermos and Frederick T Sparrow.  The traffic assignment problem for a general network.  Journal of Research of the National Bureau of Standards B, 73(2):91–118, 1969.
  * [13] M Daoud, JP Cotton, B Farnoux, Gerard Jannink, G Sarma, H Benoit, C Duplessix, Cl Picot, and PG De Gennes.  Solutions of flexible polymers. neutron experiments and interpretation.  Macromolecules, 8(6):804–818, 1975.
  * [14] Robert B Dial.  A path-based user-equilibrium traffic assignment algorithm that obviates path storage and enumeration.  Transportation Research Part B: Methodological, 40(10):917–936, 2006.
  * [15] Edsger W Dijkstra.  A note on two problems in connexion with graphs.  In Edsger Wybe Dijkstra: His Life, Work, and Legacy, pages 287–290. 2022.
  * [16] Sven Erlander.  On the relationship between the discrete and continuous models for combined distribution and assignment.  Transportation Research Part B: Methodological, 22(5):371–382, 1988.
  * [17] Shimon Even, Alon Itai, and Adi Shamir.  On the complexity of time table and multi-commodity flow problems.  In 16th annual symposium on foundations of computer science (sfcs 1975), pages 184–193. IEEE, 1975.
  * [18] Alex Fabrikant, Christos Papadimitriou, and Kunal Talwar.  The complexity of pure nash equilibria.  In Proceedings of the thirty-sixth annual ACM symposium on Theory of computing, pages 604–612, 2004.
  * [19] Angelo Fanelli, Michele Flammini, and Luca Moscardelli.  The speed of convergence in congestion games under best-response dynamics.  In International Colloquium on Automata, Languages, and Programming, pages 796–807. Springer, 2008.
  * [20] Transportation Networks for Research Core Team.  Transportation networks for research.  https://github.com/bstabler/TransportationNetworks.
  * [21] Marguerite Frank, Philip Wolfe, et al.  An algorithm for quadratic programming.  Naval research logistics quarterly, 3(1-2):95–110, 1956.
  * [22] Stuart Geman and Donald Geman.  Stochastic relaxation, gibbs distributions, and the bayesian restoration of images.  IEEE Transactions on pattern analysis and machine intelligence, PAMI-6(6):721–741, 1984.
  * [23] Kerry Hinton, Jayant Baliga, Michael Feng, Robert Ayre, and Rodney S Tucker.  Power consumption and energy efficiency in the internet.  IEEE Network, 25(2):6–12, 2011.
  * [24] Scott Kirkpatrick, C Daniel Gelatt Jr, and Mario P Vecchi.  Optimization by simulated annealing.  science, 220(4598):671–680, 1983.
  * [25] Pieter Kleer.  Sampling from the gibbs distribution in congestion games.  In Proceedings of the 22nd ACM Conference on Economics and Computation, pages 679–680, 2021.
  * [26] Pieter Kleer and Guido Schäfer.  Computation and efficiency of potential function minimizers of combinatorial congestion games.  Mathematical Programming, 190:523–560, 2021.
  * [27] Donald E Knuth.  Mathematics and computer science: Coping with finiteness: Advances in our ability to compute are bringing us substantially closer to ultimate limitations.  Science, 194(4271):1235–1242, 1976.
  * [28] Elias Koutsoupias and Christos Papadimitriou.  Worst-case equilibria.  In Annual symposium on theoretical aspects of computer science, pages 404–413. Springer, 1999.
  * [29] Neal Madras and Alan D Sokal.  The pivot algorithm: A highly efficient monte carlo method for the self-avoiding walk.  Journal of Statistical Physics, 50:109–186, 1988.
  * [30] A Meli.  Simulation and analysis of various routing algorithms for optical networks.  2004\.
  * [31] Marc Mezard and Andrea Montanari.  Information, physics, and computation.  Oxford University Press, 2009.
  * [32] Marc Mézard, Giorgio Parisi, and Miguel Angel Virasoro.  Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications, volume 9.  World Scientific Publishing Company, 1987.
  * [33] Maria Mitradjieva and Per Olov Lindberg.  The stiff is moving—conjugate direction frank-wolfe methods with applications to traffic assignment.  Transportation Science, 47(2):280–293, 2013.
  * [34] Asuman E Ozdaglar and Dimitri P Bertsekas.  Optimal solution of integer multicommodity flow problems with application in optical networks.  In Frontiers in global optimization, pages 411–435. Springer, 2004.
  * [35] Michael Patriksson.  The traffic assignment problem: models and methods.  Courier Dover Publications, 2015.
  * [36] Judea Pearl.  Reverend bayes on inference engines: a distributed hierarchical approach.  In Proceedings of the Second AAAI Conference on Artificial Intelligence, AAAI’82, page 133–136. AAAI Press, 1982.
  * [37] Rudolf Peierls.  Statistical theory of superlattices with unequal concentrations of the components.  Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences, 154(881):207–222, 1936.
  * [38] Olga Perederieieva, Matthias Ehrgott, Andrea Raith, and Judith YT Wang.  A framework for and empirical study of algorithms for traffic assignment.  Computers & Operations Research, 54:90–107, 2015.
  * [39] Arthur Pigou.  The economics of welfare.  Macmillan, 1920.
  * [40] Christian P Robert, George Casella, and George Casella.  Monte Carlo statistical methods, volume 2.  Springer, 1999.
  * [41] Robert W Rosenthal.  A class of games possessing pure-strategy nash equilibria.  International Journal of Game Theory, 2:65–67, 1973.
  * [42] Tim Roughgarden.  Selfish routing and the price of anarchy.  MIT press, 2005.
  * [43] Tim Roughgarden.  Algorithmic game theory.  Communications of the ACM, 53(7):78–86, 2010.
  * [44] Khodakaram Salimifard and Sara Bigharaz.  The multicommodity network flow problem: state of the art classification, applications, and solution methods.  Operational Research, pages 1–47, 2022.
  * [45] Farhad Shahrokhi and David W Matula.  The maximum concurrent flow problem.  Journal of the ACM (JACM), 37(2):318–334, 1990.
  * [46] Mandana Tabrizi.  Rosenbluth algorithm studies of self-avoiding walks.  2015\.
  * [47] David J Thouless, Philip W Anderson, and Robert G Palmer.  Solution of’solvable model of a spin glass’.  Philosophical Magazine, 35(3):593–601, 1977.
  * [48] Ido Tishby, Ofer Biham, Reimer Kühn, and Eytan Katzav.  The mean and variance of the distribution of shortest path lengths of random regular graphs.  Journal of Physics A: Mathematical and Theoretical, 55(26):265005, 2022.
  * [49] EJ Van Rensburg.  The statistical mechanics of interacting walks, polygons, animals and vesicles (2nd edn).  The Statistical Mechanics of Interacting Walks, 2015.
  * [50] J G WARDROP.  Road paper. some theoretical aspects of road traffic research.  Proceedings of the Institution of Civil Engineers, 1(3):325–362, 1952.
  * [51] Wei Wei, Longxian Hu, Qiuwei Wu, and Tao Ding.  Efficient computation of user optimal traffic assignment via second-order cone and linear programming techniques.  IEEE Access, 7:137010–137019, 2019.
  * [52] Andrés Weintraub, Carmen Ortiz, and Jaime González.  Accelerating convergence of the frank-wolfe algorithm.  Transportation Research Part B: Methodological, 19(2):113–122, 1985.
  * [53] Jin Y Yen.  Finding the k shortest loopless paths in a network.  management Science, 17(11):712–716, 1971.
  * [54] Chi Ho Yeung, David Saad, and KY Michael Wong.  From the physics of interacting polymers to optimizing routes on the london underground.  Proceedings of the National Academy of Sciences, 110(34):13717–13722, 2013.
  * [55] Lenka Zdeborová and Florent Krzakala.  Statistical physics of inference: Thresholds and algorithms.  Advances in Physics, 65(5):453–552, 2016.

##  Appendix A Pseudocode of the greedy algorithm

The following is a pseudocode of an efficient implementation of the greedy
algorithm.

Algorithm 2 Greedy algorithm

1:Input:  $G=(V,E)$ graph over which paths are sampled,
$\\{s^{\mu}\\}_{\mu\in[M]},\\{t^{\mu}\\}_{\mu\in[M]}$ respectively the origins
and destinations of each path, $\\{\pi^{\mu}_{0}\\}_{\mu\in[M]}$ set of
initial paths.

2:Output: A set of paths $\\{\pi^{\mu}_{0}\\}_{\mu\in[M]}$, that correspond to
a local minimum of $H$.

3:Gflow$\leftarrow$ copy($G$) $\triangleright$ the weight on edge $e$ of Gflow
is $I_{e}$

4:Gcost $\leftarrow$ copy($G$) $\triangleright$ the weight on edge $e$ of
Gcost is $\phi(I_{e})$

5:for $e\in G.$edges do

6: Gflow.weight($e$)$\leftarrow$0

7:end for

8:for $\mu=1,\dots,M$ do

9: $\pi^{\mu}\leftarrow\pi_{0}^{\mu}$

10: for $k=1,\dots,$ length($\pi^{\mu})-1$ do

11: $e\leftarrow[\pi^{\nu}[k],\pi^{\nu}[k+1]]$

12: Gflow.weight($e$)$\leftarrow$ Gflow.weight($e$)+1 $\triangleright$
populating the flow graph with paths

13: end for

14:end for

15:H$\leftarrow 0$

16:for $e\in G.$edges do $\triangleright$ loop to compute current value of
$H(\pmb{I})$

17:
Gcost.weight($e$)$\leftarrow\phi(\text{Gflow.weight($e$)}+1)-\phi(\text{Gflow.weight($e$)})$

18: H$\leftarrow$ H+Gcost.weight($e$)

19:end for

20:flagconv $\leftarrow$ False$\triangleright$ flag for convergence

21:while not flagconv do$\triangleright$ main loop

22: for $\nu=1,\dots,M$ do$\triangleright$ loop over all paths

23: for $k=1,\dots,\text{length($\pi^{\nu}$)}-1$ do

24: $e\leftarrow[\pi^{\nu}[k],\pi^{\nu}[k+1]]$

25: Gflow.weight($e$)$\leftarrow$ Gflow.weight($e)-1$$\triangleright$ removing
path $\nu$ from the graph so that one has $I^{\backslash\nu}$

26:
Gcost.weight($e$)$\leftarrow\phi(\text{Gflow.weight($e$)}+1)-\phi(\text{Gflow.weight($e$)})$

27: end for

28: $\pi^{\nu}\leftarrow$ shortestPath(Gcost,
$s^{\nu},t^{\nu}$)$\triangleright$ updating path $\nu$

29: for $k=1,\dots,\text{length($\pi^{\nu}$)}-1$ do

30: $e\leftarrow[\pi^{\nu}[k],\pi^{\nu}[k+1]]$

31: Gflow.weight($e$)$\leftarrow$ Gflow.weight($e$)+1$\triangleright$ adding
back path $\nu$ to the graph

32:
Gcost.weight($e$)$\leftarrow\phi(\text{Gflow.weight($e$)}+1)-\phi(\text{Gflow.weight($e$)})$

33: end for

34: end for

35: newH$\leftarrow 0$ $\triangleright$ measuring energy after all paths are
updated

36: for $e\in G.$edges do

37: newH$\leftarrow$newH+Gcost.weight($e$)

38: end for

39: if H == newH then

40: flagconv$\leftarrow$True $\triangleright$ If $H(\pmb{I})$ does not change
once all the paths are updated, then the algorithm has converged

41: end if

42: H $\leftarrow$ newH

43:end while

The method ’.weight($e$)’ accesses the weight of edge $e$. ’$G$.edges’ is a
list of all the edges in $G$. Each path is represented as a list of nodes,
i.e., $\pi=[\pi[1],\pi[2],\dots,\pi_{\text{length}(\pi)}]$. Similarly an edge
is represented as a list of two nodes. The function ’shortestPath($G,s,t$)’
instead returns she shortest path from $s$ to $t$ in the weighted graph $G$.

##  Appendix B Derivation of CBP

The conditional probability distribution (19) can be represented by a factor
graph, where variable nodes are the at the edges and factor nodes at the nodes
of the original network. Let $(ij)$ be an (undirected) edge on the graph $G$.
We place two variables on this edge $\sigma_{i\to j}^{\nu}=\mathbb{I}[i\to
j\in\pi^{\nu}]$ and $\sigma_{j\to i}^{\nu}=\mathbb{I}[j\to i\in\pi^{\nu}]$.
The probability distribution (II.2)can be decomposed into factors as
$P\left(\\{\sigma^{\nu}_{i\to j}\\}_{(ij)\in E}|\\{H^{\nu}_{ij}\\}_{(ij)\in
E},s^{\nu},t^{\nu}\right)\propto
f_{s^{\nu}}f_{t^{\nu}}\prod_{i\in[N]\backslash s^{\nu},t^{\nu}}f_{i}$. Let us
look at these factors:

  * •

$f_{i}=\mathbb{I}[\sum_{\begin{subarray}{c}k,l\in\partial i\\\ k\neq
l\end{subarray}}\sigma^{\nu}_{k\to i}+\sigma^{\nu}_{i\to l}=(0\text{ OR
}2)]\exp{(-\beta\sum_{k\in\partial i}H^{\nu}_{ik}\sigma^{\nu}_{i\rightarrow
k})}]$. This constraint ensures that $\pi^{\nu}$ passes through node $i$ no
more than once. It also ensures that the $\\{\sigma_{i\to j}\\}$ form a valid
path and it contains the weighting factor $e^{-\beta H^{\nu}_{ik}}$ that
represents the cost of passing through edge $i,k$, independently of direction.
This constraint is placed at every node in the graph except for
$s^{\nu},t^{\nu}$

  * •

$f_{t^{\nu}}=\mathbb{I}\left[\sum_{k\in\partial t^{\nu}}\sigma^{\nu}_{k\to
t^{\nu}}=1\right]\mathbb{I}\left[\sum_{k\in\partial
t^{\nu}}\sigma^{\nu}_{t^{\nu}\to k}=0\right]$ imposes that the path ends at
$t^{\nu}$.

  * •

$f_{s^{\nu}}=\mathbb{I}\left[\sum_{k\in\partial
s^{\nu}}\sigma^{\nu}_{s^{\nu}\to k}=1\right]\mathbb{I}\left[\sum_{k\in\partial
s^{\nu}}\sigma^{\nu}_{k\to s^{\nu}}=0\right]\prod_{k\in\partial
s^{\nu}}(\exp{(-\beta H^{\nu}_{s^{\nu}k}\sigma^{\nu}_{s^{\nu}\rightarrow
k})})$ imposes that the path starts at $s^{\nu}$.

Having defined the quantity $f_{i}$, one can draw the factor graph as shown in
Figure 9.

Figure 9: Factor graph associated to the probability conditional probability
distribution (II.2).

From the construction of the graph, one can deduce that if the original
network is tree like, the factor graph will remain tree like given our choice
of variables. Given this property, BP equations are expected to converge to a
good solution if the graph does not have too many loops.

We can now derive the BP equations for this factor graph; this gives two sets
of messages. Since each variable node has exactly two neighbouring factor
nodes, simplifications of the BP equations leads to the iteration of only one
set of messages. We denote by $\chi^{\nu}_{i\rightarrow j}$ the message going
from node $i$ to node $j$. The message $\chi^{\nu}_{i\rightarrow
j}(\sigma^{\nu}_{i\to j},\sigma^{\nu}_{j\to i})$ depends on the value of the
variable, that can take four different values. First notice that $\nu_{i\to
j}(1,1)=0$, since the path cannot backtrack on one edge. For the other three
values ($\sigma^{\nu}_{i\to j},\sigma^{\nu}_{j\to i})\in$ {(0,0),(1,0),(0,1)}
(for a node $i\neq s^{\nu},e^{\nu}$), update equations are derived:

$\displaystyle\chi^{\nu}_{i\to j}(0,0)$ $\displaystyle=\frac{1}{Z_{i\to
j}}\left[\prod_{k\in\partial i\backslash j}\chi^{\nu}_{k\to
i}(0,0)+\sum_{\begin{subarray}{c}k,l\in\partial i\backslash j\\\ k\neq
l\end{subarray}}e^{-\beta H^{\nu}_{il}}\chi^{\nu}_{k\to
i}(1,0)\chi^{\nu}_{l\to i}(0,1)\prod_{r\in\partial i\backslash
j,k,l}\chi^{\nu}_{r\to i}(0,0)\right]$ $\displaystyle\chi^{\nu}_{i\to j}(1,0)$
$\displaystyle=\frac{e^{-\beta H^{\nu}_{ij}}}{Z_{i\to j}}\sum_{k\in\partial
i\backslash j}\chi^{\nu}_{k\to i}(1,0)\prod_{l\in\partial i\backslash
j,k}\chi^{\nu}_{l\to i}(0,0)$ (25) $\displaystyle\chi^{\nu}_{i\to j}(0,1)$
$\displaystyle=\frac{1}{Z_{i\to j}}\sum_{k\in\partial i\backslash j}e^{-\beta
H^{\nu}_{ik}}\chi^{\nu}_{k\to i}(0,1)\prod_{l\in\partial i\backslash
j,k}\chi^{\nu}_{l\to i}(0,0)$

Equations B are intuitive. For example, to derive the expression of
$\chi^{\nu}_{i\to j}(0,0)$, two cases need to be considered: the case where
none of the neighboring edges to $i$ is occupied, and the case where a path
passes through node $i$, but without passing through node $j$ (and therefore
passes through $k\to i\to l$ with $k,l\neq j$). Same reasoning can be applied
for the two other messages.

Special equations apply for the source node $s^{\nu}$:

$\displaystyle\chi^{\nu}_{s^{\nu}\to j}(0,0)$
$\displaystyle=\frac{1}{Z_{s^{\nu}\to j}}\sum_{k\in\partial s^{\nu}\backslash
j}e^{-\beta H^{\nu}_{s^{\nu}k}}\chi^{\nu}_{k\to
s^{\nu}}(0,1)\prod_{l\in\partial s^{\nu}\backslash j,k}\chi^{\nu}_{l\to
s^{\nu}}(0,0)$ $\displaystyle\chi^{\nu}_{s^{\nu}\to j}(1,0)$
$\displaystyle=e^{-\beta H^{\nu}_{s^{\nu}j}}\frac{1}{Z_{s^{\nu}\to
j}}\prod_{k\in\partial s^{\nu}\backslash j}\chi^{\nu}_{k\to s^{\nu}}(0,0)$
(26) $\displaystyle\chi^{\nu}_{s^{\nu}\to j}(0,1)$ $\displaystyle=0$

and for the end node $e^{\nu}$ :

$\displaystyle\chi^{\nu}_{e^{\nu}\to j}(0,0)$
$\displaystyle=\frac{1}{Z_{e^{\nu}\to j}}\sum_{k\in\partial e\backslash
j}\chi^{\nu}_{k\to e^{\nu}}(1,0)\prod_{l\in\partial e\backslash
j,k}\chi^{\nu}_{l\to e^{\nu}}(0,0)$ $\displaystyle\chi^{\nu}_{e^{\nu}\to
j}(0,1)$ $\displaystyle=\frac{1}{Z_{e^{\nu}\to j}}\prod_{k\in\partial
e\backslash j}\chi^{\nu}_{k\to e^{\nu}}(0,0)$ (27)
$\displaystyle\chi^{\nu}_{e^{\nu}\to j}(1,0)$ $\displaystyle=0$

To get rid of the normalization factor $1/Z_{i\rightarrow j}$ present in the
messages, the following change of variables can be applied:

$\displaystyle z^{\nu}_{i\rightarrow j}=\frac{\chi^{\nu}_{i\rightarrow
j}(1,0)}{\chi^{\nu}_{i\rightarrow j}(0,0)}\text{ and
}\tilde{z}^{\nu}_{i\rightarrow j}=\frac{\chi^{\nu}_{i\rightarrow
j}(0,1)}{\chi^{\nu}_{i\rightarrow j}(0,0)}e^{-\beta H^{\nu}_{ji}}$ (28)

Using equations B, one finds a closed form for the normalized messages
$z^{\nu}_{i\rightarrow j}$ and $\tilde{z}^{\nu}_{i\rightarrow j}$ :

$\displaystyle z^{\nu}_{i\to j}$
$\displaystyle=\operatorname{\mathbb{I}}(i\neq
s^{\nu})\operatorname{\mathbb{I}}(i\neq e^{\nu})\frac{e^{-\beta
H^{\nu}_{ij}}\sum_{k\in\partial i\backslash j}z^{\nu}_{k\to
i}}{1+\sum_{\begin{subarray}{c}k,l\in\partial i\backslash j\\\ k\neq
l\end{subarray}}z^{\nu}_{k\to i}\tilde{z}^{\nu}_{l\to
i}}+\operatorname{\mathbb{I}}(i=s^{\nu})\frac{e^{-\beta
H^{\nu}_{s^{\nu}j}}}{\sum_{k\in\partial s^{\nu}\backslash
j}\tilde{z}^{\nu}_{k\to s^{\nu}}}$ (29) $\displaystyle\tilde{z}^{\nu}_{i\to
j}$ $\displaystyle=\operatorname{\mathbb{I}}(i\neq
s^{\nu})\operatorname{\mathbb{I}}(i\neq e^{\nu})\frac{\sum_{k\in\partial
i\backslash j}\tilde{z}^{\nu}_{k\to
i}}{1+\sum_{\begin{subarray}{c}k,l\in\partial i\backslash j\\\ k\neq
l\end{subarray}}z^{\nu}_{k\to i}\tilde{z}^{\nu}_{l\to
i}}+\operatorname{\mathbb{I}}(i=e^{\nu})\frac{1}{\sum_{k\in\partial
e^{\nu}\backslash j}z^{\nu}_{k\to e^{\nu}}}$ (30)

Equations 29 and 30 are fixed point equations for the messages $z^{\nu}_{i\to
j}$ and $\tilde{z}^{\nu}_{i\to j}$ at finite $\beta$. In order to take the
$\beta\rightarrow\infty$ limit, it is useful to define the quantities
$a^{\nu}_{i\to j}=-\log(z^{\nu}_{i\to j})/\beta$ and $b^{\nu}_{i\to
j}=-\log(\tilde{z}^{\nu}_{i\to j})/\beta$. After taking the zero temperature
limit ($\beta\to\infty$), the belief propagation equations thus provide update
equations for the messages $a^{\nu}_{i\rightarrow j}$ 13 and
$b^{\nu}_{i\rightarrow j}$ 18.

The zero temperature belief propagation equations presented allow, if iterated
until convergence, to find the path $\nu$ minimizing the energy once all other
paths $\mu\in[M]\backslash\nu$ are fixed. Indeed, it samples from the
conditional probability distribution II.2 that concentrates around the
configurations of lowest energy when $\beta$ goes to $\infty$. As already
mentioned in section II.2, following this procedure results in one step of the
greedy procedure. In other words, BP is used as an algorithm to find the
shortest path. However, the procedure described in [54] consists in only
making one iteration of the update equations 13 and 18 before updating the
cost on edges. The cost $H_{ij}^{\nu}$ at the edge $(ij)$ being defined as
$H_{ij}^{\nu}=\phi(I^{\backslash\nu}_{ij}+1)-\phi(I^{\backslash\nu}_{ij})$, it
is necessary to compute the expected total traffic of the other paths
$I^{\backslash\nu}_{ij}$ between $i$ and $j$. We can write the expected
traffic as
$I^{\backslash\nu}_{ij}=\sum_{\mu\in[M]\backslash\nu}P(\sigma^{\mu}_{(ij)}=1)$
where $\sigma^{\mu}_{(ij)}=\operatorname{\mathbb{I}}[i\to j\in\pi^{\nu}\text{
OR }j\to i\in\pi^{\nu}]$. The expected contribution of path $\mu$ to the
traffic at edge $(ij)$ is equal to the probability $P(\sigma^{\mu}_{(ij)}=1)$.
Naturally, we have that $P(\sigma^{\mu}_{(ij)}=1)=P(\sigma^{\mu}_{i\rightarrow
j}=1)+P(\sigma^{\mu}_{j\rightarrow i}=1)$, where $P(\sigma^{\mu}_{i\rightarrow
j}=1)$ can be expressed using the messages $a^{\mu}_{i\rightarrow j}$ and
$b^{\mu}_{i\rightarrow j}$. In order to derive the expression, one can start
by deducing the marginal from the finite temperature messages
$\chi^{\mu}_{i\rightarrow j}$, then operate the same change of variables 28 as
follows :

$\displaystyle P(\sigma^{\mu}_{i\rightarrow j}=1)$
$\displaystyle=\frac{{\chi^{\mu}_{i\rightarrow j}(1,0)}\sum_{k\in\partial
j\backslash i}e^{-\beta H^{\mu}_{jk}}\chi^{\mu}_{k\rightarrow
j}(0,1)\prod_{l\in\partial j\backslash i,k}\chi^{\mu}_{l\rightarrow
j}(0,0)}{\prod_{k\in\partial j}\chi^{\mu}_{k\rightarrow
j}(0,0)+\sum_{k,l\in\partial j,k\neq l}e^{-\beta
H^{\mu}_{jl}}\chi^{\mu}_{k\rightarrow j}(1,0)\chi^{\mu}_{l\rightarrow
j}(0,1)\prod_{r\in\partial j\backslash l,k}\chi^{\mu}_{r\rightarrow j}(0,0)}$
$\displaystyle=\frac{z^{\mu}_{i\rightarrow j}\sum_{k\in\partial j\backslash
i}\tilde{z}^{\mu}_{k\rightarrow j}}{1+\sum_{k,l\in\partial j,k\neq
l}z^{\mu}_{k\rightarrow j}\tilde{z}^{\mu}_{l\rightarrow j}}$ (31)

Taking the limit $\beta\rightarrow\infty$, the probability
$P(\sigma^{\mu}_{i\rightarrow j}=1)$ concentrates at one value (either 0 or 1)
depending on the values of the messages. Explicitly, the value of
$\sigma^{\mu}_{i\rightarrow j}$ (including the case $i=s^{\mu}$ and
$i=t^{\mu}$) will follow the expression:

$\displaystyle\sigma^{\mu}_{i\rightarrow
j}=\operatorname{\mathbb{I}}(i=s^{\mu})$
$\displaystyle\Theta\bigg{(}\min_{l\in\partial j\backslash
i}(b^{\mu}_{l\rightarrow j})-b^{\mu}_{i\rightarrow j}\bigg{)}$
$\displaystyle+\operatorname{\mathbb{I}}(i\neq
s^{\mu})\Theta\bigg{(}\min{\bigg{(}0,\min_{\begin{subarray}{c}l\in\partial
j\backslash i\end{subarray}}{(a^{\mu}_{l\rightarrow
j}+\min_{\begin{subarray}{c}k\in\partial j\backslash
l\end{subarray}}b^{\mu}_{k\rightarrow j})}\bigg{)}}-(a^{\mu}_{i\rightarrow
j}+\min_{\begin{subarray}{c}k\in\partial j\backslash
i\end{subarray}}b^{\mu}_{k\rightarrow j})\bigg{)}$ (32)

The results derived in this appendix can be compared to the ones exposed in
the supplementary material of [54] and to the pseudocode of the CBP algorithm
(also taken from [54]). The update equations for $a^{\nu}_{i\rightarrow j}$
(13) and $b^{\nu}_{i\rightarrow j}$ (18) can be compared to equations
(35),(36) which are the ones originally derived in [54]. The two sets of
equations, are almost identical, we shall now focus on the differences. The
main difference is the presence of $H_{ij}^{\nu}$ in place of
$\phi^{\prime}(\eta_{ji}^{\nu*})$. This is due to the fact that the authors
use an approximated form of the energy obtained by Taylor expanding
$H_{ij}^{\nu}$. This approximation is only valid if $I_{i\rightarrow
j}^{\backslash\nu}\gg 1$ which is not the case in the regime $\rho=O(1)$,
which is the one considered in [54]. Another difference lies in the update
equation of message $b^{\nu}_{i\rightarrow j}$ in the case $i=e^{\nu}$.
Indeed, the term $-\min_{k\in\partial e^{\nu}\backslash
j}[a^{\nu}_{k\rightarrow e^{\nu}}]$ is present in equation (18) and is missing
from equation (36). Additionally, some algebraic computations reveal that
equation (19) exactly corresponds to equation (B.1), originally derived in
[54].

Despite the presence of these differences, in our experiments we use the code
from [54]. We further tried using the exact formula for $H_{ij}^{\nu}$ in
place of the Taylor expansion and this did not improve the energy reached by
CBP.

###  B.1 CBP pseudocode

The pseudo-code corresponding to the implementation made by the authors of
[54] is given in 3. We report below the formulas to update the binary
variables of each path based on the values of the messages. We say
$\sigma_{ji}^{\nu}$ is a binary variable indicating if path $\nu$ goes through
the undirected edge $(ji)$ and $\eta_{ji}^{\nu*}$ the expected traffic through
edge $(ji)$ taking into account that path $\nu$ passes through it (thus the
+1). We can determine these quantities from the messages as

$\displaystyle\sigma_{ji}^{\nu}$
$\displaystyle=\operatorname{\mathbb{I}}(i=s^{\nu})\Theta\left(\min_{l\in\partial
j\backslash i}\left[b_{l\rightarrow j}^{\nu}\right]-b_{i\rightarrow
j}^{\nu}\right)+\operatorname{\mathbb{I}}(i=t^{\nu})\Theta\left(\min_{l\in\partial
j\backslash i}\left[a_{l\rightarrow j}^{\nu}\right]-a_{i\rightarrow
j}^{\nu}\right)$ $\displaystyle+\operatorname{\mathbb{I}}(i\neq
s^{\nu})\operatorname{\mathbb{I}}(i\neq
t^{\nu})\Theta\left(\min\left[0,\min_{\begin{subarray}{c}l,r\in\partial
j\backslash i\\\ l\neq r\end{subarray}}\left[a_{l\rightarrow
j}^{\nu}+b_{r\rightarrow j}^{\nu}\right]\right]\right.$
$\displaystyle\left.-\min\left[a_{i\rightarrow j}^{\nu}+\min_{l\in\partial
j\backslash i}\left[b_{l\rightarrow j}^{\nu}\right],b_{i\rightarrow
j}^{\nu}+\min_{l\in\partial j\backslash i}\left[a_{l\rightarrow
j}^{\nu}\right]\right]\right),$ (33)
$\eta_{ji}^{\nu*}=1+\sum_{\mu\neq\nu}\sigma_{ji}^{\mu},$ (34)

The last two equations are respectively referenced to in the pseudo-code by
formulas $F_{1}(\\{a_{r\rightarrow j}^{\nu}\\}_{r\in\partial
j},\\{b_{s\rightarrow j}^{\nu}\\}_{s\in\partial j})$ and
$F_{2}(\\{\sigma^{\mu}_{(ji)}\\}_{\mu\in[M]\backslash\nu})$.

We now give the update equations for the sets of messages $a$ and $b$.

$\displaystyle a_{j\rightarrow i}^{\nu}=\begin{cases}\min_{l\in\partial
j\backslash i}\left[a_{l\rightarrow
j}^{\nu}\right]+\phi^{\prime}\left(\eta_{ji}^{\nu*}\right)-\min\left[0,\min_{\begin{subarray}{c}l,r\in\partial
j\backslash i\\\ l\neq r\end{subarray}}\left[a_{l\rightarrow
j}^{\nu}+b_{r\rightarrow j}^{\nu}\right]\right],&\text{if
}j\not\in\\{s^{\nu},e^{\nu}\\}\\\ -\min_{l\in\partial j\backslash
i}\left[b_{l\rightarrow
j}^{\nu}\right]+\phi^{\prime}\left(\eta_{ji}^{\nu*}\right),&\text{if
}j=s^{\nu}\\\ \infty,&\text{if }j=e^{\nu}\end{cases}$ (35) $\displaystyle
b_{j\rightarrow i}^{\nu}=\begin{cases}\min_{l\in\partial j\backslash
i}\left[b_{l\rightarrow
j}^{\nu}\right]+\phi^{\prime}\left(\eta_{ji}^{\nu*}\right)-\min\left[0,\min_{\begin{subarray}{c}l,r\in\partial
j\backslash i\\\ l\neq r\end{subarray}}\left[a_{l\rightarrow
j}^{\nu}+b_{r\rightarrow j}^{\nu}\right]\right],&\text{if
}j\not\in\\{s^{\nu},e^{\nu}\\}\\\ \infty,&\text{if }j=s^{\nu}\\\
\phi^{\prime}\left(\eta_{ji}^{\nu*}\right),&\text{if }j=e^{\nu}\end{cases}$
(36)

where the term $\phi^{\prime}\left(\eta_{ji}^{\nu*}\right)$ quantifies the
contribution for the total energy of path $\nu$ passing through the undirected
edge $(ji)$, using the Taylor expansion as explained previously. Formula
$F_{3}(\\{a_{r\rightarrow i}^{\nu}\\}_{r\in\partial i\backslash
j},\\{b_{s\rightarrow i}^{\nu}\\}_{s\in\partial i\backslash
j},\phi^{\prime}(\eta_{ji}^{\nu*}))$ allows to update the messages and is
reported in the equations (35) and (36). The initialization of the messages is
chosen in order to bias the algorithm to converge towards a configurations
with consistent paths and this is done in the following way. Let us consider
the path $\nu$. While the messages related to path $\nu$ involving all nodes
except $s^{\nu}$ and $e^{\nu}$ are initialized to a positive random value
between 0 and 1. Instead the messages $\\{a_{s^{\nu}\rightarrow
j}^{\nu}\\}_{j\in\partial s^{\nu}}$ exiting source node are set to a random
negative value between -1 and 0 (the probability for the edge to be occupied
is higher). Also, given that path $\nu$ is not allowed to exit the end node
$t^{\nu}$, it is natural to initialize the messages $\\{a_{t^{\nu}\rightarrow
j}^{\nu}\\}_{j\in\partial t^{\nu}}$ to infinity (corresponding to a
probability 0). In the same way, all messages in $\\{a_{s^{\nu}\rightarrow
j}^{\nu}\\}_{j\in\partial s^{\nu}}$ are initialized to infinity as no value as
edges entering node $s^{\nu}$ should not be occupied by the path $\nu$. Some
other implementation choices have been made in order to improve the
convergence of the algorithm. For example, an additive noise has been
introduced in the computation on the costs on the edges (term
$\phi^{\prime}(\eta_{ij}^{\nu*})$ line 13 of pseudo-code 3) in order to break
the degeneracy (cases where two path configurations have the same energy).

Algorithm 3 CBP algorithm for routing

1:Input:  $G=(V,E)$ graph over which paths are computed. Sources
$\\{s^{\mu}\\}_{\mu\in[M]}$ and destinations $\\{s^{\mu}\\}_{\mu\in[M]}$ of
each path. The maximum number of iterations $MaxIter$. Equilibrium time
$EqTime$ (stopping condition). A tolerance level $\epsilon$ to determine if
the messages have converged.

2:Output: $\\{\sigma_{(ij)}^{\nu}\\}$ are binary variables representing the
occupancy of each edge by each path (=1 if undirected edge $(ij)$ is occupied
by path $\nu$ and 0 otherwise).

3:For all $(ij)\in E,\nu\in[M]$ initialize messages $a^{\nu}_{i\rightarrow j}$
and $b^{\nu}_{i\rightarrow j}$ as :  $\displaystyle\begin{aligned}
a^{\nu}_{i\to j}&=\left\\{\begin{array}[]{lr}-U[0,1]\text{ \ \ if $i=s^{\nu}$
},\\\ \infty\text{ \ \ \ \ \ \ \ \ \ \ if $i=t^{\nu}$ },\\\ U[0,1]\text{ \ \ \
\ \ otherwise.}\end{array}\right.\end{aligned}$ $\displaystyle\begin{aligned}
b^{\nu}_{i\to j}&=\left\\{\begin{array}[]{lr}\infty\text{ \ \ \ \ \ \ \ \ if
$i=s^{\nu}$ },\\\ 0\text{\ \ \ \ \ \ \ \ \ \ \ if $i=t^{\nu}$ },\\\
U[0,1]\text{ \ \ \ otherwise.}\end{array}\right.\end{aligned}$

4:$t$ $\leftarrow$ 0 $\triangleright$ counter on number of iterations.

5:$EqCount$ $\leftarrow$ 0 $\triangleright$ counter used to condition on the
convergence.

6:repeat

7: Select a random node $i\in V$ and a random neighbour $j\in\partial i$

8: $pre_{a}^{\nu}$ $\leftarrow$ $a_{i\rightarrow j}^{\nu}\text{\ \ , \ \
}\forall\nu\in[M]$

9: $pre_{b}^{\nu}$ $\leftarrow$ $b_{i\rightarrow j}^{\nu}\text{\ \ , \ \
}\forall\nu\in[M]$ $\triangleright$ Store old messages.

10: for each path $\nu\in[M]$ do

11: $\sigma_{(ij)}^{\nu}\leftarrow$
$F_{1}(pre_{a}^{\nu},pre_{b}^{\nu},\\{a_{r\rightarrow
j}^{\nu}\\}_{r\in\partial i\backslash j},\\{b_{s\rightarrow
j}^{\nu}\\}_{s\in\partial i\backslash j})$ $\triangleright$ assign spins,
based on old messages.

12: $\eta_{ji}^{\nu*}\leftarrow$
$F_{2}(\\{\sigma^{\mu}_{(ji)}\\}_{\mu\in[M]\backslash\nu})$ $\triangleright$
Compute the undirected traffic at edge $(ji)$ ignoring path $\nu$.

13: $a_{i\rightarrow j}^{\nu},b_{i\rightarrow j}^{\nu}\leftarrow$
$F_{3}(\\{a_{r\rightarrow i}^{\nu}\\}_{r\in\partial i\backslash
j},\\{b_{s\rightarrow i}^{\nu}\\}_{s\in\partial i\backslash
j},\phi^{\prime}(\eta_{ji}^{\nu*}))$ $\triangleright$ Update of the messages.

14: end for

15: if $\sum_{\mu=0}^{N}|a_{j\rightarrow
i}^{\mu}-pre_{a}^{\mu}|+|b_{j\rightarrow i}^{\mu}-pre_{b}^{\mu}|<\epsilon$
then $\triangleright$ Check for convergence

16: $EqCount\leftarrow EqCount+1$

17: else

18: $EqCount\leftarrow 0$

19: end if

20: $t\leftarrow t+1$

21:until ($EqCount$ $\geq$ $EqTime$ or $t\geq$ $MaxIter*$card($V$))

22:for each edge $(ij)\in E$ do

23: for each path $\nu\in[M]$ do

24: $\sigma_{(ij)}^{\nu}\leftarrow$ $F_{1}(\\{a_{r\rightarrow
j}^{\nu}\\}_{r\in\partial j},\\{b_{q\rightarrow j}^{\nu}\\}_{q\in\partial j})$

25: end for

26:end for

##  Appendix C Derivation of the SAW sampling algorithm

This appendix we detail the derivation of the MCMC to sample SAWs on a
weighted graph. Let us first set up some notation. We will indicate with
$\pi=(\pi_{1},\dots,\pi_{L(\pi)})$ a path going from $\pi_{1}$ to
$\pi_{L(\pi)}$. $L(\pi)$ indicates the number of nodes in the path. Consider a
graph with edges weighted with the matrix
$W\in\operatorname{\mathbb{R}}^{N\times N}$. $W_{ij}$ represents the cost for
the walk to traverse edge $i\to j$. If the edge $i\to j$ is not present we set
$W_{ij}=\infty$. We suppose $W_{ij}\geq 0$ for all $i,j$. When referring to
shortest paths, we will mean paths with minimum weight with respect to $W$. We
wish to sample from the following probability measure, over the SAWs going
from node $s$ to node $t$.

$P(\pi|s,t)=\frac{1}{Z(s,t)}e^{-\beta
E_{W}(\pi)}\operatorname{\mathbb{I}}[\pi_{1}=s]\operatorname{\mathbb{I}}[\pi_{L(\pi)}=e]\operatorname{\mathbb{I}}[\pi\text{
is a SAW}],\quad
E_{W}(\pi)\coloneqq\sum_{k=1}^{L(\pi)-1}W_{\pi_{k},\pi_{k+1}}$ (37)

When $\beta$ is large, the probability is concentrated around the shortest
path (according to $E$). Instead when $\beta\to 0$ the measure is uniform over
all SAWs. To sample from this measure, we design a Metropolis based sampler.
Let $\pi$ be the current path in the MCMC. We first propose a new path
$\tilde{\pi}$ by sampling it from a proposal probability $Q(\tilde{\pi}|s,t)$
and then we accept the move with probability
$p_{acc}=\frac{P(\tilde{\pi}|s,t)Q(\pi|s,t)}{P(\pi|s,t)Q(\tilde{\pi}|s,t)}$.
Notice that differently from the usual Metropolis algorithm, the proposed path
is independent from the current path. With use the notation $\pi_{:k}$ and
$\pi_{k:}$ we indicate respectively the first $k$ nodes of path $\pi$, and the
nodes from $\pi_{k}$ (included) to the end of $\pi$. As a proposal probability
we use an autoregressive process of the following form

$Q(\tilde{\pi}|s,t)=\prod_{k=1}^{L(\pi)-1}q(\tilde{\pi}_{k+1}|\tilde{\pi}_{:k},s,t).$
(38)

Thanks to the autoregressive property, $Q$ is easy to sample: it is sufficient
to draw iteratively from the distribution $q(\cdot|\cdot)$ until the node $e$
is reached. Given a graph $G=(V,E)$ and a path (a sequence of nodes connected
by edges) $\pi$ in $G$, we denote by $G\backslash\pi$ the graph where we have
removed the nodes in $\pi$ and the edges connected to them. Define also
$V_{e}(y,\pi)$ to be the cost of the shortest path from $y$ to $e$ in the
graph $G\backslash\pi$. If such path does not exist, we define $V_{e}(y,\pi)$
to be $\infty$. we are now ready to define $q(\cdot|\cdot)$.

$q(\tilde{\pi}_{k+1}=y|\tilde{\pi}_{:k},s,t)=\frac{1}{\mathcal{Z}(x_{:k},s,t)}\operatorname{\mathbb{I}}[y\not\in\tilde{\pi}_{:k}]\exp\left[-\beta(W_{\tilde{\pi}_{k}y}+V_{e}(y,\tilde{\pi}_{:k}))\right].$
(39)

In words, supposing we have sampled the first $k$ elements of the path and we
must pick the $k+1$, for each neighbor $y$ of $\tilde{\pi}_{k}$,
$W_{\tilde{\pi}_{k}y}+V_{e}(y,\tilde{\pi}_{:k})$ represents the cost of going
from $x_{k}$ to $y$ and from there continuing along the shortest SAW from $y$
to $e$. As explained in section C.1, this can be seen as a high $\beta$
expansion where one only keeps track of the first order term in $\beta$. The
pseudocode of this algorithm is presented in 4.

###  C.1 SAW sampler as high $\beta$ expansion

In this section we show how the autoregressive process to propose SAWs
described by (39) corresponds to a high $\beta$ expansion. Let $P$ be the
probability measure on SAWs with constrained ends defined in (37). Consider
the following way of writing the conditional probability for $k+1$th node in
the path given the first $k$ nodes, under $P$.

$\displaystyle P(\pi_{k+1}=y|\pi_{:k},s,t)=$ (40)
$\displaystyle=\sum_{\pi_{k+2},\dots,\pi_{L(\pi)}}P(\pi_{k+1:}=[y,\pi_{k+2:}]|\pi_{:k},s,t)=\frac{1}{Z(\pi_{:k},s,t)}\sum_{\pi_{k+2},\dots,\pi_{L(\pi)}}\operatorname{\mathbb{I}}[\pi\text{
is a SAW}]\,e^{-\beta(W_{\pi_{k}y}+E_{W}(\pi_{k+2:}))}$ (41)
$\displaystyle=\frac{1}{Z(\pi_{:k},s,t)}e^{-\beta(W_{\pi_{k}y}+V_{e}(y,\pi_{:k}))}\sum_{\pi_{k+2},\dots,\pi_{L(\pi)}\in
V}\operatorname{\mathbb{I}}[\pi\text{ is a
SAW}]\,e^{-\beta(E_{W}(\pi_{k+2:})-V_{e}(y,\pi_{:k}))}.$ (42)

$Z(\pi_{:k,s,t})$ is as always a normalizing constant. In the second
expression we sum over all the possible continuations of $\pi_{:k+1}$ and
$[y,\pi_{k+2:}]$ indicates the path obtained by concatenating $y$ and
$\pi_{k+2:}$, similarly in the following steps $\pi$ indicates the
concatenation of $\pi_{:k},y,\pi_{k+2:}$. In the last passage we have brought
out of the sum the term $V_{e}(y,\pi_{:k}))$ corresponding to the shortest
path. Taking the limit $\beta\to\infty$ we obtain

$\lim_{\beta\to\infty}\sum_{\pi_{k+2},\dots,\pi_{L(\pi)}\in
V}\operatorname{\mathbb{I}}[\pi\text{ is a
SAW}]\,e^{-\beta(E_{W}(\pi_{k+2:})-V_{e}(y,\pi_{:k}))}=\sum_{\pi_{k+2},\dots,\pi_{L(\pi)}\in
V}\operatorname{\mathbb{I}}[\pi\text{ is a
SAW}]\operatorname{\mathbb{I}}[E_{W}(\pi_{k+2:})=V_{e}(y,\pi_{:k}))]$ (43)

This is the number of shortest paths from $y$ to $e$ in the graph
$G\backslash\pi_{:k}$. If one assumes $W$ to be a random matrix with
continuous entries, then the shortest path is almost surely not degenerate and
thus limit is one. Then it is clear that in the limit $\beta\to\infty$,
$\frac{1}{\beta}\log P(\pi_{k+1}=y|\pi_{:k},s,t)\to\frac{1}{\beta}\log
q(\pi_{k+1}=y|\pi_{:k},s,t)$, and therefore $q(\cdot|\cdot)$ captures the
leading order in $\beta$. In some special cases cases there can be a huge
degeneracy in the number of shortest paths: consider for example the case of
$W$ being the adjacency matrix of a grid.

Algorithm 4 Metropolis SAW sampler for weighted graph

1:Input:  Weight matrix $W\in\mathbb{R}^{N\times N}_{+}$, initial path
$\pi_{0}$, inverse temperature $\beta$, $G$ graph over which the path is
sampled, $s,t$ respectively the origin and destination of the path,
$t_{\text{max}}$ number of MCMC steps.

2:Output: a list $S$ of paths sampled from the Gibbs distribution (37)

3:$S=[\pi_{0}]$

4:$\pi\leftarrow\pi_{0}$

5:$K\leftarrow\log Q(\pi|s,t)$

6:for $k=1,\dots,t_{\text{max}}$ do

7: $\tilde{\pi}\sim Q(\cdot|s,t)$ $\triangleright$ Sampling proposed path
using 5. $Q$ depends implicitly on $W,G,\beta$.

8: $\tilde{K}\leftarrow\log Q(\tilde{\pi}|s,t)$

9:
$p_{\text{acc}}\leftarrow\exp\left[K-\tilde{K}-\beta(E_{W}(\tilde{\pi})-E_{W}(\pi))\right]$

10: $U\sim\text{Uniform}([0,1])$

11: if $U<p_{\text{acc}}$ then $\triangleright$ Accept-reject step

12: $\pi\leftarrow\tilde{\pi}$

13: $K\leftarrow\tilde{K}$

14: end if

15: $S.\text{append}(\pi)$

16:end for

Algorithm 5 Proposal algorithm for SAW sampler

1:Input:  Weight matrix $W\in\mathbb{R}^{N\times N}_{+}$, inverse temperature
$\beta$, $G$ directed weighted graph over which paths are sampled, $(s,t)$
respectively the origin and destination of the path.

2:Output: a SAW $\pi$ sampled from the proposal distribution (38)

3:$\pi\leftarrow[s]$

4:$x\leftarrow s$ $\triangleright$ $x$ is the current endpoint of the path

5:while $x\neq t$ do $\triangleright$ loop until the end node is reached

6: for $z\in\partial_{in}x$ do$\triangleright$ $\partial_{in}x$ is the set of
nodes $z$ such that the edge $z\to x$ exists in $G$

7: $G$.removeEdge($z,x$) $\triangleright$ Removing all incoming edges into
$x$. This guarantees that the path is self avoiding

8: end for

9: probnextnode $\leftarrow[\;]$

10: nextnodes $\leftarrow[\;]$ $\triangleright$ list of possible next nodes

11: for $y\in\partial_{out}x$ do $\triangleright$ $\partial_{out}x$ is the set
of nodes $z$ such that the edge $x\to z$ exists in $G$

12: nextnodes.append($y$)

13: $V\leftarrow$ shortestPath($G,W,y,t$)

14: probnextnode.append($e^{-\beta(V+W_{xy})}$)

15: end for

16: probnextnode $\leftarrow$ probnextnode/sum(probnextnode) $\triangleright$
normalizing the probability

17: $i\sim$ probnextnode $\triangleright$ sampling a index according to
probnextnode

18: $x\leftarrow$ nextnodes[i]

19: $\pi.$append(x)

20:end while

The function shortestPath($G,W,y,t$) returns
$\min_{\pi:\pi_{1}=y,\pi_{L(\pi)=t}}E_{W}(\pi)$, where the minimum is taken
over all the paths in graph $G$. Numerically the normalization step is tricky,
since all the exponentials can be very negative, this might result in a $0/0$
division. To avoid this it is better to compute probabilities as using the
formula
$e^{-x_{i}}/\sum_{j}e^{-x_{j}}=e^{-(x_{i}-x_{1})}/(1+\sum_{j=2}e^{-(x_{j}-x_{1})})$,
where we supposed $x_{1}=\min_{i}x_{i}$. In this way the expression is
regularized. Notice that with a small variation, algorithm 5 can return also
$Q(\pi|s,t)$, which is needed in the Metropolis step. One initializes $\log
Q=0$ and every time a node is added one takes $\log Q=\log
Q+\log($probnextnode[i]$)$, where is the index sampled from probnextnode the
index of the next node in the path. Notice that in 5 we always assume the
graph to be directed. If we are dealing with an undirected graph, then we just
transform it into a directed one by replacing each undirected edge with two
opposite directed ones.

##  Appendix D Derivation of the Simulated annealing algorithm

In this appendix we give the pseudocode of the simulated annealing algorithm
used to find a solution to ITAP. As annealing schedule we use

$\beta(t)=\begin{cases*}\beta_{\text{min}}\frac{T_{\text{ann}}}{T_{\text{ann}}-t}&if
$t<T_{\text{ann}}$\\\ \infty&if $t\geq T_{\text{ann}}$\end{cases*}$ (44)

In words, we start from $\beta_{\text{min}}$, then increase $\beta$, until it
diverges after $T_{\text{ann}}$ steps, after that we keep $\beta=\infty$ until
the iterations converge.

Algorithm 6 Simulated annealing for ITAP

1:Input:  Initial paths $\\{\pi_{0}^{\mu}\\}_{\mu\in[M]}$, OD pairs
$\\{s^{\mu},t^{\mu}\\}_{\mu\in[M]}$ inverse temperature schedule $\beta(t)$,
$G$ graph over which paths are sampled, $t_{\text{singlepath}}$ number of
Metropolis steps to make when sampling from (5)

2:Output: $\\{\pi^{\mu}\\}_{\mu\in[M]}$ a set of paths to which the algorithm
converges

3:Gflow$\leftarrow$ copy($G$) $\triangleright$ the weight on edge $e$ of Gflow
is $I_{e}$

4:Gcost $\leftarrow$ copy($G$) $\triangleright$ the weight on edge $e$ of
Gcost is $\phi(I_{e})$

5:for $e\in G.$edges do

6: Gflow.weight($e$)$\leftarrow$0

7:end for

8:for $\mu=1,\dots,M$ do

9: $\pi^{\mu}\leftarrow\pi_{0}^{\mu}$

10: for $k=1,\dots,$ length($\pi^{\mu})-1$ do

11: $e\leftarrow[\pi^{\nu}[k],\pi^{\nu}[k+1]]$

12: Gflow.weight($e$)$\leftarrow$ Gflow.weight($e$)+1 $\triangleright$
populating the flow graph with paths

13: end for

14:end for

15:H$\leftarrow 0$ $\triangleright$ current value of $H(\pmb{I})$

16:for $e\in G.$edges do

17:
Gcost.weight($e$)$\leftarrow\phi(\text{Gflow.weight($e$)}+1)-\phi(\text{Gflow.weight($e$)})$

18: H$\leftarrow$ H+Gcost.weight($e$)

19:end for

20:flagconv $\leftarrow$ False

21:while not flagconv do$\triangleright$ main loop

22: for $\nu=1,\dots,M$ do$\triangleright$ loop over all paths

23: for $k=1,\dots,\text{length($\pi^{\nu}$)}-1$ do

24: $e\leftarrow[\pi^{\nu}[k],\pi^{\nu}[k+1]]$

25: Gflow.weight($e$)$\leftarrow$ Gflow.weight($e$)-1$\triangleright$ removing
path $\nu$ from the graph so that one has $I^{\backslash\nu}$

26:
Gcost.weight($e$)$\leftarrow\phi(\text{Gflow.weight($e$)}+1)-\phi(\text{Gflow.weight($e$)})$

27: end for

28: $W\leftarrow$ weightMatrix(Gcost)

29: $\pi^{\nu}\leftarrow$
SAWsampler($W,\pi^{\nu},\beta(t),G,s^{\nu},t^{\nu},t_{\text{singlepath}}$)$\triangleright$
Run algorithm 4 for $t_{\text{singlepath}}$ steps to sample from the
conditional (5)

30: for $k=1,\dots,\text{length($\pi^{\nu}$)}-1$ do

31: $e\leftarrow[\pi^{\nu}[k],\pi^{\nu}[k+1]]$

32: Gflow.weight($e$)$\leftarrow$ Gflow.weight($e$)+1$\triangleright$ adding
back path $\nu$ to the graph

33:
Gcost.weight($e$)$\leftarrow\phi(\text{Gflow.weight($e$)}+1)-\phi(\text{Gflow.weight($e$)})$

34: end for

35: end for

36: newH$\leftarrow 0$ $\triangleright$ measuring energy after all paths are
updated

37: for $e\in G.$edges do

38: newH$\leftarrow$newH+Gcost.weight($e$)

39: end for

40: if H == newH and $\beta==\infty$ then

41: flagconv$\leftarrow$True $\triangleright$ If $\beta=\infty$ and
$H(\pmb{I})$ does not change once all the paths are updated, then the
algorithm has converged

42: end if

43: H $\leftarrow$ newH

44:end while

The function ’weightMatrix(Gcost)’ returns an $N\times N$ (with $N$ being the
number of vertices) matrix $W$ where $W_{ij}=$Gcost.weight($[i,j]$) if edge
$i\to j$ exists in $G$ and $W_{ij}=\infty$ otherwise. Notice that the
pseudocode of this algorithm is almost identical to that of the greedy
algorithm 2. The only difference is in lines 28, 29. Erasing line 28 and
replacing 29 with ’$\pi^{\nu}\leftarrow$ shortestPath(Gcost,
$s^{\nu},t^{\nu}$)’ (line 28 in 2) one obtains again the greedy algorithm.
This relation sheds some light onto the relation between the two algorithms.
While the greedy algorithm always picks the shortest path (weighted with
$\Delta H$), the simulated annealing is noisy and samples a path from $\propto
e^{-\beta\Delta H}$.

##  Appendix E RITAP details

The RITAP algorithm 1 is composed of two parts:

  * •

An algorithm named $\verb|solveTAP|(G,D)$ that outputs a solution to TAP in
the form of a collection
$\operatorname{\mathcal{P}}=\\{\operatorname{\mathcal{P}}_{xy}\\}$ with
$\operatorname{\mathcal{P}}_{xy}$ defined in (20).

  * •

A projection algorithm that takes as input a TAP solution and outputs an ITAP
solution, with the same $\operatorname{\mathcal{P}}^{\text{ITAP}}$.

The next two sections focus respectively on the implementation of
$\verb|solveTAP|(G,D)$, and on the integer projection step.

###  E.1 Frank-Wolfe as TAP solver

To solve TAP we use Frank-Wolfe (FW) algorithm. FW is an iterative algorithm
that produces a sequence of edge flows $\\{\pmb{I}(t)\\}$ that approaches the
optimal edge flows $\pmb{I}^{\star}$. Each step of FW consists of two sub
steps:

  1. 1.

Linearizing $H(\pmb{I}(t))$ near the current point and minimizing the
linearized objective. Denoting with $\tilde{H}$ this linear approximation we
have $\tilde{H}(\pmb{I};\pmb{I}(t))=H(\pmb{I}(t))+\sum_{e\in
E}\phi^{\prime}(I_{e}(t))(I_{e}-I_{e}(t))$. $\tilde{H}(\pmb{I};\pmb{I}(t))$ is
minimized with respect to $\pmb{I}$ under the TAP constraints. Since
$\tilde{H}$ is linear in $\pmb{I}$, at its minimum for each pair of nodes
$x,y$, all the flow from $x$ to $y$ is routed through a shortest path in the
graph whose edges are weighted with $\phi^{\prime}(I_{e}(t))$. The new flow
$\tilde{\pmb{I}}(t)$ obtained from this minimization is called the all-or-
nothing flow, since between every pair of nodes all the flow is routed through
a single path.

  2. 2.

Step towards the all-or-nothing solution: We move from $\pmb{I}(t)$ towards
the all-or-nothing solution with a step size $\eta(t)$.
$\pmb{I}(t+1)=(1-\eta(t))\pmb{I}(t)+\eta(t)\pmb{\tilde{\pmb{I}}}(t)$. The step
size is taken to be decreasing in $t$ with an appropriate schedule that
guarantees convergence.

When $H(\pmb{I})$ is convex, FW gives, at each step, a lower bound
$H_{\text{LB}}(t)$ to the optimal energy. This allows to assess convergence of
the algorithm by measuring $H(\pmb{I}(t))-H_{\text{LB}}(t)$.

The following is a one line proof of the fact that
$H_{\text{LB}}(t)=\min_{\tilde{\pmb{I}}}\tilde{H}(\pmb{I},\pmb{I}(t))$
999where the minimization over $\tilde{\pmb{I}}$ is carried out under the TAP
constraints is a lower bound to $H(\pmb{I}^{\star})$

$H(\pmb{I}^{\star})\geq
H(\pmb{I}(t))+\sum_{e}\phi^{\prime}(I_{e})(I_{e}^{\star}-I_{e}(t))\geq\min_{\tilde{\pmb{I}}}H(\pmb{I}(t))+\sum_{e}\phi^{\prime}(I_{e})(\tilde{I}_{e}-I_{e}(t))=\min_{\tilde{\pmb{I}}}\tilde{H}(\pmb{I},\pmb{I}(t)),$
(45)

where the minimization over $\tilde{\pmb{I}}$ is carried out under the TAP
constraints, and the first inequality holds by convexity. Notice that
$\tilde{H}(\pmb{I},\pmb{I}(t))$ is minimized at every time step, therefore the
bound is readily available.

In our experiments we run FW with the following specifications:

  * •

The edge flows at initialization $\pmb{I}(t=0)$ are obtained by computing the
shortest paths (in the topological distance on $G$) between every pair of
nodes $x,y$ with positive demand. Denoting with $\pi_{xy}$ the a topological
shortest path between $x$ and $y$, we have $I(t=0)_{e}=\sum_{x,y\in(V\times
V)}D_{xy}\mathbb{I}[e\in\pi_{xy}]$.

  * •

Supposing we run FW for $t_{\text{max}}$ steps, our step size schedule is
$\eta(t)=t^{-\alpha(t)}$ with

$\alpha(t)=\begin{cases}\frac{1}{2}(1+\frac{t}{t_{\text{switch}}})&\text{if
$t\leq t_{\text{switch}}$}\\\ 1.35&\text{if
$t>t_{\text{switch}}$},\end{cases}$ (46)

with $t_{\text{switch}}=0.9\,t_{\text{max}}$. We picked this schedule since we
observed it works well in practice. The idea behind it is to start with a step
size $\eta(t)=1/\sqrt{t}$ and gradually move towards a rate $\eta(t)=1/t$. The
last 10% of iterations we decrease to $t^{-1.35}$, this has the effect of
obtaining a better lower bound to the TAP optimal energy. Without worrying
about the lower bound the schedule
$\alpha(t)=\frac{1}{2}(1+\frac{t}{t_{\text{max}}})$ is found to work very well
and leads to $H(\pmb{I}(t))/H(\pmb{I}^{\star})-1\approx 1/t^{3}$.

  * •

The classic FW algorithm outputs only the optimal edge flows, and does not
return the optimal paths. We modify FW to keep track, throughout the
iterations, of the paths on which flow is allocated. The paths together with
the corresponding amount of flow routed over them are returned.

Notice that the FW algorithm does not depend explicitly on $M$, the number of
OD pairs. $M$ is instead reflected in the magnitude of the entries of $D$.
Thanks to this independence FW is particularly efficient for very large $M$.

###  E.2 ITAP projection

Once we have a TAP solution we must project it onto a set of path flows
satisfying the integrality constraint. RITAP as presented in 1 implements one
such procedure outputting the paths $\\{\pi^{\mu}\\}_{\mu\in[M]}$. Crucially
this procedure complexity is linear in $M$. Below in 7 we provide a different
algorithm that achieves the same result, but whose dependence on the number of
paths goes like $\min(M,N^{2})$.

Algorithm 7 Relaxed ITAP solver (RITAP)

1:Input: $G$ graph over which paths are sampled, $D$ demand matrix.

2:Output: $\operatorname{\mathcal{P}}^{\text{ITAP}}$ collection of sets in the
form 20 containing the paths flows satisfying ITAP constraints.

3:$\operatorname{\mathcal{P}}\leftarrow\verb|solveTAP|(G,D)$. $\triangleright$
Solving TAP

4:$\operatorname{\mathcal{P}}^{\text{ITAP}}\leftarrow\text{copy}(\operatorname{\mathcal{P}})$

5:$\operatorname{\mathcal{P}}^{\text{tmp}}\leftarrow\text{copy}(\operatorname{\mathcal{P}})$

6:for $x\in G.$nodes do

7: for $y\in G.$nodes do

8: if $D_{xy}>0$ then

9: allocatedFlow$\leftarrow 0$

10: for $a\in[1,2,\dots,R_{xy}]$ do

11:
$\operatorname{\mathcal{P}}^{\text{ITAP}}_{xy}.h[a]\leftarrow\lfloor\operatorname{\mathcal{P}}^{\text{ITAP}}_{xy}.h[a]\rfloor$$\triangleright$
Keeping only the integer part of the flow on each path

12: allocatedFlow$\leftarrow$ allocatedFlow
$+\operatorname{\mathcal{P}}^{\text{ITAP}}_{xy}.h[a]$

13:
$\operatorname{\mathcal{P}}^{\text{tmp}}_{xy}.h[a]\leftarrow\operatorname{\mathcal{P}}^{\text{tmp}}_{xy}.h[a]$
mod 1 $\triangleright$ Keeping only the fractional part

14: end for

15: sortedFracFlows, permIdx $\leftarrow$
sortTogether($\operatorname{\mathcal{P}}^{\text{tmp}}_{xy}.h,[1,2,\dots,R_{xy}]$)
$\triangleright$ Sort by fractional part

16: for $k\in[0,1,\dots,D_{xy}-\text{allocatedFlow}]$ do $\triangleright$
Allocating the remaining flow in decreasing order of fractional flow

17:
$\operatorname{\mathcal{P}}^{\text{ITAP}}_{xy}.h[\text{permIdx[$k$]}]\leftarrow\operatorname{\mathcal{P}}^{\text{ITAP}}_{xy}.h[\text{permIdx[$k$]}]+1$

18: end for

19: end if

20: end for

21:end for

22:return  $\operatorname{\mathcal{P}}^{\text{ITAP}}$

the function sortTogether takes two lists as arguments, it sorts the first
list in descending order, and applies the same permutation to the second list.
For example SortTogether($[3,4,1,5,2],[a,b,c,d,e]$)=$[5,4,3,2,1],[d,b,a,e,c]$.
Basically this algorithm builds $\operatorname{\mathcal{P}}^{\text{ITAP}}$ by
first keeping the integer part of the TAP path flows and then adding the
remaining integer units of flow, on the paths with the highest fractional
part.

##  Appendix F Additional numerical experiments in the regime
$M=\Theta(N^{2})$

In this section we present some complementary results to those shown in
section III.2 and IV.2. We consider the setting of section III.2, consisting
of undirected RRG with $d=3$, OD pairs sampled uniformly at random and
$\phi(x)=x^{2}$. We run experiments each time changing one of these factors
and leaving the others unaltered. Below is a list of the modified settings.

  * •

Considering RRG with higher degrees, instead of $d=3$

  * •

Considering real world traffic networks instead of RRGs

  * •

Considering a realistic Bureau of Public Roads (BPR) nonlinearity instead of
$\phi(x)=x^{\gamma}$

###  F.1 Varying the degree

We start by presenting data about RRGs with degree larger than $3$ and fixed
size $N=128$. Increasing the degree alters some qualitative characteristics
previously observed for $d=3$ in figure 7. However crucially the asymptotic
behavior is robust to changes in $d$. Panel (a) and (b) show that the gap
between the TAP optimal energy and the one outputted by ITAP goes to zero as
$\eta$ increases. In panel (c) we observe that the onset of the power law
behavior for $1-F_{\text{int}}$ is displaced when changing $d$, this is
potentially because the bound (24) becomes non-vacuous (and hence the $1/\eta$
is expected to arise) when $\eta\geq(K-1)$. Since $K-1$ varies from 0.1 for
$d=3$ to 10 for $d=48$, the knee in panel (c) moves accordingly. Ultimately
for large $\eta$ all curves display the predicted $1/\eta$ rate. Similarly,
for all degrees, the average degeneracy $K$ (plotted in (d)) plateaus when
$\eta$ increases. The two-scale (i.e. $\rho,\eta$) behavior also persists when
changing $d$. For $\rho=\Theta(1)$ the curves collapse when plotted against
$\rho$ and for $\eta\geq 1$ (approximately) one enters in the asymptotic
regime. In summary, all of the claims we made about the model still hold when
looking as experiments with higher $d$.

Figure 10: Experiments for $N=128,\,\gamma=2$, on RRG with varying degree.
From darker to lighter the degree is $2,6,12,24,48$. Results are averaged over
20 independent realizations of the RRG and of the OD pairs. Shaded areas
represent a $\pm 1$ standard deviation interval around the mean. Insets show
the same quantity plotted in the main panel but as a function of $\rho$. (a)
The upper and lower curves correspond respectively to TAP optimal solution and
the RITAP solution to ITAP. (b) Ratio of the energies reached respectively by
RITAP and TAP. (c) Fraction of integer paths as a function of $\eta$. The
dashed line has a slope of $1/\eta$.(d) Average degeneracy as a function of
$\eta$

###  F.2 Varying the nonlinearity

In traffic the assignment literature the time a user takes to traverse an edge
with total traffic $I$ is usually modeled using the Bureau of Public Roads
(BPR) functions.

$f(I)=t_{0}\left(1+\alpha\left(\frac{I}{c}\right)^{\beta}\right),$ (47)

where $t_{0}$ is the crossing time in absence of traffic, $c$ is the capacity
of the edge; when the traffic exceeds the capacity the cost quickly rises. In
our experiments we take $t_{0}=1,\alpha=3.8\times 10^{-2},\beta=4,c=1$.

Since we are interested in minimizing the average cost we set
$\phi(x)=xf(x)$101010recall that $\phi$ represents the total cost of one edge.
Therefore it is obtained by multiplying the individual cost $f$ times the flow
on the edge.. Results of thee experiments are presented in figure 11. Each
line is the result of averaging over 20 independent realizations of the OD
pairs and the RRG. We notice an interesting difference with the
$\phi(x)=x^{\gamma}$ case: in panel (a), for small $\eta$ (small $M$), the TAP
and RITAP curves are close. This is in contrast with the $\phi(x)=x^{\gamma}$
case, where the gap between TAP and RITAP was maximum for small $M$. The
reason behind this difference lies in the fact that the $xf(x)$ has a linear
behavior near zero, therefore when the flow per edge is small, the flow is
routed on shortest paths in the topological distance, instead of being split
across many paths. In other words, having a nonlinearity that has a nonzero
derivative near $I_{e}=0$ results in TAP support paths that are shortest paths
in the topological distance when $M$ is small. Thanks to this fact the
degeneracy at small $M$ is quite small compared to the case $\gamma=2$ (see
figure 7 (d)). Apart from this difference in the small $M$ regime, the
asymptotic behavior is unchanged. As in the case $\gamma=2$, when $\eta\gg 1$,
$H_{\text{ITAP}}/H_{\text{TAP}}\to 1$, $1-F_{\text{int}}\to 0$ as $1/\eta$,
$K$ approaches a constant value. We conclude that the asymptotic phenomenology
is robust to changes in the nonlinearity. One caveat to this statement
consists in the choice of the capacity parameter. Our experiments suggest that
if $c$ is taken to be large (e.g. a few hundreds), then this introduces an
additional scale into the problem. The new scale alters the two regimes
($\rho,\eta$) phenomenology we so far encountered. Finally we remark that also