# Towards Inference Delivery Networks:
Distributing Machine Learning
with Optimality Guarantees

Tareq Si Salem1, Gabriele Castellano12, Giovanni Neglia1, Fabio Pianese2,
Andrea Araldo3
1Inria, Université Côte d’Azur, France, {tareq.si-salem, gabriele.castellano,
<EMAIL_ADDRESS> 2Nokia Bell Labs, France, {fabio.pianese,
<EMAIL_ADDRESS> 3Télécom SudParis - Institut
Polytechnique de Paris, France<EMAIL_ADDRESS>
###### 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. (1)

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

  2. (2)

We frame the allocation of ML model in IDNs as an (NP-hard) optimization
problem that captures the trade-off between latency and accuracy (Sec. III).

  3. (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 (Sec. V).

  4. (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, 9, 10, 11, 12, 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, 18, 19, 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 auhtors 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. The
similarities between this problem and inference provisioning inspired us in
the design of Inference Delivery Networks. However, the two problems feature
some crucial differences. For instance, 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]; nonetheless the problem considered in this paper differs from the above,
as inference services can run under different configurations that lead to
variable resource consumption and inference accuracy. In some sense,
traditional services provide a binary answer to any user request: the request
can be satisfied if and only if the service is deployed at the node and any
two nodes deploying the service are equivalent. In ML inference, however,
several models can provide an answer but the accuracy of the answer can be
different [36].

A similar trade-off between resource usage and perceived quality typically
emerges in the context of video caching [37, 38, 39, 40, 41, 42, 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 [42, 41]. 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 optimize inference quality
and latency and provides strong performance guarantees without requiring any
knowledge about the request process thanks to our adversarial setting (Sec.
V).

Finally, ML model allocation in an IDN can also be considered as a particular
instance of _similarity caching_ [44], 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 [45],
contextual advertising [46], object recognition [47], and recommender systems
[48]. To the best of our knowledge, the literature on similarity caching has
restricted itself to (i) a single cache (with the exception of [48, 49, 50]),
and (ii) homogeneous items with identical resource 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_.111 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.  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 LABEL:tab:notations_definitions.

Figure 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.

###  III-A Compute Nodes and Models

(a) Keras pre-trained models (b) Pytorch pre-trained models

Figure 2: Example of pre-trained model catalog for the image classification task. Data from [51]. TABLE I: Notation Summary.  |  Inference Delivery Networks
---|---
$G(\mathcal{V},\mathcal{E})$ |  Undirected weighted graph, with nodes $\mathcal{V}$ and edges $\mathcal{E}$
$w_{u,v}$ |  Weight of edge $(u,v)\in\mathcal{E}$
$\mathcal{N}$ / $\mathcal{M}$ / $\mathcal{M}_{i}$ |  Tasks / models catalog / models catalog for task $i$
$s^{v}_{m}$ |  Size of model $m\in\mathcal{M}$ on node $v\in\mathcal{V}$
$a_{m}$ |  Prediction accuracy of model $m$
$d^{v}_{m}$ |  Average inference delay of model $m$ at node $v$
$b^{v}$ |  Allocation budget constraint at node $v$
$x^{v}_{m}$ |  0-1 indicator variable set to 1 if model m at node $v$ is allocated
$\omega^{v}_{m}$ |  0-1 indicator constant set to 1 if model m at node $v$ is a permanent repository model
$\boldsymbol{{\omega}}$ / $\boldsymbol{{x}}$ |  Minimal allocation vector (Sec. III-D) / allocation vector (Sec. III-A)
$\boldsymbol{{x}}^{v}$ / $\boldsymbol{{x}}$ |  Allocation vector at node $v$ / global allocation vector
$\boldsymbol{{s}}^{v}$ |  Vector of model sizes at node $v$
$\boldsymbol{{p}}$ |  Routing path $\\{p_{1},\dots,p_{J}\\}$ of connected nodes $p_{j}\in\mathcal{V}$
$\nu(\boldsymbol{{p}})$ |  Repository node associated to the path $\boldsymbol{{p}}$
$C^{p_{j}}_{\boldsymbol{{p}},m}$ |  Cost of serving on $p_{j}$ along path $\boldsymbol{{p}}$ using model $m$
$\rho$ |  Request type $\rho=(i,\boldsymbol{{p}})$, i.e., the requested task $i$ and the request’s routing path $\boldsymbol{{p}}$
$\mathcal{R}/\mathcal{R}_{i}$ |  Set of all the possible request types / request types for task $i$
$R$ |  Total number of request types
$r^{t}_{\rho}$ |  Number of times $\rho$ is requested during time slot $t$
$\mathrm{load}^{t,v}_{m}(\rho^{\prime})$ |  Number of type-$\rho$ requests served by model $m$ at node $v$ during the $t$-slot
$L^{v}_{m}$ |  Maximum capacity of model $m$ on node $v$
$l^{t,v}_{\rho,m}$ |  Potential available capacity of $m$ on node $v$ for request type $\rho$ at time $t$
$\gamma^{k}_{\rho}$ |  $k$-th smallest cost for request type $\rho$ along its path (Sec. III-E).
$\boldsymbol{{r}}_{t}$ / $\boldsymbol{{l}}_{t}$ |  Request batch vector / potential available capacities vector
$\lambda^{k}_{\rho}{(\boldsymbol{{l}}_{t})}$ / $z^{k}_{\rho}{(\boldsymbol{{l}}_{t},\boldsymbol{{x}})}$ |  Potential / effective available capacity of the model serving $\rho$ with cost $\gamma^{k}_{\rho}$ (Sec. III-E).
$Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}})$ |  Number of requests of type $\rho$ that can be served by the $k$-th smallest cost models found along path $\boldsymbol{{p}}$ under allocation $\boldsymbol{{x}}$ at time slot $t$ (15)
$\kappa_{\rho}(v,m)$ |  Rank of model $m$ allocated to node $v$, among all the models that can potentially serve requests of type $\rho$
$K_{\rho}$ |  Maximum number of models that request $\rho$ may encounter along its serving path
$\mathcal{X}^{v}$ / $\mathcal{X}$ |  Set of valid integral allocations at node $v$ / at all nodes
$C$ / $G$ |  The overall system cost (12) / allocation gain (13)
|  INFIDA
$\Phi^{v}$ |  Weighted negative entropy map $\Phi^{v}:\mathcal{D}^{v}\to\mathbb{R}$ given as $\Phi^{v}(\boldsymbol{{y}}^{v})=\sum_{m\in\mathcal{M}}s^{v}_{m}y^{v}_{m}\log(y^{v}_{m})$
$\Phi$ |  Global mirror map $\Phi(\boldsymbol{{y}})=\sum_{v\in\mathcal{V}}\Phi^{v}(\boldsymbol{{y}}^{v})$
$y^{v}_{m}$ |  Fractional allocation (allocation probability) of model $m$ is at node $v$
$\boldsymbol{{y}}^{v}$ / $\boldsymbol{{y}}$ |  Fractional allocation vector at node $v$ / global fractional allocation vector
$\mathcal{Y}^{v}$ / $\mathcal{Y}$ |  Set of valid fractional allocations at node $v$ / at all nodes
$\vec{g}_{t}$ |  Subgradient vector of $G$ over $\mathcal{Y}$ at point $\boldsymbol{{y}}_{t}$
$g^{v}_{t,m}$ |  Component $(v,m)$ of the subgradient vector $\boldsymbol{{g}}_{t}$
$\mathcal{P}_{\mathcal{Y}^{v}\cap\mathcal{D}^{v}}^{\Phi^{v}}(\,\cdot\,)$ |  Projection operator onto $\mathcal{Y}^{v}\cap\mathcal{D}^{v}$
$\eta$ |  learning rate
$B$ |  Refresh period
$T$ |  Time horizon
$t$ |  Time slot $t\in[T]$
$A$ |  Regret constant
$L_{\max}$ |  Upper bound on model capacities
$\Delta_{C}$ |  Upper bound on maximum serving cost difference
$\psi$ |  Regret discount factor equal to $1-\frac{1}{e}$

We represent the inference delivery network (IDN) as a weighted graph
$G(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ is the set of compute nodes,
and $\mathcal{E}$ represents their interconnections. Each node
$v\in\mathcal{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,\dots,|\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
$\mathcal{M}_{i}$. Each task is served by a separate set of models, i.e.,
$\mathcal{M}_{i}{\cap}\mathcal{M}_{i^{\prime}}{=}\emptyset,\forall
i,i^{\prime}\in\mathcal{N},i\neq i^{\prime}$. Catalog $\mathcal{M}_{i}$ may
encompass, for instance, independently trained models or shrunk versions of a
high quality model generated through distillation [52, 53]. We denote by
$\mathcal{M}{=}\cup_{i\in\mathcal{N}}\mathcal{M}_{i}=\\{1,2,\dots,|\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\mathcal{V}$, we denote by

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

the decision variable that indicates if model $m\in\mathcal{M}$ is deployed on
node $v$.222 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.  Therefore,
$\boldsymbol{{x}}^{v}=[x^{v}_{m}]_{m\in\mathcal{M}}$ is the allocation vector
on node $v$, and $\boldsymbol{{x}}=[\boldsymbol{{x}}^{v}]_{v\in\mathcal{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\mathcal{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$.333Note 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).
Therefore, budget constraints are expressed as

$\sum_{m\in\mathcal{M}}x^{v}_{m}s_{m}^{v}\leq b^{v},\forall v\in V.$ (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
$\boldsymbol{{\omega}}=[\omega^{v}_{m}]_{(v,m)\in\mathcal{V}\times\mathcal{M}}$
the _minimal allocation_. Note that the presence of repositories introduce the
following constraints to the allocation vector:

$\displaystyle x^{v}_{m}\geq\omega^{v}_{m},\forall v\in\mathcal{V},\forall
m\in\mathcal{M}.$ (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.,

$\displaystyle\mathcal{X}^{v}\triangleq\left\\{\boldsymbol{{x}}^{v}\in\\{0,1\\}^{\mathcal{M}}:\boldsymbol{{x}}^{v}\text{
satisfies Eqs.~{}\eqref{eq:decision-
variable}--\eqref{eq:repo_ineq}}\right\\}.$ (4)

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

###  III-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}$. The routing is therefore
predetermined, and our decisions only concern placement of models (i.e.,
variables $x^{v}_{m}$).

A routing path $\boldsymbol{{p}}$ of length $|\boldsymbol{{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 [54], we assume that paths are simple,
i.e., they do not contain repeated nodes. A request is therefore characterized
by the pair $(i,\boldsymbol{{p}})$, where $i$ is the task requested and
$\boldsymbol{{p}}$ is the routing path to be traversed. We call the pair
$(i,\boldsymbol{{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(\boldsymbol{{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.

###  III-C Cost Model

When serving a request of type $\rho{=}(i,\boldsymbol{{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. Additionally, the system experiences a network cost, due to
using the path between $p_{1}$ and $p_{j}$. Similarly to previous work [55],
we can write the total cost of serving a request as

$C^{p_{j}}_{\boldsymbol{{p}},m}=f((p_{1},\dots,p_{j}),m).$ (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^{p_{j}}_{\boldsymbol{{p}},m}=\sum_{j^{\prime}=1}^{j-1}w_{p_{j^{\prime}},p_{j^{\prime}+1}}+d^{p_{j}}_{m}+\alpha(1{-}a_{m}),$
(6)

where $a_{m}$ and $d^{p_{j}}_{m}$ 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 [56]. Parameter
$w_{v,v^{\prime}}\in\mathbb{R}_{+}$ is the (round-trip) latency of edge
$(v,v^{\prime})\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^{p_{j}}_{m}+\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).

###  III-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
$\boldsymbol{{r}}_{t}=[r^{t}_{\rho}]_{\rho\in\mathcal{R}}$, where
$r^{t}_{\rho}\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^{v}_{m}\in\mathbb{N}$ when
deployed at node $v\in\mathcal{V}$, i.e., it can serve at most $L^{v}_{m}$
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^{t,v}_{\rho,m}\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 $\boldsymbol{{r}}_{t}$ and allocation vector $\boldsymbol{{x}}_{t}^{v}$.
Formally, let $\mathrm{load}^{t,v}_{m}(\rho)$ denote the number of type-$\rho$
requests served by model $m$ at node $v$ during the $t$-slot, then

$\displaystyle
l^{t,v}_{\rho,m}\triangleq\min\left\\{L^{v}_{m}-\sum_{\rho^{\prime}\in\mathcal{R}\setminus\\{\rho\\}}\mathrm{load}^{t,v}_{m}(\rho^{\prime}),\,r^{t}_{\rho}\right\\}.$
(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,\boldsymbol{{p}})$ and
$\rho^{\prime}=(i,\boldsymbol{{p}}^{\prime})$ arrive at node $v$. The arrival
order and the node scheduling discipline may determine that many requests of
type $\rho^{\prime}$ be served, which would leave a small $l^{t,v}_{\rho,m}$
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^{t,v}_{\rho,m}\triangleq\min\\{L^{v}_{m},r^{t}_{\rho}\\}$. The _effective
available capacity_ is then equal to $l^{t,v}_{\rho,m}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 guarantees in the adversarial setting. In what follows, we can then
consider that values $l^{t,v}_{\rho,m}$ are exogeneously determined. The
vector of potential available capacities at time $t\in[T]$ is denoted by

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

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

$\displaystyle\sum_{\rho\in\mathcal{R}_{i}}r^{t}_{\rho}\leq\sum_{m\in\mathcal{M}_{i}}\omega^{\nu(\boldsymbol{{p}})}_{m}L^{\nu(\boldsymbol{{p}})}_{m}\ignorespaces,\forall
i\in\mathcal{N}.$ (9)

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

$\displaystyle\mathcal{A}\triangleq\bigg{\\{}$
$\displaystyle(\boldsymbol{{r}},\boldsymbol{{l}})\in(\mathbb{N}\cup\\{0\\})^{\mathcal{R}}\times(\mathbb{N}\cup\\{0\\})^{\bigcup_{i\in\mathcal{N}}\mathcal{R}_{i}\times\mathcal{M}_{i}\times\mathcal{V}}:$
$\displaystyle\sum_{\rho\in\mathcal{R}_{i}}r^{t}_{\rho}\leq\sum_{m\in\mathcal{M}_{i}}\omega^{\nu(\boldsymbol{{p}})}_{m}L^{\nu(\boldsymbol{{p}})}_{m},l^{v}_{\rho,m}\leq\min\\{L^{v}_{m},r_{\rho}\\},$
$\displaystyle\forall
i\in\mathcal{N},v\in\mathcal{V},m\in\mathcal{M},\rho\in\mathcal{R}\bigg{\\}}.$
(10)

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

###  III-E Serving Model

Figure 3: Necessity of partial synchronization in IDN among close-by computing
nodes under the cost model in Eq. (6).

Given request $\rho{=}(i,\boldsymbol{{p}}){\in}\mathcal{R}$, let
$K_{\rho}=|\boldsymbol{{p}}||\mathcal{M}_{i}|$ denote the maximum number of
models that may encounter along its serving path $\boldsymbol{{p}}$. We order
the corresponding costs $\\{C^{p_{j}}_{\boldsymbol{{p}},m},\forall
m\in\mathcal{M}_{i},\forall p_{j}\in\boldsymbol{{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.444 Note that we do not
consider only the models deployed in the network, but all the possible node-
model pairs.  If $v\notin\boldsymbol{{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^{k}_{\rho}$,
$\lambda^{k}_{\rho}(\boldsymbol{{l}}_{t})$, and
$z^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})$, respectively:

$\displaystyle\gamma^{k}_{\rho}=C^{v}_{\boldsymbol{{p}},m},\quad\lambda^{k}_{\rho}(\boldsymbol{{l}}_{t})=l^{t,v}_{\rho,m},\quad
z^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})=x^{v}_{m}l^{t,v}_{\rho,m}.$
(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,\boldsymbol{{p}})$, it only needs information about those models that (i)
are located upstream on the serving path (i.e., on nodes
$p_{l}\in\boldsymbol{{p}}$ with $l>h$), and (ii) provide a cost smaller than
$\gamma^{k}_{\rho}$. Since the cost increases with the network latency (as
illustrated in Fig. 3), the number of models satisfying these criteria is
small in practice.555 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$).  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^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})<r^{t}_{\rho}$).
If this is the case, model $m$ will serve with cost $\gamma_{\rho}^{k}$ at
most $z^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})$ requests (its
effective available capacity) out of the
$r^{t}_{\rho}-\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})$
requests still to be satisfied. The aggregate cost incurred by the system at
time slot $t$ is then given by

$\displaystyle
C(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}})\;{=}\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}}\gamma^{k}_{\rho}$
$\displaystyle\cdot\min\biggl{\\{}r^{t}_{\rho}\;{-}\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}}),\hskip
1.00006ptz^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\biggr{\\}}$
$\displaystyle\cdot\mathds{1}_{\left\\{\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})<r^{t}_{\rho}\right\\}}.$
(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^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})$. We add the
indicator function $\mathds{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.

###  III-F Allocation Gain and Static Optimal Allocations

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

$G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}})=C(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})-C(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}).$
(13)

The first term
$C(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})$ on the
right hand side is the service cost when only repository models are present in
the network. Since intermediate nodes can help serving the requests at a
reduced cost,
$C(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})$ is an
upper bound on the aggregate serving cost, and the allocation gain captures
the cost reduction achieved by model allocation $\boldsymbol{{x}}$.

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

$\displaystyle\boldsymbol{{x}}_{*}=\underset{\boldsymbol{{x}}\in\mathcal{X}}{\operatorname*{\arg\max}}\left(G_{T}(\boldsymbol{{x}})\triangleq\frac{1}{T}\sum^{T}_{t=1}G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right).$
(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 . It is
known that submodular maximization problems are NP-hard and cannot be
approximated with a ratio better than $(1-1/e)$ even under simpler cardinality
constraints [57]. 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 [58].

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

$\displaystyle
Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}})\triangleq\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}{z}^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}.$
(15)

The $\min\\{\,\cdot\,,\,\cdot\,\\}$ operator denotes that we can never serve
more than the number of requests $r^{t}_{\rho}$ issued by users. Observe that,
being the minimal allocation $\boldsymbol{{\omega}}$ an input parameter not
dependent on our decisions,
$Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})$
is a constant. Additionally, since the models allocated in $\boldsymbol{{x}}$
always include those allocated in $\boldsymbol{{\omega}}$, we have
$Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}})\geq
Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\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:

$\displaystyle G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}})=$
$\displaystyle\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\underset{\text{cost
saving}}{\underbrace{\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)}}$
$\displaystyle\underset{\text{additional
requests}}{\underbrace{\left(Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}})-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})\right).}}$
(16)

We prove this lemma in Appendix B. This result tells us that the gain of a
certain allocation $\boldsymbol{{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 the gain with a component (i)
proportional to its cost saving $\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}$ 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 $\boldsymbol{{x}}$ can serve with respect to the minimal allocation
$\boldsymbol{{\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.

###  IV-A Algorithm Overview

1:procedure
INFIDA($\boldsymbol{{y}}^{v}_{1}{=}\underset{\boldsymbol{{y}}^{v}\in\mathcal{Y}^{v}\cap\mathcal{D}^{v}}{\arg\min}\,\Phi^{v}(\boldsymbol{{y}}^{v})$,
$\boldsymbol{{x}}^{v}_{1}{=}\textsc{DepRound}(\boldsymbol{{y}}^{v}_{1})$,
$\eta{\in}\mathbb{R}_{+}$)

2: for $t=1,2,\dots,T$ do

3: Compute
$\boldsymbol{{g}}^{v}_{t}\in\partial_{\boldsymbol{{y}}^{v}}G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})$
through (18).

4:
$\hat{\boldsymbol{{y}}}^{v}_{t}\leftarrow\nabla\Phi^{v}(\boldsymbol{{y}}^{v}_{t})$$\triangleright$
Map state to the dual space

5:
$\hat{\boldsymbol{{h}}}^{v}_{t+1}\leftarrow\hat{\boldsymbol{{y}}}^{v}_{t}+\eta\boldsymbol{{g}}^{v}_{t}$$\triangleright$
Take gradient step in the dual space

6:
$\boldsymbol{{h}}^{v}_{t+1}\leftarrow\left(\nabla\Phi^{v}\right)^{-1}(\hat{\boldsymbol{{h}}}^{v}_{t+1})$$\triangleright$
Map dual state back to the primal space

7:
$\boldsymbol{{y}}^{v}_{t+1}\leftarrow\mathcal{P}_{\mathcal{Y}^{v}\cap\mathcal{D}^{v}}^{\Phi^{v}}(\boldsymbol{{h}}^{v}_{t+1})$$\triangleright$
Project new state onto the feasible region using Algorithm 2

8:
$\boldsymbol{{x}}^{v}_{t+1}\leftarrow\textsc{DepRound}(\boldsymbol{{y}}^{v}_{t+1})$
$\triangleright$ Sample a discrete allocation

Algorithm 1 INFIDA distributed allocation on node $v$

On every node $v\in\mathcal{V}$, INFIDA updates the allocation
$\boldsymbol{{x}}^{v}\in\mathcal{X}^{v}{\subset}\\{0,1\\}^{|\mathcal{M}|}$, by
operating on a correspondent fractional state
$\boldsymbol{{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 $\norm{\boldsymbol{{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
$\norm{\boldsymbol{{s}}^{v}}_{1}\geq b^{v}$ then

$\displaystyle\mathcal{Y}^{v}\triangleq\left\\{\boldsymbol{{y}}^{v}\in[0,1]^{\mathcal{M}}:\sum_{m\in\mathcal{M}}y^{v}_{m}s_{m}^{v}=b^{v},\forall
v\in V\right\\};$ (17)

otherwise, for the corner case $\norm{\boldsymbol{{s}}^{v}}_{1}<b^{v}$, we
have $\mathcal{Y}^{v}\triangleq\left\\{[1]^{\mathcal{M}}\right\\}$.

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

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

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 $\boldsymbol{{y}}^{v}$, and (ii) updates
its local allocation $\boldsymbol{{x}}^{v}$ via randomized rounding (Sec.
IV-C). INFIDA is summarized in Algorithm 1 and detailed below.

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

In our work, we use a generalized version of the gradient method called Online
Mirror Ascent (OMA) [59, Ch. 4]. OMA uses a function
$\Phi^{v}:\mathcal{D}^{v}\to\mathbb{R}_{+}$ (mirror map) to map
$\boldsymbol{{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}(\boldsymbol{{y}}^{v})=\sum_{m\in\mathcal{M}}s^{v}_{m}y^{v}_{m}\log(y^{v}_{m})$,
which is known to achieve better convergence rate in high dimensional spaces
when each subgradient component is bounded.666Technically, 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 [59, Sec. 4.3]. To
compute a feasible fractional state $\boldsymbol{{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 [60] to obtain a negative entropy
projection
$\mathcal{P}^{\Phi^{v}}_{\mathcal{Y}^{v}\cap\mathcal{D}^{v}}(\,\cdot\,)$. Our
adaptation is described in Appendix C.

Allocation update. Once the fractional state $\boldsymbol{{y}}^{v}$ has been
updated, the final step of INFIDA is to determine a new random discrete
allocation $\boldsymbol{{x}}^{v}$ and update the local models accordingly. The
sampled allocation $\boldsymbol{{x}}^{v}$ should (i) comply with the budget
constraint (2) on node $v$ and (ii) be consistent with the fractional state,
i.e., $\mathds{E}[x^{v}_{m}]=y^{v}_{m}\;\forall m\in\mathcal{M}$. To this
purpose, we use the DepRound [61] 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.

###  IV-B Subgradient Computation

At the end of every time slot $t$, a subgradient $\boldsymbol{{g}}_{t}$ of the
gain function in Eq. (16) at point $\boldsymbol{{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 D):

$\displaystyle g^{v}_{t,m}=\sum_{\rho\in\mathcal{R}}l^{t,v}_{\rho,m}$
$\displaystyle\cdot\left(\gamma^{K^{*}_{\rho}(\boldsymbol{{y}}_{t})}_{\rho}-C^{v}_{\boldsymbol{{p}},m}\right)\cdot\mathds{1}_{\\{\kappa_{\rho}(v,m)<K^{*}_{\rho}(\boldsymbol{{y}}_{t})\\}},$
(18)

where $K^{*}_{\rho}(\boldsymbol{{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^{t}_{\rho}$ requests in the batch given the fractional state
$\boldsymbol{{y}}_{t}$. Formally,
$K^{*}_{\rho}(\boldsymbol{{y}}_{t})\triangleq\min\big{\\{}k\in[K_{\rho}-1]:\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})\geq
r^{t}_{\rho}\big{\\}}$.

For the sake of clarity assume that the mirror map is Euclidean, and then the
dual and primal spaces collapse and
$\hat{\boldsymbol{{y}}}_{t}=\boldsymbol{{y}}_{t}$. At each iteration, each
component $y^{t,v}_{m}$ 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^{v}_{m,t}$ 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^{v}_{m,t}$. This contribution is
proportional to the potential available capacity of model $m$ on node $v$ and
to the relative gain
$\left(\gamma^{K^{*}_{\rho}(\boldsymbol{{y}}_{t})}_{\rho}-C^{v}_{\boldsymbol{{p}},m}\right)$,
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. 1.

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

  2. 2.

As the control message travels upstream, intermediate nodes add to $Z$ the
local values $z^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{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^{t}_{\rho}$, that is until the message reaches the
$K^{*}_{\rho}(\boldsymbol{{y}}_{t})$-th model.

  3. 3.

Once the $K^{*}_{\rho}(\boldsymbol{{y}}_{t})$-th model is detected, a control
message is sent down in the opposite direction, containing the cost
$\gamma^{K^{*}_{\rho}(\boldsymbol{{y}}_{t})}_{\rho}$ 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

$\displaystyle
h^{v}_{m}=l^{t,v}_{\rho,m}\cdot\left(\gamma^{K^{*}_{\rho}(\boldsymbol{{y}}_{t})}_{\rho}-C^{v}_{\boldsymbol{{p}},m}\right).$
(19)

  4. 4.

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

$\displaystyle g^{v}_{t,m}=\sum_{m\in\mathcal{M}_{i}}h^{v}_{m}.$

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^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t}),\gamma^{k}_{\rho})$
to the message, and lets upstream nodes to apply any pending update in the
correct order.

###  IV-C State Rounding

Once the new fractional state $\boldsymbol{{y}}_{t+1}$ is computed, each node
$v$ independently draws a random set of models to store locally in such a way
that $\mathds{E}[\boldsymbol{{x}}^{v}_{t+1}]=\boldsymbol{{y}}^{v}_{t+1}$. This
sampling guarantees that the final allocation $\boldsymbol{{x}}^{v}_{t+1}$
satisfies constraint (2) in expectation. A naive approach is to draw each
variable $x^{v,t+1}_{m}$ independently, but it leads to a large variance of
the 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 [61].
The procedure modifies the fractional state $\boldsymbol{{y}}^{v}_{t+1}$
iteratively: at each iteration, DepRound operates on two fractional variables
$y^{v,t+1}_{m},y^{v,t+1}_{m^{\prime}}$ so that at least one of them becomes
integral and the aggregate size of the corresponding models
$s^{v}_{m}y^{v,t+1}_{m}+s^{v}_{m^{\prime}}y^{v,t+1}_{m^{\prime}}$ 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 $\mathcal{O}(|\mathcal{M}|)$ steps.

Note that, to satisfy
$\mathds{E}[\boldsymbol{{x}}^{v}_{t+1}]=\boldsymbol{{y}}^{v}_{t+1}$, the
residual fractional variable, say it $y^{v,t+1}_{\bar{m}}$, needs to be
rounded. At this point $x^{v,t+1}_{\bar{m}}$ 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 [62]. In our scenario, the $\psi$-regret is defined as the gain
loss in comparison to the best static allocation in hindsight, i.e.,
$\boldsymbol{{x}}_{*}\in\operatorname*{\arg\max}_{\boldsymbol{{x}}\in\mathcal{X}}\sum^{T}_{t=1}G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}})$,
discounted by a factor $\psi\in(0,1]$. Formally,

$\displaystyle\psi\text{-}\mathrm{Regret}_{T,\mathcal{X}}\triangleq$ (20)
$\displaystyle\underset{{\scriptscriptstyle{\\{\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t}\\}_{t=1}^{T}\in\mathcal{A}^{T}}}}{\sup}\left\\{\psi\sum^{T}_{t=1}G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}_{*}){-}\mathbb{E}\left[\sum^{T}_{t=1}G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}_{t})\right]\right\\},$

where allocations $\boldsymbol{{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 $\boldsymbol{{r}}_{t}$ and
$\boldsymbol{{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 $\boldsymbol{{r}}_{t}$) or difficult to model system interactions
(e.g., the available capacities $\boldsymbol{{l}}_{t}$). This technique has
been recently successfully used to model caching problems (e.g., in [63, 60]).
Our main result is the following (the full proof is in Appendix F):

######  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:

$\displaystyle\text{$\left(1-1/e\right)$-}\mathrm{Regret}_{T,\mathcal{X}}\leq
A\sqrt{T},$ (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 $\boldsymbol{{x}}_{t}$ is a $(1{-}1/e)$-approximation of the
fractional gain. Then, we use online learning results [59] 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 F.

∎

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

######  Proposition V.1.1.

(offline solution) Replace in INFIDA the allocation gain
$G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}})$ by
$G_{T}(\boldsymbol{{y}})$ (defined in (14)). After $\tilde{T}$ iterations, let
$\bar{\boldsymbol{{y}}}$ be the average fractional allocation
$\bar{\boldsymbol{{y}}}=\frac{1}{\tilde{T}}\sum^{\tilde{T}}_{\tilde{t}=1}\boldsymbol{{y}}_{t}$,
and $\bar{\boldsymbol{{x}}}$ the random state sampled from
$\bar{\boldsymbol{{y}}}$ using DepRound. $\forall\epsilon>0$, for $\bar{T}$
large enough, $\bar{\boldsymbol{{x}}}$ satisfies

$\displaystyle\mathbb{E}\left[G_{T}(\bar{\boldsymbol{{x}}})\right]\geq\left(1-\frac{1}{e}-\epsilon\right)G_{T}({\boldsymbol{{x}}}_{*}),$
(22)

where
${\boldsymbol{{x}}}_{*}=\operatorname*{\arg\max}_{\boldsymbol{{x}}\in\mathcal{X}}G_{T}(\boldsymbol{{x}})$.

The proof is given in Appendix. G.

##  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 [62].

Topology. We simulate a hierarchical topology similar to [64] 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 [65, 66], profiled for each simulated
processing unit, for a total of 10 models (Table II). We consider a task
catalog with $|\mathcal{N}|=20$ different object detection tasks. We allow $3$
duplicates per model; this gives $|\mathcal{M}_{i}|=30$ alternative models per
task $i\in\mathcal{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. 4(a)), a request is for task $i$ with
constant probability
$p(i)=\frac{(i+1)^{-1.2}}{\sum_{i^{\prime}\in\mathcal{N}}(i^{\prime}+1)^{-1.2}}$
(a Zipf distribution with exponent $1.2$), while _(ii)_ in the _Sliding
Popularity Profile_ (Fig. 4(b)), the $l$-th consecutive request is for task
$i$ with probability $\tilde{p}(i,l)=p\left(\left(i+5\lfloor
l/W\rfloor\right)\mod 20\right)$, 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}$).

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

(a) _Fixed Popularity Profile_ (b) _Sliding Popularity Profile_

Figure 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.

Static greedy. We adapt the static greedy (SG) heuristic from the cost-benefit
greedy in [62]. 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 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^{v}_{m,\rho}$ 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\mathcal{M}$ with positive gain that we
denote by $q^{v}_{m,\rho}$. For every model $m$, an importance weight is
computed as
$w^{v}_{m}{=}\frac{1}{s^{v}_{m}}\frac{1}{|\mathcal{R}|}\sum_{\rho\in\mathcal{R}}q^{v}_{m,\rho}\min\\{\phi^{v}_{m,\rho},L^{v}_{m}\\}$,
where $s^{v}_{m}$ is the size of model $m$ and
$\min\\{\phi^{v}_{m,\rho},L^{v}_{m}\\}$ 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^{v}_{m,\rho},L^{v}_{m}\\}$ from $\phi^{v}_{m_{*},\rho}$ and from
all the $\phi^{v}_{m^{\prime},\rho}\colon
q^{v}_{m^{\prime},\rho}<q^{v}_{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 $\textsc{INFIDA}_{\textsc{Offline}}$, 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 $\mathcal{P}$ with the
associated sequence of allocation decisions
$\\{\boldsymbol{{x}}_{t}\\}^{T}_{t=1}$ is evaluated in terms of the time-
averaged gain normalized to the number of requests per time slot (NTAG):

$\displaystyle\mathrm{NTAG}(\mathcal{P})=\sum^{T}_{t=1}\frac{1}{T\norm{\boldsymbol{{r}}_{t}}_{1}}G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}_{t}).$
(23)

Moreover, we evaluate the update cost of a policy $\mathcal{P}$ with the
associated sequence of allocation decisions
$\\{\boldsymbol{{x}}_{t}\\}^{T}_{t=1}$ 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:

$\displaystyle\mathrm{MU}(\mathcal{P})\triangleq\frac{1}{T}\sum^{T}_{t=2}\sum_{(v,m)\in\mathcal{V}\times\mathcal{M}}s^{v}_{m}\max\\{0,x^{v}_{t,m}-x^{v}_{t-1,m}\\}.$
(24)

###  VI-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$.777The inaccuracy cost is taken in 0–100, then $\alpha$ picked here
corresponds to $100\times$ scaling of the parameter defined in Eq. (6).

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. 5(a)), 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. 5(b)). 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. 5(c)). Since higher quality models feature, in
general, a lower serving capacity (Table II), Fig. 5(c) suggests that a
significant number of requests is served in the cloud (Tier 0) for this value
of $\alpha$.

(a) $\alpha=3$ (b) $\alpha=4$ (c) $\alpha=5$

Figure 5: Fractional allocation decisions $y^{v}_{m}$ 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.  Figure 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_.

(a) *

(b) (c)

Figure 7: NTAG of the different policies under _Sliding Popularity Profile_
and network topologies: (a) _Network Topology I_ , and (b) _Network Topology
II_.

(a) *

(b) Models Updates (MU) (c) NTAG

Figure 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_.

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 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 $\textsc{INFIDA}_{\textsc{Offline}}$ for different values of
$\alpha$ under the _Sliding Popularity Profile_. Results are shown both for
_Network Topology I_ (Fig. 7(b)) and for _Network Topology II_ (Fig. 7(c)).

The plot shows that the gain decreases by increasing $\alpha$. This is
expected since the gain (13) is defined as the improvement 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
$\textsc{INFIDA}_{\textsc{Offline}}$ 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
$\textsc{INFIDA}_{\textsc{Offline}}$ has better performance than SG: although
both policies are offline, $\textsc{INFIDA}_{\textsc{Offline}}$ manages to
provide a better allocation.

###  VI-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 stability of the system. Indeed,
frequent updates 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 = 1H. We run this
experiment under _Network Topology I_ and _Sliding Popularity Profile_. We set
the trade-off parameter $\alpha=1$.

In particular, Figure 8(b) shows the update cost (MU) for different refresh
periods, while Figure 8(c) 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. 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. 8(c): 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}}$.

###  VI-C Scalability on Requests Load

(a) *

(b) _Fixed Popularity Profile_ (c) _Sliding Popularity Profile_

Figure 9: NTAG of the different policies for different request rates under
_Network Topology I_.

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 $\textsc{INFIDA}_{\textsc{Offline}}$ 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 $\textsc{INFIDA}_{\textsc{Offline}}$ ($\approx$8%) is
considerably limited compared to SG ($\approx$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 $\approx$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. 9(c)). 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
$\textsc{INFIDA}_{\textsc{Offline}}$ consistently provide the lowest average
inaccuracy and latency under both high request load (10,000 rps) and default
request load (7,500 rps). $\textsc{INFIDA}_{\textsc{Offline}}$ 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. 9(b)).

Figure 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_.

##  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. G. Neglia acknowledges support from 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: http://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] G. Goel, Y. Lin, H. Sun _et al._ , “Beyond Online Balanced Descent: An Optimal Algorithm for Smoothed Online Optimization,” _Advances in Neural Information Processing Systems_ , 2019.
  * [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,” https://cloud.google.com/vertex-ai.
  * [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 provisioning 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 Architecture 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 Accelerator for Low-Power Deep Learning Inference on Mobile Devices,” in _2016 15th ACM/IEEE International Conference on Information Processing 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: Processing 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] M. Garetto, E. Leonardi, and G. Neglia, “Similarity Caching: Theory and Algorithms,” in _IEEE INFOCOM 2020-IEEE Conference on Computer Communications_ , 2020.
  * [45] 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.
  * [46] 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.
  * [47] 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.
  * [48] 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.
  * [49] 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.
  * [50] M. Garetto, E. Leonardi, and G. Neglia, “Content Placement in Networks of Similarity Caches,” _arXiv preprint arXiv:2102.04974_ , 2021.
  * [51] D. Blalock, J. J. G. Ortiz, J. Frankle _et al._ , “What is the State of Neural Network Pruning?” _arXiv preprint arXiv:2003.03033_ , 2020.
  * [52] G. Hinton _et al._ , “Distilling the Knowledge in a Neural Network,” _preprint arXiv:1503.02531_ , 2015.
  * [53] S. Ravi, “Custom On-Device ML Models with Learn2Compress,” 2018.
  * [54] S. Ioannidis and E. Yeh, “Adaptive caching networks with optimality guarantees,” _ACM SIGMETRICS Performance Evaluation Review_ , vol. 44, no. 1, pp. 113–124, 2016.
  * [55] 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.
  * [56] Y. Wang, G. Y. Wei, and D. Brooks, “Benchmarking TPU, GPU, and CPU platforms for deep learning,” _arXiv preprint arXiv:1907.10701_ , 2019.
  * [57] G. L. Nemhauser and L. A. Wolsey, “Best Algorithms for Approximating the Maximum of a Submodular Set Function,” _Mathematics of operations research_ , 1978.
  * [58] Y. Fairstein, A. Kulik, J. S. Naor _et al._ , “A (1-1/e-$\varepsilon$)-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.
  * [59] S. Bubeck, “Convex Optimization: Algorithms and Complexity,” _Foundations and Trends® in Machine Learning_ , vol. 8, no. 3-4, pp. 231–357, Nov. 2015.
  * [60] T. Si Salem, G. Neglia, and S. Ioannidis, “No-Regret Caching via Online Mirror Descent,” in _IEEE International Conference on Communications (ICC)_ , 2021.
  * [61] 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.
  * [62] A. Krause and D. Golovin, “Submodular Function Maximization,” _Tractability_ , vol. 3, pp. 71–104, 2014.
  * [63] G. S. Paschos, A. Destounis, L. Vigneri _et al._ , “Learning to Cache With No Regrets,” in _IEEE INFOCOM_ , 2019.
  * [64] 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.
  * [65] A. Bochkovskiy, C.-Y. Wang, and H.-Y. M. Liao, “YOLOv4: Optimal Speed and Accuracy of Object Detection,” _preprint arXiv:2004.10934_ , 2020.
  * [66] Y. Cai _et al._ , “YOLObile: Real-Time Object Detection on Mobile Devices via Compression-Compilation Co-Design,” _preprint arXiv:2009.05697_ , 2020\.
  * [67] B. S. Mordukhovich and N. M. Nam, “Geometric Approach to Convex Subdifferential Calculus,” _Optimization_ , vol. 66, no. 6, pp. 839–873, 2017.
  * [68] R. T. Rockafellar, _Convex Analysis_. Princeton university press, 2015.
  * [69] S. Shalev-Shwartz, “Online Learning: Theory, Algorithms, and Applications,” Ph.D. dissertation, The Hebrew University of Jerusalem, 2007.
  * [70] 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.
  * [71] 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 $\boldsymbol{{x}}(S)$ with $x_{m}^{v}(S)=1$ if
$((v,m)\in S\lor\omega^{v}_{m}=1)$, $x^{v}_{m}(S)=0$ otherwise. Now we can
define the set function $f_{t}\colon
2^{\mathcal{V}\times\mathcal{M}}\to\mathbb{R}$:

$\displaystyle f_{t}(S)\triangleq
G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}(S)).$ (25)

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

$\displaystyle
F_{T}(S)\triangleq\frac{1}{T}\sum^{T}_{t=1}f_{t}(S)=\frac{1}{T}\sum^{T}_{t=1}G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}(S)).$
(26)

We first start proving that $f_{t}$ is submodular in the following lemma.

######  Lemma A.1.

The set function $f_{t}\colon 2^{\mathcal{V}\times\mathcal{M}}\to\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

$\displaystyle
f_{t}(\emptyset)=G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}(\emptyset))=G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})=C(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})-C(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})=0.$
(27)

Submodularity. A function $f_{t}\colon
2^{\mathcal{V}\times\mathcal{M}}\to\mathbb{R}$ is submodular [62] if for every
$S^{\prime}\subset S^{\prime\prime}\subset\mathcal{V}\times\mathcal{M}$ and
$(\bar{v},\bar{m})\in(\mathcal{V}\times\mathcal{M})\setminus S^{\prime\prime}$
it holds that

$f_{t}(S^{\prime\prime}\cup\\{(\bar{v},\bar{m})\\})-f_{t}(S^{\prime\prime})\leq
f_{t}(S^{\prime}\cup\\{(\bar{v},\bar{m})\\})-f_{t}(S^{\prime}).$

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

Since $S^{\prime}\subset S^{\prime\prime}$, we have that if
${x^{\prime}}^{v}_{m}=1\implies{x^{\prime\prime}}^{v}_{m}=1$, and thus
$z^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime})\leq
z^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})$, for any
$k$ (see Eq. (11)). Therefore

$\displaystyle
Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime})\leq
Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime}),\forall
k\in[K_{\rho}-1],\forall\rho\in\mathcal{R}.$ (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:

$\displaystyle\min\\{r_{\rho}^{t}-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime}),\lambda^{{\bar{k}_{\rho}}}_{\rho}{(\boldsymbol{{l}}_{t})}\\}$
$\displaystyle\leq\min\\{r_{\rho}^{t}-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime}),\lambda^{{\bar{k}_{\rho}}}_{\rho}{(\boldsymbol{{l}}_{t})}\\},$
(29)

or equivalently:

$\displaystyle\min\\{r_{\rho}^{t},Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})+\lambda^{{\bar{k}_{\rho}}}_{\rho}{(\boldsymbol{{l}}_{t})}\\}-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})$
$\displaystyle\leq\min\\{r_{\rho}^{t},Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime})+\lambda^{{\bar{k}_{\rho}}}_{\rho}{(\boldsymbol{{l}}_{t})}\\}-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime}).$
(30)

Observe that

${z}^{k}_{\rho}(\boldsymbol{{l}}_{t},\bar{\boldsymbol{{x}}}^{\prime\prime})=\begin{cases}{z}^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})&\text{if
}k\neq{\bar{k}_{\rho}},\\\
{z}^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})+\underset{=1}{\underbrace{x^{\bar{v}}_{\bar{m}}}}\cdot
l^{t,\bar{v}}_{\rho,\bar{m}}\stackrel{{\scriptstyle\text{\eqref{eq:several-
definitions}}}}{{=}}{z}^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})+\lambda^{{\bar{k}_{\rho}}}_{\rho}{(\boldsymbol{{l}}_{t})}&\text{if
}k={\bar{k}_{\rho}},\vspace{-.5em}\end{cases}$

and thus

$\displaystyle
Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\bar{\boldsymbol{{x}}}^{\prime\prime})=\begin{cases}Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})\quad\text{if}\,\,k<{\bar{k}_{\rho}}\\\
\min\left\\{\textstyle
r_{\rho}^{t},\sum^{k}_{k^{\prime}=1}{z}^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})+\lambda^{{\bar{k}_{\rho}}}_{\rho}{(\boldsymbol{{l}}_{t})}\right\\}=\min\\{r_{\rho}^{t},Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})+\lambda^{{\bar{k}_{\rho}}}_{\rho}{(\boldsymbol{{l}}_{t})}\\}\quad\text{if}\,\,k\geq{\bar{k}_{\rho}}.\end{cases}$
(31)

Note that the same equality holds between $\bar{\boldsymbol{{x}}}^{\prime}$
and $\boldsymbol{{x}}^{\prime}$ since $(\bar{v},\bar{m})\not\in S^{\prime}$.

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

$\displaystyle
f_{t}(S^{\prime\prime}\cup\\{(\bar{v},\bar{m})\\})-f_{t}(S^{\prime\prime})=G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\bar{\boldsymbol{{x}}}^{\prime\prime})-G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:gain-
compact}}}{{=}}\sum_{\rho\in\mathcal{R}}\Bigg{[}\sum_{k=1}^{{\bar{k}_{\rho}}-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\left(Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})\right)$
$\displaystyle\qquad\qquad+\sum_{k={\bar{k}_{\rho}}}^{K_{\rho}-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\left(\min\left\\{r_{\rho}^{t},Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})+\lambda^{{\bar{k}_{\rho}}}_{\rho}{(\boldsymbol{{l}}_{t})}\right\\}-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})\right)$
$\displaystyle\qquad\qquad-\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\left(Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})\right)\Bigg{]}$
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k={\bar{k}_{\rho}}}^{K_{\rho}-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\Bigg{[}\left(\min\left\\{r_{\rho}^{t},Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})+\lambda^{{\bar{k}_{\rho}}}_{\rho}{(\boldsymbol{{l}}_{t})}\right\\}-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})\right)-\left(Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})\right)\Bigg{]}$
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k={\bar{k}_{\rho}}}^{K_{\rho}-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\left(\min\left\\{r_{\rho}^{t},Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})+\lambda^{{\bar{k}_{\rho}}}_{\rho}{(\boldsymbol{{l}}_{t})}\right\\}-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})\right).$
(32)

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

$f_{t}(S^{\prime\prime}\cup\\{(\bar{v},\bar{m})\\})-f_{t}(S^{\prime\prime})\leq
f_{t}(S^{\prime}\cup\\{(\bar{v},\bar{m})\\})-f_{t}(S^{\prime}).$

We conclude that $f_{t}$ is a submodular set function.

Monotonicity. The function $f_{t}\colon
2^{\mathcal{V}\times\mathcal{M}}\to\mathbb{R}$ is monotone [62] if for every
$S^{\prime}\subset S^{\prime\prime}\subset\mathcal{V}\times\mathcal{M}$ it
holds that

$f_{t}(S^{\prime\prime})\geq f_{t}(S^{\prime}).$ (33)

We have

$\displaystyle f_{t}(S^{\prime\prime})-f_{t}(S^{\prime})$
$\displaystyle=G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},{\boldsymbol{{x}}}^{\prime\prime})-G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime})$
(34)
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\left(Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime\prime})-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}}^{\prime})\right)$
(35) $\displaystyle\geq 0.$ (36)

The last inequality is obtained using Eq. (28). We conclude that $f_{t}$ is
monotone. ∎

######  Lemma A.2.

The set function $F_{T}:2^{\mathcal{V}\times\mathcal{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 [62]. Moreover, each $f_{t}$ is normalized, then it follows that
$F_{T}$ is normalized.

∎

##  Appendix B Equivalent Expression of the Gain Function

######  Lemma B.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 $\boldsymbol{{l}}_{t}$ and allocation vector $\boldsymbol{{x}}$. For
brevity, let us denote
$z^{k^{\prime}}_{\rho}=z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})$.
The following formula holds:

$\displaystyle\min\left\\{c,\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}\right\\}-\min\left\\{c,\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}\right\\}=\min\left\\{c-\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho},z^{k}_{\rho}\right\\}\cdot\mathds{1}\left\\{\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}<c\right\\}.$
(37)

###### Proof.

We distinguish two cases:

  1. (I)

When the $k-1$ less costly models have at least $c$ effective capacity, i.e.,
$\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}\geq c$, we obtain:

$\displaystyle\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}=\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}+z^{k}_{\rho}\geq
c+z^{k}_{\rho}\geq c.$ (38)

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

  2. (II)

When $\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}<c$, we obtain:

$\displaystyle\min\left\\{c,\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}\right\\}-\min\left\\{c,\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}\right\\}$
$\displaystyle=\min\left\\{c,\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}\right\\}-\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}=\min\left\\{c-\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho},z^{k}_{\rho}\right\\}.$
(39)

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

By combining Eq. (38) and Eq. (39) we obtain Eq. (37). ∎

######  Lemma B.2.

The cost function given by Eq. (12) can be expressed as:

$\displaystyle
C(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}})=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k}_{\rho}-\gamma^{k+1}_{\rho}\right)\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}+\gamma^{K_{\rho}}_{\rho}r^{t}_{\rho}.$
(40)

###### Proof.

The sum
$\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{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}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\geq
r^{t}_{\rho}$ (see Eq. (9)) and

$\displaystyle\gamma^{K_{\rho}}_{\rho}\min\left\\{r^{t}_{\rho},\sum^{K_{\rho}}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}=\gamma^{K_{\rho}}_{\rho}r^{t}_{\rho}.$
(41)

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

$\displaystyle
C(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}})\;$
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}}\gamma^{k}_{\rho}\cdot\min\left\\{r^{t}_{\rho}\;{-}\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}}),\hskip
1.00006ptz^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}\cdot\mathds{1}_{\left\\{\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})<r^{t}_{\rho}\right\\}}$
(42)
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}}\gamma^{k}_{\rho}\left(\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}-\min\left\\{r^{t}_{\rho},\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}\right)$
(43)
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:repo_inclusion}}}{{=}}\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\gamma^{k}_{\rho}\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}-\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}}\gamma^{k}_{\rho}\min\left\\{r^{t}_{\rho},\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}+\gamma^{K_{\rho}}_{\rho}r^{t}_{\rho}$
(44)
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\gamma^{k}_{\rho}\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}-\sum_{\rho\in\mathcal{R}}\sum_{k=2}^{K_{\rho}}\gamma^{k}_{\rho}\min\left\\{r^{t}_{\rho},\sum^{k-1}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}+\gamma^{K_{\rho}}_{\rho}r^{t}_{\rho}$
(45)
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\gamma^{k}_{\rho}\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}-\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\gamma^{k+1}_{\rho}\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}+\gamma^{K_{\rho}}_{\rho}r^{t}_{\rho}$
(46)
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k}_{\rho}-\gamma^{k+1}_{\rho}\right)\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}+\gamma^{K_{\rho}}_{\rho}r^{t}_{\rho}.$
(47)

∎

Proof of Lemma III.1.

###### Proof.

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

$\displaystyle
G(\boldsymbol{{r}},\boldsymbol{{l}}_{t},\boldsymbol{{x}})=C(\boldsymbol{{r}},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})-C(\boldsymbol{{r}},\boldsymbol{{l}}_{t},\boldsymbol{{x}})$
(48)
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k}_{\rho}-\gamma^{k+1}_{\rho}\right)\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}{z}^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})\right\\}-\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k}_{\rho}-\gamma^{k+1}_{\rho}\right)\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}$
(49)
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k}_{\rho}-\gamma^{k+1}_{\rho}\right)\cdot\left\\{\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})\right\\}-\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}{z}^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}\right\\}$
(50)
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\cdot\left\\{\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}{z}^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{x}})\right\\}-\min\left\\{r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})\right\\}\right\\}$
(51)
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right){\left(Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{x}})-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})\right)}.$
(52)

∎

##  Appendix C Projection Algorithm

In order to project a fractional allocation $\boldsymbol{{y}}^{\prime}$ lying
outside the constraint set $\mathcal{Y}$ to a feasible allocation
$\boldsymbol{{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

$\displaystyle
D_{\Phi}(\boldsymbol{{y}},\boldsymbol{{y}}^{\prime})=\Phi(\boldsymbol{{y}})-\Phi(\boldsymbol{{y}}^{\prime})-{\nabla\Phi(\boldsymbol{{y}}^{\prime})}^{T}(\boldsymbol{{y}}-\boldsymbol{{y}}^{\prime}),$
(53)

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

$\displaystyle D_{\Phi^{v}}(\boldsymbol{{y}}^{v},\boldsymbol{{y}}^{\prime
v})={\Phi^{v}}(\boldsymbol{{y}}^{v})-{\Phi^{v}}(\boldsymbol{{y}}^{\prime
v})-{\nabla{\Phi^{v}}(\boldsymbol{{y}}^{\prime
v})}^{T}(\boldsymbol{{y}}^{v}-\boldsymbol{{y}}^{\prime v}).$ (54)

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

$\displaystyle\boldsymbol{{y}}=\textstyle\prod_{\mathcal{Y}\cap\mathcal{D}}^{\Phi}(\boldsymbol{{y}}^{\prime})$
$\displaystyle=\underset{\boldsymbol{{y}}\in\mathcal{Y}\cap\mathcal{D}}{\mathrm{argmin}}\,D_{\Phi}(\boldsymbol{{y}},\boldsymbol{{y}}^{\prime}).$
(55)

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.,

$\displaystyle
D_{\Phi}(\boldsymbol{{y}},\boldsymbol{{y}}^{\prime})=\sum_{v\in\mathcal{V}}D_{\Phi^{v}}(y^{v},y^{\prime
v})$ (56)

where
$\boldsymbol{{y}}\in\mathcal{Y}=\bigtimes_{v\in\mathcal{V}}\mathcal{Y}^{v}$,
and $\boldsymbol{{y}}^{v}\in\mathcal{Y}^{v},\forall v\in\mathcal{V}$. In order
to minimize the value $D_{\Phi}(\boldsymbol{{y}},\boldsymbol{{y}}^{\prime})$
for $\boldsymbol{{y}}\in\mathcal{Y}\cap\mathcal{D}$, we can independently
minimize the values $D_{\Phi^{v}}(y^{v},y^{\prime v})$ for
$\boldsymbol{{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

$\displaystyle\boldsymbol{{y}}^{v}=\textstyle\prod_{\mathcal{Y}^{v}\cap\mathcal{D}^{v}}^{\Phi^{v}}(\boldsymbol{{y}}^{\prime
v})$
$\displaystyle=\underset{\boldsymbol{{y}}^{v}\in\mathcal{Y}^{v}\cap\mathcal{D}^{v}}{\mathrm{argmin}}\,D_{\Phi^{v}}(\boldsymbol{{y}}^{v},\boldsymbol{{y}}^{\prime
v}).$ (57)

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

1:$|\mathcal{M}|$; $b^{v}$; $\boldsymbol{{s}}^{v}$; Sorted
$\boldsymbol{{y}}^{\prime v}$ where ${{y_{|\mathcal{M}|}^{\prime
v}}}\geq\dots\geq{{y_{1}^{\prime v}}}$

2:${y_{|\mathcal{M}|+1}^{\prime v}}\leftarrow+\infty$

3:for $k\in\\{|\mathcal{M}|,|\mathcal{M}|-1,\dots,1\\}$ do

4:
$m_{k}\leftarrow\frac{b^{v}-\sum^{|\mathcal{M}|}_{m=k+1}s^{v}_{m}}{\sum_{m=1}^{k}s^{v}_{m}{y_{m}^{\prime
v}}}$

5: if ${y_{k}^{\prime v}}m_{k}<1\leq{y_{k+1}^{\prime v}}m_{k}$
then$\triangleright$ Appropriate $k$ is found

6: for $k^{\prime}\in\\{1,2,\dots,k\\}$ do

7: $y^{v}_{k^{\prime}}\leftarrow m_{k}y^{\prime v}_{k^{\prime}}$
$\triangleright$ Scale the variable’s components

8: for $k^{\prime}\in\\{k+1,k+2,\dots,|\mathcal{M}|\\}$ do

9: $y^{v}_{k^{\prime}}\leftarrow 1$$\triangleright$ Cap the variable’s
components to 1

10: return $\boldsymbol{{y}}^{v}$$\triangleright$ $\boldsymbol{{y}}^{v}$ is
the result of the projection

######  Theorem C.1.

Algorithm 2 when executed at node $v\in\mathcal{V}$ returns
$\Pi^{\Phi^{v}}_{\mathcal{Y}^{v}\cap\mathcal{D}^{v}}({\boldsymbol{{y}}^{\prime
v}})$, i.e., the projection of the vector $\boldsymbol{{y}}^{\prime v}$ onto
the weighted capped simplex $\mathcal{Y}^{v}\cap\mathcal{D}^{v}$ under the
weighted negative entropy
$\Phi^{v}(\boldsymbol{{y}}^{v})=\sum_{m\in\mathcal{M}}s^{v}_{m}y^{v}_{m}\log(y^{v}_{m})$.
The time complexity of the projection is
$\mathcal{O}\left(|\mathcal{M}|\log(|\mathcal{M}|)\right)$.

###### Proof.

$\displaystyle\textstyle\prod_{\mathcal{Y}^{v}\cap\mathcal{D}^{v}}^{\Phi^{v}}(\boldsymbol{{y}}^{\prime
v})$
$\displaystyle=\underset{\boldsymbol{{y}}^{v}\in\mathcal{Y}^{v}\cap\mathcal{D}^{v}}{\mathrm{argmin}}\,D_{{\Phi^{v}}}(\boldsymbol{{y}}^{v},\boldsymbol{{y}}^{\prime
v})$ (58)
$\displaystyle=\underset{\boldsymbol{{y}}^{v}\in\mathcal{Y}^{v}\cap\mathcal{D}^{v}}{\mathrm{argmin}}\sum_{m\in\mathcal{M}}s^{v}_{m}\left({y^{v}_{m}}\mathrm{log}\left(\frac{{y^{v}_{m}}}{{y^{\prime
v}_{m}}}\right)-{y^{v}_{m}}+{y^{\prime v}_{m}}\right).$ (59)

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

$\displaystyle\mathcal{J}(\boldsymbol{{y}}^{v},\mathbf{\beta},\tau)=\sum_{m\in\mathcal{M}}{s^{v}_{m}}\left({y^{v}_{m}}\mathrm{log}\left(\frac{{y^{v}_{m}}}{{y^{\prime
v}_{m}}}\right)-{y^{v}_{m}}+{y^{\prime
v}_{m}}\right)-\sum_{m\in\mathcal{M}}\beta_{m}\left(1-{y^{v}_{m}}\right)-\tau\left(\sum_{m\in\mathcal{M}}{s^{v}_{m}}{y^{v}_{m}}-b^{v}\right).$
(60)

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

$\displaystyle{s^{v}_{m}}\mathrm{log}({\hat{y}^{v}_{m}})-{s^{v}_{m}}\mathrm{log}({y^{\prime
v}_{m}})+\beta_{m}-{s^{v}_{m}}\tau=0,$ (61a)
$\displaystyle{\hat{y}^{v}_{m}}\leq 1,$ (61b) $\displaystyle\beta_{m}\geq 0,$
(61c) $\displaystyle\sum_{m\in\mathcal{M}}{s^{v}_{m}}{\hat{y}^{v}_{m}}=b^{v},$
(61d) $\displaystyle\beta_{m}(1-{\hat{y}^{v}_{m}})=0.$ (61e)

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

$\displaystyle{\hat{y}}^{v}_{1}\leq...\leq{\hat{y}}^{v}_{k}<{\hat{y}}^{v}_{k+1}=...={\hat{y}}^{v}_{|\mathcal{M}|}=1\,\,\,\text{if}\,k<{|\mathcal{M}|},$
(62)
$\displaystyle{\hat{y}}^{v}_{1}\leq{\hat{y}}^{v}_{2}\leq...\leq{\hat{y}}^{v}_{|\mathcal{M}|}<1\,\,\,\text{if}\,k={|\mathcal{M}|}.$
(63)

The goal here is to identify a valid value for $k$ (number of components of
$\boldsymbol{{\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,...,k$, we have from (61e) that $\beta_{m_{L}}=0$, and then from
(61a), $s^{v}_{m_{L}}\mathrm{log}(y^{\prime
v}_{m_{L}})+s^{v}_{m_{L}}\tau=s^{v}_{m_{L}}\mathrm{log}({\hat{y}}^{v}_{m_{L}})<s^{v}_{m_{L}}\mathrm{log}(1)=0$,
and can be simplified to

$\displaystyle y^{\prime
v}_{m_{L}}e^{\tau}<1,\forall{m_{L}}\in\\{1,\dots,k\\}.$ (64)

  * •

For ${m_{U}}=k+1,...,{|\mathcal{M}|}$: as $\beta_{m_{U}}\geq 0$ from (61c), we
get
$0=s^{v}_{m_{U}}\mathrm{log}({\hat{y}}^{v}_{m_{U}})=s^{v}_{m_{U}}\mathrm{log}(y^{\prime
v}_{m_{U}})-\beta_{m_{U}}+s^{v}_{m_{U}}\tau\leq
s^{v}_{m_{U}}\mathrm{log}(y^{\prime v}_{m_{U}})+s^{v}_{m_{U}}\tau$, and can be
simplified to

$\displaystyle y^{\prime v}_{m_{U}}e^{\tau}\geq
1,\forall{m_{U}}\in\\{k+1,\cdots,{|\mathcal{M}|}\\}.$ (65)

Consider Eqs. (62) and (63), and since for $m_{L}\in\\{1,\dots,k\\}$ we have
$y^{\prime v}_{m_{L}}e^{\tau}=\hat{y}^{v}_{m_{L}}$ (the order is preserved),
then the conditions in Eq. (64) are

$\displaystyle y^{\prime v}_{1}e^{\tau}\leq...\leq y^{\prime
v}_{k}e^{\tau}<1.$ (66)

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

$y^{\prime v}_{k}e^{\tau}<1\leq y^{\prime v}_{k+1}e^{\tau}.$

Note that the r.h.s inequality is ignored when $k=|\mathcal{M}|$ by
construction ($y^{\prime v}_{|\mathcal{M}|+1}=+\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. (61d):

$b^{v}=\sum^{|\mathcal{M}|}_{m=1}{{s^{v}_{m}}{\hat{y}^{v}_{m}}}=\sum^{|\mathcal{M}|}_{m=k+1}{s^{v}_{m}}+e^{\tau}\sum_{m=1}^{k}{s^{v}_{m}}{y^{\prime
v}_{m}}.$

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

$\displaystyle m_{k}\triangleq e^{\tau}$
$\displaystyle=\frac{b^{v}-\sum^{|\mathcal{M}|}_{m=k+1}{s^{v}_{m}}}{\sum_{m=1}^{k}s^{v}_{m}{y^{\prime
v}_{m}}}.$ (67)

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

$\displaystyle y^{\prime v}_{k}m_{k}<1\leq y^{\prime v}_{k+1}m_{k}.$ (68)

The appropriate $k$ satisfying the KKT conditions is contained in
$\mathcal{M}$, and due to the sorting operation this gives total time
complexity of $\mathcal{O}\left(|\mathcal{M}|\log(|\mathcal{M}|)\right)$ per
iteration. In practice, the online mirror ascent method quickly sets
irrelevant items in the fractional allocation vector $\boldsymbol{{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. ∎

##  Appendix D Subgradient Expression

######  Lemma D.1.

The gain function in Eq. (16) has a subgradient $\boldsymbol{{g}}_{t}$ at
point $\boldsymbol{{y}}_{t}\in\mathcal{Y}$ given by

$\displaystyle\boldsymbol{{g}}_{t}=\left[\sum_{\rho\in\mathcal{R}}l^{t,v}_{\rho,m}\left(\gamma^{K^{*}_{\rho}(\boldsymbol{{y}}_{t})}_{\rho}-C^{v}_{\boldsymbol{{p}},m}\right)\mathds{1}_{\\{\kappa_{\rho}(v,m)<K^{*}_{\rho}(\boldsymbol{{y}}_{t})\\}}\right]_{(v,m)\in\mathcal{V}\times\mathcal{M}},$
(69)

where
$K^{*}_{\rho}(\boldsymbol{{y}}_{t})=\min\big{\\{}k\in[K_{\rho}-1]:\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})\geq
r^{t}_{\rho}\big{\\}}$.

###### Proof.

The function given by
$Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})=\min\left\\{\textstyle
r^{t}_{\rho},\sum^{k}_{k^{\prime}=1}{z}^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})\right\\}$
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 [67, Theorem 8.2], at point
$\boldsymbol{{y}}_{t}\in\mathcal{Y}$ as

$\displaystyle\textstyle\partial
Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})=\begin{cases}\left\\{\nabla(\sum^{k}_{k^{\prime}=1}{z}^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t}))\right\\}&\mathrm{if}\,\,Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})<r^{t}_{\rho}\qquad(\text{r.h.s.
argument of the min is active}),\\\
\mathrm{conv}\left(\left\\{\boldsymbol{{0}},\nabla(\sum^{k}_{k^{\prime}=1}{z}^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t}))\right\\}\right)&\mathrm{if}\,\,Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})=r^{t}_{\rho}\qquad(\text{both
arguments of the min are active}),\\\
\\{\boldsymbol{{0}}\\}&\mathrm{otherwise}\ \qquad\qquad\qquad(\text{l.h.s.
argument of the min is active}),\end{cases}$ (70)

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

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

$\displaystyle\tilde{\boldsymbol{{g}}}^{k}_{\rho,t}=\begin{cases}\nabla(\sum^{k}_{k^{\prime}=1}{z}^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t}))&\mathrm{if}\,\,Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})<r^{t}_{\rho},\\\
\boldsymbol{{0}},&\mathrm{otherwise}.\end{cases}$ (71)

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

$\displaystyle\frac{\partial\phantom{y^{v}_{m}}}{\partial
y^{v}_{m}}\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})$
$\displaystyle=\sum^{k}_{k^{\prime}=1}\frac{\partial\phantom{y^{v}_{m}}}{\partial
y^{v}_{m}}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})\stackrel{{\scriptstyle\eqref{eq:several-
definitions}}}{{=}}l^{t,v}_{\rho,m}\cdot\mathds{1}_{\\{\kappa_{\rho}(v,m)\leq
k\\}}.$

The indicator variable $\mathds{1}_{\\{\kappa_{\rho}(v,m)\leq k\\}}$ is
introduced since the partial derivative of
$\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})$
w.r.t. $y^{v}_{m}$ 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^{v}_{m}$ appears once in the summation).

We obtain from Eq. (71)

$\displaystyle\tilde{g}^{k,v}_{\rho,t,m}$
$\displaystyle=\begin{cases}l^{t,v}_{\rho,m}\cdot\mathds{1}_{\\{\kappa_{\rho}(v,m)\leq
k\\}}&\mathrm{if}\,\,Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})<r^{t}_{\rho},\\\
0&\mathrm{otherwise.}\end{cases}$ (72)
$\displaystyle=l^{t,v}_{\rho,m}\cdot\mathds{1}_{\\{\kappa_{\rho}(v,m)\leq
k\,\,\land\,\,Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})<r^{t}_{\rho}\\}},\forall(v,m)\in\mathcal{V}\times\mathcal{M}.$
(73)

By considering the subdifferential

$\displaystyle\partial
G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})=\partial\left(\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\left(Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}})-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})\right)\right),$
(74)

and using [68, Theorem 23.6], we get

$\displaystyle\partial
G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\partial\left(\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\left(Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}})-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})\right)\right).$
(75)

The constant factors $\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)$ are
non-negative, so we can multiply both sides of the subgradient inequality by a
non-negative constant [68, Sec. 23]; furthermore, the subgradient of the
constants
$Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})$
is $\boldsymbol{{0}}$. We get

$\displaystyle\partial
G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\partial\left(Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}})-Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{\omega}})\right)=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\partial
Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}).$
(76)

Then, a subgradient $\boldsymbol{{g}}_{t}\in\partial
G(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})$ at point
$\boldsymbol{{y}}_{t}\in\mathcal{Y}$ is given by

$\displaystyle\boldsymbol{{g}}_{t}=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\tilde{\boldsymbol{{g}}}^{k}_{\rho,t}.$
(77)

The $(v,m)$-th component of the subgradient $\boldsymbol{{g}}_{t}$ is

$\displaystyle g^{v}_{t,m}$
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\cdot\tilde{g}^{k,v}_{\rho,t,m}$
(78)
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=1}^{K_{\rho}-1}l^{t,v}_{\rho,m}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\cdot\mathds{1}_{\\{\kappa_{\rho}(v,m)\leq
k\,\,\land\,\,Z^{k}_{\rho}(\boldsymbol{{r}}_{t},\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})<r^{t}_{\rho}\\}}$
(79)
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=\kappa_{\rho}(v,m)}^{K_{\rho}-1}l^{t,v}_{\rho,m}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)\cdot\mathds{1}_{\\{\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})<r^{t}_{\rho}\\}}$
(80)
$\displaystyle=\sum_{\rho\in\mathcal{R}}\sum_{k=\kappa_{\rho}(v,m)}^{K^{*}_{\rho}(\boldsymbol{{y}}_{t})-1}l^{t,v}_{\rho,m}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)$
(81)
$\displaystyle=\sum_{\rho\in\mathcal{R}}l^{t,v}_{\rho,m}\cdot\sum_{k=\kappa_{\rho}(v,m)}^{K^{*}_{\rho}(\boldsymbol{{y}}_{t})-1}\left(\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}\right)$
(82)
$\displaystyle=\sum_{\rho\in\mathcal{R}}l^{t,v}_{\rho,m}\left(\gamma^{K^{*}_{\rho}(\boldsymbol{{y}}_{t})}_{\rho}-\gamma^{\kappa_{\rho}(v,m)}_{\rho}\right)\cdot\mathds{1}_{\\{\kappa_{\rho}(v,m)<K^{*}_{\rho}(\boldsymbol{{y}}_{t})\\}}$
(83) $\displaystyle\stackrel{{\scriptstyle\eqref{eq:several-
definitions}}}{{=}}\sum_{\rho\in\mathcal{R}}l^{t,v}_{\rho,m}\left(\gamma^{K^{*}_{\rho}(\boldsymbol{{y}}_{t})}_{\rho}-C^{v}_{\boldsymbol{{p}},m}\right)\cdot\mathds{1}_{\\{\kappa_{\rho}(v,m)<K^{*}_{\rho}(\boldsymbol{{y}}_{t})\\}},\forall(v,m)\in\mathcal{V}\times\mathcal{M},$
(84)

where
$K^{*}_{\rho}(\boldsymbol{{y}}_{t})=\min\big{\\{}k\in[K_{\rho}-1]:\sum^{k}_{k^{\prime}=1}z^{k^{\prime}}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}}_{t})\geq
r^{t}_{\rho}\big{\\}}$.

∎

##  Appendix E Supporting Lemmas for the Proof of Theorem V.1

###  E-A Concavity of the Gain Function

######  Lemma E.1.

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

###### Proof.

Since $\lambda^{k}_{\rho}$ is defined to be the $k$-th smallest cost for any
$k\in[K_{\rho}]$ (see Eq. (11)), then the factors
$\gamma^{k+1}_{\rho}-\gamma^{k}_{\rho}$ are always non-negative. Moreover
$z^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}})=y^{v}_{m}l^{t,v}_{\rho,m}$,
where $v,m$ are such that $\kappa_{\rho}(v,m)=k$. Therefore,
$Z^{k}_{\rho}(\boldsymbol{{l}}_{t},\boldsymbol{{y}})$ is the minimum between a