Towards Inference Delivery Networks:
Distributing Machine Learning
with Optimality Guarantees

Tareq Si Salem*, Gabriele Castellano*†, Giovanni Neglia*, Fabio Pianese†, Andrea Araldo‡

*Inria, Université Côte d'Azur, France, {tareq.si-salem, gabriele.castellano, giovanni.neglia}@inria.fr,

†Nokia Bell Labs, France, {fabio.pianese, gabriele.castellano.ext}@nokia.com,

‡Samovar, Télécom SudParis, Institut Polytechnique de Paris, France, andrea.araldo@telecom-sudparis.eu

**Abstract**—An increasing number of applications rely on complex inference tasks that are based on machine learning (ML). Currently, there are two options to run such tasks: either they are served directly by the end device (e.g., smartphones, IoT equipment, smart vehicles), or offloaded to a remote cloud. Both options may be unsatisfactory for many applications: local models may have inadequate accuracy, while the cloud may fail to meet delay constraints. In this paper, we present the novel idea of inference delivery networks (IDNs), networks of computing nodes that coordinate to satisfy ML inference requests achieving the best trade-off between latency and accuracy. IDNs bridge the dichotomy between device and cloud execution by integrating inference delivery at the various tiers of the infrastructure continuum (access, edge, regional data center, cloud). We propose a distributed dynamic policy for ML model allocation in an IDN by which each node dynamically updates its local set of inference models based on requests observed during the recent past plus limited information exchange with its neighboring nodes. Our policy offers strong performance guarantees in an adversarial setting and shows improvements over greedy heuristics with similar complexity in realistic scenarios.

# I. INTRODUCTION

Machine learning (ML) models are often trained to perform inference, that is to elaborate predictions based on input data. ML model training is a computationally and I/O intensive operation and its streamlining is the object of much research effort. Although inference does not involve complex iterative algorithms and is therefore generally assumed to be easy, it also presents fundamental challenges that are likely to become dominant as ML adoption increases [1]. In a future where AI systems are ubiquitously deployed and need to make timely and safe decisions in unpredictable environments, inference requests will have to be served in real-time and the aggregate rate of predictions needed to support a pervasive ecosystem of sensing devices will become overwhelming.

Today, two deployment options for ML models are common: inferences can be served by the end devices (smartphones, IoT equipment, smart vehicles, etc.), where only simple models can run, or by a remote cloud infrastructure, where powerful “machine learning as a service” (MLaaS) solutions rely on sophisticated models and provide inferences at extremely high throughput.

However, there exist applications for which both options may be unsuitable: local models may have inadequate accuracy, while the cloud may fail to meet delay constraints.

As an example, popular applications such as recommendation systems, voice assistants, and ad-targeting, need to serve predictions from ML models in less than 200 ms. Future wireless services, such as connected and autonomous cars, industrial robotics, mobile gaming, augmented/virtual reality, have even stricter latency requirements, often below 10 ms and in the order of 1 ms for what is known as the tactile Internet [2]. In enabling such strict latency requirements, the advent of Edge Computing plays a key role, as it deployes computational resources at the edge of the network (base stations, access points, ad-hoc servers). However, edge resources have limited capacity in comparison to the cloud and need to be wisely used. Therefore, integrating ML inference in the continuum between end devices and the cloud—passing through edge servers and regional micro data-centers—will require complex resource orchestration.

We believe that, to allocate resources properly, it will be crucial to study the trade-offs between accuracy, latency and resource-utilization, adapted to the requirements of the specific application. In fact, inference accuracy and, in general, resource efficiency increase toward the cloud, but so does communication latency. In this paper, we present the novel idea of inference delivery networks (IDN): networks of computing nodes that coordinate to satisfy inference requests achieving the best trade-off. An IDN may be deployed directly by the ML application provider, or by new IDN operators that offer their service to different ML applications, similarly to what happens for content delivery networks. The same inference task can be served by a set of heterogeneous models featuring diverse performance and resource requirements (e.g., different model architectures [3], multiple downsized versions of the same pre-trained model [4], different configurations and execution setups). Therefore, we study the novel problem of how to deploy the available ML models on the available IDN nodes, where a deployment strategy consists in two coupled decisions: (i) where to place models for serving a certain task and (ii) how to select their size/complexity among the available alternatives.

In this paper, we first define a specific optimization problem for ML model allocation in IDNs. We characterize the complexity of such problem and then introduce INFIDA (INFerence Intelligent Distributed Allocation), a distributed dynamic allocation policy. Following this policy, each IDN node periodically updates its local allocation of inference
models on the basis of the requests observed during the recent past and limited information exchange with its neighbors. The policy offers strong performance guarantees in an adversarial setting [5], that is a worst case scenario where the environment evolves in the most unfavorable way. Numerical experiments in realistic settings show that our policy outperforms heuristics with similar complexity. Our contributions are as follows:

(1) We present the novel idea of inference delivery networks (IDNs).

(2) We frame the allocation of ML models in IDNs as an (NP-hard) optimization problem that captures the trade-off between latency and accuracy, and study how this problem diverges from settings considered in previous works (Sec. III).

(3) We propose INFIDA, a distributed and dynamic allocation algorithm for IDNs (Sec. IV), and we show it provides strong guarantees in the adversarial setting, providing novel theoretical results in approximating budget-additive (submodular) set functions (Sec. V).

(4) We evaluate INFIDA in a realistic simulation scenario and compare its performance both with an offline greedy heuristic and with its online variant under different topologies and trade-off settings (Sec. VI).

## II. RELATED WORK

The problem of machine learning is often reduced to the training task, i.e., producing statistical models that can map input data to certain predictions. A considerable amount of existing works addresses the problem of model training: production systems such as Hadoop [6] and Spark [7] provide scalable platforms for analyzing large amount of data on centralized systems, and even the problem of distributing the training task over the Internet has been largely addressed recently by many works on federated learning [8]–[13]. However, there is surprisingly less research on how to manage the deployment of ML models once they have been trained (inference provisioning).

Most of the existing solutions on inference provisioning (e.g., Tensorflow Serving [14], Azure ML [15], and Cloud ML [16]) address the scenario where inference queries are served by a data center. Recent works [17]–[20] propose improvements on performance and usability of such cloud inference systems. Clipper [17] provides a generalization of TensorFlow Serving [14] to enable the usage of different ML frameworks, such as Apache Spark MLLib [21], Scikit-Learn [22], and Caffe [23]. The authors of [18] propose a reinforcement learning scheduler to improve the system throughput. INFaaS [19] provides a real-time scheduling of incoming queries on available model variants, and scales deployed models based on load thresholds. Last, InferLine [20] extends Clipper to minimize the end-to-end latency of a processing pipeline, both periodically adjusting the models allocation and constantly monitoring and handling unexpected query spikes; the solution can be applied to any inference serving system that features a centralized queue of queries. All these solutions address the problem of inference provisioning in the scenario where the requests are served within a data

center and are not suitable for a geographically distributed infrastructure where resources are grouped in small clusters and network latency is crucial (e.g., Edge Computing). For instance, none of the previous works consider the network delay between different compute nodes, it being negligible in a data center.

For what concerns inference provisioning in constrained environments, fewer works exist. Some solutions attempt to adapt inference to the capabilities of mobile hardware platforms through the principle of model splitting, a technique that distributes a ML model by partitioning its execution across multiple discrete computing units. Model splitting was applied to accommodate the hardware constraints in multi-processor mobile devices [24], to share a workload among mobile devices attached at the same network edge [25], and to partially offload inferences to a remote cloud infrastructure [26], possibly coupled with early exit strategies [27] and conditional hierarchical distributed deployment [28]. Model splitting is orthogonal to our concerns and could be accounted for in an enhanced IDN scheme.

There has been some work on ML model placement at the edge in the framework of what is called “AI on Edge” [29], but it considers a single intermediate tier between the edge device and the cloud, while we study general networks with nodes in the entire cloud-to-the-edge continuum. Our dynamic placement INFIDA algorithm could be applied also in this more particular setting, for example in the MODI platform [30]. The work closest to ours is [31], which proposes an online learning policy, with the premise of load balancing over a pool of edge devices while maximizing the overall accuracy. INFIDA has more general applicability, as we do not make any assumption on the network topology and also employ a more flexible cost model that takes into account network delay. Another related work in this framework is VideoEdge [32], which studies how to split the analytics pipeline across different computational clusters to maximize the average inference accuracy. Beside the focus on the specific video application, the paper does not propose any dynamic allocation placement algorithm.

Even if the problem of inference provisioning is currently overlooked in the context of distributed systems, there exists a vast literature on the problem of content placement [33] where objects can be stored (cached) into different nodes in order to reduce the operational cost of content delivery. Content placement has been extended to the case of service caching (or placement), where an entire service can be offloaded onto nodes co-located with base-stations or mobile-micro clouds, engaging not only storage but also computational resources and energy [34], [35]. The similarities between this problem and inference provisioning inspired us in the design of Inference Delivery Networks. However, the two problems feature crucial differences. First, in a content delivery network a request for a given item may only be served by a server storing that specific item. Whereas, in IDNs, several models can provide an answer but the accuracy of the answer can be different [36]. Second, a key property of content placement is that the service cost always increases together with the path length, i.e., the distance between the request source and the node serving the request. This is not the case for
inference delivery networks, as 1) upstream models may be more accurate and 2) the same model at may feature different processing delays based on the serving node properties. This leads to a more complex cost function, as the first node receiving the request may not be the optimal one to serve it. Note that this difference was crucial in the design of our algorithm (see Fig. 3). Finally, multiple requests can simultaneously be processed by a given model, leading to additional considerations about requests load and serving capacities.

A similar trade-off between resource usage and perceived quality typically emerges in the context of video caching [37]–[43], where the same video can be cached into multiple network nodes at different qualities (or resolutions): the operator optimizes the user experience by jointly deciding the placement of videos and their quality. These works either maximize the video quality perceived by the user [37], [40], [41], minimize the download time [39], [43] and the backhaul traffic [38], or minimize a combined cost [42]. Although some of the models in these papers may be adapted to inference provisioning in IDNs, these works in general study static optimization problems under a known request process and consider simple network topologies: a single cache [38], [39], a pool of parallel caches [37], [43], bipartite networks [41], [42]. The only exception is [40], which considers an arbitrary topology and provides some online heuristics, but ignores the service latency, which is of paramount importance when placing interactive ML models (e.g., for applications like augmented reality or autonomous driving). Instead, we propose a dynamic policy that jointly optimizes inference quality and latency and provides strong performance guarantees without requiring any knowledge about the request process thanks to our adversarial setting (Sec. V).

Adversarial analysis is typically studied through the lens of online convex optimization (OCO) [5]. OCO models can be tackled with well-understood learning algorithms [5], [44], [45]. However, the problem of optimizing the allocation of ML models in IDNs diverges from the template of OCO. In particular, the decision set and the cost functions are non-convex because the allocation decisions are not continuous; moreover, as we show in Appendix B, computing the optimal allocation is NP-hard contrarily to the OCO setting. In our work, we generalize findings from [46]–[48], and provide novel results to approximate budget-additive (submodular) set functions [49], which are of independent interest beyond this work (e.g., for online advertising [50], [51], market equilibrium [52], [53]).

Finally, ML model allocation in an IDN can also be considered as a particular instance of similarity caching [54], a general model where items and requests can be thought as embedded in a metric space: edge nodes can store a set of items, and the distance between a request and an item determines the quality of the matching between the two. Similarity caching was applied to a number of applications including content-based image retrieval [55], contextual advertising [56], object recognition [57], and recommender systems [58]. To the best of our knowledge, the literature on similarity caching has restricted itself to (i) a single cache (with the exception of [58]–[60]), and (ii) homogeneous items with identical resource

<figure><img src="image_1.png" /><figcaption>Fig. 1: System overview: a network of compute nodes serves inference requests along predefined routing paths. A repository node at the end of each path ensures that requests are satisfied even when there are no suitable models on intermediate nodes.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Fig. 2: Example of pre-trained model catalog for the image classification task. Data from [61].</figcaption></figure>

requirements. A consequence is that in our setting similarity caching policies would only allocate models based on their accuracy, ignoring the trade-offs imposed by their resource requirements. Moreover, the literature on similarity caching ignores system throughput constraints, while we explicitly take into account that each model can only serve a bounded number of requests per second, according to its capacity.

### III. INFERENCE SYSTEM DESIGN

We consider a network of compute nodes, each capable of hosting some pre-trained ML models depending on its capabilities. Such ML models are used to serve inference requests for different classification or regression tasks.<fnref n="1" /> As shown in Fig. 1, requests are generated by end devices and routed over given serving paths (e.g., from edge to cloud nodes). The goal of the system is to optimize the allocation of ML models across the network so that the aggregate serving cost is minimized. Our system model is detailed below, and the notation used across the paper is summarized in Table I.

#### A. Compute Nodes and Models

We represent the inference delivery network (IDN) as a weighted graph $G(\mathcal{V}, \mathcal{E})$, where $\mathcal{V}$ is the set of compute

<footnote n="1">¹We are using the term task according to its meaning in the ML community, e.g., a task could be to detect objects in an image, to predict its future position, to recognize vocal commands.</footnote>
<table><caption>TABLE I: Notation Summary.</caption><thead><tr><th colspan="2">Inference Delivery Networks</th></tr></thead><tbody><tr><td>G(V, E)</td><td>Undirected weighted graph, with nodes V and edges E</td></tr><tr><td>w<sub>u,v</sub></td><td>Weight of edge (u, v) ∈ E</td></tr><tr><td>N | M | M<sub>i</sub></td><td>Tasks / models catalog / models catalog for task i</td></tr><tr><td>s<sup>v</sup><sub>m</sub></td><td>Size of model m ∈ M on node v ∈ V</td></tr><tr><td>a<sub>m</sub></td><td>Prediction accuracy of model m</td></tr><tr><td>d<sup>v</sup><sub>m</sub></td><td>Average inference delay of model m at node v</td></tr><tr><td>b<sup>v</sup></td><td>Allocation budget constraint at node v</td></tr><tr><td>x<sup>v</sup><sub>m</sub></td><td>0-1 indicator variable set to 1 if model m at node v is allocated</td></tr><tr><td>ω<sup>v</sup><sub>m</sub></td><td>0-1 indicator constant set to 1 if model m at node v is a permanent repository model</td></tr><tr><td>ω / x</td><td>Minimal allocation vector (Sec. III-D) / allocation vector (Sec. III-A)</td></tr><tr><td>x<sup>v</sup> / x</td><td>Allocation vector at node v / global allocation vector</td></tr><tr><td>s<sup>v</sup></td><td>Vector of model sizes at node v</td></tr><tr><td>p</td><td>Routing path {p<sub>1</sub>, ..., p<sub>J</sub>} of connected nodes p<sub>j</sub> ∈ V</td></tr><tr><td>ν(p)</td><td>Repository node associated to the path p</td></tr><tr><td>C<sup>p<sub>j</sub></sup><sub>p,m</sub></td><td>Cost of serving at node p<sub>j</sub> along path p using model m</td></tr><tr><td>ρ</td><td>Request type ρ = (i, p), i.e., the requested task i and the request's routing path p</td></tr><tr><td>R/R<sub>i</sub></td><td>Set of all the possible request types / request types for task i</td></tr><tr><td>R</td><td>Total number of request types</td></tr><tr><td>r<sup>t</sup><sub>ρ</sub></td><td>Number of times ρ is requested during time slot t</td></tr><tr><td>load<sup>t,v</sup><sub>m</sub>(ρ')</td><td>Number of type-ρ requests served by model m at node v during the t-slot</td></tr><tr><td>L<sup>v</sup><sub>m</sub></td><td>Maximum capacity of model m on node v</td></tr><tr><td>l<sup>t,v</sup><sub>ρ,m</sub></td><td>Potential available capacity of m on node v for request type ρ at time t</td></tr><tr><td>γ<sup>k</sup><sub>ρ</sub></td><td>k-th smallest cost for request type ρ along its path (Sec. III-E).</td></tr><tr><td>r<sub>t</sub> / l<sub>t</sub></td><td>Request batch vector / potential available capacities vector</td></tr><tr><td>λ<sup>k</sup><sub>ρ</sub>(l<sub>t</sub>) / z<sup>k</sup><sub>ρ</sub>(l<sub>t</sub>, x)</td><td>Potential / effective available capacity of the model serving ρ with cost γ<sup>k</sup><sub>ρ</sub> (Sec. III-E).</td></tr><tr><td>Z<sup>k</sup><sub>ρ</sub>(r<sub>t</sub>, l<sub>t</sub>, x)</td><td>Number of requests of type ρ that can be served by the k-th smallest cost models found along path p under allocation x at time slot t (15)</td></tr><tr><td>κ<sub>ρ</sub>(v, m)</td><td>Rank of model m allocated to node v, among all the models that can potentially serve requests of type ρ</td></tr><tr><td>K<sub>ρ</sub></td><td>Maximum number of models that request ρ may encounter along its serving path</td></tr><tr><td>X<sup>v</sup> / X</td><td>Set of valid integral allocations at node v / at all nodes</td></tr><tr><td>C / G / G<sub>T</sub></td><td>The overall system cost (12) / allocation gain (13) / static allocation gain (14)</td></tr></tbody></table><table><thead><tr><th colspan="2">INFIDA</th></tr></thead><tbody><tr><td>Φ<sup>v</sup></td><td>Weighted negative entropy map Φ<sup>v</sup> : D<sup>v</sup> → R given as Φ<sup>v</sup>(y<sup>v</sup>) = &sum;<sub>m&isin;M</sub> s<sup>v</sup><sub>m</sub>y<sup>v</sup><sub>m</sub> log(y<sup>v</sup><sub>m</sub>)</td></tr><tr><td>Φ</td><td>Global mirror map Φ(y) = &sum;<sub>v&isin;V</sub> Φ<sup>v</sup>(y<sup>v</sup>)</td></tr><tr><td>y<sup>v</sup><sub>m</sub></td><td>Fractional allocation (allocation probability) of model m at node v</td></tr><tr><td>y<sup>v</sup> / y</td><td>Fractional allocation vector at node v / global fractional allocation vector</td></tr><tr><td>y<sup>v</sup> / Y</td><td>Set of valid fractional allocations at node v / at all nodes</td></tr><tr><td>g<sub>t</sub></td><td>Subgradient vector of G over Y at point y<sub>t</sub></td></tr><tr><td>g<sup>v</sup><sub>t,m</sub></td><td>Component (v, m) of the subgradient vector g<sub>t</sub></td></tr><tr><td>P<sup>Φ<sup>v</sup></sup><sub>Y<sup>v</sup>&cap;D<sup>v</sup></sub>(·)</td><td>Projection operator onto Y<sup>v</sup> &cap; D<sup>v</sup></td></tr><tr><td>&eta;</td><td>learning rate</td></tr><tr><td>B</td><td>Refresh period</td></tr><tr><td>T</td><td>Time horizon</td></tr><tr><td>t</td><td>Time slot t &isin; [T]</td></tr><tr><td>A</td><td>Regret constant</td></tr><tr><td>L<sub>max</sub></td><td>Upper bound on model capacities</td></tr><tr><td>&Delta;<sub>C</sub></td><td>Upper bound on maximum serving cost difference</td></tr><tr><td>&psi;</td><td>Regret discount factor equal to 1 &minus; 1/c</td></tr></tbody></table>

nodes, and $\mathcal{E}$ represents their interconnections. Each node $v \in V$ is capable of serving inference tasks that are requested from anywhere in the network (e.g., from end-users, vehicles, IoT devices). We denote by $\mathcal{N} = \{1, 2, ..., |\mathcal{N}|\}$ the set of tasks the system can serve (e.g., object detection, speech recognition, classification), and assume that each task $i \in \mathcal{N}$ can be served with different quality levels (e.g., different accuracy as illustrated in Fig. 2) and different resources' requirements by a set of suitable models $M_i$. Each task is served by a separate set of models, i.e., $M_i \cap M_{i'} = \emptyset, \forall i, i' \in \mathcal{N}, i \neq i'$. Catalog $M_i$ may encompass, for instance, independently trained models or shrunk versions of a high quality model generated through distillation [62], [63]. We denote by $\mathcal{M} = \cup_{i \in \mathcal{N}} M_i = \{1, 2, ..., |\mathcal{M}|\}$ the catalog of all the available models.

Finally, every model of the catalog may provide a different throughput (i.e., number of requests it can serve in a given time period), and therefore, support a different load (we formalize this in Sec. III-D).

For each compute node $v \in V$, we denote by

$$x_m^v \in \{0, 1\}, \text{ for } m \in \mathcal{M}, \quad (1)$$

the decision variable that indicates if model $m \in \mathcal{M}$ is deployed on node $v$.<fnref n="2" /> Therefore, $x^v = [x_m^v]_{m \in \mathcal{M}}$ is the allocation vector on node $v$, and $x = [x^v]_{v \in V}$ denotes the global allocation decision.

We assume that the allocation of ML models at each node is constrained by a single resource dimension, potentially different at each node. A node could be, for instance, severely limited by the amount of available GPU memory, another by the maximum throughput in terms of instructions per second. The limiting resource determines the allocation budget $b^v \in \mathbb{R}_+$ at node $v \in V$. We also denote with $s_m^v \in \mathbb{R}_+$ the size of model $m \in \mathcal{M}$, i.e., the consumed amount of the limiting resource at node $v$.<fnref n="3" /> Therefore, budget constraints are expressed as

$$\sum_{m \in \mathcal{M}} x_m^v s_m^v \leq b^v, \forall v \in V. \quad (2)$$

To every task $i \in \mathcal{N}$, we associate a fixed set of repository nodes that always run one model capable of serving all the requests for task $i$ (e.g., high-performance models deployed in large data centers). We call these models repository models and they are statically allocated. Repository models ensure requests are satisfied even when the rest of the network is not hosting any additional model.

We discern repository models through constants $\omega_m^v \in \{0, 1\}$, each indicating if model $m$ is permanently deployed on node $v$. We assume that values $\omega_m^v$ are given as input. We call the vector $\omega = [\omega_m^v]_{(v,m) \in V \times \mathcal{M}}$ the minimal allocation. Note that the presence of repositories introduce the following

<footnote n="2">²Our formulation allows each node to host multiple copies of the same model to satisfy a larger number of request. For example, two copies of the same model can be represented as two distinct models with identical performance and requirements.</footnote>

<footnote n="3">³Note that, even when the limiting resource is the same, say computing, the budget consumed by a model may be different across nodes, as they may have different hardware (e.g., GPUs, CPUs, or TPUs).</footnote>
constraints to the allocation vector:

$$x_m^v \ge \omega_m^v, \forall v \in \mathcal{V}, \forall m \in \mathcal{M}. \quad (3)$$

The set of possible allocations at node $v \in \mathcal{V}$ is determined by the integrality constraints (1), budget constraints (2), and repository constraints (3), i.e.,

$$\mathcal{X}^v \triangleq \{\mathbf{x}^v \in \{0,1\}^M : \mathbf{x}^v \text{ satisfies Eqs. (1)-(3)}\}. \quad (4)$$

The set of possible global allocations is given as $\mathcal{X} \triangleq \mathcal{X}_{v \in \mathcal{V}}^v$.

## B. Inference Requests

We assume that every node has a predefined routing path towards a suitable repository node for each task $i \in \mathcal{N}$. Therefore, for a given request for task $i$, the routing path is a set of network nodes towards a repository node able to serve task $i$. Since we assume repository nodes are predefined, the routing path does not depend on the placement decisions (i.e., on the variables $x_m^v$). Hence, a request always follows its predetermined path, but intermediate nodes that host suitable models can serve it directly instead of forwarding it all the way to the repository node. In such cases, the request would traverse just a portion of the path. A routing path $\mathbf{p}$ of length $|\mathbf{p}| = J$ is a sequence $\{p_1, p_2, \dots, p_J\}$ of nodes $p_j \in \mathcal{V}$ such that edge $(p_j, p_{j+1}) \in \mathcal{E}$ for every $j \in \{1, 2, \dots, J-1\}$. As in [46], we assume that paths are simple, i.e., they do not contain repeated nodes. A request is therefore characterized by the pair $(i, \mathbf{p})$, where $i$ is the task requested and $\mathbf{p}$ is the routing path to be traversed. We call the pair $(i, \mathbf{p})$ the request type. We denote by $\mathcal{R}$ the set of all possible request types, and by $\mathcal{R}_i$ all possible request types for tasks $i$. When a request for task $i$ is propagated from node $p_1$ toward the associated repository node $\nu(\mathbf{p}) \triangleq p_j$, any intermediate node along the path that hosts a suitable model $m \in \mathcal{M}_i$ can serve it. The actual serving strategy is described in Sec. III-E.

## C. Cost Model

When serving a request of type $\rho=(i, \mathbf{p}) \in \mathcal{R}$ on node $p_j$ using model $m$, the system experiences an inference cost that depends on the quality of the model (i.e., on inference inaccuracy) and the inference time.<fnref n="4" /> Additionally, the system experiences a network cost, due to using the path between $p_1$ and $p_j$. Similarly to previous work [64], we can write the total cost of serving a request as

$$C_{\mathbf{p},m}^{p_j} = f((p_1, \dots, p_j), m). \quad (5)$$

While our theoretical results hold under this very general cost model, in what follows—for the sake of concreteness—we refer to the following simpler model:

$$C_{\mathbf{p},m}^{p_j} = \sum_{j'=1}^{j-1} w_{p_{j'}, p_{j'+1}} + d_m^{p_j} + \alpha(1-a_m), \quad (6)$$

where $a_m$ and $d_m^{p_j}$ are respectively the prediction accuracy (in a scale from 0 to 1) and the average inference delay of model $m$ on node $p_j$. Indeed, the same model may provide different inference delays, depending on the hardware capabilities of the node on which it is deployed, e.g., the type of GPU or TPU [65]. Parameter $w_{v,v'} \in \mathbb{R}_+$ is the (round-trip) latency of edge $(v, v') \in \mathcal{E}$. Parameter $\alpha$ weights the importance of accuracy w.r.t. the overall latency and can be set depending on the application. Note that seeking cost minimization along a serving path usually leads to a trade-off: while the network cost always increases with $j$, in a typical network the service cost $d_m^{p_j} + \alpha(1-a_m)$ tends to decrease, as farther nodes (e.g., data centers) are better equipped and can run more accurate models (Fig. 2). We remark that models' sizes determine which allocations are feasible, but do not affect directly the service costs. As a consequence, even if we assume different limiting resources on different nodes (e.g., GPU, memory), we do not need to convert amounts of different resources to a common unit (e.g., a monetary cost).

## D. Request Load and Serving Capacity

Let us assume that time is split in slots of equal duration. We consider a time horizon equal to $T$ slots. At the beginning of a slot $t$, the system receives a batch of requests $\mathbf{r}_t = [r_{\rho}^t]_{\rho \in \mathcal{R}}$, where $r_{\rho}^t \in \mathbb{N} \cup \{0\}$ denotes the number of requests of type $\rho \in \mathcal{R}$.

Model $m \in \mathcal{M}$ has maximum capacity $L_m^v \in \mathbb{N}$ when deployed at node $v \in \mathcal{V}$, i.e., it can serve at most $L_m^v$ requests during one time slot $t \in [T]$, in absence of other requests for other models. We do not make specific assumptions on the time required to serve a request.

We denote by $l_{\rho,m}^{t,v} \in \mathbb{N} \cup \{0\}$ the potential available capacity, defined as the maximum number of type-$\rho$ requests node $v$ can serve at time $t$ through model $m$, under the current request load $\mathbf{r}_t$ and allocation vector $\mathbf{x}_t^v$. Formally, let $\text{load}_m^{t,v}(\rho)$ denote the number of type-$\rho$ requests served by model $m$ at node $v$ during the $t$-slot, then

$$l_{\rho,m}^{t,v} \triangleq \min \left\{ L_m^v - \sum_{\rho' \in \mathcal{R} \setminus \{\rho\}} \text{load}_m^{t,v}(\rho'), r_{\rho}^t \right\}. \quad (7)$$

The potential available capacity depends on the request arrival order and the scheduling discipline at node $v$. For instance, suppose that in time slot $t$, requests of two types $\rho = (i, \mathbf{p})$ and $\rho' = (i, \mathbf{p}')$ arrive at node $v$. The arrival order and the node scheduling discipline may determine that many requests of type $\rho'$ be served, which would leave a small $l_{\rho,m}^{t,v}$ available for requests of type $\rho$. Or the opposite may happen. It is useful to define the potential available capacity also for models that are not currently deployed at the node, as $l_{\rho,m}^{t,v} \triangleq \min\{L_m^v, r_{\rho}^t\}$. The effective available capacity is then equal to $l_{\rho,m}^{t,v}x_m^v$.

Our analysis in Sec. V considers a “pessimistic” scenario where an adversary selects both requests and available capacities for all models but the repository ones. This approach relieves us from the need to model system detailed operations, while our proposed algorithm (Sec. IV) benefits from strong

<footnote n="4">⁴Note that deployed models may need to be re-trained from time to time, and we do not consider the corresponding costs. Moreover, to streamline the presentation, we assume the inference costs to be static over time; nonetheless, one could easily extend our model to the case in which these costs are time-varying.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 3: Necessity of partial synchronization in IDN among close-by computing nodes under the cost model in Eq. (6).</figcaption></figure>

guarantees in the adversarial setting. In what follows, we can then consider that values $l_{\rho,m}^{t,v}$ are exogeneously determined. The vector of potential available capacities at time $t \in [T]$ is denoted by

$$ \mathbf{l}_t = [l_{\rho,m}^{t,v}]_{(\rho,m,v) \in \bigcup_{i \in \mathcal{N}} \mathcal{R}_i \times \mathcal{M}_i \times \mathcal{V}}. \quad (8) $$

As we mentioned in Sec. III-A, any request of type $\rho = (i, \mathbf{p}) \in \mathcal{R}$ can always be served by the associated repository model at node $\nu(\mathbf{p})$. This requirement can be expressed as follows:

$$ \sum_{\rho \in \mathcal{R}_i} r_{\rho}^{t} \le \sum_{m \in \mathcal{M}_i} \omega_m^{\nu(\mathbf{p})} L_m^{\nu(\mathbf{p})}, \forall i \in \mathcal{N}. \quad (9) $$

Thus, at any time $t \in [T]$ the adversary can select a request batch $\mathbf{r}_t$ and potential available capacity $\mathbf{l}_t$ from the set

$$ \mathcal{A} \triangleq \left\{ \begin{array}{l} (\mathbf{r}, \mathbf{l}) \in (\mathbb{N} \cup \{0\})^{\mathcal{R}} \times (\mathbb{N} \cup \{0\})^{\mathcal{M}_i \times \mathcal{V}} \cup_{i \in \mathcal{N}} \mathcal{R}_i \times \mathcal{M}_i \times \mathcal{V}: \\ \sum_{\rho \in \mathcal{R}_i} r_{\rho}^{t} \le \sum_{m \in \mathcal{M}_i} \omega_m^{\nu(\mathbf{p})} L_m^{\nu(\mathbf{p})}, l_{\rho,m}^{v} \le \min\{L_m^{v}, r_{\rho}\}, \\ \forall i \in \mathcal{N}, v \in \mathcal{V}, m \in \mathcal{M}, \rho \in \mathcal{R} \end{array} \right\}. \quad (10) $$

Note the constraint on potential available capacities is looser than the definition in (7) corresponding to a more powerful adversary.

## E. Serving Model

Given request $\rho=(i, \mathbf{p}) \in \mathcal{R}$, let $K_{\rho} = |\mathbf{p}||\mathcal{M}_i|$ denote the maximum number of models that may encounter along its serving path $\mathbf{p}$. We order the corresponding costs $\{C_{\mathbf{p},m}^{p_j}, \forall m \in \mathcal{M}_i, \forall p_j \in \mathbf{p}\}$ in increasing order and we denote by $\kappa_{\rho}(v, m)$ the rank of model $m \in \mathcal{M}_i$ allocated at node $v$ within the order defined above.<fnref n="5" /> If $v \notin \mathbf{p}$ we have $\kappa_{\rho}(v, m) = \infty$.

If $\kappa_{\rho}(v, m) = k$, then model $m$ at node $v$ has the $k$-th smallest cost to serve request $\rho$. We denote the model service cost, its potential available capacity, and its effective capacity as $\gamma_{\rho}^k$, $\lambda_{\rho}^k(l_t)$, and $z_{\rho}^k(l_t, \mathbf{x})$, respectively:

$$ \gamma_{\rho}^{k} = C_{\mathbf{p},m}^{v}, \quad \lambda_{\rho}^{k}(l_t) = l_{\rho,m}^{t,v}, \quad z_{\rho}^{k}(l_t, \mathbf{x}) = x_{m}^{v} l_{\rho,m}^{t,v}. \quad (11) $$

We assume the IDN serves requests as follows. Each request is forwarded along its serving path and served when it

encounters a model with the smallest serving cost among those that are not yet saturated, i.e., that may still serve requests.

Since models do not necessarily provide increasing costs along the path, this serving strategy requires that a node that runs a model $m \in \mathcal{M}_i$ and receives a request for task $i$, knows whether there are better alternatives for serving task $i$ upstream or not. In the first case, it will forward the request along the path, otherwise it will serve it locally. We argue that, in a real system, this partial knowledge can be achieved with a limited number of control messages. In fact, if node $v = p_h$ hosts the model with the $k$-th cost for request $(i, \mathbf{p})$, it only needs information about those models that $(i)$ are located upstream on the serving path (i.e., on nodes $p_l \in \mathbf{p}$ with $l > h$), and $(ii)$ provide a cost smaller than $\gamma_{\rho}^k$. Since the cost increases with the network latency (as illustrated in Fig. 3), the number of models satisfying these criteria is small in practice.<fnref n="6" /> A node needs to propagate downstream a control message with the information about the requests it can serve and the corresponding costs. Nodes forwarding the control message progressively remove the information about the tasks they can serve with a smaller cost, until the control message payload is empty and the message can be dropped. Every node $v \in \mathcal{V}$ generates this control message whenever the available capacity of any of the local models in $v$ changes.

According to the presented serving strategy, the requests load is split among the currently available models giving priority to those that provide the smallest serving costs up to their saturation. In particular, model $m$ with the $k$-th smallest cost will serve some requests of type $\rho$ only if the less costly models have not been able to satisfy all of such requests (i.e., if $\sum_{k'=1}^{k-1} z_{\rho}^{k'}(l_t, \mathbf{x}) < r_{\rho}^t$). If this is the case, model $m$ will serve with cost $\gamma_{\rho}^k$ at most $z_{\rho}^k(l_t, \mathbf{x})$ requests (its effective available capacity) out of the $r_{\rho}^t - \sum_{k'=1}^{k-1} z_{\rho}^{k'}(l_t, \mathbf{x})$ requests still to be satisfied. The aggregate cost incurred by the system at time slot $t$ is then given by

$$ C(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}) = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}} \gamma_{\rho}^{k} \cdot \min \left\{ r_{\rho}^{t} - \sum_{k'=1}^{k-1} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}), z_{\rho}^{k}(\mathbf{l}_t, \mathbf{x}) \right\} \cdot \mathbf{1}_{\left\{ \sum_{k'=1}^{k-1} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) < r_{\rho}^{t} \right\}}. \quad (12) $$

Note that we introduce the $\min\{\cdot, \cdot\}$ operator, since the number of requests served by the $k$-th best model cannot exceed its effective capacity $z_{\rho}^k(l_t, \mathbf{x})$. We add the indicator function $\mathbf{1}_{\{\cdot\}}$ to indicate that the $k$-th best model does not serve any requests, in case better models (ranked from 1 to $k-1$) are able to satisfy all of them.

## F. Allocation Gain and Static Optimal Allocations

We are interested in model allocations $\mathbf{x}$ that minimize the aggregate cost (12), or, equivalently, that maximize the allocation gain defined as

$$ G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}) = C(\mathbf{r}_t, \mathbf{l}_t, \omega) - C(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}). \quad (13) $$

The first term $C(\mathbf{r}_t, \mathbf{l}_t, \omega)$ on the right hand side is the service cost when only repository models are present in

<footnote n="5">⁵Note that we do not consider only the models deployed in the network, but all the possible node-model pairs.</footnote>

<footnote n="6">⁶In realistic settings (Sec. VI), we experienced that each deployed model has at most 6 better alternatives on upstream nodes (worst case with $\alpha=1$).</footnote>
the network. Since intermediate nodes can help serving the requests at a reduced cost, $C(\mathbf{r}_t, \mathbf{l}_t, \omega)$ is an upper bound on the aggregate serving cost, and the allocation gain captures the cost reduction achieved by model allocation $\mathbf{x}$.

The static model allocation problem can then be formulated as finding the model allocation $\mathbf{x}^*$ that maximizes the time-averaged allocation gain over the time horizon $T$, i.e.,

$$ \mathbf{x}^* = \arg\max_{\mathbf{x} \in \mathcal{X}} \left( G_T(\mathbf{x}) \triangleq \frac{1}{T} \sum_{t=1}^{T} G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}) \right). \quad (14) $$

This is a submodular maximization problem under multiple knapsack constraints. In our context, this intuitively means that the problem is characterized by a diminishing return property: adding a model $m$ to any node $v$ gives us a marginal gain that depends on the current allocation: the more the models already deployed in the current allocation, the less the marginal gain we get by the new $m$. We prove submodularity in Lemma A.2 in Appendix A. In Appendix B Theorem B.1, we prove that this problem is NP-hard even under cardinality constraints (i.e., the models have equal size) and a two nodes scenario. We demonstrate the hardness of the problem by a reduction of the similarity caching problem [54], [66], which is NP-hard; a result that follows from a reduction of the dominating set problem. It is known that submodular maximization problems cannot be approximated with a ratio better than $(1 - 1/e)$ even under simpler cardinality constraints [67]. Under the multi-knapsack constraint, it is possible to solve the offline problem achieving a $(1 - 1/e - \epsilon)$-approximation through a recent algorithm proposed in [68].

Let us consider a model allocation $\mathbf{x}$. Within time slot $t$, the $k$ smallest cost models along a path $\mathbf{p}$ that are suitable for request type $\rho = (i, \mathbf{p})$ can serve up to $Z_{\rho}^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x})$ requests, where $Z_{\rho}^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x})$ is defined as

$$ Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}) \triangleq \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_{t}, \mathbf{x}) \right\}. \quad (15) $$

The min{·, ·} operator denotes that we can never serve more than the number of requests $r_{\rho}^t$ issued by users. Observe that, being the minimal allocation $\omega$ an input parameter not dependent on our decisions, $Z_{\rho}^k(\mathbf{r}_t, \mathbf{l}_t, \omega)$ is a constant. Additionally, since the models allocated in $\mathbf{x}$ always include those allocated in $\omega$, we have $Z_{\rho}^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}) \ge Z_{\rho}^k(\mathbf{r}_t, \mathbf{l}_t, \omega)$.

Using (15), we provide the following alternative formulation of the allocation gain.

**Lemma III.1.** The allocation gain (13) has the following equivalent expression:

$$ G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}) = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} \underbrace{(\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k})}_{\text{cost saving}} \underbrace{(Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}) - Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \omega))}_{\text{additional requests}}. \quad (16) $$

We prove this lemma in Appendix C. This result tells us that the gain of a certain allocation $\mathbf{x}$ can be expressed as a sum of several components. In particular, for each request type $\rho$, the $k$-th smallest cost model along the path contributes to

**Algorithm 1 INFIDA distributed allocation on node v**

1: **procedure** INFIDA($\mathbf{y}_1^v = \arg \min_{\mathbf{y}^v \in \mathcal{Y}^v \cap \mathcal{D}^v} \Phi^v(\mathbf{y}^v)$), $\mathbf{x}_1^v = \text{DEPROUND}(\mathbf{y}_1^v)$,
   $\eta \in \mathbb{R}_+$
2: **for** $t = 1, 2, \dots, T$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Compute $\mathbf{g}_t^v \in \partial \mathbf{y}^v G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t)$ through (18).
4: &nbsp;&nbsp;&nbsp;&nbsp;$\hat{\mathbf{y}}_t^v \leftarrow \nabla \Phi^v(\mathbf{y}_t^v)$ ▷ Map state to the dual space
5: &nbsp;&nbsp;&nbsp;&nbsp;$\hat{\mathbf{h}}_{t+1} \leftarrow \hat{\mathbf{y}}_t^v + \eta \mathbf{g}_t^v$ ▷ Take gradient step in the dual space
6: &nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{h}_{t+1}^v \leftarrow (\nabla \Phi^v)^{-1}(\hat{\mathbf{h}}_{t+1})$ ▷ Map dual state back to the primal space
7: &nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{y}_{t+1}^v \leftarrow \mathcal{P}_{\mathcal{Y}^v \cap \mathcal{D}^v}^{\Phi^v}(\mathbf{h}_{t+1}^v)$ ▷ Project new state onto the feasible region using Algorithm 2
8: &nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{x}_{t+1}^v \leftarrow \text{DEPROUND}(\mathbf{y}_{t+1}^v)$ ▷ Sample a discrete allocation

the gain with a component (i) proportional to its cost saving $\gamma_{\rho}^{k+1} - \gamma_{\rho}^k$ with respect to the $(k+1)$-th smallest cost model and (ii) proportional to the amount of additional requests that the $k$-th smallest cost models in allocation $\mathbf{x}$ can serve with respect to the minimal allocation $\omega$.

## IV. INFIDA ALGORITHM

In this section, we propose INFIDA, an online algorithm that can operate in a distributed fashion without requiring global knowledge of the allocation state and requests arrival. In Sec. V, we show that INFIDA generates dynamically allocations experiencing average costs that converge to a $(1 - 1/e - \epsilon)$-approximation of the optimum, which matches the best approximation ratio achievable in polynomial time even in this online setting.

### A. Algorithm Overview

On every node $v \in \mathcal{V}$, INFIDA updates the allocation $\mathbf{x}^v \in \mathcal{X}^v \subset \{0, 1\}^{\mathcal{M}}$, by operating on a correspondent fractional state $\mathbf{y}^v \in \mathcal{Y}^v \subset [0, 1]^{\mathcal{M}}$, and the fractional allocations satisfy the budget constraint in Eq. (2). Note that, if $\|\mathbf{s}^v\|_1 < b^v$ for a node $v \in \mathcal{V}$, we can always consider fractional allocations that consume entirely the allowed budget; otherwise, all the allocations are set to 1 (node $v$ can store the whole catalog of models). Formally, if $\|\mathbf{s}^v\|_1 \ge b^v$ then

$$ \mathcal{Y}^v \triangleq \left\{ \mathbf{y}^v \in [0, 1]^M : \sum_{m \in \mathcal{M}} y_m^v s_m^v = b^v, \forall v \in V \right\}; \quad (17) $$

otherwise, for the corner case $\|\mathbf{s}^v\|_1 < b^v$, we have $\mathcal{Y}^v \triangleq \{[1]^M\}$.

Each variable $y_m^v$ can be interpreted as the probability of hosting model $m$ on node $v$, i.e., $y_m^v = \mathbb{P}[x_m^v = 1] = \mathbb{E}[x_m^v]$.

We define $G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y})$ as in (13), replacing $\mathbf{x}$ with $\mathbf{y}$. Note that $G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y})$ is a concave function of variable $\mathbf{y} \in \mathcal{Y} = \times_{v \in \mathcal{V}} \mathcal{Y}^v$ (see Lemma F.1 in Appendix F).

Within a time slot $t$, node $v$ collects measurements from messages that have been routed through it (Sec. IV-B). At the end of every time slot, the node (i) computes its new fractional state $\mathbf{y}^v$, and (ii) updates its local allocation $\mathbf{x}^v$ via randomized rounding (Sec. IV-C). INFIDA is summarized in Algorithm 1 and detailed below.

**State computation.** The fractional state $\mathbf{y}^v$ is updated through an iterative procedure aiming to maximize $G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y})$. This could be the standard gradient ascent method, which updates
the fractional state at each node as $\mathbf{y}_{t+1}^v = \mathbf{y}_t^v + \eta \mathbf{g}_t^v$, where $\eta \in \mathbb{R}_+$ is the step size and $\mathbf{g}_t^v$ is a subgradient of $G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y})$ with respect to $\mathbf{y}^v$.

In our work, we use a generalized version of the gradient method called Online Mirror Ascent (OMA) [69, Ch. 4]. OMA uses a function $\Phi^v: D^v \to \mathbb{R}_+$ (mirror map) to map $\mathbf{y}$ to a dual space before applying the gradient ascent method; then the obtained state is mapped back to the primal space (lines 3–5 of Algorithm 1). OMA reduces to the classic gradient ascent method if $\Phi^v$ is the squared Euclidean norm (in this case the primal space coincides with the dual one). Instead, we use the weighted negative entropy map $\Phi^v(\mathbf{y}^v) = \sum_{m \in \mathcal{M}} s_m^v y_m^v \log(y_m^v)$, which is known to achieve better convergence rate in high dimensional spaces when each subgradient component is bounded.<fnref n="7" /> To compute a feasible fractional state $\mathbf{y}^v$, we then perform a projection to the set $\mathcal{Y}^v$ on node $v$ (line 6 of Algorithm 1). We adapt the projection algorithm from [70] to obtain a negative entropy projection $\mathcal{P}_{\mathcal{Y}^v \cap D^v}(\cdot)$. Our adaptation is described in Appendix D.

**Allocation update.** Once the fractional state $\mathbf{y}^v$ has been updated, the final step of INFIDA is to determine a new random discrete allocation $\mathbf{x}^v$ and update the local models accordingly. The sampled allocation $\mathbf{x}^v$ should (i) comply with the budget constraint (2) on node $v$ and (ii) be consistent with the fractional state, i.e., $\mathbb{E}[x_m^v] = y_m^v \forall m \in \mathcal{M}$. To this purpose, we use the DepRound [71] subroutine (line 7 of Algorithm 1).

In the remainder of this section we detail how each node computes its contribution to the global subgradient, and the rounding strategy used to determine the discrete allocation.

## B. Subgradient Computation

At the end of every time slot $t$, a subgradient $\mathbf{g}_t$ of the gain function in Eq. (16) at point $\mathbf{y}_t \in \mathcal{Y}$ is computed in a distributed fashion: each node $v$ evaluates the $(v, m)$-th component of the subgradient for any $m \in \mathcal{M}$ as follows (see Appendix E):

$$g_{t,m}^v = \sum_{\rho \in \mathcal{R}} l_{\rho,m}^{t,v} \cdot (\gamma_{\rho}^{K_{\rho}^*(\mathbf{y}_t)} - C_{\mathbf{p},m}^v) \cdot \mathbf{1}_{\{\kappa_{\rho}(v,m) < K_{\rho}^*(\mathbf{y}_t)\}}, \quad (18)$$

where $K_{\rho}^*(\mathbf{y}_t)$ is the order of the worst needed model, i.e., the model with the highest cost that is needed to serve all the $r_{\rho}^t$ requests in the batch given the fractional state $\mathbf{y}_t$. Formally, $K_{\rho}^*(\mathbf{y}_t) \triangleq \min \{k \in [K_{\rho} - 1] : \sum_{k'=1}^k z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{y}_t) \geq r_{\rho}^t\}$.

For the sake of clarity assume that the mirror map is Euclidean, and then the dual and primal spaces collapse and $\hat{\mathbf{y}}_t = \mathbf{y}_t$. At each iteration, each component $y_m^{t,v}$ of the fractional allocation vector is updated by adding a mass equal to the product of $\eta > 0$ and the corresponding component of the subgradient (Algorithm 1, line 5). Observe that $g_{m,t}^v$ is the sum of different contributions, one per each request type $\rho$. Thanks to the indicator function, only the terms of the

request types that are served by model $m$ on $v$ contribute to $g_{m,t}^v$. This contribution is proportional to the potential available capacity of model $m$ on node $v$ and to the relative gain $(\gamma_{\rho}^{K_{\rho}^*(\mathbf{y}_t)} - C_{\mathbf{p},m}^v)$, i.e., the cost reduction achieved when serving request type $\rho$ with model $m$ on $v$, rather than with the worst needed model. Then, gradient updates add more mass to the models that can contribute more to increase the gain. On the contrary, the projection step tends to remove the added mass from all components to satisfy the constraints. The overall effect is that fractional allocations of more (resp. less) useful models tend to increase (resp. decrease).

The subgradient in Eq. (18) can be computed at each node using only information from the control messages collected at the end of the time slot $t$. The steps needed to compute the subgradient are as follows.

1. At the end of the time slot, each node generates a control message for every received request type $\rho = (i, \mathbf{p})$ that is propagated along $\mathbf{p}$. The control message contains the quantity $r_{\rho}^t \geq 1$ (the multiplicity of the request), and a cumulative counter $Z$ initialized to zero.

2. As the control message travels upstream, intermediate nodes add to $Z$ the local values $z_{\rho}^k(l_t, \mathbf{y}_t)$ (fractional effective capacity in Eq. (11)). These values are added following increasing values of cost. This message is propagated until $Z \geq r_{\rho}^t$, that is until the message reaches the $K_{\rho}^*(\mathbf{y}_t)$-th model.

3. Once the $K_{\rho}^*(\mathbf{y}_t)$-th model is detected, a control message is sent down in the opposite direction, containing the cost $\gamma_{\rho}^{K_{\rho}^*(\mathbf{y}_t)}$ of the last checked model. Every node $v$ in the reverse direction reads the cost value from the control message and, for each model $m \in \mathcal{M}_i$, computes the quantity

$$h_m^v = l_{\rho,m}^{t,v} \cdot (\gamma_{\rho}^{K_{\rho}^*(\mathbf{y}_t)} - C_{\mathbf{p},m}^v). \quad (19)$$

4. Node $v$ can then compute $g_{v,m}^v$ in Eq. (18) as follows

$$g_{t,m}^v = \sum_{m \in \mathcal{M}_i} h_m^v.$$

Note that the cost in Eq. (6) does not necessarily increase along the path. Therefore, a traversed node is not able to update directly the variable $Z$ when there exist upstream nodes with lower cost. In this case, the node simply appends the information $(z_{\rho}^k(l_t, \mathbf{y}_t), \gamma_{\rho}^k)$ to the message, and lets upstream nodes to apply any pending update in the correct order. In our work, we assume to operate on a reliable communication channel. Nonetheless, we note that INFIDA is robust to noise $\xi_t$ affecting the sub-gradient, as long as such noise is not biased, i.e., $\mathbb{E}[\xi_t] = 0$ for every timeslot $t$ [5, Theorem 3.4].

## C. State Rounding

Once the new fractional state $\mathbf{y}_{t+1}$ is computed, each node $v$ independently draws a random set of models to store locally in such a way that $\mathbb{E}[x_{t+1}^v] = \mathbf{y}_{t+1}^v$. This sampling guarantees that the final allocation $x_{t+1}^v$ satisfies constraint (2) in expectation. A naive approach is to draw each variable $x_m^{v,t+1}$ independently, but it leads to a large variance of the

<footnote n="7">⁷Technically, the advantage in this setting derives from the infinite norm of the subgradient being independent from the space dimension, while the Euclidean norm grows proportionally to the squared root of the space dimension [69, Sec. 4.3].</footnote>
total size of the models selected, potentially exceeding by far
the allocation budget at node $v$.

To construct a suitable allocation we adopt the DepRound procedure from [71]. The procedure modifies the fractional state $\mathbf{y}_{t+1}^v$ iteratively: at each iteration, DepRound operates on two fractional variables $y_m^{v,t+1}$, $y_{m'}^{v,t+1}$ so that at least one of them becomes integral and the aggregate size of the corresponding models $s_m^v y_m^{v,t+1} + s_{m'}^v y_{m'}^{v,t+1}$ does not change. This operation is iterated until all variables related to node $v$ are rounded except (at most) one, which we call residual fractional variable. This is done in $O(|\mathcal{M}|)$ steps.

Note that, to satisfy $\mathbb{E}[\mathbf{x}_{t+1}^v] = \mathbf{y}_{t+1}^v$, the residual fractional variable, say it $y_m^{v,t+1}$, needs to be rounded. At this point $x_m^{v,t+1}$ can be randomly drawn. Now the final allocation can exceed the budget bound $b^v$ by at most $s_{\bar{m}}$. These (slight) occasional violations of the constraint may not be a problem, e.g., at an edge server running multiple applications, where resources may be partially redistributed across different applications; they may be explicitly accounted for in the service level agreements. If the budget bound cannot be exceeded even temporarily, the node is not able to store the model $\bar{m}$, but it may still exploits the residual free resources to deploy the model that provides the best marginal gain among those that fit the available budget. In practice, we expect the corresponding gain decrease to be negligible.

V. THEORETICAL GUARANTEES

We provide the optimality guarantees of our INFIDA algorithm in terms of the $\psi$-regret [72]. In our scenario, the $\psi$-regret is defined as the gain loss in comparison to the best static allocation in hindsight, i.e., $\mathbf{x}_* \in \arg\max_{\mathbf{x} \in \mathcal{X}} \sum_{t=1}^T G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x})$, discounted by a factor $\psi \in (0, 1]$. Formally,

$$ \psi\text{-Regret}_{T, \mathcal{X}} \triangleq \sup_{\{\mathbf{r}_t, \mathbf{l}_t\}_{t=1}^T \in \mathcal{A}^T} \left\{ \psi \sum_{t=1}^T G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}_*) - \mathbb{E} \left[ \sum_{t=1}^T G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}_t) \right] \right\}, \quad (20) $$

where allocations $\mathbf{x}_t$ are computed using INFIDA and the expectation is over the randomized choices of DEPROUND. Note that, by taking the supremum over all request sequences and potential available capacities, we measure regret in an adversarial setting, i.e., against an adversary that selects, for every $t \in [T]$, vectors $\mathbf{r}_t$ and $\mathbf{l}_t$ to jeopardize the performance of our algorithm. Obviously, we do not expect such an adversary would exist in reality, but the adversarial analysis provides bounds on the behavior of INFIDA in the worst case.

The adversarial analysis is a modeling technique to characterize system performance under highly volatile external parameters (e.g., the sequence of requests $\mathbf{r}_t$) or difficult to model system interactions (e.g., the available capacities $\mathbf{l}_t$). This technique has been recently successfully used to model caching problems (e.g., in [70], [73]). Our main result is the following (the full proof is in Appendix G):

**Theorem V.1.** INFIDA has a sublinear $(1 - 1/e)$-regret w.r.t. the time horizon $T$, i.e., there exists a constant $A$ such that:

$$ (1 - 1/e)\text{-Regret}_{T, \mathcal{X}} \le A\sqrt{T}, \quad (21) $$

where $A \propto RL_{\max}\Delta_C$. $R$, $L_{\max}$, and $\Delta_C$ are upper bounds, respectively, on the total number of request types at any time slot, on the model capacities, and on the largest serving cost difference between serving at a repository node and at any other node.

*Proof.* (sketch) We first prove that the expected gain of the randomly sampled allocations $\mathbf{x}_t$ is a $(1-1/e)$-approximation of the fractional gain. Then, we use online learning results [69] to bound the regret of Online Mirror Ascent schemes operating on a convex decision space and against concave gain functions picked by an adversary. The two results are combined to obtain an upper bound on the $(1-1/e)$-regret. We fully characterize the regret constant $A$ in Appendix G. $\square$

Note that this result holds over the integral domain (see Appendix F, Lemmas F.7–F.11), thus generalizing the approximation techniques in [46]–[48] and providing a novel result in approximating budget-additive (submodular) set functions [49].

We observe that the regret bound depends crucially on the maximum number of request types $R$, maximum model capacity $L_{\max}$ and maximum serving cost difference $\Delta_C$. When considering the cost model in Eq. (6), we can consider for $\Delta_C$ the sum of the total latency of the heaviest path, the parameter $\alpha$, and the largest inference delay. This result is intuitive: when these values are bigger, the adversary has a larger room to select values that can harm the performance of the system.

As a direct consequence of Theorem V.1, the expected time averaged $(1 - 1/e)$-regret of INFIDA can get arbitrarily close to zero for large time horizon. Hence, INFIDA achieves a time averaged expected gain that is a $(1 - 1/e - \epsilon)$-approximation of the optimal time averaged static gain, for arbitrarily small $\epsilon$.

Observe that INFIDA computes a different $\mathbf{x}_t$ at every time slot. Intuitively, this allows it to “run after” the exogenous variation of the adversarial input $\{\mathbf{r}_t, \mathbf{l}_t\}_{t=1}^T \in \mathcal{A}^T$. An alternative goal that can be achieved by INFIDA is to find a static allocation $\bar{\mathbf{y}}$. In order to do so, we need to (i) run INFIDA for $\tilde{T}$ time-slots, (ii) based on the $\{\mathbf{y}_t\}_{t=1}^T$ computed by INFIDA, calculate $\bar{\mathbf{x}}$ (the exact calculation is in Proposition V.1.1), (iii) deploy in the IDN the allocation $\bar{\mathbf{x}}$ and keep it static, in order to avoid switches. Obviously, we would like the quality of $\bar{\mathbf{x}}$ to be close to the best $\mathbf{x}_*$, defined in (14). The following proposition shows that the gain achieved with our $\bar{\mathbf{x}}$ is boundedly close to the optimum. Moreover, since (14) is NP-hard, there cannot exist better bounds than the one we achieve, assuming $P \neq \text{NP}$ [67].

**Proposition V.1.1.** (offline solution) Replace in INFIDA the allocation gain $G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y})$ by $G_T(\mathbf{y})$ (defined in (14)). After $\tilde{T}$ iterations, let $\bar{\mathbf{y}}$ be the average fractional allocation $\bar{\mathbf{y}} = \frac{1}{\tilde{T}} \sum_{t=1}^T \mathbf{y}_t$, and $\bar{\mathbf{x}}$ the random state sampled from $\bar{\mathbf{y}}$ using DEPROUND. $\forall \epsilon > 0$, for $\tilde{T}$ large enough, $\bar{\mathbf{x}}$ satisfies

$$ \mathbb{E}[G_T(\bar{\mathbf{x}})] \ge \left(1 - \frac{1}{e} - \epsilon\right) G_T(\mathbf{x}_*), \quad (22) $$

where $\mathbf{x}_* = \arg\max_{\mathbf{x} \in \mathcal{X}} G_T(\mathbf{x})$.

The proof is given in Appendix H.
## VI. EXPERIMENTAL RESULTS

We evaluate INFIDA by simulating a realistic scenario based on the typical structure of ISP networks. We compare our solution with a greedy heuristic and its online variant (described below), as the greedy heuristic is known to achieve good performance in practice for submodular optimization [72].

**Topology.** We simulate a hierarchical topology similar to [74] that spans between edge and cloud, with different capacities at each tier. We consider 5 tiers: base stations (tier 4), central offices (tiers 3, 2), ISP data center (tier 1), a remote cloud (tier 0). We assume a hierarchical geographic distribution similar to LTE. We take the Round-Trip Time (RTT) across the different tiers as follows: tier 4 to tier 3 takes 6 ms, tier 3 to tier 2 takes 6 ms, tier 2 to tier 1 takes 15 ms, and tier 1 to tier 0 takes 40 ms. We execute our experiments at two different scales: *Network Topology I* counts 24 base stations and 36 nodes in total, while *Network Topology II* is a simpler 5-node scenario with 2 base stations.

**Processing Units.** We take GPU memory of the computing nodes as the limiting budget. The node at tier 0 can store the entire models catalog. We simulate the performance of two different processing units: the computing nodes at tiers 0 and 1 are equipped with high-end GPUs (Titan RTX), and the remaining tiers 2–4 have mid-tier GPUs (GeForce GTX 980). The budget of each computing tier is given as follows: a tier-1 node has 16GB GPUs, a tier-2 node has 12GB GPUs, a tier-3 node has 8GB GPUs, and a tier-4 node has 4GB GPUs.

**Catalog and requests.** We simulate performance based on state-of-the-art pre-trained models and their pruned versions [75], [76], profiled for each simulated processing unit, for a total of 10 models (Table II). We consider a task catalog with $|N| = 20$ different object detection tasks. We allow 3 duplicates per model; this gives $|M_i| = 30$ alternative models per task $i \in N$. Note how, as model complexity decreases, the number of frames a GPU can process per second increases, and consequently the average inference delay decreases.

The time slot duration is set to 1 minute and requests arrive at a constant rate of 7,500 requests per second (rps), unless otherwise said. Each request type is assigned randomly to two base stations in tier 4. The corresponding task is selected according to two different popularity profiles: (i) in the *Fixed Popularity Profile* (Fig. 4a), a request is for task *i* with constant probability $p(i) = \frac{(i+1)^{-1.2}}{\sum_{i' \in N} (i'+1)^{-1.2}}$ (a Zipf distribution with exponent 1.2), while (ii) in the *Sliding Popularity Profile* (Fig. 4b), the *l*-th consecutive request is for task *i* with probability $\tilde{p}(i, l) = p((i+5|l/W]) \mod 20)$, that is, the popularity of the tasks changes through a cyclic shift of 5 tasks every 1 hour for a request rate of 7,500 rps ($W = 2.7 \times 10^7$). **Static greedy.** We adapt the static greedy (SG) heuristic from the cost-benefit greedy in [72]. SG operates in hindsight seeking maximization of the time averaged allocation gain over the whole time horizon *T*, as in Eq. (14). Starting from an empty allocation, this policy progressively allocates the model that provides the highest marginal gain normalized by size, among those that meet the budget constraints. This process is repeated until either the intermediate allocation is capable of

<table><caption>TABLE II: Catalog for variants of YOLOv4 [75] profiled on two different Processing Units. Accuracy is for the MS COCO dataset. Values for the pruned variants are adapted from [76].</caption><thead><tr><th rowspan="2">variants of yolov4</th><th rowspan="2">accuracy (mAP@0.5)</th><th rowspan="2">memory (MB)</th><th colspan="2">frames per second</th></tr><tr><th>Titan RTX</th><th>GTX 980</th></tr></thead><tbody><tr><td>608p</td><td>65.7</td><td>1577</td><td>41.7</td><td>14.2</td></tr><tr><td>512p</td><td>64.9</td><td>1185</td><td>55.5</td><td>18.9</td></tr><tr><td>416p</td><td>62.8</td><td>1009</td><td>73.8</td><td>25.1</td></tr><tr><td>320p</td><td>57.3</td><td>805</td><td>100</td><td>34.1</td></tr><tr><td>3.99pruned</td><td>55.1</td><td>395</td><td>209</td><td>71.0</td></tr><tr><td>8.09pruned</td><td>51.4</td><td>195</td><td>329</td><td>112</td></tr><tr><td>10.10pruned</td><td>50.9</td><td>156</td><td>371</td><td>126</td></tr><tr><td>14.02pruned</td><td>49.0</td><td>112</td><td>488</td><td>166</td></tr><tr><td>tiny-416p</td><td>38.7</td><td>187</td><td>888</td><td>302</td></tr><tr><td>tiny-288p</td><td>34.4</td><td>160</td><td>1272</td><td>433</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 4: Popularity profiles of inference tasks for request rate 7,500 rps. Each dot represents a request for a given task at given time. Therefore, popular tasks correspond to denser lines. In (b), popularity changes at fixed time intervals through a cyclic shift. At any moment the requests are i.i.d. and sampled from a Zipf distribution with exponent 1.2. The figure shows a random down-sample of 5,000 requests to emphasize the density difference across the different tasks.</figcaption></figure>

serving all requests or none of the remaining valid allocations introduces a positive marginal gain.

**Online load-aware greedy heuristic.** As INFIDA is the first online policy for ML models’ allocation in IDNs, there is no clear baseline to compare it with. We then propose an online heuristic based on SG, which we call online load-aware greedy (OLAG). A node *v* uses counters $\phi_{m,\rho}^v$ to keep track of the number of times a request $\rho \in \mathcal{R}$ is forwarded upstream but could have been served locally at a lower cost compared to the repository, i.e., using a model $m \in M$ with positive gain that we denote by $q_{m,\rho}^v$. For every model *m*, an importance weight is computed as $w_m^v = \frac{1}{s_m^v |\mathcal{R}|} \sum_{\rho \in \mathcal{R}} q_{m,\rho}^v \min\{\phi_{m,\rho}^v, L_m^v\}$, where $s_m^v$ is the size of model *m* and $\min\{\phi_{m,\rho}^v, L_m^v\}$ is the number of requests that could have been improved by *m*. At the end of a time slot, the node selects the model $m_*$ with the highest importance while respecting the resource budget constraint, then subtracts the quantity $\min\{\phi_{m,\rho}^v, L_m^v\}$ from $\phi_{m_*,\rho}^v$ and from all the $\phi_{m',\rho}^v: q_{m',\rho}^v < q_{m_*,\rho}$, i.e., models that provide a gain lower than $m_*$. This procedure is repeated until the resource budget of the node is consumed.
**Offline INFIDA.** Motivated by Proposition V.1.1, we implemented also an offline version of INFIDA that we call INFIDA<sub>OFFLINE</sub>, which maximizes the time-averaged gain (14) over the whole time horizon *T*. The potential available capacities are determined at runtime from the current allocations and request batches (rather than by an adversary).

**Performance Metrics.** The performance of a policy *P* with the associated sequence of allocation decisions {*x*<sub>*t*</sub>}<sub>*t*<sup>*T*</sup></sub> is evaluated in terms of the time-averaged gain normalized to the number of requests per time slot (NTAG):

$$ \text{NTAG}(\mathcal{P}) = \sum_{t=1}^{T} \frac{1}{T \|\mathbf{r}_t\|_1} G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}_t). \quad (23) $$

Moreover, we evaluate the update cost of a policy *P* with the associated sequence of allocation decisions {*x*<sub>*t*</sub>}<sub>*t*<sup>*T*</sup></sub> by quantifying the total size of fetched models over *T* time slots. The update cost is reflected by the Time-Averaged Model Updates (MU) metric defined as:

$$ \text{MU}(\mathcal{P}) \triangleq \frac{1}{T} \sum_{t=2}^{T} \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}} s_m^v \max\{0, x_{t,m}^v - x_{t-1,m}^v\}. \quad (24) $$

### A. Trade-off between Latency and Accuracy

We first evaluate how INFIDA adapts to different trade-offs between end-to-end latency and inference accuracy by varying the trade-off parameter $\alpha$.<fnref n="8" />

Figure 5 shows the fractional allocation decision at each tier of the network topology for different values of $\alpha$ (remember that the smaller $\alpha$ the more importance is given to the latency rather than to inaccuracy, see Eq. (6)). Models are ordered horizontally by increasing accuracy with 3 potential replicas for each model, and only the models able to serve the most popular request are shown. Note that the tier-0 node acts as a repository and its allocation is fixed; moreover, in Fig 5 the repository node picks the second most accurate model because it provides the smallest combined cost in Eq. (6).

For $\alpha = 3$ (Fig. 5a), INFIDA allocates a considerable amount of small models (which provide low accuracy) near the edge (tiers 1–3 and model IDs 0–18), as they can serve a larger number of requests compared to higher quality models with low inference delay. By giving more importance to the accuracy ($\alpha = 4$) the system tends to deploy more accurate models and rarely allocates small models (Fig. 5b). For $\alpha = 5$, the number of models deployed on lower tiers decreases, as the system allocates no small models in practice (model IDs 0–20) and selects instead multiple replicas of the most accurate models (Fig. 5c). Since higher quality models feature, in general, a lower serving capacity (Table II), Fig. 5c suggests that a significant number of requests is served in the cloud (Tier 0) for this value of $\alpha$.

Figure 6 shows the average experienced inaccuracy (inaccuracy is given by 100 – mAP and mAP is the mean average precision) and latency for different values of $\alpha$ under Network

<figure><img src="image_5.png" /><figcaption>Fig. 5: Fractional allocation decisions $y_m^v$ of INFIDA on the various tiers of Network Topology I under Fixed Popularity Profile. We only show the allocations corresponding to the models capable to serve the most popular request. The model IDs are sorted by increasing accuracy.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 6: Average latency (dashed line) and inaccuracy (solid line) costs experienced with INFIDA for different values of $\alpha$ under Network Topology I and Fixed Popularity Profile.</figcaption></figure>

Topology I and Fixed Popularity Profile. When accuracy is not important (i.e., $\alpha \approx 0$), INFIDA effectively achieves very low end-to-end latency (few milliseconds) by prioritizing the deployment of small and inaccurate models near to the edge nodes. Noticeably, the trend in both curves (decreasing inaccuracy and increasing latency) suggests that, when higher accuracy is required, the system starts to prefer models deployed close to the cloud, leading to a sudden change in the trade-off and to a significant increase in latency.

In Fig. 7 we show the normalized time-averaged gain of INFIDA compared to OLAG, SG, and INFIDA<sub>OFFLINE</sub> for different values of $\alpha$ under the Sliding Popularity Profile. Results are shown both for Network Topology I (Fig. 7a) and for Network Topology II (Fig. 7b).

The plot shows that the gain decreases by increasing $\alpha$. This is expected since the gain (13) is defined as the improvement

<footnote n="8">⁸The inaccuracy cost is taken in 0–100, then $\alpha$ picked here corresponds to 100× scaling of the parameter defined in Eq. (6).</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 7: NTAG of the different policies under Sliding Popularity Profile and network topologies: (a) Network Topology I, and (b) Network Topology II.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 8: (a) Models Updates (MU) and (b) NTAG of OLAG and INFIDA for different values of refresh period $B \in \{4, 8, 16\}$, and for a dynamic refresh period with initial value $B_{\text{init}} = 1$, target value $B_{\text{target}} = 32$ and stretching duration $\Delta t = 60$ (1H). The experiment is run under Network Topology I and Sliding Popularity Profile.</figcaption></figure>

w.r.t. the repository allocation (tier 0). Therefore, when the latency is not important, high accuracy models at tier 0 are preferred, and there is no much room for improvement (the optimal gain eventually tends to zero for $\alpha \to +\infty$). Note that, in general, SG and INFIDA<sub>OFFLINE</sub> policies perform worse than their offline counterparts, as they pick a single allocation that is the best w.r.t. the whole sequence of requests. However, in the *Sliding Popularity Profile* (Fig. 4) the best decision changes periodically, and only the online policies have the ability to adapt to such change. Moreover, we observe that consistently INFIDA<sub>OFFLINE</sub> has better performance than SG: although both policies are offline, INFIDA<sub>OFFLINE</sub> manages to provide a better allocation.

## B. Trade-off between model updates and service cost.

In this set of experiments, we evaluate how the frequency at which INFIDA updates the model allocation affects the update cost incurred by the system. Indeed, frequent updates

<figure><img src="image_9.png" /><figcaption>Fig. 9: NTAG of the different policies for different request rates under Network Topology I.</figcaption></figure>

could lead to massive migrations with an overhead on network bandwidth. As an evaluation metric, we measure the total size of fetched models averaged over time (see the performance metric in Eq. (24)). We introduce *B* that we call the *refresh period*, and we restrict INFIDA to only sample a physical allocation every $B \in \{4, 8, 32\}$ time slots (line 8 in Algorithm 1). Additionally, we experiment linear stretching of the refresh period *B* with initial period $B_{\text{init}} = 1$ and target period $B_{\text{target}} = 32$ in a stretching duration of $\Delta t = 1$H. We run this experiment under *Network Topology I* and *Sliding Popularity Profile*. We set the trade-off parameter $\alpha = 1$.

In particular, Figure 8a shows the update cost (MU) for different refresh periods, while Figure 8b shows the NTAG. Both plots include the performance of the OLAG heuristic. We observe that, by increasing the refresh period *B*, the system fetches a smaller number of models, and therefore the update cost decreases, at the expense of reactivity. This tradeoff was characterized formally in [66], wherein the regret is sublinear for $B = \Theta(T^\beta)$ for $\beta \in [0, 1)$, and the update costs are sublinear for $\beta \in (0, 1)$. Nevertheless, even for large values of *B* INFIDA eventually exceeds OLAG in performance: this result is expected since the algorithm continues to learn on the fractional (virtual) states and only the physical allocations are delayed and eventually catch-up for a large time horizon. On the other hand, we observe that OLAG is relatively conservative in updating its allocation, as it quickly picks a sub-optimal allocation and rarely updates it.

The previous observation motivates the use of a dynamic refresh period. By refreshing more frequently at the start we allow the physical allocation generated by INFIDA to catch up quickly with the fractional states as shown in Fig. 8b: a dynamic refresh period that stretches from $B_{\text{init}} = 1$ to $B_{\text{target}} = 32$ attains much faster and more precise convergence. This is achieved at the expense of a high update cost at the start, which is, however, quickly dampened until it matches the same update cost of fixing the refresh period to $B_{\text{target}}$.
<figure><img src="image_10.png" /><figcaption>Fig. 10: Average Latency vs. Average Inaccuracy obtained for different values of $\alpha \in \{0.5, 1, 2, 3, 4, 5, 6\}$ under Fixed Popularity Profile and Network Topology II.</figcaption></figure>

### C. Scalability on Requests Load

We show how the system performs under different requests loads. For this set of experiments, we set $\alpha = 1$. Figure 9 compares the results for the different allocation policies.

We notice that, being INFIDA<sub>OFFLINE</sub> and SG offline policies, they perform well when the popularity profile is static (Fig. 9a), but deteriorate under *Sliding Popularity Profile*. Notably, the performance degradation of INFIDA<sub>OFFLINE</sub> (≈8%) is considerably limited compared to SG (≈30%), which even gets worse when increasing the requests load.

Figure 9 shows that, in general, INFIDA provides a higher gain compared to the OLAG heuristic. In particular, under *Fixed Popularity Profile* INFIDA manages to converge to the same NTAG provided by its offline counterpart (Fig. 9a), which is ≈10% better than the one provided by OLAG when the load is 7,083 rps. Additionally, OLAG's performance deteriorates when the requests load increases from 7,083 rps to 10,000 rps. It is also noteworthy that OLAG visibly suffers from perturbed performance when the popularity of the tasks changes over time (Fig. 9b). On the other hand, results show the robustness of INFIDA against changing request loads and popularity: the algorithm preserves its performance in terms of normalized time-averaged gain for the analyzed request loads and under both *Fixed Popularity Profile* and *Sliding Popularity Profile*, always converging to the highest NTAG.

Last, in Fig. 10 we evaluate separately the average latency and inaccuracy attained by the different policies using different values of $\alpha \in \{0.5, 1, 2, 3, 4, 5, 6\}$ under *Fixed Popularity Profile* and *Network Topology II*. We observe that INFIDA and its offline counterpart INFIDA<sub>OFFLINE</sub> consistently provide the lowest average inaccuracy and latency under both high request load (10,000 rps) and default request load (7,500 rps). INFIDA<sub>OFFLINE</sub> is run with hindsight and serves as a lower bound on the achievable latency and inaccuracy under a fixed popularity request process (as in Fig. 9a).

## VII. CONCLUSIONS

In this paper, we introduced the idea of inference delivery networks (IDNs), networks of computing nodes that coordinate to satisfy inference requests in the continuum between Edge and Cloud. IDN nodes can serve inference requests with different levels of accuracy and end-to-end latency, based on their geographic location and processing capabilities. We formalized the NP-hard problem of allocating ML models on IDN nodes, capturing the trade-off between latency and accuracy. We proposed INFIDA, a dynamic ML model allocation algorithm that operates in a distributed fashion and provides strong guarantees in an adversarial setting. We evaluated INFIDA simulating the realistic scenario of an ISP network, and compared its performance under two different topologies with both an offline greedy heuristic and its online variant. Our results show that INFIDA adapts to different latency/accuracy trade-offs and scales well with the number of requests, outperforming the greedy policies in all the analyzed settings.

## VIII. ACKNOWLEDGEMENT

This work has been carried out in the framework of a common lab agreement between Inria and Nokia Bell Labs. This research was supported in part by the French Government through the "Plan de Relance" and "Programme d'investissements d'avenir" and by Inria under the exploratory action MAMMALS.

## REFERENCES

[1] I. Stoica, D. Song, R. A. Popa et al., "A Berkeley View of Systems Challenges for AI," EECS Department, University of California, Berkeley, Tech. Rep. UCB/EECS-2017-159, Oct 2017. [Online]. Available: https://www2.eecs.berkeley.edu/Pubs/TechRpts/2017/EECS-2017-159.html

[2] M. Simsek, A. Aijaz, M. Dohler et al., "5G-Enabled Tactile Internet," *IEEE Journal on Selected Areas in Communications*, vol. 34, no. 3, pp. 460-473, 2016.

[3] A. G. Howard, M. Zhu, B. Chen et al., "MobileNets: Efficient Convolutional Neural Networks for Mobile Vision Applications," preprint arXiv:1704.04861, 2017.

[4] L. Deng, G. Li, S. Han et al., "Model Compression and Hardware Acceleration for Neural Networks: A Comprehensive Survey," *Proceedings of the IEEE*, 2020.

[5] E. Hazan et al., "Introduction to Online Convex Optimization," *Foundations and Trends in Optimization*, vol. 2, no. 3-4, pp. 157-325, 2016.

[6] S. Babu and H. Herodotou, "Massively Parallel Databases and MapReduce Systems," *Foundations and Trends in Databases*, 2013.

[7] M. Zaharia, M. Chowdhury, T. Das et al., "Resilient Distributed Datasets: A Fault-Tolerant Abstraction for In-Memory Cluster Computing," in *9th USENIX Symposium on Networked Systems Design and Implementation (NSDI 12)*, 2012, pp. 15-28.

[8] P. Kairouz, H. B. McMahan, B. Avent et al., "Advances and Open Problems in Federated Learning," arXiv preprint arXiv:1912.04977, 2021.

[9] T. Li, A. K. Sahu, A. Talwalkar et al., "Federated Learning: Challenges, Methods, and Future Directions," *IEEE Signal Processing Magazine*, vol. 37, no. 3, p. 50-60, May 2020.

[10] J. Konečný, H. B. McMahan, D. Ramage et al., "Federated Optimization: Distributed Machine Learning for On-Device Intelligence," arXiv preprint arXiv:1610.02527, 2016.

[11] J. Konečný, H. B. McMahan, F. X. Yu et al., "Federated Learning: Strategies for Improving Communication Efficiency," arXiv preprint arXiv:1610.05492, 2016.

[12] W. Wu, L. He, W. Lin et al., "Accelerating Federated Learning over Reliability-Agnostic Clients in Mobile Edge Computing Systems," *IEEE Transactions on Parallel and Distributed Systems*, 2020.
[13] G. Neglia, G. Calbi, D. Towsley et al., "The Role of Network Topology for Distributed Machine Learning," in *IEEE Conference on Computer Communications (INFOCOM)*. IEEE, 2019, pp. 2350-2358.

[14] C. Olston, N. Fiedel, K. Gorovoy et al., "TensorFlow-Serving: Flexible, High-Performance ML Serving," preprint arXiv:1712.06139, 2017.

[15] D. Chappell, "Introducing Azure Machine Learning," A guide for technical professionals, sponsored by microsoft corporation, 2015.

[16] "Vertex AI | google cloud."

[17] D. Crankshaw, X. Wang, G. Zhou et al., "Clipper: A Low-Latency Online Prediction Serving System," in *14th USENIX Symposium on Networked Systems Design and Implementation (NSDI 17)*, 2017, pp. 613-627.

[18] H. Mao, M. Schwarzkopf, S. B. Venkatakrishnan et al., "Learning Scheduling Algorithms for Data Processing Clusters," in *Proceedings of the ACM Special Interest Group on Data Communication*, 2019, pp. 270-288.

[19] F. Romero, Q. Li, N. J. Yadwadkar et al., "INFaaS: Managed and Model- less Inference Serving," preprint arXiv:1905.13348, 2019.

[20] D. Crankshaw, G.-E. Sela, X. Mo et al., "InferLine: Latency-aware Pro- visioning and Scaling for Prediction Serving Pipelines," in *Proceedings of the 11th ACM Symposium on Cloud Computing*, 2020, pp. 477-491.

[21] X. Meng, J. Bradley, B. Yavuz et al., "MLlib: Machine Learning in Apache Spark," *The Journal of Machine Learning Research*, vol. 17, no. 1, pp. 1235-1241, 2016.

[22] Y. Jia, E. Shelhamer, J. Donahue et al., "Caffe: Convolutional Archi- tecture for Fast Feature Embedding," in *Proceedings of the 22nd ACM international conference on Multimedia*, 2014, pp. 675-678.

[23] F. Pedregosa, G. Varoquaux, A. Gramfort et al., "Scikit-learn: Machine Learning in Python," *the Journal of machine Learning research*, vol. 12, pp. 2825-2830, 2011.

[24] N. D. Lane, S. Bhattacharya, P. Georgiev et al., "DeepX: A Software Ac- celerator for Low-Power Deep Learning Inference on Mobile Devices," in *2016 15th ACM/IEEE International Conference on Information Pro- cessing in Sensor Networks (IPSN)*, 2016.

[25] N. Fernando, S. W. Loke, and W. Rahayu, "Computing with Nearby Mobile Devices: A Work Sharing Algorithm for Mobile Edge-Clouds," *IEEE Transactions on Cloud Computing*, vol. 7, no. 2, pp. 329-343, 2019.

[26] Y. Kang, J. Hauswald, C. Gao et al., "Neurosurgeon: Collaborative Intelligence Between the Cloud and Mobile Edge," *ACM SIGARCH Computer Architecture News*, 2017.

[27] S. Teerapittayanon, B. McDanel, and H. T. Kung, "BranchyNet: Fast Inference via Early Exiting from Deep Neural Networks," in *2016 23rd International Conference on Pattern Recognition (ICPR)*, 2016.

[28] S. Teerapittayanon, B. McDanel, and H.-T. Kung, "Distributed Deep Neural Networks Over the Cloud, the Edge and End Devices," in *2017 IEEE 37th International Conference on Distributed Computing Systems (ICDCS)*, 2017.

[29] S. Deng, H. Zhao, J. Yin et al., "Edge Intelligence: The Confluence of Edge Computing and Artificial Intelligence," *IEEE Internet of Things Journal*, 2020.

[30] S. S. Ogden and T. Guo, "MODI: Mobile Deep Inference Made Efficient by Edge Computing," in *USENIX Workshop on Hot Topics in Edge Computing (HotEdge 18)*, 2018.

[31] Y. Jin, L. Jiao, Z. Qian et al., "Provisioning Edge Inference as a Service via Online Learning," in *2020 17th Annual IEEE International Conference on Sensing, Communication, and Networking (SECON)*, 2020.

[32] C.-C. Hung, G. Ananthanarayanan, P. Bodik et al., "VideoEdge: Pro- cessing Camera Streams using Hierarchical Clusters," in *IEEE/ACM Symposium on Edge Computing*, 2018.

[33] X. Qiu, H. Li, C. Wu et al., "Cost-Minimizing Dynamic Migration of Content Distribution Services into Hybrid Clouds," *IEEE Transactions on Parallel and Distributed Systems*, 2015.

[34] J. Xu, L. Chen, and P. Zhou, "Joint Service Caching and Task Offloading for Mobile Edge Computing in Dense Networks," in *IEEE INFOCOM 2018 - IEEE Conference on Computer Communications*, 2018, pp. 207-215.

[35] S. Wang, R. Urgaonkar, T. He et al., "Dynamic Service Placement for Mobile Micro-Clouds with Predicted Future Costs," *IEEE Transactions on Parallel and Distributed Systems*, 2017.

[36] A. Ben Ameur, A. Araldo et al., "On the Deployability of Augmented Reality Using Embedded Edge Devices," in *2021 IEEE 18th Annual Consumer Communications & Networking Conference (CCNC)*, 2021.

[37] M. Choi, J. Kim, and J. Moon, "Wireless Video Caching and Dynamic Streming Under Differentiated Quality Requirements," *IEEE Journal on Selected Areas in Communications*, vol. 36, 2018.

[38] Z. Ye, F. De Pellegrini, R. El-Azouzi et al., "Quality-Aware DASH Video Caching Schemes at Mobile Edge," in *2017 29th International Teletraffic Congress (ITC 29)*, vol. 1. IEEE, 2017, pp. 205-213.

[39] C. Zhan and Z. Wen, "Content Cache Placement for Scalable Video in Heterogeneous Wireless Network," *IEEE Communications Letters*, vol. 21, no. 12, pp. 2714-2717, 2017.

[40] A. Araldo, F. Martignon, and D. Rossi, "Representation Selection Problem: Optimizing Video Delivery through Caching," in *2016 IFIP Networking Conference (IFIP Networking) and Workshops*. IEEE, 2016, pp. 323-331.

[41] Z. Qu, B. Ye, B. Tang et al., "Cooperative Caching for Multiple Bitrate Videos in Small Cell Edges," *IEEE Transactions on Mobile Computing*, vol. 19, no. 2, pp. 288-299, 2020.

[42] K. Poularakis, G. Iosifidis, A. Argyriou et al., "Video delivery over heterogeneous cellular networks: Optimizing cost and performance," in *IEEE INFOCOM 2014-IEEE Conference on Computer Communications*. IEEE, 2014, pp. 1078-1086.

[43] ——, "Caching and Operator Cooperation Policies for Layered Video Content Delivery," in *IEEE INFOCOM 2016-The 35th Annual IEEE International Conference on Computer Communications*. IEEE, 2016, pp. 1-9.

[44] S. Shalev-Shwartz, "Online Learning and Online Convex Optimization," *Found. Trends Mach. Learn.*, vol. 4, no. 2, p. 107-194, Feb. 2012.

[45] H. B. McMahan, "A survey of algorithms and analysis for adaptive online learning," *The Journal of Machine Learning Research*, vol. 18, no. 1, pp. 3117-3166, 2017.

[46] S. Ioannidis and E. Yeh, "Adaptive Caching nNetworks with Optimal- ity Guarantees," *ACM SIGMETRICS Performance Evaluation Review*, vol. 44, no. 1, pp. 113-124, 2016.

[47] D. Paria and A. Sinha, "LeadCache: Regret-Optimal Caching in Networks," in *Advances in Neural Information Processing Systems*, M. Ranzato, A. Beygelzimer, Y. Dauphin et al., Eds., vol. 34. Curran Associates, Inc., 2021, pp. 4435-4447.

[48] K. Shanmugam, N. Golrezaei, A. G. Dimakis et al., "FemtoCaching: Wireless Content Delivery Through Distributed Caching Helpers," *IEEE Transactions on Information Theory*, vol. 59, no. 12, pp. 8402-8413, 2013.

[49] J. Garg, M. Hoefer, and K. Mehlhorn, *Approximating the Nash Social Welfare with Budget-Additive Valuations*, pp. 2326-2340.

[50] A. Mehta, "Online Matching and Ad Allocation," *Foundations and Trends® in Theoretical Computer Science*, vol. 8, no. 4, pp. 265-368, 2013.

[51] A. Mehta, A. Saberi, U. Vazirani et al., "AdWords and Generalized Online Matching," *J. ACM*, vol. 54, no. 5, p. 22-es, oct 2007.

[52] M. Feldman, N. Gravin, and B. Lucier, "Combinatorial Walrasian Equilibrium," *SIAM Journal on Computing*, vol. 45, no. 1, pp. 29-48, 2016.

[53] T. Roughgarden and I. Talgam-Cohen, "Why Prices Need Algorithms," in *Proceedings of the Sixteenth ACM Conference on Economics and Computation*, ser. EC '15. New York, NY, USA: Association for Computing Machinery, 2015, p. 19-36.

[54] M. Garetto, E. Leonardi, and G. Neglia, "Similarity Caching: Theory and Algorithms," in *IEEE INFOCOM 2020-IEEE Conference on Computer Communications*, 2020.

[55] F. Falchi, C. Lucchese, S. Orlando et al., "A Metric Cache for Similarity Search," in *Proceedings of the 2008 ACM workshop on Large-Scale distributed systems for information retrieval*, 2008, pp. 43-50.

[56] S. Pandey, A. Broder, F. Chierichetti et al., "Nearest-Neighbor Caching for Content-Match Applications," in *Proceedings of the 18th international conference on World wide web*, 2009, pp. 441-450.

[57] U. Drolia, K. Guo, J. Tan et al., "Cachier: Edge-Caching for Recognition Applications," in *2017 IEEE 37th International Conference on Distributed Computing Systems (ICDCS)*. IEEE, 2017, pp. 276-286.

[58] P. Sermpezis, T. Giannakas, T. Spyropoulos et al., "Soft Cache Hits: Improving Performance Through Recommendation and Delivery of Related Content," *IEEE Journal on Selected Areas in Communications*, 2018.

[59] J. Zhou, O. Simeone, X. Zhang et al., "Adaptive Offline and Online Similarity-Based Caching," *IEEE Networking Letters*, vol. 2, no. 4, pp. 175-179, 2020.

[60] M. Garetto, E. Leonardi, and G. Neglia, "Content Placement in Networks of Similarity Caches," *Computer Networks*, vol. 201, p. 108570, 2021.

[61] D. Blalock, J. J. G. Ortiz, J. Frankle et al., "What is the State of Neural Network Pruning?" *arXiv preprint arXiv:2003.03033*, 2020.

[62] G. Hinton et al., "Distilling the Knowledge in a Neural Network," preprint arXiv:1503.02531, 2015.

[63] S. Ravi, "Custom On-Device ML Models with Learn2Compress," 2018.
[64] Z. Fang, T. Yu, O. J. Mengshoel et al., "QoS-Aware Scheduling of Heterogeneous Servers for Inference in Deep Neural Networks," in *ACM Conference on Information and Knowledge Management (CIKM)*, vol. Part F1318, 2017, pp. 2067–2070.

[65] Y. Wang, G. Y. Wei, and D. Brooks, "Benchmarking TPU, GPU, and CPU platforms for deep learning," *arXiv preprint arXiv:1907.10701*, 2019.

[66] T. Si Salem, G. Neglia, and D. Carra, "Ascent Similarity Caching With Approximate Indexes," *IEEE/ACM Transactions on Networking*, vol. 31, no. 3, pp. 1173–1186, 2023.

[67] G. L. Nemhauser and L. A. Wolsey, "Best Algorithms for Approximating the Maximum of a Submodular Set Function," *Mathematics of operations research*, 1978.

[68] Y. Fairstein, A. Kulik, J. S. Naor et al., "A $(1-1/e-\epsilon)$-Approximation for the Monotone Submodular Multiple Knapsack Problem," in *28th Annual European Symposium on Algorithms (ESA 2020)*. Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2020.

[69] S. Bubeck, "Convex Optimization: Algorithms and Complexity," *Foundations and Trends in Machine Learning*, vol. 8, no. 3-4, pp. 231–357, Nov. 2015.

[70] T. Si Salem, G. Neglia, and S. Ioannidis, "No-Regret Caching via Online Mirror Descent," *ACM Trans. Model. Perform. Eval. Comput. Syst.*, vol. 8, no. 4, aug 2023.

[71] J. Byrka, T. Pensyl, B. Rybicki et al., "An Improved Approximation for $k$-median, and Positive Correlation in Budgeted Optimization - Extended Arxiv version," in *ACM-SIAM symposium on Discrete algorithms*, 2014.

[72] A. Krause and D. Golovin, "Submodular Function Maximization," *Tractability*, vol. 3, pp. 71–104, 2014.

[73] G. S. Paschos, A. Destounis, L. Vigneri et al., "Learning to Cache With No Regrets," in *IEEE INFOCOM*, 2019.

[74] A. Ceselli, M. Premoli, and S. Secci, "Mobile Edge Cloud Network Design Optimization," *IEEE/ACM Transactions on Networking*, vol. 25, no. 3, pp. 1818–1831, 2017.

[75] A. Bochkovskiy, C.-Y. Wang, and H.-Y. M. Liao, "YOLOv4: Optimal Speed and Accuracy of Object Detection," *preprint arXiv:2004.10934*, 2020.

[76] Y. Cai et al., "YOLObile: Real-Time Object Detection on Mobile Devices via Compression-Compilation Co-Design," *preprint arXiv:2009.05697*, 2020.

[77] G. Neglia, M. Garetto, and E. Leonardi, "Similarity Caching: Theory and Algorithms," *IEEE/ACM Trans. Netw.*, vol. 30, no. 2, p. 475–486, dec 2021.

[78] B. S. Mordukhovich and N. M. Nam, "Geometric Approach to Convex Subdifferential Calculus," *Optimization*, vol. 66, no. 6, pp. 839–873, 2017.

[79] R. T. Rockafellar, *Convex Analysis*. Princeton University Press, 2015.

[80] S. Shalev-Shwartz, "Online Learning: Theory, Algorithms, and Applications," Ph.D. dissertation, The Hebrew University of Jerusalem, 2007.

[81] M. X. Goemans and D. P. Williamson, "New $\frac{3}{4}$-Approximation Algorithms for the Maximum Satisfiability Problem," *SIAM Journal on Discrete Mathematics*, vol. 7, no. 4, pp. 656–666, 1994.

[82] H. S. Wilf, "Some Applications of the Inequality of Arithmetic and Geometric Means to Polynomial Equations," in *Proceedings of the American Mathematical Society*, vol. 14, no. 2. JSTOR, 1963, pp. 263–265.
Supplementary Material for Paper:

Towards Inference Delivery Networks:
Distributing Machine Learning with Optimality Guarantees

APPENDIX A

SUBMODULARITY OF THE GAIN FUNCTION

As submodularity is defined for set functions, let us associate to the gain function $G(\cdot)$ an opportune set function as follows. Given a set $S \subseteq \mathcal{V} \times \mathcal{M}$ of pairs $(v, m) \in \mathcal{V} \times \mathcal{M}$ (nodes and models), we define the corresponding associated vector $\mathbf{x}(S)$ with $x_m^v(S) = 1$ if $((v, m) \in S \lor \omega_m^v = 1)$, $x_m^v(S) = 0$ otherwise. Now we can define the set function $f_t: 2^{\mathcal{V} \times \mathcal{M}} \rightarrow \mathbb{R}$:

$$f_t(S) \triangleq G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}(S)). \qquad (25)$$

We can also define the set function $F_T: 2^{\mathcal{V} \times \mathcal{V}} \rightarrow \mathbb{R}$ associated to the time-averaged gain in Eq. (14) as

$$F_T(S) \triangleq \frac{1}{T} \sum_{t=1}^{T} f_t(S) = \frac{1}{T} \sum_{t=1}^{T} G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}(S)). \qquad (26)$$

We first start proving that $f_t$ is submodular in the following lemma.

**Lemma A.1.** The set function $f_t: 2^{\mathcal{V} \times \mathcal{M}} \rightarrow \mathbb{R}$ in Eq. (25) is normalized (i.e., $f_t(\emptyset) = 0$), submodular, and monotone.

*Proof.*

**Normalization.** The constructed set function is normalized (as in $f_t(\emptyset) = 0$), we have

$$f_t(\emptyset) = G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}(\emptyset)) = G(\mathbf{r}_t, \mathbf{l}_t, \omega) = C(\mathbf{r}_t, \mathbf{l}_t, \omega) - C(\mathbf{r}_t, \mathbf{l}_t, \omega) = 0. \qquad (27)$$

**Submodularity.** A function $f_t: 2^{\mathcal{V} \times \mathcal{M}} \rightarrow \mathbb{R}$ is submodular [72] if for every $S' \subset S'' \subset \mathcal{V} \times \mathcal{M}$ and $(\bar{v}, \bar{m}) \in (\mathcal{V} \times \mathcal{M}) \setminus S''$ it holds that

$$f_t(S'' \cup \{(\bar{v}, \bar{m})\}) - f_t(S'') \le f_t(S' \cup \{(\bar{v}, \bar{m})\}) - f_t(S'').$$

Let us consider $(\bar{v}, \bar{m}) \in (\mathcal{V} \times \mathcal{M}) \setminus S''$. We take $\mathbf{x}', \mathbf{x}'', \bar{\mathbf{x}}',$ and $\bar{\mathbf{x}}''$ as short hand notation for $\mathbf{x}(S')$, $\mathbf{x}(S'')$, $\mathbf{x}(S' \cup \{(\bar{v}, \bar{m})\})$, and $\mathbf{x}(S'' \cup \{(\bar{v}, \bar{m})\})$, respectively.

Since $S' \subset S''$, we have that if $x_m^v = 1 \Rightarrow x_m'' = 1$, and thus $z_\rho^k(l_t, \mathbf{x}') \le z_\rho^k(l_t, \mathbf{x}'')$, for any $k$ (see Eq. (11)). Therefore

$$Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}') \leq Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}''), \forall k \in [K_{\rho} - 1], \forall \rho \in \mathcal{R}. \qquad (28)$$

Let us denote $\bar{k}_\rho \triangleq \kappa_\rho(\bar{v}, \bar{m})$. Due to Eq. (28), $\forall k \in [K_\rho - 1], \forall \rho \in \mathcal{R}$ the following inequality holds:

$$\min\{r_{\rho}^{t} - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}''), \lambda_{\rho}^{\bar{k}_{\rho}}(\mathbf{l}_{t})\} \le \min\{r_{\rho}^{t} - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}'), \lambda_{\rho}^{\bar{k}_{\rho}}(\mathbf{l}_{t})\}, \qquad (29)$$

or equivalently:

$$\min\{r_{\rho}^{t}, Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}^{''}) + \lambda_{\rho}^{\bar{k}_{\rho}}(\mathbf{l}_{t})\} - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}^{''}) \le \min\{r_{\rho}^{t}, Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}^{''}) + \lambda_{\rho}^{\bar{k}_{\rho}}(\mathbf{l}_{t})\} - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}^{''}). \qquad (30)$$

Observe that

$$z_{\rho}^{k}(\mathbf{l}_{t}, \bar{\mathbf{x}}'') = \begin{cases} z_{\rho}^{k}(\mathbf{l}_{t}, \mathbf{x}'') & \text{if } k \neq \bar{k}_{\rho}, \\ z_{\rho}^{k}(\mathbf{l}_{t}, \mathbf{x}'') + \underbrace{x_{\bar{m}}^{\bar{v}}}_{=1} \cdot l_{\rho, \bar{m}}^{t, \bar{v}} \stackrel{(11)}{=} z_{\rho}^{k}(\mathbf{l}_{t}, \mathbf{x}'') + \lambda_{\rho}^{\bar{k}_{\rho}}(\mathbf{l}_{t}) & \text{if } k = \bar{k}_{\rho}, \end{cases}$$

and thus

$$Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \bar{\mathbf{x}}'') = \begin{cases} Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}'') & \text{if } k < \bar{k}_{\rho}, \\ \min\{r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_{t}, \mathbf{x}'') + \lambda_{\rho}^{\bar{k}_{\rho}}(\mathbf{l}_{t})\} & \text{if } k \ge \bar{k}_{\rho}, \end{cases} \qquad (31)$$

Note that the same equality holds between $\bar{\mathbf{x}}'$ and $\mathbf{x}'$ since $(\bar{v}, \bar{m}) \notin S'$.

The marginal gain of adding pair $(\bar{v}, \bar{m}) \in \mathcal{V} \times \mathcal{M}$ to the allocation set $S''$ is

$$f_t(S'' \cup \{(\bar{v}, \bar{m})\}) - f_t(S'') = G(\mathbf{r}_t, \mathbf{l}_t, \bar{\mathbf{x}}'') - G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}'')$$

$$\stackrel{(16)}{=} \sum_{\rho \in \mathcal{R}} \left[ \sum_{k=1}^{\bar{k}_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) (Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}'') - Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \omega)) \right]$$
$$
\begin{align}
&+ \sum_{k=\bar{k}_{\rho}}^{K_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) \left( \min \left\{ r_{\rho}^{t}, Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}^{\prime\prime}) + \lambda_{\rho}^{\bar{k}_{\rho}}(\mathbf{l}_{t}) \right\} - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{\omega}) \right) \nonumber \\
&- \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) \left( Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}^{\prime\prime}) - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{\omega}) \right) \nonumber \\
&= \sum_{\rho \in \mathbb{R}} \sum_{k=\bar{k}_{\rho}}^{K_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) \left[ \left( \min \left\{ r_{\rho}^{t}, Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}^{\prime\prime}) + \lambda_{\rho}^{\bar{k}_{\rho}}(\mathbf{l}_{t}) \right\} - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{\omega}) \right) - (Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}^{\prime\prime}) - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{\omega})) \right] \nonumber \\
&= \sum_{\rho \in \mathbb{R}} \sum_{k=\bar{k}_{\rho}}^{K_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) \left( \min \left\{ r_{\rho}^{t}, Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}^{\prime\prime}) + \lambda_{\rho}^{\bar{k}_{\rho}}(\mathbf{l}_{t}) \right\} - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}^{\prime\prime}) \right). \tag{32}
\end{align}
$$

By bounding up each term of the marginal gain (32) as in (30) we get the following:

$$
f_t(S'' \cup \{(\bar{v}, \bar{m})\}) - f_t(S'') \le f_t(S' \cup \{(\bar{v}, \bar{m})\}) - f_t(S').
$$

We conclude that $f_t$ is a submodular set function.

**Monotonicity.** The function $f_t: 2^{V \times M} \to \mathbb{R}$ is monotone [72] if for every $S' \subset S'' \subset V \times M$ it holds that

$$
f_t(S'') \geq f_t(S'). \tag{33}
$$

We have

$$
f_t(S'') - f_t(S') = G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}'') - G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}') \quad (34)
$$

$$
= \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) (Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}^{\prime\prime}) - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}^{\prime})) \quad (35)
$$

$$
\ge 0. \tag{36}
$$

The last inequality is obtained using Eq. (28). We conclude that $f_t$ is monotone. $\square$

**Lemma A.2.** The set function $F_T : 2^{V \times M} \to \mathbb{R}$ in Eq. (26) is normalized (i.e., $f_T(\emptyset) = 0$), submodular, and monotone.

*Proof.* The set function $F_T$ is a nonnegative linear combination of submodular and monotone functions $f_t$ (see Lemma A.1), then $F_T$ is also submodular and monotone [72]. Moreover, each $f_t$ is normalized, then it follows that $F_T$ is normalized. $\square$

<figure><img src="image_11.png" /></figure>

APPENDIX B

NP-HARDNESS

**Theorem B.1.** *The problem of maximizing the time-averaged allocation gain (14) is NP-hard.*

*Proof.* We demonstrate the hardness of the problem by a reduction of the similarity caching problem, which is NP-hard [77] (a result that follows a reduction of the dominating set problem). We first define the similarity caching problem and characterize its inputs and variables.

**Similarity Caching Problem.** Consider the network topology comprising of a repository node and a single cache node. The repository node stores a catalog of files $\mathcal{N} = \{1, 2, \dots, \hat{N}\}$. A cache node can store a subset $\hat{\mathcal{S}} \subset \mathcal{N}$ of $\hat{k}$ files (i.e., $|\hat{\mathcal{S}}| = \hat{k}$), where $\hat{k} \in \{1, 2, \dots, |\hat{\mathcal{N}}| - 1\}$. The system incurs an approximation cost $C_a(\hat{r}, \hat{s}) \in \mathbb{R} \cup \{-\infty, +\infty\}$ when a request for file $\hat{r} \in \hat{\mathcal{N}}$ is satisfied by the cache serving file $\hat{s} \in \hat{\mathcal{N}}$. The approximation cost is null if $\hat{s} = \hat{r}$ (i.e., $C_a(\hat{r}, \hat{r}) = 0$ for all $\hat{r} \in \hat{\mathcal{N}}$). The system also incurs an additional retrieval cost $C_r \in \mathbb{R} \cup \{-\infty, +\infty\}$ to serve the request with a file stored at the repository node. For a given request $\hat{r} \in \hat{\mathcal{N}}$ and cache state $\hat{\mathcal{S}}$, the system incurs the following overall cost.

$$
C_{\text{overall}}(\hat{r}, \hat{\mathcal{S}}) \triangleq \min \left\{ C_r, \min \left\{ C_a(\hat{r}, \hat{s}) : \hat{s} \in \hat{\mathcal{S}} \right\} \right\}, \quad (37)
$$

where the term $\min \{C_a(\hat{r}, \hat{s}) : \hat{s} \in \hat{\mathcal{S}}\}$ signifies that the best approximating file stored at the cache is selected as a candidate to serve the request $\hat{r}$, and the term $\min \{C_r, \cdot\}$ signifies that when the cost of approximating the request with the candidate file exceeds the retrieval cost, the system serves the request by fetching an identical file from the repository and incurs a retrieval cost $C_r$. At timeslot $t \in \{1, 2, \dots, T\}$, the system receives a request $\hat{r}_t \in \hat{\mathcal{N}}$.

The static offline problem is formulated as follows.

$$
\hat{\mathcal{S}}_{\star} \in \arg \min_{\hat{\mathcal{S}} \subset \mathcal{N}, |\hat{\mathcal{S}}|=\hat{k}} \sum_{t=1}^{T} C_{\text {overall }}(\hat{r}_{t}, \hat{\mathcal{S}}) .
\qquad(38)
$$
Inputs. The similarity caching problem (38) takes the following inputs: $\mathcal{N} \in \{1, 2, \dots, \hat{N}\}$, $\hat{k} \in \{1, 2, \dots, |\mathcal{N} - 1|\}$, $C_r \in \mathbb{R} \cup \{-\infty + \infty\}$, $C_a : \mathcal{N} \times \mathcal{N} \to \mathbb{R} \cup \{-\infty + \infty\}$, $\{\hat{r}_1, \hat{r}_2, \dots, \hat{r}_T\} \in \mathcal{N}^T$.

Decision variable. The similarity caching problem (38) seeks a cache allocation $\hat{S} \subset \mathcal{N}$ such that $|\hat{S}| = \hat{k}$.

**Reduction of similarity caching problem.** We show that the similarity caching problem (38) can be reduced to the static model allocation problem in Eq. (14). We start observing that problem (14) is equivalent to

$$ \mathbf{x}_{*} \in \arg \min_{\mathbf{x} \in \mathcal{X}} \sum_{t=1}^{T} C(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}). \qquad (39) $$

The set of nodes is $\mathcal{V} = \{1, 2\}$, and the set of edges is $\mathcal{E} = \{(1, 2)\}$ (a topology with a single repository and a single cache, which we assume to be node 2 and node 1, respectively). The set of models is $\mathcal{M} = \mathcal{N}$ and the set of tasks is $\mathcal{N} = \mathcal{N}^T$. The set of request types is $\mathcal{R} = \mathcal{N} \times \{(1, 2)\}$. The maximum capacity is $L_m^v = 1$, and the potential available capacity is $l_{\rho,m}^{t,v} = 1$ for every node $v \in \mathcal{V}$, timeslot $t \in \{1, 2, \dots, T\}$, $\rho \in \mathcal{R}$, and model $m \in \mathcal{M}$. Given a request type $\rho = (i, \mathbf{p}) \in \mathcal{R}$, we define a cost $C_{\mathbf{p},m}^{p_j} \triangleq C_a(i, m)$ for node $p_j = 1$ and $C_{\mathbf{p},m}^{p_j} = C_r$ for node $p_j = 2$ for $m \in \mathcal{M}$. The allocation vector at node 1 is $\mathbf{x}^1 = (1_{\{m \in \mathcal{S}\}})_{m \in \mathcal{M}}$, and at node 2 is $\mathbf{x}^2 = (1)_{m \in \mathcal{M}}$. The model size is $s_m^v = 1$ for every $v \in \mathcal{V}$ and $m \in \mathcal{N}$. The allocation budget is $b_m^1 = \hat{k}$, and $b_m^2 = |\mathcal{M}|$. At timeslot $t$, a single request constitutes the request batch $\mathbf{r}_t = (1_{\{i=\hat{r}_t\}})_{(i, \mathbf{p}) \in \mathcal{R}}$.

The aggregate cost in Eq. (12) incurred by the system at time slot $t$ is given by

$$ C(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}) = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}} \gamma_{\rho}^{k} \min \left\{ r_{\rho}^{t} - \sum_{k'=1}^{k-1} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}), z_{\rho}^{k}(\mathbf{l}_t, \mathbf{x}) \right\} \cdot 1\{\sum_{k'=1}^{k-1} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) < r_{\rho}^{t}\} \quad (40) $$

$$ = \sum_{k=1}^{2\hat{N}} \gamma_{(\hat{r}_t, \mathbf{p})}^{k} \min \left\{ 1 - \sum_{k'=1}^{k-1} z_{(\hat{r}_t, \mathbf{p})}^{k'}(\mathbf{l}_t, \mathbf{x}), z_{(\hat{r}_t, \mathbf{p})}^{k}(\mathbf{l}_t, \mathbf{x}) \right\} \cdot 1\{\sum_{k'=1}^{k-1} z_{(\hat{r}_t, \mathbf{p})}^{k'}(\mathbf{l}_t, \mathbf{x}) < 1\} \quad (41) $$

$$ = \sum_{k=1}^{2\hat{N}} \gamma_{(\hat{r}_t, \mathbf{p})}^{k} \min \left\{ 1 - 0, z_{(\hat{r}_t, \mathbf{p})}^{k}(\mathbf{l}_t, \mathbf{x}) \right\} \cdot 1\{\sum_{k'=1}^{k-1} z_{(\hat{r}_t, \mathbf{p})}^{k'}(\mathbf{l}_t, \mathbf{x}) < 1\} \quad (42) $$

$$ = \sum_{k=1}^{2\hat{N}} \gamma_{(\hat{r}_t, \mathbf{p})}^{k} z_{(\hat{r}_t, \mathbf{p})}^{k}(\mathbf{l}_t, \mathbf{x}) \cdot 1\{\sum_{k'=1}^{k-1} z_{(\hat{r}_t, \mathbf{p})}^{k'}(\mathbf{l}_t, \mathbf{x}) < 1\} \quad (43) $$

$$ = \min \left\{ \gamma_{(\hat{r}_t, \mathbf{p})}^{k} : z_{(\hat{r}_t, \mathbf{p})}^{k}(\mathbf{l}_t, \mathbf{x}) = 1, k \in \{1, 2, \dots, 2\hat{N}\} \right\} \quad (44) $$

$$ = \min(\{C_r\} \cup \{C_a(\hat{r}_t, m) : m \in \hat{S}\}) \quad (45) $$

$$ = C_{\text{overall}}(\hat{r}_t, \hat{S}). \quad (46) $$

Equation (40) is the definition of $C(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x})$ that we provide in Eq. (12). The batch of requests consists only of a single request $\mathbf{r}_t = (1_{\{i=\hat{r}_t\}})_{(i, \mathbf{p}) \in \mathcal{R}}$; this gives Equation (41). Note that $K_{\rho} = 2\hat{N}$ since a request can be served with $\hat{N}$ different models at node 2 and node 1. Equation (42) can be understood as follows. Consider a model rank $k$ that makes the indicator function $1\{\sum_{k'=1}^{k-1} z_{(\hat{r}_t, \mathbf{p})}^{k'}(\mathbf{l}_t, \mathbf{x})<1\}$ true. This means that no model with lower rank $k' < k$ is able to satisfy the single request we are considering, which implies that the effective capacity is always $z_{(\hat{r}_t, \mathbf{p})}^{k'}(\mathbf{l}_t, \mathbf{x}) = 0$ for any $k' < k$. Equation (43) follows from $z_{(\hat{r}_t, \mathbf{p})}^{k'}(\mathbf{l}_t, \mathbf{x}) \in \{0, 1\}$ for $k' \in \{1, \dots, 2\hat{N}\}$. Equation (44) follows from the fact that 1) for each request $\rho$, model service costs are ordered in non-decreasing order (i.e., $\gamma_{(\hat{r}_t, \mathbf{p})}^{k} \le \gamma_{(\hat{r}_t, \mathbf{p})}^{k+1}$) and 2) the term $z_{(\hat{r}_t, \mathbf{p})}^{k}(\mathbf{l}_t, \mathbf{x})1\{\sum_{k'=1}^{k-1} z_{(\hat{r}_t, \mathbf{p})}^{k'}(\mathbf{l}_t, \mathbf{x})<1\}$ can only be non-zero for a single value of $k$; indeed, as we only need to satisfy one request, only one model will be involved in serving that request. Equation (45) follows from the definition of the costs $C_{\mathbf{p},m}^{p_j} = C_a(i, m)$ for $p_j = 1$ and $C_{\mathbf{p},m}^{p_j} = C_r$ for $p_j = 2$ for $m \in \mathcal{M}$. Equation (46) is a direct result from the definition of $C_{\text{overall}}$ in Eq. (37).

Hence, finding the optimal allocation $\mathbf{x}_*$ that minimizes the total allocation cost in (39) is equivalent to solve problem (38), which is NP-hard.

□

# APPENDIX C

## EQUIVALENT EXPRESSION OF THE GAIN FUNCTION

**Lemma C.1.** Let us fix the threshold $c \in \mathbb{N} \cup \{0\}$, request type $\rho \in \mathcal{R}$, model rank $k \in [K_{\rho}]$, time slot $t$, load vector $\mathbf{l}_t$ and allocation vector $\mathbf{x}$. For brevity, let us denote $z_{\rho}^{k'} = z_{\rho}^{k}(\mathbf{l}_t, \mathbf{x})$. The following formula holds:

$$ \min \left\{ c, \sum_{k'=1}^{k} z_{\rho}^{k'} \right\} - \min \left\{ c, \sum_{k'=1}^{k-1} z_{\rho}^{k'} \right\} = \min \left\{ c - \sum_{k'=1}^{k-1} z_{\rho}^{k'}, z_{\rho}^{k} \right\} \cdot 1 \left\{ \sum_{k'=1}^{k-1} z_{\rho}^{k'} < c \right\}. \quad (47) $$
Proof. We distinguish two cases:

(I) When the $k-1$ less costly models have at least $c$ effective capacity, i.e., $\sum_{k'=1}^{k-1} z_{\rho}^{k'} \ge c$, we obtain:

$$ \sum_{k'=1}^{k} z_{\rho}^{k'} = \sum_{k'=1}^{k-1} z_{\rho}^{k'} + z_{\rho}^{k} \ge c + z_{\rho}^{k} \ge c. \quad (48) $$

The last inequality is obtained using $z_{\rho}^{k} \ge 0$. Therefore, the left term of Eq. (47) becomes $c-c=0$, and the indicator function of the right term becomes zero. Hence, Eq. (47) is verified in this case.

(II) When $\sum_{k'=1}^{k-1} z_{\rho}^{k'} < c$, we obtain:

$$ \min \left\{ c, \sum_{k'=1}^{k} z_{\rho}^{k'} \right\} - \min \left\{ c, \sum_{k'=1}^{k-1} z_{\rho}^{k'} \right\} = \min \left\{ c, \sum_{k'=1}^{k} z_{\rho}^{k'} \right\} - \sum_{k'=1}^{k-1} z_{\rho}^{k'} = \min \left\{ c - \sum_{k'=1}^{k-1} z_{\rho}^{k'}, z_{\rho}^{k} \right\}. \quad (49) $$

Hence Eq. (47) is verified, being the indicator function equal to 1 in this case.

By combining Eq. (48) and Eq. (49) we obtain Eq. (47). □

**Lemma C.2.** The cost function given by Eq. (12) can be expressed as:

$$ C(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}) = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k} - \gamma_{\rho}^{k+1}) \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) \right\} + \gamma_{\rho}^{K_{\rho}} r_{\rho}^{t}. \quad (50) $$

Proof. The sum $\sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x})$ for $k = K_{\rho}$ surely includes a repository model as it sums all the models along the path of request type $\rho$; thus, we have $\sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) \ge r_{\rho}^{t}$ (see Eq. (9)) and

$$ \gamma_{\rho}^{K_{\rho}} \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{K_{\rho}} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) \right\} = \gamma_{\rho}^{K_{\rho}} r_{\rho}^{t}. \quad (51) $$

Now we use Lemma C.1 to express the cost function in Eq. (12) as a sum of the difference of min functions.

$$ C(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}) = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}} \gamma_{\rho}^{k} \cdot \min \left\{ r_{\rho}^{t} - \sum_{k'=1}^{k-1} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}), z_{\rho}^{k}(\mathbf{l}_t, \mathbf{x}) \right\} \cdot 1\{\sum_{k'=1}^{k-1} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) < r_{\rho}^{t}\} \quad (52) $$

$$ = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}} \gamma_{\rho}^{k} \left( \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) \right\} - \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k-1} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) \right\} \right) \quad (53) $$

$$ \stackrel{(51)}{=} \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} \gamma_{\rho}^{k} \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) \right\} - \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}} \gamma_{\rho}^{k} \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k-1} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) \right\} + \gamma_{\rho}^{K_{\rho}} r_{\rho}^{t} \quad (54) $$

$$ = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} \gamma_{\rho}^{k} \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) \right\} - \sum_{\rho \in \mathcal{R}} \sum_{k=2}^{K_{\rho}} \gamma_{\rho}^{k} \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k-1} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) \right\} + \gamma_{\rho}^{K_{\rho}} r_{\rho}^{t} \quad (55) $$

$$ = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} \gamma_{\rho}^{k} \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) \right\} - \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} \gamma_{\rho}^{k+1} \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) \right\} + \gamma_{\rho}^{K_{\rho}} r_{\rho}^{t} \quad (56) $$

$$ = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k} - \gamma_{\rho}^{k+1}) \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{x}) \right\} + \gamma_{\rho}^{K_{\rho}} r_{\rho}^{t}. \quad (57) $$

□

**Proof of Lemma III.1.**

Proof. By using the expression Eq. (50) for a generic allocation vector $\mathbf{x}$ and for $\omega$, we obtain:

$$ G(\mathbf{r}, \mathbf{l}_t, \mathbf{x}) = C(\mathbf{r}, \mathbf{l}_t, \omega) - C(\mathbf{r}, \mathbf{l}_t, \mathbf{x}) \quad (58) $$

$$ = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k} - \gamma_{\rho}^{k+1}) \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \omega) \right\} - \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k} - \gamma_{\rho}^{k+1}) \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \omega) \right\} \quad (59) $$

$$ = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k} - \gamma_{\rho}^{k+1}) \cdot \left\{ \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \omega) \right\} - \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \omega) \right\} \right\} \quad (60) $$
$$ = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) \cdot \left\{ \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_{t}, \mathbf{x}) \right\} - \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_{t}, \mathbf{\omega}) \right\} \right\} \quad (61) $$

$$ = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) (Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{x}) - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{\omega})). \quad (62) $$

□
**Algorithm 2** Weighted negative entropy Bregman projection onto the weighted capped simplex

**Input:** $|\mathcal{M}|; b^v; s^v$; Sorted $\mathbf{y}^v$ where $y_{|\mathcal{M}|+1}^v \leftarrow +\infty$
1: $\mathbf{y}_{|\mathcal{M}|+1}^v \leftarrow +\infty$
2: **for** $k \in \{|\mathcal{M}|, |\mathcal{M}| - 1, \dots, 1\}$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;$m_k \leftarrow \frac{b^v - \sum_{m=k+1}^{|\mathcal{M}|} s_m^v}{\sum_{m=1}^k s_m^v y_{k+1}^{v'}}$
4: &nbsp;&nbsp;&nbsp;&nbsp;**if** $y_k^{v'} m_k < 1 \le y_{k+1}^{v'} m_k$ **then**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $k' \in \{1, 2, \dots, k\}$ **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$y_{k'}^{v'} \leftarrow m_k y_{k'}^{v'}$
7: &nbsp;&nbsp;&nbsp;&nbsp;**for** $k' \in \{k+1, k+2, \dots, |\mathcal{M}|\}$ **do**
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$y_{k'}^{v'} \leftarrow 1$
9: &nbsp;&nbsp;&nbsp;&nbsp;**return** $\mathbf{y}^v$

$\triangleright$ Appropriate $k$ is found

$\triangleright$ Scale the variable's components

$\triangleright$ Cap the variable's components to 1

$\triangleright$ $\mathbf{y}^v$ is the result of the projection

## APPENDIX D
PROJECTION ALGORITHM

In order to project a fractional allocation $\mathbf{y}'$ lying outside the constraint set $\mathcal{Y}$ to a feasible allocation $\mathbf{y}$, we perform a Bregman projection associated to the global mirror map $\Phi: \mathcal{D} \to \mathbb{R}$, where the Bregman divergence associated to the mirror map $\Phi: \mathcal{D} \to \mathbb{R}$ is given by

$$D_{\Phi}(\mathbf{y}, \mathbf{y}') = \Phi(\mathbf{y}) - \Phi(\mathbf{y}') - \nabla\Phi(\mathbf{y}')^T (\mathbf{y} - \mathbf{y}'), \quad (63)$$

and the Bregman divergences associated to the mirror maps $\Phi^v: \mathcal{D}^v \to \mathbb{R}$ are also given by

$$D_{\Phi^v}(\mathbf{y}^v, \mathbf{y}'^v) = \Phi^v(\mathbf{y}^v) - \Phi^v(\mathbf{y}'^v) - \nabla\Phi^v(\mathbf{y}'^v)^T (\mathbf{y}^v - \mathbf{y}'^v). \quad (64)$$

The projection operation yields a constrained minimization problem, i.e.,

$$\mathbf{y} = \prod_{\mathcal{Y} \cap \mathcal{D}}^{\Phi} (\mathbf{y}') = \underset{\mathbf{y} \in \mathcal{Y} \cap \mathcal{D}}{\operatorname{argmin}} D_{\Phi}(\mathbf{y}, \mathbf{y}'). \quad (65)$$

The global mirror map $\Phi: \mathcal{D} \to \mathbb{R}$ is defined as the sum of the weighted negative entropy maps $\Phi^v: \mathcal{D}^v \to \mathbb{R}$, where $\mathcal{D} = \mathbb{R}_+^{\mathcal{V} \times \mathcal{M}}$ is the domain of $\Phi$, and the set $\mathcal{D}^v = \mathbb{R}_+^\mathcal{V}$ is the domain of $\Phi^v$ for all $v \in \mathcal{V}$. Thus, it follows that the global Bregman divergence is the sum of the Bregman divergences local to each node $v \in \mathcal{V}$, i.e.,

$$D_{\Phi}(\mathbf{y}, \mathbf{y}') = \sum_{v \in \mathcal{V}} D_{\Phi^v}(y^v, y'^v) \quad (66)$$

where $\mathbf{y} \in \mathcal{Y} = \mathcal{X}_{v \in \mathcal{V}} \mathcal{Y}^v$, and $\mathbf{y}^v \in \mathcal{Y}^v, \forall v \in \mathcal{V}$. In order to minimize the value $D_{\Phi}(\mathbf{y}, \mathbf{y}')$ for $\mathbf{y} \in \mathcal{Y} \cap \mathcal{D}$, we can independently minimize the values $D_{\Phi^v}(y^v, y'^v)$ for $\mathbf{y}^v \in \mathcal{Y}^v \cap \mathcal{D}^v$ giving $|\mathcal{V}|$ subproblems; for every $v \in \mathcal{V}$ we perform the following projection

$$\mathbf{y}^v = \prod_{\mathcal{Y}^v \cap \mathcal{D}^v}^{\Phi^v}(\mathbf{y}'^v) = \underset{\mathbf{y}^v \in \mathcal{Y}^v \cap \mathcal{D}^v}{\operatorname{argmin}} D_{\Phi^v}(\mathbf{y}^v, \mathbf{y}'^v). \quad (67)$$

**Theorem D.1.** Algorithm 2 when executed at node $v \in \mathcal{V}$ returns $\prod_{\mathcal{Y}^v \cap \mathcal{D}^v}^{\Phi^v}(\mathbf{y}^v)$, i.e., the projection of the vector $\mathbf{y}^v$ onto the weighted capped simplex $\mathcal{Y}^v \cap \mathcal{D}^v$ under the weighted negative entropy $\Phi^v(\mathbf{y}^v) = \sum_{m \in \mathcal{M}} s_m^v y_m^v \log(y_m^v)$. The time complexity of the projection is $O(|\mathcal{M}| \log(|\mathcal{M}|))$.

*Proof.*

$$\prod_{\mathcal{Y}^v \cap \mathcal{D}^v}^{\Phi^v}(\mathbf{y}'^v) = \underset{\mathbf{y}^v \in \mathcal{Y}^v \cap \mathcal{D}^v}{\operatorname{argmin}} D_{\Phi^v}(\mathbf{y}^v, \mathbf{y}'^v) \quad (68)$$

$$= \underset{\mathbf{y}^v \in \mathcal{Y}^v \cap \mathcal{D}^v}{\operatorname{argmin}} \sum_{m \in \mathcal{M}} s_m^v \left( y_m^v \log \left( \frac{y_m^v}{y_m^{v'}} \right) - y_m^v + y_m^{v'} \right). \quad (69)$$

We adapt the negative entropy projection algorithm in [70]. The constraints $y_m^v > 0, \forall m \in \mathcal{M}$ are implicitly enforced by the negentropy mirror map $\Phi^v$ and $D_{\Phi^v}(\mathbf{y}^v, \mathbf{y}'^v)$ is convex in $\mathbf{y}^v$. The Lagrangian function of the above problem:

$$\mathcal{J}(\mathbf{y}^v, \beta, \tau) = \sum_{m \in \mathcal{M}} s_m^v \left( y_m^v \log \left( \frac{y_m^v}{y_m^{v'}} \right) - y_m^v + y_m^{v'} \right) - \sum_{m \in \mathcal{M}} \beta_m (1 - y_m^v) - \tau \left( \sum_{m \in \mathcal{M}} s_m^v y_m^v - b^v \right). \quad (70)$$

At optimal point $\hat{\mathbf{y}}^v$ the following KKT conditions hold:

$$s_m^v \log(\hat{y}_m^v) - s_m^v \log(y_m^{v'}) + \beta_m - s_m^v \tau = 0, \quad (71a)$$
$$ \hat{y}_m^v \le 1, \qquad (71b) $$

$$ \beta_m \ge 0, \qquad (71c) $$

$$ \sum_{m \in \mathcal{M}} s_m^v \hat{y}_m^v = b^v, \qquad (71d) $$

$$ \beta_m(1 - \hat{y}_m^v) = 0. \qquad (71e) $$

Without loss of generality, assume the components of $\hat{y}^v$ are in non-decreasing order. Let $k$ be the index of the largest component of $\hat{y}^v$ strictly smaller than , i.e.,

$$ \hat{y}_1^v \le \dots \le \hat{y}_k^v < \hat{y}_{k+1}^v = \dots = \hat{y}_{|\mathcal{M}|}^v = 1 \quad \text{if } k < |\mathcal{M}|, \qquad (72) $$

$$ \hat{y}_1^v \le \hat{y}_2^v \le \dots \le \hat{y}_{|\mathcal{M}|}^v < 1 \quad \text{if } k = |\mathcal{M}|. \qquad (73) $$

The goal here is to identify a valid value for $k$ (number of components of $\hat{y}^v$ different from 1) and $\tau$. For now assume that $\tau$ is known, so a valid $k \in \mathcal{M}$ should satisfy the following:

* For $m_L = 1, \dots, k$, we have from (71e) that $\beta_{m_L} = 0$, and then from (71a), $s_{m_L}^v \log(y_{m_L}^{\prime v}) + s_{m_L}^v \tau = s_{m_L}^v \log(\hat{y}_{m_L}^v) < s_{m_L}^v \log(1) = 0$, and can be simplified to

$$ y_{m_L}^{\prime v} e^{\tau} < 1, \forall m_L \in \{1, \dots, k\}. \qquad (74) $$

* For $m_U = k + 1, \dots, |\mathcal{M}|$: as $\beta_{m_U} \ge 0$ from (71c), we get $0 = s_{m_U}^v \log(\hat{y}_{m_U}^v) = s_{m_U}^v \log(y_{m_U}^{\prime v}) - \beta_{m_U} + s_{m_U}^v \tau \le s_{m_U}^v \log(y_{m_U}^{\prime v}) + s_{m_U}^v \tau$, and can be simplified to

$$ y_{m_U}^{\prime v} e^{\tau} \ge 1, \forall m_U \in \{k+1, \dots, |\mathcal{M}|\}. \qquad (75) $$

Consider Eqs. (72) and (73), and since for $m_L \in \{1, \dots, k\}$ we have $y_{m_L}^{\prime v} e^{\tau} = \hat{y}_{m_L}^v$ (the order is preserved), then the conditions in Eq. (74) are

$$ y_1^{\prime v} e^{\tau} \le \dots \le y_k^{\prime v} e^{\tau} < 1. \qquad (76) $$

If the components of $\mathbf{y}^{\prime v}$ are ordered in ascending order, then it is enough to check if $y_k^{\prime v} e^{\tau} < 1$ holds for Eq. (76) to be true. Moreover, for $m_U \in \{k+1, \dots, |\mathcal{M}|\}$, we have $\hat{y}_{m_U}^v = 1$ and $y_{m_U}^{\prime v} e^{\tau} \ge 1$. Then, by taking $y_{|\mathcal{M}|+1}^{\prime v} \triangleq +\infty$ ($k$ can be equal to $|\mathcal{M}|$ as in Eq. (73)) it is enough to check with the smallest $y_{m_U}^{\prime v}$ to summarize all the conditions in Eq. (75). Thus, all the needed conditions can be further simplified to:

$$ y_k^{\prime v} e^{\tau} < 1 \le y_{k+1}^{\prime v} e^{\tau}. $$

Note that the r.h.s inequality is ignored when $k = |\mathcal{M}|$ by construction ($y_{|\mathcal{M}|+1}^{\prime v} = +\infty$).

Now we established how to verify if a given $k \in \mathcal{M}$ is valid, what remains is to give the expression of $\tau$ using the knapsack constraint in Eq. (71d):

$$ b^v = \sum_{m=1}^{|\mathcal{M}|} s_m^v \hat{y}_m^v = \sum_{m=k+1}^{|\mathcal{M}|} s_m^v + e^{\tau} \sum_{m=1}^k s_m^v y_m^{\prime v}. $$

For a given $k \in \mathcal{M}$, we define

$$ m_k \triangleq e^{\tau} = \frac{b^v - \sum_{m=k+1}^{|\mathcal{M}|} s_m^v}{\sum_{m=1}^k s_m^v y_m^{\prime v}}. \qquad (77) $$

Thus, a valid $k$ is the value satisfying the following inequalities (line 7 of Algorithm 2):

$$ y_k^{\prime v} m_k < 1 \le y_{k+1}^{\prime v} m_k. \qquad (78) $$

The appropriate $k$ satisfying the KKT conditions is contained in $\mathcal{M}$, and due to the sorting operation this gives total time complexity of $O(|\mathcal{M}| \log(|\mathcal{M}|))$ per iteration. In practice, the online mirror ascent method quickly sets irrelevant items in the fractional allocation vector $\mathbf{y}^{\prime v}$ very close to 0. Therefore, we can keep track only of items with a fractional value above a threshold $\epsilon > 0$, and the size of this subset is practically $\ll |\mathcal{M}|$. Therefore, the projection can be very efficient in practice. $\square$

## APPENDIX E
SUBGRADIENT EXPRESSION

**Lemma E.1.** The gain function in Eq. (16) has a subgradient $g_t$ at point $\mathbf{y}_t \in \mathcal{Y}$ given by

$$ g_t = \left[ \sum_{\rho \in \mathcal{R}} l_{\rho,m}^{t,v} \left( \gamma_{\rho}^{K_{\rho}^{*}(\mathbf{y}_t)} - C_{\mathbf{p},m}^{v} \right) \mathbf{1}_{\{\kappa_{\rho}(v,m) < K_{\rho}^{*}(\mathbf{y}_t)\}} \right]_{(v,m) \in \mathcal{V} \times \mathcal{M}} , \qquad (79) $$
where $K_{\rho}^{*}(\mathbf{y}_{t}) = \min \{k \in [K_{\rho} - 1] : \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_{t}, \mathbf{y}_{t}) \geq r_{\rho}^{t}\}$.

*Proof.* The function given by $Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{y}_{t}) = \min \{r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_{t}, \mathbf{y}_{t})\}$ is a minimum of two concave differentiable functions (a constant, and a linear function). We can characterize its subdifferential (set of all possible subgradients), using [78, Theorem 8.2], at point $\mathbf{y}_{t} \in \mathcal{Y}$ as

$$ \partial Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{y}_{t}) = \begin{cases} \left\{ \nabla\left(\sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_{t}, \mathbf{y}_{t})\right) \right\} & \text{if } Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{y}_{t}) < r_{\rho}^{t} & (\text{r.h.s. argument of the min is active}), \\ \operatorname{conv}\left(\left\{\mathbf{0}, \nabla\left(\sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_{t}, \mathbf{y}_{t})\right)\right\}\right) & \text{if } Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{y}_{t}) = r_{\rho}^{t} & (\text{both arguments of the min are active}), \\ \{\mathbf{0}\} & \text{otherwise} & (\text{l.h.s. argument of the min is active}), \end{cases} \tag{80} $$

where $\operatorname{conv}(\cdot)$ is the convex hull of a set, and the gradient $\nabla$ is given by $\nabla(\cdot) = [\frac{\partial}{\partial y_m^v}(\cdot)]_{(v,m) \in \mathcal{V} \times \mathcal{M}}$. The operator $\frac{\partial}{\partial y_m^v}(\cdot)$ is the partial derivative w.r.t. $y_m^v$ (not to be confused with the subdifferential notation).

We restrict ourselves to the valid subgradient $\tilde{g}_{\rho,t}^k \in \partial Z_{\rho}^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t)$ given by

$$ \tilde{g}_{\rho,t}^{k} = \begin{cases} \nabla(\sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{y}_t)) & \text{if } Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t) < r_{\rho}^{t}, \\ \mathbf{0}, & \text{otherwise.} \end{cases} \tag{81} $$

Note that for every $(v, m) \in \mathcal{V} \times \mathcal{M}$ we have

$$ \frac{\partial}{\partial y_m^v} \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{y}_t) = \sum_{k'=1}^{k} \frac{\partial}{\partial y_m^v} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{y}_t) \stackrel{(11)}{=} l_{\rho,m}^{t,v} \cdot \mathbf{1}_{\{\kappa_{\rho}(v,m) \le k\}}. $$

The indicator variable $\mathbf{1}_{\{\kappa_{\rho}(v,m) \le k\}}$ is introduced since the partial derivative of $\sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{y}_t)$ w.r.t. $y_m^v$ is non-zero only if model $m$ at node $v$ is among the $k$ best models to serve requests of type $\rho$ (in this case, the variable $y_m^v$ appears once in the summation).

We obtain from Eq. (81)

$$ \tilde{g}_{\rho,t,m}^{k,v} = \begin{cases} \mathbf{1}_{\{\kappa_{\rho}(v,m) \le k\}} \cdot \mathbf{1}_{\{\kappa_{\rho}(v,m) \le k\}} & \text{if } Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t) < r_{\rho}^{t}, \\ 0 & \text{otherwise.} \end{cases} \tag{82} $$

$$ = l_{\rho,m}^{t,v} \cdot \mathbf{1}_{\{\kappa_{\rho}(v,m) \le k \land Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t) < r_{\rho}^{t}\}}, \forall (v,m) \in \mathcal{V} \times \mathcal{M}. \tag{83} $$

By considering the subdifferential

$$ \partial G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t) = \partial \left( \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) (Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t) - Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \omega)) \right), \tag{84} $$

and using [79, Theorem 23.6], we get

$$ \partial G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t) = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} \partial \left( (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) (Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t) - Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \omega)) \right). \tag{85} $$

The constant factors $(\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k})$ are non-negative, so we can multiply both sides of the subgradient inequality by a non-negative constant [79, Sec. 23]; furthermore, the subgradient of the constants $Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \omega)$ is $\mathbf{0}$. We get

$$ \partial G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t) = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) \partial (Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t) - Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \omega)) = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) \partial Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t). \tag{86} $$

Then, a subgradient $g_t \in \partial G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t)$ at point $\mathbf{y}_t \in \mathcal{Y}$ is given by

$$ g_t = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) \tilde{g}_{\rho,t}^{k,v}. \tag{87} $$

The $(v, m)$-th component of the subgradient $g_t$ is

$$ g_{t,m}^{v} = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) \cdot \tilde{g}_{\rho,t,m}^{k,v} \tag{88} $$
$$ = \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_{\rho}-1} l_{\rho,m}^{t,v} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) \cdot \mathbb{1}_{\{\kappa_{\rho}(v,m) \le k \wedge Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t) < r_{\rho}^{t}\}} \quad (89) $$

$$ = \sum_{\rho \in \mathcal{R}} \sum_{k=\kappa_{\rho}(v,m)}^{K_{\rho}-1} l_{\rho,m}^{t,v} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) \cdot \mathbb{1}_{\{\sum_{k'=1}^{k} z_{\rho}^{k'} (\mathbf{l}_t, \mathbf{y}_t) < r_{\rho}^{t}\}} \quad (90) $$

$$ = \sum_{\rho \in \mathcal{R}} \sum_{k=\kappa_{\rho}(v,m)}^{K_{\rho}^{*}(\mathbf{y}_{t})-1} l_{\rho,m}^{t,v} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) \quad (91) $$

$$ = \sum_{\rho \in \mathcal{R}} l_{\rho,m}^{t,v} \cdot \sum_{k=\kappa_{\rho}(v,m)}^{K_{\rho}^{*}(\mathbf{y}_{t})-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}) \quad (92) $$

$$ = \sum_{\rho \in \mathcal{R}} l_{\rho,m}^{t,v} (\gamma_{\rho}^{K_{\rho}^{*}(\mathbf{y}_{t})} - \gamma_{\rho}^{\kappa_{\rho}(v,m)}) \cdot \mathbb{1}_{\{\kappa_{\rho}(v,m) < K_{\rho}^{*}(\mathbf{y}_{t})\}} \quad (93) $$

$$ \stackrel{(11)}{=} \sum_{\rho \in \mathcal{R}} l_{\rho,m}^{t,v} (\gamma_{\rho}^{K_{\rho}^{*}(\mathbf{y}_{t})} - C_{\mathbf{p},m}^{v}) \cdot \mathbb{1}_{\{\kappa_{\rho}(v,m) < K_{\rho}^{*}(\mathbf{y}_{t})\}}, \forall (v,m) \in \mathcal{V} \times \mathcal{M}, \quad (94) $$

where $K_{\rho}^{*}(\mathbf{y}_{t}) = \min \{k \in [K_{\rho}-1] : \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_{t}, \mathbf{y}_{t}) \ge r_{\rho}^{t}\}$.

□

# APPENDIX F
## SUPPORTING LEMMAS FOR THE PROOF OF THEOREM V.1

### A. Concavity of the Gain Function

**Lemma F.1.** The gain function given by Eq. (16) is concave over its domain $\mathcal{Y}$ of possible fractional allocations.

*Proof.* Since $\lambda_{\rho}^{k}$ is defined to be the k-th smallest cost for any $k \in [K_{\rho}]$ (see Eq. (11)), then the factors $\gamma_{\rho}^{k+1} - \gamma_{\rho}^{k}$ are always non-negative. Moreover $z_{\rho}^{k}(\mathbf{l}_{t}, \mathbf{y}) = y_{m}^{v} l_{\rho,m}^{t,v}$, where $v, m$ are such that $\kappa_{\rho}(v, m) = k$. Therefore, $Z_{\rho}^{k}(\mathbf{l}_{t}, \mathbf{y})$ is the minimum between a constant $r_{\rho}^{t}$ and a sum $\sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_{t}, \mathbf{x})$ of linear functions of $\mathbf{y}$. Such minimum is thus a concave function of $\mathbf{y}$. Therefore, the gain in Eq. (16) is a weighted sum with positive weights of concave functions in $\mathbf{y}$, which is concave. □

### B. Strong convexity of the Mirror Map

**Lemma F.2.** The global mirror map $\Phi(\mathbf{y}) = \sum_{v \in \mathcal{V}} \Phi(y^{v}) = \sum_{v \in \mathcal{V}} \sum_{m \in \mathcal{M}} s_{m}^{v} y_{m}^{v} \log(y_{m}^{v})$ defined over the domain $\mathcal{D} = \mathbb{R}_{>0}^{|\mathcal{M}| \times |\mathcal{V}|}$ is $\theta$-strongly convex w.r.t. the norm $\| \cdot \|_{l_{1}(\mathbf{s})}$ over $\mathcal{Y} \cap \mathcal{D}$, where

$$ \theta \triangleq \frac{1}{s_{\max}|\mathcal{V}|\|\mathcal{M}\|}, \quad (95) $$

$$ \|\mathbf{y}\|_{l_1(\mathbf{s})} \triangleq \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}} s_m^v |y_m^v| \quad (\text{weighted } l_1 \text{ norm}), \quad (96) $$

and $s_{\max} \triangleq \{s_m^v : (v, m) \in \mathcal{V} \times \mathcal{M}\}$ is the maximum model size. In words, this means that the mirror map $\Phi$'s growth is lower bounded by a quadratic with curvature $\frac{1}{s_{\max}|\mathcal{V}|\|\mathcal{M}\|}$.

*Proof.* We extend the proof of the strong convexity of the negative entropy w.r.t. to the $l_{1}$ norm over the simplex given in [80, Lemma 16]. The map $\Phi(\mathbf{y})$ is differentiable over $\mathcal{Y} \cap \mathcal{D}$, so a sufficient (and also necessary) condition for $\Phi(\mathbf{y})$ to be $\theta$-strongly convex w.r.t. $\| \cdot \|_{l_{1}(\mathbf{s})}$ is:

$$ (\nabla \Phi(\mathbf{y}') - \nabla \Phi(\mathbf{y}))^T (\mathbf{y}' - \mathbf{y}) \ge \theta \|\mathbf{y}' - \mathbf{y}\|_{l_1(\mathbf{s})}^2, \quad \forall \mathbf{y}', \mathbf{y} \in \mathcal{Y} \cap \mathcal{D}. \quad (97) $$

We have

$$ (\nabla \Phi(\mathbf{y}') - \nabla \Phi(\mathbf{y}))^T (\mathbf{y}' - \mathbf{y}) = \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}} s_m^v (\log(y_m^{v}) - \log(y_m^v))(y_m^{v} - y_m^v). \quad (98) $$

Take $\mu_{m}^{v} \triangleq s_{m}^{v}(\log(y_{m}^{v}) - \log(y_{m}^{v}))(y_{m}^{v} - y_{m}^{v})$, and note that $\mu_{m}^{v} \ge 0$ (because log is an increasing function).

$$ \|\mathbf{y}' - \mathbf{y}\|_{l_1(\mathbf{s})}^2 = \left( \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}} s_m^v |y_m^{v} - y_m^v| \right)^2 = \left( \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}: \mu_m^v \neq 0} \sqrt{\mu_m^v \frac{s_m^v |y_m^{v} - y_m^v|}{\sqrt{\mu_m^v}}} \right)^2 $$
$$
\begin{aligned}
&\le \left( \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}: \mu_m^v \ne 0} \mu_m^v \right) \left( \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}: \mu_m^v \ne 0} (s_m^v)^2 \frac{(y_m^v - y_m^v)^2}{\mu_m^v} \right) \\
&= \left( \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}: \mu_m^v \ne 0} s_m^v (\log(y_m^v) - \log(y_m^v))(y_m^v - y_m^v) \right) \left( \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}: \mu_m^v \ne 0} s_m^v \frac{y_m^v - y_m^v}{\log(y_m^v) - \log(y_m^v)} \right).
\end{aligned}
$$

The inequality is obtained using Cauchy–Schwarz inequality. Take $s_m^v \le s_{\max}$, $\forall (v, m) \in \mathcal{V} \times \mathcal{M}$, we obtain:

$$
\begin{aligned}
\sum_{(v,m) \in \mathcal{V} \times \mathcal{M}: \mu_m^v \neq 0} s_m^v \frac{y_m^v - y_m^v}{\log(y_m^v) - \log(y_m^v)} &\le s_{\max} \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}: \mu_m^v \neq 0} \frac{y_m^v - y_m^v}{\log(y_m^v) - \log(y_m^v)} \\
&= s_{\max} \sum_{v \in \mathcal{V}} \sum_{m \in \mathcal{M}: \mu_m^v \neq 0} \frac{y_m^v - y_m^v}{\log(y_m^v) - \log(y_m^v)} \\
&\le s_{\max} \sum_{v \in \mathcal{V}} \sum_{m \in \mathcal{M}} \frac{y_m^v + y_m^v}{2} \\
&\le s_{\max} |\mathcal{V}| |\mathcal{M}|.
\end{aligned}
$$

The second inequality is shown in [80, Eq. (A.16)]. We find that $\forall \mathbf{y}', \mathbf{y} \in \mathcal{Y} \cap \mathcal{D}$:

$$ \frac{1}{s_{\max}|\mathcal{V}|U} \|\mathbf{y}' - \mathbf{y}\|_{l_1(s)}^2 \le \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}: \mu_m^v \ne 0} s_m^v \log(y_m^v) - \log(y_m^v) = (\nabla \Phi(\mathbf{y}') - \nabla \Phi(\mathbf{y}))^T (\mathbf{y}' - \mathbf{y}). \quad (99) $$

The strong convexity constant $\theta$ is $\frac{1}{s_{\max}|\mathcal{V}| |\mathcal{M}|}$. $\square$

### C. Subgradient Bound

**Lemma F.3.** For any $(\mathbf{r}_t, \mathbf{l}_t) \in \mathcal{A}$, the subgradients $\mathbf{g}_t$ of the gain function in Eq. (16) at point $\mathbf{y}_t \in \mathcal{Y}$ are bounded under the norm $\|\cdot\|_{l_\infty(\frac{1}{s})}$ by $\sigma = \frac{RL_{\max}\Delta_C}{s_{\min}}$, where $s_{\min} \triangleq \min\{s_m^v : \forall (v, m) \in \mathcal{V} \times \mathcal{M}\}$, $L_{\max} \triangleq \max\{L_m^v : \forall (v, m) \in \mathcal{V} \times \mathcal{M}\}$, $R = |\mathcal{R}|$, and $\Delta_C \triangleq \max\{(\sum_{m \in \mathcal{M}} \omega_{m'}^{\nu(\mathbf{p})} C_{\mathbf{p},m'}^{\nu(\mathbf{p})}) - C_{\mathbf{p},m}^v : \forall (i, \mathbf{p}) \in \mathcal{R}, (v, m) \in \mathbf{p} \times \mathcal{M}\}$ is the maximum serving cost difference between serving at a repository node $\nu(\mathbf{p})$ and at any other node $v \in \mathbf{p}$. The norm $\|\cdot\|_{l_\infty(\frac{1}{s})}$ is defined as

$$ \|\mathbf{y}\|_{l_{\infty}(\frac{1}{s})} \triangleq \max \left\{ \frac{|y_m^v|}{s_m^v} : (v, m) \in \mathcal{V} \times \mathcal{M} \right\}. \quad (100) $$

*Proof.* We have for any $t \in [T]$

$$ \|\mathbf{g}_t\|_{l_\infty(\frac{1}{s})} = \max \left\{ \frac{|g_{t,m}^v|}{s_m^v}, \forall (v, m) \in \mathcal{V} \times \mathcal{M} \right\} \le \max \left\{ \frac{|g_{t,m}^v|}{s_{\min}}, \forall (v, m) \in \mathcal{V} \times \mathcal{M} \right\} \quad (101) $$

$$ \le \frac{L_{\max}}{s_{\min}} \max \left\{ \sum_{\rho \in \mathcal{R}} \left( \gamma_{\rho}^{K_{\rho}(\mathbf{y}_t)} - C_{\mathbf{p},m}^{v} \right) \cdot \mathbf{1}_{\{\kappa_{\rho}(v,m) < K_{\rho}^{*}(\mathbf{y}_t)\}}, \forall (v,m) \in \mathcal{V} \times \mathcal{M} \right\} \quad (102) $$

$$ \le \frac{L_{\max}R}{s_{\min}} \max\{\gamma_{\rho}^{K_{\rho}} - \gamma_{\rho}^{1}, \forall \rho \in \mathcal{R}\} \le \frac{L_{\max}R\Delta_C}{s_{\min}} = \sigma. \quad (103) $$

$\square$

### D. Dual Norm

**Lemma F.4.** $\|\cdot\|_{l_\infty(\frac{1}{s})}$ is the dual norm of $\|\cdot\|_{l_1(s)}$ defined in (100) and (96), respectively.

*Proof.* The dual norm $\|\cdot\|_*$ of $\|\cdot\|_{l_1(s)}$ is defined as (e.g., [69])

$$ \|\mathbf{z}\|_* \triangleq \sup_{\mathbf{y} \in \mathbb{R}^{\mathcal{V} \times \mathcal{M}}} \left\{ \mathbf{z}^T \mathbf{y} : \|\mathbf{y}\|_{l_1(s)} \le 1 \right\}, \forall \mathbf{z} \in \mathbb{R}^{\mathcal{V} \times \mathcal{M}}. \quad (104) $$

We thus need to show that $\|\mathbf{z}\|_{l_\infty(\frac{1}{s})} = \sup_{\mathbf{y} \in \mathbb{R}^{\mathcal{V} \times \mathcal{M}}} \{\mathbf{z}^T \mathbf{y} : \|\mathbf{y}\|_{l_1(s)} \le 1\}$, $\forall \mathbf{z} \in \mathbb{R}^{\mathcal{V} \times \mathcal{M}}$. Consider any two vectors $\mathbf{y}$ and $\mathbf{z}$ in $\mathbb{R}^{\mathcal{V} \times \mathcal{M}}$. We have

$$ \mathbf{z}^T \mathbf{y} = \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}} y_m^v z_m^v = \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}} (s_m^v y_m^v) \left( \frac{z_m^v}{s_m^v} \right) \le \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}} (s_m^v \cdot |y_m^v|) \left( \frac{|z_m^v|}{s_m^v} \right) \quad (105) $$
$$ \leq \left( \sum_{(v,m) \in \mathcal{V} \times \mathcal{M}} s_m^v |y_m^v| \right) \max \left\{ \frac{|z_m^v|}{s_m^v} : (v,m) \in \mathcal{V} \times \mathcal{M} \right\} = \|\mathbf{y}\|_{l_1(\mathbf{s})} \|\mathbf{z}\|_{l_\infty(\frac{1}{\mathbf{s}})}. \quad (106) $$

Observe that

$$ \mathbf{y}^T \mathbf{z} \leq \|\mathbf{z}\|_{l_{\infty}(\frac{1}{\mathbf{s}})}, \quad \forall \mathbf{y} : \|\mathbf{y}\|_{l_1(\mathbf{s})} \leq 1. \quad (107) $$

Let $(v_*, m_*) = \arg \max_{(v,m) \in \mathcal{V} \times \mathcal{M}} \{\frac{|z_m^v|}{s_m^v}\}$. The equality is achieved in (107) when $\mathbf{y}_* = \left[ \frac{\operatorname{sign}(z_m^v)}{s_m^v} \mathbf{1}_{\{(v,m)=(v_*,m_*)\}} \right]_{(v,m) \in \mathcal{V} \times \mathcal{M}}$. Note that $\|\mathbf{y}_*\|_{l_1(\mathbf{s})} = 1 \le 1$, then the supremum in (104) is attained for $\mathbf{y} = \mathbf{y}_*$ and has value $\|\mathbf{z}\|_{l_{\infty}(\frac{1}{\mathbf{s}})}$; therefore, $\|\cdot\|_{l_{\infty}(\frac{1}{\mathbf{s}})}$ is the dual norm of $\|\cdot\|_{l_1(\mathbf{s})}$. $\square$

## E. Bregman Divergence Bound

**Lemma F.5.** The value of the Bregman divergence $D_{\Phi}(\mathbf{y}, \mathbf{y}_1)$ in Eq. (63) associated with the mirror map $\Phi(\mathbf{y}) = \sum_{v \in \mathcal{V}} \Phi^v(\mathbf{y}^v) = \sum_{v \in \mathcal{V}} \sum_{m \in \mathcal{M}} s_m^v y_m^v \log(y_m^v)$ is upper bounded by a constant

$$ D_{\max} \triangleq \sum_{v \in \mathcal{V}} \min\{b^v, \|\mathbf{s}^v\|_1\} \log \left( \frac{\|\mathbf{s}^v\|_1}{\min\{b^v, \|\mathbf{s}^v\|_1\}} \right) \geq D_{\Phi}(\mathbf{y}, \mathbf{y}_1). \quad (108) $$

where $y_{1,m}^v = \frac{\min\{b^v, \|\mathbf{s}^v\|_1\}}{\|\mathbf{s}^v\|_1}$, $\forall (v,m) \in \mathcal{V} \times \mathcal{M}$ and $\mathbf{s}^v = [s_m^v]_{m \in \mathcal{M}}$ for every $v \in \mathcal{V}$.

*Proof.* We prove that $\mathbf{y}_1^v$ is the minimizer of $\Phi^v$ over $\mathcal{Y}^v$. As $\Phi^v$ is convex over $\mathcal{Y}^v$ and differentiable in $\mathbf{y}_1^v$, $\mathbf{y}_1^v$ is a minimizer if and only if $\nabla\Phi^v(\mathbf{y}_1^v)^T(\mathbf{y}_1^v - \mathbf{y}) \le 0$, $\forall \mathbf{y}^v \in \mathcal{Y}^v$ [69, Proposition 1.3] (first order optimality condition). Note that from the definition of $\mathcal{Y}^v$ (see Sec. IV) we have for any $\mathbf{y}^v \in \mathcal{Y}^v$

$$ \sum_{m \in \mathcal{M}} s_m^v y_m^v = \min\{b^v, \|\mathbf{s}^v\|_1\}. \quad (109) $$

Let $c = \frac{\min\{b^v, \|\mathbf{s}^v\|_1\}}{\|\mathbf{s}^v\|_1}$, we get

$$ \nabla \Phi^v (\mathbf{y}_1^v)^T (\mathbf{y}_1^v - \mathbf{y}) = \sum_{m \in \mathcal{M}} s_m^v (\log(c) + 1)(c - y_m^v) = (\log(c) + 1)(c \|\mathbf{s}^v\|_1 - \sum_{m \in \mathcal{M}} s_m^v y_m^v) \quad (110) $$

$$ \begin{align} & \stackrel{(109)}{=} (\log(c) + 1) \left( c \|\mathbf{s}^v\|_1 - \sum_{m \in \mathcal{M}} s_m^v y_m^v \right) = (\log(c) + 1) \left( \min\{b^v, \|\mathbf{s}^v\|_1\} - \min\{b^v, \|\mathbf{s}^v\|_1\} \right) \tag{111} \\ & = 0. \tag{112} \end{align} $$

We confirmed that $\mathbf{y}_1^v$ is a minimizer of $\Phi^v$ over $\mathcal{Y}^v$. We have $\Phi^v(\mathbf{y}^v) \le 0, \forall \mathbf{y}^v \in \mathcal{Y}^v$, and using the first order optimality condition we obtain

$$ D_{\Phi^v}(\mathbf{y}^v, \mathbf{y}_1^v) \stackrel{(64)}{=} \Phi^v(\mathbf{y}^v) - \Phi^v(\mathbf{y}_1^v) + \nabla\Phi^v(\mathbf{y}_1^v)^T(\mathbf{y}_1^v - \mathbf{y}^v) \le \Phi^v(\mathbf{y}^v) - \Phi^v(\mathbf{y}_1^v) \le -\Phi^v(\mathbf{y}_1^v) \quad (113) $$

$$ = \min\{b^v, \|\mathbf{s}^v\|_1\} \log \left( \frac{\|\mathbf{s}^v\|_1}{\min\{b^v, \|\mathbf{s}^v\|_1\}} \right). \quad (114) $$

Thus, we obtain

$$ \sum_{v \in \mathcal{V}} D_{\Phi^v}(\mathbf{y}^v, \mathbf{y}_1^v) \stackrel{(66)}{=} D_{\Phi}(\mathbf{y}, \mathbf{y}_1) \le \sum_{v \in \mathcal{V}} \min\{b^v, \|\mathbf{s}^v\|_1\} \log \left( \frac{\|\mathbf{s}^v\|_1}{\min\{b^v, \|\mathbf{s}^v\|_1\}} \right). \quad (115) $$

## F. Bounds on the Gain Function

Upper and lower bounds on the gain function in Eq. (16) will be established using the following bounding function

$$ \Lambda(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}) \triangleq \sum_{\rho \in \text{supp}(\mathbf{r}^t)} \sum_{k=1}^{K_\rho-1} (\gamma_\rho^{k+1} - \gamma_\rho^k) r_\rho^t \left( 1 - \prod_{k'=1}^k \left(1 - z_\rho^{k'} (\mathbf{l}_t, \mathbf{y})/r_\rho^t \right) \right) \mathbf{1}_{\{Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega})=0\}}, \forall \mathbf{y} \in \mathcal{X} \cup \mathcal{Y}, \quad (116) $$

where

$$ \text{supp}(\mathbf{r}^t) \triangleq \{\rho \in \mathcal{R} : r_\rho^t \neq 0\} \quad (117) $$

is the set of request types for which there is a non-zero number of requests in the request batch $\mathbf{r}^t$.
**Lemma F.6.** The gain function in Eq. (16) can be equivalently expressed as

$$G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}) = \sum_{\rho \in \text{supp}(\mathbf{r}^t)} \sum_{k=1}^{K_\rho-1} (\gamma_\rho^{k+1} - \gamma_\rho^k) \min \left\{ r_\rho^t, \sum_{k'=1}^k z_\rho^{k'}(\mathbf{l}_t, \mathbf{y}) \right\} \mathbf{1}_{\{Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega})=0\}, \forall \mathbf{y} \in \mathcal{X} \cup \mathcal{Y}. \quad (118)$$

*Proof.* Remember from the definition in Eq. (15) that $Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}) = \min\{r_\rho^t, \sum_{k'=1}^k z_\rho^{k'}(\mathbf{l}_t, \mathbf{y})\}$. We observe that $Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega})$ is not a function of $\mathbf{y}$ and it is equal to 0, when there is no repository with model's rank smaller or equal to $k$, and to $r_\rho^t$, otherwise; therefore, $Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega}) \in \{0, r_\rho^t\}$.

When $Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega}) \neq 0$, and thus $Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega}) = r_\rho^t$, the following holds

$$Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}) - Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega}) = \min \left\{ r_\rho^t, \sum_{k'=1}^k z_\rho^{k'}(\mathbf{l}_t, \mathbf{y}) \right\} - Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega}) \quad (119)$$

$$= \min \left\{ 0, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_t, \mathbf{y}) - Z_{\rho}^{k}(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega}) \right\} = 0. \quad (120)$$

The last equality holds because $\sum_{k'=1}^k z_\rho^{k'}(\mathbf{l}_t, \mathbf{y}) - Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega}) \ge 0, \forall \mathbf{y} \in \mathcal{X} \cup \mathcal{Y}$ from Eq. (3). Otherwise, when $Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega}) = 0$, we have $Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}) - Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega}) = Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}) \stackrel{(15)}{=} \min \{r_\rho^t, \sum_{k'=1}^k z_\rho^{k'}(\mathbf{l}_t, \mathbf{y})\}$.

Hence, we can succinctly write, for any value of $Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega})$, that

$$Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{y}) - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{\omega}) = \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_{t}, \mathbf{y}) \right\} \mathbf{1}_{\{Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{\omega})=0\}}. \quad (121)$$

Since $\sum_{k'=1}^k z_\rho^{k'}(\mathbf{l}_t, \mathbf{y})$ is non negative, the previous formula implies that

$$Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{y}) - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{\omega}) = 0, \quad \text{if } r_{\rho}^{t} = 0. \quad (122)$$

By combining Eq. (121) and (122) that

$$Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{y}) - Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{\omega}) = \min \left\{ r_{\rho}^{t}, \sum_{k'=1}^{k} z_{\rho}^{k'}(\mathbf{l}_{t}, \mathbf{y}) \right\} \mathbf{1}_{\{Z_{\rho}^{k}(\mathbf{r}_{t}, \mathbf{l}_{t}, \mathbf{\omega})=0 \land r_{\rho}^{t} \neq 0\}}. \quad (123)$$

Hence, applying the above equalities on the gain expression we get

$$G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}) \stackrel{(16)}{=} \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_\rho-1} (\gamma_\rho^{k+1} - \gamma_\rho^k) (Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}) - Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega})) \quad (124)$$

$$\stackrel{(123)}{=} \sum_{\rho \in \mathcal{R}} \sum_{k=1}^{K_\rho-1} (\gamma_\rho^{k+1} - \gamma_\rho^k) \min \left\{ r_\rho^t, \sum_{k'=1}^k z_\rho^{k'}(\mathbf{l}_t, \mathbf{y}) \right\} \mathbf{1}_{\{Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega})=0 \land r_\rho^t \neq 0\}} \quad (125)$$

$$\stackrel{(117)}{=} \sum_{\rho \in \text{supp}(\mathbf{r}^t)} \sum_{k=1}^{K_\rho-1} (\gamma_\rho^{k+1} - \gamma_\rho^k) \min \left\{ r_\rho^t, \sum_{k'=1}^k z_\rho^{k'}(\mathbf{l}_t, \mathbf{y}) \right\} \mathbf{1}_{\{Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{\omega})=0\}}. \quad (126)$$

□

**Lemma F.7.** Consider $n \in \mathbb{N}$, $\mathbf{y} \in [0, 1]^n$, $\mathbf{q} \in \mathbb{N}^n$, and $c \in \mathbb{N}$. We assume that $q_i \le c, \forall i \in [n]$. The following holds

$$\min \left\{ c, \sum_{i \in [n]} y_i q_i \right\} \ge c - c \prod_{i \in [n]} (1 - y_i q_i / c). \quad (127)$$

*Proof.* We define $a_n \triangleq c - c \prod_{i \in [n]} (1 - y_i q_i / c)$ and $b_n \triangleq \min\{c, \sum_{i \in [n]} y_i q_i\}$.

We first show by induction that, if $a_n \le b_n$, then this inequality holds also for $n + 1$. Base case ($n = 1$).

$$a_1 = c - c + y_1 q_1 = y_1 q_1 = \min\{c, q_1 y_1\} = b_1. \quad (128)$$

Induction step.

$$a_{n+1} = c - c \prod_{i \in [n+1]} (1 - y_i q_i / c) \quad (129)$$
$$= c - c \prod_{i \in [n]} (1 - y_i q_i / c) (1 - y_{n+1} q_{n+1} / c) \qquad (130)$$

$$= c - c \prod_{i \in [n]} (1 - y_i q_i / c) + (c y_{n+1} q_{n+1} / c) \prod_{i \in [n]} (1 - y_i q_i / c) \qquad (131)$$

$$= a_n + y_{n+1} q_{n+1} \prod_{i \in [n]} (1 - y_i q_i / c) \qquad (132)$$

$$\le a_n + y_{n+1} q_{n+1}. \qquad (133)$$

The last inequality holds since by construction $q_i \le c$ and thus $0 \le y_i q_i / c \le 1$, and $0 \le \prod_{i \in [n]} (1 - y_i q_i / c) \le 1$. For the same reason, $0 \le \prod_{i \in [n+1]} (1 - y_i q_i / c) \le 1$ and thus, by (129), we have $a_{n+1} \le c$. Moreover, note that if $b_n = c$ then $b_{n+1} = c$. Therefore:

$$a_{n+1} \le \min \{c, a_n + y_{n+1} q_{n+1}\} \le \min \{c, b_n + y_{n+1} q_{n+1}\} \qquad (134)$$

$$= \begin{cases} \min \{c, \sum_{i=1}^{n+1} y_i q_i\} = b_{n+1}, & \text{if } b_n \le c, \\ \min \{c, c + y_{n+1} q_{n+1}\} = c = b_{n+1}, & \text{if } b_n = c, \end{cases} \qquad (135)$$

and the proof by induction is completed. □

**Lemma F.8.** Consider $\mathbf{y} \in [0, 1]^n$, $\mathbf{q} \in \mathbb{N}^n$, $c \in \mathbb{N}$ and $n \in \mathbb{N}$. We assume that $q_i \le c, \forall i \in [n]$. The following holds

$$c - c \prod_{i \in [n]} (1 - y_i q_i / c) \ge (1 - 1/e) \min \left\{ c, \sum_{i \in [n]} y_i q_i \right\}. \qquad (136)$$

*Proof.* Our proof follows the same lines of the proof of [81, Lemma 3.1]. We use the arithmetic/geometric mean inequality [82] on the non-negative variables $1 - y_i q_i / c, i \in [n]$ to obtain:

$$\frac{1}{n} \sum_{i \in [n]} (1 - y_i q_i / c) \ge \left( \prod_{i \in [n]} (1 - y_i q_i / c) \right)^{\frac{1}{n}}. \qquad (137)$$

We reformulate the above as:

$$1 - \prod_{i \in [n]} (1 - y_i q_i / c) \ge 1 - \left( \frac{1}{n} \sum_{i \in [n]} (1 - y_i q_i / c) \right)^n = 1 - \left( 1 - \frac{1}{n} \sum_{i \in [n]} y_i q_i / c \right)^n \qquad (138)$$

$$\ge 1 - \left( 1 - \frac{1}{n} \min \left\{ 1, \sum_{i \in [n]} y_i q_i / c \right\} \right)^n . \qquad (139)$$

To obtain the last inequality, consider that, for any number $z$, we have $z \ge \min\{1, z\}$, and thus $\sum_{i \in [n]} y_i q_i / c \ge \min\{1, \sum_{i \in [n]} y_i q_i / c\}$.

The function $f(z) = 1 - (1 - z/n)^n$ is concave for $z \in [0, 1]$, then, for $z \in [0, 1]$, $f(z) \ge f(0) + z \frac{f(1) - f(0)}{1 - 0} = zf(1)$, as $f(0) = 0$. Setting $z = \min\{1, \sum_{i \in [n]} y_i q_i / c\}$, we obtain the following:

$$1 - \prod_{i \in [n]} (1 - y_i q_i / c)^{(139)} \ge 1 - \left(1 - \frac{1}{n}z\right)^n \ge (1 - (1 - 1/n)^n)z \ge (1 - 1/e)z. \qquad (140)$$

The last inequality is obtained since $1 - (1 - 1/n)^n$ decreases in $n$, and it is lower bounded by $1 - 1/e$. By multiplying both sides of the above inequality by $c \in \mathbb{N}$, and replacing $z$ with its value we conclude the proof. □

**Lemma F.9.** for any request batch $\mathbf{r}_t$ and potential available capacity $\mathbf{l}_t$ such that $(\mathbf{r}_t, \mathbf{l}_t) \in \mathcal{A}$, the allocation gain $G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y})$ has the following lower and upper bounds

$$(1 - \frac{1}{e})^{-1} \Lambda(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}) \ge G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}) \ge \Lambda(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}), \qquad \forall \mathbf{y} \in \mathcal{X} \cup \mathcal{Y}. \qquad (141)$$

*Proof.* We have the following

$$G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}) \stackrel{(118)}{=} \sum_{\rho \in \text{supp}(\mathbf{r}^t)} \sum_{k=1}^{K_{\rho}-1} (\gamma_{\rho}^{k+1} - \gamma_{\rho}^k) \min \left\{ r_{\rho}^t, \sum_{k'=1}^{k} z_{\rho}^{k'} (\mathbf{l}_t, \mathbf{y}) \right\} \mathbf{1}_{\{Z_{\rho}^k(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y})=0\}} \qquad (142)$$
$$ \overset{(127)}{\ge} \sum_{\rho \in \text{supp}(\boldsymbol{r}^t)} \sum_{k=1}^{K_\rho-1} (\gamma_\rho^{k+1} - \gamma_\rho^k) r_\rho^t \left( 1 - \prod_{k'=1}^k \left(1 - \frac{z_{\rho,k'}(\boldsymbol{l}_t, \boldsymbol{y})}{r_\rho^t}\right) \right) \mathbb{1}_{\{Z_\rho^k(\boldsymbol{r}_t, \boldsymbol{l}_t, \boldsymbol{\omega})=0\}} \quad (143) $$

$$ \overset{(116)}{=} \Lambda(\boldsymbol{r}_t, \boldsymbol{l}_t, \boldsymbol{y}), \forall \boldsymbol{y} \in \mathcal{X} \cup \mathcal{Y}, \quad (144) $$

and

$$ (1 - 1/e)G(\boldsymbol{r}_t, \boldsymbol{l}_t, \boldsymbol{y}) \overset{(118)}{=} (1 - 1/e) \sum_{\rho \in \text{supp}(\boldsymbol{r}^t)} \sum_{k=1}^{K_\rho-1} (\gamma_\rho^{k+1} - \gamma_\rho^k) \min \left\{ r_\rho^t, \sum_{k'=1}^k z_\rho^{k'}(\boldsymbol{l}_t, \boldsymbol{y}) \right\} \mathbb{1}_{\{Z_\rho^k(\boldsymbol{r}_t, \boldsymbol{l}_t, \boldsymbol{\omega})=0\}} \quad (145) $$

$$ \overset{(136)}{\le} \sum_{\rho \in \text{supp}(\boldsymbol{r}^t)} \sum_{k=1}^{K_\rho-1} (\gamma_\rho^{k+1} - \gamma_\rho^k) r_\rho^t \left( 1 - \prod_{k'=1}^k \left(1 - \frac{z_{\rho,k'}(\boldsymbol{l}_t, \boldsymbol{y})}{r_\rho^t}\right) \right) \mathbb{1}_{\{Z_\rho^k(\boldsymbol{r}_t, \boldsymbol{l}_t, \boldsymbol{\omega})=0\}} \quad (146) $$

$$ \overset{(116)}{=} \Lambda(\boldsymbol{r}_t, \boldsymbol{l}_t, \boldsymbol{y}), \forall \boldsymbol{y} \in \mathcal{X} \cup \mathcal{Y}. \quad (147) $$

Inequalities in Eq. (143) and Eq. (146) follow from Eq. (127) and Eq. (136), respectively, by replacing $c = r_\rho^t$, $q_{k'} = \lambda_\rho^{k'}(\boldsymbol{l}_t)$, $x_{k'} = z_{\rho,k'}(\boldsymbol{l}_t, \boldsymbol{y})/\lambda_\rho^{k'}(\boldsymbol{l}_t)$, and $n = k$. $\square$

**Lemma F.10.** Let the allocation $\boldsymbol{x}^v$ be the random output of DEPROUND on node $v \in \mathcal{V}$ given the fractional allocation $\boldsymbol{y}^v \in \mathcal{Y}^v$. For any subset of the model catalog $S \subset \mathcal{M}$ and any number $c_m \in [0, 1], \forall m \in S$, DEPROUND satisfies the following:

$$ \mathbb{E} \left[ \prod_{m \in S} (1 - x_m^v c_m) \right] \le \prod_{m \in S} (1 - y_m^v c_m). \quad (148) $$

*Proof.* DEPROUND uses a subroutine SIMPLIFY, which, given input variables $y_m, y_{m'} \in (0, 1)$, outputs $x_m, x_{m'} \in [0, 1]$ with at least one of them being integral (0 or 1). Note that the input to SIMPLIFY is never integral since it is only called on fractional and yet unrounded variables. The property (B3) in [71, Lemma 2.1] implies that the output variables $x_m$ and $x_{m'}$ satisfy the following inequality:

$$ \mathbb{E}[x_m x_{m'}] \le y_m y_{m'}. \quad (149) $$

We have for any $c_m, c_{m'} \in [0, 1]$:

$$ \mathbb{E}[(1 - x_m c_m)(1 - x_{m'} c_{m'})] = \mathbb{E}[1 - x_m c_m - x_{m'} c_{m'} + x_m x_{m'} c_m c_{m'}] \quad (150) $$

$$ = 1 - y_m c_m - y_{m'} c_{m'} + \mathbb{E}[x_m x_{m'}] c_m c_{m'} \quad (151) $$

$$ \le 1 - y_m c_m - y_{m'} c_{m'} + y_m y_{m'} c_m c_{m'} \quad (152) $$

$$ = (1 - y_m c_m)(1 - y_{m'} c_{m'}). \quad (152) $$

where the second equality is obtained recalling that, by construction, $\mathbb{E}[x_m^v] = y_m, \forall m \in \mathcal{M}$ (Sec. IV-C).

Thus, the two fractional inputs $y_m$ and $y_{m'}$ to the SIMPLIFY subroutine, return $x_m$ and $x_{m'}$ satisfying the property (152). By induction as in the proof in [71, Lemma 2.2], we obtain for any $S \subset \mathcal{M}$:

$$ \mathbb{E} \left[ \prod_{m \in S} (1 - x_m c_m) \right] \le \prod_{m \in S} (1 - y_m c_m). \quad (153) $$

Note that the above property is satisfied with equality if the components of $\boldsymbol{x} \in \{0, 1\}^{|\mathcal{M}|}$ are sampled independently with $\mathbb{E}[x_m] = y_m$. $\square$

**Lemma F.11.** Let the allocation $\boldsymbol{x}^v$ be the random output of DEPROUND on node $v \in \mathcal{V}$ given the fractional allocation $\boldsymbol{y}^v \in \mathcal{Y}^v$. The following holds

$$ \mathbb{E}[\Lambda(\boldsymbol{r}_t, \boldsymbol{l}_t, \boldsymbol{x}_t)] \ge \Lambda(\boldsymbol{r}_t, \boldsymbol{l}_t, \boldsymbol{y}_t). \quad (154) $$

*Proof.* For any $k'$, assume $m \in \mathcal{M}$ and $v \in \mathcal{V}$ are such that $\kappa_\rho(v, m) = k'$ (Sec. III-E). Since $z_\rho^{k'}(\boldsymbol{l}_t, \boldsymbol{y}) = y_m^v l_{\rho,m}^{t,v}$ for a given $(m, v) \in \mathcal{M} \times \mathcal{V}$, then $z_{\rho,k'}(\boldsymbol{l}_t, \boldsymbol{y})/r_\rho^t$ can be written as:

$$ z_{\rho,k'}(\boldsymbol{l}_t, \boldsymbol{y})/r_\rho^t = \frac{l_{\rho,m}^{t,v}}{r_\rho^t} y_m^v. \quad (155) $$

where $l_{\rho,m}^{t,v}$ is a constant in $[0, 1]$ ($l_{\rho,m}^{t,v} \le \min\{r_\rho^t, L_m^v\}$ - see Sec. III-D) that scales variable $y_m^v$; therefore, by applying Lemma F.10, we obtain the following upper bound on the bounding function. Consider for all $v \in \mathcal{V}$ and $t \in [T]$ that $\boldsymbol{x}_t^v$ is
the random allocation obtained by running DEPROUND on the fractional allocation $y_t^v$, then

$$ \mathbb{E}[\Lambda(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}_t)] \stackrel{(116)}{=} \mathbb{E} \left[ \sum_{\rho \in \text{supp}(\mathbf{r}^t)} \sum_{k=1}^{K_\rho-1} (\gamma_\rho^{k+1} - \gamma_\rho^k) r_\rho^t \left( 1 - \prod_{k'=1}^k (1 - z_{\rho,k'}(\mathbf{l}_t, \mathbf{x}_t)/r_\rho^t) \right) \mathbb{1}_{\{Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \omega)=0\}} \right] \quad (156) $$

$$ = \sum_{\rho \in \text{supp}(\mathbf{r}^t)} \sum_{k=1}^{K_\rho-1} (\gamma_\rho^{k+1} - \gamma_\rho^k) r_\rho^t \left( 1 - \mathbb{E} \left[ \prod_{k'=1}^k (1 - z_{\rho,k'}(\mathbf{l}_t, \mathbf{x}_t)/r_\rho^t) \right] \right) \mathbb{1}_{\{Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \omega)=0\}} \quad (157) $$

$$ \geq \sum_{\rho \in \text{supp}(\mathbf{r}^t)} \sum_{k=1}^{K_\rho-1} (\gamma_\rho^{k+1} - \gamma_\rho^k) r_\rho^t \left( 1 - \prod_{k'=1}^k (1 - z_{\rho,k'}(\mathbf{l}_t, \mathbf{y}_t)/r_\rho^t) \right) \mathbb{1}_{\{Z_\rho^k(\mathbf{r}_t, \mathbf{l}_t, \omega)=0\}} \quad (158) $$

$$ = \Lambda(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t). \quad (159) $$

The equality is obtained using the linearity of the expectation, and the inequality is obtained by applying directly Lemma F.10. □

<figure><img src="image_12.png" /></figure>

# APPENDIX G

## PROOF OF THEOREM V.1

*Proof.* To prove the $\psi$-regret guarantee: (i) we first establish an upper bound on the regret of the INFIDA policy over its fractional allocations domain $\mathcal{Y}$ against a fractional optimum, then (ii) we use it to derive a corresponding $\psi$-regret guarantee over the integral allocations domain $\mathcal{X}$.

**Fractional domain regret guarantee.** To establish the regret guarantee of running Algorithm 1 at the level of each computing node $v \in \mathcal{V}$, we showed that the following properties hold:

1) The function $G$ is concave over its domain $\mathcal{Y}$ (Lemma F.1).

2) The mirror map $\Phi: \mathcal{D} \to \mathbb{R}$ is $\theta$-strongly convex w.r.t. the norm $\|\cdot\|_{l_1(\mathbf{s})}$ over $\mathcal{Y} \cap \mathcal{D}$, where $\theta$ is equal to Eq. (95) (Lemma F.2).

3) The gain function $G: \mathcal{Y} \to \mathbb{R}$ is $\sigma$-Lipchitz w.r.t $\|\cdot\|_{l_1(\mathbf{s})}$: the subgradients are bounded under the norm $\|\cdot\|_{l_\infty(\frac{1}{\sigma})}$ by $\sigma$, i.e., the subgradient of $G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y})$ at point $\mathbf{y}_t \in \mathcal{Y}$ is upper bounded ($\|\mathbf{g}_t\|_{l_\infty(\frac{1}{\sigma})} \le \sigma$) for any $(\mathbf{r}_t, \mathbf{l}_t) \in \mathcal{A}$ (Lemma F.3).

4) $\|\cdot\|_{l_\infty(\frac{1}{\sigma})}$ is the dual norm of $\|\cdot\|_{l_1(\mathbf{s})}$ (Lemma F.4).

5) The Bregman divergence $D_{\Phi}(\mathbf{y}_*, \mathbf{y}_1)$ in Eq. (63) is upper bounded by a constant $D_{\max}$ where $\mathbf{y}_* = \arg\max_{\mathbf{y} \in \mathcal{Y}} \sum_{t=1}^T G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y})$ and $\mathbf{y}_1 = \arg\min_{\mathbf{y} \in \mathcal{Y} \cap \mathcal{D}} \Phi(\mathbf{y})$ is the initial allocation (Lemma F.5).

Because of properties 1–5 above, the following bound holds for the regret of INFIDA over its fractional domain $\mathcal{Y}$ (vector field point of view of Mirror Descent in [69, Sec. 4.2] combined with [69, Theorem 4.2]):

$$ \text{Regret}_{T, \mathcal{Y}} = \sup_{\{(\mathbf{r}_t, \mathbf{l}_t)\}_{t=1}^T \in \mathcal{A}^T} \left\{ \sum_{t=1}^T G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_*) - \sum_{t=1}^T G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t) \right\} \quad (160) $$

$$ \leq \frac{D_{\Phi}(\mathbf{y}_*, \mathbf{y}_1)}{\eta} + \frac{\eta}{2\theta} \sum_{t=1}^{T} \|\mathbf{g}_t\|_{l_{\infty}(\frac{1}{\sigma})}^2 \leq \frac{D_{\max}}{\eta} + \frac{\eta\sigma^2 T}{2\theta}. \quad (161) $$

where $\eta$ is the learning rate of INFIDA (Algorithm 1, line 6). By selecting the learning rate $\eta = \frac{1}{\sigma} \sqrt{\frac{2\theta D_{\max}}{T}}$ giving the tightest upper bound we obtain

$$ \text{Regret}_{T, \mathcal{Y}} \le \sigma \sqrt{\frac{2D_{\max}}{\theta} T}. \quad (162) $$

**Integral domain regret guarantee.** Note that, by restricting the maximization to the subset of integral allocations $\mathbf{x} \in \mathcal{X}$, the optimal allocation $\mathbf{x}_* = \arg\max_{\mathbf{x} \in \mathcal{X}} \sum_{t=1}^T G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x})$ can only lead to a lower gain, i.e.,

$$ \sum_{t=1}^{T} G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}_*) \le \sum_{t=1}^{T} G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_*). \quad (163) $$

By taking $\psi = 1 - \frac{1}{e}$, and using the bounding function $\Lambda$ defined in Eq. (116), with the expectation taken over the random choices of the policy (DEPROUND at line 8 in Algorithm 1) we obtain

$$ \begin{aligned} \mathbb{E} \left[ \sum_{t=1}^{T} G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}_t) \right] & \stackrel{(141)}{\ge} \mathbb{E} \left[ \sum_{t=1}^{T} \Lambda(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}_t) \right] \stackrel{(154)}{\ge} \sum_{t=1}^{T} \Lambda(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t) \stackrel{(141)}{\ge} \psi \sum_{t=1}^{T} G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_t) \\ & \stackrel{(162)}{\ge} \psi \sum_{t=1}^{T} G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y}_*) - \psi \sigma \sqrt{\frac{2D_{\max}}{\theta} T} \stackrel{(163)}{\ge} \psi \sum_{t=1}^{T} G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}_*) - \psi \sigma \sqrt{\frac{2D_{\max}}{\theta} T}. \end{aligned} \quad (164) $$
Thus, we have

$$ \psi \sum_{t=1}^{T} G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}_*) - \mathbb{E} \left[ \sum_{t=1}^{T} G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}_t) \right] \leq \psi \sigma \sqrt{\frac{2D_{\max}}{\theta}} T. \quad (165) $$

The above inequality holds for any sequence $\{(\mathbf{r}_t, \mathbf{l}_t)\}_{t=1}^T \in \mathcal{A}^T$. Thus, the $\psi$-regret is given by

$$ \psi\text{-Regret}_{T, \mathcal{A}} \stackrel{(20)}{=} \sup_{\{(\mathbf{r}_t, \mathbf{l}_t)\}_{t=1}^T \in \mathcal{A}^T} \left\{ \psi \sum_{t=1}^T G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}_*) - \mathbb{E} \left[ \sum_{t=1}^T G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{x}_t) \right] \right\} \leq A\sqrt{T}, \quad (166) $$

where

$$ A = \psi \sigma \sqrt{\frac{2D_{\max}}{\theta}} = \psi \frac{RL_{\max} \Delta_C}{s_{\min}} \sqrt{s_{\max} |\mathcal{V}| |\mathcal{M}|} \sqrt{2 \sum_{v \in \mathcal{V}} \min\{b^v, \|\mathbf{s}^v\|_1\} \log \left( \frac{\|\mathbf{s}^v\|_1}{\min\{b^v, \|\mathbf{s}^v\|_1\}} \right)}, $$

using the upper bounds on $\theta$, $\sigma$, and $D_{\max}$ determined in Lemmas F.2, F.3, and F.5, respectively.

This proves Theorem V.1. $\square$

<figure><img src="image_13.png" /></figure>

# APPENDIX H

## PROOF OF PROPOSITION V.1.1

*Proof.* Let $\bar{\mathbf{y}}$ be the average fractional allocation $\bar{\mathbf{y}} = \frac{1}{T} \sum_{t=1}^{\tilde{T}} \mathbf{y}_t$ of INFIDA, and $\bar{\mathbf{x}}$ the random state sampled from $\bar{\mathbf{y}}$ using DEPROUND. We take $G_T(\mathbf{y}) = \frac{1}{T} \sum_{t=1}^{\tilde{T}} G(\mathbf{r}_t, \mathbf{l}_t, \mathbf{y})$, $\forall \mathbf{y} \in \mathcal{Y}$. We have

$$ \mathbb{E}[G_T(\bar{\mathbf{x}})] \stackrel{(141)}{\geq} \mathbb{E} \left[ \frac{1}{T} \sum_{t=1}^{\tilde{T}} \Lambda(\mathbf{r}_t, \mathbf{l}_t, \bar{\mathbf{x}}) \right] \stackrel{(154)}{\geq} \frac{1}{T} \sum_{t=1}^{\tilde{T}} \Lambda(\mathbf{r}_t, \mathbf{l}_t, \bar{\mathbf{y}}) \stackrel{(141)}{\geq} \psi G_T(\bar{\mathbf{y}}). \quad (167) $$

Using Jensen's inequality we get

$$ \psi G_T(\bar{\mathbf{y}}) \geq \psi \frac{1}{\tilde{T}} \sum_{t=1}^{\tilde{T}} G_T(\mathbf{y}_t). \quad (168) $$

<figure><img src="image_14.png" /></figure>

It straightforward to check that $G_T$ satisfies the same properties 1 (concavity) and 3 (subgradient boundedness) as $G$ and the remaining properties are preserved under the same mirror map and convex decision set. With properties 1–5 satisfied, we can apply [69, Theorem 4.2] to obtain

$$ \sum_{t=1}^{\tilde{T}} G_T(\mathbf{y}_*) - \sum_{t=1}^{\tilde{T}} G_T(\mathbf{y}_t) = \tilde{T} G_T(\mathbf{y}_*) - \sum_{t=1}^{\tilde{T}} G_T(\mathbf{y}_t) \leq \sigma \sqrt{\frac{2D_{\max}}{\theta}} \tilde{T}. \quad (169) $$

Dividing both sides of the above inequality by $\tilde{T}$ gives

$$ \frac{1}{\tilde{T}} \sum_{t=1}^{\tilde{T}} G_T(\mathbf{y}_t) \geq G_T(\mathbf{y}_*) - \sigma \sqrt{\frac{2D_{\max}}{\theta \tilde{T}}}. \quad (170) $$

Using the same argument to obtain Eq. (163), i.e., restricting the maximization to the integral domain gives a lower value, we get

$$ \frac{1}{\tilde{T}} \sum_{t=1}^{\tilde{T}} G_T(\mathbf{y}_t) \geq G_T(\mathbf{x}_*) - \sigma \sqrt{\frac{2D_{\max}}{\theta \tilde{T}}}. \quad (171) $$

Using Eq. (167), and Eq. (171) we obtain

$$ \mathbb{E}[G_T(\bar{\mathbf{x}})] \geq \psi G_T(\mathbf{x}_*) - \psi \sigma \sqrt{\frac{2D}{\theta \tilde{T}}}. \quad (172) $$

Thus, $\forall \epsilon > 0$ and over a sufficiently large running time $\tilde{T}$ for INFIDA, $\bar{\mathbf{x}}$ satisfies

$$ \mathbb{E}[G_T(\bar{\mathbf{x}})] \geq \left(1 - \frac{1}{e} - \epsilon\right) G_T(\mathbf{x}_*). \quad (173) $$