# Routing in an Uncertain World:
Adaptivity, Efficiency, and Equilibrium

Dong Quan Vu∗ ∗ Univ. Grenoble Alpes, CNRS, Inria, LIG, 38000, Grenoble,
France<EMAIL_ADDRESS>,  Kimon Antonakopoulos∗
<EMAIL_ADDRESS>and  Panayotis Mertikopoulos∗,♯ ♯ Criteo AI Lab.
<EMAIL_ADDRESS>
###### Abstract.

We consider the traffic assignment problem in nonatomic routing games where
the players’ cost functions may be subject to random fluctuations (e.g.,
weather disturbances, perturbations in the underlying network, etc.). We
tackle this problem from the viewpoint of a control interface that makes
routing recommendations based solely on observed costs and without any further
knowledge of the system’s governing dynamics – such as the network’s cost
functions, the distribution of any random events affecting the network, etc.
In this online setting, learning methods based on the popular exponential
weights algorithm converge to equilibrium at an
$\operatorname{\mathcal{O}}(1/\sqrt{T})$ rate: this rate is known to be order-
optimal in stochastic networks, but it is otherwise suboptimal in static
networks. In the latter case, it is possible to achieve an
$\operatorname{\mathcal{O}}(1/T^{2})$ equilibrium convergence rate via the use
of finely tuned accelerated algorithms; on the other hand, these accelerated
algorithms fail to converge altogether in the presence of persistent
randomness, so it is not clear how to achieve the “best of both worlds” in
terms of convergence speed. Our paper seeks to fill this gap by proposing an
adaptive routing algortihm with the following desirable properties:  (it n) it
seamlessly interpolates between the $\operatorname{\mathcal{O}}(1/T^{2})$ and
$\operatorname{\mathcal{O}}(1/\sqrt{T})$ rates for static and stochastic
environments respectively;  (it n) its convergence speed is polylogarithmic in
the number of paths in the network;  (it n) the method’s per-iteration
complexity and memory requirements are both linear in the number of nodes and
edges in the network; and  (it n) it does not require any prior knowledge of
the problem’s parameters.

###### Key words and phrases:

Nonatomic congestion games; traffic equilibrium; adaptive methods; exponential
weights; uncertainty

###### 2010 Mathematics Subject Classification:

Primary 91A10, 91A26; secondary 68Q32, 68T02.

The authors are grateful for financial support by the French National Research
Agency (ANR) in the framework of the “Investissements d’avenir” program
(ANR-15-IDEX-02), the LabEx PERSYVAL (ANR-11-LABX-0025-01), MIAI@Grenoble
Alpes (ANR-19-P3IA-0003), and the grant ALIAS (ANR-19-CE48-0018-01). The
authors’ research was also supported by the COST Action CA16228 “European
Network for Game Theory” (GAMENET)

##  1\. Introduction

Transportation networks in major metropolitan areas carry several million car
trips and commutes per day, giving rise to a chaotic and highly volatile
environment for the average commuter. As a result, navigation apps like Google
Maps, Waze and MapQuest have seen an explosive growth in their user base,
routinely receiving upwards of $10^{4}$ routing requests per second during
rush hour – and going up to $10^{5}$ queries/second in the largest cities in
the US and China [10]. This vast number of users must be routed efficiently,
in real-time, and without causing any “ex-post” regret at the user end;
otherwise, if a user could have experienced better travel times along a non-
recommended route, they would have no incentive to follow the app
recommendation in the first place. In the language of congestion games [35],
this requirement is known as a “Wardrop equilibrium”, and it is typically
represented as a high-dimensional vector describing the traffic flow along
each path in the network [44].

Ideally, this traffic equilibrium should be computed _before_ making a routing
recommendation in order to minimize the number of disgruntled users in the
system. In practice however, this is rarely possible: the state of the network
typically depends on random, unpredictable factors that may vary considerably
from one epoch to the next – e.g., due to road incidents, rain, fog and/or
other weather conditions – so it is generally unrealistic to expect that such
a recommendation can be made in advance. Instead, it is more apt to consider
an online recommendation paradigm that unfolds as follows:

  1. i.5pt)

At each time slot $t=1,2,\dotsc\,,T$, a control interface – such as Google
Maps – determines a traffic assignment profile and routes all demands received
within this slot according to this profile.

  2. ii.5pt)

The interface observes and records the travel times of the network’s users
within the time slot in question; based on this feedback, it updates its
candidate profile for the next epoch and the process repeats.

Of course, the crux of this paradigm is the optimization algorithm used to
update traffic flow profiles from one epoch to the next. Any such algorithm
would have to satisfy the following _sine qua non_ requirements:

  1. i.5pt)

_Universal convergence rate guarantees in terms of the number of epochs:_ The
distance from equilibrium of the candidate traffic assignment after $T$ epochs
should be as small as possible in terms of $T$. However, learning methods that
are well-suited for rapidly fluctuating environments may be too slow in static
environments (i.e., when cost functions do not vary over time); conversely,
methods that are optimized for static environments may fail to converge in the
presence of randomness. As a result, it is crucial to employ _universal_
algorithms that achieve the “best of both worlds” in terms of convergence
speed.

  2. ii.5pt)

_Fast convergence in terms of the size of the network:_ The algorithm’s
convergence speed must be at most polynomial in the size $G$ of the underlying
graph (i.e., its number of edges $E$ plus the number of vertices $V$). In
turn, since the number of paths $P$ in the network is typically exponential in
$G$, the algorithm’s convergence speed must be at most polylogarithmic in $P$.

  3. iii.5pt)

_Scalable per-iteration complexity:_ An algorithm can be implemented
efficiently only if the number of arithmetic operations and total memory
required at each iteration remains scalable as the network grows in size. In
practice, this means that that the algorithm’s per-itration complexity must
not exceed $\operatorname{\mathcal{O}}(G)$.

  4. iv.5pt)

_Parameter-freeness:_ In most practical situations, the parameters of the
model – e.g., the distribution of random events or the smoothness modulus of
the network’s cost functions – cannot be assumed known, so any routing
recommendation algorithm must likewise not require such parameters as input.

To put these desiderata in context, we begin below by reviewing the most
relevant works in this direction.

### Related work.

The _static_ regime of our model matches the framework of Blum et al. [8] who
showed that a variant of the _exponential weights_ (ExpWeight) algorithm [43,
28, 3, 5] converges to equilibrium at an $\operatorname{\mathcal{O}}(\log
P\big{/}\sqrt{T})$ rate (in the sense of time averages). This result was
subsequently extended to routing games with stochastic cost functions by
Krichene et al. [22, 21], who showed that exponential weights (ExpWeight) also
enjoys an $\operatorname{\mathcal{O}}(\log P\big{/}\sqrt{T})$ convergence rate
to _mean_ Wardrop equilibria (again, in the Cesàro sense). However, if the
learning rate of the ExpWeight algorithm is not chosen appropriately in terms
of $T$, the method may lead to non-convergent, chaotic behavior, even in
symmetric congestion games over a $2$-link Pigou network [37].

From an optimization viewpoint, nonatomic congestion games with stochastic
cost functions correspond to convex minimization problems with stochastic
first-order oracle feedback. In this setting, the
$\operatorname{\mathcal{O}}(1/\sqrt{T})$ rate is, in general, unimprovable
[31, 9], so the guarantee of Krichene et al. [22] is order-optimal. On the
other hand, in the static regime, the $\operatorname{\mathcal{O}}(1/\sqrt{T})$
convergence speed of Blum et al. [8] is _not_ optimal: since nonatomic
congestion games admit a smooth convex potential – known as the
_Beckmann–McGuire–Winsten_ (BMW) potential [6] – this rate can be improved to
$\operatorname{\mathcal{O}}(1/T^{2})$ via the seminal “accelerated gradient”
algorithm of Nesterov [30]. On the downside, if applied directly to our
problem, the algorithm of Nesterov [30] has a catastrophic $\Theta(P)$
dependence on the number of paths; however, by coupling it with a “mirror
descent” template in the spirit of [33, 45], Krichene et al. [22] proposed an
accelerated method with an exponential projection step that is particularly
well-suited for congestion problems. In fact, as we show in the sequel, it is
possible to design an accelerated exponential weights method – which we call
AcceleWeight – that achieves an $\operatorname{\mathcal{O}}(\log
P\big{/}T^{2})$ rate in static environments.

Importantly, despite its optimality in the static regime, accelerated
exponential weights (AcceleWeight) fails to converge altogether in stochastic
problems; moreover, the method’s step-size must also be tuned with prior
knowledge of the problem’s smoothness parameters (which are not readily
available to the optimizer). The _universal primal gradient descent_ (UPGD)
algorithm of Nesterov [34] provides a work-around to resolve the latter issue,
but it relies on a line-search mechanism that cannot be applied to stochastic
problems, so it does not resolve the former. Instead, a partial solution for
the stochastic case is achieved by the accelerated stochastic approximation
(AC-SA) algorithm of Lan [25] which achieves order-optimal rates in both the
static and stochastic regimes; however:  itn) the running time $T$ of the AC-
SA algorithm must be fixed in advance as a function of the accuracy threshold
required; and  itn) AC-SAfurther assumes full knowledge of the smoothness
modulus of the game’s cost functions (which cannot be computed ahead of time).

To the best of our knowledge, the first parameter-free algorithm with optimal
rate interpolation guarantees is the AcceleGrad algorithm of Levy et al. [27],
which was originally developed for _unconstrained_ problems (and assumes
knowledge of a compact set containing a solution of the problem). The UnixGrad
proposal of Kavis et al. [20] subsequently achieved the desired adaptation in
constrained problems, but under the requirement of a bounded Bregman diameter.
This requirement rules out the ExpWeight template (the simplex has infinite
diameter under the entropy regularizer that generates the ExpWeight
algorithm), so the convergence speed of UnixGrad ends up being _polynomial_ in
the number of paths in the network – and since the latter scales exponentially
with the size of the network, UnixGrad is unsuitable for networks with more
than $30$ or so nodes.

Finally, another major challenge in terms of scalability is the algorithm’s
_per-iteration complexity_ , i.e., the memory and processing requirements for
making a single update. In regard to this point, the per-iteration complexity
of the universal primal gradient descent (UPGD), AC-SA, AcceleGrad and
UnixGrad algorithms is linear in the number of paths $P$ in the network, and
hence _exponential_ in the size of the network. The only exception to this
list is the (non-accelerated) ExpWeight algorithm, which can be implemented
efficiently with $\operatorname{\mathcal{O}}(G)$ operations per iteration via
a technique known as “weight-pushing” [42, 18]. However, the weight-pushing
technique is by no means universal: it was specifically designed for the
ExpWeight algorithm, and it is not applicable to any of the universal methods
discussed above (precisely because of the acceleration mechanism involved).

### Our contributions.

In summary, all existing traffic assingment algorithms and methods fail in at
least two of the axes mentioned above: either they are slow / non-convergent
outside the specific regime for which they were designed, or they cannot be
implemented efficiently (see also Table 1 for an overview). Consequently, our
paper focuses on the following question:

Is it possible to design a scalable, paramater-free traffic assignment
algorithm
which is simultaneously order-optimal in both static and stochastic
environments?

To provide a positive answer to this question, we take a two-step approach.
Our first contribution is to design an adaptive exponential weights algorithm
– dubbed AdaWeight – which is simultaneously order-optimal, in both $T$ and
$P$, and in both static and stochastic environments. Informally, we have:

######  Informal Version of Theorem 0.

The adaptive exponential weights (AdaWeight) algorithm enjoys the following
equilibrium convergence guarantees after $T$ epochs:

  1. i.5pt)

In static networks, AdaWeight converges to a Wardrop equilibrium at a rate of
$\operatorname{\mathcal{O}}\big{(}(\log P)^{3/2}\big{/}T^{2}\big{)}$.

  2. ii.5pt)

In stochastic networks, it converges to a mean Wardrop equilibrium at a rate
of $\operatorname{\mathcal{O}}\big{(}(\log P)^{3/2}\big{/}\sqrt{T}\big{)}$.

In the above, the logarithmic dependency on $P$ ensures that the convergence
speed of AdaWeight is polynomial – and in fact, _subquadratic_ – in $G$. In
this way, AdaWeight successfully solves the challenge of achieving optimal
convergence rates in both static and stochastic environments, while remaining
parameter-free and “anytime” (i.e., there is no need to tune the algorithm in
terms of $T$). To the best of our knowledge, this is the first method that
simultaneously achieves these desiderata.

On the other hand, a crucial drawback of AdaWeight is that it operates at the
path level, i.e., it updates at each stage a state variable of dimension $P$;
as a result, the algorithm’s per-iteration complexity is typically
_exponential_ in $G$. To overcome this, we propose an improved algorithm,
which we call AdaLight (short for _adaptive local weights_ algorithm), and
which emulates the AdaWeight template with two fundamental differences: First,
instead of maintaining a state variable _per path_ , adaptive local weights
(AdaLight) maintains a state vector _per node_ , with dimension equal to the
node’s out-degree. Second, the recommended path for any particular user of the
navigation interface is constructed “bottom-up” via a sequence of routing
probabilities updated at each node, so there is no need to ever keep in memory
(or update) a path variable.

This process looks similar to the weight-pushing technique of [42, 18] used to
scale down the per-iteration complexity of the ExpWeight algorithm. However,
the acceleration mechanism of AdaWeight involves several moving averages that
_cannot_ be implemented via weight-pushing; and since these averaging steps
are downright essential for the universality of AdaWeight, they cannot be
omitted from the algorithm. Instead, our solution for executing these
averaging steps is based on the following observation: although the averaged
recommendations profiles cannot be “weight-pushed”, the induced edge-load
profiles (i.e., the mass of traffic induced on each edge) can. This
observation allows us to introduce a novel dynamic programming subroutine –
called PushPullMatch – which computes efficiently these averaged loads and
then uses them to reconstruct a set of routing recommendations that are
consistent with these loads.

In this way, by combining the AdaWeight blueprint with the PushPullMatch
subroutine, AdaLight only requires $\operatorname{\mathcal{O}}(G)$ memory and
processing power per update, all the while retaining the sharp convergence
properties of the AdaWeight mother scheme. In more detail, we have:

######  Informal Version of Theorem 0.

The AdaLight algorithm converges to equilibrium at a rate of
$\operatorname{\mathcal{O}}\big{(}(\log P)^{3/2}\big{/}T^{2}\big{)}$ in static
environments and $\operatorname{\mathcal{O}}\big{(}(\log
P)^{3/2}\big{/}\sqrt{T}\big{)}$ in stochastic environments; moreover, the
total amount of arithmetic operations required per epoch is
$\operatorname{\mathcal{O}}(G)$ where $G$ is the size of the network.

This result shows that AdaLight successfully meets the desiderata stated
earlier, so AdaWeight might appear redundant. However, it is not possible to
derive the equilibrium convergence properties of AdaLight without first going
through AdaWeight, so, to better convey the ideas involved, we also describe
the AdaWeight algorithm in detail (even though it is not scalable per se).

| ExpWeight | AcceleWeight | UnixGrad | UPGD | AC-SA | AdaLight
---|---|---|---|---|---|---
Static | ${\log P\big{/}\sqrt{T}}$ | ${\log P\big{/}T^{2}}$ | ${P\big{/}T^{2}}$ | ${\log P\big{/}T^{2}}$ | ${\log P\big{/}T^{2}}$ | ${(\log P)^{3/2}/T^{2}}$
Stochastic | ${\log P\big{/}\sqrt{T}}$ | ✕ | ${P\big{/}\sqrt{T}}$ | ✕ | ${\log P\big{/}\sqrt{T}}$ | ${(\log P)^{3/2}/\sqrt{T}}$
Anytime | partially | ✓ | ✓ | ✕ | ✕ | ✓
Param.-Agn. | ✕ | ✕ | ✓ | ✓ | ✕ | ✓
Compl./Iter. | ${G}$ | ${P}$ | ${P}$ | ${P}$ | ${P}$ | ${G}$
Table 1. Overview of related work in comparison to the AdaLight algorithm
(this paper). For the purposes of this table, $G$ refers to the size of the
underlying graph while $P$ refers to the number of relevant paths in the
network (so $P$ is typically exponential in $G$). The “anytime” property
refers to whether the number of iterations $T$ must be fixed at the outset and
included as a parameter in the algorithm; if not, the algorithm is labeled
“anytime”. Finally, the “parameter-agnostic” property refers to whether any
other parameters – such as the Lipschitz modulus of the network’s cost
functions – need to be known beforehand or not. All estimates are reported in
the $\operatorname{\mathcal{O}}(\cdot)$ sense.

### Paper outline.

Our paper is structured as follows. In Section 2, we formally define our
congestion game setup, the relevant equilibrium notions, and our learning
model. Subsequently, to set the stage for our main results, we present in
Section 3 the classic exponential weights algorithm as well as the
AcceleWeight variant which achieves an $\operatorname{\mathcal{O}}(\log
P/T^{2})$ convergence rate in static environments; both algorithms are non-
adaptive, and they are used as a baseline for our adaptive results. Our
analysis proper begins in Section 4, where we present the AdaWeight algorithm
and its convergence analysis. Then, in Section 5, we present the local flow
setup used to construct the AdaLight algorithm and prove its convergence
guarantees. Finally, in Section 6, we report a series of numerical experiments
validating our theoretical results in real-life transport networks.

##  2\. Problem setup

We begin in this section by introducing the basic elements of our model.

###  2.1. The game.

Building on the classic congestion framework of Beckmann et al. [6], we
consider a class of nonatomic routing games defined by the following three
primitives:  (it n) an underlying _network structure_ ;  (it n) the associated
set of _traffic demands_ ; and  (it n) the network’s _cost functions_.  The
formal definition of each of these primitives is as follows:

  1. 1.

Network structure: Consider a multi-graph
$\mathcal{G}\equiv\mathcal{G}(\mathcal{V},\mathcal{E})$ with vertex set
$\mathcal{V}$ and edge set $\mathcal{E}$. The focal point of interest is a set
of distinct _origin-destination_ (O/D) pairs
$(O^{i},D^{i})\in\mathcal{V}\times\mathcal{V}$ indexed by
$i\in\mathcal{N}=\\{1,\dotsc,N\\}$. For each $i\in\mathcal{N}$, we assume
given a directed acyclic subgraph
$\mathcal{G}^{i}\equiv(\mathcal{V}^{i},\mathcal{E}^{i})$ of $\mathcal{G}$ that
determines the set of routing paths from $O^{i}$ to $D^{i}$, and we write
$\mathcal{P}^{i}$ for the corresponding set of paths joining $O^{i}$ to
$D^{i}$ in $\mathcal{G}^{i}$. For posterity, we will also write
$\mathcal{P}=\operatorname*{\bigcup}_{i\in\mathcal{N}}\mathcal{P}^{i}$ for the
set of all routing paths in the network, and $P=\lvert\mathcal{P}\rvert$ and
$P^{i}=\lvert\mathcal{P}^{i}\rvert$ for the respective cardinalities;
likewise, we will write
$G^{i}=\lvert\mathcal{V}^{i}\rvert+\lvert\mathcal{E}^{i}\rvert$ for the size
of $\mathcal{G}^{i}$, and $G=\sum_{i\in\mathcal{N}}G^{i}$ for the total size
of the network.

  2. 2.

Traffic demands and flows: Each pair $i\in\mathcal{N}$ is associated to a
_traffic demand_ $M^{i}>0$ that is to be routed from $O^{i}$ to $D^{i}$ via
$\mathcal{P}^{i}$; we also write $M_{\mathrm{tot}}=\lVert
M\rVert_{1}=\sum_{i\in\mathcal{N}}M^{i}$ and $M_{\max}=\lVert
M\rVert_{\infty}=\max_{i\in\mathcal{N}}M^{i}$ for the total and maximum
traffic demand associated to the network’s origin-destination (O/D) pairs
respectively.

Now, to route this traffic, the set of feasible _traffic assignment profiles_
– or _flows_ – is defined as

$\textstyle\mathcal{F}\coloneqq\big{\\{}f\in\mathbb{R}_{+}^{\mathcal{P}}:\sum_{p\in\mathcal{P}^{i}}f_{p}=M^{i},i=1,\dotsc,N\big{\\}}$
(1)

i.e., as the product of scaled simplices
$\mathcal{F}=\prod_{i}M^{i}\operatorname{\Delta}(\mathcal{P}^{i})$. In turn,
each feasible flow profile $f\in\mathcal{F}$ induces on each edge
$e\in\mathcal{E}$ the corresponding _traffic load_

$\ell_{e}(f)=\sum_{p\in\mathcal{P}}\operatorname{\mathds{1}}_{\left\\{e\in
p\right\\}}f_{p}$ (2)

i.e., the accumulated mass of all traffic going through $e$. It is also worth
noting here that the path index $p\in\mathcal{P}$ completely characterizes the
O/D pair $i\in\mathcal{N}$ to which it belongs; when we want to make this
relation explicit, we will write $f_{p}^{i}$ instead of $f_{p}$.

  3. 3.

Congestion costs: The traffic routed through a given edge $e\in\mathcal{E}$
incurs a _congestion cost_ depending on the total traffic on the edge and/or
any other exogenous factors. Formally, we will collectively encode all
exogenous factors in a _state variable_ $\omega\in\Omega$ – the “ _state of
the world_ ” – that takes values in some ambient probability space
$(\Omega,\mathcal{F},\operatorname{\mathbb{P}})$. In addition, we will assume
that each edge $e\in\mathcal{E}$ is endowed with a _cost function_
$c_{e}\colon\mathbb{R}_{+}\times\Omega\to\mathbb{R}_{+}$ which determines the
cost $c_{e}(\ell_{e}(f);\omega)$ of traversing $e\in\mathcal{E}$ when the
network is at state $\omega\in\Omega$ and traffic is assigned according to the
flow profile $f\in\mathcal{F}$. Analogously, the cost to traverse a path
$p\in\mathcal{P}$ will be given by the induced _path-cost function_
$c_{p}\colon\mathcal{F}\times\Omega\to\mathbb{R}_{+}$ defined as

$c_{p}(f;\omega)=\sum_{e\in p}c_{e}(\ell_{e}(f);\omega).$ (3)

In this general setting, the only assumption that we will make for the game’s
cost functions is as follows:

######  Assumption 1.

Each cost function $c_{e}(;\omega)$, $e\in\mathcal{E}$, is measurable in
$\omega$ and non-decreasing, bounded and Lipschitz continuous in .
Specifically, there exist $H>0$ and $L>0$ such that $c_{e}(;\omega)\leq H$ and
$\lvert c_{e}(;\omega)-c_{e}(^{\prime};\omega)\rvert\leq
L\lvert^{\prime}-\rvert$ for all $,^{\prime}\in[0,M_{\mathrm{tot}}]$, all
$e\in\mathcal{E}$, and $\operatorname{\mathbb{P}}$-almost all
$\omega\in\Omega$.

Assumption 1 represents a very mild regularity requirement that is satisfied
by most congestion models that occur in practice – including BPR, polynomial,
or regularly-varying latency functions, cf. [35, 8, 26, 12, 38, 36] and
references therein. For this reason, we will treat Assumption 1 as a standing,
blanket assumption and we will not mention it explicitly in the sequel.

###  2.2. Regimes of uncertainty and examples.

The advent of uncertainty in the game’s cost functions (as modeled by
$\omega\in\Omega$) is an important element of our congestion game framework.
To elaborate further on this, it will be convenient to define the _mean cost
function_ of edge $e\in\mathcal{E}$ as

$\displaystyle\bar{c}_{e}(f)$
$\displaystyle=\operatorname{\mathbb{E}}_{\omega}\left[c_{e}(\ell_{e}(f);\omega)\right]$
(4) and consider the corresponding (random) fluctuation process $\displaystyle
U_{e}(f;\omega)$ $\displaystyle=c_{e}(\ell_{e}(f);\omega)-\bar{c}_{e}(f)$ (5)

i.e., the deviation of the network’s cost functions at state $\omega\in\Omega$
from their mean value. The magnitude of these deviations may then be
quantified by the randomness parameter

$\sigma=\operatorname{ess}\sup\nolimits_{\omega\in\Omega}\max\nolimits_{f\in\mathcal{F},e\in\mathcal{E}}\lvert
U_{e}(f;\omega)\rvert$ (6)

where $\operatorname{ess}\sup$ denotes the essential supremum over $\Omega$
with respect to $\operatorname{\mathbb{P}}$ (so $\sigma\leq 2H$ by Assumption
1). Informally, larger values of $\sigma$ indicate a higher degree of
randomness, implying in turn that the traffic assignment problem becomes more
difficult to solve; on the other hand, if $\sigma=0$, one would expect
algorithmic methods to achieve better results. This distinction plays a key
role in the sequel so we formalize it as follows:

######  Definition 1 (Static and Stochastic Environments).

When $\sigma=0$, we say that the environment is _static_ ; otherwise, if
$\sigma>0$, we say that the environment is _stochastic_.

The dependence of the network’s cost functions on exogenous random factors is
the main difference of our setup with standard congestion models in the spirit
of Beckmann et al. [6] and Nisan et al. [35]. For concreteness, we mention
below two existing models that can be seen as special cases of our framework:

######  Example 2.1 (Deterministic routing games).

The _static regime_ described above matches the deterministic routing models
of Fischer and Vöcking [16], Blum et al. [8] and Krichene et al. [23, 24]. In
these models, the network only has a single (deterministic) state so,
trivially, $\sigma=0$. §

######  Example 2.2 (Noisy cost measurements).

To model uncertainty in the cost measurement process, Krichene et al. [24, 21]
considered the case where the network’s cost functions are fixed, but cost
measurements are only accurate up to a random, zero-mean error. In our
framework, this can be modeled by simply assuming a familly of cost functions
of the form $c_{e}(f;\omega)=\bar{c}_{e}(f)+\omega_{e}$ with
$\operatorname{\mathbb{E}}[\omega_{e}]=0$ for all $e\in\mathcal{E}$. §

Besides these standard examples, our model is sufficiently flexible to capture
other random factors such as weather conditions, traffic accidents, road
incidents, etc. For instance, to model the difference between dry and wet
weather, one can take
$\Omega=\\{\omega_{\textrm{dry}},\omega_{\textrm{wet}}\\}$ and consider a set
of cost functions with
$c_{e}(\ell_{e}(f);\omega_{\textrm{dry}})<c_{e}(\ell_{e}(f);\omega_{\textrm{wet}})$
to reflect the fact that the congestion costs are higher when the roads are
wet. However, to maintain the generality of our model, we will not focus on
any particular application.

###  2.3. Notions of equilibrium.

In the above framework, each state variable $\omega\in\Omega$ determines an
instance of a routing game, defined formally as a tuple
$\Gamma_{\omega}\equiv\Gamma_{\omega}(\mathcal{G},\mathcal{N},\mathcal{P},c_{\omega})$
where $c_{\omega}$ is shorthand for the network’s cost functions
$\\{c_{e}(\cdot;\omega)\\}_{e\in\mathcal{E}}$ instantiated at $\omega$. Of
course, in analyzing the game, each individual instance $\Gamma_{\omega}$ is
meaningless by itself unless $\operatorname{\mathbb{P}}$ assigns positive
probability only to a _single_ $\omega$. For this reason, we will instead
focus on the _mean game_
$\Gamma\equiv\Gamma(\mathcal{G},\mathcal{N},\mathcal{P},\bar{c})$ which has
the same network and routing structure as every $\Gamma_{\omega}$,
$\omega\in\Omega$, but with congestion costs given by the _mean cost
functions_

$\bar{c}_{p}(f)=\operatorname{\mathbb{E}}[c_{p}(f;\omega)]=\sum\nolimits_{e\in
p}\operatorname{\mathbb{E}}[c_{e}(\ell_{e}(f);\omega)].$ (7)

Motivated by the route recommendation problem described in the introduction,
we will focus on traffic assignment profiles where Wardrop’s unilateral
optimality principle [44] holds on average, i.e., _all traffic is routed along
a path with minimal mean cost._ Formally, we have:

######  Definition 2 (Mean equilibrium flows).

We say that $f^{\ast}\in\mathcal{F}$ is a _mean equilibrium flow_ if and only
if $\bar{c}_{p}(f^{\ast})\leq\bar{c}_{q}(f^{\ast})$ for all $i\in\mathcal{N}$
and all $p,q\in\mathcal{P}^{i}$ such that $(f^{\ast})_{p}^{i}>0$.

######  Remark 1.

Definition 2 means that, on average, no user has an incentive to deviate from
the recommended route; obviously, when the support of
$\operatorname{\mathbb{P}}$ is a singleton, we recover the usual definition of
a _Wardrop equilibrium_ [44, 6, 35]. This special case will be particularly
important and we discuss it in detail in the next section.

Importantly, the problem of finding an equilibrium flow of a (fixed) routing
game $\Gamma_{\omega}$ admits a _potential function_ – often referred to as
the Beckmann–McGuire–Winsten (BMW) potential [6, 14]. Specifically, for a
given instance $\omega\in\Omega$, the BMW potential is defined as

$\Phi_{\omega}(f)=\sum\nolimits_{e\in\mathcal{E}}\int_{0}^{\ell_{e}(f)}c_{e}(;\omega)\>d\quad\text{for
all $f\in\mathcal{F}$},$ (BMW)

and it has the fundamental property that
$\operatorname*{arg\,min}\Phi_{\omega}$ coincides with the set
$\operatorname{Eq}(\Gamma_{\omega})$ of equilibria of $\Gamma_{\omega}$. In
our stochastic context, a natural question that arises is whether this
property can be extended to the mean game $\Gamma$ where the randomness is
“averaged over”. Clearly, the most direct candidate for a potential function
in this case is the _mean BMW potential_

$\Phi(f)=\operatorname{\mathbb{E}}[\Phi_{\omega}(f)]=\operatorname{\mathbb{E}}\left[\sum_{e\in\mathcal{E}}\int_{0}^{\ell_{e}(f)}c_{e}(;\omega)\>d\right].$
(8)

Since $c_{e}$ is continuous and non-decreasing for all $e\in\mathcal{E}$, the
mean BMW potential $\Phi$ is, in turn, continuously differentiable and convex
on $\mathcal{F}$. Of course, since the probability law
$\operatorname{\mathbb{P}}$ is not known, we will _not_ assume that $\Phi$
(and/or its gradients) can be explicitly computed in general. Nevertheless,
building on the deterministic characterization of Beckmann et al. [6], we have
the following equivalence between mean equilibria and minimizers of $\Phi$:

######  Proposition 1.

$\Phi$ is a potential function for the mean game
$\Gamma\equiv\Gamma(\mathcal{G},\mathcal{N},\mathcal{P},\bar{c})$; in
particular, we have

$\frac{\partial\Phi}{\partial f_{p}}=\bar{c}_{p}(f)\quad\text{for all
$f\in\mathcal{F}$ and all $p\in\mathcal{P}$}.$ (9)

Accordingly, a flow profile $f^{\ast}\in\mathcal{F}$ is an equilibrium of the
mean game $\Gamma$ if and only if it is a minimizer of $\Phi$ over
$\mathcal{F}$; more succinctly, the equilibrium set
$\operatorname{Eq}(\Gamma)$ of $\Gamma$ satisfies
$\operatorname{Eq}(\Gamma)=\operatorname*{arg\,min}_{f\in\mathcal{F}}\Phi(f)$.

###### Proof.

First, apply the dominated convergence theorem (saying that gradients and
expectations commute), we have $\frac{\partial\Phi}{\partial
f_{p}}=\operatorname{\mathbb{E}}\left[\frac{\partial\Phi_{\omega}}{\partial
f_{p}}\right]=\operatorname{\mathbb{E}}\left[c_{p}(f;\omega)\right]=\bar{c}_{p}(f)$.
The rest of the proof of Proposition 1 can be obtained by following [6, 40]
for the mean game $\Gamma$. ∎

In view of Proposition 1, $\Phi$ provides a natural merit function for
examining how close a given flow profile $f\in\mathcal{F}$ is to being an
equilibrium. Formally, given a sequence of candidate flow profiles
$f^{t}\in\mathcal{F}$, $t=1,2,\dotsc\,$, we define the associated _equilibrium
gap_ after $T$ epochs as

$\operatorname{Gap}(T)=\Phi(f^{T})-\min\Phi.$ (10)

Clearly, the sequence in question converges to $\operatorname{Eq}(\Gamma)$ if
and only if $\operatorname{Gap}(T)\to 0$. Moreover, the speed of this
convergence is also captured by $\operatorname{Gap}(T)$, so all our
convergence certificates will be stated in terms of $\operatorname{Gap}(T)$.

###  2.4. Sequence of events.

With all this in hand, the sequence of events in our model unfolds as follows:

  1. i.5pt)

At each time slot $t=1,2,\dotsc\,$, the navigation interface selects a traffic
assignment profile $f^{t}\in\mathcal{F}$ and it routes all demands received
within this time slot according to $f^{t}$.

  2. ii.5pt)

Concurrently, the state $\omega^{t}$ of the network is drawn i.i.d. from
$\Omega$ according to $\operatorname{\mathbb{P}}$.

  3. iii.5pt)

The interface observes the realized congestion costs
$C^{t}_{e}=c_{e}(\ell_{e}(f^{t});\omega^{t})$ along each edge
$e\in\mathcal{E}$; subsequently, it uses this information to update the
routing profile $f^{t}$, and the process repeats.

In what follows, we will analyze different learning algorithms for updating
$f^{t}$ in the above sequence of events, and we will examine their equilibrium
convergence properties in terms of the equilibrium gap function (10). For
posterity, we emphasize again that the optimizer has no prior knowledge of the
enviroment, the distribution of random events, the network’s cost functions,
etc., so we will pay particular attention to the dependence of each
algorithm’s guarantees on these (otherwise unknown) parameters.

##  3\. Non-adaptive methods

In the rest of our paper, we discuss a series of routing algorithms and their
equilibrium convergence guarantees. To set the stage for our main
contributions, we begin by presenting two learning methods that are order-
optimal in the two basic regimes described in the previous section:

  1. i.5pt)

The “vanilla” exponential weights (ExpWeight) method of Blum et al. [8] for
_stochastic_ environments.

  2. ii.5pt)

An accelerated exponential weights algorithm – which we call AcceleWeight –
for _static_ environments.

Both methods rely crucially on the (rescaled) _logit choice map_
$\operatorname{\Lambda}\colon\mathbb{R}^{\mathcal{P}}\to\mathcal{F}$, given in
components as

$\operatorname{\Lambda}_{p}(w)=\frac{M^{i}\exp(w_{p})}{\sum_{q\in\mathcal{P}^{i}}\exp(w_{q})}\qquad\text{for
all $w\in\mathbb{R}^{\mathcal{P}}$, $p\in\mathcal{P}^{i}$ and
$i\in\mathcal{N}$}.$ (11)

For concision, we will also write $C^{t}_{p}=\sum_{e\in p}C^{t}_{e}$ for the
cost of path $p\in\mathcal{P}$ at stage $t$, and
$C^{t}=(C^{t}_{p})_{p\in\mathcal{P}}$ for the profile thereof.

###  3.1. EW in stochastic environments.

We begin by presenting the exponential weights algorithm of Blum et al. [8] in
standard pseudocode form:

1

2Initialize $w^{1}\leftarrow 0$

3 for  _$t=1,2,\dotsc\,$_ do

set $f^{t}\leftarrow\operatorname{\Lambda}(w^{t})$ and get $C^{t}\leftarrow
c(f^{t},\omega^{t})$

// route and measure costs

set $w^{t+1}=w^{t}-\gamma^{t}C^{t}$

// update path scores

4

5

Algorithm 1 EW

We then have the following equilibrium convergence result:

######  Theorem 1.

Suppose that ExpWeight (Algorithm 1) is run for $T$ epochs with a fixed
learning rate
$\gamma^{t}=\sqrt{\log\left({M_{\max}P}\big{/}{M_{\mathrm{tot}}}\right)}/(H\sqrt{T})$
for any $t$. Then the time-averaged flow profile
$\bar{f}^{T}=(1/T)\sum_{t=1}^{T}f^{t}$ enjoys the equilibrium convergence
rate:

$\operatorname{\mathbb{E}}\left[\operatorname{Gap}\left(T\right)\right]\leq\frac{2M_{\mathrm{tot}}H\sqrt{\log\left({M_{\max}P}\big{/}{M_{\mathrm{tot}}}\right)}}{\sqrt{T}}=\operatorname{\mathcal{O}}\left(\frac{(\log
P)^{1/2}}{\sqrt{T}}\right).$ (12)

The proof of Theorem 1 relies on standard techniques, so we omit it; for a
series of closely related results, we refer the reader to Krichene et al. [21]
and Shalev-Shwartz [41, Corollary 2.14 and Theorem 4.1]; What is more
important for our purposes is that Theorem 1 confirms that ExpWeight achieves
an $\operatorname{\mathcal{O}}\big{(}1\big{/}\sqrt{T}\big{)}$ convergence
speed in the number of epochs $T$, so it is order-optimal in stochastic
environments [31, 9]. Moreover, the algorithm’s convergence speed is
logarithmic in $P$, and hence (at most) linear in the size $G$ of the
underlying network.

On the downside, the convergence rate obtained for Algorithm 1 concerns the
time-averaged flow $\bar{f}^{T}$, not the actual recommendation, a distinction
which is important for practical applications. Second, Theorem 1 is not an
“anytime” convergence result in the sense that the algorithm’s learning rate
$\gamma^{t}$ must be tuned relative to $T$ (the rate is not valid otherwise).
This issue can be partially resolved either via a doubling trick [41] (at the
cost of restarting the algorithm ever so often) or by using a variable
learning rate of the form $\gamma^{t}\propto 1/\sqrt{t}$ (at the cost of
introducing an additional $\log T$ factor in the algorithm’s convergence
speed). However, in either case, the rate (12) only applies to the sequence of
time averages, not the actual recomendations.

###  3.2. XLEW in static environments.

We now turn our attention to the static regime, i.e., when there are no
exogenous variations in the network’s cost functions ($\sigma=0$); in this
case, it is reasonable to expect that the
$\operatorname{\mathcal{O}}(1/\sqrt{T})$ convergence speed of ExpWeight can be
improved. Indeed, as we show below, the BMW potential of the game is
_Lipschitz smooth_ , so the optimal convergence speed in static environments
is $\operatorname{\mathcal{O}}(1/T^{2})$ [31, 9].

######  Proposition 2.

The BMW potential $\Phi$ is Lipschitz smooth relative to the $L^{1}$ norm on
$\mathcal{F}$. In particular, its smoothness modulus is $\beta=KL$, where
$K=\max_{p\in\mathcal{P}}\sum_{e\in p}\operatorname{\mathds{1}}_{\left\\{e\in
p\right\\}}$ is the length of the longest path in $\mathcal{P}$.

###### Proof of Proposition 2.

It suffices to show that the gradient $\nabla\Phi$ of $\Phi$ is Lipschitz
continuous with respect to the (primal) $L^{1}$ norm on $\mathcal{F}$ and the
(dual) $L^{\infty}$ norm on $\mathbb{R}^{\mathcal{P}}$. Indeed, for all
$f,f^{\prime}\in\mathcal{F}$, combining Propositions 1 and 1 yields

$\displaystyle\lVert\nabla\Phi(f)-\nabla\Phi(f^{\prime})\rVert_{\infty}$
$\displaystyle=\max_{p\in\mathcal{P}}\lvert\bar{c}_{p}(f)-\bar{c}_{p}(f^{\prime})\rvert$
[by Proposition 1] $\displaystyle\leq\max_{p\in\mathcal{P}}\sum_{e\in
p}L\left\lvert\ell_{e}(f)-\ell_{e}(f^{\prime})\right\rvert$ [by Assumption 1]
$\displaystyle\leq
L\cdot\max_{p\in\mathcal{P}}\sum_{e\in\mathcal{E}}\operatorname{\mathds{1}}_{\left\\{e\in
p\right\\}}\sum_{q\in\mathcal{P}}\operatorname{\mathds{1}}_{\left\\{e\in
q\right\\}}\lvert f_{q}-f^{\prime}_{q}\rvert$ [definition of $\ell_{e}$]
$\displaystyle\leq
L\cdot\sum_{e\in\mathcal{E}}\sum_{q\in\mathcal{P}}\operatorname{\mathds{1}}_{\left\\{e\in
q\right\\}}\lvert f_{q}-f^{\prime}_{q}\rvert$ [since
$\operatorname{\mathds{1}}_{\left\\{e\in p\right\\}}\leq 1$]
$\displaystyle\leq
L\cdot\max_{p\in\mathcal{P}}\sum_{e\in\mathcal{E}}\operatorname{\mathds{1}}_{\left\\{e\in
p\right\\}}\cdot\sum_{q\in\mathcal{P}}\lvert
f_{q}-f^{\prime}_{q}\rvert=KL\lVert f-f^{\prime}\rVert_{1}$ [by Cauchy-
Schwarz]

so our claim follows. ∎

Given that $\Phi$ is smooth, the iconic $\operatorname{\mathcal{O}}(1/T^{2})$
rate mentioned above can be achieved by the accelerated gradient algorithm of
Nesterov [30]. However, if Nesterov’s algorithm were applied “off the shelf”,
the constants involved would be linear in $P$ (the problem’s dimensionality),
and hence exponential in the size of the network. Instead, building on ideas
by Allen-Zhu and Orecchia [1] and Krichene et al. [22], this problem can be
overcome by combining the acceleration mechanism of Nesterov [30] with the
dimensionally-efficient ExpWeight template. The resulting _accelerated
exponential weights_ (AcceleWeight) method proceeds in pseudocode form as
follows:

1

Input: Smoothness parameter $\beta=KL$

2

3Initialize $w^{1}\leftarrow 0$; $f^{0}\leftarrow 0$; $\alpha^{0}\leftarrow
0$; $\gamma^{0}\leftarrow(NM_{\max}\beta)^{-1}$

4

5for  _$t=1,2,\dotsc\,$_ do

set $z^{t}\leftarrow\operatorname{\Lambda}(w^{t})$

// exploratory flow obtained from path scores

set $f^{t}\leftarrow\alpha^{t-1}f^{t-1}+(1-\alpha^{t-1})z^{t}$

// average with previous state

set
$\gamma^{t}\leftarrow\gamma^{t-1}+\gamma^{0}/2+\sqrt{\gamma^{t-1}\gamma^{0}+(\gamma^{0}/2)^{2}}$

// update step-size

set $\alpha^{t}\leftarrow{\gamma^{t-1}}\big{/}{\gamma^{t}}$

// update averaging weight

set $\bar{f}^{t}\leftarrow\alpha^{t}f^{t}+(1-\alpha^{t})z^{t}$ and get
$\bar{C}^{t}\leftarrow c\left(\bar{f}^{t},\omega^{t}\right)$

// route and measure costs

set $w^{t+1}\leftarrow w^{t}-(1-\alpha^{t})\gamma^{t}\bar{C}^{t}$

// update path scores

6

7

Algorithm 2 XLEW

At a high level, AcceleWeight follows the acceleration template of Nesterov
[30], but replaces Euclidean projections with the logit map (11). With this in
mind, we obtain the following convergence guarantee:

######  Theorem 2.

The sequence of flow profiles $f^{1},f^{2},\dots$ generated by AcceleWeight
(Algorithm 2) in a static environment ($\sigma=0$) enjoys the equilibrium
convergence rate

$\operatorname{Gap}\left({T}\right)\leq\frac{4\beta
N^{2}M_{\max}^{2}\log\left(M_{\max}P/M_{\mathrm{tot}}\right)}{(T-1)^{2}}=\operatorname{\mathcal{O}}\left(\frac{\log
P}{T^{2}}\right).$ (13)

The proof of Theorem 2 is based on techniques that are widely used in the
analysis of accelerated methods so, to streamline our presentation, we
relegate it to Appendix A.

Theorem 2 confirms that, in the static regime, AcceleWeight converges to
equilibrium with an optimal
$\operatorname{\mathcal{O}}\big{(}1\big{/}T^{2}\big{)}$ convergence speed, as
desired. For our purposes, it is more important to note that Theorem 2
confirms that, in static environments, accelerated exponential weights
achieves a much faster convergence rate than ExpWeight. However, this
improvement comes with several limitations: First and foremost, as we discuss
in Section 6, AcceleWeight fails to converge altogether in the stochastic
regime, so it is not a universal method (i.e., an algorithm that is
simultaneously order-optimal in both static and stochastic environments).
Second, to obtain the accelerated rate (13), the step-size of AcceleWeight
must be tuned with prior knowledge of the problem’s smoothness modulus $\beta$
(which is not realistically available to the learner), so the method is not
adaptive. These are both crucial limitations that we seek to overcome in
Sections 4 and 5.

###  3.3. Per-iteration complexity of each method.

We close this section with a short discussion on the per-iteration complexity
of ExpWeight and AcceleWeight. The key point of note is that Algorithms 1 and
2 both require $\operatorname{\mathcal{O}}(P)$ time and space per iteration
(i.e., in terms of arithmetic operations and memory storage). This is
inefficient in large-scale networks where $P$ is exponentially large in $G$
(the size of the network), so it is not clear how to implement either of these
algorithms in practice.

As discussed in Section 1, there exist variants of the exponential weights
algorithm that _can_ be implemented efficiently in
$\operatorname{\mathcal{O}}(G)$ time and space per iteration. This is
accomplished via a dynamic programming procedure known as “weight-pushing”
[42, 18]; however, this technique is limited to boosting the computation step
involving the logit mapping (Line 3 of Algorithm 1 and Algorithm 2 of
Algorithm 2) and cannot efficiently execute other steps of AcceleWeight that
require $\operatorname{\mathcal{O}}(P)$ operations or storage bits (like the
averaging steps involved in the acceleration mechanism). To maintain the flow
of our discussion, we will revisit this issue in detail in Section 5.

##  4\. AEW: Adaptive learning in the presence of uncertainty

###  4.1. Adaptivity of AdaWeight in optimizing equilibrium-convergence rates

To summarize the situation so far, we have seen that ExpWeight attains an
$\operatorname{\mathcal{O}}(\log(P)/\sqrt{T})$ rate, which is order-optimal in
the stochastic case but suboptimal in static environments; by contrast,
AcceleWeight attains an $\operatorname{\mathcal{O}}(\log(P)/T^{2})$ rate in
static environment, but has no convergence guarantees in the presence of
randomness and uncertainty. Consequently, neither of these algorithms meets
our stated objective to concurrently achieve order-optimal guarantees in both
the static and stochastic cases (and without requiring prior knowledge of the
problem’s smoothness modulus).

To resolve this gap, we propose below an _adaptive exponential weights_ method
– AdaWeight for short – which achieves these objectives by mixing the
acceleration template of AcceleWeight with the dual extrapolation method of
Nesterov [32]. We present the pseudocode of AdaWeight below:

1 Initialize $\alpha^{0}\leftarrow 0$, $z^{0}\leftarrow 0$,
$\eta^{1}\leftarrow 1$ and $w^{1}\leftarrow 0$

2 for  _$t=1,2,\ldots$_ do

set $r^{t}\leftarrow\sum_{s=0}^{t-1}\alpha^{s}z^{s}$

// set an anchor primal point

// the test phase:

3 set
$\tilde{z}^{t}\leftarrow\operatorname{\Lambda}\left(\eta^{t}w^{t}\right)$

set
$\tilde{f}^{t}\leftarrow\left(\alpha^{t}\tilde{z}^{t}+r^{t}\right)\big{/}\sum_{s=0}^{t}\alpha^{s}$
and get $\tilde{C}^{t}\leftarrow c(\tilde{f}^{t},\omega^{t})$

// reweigh with anchor + query a test point

set $\tilde{w}^{t}\leftarrow w^{t}-\alpha^{t}\tilde{C}^{t}$

// exploratory score update

// the recommendation phase:

4 set
$z^{t}\leftarrow\operatorname{\Lambda}\left(\eta^{t}\tilde{w}^{t}\right)$

set
$f^{t}\leftarrow\left(\alpha^{t}z^{t}+r^{t}\right)\big{/}\sum_{s=0}^{t}\alpha^{s}$
and get $C^{t}\leftarrow c(f^{t},\omega^{t})$

// reweigh with anchor + route and measure costs

set $w^{t+1}\leftarrow w^{t}-\alpha^{t}C^{t}$

// update path scores

set $\eta^{t+1}\leftarrow
1\big{/}\sqrt{1+\sum_{s=0}^{t}\left\lVert\alpha^{s}(C^{s}-\tilde{C}^{s})\right\rVert_{\infty}^{2}}$

// update learning rate

5

Algorithm 3 AEW

The main novelty in the definition of the AdaWeight algorithm is the
introduction of two “extrapolation” sequences, $z^{t}$ and $\tilde{w}^{t}$,
that venture outside the convex hull of the generated primal (flow) and dual
(score) variables respectively. These leading states are subsequently
averaged, and the method proceeds with an adaptive step-size rule. In more
details, AdaWeight relies on three key components:

  1. i.5pt)

A dual extrapolation mechanism for generating the leading sequences in Lines 3
and 3; these sequences are central for anticipating the loss landscape of the
problem.

  2. ii.5pt)

An acceleration mechanism obtained from the ($\alpha^{t}$)-_weighted average_
steps in Lines 3 and 3; in the analysis, $\alpha^{t}$ will grow as $t$, so
almost all the weight will be attributed to the state closest to the current
one.

  3. iii.5pt)

An _adaptive sequence of learning rates_ (cf. Line 3) in the spirit of Rakhlin
and Sridharan [39], Kavis et al. [20] and Antonakopoulos et al. [2]. This
choice is based on the ansatz that, if the algorithm encounters coherent
gradient updates (which can only occur in static environments), it will
eventually stabilize to a strictly positive value; otherwise, it will decrease
to zero at a $\Theta(1/\sqrt{t})$ rate. This property is crucial to
interpolate between the stochastic and static regimes.

The combination of the weighted average iterates and adaptive learning rate in
AdaWeight is shared by the UnixGrad algorithm proposed by Kavis et al. [20],
which also provides rate interpolation in constrained problems. However,
UnixGrad requires the problem’s domain to have a finite Bregman diameter –
and, albeit compact, the set of feasible flows $\mathcal{F}$ has an _infinite_
diameter under the entropic regularizer that generates the ExpWeight template.
Therefore, UnixGrad is not applicable to our routing games. This is the reason
for switching gears to the “primal-dual” approach offered by the dual
extrapolation template; this primal-dual interplay provides the missing link
that allows AdaWeight to simultaneously enjoy order-optimal convergence
guarantees in both settings while maintaining the desired polynomial
dependency on the problem’s dimension.

In light of the above, our main convergence result for AdaWeight is as
follows:

######  Theorem 3.

Let $f^{1},f^{2},\ldots,$ be the sequence of flows recommended by AdaWeight
(Algorithm 3) running with $\alpha^{t}=t$ for any $t=1,2,\dotsc\,$, then it
enjoys the following equilibrium convergence rate:

$\displaystyle\operatorname{\mathbb{E}}\left[\operatorname{Gap}\left(T\right)\right]\leq\frac{2\sqrt{2}{A}\sigma}{\sqrt{T}}+\frac{16\beta\sqrt{NM_{\max}}{A}^{{3}/{2}}+{B}}{T^{2}}=\operatorname{\mathcal{O}}\left(\left(\log
P\right)^{{3}/{2}}\left(\frac{\sigma}{\sqrt{T}}+\frac{1}{T^{2}}\right)\right).$
(14)

Specifically, in the static case, AdaWeight enjoys the sharper rate:
$\operatorname{Gap}\left(T\right)\leq\operatorname{\mathcal{O}}\left({\left(\log
P\right)^{{3}/{2}}}\big{/}{T^{2}}\right)$. In the above expressions, $\sigma$
is as defined in Equation 6 while ${A}$ and ${B}$ are positive constants given
by ${A}\coloneqq
NM_{\max}\left[2\log\left({PM_{\max}}/{M_{\mathrm{tot}}}\right)+13\right]=\operatorname{\mathcal{O}}(\log{P})$
and ${B}\coloneqq
M_{\mathrm{tot}}\log(PM_{\max}/M_{\mathrm{tot}})=\operatorname{\mathcal{O}}\left({\log
P}\right)$.

Theorem 3 characterizes the convergence speed of AdaWeight according to the
number of learning iterations $T$, the number of paths $P$ and the level of
randomness $\sigma$ of the network’s states (as defined in (6)). In the
stochastic environment, $\sigma\geq 0$ and hence, AdaWeight enjoys a
convergence rate of order $\operatorname{\mathcal{O}}\left((\log
P)^{3/2}\sigma\big{/}\sqrt{T}\right)$. In the static case, $\sigma=0$ and the
convergence rate is accelerated to $\operatorname{\mathcal{O}}\left((\log
P)^{3/2}\big{/}T^{2}\right)$. In instances where the cost observations become
more accurate over time, Theorem 3 also allows us to achieve a smooth trade-
off in the convergence rate. For example, if
$\sigma=\operatorname{\mathcal{O}}\left(1/T\right)$ for some suitable ,
AdaWeight ’s rate of convergence carries a dependence of the order of
$\operatorname{\mathcal{O}}\left(T^{\max\\{--1/2,-2\\}}\right)$.

In view of these result, Theorem 3 confirms that AdaWeight satisfies the
following desiderata:  (itn) it achieves _simultaneously optimal guarantees_
in both stochastic and static environments in the number of learning
iterations (i.e., $\operatorname{\mathcal{O}}(1/\sqrt{T})$ and
$\operatorname{\mathcal{O}}(1/T^{2})$ respectively);  (itn) the derived rates
maintain a _polynomial dependency_ in terms of the network’ s combinatorial
primitives; and  (itn) it requires _no prior tuning_ by the learner.
Moreover, unlike ExpWeight, the convergence of AdaWeight corresponds to an
actual traffic flow profile that is implemented in epoch $t$ and not the
average flow.

##### The per-iteration complexities of AdaWeight .

The AdaWeight method, as presented in Algorithm 3, while having a simple
presentation, takes $\mathcal{O}(P)$ time (and space) to complete each
iteration. Therefore, implementing AdaWeight is inefficient. Although
AdaWeight shares certain elements of the ExpWeight template, the naive
application of weight-pushing _fails_ to obtain an implementation of AdaWeight
having a polynomial per-iteration complexity in the network’s size. This is
due to the complicated averaging steps in AdaWeight (Lines 3 and 3 of
Algorithm 3). In Section 5, we re-discuss this in more details and we propose
another algorithm, called AdaLight, that maintains all desired features of
AdaWeight and achieves an efficient per-iteration complexity (i.e., sub-
quadratic in the size of the underlying graph).

###  4.2. Proof of Theorem 3

AdaWeight is _not_ merely a “convex combination” of the two non-adaptive
optimal algorithms (AcceleWeight and ExpWeight). For this reason, the proof of
Theorem 3 is technically involved. Particularly, the primal-dual averaging
method that we use in AdaWeight create sequences of filtration-dependent step-
sizes; while this is the key element allowing AdaWeight to achieve the best-
of-both-worlds, previously known approaches in deriving the convergence
analyses (as in non-adaptive algorithmic schemes – ExpWeight and AcceleWeight
– and/or in UnixGrad) are not applicable in this case. To handle this new
challenge, we propose a completely new way to treat the learning rate in order
to make the derived bounds summable. To ensure the comprehensibility, in this
section, we first present a sketch of proof explaining the high-level idea
before presenting the complete proof with all technical details.

####  4.2.1. Sketch of proof (Theorem 3)

Let
$\operatorname{\overline{\operatorname{\mathcal{R}}}}_{T}(f^{\ast})\coloneqq\sum_{t=1}^{T}{\alpha^{t}}\langle\nabla\Phi(f^{t}),z^{t}-f^{\ast}\rangle$,
the starting point of our proof is the following result:

$\operatorname{\mathbb{E}}\left[\operatorname{Gap}({T})\right]=\operatorname{\mathbb{E}}\left[\Phi(f^{T})-\Phi(f^{\ast})\right]\leq{2\operatorname{\mathbb{E}}\left[\operatorname{\overline{\operatorname{\mathcal{R}}}}_{T}\left(f^{\ast}\right)\right]}\big{/}{T^{2}}\textrm{
for any }T.$ (15)

(15) is a standard result in working with the $\alpha^{t}$-weighted averaging
technique that also appears in several previous works [13, 20]. For the sake
of completeness, we provide the proof of (15) in B.

Recall that as AdaWeight is run, $\operatorname{\nabla}\Phi(f^{t})$ is _not_
observable and only the cost $C^{t}$ is observed and used. Therefore, instead
of $\operatorname{\overline{\operatorname{\mathcal{R}}}}_{T}(f^{\ast})$, we
now _focus on the term_
$\operatorname{\mathcal{R}}_{T}(f^{\ast})\coloneqq\sum_{t=1}^{T}{\alpha^{t}}\left\langle
C^{t},z^{t}-f^{\ast}\right\rangle$. In the stochastic case, we can prove that
$\operatorname{\mathcal{R}}_{T}(f^{\ast})$ is an unbiased estimation of
$\operatorname{\overline{\operatorname{\mathcal{R}}}}_{T}(f^{\ast})$ (and in
the static case, they coincide). Therefore, any upper-bound of
$\operatorname{\mathcal{R}}_{T}(f^{\ast})$ can be translated via (15) into an
upper-bound of the left-hand-sides of (14). The key question becomes _“Which
upper-bound of $\operatorname{\mathcal{R}}_{T}(f^{\ast})$ can be derived to
guarantee the convergence rates in (14)?”_

To answer the question above, we note that AdaWeight is built on the dual
extrapolation template with two phases: the test phase and the recommendation
phase. This allows us to upper-bound
$\operatorname{\mathcal{R}}_{T}(f^{\ast})$ in terms of the “distance” between
the pivot primal points in these phases (i.e., $\tilde{z}^{t}$ and $z^{t}$)
and the difference between the costs measured in these phases (i.e.,
$\tilde{C}^{t}$ and $C^{t}$). Particularly, we have

$\operatorname{\mathcal{R}}_{T}(f^{\ast})\leq\sum_{t=1}^{T}\textsl{g}_{\textrm{primal}}(\eta^{t+1}){\lVert
z^{t}-\tilde{z}^{t}\rVert_{1}^{2}}+\sum_{t=1}^{T}\textsl{g}_{\textrm{dual}}(\eta^{t+1}){(\alpha^{t})^{2}\left\lVert
C^{t}-\tilde{C}^{t}\right\rVert_{\infty}^{2}},$ (16)

where $\textsl{g}_{\textrm{primal}}(\eta^{t+1})$ and
$\textsl{g}_{\textrm{dual}}(\eta^{t+1})$ are certain functions that also
depend on other parameters of the game. Importantly, the terms in the right-
hand-size of (16) are actually summable: they are bounded by the summation of
two terms $\sum_{t=1}^{T}\textsl{g}(\eta^{t})(\alpha^{t})^{2}\lVert
C^{t}-\tilde{C}^{t}\rVert_{\infty}^{2}$ and
$\sqrt{\sum_{t=1}^{T}(\alpha^{t})^{2}\lVert[C^{t}-\tilde{C}^{t}]-[\operatorname{\nabla}\Phi(f^{t})-\operatorname{\nabla}\Phi(\tilde{f}^{t})]\rVert_{\infty}^{2}}$.
Here, the first term arises when we bound $\lVert z^{t}-\tilde{z}^{t}\rVert$
using the smoothness of the BMW potential $\Phi$. Moreover, by our _special
choice of adaptive sequence of learning rate_ $\eta^{t}$, there exists
$T_{0}\ll T$ such that only the first $T_{0}$ components in this summation are
positive and hence, this summation is actually of order
$\operatorname{\mathcal{O}}(1)$ (i.e., it does not depend on $T$). The second
term quantifies the error of the observed costs with respect to the actual
gradients of the potential. This term is of order
$\operatorname{\mathcal{O}}(\sigma T^{3/2})$ when we choose $\alpha^{t}=t$.
Combine these results with (15), we obtain the convergence rates indicated in
(14).

####  4.2.2. Proof with technical details (Theorem 3)

In this section, we work with the _entropy regularizer_
$h(f)\coloneqq\sum_{i\in\mathcal{N}}\sum_{p\in\mathcal{P}^{i}}f^{i}_{p}\log(f^{i}_{p}),\forall
f\in\mathcal{F}$ which can be used to define the exponential weights template.
Trivially, $h$ is $1/\kappa$-_strongly convex_ on $\mathcal{F}$ w.r.t the
$\lVert\cdot\rVert_{1}$ norm, where we define $\kappa\coloneqq M_{\max}N$. We
will also work with the _Fenchel conjugate_ of $h$, defined as
$h^{\ast}(w)\coloneqq\max_{f\in\mathcal{F}}\left\langle
f,w\right\rangle-h(f)$, and its _Fenchel coupling_ , defined as
$F(f,w)=h\left(f\right)+{h^{*}}\left(w\right)-\left\langle f,w\right\rangle$
for any $f\in\mathcal{F}$ and any $w\in\mathbb{R}^{P}$. For the sake of
brevity, we will also define $\min
h\coloneqq\min_{f^{\prime}\in\mathcal{F}}h(f)$ and denote the
$\lVert\cdot\rVert_{1}$-diameter of $\mathcal{F}$ by
$D\coloneqq\max_{f,f^{\prime}\in\mathcal{F}}\lVert f-f^{\prime}\rVert_{1}$. We
also denote the Kullback–Leibler divergence between two flows $f$ and
$f^{\prime}$ by
$\mathrm{KL}(f\nonscript\,\|\nonscript\,\mathopen{}f^{\prime})$. Note also
that throughout this section, when we mention Algorithm 3, we understand that
it is run with $\alpha^{t}\coloneqq t$ as chosen in Theorem 3.

First, we prove the following proposition (corresponding to Equation (16) in
the proof sketch):

######  Proposition 3.

Run Algorithm 3, for any $f\in\mathcal{F}$, we have:

$\operatorname{\mathcal{R}}_{T}(f^{\ast})\leq h(f^{\ast})-\min
h+{A}_{\texttt{dual}}\sum_{t=1}^{T}(\alpha^{t})^{2}\eta^{t+1}\left\lVert
C^{t}-\tilde{C}^{t}\right\rVert_{\infty}^{2}-\sum_{t=1}^{T}\frac{1}{2\kappa\eta^{t+1}}\lVert
z^{t}-\tilde{z}^{t}\rVert_{1}^{2},$ (17)

where ${A}_{\texttt{dual}}\coloneqq{h(f^{\ast})-\min
h+\frac{\kappa^{2}+2D^{2}}{2\kappa}}$.

##### Proof of Proposition 3.

Consider an “intermediate” point
$y^{t}\coloneqq\operatorname{\Lambda}\left(\eta^{t}w^{t+1}\right)$, we focus
on the terms $\alpha^{t}\langle C^{t},y^{t}-f^{\ast}\rangle$ and
$\alpha^{t}\langle C^{t},z^{t}-y^{t}\rangle$. These terms sum up to
$\alpha^{t}\langle C^{t},z^{t}-f^{\ast}\rangle$ which defines
$\operatorname{\mathcal{R}}_{T}(f^{\ast})$. First, from the update rule of
$w^{t}$ in Algorithm 3, we have

$\displaystyle{\alpha^{t}}\left\langle
C^{t},y^{t}\\!-\\!f^{\ast}\right\rangle\\!=\\!\frac{1}{\eta^{t}}\langle\eta^{t}w^{t}\\!-\\!\eta^{t}w^{t+1}\mathopen{},\mathopen{}y^{t}-f^{\ast}\rangle$
$\displaystyle=\frac{1}{\eta^{t}}\left[-h(y^{t})\\!-\\!h^{\ast}\left(\eta^{t}w^{t+1}\right)+\left\langle
f^{\ast},\eta^{t}w^{t+1}-\eta^{t}w^{t}\right\rangle+\left\langle
y^{t},\eta^{t}w^{t}\right\rangle\right]$
$\displaystyle\leq\frac{1}{\eta^{t}}\left[F(f^{\ast},\eta^{t}w^{t})\\!-\\!F(f^{\ast},\eta^{t}w^{t+1})\\!-\\!\mathrm{KL}(y^{t}\nonscript\,\|\nonscript\,\mathopen{}\tilde{z}^{t})\right].$

Here, the last inequality comes from the definition of Fenchel coupling and
the fact that $\tilde{z}^{t}=\operatorname{\Lambda}(\eta^{t}w^{t})$. Now,
apply the three-point inequality with Bregman divergence (Lemma 3.1 of [11])
to the KL-divergence,

$\frac{1}{\eta^{t}}\mathrm{KL}(y^{t}\nonscript\,\|\nonscript\,\mathopen{}\tilde{z}^{t})-\frac{1}{\eta^{t}}\mathrm{KL}(z^{t}\nonscript\,\|\nonscript\,\mathopen{}\tilde{z}^{t})-\frac{1}{\eta^{t}}\mathrm{KL}(y^{t}\nonscript\,\|\nonscript\,\mathopen{}z^{t})=\langle\operatorname{\nabla}h(\tilde{z})-\operatorname{\nabla}h(z),z^{t}-y^{t}\rangle\geq\alpha^{t}\left\langle\tilde{C}^{t},z^{t}-y^{t}\right\rangle.$

Combining the two inequalities derived above, we have:

$\displaystyle\operatorname{\mathcal{R}}_{T}(f^{\ast})=\sum_{t=1}^{T}\alpha^{t}\langle
C^{t},y^{t}-f^{\ast}\rangle+t\langle C^{t},z^{t}-y^{t}\rangle$
$\displaystyle\leq$
$\displaystyle\underbrace{\sum_{t=1}^{T}\frac{1}{\eta^{t}}\left[F(f^{\ast},\eta^{t}w^{t})\\!-\\!F(f^{\ast},\eta^{t}w^{t+1})\right]}_{\textstyle\coloneqq{A}_{1}}+\underbrace{\left[\sum_{t=1}^{T}{\alpha^{t}}\left\langle
C^{t}-\tilde{C}^{t},z^{t}\\!-\\!y^{t}\right\rangle\\!-\\!\sum_{t=1}^{T}\frac{1}{\eta^{t}}\mathrm{KL}(y^{t}\nonscript\,\|\nonscript\,\mathopen{}z^{t})\right]}_{\textstyle\coloneqq{A}_{2}}-\underbrace{\sum_{t=1}^{T}\frac{1}{\eta^{t}}\mathrm{KL}(z^{t}\nonscript\,\|\nonscript\,\mathopen{}\tilde{z}^{t})}_{\textstyle\coloneqq{A}_{3}}.$
(18)

Now, we look for upper-bounds of the three terms in the right-hand-side of
(18). First, we trivially have

$\displaystyle{A}_{1}\leq\frac{1}{\eta^{1}}F(f^{\ast},\eta^{1}w^{1})+\left(\frac{1}{\eta^{T+1}}-\frac{1}{\eta^{1}}\right)\left[h(f^{\ast})-\min
h\right]=\frac{1}{\eta^{T+1}}\left[h(f^{\ast})-\min h\right].$ (19)

Second, for any $T$, from the Cauchy-Schwarz inequality and the fact that
$\left\lVert f-f^{\prime}\right\rVert_{1}\left\lVert
w-w^{\prime}\right\rVert_{\infty}=\min_{>0}\left\\{\frac{missing}{2}\left\lVert
f-f^{\prime}\right\rVert_{1}^{2}+\frac{1}{2missing}\left\lVert
w-w^{\prime}\right\rVert_{\infty}^{2}\right\\}$ for any
$f,f^{\prime},w,w^{\prime}\in\mathbb{R}^{d}$,111Consider the function
$\psi()=\frac{missing}{2}\left\lVert
f-f^{\prime}\right\rVert_{1}^{2}+\frac{1}{2missing}\left\lVert
w-w^{\prime}\right\rVert_{\infty}^{2}$; then ${}^{*}=\left\lVert
f-f^{\prime}\right\rVert_{1}\left\lVert w-w^{\prime}\right\rVert_{\infty}$ is
the minimizer of $\psi$. we have

$\displaystyle{\alpha^{t}}\langle
C^{t}\\!-\\!\tilde{C}^{t}\mathopen{},\mathopen{}y^{t}\\!-\\!z^{t}\rangle\leq\left\lVert\alpha^{t}(C^{t}\\!-\\!\tilde{C}^{t})\right\rVert_{\infty}\left\lVert
y^{t}\\!-\\!z^{t}\right\rVert_{1}\leq{\frac{(\alpha^{t})^{2}\kappa\eta^{t+1}}{2}\left\lVert
C^{t}-\tilde{C}^{t}\right\rVert_{\infty}^{2}+\frac{1}{2\kappa\eta^{t+1}}\left\lVert
y^{t}-z^{t}\right\rVert_{1}^{2}}.$ (20)

Combine this with the strong convexity of $h$, we obtain:

$\displaystyle{A}_{2}\leq\frac{\kappa}{2}\sum_{t=1}^{T}{(\alpha^{t})}^{2}\eta^{t+1}\left\lVert
C^{t}\\!-\\!\tilde{C}^{t}\right\rVert_{\infty}^{2}+\frac{1}{2\kappa}\sum_{t=1}^{T}\left(\frac{1}{\eta^{t+1}}-\frac{1}{\eta^{t}}\right)\lVert
y^{t}-\tilde{z}^{t}\rVert_{1}^{2}\leq\frac{\kappa}{2}\sum_{t=1}^{T}{(\alpha^{t})}^{2}\eta^{t+1}\left\lVert
C^{t}-\tilde{C}^{t}\right\rVert_{\infty}^{2}+\frac{D^{2}}{2\kappa}\left(\frac{1}{\eta^{T+1}}-1\right).$
(21)

Third, recall the notion $D$ denoting the $\lVert\rVert_{1}$-diameter of
$\mathcal{F}$, we have

${A}_{3}\geq\sum_{t=1}^{T}\frac{1}{2\kappa\eta^{t}}\lVert
z^{t}-\tilde{z}^{t}\rVert_{1}^{2}=\sum_{t=1}^{T}\frac{1}{2\kappa\eta^{t+1}}\lVert
z^{t}-\tilde{z}^{t}\rVert_{1}^{2}-\frac{D^{2}}{2\kappa}\left(\frac{1}{\eta^{T+1}}-1\right).$
(22)

Moreover, from the update rule of $\eta^{t+1}$ and Lemma 2 of [20] (also
presented in [27, 29]), we have:

$\displaystyle\frac{1}{\eta^{T+1}}=\sqrt{1+\sum_{t=1}^{T}{(\alpha^{t})}^{2}\lVert
C^{t}-\tilde{C}^{t}\rVert_{\infty}^{2}}\leq\sum_{t=1}^{T}\frac{{(\alpha^{t})}^{2}\lVert
C^{t}-\tilde{C}^{t}\rVert_{\infty}^{2}}{\sqrt{1+\sum_{s=1}^{t}{(\alpha^{s})}^{2}\lVert
C^{s}-\tilde{C}^{s}\rVert_{\infty}^{2}}}+1=\sum_{t=1}^{T}{(\alpha^{t})}^{2}\eta^{t+1}\lVert
C^{t}-\tilde{C}^{t}\rVert_{\infty}^{2}+1.$

Apply this inequality to (22) then apply it with (19), (21) into (18), we
obtain precisely (17). ∎

Now, let us define the noises of the observed costs $\tilde{C}^{t}$ and
$C^{t}$, induced by the flow profiles $\tilde{f}^{t}$ and $f^{t}$, in
comparison with the corresponding gradients of the BMW potential $\Phi$ as
follows:

$\tilde{U}^{t}\coloneqq\tilde{C}^{t}-\nabla\Phi\left(\tilde{f}^{t}\right)\textrm{
and }U^{t}\coloneqq C^{t}-\nabla\Phi\left(f^{t}\right).$

By the definition of $\sigma$ (cf. Equation 6), we have
$\sup\left\\{\operatorname{\mathbb{E}}\left[\lVert\tilde{U}^{t}\rVert_{\infty}^{2}\right],\operatorname{\mathbb{E}}\left[\lVert
U^{t}\rVert_{\infty}^{2}\right]\right\\}\leq\sigma^{2}$. From these
definitions, we also deduce that
$\operatorname{\mathbb{E}}\left[\tilde{U}^{t}\middle|\mathcal{H}^{t-1}\right]=\operatorname{\mathbb{E}}\left[U^{t}\middle|\mathcal{H}^{t-1}\right]=0$
where
$\mathcal{H}^{t-1}=\left\\{f^{t-1},\tilde{f}^{t-1},\omega^{t-1},\ldots,f^{1},\tilde{f}^{1},\omega^{1}\right\\}$
is the filtration up to time epoch $t-1$.

Equation 17 involves the difference of the actual gradients of the BMW
potential $\nabla\Phi(f^{t})-\nabla\Phi(\tilde{f}^{t})$ and the difference of
costs at $f$ and $\tilde{f}$, i.e., $C^{t}-\tilde{C}^{t}$ (in the static case,
these differences are equal). For the sake of brevity, we define the following
terms in order to analyze the gap between these differences:

$\kappa^{t}=\min\left\\{\lVert\nabla\Phi(f^{t})-\nabla\Phi(\tilde{f}^{t})\rVert_{\infty}^{2},\lVert
C^{t}-\tilde{C}^{t}\rVert_{\infty}^{2}\right\\}\textrm{ and
}\xi^{t}=\left[C^{t}-\tilde{C}^{t}\right]-\left[\nabla\Phi(f^{t})-\nabla\Phi(\tilde{f}^{t})\right].$

We aim to construct the upper-bounds the last two terms in the right-hand-side
of (17) in terms of $\kappa^{t}$ and $\xi^{t}$. Particularly, from (17), we
can prove the following proposition:

######  Proposition 4.

Run Algorithm 3 and define
$\tilde{\eta}^{t}\coloneqq{1}\big{/}{\sqrt{1+2\sum_{s=1}^{t-1}(\alpha^{s})^{2}\kappa^{s}}}$,
we have

$\operatorname{\mathcal{R}}_{T}(f^{\ast})\leq h(f^{\ast})-\min
h+{A}_{\texttt{stoch}}2\sqrt{2}\left(\sum_{t=1}^{T}(\alpha^{t})^{2}\lVert\xi_{t}\rVert_{\infty}^{2}\right)^{1/2}+\sum_{t=1}^{T}\textsl{g}(\tilde{\eta}^{t+1})(\alpha^{t})^{2}\kappa^{t},$
(23)

where ${A}_{\texttt{stoch}}\coloneqq h(f^{\ast})-\min
h+\frac{\kappa^{2}+3D^{2}}{2\kappa}$ and
$\textsl{g}(\tilde{\eta}^{t+1})\coloneqq
4{A}_{\texttt{stoch}}\tilde{\eta}^{t+1}-\frac{1}{8\beta^{2}\kappa\tilde{\eta}^{t+1}}$.

##### Proof of Proposition 4.

First, by definitions of $\kappa^{t}$ and $\xi^{t}$, we have $\lVert
C^{t}-\tilde{C}^{t}\rVert_{\infty}^{2}\leq
2\kappa^{t}+2\lVert\xi^{t}\rVert_{\infty}^{2}$. Apply this and Lemma 2 of
[20], we have

$\displaystyle\sum_{t=1}^{T}(\alpha^{t})^{2}\eta^{t+1}\lVert
C^{t}-\tilde{C}^{t}\rVert_{\infty}^{2}$
$\displaystyle=\sum_{t=1}^{T}\frac{(\alpha^{t})^{2}\lVert
C^{t}-\tilde{C}^{t}\rVert_{\infty}^{2}}{\sqrt{1+\sum_{s=1}^{t}\alpha^{s}\lVert
C^{s}-\tilde{C}^{s}\rVert_{\infty}^{2}}}$ $\displaystyle\leq
2\sqrt{1+\sum_{t=1}^{T}(\alpha^{t})^{2}\lVert
C^{t}-\tilde{C}^{t}\rVert_{\infty}^{2}}$
$\displaystyle\leq\sqrt{1+2\sum_{t=1}^{T}(\alpha^{t})^{2}\kappa^{t}}+\sqrt{2\sum_{t=1}^{T}(\alpha^{t})^{2}\lVert\xi^{t}\rVert_{\infty}^{2}}$
$\displaystyle\leq
4\sum_{t=1}^{T}\frac{(\alpha^{t})^{2}\kappa^{t}}{\sqrt{1+2\sum_{s=1}^{t}(\alpha^{s})^{2}\kappa^{t}}}+\sqrt{2\sum_{t=1}^{T}(\alpha^{t})^{2}\lVert\xi^{t}\rVert_{\infty}^{2}}.$
(24)

On the other hand, from the definition of $\tilde{\eta}^{t}$, we have
$\frac{1}{\tilde{\eta}^{t}}\leq\frac{1}{\eta^{t}},\forall t$ and hence,

$\displaystyle\sum_{t=1}^{T}\frac{1}{\tilde{\eta}^{t+1}}\left\lVert\tilde{z}^{t}-z^{t}\right\rVert_{1}^{2}\leq\sum_{t=1}^{T}\left(\frac{1}{\eta^{t+1}}-\frac{1}{\eta^{t}}\right)\left\lVert
z^{t}-\tilde{z}^{t}\right\rVert_{1}^{2}+\frac{1}{\eta^{t}}\left\lVert
z^{t}-\tilde{z}^{t}\right\rVert_{1}^{2}\leq\left(\frac{1}{\eta^{T+1}}-1\right)D^{2}+\sum_{t=1}^{T}\frac{1}{\eta^{t}}\left\lVert
z^{t}-\tilde{z}^{t}\right\rVert_{1}^{2}.$

Combine this with the fact that
$z^{t}-\tilde{z}^{t}=\frac{\alpha^{t+1}}{2}\left(f^{t}-\tilde{f}^{t}\right)$
(from Lines 3 and 3 of Algorithm 3) and choose an increasing sequence of
$\alpha^{t}$ (such as $\alpha^{t}=t$), we have

$\displaystyle-\sum_{t=1}^{T}\frac{1}{2\kappa\eta^{t+1}}\lVert
z^{t}-\tilde{z}^{t}\rVert_{1}^{2}\leq$
$\displaystyle-\frac{1}{2\kappa}\sum_{t=1}^{T}\frac{1}{\tilde{\eta}^{t+1}}\left\lVert
z^{t}-\tilde{z}^{t}\right\rVert_{1}^{2}+\frac{D^{2}}{2\kappa}\left(\frac{1}{\eta^{T+1}}-1\right)$
$\displaystyle\leq$
$\displaystyle-\frac{1}{2\kappa}\sum_{t=1}^{T}\frac{\alpha^{t+1}}{2\tilde{\eta}^{t+1}}\left\lVert
f^{t}-\tilde{f}^{t}\right\rVert_{1}^{2}+\frac{D^{2}}{2\kappa}\left(\frac{1}{\eta^{T+1}}-1\right)$
$\displaystyle\leq$
$\displaystyle-\frac{1}{2\kappa}\sum_{t=1}^{T}\frac{1}{\tilde{\eta}^{t+1}}\frac{(\alpha^{t+1})^{2}\left\lVert\nabla\Phi(f^{t})-\nabla\Phi(\tilde{f}^{t})\right\rVert_{\infty}^{2}}{4\beta^{2}}+\frac{D^{2}}{2\kappa}\left(\frac{1}{\eta^{T+1}}-1\right)$
$\displaystyle\leq$
$\displaystyle-\frac{1}{8\beta^{2}\kappa}\sum_{t=1}^{T}\frac{(\alpha^{t})^{2}\kappa^{t}}{\tilde{\eta}^{t+1}}+\frac{D^{2}}{2\kappa}\left(\frac{1}{\eta^{T+1}}-1\right).$
(25)

Apply (24) and (25) into (17), we obtain (23) and finish the proof of
Proposition 4. ∎

We now focus in the last term of (23). We denote
$T_{0}\coloneqq\max\left\\{1\leq t\leq
T:\tilde{\eta}^{t+1}\geq\left[32\beta^{2}\kappa{A}_{\texttt{stoch}}\right]^{-1/2}\right\\}$.
Then for any $t\geq T_{0}$, we have $\textsl{g}(\tilde{\eta}^{t+1})<0$. In
other words, in
$\sum_{t=1}^{T}\textsl{g}(\tilde{\eta}^{t+1})(\alpha^{t})^{2}\kappa^{t}$, only
the first $T_{0}$ components are positive. Therefore,

$\displaystyle\sum_{t=1}^{T}\textsl{g}(\tilde{\eta}^{t+1})(\alpha^{t})^{2}\kappa^{t}\leq\sum_{t=1}^{T_{0}}\textsl{g}(\tilde{\eta}^{t+1})(\alpha^{t})^{2}\kappa^{t}\leq$
$\displaystyle\sum_{t=1}^{T_{0}}4{A}_{\texttt{stoch}}\tilde{\eta}^{t+1}(\alpha^{t})^{2}\kappa^{t}$
$\displaystyle=$ $\displaystyle
4{A}_{\texttt{stoch}}\cdot\sum_{t=1}^{T_{0}}\frac{(\alpha^{t})^{2}\kappa^{t}}{\sqrt{1+2\sum_{s=1}^{t}(\alpha^{s})^{2}\kappa^{s}}}$
$\displaystyle\leq$ $\displaystyle 4{A}_{\texttt{stoch}}\cdot
2\sqrt{1+2\sum_{t=1}^{T_{0}}(\alpha^{t})^{2}\kappa^{t}}$ $\displaystyle=$
$\displaystyle\frac{8}{\tilde{\eta}^{T_{0}}}{A}_{\texttt{stoch}}$
$\displaystyle\leq$ $\displaystyle
32\beta\sqrt{2\kappa}({A}_{\texttt{stoch}})^{3/2}.$ (26)

Combine (26), (23) and the fact that
$\operatorname{\mathbb{E}}\left[\lVert\xi_{t}\rVert_{\infty}^{2}\right]\leq
2\operatorname{\mathbb{E}}\left[\lVert
U^{t}\rVert_{\infty}^{2}\right]+2\operatorname{\mathbb{E}}\left[\lVert\tilde{U}^{t}\rVert_{\infty}^{2}\right]\leq
4\sigma^{2}$, we have

$\displaystyle\operatorname{\mathbb{E}}\left[\operatorname{\mathcal{R}}_{T}(f^{\ast})\right]\leq
h(f^{\ast})-\min
h+\sigma{A}_{\texttt{stoch}}4\sqrt{2\sum\nolimits_{t=1}^{T}(\alpha^{t})^{2}}+32\beta\sqrt{2\kappa}({A}_{\texttt{stoch}})^{3/2}.$
(27)

Finally, in order to apply (15), we need to make the connection between
$\operatorname{\mathbb{E}}\left[\operatorname{\overline{\operatorname{\mathcal{R}}}}_{T}(f^{\ast})\right]$
and
$\operatorname{\mathbb{E}}\left[\operatorname{\mathcal{R}}_{T}(f^{\ast})\right]$.
Particularly, we have:

$\displaystyle\operatorname{\mathbb{E}}\left[\operatorname{\overline{\operatorname{\mathcal{R}}}}_{T}(f^{\ast})\right]=\operatorname{\mathbb{E}}\left[\operatorname{\mathcal{R}}_{T}(f^{\ast})\right]-\operatorname{\mathbb{E}}\left[\sum_{t=1}^{T}{\alpha^{t}}\left\langle
C^{t}-\operatorname{\nabla}\Phi(f^{t}),z^{t}-f^{\ast}\right\rangle\right]$
$\displaystyle=\operatorname{\mathbb{E}}\left[\operatorname{\mathcal{R}}_{T}(f^{\ast})\right]-\operatorname{\mathbb{E}}\left[\sum_{t=1}^{T}{\alpha^{t}}\left\langle
U^{t},z^{t}-f^{\ast}\right\rangle\right]$
$\displaystyle=\operatorname{\mathbb{E}}\left[\operatorname{\mathcal{R}}_{T}(f^{\ast})\right].$
(28)

Here, the last equality comes from the fact that
$\operatorname{\mathbb{E}}\left[{\alpha^{t}}\left\langle
U^{t},z^{t}-f^{\ast}\right\rangle\right]=\operatorname{\mathbb{E}}\left[{\alpha^{t}}\left\langle\operatorname{\mathbb{E}}\left[U^{t}|\mathcal{H}^{t-1}\right],z^{t}-f^{\ast}\right\rangle\right]=0$
(by law of total expectation). Combine (28) and (27), then apply (15) with the
choice of step-size $\alpha^{t}=t,\forall t$, we have:

$\displaystyle\operatorname{\mathbb{E}}\left[\Phi(f^{T})-\Phi(f^{\ast})\right]\leq\frac{1}{T^{2}}\left(h(f^{\ast})-\min
h+\sigma{A}_{\texttt{stoch}}4\sqrt{2}T^{3/2}+32\beta\sqrt{2\kappa}({A}_{\texttt{stoch}})^{3/2}\right).$
(29)

Finally, recall that $h(f^{\ast})\leq M_{\mathrm{tot}}\log(M_{\max})$, $-\min
h\leq M_{\mathrm{tot}}\log(P/M_{\mathrm{tot}})$, $\kappa\coloneqq NM_{\max}$,
$D\leq 2M_{\mathrm{tot}}$, $\sigma\leq 2H$ and the definition of
${A}_{\texttt{stoch}}$ (cf. Proposition 4), we have
${A}_{\texttt{stoch}}\leq\frac{1}{2}NM_{\max}\left[2\log\left(\frac{PM_{\max}}{M_{\mathrm{tot}}}\right)+13\right]\coloneqq\frac{1}{2}{A}$
for the constant ${A}$ defined in Theorem 3. Plug these results into (29), we
obtain (14) and conclude the proof.

##  5\. AdaLight: Adaptive learning with efficient per-iteration complexity

In this section, our main focus is to design an algorithm that not only has
the adaptive optimality of AdaWeight but also has a per-iteration complexity
that is scalable in the network’s size. To do this, in Section 5.1, we first
introduce an alternative routing paradigm, called local flows. In Section 5.2
and Section 5.3, we then propose the AdaLight algorithm and analyze its
convergence properties and its per-iteration complexities.

In the sequel, we will additionally use the following set of notation. For any
(directed) edge $e\in\mathcal{E}$, we denote by $v^{+}_{e}$ and $v^{-}_{e}$
the tail vertex and head vertex of $e$ respectively, i.e., $e$ goes from
$v^{+}_{e}$ to $v^{-}_{e}$. For any vertex $v$ and any sub-graph
$\mathcal{G}^{i}$, let ${\text{In}^{i}_{v}}$ and ${\text{Out}^{i}_{v}}$
respectively denote the set of incoming edges to $v$ and the set of outgoing
edges from $v$, and let $\text{Child}^{i}_{v}$ and $\text{Parent}^{i}_{v}$
denote the set of its direct successors and the set of its direct
predecessors. For the sake of brevity, we also write
$\dim^{i}_{v}\coloneqq\left\lvert{\text{Out}^{i}_{v}}\right\rvert$ which is
the out degree of vertex $v$ in $\mathcal{G}^{i}$ and
$\dim_{v}\coloneqq\sum_{i\in\mathcal{N}}\dim^{i}_{v}$.

###  5.1. dis routing paradigm

On the defined network structure (cf. Section 2), let us consider a novel
route-recommendation paradigm that unfolds as follows:

  1. i.5pt)

The control interface determines at each vertex $v$ a profile
$x^{i}_{v}\in\Delta({\text{Out}^{i}_{v}})$ for each O/D pair $i$.

  2. ii.5pt)

It then makes a route recommendation for each $i$-type traffic by choosing the
vertices (in $\mathcal{G}^{i}$) one by one such that if the vertices
$O^{i},v_{1},\ldots,v_{k}$ is chosen, then a vertex
$v_{k+1}\in{\text{Out}^{i}_{v_{k}}}$ will be chosen with probability
$x^{i}_{v_{k},e}$ where $e$ is the edge going from $v_{k}$ to $v_{k+1}$.

In this paradigm, the recommendations are made based on “local” decisions at
each vertex. In average, among the mass of $i$-type traffic arriving at a
vertex $v$, the proportion that is routed on edge $e\in{\text{Out}^{i}_{v}}$
is precisely $x^{i}_{v,e}$. We formalize this idea via the notion of local
flow defined as follows.

#####  dis.

For each pair $i$ and vertex $v$, we define the set
${\mathcal{X}}_{v}^{i}\coloneqq\big{\\{}x_{v}^{i}\in[0,1]^{\dim^{i}_{v}}:\sum\nolimits_{e\in{\text{Out}^{i}_{v}}}x^{i}_{v,e}=1\big{\\}}$.
We call the elements of ${\mathcal{X}}_{v}^{i}$ the local flows corresponding
to $i$ and $v$. We also use the notations
${\mathcal{X}}_{v}\coloneqq\prod_{i\in\mathcal{N}}{\mathcal{X}}^{i}_{v}$ for
the set of local flow profiles chosen at $v$ across all O/D pairs,
${\mathcal{X}}^{i}\coloneqq\prod_{v\in\mathcal{V}}{\mathcal{X}}^{i}_{v}$ for
the set of local flow profiles of $i$-type traffic across the vertices and
${\mathcal{X}}\coloneqq\prod_{v\in\mathcal{V}}{\mathcal{X}}_{v}=\prod_{i\in\mathcal{N}}{\mathcal{X}}^{i}$
for the set of all local flow profiles.

##### Traffic loads induced by local flows.

When the network’s traffic is routed according to a local flow profile
$x\in{\mathcal{X}}$, we denote the _mass_ of $i$-type traffic arriving at $v$
by $\mu^{i}_{v}(x)$ and denote the _load_ of $i$-type traffic routed on $e$ by
${\ell}^{i}_{e}(x)$. Formally, they can be computed (recursively) as follows:

${\ell}^{i}_{e}(x)=x^{i}_{v,e}\cdot\mu^{i}_{v}(x)\textrm{,~{}and~{}}\mu^{i}_{v}(x)\coloneqq\left\\{\begin{array}[]{ll}M^{i}&\textrm{,
if }v=O^{i}\\\
\sum_{e^{\prime}\in{\text{In}^{i}_{v}}}{\ell}^{i}_{e^{\prime}}(x)&\textrm{, if
}v\neq O^{i}\\\ \end{array}\right.$ (30)

We define the total load induced on $e$ as
${\ell}_{e}\left(x\right)=\sum_{i\in\mathcal{N}}{\ell}^{i}_{e}(x)$.

The congestion cost on each edge $e$ is determined by the cost function
$c_{e}$ as defined in Section 2. Particularly,
$c_{e}\left({\ell}_{e}(x),\omega\right)$ is the cost on $e$ at state
$\omega\in\Omega$ when the local flow profile $x\in{\mathcal{X}}$ is
implemented.

Finally, corresponding to the local routing paradigm described above, we re-
adjust the learning model at each time epoch $t=1,2,\dotsc\,$ as follows:

  1. i.5pt)

The navigation chooses a local flow profile $x^{t}\in{\mathcal{X}}$ and routes
the traffic accordingly.

  2. ii.5pt)

Concurrently, the state $\omega^{t}$ of the network is drawn from $\Omega$
(i.i.d. relative to $\operatorname{\mathbb{P}}$).

  3. iii.5pt)

The congestion costs $c_{e}\left({\ell}_{e}(x^{t}),\omega^{t}\right)$ on each
edge $e$ is observed.

###  5.2. AdaLight algorithm

In this section, we propose a new equilibrium learning method called AdaLight.
This method will be implemented in the local flow learning model described in
the previous section.

To build AdaLight, our starting point is the AdaWeight method. In Algorithm 3,
there are two major types of flows updating: the flows $\tilde{z}^{t}$ and
$z^{t}$ outputting from the logit mapping $\operatorname{\Lambda}$, and the
flows $\tilde{f}^{t}$ and $f^{t}$ obtained from averaging between two pre-
computed flows. In Algorithm 3, it takes $\operatorname{\mathcal{O}}(P)$
computations to update each of these flow profiles in each iteration. To
improve the per-iteration complexity, we can leverage the weight-pushing
technique [42, 18]: by assigning weights on edges and use dynamic programming
principles, we can find the local flow profiles inducing the same loads (and
costs) as the logit-mapped flows $\tilde{z}^{t}$ and $z^{t}$. However, the
weight-pushing technique fails to derive local flow profiles matching the
averaged flows $\tilde{f}^{t}$ and $f^{t}$ of Algorithm 3. Particularly,
weight-pushing two flows and then taking the average will not incur the same
costs as implementing the average of the two flows. Recall that these
averaging steps are the key elements allowing the adaptability of AdaWeight,
the key challenge now is to implement efficiently these steps.

Facing up this challenge, we make the following observation: _while the
averaged flows of AdaWeight are not “weight-pushable”, their induced loads
are_. In this perspective, there arise two new challenges. First, we need to
efficiently compute the load profiles induced by the averaged flows (without
explicitly computing these flows). To do this, we introduce a sub-routine
called pulling-forward – reflecting the fact that we start from the origin and
work forward (unlike the classical weight-pushing that starts from the
destination and works backward). Second, we need to derive local flow profiles
that matches these averaged loads. We refer to this as matching-load. Note
that _pulling-forward and matching-load are completely novel contributions_.

Taking an overall view, AdaLight is a combination of the weight-pushing
technique (to handle the logit-mapped flows) and the pulling-forward and
matching-load procedures (to handle the averaged flows) into the AdaWeight
template. For the sake of conciseness, we first present a sub-algorithm in
Section 5.2.1 and present a pseudo-code form of AdaLight in Section
5.2.2.222As a side note on AcceleWeight (cf. Section 3.2), a direct
application of weight-pushing also fails to achieve an efficient
implementation for the same reasons of that of AdaWeight. Note that we can use
techniques that we introduce in AdaLight to improve AcceleWeight and achieve
an efficient version. As AcceleWeight is not the focus of our work, we omit
the details.

####  5.2.1. PushPullMatch sub-algorithm

In this section, we present a sub-algorithm, called PushPullMatch, that
combines the three auxiliary routines: pushing-backward, pulling-forward and
matching-load. In Algorithm 4, we give a pseudo-code form of PushPullMatch. It
takes four inputs: an “anchor local-load” profile
$\textsf{r}=(\textsf{r}^{i}_{e})_{i\in\mathcal{N},e\in\mathcal{E}^{i}}$, a
“local-weight” profile
$\textsf{w}=(\textsf{w}^{i}_{e})_{i\in\mathcal{N},e\in\mathcal{E}^{i}}$; and
two scalar numbers $\alpha^{t}$ and $\sum_{s=0}^{t}\alpha^{s}$ (we use these
notations to facilitate the presentation of AdaLight in Section 5.2). When it
finishes, Algorithm 4 outputs a local flow profile and updates the anchor
local-load profile.

1

2

Input: $\textsf{r},\textsf{w}\in\mathbb{R}^{\sum_{i}\sum_{v}\dim^{i}_{v}}$ and
$\alpha^{t},\sum_{s=0}^{t}\alpha^{s}>0$

Output: $x\in{\mathcal{X}}$ and update r

3 In parallel for all  _$i\in\mathcal{N}$_ do

4 Fix an topological order
$v_{0},v_{1},\ldots,v_{\lvert\mathcal{V}^{i}\rvert}$ of the graph
$\mathcal{G}^{i}$ (such that $v_{0}\coloneqq O^{i}$ and
$v_{\lvert\mathcal{V}^{i}\rvert}\coloneqq D^{i}$)

// Pushing-backward phase

5 for  _$v=v_{\lvert\mathcal{V}^{i}\rvert},\ldots,v_{0}$_ do

6 if  _$v=D^{i}$_ then set $\texttt{BWscore}_{v}\leftarrow 1$

7 // initialize BWscore at destination else set
$\texttt{BWscore}_{v}\leftarrow\sum_{e\in{\text{Out}^{i}_{v}}}{\texttt{BWscore}_{v^{-}_{e}}\cdot\exp\left(\textsf{w}^{i}_{e}\right)}$

// compute BWscore based on BWscore of children set
$z^{i}_{v,e}\leftarrow\exp\left(\textsf{w}^{i}_{e}\right){\texttt{BWscore}_{v^{-}_{e}}}\big{/}{\texttt{BWscore}_{v}}$
for $e\in{\text{Out}^{i}_{v}}$

// compute pivot outgoing flow

8

// Pulling-forward phase

9 for  _$v=v_{0},\ldots,v_{\lvert\mathcal{V}^{i}\rvert}$_ do

10 if  _$v=O^{i}$_ then set $\texttt{FWscore}_{v}\leftarrow M^{i}$

11 // initialize FWscore at origin else set
$\texttt{FWscore}_{v}\leftarrow\sum_{e^{\prime}\in{\text{In}^{i}_{v}}}{{\ell}^{i}_{e^{\prime}}(z)}$

12 // compute FWscore from loads of incoming edges for
_$e\in{\text{Out}^{i}_{v}}$_ do

set ${\ell}^{i}_{e}(z)\leftarrow\texttt{FWscore}_{v}\cdot z^{i}_{v,e}$

// pull the loads forward

set
${\ell}^{i}_{e}(x)\leftarrow\left[\alpha^{t}{\ell}^{i}_{e}\left(z\right)+\textsf{r}^{i}_{e}\right]\big{/}\sum_{s=0}^{t}\alpha^{s}$

// average $\&$ pull forward

set
$\textsf{r}^{i}_{e}\leftarrow\textsf{r}^{i}_{e}+\alpha^{t}\cdot{\ell}^{i}_{e}(z)$

// update the anchor

13

14

15

// Matching-load phase

16 for  _$v=v_{0},\ldots,v_{\lvert\mathcal{V}^{i}\rvert}$_ do

17 if  _$v=O^{i}$_ then set $\mu^{i}_{v}({x})\leftarrow M^{i}$

18 // initialize mass arriving at origin else set
$\mu^{i}_{v}({x})\leftarrow\sum_{e^{\prime}\in{\text{In}^{i}_{v}}}{{\ell}^{i}_{e^{\prime}}(x)}$

// compute mass that would arrive at $v$ if $x$ is implemented set
$x^{i}_{v,e}\leftarrow{\ell}^{i}_{e}(x)\big{/}\mu^{i}_{v}({x})$ for
$e\in{\text{Out}^{i}_{v}}$

// compute outgoing flow matching ${\ell}^{i}_{e}(x)$

19

20

21

Algorithm 4
$\textsc{Push\\-Pull\\-Match}(\textsf{r},\textsf{w},\alpha^{t},\sum_{s=0}^{t}\alpha^{s})$

Particularly, when we run PushPullMatch, it executes the following three
phases for each O/D pair $i$:

  1. i.5pt)

The pushing-backward phase: in this phase, we consider vertices in
$\mathcal{G}^{i}$ one by one in a _reversed topological order_ , i.e., we work
backwardly from the destination to the origin.333Note that since
$\mathcal{G}^{i}$ are DAGs, such an topological order exists and can be found
in $\operatorname{\mathcal{O}}{\lvert\mathcal{E}^{i}\rvert}$ time. For each
vertex $v$, we assign a score – called the _backward score_ (denoted by
BWscore in Algorithm 4) – depending on the weights of $v$’s outgoing edges
(i.e., $\textsf{w}^{i}_{e},\forall e\in{\text{Out}^{i}_{v}}$) and its
children’s scores. Then, we “pushes” the backward score of $v$ to its parents
so that, in their turns, we can compute their scores. Finally, based on the
backward score of $v$, we compute a local flow
$z^{i}_{v}\in{\mathcal{X}}^{i}_{v}$.

  2. ii.5pt)

The pulling-forward phase: in this phase, we consider vertices in
$\mathcal{G}^{i}$ one by one in a _topological order_ , i.e., we work
forwardly from the origin to the destination. On each vertex $v$, we assign a
score – called the _forward score_ (denoted by FWscore in Algorithm 4) –
computed from the loads induced by $z$ on the incoming edges of $v$ that are
“pulled" from $v$’s parents. Then, for any $e\in{\text{Out}^{i}_{v}}$, we
derive $\ell^{i}_{e}(z)$ and compute another term – called $\ell^{i}_{e}(x)$ –
by taking the average of ${\ell}^{i}_{e}(z)$ and the anchor load
$\textsf{r}^{i}_{e}$. Finally, we update the anchor load-profile r (that will
be used later in AdaLight).

  3. iii.5pt)

The matching-load phase: in this phase, we once again consider vertices in
$\mathcal{G}^{i}$ in a topological order. For each $v$, we compute a local
flow $x^{i}_{v}\in{\mathcal{X}}^{i}_{v}$ that matches precisely the load
$\ell^{i}_{e}(x),\forall e\in{\text{Out}^{i}_{v}}$ obtained in the pulling-
forward phase. To do this, we need to compute $\mu^{i}_{v}(x)$ – the mass of
$i$-type traffic arriving at $v$ induced by $x$ – this can be done by pulling
the loads corresponding to $x^{i}$ from the parents of $v$. The local flow
profiles $x^{i},i\in\mathcal{N}$ constitute the output of PushPullMatch.

####  5.2.2. Pseudo-code of AdaLight.

Having all the necessary preparations, we now present the AdaLight method. The
key difference that makes AdaLight stand out from AdaWeight is that it uses
the sub-algorithm PushPullMatch to efficiently compute the corresponding local
flows instead of working with the $P$-dimensional flow profiles as in
AdaWeight. We present a pseudo-code form of AdaLight in Algorithm 5. Similar
to AdaWeight, AdaLight follows two phases: the test phase and the
recommendation phase:

  1. i.5pt)

In the _test phase_ , we run PushPullMatch with the anchor load profile
$\textsf{r}^{t}$ and the weight profile $\eta^{t}\textsf{w}^{t}$ as inputs
(here, $\eta^{t}$ is a learning rate). At the end of the test phase, the cost
$\tilde{C}^{t}$ at a “test” local flow is derived and used to update the test-
weight profile $\tilde{\textsf{w}}^{t}$.

  2. ii.5pt)

In the _recommendation phase_ , we re-run PushPullMatch but this time, with
the test weight $\eta^{t}\tilde{\textsf{w}}^{t}$ obtained previously in the
test phase as input. Moreover, we will also use the output of PushPullMatch to
update the anchor load-profile $\textsf{r}^{t}$. At the end of the
recommendation phase, we obtain the local flow profile $x^{t}$, route the
traffic accordingly, then update the weight profile $\textsf{w}^{t}$ by using
the incurred costs.

Finally, we observe that in each iteration of Algorithm 5, it is required to
compute the term $\max_{p\in\mathcal{P}}\left(\sum_{e\in p}\lvert
C^{t}_{e}-\tilde{C}^{t}_{e}\rvert\right)$ for updating the learning rate
$\eta^{t}$. This can also be done in efficiently: for any $i\in\mathcal{N}$,
we solve a weighted-shortest-path problem on $\mathcal{G}^{i}$ where each edge
$e$ is assigned with a weight ${}_{e}\coloneqq-\left\lvert
C^{t}_{e}-\tilde{C}^{t}_{e}\right\rvert$. For instance, this step can be done
by the classical Bellman-Ford algorithm [7, 17] that takes only
$\operatorname{\mathcal{O}}(\lvert\mathcal{V}\rvert)$ rounds of
computations.444There exist more complicated shortest-path algorithms that
have better complexities (e.g., [4, 19, 15]); however, the complication of
these algorithms are beyond the purpose of our work and they do not improve
the complexity of AdaLight in general. In this work, we only analyze AdaLight
with the simple Bellman-Ford algorithm.

1

2

3Initialize $\textsf{w}^{1}\leftarrow 0$, $\textsf{r}^{1}=0$,
$\alpha^{0}\leftarrow 0$ and $\eta^{1}\leftarrow 1$

4 for  _$t=1,2,\dotsc\,$_ do

// Test phase

set $\tilde{x}^{t}\leftarrow$ output of
$\textsc{Push\\-Pull\\-Match}(\textsf{r}^{t},\eta^{t}\textsf{w}^{t},\alpha^{t},\sum_{s=0}^{t}\alpha^{s})$

// compute a test local flow

get
$\tilde{C}^{t}_{e}\\!\leftarrow\\!c_{e}\left({\ell}_{e}\left(\tilde{x}^{t}\right),\omega^{t}\right)$
for $e\in\mathcal{E}$

// querry the test local flow

set
$\tilde{\textsf{w}}^{t,i}_{e}\leftarrow\textsf{w}^{t,i}_{e}-\alpha^{t}\tilde{C}^{t}_{e}$
for any $i\in\mathcal{N}$ and $e\in\mathcal{E}^{i}$

// update test weight

5

// Recommendation phase

set $x^{t},\textsf{r}^{t+1}\leftarrow$ outputs of
$\textsc{Push\\-Pull\\-Match}(\textsf{r}^{t},\eta^{t}\tilde{\textsf{w}}^{t},\alpha^{t},\sum_{s=0}^{t}\alpha^{s})$

// compute local flow $\&$ update anchor

Route according to $x^{t}$ and get
$C^{t}_{e}\\!\leftarrow\\!c_{e}\left({\ell}_{e}\left(x^{t}\right),\omega^{t}\right)$
for $e\in\mathcal{E}$

// route, measure costs

set $\textsf{w}^{t+1,i}_{e}\leftarrow\textsf{w}^{t,i}_{e}-\alpha^{t}C^{t}_{e}$
for any $i\in\mathcal{N}$ and $e\in\mathcal{E}^{i}$

// update weight

// Update the learning rate

6 set $\eta^{t+1}\leftarrow
1\big{/}\sqrt{1+\sum_{s=0}^{t}\alpha^{s}\max_{p\in\mathcal{P}}\left(\sum_{e\in
p}\lvert C^{s}_{e}-\tilde{C}^{s}_{e}\rvert\right)^{2}}$

Algorithm 5 ALW (AdaLight)

###  5.3. AdaLight: convergence results and per-iteration complexities

The convergence properties of AdaLight is formally presented in the following
theorem:

######  Theorem 4.

Let $x^{1},x^{2},\ldots$ be the sequence of local flow profiles recommended by
AdaLight (i.e., Algorithm 5) running with $\alpha^{t}=t,\forall
t=1,2,\dotsc\,$; then the following results hold:

  1. i.5pt)

The sequence of flow profiles $f^{t},\forall t=1,2,\dotsc\,$ where
$f^{t,i}_{p}\coloneqq M^{i}\prod_{e\in p}x^{t,i}_{v^{+}_{e},e}$ for any
$i\in\mathcal{N}$ and $p\in\mathcal{P}$ enjoys the following equilibrium
convergence rate:

$\operatorname{\mathbb{E}}\left[\operatorname{Gap}\left(T\right)\right]\leq\operatorname{\mathcal{O}}\left(\left(\log
P\right)^{{3}/{2}}\left(\frac{\sigma}{\sqrt{T}}+\frac{1}{T^{2}}\right)\right),$

and specifically, in the static case, it enjoys the convergence rate
${\operatorname{Gap}\left(T\right)}\leq\operatorname{\mathcal{O}}\left(\log(P)^{3/2}/T^{2}\right)$.

  2. ii.5pt)

Moreover, each iteration of AdaLight requires only an
$\operatorname{\mathcal{O}}\left(\lvert\mathcal{N}\rvert\lvert\mathcal{V}\rvert\lvert\mathcal{E}\rvert\right)$
number of computations.

Result (i) of Theorem 4 shows that AdaLight also converges toward an
equilibrium with the same rate as AdaWeight: an
$\operatorname{\mathcal{O}}((\log P)^{3/2}/\sqrt{T})$ rate in the stochastic
regime and an $\operatorname{\mathcal{O}}((\log P)^{3/2}/{T}^{2})$ rate in the
static regime. More precisely, we shall see in the proof of Theorem 4 that
AdaLight recommends a sequence of local flows that induce the same BMW-
potential values as the sequence of flows recommended by AdaWeight which, in
turn, approaches the potential value at an equilibrium (see Proposition 5
below). Note that although the flow $f^{t}$ is mentioned in Theorem 4, it is
never computed by AdaLight and it is not needed for routing the traffic in
practice (only the local flow $x^{t}$ is needed).

Moreover, Result (ii) of Theorem 4 shows the main difference between using
AdaLight and using AdaWeight: the per-iteration (space and time) complexities
of AdaLight are polynomial in terms of the network’s primitive parameters
(numbers of O/D pairs, numbers of vertices and numbers of edges). Therefore,
unlike AdaWeight, AdaLight can run efficiently even in large networks.

###  5.4. Proof of Theorem 4

##### First, we prove that $x^{T}$ and $f^{T}$ – as defined in Theorem 4 –
induce the same costs.

In fact, we can prove a stronger result as follows:

######  Proposition 5.

Any local flow profile $x\in{\mathcal{X}}$ and the flow profile
$f\in\mathcal{F}$ such that $f^{i}_{p}\coloneqq M^{i}\prod_{e\in p}x^{i}_{e}$
induce the same load profiles, i.e., $\ell_{e}(x)=\ell_{e}(f)$ for any $e$.

###### Proof.

_Proof of Proposition 5._ Fix an O/D pair $i\in\mathcal{N}$. Let
$\mathcal{P}^{O^{i}}_{v}$ and $\mathcal{P}^{v}_{D^{i}}$ denotes the set of
paths in $\mathcal{G}^{i}$ going from $O^{i}$ to $v$ and the set of paths
going from $v$ to $D^{i}$ respectively. Let us denote
$\texttt{H}^{i}_{v}\coloneqq\sum_{p\in\mathcal{P}^{v}_{D^{i}}}\prod_{e\in
p}x^{i}_{v^{+}_{e},e}$ for any $v\in\mathcal{V}^{i}\backslash\\{D^{i}\\}$ and
conventionally set $\texttt{H}^{i}_{D^{i}}=1$. By induction (following any
topological order of vertices), we can prove that $\texttt{H}^{i}_{v}=1$ for
any $v\in\mathcal{V}^{i}$. From this result, the definition of the load
profile in (30) and the definition of $\texttt{H}^{i}_{v}$, the following
equality holds for any $v$ and $e\in{\text{Out}^{i}_{v}}$:

$\displaystyle\ell^{i}_{e}(x)=\mu^{i}_{v}(x)\cdot
x^{i}_{e}\cdot\texttt{H}^{i}_{v^{-}_{e}}$
$\displaystyle=M^{i}\sum_{q\in\mathcal{P}^{O^{i}}_{v}}\prod_{e^{\prime}\in
q}x_{e^{\prime}}\cdot
x_{e}\cdot\sum_{p\in\mathcal{P}^{v^{-}_{e}}_{D^{i}}}\prod_{e^{\prime\prime}\in
p}x_{e^{\prime\prime}}=M^{i}\sum_{p\in\mathcal{P}_{i}\atop p\ni
e}\prod_{e^{\prime}\in p}x_{e^{\prime}}=\sum_{p\in\mathcal{P}_{i}\atop p\ni
e}f^{i}_{p}=\ell^{i}_{e}(f).$

∎∎

Apply Proposition 5, since $x^{T}$ and its corresponding flow $f^{T}$ induce
the same load profiles, by definition of the cost functions, they induce the
same costs, i.e., $c(f^{T},\omega)=c(x^{T},\omega)$ for any state
$\omega\in\Omega$.

##### Second, we prove that the recommendations of AdaLight and AdaWeight
coincide.

Let us first focus on the recommendation phase of Algorithm 5. Particularly,
assume that at time epoch $t$, the anchor local-load profile $\textsf{r}^{t}$
in Algorithm 5 matches the load induced by the anchor flow $r^{t}$ defined at
Line 3 of Algorithm 3 and also assume that the local-weight
$\tilde{\textsf{w}}^{t}$ used in Algorithm 5 and the weight $\tilde{w}^{t}$
defined at Line 3 of Algorithm 3 satisfy that $\sum_{e\in
p}\tilde{\textsf{w}}^{t,i}_{e}=\tilde{w}^{t}_{p}$ for any $p\in\mathcal{P}$
and $i\in\mathcal{N}$. We will prove that the local flow profile $x^{t}$
output from
$\textsc{Push\\-Pull\\-Match}(\textsf{r}^{t},\eta^{t}\tilde{\textsf{w}}^{t},\alpha^{t},\sum_{s=0}^{t}\alpha^{s})$
induces the same load profiles (and the same costs) as the flow profile
$f^{t}$ in Algorithm 3.

To do this, we observe that the following equality holds true for any
$i\in\mathcal{N}$ and $p\in\mathcal{P}^{i}$ and with the local flow profile
$z$ computed in the pushing-backward phase of
$\textsc{Push\\-Pull\\-Match}(\textsf{r}^{t},\eta^{t}\tilde{\textsf{w}}^{t},\alpha^{t},\sum_{s=0}^{t}\alpha^{s})$:

$M^{i}\prod_{e\in p}z^{i}_{v^{+}_{e},e}=M^{i}\prod_{e\in
p}\exp(\tilde{\textsf{w}}^{t,i}_{e})\frac{\texttt{BWscore}_{v^{-}_{e}}}{\texttt{BWscore}_{v^{+}_{e}}}=\frac{\exp\left(\sum_{e\in
p}\tilde{\textsf{w}}^{t,i}_{e}\right)}{\texttt{BWscore}_{O^{i}}}=\frac{\exp(\tilde{w}^{t}_{p})}{\sum_{q\in\mathcal{P}^{i}}\exp(\tilde{w}^{t}_{q})}=z^{t}_{p}.$
(31)

Here, the last equality comes directly from the update rule of $z^{t}$ in Line
3 of Algorithm 3.

Now, combine (31) with Proposition 5, the loads induced by $z$ of
$\textsc{Push\\-Pull\\-Match}(\textsf{r}^{t},\eta^{t}\tilde{\textsf{w}}^{t},\alpha^{t},\sum_{s=0}^{t}\alpha^{s})$
are precisely the loads induced by $z^{t}$ in Algorithm 3. Moreover, by
definition of local-load (cf. Section 5.1), we have that the term
${\ell}^{i}_{e}(z)$ computed at Line 4 of PushPullMatch are precisely these
loads.

Due to the arguments above and the assumption on the relation of
$\textsf{r}^{t}$ and $r^{t}$, the term ${\ell}^{i}_{e}(x)$ computed at Line 4
of PushPullMatch is precisely the load induced by the averaged flow $f^{t}$
computed at Line 3 of Algorithm 3. Finally, from the updating rule of the
local flow profile $x$ at Line 4 of PushPullMatch, we deduce that these loads
match precisely with the loads induced by $x$ (which is the output of
PushPullMatch).

As a conclusion, since we set $x^{t}$ as the output of
$\textsc{Push\\-Pull\\-Match}(\textsf{r}^{t},\eta^{t}\tilde{\textsf{w}}^{t},\alpha^{t},\sum_{s=0}^{t}\alpha^{s})$,
the loads (and hence, the costs) induced by $x^{t}$ are precisely that of the
flows $f^{t}$ recommended by Algorithm 3. Using a similar line of arguments,
we can also deduce that the local flow profile $\tilde{x}^{t}$ computed in
Algorithm 5 also matches the flow profile $\tilde{f}^{t}$ computed in
AdaWeight. For this reason, the assumptions we made on
$\tilde{\textsf{w}}^{t}$ and $\textsf{r}^{t}$ in the above proof actually hold
true.

To sum up, we have proved that the flow profiles defined in Theorem 4 have the
same loads as the local flow profiles recommended by AdaLight that, in turn,
coincide with the loads incurred by the recommended flows of AdaWeight. This
leads to the fact that AdaLight (Algorithm 5) inherits the convergence rates
of AdaWeight (Algorithm 3). This concludes the proof of Result (i) of Theorem
3.

##### Finally, we justify the per-iteration complexity of AdaLight

. It is trivial to see that when PushPullMatch is run, at any vertex $v$ and
O/D pair $i$, each of its phases only takes
$\operatorname{\mathcal{O}}(\max\\{\texttt{In}^{i}_{v},\texttt{Out}^{i}_{v}\\})$
rounds of computations where $\texttt{In}^{i}_{v}$ and $\texttt{Out}^{i}_{v}$
are the in-degree and out-degree of $v$ in $\mathcal{G}^{i}$. As mentioned
above, the learning rate update step can be done in
$\operatorname{\mathcal{O}}(\lvert\mathcal{V}\rvert)$ time. We conclude that
each iteration of AdaLight requires only an
$\operatorname{\mathcal{O}}\left(\lvert\mathcal{N}\rvert\lvert\mathcal{V}\rvert\lvert\mathcal{E}\rvert\right)$
number of computations at each vertex. ∎

##  6\. Numerical Experiments

In this section, we report the results of several numerical experiments that
we conducted to justify the theoretical convergence results of AdaWeight and
AdaLight. Particularly, in Section 6.1, we present the experiments on a toy
example to highlight the advantages of AdaWeight / AdaLight over several
benchmark algorithms. Then, in Section 6.2, we show the superiority in
performance of AdaLight in real-word networks. The experiments with real-world
data pose several additional computational challenges; we also discuss these
challenges and provide quick-fix solutions. The codes of our experiments are
available at dongquan-vu/Adaptive_Distributed_Routing.

###  6.1. Experiments on small-size networks

We first consider a toy-example of on a network with 4 vertices and 4 edges.
Specifically, the edges in this network are arranged to form 2 parallel paths
going from an origin vertex, namely $O$, to a destination vertex, namely $D$.
At each time epoch, a traffic demand (i.e., inflow) of size $M=10$ is sent
from $O$ to $D$.

In the experiments presented below, we run AdaWeight, ExpWeight and
AcceleWeight in $T_{\max}=10\textrm{e}5$ epochs. We record the values of the
BMW potential $c$ at outputs of each algorithm at each time
$T=1,2,\ldots,T_{\max}$. We then identify the value $c(f^{\ast})$ which is the
minimum among all $c$-value of flow profiles computed by these algorithms in
all iterations; thus, $f^{\ast}$ represents the equilibrium flow. In order to
track down the convergence properties of these algorithms, we will compute and
plot out the evolution (when $T=1,2,\ldots,T_{\max}$) of the _gaps_ (cf.
Equation 10) of the flows derived from these algorithms that are analyzed in
Theorem 1, Theorem 2 and Theorem 3 (we denote these gaps by
$\operatorname{Gap}_{\textrm{\acs{AEW}}}(T)$,
$\operatorname{Gap}_{\acs{EW}}(T)$, and
$\operatorname{Gap}_{\textrm{\acs{XLEW}}}(T)$ for short).

${T}$${\operatorname{Gap}(T)}$ (a) Log-log plot of $\operatorname{Gap}(T)$ in
static setting

${T}$${\operatorname{Gap}(T)}$ (b) Log-log plot of $\operatorname{Gap}(T)$ in
stochastic setting

Figure 1. Convergence speed of AdaWeight, ExpWeight, AcceleWeight in a
parallel network.

First, we consider a _static environment_. Particularly, the costs on edges
are determined via fixed linear cost functions such that when being routed
with the same load, one path has a higher cost than the other. We plot out the
evolution of $\operatorname{Gap}_{\textrm{\acs{AEW}}}(T)$,
$\operatorname{Gap}_{\acs{EW}}(T)$, and
$\operatorname{Gap}_{\textrm{\acs{XLEW}}}(T)$ in Fig. 1(a). We observe that in
this static setting, all these terms converge toward 0; in other words, all
three algorithms converge towards equilibria of the game. Importantly, we
observe that AdaWeight and AcceleWeight enjoy an accelerated convergence speed
that is much faster than that of ExpWeight; this coincides with theoretical
results in the previous sections. Note that if AcceleWeight is run with a
badly-tuned initial step-size $\gamma^{0}$ that is too small or too large
(e.g., when $\gamma^{0}=1$ or $\gamma^{0}=1\textrm{e}-8$), it requires a long
warming-up phase and might converge slowly. We also observe that in Fig. 1,
$\operatorname{Gap}_{\acs{EW}}(T)$ does not seem to fluctuate as much as
$\operatorname{Gap}_{\textrm{\acs{AEW}}}(T)$ and
$\operatorname{Gap}_{\textrm{\acs{XLEW}}}(T)$, this comes from the fact that
$\operatorname{Gap}_{\acs{EW}}(T)$ is computed from the time-averaged flow
profiles of outputs of ExpWeight. We recall that these time-averaged flows are
not the flows recommended by ExpWeight.

Second, we consider a _stochastic_ setting. Particularly, the costs on edges
are altered by adding noises generated randomly from a zero-mean normal
distribution. To have a better representation of these uncertainties, we run
each algorithm in 5 different instances of the noises’ layouts. We take the
averaged results among these instances and plot out the evolution of
$\operatorname{Gap}_{\textrm{\acs{AEW}}}(T)$,
$\operatorname{Gap}_{\acs{EW}}(T)$ and
$\operatorname{Gap}_{\textrm{\acs{XLEW}}}(T)$ in Fig. 1(b). In this setting,
we observe that while AdaWeight and ExpWeight converge toward equilibria,
AcceleWeight fails to do so. Moreover, in this stochastic environment, the
accelerated rate is no longer obtainable: AdaWeight and ExpWeight converge
with the same rate (of order $\operatorname{\mathcal{O}}(1/\sqrt{T})$). This
confirms the theoretical results of AdaWeight, ExpWeight and AcceleWeight.

###  6.2. Experiments on real-world datasets

In this section, we present several experiments using a real-world dataset
collected and provided in [dataset2021] (with free license). This dataset
contains the road networks of different cities in the world. The congestion
model in this dataset is assumed to follow the BPR cost functions whose
coefficients are estimated a priori. The purpose of our experiments is to
measure the equilibrium convergence of our proposed methods with the presence
of uncertainties and no knowledge on the cost functions.

Due to their per-iteration complexity, the naive implementations of ExpWeight,
AcceleWeight and AdaWeight in Algorithms 1, 2, 3 take an extremely long
running time when being run in these real-world networks. Therefore, to deal
with these large-scale networks, we consider the distributed model (presented
in Section 5.1) and run AdaLight in place of AdaWeight. As a benchmark, we
consider a variant of ExpWeight implemented with the classical weight-pushing
technique. These implementations enjoy a per-iteration complexity that is
polynomial in terms of the networks’ sizes.555In this section, to improve the
visibility of numerical results, we choose to only use ExpWeight as the
benchmark and do not report the performance of AcceleWeight that does not
guarantee an equilibrium convergence in stochastic settings.

Before presenting the obtained experimental results (in Section 6.2.2), we
first address a computational issue in using weight-pushing ideas in large-
scale networks. As far as we know, this issue has not been reported in any
previous works (prior to this work, weight-pushing is mostly analyzed
theoretically and real-world implementations have not been provided). We
formally address this computational issue and introduce a quick-fix solution
in Section 6.2.1 (readers who are eager to see the numerical experiments can
skip this section).

####  6.2.1. Computational issues of weight-pushing with large weights.

The PushPullMatch sub-algorithm that we constructed in Section 5.2 involves
the pushing-backward phase where a backward score (denoted by BWscore) is
assigned on each vertex in the network. In theory, this score can be computed
efficiently. In practice though, when the magnitude of costs values are large,
the input weights of PushPullMatch is proportional to the negative
accumulation of costs (that decreases quickly to $-\infty$) and hence, the
computation of BWscore involves a division of two infinitesimally small
numbers; this is often unsolvable by computers.

To resolve this issue, instead of keeping track of BWscore, we keep track of
its logarithm. Particularly, fix an O/D pair $i$, for each vertex $v$ and each
outgoing edge $e\in{\text{Out}^{i}_{v}}$, let us denote
$\texttt{logscore}_{v}\coloneqq\log(\texttt{BWscore}_{v})$ and
$\texttt{edgescore}_{e}=\texttt{logscore}_{v}+\textsf{w}^{i}_{e}$, then from
Line 4 of Algorithm 4, it can be computed by

$\displaystyle\texttt{logscore}_{v}$
$\displaystyle=\log\left[\sum_{e\in{\text{Out}^{i}_{v}}}\exp\left(\texttt{logscore}_{v}+\textsf{w}^{i}_{e}\right)\right]$
$\displaystyle=\max_{e^{\prime}\in{\text{Out}^{i}_{v}}}(\texttt{edgescore}_{e^{\prime}})+\log\left[\sum_{e\in{\text{Out}^{i}_{v}}}\exp\left(\texttt{edgescore}_{e}-\max_{e^{\prime}\in{\text{Out}^{i}_{v}}}(\texttt{edgescore}_{e^{\prime}})\right)\right].$
(32)

The expression in (32) allows computers to compute $\texttt{logscore}_{v}$
without any issues even when the magnitude of costs and $\textsf{w}^{i}_{e}$
is large. Finally, the local flow profile at Line 4 of the pushing-backward
phase of PushPullMatch can be computed as
$z^{i}_{v,e}={1}\big{/}{\exp(\texttt{logscore}_{v}-\texttt{edgescore}_{e})}$.

####  6.2.2. Experimental results.

We consider several instances in the dataset [dataset2021] representing the
urban traffic networks of different cities. The primitive parameters of these
networks are summarized in Table 2. More information on these networks are
given at https://github.com/bstabler/TransportationNetworks. Note that in
previous work, this dataset is used mostly for social-costs optimization; in
our knowledge, no work has used this dataset for equilibrium searching
problems.

Table 2. Networks’ sizes of several real-world instances Network name |  $\\#$ of vertices |  $\\#$ of edges |  $\\#$ of O/D pairs |  $\min$ demands |  $\max$ demands
---|---|---|---|---|---
SiouxFalls | 24 | 76 | 528 | 100 | 4400
Eastern-Massachusetts | 74 | 258 | 1113 | 0.5 | 957.7
Berlin-Friedrichshain | 224 | 523 | 506 | 1.05 | 107.53
Anaheim | 416 | 914 | 1406 | 1 | 2106.5

For our experiments, in each network, we run AdaLight and ExpWeight in
$T_{\max}=10\textrm{e}4$ time epochs.

##### Static environment.

First, we consider the static setting where in each time epoch, for each
network, the costs on edges are determined by a BPR function with coefficients
given by the dataset. The results are reported in Fig. 2, particularly as
follows:

  1. i.5pt)

In Fig. 2(a), Fig. 2(c), Fig. 2(e) and Fig. 2(g), we plot out the evolution of
$\operatorname{Gap}_{\acs{ALW}}(T)$ and $\operatorname{Gap}_{\acs{EW}}(T)$ in
different networks. Again, both algorithms converge toward equilibria.
Importantly, the results in these large-scale networks show the superiority of
AdaLight (and AdaWeight) in the convergence speed in comparison with
ExpWeight. We also make several side-comments as follows. First, the evolution
of $\operatorname{Gap}_{\textrm{\acs{ALW}}}(T)$ appears to involve less
fluctuations than that in the small-size network considered in Section 6.1.
This phenomenon is trivially explicable: in the toy-example with only two
paths, any small changes in the implemented flow profile might lead to a large
impact on the costs; in real-world instances, moving a mass of traffic from
one path to the others does not make so much differences in the total induced
costs accumulated from a large number of paths. Second, as we consider
networks with larger sizes, both AdaLight and ExpWeight require a longer
warming-up phase (where $\operatorname{Gap}(T)$ does not decrease
significantly) before picking up good recommendations and eventually converge.
Albeit the case, this warming-up phase is still reasonably small and it does
not restrict the applicability of our proposed algorithms.

  2. ii.5pt)

To validate the rate of convergence (in terms of $T$), we plot out the
evolution of the terms $T^{2}\cdot\operatorname{Gap}_{\acs{ALW}}(T)$ and
$T^{2}\cdot\operatorname{Gap}_{\acs{EW}}(T)$ in Fig. 2(b), Fig. 2(d), Fig.
2(f) and Fig. 2(h). In these plots,
$T^{2}\cdot\operatorname{Gap}_{\acs{ALW}}(T)$ approach a horizontal line; this
confirms that the convergence speed of AdaLight (also of AdaWeight) are
precisely in order $\operatorname{\mathcal{O}}(T^{2})$. On the other hand,
ExpWeight does not reach this convergence order (and hence,
$T^{2}\cdot\operatorname{Gap}_{\acs{EW}}(T)$ continue increasing).

##### Stochastic environment.

In this setting, the cost of edges are altered by adding random noises
generated from zero-mean normal distributions. For validation purposes, for
each network, we run AdaLight and ExpWeight in 5 instances (of noises’ layout)
and report the averaged results across these instances in Fig. 3.

  1. i.5pt)

In Fig. 3(a), Fig. 3(c), Fig. 3(e) and Fig. 3(g), we plot out the evolution of
$\operatorname{Gap}_{\acs{EW}}(T)$ and
$\operatorname{Gap}_{\textrm{\acs{ALW}}}(T)$. These terms tend to zero as $T$
increases; this confirms that ExpWeight and AdaLight converge toward
equilibria in this stochastic setting. Moreover, in this setting, the
convergence rate of AdaLight (also of AdaWeight) and ExpWeight are of the same
order.

  2. ii.5pt)

To justify the convergence speed of the considered algorithms, in Fig. 3(b),
Fig. 3(d), Fig. 3(f) and Fig. 3(h), we plot our the terms
$\sqrt{T}\cdot\operatorname{Gap}_{\textrm{\acs{ALW}}}(T)$ and
$\sqrt{T}\cdot\operatorname{Gap}_{\acs{EW}}(T)$. These terms approach
horizontal lines as $T$ increases. This reaffirms the fact that the speed of
convergence of ExpWeight and AdaWeight are
$\operatorname{\mathcal{O}}(1/\sqrt{T})$ in the stochastic regime. This result
is consistent with our theoretical results.

##### On the elapsed time of AdaLight.

We end this section with an interesting remark. Even in the largest network
instance in our experiments (Anaheim), it only takes AdaLight around _4
seconds to finish one round of learning iteration_ and to output a route
recommendation (the computations are conducted on a machine with the following
specs: Intel Core(TM) i7-9750H CPU 2.60GHz and 8GB RAM). For most of the
applications in urban traffic routing, the scale of time between fluctuations
of networks’ states is often much larger than this elapsed time (it might take
hours or even days for a significant change to happen). This highlights the
implementability and practicality of our proposed methods, even in networks
with much larger sizes.

${T}$${\operatorname{Gap}(T)}$ (a) SiouxFalls: Log-log plot of
$\operatorname{Gap}(T)$

${T}$$T^{2}\cdot{\operatorname{Gap}(T)}$ (b) SiouxFalls: Log-log plot of
$T^{2}\cdot\operatorname{Gap}(T)$

${T}$${\operatorname{Gap}(T)}$ (c) Eastern-Massachusetts: Log-log plot of
$\operatorname{Gap}(T)$

${T}$$T^{2}\cdot{\operatorname{Gap}(T)}$ (d) Eastern-Massachusetts: Log-log
plot of $T^{2}\cdot\operatorname{Gap}(T)$

${T}$${\operatorname{Gap}(T)}$ (e) Berlin-Friedrichshain: Log-log plot of
$\operatorname{Gap}(T)$

${T}$$T^{2}\cdot{\operatorname{Gap}(T)}$ (f) Berlin-Friedrichshain: Log-log
plot of $T^{2}\cdot\operatorname{Gap}(T)$

${T}$${\operatorname{Gap}(T)}$ (g) Anaheim: Log-log plot of
$\operatorname{Gap}(T)$

${T}$$T^{2}\cdot{\operatorname{Gap}(T)}$ (h) Anaheim: Log-log plot of
$T^{2}\cdot\operatorname{Gap}(T)$

Figure 2. Convergence speed of AdaLight (AdaWeight) and ExpWeight in static
environments.

${T}$${\operatorname{Gap}(T)}$ (a) SiouxFalls: Log-log plot of
$\operatorname{Gap}(T)$

${T}$$\sqrt{T}\cdot{\operatorname{Gap}(T)}$ (b) SiouxFalls: Log-log plot of
$\sqrt{T}\cdot\operatorname{Gap}(T)$

${T}$${\operatorname{Gap}(T)}$ (c) Eastern-Massachusetts: Log-log plot of
$\operatorname{Gap}(T)$

${T}$$\sqrt{T}\cdot{\operatorname{Gap}(T)}$ (d) Eastern-Massachusetts: Log-log
plot of $\sqrt{T}\cdot\operatorname{Gap}(T)$

${T}$${\operatorname{Gap}(T)}$ (e) Berlin-Friedrichshain: Log-log plot of
$\operatorname{Gap}(T)$

${T}$$\sqrt{T}\cdot{\operatorname{Gap}(T)}$ (f) Berlin-Friedrichshain: Log-log
plot of $\sqrt{T}\cdot\operatorname{Gap}(T)$

${T}$${\operatorname{Gap}(T)}$ (g) Anaheim: Log-log plot of
$\operatorname{Gap}(T)$

${T}$$\sqrt{T}\cdot{\operatorname{Gap}(T)}$ (h) Anaheim: Log-log plot of
$\sqrt{T}\cdot\operatorname{Gap}(T)$

Figure 3. Convergence speed of AdaLight (AdaWeight) and ExpWeight in
stochastic environments.

##  Appendix A Proof of Theorem 2

Our goal in this appendix is to prove the $\operatorname{\mathcal{O}}(\log
P\big{/}T^{2})$ equilibrium convergence rate of AcceleWeight in static
environments. In the following, we use the entropy regularizer $h$ as defined
in Section 4.2 and let $f^{t}$, $z^{t}$ and $\gamma^{t}$ be defined as per
Algorithm 2. With all this in hand, our proof of Theorem 2 will proceed in two
basic steps:

##### Step 1: Establish an energy function.

Building on the analysis of accelerated mirror descent algorithms [1, 22], we
will consider the energy fuction

$\Delta^{t}=\gamma^{t-1}\operatorname{Gap}(t)+\mathrm{KL}(f^{\ast}\nonscript\,\|\nonscript\,\mathopen{}z^{t})$
(A.1)

where $\operatorname{Gap}(t)=\Phi(f^{t})-\min\Phi$ and
$\mathrm{KL}(f^{\ast}\nonscript\,\|\nonscript\,\mathopen{}z^{t})$ denotes the
Kullback–Leibler divergence between $f^{\ast}$ and $z^{t}$.

Our goal in the sequel will be to prove that $\Delta^{t}$ is decreasing in
$t$. Indeed, from Proposition 2 and the strong convexity of $h$, we have

$\displaystyle\Phi(f^{t+1})$
$\displaystyle\leq\Phi(\bar{z}^{t})+\langle\nabla\Phi(\bar{z}^{t})\mathopen{},\mathopen{}f^{t+1}-\bar{z}^{t}\rangle+\frac{\beta}{2}\lVert
f^{t+1}-\bar{z}^{t}\rVert_{1}^{2}$
$\displaystyle=\Phi(\bar{z}^{t})+\langle\nabla\Phi(\bar{z}^{t})\mathopen{},\mathopen{}f^{t+1}-\bar{z}^{t}\rangle+\frac{\beta}{2}(1-\alpha^{t})^{2}\lVert
z^{t+1}-z^{t}\rVert_{1}^{2}$
$\displaystyle\leq\Phi(\bar{z}^{t})+\langle\nabla\Phi(\bar{z}^{t})\mathopen{},\mathopen{}f^{t+1}-\bar{z}^{t}\rangle+\beta\kappa(1-\alpha^{t})^{2}\mathrm{KL}(z^{t+1}\nonscript\,\|\nonscript\,\mathopen{}z^{t}).$
(A.2)

Moreover, from (A.2) and the convexity of $\Phi$, for all
$f^{\ast}\in\mathcal{F}^{\ast}\coloneqq\operatorname*{arg\,min}\Phi$, we have

$\displaystyle\operatorname{Gap}(t+1)-\alpha^{t}\operatorname{Gap}(t)$
$\displaystyle=\Phi(f^{t+1})-[\alpha^{t}\Phi(f^{t})+(1-\alpha^{t})\Phi(f^{\ast})]$
$\displaystyle\leq\Phi(f^{t+1})-\Phi(\bar{z}^{t})+\langle\nabla\Phi(\bar{z}^{t})\mathopen{},\mathopen{}\bar{z}^{t}-\alpha^{t}f^{t}-(1-\alpha^{t})f^{\ast}\rangle$
$\displaystyle\leq\langle\nabla\Phi(\bar{z}^{t})\mathopen{},\mathopen{}f^{t+1}-\bar{z}^{t}\rangle+\beta\kappa(1-\alpha^{t})^{2}\mathrm{KL}(z^{t+1}\nonscript\,\|\nonscript\,\mathopen{}z^{t})+\langle\nabla\Phi(\bar{z}^{t})\mathopen{},\mathopen{}\bar{z}^{t}-\alpha^{t}f^{t}-(1-\alpha^{t})f^{\ast}\rangle$
$\displaystyle=(1-\alpha^{t})\langle\nabla\Phi(\bar{z}^{t})\mathopen{},\mathopen{}z^{t+1}-f^{\ast}\rangle+\beta\kappa(1-\alpha^{t})^{2}\mathrm{KL}(z^{t+1}\nonscript\,\|\nonscript\,\mathopen{}z^{t}).$
(A.3)

To proceed, let ${W}^{t}=(1-\alpha^{t})\gamma^{t}\nabla\Phi(\bar{z}^{t})$.
Then, from Proposition 1 and the update structure of Algorithm 2, we have:

$\displaystyle\langle\nabla h(z^{t})-\nabla
h(z^{t+1})\mathopen{},\mathopen{}z^{t+1}-f^{\ast}\rangle$ $\displaystyle=$
$\displaystyle\sum_{i\in\mathcal{N}}\log\left(\frac{\sum_{q\in\mathcal{P}^{i}}z^{t}_{q}\exp(-{W}^{t}_{q})}{M^{i}}\right)\sum_{p\in\mathcal{P}^{i}}\big{(}z^{t+1}_{p}-f^{\ast}_{p}\big{)}+\sum_{i\in\mathcal{N}}\sum_{p\in\mathcal{P}^{i}}{W}^{t}_{p}\big{(}z^{t+1}_{p}-f^{\ast}_{p}\big{)}\hskip
120.00018pt$ $\displaystyle=$ $\displaystyle
0+\langle{W}^{t}\mathopen{},\mathopen{}z^{t+1}-f^{\ast}\rangle$ [since
$\sum_{p\in\mathcal{P}^{i}}z^{t+1}_{p}=\sum_{p\in\mathcal{P}^{i}}f^{\ast}_{p}=M^{i}$]
$\displaystyle=$
$\displaystyle(1-\alpha^{t})\gamma^{t}\langle\nabla\Phi(\bar{z}^{t})\mathopen{},\mathopen{}z^{t+1}-f^{\ast}\rangle.$
(A.4)

Thus, multiplying both sides of (A) by $\gamma^{t}$ and combining them with
(A.4), we obtain

$\displaystyle\gamma^{t}\operatorname{Gap}(t+1)-\alpha^{t}\gamma^{t}\operatorname{Gap}(t)$
$\displaystyle\leq\langle\nabla h(z^{t})-\nabla
h(z^{t+1})\mathopen{},\mathopen{}z^{t+1}-f^{\ast}\rangle+\gamma^{t}\kappa\beta(1-\alpha^{t})^{2}\mathrm{KL}(z^{t+1}\nonscript\,\|\nonscript\,\mathopen{}z^{t})$
$\displaystyle=\mathrm{KL}(f^{\ast}\nonscript\,\|\nonscript\,\mathopen{}z^{t})-\mathrm{KL}(f^{\ast}\nonscript\,\|\nonscript\,\mathopen{}z^{t+1})+\big{(}\gamma^{t}\kappa\beta(1-\alpha^{t})^{2}-1\big{)}\mathrm{KL}(z^{t+1}\nonscript\,\|\nonscript\,\mathopen{}z^{t}).$
(A.5)

Now, by the update rule of $\gamma^{t}$ in Algorithm 2 of Algorithm 2, the
choice of $\gamma^{0}$ in Algorithm 2 and the update rule of $\alpha^{t}$ in
Algorithm 2, we get

$\gamma^{t}\kappa\beta(1-\alpha^{t})^{2}=\frac{\gamma^{t}}{\gamma^{0}}\left(1-\frac{\gamma^{t-1}}{\gamma^{t}}\right)^{2}=\frac{\left(\gamma^{t-1}-\gamma^{t}\right)^{2}}{\gamma^{t}\gamma^{0}}=1.$
(A.6)

Therefore, the last term in (A.5) vanishes and we can rewrite (A.5) as

$\gamma^{t}\operatorname{Gap}(t+1)+\mathrm{KL}(f^{\ast}\nonscript\,\|\nonscript\,\mathopen{}z^{t+1})\leq\gamma^{t-1}\operatorname{Gap}(t)+\mathrm{KL}(f^{\ast}\nonscript\,\|\nonscript\,\mathopen{}z^{t}).$
(A.7)

This shows that $\Delta^{t}\leq\Delta^{t-1}\leq\dotsm\leq\Delta^{1}$ (by
convention, we set $\gamma^{-1}=0$), i.e., $\Delta^{t}$ is decreasing in $t$,
as claimed. ∎

##### Step 2: Upper-bounding the equilibrium gap.

By iterating (A.7), we readily obtain

$\gamma^{T-1}\operatorname{Gap}(T)\leq\Delta^{T}\leq\dotsm\leq\Delta^{1}=\mathrm{KL}(f^{\ast}\nonscript\,\|\nonscript\,\mathopen{}z^{1})$
(A.8)

and hence

$\operatorname{Gap}(T)\leq{\mathrm{KL}(f^{\ast}\nonscript\,\|\nonscript\,\mathopen{}z^{1})}\big{/}{\gamma^{T-1}}$
(A.9)

To finish the proof, we need to bound $\gamma^{T-1}$ and
$\mathrm{KL}(f^{\ast}\nonscript\,\|\nonscript\,\mathopen{}f^{1})$ from above.

We begin by noting that

$\sqrt{\kappa\beta\gamma^{t-1}}=\sqrt{\kappa\beta\gamma^{t}}\sqrt{1-\frac{1}{\sqrt{\gamma^{t}\kappa\beta}}}\leq\sqrt{\kappa\beta\gamma^{t}}\cdot\left(1-\frac{1}{2\sqrt{\kappa\beta\gamma^{t}}}\right)=\sqrt{\kappa\beta\gamma^{t}}-\frac{1}{2}$
(A.10)

so $\sqrt{\kappa\beta\gamma^{t}}\geq\sqrt{\kappa\beta\gamma^{t-1}}+1/2$.
Therefore, telescoping this last bound, we get

$\displaystyle\sqrt{\kappa\beta\gamma^{T-1}}\geq\sqrt{\kappa\beta\gamma^{0}}+\frac{T-1}{2}=\frac{T+1}{2}\geq\frac{T}{2}$

and hence

$\gamma^{T-1}>{T^{2}}\big{/}(4\kappa\beta)>0.$ (A.11)

Second, from the choice of $w^{0}$ and $\alpha^{0}$ in Algorithm 2 of
Algorithm 2, it follows that $f^{1}_{p}=z^{1}_{p}=M^{i}\big{/}P_{i}$ for all
$p\in\mathcal{P}^{i}$, $i\in\mathcal{N}$. We thus get $\langle\nabla
h(f^{1})\mathopen{},\mathopen{}f^{1}-f^{\ast}\rangle=0$ so
$\mathrm{KL}(f^{\ast}\nonscript\,\|\nonscript\,\mathopen{}f^{1})=h(f^{\ast})-h(f^{1})\leq\max
h-\min h$. Moreover, by a straightforward calculation, we get $\max
h=M_{\mathrm{tot}}\log(M_{\max})$ and $\min
h=M_{\mathrm{tot}}\log(P/M_{\mathrm{tot}})$, so

$\displaystyle\mathrm{KL}(f^{\ast}\nonscript\,\|\nonscript\,\mathopen{}f^{1})\leq
M_{\mathrm{tot}}\log(PM_{\max}/M_{\mathrm{tot}}).$ (A.12)

Thus, combining Equations A.12, A.11 and A.8 and recalling that
$\kappa\coloneqq M_{\max}N$ (i.e., the strongly-convexity constant of $h$) and
the fact that $M_{\mathrm{tot}}=\sum_{i\in\mathcal{N}}M^{i}\leq NM_{\max}$ we
finally obtain

$\operatorname{Gap}(T)\leq\frac{4\beta
NM_{\max}M_{\mathrm{tot}}\log\left(PM_{\max}\big{/}M_{\mathrm{tot}}\right)}{\left(T-1\right)^{2}}\leq\frac{4\beta
N^{2}M_{\max}^{2}\log\left(PM_{\max}\big{/}M_{\mathrm{tot}}\right)}{\left(T-1\right)^{2}}$
(A.13)

and our proof of Theorem 2 is complete. ∎

##  Appendix B Proof of Equation 15

Let us denote $R^{t}:=\sum_{s=1}^{t}\alpha^{s}$. By Line 3 of Algorithm 3, for
any $t$, we have:
$z^{t}=\frac{R^{t}}{\alpha^{t}}f^{t}-\frac{R^{t-1}}{\alpha^{t}}f^{t-1}$. As a
consequence,

$\displaystyle\sum_{t=1}^{T}{\alpha^{t}}\langle
z^{t}-f^{\ast},\nabla\Phi(f^{t})\rangle=$
$\displaystyle\sum_{t=1}^{T}{\alpha^{t}}\left\langle\frac{R^{t}}{\alpha^{t}}f^{t}-\frac{R^{t-1}}{\alpha^{t}}f^{t-1}-f^{\ast},\nabla\Phi(f^{t})\right\rangle$
$\displaystyle=$ $\displaystyle\sum_{t=1}^{T}\left[R^{t-1}\left\langle
f^{t}-f^{t-1},\nabla\Phi(f^{t})\right\rangle+\alpha^{t}\left\langle
f^{t}-f^{\ast},\nabla\Phi(f^{t})\right\rangle\right]$ $\displaystyle\geq$
$\displaystyle\sum_{t=1}^{T}R^{t-1}\left[\Phi(f^{t})-\Phi(f^{t-1})\right]+\sum_{t=1}^{T}\alpha^{t}\left[\Phi(f^{t})-\Phi(f^{\ast})\right]$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{T}\alpha^{t}\left[\Phi(f^{T})-\Phi(f^{\ast})\right].$
(B.1)

Here, the last equality is achieved via telescopic sum. Now, when we choose
$\alpha^{t}=t$ (as indicated in Theorem 3), we notice that
$R^{t}>\frac{T^{2}}{2}$. Divide two sides of (B.1) by $R^{t}$ and taking
expectations, we obtain precisely (15). ∎

## References

  * Allen-Zhu and Orecchia [2017] Allen-Zhu Z, Orecchia L (2017) Linear coupling: An ultimate unification of gradient and mirror descent. _ITCS ’17: Proceedings of the 8th Conference on Innovations in Theoretical Computer Science_.
  * Antonakopoulos et al. [2021] Antonakopoulos K, Belmega EV, Mertikopoulos P (2021) Adaptive extra-gradient methods for min-max optimization and games. _ICLR ’21: Proceedings of the 2021 International Conference on Learning Representations_.
  * Auer et al. [1995] Auer P, Cesa-Bianchi N, Freund Y, Schapire RE (1995) Gambling in a rigged casino: The adversarial multi-armed bandit problem. _Proceedings of the 36th Annual Symposium on Foundations of Computer Science_.
  * Awerbuch [1989] Awerbuch B (1989) Randomized distributed shortest paths algorithms. _Proceedings of the twenty-first Annual ACM symposium on Theory of Computing_ , 490–500.
  * Beck and Teboulle [2003] Beck A, Teboulle M (2003) Mirror descent and nonlinear projected subgradient methods for convex optimization. _Operations Research Letters_ 31(3):167–175.
  * Beckmann et al. [1956] Beckmann M, McGuire CB, Winsten C (1956) _Studies in the Economics of Transportation_ (Yale University Press).
  * Bellman [1958] Bellman R (1958) On a routing problem. _Quarterly of applied mathematics_ 16(1):87–90.
  * Blum et al. [2006] Blum A, Even-Dar E, Ligett K (2006) Routing without regret: on convergence to Nash equilibria of regret-minimizing in routing games. _PODC ’06: Proceedings of the 25th annual ACM SIGACT-SIGOPS symposium on Principles of Distributed Computing_ , 45–52.
  * Bubeck [2015] Bubeck S (2015) Convex optimization: Algorithms and complexity. _Foundations and Trends in Machine Learning_ 8(3-4):231–358.
  * Cabannes et al. [2019] Cabannes T, Sangiovanni M, Keimer A, Bayen AM (2019) Regrets in routing networks: Measuring the impact of routing apps in traffic. _ACM Transactions on Spatial Algorithms and Systems (TSAS)_ 5(2):1–19.
  * Chen and Teboulle [1993] Chen G, Teboulle M (1993) Convergence analysis of a proximal-like minimization algorithm using Bregman functions. _SIAM Journal on Optimization_ 3(3):538–543.
  * Colini-Baldeschi et al. [2020] Colini-Baldeschi R, Cominetti R, Mertikopoulos P, Scarsini M (2020) When is selfish routing bad? The price of anarchy in light and heavy traffic. _Operations Research_ 68(2):411–434.
  * Cutkosky [2019] Cutkosky A (2019) Anytime online-to-batch, optimism and acceleration. _International Conference on Machine Learning_ , 1446–1454 (PMLR).
  * Dafermos and Sparrow [1969] Dafermos SC, Sparrow FT (1969) The traffic assignment problem for a general network. _Journal of Research of the National Bureau of Standards_ 73B(2):91–118.
  * Elkin [2020] Elkin M (2020) Distributed exact shortest paths in sublinear time. _Journal of the ACM (JACM)_ 67(3):1–36.
  * Fischer and Vöcking [2004] Fischer S, Vöcking B (2004) On the evolution of selfish routing. _Proceedings of the 12th European Symposium on Algorithms_.
  * Ford Jr [1956] Ford Jr LR (1956) Network flow theory. Technical report, Rand Corp Santa Monica Ca.
  * György et al. [2007] György A, Linder T, Lugosi G, Ottucsák G (2007) The online shortest path problem under partial monitoring. _Journal of Machine Learning Research_ 8:2369–2403.
  * Huang et al. [2017] Huang CC, Nanongkai D, Saranurak T (2017) Distributed exact weighted all-pairs shortest paths in õ (n^$\\{$5/4$\\}$) rounds. _2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS)_ , 168–179 (IEEE).
  * Kavis et al. [2019] Kavis A, Levy KY, Bach F, Cevher V (2019) UnixGrad: A universal, adaptive algorithm with optimal guarantees for constrained optimization. _NeurIPS ’19: Proceedings of the 33rd International Conference on Neural Information Processing Systems_.
  * Krichene et al. [2015a] Krichene S, Krichene W, Dong R, Bayen A (2015a) Convergence of heterogeneous distributed learning in stochastic routing games. _Allerton ’15: Proceedings of the 52nd Annual Allerton Conference on Communication, Control, and Computing_.
  * Krichene et al. [2015b] Krichene W, Bayen A, Bartlett PL (2015b) Accelerated mirror descent in continuous and discrete time. _NIPS ’15: Proceedings of the 29th International Conference on Neural Information Processing Systems_.
  * Krichene et al. [2014] Krichene W, Drighès B, Bayen A (2014) On the convergence of no-regret learning in selfish routing. _ICML ’14: Proceedings of the 31st International Conference on Machine Learning_.
  * Krichene et al. [2015c] Krichene W, Drighès B, Bayen AM (2015c) Online learning of Nash equilibria in congestion games. _SIAM Journal on Control and Optimization_ 53(2):1056–1081.
  * Lan [2012] Lan G (2012) An optimal method for stochastic composite optimization. _Mathematical Programming_ 133(1-2):365–397.
  * LeBlanc et al. [1975] LeBlanc LJ, Morlok EK, Pierskalla WP (1975) An efficient approach to solving the road network equilibrium traffic assignment problem. _Transportation Research_ 9:309–318.
  * Levy et al. [2018] Levy KY, Yurtsever A, Cevher V (2018) Online adaptive methods, universality and acceleration. _NeurIPS ’18: Proceedings of the 32nd International Conference of Neural Information Processing Systems_.
  * Littlestone and Warmuth [1994] Littlestone N, Warmuth MK (1994) The weighted majority algorithm. _Information and Computation_ 108(2):212–261.
  * McMahan and Streeter [2010] McMahan HB, Streeter M (2010) Adaptive bound optimization for online convex optimization. _COLT ’10: Proceedings of the 23rd Annual Conference on Learning Theory_.
  * Nesterov [1983] Nesterov Y (1983) A method for unconstrained convex minimization problem with the rate of convergence $\mathcal{O}(1/k^{2})$. _Proceedings of the USSR Academy of Sciences_ 269(543-547).
  * Nesterov [2004] Nesterov Y (2004) _Introductory Lectures on Convex Optimization: A Basic Course_. Number 87 in Applied Optimization (Kluwer Academic Publishers).
  * Nesterov [2007] Nesterov Y (2007) Dual extrapolation and its applications to solving variational inequalities and related problems. _Mathematical Programming_ 109(2):319–344.
  * Nesterov [2009] Nesterov Y (2009) Primal-dual subgradient methods for convex problems. _Mathematical Programming_ 120(1):221–259.