Mario Bravo
email: mario.bravo.g@usach.cl

David Stuart Leslie
email: d.leslie@lancaster.ac.uk

AND Panayotis Mertikopoulos
email: panayotis.mertikopoulos@imag.fr.2010

BANDIT LEARNING IN CONCAVE N -PERSON GAMES

Keywords: Mathematics Subject Classification. Primary 91A10, 91A26, secondary 68Q32, 68T02 Bandit feedback, concave games, Nash equilibrium, mirror descent

Introduction

The bane of decision-making in an unknown environment is regret: noone wants to realize in hindsight that the decision policy they employed was strictly inferior to a plain policy prescribing the same action throughout. For obvious reasons, this issue becomes considerably more intricate when the decision-maker is subject to situational uncertainty and the "fog of war": when the only information at the optimizer's disposal is the reward obtained from a given action (the so-called "bandit" framework), is it even possible to design a no-regret policy? Especially in the context of online convex optimization (repeated decision problems with continuous action sets and convex costs), this problem becomes even more challenging because the decision-maker typically needs to infer gradient information from the observation of a single scalar. Nonetheless, despite this extra degree of difficulty, this question has been shown to admit a positive answer: regret minimization is possible, even with bandit feedback [START_REF] Flaxman | Online convex optimization in the bandit setting: gradient descent without a gradient[END_REF][START_REF] Kleinberg | Nearly tight bounds for the continuum-armed bandit problem[END_REF].

In this paper, we consider a multi-agent extension of this framework where, at each stage n = 1, 2, . . . , of a repeated decision process, the reward of an agent is determined by the actions of all agents via a fixed mechanism: a non-cooperative N -person game. In general, the agents -or players -might be completely oblivious to this mechanism, perhaps even ignoring its existence: for instance, when choosing how much to bid for a good in an online auction, an agent is typically unaware of who the other bidders are, what are their specific valuations, etc. Hence, lacking any knowledge about the game, it is only natural to assume that agents will at least seek to achieve a minimal worst-case guarantee and minimize their regret. As a result, a fundamental question that arises is a) whether the agents' sequence of actions stabilizes to a rationally admissible state under no-regret learning; and b) if it does, whether convergence is affected by the information available to the agents.

Related work.

In finite games, no-regret learning guarantees that the players' timeaveraged, empirical frequency of play converges to the game's set of coarse correlated equilibria (CCE), and the rate of this convergence is O(1/n) for (λ, µ)-smooth games [START_REF] Foster | Learning in games: Robustness of fast convergence[END_REF]Syrgkanis et al., 2015). In general however, this set might contain highly subpar, rationally inadmissible strategies: for instance, Viossat and Zapechelnyuk (2013) provide examples of CCE that assign positive selection probability only to strictly dominated strategies. In the class of potential games, [START_REF] Cohen | Learning with bandit feedback in potential games[END_REF] recently showed that the actual sequence of play (i.e., the sequence of actions that determine the agents' rewards at each stage) converges under no-regret learning, even with bandit feedback. Outside this class however, the players' chosen actions may cycle in perpetuity, even in simple, two-player zero-sum games with full information (Mertikopoulos et al., 2018a,b); in fact, depending on the parameters of the players' learning process, agents could even exhibit a fully unpredictable, aperiodic and chaotic behavior [START_REF] Palaiopanos | Multiplicative weights update with constant step-size in congestion games: Convergence, limit cycles and chaos[END_REF]. As such, without further assumptions in place, no-regret learning in a multi-agent setting does not necessarily imply convergence to a unilaterally stable, equilibrium state.

In the broader context of games with continuous action sets (the focal point of this paper), the long-run behavior of no-regret learning is significantly more challenging to analyze. In the case of mixed-strategy learning, [START_REF] Perkins | Stochastic fictitious play with continuous action sets[END_REF] and [START_REF] Perkins | Mixed-strategy learning with continuous action sets[END_REF] showed that mixed-stratgy learning based on stochastic fictitious play converges to an ε-perturbed Nash equilibrium in potential games (but may lead to as much as O(εn) regret in the process). More relevant for our purposes is the analysis of [START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF] who showed that the time-averaged sequence of play induced by a no-regret dual averaging (DA) process with noisy gradient feedback converges to Nash equilibrium in monotone games (a class which, in turn, contains all concave potential games).

The closest antecedent to our approach is the recent work of [START_REF] Mertikopoulos | Learning in games with continuous action sets and unknown payoff functions[END_REF] who showed that the actual sequence of play generated by dual averaging converges to Nash equilibrium in the class of variationally stable games (which includes all monotone games). To do so, the authors first showed that a naturally associated continuous-time dynamical system converges, and then used the so-called asymptotic pseudotrajectory (APT) framework of [START_REF] Benaïm | Dynamics of stochastic approximation algorithms[END_REF] to translate this result to discrete time. Similar APT techniques were also used in a very recent preprint by [START_REF] Bervoets | Learning with minimal information in continuous games[END_REF] to establish the convergence of a payoffbased learning algorithm in two classes of one-dimensional concave games: games with strategic complements, and ordinal potential games with isolated equilibria. The algorithm of [START_REF] Bervoets | Learning with minimal information in continuous games[END_REF] can be seen as a special case of mirror descent coupled with a two-point gradient estimation process, suggesting several interesting links with our paper.

Our contributions. In this paper, we drop all feedback assumptions and we focus on the bandit framework where the only information at the players' disposal is the payoffs they receive at each stage. As we discussed above, this lack of information complicates matters considerably because players must now estimate their payoff gradients from their observed rewards. What makes matters even worse is that an agent may introduce a significant bias in the (concurrent) estimation process of another, so traditional, multiple-point estimation techniques for derivative-free optimization cannot be applied (at least, not without significant communication overhead between players).

To do away with player coordination requirements, we focus on learning processes which could be sensibly deployed in a single-agent setting and we show that, in monotone games, the sequence of play induced by a wide class of no-regret learning policies converges to Nash equilibrium with probability 1. Furthermore, by specializing to the class of strongly monotone games, we show that the rate of convergence is O(n -1/3 ), i.e., it is nearly optimal with respect to the attainable O(n -1/2 ) rate for bandit, single-agent stochastic optimization with strongly convex and smooth objectives [START_REF] Agarwal | Optimal algorithms for online convex optimization with multi-point bandit feedback[END_REF][START_REF] Shamir | On the complexity of bandit and derivative-free stochastic convex optimization[END_REF].

We are not aware of a similar Nash equilibrium convergence result for concave games with general convex action spaces and bandit feedback: the analysis of [START_REF] Mertikopoulos | Learning in games with continuous action sets and unknown payoff functions[END_REF] requires first-order feedback, while the analysis of [START_REF] Bervoets | Learning with minimal information in continuous games[END_REF] only applies to one-dimensional games. We find this outcome particularly appealing for practical applications of game theory (e.g., in network routing) because it shows that in a wide class of (possibly very complicated) nonlinear games, the Nash equilibrium prediction does not require full rationality, common knowledge of rationality, flawless execution, or even the knowledge that a game is being played: a commonly-used, individual no-regret algorithm suffices.

Problem setup and preliminaries

2.1. Concave games. Throughout this paper, we will focus on games with a finite number of players i ∈ N = {1, . . . , N } and continuous action sets. During play, every player i ∈ N selects an action x i from a compact convex subset X i of a d i -dimensional normed space V i ; subsequently, based on each player's individual objective and the action profile x = (x i ; x -i ) ≡ (x 1 , . . . , x N ) of all players' actions, every player receives a reward, and the process repeats. In more detail, writing X ≡ i X i for the game's action space, we assume that each player's reward is determined by an associated payoff (or utility) function u i : X → R. Since players are not assumed to "know the game" (or even that they are involved in one) these payoff functions might be a priori unknown, especially with respect to the dependence on the actions of other players. Our only structural assumption for u i will be that

u i (x i ; x -i ) is concave in x i for all x -i ∈ X -i ≡ j =i X j , i ∈ N .
With all this in hand, a concave game will be a tuple G ≡ G(N , X , u) with players, action spaces and payoffs defined as above. Below, we briefly discuss some examples thereof:

Example 2.1 (Cournot competition). In the standard Cournot oligopoly model, there is a finite set of firms indexed by i = 1, . . . , N , each supplying the market with a quantity x i ∈ [0, C i ] of some good (or service), up to the firm's production capacity C i . By the law of supply and demand, the good is priced as a decreasing function P (x tot ) of the total amount x tot = N i=1 x i supplied to the market, typically following a linear model of the form P (x tot ) = a-bx tot for positive constants a, b > 0. The utility of firm i is then given by

u i (x i ; x -i ) = x i P (x tot ) -c i x i , (2.1)
i.e., it comprises the total revenue from producing x i units of the good in question minus the associated production cost (in the above, c i > 0 represents the marginal production cost of firm i).

Example 2.2 (Resource allocation auctions). Consider a service provider with a number of splittable resources s ∈ S = {1, . . . , S} (bandwidth, server time, GPU cores, etc.). These resources can be leased to a set of N bidders (players) who can place monetary bids x is ≥ 0 for the utilization of each resource s ∈ S up to each player's total budget b i , i.e., s∈S x is ≤ b i . Once all bids are in, resources are allocated proportionally to each player's bid, i.e., the i-th player gets ρ is = (q s x is ) (c s + j∈N x js ) units of the s-th resource (where q s denotes the available units of said resource and c s ≥ 0 is the "entry barrier" for bidding on it). A simple model for the utility of player i is then given by

u i (x i ; x -i ) = s∈S [g i ρ is -x is ],
(2.2) with g i denoting the marginal gain of player i from acquiring a unit slice of resources.

2.2. Nash equilibrium and monotone games. The most widely used solution concept for non-cooperative games is that of a Nash equilibrium (NE), defined here as any action profile x * ∈ X that is resilient to unilateral deviations, viz.

u i (x * i ; x * -i ) ≥ u i (x i ; x * -i ) for all x i ∈ X i , i ∈ N . (NE)
By the classical existence theorem of [START_REF] Debreu | A social equilibrium existence theorem[END_REF], every concave game admits a Nash equilibrium. Moreover, thanks to the individual concavity of the game's payoff functions, Nash equilibria can also be characterized via the first-order optimality condition v i (x * ), x i -x * i ≤ 0 for all x i ∈ X i , (2.3) where v i (x) denotes the individual payoff gradient of the i-th player, i.e.,

v i (x) = ∇ i u i (x i ; x -i ),
(2.4) with ∇ i denoting differentiation with respect to x i . 1 In terms of regularity, it will be convenient to assume that each v i is Lipschitz continuous; to streamline our presentation, this will be our standing assumption in what follows.

Starting with the seminal work of [START_REF] Rosen | Existence and uniqueness of equilibrium points for concave N -person games[END_REF], much of the literature on continuous games and their applications has focused on games that satisfy a condition known as diagonal strict concavity (DSC). In its simplest form, this condition posits that there exist positive constants λ i > 0 such that

i∈N λ i v i (x ) -v i (x), x i -x i < 0 for all x, x ∈ X , x = x .
(DSC)

Owing to the formal similarity between (DSC) and the various operator monotonicity conditions in optimization (see e.g., [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF], games that satisfy (DSC) are commonly referred to as (strictly) monotone. As was shown by Rosen (1965, Theorem 2), monotone games admit a unique Nash equilibrium x * ∈ X ,

1
We adopt here the standard convention of treating v i (x) as an element of the dual space

Y i ≡ V * i of V i , with y i , x i denoting the duality pairing between y i ∈ Y i and x i ∈ X i ⊆ V i .
which, in view of (DSC) and (NE), is also the unique solution of the (weighted) variational inequality

i∈N λ i v i (x), x i -x * i < 0 for all x = x * . (VI)
This property of Nash equilibria of monotone games will play a crucial role in our analysis and we will use it freely in the rest of our paper.

In terms of applications, monotonicity gives rise to a very rich class of games. As we show in the paper's supplement, Examples 2.1 and 2.2 both satisfy diagonal strict concavity (with a nontrivial choice of weights for the latter), as do atomic splittable congestion games in networks with parallel links [START_REF] Orda | Competitive routing in multi-user communication networks[END_REF]Sorin and Wan, 2016), multi-user covariance matrix optimization problems in multiple-input and multiple-output (MIMO) systems [START_REF] Mertikopoulos | Distributed stochastic optimization via matrix exponential learning[END_REF], and many other problems where online decision-making is the norm. Namely, the class of monotone games contains all strictly convex-concave zero-sum games and all games that admit a (strictly) concave potential, i.e., a function f :

X → R such that v i (x) = ∇ i f (x) for all x ∈ X , i ∈ N .
In view of all this (and unless explicitly stated otherwise), we will focus throughout on monotone games; for completeness, we also include in the supplement a straightforward second-order test for monotonicity.

Regularized no-regret learning

We now turn to the learning methods that players could employ to increase their individual rewards in an online manner. Building on Zinkevich's ( 2003) online gradient descent policy, the most widely used algorithmic schemes for no-regret learning in the context of online convex optimization invariably revolve around the idea of regularization. To name but the most well-known paradigms, "following the regularized leader" (FTRL) explicitly relies on best-responding to a regularized aggregate of the reward functions revealed up to a given stage, while online mirror descent (OMD) and its variants use a linear surrogate thereof. All these no-regret policies fall under the general umbrella of "regularized learning" and their origins can be traced back to the seminal mirror descent (MD) algorithm of [START_REF] Nemirovski | Problem Complexity and Method Efficiency in Optimization[END_REF]. 2 The basic idea of mirror descent is to generate a new feasible point x + by taking a so-called "mirror step" from a starting point x along the direction of an "approximate gradient" vector y (which we treat here as an element of the dual space

Y ≡ i Y i of V ≡ i V i ). 3 To do so, let h i : X i → R be a continuous and K i -strongly convex distance-generating (or regularizer ) function, i.e., h i (tx i + (1 -t)x i ) ≤ th i (x i ) + (1 -t)h i (x i ) -1 2 K i t(1 -t) x i -x i 2 , (3.1)
for all x i , x i ∈ X i and all t ∈ [0, 1]. In terms of smoothness (and in a slight abuse of notation) we also assume that the subdifferential of h i admits a continuous selection, i.e., a continuous function

∇h i : dom ∂h i → Y i such that ∇h i (x i ) ∈ ∂h i (x i ) for all
2 In a utility maximization setting, mirror descent should be called mirror ascent because players seek to maximize their rewards (as opposed to minimizing their losses). Nonetheless, we keep the term "descent" throughout because, despite the role reversal, it is the standard name associated with the method.

3 For concreteness (and in a slight abuse of notation), we assume in what follows that V is equipped with the product norm x 2 = i x i 2 and Y with the dual norm y * = max{ y, x : x ≤ 1}.

x i ∈ dom ∂h i . 4 Then, letting h(x) = i h i (x i ) for x ∈ X (so h is strongly convex with modulus K = min i K i ), we get a pseudo-distance on X via the relation

D(p, x) = h(p) -h(x) -∇h(x), p -x , (3.2)
for all p ∈ X , x ∈ dom ∂h. This pseudo-distance is known as the Bregman divergence and we have D(p, x) ≥ 0 with equality if and only if x = p; on the other hand, D may fail to be symmetric and/or satisfy the triangle inequality so, in general, it is not a bona fide distance function on X . Nevertheless, we also have D(p, x) ≥ 1 2 K x -p 2 (see the paper's supplement), so the convergence of a sequence X n to p can be checked by showing that D(p, X n ) → 0. For technical reasons, it will be convenient to also assume the converse, i.e., that D(p, X n ) → 0 when X n → p. This condition is known in the literature as "Bregman reciprocity" [START_REF] Chen | Convergence analysis of a proximal-like minimization algorithm using Bregman functions[END_REF], and it will be our blanket assumption in what follows (note that it is trivially satisfied by Examples 3.1 and 3.2 below). Now, as with true Euclidean distances, D(p, x) induces a prox-mapping given by

P x (y) = arg min x ∈X { y, x -x + D(x , x)} (3.3)
for all x ∈ dom ∂h and all y ∈ Y. Just like its Euclidean counterpart below, the prox-mapping (3.3) starts with a point x ∈ dom ∂h and steps along the dual (gradient-like) vector y ∈ Y to produce a new feasible point x + = P x (y). Standard examples of this process are:

Example 3.1 (Euclidean projections). Let h(x) = 1 2 x 2 2 denote the Euclidean squared norm. Then, the induced prox-mapping is

P x (y) = Π(x + y), (3.4)
with Π(x) = arg min x ∈X x -x 2 denoting the standard Euclidean projection onto X . Hence, the update rule x + = P x (y) boils down to a "vanilla", Euclidean projection step along y.

Example 3.2 (Entropic regularization and multiplicative weights). Suppressing the player index for simplicity, let X be a d-dimensional simplex and consider the entropic regularizer h(x) = d j=1 x j log x j . The induced pseudo-distance is the so-called Kullback-Leibler (KL) divergence D KL (p, x) = d j=1 p j log(p j /x j ), which gives rise to the prox-mapping

P x (y) = (x j exp(y j )) d j=1 d j=1 x j exp(y j ) (3.5)
for all x ∈ X • , y ∈ Y. The update rule x + = P x (y) is widely known as the multiplicative weights (MW) algorithm and plays a central role for learning in multiarmed bandit problems and finite games [START_REF] Arora | The multiplicative weights update method: A meta-algorithm and applications[END_REF][START_REF] Auer | Gambling in a rigged casino: The adversarial multi-armed bandit problem[END_REF][START_REF] Freund | Adaptive game playing using multiplicative weights[END_REF]. 4 Recall here that the subdifferential of

h i at x i ∈ X i is defined as ∂h i (x i ) ≡ {y i ∈ Y i : h i (x i ) ≥ h i (x i )+ y i , x i -x i for all x i ∈ V i }, with the standard convention that h i (x i ) = +∞ if x i ∈ V i \X i .
By standard results, the domain of subdifferentiability ∂h i ≡ {x i ∈ X i :

∂h i = ∅} of h i satisfies X • i ⊆ dom ∂h i ⊆ X i .
With all this in hand, the multi-agent mirror descent (MD) algorithm is given by the recursion

X n+1 = P Xn (γ n vn ), (MD)
where γ n is a variable step-size sequence and vn = (v i,n ) i∈N is a generic feedback sequence of estimated gradients. In the next section, we detail how this sequence is generated with first-or zeroth-order (bandit) feedback.

4. First-order vs. bandit feedback 4.1. First-order feedback. A common assumption in the literature is that players are able to obtain gradient information by querying a first-order oracle [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF]. i.e., a "black-box" feedback mechanism that outputs an estimate vi of the individual payoff gradient v i (x) of the i-th player at the current action profile x = (x i ; x -i ) ∈ X . This estimate could be either perfect, giving vi = v i (x) for all i ∈ N , or imperfect, returning noisy information of the form vi = v i (x) + U i where U i denotes the oracle's error (random, systematic, or otherwise).

Having access to a perfect oracle is usually a tall order, either because payoff gradients are difficult to compute directly (especially without global knowledge), because they involve an expectation over a possibly unknown probability law, or for any other number of reasons. It is therefore more common to assume that each player has access to a stochastic oracle which, when called against a sequence of actions X n ∈ X , produces a sequence of gradient estimates vn = (v i,n ) i∈N that satisfies the following statistical assumptions:

a) Unbiasedness: E[v n | F n ] = v(X n ). b) Finite mean square: E[ vn 2 * | F n ] ≤ V 2 for some finite V ≥ 0.
(4.1) In terms of measurability, the expectation in (4.1) is conditioned on the history F n of X n up to stage n; in particular, since vn is generated randomly from X n , it is not F n -measurable (and hence not adapted). To make this more transparent, we will write vn = v(X n ) + U n+1 where U n is an adapted martingale difference sequence with

E[ U n+1 2 * | F n ] ≤ σ 2 for some finite σ ≥ 0.
4.2. Bandit feedback. Now, if players don't have access to a first-order oracle -the so-called bandit or payoff-based framework -they will need to derive an individual gradient estimate from the only information at their disposal: the actual payoffs they receive at each stage. When a function can be queried at multiple points (as few as two in practice), there are efficient ways to estimate its gradient via directional sampling techniques as in [START_REF] Agarwal | Optimal algorithms for online convex optimization with multi-point bandit feedback[END_REF]. In a game-theoretic setting however, multiple-point estimation techniques do not apply because, in general, a player's payoff function depends on the actions of all players. Thus, when a player attempts to get a second query of their payoff function, this function may have already changed due to the query of another player -i.e., instead of sampling u i (•; x -i ), the i-th player would be sampling u i (•; x -i ) for some x -i = x -i .

Following Spall (1997) and [START_REF] Flaxman | Online convex optimization in the bandit setting: gradient descent without a gradient[END_REF], we posit instead that players rely on a simultaneous perturbation stochastic approximation (SPSA) approach that allows them to estimate their individual payoff gradients v i based off a single function evaluation. In detail, the key steps of this one-shot estimation process for each player i ∈ N are:

(0) Fix a query radius δ > 0. 5 (1) Pick a pivot point x i ∈ X i where player i seeks to estimate their payoff gradient.

(2) Draw a vector z i from the unit sphere

S i ≡ S di of V i ≡ R di and play xi = x i + δz i . 6 (3) Receive ûi = u i (x i ; x-i ) and set vi = d i δ ûi z i . (4.2)
By adapting a standard argument based on Stokes' theorem (detailed in the supplement), it can be shown that vi is an unbiased estimator of the individual gradient of the δ-smoothed payoff function

u δ i (x) = 1 vol(δB i ) j =i vol(δS j ) δBi j =i δSj u i (x i +w i ; x -i +z -i ) dz 1 • • • dw i • • • dz N (4.3) with B i ≡ B di denoting the unit ball of V i . The Lipschitz continuity of v i guarantees that ∇ i u i -∇ i u δ i ∞ = O(δ)
, so this estimate becomes more and more accurate as δ → 0 + . On the other hand, the second moment of vi grows as O(1/δ 2 ), implying in turn that the variability of vi grows unbounded as δ → 0 + . This manifestation of the bias-variance dilemma plays a crucial role in designing no-regret policies with bandit feedback [START_REF] Flaxman | Online convex optimization in the bandit setting: gradient descent without a gradient[END_REF][START_REF] Kleinberg | Nearly tight bounds for the continuum-armed bandit problem[END_REF], so δ must be chosen with care.

Before dealing with this choice though, it is important to highlight two feasibility issues that arise with the single-shot SPSA estimate (4.2). The first has to do with the fact that the perturbation direction z i is chosen from the unit sphere S i so it may fail to be tangent to X i , even when x i is interior. To iron out this wrinkle, it suffices to sample z i from the intersection of S i with the affine hull of X i in V i ; on that account (and without loss of generality), we will simply assume in what follows that each X i is a convex body of V i , i.e., it has nonempty topological interior.

The second feasibility issue concerns the size of the perturbation step: even if z i is a feasible direction of motion, the query point xi = x i + δz i may be unfeasible if x i is too close to the boundary of X i . For this reason, we will introduce a "safety net" in the spirit of [START_REF] Agarwal | Optimal algorithms for online convex optimization with multi-point bandit feedback[END_REF], and we will constrain the set of possible pivot points x i to lie within a suitably shrunk zone of X .

In detail, let B ri (p i ) be an r i -ball centered at p i ∈ X i so that B ri (p i ) ⊆ X i . Then, instead of perturbing x i by z i , we consider the feasibility adjustment

w i = z i -r -1 i (x i -p i ), (4.4)
and each player plays xi = x i + δw i instead of x i + δz i . In other words, this adjustment moves each pivot to

x δ i = x i -r -1 i δ(x i -p i ), i.e.
, O(δ)-closer to the interior base point p i , and then perturbs x δ i by δz i . Feasibility of the query point is then ensured by noting that

xi = x δ i + δz i = (1 -r -1 i δ)x i + r -1 i δ(p i + r i z i ), (4.5) so xi ∈ X i if δ/r i < 1 (since p i + r i z i ∈ B ri (p i ) ⊆ X i ).
5 For simplicity, we take δ equal for all players; the extension to player-specific δ is straightforward, so we omit it. 6 We tacitly assume here that the query directions z i ∈ S d i are drawn independently across players.

Algorithm 1: Multi-agent mirror descent with bandit feedback (player indices suppressed)

Require: step-size γn > 0, query radius δn > 0, safety ball Br(p) ⊆ X 1: choose X ∈ dom ∂h

# initialization

2: repeat at each stage n = 1, 2, . . . The difference between this estimator and the oracle framework we discussed above is twofold. First, each player's realized action is xi = x i + δw i , not x i , so there is a disparity between the point at which payoffs are queried and the action profile where the oracle is called. Second, the resulting estimator v is not unbiased, so the statistical assumptions (4.1) for a stochastic oracle do not hold. In particular, given the feasibility adjustment (4.4), the estimate (4.2) with x given by (4.5) satisfies

E[v i ] = ∇ i u δ i (x δ i ; x δ -i ), (4.6)
so there are two sources of systematic error: an O(δ) perturbation in the function, and an O(δ) perturbation of each player's pivot point from x i to x δ i . Hence, to capture both sources of bias and separate them from the random noise, we will write

vi = v i (x) + U i + b i (4.7)
where

U i = vi -E[v i ] and b i = ∇ i u δ i (x δ ) -∇ i u i (x).
We are thus led to the following manifestation of the bias-variance dilemma: the bias term b in (4.7) is O(δ), but the second moment of the noise term U is O(1/δ 2 ); as such, an increase in accuracy (small bias) would result in a commensurate loss of precision (large noise variance). Balancing these two factors will be a key component of our analysis.

Convergence analysis and results

Combining the learning framework of Section 3 with the single-shot gradient estimation machinery of Section 4, we obtain the following variant of (MD) with payoff-based, bandit feedback:

Xn = X n + δ n W n , X n+1 = P Xn (γ n vn ). (MD-b)
In the above, the perturbations W n and the estimates vn are given respectively by (4.4) and (4.2), i.e.,

W i,n = Z i,n -r -1 i (X i,n -p i ) vi,n = (d i /δ n )u i ( Xn ) Z i,n (5.1)
and Z i,n is drawn independently and uniformly across players at each stage n (see also Algorithm 1 for a pseudocode implementation and Fig. 1 for a schematic representation).

In the rest of this paper, our goal will be to determine the equilibrium convergence properties of this scheme in concave N -person games. Our first asymptotic result below shows that, under (MD-b), the players' learning process converges to Nash equilibrium in monotone games: Theorem 5.1. Suppose that the players of a monotone game G ≡ G(N , X , u) follow (MD-b) with step-size γ n and query radius δ n such that

X X 1 . X1 . γ(d/δ) û1 z 1 δz 1 . X 2 Π . X2 . X 3 γ(d/δ) û2 z 2 δz 2
lim n→∞ γ n = lim n→∞ δ n = 0, ∞ n=1 γ n = ∞, ∞ n=1 γ n δ n < ∞, and ∞ n=1 γ 2 n δ 2 n < ∞.
(5.2) Then, the sequence of realized actions Xn converges to Nash equilibrium with probability 1.

Even though the setting is different, the conditions (5.2) for the tuning of the algorithm's parameters are akin to those encountered in Kiefer-Wolfowitz stochastic approximation schemes and serve a similar purpose. First, the conditions lim n→∞ γ n = 0 and ∞ n=1 γ n = ∞ respectively mitigate the method's inherent randomness and ensure a horizon of sufficient length. The requirement lim n→∞ δ n = 0 is also straightforward to explain: as players accrue more information, they need to decrease the sampling bias in order to have any hope of converging. However, as we discussed in Section 4, decreasing δ also increases the variance of the players' gradient estimates, which might grow to infinity as δ → 0. The crucial observation here is that new gradients enter the algorithm with a weight of γ n so the aggregate bias after n stages is of the order of O(

n k=1 γ k δ k ) and its variance is O( n k=1 γ 2 k /δ 2 k )
. If these error terms can be controlled, there is an underlying drift that emerges over time and which steers the process to equilibrium. We make this precise in the supplement by using a suitably adjusted variant of the Bregman divergence as a quasi-Féjér energy function for (MD-b) and relying on a series of (sub)martingale convergence arguments to establish the convergence of Xn (first as a subsequence, then with probability 1).

Of course, since Theorem 5.1 is asymptotic in nature, it is not clear how to choose γ n and δ n so as to optimize the method's convergence rate. Heuristically, if we take schedules of the form γ n = γ/n p and δ n = δ/n q with γ, δ > 0 and 0 < p, q ≤ 1, the only conditions imposed by (5.2) are p + q > 1 and p -q > 1/2. However, as we discussed above, the aggregate bias in the algorithm after n stages is O(

n k=1 γ n δ n ) = O(1/n p+q-1 ) and its variance is O( n k=1 γ 2 k /δ 2 k ) = O(1/n 2p-2q-1
): if the conditions (5.2) are satisfied, both error terms vanish, but they might do so at very different rates. By equating these exponents in order to bridge this gap, we obtain q = p/3; moreover, since the single-shot SPSA estimator (4.2) introduces a Θ(δ n ) random perturbation, q should be taken as large as possible to ensure that this perturbation vanishes at the fastest possible rate. As a result, the most suitable choice for p and q seems to be p = 1, q = 1/3, leading to an error bound of O(1/n 1/3 ).

We show below that this bound is indeed attainable for games that are strongly monotone, i.e., they satisfy the following stronger variant of diagonal strict concavity:

i∈N λ i v i (x ) -v i (x), x i -x i ≤ - β 2 x -x 2 (β-DSC)
for some λ i , β > 0 and for all x, x ∈ X . Focusing for expository reasons on the most widely used, Euclidean incarnation of the method (Example 3.1), we have:

Theorem 5.2. Let x * be the (necessarily unique) Nash equilibrium of a β-strongly monotone game. If the players follow (MD-b) with Euclidean projections and parameters γ n = γ/n and δ n = δ/n 1/3 with γ > 1/(3β) and δ > 0, we have

E[ Xn -x * 2 ] = O(n -1/3
).

(5.3) Theorem 5.2 is our main finite-time analysis result, so some remarks are in order. First, the step-size schedule γ n ∝ 1/n is not required to obtain an O(n -1/3 ) convergence rate: as we show in the paper's supplement, more general schedules of the form γ n ∝ 1/n p and δ n ∝ 1/n q with p > 3/4 and q = p/3 > 1/4, still guarantee an O(n -1/3 ) rate of convergence for (MD-b). To put things in perspective, we also show in the supplement that if (MD) is run with first-order oracle feedback satisfying the statistical assumptions (4.1), the rate of convergence becomes O(1/n). Viewed in this light, the price for not having access to gradient information is no higher than O(n -2/3 ) in terms of the players' equilibration rate.

Finally, it is also worth comparing the bound (D.2) to the attainable rates for stochastic convex optimization (the single-player case). For problems with objectives that are both strongly convex and smooth, [START_REF] Agarwal | Optimal algorithms for online convex optimization with multi-point bandit feedback[END_REF] attained an O(n -1/2 ) convergence rate with bandit feedback, which [START_REF] Shamir | On the complexity of bandit and derivative-free stochastic convex optimization[END_REF] showed is unimprovable. Thus, in the single-player case, the bound (D.2) is off by n 1/6 and coincides with the bound of [START_REF] Agarwal | Optimal algorithms for online convex optimization with multi-point bandit feedback[END_REF] for strongly convex functions that are not necessarily smooth. One reason for this gap is that the Θ(n -1/2 ) bound of Shamir ( 2013) concerns the smoothed-out time average Xn = n -1 n k=1 X k , while our analysis concerns the sequence of realized actions Xn . This difference is semantically significant: In optimization, the query sequence is just a means to an end, and only the algorithm's output matters (i.e., Xn ). In a game-theoretic setting however, it is the players' realized actions that determine their rewards at each stage, so the figure of merit is the actual sequence of play Xn . This sequence is more difficult to control, so this disparity is, perhaps, not too surprising; nevertheless, we believe that this gap can be closed by using a more sophisticated single-shot estimate, e.g., as in [START_REF] Ghadimi | Stochastic first-and zeroth-order methods for nonconvex stochastic programming[END_REF]. We defer this analysis to the future.

Concluding remarks

The most sensible choice for agents who are oblivious to the presence of each other (or who are simply conservative), is to deploy a no-regret learning algorithm. With this in mind, we studied the long-run behavior of individual regularized noregret learning policies and we showed that, in monotone games, play converges to equilibrium with probability 1, and the rate of convergence almost matches the optimal rates of single-agent, stochastic convex optimization. Nevertheless, several questions remain open: whether there is an intrinsic information-theoretic obstacle to bridging this gap; whether our convergence rate estimates hold with high probability (and not just in expectation); and whether our analysis extends to a fully decentralized setting where the players' updates need not be synchronous. We intend to address these questions in future work.

Appendix A. Monotone games

Our aim in this appendix is to show that the game-theoretic examples of Section 2 are both monotone. Before studying them in detail, it will be convenient to introduce a straightforward second-order test for monotonicity based on the game's Hessian matrix.

Specifically, extending the notion of the Hessian of an ordinary (scalar) function, the (λ-weighted ) Hessian of a game G is defined as the block matrix

H G (x; λ) = (H ij (x; λ)) i,j∈N with blocks H ij (x; λ) = λ i 2 ∇ j ∇ i u i (x) + λ j 2 (∇ i ∇ j u j (x)) . (A.1)
As was shown by Rosen (1965, Theorem 6), G satisifes (DSC) with weight vector λ whenever z H G (x; λ)z < 0 for all x ∈ X and all nonzero z ∈ V ≡ i V i that are tangent to X at x.7 It is thus common to check for monotonicity by taking λ i = 1 for all i ∈ N and verifying whether the unweighted Hessian of G is negative-definite on the affine hull of X .

A.1. Cournot competition (Example 2.1). In the standard Cournot oligopoly model described in the main body of the paper, the players' payoff functions are given by

u i (x) = x i a -b j x j -c i x i . (A.2)
Consequently, a simple differentiation yields

H ij (x) = 1 2 ∂ 2 u i ∂x i ∂x j + 1 2 ∂ 2 u j ∂x j ∂x i = -b(1 + δ ij ), (A.3)
where δ ij = 1{i = j} is the Kronecker delta. This matrix is clearly negative-definite, so the game is monotone.

A.2. Resource allocation auctions (Example 2.2). In our auction-theoretic example, the players' payoff functions are given by

u i (x i ; x -i ) = s∈S g i q s x is c s + j∈N x js -x is (A.4)
To prove monotonicity in this example, we will consider the following criterion due to [START_REF] Goodman | Note on existence and uniqueness of equilibrium points for concave N -person games[END_REF]: a game G satisfies (DSC) with weights λ i , i ∈ N , if: a) Each payoff function u i is strictly concave in x i and convex in x -i .

b) The function i∈N λ i u i (x) is concave in x. Since the function φ(x) = x/(c + x) is strictly concave in x for all c > 0, the first condition above is trivial to verify. For the second, letting

λ i = 1/g i gives i∈N λ i u i (x) = i∈N s∈S q s x is c s + j∈N x js - i∈N s∈S x is = s∈S q s i∈N x is c s + i∈N x is - i∈N s∈S x is . (A.5)
Since the summands above are all concave in their respective arguments, our claim follows.

Proof of Lemma B.1.

if x + = Q(∇h(x) + y).
For the inequality (B.5), it suffices to show it holds for interior p ∈ X • (by continuity). To do so, let

φ(t) = h(x + t(p -x)) -[h(x) + y, x + t(p -x) ]. (B.6)
Since h is strongly convex and y ∈ ∂h(x) by (B.4a), it follows that φ(t) ≥ 0 with equality if and only if t = 0. Moreover, note that ψ(t) = ∇h(x + t(p -x)) -y, p -x is a continuous selection of subgradients of φ. Given that φ and ψ are both continuous on [0, 1], it follows that φ is continuously differentiable and φ = ψ on [0, 1]. Thus, with φ convex and φ(t) ≥ 0 = φ(0) for all t ∈ [0, 1], we conclude that φ (0) = ∇h(x) -y, p -x ≥ 0, from which our claim follows.

We continue with some basic relations connecting the Bregman divergence relative to a target point before and after a prox step. The basic ingredient for this is a generalization of the law of cosines which is known in the literature as the "three-point identity" [START_REF] Chen | Convergence analysis of a proximal-like minimization algorithm using Bregman functions[END_REF]:

Lemma B.2.
Let h be a regularizer on X . Then, for all p ∈ X and all x, x ∈ dom ∂h, we have

D(p, x ) = D(p, x) + D(x, x ) + ∇h(x ) -∇h(x), x -p . (B.7)
Proof. By definition, we get:

D(p, x ) = h(p) -h(x ) -∇h(x ), p -x D(p, x) = h(p) -h(x) -∇h(x), p -x D(x, x ) = h(x) -h(x ) -∇h(x ), x -x . (B.8)
The lemma then follows by adding the two last lines and subtracting the first.

With all this at hand, we have the following upper and lower bounds:

Proposition B.3. Let h be a K-strongly convex regularizer on X , fix some p ∈ X , and let x + = P x (y) for x ∈ dom ∂h, y ∈ Y. Then, we have:

D(p, x) ≥ K 2 x -p 2 . (B.9a) D(p, x + ) ≤ D(p, x) -D(x + , x) + y, x + -p (B.9b) ≤ D(p, x) + y, x -p + 1 2K y 2 * (B.9c)
Proof of (B.9a). By the strong convexity of h, we get where, in the last line, we used the identity

h(p) ≥ h(x) + ∇h(x), p -x + K 2 p -x 2 (B.
∇ δB f (x + w) dw = δS f (x + z) z z dz (C.3)
which, in turn, follows from Stokes' theorem [START_REF] Flaxman | Online convex optimization in the bandit setting: gradient descent without a gradient[END_REF][START_REF] Lee | Introduction to Smooth Manifolds[END_REF]. Since vol(δB i ) = (δ/d i ) vol(δS i ), the above yields E[v i ] = ∇ i u δ i with u δ i given by (4.3). For the second part of the lemma, let L i denote the Lipschitz constant of v i , i.e., v i (x ) -v i (x) * ≤ L i x -x for all x, x ∈ X . Then, for all w i ∈ δB i and all z j ∈ δS j , j = i, we have

∇ i u i (x i + w i ; x -i + z -i ) -∇ i u i (x) ≤ L i w i 2 + j =i z j 2 ≤ L i √ N δ. (C.4)
Our assertion then follows by integrating and differentiating under the integral sign.

With this basic estimate at hand, we proceed to establish the convergence of the Bregman divergence relative to the game's Nash equilibria:

Proposition C.2. Let x * be a Nash equilibrium of G. Then, with assumptions as in Theorem 5.1, the Bregman divergence D(x * , X n ) converges (a.s.) to a finite random variable D ∞ .

Remark. For expository reasons, we tacitly assume above (and in what follows) that G satisfies (DSC) with weights λ i = 1 for all i ∈ N . If this is not the case, the Bregman divergence D(p, x) should be replaced by the weight-adjusted variant

D λ (p, x) = i∈N λ i D(p i , x i ).
(C.5)

Since this adjustment would force us to carry around all player indices, the presentation would become significantly more cumbersome; to avoid this, we stick with the simpler, unweighted case.

Proof. Let D n = D(x * , X n ) for some Nash equilibrium x * of G and write

vn = v(X n ) + U n+1 + b n , (C.6)
where, recalling the setup of Section 4 in the main body of the paper, the noise process

U n+1 = vn -E[v n | F n ] is an F n -adapted martingale difference sequence and b n = v δn (X δn n ) -v(X n )
denotes the systematic bias of the estimator vn . 8 Then, by Proposition B.3, we have

D n+1 = D(x * , P Xn (γ n vn )) ≤ D(x * , X n ) + γ n vn , X n -x * + γ 2 n 2K vn 2 * = D n + γ n v(X n ) + U n+1 + b n , X n -x * + γ 2 n 2K vn 2 * ≤ D n + γ n ξ n+1 + γ n r n + γ 2 n 2K vn 2 * , (C.7)
8 Recall here that X δ i , i ∈ N , denotes the δ-adjusted pivot X δ i = X i + r -1 i δ(X i -p i ), i.e., including the feasibility adjustment r -1 i (X i -p i ).

where, in the last line, we set ξ n+1 = U n+1 , X n -x * , r n = b n , X n -x * , and we used the variational characterization (VI) of Nash equilibria of monotone games. Thus, conditioning on F n and taking expectations, we get

E[D n+1 | F n ] ≤ D n + E[ξ n+1 | F n ] + γ n E[r n | F n ] + γ 2 n 2K E[ vn 2 * | F n ] ≤ D n + γ n E[r n | F n ] + V 2 2K γ 2 n δ 2 n . (C.8)
where we set

V 2 = i d 2 i max x∈X |u i (x)
| 2 and we used the fact that

X n is F n - measurable, so E[ξ n+1 | F n ] = E[U n+1 | F n ], X n -x * = 0.
(C.9) Finally, by Lemma C.1, we have

b n * = v δn (X δn n ) -v(X n ) * ≤ v δn (X δn n ) -v(X δn n ) * + v(X δn n ) -v(X n ) * = O(δ n ), (C.10)
where we used the fact that v is Lipschitz continuous and v δ -v ∞ = O(δ). This shows that there exists some B > 0 such that r n ≤ Bδ n ; as a consequence, we obtain (C.12) i.e., R n is an F n -adapted supermartingale.9 Since the series

E[D n+1 | F n ] ≤ D n + Bγ n δ n + V 2 2K γ 2 n δ 2 n . (C.11) Now, letting R n = D n + ∞ k=n [Bγ k δ k + (2K) -1 V 2 γ 2 k /δ 2 k ], the estimate (C.7) gives E[R n+1 | F n ] = E[D n+1 | F n ] + ∞ k=n+1 Bγ k δ k + V 2 2K γ 2 k δ 2 k ≤ D n + Bγ n δ n + V 2 2K γ 2 n δ 2 n + ∞ k=n+1 Bγ k δ k + V 2 2K γ 2 k δ 2 k ≤ D n + ∞ k=n Bγ k δ k + V 2 2K γ 2 k δ 2 k = R n ,
∞ n=1 γ n δ n and ∞ n=1 γ 2 n /δ 2 n are both summable, it follows that E[R n ] = E[E[R n | F n-1 ]] ≤ E[R n-1 ] ≤ • • • ≤ E[R 1 ] ≤ E[D 1 ] + ∞ n=1 Bγ n δ n + V 2 2K γ 2 n δ 2 n < ∞ (C.13)
i.e., R n is uniformly bounded in L 1 . Thus, by Doob's convergence theorem for supermartingales (Hall and Heyde, 1980, Theorem 2.5), it follows that R n converges (a.s.) to some finite random variable R ∞ . In turn, by inverting the definition of R n , it follows that D n converges (a.s.) to some random variable D ∞ , as claimed.

Proposition C.3. Suppose that the assumptions of Theorem 5.1 hold. Then, with probability 1, there exists a (random) subsequence X n k of (MD-b) which converges to Nash equilibrium.

Proof. We begin with the technical observation that the set X * of Nash equilibria of G is closed (and hence, compact). Indeed, let x * n , n = 1, 2, . . . , be a sequence of Nash equilibria converging to some limit point x * ∈ X ; to show that X * is closed, it suffices to show that x * ∈ X . However, since Nash equilibria of G satisfy the variational characterization (VI), we also have

v(x), x -x * n ≤ 0 for all x ∈ X . Hence, with x * n → x * as n → ∞, it follows that v(x), x -x * = lim n→∞ v(x), x -x * n ≤ 0 for all x ∈ X , (C.14)
i.e., x * satisfies (VI). Since G is monotone, we conclude that x * is a Nash equilibrium, as claimed. Suppose now ad absurdum that, with positive probability, the pivot sequence X n generated by (MD-b) admits no limit points in X * . 10 Conditioning on this event, and given that X * is compact, there exists a (nonempty) compact set C ⊂ X such that C ∩ X * = ∅ and X n ∈ C for all sufficiently large n. Moreover, by (VI), we have v(x), x -x * < 0 whenever x ∈ C and x * ∈ X * . Therefore, by the continuity of v and the compactness of X * and C, there exists some c > 0 such that

v(x), x -x * ≤ -c for all x ∈ C, x * ∈ X . (C.15)
To proceed, fix some x * ∈ X * and let D n = D(x * , X n ) as in the proof of Proposition C.2. Then, telescoping (C.7) yields the estimate

D n+1 ≤ D 1 + n k=1 γ k v(X n ), X n -x * + n k=1 γ k ξ k+1 + n k=1 γ k r k + n k=1 γ 2 k 2K vn 2 * , (C.16)
where, as in the proof of Proposition C.2, we set

ξ n+1 = U n+1 , X n -x * (C.17
D n+1 ≤ D 1 -τ n c - n k=1 γ k ξ k+1 τ n - n k=1 γ k r k τ n - (2K) -1 n k=1 γ 2 k vk 2 * τ n .
(C.19) Since U n is a martingale difference sequence with respect to F n , we have E[ξ n+1 | F n ] = 0 (recall that X n is F n -measurable by construction). Moreover, by construction, there exists some constant σ > 0 such that

U n+1 2 * ≤ σ 2 δ 2 n , (C.20)
and hence:

∞ n=1 γ 2 n E[ξ 2 n+1 | F n ] ≤ ∞ n=1 γ 2 n X n -x * 2 E[ U n+1 2 * | F n ] ≤ diam(X ) 2 σ 2 ∞ n=1 γ 2 n δ 2 n < ∞. (C.21)
Therefore, by the law of large numbers for martingale difference sequences (Hall and Heyde, 1980, Theorem 2.18), we conclude that τ -1 n n k=1 γ k ξ k+1 converges to 0 with probability 1.

For the third term in the brackets of (C.19) we have r n → 0 as n → ∞ (a.s.). Since

∞ n=1 γ n = ∞, it follows n k=1 γ k r k n k=1 γ k → 0.
Finally, for the last term in the brackets of (C.19), let

S n+1 = n k=1 γ 2 k vk 2 * . Since vk is F n -measurable for all k = 1, 2, . . . , n -1, we have E[S n+1 | F n ] = E n-1 k=1 γ 2 k vk 2 * + γ 2 n vn 2 * F n = S n + γ 2 n E[ vn 2 * | F n ] ≥ S n , (
C.22) i.e., S n is a submartingale with respect to F n . Furthermore, by the law of total expectation, we also have

E[S n+1 ] = E[E[S n+1 | F n ]] ≤ V 2 n k=1 γ 2 k δ 2 k ≤ V 2 ∞ k=1 γ 2 k δ 2 k < ∞, (C.23)
implying in turn that S n is uniformly bounded in L 1 . Hence, by Doob's submartingale convergence theorem [START_REF] Hall | Martingale Limit Theory and Its Application[END_REF], Theorem 2.5), we conclude that S n converges to some (almost surely finite) random variable S ∞ with E[S ∞ ] < ∞. Consequently, we have lim n→∞ S n+1 /τ n = 0 with probability 1. Applying all of the above to the estimate (C.19), we get D n+1 ≤ D 1 -cτ n /2 for sufficiently large n, and hence, D(x * , X n ) → -∞, a contradiction. Going back to our original assumption, this shows that at least one of the limit points of X n must lie in X * , so our proof is complete.

We are finally in a position to prove Theorem 5.1 regarding the convergence of (MD-b):

Proof of Theorem 5.1. By Proposition C.3, there exists a (possibly random) Nash equilibrium x * of G such that X n k -x * → 0 for some (random) subsequence X n k . By the assumed reciprocity of the Bregman divergence, this implies that lim inf n→∞ D(x * , X n ) = 0 (a.s.). Since lim n→∞ D(x * , X n ) exists with probability 1 (by Proposition C.2), it follows that where we set Q n = Q(1+q n /n), so Q n → Q as n → ∞. Then, under the assumption that P > q when p = 1, (D.5) can be rewritten as (D.6) for some sequence R n with R n → R as n → ∞. Now, fix some small enough ε > 0. From (D.6), we readily get

ãn+1 ≤ ãn 1 - R n n p + Q n n p ,
ãn+1 ≤ ãn - R n ãn -Q n n p . (D.7)
Since R n → R and Q n → Q as n → ∞, we will have R n > R -ε and Q n < Q + ε for all n greater than some n ε . Thus, if n ≥ n ε and (R -ε)ã n -(Q + ε) > ε, we will also have

ãn+1 ≤ ãn - R n ãn -Q n n p ≤ ãn - (R -ε)ã n -(Q + ε) n p ≤ ãn - ε n p . (D.8)
The above shows that, as long as ãn > (Q + 2ε)/(R -ε), ãn will decrease at least by ε/n p at each step. In turn, since ∞ n=1 (1/n p ) = ∞, it follows by telescoping that lim sup n→∞ ãn ≤ (Q + 2ε)/(R -ε). Hence, with ε arbitrary, we conclude that lim sup n→∞ a n n q ≤ Q/R, as claimed.

Proof of Theorem D.1. We begin with the second part of the theorem; the first part will follow by setting some estimates equal to zero, so the analysis is more streamlined that way. Also, as in the previous section, we tacitly assume that (β-DSC) holds with weights λ i = 1 for all i ∈ N . If this is not the case, the Bregman divergence D(p, x) should be replaced by the weight-adjusted variant (C.5), but this would only make the presentation more difficult to follow, so we omit the details.

The main component of our proof is the estimate (C.7), which, for convenience (and with notation as in the previous section), we also reproduce below:

D n+1 ≤ D n + γ n v(X n ), X n -x * + γ n ξ n+1 + γ n r n + γ 2 n 2K vn 2 * .
(D.9)

In the above, since the algorithm is run with Euclidean projections, D n = 1 2 X nx * 2 ; other than that, ξ n and r n are defined as in (C.17) and (C.18) respectively. Since the game is β-strongly monotone and x * is a Nash equilibrium, we further have with B and V defined as in the proof of Theorem 5.1 in the previous section. Now, substituting γ n = γ/n p and δ n = δ/n q in (D.12) readily yields Dn+1 ≤ 1 -βγ n p Dn + Bγδ n p+q + V 2 γ 2 δ 2 2Kn 2(p-q) . (D.13) Hence, taking p = 1 and q = 1/3, the last two exponents are equated, leading to the estimate Dn+1 ≤ 1 -βγ n Dn + C n 4/3 , (D.14)

v(X n ), X n -x * ≤ v(X n ) -v(x * ), X n -x * ≤ - β 2 X n -x * 2 = -
with C = γδB + (2K) -1 γ 2 δ 2 V 2 . Thus, with βγ > 1/3, applying Lemma D.2 with p = 1 and q = 1/3, we finally obtain Dn = O(1/n 1/3 ).

The proof for the oracle case is similar: the key observation is that the bound (D.12) becomes (D.15) with V defined as in (4.1). Hence, taking γ n = γ/n with βγ > 1 and applying again Lemma D.2 with p = q = 1, we obtain Dn = O(1/n) and our proof is complete.

Dn+1 ≤ (1 -βγ n ) Dn + V 2 2K γ 2 n ,
To conclude, we note that the O(1/n 1/3 ) bound of Theorem D.1 cannot be readily improved by choosing a different step-size schedule of the form γ n ∝ 1/n p for some p < 1. Indeed, applying Lemma D.2 to the estimate (D.13) yields a bound which is either O(1/n q ) or O(1/n p-2q ), depending on which exponent is larger. Equating the two exponents (otherwise, one term would be slower than the other), we get q = p/3, leading again to a O(1/n 1/3 ) bound. Unless one has finer control on the bias/variance of the SPSA gradient estimator used in (MD-b), we do not see a way of improving this bound in the current context.

Figure 1 :

 1 Figure 1: Schematic representation of Algorithm 1 with ordinary, Euclidean projections. To reduce visual clutter, we did not include the feasibility adjustment r -1 (x -p) in the action selection step Xn → Xn.

  10) so (B.9a) follows by gathering all terms involving h and recalling the definition of D(p, x). Proof of (B.9b) and (B.9c). By the three-point identity (B.7), we readily obtain D(p, x) = D(p, x + ) + D(x + , x) + ∇h(x) -∇h(x + ), x + -p , (B.11) and hence:

  ) and r n = b n , X n -x * . (C.18) Subsequently, letting τ n = n k=1 γ k and using (C.15), we obtain

  * , X n ) = lim inf n→∞ D(x * , X n ) = 0, (C.24)i.e., X n converges to x * by the first part of Proposition B.3. Since δ n → 0 and Xn -X n = δ n W n = O(δ n ), our claim follows.

  βD n , (D.10) so (D.9) becomesD n+1 ≤ (1 -βγ n )D n + γ n ξ n+1 + γ n r n + Dn = E[D n ]and taking expectations, we obtainDn+1 ≤ (1 -βγ n ) Dn + Bγ n δ n +

  To prove (B.4a), note that x solves (B.2) if and only if y -∂h(x) 0, i.e., if and only if y ∈ ∂h(x). Similarly, for (B.4b), comparing (3.3) and (B.1), we see that x + solves (3.3) if and only if ∇h(x) + y ∈ ∂h(x + ), i.e., if and only

By "tangent" we mean here that z belongs to the tangent cone TC(x) to X at x, i.e., the intersection of all supporting (closed) half-spaces of X at x.

In particular, this shows that E[Dn | F n-1 ] is quasi-Fejér in the sense of[START_REF] Combettes | Quasi-Fejérian analysis of some optimization algorithms[END_REF].

M. Bravo gratefully acknowledges the support provided by FONDECYT grant 11151003. P. Mertikopoulos was partially supported by the Huawei HIRP flagship grant ULTRON, and the French National Research Agency (ANR) grant ORACLESS (ANR-16-CE33-0004-01). Part of this work was carried out with financial support by the ECOS project C15E03.

1

Appendix B. Properties of Bregman proximal mappings

In this appendix, we provide some auxiliary results and estimates that are used throughout the convergence analysis of Appendix C. Some of the results we present here are not new (see e.g., [START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF]; however, the set of hypotheses used to obtain them varies widely in the literature, so we provide all proofs for completeness.

In what follows, we will make frequent use of the convex conjugate h * : Y → R of h, defined here as h * (y) = max By standard results in convex analysis (Rockafellar, 1970, Chap. 26), h * is differentiable on Y and its gradient satisfies the identity

For notational convenience, we will also write

3) and we will refer to Q : Y → X as the mirror map generated by h.

Together with the prox-mapping induced by h, all these notions are related as follows:

Lemma B.1. Let h be a regularizer on X . Then, for all x ∈ dom ∂h, y ∈ Y, we have:

Finally, if x = Q(y) and p ∈ X , we have

Remark. Note that (B.4b) directly implies that ∂h(x + ) = ∅, i.e., x + ∈ dom ∂h. An immediate consequence of this is that the update rule x ← P x (y) is well-posed, i.e., it can be iterated in perpetuity.

where, in the last step, we used (B.5) and the fact that x + = Q(∇h(x) + y), by (B.4b), since x + = P x (y). The above is just (B.9b), so the first part of our proof is complete.

To proceed with the proof of (B.9c), note that (B.12) gives

By Young's inequality [START_REF] Rockafellar | Convex Analysis[END_REF], we also have

and hence

with the last step following from Lemma B.1 after plugging in x in place of p.

Appendix C. Asymptotic convergence analysis

Our goal in this appendix is to prove Theorem 5.1. Our proof strategy will be based on a two-pronged approach. First, we will show that the pivot sequence X n satisfies a "quasi-Fejér" property [START_REF] Combettes | Quasi-Fejérian analysis of some optimization algorithms[END_REF][START_REF] Combettes | Stochastic quasi-Fejér block-coordinate fixed point iterations with random sweeping[END_REF] with respect to the Bregman divergence. This quasi-Fejér property allows us to show that the Bregman divergence D(x * , X n ) with respect to a Nash equilibrium x * of G converges. To show that this limit is actually zero for some Nash equilibrium, we prove that, with probability 1, the sequence X n admits a (random) subsequence that converges to a Nash equilibrium. The theorem then follows by combining these two results.

To carry all this out, we begin with an auxiliary lemma for the SPSA estimation process of Section 4:

with u δ i as in (4.3). Moreover, we have

By the independence of the sampling directions z i , i ∈ N , we have

Appendix D. Finite-time analysis and rates of convergence

We now turn to the finite-time analysis of (MD-b). To begin, we briefly recall that a game G is β-strongly monotone if it satisfies the condition

for some λ i , β > 0 and for all x, x ∈ X . Our aim in what follows will be to prove the following convergence rate estimate for multi-agent mirror descent in strongly monotone games:

Theorem D.1. Let x * be the (unique) Nash equilibrium of a β-strongly monotone game. Then: a) If the players have access to a gradient oracle satisfying (4.1) and they follow (MD) with Euclidean projections and step-size sequence γ n = γ/n for some γ > 1/β, we have

b) If the players only have bandit feedback and they follow (MD-b) with Euclidean projections and parameters γ n = γ/n and δ n = δ/n 1/3 with γ > 1/(3β) and δ > 0, we have

Remark. Theorem 5.2 is recovered by the second part of Theorem D.1 above; the first part (which was alluded to in the main paper) serves as a benchmark to quantify the gap between bandit and oracle feedback.

For the proof of Theorem D.1 we will need the following lemma on numerical sequences, a version of which is often attributed to [START_REF] Chung | On a stochastic approximation method[END_REF]:

Lemma D.2. Let a n , n = 1, 2, . . . , be a non-negative sequence such that

where 0 < p ≤ 1, q > 0, and P, Q > 0. Then, assuming P > q if p = 1, we have

Proof. Clearly, it suffices to show that lim sup n→∞ n q a n ≤ Q/R. To that end, write q n = n[(1 + 1/n) q -1], so (1 + 1/n) q = 1 + q n /n and q n → q as n → ∞. Then, multiplying both sides of (D.3) by (n + 1) q and letting ãn = a n n q , we get


Andrew Ward

Zhengyuan Zhou

Panayotis Mertikopoulos

Nicholas Bambos

Power Control with Random Delays: Robust Feedback Averaging

Distributed power control schemes in wireless networks have been well-examined, but standard methods rarely consider the effect of potentially random delays, which occur in almost every real-world network. We present Robust Feedback Averaging, a novel power control algorithm that is capable of operating in delay-ridden and noisy environments. We prove optimal convergence of this algorithm in the presence of random, time-varying delays, and present numerical simulations that indicate that Robust Feedback Averaging outperforms the ubiquitous Foschini-Miljanic algorithm in several regimes.

I. INTRODUCTION

Proper management of wireless communication networks relies on efficient and robust power control strategies. Such strategies may be broadly categorized into either centralized or distributed paradigms. While centralized methods may allow for more globally optimal theoretical results, this is counterbalanced by an increased communication requirement for implementation, and a reduction in system robustness due to the presence of a single total failure point. In contrast, divesting the control into a distributed scheme alleviates these concerns and opens up the study of both cooperative and competitive optimal decentralized control strategies [START_REF] Goldsmith | Wireless communications[END_REF], [START_REF] Rappaport | Wireless communications: principles and practice[END_REF]. Such strategies have been well-examined, and solutions for the fundamental problems have been established.

One foremost power control method in such a distributed wireless network was developed by Foschini and Miljanic (FM) [START_REF] Foschini | A simple distributed autonomous power control algorithm and its convergence[END_REF]. Since its development, the FM algorithm has proved applicable and stable in far more general environments than originally considered [START_REF] Zhou | A stochastic stability characterization of the foschini-miljanic algorithm in random wireless networks[END_REF]. As a clean and polished control scheme which possesses several strong and useful convergence properties, the FM algorithm has heavily inspired subsequent research in the field [START_REF] Yates | A framework for uplink power control in cellular radio systems[END_REF], [START_REF] Zhou | Repeated games for power control in wireless communications: Equilibrium and regret[END_REF], [START_REF] Ulukus | Stochastic power control for cellular radio systems[END_REF]. Extensive studies have examined the power control problem from game-theoretic and utility based perspectives, allowing the consideration of cooperative and competitive network participants [START_REF] Han | Game theory in wireless and communication networks: theory, models, and applications[END_REF], [START_REF] Menache | Network games: Theory, models, and dynamics[END_REF], [START_REF] Zhou | Wireless communications games in fixed and random environments[END_REF], [START_REF] Zhou | Deterministic and stochastic wireless networks games: Equilibrium, dynamics and price of anarchy[END_REF]. These control methods generally focus on "memoryless" implementations, which only consider the current power signal vector when deciding the optimal power for a subsequent iteration. This restriction results in well-known stability issues.

In a 2017 study which addresses these stability issues, Zhou et al. presented a power control algorithm known as dual averaging (DA) [START_REF] Zhou | Stable power control in wireless networks via dual averaging[END_REF]. This algorithm uses historical A. Ward, Z. Zhou, and N. Bambos are with the Department of Electrical Engineering, Stanford University, Stanford, CA, 94305, USA. {atward, zyzhou, bambos}@stanford.edu P. Mertikopoulos is with the French National Center for Scientific Research (CNRS) and with Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LIG, F-38000, Grenoble, France. A. Ward is funded by the Department of Defense through an NDSEG fellowship information in a way that prevents the current transmitted power from being dominated by any particular previously received signal. The algorithm strictly enforces fully distributed control and eliminates extraneous communication requirements. Furthermore, the total historical information is not stored explicitly, but is carried forward implicitly and efficiently in rolling summaries to reduce memory requirements.

A further concern beyond memory requirements is the delays inherent in any physically realizable system. These delays can result from either signal propagation time or from the computational requirements to implement such power control algorithms in any real-world device. Previous literature has studied the effect of relatively simple constant and time-varying delays on the FM algorithm and adaptions thereof [START_REF] Charalambous | On the stability of the foschini-miljanic algorithm with time-delays[END_REF], [START_REF] Charalambous | On the stability of a power control algorithm for wireless networks in the presence of time-varying delays[END_REF], [START_REF] Campos-Delgado | Characterization of stability margins of the foschini-miljanic power allocation strategy under constant and time-varying delays[END_REF], [START_REF] Mitra | An asynchronous distributed algorithm for power control in cellular radio systems[END_REF]. However, modern mobile and mesh wireless networks may consist of a wide range of distinct devices, each attempting to send different types of information. This may lead to a spread of delays culminating in either no power signals in a time slot, or multiple signals being delayed to the same subsequent time slot. Furthermore, it could also lead to power signals arriving out of order. Since the presence of delay is so ubiquitous in real-world systems, any new power control algorithm should be shown to handle the concerns that delayed information introduces. Several previous studies have presented neural-network techniques to combat time-varying delay [START_REF] Wei | Fixed-order piecewise-affine output feedback controller for fuzzy-affine-model-based nonlinear systems with time-varying delay[END_REF], [START_REF] Qiu | Data-based optimal control for networked double-layer industrial processes[END_REF], but these can be computationally intensive and require parameter tuning.

Our goal is to present a new power control algorithm that maintains the desirable stability and efficiency of the DA algorithm [START_REF] Zhou | Stable power control in wireless networks via dual averaging[END_REF], and is capable of handling arbitrary delays. We must first address the situation in which multiple power iterates are delayed such that they arrive within a single time step. In this case, the control algorithm must either choose from or combine all available data without prior knowledge of the order in which the iterates were generated. We must also allow for the situation in which power iterates are received out of order. These considerations allow our new algorithm to address a wide range of both deterministic and random delays, a highly desirable capability which greatly increases the robustness of the system.

Delays have been studied in wireless networks in the context of the FM algorithm [START_REF] Charalambous | On the stability of the foschini-miljanic algorithm with time-delays[END_REF], [START_REF] Charalambous | On the stability of a power control algorithm for wireless networks in the presence of time-varying delays[END_REF], [START_REF] Campos-Delgado | Characterization of stability margins of the foschini-miljanic power allocation strategy under constant and time-varying delays[END_REF], [START_REF] Mitra | An asynchronous distributed algorithm for power control in cellular radio systems[END_REF], though none have addressed the case in which multiple SINR measurements arrive simultaneously or the case where SINR measurements arrive out of order. So far, properties of the DA algorithm in a network with delays have not been studied.

A. Our Contributions

Our contributions are threefold: we first present Robust Feedback Averaging (RFA), a novel algorithm which utilizes the core idea of the dual averaging distributed power control algorithm. This new algorithm expands the DA algorithm's capabilities, and allows it to address a wide range of delay types. RFA is a fully asynchronous control algorithm, and this asynchronicity is central in RFA's ability to handle arbitrary delays. This property relaxes the need for update synchronization, and even allows different links to perform power updates at different rates. Furthermore, the RFA algorithm maintains both the constant memory property and the baseline power convergence guarantee of the original DA algorithm in this more general environment.

Second, we prove that the RFA algorithm maintains the original performance guarantees: if the channel is feasible, the RFA algorithm converges to the optimal power vector. In addition, we prove a new guarantee for the RFA algorithm: even subject to random delays, the algorithm converges almost surely to the optimal power vector. These guarantees are not matched by the FM algorithm, or by parallel adaptations of such. While the FM algorithm has been shown to converge under certain time-varying delays, it may fail to converge when said delays are arbitrarily random. Additionally, the FM algorithm as it stands cannot address the situation of multiple iterates within a time slot, and even if augmented in a manner similar to the RFA algorithm, does not maintain its convergence guarantee.

Thirdly, we demonstrate the robustness of the RFA algorithm in the presence of both random channels and delays via a series of numerical simulations. These simulations demonstrate that the RFA algorithm achieves convergence in a variety of scenarios, both ones that we prove will converge and ones where the RFA algorithm is not proven to converge. We also show that our algorithm empirically outperforms a corresponding augmentation of the FM algorithm.

II. BACKGROUND, MODEL, MOTIVATION A. Power Control in Wireless Networks

We consider a typical wireless network model consisting of N communication links, each comprised of a transmitter and an intended receiver. The transmitted power vector is denoted by p = (p 1 , . . . , p N ) where p i is the power used by the transmitter of link i. For this paper, we consider only non-negative real transmission powers. We also enforce a maximum power p max i for each link, which is a standard practical constraint [12] [19]. The feasible support set for the power vector, then, is

P = N i=1 [0, p max i ].
The most common service quality measure used in power control systems is signal to interference and noise ratio (SINR). Given a power vector p, link i's SINR r i (p) is given by the following ratio:

r i (p) = G ii p i j =i G ij p j + η i (1)
where η i is the thermal noise associated with the receiver of link i and G ij ≥ 0 is the power gain between transmitter j and receiver i; for i = j, G ij represents the interference receiver i experiences due to transmitter j per unit transmission power used. Here, we represent all the power gains G ij with the gain matrix G and all the thermal noises with the vector η.

Additionally, each link i has a threshold r * i > 0, which denotes the minimum SINR necessary to achieve acceptable service quality for link i. The standard power control problem [START_REF] Weeraddana | Weighted sum-rate maximization in wireless networks: A review[END_REF], [START_REF] Tan | Wireless network optimization by perron-frobenius theory[END_REF] is to find a power vector p such that the quality of service constraints hold:

r i (p) ≥ r * i , ∀ i (2)
where inequality is component-wise. Finding a solution to this problem assumes that such a power vector exists; to characterize the existence of such power vectors, we shall define the notion of channel feasibility:

Definition 1.
A channel, given by G, η is said to be feasible with respect to a target SINR vector r * = (r * 1 , . . . , r * N ) if there exists a p satisfying (2). The channel is otherwise said to be infeasible.

We can easily check for channel feasibility by constructing the re-weighted noise vector γ and the re-weighted gain matrix W:

γ i = r * i η i G ij W ij = 0 , i = j r * i Gii Gij , i = j. (3)
Based on [START_REF] Foschini | A simple distributed autonomous power control algorithm and its convergence[END_REF], we have the following definition, which gives us a convenient way to check channel feasibility: Definition 2. A channel {G, η} is feasible with respect to a target SINR vector r * if and only if the largest eigenvalue λ max (W) < 1.

We also define the optimal power vector p * ∈ P to be the smallest component-wise p that satisfies [START_REF] Rappaport | Wireless communications: principles and practice[END_REF] and lies in the feasible support set.

B. Power Control Algorithms

The standard algorithm used in distributed power control systems is the Foschini Miljanic (FM) algorithm [START_REF] Foschini | A simple distributed autonomous power control algorithm and its convergence[END_REF]. The FM power update rule, modified slightly to accommodate the maximum power bounds, is

p t+1 i = min p t i r * i r i (p t ) , p max i (4)
A novel alternative update to the FM algorithm is the dual averaging (DA) power control algorithm [START_REF] Zhou | Stable power control in wireless networks via dual averaging[END_REF], which was inspired by the dual averaging/mirror descent optimization algorithm [START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF], [START_REF] Zhou | Stochastic mirror descent in variationally coherent optimization problems[END_REF], [START_REF] Mertikopoulos | Learning in games with continuous action sets and unknown payoff functions[END_REF] . This algorithm performs updates on a dual power vector, y, and then transforms these dual power vectors into the transmitted power vector p by finding the closest p to y that satisfies the constraints given in [START_REF] Rappaport | Wireless communications: principles and practice[END_REF] and is feasible. The dual averaging update equation is:

y t+1 i = y t i - 1 t + 1 G ii p t i -r * i ( j =i G ij p t j + η i ) , (5)
and these dual power vectors are simply clipped to ensure feasibility, i.e. p t i = max(0, min(y t i , p max i )).

C. Delay in Wireless Networks

While both the FM algorithm and the DA algorithm enjoy optimal convergence properties in a delay-free network, realworld networks include multiple sources of latency. As they stand, neither the FM update presented in Equation 4nor the DA update shown in Equation 5 is able to accommodate delays; both of these update equations assume the SINR r i (p t-1 ) from the previous iteration, measured at receiver i, is immediately available to the transmitter at time step t. Furthermore, the power updates are assumed to be synchronized across all links in the system.

We aim to explore the properties of the dual averaging algorithm in the presence of delayed SINR measurements. However, before we present our Robust Feedback Averaging algorithm, we must formally introduce the notion of delay in our model.

We now consider a channel in which, for each time step and for each link, the SINR signal is delayed by D t i ∈ Z + units. Notice that these delay variables are distinct for every transmitter i and every time step t. Therefore, if transmitter i sends a power vector (p t i ) at time t , that transmitter will receive the SINR measurement r i (p t ) from receiver i at time t + D t i . This is a more comprehensive model of delay than those previously considered in the context of the FM algorithm [START_REF] Charalambous | On the stability of the foschini-miljanic algorithm with time-delays[END_REF], [START_REF] Charalambous | On the stability of a power control algorithm for wireless networks in the presence of time-varying delays[END_REF], [START_REF] Campos-Delgado | Characterization of stability margins of the foschini-miljanic power allocation strategy under constant and time-varying delays[END_REF], [START_REF] Mitra | An asynchronous distributed algorithm for power control in cellular radio systems[END_REF]. In these works, the delay was either constant or varying with respect to time, but not also varying with respect to each receiver. These conveniently avoid the aforementioned issues we have highlighted in this paper: simultaneous SINRs and out-of-order SINRs.

The notation presented here allows for a general model for delay, and can account for the both of these cases. Additionally, by making D t i a random variable, we can explore the effect of random delays on our algorithm.

III. ROBUST FEEDBACK AVERAGING WITH DELAYED

AND ASYNCHRONOUS UPDATES We now discuss the Robust Feedback Averaging algorithm. In the traditional formulation of both the FM algorithm and the DA algorithm, the SINR r i (p t ), is assumed to be immediately available to the transmitter at every time step, and updates are assumed to be synchronized across all links in the system. By adjusting the DA algorithm for delayed, asynchronous updates, we can relax these two assumptions.

The RFA algorithm, shown in Algorithm 1, is able to accommodate these delayed SINR measurements. Since the delays D t i can all be different, it is possible for a link to receive no signals at all on a given time step. Conversely, it is also possible for a link to receive multiple SINR measurements in a given time step. When no new SINR signals are received by transmitter i, the value of y i (and the transmitted power p i ) do not change. When multiple SINR measurements are received, they are simply averaged.

To perform these two actions, we calculate n, which represents the number of SINR measurements received by a transmitter in a given time step. If the transmitter has received any measurements, then an average is taken and the dual power vector is updated as shown in Line 7 of Algorithm 1. Noticeably, the algorithm doesn't know which SINRs were received, nor if the SINRs arrive out of order. Because of the running average implicit in the dual power vectors y i , the algorithm can dynamically adapt to whatever SINR measurements it receives, regardless of the order. This implicit storing of all historical power and update information is a necessary practical detail in the context of memory-constrained wireless transmitters.

In addition to the minimal memory requirements, the RFA algorithm as presented satisfies a fundamental requirement: it converges to the optimal power vector p * in the presence of random delays.

Theorem 1. Let p t be the iterates generated by Algorithm 1. If the random delays are all bounded: D t i < ∞, ∀i, t, then p t → p * , a.s..

The proof is omitted due to space limitation. The property of almost sure convergence to the optimal vector makes the RFA algorithm a powerful tool in delay-ridden networks. Furthermore, the only condition required of the proof is that the delays be finite. In real-world systems, this is always true, so this is an extremely weak assumption. And, in the presence of no delay, RFA reduces to the dual averaging algorithm and enjoys almost sure convergence in both deterministic and stochastic time-varying channels.

IV. SIMULATION RESULTS

To validate the convergence and performance guarantees of Robust Feedback Averaging, we provide several numerical experiments that illustrate the algorithm's various properties. We consider two sets of experiments in this section: first, we introduce two different forms of delay into a deterministic, time invariant channel; second, we examine the effect of adding delays in a more complex stochastic and time-varying channel.

Algorithm 1 Robust Feedback Averaging Algorithm 1: Each link i arbitrarily chooses an initial y 0 i . 2: for t = 0, 1, 2, . . . do

3:

for i = 1, . . . , N do 4:

p t i = max(0, min(y t i , p max i )) 5: n = t t =0 I t + D t i = t 6:
if n > 0 then 7: end for 12: end for

y t+1 i =y t i -1 n(t+1) t t =0 G ii p t i - r * i j =i G ij p t j + η i I t + D t i = t 8:

A. Foschini-Miljanic Algorithm with Asynchronous Delay

As an algorithmic comparison, we also alter the Foschini-Miljanic algorithm to account for transmitter-distinct, timevarying delays. We further modify the FM update presented in [START_REF] Charalambous | On the stability of the foschini-miljanic algorithm with time-delays[END_REF] to allow the FM algorithm to handle delays which are specific for each link. Our modification is similar to the modification in Algorithm 1: all of the SINR values returned in a given time step are averaged, and the average value is used for the FM update. The modified FM update rule is:

p t+1 i = 1 n t t =0 p t i r * i r i (p t ) I(t + D t i = t) (6)
where n = t t =0 I(t + D t i = t) as before. As expected, our modified FM algorithm reduces to the original FM algorithm when D t i = 0 ∀ i, t. We compare the performance of this modified FM algorithm and our RFA algorithm in the following experiments.

B. Deterministic and Time-Invariant Channel

We begin with a deterministic, time-invariant channel with 4 links. The channel properties G, η and optimal power vector p * are

G =     3 1 2 1 1 3 1 2 2 1 3 1 1 2 1 3     η =     0.1 0.2 0.3 0.1     p * =     0.079167 0.089583 0.104167 0.077083    
For these simulations, we select the target SINR r * i = 0.5 and the maximum power p max i = 1.0 for all links in the channel. We chose these parameters to present a relatively challenging environment for the algorithms; as the G matrix indicates, each link experiences a lot of interference from every other link in the system. Due to this challenging environment, we only require the links to achieve a SINR of 0.5. The maximum power p max i = 1.0 was chosen by considering the optimal power vector p * ; any algorithm which finds the optimal transmitter powers will not be affected by this upper bound in the long run.

We can calculate the re-weighted gain matrix and reweighted noise vector W, γ, and verify that λ max (W) = 0.67 < 1. Therefore, based on Definition 2, the channel is stable. This means that, for both the FM algorithm and the RFA algorithm, the transmitted link powers will converge to the optimal values in the absence of delay.

We will examine algorithm performance on this channel in several different delay regimes. For the simulations shown, all powers p 0 i (FM) and dual powers y 0 i (RFA) were initialized to zero.

Experiment 1: We first examine the algorithmic performance under the regime of constant delay. However, we make the delay distinct for each channel. For this case, we set D = [D 1 D 2 D 3 D 4 ] T = [25 50 75 100] T ∀t. This is a reasonable experiment, as in a practical setting, each transmitter-receiver pair could have different computational capabilities or could be separated by a greater physical distance, causing the delays for each link to be distinct. Both algorithms converge to the optimal power vector, and RFA uses less overall power.

Even this simple scenario is not covered by the works which studied the FM algorithm in the presence of delay [START_REF] Charalambous | On the stability of the foschini-miljanic algorithm with time-delays[END_REF], [START_REF] Charalambous | On the stability of a power control algorithm for wireless networks in the presence of time-varying delays[END_REF], [START_REF] Campos-Delgado | Characterization of stability margins of the foschini-miljanic power allocation strategy under constant and time-varying delays[END_REF], [START_REF] Mitra | An asynchronous distributed algorithm for power control in cellular radio systems[END_REF], so the FM algorithm is not guaranteed to converge. However, as we proved in Theorem 1, RFA will eventually return the optimal power vector p * . Figure 1 shows the power for each link with this constant delay. Specifically, Figure 1(a) shows the transmitted power for each link plotted against time when using the RFA algorithm, and Figure 1(b) does so using the modified FM algorithm. As expected, neither algorithm updates until the 25 th time step, when the first SINR is received. At t = 25, only Link 1 updates, as it is the only one to receive a SINR measurement; 25 time steps later, Link 2 begins updating, then Link 3, and finally Link 4.

The FM algorithm exhibits interesting, step-wise behavior in Figure 1(b). This is due to the fact that the FM algorithm only depends on the most recent SINR measurement to set its next transmission power. Because the transmission powers at time t = [0, 25) were all the same, the SINRs received during t = [25, 50) were all the same, so the transmission powers in that interval remained constant.

Regardless of the stepwise behavior, the power appears to converge to the optimal values both while using the FM algorithm and while using the RFA algorithm. however, the transient behavior is quite different. Noticeably, the FM algorithm encourages all of the transmitted powers to rise to the maximum power at the beginning of the simulation. This is an effect of the FM algorithm only considering the most recent SINR measurement when performing an update; at the beginning, all SINRs are below the target values, so they are increased drastically. It then takes a considerable amount of time (200 time steps) for the values to eventually stabilize, and it takes another 400 time steps for all of the powers to simultaneously decrease to eventually settle on the optimal powers.

Meanwhile, the RFA algorithm encourages a more gradual approach. Because of the role of averaging, when a transmitter receives a SINR measurement that is far below the target SINR, it slowly increases its transmitted power. With all transmitters doing this in tandem, this causes the trajectories of the transmitted powers to approach the optimal power vector from below, rather than from above. This is significant, as the area under each curve represents the total power used by each transmitter. As seen in Figure 1, there is about an order of magnitude difference in the area under the curves in Figure 1(a) and Figure 1(b), which represents an order of magnitude power savings caused by RFA.

Experiment 2: The advantage of RFA is even more apparent when we allow the delays to be random. In our second simulation, we sample each delay D t i individually from a uniform distribution over the integers, i.e. D t i ∼ U (0, 100). This is a more challenging environment than in the previous experiment, because the delay in the previous experiment, though it varied across channels, was constant with respect to a single channel. That meant that (after t = 100) every transmitter received a single SINR measurement. In this experiment, we examine the impact of multiple SINR measurements per time step, and out-of-order SINR measurements.

Figure 2 shows the transmission power of each transmitter vs time. In Figure 2(a), we see a much more noisy convergence curve for the FM algorithm, with frequent spikes in the transmitted powers. Each spike in the plot after the initial transient behavior is due to an out-of-order SINR measurement; here, again, the "memoryless" property of the FM algorithm hinders its performance. The initial spike in power, which is similar to the behavior observed in Experiment 1, is again an unnecessary waste of power.

The RFA algorithm, on the other hand, does not cause any of the transmitters to send spikes in power when an out-of-order SINR is received. Similarly to Experiment 1, it converges from below to the optimal power vector in an energy-efficient manner.

C. Stochastic and Time-Varying Channel

We now turn to the case of a stochastic, time-varying channel with delay. We define a stochastic channel as one in which the gain matrices and noise evolve stochastically over time: {G t } ∞ t=0 , {η t } ∞ t=0 . For these experiments, we consider a stochastic channel that is drawn i.i.d. from a set of 2 deterministic channels, each with equal probability:

(G t , η t ) = (G 1 , η 1 ), w.p. 0.5 (G 2 , η 2 ), w.p. 0.5, G 1 = 2 5 6 3 η 1 = 0.1 0.2 G 2 = 5 3 3 3 η 2 = 0.15 0.05
Again, we chose these matrices to present a challenging environment for the algorithms; even though there are only two channels, the G matrices dictate that the signal in each link will heavily interfere with the signals from the other link.

Experiment 3: First, we set r * i = 0.5 and p max = 1.0 as before. Using this r * i , we check channel feasibility by constructing W 1 , γ 1 , W 2 , and γ 2 . We see that λ max (W 1 ) = 1.34 and λ max (W 2 ) = 0.46. This indicates that the channel G 2 , η 2 is feasible (by Definition 2, but the channel G 1 , η 1 is infeasible. This means that, half of the time, the algorithms will be operating on an infeasible channel; the FM algorithm will not converge to the optimal power vector in this scenario.

However, the average channel, G = 3.5 4 4.5 3 η = 0.125 0.125 , is feasible (with λ max (W) = 0.79). So, by the results shown in [START_REF] Zhou | Stable power control in wireless networks via dual averaging[END_REF], the RFA algorithm allows all the link powers in this channel to converge to their optimal value (with no delay). Now, we add delay to this stochastic channel. We again add uniform independent delay for each channel for each time step, i.e. D t i ∼ U (0, 100). Figure 3 shows the results of this simulation, with Figures 3(a) and 3(b) illustrating the performance of RFA and FM, respectively.

We see that, again, the RFA algorithm converges smoothly to the optimal power vector from below, just as in previous experiments. However, the FM algorithm, as expected, does not converge. Since FM only relies on the most recently returned SINR measurement, whenever the channel is infeasible (which occurs in this experiment with probability 0.5), the transmitted powers will diverge. Of course, when the channel is feasible, the power vector decreases, so the power will only reach the maximum value if a long sequence of SINR measurements from the infeasible channel were returned sequentially, which will only happen with very low probability.

Experiment 4: For this experiment, we did not examine the performance of FM, but instead focused on the performance of RFA. We use the same stochastic channel as in Experiment 3, but here we vary r * i between 0.5 and 0.9 (where we always set r * 1 = r * 2 ). As in Experiment 3, setting r * i = 0.5 causes the average channel, G = [ 3. 5 4 4.5 3 ] η = [ 0.125 0.125 ] to be feasible (with λ max (W) = 0.79). However, when r * i = 0.9, we reconstruct W i and γ i and see that λ max (W) = 1.31 in this case. Therefore, the channel is on average infeasible, and we do not expect the link powers to converge when using RFA.

Figure 4 shows the results of this experiment. Unlike in the other figures, here Figure 4 i . As we expect, when the convergence condition in Definition 2 holds, then the power vector converges to the optimal p * . But, when λ max (W) > 1, the powers diverge, and the powers diverge faster with larger λ max (W).

Fig. 1 :

 1 Fig. 1: Results of Experiment 1: power marginal convergence behavior of the RFA algorithm (a) and the FM algorithm (b) with a fixed, constant delay for each channel: D = [25 50 75 100] T ∀t for channels 1-4, respectively.Both algorithms converge to the optimal power vector, and RFA uses less overall power.

Fig. 2 :

 2 Fig. 2: Results of Experiment 2: Power marginal convergence behavior of the RFA algorithm (a) and FM algorithm (b) with a uniform random delay of [0, 100] time steps. Due to the out-of-order nature of the SINR measurements, the FM algorithm frequently unnecessarily increases the transmission power. Both algorithms appear to converge.

Fig. 3 :

 3 Fig.3: Results of Experiment 3: power marginal convergence behavior in the presence of a stochastic, timevarying channel and uniformly random delay. The delays are sampled from uniform distribution: D t i ∼ U (0, 100). Again, RFA (a) converges to the optimal average power, but the FM algorithm (b) does not converge, since the channel is infeasible half of the time.

Fig. 4 :

 4 Fig. 4: Results of Experiment 4: power marginal convergence behavior of the RFA algorithm in the presence of a stochastic, time-varying channel and uniformly random delay (D t i ∼ U (0, 100)). Five different values of r * i are plotted here, and the corresponding λ max (W) for each is shown. The two links in the channel are plotted separately, with the transmitted power of Link 1 (a) above, and the transmitted power of Link 2 (b) below.

  (a) shows the power transmitted by Link 1, and Figure 4(b) shows the power transmitted by Link 2. Each individual trace represents a different target SINR r *