Asynchronous Incremental Stochastic Dual Descent Algorithm for Network Resource Allocation
Amrit S. Bedi Ketan Rajawat

	
	
	
	
	
	
	
	
	
	
	
Asynchronous Incremental Stochastic  Dual Descent Algorithm for Network Resource Allocation
    Amrit Singh Bedi, Student Member, IEEE and Ketan Rajawat, Member, IEEE The	authors are with the Department of Electrical Engineering, IIT Kanpur, Kanpur
			(UP), India 208016 (email: amritbd, ketan@iitk.ac.in).

    
=======================================================================================================================================================================================================================


	
	
Stochastic network optimization problems entail finding resource allocation policies that are optimum on an average but must be designed in an online fashion. Such problems are ubiquitous in communication networks, where resources such as energy and bandwidth are divided among nodes to satisfy certain long-term objectives. This paper proposes an asynchronous incremental dual decent resource allocation algorithm that utilizes delayed stochastic gradients for carrying out its updates. The proposed algorithm is well-suited to heterogeneous networks as it allows the computationally-challenged or energy-starved nodes to, at times,  postpone the updates. The asymptotic analysis of the proposed algorithm is carried out, establishing dual convergence under both, constant and diminishing step sizes. It is also shown that with constant step size, the proposed resource allocation policy is asymptotically near-optimal. An application involving multi-cell coordinated beamforming is detailed, demonstrating the usefulness of the proposed algorithm. 
	
	
	
		Stochastic subgradient, resource allocation, asynchronous algorithm, incremental algorithm. 
 
	
	



§ INTRODUCTION

The recent years have witnessed an unprecedented growth in the complexity and bandwidth requirements of network services. The resulting stress on the network infrastructure has motivated the network designers to move away from simpler or modular architectures and towards optimum ones. To make sure that resources such as bandwidth and energy are allocated efficiently, optimum designs advocate cooperation between the network nodes <cit.>. This paper considers the problem of cooperative network resource allocation that arises in wireless communication networks <cit.>, smart grid systems <cit.>, and in the context of scheduling <cit.>. Of particular interest is the stochastic resource allocation problem, where the goal is to find an allocation policy that is asymptotically optimal <cit.>. Although such problems are infinite dimensional in nature, they can be solved in an online fashion via stochastic dual descent methods, allowing real-time resource allocation that is also asymptotically near-optimal <cit.>. 

Heterogeneous networks are common to a number of applications where the energy availability, computational capability, or the mode of operation of the nodes is not the same across the network. Key requirements for heterogeneous network protocols include scalability, robustness, and tolerance to delays and packet losses. Towards this end, a number of distributed algorithms have been proposed in the literature <cit.>. By eliminating the need for a fusion center, the distributed algorithms operate with reduced communication overhead, and render the network resilient to single-point failures. 

Most distributed algorithms still place stringent communication and computational requirements on the network nodes. For instance, the dual stochastic gradient methods entail multiple updates and message exchanges per time slot, and cannot handle missed or delayed updates. In heterogeneous networks, such delays are often unavoidable, arising due to poor channel conditions, traffic congestion, or limited processing power at certain nodes. This paper proposes a distributed asynchronous stochastic resource allocation algorithm that tolerates such delays. The next subsection outlines the main contributions of this paper.
 




 §.§ Contributions and organization

The stochastic resource allocation problem is formulated as a constrained optimization problem where the goal is to maximize a network-wide utility function. The allocated resources at the different nodes in the network are coupled through constraint functions that involve expectations with respect to a random network state. Specifically, the aim is to find an allocation policy that satisfies the constraints on an average. The distribution of the state variables is not known, so that the optimization problem does not admit an offline solution. Instead, the idea is to observe the instances of the state variables over time, and allocate resources in an online manner. It is well-known that stochastic dual descent algorithms yield viable online algorithms for such problems <cit.>. 

Within the heterogeneous network setting considered here, the focus is on distributed algorithms that can tolerate communication and processing delays <cit.>. Different from the state-of-the-art algorithms that utilize the standard stochastic  gradient methods <cit.>, we develop two variants of the asynchronous dual descent algorithm that allow some of the nodes in the network to temporarily “fall back,” in the event of low energy availability, unusually large processing delay, node shutdown, or channel impairments. The first asynchronous variant utilizes a fusion center to collect the possibly delayed  gradients from various nodes and carry out the updates (cf. <ref> ). The second variant eliminates the need for the fusion center, and instead utilizes the fully distributed and incremental stochastic  gradient descent algorithm, where the nodes carry out updates in a round-robin fashion and pass messages along a cycle (cf. <ref> ). As earlier, the use of stale  gradients for primal and dual updates, allows the algorithm to be run on two different clocks, one corresponding to the local resource allocation and tuned to the changing random network state, while the other dictated by the message passing protocol. The key feature of the proposed algorithm is the possibility for the second clock to slow down temporarily and wait for slower nodes to catch up. The proposed algorithm thus allows timely resource allocation, while tolerating occasional delays in message passing.

The asymptotic performance of the proposed algorithm is studied under certain regularity conditions on the problem structure and bounded delays. In particular, the asymptotic performance of the asynchronous incremental stochastic  subgradient descent (AIS-SD) algorithm is characterized under both, diminishing and constant step-sizes. The overall structure of the proof is based on the convergence results in the incremental stochastic  subgradient descent algorithm of <cit.>  and the asynchronous incremental subgradient method of <cit.>. Specific to the resource allocation problem at hand, the asymptotic near-optimality and almost sure feasibility of the primal allocation policy is established for the case of constant step sizes. When applied to resource allocation problems, the proposed algorithm is called asynchronous incremental stochastic dual descent (AIS-DD). It is remarked that since the proposed algorithms utilize stochastic  subgradient descent, their computational complexity is also comparable to other distributed stochastic algorithms <cit.>. The calculation of the  subgradient is the most computationally expensive step, and like other first-order algorithms, must be carried out at every time slot.

Finally, the stochastic coordinated multi-cell beamforming problem is formulated and solved via the proposed algorithm. Detailed simulations are carried out to demonstrate the usefulness of the proposed algorithm in delay-prone and distributed environments. Summarizing, the main contributions of the paper include (a) the AIS-SD algorithm and its convergence (b) primal near-optimality and feasibility results for the allocated resources using AIS-DD; and (c) demonstration of the proposed algorithm on a practical stochastic coordinated multi-cell beamforming problem. 
	 
	The rest of the paper is organized as follows. Sec. <ref> provides an outline of the related literature. Sec. <ref> describes the problem formulation and recapitulates the known results. Sec. <ref> details  the proposed algorithm. Sec. <ref> lists the required assumptions, and provides the primal and dual convergence results. Sec. <ref> formulates the stochastic version of the coordinated beamforming problem along with the relevant simulation results. Finally, Sec. <ref> concludes the paper. 

	


 §.§ Related work
 
Resource allocation problems have been well-studied in the context of cross-layer optimization in networks <cit.>. Popular tools for solving stochastic resource allocation problems include the backpressure algorithm <cit.> and variants of the stochastic dual descent method <cit.>. However, most of these works only consider synchronous algorithms, and the effect of communication delays has not been examined in detail. An exception is the asynchronous subgradient method proposed in <cit.>, where delayed subgradients were utilized for resource allocation. The present work extends the algorithm in <cit.> by allowing delayed stochastic subgradients. Additionally, the proposed algorithm is also incremental, and is therefore applicable to a wider variety of problems. 

Depending on the mode of communication among the nodes, distributed algorithms can be broadly classified into three categories, namely, diffusion, consensus, and incremental <cit.>. Of these, the incremental update rule generally incurs the least amount of message passing overhead <cit.>, and is of interest in the present context. The incremental subgradient descent and its variants have been widely applied to large-scale problems, and generally exhibit faster convergence than the traditional steepest descent algorithm and its variants <cit.>.

The stochastic gradient and subgradient algorithms are well-known within the machine learning and signal processing communities <cit.>. The incremental stochastic subgradient method, with cyclic, random, and Markov incremental variants, was first proposed in <cit.>. The asymptotic analysis of dual problem in the present work follows the same general outline as that of the cyclic incremental algorithm in <cit.>, with additional modifications introduced to handle asynchrony. It is emphasized that these modifications are not straightforward, since the delayed stochastic subgradient is not generally a descent direction on an average. The present work also allows delays in both, primal and dual update steps, and establishes asymptotic near-optimality and feasibility of the primal allocation policies. Finally, saddle point algorithms have recently been applied to unconstrained <cit.> or proximity-constrained <cit.> network optimization problems, but do not readily generalize to the general form constrained optimization problem considered here.

Asynchronous algorithms have also been considered within the Markov decision process framework <cit.>, though the setup there is quite different and does not apply to the problem at hand. On the other hand, asynchronous first order methods have attracted a significant interest from the machine learning community <cit.>. For problems where the exact subgradient is available at each node, the asynchronous alternating directions method of multiplier (ADMM) has been well-studied <cit.>. The present work considers stochastic algorithms, and thus differs considerably in terms of both analysis and the final results. Even among algorithms utilizing stochastic subgradients, the definition of asynchrony varies across different works. One way to model asynchrony is to allow each node to carry out its update according to a local Poisson clock. This approach is followed in <cit.>, all of which consider various consensus-based distributed subgradient algorithms. The asynchronous adaptive algorithms in <cit.> also subscribe to the same philosophy, with decoupled node-updates due to communication errors, changing topology, and node failures. The incremental algorithm considered here is very different in terms of operation and analysis.

On the other hand, asynchronous operation can be modeled via delayed gradients or subgradients utilized for the updates. A  consensus-based stochastic algorithm proposed in <cit.>, and utilizes randomly delayed stochastic gradients. Along similar lines, asynchronous saddle point algorithms for network problems with edge-based constraints have recently been proposed <cit.>. Finally, for the unconstrained variants of the problem, a non-parametric approach has been proposed in <cit.>. Different from these works, the network resource allocation framework considered here allows generic convex constraints. Further, the incremental algorithm developed here handles stale subgradients while incurring significantly lower communication overheads. Asynchronous variants of the classical or averaged stochastic gradient methods have been proposed in <cit.>. The generic problem of interest here is that of the minimization of a sum of private functions at various nodes. Further, a network with star topology is considered, with updates being carried out using delayed  gradients collected at the fusion center. Different from these works, the proposed algorithm is incremental, does not require a fusion center, and is therefore more relevant to the network resource allocation problem at hand. Unlike these works, the present work also avoids making any assumptions on the compactness of the domain of the dual optimization problem. Before concluding, it is remarked that this work develops convergence results that hold on an average. Stronger results, where convergence is established in an almost sure sense, require a more involved analysis, and are not pursued here. 



The notation used in this paper is as follows. Scalars are represented by small letters, vectors by small boldface letters, and constants by capital letters. The index t is used for the time or iteration index. The inner product between vectors a and b is denoted by a,b. For a vector , projection onto the non-negative orthant is denoted by []^+. The expectation operation is denoted by . The notation ∇ is used for gradient and ∂ is used for the subgradient. The Euclidean norm is denoted by ·.


		


§ PROBLEM FORMULATION




 §.§ Problem statement

This section details the stochastic resource allocation problem at hand for a network with K nodes. The stochastic component of the problem is captured through the random network state, comprising of the random vectors ^i ∈^q for each node i ∈{1, …, K}, with unknown distributions. The overall problem is formulated as follows. 

	
		
    :=  max      f^i(^i) 
    s.t.      ^̆i(^i) + [^̌i(^i,^i_^i)]≽ 0
       ^i∈^i,  ^i∈^i

	
where the expectation in (<ref>) is with respect to the random vector ^i and 𝖯 is finite. The optimization variables in (<ref>) include the resource allocation variables {^i∈ℝ^n}_i=1^K and the policy functions {^i:^q→^p}_i=1^K, under the constraints (<ref>)-(<ref>).  Note that, the constraints in (<ref>) are required to be satisfied on an average, whereas those in (<ref>) are needed to be satisfied instantaneously. The functions f^i:^n → are assumed to be concave, and the sets ^i ⊆^n, convex and compact. The constraint function at node i is vector-valued, and is given by ^̆i(^i):=[u^i_1(^i) ⋯ u^i_d(^i)]^T, where {u_k^i(^i):^n →}_k=1^d are concave functions. On the other hand, no such restriction is imposed upon the vector-valued function ^̌i:^p×^q→^d and the compact set of functions {𝒫^i}_i=1^K. Of course, the overall problem still needs to adhere to certain regularity conditions (see Sec. <ref>), such as the Slater's constraint qualification and Lipschitz continuity of the  gradient function; see (A1)-(A7). 
	
Since the distribution of ^i is also not known in advance, it is generally not possible to solve for  in an offline manner. Therefore, an online algorithm is sought to solve problem `on the fly' as the independent identically distributed (i.i.d.) random variables {_t^i}_t∈ are realized and observed. For brevity, we denote _t^i:= _^i_t^i and _t^i(_t^i,^i):=^̆i(^i) + ^̌i(_t^i,^i_^i_t). Therefore, it is possible to write (<ref>) equivalently as [_t^i(_t^i,^i)]≽ 0. The algorithm outputs a sequence of vector pairs {_t^i, _t^i}_t, that are used for allocating resources in a timely manner. Towards this end, the stochastic dual descent algorithm has been proposed in <cit.>, which yields allocations that are almost surely near-optimal and provably convergent.


	
In the present paper, the focus is on networked systems where both, allocations (^i,^i) and the functions f^i and ^i_t are private to each node i. Likewise, the random variable _t^i is also observed and estimated locally at each node i. In other words, while the nodes can exchange dual variables and numerical values of the  gradients, they may not be willing to reveal the full functional form of the objective or constraint functions and other locally estimated quantities, owing to privacy and security concerns. Such privacy-preserving cooperation is common for many secure multi-agent systems <cit.>. To this end, the nodes may be arranged in a star topology, and utilize a centralized controller for collecting and distributing various algorithm iterates. Alternatively, ring topology may be used, allowing a fully distributed implementation, where the exchanges occur only between two immediate neighbors.
 In order to clarify the problem formulation considered in (<ref>), the following simple example is considered. 

 
 Example 1. Consider the problem of network utility maximization over a wireless network consisting of K nodes. The aim is to maximize the network-wide utility given by
 
    ∑_i=1^KU(r^i)
 
 where U(·) is a concave function that quantifies the utility obtained by the node i upon achieving a rate r^i ∈[r_min,r_max]. The channel is assumed to be time-varying, and for each channel realization h^i, node i allocates the power p^i_h^i, achieving the instantaneous rate of log(1+h^ip^i_h^i), where the noise power is assumed to be one. The goal is to maximize the utility in (<ref>) subject to constraints on the average rate and the average power consumption, and the full problem can be written as (cf. (<ref>)):
 

 
 	
    max_r^i,p^i   ∑_i=1^KU(r^i)   
    s.t. ∑_i=1^K(1/2log(1+h^ip^i_h^i))   ≥∑_i=1^Kr^i
    ∑_i=1^Kp^i_h^i   ≤ P_max
    r^i∈[r_min,r_max],  p^i∈𝒫^i

 
 It is remarked that 𝒫^i is a set of functions p^i:→, while p^i_h^i is a random variable that depends on h^i. That is, the optimization variables in (<ref>) include the rates r^i and the power allocation functions p^i.	

	


 §.§ Existing approaches and challenges

 We begin with explicating the desirable features of an algorithm that seeks to solve (<ref>). Specifically, it is required that any such algorithm meets the following requirements.
 
 	
  F1. The algorithm should allow nodes to “fall behind” temporarily, e.g., under poor channel conditions and intermittent transmission failures.
 	
  F2. The algorithm should allow a distributed implementation, that is, without requiring a star-topology or an FC.
 

 These features are particularly important for large and heterogeneous networks where delays may be unavoidable and designating an FC may be impractical. Put differently, (F1) requires the algorithm to handle the inevitable delays that may occur due to temporarily poor channel conditions or noise. Complementarily, (F2) is an  architectural requirement that must be kept in mind when choosing or designing the algorithm.

	Since the number of constraints in (<ref>) are finite, the problem is more tractable in the dual domain. To this end, introducing a dual variable ∈ℝ^d_+ corresponding to constraint in (<ref>), the stochastic (sub-)gradient descent method was proposed for solving such problems in <cit.>. The Lagrangian of (<ref>) is given by
	
    L(,,)    ={f^i(^i)+,[ _t^i(_t^i,^i)]}

	where  and  collect the primal optimization variables {^i}_i=1^K and {^i}_i=1^K respectively. Next, the dual function is obtained by maximizing L with respect to  and . Since the Lagrangian is expressed as a sum of K terms, each depending on a different set of variables, the maximization operation is separable and the dual function takes the following form:
	
    D()    = ∑_i=1^K max_^i∈^i, ^i∈^i[f^i(^i)+,([ _t^i(_t^i,^i)])] 
       =: ∑_i=1^K D^i().

	The dual problem is given by
		

	
    = min_∈ℝ^d_+ D^i().

While for general problems, it only holds that 𝖣≥𝖯, that the stochastic resource allocation problem considered here has a zero duality gap, i.e., 𝖯 = 𝖣 <cit.>. The result utilizes the Lyapunov's convexity theorem and holds under strict feasibility (Slater's condition), bounded subgradients, and continuous cumulative distribution function of ^i for each i. It is remarked that similar results are well-known in economics <cit.>, wireless communications <cit.>, and control theory <cit.>.

	The result on zero duality gap legitimizes the dual descent approach, since the dual problem is always convex, and the resultant dual solution can be used for primal recovery. To this end, similar problems in various contexts have been solved via the classical dual descent algorithm <cit.>, wherein the primal updates utilize various sampling techniques. It is remarked however that from a practical perspective, solving the dual problem alone is not sufficient, since online allocation of power or rate variables necessitates determining the primal optimum variables {^i^⋆, ^i^⋆}. In the present case, since ^i is infinite dimensional, primal recovery and consequently, online resource allocation, is not straightforward.   
	




Since the distribution of ^i is not known in advance, solving (<ref>) via classical first or second order descent methods requires a costly Monte Carlo sampling step <cit.>. Instead, the use of stochastic subgradient descent has been proposed in <cit.>, which takes the following form for t≥ 1,
	
		
  D1. Primal updates: At time t, node i observes or estimates _t^i, and allocates the resources in accordance with:
		
    {^i_t(_t), ^i_t(_t)}   = _∈^i, ∈Π^i_t  f^i()+_t, _t^i(,)

		
  D2. Dual update: The dual updates at time t take the form:
		
    _t+1=[_t-_t^i(_t^i(_t),^i_t (_t))]^+.

	
	Here, Π_t^i :={^i_^i_t∈^p | ^i ∈^i} is the set of all legitimate values of the vector ^i_^i_t. The term _t^i(_t^i(_t),^i_t (_t)) is a stochastic gradient of the dual function D^i() at =_t. Further for notational brevity, _t^i():=_t^i(_t^i(),^i_t ()) is used throughout the paper. Recall that for a given , _t^i() is stochastic and depends on the random variable ^i_t, as discussed in Sec. <ref>. The algorithm is initialized with an arbitrary _1 and the resulting allocations are asymptotically near optimal and feasible. A constant step-size stochastic gradient descent algorithm is utilized in the dual domain, which not only allows recovery of optimal primal variables via averaging, but also bestows it the ability to handle small changes in the network topology or other problem parameters. The algorithm can be implemented in a distributed fashion in a network with star-topology, with the help of a fusion center (FC). Within the FC-based implementation, the primal iterates are calculated and used locally at each node i. At the end of each time slot, the node i communicates the gradient component ^i_t(_t) to the FC, which carries out the dual update (<ref>) and broadcasts _t+1 to all the nodes in the network. Summarizing, the stochastic algorithm is preferred over its deterministic counterpart since it does not require Monte Carlo iterations, yields asymptotically near-optimal resource allocations, and is provably convergent if the stochastic process {^i_t} is stationary. 
	
It is remarked that since (<ref>) is infinite dimensional, full primal recovery is generally not possible using such dual methods. Existing algorithms only allow partial primal recovery, as will also be possible via Theorem 2. Specifically, it is well-known that while the running average 1/T∑_t = 1^T ^i_t can be viewed as the approximate version of the primal optimum ^i^⋆, no such interpretation exists for the infinite-dimensional variable ^i. For instance, the running average of ^i_t cannot be meaningfully related to the corresponding optimum ^i^⋆ <cit.>. Nevertheless, the resource allocation carried out using the primal iterates {^i_t(_t),^i_t(_t)} still ensures near-optimality and asymptotic feasibility (cf. Theorem 2). 	

In view of the desiderata (F1)-(F2), observe that a network implementation of (<ref>)-(<ref>) is still impractical since it is synchronous and FC-based, and thus has relatively stringent communication requirements. In particular, the algorithm necessitates that each node exchanges messages (i.e. ^i_t(_t) & _t) with the FC at every time-slot, thereby incurring a large communications cost. Since the updates (<ref>)-(<ref>) must occur before the network state changes, the nodes must synchronize and cooperate in order to meet these deadline constraints, ultimately increasing message passing overhead and consuming more energy. Further, nodes in large networks are often heterogeneous, and may not always be able to transmit the gradients within the stipulated time. Finally, if the nodes are not deployed in a star-topology around the FC, the need for multi-hop communications further increases the delays, results in heterogeneous energy consumption, and increases protocol overhead. In all such cases, the FC must wait for the updates to arrive from all the nodes, possibly requiring all the nodes to skip resource allocation for one or more time slots, and resulting in a suboptimal asymptotic objective value.
	





	

§ PROPOSED ALGORITHM

	This section details the proposed stochastic dual descent algorithm that incorporates the features (F1)-(F2) in its design. To begin with, Sec. <ref> describes the asynchronous variant that tolerates delayed gradients still resulting in near-optimal resource allocation. Next, Sec. <ref> details the more general AIS-DD algorithm that is amenable to a distributed implementation. 

	
	

 §.§ Asynchronous stochastic dual descent

	The asynchronous stochastic dual descent algorithm addresses (F1), and proceeds as follows for all t ≥ 1:
	
		
  * Primal update: At each time t, node i solves
		
    {_t^i(_t-π_i(t))   , _t^i(_t-π_i(t))}
       :=_∈𝒳^i,∈Π_t^i  f^i()+_t-π_i(t),_t^i(,)

		for all 1 ≤ i ≤ K, and some finite delay π_i(t) ≥ 0. 
		
  * Dual update: The dual update at time t is given by
		
    _t+1=[_t-(_t-δ_i(t)^i(_t-τ_i(t)))]^+

		
the stale  gradient, evaluated at time t-δ_i(t), is given by  
    _t-δ_i(t)^i   (_t-τ_i(t))
       :=_t-δ_i(t)^i(_t-δ_i(t)^i(_t-τ_i(t)),_t-δ_i(t)^i(_t-τ_i(t)))
 where the total delay is denoted by  τ_i(t):=π_i(t)+δ_i(t) and π_i(t), δ_i(t) ≥ 0. 







			
Different from (<ref>), the resource allocation in (<ref>) utilizes an old dual variable, _t-π_i(t).	Further, the dual update is also carried out using an old gradient _t-δ_i(t)^i(_t-τ_i(t)). The two modes of asynchrony introduced in (<ref>)-(<ref>) allow the primal and dual updates to be carried out at different time scales. In other words, while the resource allocation at each node still occurs at every time slot, the rate at which the dual variables and the  gradients are exchanged may be different. In order to highlight the asynchronous nature of the algorithm, the implementation of (<ref>)-(<ref>) is now described from the perspective of the FC and that of node i, in Algorithms <ref> and <ref>, respectively.
		

		

				

						
Observe that in Algorithms 1 and 2, some steps are `optional,' which in the present case, means that they can, at times, be skipped. These steps are however still required to be carried out `often enough', so that the total delay τ_i(t) is bounded for each node i; cf. (A4) in Sec. <ref>. Nevertheless, the optional steps in these algorithms allow the dual updates to occur at a different rate. For instance, as long as each packet is correctly time-stamped, the dual updates at the FC may occur as and when the  gradients become available, instead of following a fixed schedule. 
							
						
						The ability to postpone or skip transmissions is important in the context of large heterogeneous networks. For instance, transmissions from the nodes to the FC often requires a multiple access protocol, inter-node coordination, and energy budgeting at each node. Consequently, energy-constrained nodes may extend their lifetime simply by scheduling their transmissions once every few time slots. Similarly, energy harvesting nodes may only transmit when sufficient energy is available, choosing to stay silent in times of energy paucity. The slower nodes may even skip the  gradient calculation, as long as the resources are allocated in time. 						
						Finally, the communication between the nodes and the FC may also incur delays, arising from queueing, processing, or retransmission at various layers in the protocol stack. The flexibility of carrying out updates with stale information makes the network tolerant to such delays. 
						
						
						

						

 §.§ Asynchronous Incremental Stochastic Dual Descent
 
						This subsection details an incremental version of the asynchronous algorithm introduced in Sec. <ref>, that obviates the need for an FC and is thus endowed with both (F1) and (F2). The AIS-DD algorithm allows each node to perform the partial dual update itself, while passing messages to nodes along a cycle. Specifically, for a network with a ring topology, such that node i passes dual variable _t^i to node i+1 and so on, the primal and dual updates take the following form. 
						
						
							
  * Primal update:   At time t, node i solves
							
    (^i_t   (^i-1_t-π_i(t)), _t^i(^i-1_t-π_i(t)))
       :=_∈𝒳^i, ∈Π_t^i f^i()+^i-1_t-π_i(t),_t^i(,).

							
  * Dual update: 	At time t, the dual update at node i takes the form
							
    ^i_t=   [^i-1_t-  [_t-δ_i(t)^i(^i-1_t-τ_i(t))]]^+
 
						
						 where,
 
    _t-δ_i(t)^i   (^i-1_t-τ_i(t))
    
    						 :=   _t-δ_i(t)^i(_t-δ_i(t)^i(^i-1_t-τ_i(t)),_t-δ_i(t)^i(^i-1_t-τ_i(t)))
 and _t^0 is read as _t-1^K and _t=_t^0 will be used to evaluate the performance of the asynchronous incremental  algorithms.
						A key feature of the AIS-DD algorithm is that the message passing and the dual updates occur in parallel with the resource allocation, as shown in Fig. <ref>. The full implementation details are provided in Algorithm <ref>. 
						
								
								Here, the two optional steps may be repeated as long as the received _t'^i-1 is still old, that is, t' ≤ t. As in Sec. <ref>, the nodes are allowed to halt the updates temporarily, as long as they “catch up,” eventually. In other words, the updates for time t' must be carried out before time t'+τ so as to ensure that τ_i(t) ≤τ for all t. Interestingly, although resources are allocated at every time slot, the network may or may not carry out one or more message passing rounds per time-slot. It is remarked that the update in (<ref>) must still be performed once at every node for each time index t'. Equivalently, the algorithm runs on two `clocks,' one dictating the resource allocation and synchronous with the changes in the network state, and the other governed by the rate at which messages get passed around the network. In the next section, we will establish that the such an algorithm still converges, as long as the difference between the two clocks is bounded. In summary, the AIS-DD algorithm has all the benefits of the asynchronous dual descent algorithm of (<ref>)-(<ref>), while allowing a distributed implementation. 
								
As with classical incremental algorithms, the nodes must communicate along a ring topology. Strictly speaking, the message passing overhead is minimized if the updates occur along a Hamilton cycle <cit.>. Even when the network does not admit a Hamilton cycle, an approximate cycle can be found using a random walk protocol <cit.> or the protocol described in <cit.>. It is remarked that such a route need only be found once, at the start of the algorithm.
							

								

§ CONVERGENCE RESULTS

								This section provides the convergence results for the AIS-SD and AIS-DD algorithm. We begin with developing and analyzing the convergence of AIS-SD algorithm  (cf. Theorem <ref>). It is emphasized that the AIS-SD algorithm is general-purpose, and can be used to minimize any sum of functions in an incremental and asynchronous manner.  Subsequently, the asynchronous incremental stochastic gradient descent algorithm is applied to (<ref>) in the dual domain, and a primal-averaging method is proposed that yields asymptotically near-optimal allocations (cf. Theorem <ref>). We begin with stating the assumptions and briefly reviewing some of the known results (Sec. <ref>). The results for the dual case are outlined in Sec. <ref>, while the near-optimality of the resource allocation is established in Sec. <ref>. 
								
								

								

 §.§ Assumptions and known results

								This subsection begins with the discussion of the following general optimization problem:

								
    = min_∈Ł D^i()

								where,  is the optimization variable, Ł⊆ℝ^d is a non-empty, closed, and convex set,  is finite, and the objective function separates into node-specific cost functions D^i. The goal is to solve  (<ref>) using only the stochastic  subgradients _t^i() of D^i(). It is emphasized that the general results presented in this subsection do not required D^i to be differentiable. As in (<ref>), _t^i() is stochastic due to its dependence on the random variable _t^i that is first observed at node i at time t. Besides the network resource allocation problem considered here, (<ref>) also arises in the context of machine learning <cit.> and distributed parameter estimation <cit.>. Before describing the known results related to (<ref>), the necessary assumptions are first stated. 
								
								
									
  A1.  Non-expansive projection mapping. The projection mapping  satisfies -≤- for all  and . 
									
									
  A2.  Zero-mean time-invariant error. Given , the averaged subgradient function satisfies ∂ D^i() = [_t^i()]. 
									
									
  A3.  Bounded moments.   Given ∈Ł, the second moment of _t^i() is bounded as follows:
									
    [_t^i()^2] ≤V_i^2.

								 
								
								These assumptions are not very restrictive, and hold for most real-world resource allocation problems. A stochastic incremental algorithm for solving (<ref>) was first proposed in <cit.>. Given a network with ring topology, the updates in <cit.> take the form
								
    _t^i    = _t^i-1 - ϵ_t^i(_t^i-1)

								where _t^0 is read as _t-1^K. It was shown in <cit.>, that under (A1)-(A3), the iterates ^i_t are asymptotically near optimal in the following sense

								
    t→∞liminf [D(_t)]≤𝖣+Ø().


								
								where _t= _t^0=_t-1^K. Further, for the case when the step size is diminishing, i.e. _t satisfies lim_T→∞∑_t=1^T_t = ∞ and lim_T→∞∑_t=1^T_t^2 < ∞, it holds that
								
    lim inf_t→∞D(_t)= 𝖣.

								
								This paper provides the corresponding results for the asynchronous case, where the subgradient in (<ref>) is replaced by an older copy _t-δ_i(t)^i(_t-τ_i(t)^i-1), that is, the stochastic subgradient of D^i()  that depends on the random variable _t-δ_i(t)^i and is evaluated at = _t-τ_i(t)^i-1. The delays satisfy τ_i(t) ≥δ_i(t) ≥ 0 and for the special case of no delay, the stochastic subgradient simplifies to _t^i(_t^i-1) as in (<ref>). The following additional assumption regarding the delays δ_i(t) and τ_i(t) is stated. 
								
									
  A4.  Bounded delay.  For each 1≤ i ≤ K and t≥ 1, it holds that 0≤δ_i(t) ≤τ_i(t) ≤τ < ∞.
								
								The boundedness assumption on the delay in (A4) allows us to develop convergence results that hold in the worst case, and has been widely used in the context of asynchronous algorithms <cit.>. It is remarked that an alternative assumption, made in <cit.>, allows the delays δ_i(t) and τ_i(t) to be random variables with unbounded supports but finite means, but is not pursued here. Even with bounded delays, the extension to the asynchronous case is not straightforward, since the the old stochastic subgradients are not necessarily descent directions on an average. Indeed, the resulting subgradient error at time t, defined as
								
    _t,δ_i(t)^i:=   [∂ D^i(_t^i-1)-_t-δ_i(t)^i(_t-τ_i(t)^i-1)]

								is neither zero-mean nor i.i.d. In other words, the asynchronous algorithm cannot simply be considered as a special case of the inexact subgradient method.  
								
It is worth pointing out that there is a subtle difference between the definition of the delayed stochastic gradient considered here, and those considered in <cit.>. Specifically, the delayed gradient in these works takes the form _t^i(_t-τ_i(t)^i-1) instead of the one in (<ref>). As a result, given _t-τ_i(t)^i-1, the gradient error at time t in these papers is indeed zero mean and i.i.d., an assumption that simplifies the analysis to a certain extent. It is also remarked that the definition of the delayed stochastic gradients in <cit.> is however similar to that considered here. Different from these works, the dual convergence results developed here consider subgradients instead of gradients, and are therefore applicable to a wider range of problems. 
								
								
								Within the context of network resource allocation, it is also important to study the (near-)optimality of the allocations {^i_t, ^i_t}. Towards this end, some additional assumptions are first stated.
								
									
  A5.  Non-atomic probability density function:  The random variables {_t^i}_i=1^K have non-atomic probability density functions (pdf).
									
									
  A6.  Slater's condition:  There exists strictly feasible (^i,^i), i.e., [_t^i(_t^i,^i)]> 0.
									
  A7.  Lipschitz continuous  gradients.  Given , ' ∈Ł, there exists L_i < ∞ such that 
									
    ∇ D^i()-∇ D^i(')≤ L_i-'.
 
								
								In (A5), for {_t^i}_i=1^K to have a non-atomic pdf, it should not have any point masses or delta functions. Note that this requirement is not restrictive for most applications arising in wireless communications; see e.g. <cit.>. The Slater's condition is a standard assumption that ensures that 𝖯 = 𝖣 and consequently, since 𝖯 is finite, so is 𝖣. The Lipschitz condition in (A7) is however restrictive, since it requires the dual functions D^i() to be differentiable with respect to . In other words, with (A7), _t^i() is a stochastic  gradient, not a subgradient. It is remarked however that (A7) always holds if f^i(^i) is strongly convex. Moreover, it is generally possible to enforce (A7) artificially by adding a strongly convex regularizer (such as θ^i_2^2) to the primal objective <cit.>. Note however that (A5)-(A7) will not be utilized while establishing the dual convergence results. 
								
								The incremental or asynchronous  gradient methods have thus far never been applied to the problem of network resource allocation. For the classical stochastic dual descent method (cf. (<ref>)-(<ref>)], it is known that under (A1)-(A3) and (A6), the average resource allocations ^i:=^i_t are asymptotically feasible and near-optimal <cit.>.
				

								

 §.§ Convergence of the AIS-SD algorithm
 
								This subsection provides the convergence results for the AIS-SD algorithm, applied to (<ref>). For the general case, the updates take the following form:
								
    _t^i    = _t^i-1 - _t^i_t-δ_i(t)(^i-1_t-τ_i(t))    1≤ i ≤ K

								where ϵ_t is the step-size, _t^i() is a stochastic subgradient of D^i() and ^0_t is read as ^K_t-1. Since the dual problem (<ref>) is simply a special case of (<ref>), the results developed here also apply to the iterates {_t^i} generated by Algorithm 1. In order to keep the discussion generic, the results are presented for both, diminishing and constant step sizes. 
								
								
									The following results apply to the iterates generated by (<ref>) with _t=_t^0 under (A1)-(A4).
									
										
  * Diminishing  step-size: If the positive sequence {_t} satisfies lim_T→∞∑_t=1^T_t = ∞ and lim_T→∞∑_t=1^T_t^2 < ∞, then it holds that
										
    lim inf_t→∞[∑_i=1^KD^i(_t)]=𝖣.

										
  * Error bound for constant step size: For _t =  > 0, and any arbitrary scalar η > 0, it holds that

    min_1≤ t≤ T∑_i=1^K   D^i(_t)≤𝖣 +  C(τ)+η/2

where T ≤ B^2_0/η. Here, τ is the maximum delay as defined in (A4), C(τ):=C_1+(C_2+τ C_2'), C_1=KV^2, C_2:=2KV^2K-1/2, C_2':=4K^2V^2, and B_0 is such that _1-≤ B_0.
										









									 
								
								A popular choice for the diminishing step-size parameter _t required in Theorem <ref>(a) is _t=t^-α for α∈ (1/2, 1). For this case, the objective function in (<ref>) converges exactly to the dual optimum. On the other hand, with a constant step size , the minimum objective value comes to within an O()-sized ball around the optimum as T →∞. More precisely, the result in (<ref>) provides an upper bound on the number of iterations required to come η-close to this ball. Different from the results in <cit.>, the size of the ball now depends on the maximum delay τ, quantifying the worst-case impact of using delayed subgradients.
								
								The proof of Theorem <ref> follows the same overall structure as in <cit.>, with appropriate modifications introduced to handle the asynchrony. To begin with, the following intermediate lemma splits a function related to the optimality gap in Theorem <ref> into three different terms, and develops bounds on each. The proof of the following lemma is provided in Appendix <ref>.
								
								 
									Under (A1)-(A4), the iterates generated by (<ref>) with _t=_t^0 satisfy the following bounds:
									
    ∑_t=1^T∑_i=1^K2_tD^i(_t)-𝖣   ≤B_0^2+ I_0+I_1

									where,

									
    I_0    := ∑_t(_t^2 KV^2+2τ KV^2∑_i_t_t-τ_i(t)
        + 2V^2∑_i(i-1)_t_t-τ_i(t))
    
    									I_1   :=2∑_t=1^T∑_i=1^K_t _t-δ_i(t)^i(_t-τ_i(t)^i-1), _t-τ_i(t)^i-1-_t^i-1
       ≤ 2τ KV^2∑_t=1^T∑_i=1^K_t _t-τ_i(t).

									where _1^0 - ≤ B_0. Note that _1^0=_1.
								
								Having developed the necessary bounds, the proof of Theorem <ref> is presented next. 
								[Proof of Theorem 1] 
									For the positive sequence {_t}, it holds that
									
    ∑_t=1^T∑_i=1^K2_t D^i(_t)≥(min_1≤ t≤ TD^i(_t)) ∑_t=1^T2_t .

									Substituting the bounds obtained in Lemma <ref>, and noting that it always holds that _t ≤_t-τ for all τ≥ 0, we obtain
									
									
									
									
									
    min_1≤ t≤ T   D^i(_t) -𝖣
       ≤B_0^2+C_1∑_t=1^T_t^2  +(C_2+τ C_2')∑_t=1^T^2_[t-τ]_+/2∑_t=1^T_t

									where, C_1:= KV^2, C_2':=4K^2V^2 and C_2:=2KV^2K-1/2. Note that in (<ref>), we have used the notation _[t-τ]_+:=_1 for all t ≤τ.  Next, for the case when _t is diminishing, and satisfies  lim_T→∞∑_t=1^T_t = ∞ and lim_T→∞∑_t=1^T_t^2 < ∞, the numerator of the bound on the right stays bounded, while the denominator grows to infinity. Consequently, taking the limit of T →∞ on both sides of (<ref>), the required result in (<ref>) follows.
									Observe that when the step size is constant, the bound in (<ref>) can be written as
									
    min_1≤ t≤ TD^i(_t)-𝖣   ≤B_0^2/2 T+/2 C_1+/2 (C_2+τ C'_2)
       ≤B_0^2/2 T+/2 C(τ)

									where C(τ) is as defined in Theorem <ref>. In the limit as T→∞, the bound becomes
										
    inf_t≥ 1D(_t)≤+ C(τ)/2.

									which is the asymptotic version of the result in (<ref>). 
									
									The rate result in (<ref>) builds upon a similar result from <cit.>. Intuitively, D(_t) continues to decrease as long as it is significantly larger than 𝖣. The rest of the proof characterizes the resulting decrement rigorously and subsequently invokes the monotone convergence theorem in order to establish that D(_t) must eventually come close to 𝖣. Given arbitrary η >0 and recalling that _t := _t^0, define the sequence 
	

	
    _t+1:=_t+1 ;  if  D(_t)≥+ C(τ)+η/2
    ^⋆ ;  otherwise.
	
	Alternatively, _t is same as _t until _t enters level set defined as
	

	
    L={∈Λ | D()< + C(τ)+η/2}

	and _t terminates at ^⋆. From (<ref>) and Lemma <ref>, we have for constant step size _t= that
	
    _t+1-^⋆^2≤   _t-^⋆^2-2[D(_t)-]
       +^2C(τ)

	where _t+1:=_t^K and _t=_t^0. Next define
	
    z_t:= 2[D(_t)-]-^2C(τ)  if  _t∉ L
    
    	0  if  otherwise,
 
	so that (<ref>) can be written as
	
    _t+1-^⋆^2≤_t-^⋆^2-z_t.
From the Monotone convergence theorem, we have that ∑_t=1^∞z_t<∞, implying that there exists T < ∞, such that z_t=0 for all t≥ T. Observe that for the case when _t∉ L, it holds that 
	
    z_t=   2[D(_t)-]-^2C(τ)
    ≥   2[+ C(τ)+η/2-]-^2C(τ)
    	≥η.
 
	Consequently, it follows from (<ref>) that

	
    𝔼[||_T+1-^⋆||^2]   ≤ ||_1-^⋆||^2-∑_t=1^Tz_t 
       ≤ B_0 - ∑_t=1^Tz_t.
Since the term on the left is non-negative, we have that
	B_0 ≥∑_t=1^Tz_t≥ T η, yielding the required bound on T.





 §.§  Primal near optimality and feasibility

The AIS-SD algorithm of Sec. <ref>, when applied to solve the dual problem in (<ref>), is referred to as the AIS-DD algorithm. In order to ensure that the results developed thus far continue to apply to the dual problem, assumptions (A5)-(A7) are also required. As mentioned earlier, for the primal problem, Ł is simply the non-negative orthant implying that 𝖣 is finite. This subsection establishes the average near-optimality of the AIS-DD algorithm in (<ref>)-(<ref>). Note that Theorem <ref> does not imply that the allocations {_t^i,_t^i} converge. Instead, the results will make use of the ergodic limit variable 
								
    ^i_T:=1/T∑_t=1^T^i_t
 
								for each 1 ≤ i ≤ K. The main theorem for this subsection is presented next.
								 
									Under (A1)-(A7) and for constant step size >0, the iterates generated by (<ref>)-(<ref>) follow:
									
										
  * 
    lim inf_T→∞f^i(^i_T)≥   𝖯-(C_3+τ C_4)

										where,

										
    C_3 =   (V^2K^2)/2
    
    										C_4=   K^2BLV+2K^2V^2.

										
  *  
										
    lim inf_T→∞1/T∑_t=1^T∑_i=1^K[_t-τ_i(t)^i(^i-1_t-τ_i(t))]≽ 0.

									
								
								
								Intuitively, the resource allocations in (<ref>) are near-optimal, with optimality gap depending on the step size  and the delay bound τ. Further, the allocations are almost surely asymptotically feasible, regardless of the the delay bound or the step size. As in Sec. <ref>, the proof of Theorem <ref> proceeds by first splitting the optimality gap into three terms and developing bounds on each. The required results are summarized into the following intermediate Lemmas, whose proofs are deferred to Appendices <ref> and <ref> respectively.  
								
									Under (A1)-(A6), the iterates _t^i obtained from (<ref>) are bounded on an average, i.e., there exists B < ∞ such that _t≤ B for all t ≥ 1. 
								
							 
									Under (A1)-(A7), the iterates generated by (<ref>)-(<ref>) satisfy the following bounds:
										
    f^i(_T^i)≥𝖣-	I_2 - I_3

									where,

									
    I_2 :=   1/T∑_t=1^T∑_i=1^K[D^i(_t)-D^i()]
       ≤ V^2K(K-1)/2 
    
    									I_3 :=   1/T∑_t=1^T∑_i=1^K,∇ D^i()
       ≤_1^2/2 T+ KV^2 /2 + I_4 
    
    									I_4 :=   1/T∑_t=1^T∑_i=1^K,(∇ D^i()-_t-δ_i(t)^i(_t-τ_i(t)^i-1))
       ≤τ K^2V(BL+2V).

								
								Having established the intermediate results, the proof of Theorem <ref> is now presented. 
								
								
									The primal near-optimality can be established directly from Lemma <ref>. Specifically, summing the bounds for I_2, I_3, and I_4, and taking the limit as T →∞, the bound in (<ref>) follows.
									
									In order to establish (<ref>), observe that for any t≥ 1 and 1 ≤ i ≤ K, it holds that
									
    ^i_t   = [^i-1_t-_t-δ_i(t)^i(^i-1_t-τ_i(t))]^+
       ≽^i-1_t-_t-δ_i(t)^i(^i-1_t-τ_i(t))

									where the inequality holds element-wise. Summing both sides over all 1≤ t ≤ T and 1 ≤ i ≤ K, and rearranging, it follows that
									
    1/T∑_t=1^T∑_i=1^K_t-δ_i(t)^i(^i-1_t-τ_i(t))    ≽1/ T∑_t=1^T∑_i=1^K (_t^i-1 - _t^i)
       ≽^K_1-^K_t+1/ T.

									Finally, since _1^K ≽ 0, taking expectations on both sides, it follows that
									
    1/T∑_t=1^T∑_i=1^K[_t-δ_i(t)^i(^i-1_t-τ_i(t))] ≽ -B/ T

									where (<ref>) holds due to Lemma <ref>. In other words, given any α >0, there exists t_0 ∈ℕ such that for all T ≥ t_0, 
									
    1/T∑_t=1^T∑_i=1^K[_t-δ_i(t)^i(^i-1_t-τ_i(t))] ≽ -α.

									Taking the limit as T→∞, the result in (<ref>) follows.
								
								

§ APPLICATION TO CO-ORDINATED BEAMFORMING

								This section considers the co-ordinated downlink beamforming problem in wireless communication networks. The usefulness of the proposed stochastic incremental algorithm is demonstrated by applying it to the beamforming problem and solving it in a distributed and online fashion. Simulations are carried out to confirm that the performance of the proposed algorithm is close to that of the centralized algorithm. 
												

								

 §.§ Problem formulation
  
								Consider a multi-cell multi-user wireless network with B base stations and U users. Each user j ∈{1, …, U} is associated with a single base station b(j) ∈{1, …, B}, and the set of users associated with a base station i is denoted by _i:={j | b(j) = i}. For the sake of consistency, this section will utilize indices i and m for base stations, and indices j, k, and n for users, with the additional restriction that b(j) = b(k) = i and b(n) = m. Within the downlink scenario considered here, user j can only receive data symbols s_j∈ℂ from its associated base station b(j). The signals transmitted by the base station i intended for other users k ∈_i∖{j}, as well as the signals transmitted by other base stations m ≠ i constitute, respectively, the intra-cell and inter-cell interference at user j. The base station i, equipped with N_i transmit antennas, utilizes the transmit beamforming vector _j∈ℂ^N_i × 1 for each of its associated user j ∈_i. Consequently, the received signal at user j is given by
								
    y_j = _ij^H(_js_j + ∑_k ∈_i∖{j}_ks_k) + ∑_m ≠ i 
     n ∈_m_mj^H_ns_n + e_j

								where _ij denotes the complex channel gain vector between base station i and user j, and e_j is the zero mean, complex Gaussian random variable with variance σ^2 that models the noise at user j. Assuming s_k to be independent, zero-mean, and with unit variance, the expression for the signal-to-interference-plus-noise ratio (SINR) at user j is given by
								
    SINR_j := _ij^H_j^2/∑_k ∈_i∖{j}_ij^H_k^2 + ∑_m ≠ i∑_n ∈_m_mj^H_n^2 + σ^2

								where i = b(j) is the associated base station. 
								
								Within the classical co-ordinated beamforming framework, the goal is to design the beamformers {_j}_j=1^U so as to minimize the transmit power, while meeting the SINR constraints at each user. The required optimization problem becomes <cit.>
								

								
    min_{_j}_j=1^U   ∑_j=1^U _j^2    subject to SINR_j ≥γ_j  ∀ j

								where γ_j is a pre-specified quality-of-service (QoS) threshold for user j. While the beamforming vectors resulting from (<ref>) are optimal, the centralized nature of the optimization problem renders it impractical for application to real networks. For instance, the solution proposed in <cit.> requires the estimated channel gains {_ij} to be collected at a centralized location, where (<ref>) is solved via an iterative algorithm. In practice however, the entire parameter exchange and the algorithm must complete within a fraction of the coherence time of the channel, lest the designed beamformer becomes obsolete. Such a solution is therefore difficult to implement, not robust to node or link failures, and not scalable to large networks. 
								
								Observe that the modified version of (<ref>) can be written as
								

								
									
    min_{_j, I_nj}∑_j=1^U _j^2    subject to    
       _ij^H_j^2/∑_k ∈_i∖{j}_ij^H_k^2 + I_j^2 + σ^2≥γ_j   ∀ j
       ∑_m ≠ i∑_n ∈_m_mj^H_n≤I_j     ∀ j

								
								where i = b(j). Note that constraints in (<ref>) and (<ref>) will ensure the SINR is still greater than the required threshold of γ_j. It is due to the fact that the feasible set is restricted and feasible set of (<ref>) will be subset of that of (<ref>) and solution found for (<ref>) can be used for (<ref>). Next, the use of primal or dual decomposition techniques can yield a distributed algorithm for (<ref>). Nevertheless, such distributed algorithms also suffer from the limitations mentioned earlier, since the optimum beamforming vectors are required at every time slot. 
								
								On the other hand, within the uncoordinated beamforming framework, the optimization variable I_j in (<ref>) is replaced with a pre-specified threshold ρ. This renders (<ref>) separable at each base station, allowing beamforming vectors to be designed in parallel. However, the resulting beamformers are suboptimal, and may even render the problem infeasible if ρ is too small or too large. 
								

								
									
    min_{_j, I_nj}∑_j=1^U _j^2    subject to    
       _ij^H_j^2/∑_k ∈_i∖{j}_ij^H_k^2 + ρ^2∑_m ≠ icard(_m) + σ^2≥γ_j  ;   ∀ j
       _mj^H_n(t)≤ρm≠ i, n ∈ U_m, ∀ j.

								
									
									A compromise is possible within the stochastic optimization framework by requiring the bound in (<ref>) to only be satisfied on an average. Note that this amounts to relaxing the optimization problem (<ref>) since the SINR constraint is no longer binding at every time slot. The overall stochastic optimization problem can be expressed as 

									
										
    min_{_j(t), I_nj(t)}∑_j=1^U _j(t)^2    subject to 
       _ij^H_j(t)^2/∑_k ∈_i∖{j}_ij^H_k(t)^2 + I_j^2(t) + σ^2≥γ_j    ∀ j 
       ∑_m ≠ i∑_n ∈_m_mj^H_n(t)≤I_j(t)∀ j 
       _mj^H_n(t)≤ρm≠ i, n ∈ U_m, ∀ j

									
									
									where i = b(j). Different from (<ref>) or (<ref>), the stochastic optimization problem (<ref>) involves finding policies _j(t) and I_j(t), which are not necessarily optimal for every time slot t, but only on an average. Specifically, the intercell interference is bounded on an average [cf. (<ref>)], but also instantaneously [cf. (<ref>)], so as to limit the worst case SINR. The problem in (<ref>) can be readily implemented using the proposed  distributed and asynchronous stochastic dual descent algorithm. In contrast to (<ref>), the stochastic algorithm is not required to converge at every time slot, and allows cooperation over heterogeneous nodes. 
									

													

									

 §.§ Solution to optimization problem

									The AIS-DD algorithm proposed in Sec. <ref> can now be applied to solve (<ref>). To this end, associate dual variables λ_j for all users {j ∈[1,U]}, and observe that the primal variables at node i include {_j}_j ∈_i and {I_j}_j ∈_i. Departing from the notational convention used thus far, the subscript in λ_j is used for indexing the users, while time dependence is indicated by λ_j(t). Proceeding as in Sec. <ref>, and recalling that the indices j and n are such that i = b(j) ≠ b(n) = m, the operation at node i is summarized in Algo. <ref>. Observe that such an implementation entails allocating resources prior to the dual updates, and thus results in the delay of at least one, i.e., π_i(t) ≥ 1, compared to the synchronous version. Conversely, the dual updates occur as and when they are passed around, without creating a bottleneck for the resource allocation. For the sake of simplicity, it is assumed that the dual updates occur along the route 1, 2, …, K. 
							
									
									Next, simulations are carried out demonstrate the applicability of the stochastic algorithm to the beamforming problem at hand. For the simulations, we consider a system with B = 10 and U = 10, with one users per cell. Each of the  base stations  have ten antennas (N_i = 10), while  the other algorithm parameters are =0.5, σ^2=1, ρ=1.65 γ_j = 10 dB for all  j. In order to keep the simulations realistic, we assume that the delays in the dual updates arise from random events such as node and link failures. For the centralized algorithm, a random subset of four out of ten nodes are selected to transmit their current  gradients to the FC at every time slot. Since the FC utilizes old  gradients for the other nodes, it results in an average delay of 5.4 time slots. Similarly, for the incremental algorithm, it is assumed that at every time slot, five to fifteen dual update steps (cf. (<ref>)) occur, resulting in an average delay of 5.4 time slots. For instance, if at any time slot, only 8 nodes update, it will result in a delay π_i(t) = 2 at the remaining two nodes, where the dual update will occur at the start of the next time slot. The delay may increase further if fewer than 10 nodes update for consecutive time slots and conversely, may decrease if more than 10 updates occur per time slot. Fig. <ref> shows the running average of the primal objective function as a function of time using Monte Carlo simulations. For comparison, the performance of the classical centralized stochastic gradient method [cf. (<ref>)-(<ref>)], assuming  perfect message passing, is also shown. As evident, the performance loss due to the delays in the availability of the dual variables in minimal. 
										
										
									
									In order to motivate the stochastic formulation over the deterministic one, Fig. <ref> also compares the average transmit power and SINR achieved for the various cases and for different values of the parameter ρ. As expected, the distributed deterministic algorithm performs poorly since it forces the SINR bound to be a constant that does not depend on the channel. By design, the worst case SINR is bounded below by one at every time slot in both the deterministic formulations. Interestingly, the worst case SINR achieved for the relaxed stochastic formulation is also close to one on an average. In return, the stochastic algorithm yields an average transmit power that is equal to or below that obtained by the centralized deterministic formulation. In other words, it is always possible to artificially raise γ to a value that is slightly higher than one, so as to obtain an average SINR above 10 dB, while still getting near-optimal average transmit power. 
								
					
							
Next, we study the effect of delay on the rate of convergence of the AIS-DD algorithm. For this case, a simple system with B = 10 and U=10, and constant delays at all the nodes is considered. The base stations have ten antennas each (N_i=10) and  the other algorithm parameters are =0.2, σ^2=1, and ρ=1.65. Fig. <ref> shows the evolution of the primal objective function for various delay values. As expected, the convergence is slower if both π_i(t) and δ_i(t) are consistently larger. Interestingly however, a small increase in the delays amounts to only a marginal loss in performance.

Finally, in order to demonstrate the scalability of the proposed algorithm, Fig. <ref> shows an example run for a system with B=50 nodes and U =50. The  base stations have ten antennas each (N_i= 10), while  the other algorithm parameters are =0.5, σ^2=1, and ρ=5. The delay is generated in the similar manner as for the earlier simulations. It can be observed that even when the number of nodes is large, the difference between the performance of the synchronous and asynchronous algorithms remains relatively small.





								
								

									

§ CONCLUSION

									 This paper considers a constrained stochastic resource allocation problem over a heterogeneous network. An asynchronous incremental stochastic dual descent method is proposed for solving the same. The proposed algorithm utilizes delayed gradients for carrying out the updates, resulting in an attractive feature that allows nodes to skip or postpone some updates. The convergence of the proposed algorithm is established for both constant and diminishing step sizes. Further, it is shown that the resource allocations arising from the proposed algorithm are also asymptotically near-optimal. A novel multi-cell coordinated beamforming problem is formulated within the stochastic framework considered here, and solved via the proposed algorithm. Simulation results reveal that the impact of using stale stochastic gradients is minimal. 


§ PROOF OF LEMMA <REF>

									


  §.§.§ Preliminaries
 Before deriving the required bounds, some preliminary results are first obtained. Recall that the quantity _υ^i() denotes the stochastic (sub-)gradient of D^i() at time t = υ and evaluated at . Within the context of the dual descent algorithm, we also have that _υ^i():=_υ^i(^i_υ(),_υ^i()) for υ≥ 1 and all 1≤ i≤ K. In particular the updates in (<ref>) use υ = t-δ_i(t) and = _t-τ_i(t)^i-1, where τ_i(t) = δ_i(t)+π_i(t). It can be seen that for the special case of the synchronous algorithm, we have that τ_i(t) = δ_i(t) = 0 and the stochastic (sub-)gradient is written as _t^i(_t^i-1).





For the sake of convenience, let us denote _t-δ_i(t)^i := _t-δ_i(t)^i(_t-τ_i(t)^i-1). First, we establish that the distance between the iterates _t^i and _ℓ^i is bounded by a term that is proportional to the step size. From the updates in (<ref>), it holds for all t and 1≤ i ≤ K that
									
    _t^i-_t^i-1   = _t^i-1 -_t_t-δ_i(t)^i - _t^i-1
       ≤_t_t-δ_i(t)^i≤_tV_i

									where the inequalities in (<ref>) follow from (A1) and (A3). The bound _t-δ_i(t)^i≤ V_i follows from A3 and Jensen's inequality which implies that _t-δ_i(t)^i≤√(_t-δ_i(t)^i^2). Given 1≤ i,j ≤ K and t ≥ℓ≥ 1, it follows that
									
    _t^i-_ℓ^j   ≤∑_k=1^i_t^k-_t^k-1
       +∑_s=ℓ+1^t-1∑_k=1^K_s^k-_s^k-1 +∑_k=j+1^K_ℓ^k-_ℓ^k-1
       ≤_t ∑_k=1^i V_k+ (∑_s=ℓ+1^t-1_s)(∑_k=1^K V_k)+_ℓ∑_k=j+1^K V_k
       ≤_t (iV)+ (∑_s=ℓ+1^t-1_s)(KV)+_ℓ(K-j)V
       ≤_ℓ V[i+(t-ℓ-1)K+K-j]
       ≤_ℓ V[i-j+K(t-ℓ)]
where (<ref>) is obtained by substituting V = max_i V_i and (<ref>) follows since _ℓ≥_s for all ℓ≤ s ≤ t. The result in (<ref>) holds from the inequality (i-j)≤ |i-j|.  
									
									
									
									
									
									Further, for t≥ 1 and 1≤ i≤ K, let ^i_t be the σ-algebra generated by the random variables 

									
    {_1^1, …, _1^K, _2^1 …, _t-1^K, _t^1,…,_t^i}.

									where _t^i is the random state variables observed at node i at time t. 	With this definition, it holds that 
									
    ^i_t-δ_i(t)|ℱ_t-δ_i(t)^i-1=∂ D^i(_t-τ_i(t)^i-1)

									since _t-δ_i(t)^i may depend on _s^i-1 only for s ≤t-τ_i(t).

									
									 The proof is organized into two parts. Subsection <ref> develops a bound on the optimality gap in (<ref>), in terms of I_1. Subsequently, Subsection <ref> develops the required bound on I_1. 
										


  §.§.§ Bound on the optimality gap

An upper bound on _t^K-^2 is developed by making use of the form of the updates in (<ref>) for all nodes 1≤ i ≤ K. The bound follows from the use of triangle inequality and the moment bounds in (A3). Further, the bounded delay assumption (A4) enter through the use of (<ref>). 

Observe from the updates in (<ref>) that
										
    _t^i-^2    =  _t^i-1 -_t_t-δ_i(t)^i-^2 
       ≤_t^i-1 -_t_t-δ_i(t)^i-^2 
       = _t^i-1-^2 -2_t _t-δ_i(t)^i, _t^i-1- +_t^2 _t-δ_i(t)^i^2
       = _t^i-1-^2 -2_t _t-δ_i(t)^i, _t-τ_i(t)^i-1-
       +_t^2 _t-δ_i(t)^i^2-2_t _t-δ_i(t)^i, _t^i-1-_t-τ_i(t)^i-1

										where (<ref>) follows form (A1), and the term 2_t^i_t-δ_i(t), ^i-1_t-τ_i(t) has been added and subtracted to obtain (<ref>). Taking expectations on both sides and summing over all 1≤ i ≤ K and 1≤ t≤ T, we obtain
										
    _T^K-^2   = _1^0-^2 + ∑_t=1^T∑_i=1^K_t^2_t-δ_i(t)^i^2+ I_1 
       -2 ∑_t=1^T∑_i=1^K_t _t-δ_i(t)^i, _t-τ_i(t)^i-1-
where I_1 is as defined in Lemma <ref>. Deferring the bound on I_1 to Subsection <ref>, the last term in (<ref>) is analyzed first. In particular, it holds from (<ref>) that
										
    _t-δ_i(t)^i, _t-τ_i(t)^i-1-
       = _t-δ_i(t)^i |ℱ_t-δ_i(t)^i-1,_t-τ_i(t)^i-1-
       = ∂ D^i(_t-τ_i(t)^i-1),_t-τ_i(t)^i-1-.

										Further, since the functions D^i() are convex, it holds that
										
    -∂ D^i(_t-τ_i(t)^i-1),_t-τ_i(t)^i-1-≤ D^i() - D^i(_t-τ_i(t)^i-1) 
       = D^i() - D^i(_t^0) + D^i(_t^0) - D^i(_t-τ_i(t)^i-1) 
       ≤ D^i() - D^i(_t^0) + ∂ D^i(_t^0), _t^0 -_t-τ_i(t)^i-1
       ≤ D^i() - D^i(_t^0) + V_i_t^0-_t-τ_i(t)^i-1

										where (<ref>)-(<ref>) follow from the first order convexity condition for D^i and (<ref>) follows from the use of triangle inequality, and the fact that given any ∈Ł, 
										
    ∂ D^i() =  g_t^i()≤√(g_t^i()^2)≤ V_i.

										For the last term in (<ref>), taking expectation and utilizing the result in (<ref>), it follows that
										
    _t^0-_t-τ_i(t)^i-1   ≤_t-τ_i(t)[(i-1)V+KV(τ_i(t))]

										for 1≤ i≤ K. Finally, substituting (<ref>), (<ref>) in (<ref>), and using (A3), it follows that
										
    _T^K - ^2 ≤   _1^0-^2 +∑_t=1^T∑_i=1^K_t^2V_i^2 + I_1
       - 2∑_t=1^T∑_i=1^K_t[ D^i(_t^0) - D^i()]
    										
       +2∑_t=1^T∑_i=1^K_t_t-τ_i(t)V_i[(i-1)V+KV(τ_i(t))].

										Since the left-hand side is non-negative and _1^0 - ≤ B_0, the first part of Lemma <ref> is obtained simply by rearranging the terms in (<ref>) 
										
    2∑_t=1^T∑_i=1^K_t[ D^i(_t^0) - D^i()]    ≤ B_0^2 +  I_1 
       + ∑_t=1^T∑_i=1^K(_t^2 V_i^2 + 2_t_t-τ_i(t) V_i[(i-1)V+τ KV])
       ≤B_0^2 + I_0 + I_1
 
										where the first inequality in (<ref>) follows since τ_i(t)≤τ, and _t is non-increasing sequence. Finally, the second inequality in (<ref>) follows from substituting V_i ≤ V for all 1≤ i ≤ K, and I_0 is as defined in Lemma <ref>.
										
										


  §.§.§ Bound on I_1

In order to derive a bound on I_1, we make use of the Cauchy-Schwartz inequality as follows:

    _t-δ_i(t)^i, _t-τ_i(t)^i-1-_t^i-1   ≤ V_i_t^i-1 -_t-τ_i(t)^i-1
    
    										   ≤τ_i(t)_t-τ_i(t) V_iKV

										where (<ref>) follows from (<ref>). Consequently, 
										
    I_1    ≤∑_t=1^T∑_i=1^K2_t _t-τ_i(t)τ_i(t)V_iKV≤2τ KV∑_t=1^T∑_i=1^K V_i_t_t-τ_i(t)
 
										where (<ref>) utilizes the bounds τ_i(t)≤τ. Finally, substituting V = max_i V_i, we obtain I_1 ≤ 2τ KV^2∑_t=1^T∑_i=1^K_t _t-τ_i(t)										
										which is the required bound. 
										
									



§ PROOF OF LEMMA <REF>

Here, we establish that the dual iterates always stay bounded, thanks to the Slater's condition in (A6). The proof begins with establishing an upper bound on the per-iteration increase in the value of _t^0 - ^2, and subsequently utilizes an induction argument to derive the following bound for all t ≥ 1:



    _t   ≤ 2+max{_1,θ/C[𝖣-f^i(^i)]
       + θ KV^2/2C+ 2θτV̅/C+ VK}
where, θ and C are positive constants, Ṽ= V^2K(K-1), V̅= 2K^2V^2, and {^i} is a slater point of (<ref>). Since ^⋆ is bounded, the right hand side of  (<ref>) serves as the bound on _t .
									
									
									
									 In order to prove (<ref>), we will instead establish a more general result that takes the form:
	
    _t^0 - ≤   max{_1^0-,θ/C[𝖣-f^i(^i)]
       +θ KV^2/2C+2θτV̅/C++ VK}.

	where, recall that _t^0 = _t and _1^0 = _1. The desired result in (<ref>) will follow by applying the triangle inequality to (<ref>). The proof of (<ref>) follows via induction. It can be seen that the inequality in (<ref>) holds trivially for the base case of t=1. As part of the inductive hypothesis, assume that (<ref>) holds for t where t≥1. It remains to show that it also holds for t+1. We split the argument into the following two cases. 
										
										Case 1.  D(_t^0) > 𝖣+Ṽ/2+ τV̅: In this case, it holds that _t+1^0 - ^2 ≤_t^0-^2. Consequently, the induction hypothesis for time t implies that (<ref>) also holds for time t+1.
										
										Case 2.   D(_t^0) ≤𝖣+Ṽ/2+ τV̅: Recall that the dual function in (<ref>) is defined as
										
    D(_t^0)   =max_^i∈𝒳, _t^i∈𝒫_t[f^i(^i)+_t^0,[ _t^i(_t^i,^i)]] 
       ≥[f^i(^i)+_t^0,[ _t^i(_t^i,^i)]]
where {^i,{_t^i}_t≥ 1}_i=1^K is a strictly feasible (Slater) solution to (<ref>). From (A6), such a strictly feasible solution exists and satisfies [ _t^i(_t^i,^i)] > C > 0. Substituting into (<ref>), and rearranging, we obtain

										
    1,_t^0   ≤1/C[D(_t^0)- f^i(^i)]

										Since _t^0≽ 0, it follows from equivalence of norms _t^0≤θ_t^0_1 = θ1,_t^0. Therefore, taking expectations in (<ref>) yields
										
    _t^0   ≤θ/C[D(_t^0)- f^i(^i)] 
       ≤θ/C[𝖣+Ṽ/2+ τV̅- f^i(^i)]
where the assumption for Case 2 has been used in (<ref>). Finally, the use of triangle inequality and the bound in (<ref>) yields
										
    _t+1^0-   ≤_t^0 + _t+1^0-_t^0 + 
       ≤_t^0 +  VK +

										which, together with (<ref>), yields (<ref>) for t+1. Therefore by mathematical induction, the inequality in (<ref>) holds for all t≥ 1. Finally, using (<ref>) and triangle inequality, we obtain the result in (<ref>) since _t=_t^0.
 
									
									

§ PROOF OF LEMMA <REF>
 
									
										The proof establishes a lower bound for the running average of the primal objective function, calculated at the primal iterates. The lower bound depends upon the dual optimal value, dual initialization, and the maximum delay bound τ. For ease of  exposition, the proof begins with re-arranging the optimality gap in the form required by Lemma <ref> and subsequently analyzing the resulting terms. The full proof is split into various parts that develop separate bounds on the terms I_2, I_3, and I_4. Since (A7) is required to establish Lemma <ref>, the dual function D has to be differentiable. 
										
Since the functions f^i are concave, the expected value of the primal objective can be written as

										
    f^i(^i_T)≥1/T∑_t=1^T∑_i=1^K f^i(^i_t)
       =1/T∑_t=1^T∑_i=1^K[f^i(^i_t)+,∇ D^i()-,∇ D^i()] 
       = 1/T∑_t=1^T∑_i=1^K D^i() -1/T∑_t=1^T∑_i=1^K,∇ D^i()

										Consider the following expression where we simply add subtract  D^i(_t) as follows 					
    D^i()   = [D^i()+ D^i(_t)-D^i(_t)]
       = D^i(_t)+(D^i()-D^i(_t)) 
       ≥𝖣 + (D^i()-D^i(_t))

											where (<ref>) follows since 𝖣 =  D^i() ≤ D^i() for all ∈Ł. Taking the expectation on both sides of (<ref>) and substituting the result into (<ref>), we obtain 
										
    f^i(^i_T)≥   𝖣 - I_2 - I_3

										where I_2 and I_3 are as defined in Lemma <ref>. The rest of the proof proceeds simply by developing bounds on I_2 and I_3. 
										


  §.§.§ Bound on I_2
 The bound on I_2 follows simply from the moment bounds in (A3) and the Cauchy-Schwartz inequality. We begin with the following observation
					Since the functions D^i are convex, it holds that
										
    [D^i(_t)-D^i()]    ≤∇ D^i(_t), _t-
       ≤∇ D^i(_t)_t-
       ≤ V_i(i-1)V

										where (<ref>) uses the Cauchy-Schwartz inequality, while (<ref>) uses (<ref>) and (<ref>). Therefore, substituting (<ref>) into the expression for I_2 and rearranging, we obtain
										I_2 ≤V∑_i=1^K(i-1)V_i.
										Finally, the required bound in Lemma <ref> is obtained by substituting V = max_i V_i. 

	

  §.§.§ Bound on I_3

The bound in I_3 follows from setting aside the error due to asynchrony I_4, and developing a bound on the remaining terms by telescopically summing the bounds on _t^i over all 1≤ i≤ K and 1≤ t≤ T.

Since 0∈Λ is a feasible dual solution, using the form of the updates in (<ref>) and expanding as in (<ref>), it follows that
										
    ^2   ≤^2 +_t-δ_i(t)^i^2-2, _t-δ_i(t)^i.

										Adding the term 2_t^i-1,∇ D^i(_t^i-1) on both sides, and rearranging, we obtain
										
    2,∇ D^i()≤   ^2-^2+_t-δ_i(t)^i^2
       -2,_t,δ_i(t)^i

										where _t,δ_i(t)^i is as defined in (<ref>). Summing over i=1,…,K and  t=1,⋯,T, taking expectation, and utilizing (A3), it follows that
										
    I_3    ≤_1^2/2 T + /2 V_i^2 + I_4 ≤_1^2/2 T +  KV^2 /2 + I_4

										where I_4 is as defined in Lemma <ref> and the (<ref>) uses V = max_iV_i. 



  §.§.§ Bound on I_4

The term I_4 collects the error from the terms that arise due to asynchrony. A bound on I_4 is developed from the use of the delay bound assumption in (A4).	Adding and subtracting _t^i-1, ∇ D^i(_t-τ_i(t)^i-1) to each summand of I_4, we obtain
										
    I_4=   1/T∑_t=1^T∑_i=1^K,∇ D^i()-∇ D^i(_t-τ_i(t)^i-1)
       +∑_t=1^T∑_i=1^K, ∇ D^i(_t-τ_i(t)^i-1)-_t-δ_i(t)^i.
 
										Of these, the first term in (<ref>) can be bounded using the bound in Lemma <ref> and the Cauchy-Schwartz inequality, by observing that

										
    ,∇ D^i()-∇ D^i(_t-τ_i(t)^i-1)
       ≤ B∇ D^i()-∇ D^i(_t-τ_i(t)^i-1)
       ≤B L^i-_t-τ_i(t)^i-1
       ≤ B L^i τ KV

										where (<ref>) follows from (A7) and (<ref>) from the bound developed in (<ref>).
										
										For the second term in (<ref>), recalling the definition of _t-δ_i(t)^i-1 from Appendix <ref>, observe that although _t^i |_t-δ_i(t)^i-1≠_t^i, there exists some κ_i(t) ≤ t such that  
										
    _κ_i(t)^i-1|_t-δ_i(t)^i-1 = _κ_i(t)^i-1.

										Indeed, observe that κ_i(t) ≥ t-δ_i(t) since _t-τ_i(t)^i-1 only depends on random variables contained in _t-δ_i(t)^i-1. The subsequent bounds hold for any κ_i(t) that satisfies (<ref>), including for the worst case when κ_i(t) = t-δ_i(t). It follows that 										
    ∇ D^i(_t-τ_i(t)^i-1)-_t-δ_i(t)^i,_t^i-1
        = ∇ D^i(_t-τ_i(t)^i-1)-_t-δ_i(t)^i,_t^i-1|ℱ_t-δ_i(t)^i-1
        = [[∇ D^i(_t-τ_i(t)^i-1)-_t-δ_i(t)^i,_κ_i(t)^i-1|ℱ_t-δ_i(t)^i-1]] 
          + ∇ D^i(_t-τ_i(t)^i-1)-_t-δ_i(t)^i,_t^i-1 - _κ_i(t)^i-1.

										From (<ref>) and (<ref>), it follows that the first summand in (<ref>) is zero. The second summand can be bounded by using the Cauchy-Schwartz inequality and the bounds in (A4) and (<ref>) as follows:
										
    ∇ D^i(_t-τ_i(t)^i-1)-_t-δ_i(t)^i,_t^i-1 - _κ_i(t)^i-1
       ≤∇ D^i(_t-τ_i(t)^i-1)-_t-δ_i(t)^i_t^i-1 - _κ_i(t)^i-1
       ≤2V_i(t-κ_i(t)) KV
       ≤2V_iτKV

										where the inequality in (<ref>) follows since t-τ_i(t) ≤ t-δ_i(t) ≤κ_i(t) ≤ t. Finally, substituting (<ref>) and (<ref>) into (<ref>) yields
										
    I_4    ≤∑_i=1^K τ KV(B L^i+2V_i)≤τ K^2V(BL+2V)

										which together with (<ref>) gives the desired bound. 
									 
									
									



							

										
										IEEEtran 
										
										
									