{'fraction_non_alphanumeric': 0.05016329196603527, 'fraction_numerical': 0.01410842586544742, 'mean_word_length': 4.54965927214731, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Parallel Stochastic Approximation Method for Nonconvex Multi-Agent Optimization Problems

Consider the problem of minimizing the expected value of a (possibly nonconvex) cost function parameterized by a random (vector) variable, when the expectation cannot be computed accurately (e.g., because the statistics of the random variables are unknown and/or the computational complexity is prohibitive). Classical sample stochastic gradient methods for solving this problem may empirically suffer from slow convergence. In this paper, we propose for the first time a stochastic parallel Successive Convex Approximation-based (best-response) algorithmic framework for general nonconvex stochastic sum-utility optimization problems, which arise naturally in the design of multi-agent systems. The proposed novel decomposition enables all users to update their optimization variables in parallel by solving a sequence of strongly convex subproblems, one for each user. Almost surely convergence to stationary points is proved. We then customize our algorithmic framework to solve the stochastic sum rate maximization problem over Single-Input-Single-Output (SISO) frequency-selective interference channels, multiple-input-multiple-output (MIMO) interference channels, and MIMO multiple-access channels. Numerical results show that our algorithms are much faster than state-of-the-art stochastic gradient schemes while achieving the same (or better) sum-rates.


I. INTRODUCTION
Wireless networks are composed of users that may have different objectives and generate interference when no multiplexing scheme is imposed to regulate the transmissions; examples are peer-to-peer networks, cognitive radio systems, and ad-hoc networks. A usual and convenient way of designing such multi-user systems is by optimizing the "social function", i.e., the (weighted) sum of the users' objective functions. This formulation however requires the knowledge of the system parameters, which in practice is either difficult to acquire (e.g., when the parameters are rapidly changing) or imperfect due to estimation errors. In such scenarios, it is convenient to focus on the optimization of long-term performance of the system, measured as the expected value of the social function parametrized by the random system parameters. In this paper, we consider the frequent and difficult case wherein (the expected value of) the social function is nonconvex and the expectation cannot be computed (either numerically or in closed form). Such a system design naturally falls into the class of stochastic optimization [2,3].
Gradient methods for unconstrained stochastic nonconvex optimization problems have been studied in [4,5,6], where almost sure convergence to stationary points has been established, under some technical conditions; see, e.g., [5]. The extension of these methods to constrained optimization problems is not straightforward; in fact, the descent-based convergence analysis developed for unconstrained gradient methods no longer applies to their projected counterpart (due to the presence of the projection operator). Convergence of stochastic gradient projection methods has been proved only for convex objective functions [4,7,8].
To cope with nonconvexity, gradient averaging seems to be an essential step to resemble convergence; indeed, stochastic conditional gradient methods for nonconvex constrained problems hinge on this idea [9,10,11,12]: at each iteration the new update of the variables is based on the average of the current and past gradient samples. Under some technical conditions, the average sample gradient eventually resembles the nominal (but unavailable) gradient of the (stochastic) objective function [9,13]; convergence analysis can then borrow results from deterministic nonlinear programming.
Numerical experiments for large classes of problems show that plain gradient-like methods usually converge slowly and are very sensitive to the choice of the step-size. Some acceleration techniques have been proposed in the literature [8,14], but only for strongly convex objective functions. Here we are interested in nonconvex (constrained) stochastic problems. Moreover, (proximal, accelerated) stochastic gradientbased schemes use only the first order information of the objective function (or its realizations); recently it was shown [15,16,17] that for deterministic nonconvex optimization problems exploiting the structure of the function by replacing its linearization with a "better" approximant can enhance empirical convergence speed. In this paper we aim at bringing this idea into the context of stochastic optimization problems.
Our main contribution is to develop a new, broad algorithmic framework for the computation of stationary solutions of a wide class of (stochastic) nonconvex optimization problems, encompassing many multi-agent system designs of practical interest. The essential idea underlying our approach is to decompose the original nonconvex stochastic problem into a sequence of (simpler) deterministic subproblems whereby the objective function is replaced by suitable chosen sample convex approximations; the subproblems can be then solved in a parallel and distributed fashion across the users. Other key features of our framework are: i) no knowledge of the objective function parameters (e.g., the Lipschitz constant of the gradient) is required; ii) it is very flexible in the choice of the approximant of the nonconvex objective function, which need not be necessarily its first or second order approximation (like in proximal-gradient schemes); of course it includes, among others, updates based on stochastic gradient-or Newtontype approximations; iii) it can be successfully used also to robustify distributed iterative algorithms solving deterministic social problems, when only inexact estimates of the system parameters are available; and iv) it encompasses a gamut of novel algorithms, offering a wide flexibility to control iteration complexity, communication overhead, and convergence speed, while converging under the same conditions. These desirable features make our schemes applicable to several different problems and scenarios. As illustrative examples, we customize our algorithms to some resource allocations problems in wireless communications, namely: the sum-rate maximization problems over MIMO Interference Channels (ICs) and Multiple Access Channels (MACs). The resulting algorithms outperform existing (gradient-based) methods both theoretically and numerically.
The proposed decomposition technique hinges on successive convex approximation (SCA) methods, and it is a nontrivial generalization of [15] to stochastic optimization problems. An SCA framework for stochastic optimization problems has also been proposed in a recent submission [18] 1 ; however our method differs from [18] in many features. First of all, we relax the key requirement that the convex approximation must be a tight global upper bound of the (sample) objective function, as required instead in [18]. This represents a turning point in the design of distributed stochastic SCA-based methods, enlarging substantially the class of (large scale) stochastic nonconvex problems solvable using our framework. Second, even when the aforementioned constraint can be met, it is not always guaranteed that the resulting convex (sample) subproblems are decomposable across the users, implying that a centralized implementation might be required in [18]; our schemes instead naturally lead to a parallel and distributed implementation. Third, the proposed methods converge under weaker conditions than those in [18].
Finally, within the classes of approximation-based methods for stochastic optimization problems, it is worth mentioning the so-called Sample Average Approach (SAA) [18,19,20,21]: the "true" (stochastic) objective function is approximated by an ensemble average. Then the resulting deterministic SSA optimization problem has to be solved by an appropriate numerical procedure. When the original objective function is nonconvex, however, computing the global optimal solution of the SAA at each step may not be easy, if not impossible. Therefore SSA-based methods are generally used to solve stochastic convex optimization problems only.
The rest of the paper is organized as follows. Sec. II formulates the problem along with some interesting applications. The novel stochastic decomposition framework is introduced in Sec. III; customizations of the main algorithms to sample applications are discussed in Sec. IV. Finally, Sec. V draws some conclusions.

II. PROBLEM FORMULATION
We consider the design of a multi-agent system composed of I users; each user i has his own strategy vector x i to optimize, which belongs to the feasible convex set X i ⊆ C ni . The variables of the other users are denoted by x −i (x j ) I j=1,j =i , and the joint strategy set of all users is X = X 1 × . . . × X I . 1 A preliminary version of our work appeared independently before [18] at IEEE SPAWC 2013 [1].
The stochastic social optimization problem is formulated as: where I f {1, . . . , I f }, with I f being the number of functions; each cost function f j (x, ξ) : X × D → R depends on the joint strategy vector x and a random vector ξ, whose probability distribution is defined on a set D ⊆ C m ; and the expectation is taken with respect to (w.r.t.) ξ. Note that the optimization variables can be complex-valued; in such a case, all the gradients of real-valued functions are intended to be conjugate gradients [22,23]. Assumptions: We make the following blanket assumptions: (a) Each X i is compact and convex; (b) Each f j (•, ξ) is continuously differentiable on X , for any given ξ, and the gradient is Lipschitz continuous with constant L ∇fj (ξ) . Furthermore, the gradient of U (x) is Lipschitz continuous with constant L ∇U < +∞. These assumptions are quite standard and are satisfied by a large class of problems. Note that the existence of a solution to (1) is guaranteed by Assumption (a). Since U (x) is not assumed to be jointly convex in x, (1) is generally nonconvex. Some instances of (1) satisfying the above assumptions are briefly listed next. Example #1: Consider the maximization of the ergodic sumrate over frequency-selective ICs: maximize p1,...,pI E N n=1 I i=1 log 1 + |hii,n| 2 pi,n σ 2 i,n + j =i |hij,n| 2 pj,n subject to p i ∈ P i {p i : with p i,n being the transmit power of user i on subchannel (subcarrier) n, N is the number of parallel subchannels, P i is the total power budget, h ij,n is the channel coefficient from transmitter j to receiver i on subchannel n, and σ 2 i,n is the variance of the thermal noise over subchannel n at the receiver i. The expectation is over channel coefficients (h ij,n ) i,j,n . Example #2: The following maximization of the ergodic sumrate over MIMO ICs also falls into the class of problems (1): ij is the covariance matrix of the thermal noise R N i (assumed to be full rank) plus the multi-user interference, P i is the total power budget, and the expectation in (3) is taken over the channels H (H ij ) I i,j=1 . Example #3: Another application of interest is the maximization of the ergodic sum-rate over MIMO MACs: This is a special case of (1) where the utility function is concave in Q (Q i ) I i=1 , I f = 1, I f = {1}, and the expectation in (4) is taken over the channels H (H i ) I i=1 . Example #4: The algorithmic framework that will be introduced shortly can be successfully used also to robustify distributed iterative algorithms solving deterministic (nonconvex) social problems, but in the presence of inexact estimates of the system parameters. More specifically, consider as example the following sum-cost minimization multi-agent problem: An efficient distributed algorithm converging to stationary solutions of (5) has been recently proposed in [15]: at each iteration t, given the current iterate x t , every agent i minimizes (w.r.t. x i ∈ X i ) the following convexified version of the social function: where a, a ). The evaluation of the above function requires the exact knowledge of ∇ i f j (x t ) for all j = i. In practice, however, only a noisy estimate of ∇ i f j (x t ) is available [24,25,26]. In such cases, convergence of pricingbased algorithms [15,27,28,29] is in jeopardy. We will show in Sec. IV-C that the proposed framework can be readily applied, for example, to robustify (and make convergent), e.g., pricing-based schemes, such as [15,27,28,29].
Since the class of problems (1) is in general nonconvex (possibly NP hard [30]), the focus of this paper is to design distributed solution methods for computing stationary solutions (possibly local minima) of (1). Our major goal is to devise parallel (nonlinear) best-response schemes that converge even when the expected value in (1) cannot be computed accurately and only sample values are available.

III. A NOVEL PARALLEL STOCHASTIC DECOMPOSITION
The social problem (1) faces two main issues: i) the nonconvexity of the objective functions; and ii) the impossibility to estimate accurately the expected value. To deal with these difficulties, we propose a decomposition scheme that consists in solving a sequence of parallel strongly convex subproblems (one for each user), where the objective function of user i is obtained from U (x) by replacing the expected value with a suitably chosen incremental estimate of it and linearizing the nonconvex part. More formally, at iteration t + 1, user i solves the following problem: given x t −i and ξ t , with the approximation functionf i (x i ; x t , ξ t ) defined aŝ where the pricing vector π i (x, ξ) is given by and f t i is an accumulation vector updated recursively according to with ρ t ∈ (0, 1] being a sequence to be properly chosen (ρ 0 = 1). The other symbols in (6) are defined as follows: where γ t ∈ (0, 1]. It turns our that x t is a random vector depending on F t , the past history of the algorithm up to iteration t: depends on F t as well (we omit this dependence for notational simplicity).
The subproblems (6a) have an interesting interpretation: each user solves a sample convex approximation of the original nonconvex stochastic function. The first term in (6b) preserves the convex component (or a part of it, if C t i ⊂ S t i ) of the (instantaneous) social function. The second term in (6b)−the pricing vector π i (x, ξ)−comes from the linearization of (at least) the nonconvex part. The vector f t i in the third term represents the incremental estimate of ∇ x * U (x t ) (whose value is not available), as one can readily check by substituting (6c) into (6d): Roughly speaking, the goal of this third term is to estimate onthe-fly the unknown ∇ x * U (x t ) by its samples collected over the iterations; based on (9), such an estimate is expected to become more and more accurate as t increases, provided that the sequence ρ t is properly chosen (this statement is made rigorous shortly in Theorem 1). The last quadratic term in (6b) is the proximal regularization whose numerical benefits are well-understood [31]. Given (6), we define the "best-response" mapping as Note that x(•, ξ) is well-defined for any given ξ because the objective function in (6) is strongly convex with constant τ min : Algorithm 1 : Stochastic parallel decomposition algorithm Our decomposition scheme is formally described in Algorithm 1, and its convergence properties are stated in Theorem 1, under the following standard boundedness assumptions on the instantaneous gradient errors [24,32]: The instantaneous gradient is unbiased with bounded variance in the following sense: This assumption is readily satisfied if ξ is a bounded i.i.d. random variable.

Theorem 1. Given problem (1) under Assumptions (a)-(c)
, suppose that τ min > 0 and the stepsizes {γ t } and {ρ t } are chosen so that Then, every limit point of the sequence {x t } generated by Algorithm 1 (at least one of such point exists) is a stationary point of (1) almost surely.
Proof: See Appendix A. On Assumption (c): The boundedness condition is in terms of the conditional expectation of the (random) gradient error. Compared with [18], Assumption (c) is weaker because it is required in [18] that every realization of the (random) gradient error must be bounded. On Condition (12d): The condition has the following interpretation: all increasing subsequences of j∈I f L ∇fj (ξ t ) must grow slower than 1/ρ t . We will discuss later in Sec. IV how this assumption is satisfied by specific applications. Note that if j∈I f L ∇fj (ξ) is uniformly bounded for any ξ (which is indeed the case if ξ is a bounded random vector), then (12d) is trivially satisfied. On Algorithm 1: To our best knowledge, Algorithm 1 is the first parallel best-response (e.g., nongradient-like) scheme for nonconvex stochastic social problems: all the users update in parallel their strategies (possibly with a memory) solving a sequence of decoupled (strongly) convex subproblems (6). It is expected to perform better than classical stochastic gradientbased schemes at no the cost of extra signaling, because the convexity of the objective function, if any, is better exploited. Our experiments on specific applications confirm this intuition; see Sec. IV. Moreover, it is guaranteed to converge under the weakest assumptions available in literature while offering some flexibility in the choice of the free parameters [cf. Theorem 1]. Diminishing stepsize rules: In order to have convergence, a diminishing stepsize rule satisfying (12) is necessary. An instance of (12) is, e.g., the following: Roughly speaking, (12) says that the stepsizes γ t and ρ t , while diminishing (with γ t decreasing faster than ρ t ), need not go to zero too fast. This kind of stepsize rules are of the same spirit of those used to guarantee convergence of gradient methods with error; see [33] for more details. Implementation issues: In order to compute the bestresponse, each user needs to know The signaling required to acquire this information is of course problem-dependent. If the problem under consideration does not have any specific structure, the most natural message-passing strategy is to communicate directly x t −i and π i (x t , ξ t ). However, in many specific applications much less signaling may be needed; see Sec. IV for some examples. Note that the signaling is of the same spirit of that of pricing-based algorithms proposed in the literature for the maximization of deterministic sum-utility functions [15,29]; no extra communication is required to update f t i : once the new pricing vector π i (x t , ξ t ) is available, the recursive update (6d) for the "incremental" gradient is based on a local accumulation register keeping track of the last iterate f t−1 i . Note also that, thanks to the simultaneous nature of the proposed scheme, the overall communication overhead is expected to be less than that required to implement sequential schemes, such as [29].

A. Some special cases
We customize next the proposed general algorithmic framework to specific classes of problems (1) arising naturally in many applications.
1) Stochastic proximal conditional gradient methods: Quite interestingly, the proposed decomposition technique resembles classical stochastic conditional gradient schemes [4] when one chooses in (6b) C t i = ∅, for all i and t, resulting in the following approximation function: with f t i updated according to (9). Note that, however, traditional stochastic conditional gradient methods [9] do not have the proximal regularization term in (14), which instead brings in well-understood numerical benefits. Moreover, it is worth mentioning that, for some of the applications introduced in Sec. II, it is just the presence of the proximal term that allows one to compute the best-responsex i (x t , ξ t ) resulting from the minimization of (14) in closed-form; see Sec. IV-B.
It turns out that convergence conditions of Algorithm 1 contain as special cases those of classical stochastic conditional gradient methods. But the proposed algorithmic framework is much more general and, among all, is able to better exploit the structure of the sum-utility function (if any) than just linearizing everything; it is thus expected to be faster than classical stochastic conditional gradient methods, a fact that is confirmed by our experiments; see Sec. IV.
2) Stochastic best-response algorithm for single (convex) functions: Suppose that the social function is a single func- and  f (x 1 , . . . , x I , ξ) is uniformly convex in each x i (but not necessarily jointly). Of course, this optimization problem can be interpreted as a special case of the framework (1), with is already convex separately in the variables x i 's, a natural choice for the approximantsf i is setting C t i = S t i = {1} for all t, resulting in the following: . This is the generalization of the deterministic pricing algorithms [15,29] to stochastic optimization problems. Examples of this class of problems are the ergodic sum-rate maximization problem over SISO and MIMO IC in (2) where both f j (•, ξ) and g j (•, ξ) are uniformly convex functions on X . A natural choice of the approximation functionŝ f i for (17) is linearizing the concave part of the sample social function, resulting in the following:

IV. APPLICATIONS
We customize now the proposed algorithmic framework to some of the applications introduced in Sec. II, and compare the resulting algorithms with state-of-the-art schemes proposed for the specific problems under considerations as well as more classical stochastic gradient algorithms. Numerical results provide a solid evidence of the superiority of our approach.

A. Sum-rate maximization over frequency-selective ICs
Consider the sum-rate maximization problem over frequency-selective ICs, as introduced in (2). Since the instantaneous rate of each user i, is uniformly strongly concave in p i ∈ P i , a natural choice for the approximation functionf i is the one in (16) wherein linearized. This leads to the following best-response functionŝ SINR t j,n = |h t jj,n | 2 p t j,n /MUI t j,n . Then the variable f t i is updated according to . Note that the optimization problem in (18a) has a closed-form expression [15]: where and µ ⋆ is the Lagrange multiplier such that 0 ≤ µ ⋆ ⊥ N n=1p i,n (p t , h t ) − P i ≤ 0, and it can be found by standard bisection method.
The overall stochastic pricing-based algorithm is then given by Algorithm 1 with best-response mapping defined in (19); convergence is guaranteed under conditions i)-iv) in Theorem 1. Note that the theorem is trivially satisfied using stepsizes rules as required in i)-iii) [e.g., (13)]; the only condition that needs some comment is condition iv). If lim sup t→∞ ρ t j∈I f L ∇fj (ξ t ) > 0, we can assume without loss of generality (w.l.o.g.) that the sequence of the Lipschitz constant j∈I f L ∇fj (ξ t ) is increasing monotonically at a rate no slower than 1/ρ t (we can always limit the discussion to such a subsequence). For anyh > 0, define p(h) Prob(|h ij,n | ≥h) and assume w.l.o.g. that 0 ≤ p(h) < 1. Note that the Lipschitz constant L ∇fj (ξ) is upper bounded by the maximum eigenvalue of the augmented Hessian of f j (x, ξ) [35], and the maximum eigenvalue increasing monotonically means that the channel coefficient is becoming larger and larger (this can be verified by explicitly calculating the augmented Hessian of f j (x, ξ); detailed steps are omitted due to page limit). Since Prob(|h t+1 ij,n | ≥ |h t ij,n | for all t ≥ t 0 ) ≤ Prob(|h t+1 ij,n | ≥h for all t ≥ t 0 ) = p(h) t−t0+1 → 0, we can infer that the magnitude of the channel coefficient increasing monotonically is an event of probability 0. Therefore, condition (12d) is satisfied. Numerical results. We simulated a SISO frequency selective IC under the following setting: the number of users is either five or twenty; equal power budget P i = P and white Gaussian noise variance σ 2 i = σ 2 are assumed for all users; the SNR of each user snr = P/σ 2 is set to 10dB; the instantaneous parallel subchannels h t (h t ij,n ) i,j,n are generated according to h t = h + △h t , where h (generated by MATLAB command randn) is fixed while △h t is generated at each t using δ · randn, with δ = 0.2 being the noise level. We simulated the following algorithms: i) the proposed stochastic best-response pricing algorithm (with τ i = 10 −8 for all i, γ 1 = ρ 0 = ρ 1 = 1, ρ t = 2/(t + 2) 0.6 , and γ t = 2/(t + 2) 0.61 for t ≥ 2); ii) the stochastic conditional gradient method [9] (with γ 1 = ρ 0 = ρ 1 = 1, ρ t = 1/(t + 2) 0.9 , and γ t = 1/(t + 2) 0.91 for t ≥ 2); iii) and the stochastic gradient projection method [26] (with γ 1 = 1 and γ t = γ t−1 (1 − 10 −3 γ t−1 ) for t ≥ 2). Note that the stepsizes are tuned such that all algorithms can achieve their best empirical convergence speed. We plot two merit functions, namely: i) the ergodic sum-rate, defined as E h [ N n=1 I i=1 r i (p t , h)] (with the expected value estimated by the sample mean of 1000 independent realizations); and ii) the "achievable" sum-rate, defined as 1 h m ), which represents the sum-rate that is actually achieved in practice (it is the time average of the instantaneous (random) sum-rate).
In Figure 1, we plot the two above merit functions versus the iteration index t achieved using the different algorithms. Our experiment show that for "small" systems (e.g., five active users), all algorithms perform quite well (for both the merit functions), with a gain in convergence speed for the proposed scheme. However, when the number of users increases increased (e.g., from 5 to 20), all other (gradientlike) algorithms suffer from very slow convergence. Quite interestingly, the proposed scheme seems also quite scalable: the convergence speed is not notably affected by the number of users, which makes it applicable to more realistic scenarios. The faster convergence of proposed stochastic best-response pricing algorithm comes from a better exploitation of partial convexity in the problem than what more classical gradient algorithms do, which validates the main idea of this paper.

B. Sum-rate maximization over MIMO ICs
Now we customize Algorithm 1 to solve the sum-rate maximization problem over MIMO ICs (3). Defining and following a similar approach as in the SISO case, the best-response of each user i becomes [cf. (16)]: where A, B tr(A H B); Π i (Q, H) is given by We can then apply Algorithm 1 based on the best-responsê Q(Q t , H t ) = (Q i (Q t , H t )) I i=1 whose convergence is guaranteed if the stepsizes are properly chosen [cf. Theorem 1].
Differently from the SISO case, the best-response in (20a) does not have a closed-form solution. A standard option to computeQ(Q t , H t ) is using standard solvers for strongly convex optimization problems. By exploiting the structure of problem (20), we propose next an efficient iterative algorithm converging toQ(Q t , H t ), wherein the subproblems solved at each step have a closed-form solution. Second-order dual method. To begin with, for notational simplicity, we rewrite (20a) in the following general form: whereX U H XU,Ã U H AU, andX = U HX U. We now partition D 0 in two blocks, its positive definite and zero parts, andX accordingly: where D 1 ≻ 0, andX 11 and D 1 have the same dimensions. Problem (22) can be then rewritten as: Note that, sinceX ∈ Q, by definitionX 11 must belong to Q as well. Using this fact and introducing the slack variable Now we solve (24) via dual decomposition (note that there is zero duality gap). The (partial) Lagrangian function of (24) is: denoting by Z the matrix of multipliers associated to the linear constraints Y =X 11 , The dual problem is then Problem (25) is quadratic and has a closed-form solution (see Lemma 2 below). Similarly, if Z ≺ 0, (26) can be solved in closed-form, up to a Lagrange multiplier which can be found by bisection; see, e.g., [29, Table I]. In our setting, however, Z in (26) is not necessarily negative definite. Nevertheless, the next lemma provides a closed form expression of Y(Z) [andX(Z)]. (25) and (26) in the setting above, the following hold: (25) is given bỹ
One can then solve the dual problem using standard (proximal) gradient-based methods; see, e.g., [35]. As a matter of fact, d(Z) is twice continuously differentiable, whose augmented Hessian matrix [22] is given by [35,Sec. 4.2.4]: and ∇ 2 X11X * 11 L(XY, Z) = −τ I. Since D 1 ≻ 0, it is easy to verify that ∇ 2 ZZ * d(Z) ≻ 0 and the following second-order Newton's method to update the dual variable Z is well-defined: . The convergence speed of the Newton's methods is typically very fast, and, in particular, superlinear convergence rate can be expected when Z t is close to Z ⋆ [35, Prop. 1.4.1].
As a final remark on efficient solution methods computinĝ Q i (Q t , H t ), note that one can also apply the proximal conditional gradient method as introduced in (14), which is based where µ ⋆ is the Lagrange multiplier that can be found efficiently by bisection method. Note that (29) differs from more traditional conditional stochastic gradient methods [9] by the presence of the proximal regularization, thanks to which one can solve (29) in closed form. Practical implementations. The proposed algorithm is fairly distributed: once the pricing matrix Π i is given, to compute the best-response, each user only needs to locally estimate the covariance matrix of the interference plus noise. Note that both the computation ofQ i (Q, H) and the update of F i can be implemented locally by each user. The estimation of the pricing matrix Π i requires however some signaling among nearby receivers. Quite interestingly, the pricing expression and thus the resulting signaling overhead necessary to compute it coincide with [29] (where a sequential algorithm is proposed for the deterministic maximization of the sum-rate over MIMO ICs) and the stochastic gradient projection method in [26]. We remark that the signaling to compute (20b) is lower than in [18,36], wherein signaling exchange is required twice (one in the computation of U i and another in that of A i ; see [36, Algorithm 1] for more details) in a single iteration to transmit among users the auxiliary variables which are of same dimensions as Π i . Numerical Results. We considered the same scenario as in the SISO case (cf. Sec. IV-A) with the following differences: i) there are 50 users; ii) the channels are matrices generated according to H t = H+△H t , where H is given while △H t is realization dependent and generated by δ · randn, with noise level δ = 0.2; and iii) the number of transmit and receive antennas is four. We simulate the following algorithms: i) the proposed stochastic best-response pricing algorithm (20) (with τ i = 10 −8 for all i; γ 1 = ρ 0 = ρ 1 = 1 and ρ t = 2/(t + 2) 0.6 and γ t = 2/(t + 2) 0.61 for t ≥ 2); ii) the proposed stochastic proximal gradient method (29) with τ = 0.01 and same stepsize as stochastic best-response pricing algorithm; iii) the stochastic conditional gradient method [9] (with γ 1 = ρ 0 = ρ 1 = 1 and ρ t = 1/(t + 2) 0.9 and γ t = 1/(t + 2) 0.91 for t ≥ 2); and iv) the stochastic weighted minimum meansquare-error (SWMMSE) method [36]. The best-response of the algorithm in i) is computed using Lemma 2. We observed convergence of the inner loop solving (25)- (26) in a very few iterations. Similarly to the SISO ICs case, we consider both ergodic sum-rate and achievable sum-rate. In Figure 2 we plot both objective functions versus the iteration index. It is clear from the figures that the proposed best-response pricing and proximal gradient algorithms outperform current schemes in terms of both convergence speed and achievable (ergodic or instantaneous) sum-rate. Note also that the bestresponse pricing algorithm is very scalable compared with the other algorithms. Finally, it is interesting to note that the proposed stochastic proximal gradient algorithm outperforms the conditional stochastic gradient method in terms of both convergence speed and iteration complexity. This is mainly due to the presence of the proximal regularization.

C. Sum-rate maximization over MIMO MACs
We consider now the sum-rate maximization problem over MIMO MACs, as defined in (4). Define A natural choice for the best-response of each user i is [cf. (15)]: and F t i is updated as any H, the ergodic sum-rate function is concave in Q i 's, and thus Algorithm 1 will converge (in the sense of Theorem 1) to the global optimal solution of (4). To the best of our knowledge, this is the first example of stochastic approximation algorithm based on best-response dynamics rather than gradient responses. Numerical results. We compare the proposed best-response method (30) (whose solution is computed using Method #2 in Sec. IV-B) with the stochastic conditional gradient method [9], and the stochastic gradient projection method [8]. System parameters (including the stepsize rules) are set as for the MIMO IC example in Sec. IV-B. In Figure 3 we plot both the ergodic sum-rate and the achievable sum-rate versus the iteration index. The figure clearly shows that Algorithm 1 outperforms the conditional gradient method and the gradient projection method in terms of convergence speed, and the performance gap is increasing as the number of users increases. This is because the proposed algorithm is a best-response type scheme, which thus explores the concavity of each users' rate function better than what gradient methods do. Note also the good scalability of the proposed best-response method.

D. Distributed deterministic algorithms with errors
The developed framework can be useful in the context of deterministic optimization as well to robustify best-responsebased algorithms in the presence of noisy estimates of system parameters. Consider the deterministic optimization problem introduced in (5). The main iterate of the deterministic counterpart of Algorithm 1 is still given by (7) but with each x i (x t ) defined as [15] x i (x t ) = arg min where π i (x) = j∈Ci ∇ i f j (x). In many applications (see, e.g., [24,25,26]), however, only a noisy estimate, denoted by π i (x), is available instead of π i (x). A heuristic is then to replace in (31) the exact π i (x) with its noisy estimate π i (x). The limitation of this approach, albeit natural, is that convergence of the resulting scheme is in jeopardy.
If [24,25], capitalizing on the proposed framework, we can readily deal with estimation errors while guaranteeing convergence. In particular, it is sufficient to modify (31) as follows: where f t i is updated according to . Algorithm 1 based on the best-response (32) is then guaranteed to converge to a stationary solution of (5), in the sense specified by Theorem 1.
As a case study, we consider next the maximization of the deterministic sum-rate over MIMO ICs in the presence of pricing estimation errors: subject to Q i 0, tr(Q i ) ≤ P i , i = 1, . . . , I.
Then (32) becomes: Convergence w.p.1 to a stationary point of the deterministic optimization problem (33) is guaranteed by Theorem 1. Note that if the channel matrices {H ii } are full column-rank, one can also set in (34) all τ i = 0, and compute (34) in closed form (cf. Lemma 2). Numerical results. We consider the maximization of the deterministic sum-rate (33) over a 5-user MIMO IC.  Sec. IV-B and then only its Hermitian part is kept; the noise level δ is set to 0.05. We compare the following algorithms: i) the proposed robust pricing method−Algorithm 1 based on the best-response defined in (34); and ii) the plain pricing method as proposed in [15] [cf. (31)]. We also include as a benchmark the sum-rate achieved by the plain pricing method (31) when there is no estimation noise (i.e., perfect π i (x) is available). In Figure 4 we plot the deterministic sum-rate in (33) versus the iteration index t. As expected, Figure 4 shows that the plain pricing method [15] is not robust to pricing estimation errors, whereas the proposed robustification preforms quite well. For instance, the rate achievable by the proposed method is about 50% larger than the one of [15], and is observed to reach the benchmark value (achieved by the plain pricing method when there is no estimation noise). This is due to the fact that the proposed robustification filters out the estimation noise. Note that the limit point generated by the proposed scheme is a stationary solution of the deterministic problem (33).

V. CONCLUSIONS
In this paper, we have proposed a novel best-response-based algorithmic framework converging to stationary solutions of general stochastic nonconvex optimization problems. The proposed novel decomposition enables all users to update their optimization variables in parallel by solving a sequence of strongly convex subproblems; which makes the algorithm very appealing for the distributed design of several multi-agent systems. We have then customized the general framework to solve special classes of problems as well as specific applications, including the stochastic maximization of the sum-rate over frequency-selective ICs, MIMO ICs and MACs. Extensive experiments have provided a solid evidence of the superiority in terms of both achievable sum-rate and practical convergence of the proposed schemes with respect to to state-of-the-art stochastic-based algorithms.

A. Proof of Theorem 1
We first introduce the following two preliminary results. Lemma 3. Given problem (1) under Assumptions (a)-(c), suppose that the stepsizes {γ t } and {ρ t } are chosen according to (12). Let {x t } be the sequence generated by Algorithm 1. Then, the following holds lim t→∞ f t − ∇U (x t ) = 0, w.p.1.
Proof: This lemma is a consequence of [10, Lemma 1]. To see this, we just need to verify that all the technical conditions therein are satisfied by the problem at hand. Specifically, Condition (a) of [10, Lemma 1] is satisfied because X i 's are closed and bounded in view of Assumption (a). Condition (b) of [10, Lemma 1] is exactly Assumption (c). Conditions (c)-(d) come from the stepsize rules i)-ii) in (12) of Theorem 1. Condition (e) of [10, Lemma 1] comes from the Lipschitz property of ∇U from Assumption (b) and stepsize rule iii) in (12) of Theorem 1.