Contemporary cloud data centers (DCs) are proliferating to provide major Internet services such as video distribution, and data backup . To fulfill the goal of reducing electricity cost and improving sustainability, advanced smart grid features are nowadays adopted by cloud networks . Due to their inherent stochasticity, these features challenge the overall network design, and particularly the network’s resource allocation task. Recent approaches on this topic mitigate the spatio-temporal uncertainty associated with energy prices and renewable availability [4, 5, 6, 7, 8]. Interestingly, the resultant algorithms reconfirm the critical role of the dual decomposition framework, and its renewed popularity for optimizing modern stochastic networks.
However, slow convergence and the associated network delay of existing resource allocation schemes have recently motivated improved first- and second-order optimization algorithms [9, 10, 11, 12]. Yet, historical data have not been exploited to mitigate future uncertainty. This is the key novelty of the present paper, which permeates benefits of statistical learning to stochastic resource allocation tasks.
Targeting this goal, a critical observation is that renowned network optimization algorithms (e.g., back-pressure and max-weight) are intimately linked with Lagrange’s dual theory, where the associated multipliers admit pertinent price interpretations [13, 14, 15]. We contend that learning these multipliers can benefit significantly from historical relationships and trends present in massive datasets . In this context, we revisit the stochastic network optimization problem from a machine learning vantage point, with the goal of learning the Lagrange multipliers in a fast and efficient manner. Works in this direction include  and , but the methods there are more suitable for problems with two features: 1) the network states belong to a distribution with finite support; and, 2) the feasible set is discrete with a finite number of actions. Without these two features, the involved learning procedure may become intractable, and the advantageous performance guarantee may not hold. Overall, online resource allocation capitalizing on data-driven learning schemes remains a largely uncharted territory.
Motivated by recent advances in machine learning, we systematically formulate the resource allocation problem as an online task with batch training and operational learning-and-adapting phases. In the batch training phase, we view the empirical dual problem of maximizing the sum of finite concave functions, as an empirical risk maximization (ERM) task, which is well-studied in machine learning 
. Leveraging our specific ERM problem structure, we modify the recently developed stochastic average gradient approach (SAGA) to fit our training setup. SAGA belongs to the class of fast incremental gradient (a.k.a. stochastic variance reduction) methods, which combine the merits of both stochastic gradient, and batch gradient methods. In the resource allocation setup, our offline SAGA yields empirical Lagrange multipliers with order-optimal linear convergence rate as batch gradient approach, and per-iteration complexity as low as stochastic gradient approach. Broadening the static learning setup in  and , we further introduce a dynamic resource allocation approach (that we term online SAGA) that operates in a learn-and-adapt fashion. Online SAGA fuses the benefits of stochastic approximation and statistical learning: In the learning mode, it preserves the simplicity of offline SAGA to dynamically learn from streaming data thus lowering the training error; while it also adapts by incorporating attributes of the well-appreciated stochastic dual (sub)gradient approach (SDG) [4, 5, 6, 7, 21], in order to track queue variations and thus guarantee long-term queue stability.
In a nutshell, the main contributions of this paper can be summarized as follows.
Using stochastic network management as a motivating application domain, we take a fresh look at dual solvers of constrained optimization problems as machine learning iterations involving training and operational phases.
During the training phase, we considerably broaden SAGA to efficiently compute Lagrange multiplier iterates at order-optimal linear convergence, and computational cost comparable to the stochastic gradient approach.
Outline. The rest of the paper is organized as follows. The system models are described in Section II. The motivating resource allocation setup is formulated in Section III. Section IV deals with our learn-and-adapt dual solvers for constrained optimization. Convergence analysis of the novel online SAGA is carried out in Section V. Numerical tests are provided in Section VI, followed by concluding remarks in Section VII.
Notation. denotes expectation (probability);
denotes the all-one vector; anddenotes the -norm of vector . Inequalities for vectors , are defined entry wise; ; and stands for transposition. denotes big order of , i.e., as ; and denotes small order of , i.e., as .
Ii Network Modeling Preliminaries
This section focuses on resource allocation over a sustainable DC network with mapping nodes (MNs), and DCs. MNs here can be authoritative DNS servers as used by Akamai and most content delivery networks, or HTTP ingress proxies as used by Google and Yahoo! , which collect user requests over a geographical area (e.g., a city or a state), and then forward workloads to one or more DCs distributed across a large area (e.g., a country).
Notice though, that the algorithms and their performance analysis in Sections V and IV can be applied to general resource allocation tasks, such as energy management in power systems , cross-layer rate-power allocation in communication links , and traffic control in transportation networks .
Network constraints. Suppose that interactive workloads are allocated as in e.g., , and only delay-tolerant workloads that are deferrable are to be scheduled across slots. Typical examples include system updates and data backup, which provide ample optimization opportunities for workload allocation based on the dynamic variation of energy prices, and the random availability of renewable energy sources.
Suppose that the time is indexed in discrete slots , and let denote an infinite time horizon. With reference to Fig. 1, let denote the amount of delay-tolerant workload arriving at MN in slot , and the vector collecting the workloads routed from MN to all DCs in slot . With the fraction of unserved workload buffered in corresponding queues, the queue length at each mapping node at the beginning of slot , obeys the recursion
where denotes the set of DCs that MN is connected to.
At the DC side, denotes the workload processed by DC during slot . Unserved portions of the workloads are buffered at DC queues, whose length obeys (cf. (1))
where is the queue length in DC at the beginning of slot , and denotes the set of MNs that DC is linked with. The per-slot workload is bounded by the capacity for each DC ; that is,
Since the MN-to-DC link has bounded bandwidth , it holds that
If MN and DC are not connected, then .
Operation costs. We will use the strictly convex function to denote the cost for distributing workloads from MN to DC , which depends on the distance between them.
Let denote the energy provided at the beginning of slot by the renewable generator at DC , which is constrained to lie in . The power consumption by DC is taken to be a quadratic function of the demand ; that is , where is a time-varying parameter capturing environmental factors, such as humidity and temperature . The energy transaction cost is modeled as a linear function of the power imbalance amount to capture the cost in real time power balancing; that is, , where is the buy/sell price in the local power wholesale market.111When the prices of buying and selling are not the same, the transaction cost is a piece-wise linear function of the power imbalance amount as in . Clearly, each DC should buy energy from external energy markets in slot at price if , or, sell energy to the markets with the same price, if .
The aggregate cost for the considered MN-DC network per slot is a convex function, given by
When (5) is assumed under specific operating conditions, our analysis applies to general smooth and strong-convex costs.
Iii Motivating Application: Resource Allocation
For notational brevity, collect all random variables at timein the state vector ; the optimization variables in the vector and let . Considering the scheduled as the outgoing amount of workload from DC , define the “node-incidence” matrix with entry given by
We assume that each row of has at least one entry, and each column of has at most one
entry, meaning that each node has at least one outgoing link, and each link has at most one source node. Further, collect the instantaneous workloads across MNs in the zero-paddedvector , as well as the link and DC capacities in the vector .
Hence, the sought scheduling is the solution of the following long-term network-optimization problem
where (6b) and (6c) come from concatenating (3)-(4) and (1)-(2), while (6d) ensures strong queue stability as in [21, Definition 2.7]. The objective in (6a) considers the entire time horizon, and the expectation is over all sources of randomness, namely the random vector , and the randomness of the variables and induced by the sample path of . For the problem (6), the queue dynamics in (6c) couple the optimization variables over the infinite time horizon. For practical cases where the knowledge of is causal, finding the optimal solution is generally intractable. Our approach to circumventing this obstacle is to replace (6c) with limiting average constraints, and employ dual decomposition to separate the solution across time, as elaborated next.
Iii-a Problem reformulation
Indeed, (6c) implies that after summing over and taking expectations yields . Since both and are bounded under (6d), dividing both sides by and taking , yields (7). Using (7), we can write a relaxed version of (6) as
Compared to (6), the queue variables are not present in (8), while the time-coupling constraints (6c) and (6d) are replaced with (7). As (8) is a relaxed version of (6), it follows that . Hence, if one solves (8) instead of (6), it will be prudent to derive an optimality bound on , provided that the schedule obtained by solving (8) is feasible for (6). In addition, using arguments similar to those in  and , it can be shown that if the random process is stationary, there exists a stationary control policy , which is only a function of the current ; it satisfies (6b) almost surely; and guarantees that , and . This implies that the infinite time horizon problem (8) is equivalent to the following per slot convex program
where we interchangeably used and , to emphasize the dependence of the real-time cost and decision on the random state . Note that the optimization in (9) is w.r.t. the stationary policy . Hence, there is an infinite number of variables in the primal domain. Observe though, that there is a finite number of constraints coupling the realizations (cf. (9b)). Thus, the dual problem contains a finite number of variables, hinting that the problem is perhaps tractable in the dual space [24, 12]. Furthermore, we will demonstrate in Section V that after a careful algorithmic design in Section IV, our online solution for the relaxed problem (9) is also feasible for (6).
Iii-B Lagrange dual and optimal solutions
Correspondingly, the dual problem of (8) is
where . We henceforth refer to (12) as the ensemble dual problem. Note that similar to , here and are both parameterized by .
If the optimal Lagrange multipliers were known, a sufficient condition for the optimal solution of (8) or (9) would be to minimize the Lagrangian or its instantaneous versions over the set [26, Proposition 3.3.4]. Specifically, as formalized in the ensuing proposition, the optimal routing and workload schedules in this MN-DC network can be expressed as a function of associated with (9b), and the realization of the random state .
|and the optimal instantaneous workload scheduling decisions are given by [with ]|
where and denote the inverse functions of and , respectively.
We omit the proof of Proposition 1, which can be easily derived using the KKT conditions for constrained optimization . Building on Proposition 1, it is interesting to observe that the stationary policy we are looking for in Section III-A is expressible uniquely in terms of . The intuition behind this solution is that the Lagrange multipliers act as interfaces between MN-DC and workload-power balance, capturing the availability of resources and utility information which is relevant from a resource allocation point of view.
However, to implement the optimal resource allocation in (13), the optimal multipliers must be known. To this end, we first outline the celebrated stochastic approximation-based and corresponding Lyapunov approaches to stochastic network optimization. Subsequently, we develop a novel approach in Section IV to learn the optimal multipliers in both offline and online settings.
Iii-C Stochastic dual (sub-)gradient ascent
For the ensemble dual problem (12), a standard (sub)gradient iteration involves taking the expectation over the distribution of to compute the gradient . This is challenging because the underlying distribution of
is usually unknown in practice. Even if the joint probability distribution functions were available, finding the expectations can be non-trivial especially in high-dimensional settings (and/or ).
where the stochastic (instantaneous) gradient
is an unbiased estimator of the ensemble gradient given by. The primal variables can be found by solving “on-the-fly” the instantaneous problems, one per slot
where the operator accounts for cases that the Lagrangian has multiplier minimizers. The minimization in (15) is not difficult to solve. For a number of relevant costs and utility functions in Proposition 1, closed-form solutions are available for the primal variables. Note from (14) that the iterate depends only on the probability distribution of through the stochastic gradient . Consequently, the process is Markov with transition probability that is time invariant since is stationary. In the context of Lyapunov optimization [5, 13] and , the Markovian iteration in (14) is interpreted as a virtual queue recursion; i.e., .
Thanks to their low complexity and ability to cope with non-stationary scenarios, SDG-based approaches are widely used in various research disciplines; e.g., adaptive signal processing , stochastic network optimization [21, 13, 17, 29], and energy management in power grids [5, 23]. Unfortunately, SDG iterates are known to converge slowly. Although simple, SDG does not exploit the often sizable number of historical samples. These considerations motivate a systematic design of an offline-aided-online approach, which can significantly improve online performance of SDG for constrained optimization, and can have major impact in e.g., network resource allocation tasks by utilizing streaming big data, while preserving its low complexity and fast adaptation.
Iv Learn-and-adapt Resource allocation
Before developing such a promising approach that we view as learn-and-adapt SDG scheme, we list our assumptions that are satisfied in typical network resource allocation problems.
(as1) State process is independent and identically distributed (i.i.d.), and the common probability density function (pdf) has bounded support.
is independent and identically distributed (i.i.d.), and the common probability density function (pdf) has bounded support.
(as2) The cost in (5) is non-decreasing w.r.t. ; it is a -strongly convex function222We say that a function is -strongly convex if and only if is convex for all , where .; and its gradient is Lipschitz continuous with constant for all .
(as3) There exists a stationary policy satisfying , for all , and , where is the slack constant; and
(as4) The instantaneous dual function in (12) is -strongly concave, and its gradient is -Lipschitz continuous with condition number , for all .
Although (as1) can be relaxed if ergodicity holds, it is typically adopted by stochastic resource allocation schemes for simplicity in exposition [13, 17, 30]. Under (as2), the objective function is non-decreasing and strongly convex, which is satisfied in practice with quadratic/exponential utility or cost functions . The so-called Slater’s condition in (as3) ensures the existence of a bounded Lagrange multiplier , which is necessary for the queue stability of (6); see e.g., [31, 17]. If (as3) cannot be satisfied, one should consider reducing the workload arrival rates at the MN side, or, increasing the link and facility capacities at the DC side. The -Lipschitz continuity of in (as4) directly follows from the strong-convexity of the primal function in (as2) with , where is the spectral radius of the matrix . The strong concavity in (as4) is frequently assumed in network optimization , and it is closely related to the local smoothness and the uniqueness of Lagrange multipliers assumed in [13, 17, 30]. For pessimistic cases, (as4) can be satisfied by subtracting an -regularizer from the dual function (11
), which is typically used in machine learning applications (e.g., ridge regression). We quantify its sub-optimality in Appendix A. Note that Appendix A implies that the primal solutionwill be -optimal and feasible for the regularizer . Since we are after an -optimal online solution, it suffices to set .
Iv-a Batch learning via offline SAGA based training
Consider that a training set of historical state samples is available. Using , we can find an empirical version of (11) via sample averaging as
Note that has been replaced by to differentiate training (based on historical data) from operational (a.k.a. testing or tracking) phases. Consequently, the empirical dual problem can be expressed as
Recognizing that the objective is a sum of finite concave functions, the task in (17) in the machine learning parlance is termed empirical risk maximization (ERM) , which is carried out using the batch gradient ascent iteration
where the index represents the batch learning (iteration) index, and is the stepsize that controls the learning rate. While iteration (18) exhibits a decent convergence rate, its computational complexity will be prohibitively high as the data size grows large. A typical alternative is to employ a stochastic gradient (SG) iteration, which uniformly at random selects one of the summands in (18). However, such an SG iteration relies only on a single unbiased gradient correction, which leads to a sub-linear convergence rate. Hybrids of stochastic with batch gradient methods are popular subjects recently [19, 32].333Stochastic iterations for the empirical dual problem are different from that in Section III-C, since stochasticity is introduced by the randomized algorithm itself, in oppose to the stochasticity of future states in the online setting.
Leveraging our special problem structure, we will adapt the recently developed stochastic average gradient approach (SAGA) to fit our dual space setup, with the goal of efficiently computing empirical Lagrange multipliers. Compared with the original SAGA that is developed for unconstrained optimization , here we start from the constrained optimization problem (6), and derive first a projected form of SAGA.
Per iteration , offline SAGA first evaluates at the current iterate , one gradient sample with sample index selected uniformly at random. Thus, the computational complexity of SAGA is that of SG, and markedly less than the batch gradient ascent (18), which requires such evaluations. Unlike SG however, SAGA stores a collection of the most recent gradients for all samples , with the auxiliary iteration index denoting the most recent past iteration that sample was randomly drawn; i.e., . Specifically, SAGA’s gradient combines linearly the gradient randomly selected at iteration with the stored ones to update the multipliers. The resultant gradient is the sum of the difference between the fresh gradient and the stored one at the same sample, as well as the average of all gradients in the memory, namely
|Therefore, the update of the offline SAGA can be written as|
where denotes the stepsize. The steps of offline SAGA are summarized in Algorithm 1.
To expose the merits of SAGA, recognize first that since is drawn uniformly at random from set , we have that , and thus the expectation of the corresponding gradient sample is given by
Hence, is an unbiased estimator of the empirical gradient in (18). Likewise, , which implies that the last two terms in (19a) disappear when taking the mean w.r.t. ; and thus, SAGA’s stochastic averaging gradient estimator is unbiased, as is the case with SG that only employs .
With regards to variance, SG’s gradient estimator has , which can be scaled down using decreasing stepsizes (e.g., ), to effect convergence of iterates in the mean-square sense . As can be seen from the correction term in (19a), SAGA’s gradient estimator has lower variance than . Indeed, the sum term of in (19a) is deterministic, and thus it has no effect on the variance. However, representing gradients of the same drawn sample , the first two terms are highly correlated, and their difference has variance considerably smaller than . More importantly, the variance of stochastic gradient approximation vanishes as approaches the optimal argument for (17); see e.g., . This is in oppose to SG where the variance of stochastic approximation remains even if the iterates are close to the optimal solution. This variance reduction property allows SAGA to achieve a linear convergence with constant stepsizes, which is not achievable for the SG method.
In the following theorem, we use the result in  to show that the offline SAGA method is linearly convergent.
Consider the offline SAGA iteration in (19), and assume that the conditions in (as2)-(as4) are satisfied. If denotes the unique optimal argument in (17), and the stepsize is chosen as with Lipschitz constant as in (as4), then SAGA iterates initialized with satisfy
|where with denoting the condition number in (as4), the expectation is over all choices of the sample index up to iteration , and constant is|
See Appendix -B. ∎
Since , Theorem 1 asserts that the sequence generated by SAGA converges exponentially to the empirical optimum in mean-square. Similar to in (as4), if in (16) are -smooth, then so is , and the sub-optimality gap induced by can be bounded by
Remarkably, our offline SAGA based on training data is able to obtain the order-optimal convergence rate among all first-order approaches at the cost of only a single gradient evaluation per iteration. We illustrate the convergence of offline SAGA with two SG variants in Fig. 2 for a simple example. As we observe, SAGA converges to the optimal argument linearly, while the SG method with diminishing stepsize has a sublinear convergence rate and the one with constant stepsize converges to a neighborhood of the optimal solution.
The stepsize in Theorem 1 requires knowing the Lipschitz constant of (16). In our context, it is given by , with denoting the spectral radius of the matrix , and defined in (as2). When can be properly approximated in practice, it is worth mentioning that the linear convergence rate of offline SAGA in Theorem 1 can be established under a wider range of stepsizes, with slightly different constants and [19, Theorem 1].
Iv-B Learn-and-adapt via online SAGA
Offline SAGA is clearly suited for learning from (even massive) training datasets available in batch format, and it is tempting to carry over our training-based iterations to the operational phase as well. However, different from common machine learning tasks, the training-to-testing procedure is no longer applicable here, since the online algorithm must also track system dynamics to ensure the stability of queue lengths in our network setup. Motivated by this observation, we introduce a novel data-driven approach called online SAGA, which incorporates the benefits of batch training to mitigate online stochasticity, while preserving the adaptation capability.
Unlike the offline SAGA iteration that relies on a fixed training set, here we deal with a dynamic learning task where the training set grows in each time slot . Running sufficiently many SAGA iterations for each empirical dual problem per slot could be computationally expensive. Inspired by the idea of dynamic SAGA with a fixed training set , the crux of our online SAGA in the operational phase is to learn from the dynamically growing training set at an affordable computational cost. This allows us to obtain a reasonably accurate solution at an affordable computational cost - only a few SAGA iterations.
To start, it is useful to recall a learning performance metric that uniformly bounds the difference between the empirical loss in (16) and the ensemble loss in (11) with high probability (w.h.p.), namely 
where denotes a bound on the statistical error induced by the finite size of the training set
. Under proper (so-termed mixing) conditions, the law of large numbers guarantees thatis in the order of , or, for specific function classes [16, Section 3.4]. That the statement in (23) holds with w.h.p. means that there exists a constant such that (23) holds with probability at least . In such case, the statistical accuracy depends on , but we keep that dependency implicit to simplify notation [16, 33].
Let be an -optimal solution obtained by the offline SAGA over the training samples, where will be termed optimization error emerging due to e.g., finite iterations say on average per sample; that is, , where is the optimal argument in (17). Clearly, the overall learning error, defined as the difference between the empirical loss with and the ensemble loss with , is bounded above w.h.p. as
where the bound superimposes the optimization error with the statistical error. If is relatively small, will be large, and keeping the per sample iterations small is reasonable for reducing complexity, but also because stays comparable to . With the dynamically growing training set in the operational phase, our online SAGA will aim at a “sweet-spot” between affordable complexity (controlled by ), and desirable overall learning error that is bounded by with .
The proposed online SAGA method consists of two complementary stages: offline learning and online learning and adaptation. In the offline training-based learning, it runs SAGA iteration (19) (cf. Algorithm 1) on a batch dataset with historical samples, and on average iterations per sample. In the online operational phase, initialized with the offline output ,444To differentiate offline and online iterations, let denote -th iteration in the offline learning, and denote -th iteration during slot . online SAGA continues the learning process by storing a “hot” gradient collection as in Algorithm 1, and also adapts to the queue length that plays the role of an instantaneous multiplier estimate analogous to in (14). Specifically, online SAGA keeps acquiring data and growing the training set based on which it learns online by running iterations (19) per slot . The last iterate of slot initializes the first one of slot ; that is, . The learned captures a priori information . To account for instantaneous state information as well as actual constraint violations which are valuable for adaptation, online SAGA will superimpose to the queue length (multiplier estimate) . To avoid possible bias in the steady-state, this superposition is bias-corrected, leading to an effective multiplier estimate (with short-hand notation )
where scalar tunes emphasis on past versus present state information, and the value of constant will be specified in Theorem 4. Based on the effective dual variable, the primal variables are obtained as in (13) with replaced by ; or, for general constrained optimization problems, as
The necessity of combining statistical learning and stochastic approximation in (25) can be further understood from obeying feasibility of the original problem (6); e.g., constraints (6c) and (6d). Note that in the adaptation step of online SAGA, the variable is updated according to (6c), and the primal variable is obtained as the minimizer of the Lagrangian using the effective multiplier , which accounts for the queue variable via (25). Assuming that is approaching infinity, the variable becomes large, and, therefore, the primal allocation obtained by minimizing online Lagrangian (26) with will make the penalty term negative. Therefore, it will lead to decreasing through the recursion in (6c). This mechanism ensures that the constraint (6d) will not be violated.555Without (25), for a finite time , is always near-optimal for the ensemble problem (12), and the primal variable in turn is near-feasible w.r.t. (9b) that is necessary for (6d). Since essentially accumulates online constraint violations of (9b), it will grow linearly with and eventually become unbounded. The online SAGA is listed in Algorithm 2.
Remark 2 (Connection with existing algorithms).
Online SAGA has similarities and differences with common schemes in machine learning, and stochastic network optimization.
(P1) SDG in [24, 6] can be viewed as a special case of online SAGA, which purely relies on stochastic estimates of multipliers or instantaneous queue lengths, meaning . In contrast, online SAGA further learns by stochastic averaging over (possibly big) data to improve convergence of SDG.
(P2) Several schemes have been developed for statistical learning via large-scale optimization, including SG, SAG , SAGA , dynaSAGA , and AdaNewton . However, these are generally tailored for static large-scale optimization, and are not suitable for dynamic resource allocation, because having neither tracks constraint violations nor adapts to state changes, which is critical to ensure queue stability (6d).
(P3) While online SAGA shares the learning-while-adapting attribute with the methods in  and , these approaches generally require using histogram-based pdf estimation, and solve the resultant empirical dual problem per slot exactly. However, histogram-based estimation is not tractable for our high-dimensional continuous random vector; and when the dual problem cannot be written in an explicit form, solving it exactly per iteration is computationally expensive.
Remark 3 (Computational complexity).
Considering the dimension of as , the dual update in (14) has the computational complexity of , since each row of has at most non-zero entries. As the dimension of is , the primal update in (15), in general, requires solving a convex program with the worst-case complexity of using a standard interior-point solver. Overall, the per-slot computational complexity of the Lyapunov optimization or SDG in [21, 24] is .
Online SAGA has a comparable complexity with that of SDG during online adaptation phase, where the effective multiplier update, primal update and queue update also incur a worst-case complexity in total. Although the online learning phase is composed of additional updates of the empirical multiplier through (19), we only need to compute one gradient at a randomly selected sample . Hence, per-iteration complexity is determined by that of solving (15), and it is independent of the data size . Clearly, subtracting stored gradient only incurs computations, and the last summand in (19a) can be iteratively computed using running averages. In addition, for each new per slot, finding the initial gradient requires another