Nazim Fatès
email: nazim.fates@loria.fr

Collective infotaxis with reactive amoebae: a note on a simple bio-inspired mechanism

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The purpose of this note is to study how to coordinate a team of agents to perform a collective search in a two-dimensional environment. The objective of the team is to locate a source which emits an information that can only be detected sporadically. This problem arises in various situations studied in Biology, for instance when insects emit pheromones to attract their partners.

There are cases where it is impossible to nd the partner by following a gradient because the pheromones are dispersed in a turbulent environment: they are found in patches rather than with a smooth distribution. Vergassola et al. named infotaxis a search strategy which consists in maximizing the information gain on the location of the source [START_REF] Vergassola | Infotaxis' as a strategy for searching without gradients[END_REF].

Here, we study this problem in the context of cellular automata: imagine a team of agents that are randomly dispersed on a two-dimensional grid and need to locate a particular cell, the emitting source. The agents have a very small memory; they have no map of their environment nor can they build one.

They can only perceive the state of the cells situated at short distance. At each time step, the agents may detect the presence of the source (make a 'hit'), but without knowing the direction of the source, nor its distance. The detection probability increases as one gets closer the source and this is the only "hint" that can guide the agents. At each time step, the agents can move or change the state of the cell they are located on. Under such circumstances, how can a team of agents achieve a form of coordination in order to eciently locate the source?

One may call this problem the distributed infotaxis problem with memoryless agents. In contrast, (classical) infotaxis has been mainly studied with agents that remember where the hits occurred. This allows them to build a map of their environment according to the detections they have made. Since the detections are stochastic events, the map is progressively adjusted. Interestingly, it was shown that a good search strategy does not necessarily consist in always trying to approach the source. In some cases, it may be more interesting to move orthogonally to the estimated direction of the source in order to improve the process of locating the direction of source [START_REF] Vergassola | Infotaxis' as a strategy for searching without gradients[END_REF]. However, the use of a map in virtual agents or in real robots imposes a certain degree of complication.

Dierent solutions have been proposed to simplify this scheme. In particular, Masson has proposed to use only a projection of the map of probabilities into a pre-determined form, where appears detection and non-detection terms [START_REF] Masson | Olfactory searches with limited space perception[END_REF].

Later on, this technique was extended by Zhang et al. to coordinate a team of robots [START_REF] Zhang | Multi-robot searching with sparse binary cues and limited space perception[END_REF].

We ask if we can simplify the solution further and if we can solve this problem with agents that do not remember where previous detections have occurred. In this context, the previous solutions cannot be used as the only thing that the agents can do is to use their respective positions on the grid to make better movements than they would do with a pure random walk.

The solution we propose is a variation of a previous model, the social amoebae [START_REF] Fatès | Solving the decentralised gathering problem with a reactiondiusion-chemotaxis scheme -social amoebae as a source of inspiration[END_REF], which was proposed to solve the decentralised aggregation problem: the agents need to gather and form a compact group in a totally decentralised way with only local perceptions and no memory of their previous actions. The aggregation scheme imitates the aggregation behaviour which occurs in the rst phase of grouping in the social amoeba Dictyostelim discoideum. The mechanism heavily relies on the asynchronous emission of reaction-diusion waves in the environment: each (virtual) amoebae needs to alternatively speak and listen and these dierent roles are governed by simple random events (see more details in Ref. [START_REF] Fatès | Solving the decentralised gathering problem with a reactiondiusion-chemotaxis scheme -social amoebae as a source of inspiration[END_REF]). If they all speak, (that is, emit waves) synchronously the aggregation process cannot occur.

We now modify this aggregation scheme to make the amoebae search for a hidden source. Informally, the system works as follows: each agent performs a random walk on the grid. When it makes a detection, which should be a rare event, it triggers a reaction-diusion wave in the grid. This wave will propagate and attract the other agents towards the emitting agent. The agents that get closer to the source will have a tendency to emit more often that the agents that are far from it, and the group will collectively have a tendency to move closer to the source. Little by little, the group will progress towards the source and nally reach its target. The coordination of the group is ensured by the asynchronous nature of interactions: there are no leaders to take the decisions but instead each agent has the possibility to emit waves or to "listen" to the emissions of the others according to some particular stochastic behaviour that needs to be adjusted appropriately.

It can be noted that although seducing, the success of this scenario is not guaranteed. Indeed, this indirect communication process between agents, also called stigmergy, has to occur inside a group that is neither too compact nor too sparse. Indeed, an excessive grouping of agents is detrimental because the space they explore with their random movements is more limited than if the agents were absolutely free to wander on the grid. On the contrary, if the agents stay away from each other, this means that the attractive force of the reaction-diusion waves is too weak to counter-balance the random movements of the agents. In this note, we explore this trade-o with a rst set of numerical simulations. We show that the use of the virtual amoebae search scheme can be benecial in some particular conditions. Since our experiments are limited, we will mainly aim at showing that this search strategy is worth of consideration.

We do not bring any denitive answer to the problem of distributed infotaxis but rather ask questions that will need to be explored with a deeper analysis.

Presentation of the model

The space on which our model takes is dened a nite square grid denoted by L = {1, . . . , L}×{1, . . . , L}. Our model, the we could name the amobae-infotaxis model, is composed of two interacting layers: the environment and the amoebae.

The environment is a cellular automaton composed of the cells of the grid.

The state of the cells are taken in a set Q that will be specied below. The global state of the environment is called a conguration and the set of congurations is denoted by Q L . At each time step, each cell calculates its new state according to its own state and the states of its neighbours and a possible inuence from the amoebae. In other words, the environment is a cellular automaton which is in interaction with a system of reactive agents.

The emitting source is modelled by a distinguished cell, whose location σ ∈ L is xed over the time. As mentioned earlier, in this problem, the cells can only make binary detections and the probability of detection decreases with the distance to the emitting source. This property is modelled by the introduction of a probability detection eld ρ : L → [0, 1] which associates to each cell the probability to make a detection in this given cell. Here, we work with a eld which follows an inverse square law ; it is described with: ρ(σ) = 1 and

∀c = σ, ρ(c) = K d(c, σ) 2 , ( 1
)
where d(x, y) = |x 1 -y 1 | + |x 2 -y 2 |
is the Manhattan distance between two cells of the grid x = (x 1 , x 2 ) and y = (y 1 , y 2 ), and where K < 1 is an arbitrary constant. Other authors have used an exponential decrease with the distance to the source. We prefer to use a 1/d 2 function in order to allow the existence of a small, but non neglible, probability of detection in the cells that are far from the source. Note that in our view, this eld does not make the problem non-local: the probabilities are simply some xed parameters that cannot be "read" directly by the cells or by the agents. Also note that this eld may be computed by a pure cellular automaton if one wishes to do so.

The amoebae are reactive agents: they are described by their positions on the grid and their internal state, which is binary and taken in the set I = {rec, emi} (emitter and receiver). For a set of N amoebae, we will denote by p = (p i ) ∈ L N the vector of all positions on the grid and by s = (s i ) ∈ S N the vector all internal states. The quantity p i and s i thus denote the position and the internal state of the i-th amoeba, respectively.

With these notations, we can describe our model with a dynamical system F , which operates on the space Q L × L N × S N . For an initial condition (x 0 , p 0 , s 0 ) ∈ Q L × L N × S N , its evolution is described by the recursive equation: (x t+1 , p t+1 , s t+1 ) = F (x t , p t , s t ). Let us now describe how these quantities are calculated 1 .

Evolution of the cells

The set of states of the cells is given by Q = {e, r, n}, which respectively correspond to the excited, refractory and neutral states. This environment implements a simple reaction-diusion law, called the Greenberg-Hastings model(see e.g. Ref. [START_REF] Berry | Robustness of the critical behaviour in the stochastic Greenberg-Hastings cellular automaton model[END_REF]). Informally, a cell becomes excited if it is neutral and if it has at least one excited neighbour or if it is excited externally by an amoeba. An excited cell always becomes refractory and a refractory always become neutral.

1 We make a slight abuse in notations because we use the formalism of classical dynamical systems even though our function F is stochastic.

4

To dene this law, we introduce the function S e that selects the excited cells around a given cell: S e (x, c) = {d ∈ D, c + d ∈ L, x c+d = e} where D = {(0, 1), (0, -1), (1, 0), (-1, 0)} represents the four cardinal directions. This means that our cellular automaton uses the von Neumann neighbourhood but note that it can also be easily applied to other types of neighbourhoods [START_REF] Fatès | Solving the decentralised gathering problem with a reactiondiusion-chemotaxis scheme -social amoebae as a source of inspiration[END_REF].

We then dene the modied Greenberg-Hastings law of evolution of the cells with, for all c ∈ L:

x t+1 c =          e if x t

Evolution of the amoebae

First, we need to describe the fact that some amoebae will make a hit, that is, detect the source with a given probability which depends on their location. To this end we introduce the function H, which, given a set of positions p ∈ L N , tells whether the i-th amoeba will make a 'hit': ∀p ∈ L N , ∀i ∈ {1, . . . , N }, H(p, i) = 1 with probability ρ(p i ) 0 with probability 1 -ρ(p i ).

The behaviour of the amoebae is as follows: at each time step, independently of the other amoebae, each amoeba may follow a random movement according to the four cardinal directions with probability p A , the probability of agitation.

With probability 1 -p A , the "regular" behaviour is applied. This regular behaviour consists in the following rules:

1. if an amoeba is in the state rec, if it is standing on a neutral cell and sees one or several excited cells in its neighbourhood, it chooses one of these excited cells uniformly at random for its movement ; otherwise, if it has no excited neighbour, its stays on the same cell.

2. if an amoeba is in the state emi, if it is standing on a neutral cell, it turns this cell to the excited state. It then remains in state emi with probability p R , the probability of reemission.

The motivation for introducing the probability of agitation p A is twofold:

on the one hand, we add a randomness which ensures that the amoebae are not too static, on the other hand, this allows us to model the intrinsic noise that would exist in the movement of biological organisms or robots. The use of the probability of reemission p R is not mandatory. It is simply here to allow several reaction-diusion waves to be triggered when a hit occurs. As we will see below, these waves have an attractive force that groups the amoebae; this attractive force balances the dispersive force that is created by the random movements.

To write down this rule with our notations, we have:

• for s t i = rec, s t+1 i = emi if H(p t , i) = 1 and s t i = rec rec otherwise , • for s t i = emi, s t+1 i = emi with proba. p R rec with proba. 1 -p R .
The positions evolve according to:

p t+1 i = B p t i , U(D) with proba. p A U S e (x t , p t i ) with proba. 1 -p A ,
where

• U(X) is a random variable that selects an element uniformly in a set X,

• B[c, d] is a function that equals c + d if c + d is in L and c otherwise,
• S e is the set of excited cells neighbours to a given cell (see above).

The task is completed when an amoeba moves to the emitting source σ. The search time is the smallest t which satises: ∃i, p t i = σ.

First experiments

Recall that our purpose is simply to test the relevance of doing a collective form of infotaxis with (almost) memoryless agents is feasible. For the sake of simplicity, we arbitrarily x the number of amoebae to N = 20 and the value of the constant in Eq. 1 to K = 0.5. Similarly, we x the reemission probability to p R = 0.9 in order to allow for multiple waves to be emitted when a 'hit' occurs.

A rst visual experiment will be presented for L = 20 and then statistics will be given for L = 80.

Qualitative observations

In our experiments, the amoebae are placed uniformly at random on the grid.

The emitting source is placed at coordinates (L/4, L/4) in order to make the problem more dicult. Indeed, if it were placed at the centre, the amoebae could more easily nd it since they do not need to displace the centre of gravity of the group.

First, let us observe visually how a system with only 5 amoebae behaves for a small grid of size L = 20 and a small probability of agitation p A = 0.1.

Figure 1 shows one particular random evolution of such a system. We see that various reaction-diusion waves are emitted at dierent times when the hits occur. In this evolution two hits occur in a short time around t = 20, which allows us to observe how the reaction-diusion annihilate when they meet frontally and how they merge when they meet orthogonally. This is the reason why it is important that the emission of waves do not occur simultaneously:

t = 0 t = 10 t = 15 t = 20 t = 30 t = 35 t = 100 t = 150
as the reaction-diusion waves annihilate when they meet, if all the amoebae were emitting at the same, none of them could receive the information. In other words, some amoebae need to remain in the rec state in order to be able to feel the waves emitted by the other amoebae. In this evolution, after a hundred steps the amoebae are already regrouped around the emitting source and at time t = 150, we are very close to the source detection event.

Quantitative experiments

We now make a second step by considering the statistical variations of the search time as a function of the agitation probability p A . Indeed, p A is the most important parameter in our study since it controls the average quantity of movement of the amoebae.

Taking the same settings as in the previous experiment, we observe the variation of the average detection time for a grid size of L = 80. The agitation probability p A varies from 0.1 to 1.0 with a step of 0.1 and the average search time is obtained with 1000 samples.

The results are shown on Fig. 2 p A = 0.15 to have a better observation of this zone.)

We believe that the existence of this minimum corresponds to a balance between dispersive eect of the agitation and the grouping eect the reactiondiusion waves. Indeed, if p A is too small, the amoebae have a tendency to all aggregate on the same point not the source ! and to all remain almost static close to this point. As the random movements are rare, if we imagine an amoeba that progressively gets far from the aggregation point, several emissions of waves emitted by the amoebae will call it back to this point. On the contrary, when p A is too high, the disorder introduced by the agitation is no longer benecial because the information gained by listening to the other amoebae is diluted by the random movements.

For values of p A which are even larger, that is, greater than 0.75, we observe that the detection time decreases again. This second eect is simply explained by considering the average number of movements per step: since the order of the system is already very weak (as p A is high), it becomes more ecient to mechanically increase the average speed of the amoebae by moving more often.

In order to estimate what part of the behaviour is due to the reactiondiusion waves and what part is due to the random movements, we estimated the average search time when there are no waves (we set K = 0). The corresponding curve on Fig. 2 shows no surprise: it is decreasing with p A and, more precisely, the average detection time varies as 1/p A . This is expected since the amoebae follow a pure Brownian motion and do not interact. We can observe that for small values of p A , the gain is important and for a similar average speed of movements, the search time with the amoebae-infotaxis search strategy may be divided a factor 5 with respect to the random walk strategy.

Discussion

This note presented a discrete model for performing a collective infotaxis task in an environment where only an intermittent information is available. The rst experiments show that the use of reaction-diusion waves is a simple means to perform a collective search of a hidden source. The mechanism we described relies on the probabilistic and asynchronous nature of the detections: since these events are rare, the agents of the system (the amoebae) can alternatively act as emitters or receivers and thus coordinate their movements to get closer and closer to the emitting source.

A continuation of the experimental work is needed to get more insights on how to use virtual amoebae to solve the decentralised infotaxis problem. Indeed, it is not clear yet in which situations it is more ecient to use this mechanism rather than a simple random walk. It is also important to get a more precise view on the scaling laws of the system: how does the search time depend on the number of amoebae and on the size of the grid? There are also questions which regard the possible improvement of the model. Here we presented a solution that uses only three states for the cellular automaton dynamics of the environment and only two states for the agents. The possibilities for further extending the behaviour of the agents are numerous, as for instance, making their emission auto-adjusting in order to have the correct balance between listening and speaking.

In the context of applying these elements of research in the context of distributed robotics [START_REF] Fatès | A robust scheme for aggregating quasi-blind robots in an active environment[END_REF][START_REF] Silva | Modelling synchronisation in multirobot systems with cellular automata: Analysis of update methods and topology perturbations[END_REF], it would also be interesting to test the robustness of our scheme against several perturbations such as the introduction of obstacles in the grid or a non-regular distribution of the detection probabilities.

Figure 1 :

 1 Figure 1: Evolution of the system for N = 5 amoebae on a grid of size L = 20 with K = 0.5, p A = 0.1, and p R = 0.9; the emitting source is in black (diamond shape); amoebae are shown as circles in grey (rec state) or in yellow (emi state). The colors of the grid cells are red, orange or white, which respectively correspond to the e, r and n states. [color online]

Figure 2 :

 2 Figure 2: Search time for N = 20 amoebae as a function of the agitation probability p A for L = 80 (in blue) with K = 0.5, p A = 0.1, and p R = 0.9.The curve in purple shows the detection time when the amoebae follow a pure random walk (reaction-diusion waves are suppressed).

  and, interestingly, present a rather surprising behaviour of the system. It can be noted that a local minimum on the detection time is obtained for p A ∼ 0.2. (An additional sampling point is added at

	500 1000 1500 2000 2500 avr. detection time				L=80 N=20 K=0.5		pure random walk with RD waves
	0.1 0	0.2	0.3	0.4	0.5	pA	0.6	0.7	0.8	0.9	1.0


Rodolphe Jenatton
email: jenatton@amazon.com

Jim C Huang
email: huangjim@amazon.com

Dominik Csiba
email: cdominik@gmail.com

Cedric Archambeau
email: cedrica@amazon.com

Online optimization and regret guarantees for non-additive long-term constraints

Introduction

Online optimization can be viewed as a sequential game where in each round t ∈ {1, • • • , T }, we are required to play an action, represented by a vector x t , which takes values in a set of actions {X t : X t ⊆ R d } T t=1 . We then observe a reward f t (x t ) as a function of the action we chose in round t. The goal is to generate a sequence of actions such that some measure of performance is for instance maximized over the course of T rounds.

The performance metric typically adopted in the online learning framework is the cumulative regret Cesa-Bianchi and Lugosi (2006). Moreover, online learning traditionally assumes that the objective function has an additive structure that nicely decomposes as a sum of regrets over the rounds of the optimization. This means that there is no coupling of x t 's across successive rounds when the sets X t are decoupled across time. Hence, the challenge in online learning rather lies in the fact that x t must be estimated before having access to the reward function f t , which corresponds to the so-called 0-lookahead setting [START_REF] Buchbinder | Unified algorithms for online learning and competitive analysis[END_REF]; [START_REF] Andrew | A tale of two metrics: simultaneous bounds on competitiveness and regret[END_REF]. In the particular case of dynamic regret with additive objective functions [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF]; [START_REF] Cesa-Bianchi | A new look at shifting regret[END_REF]; [START_REF] Hall | Dynamical models and tracking regret in online convex programming[END_REF]; [START_REF] Jadbabaie | Online optimization: Competing with dynamic comparators[END_REF], one seeks to analyze how, under various assumptions on f t and X t , the sequence of rewards f t (x t ) collected in an online fashion compares with the best sequence of rewards f t (x t ) collected in hindsight. The dynamic regret over T rounds is defined as follows:

1 T T t=1 f t (x t ) - 1 T T t=1 f t (x t ), (1)
with 1 T T t=1 f t (x t )=max x1∈X1,...,x T ∈X T 1 T T t=1 f t (x t ). Unlike static regret analysis [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF], where we compare against a single best action x in retrospect over all rounds played, the best dynamic comparator, as its name indicates, need not be identical for all t ∈ {1, • • • , T }.

In this work, we consider a setting similar to online learning with dynamic regret, but we do not assume that successive actions x t are decoupled. We provide novel regret bounds for the case where rewards are additive and the regrets are non-additive as the total cost over T rounds is non-decomposable. There is a small body of recent work that study non-additive regrets. For example, [START_REF] Rakhlin | Online learning: Beyond regret[END_REF] consider a wide class of non-additively decomposable objective functions in a 0-lookahead setting, covering for instance the problems of Blackwell's approachability, the calibration of forecasters and the global cost online learning game from [START_REF] Even-Dar | Online learning for global cost functions[END_REF]. More recently, [START_REF] Kar | Online and stochastic gradient methods for non-decomposable loss functions[END_REF] handle some specific form of non-decomposability in relation with the online optimization of metrics such as the precision at k. Closer in spirit to our work is the one by [START_REF] Agrawal | Fast algorithms for online stochastic convex programming[END_REF], who provide expected regret bounds in the case where non-additive costs depend only on the empirical mean of the actions.

Metrical task systems (MTS) [START_REF] Borodin | An optimal on-line algorithm for metrical task system[END_REF] provide an alternative analysis framework for online optimization. They consider movement costs that penalize variations of x t 's across time and rely on competitive analysis [START_REF] Borodin | Online computation and competitive analysis[END_REF]. Similar to online learning with dynamic regret, the performance of the online optimization algorithm is compared to the best sequence of actions. However, instead of measuring cumulative regret, competitive analysis adopts a multiplicative metric, known as the competitive ratio. We refer the interested reader to [START_REF] Buchbinder | Unified algorithms for online learning and competitive analysis[END_REF]; [START_REF] Andrew | A tale of two metrics: simultaneous bounds on competitiveness and regret[END_REF] for a detailed discussion. Moreover, the online game in MTS differs from the online learning setup in that it follows the 1-lookahead setting where the player has access to the reward function f t before estimating x t . The typical instantiation of the movement cost in MTS is a total variation penalty T -1 t=1 x tx t+1 2 , that can also be defined with non-Euclidean norms [START_REF] Bera | Fenchel duals for drifting adversaries[END_REF]. While the movement cost introduces dependencies across actions x t in successive rounds, we note that, to the best of our knowledge, previous work only considered movement costs with an additive structure. We depart from this approach by considering non-additive penalties.

The aforementioned frameworks have been successfully applied to derive and analyze online optimization algorithms in for example online routing [START_REF] Awerbuch | Online linear optimization and adaptive routing[END_REF], process migration of servers [START_REF] Borodin | Online computation and competitive analysis[END_REF] and portfolio allocation [START_REF] Helmbold | On-line portfolio selection using multiplicative updates[END_REF]. However, the assumption that regrets or movement costs have an additive structure is restrictive in practice. In particular, we consider the problem of online ad allocation, which is at the core of modern display advertising systems. The online ad allocation [START_REF] Chen | Online bidding algorithms for performance-based display ad allocation[END_REF] problem consists of sequentially allocating ad impressions (encoded by x t ) to a large number of competing ad slots across a large number of websites and mobile apps, subject to a variety of advertiser objectives and constraints. Advertisers will typically expect that in solving the ad allocation problem, we maximize a measure of ad performance (or advertiser welfare), subject to constraints on user targeting and constraints on ad delivery (e.g., spend as close to 100% of an ad's budget as possible over T rounds). The online ad allocation problem is characterized by non-stationarity in constraints (e.g., for a budget-constrained allocation, the amount of budget consumed per round varies dynamically). Moreover, in such an application, constraints are measured by the advertiser only in a long-term sense (e.g., budget consumed by the end of an advertising campaign), which requires the use of non-additive constraints.

Hence, the online ad allocation problem can be viewed as a hybrid between online learning with dynamic regret and MTS. It falls into the 1-lookahead setting and requires non-additive constraints, but there is no practically-justifiable concept of movement cost on xt that can be applied to account for non-additive constraints. Moreover, the movement cost is not a sensible penalty to capture the fact that one would like to show different ads in successive rounds. Finally, online ad allocation is concerned with satisfying long-term constraints, which means that the cumulative constraint violations resulting from the sequence of vectors {x t } T t=1 should not exceed a certain amount by the final round T . The previous approaches proposed in the online learning literature [START_REF] Mahdavi | Trading regret for efficiency: online convex optimization with long term constraints[END_REF]; [START_REF] Agrawal | Fast algorithms for online stochastic convex programming[END_REF] are not suitable for our use-case since they handle stationary long-term constraints with static regret guarantees in the 0-lookahead setting. And again, this type of constraints is difficult to encode as movement costs.

Contributions. Within the 1-lookahead setting, we study an online optimization problem with dynamic regret defined in terms of a sum of concave rewards and a penalty that does not decompose additively over the rounds of the online game. Non-additive penalties are suitable for modeling non-stationary and/or long-term constraints, which are of practical importance in display advertising. The resulting formulation is, to the best of our knowledge, novel and extends the work of [START_REF] Mahdavi | Trading regret for efficiency: online convex optimization with long term constraints[END_REF] and [START_REF] Agrawal | Fast algorithms for online stochastic convex programming[END_REF]. We further provide dynamic cumulative regret guarantees, showing that vanishing regret is driven by the convexity and the smoothness of the non-additive penalty, along with the smoothness with which the residuals vary over the rounds. Hence, the structure of our bound echoes results derived in previous work within the context of 0-lookahead dynamic regret analysis [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF]. Finally, we illustrate our methodological contribution by conducting experiments on synthetic data to validate the benefits of non-additive penalties and study the convergence of the cumulative regret on live traffic data collected by a display advertising platform in production.

Problem statement

Motivated by the practical realities of the online ad allocation problem, we focus on a class of online optimization problems in which the actions played over time must not only maximize some rewards but also minimize longterm constraints that are penalized by some non-additive error term through the function E. Formally, for any

(x 1 , . . . , x T ) ∈ X 1 × • • • × X T ,
we study the online maximization of objective functions of the form

P(x 1 , . . . , x T ) 1 T T t=1 f t (x t ) -E 1 T T t=1 A t x t -b t , (2)
where the matrices A t ∈ A ⊆ R m×d and vectors b t ∈ R m encode the non-stationary, long-term constraints. Although the sets X t 's depend on t, they are assumed to be fully available at each round t and do not take part in the online game.

Our goal is to derive an online algorithm and prove that the sequence (x 1 , . . . , xT ) ∈ X 1 × • • • × X T it generates is guaranteed to satisfy a regret bound of the form

P -P(x 1 , . . . , xT ) ≤ R T + S e + S A , (3)
where we have defined the optimal (primal) objective P max x1∈X1,...,x T ∈X T P(x 1 , . . . , x T ) = P(x 1 , . . . , x T ).

(4)

The regret bound in (3) is comprised of three terms, each capturing different aspects of the problem. The first term R T quantifies the contribution due to the application of an online projected subgradient method in the dual problem for (2), as described in Section 3.1.

The second term S e quantifies the impact of the smoothness with which the sequences of m-dimensional error vectors {e t } T t=1 evolve over time, where e t A t x tb t . In this respect, the term S e will be reminiscent of the guarantees traditionally obtained with dynamic regret analysis [START_REF] Cesa-Bianchi | A new look at shifting regret[END_REF], with the key exception that the smoothness is not based on the variables x t themselves, but rather on the error vectors e t . This term also has a natural interpretation and relevance in the online advertising setting, as discussed in Section 4.2. The derivation for S e will be the topic of Section 3.2.

The third and last term S A models the impact of having to estimate constraint matrices A t via Ât , as in practice A t is not known at the time when the action x t is played. This aspect is also motivated by practical aspects of the online ad allocation problem. As it will be made clear in Section 3.3, S A will depend on how smoothly the constraint matrices A t 's evolve over time. Given the above, our methodological contributions can thus be viewed as both a regret analysis with dynamic comparators for non-additive objective functions, and a competitive analysis evaluated through a regret criterion -or more precisely, a competitive difference criterion, as defined in [START_REF] Andrew | A tale of two metrics: simultaneous bounds on competitiveness and regret[END_REF] -where the service cost is f t and our movement cost E is convex and non-additive. We will interchangeably refer to both competitive difference and dynamic regret in the sequel.

Analysis

We describe in this section the main components of our analysis, starting with the saddle point formulation. The proofs of the results are relegated to the Appendix.

Linearization and Lagrangian formulation

Since the main challenge of our analysis lies in the fact that E breaks the additive structure of (1), a natural strategy is therefore to linearize E via Fenchel conjugacy, similar to [START_REF] Agrawal | Fast algorithms for online stochastic convex programming[END_REF]. We define the Fenchel conjugate of E as E (λ) sup z∈dom(E) {λ z -E(z)} [START_REF] Boyd | Convex Optimization[END_REF], and from now on, we assume that (A1) The function E is convex.1 (A2) The Fenchel conjugate E of E has its domain

Λ dom(E ) = {λ ∈ R m : E (λ) < +∞}
that is compact so that there exists R λ max λ∈Λ λ 2 < +∞.

(A3) For any t ∈ {1, • • • , T }, the function f t is concave, and the set X t is compact and convex with X t ⊆ dom(f t ), so there exists

R x max t∈{1,••• ,T },xt∈Xt x t 2 < +∞.
Equipped with those assumptions, we introduce for any x t ∈ X t , λ ∈ Λ the Lagrangian function for round t:

L t (x t , λ) f t (x t ) -λ (A t x t -b t ) + E (λ). (5)
It can be observed that min λ∈Λ 1 T T t=1 L t (x t , λ) = P(x 1 , . . . , x T ). In the sequel, we shall refer to λ ∈ R m as dual variables. Following previous work on online saddle point optimization [START_REF] Mahdavi | Trading regret for efficiency: online convex optimization with long term constraints[END_REF]; [START_REF] Koppel | A saddle point algorithm for networked online convex optimization[END_REF], our procedure will alternate between maximizing the Lagrangian with respect to the primal variable x t and minimizing with respect to the dual variable λ. On the primal side, we note that we are interested in the 1-lookahead setting, where we have access to L t before computing our estimate xt , whereas on the dual side, our sequence of λt 's will be generated according to the 0-lookahead setting. This follows the sequencing of events in practice for display advertising, where the dual variables can only be updated after we observe A t , which is only observed after we play xt . Before summarizing the procedure in Algorithm 1, we specify a key computational assumption with respect to the primal variables:

(A4) For any λ ∈ Λ, we can efficiently compute argmax xt∈Xt L t (x t , λ).

In this paper, we will be primarily interested in problems where the above assumption (A4) holds (e.g., online ad allocation problems consisting of maximizing linear functions over simplices), so that the complexity and tractability of our proposed algorithm will mostly depend on structural properties of f t and X t . Having laid out assumptions, we present our online algorithm in Algorithm 1, which makes use of an online projected subgradient method (OPSM) [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF]; [START_REF] Hazan | Logarithmic regret algorithms for online convex optimization[END_REF] with respect to the dual variables. The notation Π Q refers to the Euclidean projection onto the set Q. Interestingly, while the analysis of the primal objective (2) requires a comparison to T different optimal primal variables (i.e., measuring performance via dynamic regret), the online optimization with respect to the dual variables reduces to a static regret analysis that compares with an unique optimal dual variable λ . A related observation was exploited by [START_REF] Shalev-Shwartz | Convex repeated games and Fenchel duality[END_REF]; Shalev-Shwartz and Kakade ( 2009), but following an opposite route, i.e., analyzing the primal objective using static regret, while leveraging the dynamic regret structure of the dual optimization problem. Given Algorithm 1, we now study the impact of applying OPSM for our problem in Lemma 1, subject to an additional assumption about the boundedness of the subgradients of L t (x t , •):

(A5) There exists G > 0, max t∈{1,••• ,T }, xt∈Xt,λ∈Λ ∇ λ L t (x t , λ) 2 ≤ G.
Lemma 1. Let assumptions (A1)-( A3) and (A5) hold. Let κ ≥ 0 be the strong convexity parameter of L t (x t , •) with respect to the 2 norm. It holds for any λ ∈ Λ, and sequences {x t , λt } T t=1 generated by Algorithm 1 that

1 T T t=1 L t (x t , λt ) -L t (x t , λ) ≤ R T ,
Algorithm 1 Online saddle point optimization

Input: Initial dual variable λ1 , step sizes η t for t = 1 to T do Receive L t //1-lookahead setting Compute xt ∈ argmax xt∈Xt L t (x t , λt ) Update λt+1 = Π Λ [ λt -η t ∇ λ L t (x t , λt )] end for
where if L t (x t , •) is strongly convex (κ > 0), we take the step size η t = 1 κt and we have

R T = G 2 2κT log(eT ), while if L t (x t , •) is convex (κ = 0), we can choose the step size η t = 2R λ G √ T leading to R T = 2R λ G √ T .
The two claims can be found respectively in [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF][START_REF] Hazan | Logarithmic regret algorithms for online convex optimization[END_REF], a short proof is given in the Appendix for self-containedness. For the convex case, the step size η t depends on the time horizon T , and we can use the doubling trick (e.g., see Section 2.3.1 in [START_REF] Shalev-Shwartz | Online learning and online convex optimization[END_REF] to remove this dependency.

Having presented the component of the regret bound that comes from the use of OPSM, we next examine why the contribution of R T is not sufficient by itself to control the regret bound (3).

Controlling the worst sequence of dual variables

Given Assumptions (A1)-( A3), there do exist some optimal primal variables {x t } T t=1 as defined in (4). Let us denote an optimal dual variable given the primal variables {x t } T t=1 as

λ ∈ argmin λ∈Λ 1 T T t=1 L t (x t , λ),
and define λ equivalently for the sequence {x t } T t=1 where the existence of both λ and λ are guaranteed by Assumptions (A1)-(A2).

From the perspective of deriving an upper bound on P -P(x 1 , . . . , xT ), we observe that L t (x t , λ ) -L t (x t , λ) can be rewritten as

(L t (x t , λ ) -L t (x t , λt )) + (L t (x t , λt ) -L t (x t , λ)),
which in turn, using (A4), is upper-bounded by

(L t (x t , λ )-L t (x t , λt )) Gt( λt) +(L t (x t , λt )-L t (x t , λ)) C , (6)
where the second term above is readily obtained by the definition of xt ∈ argmax xt∈Xt L t (x, λt ). This observation is useful in several respects: while the term C can be directly bounded as a result of Lemma 1, the term G t ( λt ), which refers to the gap we incur as a result of the sequentially-generated ( λ1 , . . . , λT ), is not accounted for by Lemma 1, and so it needs to be controlled differently. Notably, λt → G t ( λt ) is concave over Λ, and so is ( λ1 , . . . , λT ) → 1 T T t=1 G t ( λt ) over Λ T : this implies that we can cast the problem of controlling the worst sequence of dual variables ( λ1 , . . . , λT ) as a concave maximization problem.

In order to obtain a meaningful upper bound on G t ( λt ) and 1 T T t=1 G t ( λt ), we leverage the fact that we are not dealing with any general sequence ( λ1 , . . . , λT ) ∈ Λ T , but sequences possibly output by Algorithm 1. More precisely, we make this characterization by noting that OPSM generates successive λt and λt+1 estimates whose differences λt -λt+1 2 are controlled as a function of the step size η t [START_REF] Andrew | A tale of two metrics: simultaneous bounds on competitiveness and regret[END_REF]. Following [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF], we therefore introduce the convex set

Λ T,ε (λ 1 , . . . , λ T ) ∈ Λ T : T -1 t=1 λ t -λ t+1 2 ≤ ε , (7)
which will be useful shortly in deriving an upper bound for S e . In particular, assuming (A5) holds and setting ε = G T t=1 η t , for any sequence ( λ1 , . . . , λT ) ∈ Λ T generated by Algorithm 1, we have ( λ1 , . . . , λT ) ∈ Λ T,ε and

1 T T t=1 G t ( λt ) ≤ S e max (λ1,...,λ T )∈Λ T ,ε 1 T T t=1 G t (λ t ). (8)
We next turn to some lemmas that make the expression of S e more explicit by leveraging duality arguments. To this end, we introduce some additional notation: first, we rewrite the total-variation constraint from (7) as

T -1 t=1 λ t -λ t+1 2 = Ω 1/2 ((∆ ⊗ I)λ)
where λ ∈ R T •m stands for the vector formed by stacking the m-dimensional λ t 's, ⊗ is the Kronecker product,

∆ ∈ R (T -1)×T is the discrete 1-dimensional gradient matrix, 2 while Ω 1/2 is the 1 / 2 norm with dual norm Ω ∞/2 , such that Ω 1/2 (u) T -1 t=1 u t 2 and Ω ∞/2 (v) max t∈{1,••• ,T -1} v t 2 , for any vectors u = [u 1 , . . . , u T -1 ] and v = [v 1 , . . . , v T -1 ] in R (T -1)•m
. We now show that S e can conveniently be expressed as a minimization problem.

Lemma 2. Let e t = A t x tb t . It holds that

S e = min α∈R (T -1)•m S e (α)
where we have introduced the convex function α → S e (α)

1 T T t=1 E e t -(∆ ⊗ I) α t -E 1 T T t=1 e t + ε T Ω ∞/2 (α).
The proof is given in the Appendix and relies on duality arguments. We note at this juncture that our analysis can similarly handle other penalties related to the total-variation chosen in (7), for instance (∆ ⊗ I)λ2 F (we omit the details owing to space limitations).

The expression for S e provided by Lemma 2 (the result of a minimization problem), makes it possible to obtain an upper bound for any candidate vector α ∈ R (T -1)•m . We propose below one such instantiation that highlights how S e depends on the smoothness with which the sequence {e t } T t=1 varies over time: Lemma 3. Let e t = A t x tb t be stacked in the vector e ∈ R T m . The term S e is upper-bounded by

S e ≤ T max t∈{1,...,T -1} Ψ t (e ) where Ψ t (e ) t j=1 T -t T e j - T j=t+1 t T e j 2 .
The proof of this lemma can be found in the Appendix. We can see from Lemma 3 that S e captures both (a) the cumulative variations of {e t } T t=1 through Ψ t (e ), modulated by the average of the step sizes ε/T = G T t=1 η t /T , and (b) the worst of the constraint violations that surfaces via G (see the definition of Assumption (A5)).

To intuitively understand the effect of Ψ t , we can first observe that if the residual vectors {e t } T t=1 are constant over time, that is e = c • 1 for some scalar c, then the terms Ψ t (e ) vanish for all t ∈ {1, . . . , T }, so that S e does not contribute to the regret guarantees in this case. To get a better sense of the impact of Ψ t beyond the case where the residual vectors are perfectly constant, we now assume that {e t } T t=1 are independent random (sub-Gaussian) vectors. It is important to stress the fact that our analysis and our main theorem (see Theorem 1) hold in absence of any stochastic assumptions, but we only momentarily consider random residuals in order to gain insight into how the term S e can scale in more realistic scenarios beyond the case of constant residual vectors: Lemma 4. Let e ∈ R T m be a random vector, such that for some µ ∈ R mT and σ ≥ 0

E[exp(α (x -µ))] ≤ exp( α 2 σ 2 /2) (9)
holds for every α ∈ R mT . Let Z t ∈ R + and ω µ ∈ R + be defined by

Z t Ψ t (e ) and ω µ max t∈{1,...,T -1} [Ψ t (µ)] 2 . ( 10
)
Then we have

E max t∈{1,...,T -1} Z t ≤ σ 2 mT + ω µ + σ 2 mT log(T ).
The proof of this lemma is given in the Appendix. Some comments are in order. Lemma 4 shows that when e ∈ R T m is a random (sub-Gaussian) vector-which notably covers the cases where the entries of e are independent bounded, or Gaussian, random variables-then, the term S e scales (in expectation) as

O ε T σ 2 mT + ω µ + σ 2 mT log(T ) .
Assuming that ω µ = O(σ 2 mT log(T )), we can see that it is sufficient to have ε = o( T /(σ 2 m log(T ))) in order to guarantee that S e vanishes. This condition is for instance satisfied in the setting where E has Lipschitz continuous gradients with parameter L > 0, for which we can take η t = L/t and ε = O(G log(T )). Interestingly, the constant term µ appears only through ω µ , and hence Ψ t (µ), so that, as discussed previously, a constant (even non-zero) vector µ leads to ω µ = 0. We therefore see that ω µ penalizes by how much we deviate from a constant mean vector.

Having derived the terms R T and S e in the cumulative regret bound, we now turn to the description of the last term S A of our regret guarantee (3).

Estimating the matrices {A t } T t=1

So far we have assumed that we have access to A t at each round, when in practice we only get to observe A t after we take an action xt . For example, in the online ad allocation problem, the amount of money to be charged to the advertiser for a single ad impression is revealed only after we have made a decision on which ad to allocate for the impression, i.e., once we have computed xt . We now address the cost incurred in having to estimate the constraint matrices A t before playing xt . To this end, we assume that in addition to the dual variables being bounded with radius R λ , we have (A6) The set A is convex and bounded; in particular, there exists R A max A∈A A F < +∞.

Since the choice of xt now depends on the estimate Ât , we introduce the estimated Lagrangian

Lt (x t , λ) f t (x t ) -λ ( Ât x t -b t ) + E (λ), (11)
with the direct relationship (for any x t ∈ X t , λ ∈ Λ):

Lt (x t , λ) = L t (x t , λ) -λ ( Ât -A t )x t . ( 12
)
Combining the above expression (12) along with the decomposition detailed in (6), it can be shown (see Lemma A in the Appendix) that the residual term due to the estimation of A t is given by

S A 1 T T t=1 λ t ( Ât -A t )(x t -xt ). ( 13
)
As a result of having to estimate A t , we present Algorithm 2, where we also apply an OPSM to estimate the constraint matrices A t 's. Lemma 5 then quantifies the cumulative regret incurred as a result of having to use estimated constraint matrices Ât to produce xt , where the result stems from the analysis of [START_REF] Hall | Dynamical models and tracking regret in online convex programming[END_REF]

Receive A t Update λt+1 = Π Λ [ λt -η t ∇ λ L t (x t , λt )] //G t a subgradient of A t -A F Update Ât+1 = Π A [ Ât -ν t G t ] end for
Lemma 5. Consider Algorithm 2. For OPSM on matrices Ât with step size

ν t = R A / √ t, we have 1 T T t=1 Ât -A t F ≤ 3 √ T R A + T t=1 A t -A t+1 F ,
and the term S A from the general regret bound (3), as defined in ( 13), is upper-bounded by:

S A ≤ 6R λ R x √ T R A + T t=1 A t -A t+1 F .
Having introduced the main three components of our regret bound, we next formally present our main results.

Main results and discussion

We start by stating our core theorem:

Theorem 1. Assume (A1)-(A6) hold, and let x 1 , . . . , x T be a sequence of optimal primal solutions for (2). Define the optimal error vectors e t A t x tb t stacked into the vector e ∈ R T m . Let M e max t∈{1,••• ,T -1} Ψ t (e ). Algorithm 2 generates a sequence (x 1 , . . . , xT ) ∈ T t=1 X t satisfying

P -P(x 1 , . . . , xT ) ≤ R T + S e + S A where R T = 2R λ G √ T , S e = M e 2R λ √ T and S A = 6R λ Rx √ T [R A + T t=1 A t -A t+1 F ].
If we additionally assume that (A7) E has Lipschitz continuous gradients over its domain with parameter L > 0, we can take instead R T = LG 2 2T log(eT ) and the term M e

2R λ √

T can be replaced by M e

LG log(eT ) T

.

The proof of the results can be found in the Appendix. The statement of Theorem 1 calls for some comments. Omitting the contribution of terms depending on {R λ , R x , R A , G}, Theorem 1 guarantees that when A t is assumed to be known, the regret P -P(x 1 , . . . , xT ) is upper-bounded by terms scaling with O(M e / √ T ), or O(M e log(T )/T ) respectively for the convex or strongly convex cases. As a result, the upper bound is mainly driven by how smooth the sequence {e t } T t=1 varies over time, as measured by M e . In the more challenging setting where the constraint matrices A t 's are also estimated, we pay an additional cost O

1 √ T [1 + T t=1 A t -A t+1 F ]
, so that the potential strong convexity of E * (or equivalently, the gradient Lipschitz continuity of E) plays only a secondary role compared to the leading term depending on the smoothness of the A t 's.

Possible instantiations of E: We now present possible valid instantiations of the error function E. A first example for which the set of assumptions (A1)-(A2) can be satisfied considers E(z) = r • Ω(z) where Ω refers to any norm on R m . Indeed, it is well known (e.g., see Example 3.26 in [START_REF] Boyd | Convex Optimization[END_REF]) that E corresponds in this

E(z) E (λ) Λ H R λ ,L ( z 2 ) I B R λ (λ) + 1 2L λ 2 2 B R λ H R λ ,L ( [z] + 2 ) I B R λ ∩R m + (λ) + 1 2L λ 2 2 B R λ ∩ R m + R λ • z 2 I B R λ (λ) B R λ R λ • [z] + 2 I B R λ ∩R m + (λ) B R λ ∩ R m +
Table 1: Examples of non-additive penalty functions. For simplicity, we display only expressions based on the 2 norm, but similar formula can be obtained for q , q ∈ {1, ∞}. We recall that Λ dom(E ) = {λ : E (λ) < +∞}.

We define B r as the Euclidean ball in R m with radius r, while I C is the indicator function of the convex set C. The function H r,s is defined as H r,s (t)

1 2 min{st 2 , r 2 s } + r[|t| -r s ] + .
The expressions are derived in Lemma B in the Appendix.

case to the indicator function of the ball for the dual norm Ω * with radius r. In particular, if Ω is the 2 norm, we can take R λ = r.

Tighter regret bounds can be obtained if the function E is taken to be gradient Lipschitz continuous (see Theorem 1 and Assumption (A7)), while preserving the boundedness of the domain Λ of its Fenchel conjugate. A possible choice in this case is an instance of a Huber function (e.g., see Section 10.6 in [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]), as detailed in Table 1.

In the setting of the online ad allocation problem, it is sometimes required to impose an asymmetric error function (e.g., to penalize under-and over-delivery differently). We show in Table 1 corresponding valid instantiations of E that depend only on their arguments via their positive parts.

Complexity of computing xt : Assumption (A4), about the exact computational oracle for xt , hides, and concentrates, the difficulties related to the optimization with respect to the primal variables x t . We can observe, in the light of the decomposition (6), that an approximate maximization-for instance leveraging the concavity of f t -would lead to residual terms that would be hard to control in our online setting. Although it may appear at first sight that an exact maximization is an overly strong requirement, wide classes of problems fall within the scope of this assumption; we can for instance cite the exactly-solved subproblems that are commonly encountered in the context of proximal methods [START_REF] Parikh | Proximal algorithms[END_REF] or conditional-gradient algorithms [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF], see Table 1). We conjecture that, thanks to Assumption (A4), we may weaken Assumption (A3) by, for instance, trying to drop the concavity assumption of the f t 's and the convexity of the X t 's.

Lower bounds: Based on Assumptions (A1)-( A2)-(A3), we know that the minimax equality for our saddle point problem holds [START_REF] Sion | General Minimax Theorems[END_REF]. In particular, defining the (convex) dual function

D(λ) 1 T T t=1 D t (λ) 1 T T t=1 max xt∈Xt L t (x t , λ),
we have the equality

P -P(x 1 , . . . , xT ) = min λ∈Λ D(λ) -P(x 1 , . . . , xT ).
This new relationship makes it possible to derive complementary guarantees for our problem. For instance, when we need not estimate the constraint matrices A t 's,3 it can be shown that we can lower bound P -P(x 1 , . . . , xT ) by

1 T T t=1 D t ( λt ) -P(x 1 , . . . , xT ) -R T + ,
which sets a lower limit on the best performance we could get with the sequences {x t , λt } T t=1 . Interestingly, this lower bound can not only be practically computed, but it also gives an indication about the inherent difficulty of the problem at hand (e.g., if it is large, the online strategy cannot compete efficiently with its offline counterpart).

Experiments

In this section, we conduct two sets of experiments. First, we consider a synthetic data and demonstrate the benefit of handling truly non-additive penalties over a heuristic based on an additive relaxation. Second, we report results on an online advertising data set consisting of a sample of 3.3 million bid requests gathered from a large ad serving system used in production at Amazon. Hence, these data account for external advertiser constraints encountered in practice, such as user behavioral targeting. We show that the rate at which regret vanishes matches our theory.

Additive versus non-additive modeling

We first investigate if the use of a non-additive penalty on average constraint violations yields improvements in reward and/or constraint violations compared to the additive error formulation adopted in standard online learning. The additive long-term penalty is defined as follows:

1 T T t=1 f t (x t ) -E(A t x t -b t ) . ( 14
)
In the experiments, we consider a linear reward function f t (x) = u t x t , with the simplicial constraint x t ∈ X t {x t ∈ R d + : 1 x t ≤ 1}. These choices of f t and X t are reminiscent of the online ad allocation setting we explore at greater length in Section 4.2. We simulated data with m = 25, d = 10, T = 200 and where {A t , b t , u t } were generated to be standard random Gaussian matrices and vectors normalized to unit norms. We ran our online algorithm with primal-dual updates dictated by the non-additive problem formulation in (2) and compared it to a baseline algorithm with primal-only updates following a standard online learning formulation with additive penalties on constraint violations as in ( 14). Since the exact maximization max xt∈Xt [f t (x t ) -E(A t x tb t )] may not be obtained in closed-form in this case, we used the solver CVXPY Diamond et al. (2014). We shall refer to the two approaches as respectively Non-additive and Additive.

So as to cover different types of geometries and convexity assumptions, we consider several choices of penalty functions, namely E(z) = R λ • z q for q ∈ {1, 2, ∞}, along with E(z) = R λ • H 1,1 ( z 2 ) where H 1,1 is the Huber function defined in Table 1. For each E, and for, both, Non-additive and Additive, we compute the reward over T rounds given by 1 T t u t xt and the (normalized) non-additive penalty E( 1 T t A t xtb t )/R λ for the primal variables xt generated online by the two algorithms.

Figure 1 shows the resulting reward versus constraint violation for varying values of R λ , averaged over 10 generations of {A t , b t , u t }. As can be seen from the curves, in applications where constraint violation are to be measured as a non-additive penalty, the additive relaxation ( 14) leads to significant deterioration of the performance compared to the non-additive penalty. In particular, the domain of the achievable constraint violation is very narrow when applying the additive heuristic, meaning that there is little room for making trading-offs. As expected, we also observe that points of Non-additive and Additive superimpose in the regime where R λ 1, i.e., when the two formulations focus on the optimization of the reward. Finally, we remark that since the quality of the additive relaxation ( 14) essentially hinges on the gap in the Jensen's inequality E( 1

T T t=1 A t x t -b t ) ≤ 1 T T t=1 E(A t x t -b t )
, and therefore on the distribution of the residuals {A t x tb t } T t=1 , we provide in the Appendix additional simulations where {A t , b t , u t } are generated according to different distributions, viz, Cauchy, uniform and gamma. In a nutshell, the same conclusions hold, even for the distributions that appear to make the relaxation (14) tighter.

Regret convergence on real online advertising data

In the second experiment, we focus on the online ad allocation problem in display advertising subject to long-term constraints on budget consumed per ad. More formally, in the display advertising setting, each user visit to a website or app triggers a bid request i, where each bid request has some subset of the m possible ads that can be served. Indexing ads by j ∈ {1, • • • , m}, the welfare of serving an impression for ad j (i.e., the value of showing the ad once to a user) for bid request i is given by u ij ≥ 0. For bid request i, let primal variable x ij = 1 correspond to our decision to allocate ad j an impression, with x ij = 0 otherwise. Non-additive Additive 0.00 0.02 0.04 0.06 0.08 0.10

1 T T t=1 A t xt -b t 2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 1 T T t=1 u t xt
Non-additive Additive 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10

H 1,1 ( 1 T T t=1 A t xt -b t 2 ) Non-additive Additive Figure 1: Reward 1 T t u t as a function of non-additive penalty E( 1 T t A t xt -b t )/R λ for R λ = 2 γ
with γ ∈ {-8, -7.5, . . . , 10}. Red crosses correspond to our proposed online algorithm, while blue circles stand for a baseline algorithm with additive penalties. Each subplot displays a different instantiation of E, namely

E(z) = R λ • z q for q ∈ {1, 2, ∞} and E(z) = R λ • H 1,1 ( z 2 ) where H 1,1 is defined in Table 1 (best seen in color).
To map this problem to our formulation, we will further partition the time dimension into disjoint intervals t ∈ {1, • • • , T } where the number of bid requests in each interval is equal to N . We let I t be the set of bid requests in round t with |I t | = N . Thus, for a given t, the primal variable x t will denote the flattened matrix x ij for i ∈ I t and for all ads j ∈ {1, • • • , m}, so that x t ∈ [0, 1] d with d = N • m. Since we can show at most one ad per bid request, we additionally have the simplicial constraints m j=1 x ij ≤ 1 for each bid request i, and we therefore define the set X t ⊆ [0, 1] d as the Cartesian product of those simplices for each i ∈ I t . We define the welfare vector u t ∈ R d following the same flattening operation, so that f t (x t ) = u t x t . Matrix A t ∈ R m×d is such that the j-th entry [A t x t ] j represents the amount of budget consumed for ad j as a result of the allocation vector x t under a cost per impression model. Finally, the non-additive long-term constraints are modeled as 1 T t A t x tb t . Vector b t is equal to the constant vector b of ad budgets. Hence, the penalty on 1 T t A t x tb t penalizes deviations from spending 100% of each ad's budget. We stress the fact that, although b is a constant vector, we are in a non-stationary regime since A t is not equal to some constant matrix A independent of t.

The goal of the experiment is to compare the regret behavior of our proposed online algorithm for different choices of non-additive E on the online advertising data set from live traffic. For each bid request in our data set, each eligible ad candidate comes with a pre-defined welfare and cost to the advertiser, which determine u t and A t respectively. We partitioned the 3.3 million impressions into batches of equal size, for a total of T = 10000 rounds. We computed the average cumulative regret

P -P(x 1 , • • • , xT ) as a function of T for (symmetric) E chosen to be i) E(z) = R λ • z 1 referred to as Convex and ii) E(z) = H R λ ,1 ( z 2
) referred to as Strongly convex. For all problem instantiations, we set R λ = 50000. We computed the term P(x 1 , • • • , xT ) by running our online primal-dual algorithm over 20 permutations of bid requests. The term P was obtained using an offline primal-dual method where the dual variable updates correspond to an offline subgradient (resp. gradient) method for the convex (resp. strongly convex) E * . Figure 2 shows the resulting cumulative regret as a function of T , averaged over the 20 runs (with error bars that are negligible). We observe that the cumulative regrets incurred for, both, Convex and Strongly convex decrease as T increases, with the latter decreasing at faster rate than the former as predicted by the theory. 5 Appendix: Proofs

Proof of Lemma 1

Proof. The two claims of the lemma can be found respectively in [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF]; [START_REF] Hazan | Logarithmic regret algorithms for online convex optimization[END_REF] and follows from (possibly, strong) convexity of L t (x t , λ) with respect to λ. For completeness, we briefly repeat the core arguments from [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF]; [START_REF] Hazan | Logarithmic regret algorithms for online convex optimization[END_REF]. The conclusions come from the fact that the term

T t=1 L t (x t , λt ) -L t (x t , λ)
is shown to be upper bounded by, with a t λ -λt

2 2 , 1 2 T t=1 1 η t a t -a t+1 -κa t + η t G 2
in turn upper bounded by a 1 2

1 η 1 -κ + G 2 2 T t=1 η t + 1 2 T t=2 a t 1 η t - 1 η t-1 -κ .
The choices of the step sizes lead to the advertised instantiations of R T .

Proof of Lemma 2

Proof. The result stems from an application of strong duality for the above convex program. First, notice that

1 T T t=1 G t (λ t ) = 1 T T t=1 λ t e t -E * (λ t ) -E 1 T T t=1 e t .
Introduce the equality constraint a = (∆ ⊗ I)λ-we remind that λ ∈ R T •m stands for the vector formed by stacking the m-dimensional λ t 's-and write the corresponding Lagrangian (we momentarily omit the constant terms independent of λ t that have no effects on the maximization), that is, for any

λ t ∈ Λ, a ∈ R (T -1)•m , α ∈ R (T -1)•m and β ≥ 0: L(λ, a, α, β) = 1 T T t=1 λ t e t -E * (λ t ) + α (a -(∆ ⊗ I)λ) + β(ε -Ω 1/2 (a)).
Using the fact that E * is the Fenchel conjugate of E and that the conjugate of the 1 / 2 norm Ω 1/2 is the indicator function of the unit ball for its dual ∞ / 2 norm (Example 3.26 in [START_REF] Boyd | Convex Optimization[END_REF]), we can maximize out the primal variables:

max λ,a L(λ, a, α, β) = 1 T T t=1 E e t -T (∆ ⊗ I) α t +βε
under the constraint max t=1...T -1 α t 2 ≤ β. Minimizing out β, and making the change of variable α ← T • α leads to the claimed dual minimization problem. Finally, the initial claim holds with equality, since Slater's constraint qualification applies here to guarantee strong duality (Section 5.2.3 in [START_REF] Boyd | Convex Optimization[END_REF]).

Proof of Lemma 3

Proof. The equality in Lemma 2 has a minimization over α on its right-hand-side. Therefore, we have that for every α ∈ R (T -1)•m we have S e upper bounded by S e (α). Specifically, let α [((∆∆ ) -1 ∆) ⊗ I]e * and we have that S e ≤ S e (α ). First, we compute the term inside the first sum of S e (α ) for a fixed t ∈ {1, . . . , T }

e * t -(∆ ⊗ I) [((∆∆ ) -1 ∆) ⊗ I]e * t = e * -(∆ ⊗ I) [((∆∆ ) -1 ∆) ⊗ I]e * t = e * -([(∆ (∆∆ ) -1 ∆) ⊗ I]e * t = [(I -∆ (∆∆ ) -1 ∆) ⊗ I]e * t
The matrix ∆ ∈ R (T -1)×T is defined by

∆ ij =      1 i = j, -1 i = j -1, 0 otherwise.
By direct computation one can show, that ∆∆ ∈ R (T -1)×(T -1) is the tridiagonal matrix

(∆∆ ) ij =      2 i = j, -1 |i -j| = 1, 0 otherwise.
There exists a close form solution for the inversion of an arbitrary tridiagonal matrix [START_REF] Usmani | Inversion of a tridiagonal jacobi matrix[END_REF]. Using this we get that the matrix (∆∆ ) -1 ∈ R (T -1)×(T -1) is defined by

[(∆∆ ) -1 ] ij = i(T -j) T i < j, (T -i)j T i ≥ j.
No,w we can again proceed with a direct calculation to compute (∆∆ ) -1 ∆ ∈ R (T -1)×T and get

[(∆∆ ) -1 ∆] ij = -i T i < j, T -i T i ≥ j.
Finally, computing the full product ∆ (∆∆ ) -1 ∆ ∈ R T ×T we get

[∆ (∆∆) t-1 ∆] ij = T -1 T i = j, -1 T i = j.
Observe, that the matrix I -∆(∆∆ ) -1 ∆ ∈ R T ×T is a matrix full of values 1/T . It follows that

e * t -(∆ ⊗ I) [((∆∆ ) -1 ∆) ⊗ I]e * t = 1 T T i=1
e * i , ∀t ∈ {1, . . . , T } Also, observe that the first two terms in S e (α ) cancel out and we have

S e (α ) = T Ω ∞/2 (α ) = T Ω ∞/2 [((∆∆ ) -1 ∆) ⊗ I]e *
Plugging in the matrix (∆∆) -1 ∆ computed earlier in the above expression we get the final result.

Proof of Lemma 4

Let W t ∈ R m×mT be a block matrix defined by

W t      T -t T I . . . T -t T I t - t T I . . . - t T I (T -t)      (15)
Observe that Z t = W t e * 2 . From now on t will be assumed to belong to {1, . . . , T -1} in all derivations. Let s > 0. Then using the convexity and monotonicity of the exponential together with Jensen's inequality, the fact that sum of elements is more than the maximal element, linearity of expectation, and the fact that the sum can be bounded by the number of terms times the maximal element, we get the following sequence of inequalitites Observe that the right-hand-side of ( 17) is increasing in the arguments Tr(W t t), Tr((W t W t ) 2 ), W t µ 2 , and W t W t 2 . From the definition of W t in (15) and a maximization in t we have

e sE[maxt Z 2 t ] ≤ E e s maxt Z 2 t = E max t e sZ 2 t ≤ E t e sZ 2 t = t E e sZ 2 t ≤ T max t E e sZ 2 t (16) Assume 0 ≤ s < 1/(2σ 2 W t W t 2 ).
W t 2 F (15) = mt(T -t) T ≤ mT 4
which can be further used to bound

Tr(W t W t ) = W t 2 F ≤ mT 4 Tr((W t W t ) 2 ) = W t W t 2 F ≤ W t 4 F ≤ m 2 T 2 16 W t W t 2 ≤ W t W t F ≤ W t 2 F ≤ mT 4
using standard arguments. Also, from definition (10) we have W t µ 2 ≤ ω µ . Combining these bounds with (17) and setting s s * defined by

0 <   s * 4 log(T ) √ σ 2 mT σ 2 mT + 8ω µ + 2 σ 2 mT log(T )   < 2 σ 2 mT ≤ 1 2σ 2 W t W t 2 (18)
we get

E[exp(s * W t e * 2 )] ≤ exp s * 4 σ 2 mT + σ 4 m 2 T 2 s * + 16ω µ 4 -2σ 2 mT s * (18) = exp s * 4 σ 2 mT + 4ω µ + σ 2 mT + 8ω µ σ 2 mT log(T ) . (19)
We finally have all the ingredients. Using the definition of Z t in (10), using the non-negativity of Z t , the concavity of the square root together with Jensen's inequality, the definition of the exponential, the inequality ( 16), the inequality (19) with the observation that the right-hand-side does not depend on t, the definition of s * in ( 18), the bound √ s ≤ 1 2 (s + 1), the fact that all the entries are positive, and

a + b = √ a 2 + 2ab + b 2 we get E max t Z t = E max t Z 2 t ≤ E[max t Z 2 t ] = log e s * E[maxt Z 2 t ] /s * (16) ≤ log T max t E e s * Z 2 t /s * (19) ≤ log(T ) s * + 1 4 σ 2 mT + 4ω µ + σ 2 mT + 8ω µ σ 2 mT log(T ) (18) = 1 2 σ 2 mT + 4ω µ + 2 σ 2 mT + 8ω µ σ 2 mT log(T ) + 2σ 2 mT log(T ) ≤ 1 2 2σ 2 mT + 4ω µ + 2 σ 2 mT + 8ω µ σ 2 mT log(T ) + σ 2 mT log(T ) ≤ σ 2 mT + ω µ + 2 σ 2 mT + ω µ σ 2 mT log(T ) + σ 2 mT log(T ) ≤ σ 2 mT + ω µ + σ 2 mT log(T )
which concludes the proof.

Proof of Lemma 5

Proof. We first apply the Cauchy-Schwartz inequality to get

| λ t ( Ât -A t )(x t -xt )| ≤ λt 2 x t -xt 2 ||| Ât -A t ||| 2 ≤ 2R λ R x Ât -A t F
where we have used that the operator norm is smaller than the Frobenius norm. We then leverage Theorem 4 from [START_REF] Hall | Dynamical models and tracking regret in online convex programming[END_REF]: in their notation, we have G ≤ 1 (since the sub-gradients of A → A t -A F are bounded by 1), σ = 1, M ≤ R A /2 and D max ≤ 2R 2 A along with Φ = I and the learning rate

η t = R A / √ t. This leads to T t=1 Ât -A t F - min Ã1,..., ÃT ∈A T t=1 Ãt -A t F ≤ 2R 2 A η T +1 + 2R A η T V + 1 2 T t=1 η t ≤ 2R A η T (R A + V ) + R A √ T ≤ 3 √ T (V + R A ),
where V stands for

T t=1 A t -A t+1 F .

Proof of Theorem 1

Proof. The proof consists in putting together the components introduced in the Section 3 of the core paper. We start from

P -P(x 1 , . . . , xT ) = 1 T T t=1 L t (x t , λ ) -L t (x t , λ)
which, using the decomposition (6), is equal to Recalling the definitions of R T , G t and S A respectively in Lemma 1, ( 6) and ( 13), we obtain

1 T T t=1 (L t (x t , λ ) -L t (
P -P(x 1 , . . . , xT ) ≤ 1 T T t=1 G t ( λt ) + R T + S A .
Noticing that with the choice ε/T G T t=1 η t /T , the relationship (8) holds, we finally obtain P -P(x 1 , . . . , xT ) ≤ R T + S e + S A .

The rest of the proof follows by instantiating R T , S e and S A . The value of R T is given in Lemma 1, while the upper bound for S e is described in Lemma 3.

Moreover, ε/T is equal to 2R λ / √ T and GL log(eT )/T in the convex and strongly convex cases respectively (given the η t from Lemma 1). Finally, the upper bound on S A is given by Lemma 4.

6 Appendix: Technical lemmas Lemma A. For any λ ∈ Λ and x t ∈ X t , we have L t (x t , λt ) -L t (x t , λ) ≤ L t (x t , λt ) -L t (x t , λ) + λ t ( Ât -A t )(x txt ). The problem instances are generated with a standard Cauchy distribution. Red crosses correspond to our proposed online algorithm, while blue circles stand for a baseline algorithm with additive penalties. Each subplot displays a different instantiation of E, namely E(z) = R λ • z q for q ∈ {1, 2, ∞} and E(z) = R λ • H 1,1 ( z 2 ) where H 1,1 is defined in Table 1. The problem instances are generated with an uniform(-1,1) distribution. Red crosses correspond to our proposed online algorithm, while blue circles stand for a baseline algorithm with additive penalties. Each subplot displays a different instantiation of E, namely E(z) = R λ • z q for q ∈ {1, 2, ∞} and E(z) = R λ • H 1,1 ( z 2 ) where H 1,1 is defined in Table 1. 0.09 0.10 0.11 0.12 0.13 0.14 0.15 0.16 1 T T t=1 The problem instances are generated with a gamma(2,2) distribution. Red crosses correspond to our proposed online algorithm, while blue circles stand for a baseline algorithm with additive penalties. Each subplot displays a different instantiation of E, namely E(z) = R λ • z q for q ∈ {1, 2, ∞} and E(z) = R λ • H 1,1 ( z 2 ) where H 1,1 is defined in Table 1.

  xtb t 1

Figure 2 :

 2 Figure 2: Cumulative regret, averaged over 20 runs, measured on the data collected by a real ad serving system as a function of the number of rounds T for E(z) = R λ z 1 (red) and E(z) = H R λ ,1 ( z 2 ) (blue).Figure best seen in color.

Figure 3 :

 3 Figure 3: Reward 1 T t u t xt as function of non-additive penalty E( 1T t A t xtb t )/R λ for R λ = 2 γ with γ ∈ {-8, -7.5, . . . , 10}.The problem instances are generated with a standard Cauchy distribution. Red crosses correspond to our proposed online algorithm, while blue circles stand for a baseline algorithm with additive penalties. Each subplot displays a different instantiation of E, namelyE(z) = R λ • z q for q ∈ {1, 2, ∞} and E(z) = R λ • H 1,1 ( z 2 )where H 1,1 is defined in Table1.

Figure 4 :

 4 Figure 4: Reward 1 T t t xt as a function of non-additive penalty E( 1T t A t xtb t )/R λ for R λ = 2 γ with γ ∈ {-8, -7.5, . . . , 10}.The problem instances are generated with an uniform(-1,1) distribution. Red crosses correspond to our proposed online algorithm, while blue circles stand for a baseline algorithm with additive penalties. Each subplot displays a different instantiation of E, namely E(z) = R λ • z q for q ∈ {1, 2, ∞} and E(z) = R λ • H 1,1 ( z 2 ) where H 1,1 is defined in Table1.

Figure 5 :

 5 Figure 5: Reward 1 T t t xt as a function of non-additive penalty E( 1T t A t xtb t )/R λ for R λ = 2 γ with γ ∈ {-8, -7.5, . . . , 10}.The problem instances are generated with a gamma(2,2) distribution. Red crosses correspond to our proposed online algorithm, while blue circles stand for a baseline algorithm with additive penalties. Each subplot displays a different instantiation of E, namely E(z) = R λ • z q for q ∈ {1, 2, ∞} and E(z) = R λ • H 1,1 ( z 2 ) where H 1,1 is defined in Table1.

.

  Algorithm 2 Online saddle point optimization with estimated constraint matrices Ât Input: Initial Â1 and dual variable λ1 , step sizes η t , ν t for t = 1 to T do Receive Lt //1-lookahead and Ât Compute xt ∈ argmax xt∈Xt Lt (x t , λt )

  A t xtb t ∞

			0.65
			0.55 0.60		Non-additive Additive			Non-additive Additive
	T t=1 u t xt	0.40 0.45 0.50
	1	T	0.35
			0.30
			0.25
			0.20 0.65				2.50 2.55 2.60 2.65 2.70 2.75 2.80 2.85 2.90 1 T T t=1 A t xt -b t 1
	1 T t=1 u T t xt	0.48 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60	0.50	0.52 T t=1 A t xt -b t 2 0.54 0.56 0.58 Non-additive 1 T Additive	0.60	0.48	0.50 H 1,1 ( 1 0.52 T t=1 A Non-additive 0.54 T Additive	0.56	0.58	0.60

t xtb t 2 )

When we refer to convex/concave functions, we implicitly assume closed proper convex/concave functions[START_REF] Boyd | Convex Optimization[END_REF].

The matrix ∆ ∈ R (T -1)×T contains two non-zero entries per row, with ∆t,t = 1 and ∆ t,t+1 = -1 for t ∈ {1, • • • , T -1}.

When At's are also estimated, the lower bound is more involved, and requires the introduction of an estimated dual function Dt along the lines of (11).

Proof. For any λ ∈ Λ and x t ∈ X t , L t (x t , λt ) -L t (x t , λ) is equal to

, where we have used twice the relationship Lt (x t , λ) = L t (x t , λ)λ ( Ât -A t )x t , and the fact that xt is defined as argmax xt∈Xt L t (x t , λt ).

We next provide the details of the computations related to relationship between E and E * presented in Table 1. We focus on one of the displayed instantiations, since the arguments in those other cases follow along the same lines.

where

We start by deriving the Lagrangian associated with

Maximizing out w, and recognizing the conjugate of the squared 2 norm (Example 3.27 in [START_REF] Boyd | Convex Optimization[END_REF])), we obtain the dual function

In turn, we minimize with respect to the dual variables

The optimal α can then be easily computed and is equal to

Plugging back this value into the dual function leads to the expression H r,s ( [u] + 2 ). The equality holds by invoking strong duality, which applies based on Slater's constraint qualification (Section 5.2.3 in [START_REF] Boyd | Convex Optimization[END_REF]).

7 Appendix: Additional experiments for Section 4.1

We show below additional results when {A t , b t , u t } are generated according to different distributions, namely, Cauchy, uniform and gamma. As described in the protocol of Section 4.1, we continue to normalize {A t , b t , u t } to unit norms. The same conclusions as those explained in Section 4.1 of the paper hold. We observe that some instantiations of E and distribution, e.g., 1 with gamma in Figure 5, lead to settings where the additive relaxation (12) appears as tighter, although our non-additive approach still offers better reward versus constraint violation tradeoffs.