Luigi Vigneri

Georgios Paschos

Panayotis Mertikopoulos

Large-Scale Network Utility Maximization: Countering Exponential Growth with Exponentiated Gradients

Network utility maximization (NUM) is an iconic problem in network traffic management which is at the core of many current and emerging network design paradigms -and, in particular, software-defined networks (SDNs). Thus, given the exponential growth of modern-day networks (in both size and complexity), it is crucial to develop scalable algorithmic tools that are capable of providing efficient solutions in time which is dimension-free, i.e., independent -or nearly-independent -on the size of the system. To do so, we leverage a suite of modified gradient methods known as "mirror descent", and we derive a scalable and efficient algorithm for the NUM problem based on gradient exponentiation. We show that the convergence speed of the proposed algorithm only carries a logarithmic dependence on the size of the network, so it can be implemented reliably and efficiently in massively large networks where traditional gradient methods are prohibitively slow. These theoretical results are subsequently validated by extensive numerical simulations showing an improvement of several order of magnitudes over standard gradient methods in large-scale networks.

I. Introduction

One of the most important design attributes of cloud-based networking is the ability to adjust the allocation of bandwidth and computing resources in real-time. This allows cloud-based networks to adapt to variable loads and demands "on the fly", leading in turn to dramatic performance gains in throughput, end-to-end connectivity, and reliability.

From a network control standpoint, a key obstacle to achieving these gains is that they require solving, rapidly and efficiently, complex optimization problems with several thousands -if not millions -of control variables. This is most readily seen in the reference problem of network utility maximization (NUM), a landmark framework pioneered by Kelly et al. [START_REF] Kelly | Rate control for communication networks: shadow prices, proportional fairness and stability[END_REF] to ensure fairness and stability in communication networks through the maximization of a concave utility function subject to link capacity constraints. Owing to its flexibility, this framework has found widespread applications to wireless networks [START_REF] Palomar | A tutorial on decomposition methods for network utility maximization[END_REF], sensor networks [START_REF] Deng | Maximizing network utility of rechargeable sensor networks with spatiotemporally coupled constraints[END_REF], caching systems [START_REF] Dehghan | A utility optimization approach to network cache design[END_REF], and many other fields (for a survey, see [START_REF] Shakkottai | Network optimization and control[END_REF] and references therein). However, as such networks grow in size and complexity, it is crucial to develop scalable algorithms which could provide reasonable solutions to the NUM problem in real time.

This situation is exacerbated by the massive penetration of software-defined networks (SDNs) to the access and transport layers. Here, network controllers are required to centrally manage the allocation of bandwidth resources to an exponentially large number of sessions (often numbering in the millions for IP access network domains). As such, any NUM algorithm must be dimension-free, i.e., its execution runtime must be independentor nearly-independent -of the size of the network.

A. Our contributions

In this paper, we follow the classical formulation of the NUM problem as a concave program over a network with d concurrent sessions [START_REF] Shakkottai | Network optimization and control[END_REF]. In contrast to application-driven formulations, this generalist approach provides us with a flexible theoretical springboard for the use of state-of-the-art optimization methods (which can then be ported to specific scenarios with minimal modifications). More concretely, motivated by recent breakthroughs in artificial intelligence and machine learning, we leverage an optimization technique known as "mirror descent" (or "ascent" in our case), and which was originally introduced by Nemirovski and Yudin [START_REF] Nemirovsky | Problem Complexity and Method Efficiency in Optimization[END_REF] to solve image reconstruction problems with massively large data sets. In contrast to vanilla gradient methods, mirror ascent replaces standard (projected) gradient steps by a variable-geometry step which intrinsically accelerates the method in regions of the problem's feasible space where larger step-sizes might be more beneficial. In this way, depending on the geometry of the problem, mirror ascent methods often succeed in lifting the so-called "curse of dimensionality", and they provide convergence rates that vastly outperform unmodulated gradient methods.

By coupling mirror descent methodologies with a synthetic modification of the geometry of the problem at hand, we derive a fast and scalable interior-point method based on the exponentiation of network utility gradients, which we call exp-NUM (exponentiated gradients algorithm for NUM). In terms of performance gains, the proposed exp-NUM algorithm achieves an ε-optimal solution of the NUM problem in O(log d/ε 2 ) iterations, i.e., in time which is O(log d) relative to the size of the network -and hence, effectively dimension-free. As a result, thanks to this exponential acceleration factor, exp-NUM can be implemented reliably and efficiently in massively large networks where traditional first-order methods are prohibitively slow. This theoretical analysis is subsequently validated by an extensive suite of numerical simulations that exhibit exponential improvement over standard methods in large-scale random networks, and which clearly indicate the benefits of the proposed algorithm for large-scale network utility maximization.

B. Relation to existing work

The performance of an iterative NUM algorithm can be characterized in terms of (a) the algorithm's iteration complexity, i.e., the number of iterations required to achieve a given error (typically measured in terms of the value of the problem's objective function); and (b) the algorithm's periteration complexity, i.e., the number of floating point operations required to perform an update. To the best of our knowledge, the fastest schemes available in terms of iteration complexity are based on Newton's method: if the NUM objective satisfies a self-concordance condition, the authors of [START_REF] Wei | A distributed newton method for network utility maximization; part ii: Convergence[END_REF] showed that an ε-optimal solution can be achieved in O(1/ √ ε) iterations, improving in this way the O(1/ε) bound of [START_REF] Athuraliya | Optimization flow control with newton-like algorithm[END_REF]. On the other hand, Newton's method requires the computation of a Hessian matrix at each queried point and the solution of a d × d linearsystem per update; as a result, the per-iteration complexity of Newton-based schemes is O(d 3 ), so it becomes prohibitive when d is of the order of a few thousands. The primal-dual algorithm of [START_REF] Beck | An o(1/k) gradient method for network resource allocation problems[END_REF] also achieves an iteration complexity of O(poly(d)/ε) (or O(poly(d)/ √ ε) when coupled with Nesterov acceleration), but assumes strong concavity of the NUM objective, which subsumes Hessian invertibility and does not hold in networks with several overlapping origin-destination (O/D) paths.

In terms of per-iteration complexity, the fastest known methods rely on first-order information, and, for this reason, they are the benchmark reference solutions for large-scale NUM problems. Within this class, (projected) gradient ascent [START_REF] Bertsekas | Nonlinear Programming[END_REF][START_REF] Shakkottai | Network optimization and control[END_REF] and dual (Lagrangian) ascent [START_REF] Low | Optimization flow control. i. basic algorithm and convergence[END_REF] require O(d/ε 2 ) iterations to achieve an ε-optimal solution. An appealing alternative is the distributed backpressure algorithm, where nodes route packets according to queue backlog differences in order to balance the backlogs in the network. Backpressure can be understood as a stochastic dual subgradient algorithm [START_REF] Zargham | Accelerated backpressure algorithm[END_REF] with the queue backlogs playing the role of stochastic Lagrangian multipliers, thus leading to an O(d/ε 2 ) iteration complexity.

This convergence speed can be improved to O(d/ε) by using an alternating direction method of multipliers (ADMM) approach [START_REF] Allybokus | Real-time fair resource allocation in distributed software defined networks[END_REF][START_REF] Gupta | Centralized network utility maximization over aggregate flows[END_REF]. ADMM updates are relatively cheap, but the method requires additional assumptions on the NUM objective (viz. strong concavity and strong smoothness) which may fail to hold in general. These assumptions were recently lifted by the authors of [START_REF] Yu | A simple parallel algorithm with an $o(1/t)$ convergence rate for general convex programs[END_REF] who provide a virtual queues algorithm with an O(d/ε) convergence speed; however, the proposed algorithm employs a complicated prox-mapping which makes its periteration complexity prohibitive for large networks.

By comparison, the iteration complexity of the exp-NUM algorithm is O(log d/ε 2 ), i.e., it converges to an ε-optimal solution in logarithmic time relative to the problem's dimension (in our case, the number of sessions). At the same time, the algorithm's per-iteration complexity is the cheapest among all methods, so this acceleration does not come at the expense of costlier updates. As such, of the algorithms discussed above, exp-NUM is the only one whose complexity provably scales with the size of the network. In practice, this means a runtime reduction of Θ(d), amounting to an improvement of 4 orders of magnitude in large-scale networks with O(10 5 ) sessions (and the gains only scale upwards as d becomes larger); We find this property of exp-NUM particularly appealing as it provides a very promising platform for the systematic design of scalable resource allocation protocols in massively large networks.

This dramatic runtime reduction by a factor of Õ(d) should not be interpreted as a blanket result that applies to all convex programs: it is a by-product of the simplicial structure of the convex program at hand which allows the use of a carefully crafted, geometry-aware regularization process. Thus, in addition to the novel application of mirror ascent methodologies, an important part of our contribution is the principled modification of the problem's geometry so as to be able to take full advantage of this tailor-made regularization process. We find this approach particularly appealing as it provides a springboard for the application of similar techniques to a wide range of resource allocation problems in the area of networking.

For convenience, we summarize the complexity results for first-order methods in Table I. For the purposes of validation, we also perform in Section V an extensive suite of simulations that shows the scalability gains of exp-NUM over traditional first-order methods in large-scale random networks.

II. System model and problem formulation

In this section, we describe the basic setup of the NUM problem; for our notational conventions, see Table II.

A. Network utility maximization

The core ingredients of the network utility maximization (NUM) problem are as follows:

• Network model: The underlying network is modeled as a directed multi-graph G = (N , L) where N is the set of the network's nodes and L is its set of links.

• Sessions: Each O/D pair in the network may carry a communication session. We will write S = {1, . . . , d} for the set of such sessions; we will also assume for now that each session s ∈ S is associated with a single O/D path, and transmits over said path at rate x s ≥ 0 (see Section IV-C for the multi-path case). • Utility function: We posit that, if the s-th session transmits at rate x s , its utility is given by U s (x s ), where U s : + → is assumed differentiable and concave.1 • Capacity constraints: Each link ∈ L is assumed to have a finite capacity c . In view of the above, the load on link ∈ L induced by a transmission rate profile x = (x s ) s∈S is given by

x = s x s , (1)
where the notation "s " signifies that link belongs to the path carrying the traffic of the s-th session. The network utility maximization (NUM) problem may then be formulated as

maximize V(x) ≡ s∈S U s (x s ) subject to x s ≥ 0 for all s ∈ S x ≤ c for all ∈ L (NUM 0 )
i.e., as the problem of maximizing the network's aggregate utility while respecting each link's individual capacity.

Concerning the session utility functions U s , the standard choice in the literature is the α-fair model [START_REF] Kelly | Rate control for communication networks: shadow prices, proportional fairness and stability[END_REF][START_REF] Radunović | A unified framework for max-min and min-max fairness with applications[END_REF][START_REF] Shakkottai | Network optimization and control[END_REF]]

U α (z) =        z 1-α 1-α if α ≥ 0, α 1, log z if α = 1.
(

) 2
This utility model is sufficiently flexible to account for several different optimality and fairness criteria (thus ensuring objective diversification in a network setting). In particular, for α = 0, (NUM 0 ) boils down to sum-rate maximization; for α = 1, solving (NUM 0 ) leads to proportionally fair rate allocations; for α = 2, (NUM 0 ) amounts to potential delay fairness (also linked to TCP congestion control in the Internet [START_REF] Mo | Fair end-to-end window-based congestion control[END_REF]); finally, as α → ∞, we get the standard max-min fairness formulation.

To increase the flexibility of our model, we further posit that each session is assigned a weight w s > 0 measuring its importance for the network at large (for instance, a session carrying emergency medical information should carry more weight than a routine music streaming session). With all this in mind, our core utility model will be of the form U s (x s ) = w s U α (x s ) for some fixed α (determined by the network administrator depending on the precise criterion to be optimized).

B. Penalty-based formulation

In large-scale networks, the number of sessions quickly grows massively large, rendering the rigid link capacity constraints in (NUM 0 ) effectively impossible to handle as stated. For this reason, it is very common in the literature [START_REF] Kelly | Rate control for communication networks: shadow prices, proportional fairness and stability[END_REF][START_REF] Shakkottai | Network optimization and control[END_REF] to consider instead a penalty-based formulation where the objective of (NUM 0 ) is written as

V(x) = s∈S U s (x s ) -µ ∈L g (x ), (3)
where µ > 0 is a rigidity parameter controlling the impact of the capacity penalty functions g : + → , ∈ L. Specifically, instead of maximizing the sum of utilities constrained by link capacities, the network administrator now maximizes the aggregate network utility minus a (potentially very steep) cost for overshooting the capacity of a given link. The benefit of this formulation is that the problem's constraints are now uncoupled across its control variables (x s , s ∈ S); on the other hand, this coupling now appears in the problem's objective (which was uncoupled in the original formulation), so the problem's interaction structure remains the same.

Regarding the choice of penalty functions g , the indicator function g (z) = ∞ if z > c and g (z) = 0 if z ≤ c exactly recovers (NUM 0 ) as stated above -but also exhibits the same drawbacks in terms of tractability. To reap the benefits of rigidity relaxation, it is instead more convenient to assume that each g is an increasing convex function, a typical choice being

g (x ) = [x -c ] q + , (4)
i.e., incurring a cost that grows as a q-th power (typically q = 1 or q = 2) when the load x = s x s on link exceeds the link's capacity.

Needless to say, this framework offers considerably more leeway and flexibility to the network adminstrator relative to the rigid formulation (NUM 0 ). First, from standard epiconvergence arguments [START_REF] Attouch | Variational Convergence for Functions and Operators[END_REF], the solutions of (NUM 0 ) can be recovered exactly from this relaxed formulation if g (c ) > 0 and µ is large enough. In addition, by controlling g , the administrator can over-or under-charge links depending on the situation of the network (for instance, by replacing c by c (1 ± ε ) for some ε > 0). For a detailed discussion, we refer the reader to [START_REF] Shakkottai | Network optimization and control[END_REF] and references therein.

For reasons that will become clear in the sequel, we take a slight modification of the above flexible framework and focus on the following network utility maximization problem maximize

V(x) ≡ s∈S U s (x s ) -µ ∈L g (x ) subject to x s ≥ 0 and s∈S x s ≤ C (NUM)
where C is a global capacity bound. This upper bound is an extrinsic parameter which is chosen by the network administrator: for instance, it can be chosen as the sum of link capacities c , as the maximum link capacity max c , or, more simply, as a means to implement a congestion control mechanism as in the classical paper of Kelly et al. [START_REF] Kelly | Rate control for communication networks: shadow prices, proportional fairness and stability[END_REF]. Importantly, we should also note here that it is possible to choose C in a way such that any solution of (NUM 0 ) is also a solution of (NUM).

All in all, the extra constraint in (NUM) should be treated as a synthetic geometry modification that does not alter the problem's solution set but instead guarantees that interim candidate solutions are well-behaved (in case a solution algorithm needs to be cut short of its execution runtime). Example 1 below provides some additional clarifications:

Example 1 (A toy example). Fig. 1 shows a network consisting of two links 1 and 2 with respective capacities c 1 = 1 and c 2 = 2. Let x 1 be the rate of the session that uses only 1 , x 2 that of the session using only 2 , and x 3 the rate of the session using both links, with corresponding utility functions U 1 (x 1 ) = x 1 , U 2 (x 2 ) = x 2 , and U 3 (x 3 ) = 3x 3 . The feasibility region of (NUM 0 ) is depicted by the green shaded area in Fig. 2; the larger, overlapping orange area represents the cap constraint of (NUM). In this case, it is easy to see that the solution of (NUM 0 ) is x * 1 = 0, x * 2 = 1, and x * 3 = 1. Instead of exploring the entire positive orthant as in the standard relaxation of the NUM problem [START_REF] Shakkottai | Network optimization and control[END_REF], the capped formulation introduces a cut-off which narrows the problem's feasible region while still remaining tractable and scalable from an algorithmic standpoint. Note also that, when one picks the global capacity constraint as above, the solution of (NUM 0 ) lies within the feasible region of (NUM) and can be recovered by choosing µ appropriately.

III. Gradient Ascent and Mirror Ascent

In this section, we start laying down the algorithmic groundwork for the exponentiated gradient algorithm for NUM (described in detail in Section IV). To that end, we begin with some standing notational conventions and assumptions:

1) We write Ω ≡ {x ∈ d + : x s ≥ 0 and s x s ≤ C} for the feasible region of (NUM) and V * ≡ max x∈Ω V(x) for the maximum value of V over Ω. 3) Finally, we make the blanket regularity assumption that

∂V ∂x s ≤ L (5)
for some L > 0 and all x ∈ Ω, s ∈ S. In words, this means that the derivatives of V are bounded over the problem's feasible region (i.e., V is Lipschitz continuous).

A. Projected gradient ascent

We begin by revisiting a standard benchmark for solving network utility maximization problems, that of projected gradient ascent (PGA). Following [START_REF] Bertsekas | Nonlinear Programming[END_REF], this is captured by the recursion

x k+1 = Π Ω (x k + γv k ) (PGA)
where:

1) k = 1, 2, . . . , is the algorithm's iteration counter.

2) x k ∈ Ω is the state of the algorithm at iteration k.

3) v k = ∇V(x k ) is the gradient of the NUM objective at x k . 4) Π Ω : d → Ω is the Euclidean projector Π Ω (w) = arg min w∈Ω { x -w 2 2 }. (6)
In words, the PGA algorithm takes a step along the gradient of the objective (modulated by the step-size parameter γ) and, if the result lies outside the problem's feasible region, it projects back to Ω and repeats the process.

To streamline our presentation, we defer the analysis of the algorithm until later in this section (cf. the discussion right after Theorem 1). For now, we only state that, under the blanket assumptions above, (PGA) achieves the following performance guarantee:

Proposition 1 (Convergence rate of projected gradient ascent). Suppose that (PGA) is run for n iterations with step-size γ. Then, the solution candidate x ≡ 1 n n k=1 x k enjoys the guarantee

V * -V(x) ≤ C 2 nγ + γdL 2 2 . ( 7
)
In particular, if (PGA) is run for n iterations with step-size

γ = C/L • √ 2/(nd), we get V * -V(x) ≤ LC 2d/n, (8)
i.e., x achieves an accuracy certificate V * -V(x) ≤ ε within at

most 2dL 2 C 2 /ε 2 = O(d/ε 2 ) iterations.
This is the main convergence guarantee for (PGA), so several remarks are in order. First, the step-size choice γ = C/L √ 2/(nd) has been calibrated so as to minimize the RHS of ( 7): as a result, this is the optimum parameter choice given a fixed running horizon. Hence, if we wish to guarantee an accuracy of ε in as few iterations as possible, Proposition 1 indicates that (PGA) should be run for n = 2dL2 C 2 /ε 2 iterations with step-size γ = ε/(dL 2 ).

A second issue concerns the use of averaging, i.e., taking x as a solution candidate instead of the algorithm's "last iterate", x n . The reason for this is that, without stronger assumptions, a large step-size could lead to overshooting the solution of (NUM), thus causing delays in convergence time; on the other hand, if the step-size is too small, the method will converge even slower. Averaging accelerates this convergence because the barycenter of iterates that oscillate around the solution set of (NUM) tends to be closer to a solution than any individual iterate.

B. Mirror ascent

Even though (PGA) is straightforward to implement and enjoys a clean performance guarantee, the number of iterations required to reach an accuracy certificate of ε scales with the number of sessions d (which could be of the order of millions in IP access domain networks with a few hundred nodes). Scale-up difficulties of this kind are encountered fairly often in large-scale optimization problems, and a promising way to circumvent them is to modify the algorithm's update mechanism so that it takes larger steps when the geometry of the feasible region is more constrained (meaning that there is less risk of overshooting). This is precisely the intuition behind the so-called mirror ascent (MA) method of Nemirovski and Yudin [START_REF] Nemirovsky | Problem Complexity and Method Efficiency in Optimization[END_REF] which works by replacing the Euclidean projector of (PGA) with a more geometry-aware update structure. Heuristically, instead of performing the update

x k+1 = Π Ω (x k + γv k ) = arg min x∈Ω { x k + γv k -x 2 2 } = arg min x∈Ω { γv k , x k -x + 1 2 x k -x 2 2 }, (9)
the mirror ascent algorithm employs the modified update

x k+1 = arg min x∈Ω { γv k , x k -x + D(x, x k )}, =: P Ω (x k , γv k ) (MA)
where the so-called "prox-mapping" P Ω is induced by the "Bregman divergence" term D, itself a generalization of the squared Euclidean norm x kx 2 2 . The starting point for the definition of the Bregman divergence D (and the associated prox-mapping P Ω ) is the notion of a distance-generating function (DGF), i.e., a function h : Ω → which is:

1) Continuous on Ω and continuously differentiable on the relative interior Ω • = {x ∈ d + : x s > 0, s x s < C} of Ω. 2 2) Strongly convex, i.e., there exists some σ > 0 such that

h(x ) ≥ h(x) + ∇h(x), x -x + σ 2 x -x 2 p ( 10
)
for all x, x ∈ Ω. 3The Bregman divergence associated to h is then defined as

D(p, x) = h(p) -h(x) -∇h(x), p -x , (11)
i.e., as the difference between h(p) and the best linear approximation of h(p) starting from x.

Remark. For clarity, we mention here that the above properties must be satisfied by the distance-generating function h, not the network utility maximization objective function V.

Before continuing, it is instructive to show how (PGA) can be recovered as a special case of (MA). Indeed, going back to Eq. ( 9), we see that the prox-mapping of (PGA) is simply

P Ω (x, w) = Π Ω (x + w), (12)
and is generated by the Bregman divergence

D(x, x ) = 1 2 x -x 2 2 , (13)
which is itself induced by the Euclidean distance-generating function

h(x) = 1 2 x 2 2 . (14)
The calculations needed to verify the above are straightforward, so we omit them; instead, we only note for posterity that the strong convexity modulus of the Euclidean regularizer is 1 relative to the Euclidean norm

• 2 .
With all this at hand, we are finally in a position to state and prove the main performance guarantee of the MA algorithm applied to (NUM):

Theorem 1 (Convergence rate of mirror ascent). Let h be a distance-generating function which is σ-strongly convex with respect to the p-norm on d . Then, if (MA) is run for n iterations with step-size γ, the solution candidate x ≡ 1 n n k=1 x k enjoys the guarantee

V * -V(x) ≤ Z nγ + γd 2-2/p L 2 2σ , (15)
with Z = D(x * , x 1 ) for some solution x * of (NUM). In particular, if (MA) is run for n iterations with step-size

γ = 1 L 2σ d 2-2/p Z n , (16)
we get

V * -V(x) ≤ L 2d 2-2/p σ Z n , (17)
i.e., x achieves an accuracy certificate V * -V(x) ≤ ε within at most 2d 2-2/p Z/σ • L 2 /ε 2 iterations.

This will be our main result for (MA) so, before discussing its proof, we provide some useful remarks: Remark 1. First, given that the MA algorithm subsumes the PGA algorithm, Theorem 1 subsumes in turn Proposition 1. Indeed, to obtain the bound [START_REF] Deng | Maximizing network utility of rechargeable sensor networks with spatiotemporally coupled constraints[END_REF] of Proposition 1, it suffices to recall that the Euclidean regularizer is 1-strongly convex with respect to the Euclidean norm (cf. the discussion above). It is also easy to see that Z = C 2 in this case, so the bound [START_REF] Deng | Maximizing network utility of rechargeable sensor networks with spatiotemporally coupled constraints[END_REF] follows from [START_REF] Nemirovsky | Problem Complexity and Method Efficiency in Optimization[END_REF] by setting p = 2, σ = 1 and Z = C 2 . Remark 2. Importantly, the dependence of the convergence rate of (MA) on the dimensionality of the problem is controlled by the coefficient d 2-2/p Z/σ. This observation is crucial for scalability: for the MA algorithm, the exponent of d in d 2-2/p must be minimized, i.e., we must focus on the case p = 1; concretely, this means that the most suitable distance-generating function for the network utility maximization problem should be examined relative to the Manhattan norm (p = 1). Therefore, the challenge is to construct a distance-generating function for which the ratio Z/σ grows as slowly as possible with d; this will be our main objective in the next section.

Proof of Theorem 1: By the definition of the prox-mapping and that of the Bregman divergence, (MA) satisfies the recursion

∇h(x k+1 ) -∇h(x k ) -γv k = 0. ( 18
)
Hence, by standard properties of the Bregman divergence [see e.g., 14, Lemma 2.1], we get

D k+1 -D(x * , x k ) ≤ γ x k -x * , v k + γ 2 2σ v k 2 q , (19)
where x * ∈ arg max V and q is such that 1/q + 1/p = 1. Then, letting D k = D(x * , x k ) and telescoping the above yields

n k=1 D k+1 -D k + γ x * -x k , v k ≤ γ 2 n k=1 v k 2 q 2σ ( 20
)
implying in turn that

D n+1 -D 1 + γ n k=1 x * -x k , v k ≤ γ 2 2σ n k=1 v k 2 q . (21) Then, since the Bregman divergence is nonnegative, we get n k=1 x * -x k , ∇V(x k ) ≤ D(x * , x 1 ) γ + γ n k=1 ∇V(x k ) 2 q 2σ . (22)
By concavity, we also have x * -x k , ∇V(x k ) ≥ V * -V(x k ); hence, by Jensen's inequality applied to x = n -1 n k=1 x k , Eq. ( 22) gives

V * -V(x) ≤ D(x * , x 1 ) γ + γd 2/q L 2 2σ , ( 23
)
where we used the definition of the q-norm and the standing assumption that |∂V/∂x s | ≤ L. The bound [START_REF] Nemirovsky | Problem Complexity and Method Efficiency in Optimization[END_REF] then follows by noting that 2/q = 2 -2/p (by definition). Finally, the value of the optimum step-size [START_REF] Palomar | A tutorial on decomposition methods for network utility maximization[END_REF] follows by minimizing the RHS of ( 15) with respect to γ; the rest of our claims then follow from elementary algebraic manipulations.

IV. Gradient Exponentiation

A. The exp-NUM algorithm

The analysis of the previous section shows that the most suitable distance-generating function for the NUM problem should be constructed relative to the Manhattan norm (since the value p = 1 minimizes the dependence of the coefficient d 2-2/p in [START_REF] Nemirovsky | Problem Complexity and Method Efficiency in Optimization[END_REF] with respect to d). In problems where the feasible region is the probability simplex ∆ = {z ∈ d+1 + : d s=0 z s = 1}, a popular choice is the (negative) entropy function

φ(z) = d s=0 z s log z s , (24)
which is known to be 1-strongly convex with respect to the Manhattan norm [19, p. 136]. Motivated by the above, we will focus in what follows on the modified entropy function

h(x) = s∈S x s log x s + x 0 log x 0 , (25)
where x 0 = C -d s=1 x s is a "slack variable" measuring the (Manhattan) distance between a state x ∈ Ω and the boundary face d s=1 x s of Ω. This distance-generating function will be the cornerstone of our analysis, so we begin by deriving the associated mirror ascent algorithm.

To that end, the first step is to obtain an explicit expression for the prox-mapping P Ω induced by the DGF (25). This is provided by the following lemma:

Lemma 1 (Entropic prox-mapping). With notation as above, the prox-mapping P Ω associated to (25) is

P Ω (x, w) = C (x s exp(w s )) s∈S C + s∈S x s [exp(w s ) -1] (26)
Proof: By definition, P Ω (x, w) solves the problem

P Ω (x, w) = arg min x ∈Ω { w, x -x + D(x , x)} (27)
with D(x , x) denoting the Bregman divergence

D(x , x) = s∈S x s log x s x s + x 0 log x 0 x 0 , (28)
as derived from definition [START_REF] Mo | Fair end-to-end window-based congestion control[END_REF]. To solve this problem, we will start with the ansatz that P Ω (x, w) is an interior point of Ω and treat (27) as an unconstrained problem; if the solution of the unconstrained problem ends up being an interior point of Ω, it will be a true solution of (27).

To proceed, let

Λ(x) = w, x -x + D(x , x) (29)
denote the objective function of (27). Then, by differentiating with respect to x , we obtain

∂Λ ∂x s = -w s + log x s x s + 1 -log x 0 x 0 -1, (30)
where we used the fact that ∂x 0 /∂x s = -1. Accordingly, setting these derivatives equal to zero, we get

x s x 0 = x s x 0 exp(w s ), (31)
and hence, summing over all sessions s ∈ S yields

C -x 0 x 0 = s x s exp(w s ) x 0 (32)
i.e., x 0 = C • x 0 /(x 0 + s x s exp(w s )). Our result then follows by substituting in (31) and backsolving.

With this lemma at hand, we obtain the exponentiated gradient ascent (EGA) algorithm

x k+1 s = C x k s exp(γ∂ s V(x k )) C + s∈S x k s [exp(γ∂ s V(x k )) -1] , (EGA)
which, in the present context, will be referred to as the exponentiated gradients algorithm for NUM (exp-NUM). For concreteness, a pseudocode implementation of exp-NUM can be found in Algorithm 1; for a schematic illustration of the differences between exp-NUM and the PGA algorithm of the previous section, see also Fig. 3. for k← 1 to n do 4:

v k ← ∇ V(x k ) 5: x k+1 ← C • x k exp(γv k )/(C -s x k s • (1 -exp(γv k s ))) 6:
end for

B. Convergence rate analysis of exp-NUM

We now turn to the convergence speed of the exp-NUM algorithm. Since exp-NUM is, first and foremost, a mirror ascent method, the analysis of the previous section shows that it suffices to calculate the following quantities: a) a lower bound for the strong convexity modulus of h relative to the 1 norm on Ω; b) an upper bound for Z = D(x * , x 1 ) relative to a convenient initialization of the algorithm. We begin below with the former: Lemma 2 (Strong convexity modulus). The entropic DGF (25) is (1/C)-strongly convex with respect to the 1 norm over Ω.

Proof: As stated above, the negative entropy φ of ( 24) is 1-strongly convex over the probability simplex ∆ = {x ∈ d+1 + : [19, p. 136]. To leverage this property, set z s = x s /C for s = 1, . . . , d, and z 0 = x 0 /C = 1 -d s=1 x s . Then, up to an additive constant that does not depend on x, we get

d s=0 x s = 1}
Cφ x C = C s∈S x s C • log x s C = s∈S x s log x s - s∈S x s log C = s∈S x s log x s + C - s∈S x s log C -C log C = h(x) -C • log C. ( 33
)
Since φ has strong convexity modulus 1, the function φ(x/C) will have modulus 1/C 2 . Thus,

h(x) = C • φ(x/C) + C log C is (1/C)-strongly convex over Ω.
As a last step, we need an initialization point which is both easy to implement and is as close as possible (in terms of the Bregman divergence) to all other points in Ω. In view of the above, the most natural initialization candidate is the barycenter of Ω, viz.

x c = (C/(d + 1), . . . , C/(d + 1)).

With this initialization, we obtain the following convergence speed estimate:

Theorem 2 (Convergence rate of exp-NUM). The output x of Algorithm 1 enjoys the guarantee

V * -V(x) ≤ C log(d + 1) nγ + γCL 2 2 . ( 35
)
In particular, if Algorithm 1 is run for n iterations with step-size γ = L -1 2 log(d + 1)/k, we get

V * -V(x) ≤ LC 2 log(d + 1)/n, (36)
i.e., x achieves an accuracy certificate V * -V(x) ≤ ε within at most 2 log(d + 1)

• L 2 C 2 /ε 2 = O(log d/ε 2 ) iterations.
Proof. As a first step, we will establish an upper bound for the quantity Z = D(x * , x 1 ) = D(x * , x c ). To that end, substituting in (28), we obtain for all x ∈ Ω:

D(x, x c ) =

C. Multipath scenario

We close this section by discussing the scenario where the traffic of each session can split into multiple paths. In this case, the network controller can route traffic over different paths in order to reduce the load on specific links. This multipath scenario can be formulated through the following optimization problem:

maximize s∈S U s         p∈P s x sp         -µ • ∈L g         s: ∈P s x sp         subject to s∈S p∈P s
x sp ≤ C (NUM-mult) The simplicial structure of (NUM-mult) pushes for an algorithm based on mirror ascent, for the same reasons as before. Concretely, to extend (25) to the multipath setting, consider the distance-generating function

h(x) = s,p x sp •log x sp + C - s,p x sp •log C - s,p x sp . (39)
Similarly to Lemma 1, we derive the exp-multiNUM proximal map P x (v) for the multipath scenario using the DGF above:

P Ω (x, w) = C • x exp(w) C -s,p x sp • (1 -exp(w sp )) , (40)
According to this algorithm, we are still able to obtain logarithmic dependence on the problem's dimensionality (considering the product of paths and sessions). Specifically, following the same reasoning as in the previous section, we get the convergence rate estimate

V(x * ) -V x ≤ C • L • 2 • log( s |P s | + 1) k . ( 41
)
The proof is similar to that of Theorem 2 and is omitted due to space limitations.

V. Numerical Experiments

In this section, we evaluate the performance of the proposed exp-NUM algorithm via extensive numerical simulations. Here, the main goals are to show the rate of convergence compared to other gradient-based methods, and to discuss how the problem dimensionality affects convergence and algorithm complexity per iteration. Our simulations consider a wide number of scenarios to cover complementary aspects of the NUM problem. We use two different α-fair models, namely α = 0 and α = 1, corresponding respectively to sum-rate maximization and proportional fairness. The underneath network topology used is based either on Erdős-Rényi or Barabási-Albert models. An exhaustive list of the parameters used in the simulations is provided in Table III.

A. Convergence rate

In Fig. 4 we show the main outcomes of our simulations. Specifically, we evaluate four scenarios by combining the two α-fair models (α = 0 and α = 1) with the aforementioned network models. For each scenario, we plot the objective value V(x) for the first 100.000 iterations (top), and the corresponding error (i.e., the difference with the optimal solution) in log-log scale (bottom). The plots compare exp-NUM with the PGA algorithm, with step-sizes computed according to Proposition 1 and Theorem 2 in order to reach a fixed accuracy threshold, common for both algorithms (to ensure fair testing).

The first immediate result that we observe from the plots is that exp-NUM approaches a solution of (NUM) in a relatively small number of iterations. While both algorithms eventually get close to the optimal objective value, the error for exp-NUM after 100.000 iterations is between one and two orders of magnitude smaller than gradient ascent. This observation is inline with our theoretical findings: by Theorem 2, the expected theoretical bound for exp-NUM is about 20 times tighter than that of PGA when d = 1000, and this difference only grows larger as the number of variables increases.

Another interesting observation is that the outcomes of the simulations do not show substantial differences based on the topology of the network. This is also confirmed by further analysis performed using the Watts-Strogatz model (but which we do not report due to space limitations). This is an important achievement for our algorithm which guarantees that its convergence rate is basically independent of the underlying network topology.

B. Problem dimensionality

The main reason to use first-order, gradient-based methods is their low per-iteration complexity, i.e., the fact that their updates remain relatively cheap in massively large systems. To make this precise, in Table IV, we compare the execution time per iteration of exp-NUM and PGA along with the Newton-based method of [START_REF] Athuraliya | Optimization flow control with newton-like algorithm[END_REF]. In the test machine used for our experiments, Netwon's method can be used efficiently only when d is fairly small: as the dimension of the problem increases, the need for the inversion of the Hessian matrix, which requires large processing and memory resources, makes the computation infeasible. In fact, even within the class of first-order methods, each iteration of exp-NUM is around 40% faster than a PGA update because exp-NUM does not require an projection step (which, in the case of PGA involves a sorting subroutine, and thus increases the periteration complexity by an O(log d) factor). We summarize all this in Table IV. VI. Conclusions In this paper, we examined the challenging scenario of resource allocation in large-scale networks with massive numbers of active sessions. In this setting, the dimensionality d of the system (number of sessions) represents a significant impediment to solving the network utility maximization rapidly and efficiently. To overcome this, we proposed the exp-NUM algorithm, a modified gradient algorithm that exploits the geometry of the problem to take bold (but safe) step-size decisions. We have shown that exp-NUM is effectively dimension-free and attains an ε-optimal solution within O(log d/ε 2 ) iterations -an exponential improvement over the poly(d) dimensionality dependence of first-order benchmark methods. In addition to these theoretical guarantees, we also provided a battery of simulation experiments which verify the scheme's superior scalability and fast convergence. We find this property particularly appealing as it provides a springboard for the application of similar techniques to a wide range of resource allocation problems in the area of networking.

An important question that remains is whether the algorithm can also be accelerated relative to the desired accuracy ε -e.g., to obtain an ε-optimal solution in no more than O(log d/ √ ε). We believe this should be possible by leveraging Nesterov's fast gradient method; we defer this analysis to the future.

Fig. 1 :

 1 Fig. 1: Network with three sessions sharing two links.

Fig. 2 :

 2 Fig. 2: Feasible region (green), global capacity constraint (orange) and solution profile of Example 1.

2 )

 2 We write x p = d s=1 |x s | p 1/p for the p-norm of x. When p = 2, this is the standard Euclidean norm; we will also make use of the Manhattan norm x 1 = d s=1 |x s | (corresponding to p = 1), and the sup norm x ∞ = max s=1,...,d |x s | (corresponding to p = ∞).

Algorithm 1 :

 1 exponentiated gradients algorithm for NUM Require: step-size γ, number of iterations n 1: function Main 2: x 1 ← {C/(d + 1), C/(d + 1), . . . , C/(d + 1)} 3:

Fig. 3 :

 3 Fig. 3: Projected vs. exponentiated gradient ascent.

  x s -C log C + C log(d + 1). (37) Therefore, since the maximum of the expression d s=0 x s log x s over Ω is obtained when x s = C for some s ∈ {0, . . . , d} and x s = 0 for all s s, we readily get D(x, x c ) ≤ C log C -C log C +C log(d+1) = C log(d+1). (38) This shows that D(x * , x c ) ≤ log(d + 1). The bound (35) then follows by applying Theorem 1 and substituting the parameter values Z = C log(d + 1), p = 1 and σ = 1/C in the estimated bound (15).

Fig. 4 :

 4 Fig. 4: Value convergence (in log-log scale) for exp-NUM and PGA in Erdős-Rényi and Barabási-Albert networks.

TABLE II :

 II Notation used in the paper.

		Network model
	N	Set of nodes
	L	Set of links
	S	Set of sessions
	d	Number of sessions (dimensionality)
	c	Capacity of link ∈ L
	C	Global cap
		Network utility maximization
	x	Session rates (control variable)
	U s	Utility function of session s ∈ S
	g	Penalty function for link ∈ L
	µ	Rigidity coefficient
	V	Objective function
	v	Gradient of V
	L	derivative bound
		Mirror ascent
	h	Distance-generating function
	D	Bregman divergence
	γ	Step size
	k	Iteration counter
	n	Total number of iterations
	P	Prox-mapping

TABLE III :

 III Default parameters used in the simulations.

	d	1000	|L|	100
	c l	N (0.5, 0.1)	C	max l c l
	w s	N (0.5, 0.1)	µ	0.1

TABLE IV :

 IV Execution time per iteration of different iterative algorithms.

	Algorithm	d = 100	d = 1000	d = 10.000
	exp-NUM	5,0×10 -4 s	6,9×10 -3 s	7,3×10 -2 s
	GD	7,0×10 -4 s	9,5×10 -3 s	1,1×10 -1 s
	Newton	1,1×10 -2 s	2,1×10 0 s	∞

The differentiability assumption could be lifted at the cost of using supergradients instead of gradients, but we do not do so for simplicity.

Strictly speaking, we will tacitly assume that h is continuously differentiable on the relative interior of every face of Ω, but this technicality can be safely ignored in what follows.

In the above, the value of p in the p-norm • p is a tunable parameter whose role will be discussed below.

*The work was done when the author was with Huawei Technologies. This work was partially supported by the French National Research Agency (ANR) project ORACLESS, and the Huawei Flagship project ULTRON. Algorithm