Coordinated Beamforming for Multiuser MISO Interference Channel under Rate Outage Constraints
aa r X i v : . [ c s . I T ] N ov ACCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 1
Coordinated Beamforming for Multiuser MISOInterference Channel under Rate OutageConstraints § Wei-Chiang Li ⋆ , Tsung-Hui Chang † , Che Lin ⋆ , and Chong-Yung Chi ⋆ Abstract
This paper studies the coordinated beamforming design problem for the multiple-input single-output(MISO) interference channel, assuming only channel distribution information (CDI) at the transmitters.Under a given requirement on the rate outage probability for receivers, we aim to maximize the systemutility (e.g., the weighted sum rate, weighted geometric mean rate, and the weighed harmonic meanrate) subject to the rate outage constraints and individual power constraints. The outage constraints,however, lead to a complicated, nonconvex structure for the considered beamforming design problemand make the optimization problem difficult to handle. Although this nonconvex optimization problem canbe solved in an exhaustive search manner, this brute-force approach is only feasible when the numberof transmitter-receiver pairs is small. For a system with a large number of transmitter-receiver pairs,computationally efficient alternatives are necessary. The focus of this paper is hence on the design ofsuch efficient approximation methods. In particular, by employing semidefinite relaxation (SDR) andfirst-order approximation techniques, we propose an efficient successive convex approximation (SCA)algorithm that provides high-quality approximate beamforming solutions via solving a sequence of convexapproximation problems. The solution thus obtained is further shown to be a stationary point for the SDRof the original outage constrained beamforming design problem. Furthermore, we propose a distributedSCA algorithm where each transmitter optimizes its own beamformer using local CDI and informationobtained from limited message exchange with the other transmitters. Our simulation results demonstratethat the proposed SCA algorithm and its distributed counterpart indeed converge, and near-optimalperformance can be achieved for all the considered system utilities.
Index terms − Interference channel, coordinated beamforming, outage probability, convex optimization,semidefinite relaxation.
EDICS : SAM-BEAM, MSP-APPL, MSP-CODR, SPC-APPL
Copyright (c) 2012 IEEE. Personal use of this material is permitted. However, permission to use this material for any otherpurposes must be obtained from the IEEE by sending a request to [email protected]. § The work is supported by the National Science Council, R.O.C., under Grant NSC 99-2221-E-007-052-MY3 and underGrant NSC 101-2218-E-011-043. Part of this work was presented at the IEEE ICASSP, Prague, Czech, May 22-27, 2011 [1]. † Tsung-Hui Chang is the corresponding author. Address: Department of Electronic and Computer Engineering, NationalTaiwan University of Science and Technology, Taipei, Taiwan 10607, R.O.C. E-mail: [email protected]. ⋆ Wei-Chiang Li, Che Lin and Chong-Yung Chi are with Institute of Communications Engineering & Department ofElectrical Engineering, National Tsing Hua University, Hsinchu, Taiwan 30013, R.O.C. E-mail: [email protected], { clin,cychi } @ee.nthu.edu.tw March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 2
I. I
NTRODUCTION
Inter-cell interference is known to be one of the main bottlenecks that limit the system performanceof a wireless cellular network where all transmitters share a universal frequency band. The performancedegradation caused by such interference is severe especially for the users at the cell edge and can onlybe alleviated when some sort of cooperation is available between base stations (BSs) [2]. According tothe level of cooperation, the coordinated transmission can be roughly divided into two classes: Networkmultiple-input multiple-output (MIMO) and interference coordination [3]. In network MIMO, all BSswork as a single virtual BS using all the available antennas for data transmission and reception. Each ofthe BSs requires to know all the channel state information (CSI) and data streams of users, demanding alarge amount of message exchange between BSs [4]. Interference coordination, by contrast, only needsCSI sharing between BSs; based on the shared CSI, the BSs coordinate with each other in the designof transmission strategies, e.g., coordinated beamforming [5], [6] or power allocation [7]. Our interestin this paper lies in the coordinated beamforming design. To this end, we adopt the commonly usedinterference channel (IFC) model [8]–[10]. Under this model, a Pareto optimal transmission scheme isthat the rate tuple of receivers resides on the boundary of the achievable rate region [11]. It is alwaysdesirable to have a Pareto optimal transmission scheme since, otherwise, the achievable rates of some ofthe receivers can be further improved.Consider a multiple-input single-output (MISO) IFC, where the transmitters are equipped with multipleantennas while the receivers, i.e., mobile users, have only single antenna. We assume that the receiversemploy single-user detection wherein the cross-link interference is treated as noise. Under such cir-cumstance, analyses in [12]–[14] have shown that the Pareto optimal transmission strategy is transmitbeamforming. While beamforming is a structurally simple transmission strategy, finding the optimaltransmit beamformers for the MISO IFC is intrinsically difficult. More precisely, it has been proved [15]that finding the optimal beamformers that maximize system utilities, such as the weighted sum rate,the geometric mean rate, or the harmonic mean rate, is NP-hard in general. As a result, lots of effortshave focused on characterizing the optimal beamformer structures [12], [14], [16] in order to reducethe search dimension for finding the optimal beamforming vectors, or on investigating suboptimal butcomputationally efficient beamforming algorithms [15]–[17]. Another approach to studying these resourceconflicts encountered in the IFC is to use Game theory; see [11], [18], [19] for related works.The aforementioned beamforming designs all assume that the transmitters have the complete knowledgeof CSI. To provide the transmitters with complete CSI, the receivers need to periodically send the CSI
March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 3 (e.g., for frequency division duplexing systems) or training signals (e.g., for time division duplexingsystems) back to the transmitters. In contrast to the CSI, channel distribution information (CDI) canremain unchanged for a relatively long period of time and thus the amount of feedback informationcan be significantly reduced. With CDI at the transmitters, the ergodic rate region of the K -user MISOIFC has been analyzed and the structure of the Pareto optimal beamformers has been characterized in[20]. For a two-user case, an efficient algorithm for finding the Pareto boundary of the ergodic rateregion was presented in [21]. Unlike the ergodic achievable rate where the packet delay is not taken intoconsideration, the outage constrained achievable rate is more suitable for delay-sensitive applications,such as those involving voice or video data communications. For such outage constrained achievable rateregion, the authors of [22], [23] presented a numerical method for finding the Pareto boundary; however,the complexity of this algorithm increases exponentially with the number of transmitter-receiver pairs.Developing efficient beamforming design algorithms that can approach the outage constrained Paretoboundary is therefore important. While several efficient beamforming algorithms can be found in [24],[25], a different power-minimization design criterion was considered, instead of rate utility maximization.In this paper, we investigate efficient coordinated beamforming design algorithms for maximizing thesystem utility under rate outage constraints and individual power constraints. Specifically, we assumethat the MISO channel between each transmitter and receiver is composed of zero-mean circularlysymmetric complex Gaussian fading coefficients where the corresponding covariance matrix is known tothe transmitter. We formulate an outage constrained coordinated beamforming design problem, aimingat finding the Pareto optimal beamformers that maximize the system utility (e.g., the weighted sumrate) subject to a pre-assigned rate outage probability requirement and power constraints. However,due to the complicated nonconvex outage constraints, we propose a successive convex approximation(SCA) algorithm, where the original problem is successively approximated by a convex problem and thebeamforming solution is refined in an iterative manner. The convex approximation formulation is obtainedby applying the convex optimization based semidefinite relaxation (SDR) technique [26], followed bya logarithmic change of variables and first-order approximation techniques. We analytically show thatthe proposed SCA algorithm can yield a beamforming solution that is a stationary point for the SDRof the original problem. We further propose a round-robin-fashioned distributed SCA algorithm whereeach transmitter optimizes only its beamformer using local CDI with limited communication overhead ofmessage exchange with the other transmitters. It is shown by simulations that the two proposed algorithmsyield near-optimal performance with lower complexity compared with those reported in [22], [23].The remaining part of this paper is organized as follows. The system model and the outage constrained March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 4 coordinated beamforming problem are presented in Section II. In Section III, we present the proposedSCA algorithm and analyze its convergence property. In Section IV, the distributed SCA algorithm isdeveloped and analyzed. Simulation results that demonstrate the efficacy of the proposed algorithms arepresented in Section V. Finally, the conclusions are drawn in Section VI.Notation: The n -dimensional complex vectors and complex Hermitian matrices are denoted by C n and H n , respectively. The n × n identity matrix is denoted by I n . The superscripts ‘ T ’ and ‘ H ’ represent thematrix transpose and conjugate transpose, respectively. We denote k·k as the vector Euclidean norm. A (cid:23) and a (cid:23) respectively mean that matrix A is positive semidefinite (PSD) and vector a is elementwisenonnegative. The trace and rank of matrix A are denoted as Tr( A ) and rank( A ) , respectively. We use theexpression x ∼ CN ( µ , Q ) if x is circularly symmetric complex Gaussian distributed with mean µ andcovariance matrix Q . We denote exp( · ) (or simply e ( · ) ) as the exponential function, while ln( · ) and Pr {·} represent the natural log function and the probability function, respectively. For a variable a ik , where i, k ∈ { , . . . , K } , { a ik } k denotes the set { a i , . . . , a iK } , { a ik } k = i denotes the set { a ik } k excluding a ii ,and { a ik } is defined as the set containing all possible a ik , i.e., { a , . . . , a K , a , . . . , a KK } .II. S IGNAL M ODEL AND P ROBLEM S TATEMENT
We consider the K -user MISO IFC where each transmitter is equipped with N t antennas and eachreceiver with a single antenna. It is assumed that transmitters employ transmit beamforming to commu-nicate with their respective receivers. Let s i ( t ) denote the information signal sent from transmitter i , andlet w i ∈ C N t be the corresponding beamforming vector. The received signal at receiver i is given by x i ( t ) = h Hii w i s i ( t ) + K X k =1 ,k = i h Hki w k s k ( t ) + n i ( t ) , (1)where h ki ∈ C N t denotes the channel vector from transmitter k to receiver i , and n i ( t ) ∼ CN (0 , σ i ) isthe additive white Gaussian noise at receiver i where σ i > is the noise variance. As can be seen from(1), in addition to the noise, each receiver suffers from the cross-link interference P k = i h Hki w k s k ( t ) . Weassume that all receivers employ single-user detection where the cross-link interference is simply treatedas background noise. Under Gaussian signaling, i.e., s i ( t ) ∼ CN (0 , , the instantaneous achievable rateof the i th transmitter-receiver pair is known to be r i ( { h ki } k , { w k } ) = log (cid:12)(cid:12) h Hii w i (cid:12)(cid:12) P k = i (cid:12)(cid:12) h Hki w k (cid:12)(cid:12) + σ i ! . In this paper, we assume that the channel coefficients h ki are block-faded (i.e., quasi-static), andthat the transmitters have only the statistical information of the channels, i.e., the CDI. In particular, March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 5 it is assumed that h ki ∼ CN ( , Q ki ) for all k, i = 1 , . . . , K , where Q ki (cid:23) denotes the channelcovariance matrix and is known to all the transmitters. Since the transmission rate R i cannot be adaptedwithout CSI, the communication would be in outage whenever the transmission rate R i > is higherthan the instantaneous capacity that the channel can support. For a given outage probability requirement ( ǫ , . . . , ǫ K ) , the beamforming vectors { w i } must satisfy Pr { r i ( { h ki } k , { w i } ) < R i } ≤ ǫ i . Following[23], we define the corresponding ǫ i -outage achievable rate region as follows. Definition 1 [23]
Let P i > denote the power constraint of transmitter i , for i = 1 , . . . , K . The ratetuple ( R , . . . , R K ) is said to be achievable if Pr { r i ( { h ki } k , { w k } ) < R i } ≤ ǫ i , i = 1 , . . . , K , forsome ( w , . . . , w K ) ∈ W × · · · × W K where ǫ i ∈ (0 , is the maximum tolerable outage probability ofreceiver i and W i , { w ∈ C N t | k w k ≤ P i } . The ǫ i -outage achievable rate region is given by R = [ w i ∈W i ,i =1 ,...,K { ( R , . . . , R K ) | Pr { r i ( { h ki } k , { w k } ) < R i } ≤ ǫ i , i = 1 , . . . , K } . Given an outage requirement ( ǫ , . . . , ǫ K ) and an individual power constraint ( P , . . . , P K ) , our goalis to optimize { w k } such that the predefined system utility function U ( R , . . . , R K ) is maximized. Tothis end, we consider the following outage constrained coordinated beamforming design problem max w i ∈ C Nt ,R i ≥ ,i =1 ,...,K U ( R , . . . , R K ) (2a)s.t. Pr { r i ( { h ki } k , { w k } ) < R i } ≤ ǫ i , (2b) k w i k ≤ P i , i = 1 , . . . , K. (2c)Note that, as each user would prefer a higher transmission rate, a sensible system utility function shouldbe strictly increasing with respect to the individual rate R i for i = 1 , . . . , K , such that the optimal ( R , . . . , R K ) of problem (2) would lie on the so-called Pareto boundary of R [11]. In this paper, weconsider the following system utility function which captures a tradeoff between the system throughputand user fairness [27] U β ( { R i } ) = P Ki =1 α i R − βi − β , ≤ β < ∞ , β = 1 , P Ki =1 α i ln( R i ) , β = 1 , (3)where the coefficients α i ∈ [0 , for i = 1 , . . . , K with P Ki =1 α i = 1 represent the user priority, and theparameter β ∈ [0 , ∞ ) reflects the user fairness. For example, for β being , , and infinity, U β ( { R i } ) corresponds to the weighted sum rate, weighted geometric mean rate, weighted harmonic mean rate andthe weighted minimal rate, respectively. Hence, for β being , , and infinity, maximizing U β ( { R i } ) March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 6 is respectively equivalent to achieving the maximal throughput, proportional fairness, minimal potentialdelay, and the max-min fairness of users [28]. It can be verified that U β ( · ) is concave in { R i } for all β . However, since the outage constraints (2b) have a complicated structure as will be seen later, solvingproblem (2) is still challenging.One possible approach to solving such a nonconvex problem is via exhaustive search [22]. In [22],each of the cross-link interference is discretized into M levels, and given a set of cross-link interferencelevels, the maximum achievable rate for each receiver can be computed [22]. Since there are a total of K ( K − cross-user links, one has to exhaustively search over M K ( K − rate tuples. The complexityof this method thus increases exponentially with K ( K − , making this approach only viable when K is small. For a simple example of K = 3 and M = 10 , this method requires searching over ratetuples, which is computationally prohibitive in practice.III. P ROPOSED C ONVEX A PPROXIMATION M ETHOD
Our goal in this section is to develop an efficient approximation algorithm that obtains near-optimalsolutions of problem (2) for any number of transmitter-receiver pairs, K . To begin with, we note from[24, Appendix I] that the outage probability function in (2b) can actually be expressed in closed form as Pr { r i ( { h ki } k , { w k } ) < R i } = 1 − exp (cid:18) − (2 R i − σ i w Hi Q ii w i (cid:19) Y k = i w Hi Q ii w i w Hi Q ii w i + (2 R i − w Hk Q ki w k . (4)So, problem (2) can be rewritten as max w i ∈ C Nt ,R i ≥ ,i =1 ,...,K U ( R , . . . , R K ) (5a)s.t. ρ i exp (cid:18) (2 R i − σ i w Hi Q ii w i (cid:19) Y k = i (cid:18) R i − w Hk Q ki w k w Hi Q ii w i (cid:19) ≤ , (5b) k w i k ≤ P i , i = 1 , . . . , K, (5c)where ρ i , − ǫ i . Although the outage probability can now be expressed in closed form, problem (5) isstill difficult to solve, since the constraints in (5b) are still nonconvex and complicated. In the ensuingsubsections, we present a convex approximation method to handle problem (5) efficiently. A. Convex Approximation Formulation
The proposed convex approximation method starts with applying semidefinite relaxation (SDR), aconvex optimization based approximation technique [26]. Specifically, through SDR, we approximate thequadratic terms w Hk Q ki w k = Tr( w k w Hk Q ki ) in (5b) by the linear terms Tr( W k Q ki ) , where the rank-one March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 7 matrices w k w Hk are replaced by the PSD matrices W k of arbitrary rank( W k ) ≤ N t . The approximatedproblem is thus given by max W i ∈ H Nt ,R i ≥ ,i =1 ,...,K U ( R , . . . , R K ) (6a)s.t. ρ i exp (cid:18) (2 R i − σ i Tr ( W i Q ii ) (cid:19) Y k = i (cid:18) R i − W k Q ki )Tr ( W i Q ii ) (cid:19) ≤ , (6b) Tr ( W i ) ≤ P i , (6c) W i (cid:23) , i = 1 , . . . , K. (6d) We should mention that SDR has been widely used in various beamforming design problems (see [29]for a review), where, in most cases, a convex semidefinite program (SDP) approximation formulationcan be directly obtained via SDR and thus can be efficiently solved. Problem (6), however, is still notconvex yet due to the constraints in (6b). Therefore, further approximations are needed for problem (6).In contrast to SDR that essentially results in a larger problem feasible set, the second approximationis restrictive , in the sense that the obtained solution must also be feasible to problem (6). To illustratethis restrictive approximation, let us consider the following change of variables: e x ki , Tr( W k Q ki ) , e y i , R i − , z i , R i − W i Q ii ) = e y i − x ii , (7)for i, k = 1 , . . . , K , where x ki , y i , z i are slack variables. By substituting (7) into (6), one can reformulateproblem (6) as the following problem max { W i }∈S ,R i ≥ ,x ki ,y i ,z i ∈ R ,k,i =1 ,...,K U ( R , . . . , R K ) , (8a)s.t. ρ i e σ i z i Y k = i (cid:0) e − x ii + x ki + y i (cid:1) ≤ , (8b) Tr( W k Q ki ) ≤ e x ki , (8c) Tr( W i Q ii ) ≥ e x ii , (8d) R i ≤ log (1 + e y i ) , (8e) e y i − x ii ≤ z i , ∀ k ∈ K ci , i = 1 , . . . , K, (8f) where K ci , { , . . . , K }\{ i } , and the set S is defined in (9) below. Notice that we have replaced theequalities in (7) with inequalities as in (8c) to (8f). It can be verified by the monotonicity of the objectivefunction that all the inequalities in (8b) to (8f) would hold with equalities at the optimal points. We alsonote that, for example, if the optimal solution satisfies Tr( W i Q ii ) = 0 in (8d), then the optimal x ii has March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 8 to be minus infinity which is not attainable. Similar issues occur for
Tr( W k Q ki ) and x ki . In view ofthis, in (8) we have enforced W , . . . , W K to lie in the subset S , { W , . . . , W K (cid:23) | Tr( W i ) ≤ P i , Tr( W i Q ik ) ≥ δ ∀ i, k = 1 , . . . , K } , (9)where δ > . As long as δ is set to a small number, the rate loss due to (9) would be negligible.It is interesting to see that constraint (8b) is now convex; constraints (8d) and (8f) are also convex.Constraints (8c) and (8e) are not convex; nevertheless, they are relatively easy to handle compared withthe original (6b). Let ( ¯ w ¯ w H , . . . , ¯ w K ¯ w HK , ¯ R , . . . , ¯ R K ) be a feasible point of problem (8). Define ¯ x ki , ln( ¯ w Hk Q ki ¯ w k ) , ¯ y i , ln(2 ¯ R i − , ¯ z i , e ¯ y i − ¯ x ii , (10)for i, k = 1 , . . . , K . Then, { ¯ x ki } , { ¯ y i } , { ¯ z i } together with { ¯ R i } and ¯ W i , ¯ w i ¯ w Hi , i = 1 , . . . , K , arefeasible to problem (8). Here we conservatively approximate (8c) and (8e) at the point ( { ¯ x ki } k,i = k , { ¯ y i } ) .Since both of e x ki and log (1+ e y i ) are convex, their first-order lower bounds at ¯ x ki and ¯ y i are respectivelygiven by e ¯ x ki ( x ki − ¯ x ki + 1) and log (1 + e ¯ y i ) + e ¯ y i ( y i − ¯ y i )ln 2 · (1 + e ¯ y i ) . (11)Consequently, restrictive approximations for (8c) and (8e) are given by Tr( W k Q ki ) ≤ e ¯ x ki ( x ki − ¯ x ki + 1) , k ∈ K ci , (12a) R i ≤ log (1 + e ¯ y i ) + e ¯ y i ( y i − ¯ y i )ln 2 · (1 + e ¯ y i ) . (12b)By replacing (8c) and (8e) with (12a) and (12b), we obtain the following approximation for problem (5): max { W i }∈S ,R i ≥ ,x ki ,y i ,z i ∈ R ,k,i =1 ,...,K U ( R , . . . , R K ) , (13a)s.t. ρ i e σ i z i Y k = i (cid:0) e − x ii + x ki + y i (cid:1) ≤ , (13b) Tr( W k Q ki ) ≤ e ¯ x ki ( x ki − ¯ x ki + 1) , (13c) Tr( W i Q ii ) ≥ e x ii , (13d) R i ≤ log (1 + e ¯ y i ) + e ¯ y i ( y i − ¯ y i )ln 2 · (1 + e ¯ y i ) , (13e) e y i − x ii ≤ z i , ∀ k ∈ K ci , i = 1 , . . . , K. (13f) Problem (13) is a convex optimization problem; it can be efficiently solved by standard convex solverssuch as
CVX [30].
March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 9
Let ( ˆ W , . . . , ˆ W K ) and ( ˆ R , . . . , ˆ R K ) denote the optimal beamforming matrices and achievable ratesyielded by the approximation problem (13). Since the lower bounds in (11) may not be exactly tight, itmay hold, for ( ˆ W , . . . , ˆ W K ) and ( ˆ R , . . . , ˆ R K ) and for some i = 1 , . . . , K, that ρ i exp (2 ˆ R i − σ i Tr( ˆ W i Q ii ) ! Y k = i ˆ R i − W k Q ki )Tr( ˆ W i Q ii ) ! < , (14) i.e., the SINR outage probability is strictly less than ǫ i and thus the outage constraint is over satisfied.Alternatively, one can obtain a tight rate tuple ( ˜ R , . . . , ˜ R K ) , where ˜ R i ≥ ˆ R i for all i = 1 , . . . , K , bysolving the equations ρ i exp (2 R i − σ i Tr( ˆ W i Q ii ) ! Y k = i R i − W k Q ki )Tr( ˆ W i Q ii ) ! = 1 , (15) for i = 1 , . . . , K . Note that each equation in (15) can be efficiently solved by simple line search. Theobtained ( ˆ W , . . . , ˆ W K ) and ( ˜ R , . . . , ˜ R K ) then serve as an approximate solution for problem (6).In summary, the reformulation above consists of two approximation steps: a) the rank relaxation of w k w Hk to W k (cid:23) by SDR, and b) constraint restrictions of (8c) and (8e) by (13c) and (13e). Note that ifproblem (13) yields a rank-one optimal ( W , . . . , W K ) , a rank-one beamforming solution can be readilyobtained by rank-one decomposition of W i = w i w Hi for all i = 1 , . . . , K . It is then straightforward toverify by the restrictiveness of (13c) and (13e) that this rank-one beamforming solution ( w , . . . , w K ) is also feasible to the original problem (5) [i.e., problem (2)], thereby satisfying the desired rate outagerequirement. In view of this, it is important to investigate the conditions under which problem (13) canyield rank-one optimal ( W , . . . , W K ) . The following proposition provides one such condition: Proposition 1
Assume that (13) is feasible. Then there exists an optimal ( W , . . . , W K ) satisfying rank ( W i ) ≤ , i = 1 , . . . , K, if the number of users is no larger than three, i.e., K ≤ .Proof: Let ( { ˆ W i } , { ˆ R i } , { ˆ x ik } , { ˆ y i } , { ˆ z i } ) denote an optimal solution of problem (13). Consider max W i (cid:23) Tr( W i Q ii ) (16a)s.t. δ ≤ Tr( W i Q ik ) ≤ e ¯ x ik (ˆ x ik − ¯ x ik + 1) , k ∈ K ci (16b) Tr( W i Q ii ) ≥ δ, Tr( W i ) ≤ P i , (16c)for all i = 1 , . . . , K . By (9) and (13c), ˆ W i is also feasible to the above problem (16). Moreover, by(13b), (13d), (13e), (13f) and by the monotonicity of U ( R , . . . , R K ) , one can show, by contradiction,that ˆ W i is actually optimal to problem (16), for all i = 1 , . . . , K . Let W ′ i be an optimal solution to(16), for i = 1 , . . . , K . Then, one can also verify that ( W ′ , . . . , W ′ K ) is optimal to problem (13). We March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 10 hence focus on problem (16). Firstly, since problem (16) is assumed to be feasible and the objective is tomaximize
Tr( W i Q ii ) , the constraint Tr( W i Q ii ) ≥ δ in (16c) is actually irrelevant and can be droppedwithout affecting the optimal solution. Secondly, it is easy to observe that, for each k ∈ K ci , it is either δ < Tr( W i Q ik ) = e ¯ x ik (ˆ x ik − ¯ x ik + 1) or Tr( W i Q ik ) = e ¯ x ik (ˆ x ik − ¯ x ik + 1) = δ at the optimum; thatis, for either case, it is equivalent to having one equality constraint in (16b) at the optimum for each k ∈ K ci . As a result, problem (16) equivalently has only K constraints. According to [31, Theorem 3.2],there always exists an optimal solution W i of problem (16) satisfying rank( W i ) ≤ √ K for i = 1 , . . . , K. (17)Therefore, if K ≤ , there always exists a rank-one optimal ( W , . . . , W K ) for problem (13). (cid:4) We should mention that K ≤ is only a sufficient condition but not a necessary condition. For K > ,there may exist other conditions under which a rank-one optimal solution exists for problem (13). If theoptimal ( W , . . . , W K ) is not of rank one, then one may resort to the rank-one approximation proceduressuch as Gaussian randomization [26], [29] to obtain an approximate solution to (2). Note that, in that case,the utility achieved by the randomized solution would be no larger than U ( ˜ R , . . . , ˜ R K ) . Surprisingly,in our computer simulations, we found that problem (13) is always solved with rank-one optimal { W i } .Some insightful analyses, which explain why problem (16) is often solved with rank-one optimal W i forrandomly generated problem instances, can be found in [17]. B. Successive Convex Approximation (SCA)
Formulation (13) is obtained by approximating problem (8) at the given feasible point ( { ¯ w i ¯ w Hi } , { ¯ R i } ) ,as described in (10). This approximation can be further improved by successively approximating problem(8) based on the optimal solution ( { W i } , { R i } ) obtained by solving (13) in the previous approximation.Specifically, let ( ˆ W [ n − , . . . , ˆ W K [ n − be the optimal beamforming matrices obtained in the ( n − thiteration, and, similar to (15), let ( ˜ R i [ n − , . . . , ˜ R i [ n − be the corresponding achievable rate tupleobtained by solving the following K equations ρ i exp (2 R i − σ i Tr( ˆ W i [ n − Q ii ) ! Y k = i R i − W k [ n − Q ki )Tr( ˆ W i [ n − Q ii ) ! = 1 , i = 1 , . . . , K. (18) Moreover, let ¯ x ki [ n −
1] = ln(Tr( ˆ W k [ n − Q ki )) , (19a) ¯ y i [ n −
1] = ln(2 ˜ R i [ n − − , i, k = 1 , . . . , K. (19b) March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 11
By replacing ¯ x ki and ¯ y i in (13) with ¯ x ki [ n − and ¯ y i [ n − in (19) for k ∈ K ci , i = 1 , . . . , K , we solve,in the n th iteration, the following convex optimization problem ( { ˆ W i [ n ] } , { ˆ R i [ n ] } , { ˆ x ik [ n ] } , { ˆ y i [ n ] } , { ˆ z i [ n ] } ) =arg max { W i }∈S ,R i ≥ x ik ,y i ,z i ∈ R i,k =1 ,...,K U ( R , . . . , R K ) (20a)s.t. ρ i e σ i z i Y k = i (cid:0) e − x ii + x ki + y i (cid:1) ≤ , (20b) Tr( W k Q ki ) ≤ e ¯ x ki [ n − ( x ki − ¯ x ki [ n −
1] + 1) , (20c) Tr( W i Q ii ) ≥ e x ii , (20d) R i ≤ log (1 + e ¯ y i [ n − ) + e ¯ y i [ n − ( y i − ¯ y i [ n − · (1 + e ¯ y i [ n − ) , (20e) e y i − x ii ≤ z i , ∀ k ∈ K ci , i = 1 , . . . , K. (20f) Note that the rate ˜ R i [ n − obtained by (18) is no less than ˆ R i [ n − for all i = 1 , . . . , K , and thusthe former is used to compute { ¯ y i [ n − } as the point for successive approximation. In fact, successiveapproximation ensures monotonic improvement of the utility U ( ˜ R [ n ] , . . . , ˜ R K [ n ]) . Let us define ¯ z i [ n − , e ¯ y i [ n − − ¯ x ii [ n − , i = 1 , . . . , K. (21)Then, by (18), (19) and (21), one can show that ( { ˆ W i [ n − } , { ˜ R i [ n − } , { ¯ x ik [ n − } , { ¯ y i [ n − } , { ¯ z i [ n − } ) is a feasible point of (20). As a result, we have U ( ˜ R [ n ] , . . . , ˜ R K [ n ]) ≥ U ( ˆ R [ n ] , . . . , ˆ R K [ n ]) ≥ U ( ˜ R [ n − , . . . , ˜ R K [ n − , ∀ n ≥ . (22)The proposed successive convex approximation (SCA) algorithm is summarized in Algorithm 1. C. Convergence Analysis
Convergence properties of Algorithm 1 is given below.
Theorem 1
Suppose that the utility U ( R , . . . , R K ) is differentiable and strictly increasing with respectto R i , for i = 1 , . . . , K . Then, the sequence, { U ( ˜ R [ n ] , . . . , ˜ R K [ n ]) } ∞ n =1 generated by Algorithm 1,converges, and any limit point of the sequence { ( ˆ W [ n ] , . . . , ˆ W K [ n ]) , ( ˜ R [ n ] , . . . , ˜ R K [ n ]) } ∞ n =1 is astationary point of problem (8) as well as a stationary point of problem (6) with extra constraints Tr( W i Q ik ) ≥ δ for i, k = 1 , . . . , K (see (9) ) .Proof of Theorem 1: As discussed earlier, the utility U ( ˜ R [ n ] , . . . , ˜ R K [ n ]) is nondecreasing with n .Moreover, due to the individual power constraints, the sequence { U ( ˜ R [ n ] , . . . , ˜ R K [ n ]) } ∞ n =1 is bounded,which implies the convergence of U ( ˜ R [ n ] , . . . , ˜ R K [ n ]) . March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 12
Algorithm 1
SCA algorithm for solving problem (2) Given ( ¯ w ¯ w H , . . . , ¯ w K ¯ w HK ) ∈ S and ( ¯ R , . . . , ¯ R K ) that are feasible to (6). Set ˆ W i [0] = ¯ w i ¯ w Hi and ˜ R i [0] = ¯ R i for all i = 1 , . . . , K, and set n = 0 . repeat n := n + 1 . Obtain { ¯ x ki [ n − } and { ¯ y i [ n − } by (19), and solve problem (20) to obtain the optimal solution ˆ u [ n ] , ( { ˆ W i [ n ] } , { ˆ R i [ n ] } , { ˆ x ik [ n ] } , { ˆ y i [ n ] } , { ˆ z i [ n ] } ) . Compute ( ˜ R [ n ] , . . . , ˜ R K [ n ]) by solving (18). until the stopping criterion is met. Obtain w ⋆i by decomposition of ˆ W i [ n ] = w ⋆i ( w ⋆i ) H for all i , if ˆ W i [ n ] are all of rank one; otherwiseperform Gaussian randomization [29] to obtain a rank-one approximate solution of (2).Let ˆ u [ n ] , ( { ˆ W i [ n ] } , { ˆ R i [ n ] } , { ˆ x ik [ n ] } , { ˆ y i [ n ] } , { ˆ z i [ n ] } ) , denote the optimal solution of (20). Toprove that any limit point of ˆ u [ n ] is a stationary point of (8), two key observations are needed. Firstly,we note that the proposed SCA algorithm is in fact an inner approximation algorithm in the nonconvexoptimization literature [32]. In particular, the nonconvex constraints in (8c) and (8e), i.e., Ψ ki ( W k , x ki ) , Tr( W k Q ki ) − e x ki ≤ , k ∈ K ci , Φ i ( R i , y i ) , R i − log (1 + e y i ) ≤ , i = 1 , . . . , K, are respectively replaced by ¯Ψ ki ( W k , x ki | ¯ x ki [ n − , Tr( W k Q ki ) − e ¯ x ki [ n − ( x ki − ¯ x ki [ n −
1] + 1) ≤ , (23) ¯Φ i ( R i , y i | ¯ y i [ n − , R i − log (1 + e ¯ y i [ n − ) + e ¯ y i [ n − ( y i − ¯ y i [ n − · (1 + e ¯ y i [ n − ) ≤ , (24)for all k ∈ K ci , i = 1 , . . . , K . One can verify that ¯Ψ ki ( W k , x ki | ¯ x ki [ n − and ¯Φ i ( R i , y i | ¯ y i [ n − satisfy Ψ ki ( ˆ W k [ n − , ¯ x ki [ n − ki ( ˆ W k [ n − , ¯ x ki [ n − | ¯ x ki [ n − (25) ∂ Ψ ki ( W k , x ki ) ∂ W k = ∂ ¯Ψ ki ( W k , x ki | ¯ x ki [ n − ∂ W k (26) ∂ Ψ ki ( W k , x ki ) ∂x ki (cid:12)(cid:12)(cid:12)(cid:12) x ki =¯ x ki [ n − = ∂ ¯Ψ ki ( W k , x ki | ¯ x ki [ n − ∂x ki (cid:12)(cid:12)(cid:12)(cid:12) x ki =¯ x ki [ n − (27) Φ i ( ˆ R i [ n − , ¯ y i [ n − i ( ˆ R i [ n − , ¯ y i [ n − | ¯ y i [ n − (28) ∂ Φ i ( R i , y i ) ∂R i = ∂ ¯Φ i ( R i , y i | ¯ y i [ n − ∂R i (29) March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 13 ∂ Φ i ( R i , y i ) ∂y i (cid:12)(cid:12)(cid:12)(cid:12) y i =¯ y i [ n − = ∂ ¯Φ i ( R i , y i | ¯ y i [ n − ∂y i (cid:12)(cid:12)(cid:12)(cid:12) y i =¯ y i [ n − , (30)for all k ∈ K ci and i = 1 , . . . , K .Secondly, the restrictive approximations made in (20c) and (20e) are asymptotically tight as n → ∞ : Claim 1
It holds true that lim n →∞ | ˆ x ik [ n ] − ¯ x ik [ n − | = 0 , (31) lim n →∞ | ˆ y i [ n ] − ¯ y i [ n − | = 0 , lim n →∞ ( ˜ R i [ n ] − ˆ R i [ n ]) = 0 , (32) for all k ∈ K ci , i = 1 , . . . , K . Claim 1 is proved in Appendix A. Moreover, by the monotonicity of U ( R , . . . , R K ) and due to (9), itis not difficult to verify that: Claim 2
The sequence { ˆ u [ n ] } ∞ n =0 generated by Algorithm 1 is bounded. Now let us consider the KKT conditions of (20). Denote L ( ˆ u [ n ] , λ [ n ] |{{ ¯ x ki [ n − } k = i } i , { ¯ y i [ n − } ) asthe Lagrangian of (20). For ease of explanation, let Θ i ( { x ki } k , y i , z i ) , ρ i e σ i z i Q k = i (1+ e − x ii + x ki + y i ) − denote the constraint function in (20b), and consider the following Lagrangian-stationarity condition: ∂ L ( ˆ u [ n ] , λ [ n ] |{{ ¯ x ki [ n − } k = i } i , { ¯ y i [ n − } ) ∂x ki = λ b i [ n ] ∂ Θ i ( { ˆ x ji [ n ] } j , ˆ y i [ n ] , ˆ z i [ n ]) ∂x ki + λ c ki [ n ] ∂ ¯Ψ ki ( ˆ W k [ n ] , ˆ x ki [ n ] | ¯ x ki [ n − ∂x ki = 0 , ∀ k = i, (33) where λ [ n ] , ( { λ b i [ n ] } , {{ λ c ik [ n ] } k = i } i , { λ d i [ n ] } , { λ e i [ n ] } , { λ f i [ n ] } , { λ Pi [ n ] } , { λ δik [ n ] } ) are dual variablesassociated with constraints (20b)-(20f), the transmit power constraint and Tr( W i Q ik ) ≥ δ . Since problem(20) satisfies the Slater’s condition, the dual variables are bounded [33]. Moreover, ˆ u [ n ] is bounded as wellby Claim 2. Therefore, there exists a subsequence { n , . . . , n ℓ , . . . } ⊆ { , . . . , n, . . . } and a primal-duallimit point, denoted by ˆ u ⋆ , ( { ˆ W ⋆i } , { ˆ R ⋆i } , { ˆ x ⋆ik } , { ˆ y ⋆i } , { ˆ z ⋆i } ) and λ ⋆ , ( { λ b ⋆i } , {{ λ c ⋆ik } k = i } i , { λ d ⋆i } , { λ e ⋆i } , { λ f ⋆i } , { λ P ⋆i } , { λ δ⋆ik } ) such that lim ℓ →∞ ˆ u [ n ℓ ] = ˆ u ⋆ , lim ℓ →∞ λ [ n ℓ ] = λ ⋆ , (34) where ( ˆ u ⋆ , λ ⋆ ) is primal-dual feasible to (20). Consider (33) over the subsequence { n , . . . , n ℓ , . . . } .By taking ℓ → ∞ in (33), and by (27), (31) and (34), we obtain λ b ⋆i ∂ Θ i ( { ˆ x ⋆ji } j , ˆ y ⋆i , ˆ z ⋆i ) ∂x ki + λ c ⋆ki ∂ Ψ ki ( ˆ W ⋆k , ˆ x ⋆ki ) ∂x ki = 0 , which is the Lagrangian-stationarity condition of (8) corresponding to x ki . By applying similar argumentsabove to all the other KKT conditions of (20) and by Claims 1 and 2, we end up with the conclusionthat ˆ u ⋆ satisfies the KKT conditions of problem (8) and thus is a stationary point. March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 14
What remains is to show that any stationary point of (8) is also a stationary point of (6) if the constraintset (9) is added in (6). This can be proved by carefully examining the equivalence of the KKT conditionsof the two problems. Due to the space limitation, the detailed derivations are omitted here. (cid:4)
As the SCA algorithm only guarantees to provide a stationary point, the approximation accuracydepends on the initial point ( { ˆ W i [0] } , { ˜ R i [0] } ) . A possible choice is to initialize Algorithm 1 via someheuristic transmission strategies. For example, one can obtain an initial point ( { ¯ w i } , { ¯ R i } ) of problem (5)through the simple maximum-ratio transmission (MRT) strategy. In this strategy, the beamforming vectors { ¯ w i } are simply set to ¯ w i = √ P i q i where q i ∈ C N t is the principal eigenvector of Q ii for i = 1 , . . . , K ,with k q i k = 1 . The corresponding rate ˜ R i can be obtained by solving (15) with { ˆ W i } = { ¯ w i ¯ w Hi } .Analogously, one can also obtain an initial point by the zero-forcing (ZF) transmission strategy, providedthat the column space of Q ii is not subsumed by the column space of P Kk = i Q ik , for all i = 1 , . . . , K .IV. D ISTRIBUTED I MPLEMENTATION
For Algorithm 1, we have implicitly assumed that there exists a control center in the network, collectingall the CDI of users and computing the beamforming solution in a centralized manner. In this section,we propose a distributed version for Algorithm 1, where each transmitter i only needs to optimize itsown beamformer, using only its local CDI, i.e., { Q ik } k , and some information obtained from the othertransmitters. Since each of the subproblems involved has a much smaller problem size than the originalproblem (8), even for a centralized implementation, the proposed distributed optimization algorithm canbe used to reduce the computation overhead of the control center.The idea of the proposed distributed algorithm is to solve problem (8) from one transmitter to another, ina round-robin fashion (i.e., the Gauss-Seidel fashion). Suppose that transmitter 1 optimizes its beamformerfirst, followed by transmitter 2 and so on, and let n denote the index of the current round. Then, in the n th round, transmitter i solves the following problem ( ˆ w i [ n ] , ˆ R [ n, i ] , . . . , ˆ R K [ n, i ]) = arg max k w i k ≤ P i R ,...,R K ≥ U ( R , . . . , R K ) (35a)s.t. ρ i exp (cid:18) (2 R i − σ i w Hi Q ii w i (cid:19) Y k = i (cid:0) R i − (cid:1) e ¯ x ki [ n − u ki ] w Hi Q ii w i ! ≤ , (35b) ρ j exp (2 R j − σ j e ¯ x jj [ n − u ji ] ! (cid:0) R j − (cid:1) w Hi Q ij w i e ¯ x jj [ n − u ji ] ! × Y k = jk = i (cid:0) R j − (cid:1) e ¯ x kj [ n − u ki ] e ¯ x jj [ n − u ji ] ! ≤ , j ∈ K ci , (35c) March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 15 where ¯ x kj [ n − u ki ] = ln( ˆ w Hk [ n − u ki ] Q kj ˆ w k [ n − u ki ]) , and u ki is equal to one if k > i and zero otherwise.Note that for (35), only ˆ w i [ n ] is optimized while { ¯ x kj [ n − u kj ] } k = i,j are given. Once the beamformingsolution of (35) is obtained, { ¯ x ik [ n ] } k are updated according to the optimal ˆ w i [ n ] and then passed to allthe other transmitters for their subsequent beamforming optimization . There are two interesting pointsto note here. Firstly, as can be seen from (35b) and (35c), transmitter i not only optimizes its rate R i , butalso takes into account the rate outage constraints for all the other users. The constraints in (35c) indicatethat transmitter i needs to regulate its own transmission in order not to violate the outage requirementof the other users. Secondly, to solve (35), transmitter i only needs the local CDI, i.e., { Q ik } k .Similar difficulties arise here as in problem (5) since problem (35) is not convex. We hence applythe same approximation techniques in Section III-A to approximate (35). In particular, we first applySDR, followed by the reformulation as described by (7), and the first-order approximations in (12). Theresulting convex optimization problem reads ( ˆ W i [ n ] , { ˆ R k [ n, i ] } , { ˆ x ik [ n ] } k , { ˆ y k [ n, i ] } , { ˆ z k [ n, i ] } ) = arg max W i ∈S i ,R k ,x ik ,y k ,z k k =1 ,...,K U ( R , . . . , R K ) (36a)s.t. ρ i e σ i z i Y k = i (cid:16) e − x ii +¯ x ki [ n − u ki ]+ y i (cid:17) ≤ , (36b) ρ j e σ j z j (cid:16) e − ¯ x jj [ n − u ji ]+ x ij + y j (cid:17) Y k = j,k = i (cid:16) e − ¯ x jj [ n − u ji ]+¯ x kj [ n − u ki ]+ y j (cid:17) ≤ , j ∈ K ci , (36c) Tr ( W i Q ii ) ≥ e x ii , (36d) Tr ( W i Q ik ) ≤ e ¯ x ik [ n − ( x ik − ¯ x ik [ n −
1] + 1) , k ∈ K ci , (36e) R j ≤ (cid:18) ln(1 + e ¯ y j [ n,i − ) + e ¯ y j [ n,i − e ¯ y j [ n,i − ( y j − ¯ y j [ n, i − (cid:19) , j = 1 , . . . , K, (36f) e y i − x ii ≤ z i , e y j − ¯ x jj [ n − u ji ] ≤ z j , j ∈ K ci , (36g)where S i , { W i (cid:23) | Tr( W i ) ≤ P i , Tr( W i Q ik ) ≥ δ, k = 1 , . . . , K } , ¯ x kj [ n − u ki ] = ln (cid:16) Tr( ˆ W k [ n − u ki ] Q kj ) (cid:17) , (37) ¯ y j [ n, i −
1] = ln (cid:16) ˜ R j [ n,i − − (cid:17) , (38)for j, k = 1 , . . . , K , and, similar to (18), ˜ R j [ n, i ] ≥ ˆ R j [ n, i ] is obtained by solving the following equations ρ j exp (2 R j − σ j e ¯ x jj [ n − u ji ] ! Y k = j (cid:16) R j − e ¯ x kj [ n − u ki ] − ¯ x jj [ n − u ji ] (cid:17) = 1 , (39) In this paper, we assume that the communication between transmitters for message exchange is error-free.
March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 16 for j = 1 , . . . , K. It is worth mentioning that problem (36) is only solved once and successive approxima-tion is not performed as in Algorithm 1. As long as problem (36) is solved by transmitter i , the algorithmdirectly goes to the next optimization problem performed by transmitter i + 1 . Successive approximationis now performed implicitly from one transmitter to another in a round-robin fashion. We summarize theproposed distributed SCA algorithm in Algorithm 2. Algorithm 2
Distributed SCA algorithm for solving problem (2) Given an initial beamforming matrix ˆ W i [0] at transmitter i , for i = 1 , . . . , K . For all i = 1 , . . . , K , transmitter i computes { ¯ x ik [0] } k by (37), and pass them to the other transmitters. Set n = 0 repeat n = n + 1 for i = 1 , . . . , K do User i solves (39) to obtain { ˜ R j [ n, i − } j and compute { ¯ y j [ n, i − } j , followed by solving(36) to obtain the solution ( ˆ W i [ n ] , { ˆ R k [ n, i ] } , { ˆ x ik [ n ] } k , { ˆ y k [ n, i ] } , { ˆ z k [ n, i ] } ) . User i computes { ¯ x ik [ n ] } k by (37) and passes them to all the other transmitters. end for until the predefined stopping criterion is met. For i = 1 , . . . , K , each transmitter i decomposes ˆ W i [ n ] = ˆ w i ˆ w Hi , if ˆ W i [ n ] is of rank one; otherwiseperform Gaussian randomization to obtain a rank-one approximate solution.Analogous to Algorithm 1, we can show that Algorithm 2 generates a stationary point of problem (8)as stated in the following theorem. Theorem 2
Suppose that U ( R , . . . , R K ) is differentiable and is strictly increasing with respect to R i , for i = 1 , . . . , K . Then, the sequence { U ( ˜ R [ n, i ] , . . . , ˜ R K [ n, i ]) } ∞ n =1 generated by Algorithm 2converges to a common value for all i = 1 , . . . , K . Moreover, for all i , any limit point of the sequence { ( ˆ W [ n ] , . . . , ˆ W k [ n ] , ˜ R [ n, i ] , . . . , ˜ R [ n, i ]) } ∞ n =1 is a stationary point of problem (8) as well as a sta-tionary point of problem (6) (with the extra constraints in (9) ) . Different from the proof for Theorem 1, the proof for Theorem 2 is more involved, since the beam-forming vectors of transmitters are not simultaneously optimized as in Algorithm 1 but are individuallyoptimized in a round-robin manner. The detailed proof of Theorem 2 is presented in Appendix B.
March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 17
Remark 1
An important issue concerning distributed optimization algorithms is the communicationoverhead introduced by message exchange between transmitters. To address this, we compare the com-munication overhead of the proposed Algorithm 2 with the following two schemes. Scheme 1: All thetransmitters directly exchange their CDI so that the design problem (2) can be handled independentlyby each transmitter. Scheme 2: A control center gathers the CDI from all transmitters, optimizes thebeamforming vectors, and distributes the beamforming solutions to the transmitters. We consider a cellularsystem where all the transmitters (i.e., BSs) are connected by dedicated backhaul links (e.g., optical fibers)and the BSs exchange messages in a point-to-point fashion. Since, in Algorithm 2, transmitter i needs toinform { ¯ x ik [ n ] } k ( K real values) to all the other K − transmitters in each round, the communicationoverhead due to transmitter i is quantified by the amount of K ( K − real values. Hence, the totalcommunication overhead of Algorithm 2 is N × K × K ( K −
1) = K ( K − N real values, where N is the number of rounds run by Algorithm 2. For scheme 1, each transmitter needs to send K covariance matrices (which contain KN t real values) to all the other K − transmitters. Therefore,the associated total communication overhead is given by K × ( K − × KN t = K ( K − N t realvalues. Therefore, for scheme 1, if N < N t , then the proposed Algorithm 2 has a smaller amount ofcommunication overhead. For scheme 2, there are K covariance matrices sent from the transmitters tothe control center, and the optimal solution { w ⋆i , R ⋆i } passed from the control center to transmitter i for i = 1 , . . . , K , respectively. Hence, the communication overhead is K N t + K (2 N t + 1) real values.Therefore, for scheme 2, the proposed Algorithm 2 has a smaller amount of communication overhead if N < N t / ( K −
1) + (2 N t + 1) / ( K − K ) . As we show in the simulation section, Algorithm 2 in generalcan converge in less than 15 rounds for a system with K ≤ and N t = 8 .We should mention here that, while in general the proposed distributed algorithm is more efficient interms of computation and communication overhead, it may result in larger transmission delays (due tothe iterative optimizations between transmitters) compared with the centralized schemes. Remark 2
We should emphasize that the proposed beamforming design is based on the users’ statisticalchannel information, which usually changes much more slowly compared to the instantaneous CSI, sobeamforming optimization need not be performed frequently. As a result, the throughput loss induced bythe round-robin optimization in Algorithm 2 should not be a serious concern.V. S
IMULATION R ESULTS
In the section, we demonstrate the performance of the proposed Algorithm 1 and Algorithm 2 forsolving the outage constrained coordinated beamforming problem in (2). In the simulations, we consider
March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 18 η A v e r a g e W e i g h t e d S u m R a t e ( b i t s / s e c / H z ) Exhaustive SearchSCA Algorithm1/ σ =20dB1/ σ =0dB1/ σ =10dB (a) η A v e r a g e W e i g h t e d H a r m o n i c M e a n R a t e ( b i t s / s e c / H z ) Exhaustive SearchSCA Algorithm1/ σ =10dB1/ σ =20dB 1/ σ =0dB (b) Fig. 1: Simulation results of the proposed SCA algorithm (Algorithm 1), for K = 2 , N t = 4 , and ( α , α ) = ( , ) ; (a) weighted sum rate versus η , (b) weighted harmonic mean rate versus η . Each ofthe results is obtained by averaging over 500 realizations of { Q ki } . β = 0 , β = 1 , and β = 2 for the objective function U β ( R , . . . , R K ) , corresponding to maximizationof the weighted sum rate, the weighted geometric mean rate, and the weighted harmonic mean rate,respectively. All receivers are assumed to have the same noise power, i.e., σ = · · · = σ K , σ ,and all power constraints are set to one, i.e., P = · · · = P K = 1 . The parameter δ in (9) is setto − . The channel covariance matrices Q ki are randomly generated. We normalize the maximumeigenvalue of Q ii , i.e., λ max ( Q ii ) , to one for all i , and normalize λ max ( Q ki ) to a value η ∈ (0 , forall k ∈ K ci , i = 1 , . . . , K . The parameter η , thereby, represents the relative cross-link interference level.If not mentioned specifically, all Q ki are of full rank, and the outage probability requirements are set tothe same value, i.e., ǫ = · · · = ǫ K = 0 . , indicating a outage probability. The stopping criterion ofAlgorithm 1 is | U ( ˜ R [ n ] , . . . , ˜ R K [ n ]) − U ( ˜ R [ n − , . . . , ˜ R K [ n − | U ( ˜ R [ n − , . . . , ˜ R K [ n − < . . That is, Algorithm 1 stops if the improvement in system utility is less than of the system utilityachieved in the previous iteration. The simple MRT solution is used to initialize both Algorithm 1 andAlgorithm 2. The convex solver CVX [30] is used to solve the convex problems (20) and (36).
Example 1:
We first examine the approximation performance of the proposed SCA algorithm, bycomparing it with the exhaustive search method in [22]. In view of the tremendous complexity overheadsof this exhaustive search method, we consider a simple case where only two transmitter-receiver pairsare present, i.e. K = 2 , and set N t = 4 . Figure 1(a) shows the simulation results for the comparisonof the achievable weighted sum rate between the proposed SCA algorithm and the exhaustive search March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 19 R (bits/sec/Hz) R ( b i t s / s e c / H z ) Pareto BoundaryWeighted Sum RateWeighted GeometricMean RateWeighted HarmonicMean RateOptimal Rate Tuples (a) R (bits/sec/Hz) R ( b i t s / s e c / H z ) Pareto BoundaryWeighted Sum RateWeighted GeometricMean RateWeighted HarmonicMean RateOptimal Rate Tuples (b)
Fig. 2: Converge trajectories of the proposed SCA algorithm. K = 2 , N t = 4 , η = 0 . ; (a) ( α , α ) =( , ) , (b) ( α , α ) = ( , ) . The results are obtained using a typical set of randomly generated { Q ki } .method against the cross-link interference level η , where the weights are given by ( α , α ) = ( , ) .Each simulation curve is obtained by averaging over 500 realizations of randomly generated { Q ki } . Fromthis figure, we can observe that, for /σ = 0 dB and /σ = 10 dB, the proposed SCA algorithm canattain almost the same average sum rate performance as the exhaustive search method, indicating that theproposed SCA algorithm yields near-optimal solutions for the outage constrained beamforming designproblem (2). For /σ = 20 dB, it can be observed that there is a small gap between the rate achieved bythe proposed SCA algorithm and that by the exhaustive search method. Nonetheless, this gap is relativelysmall and is within of the sum rate achieved by the exhaustive search method. Figure 1(b) displaysthe simulation results under the same setting as in Figure 1(a) except that the objective function is nowthe average harmonic mean rate. As the mean rate performance of SCA algorithm is almost the same asthat of the exhaustive search method, its solution is nearly optimal for problem (2).To examine how the proposed SCA algorithm converges, we illustrate in Figure 2(a) the trajectoriesof the optimal rate tuple of problem (20) in each iteration of Algorithm 1, where the weighted sum rate,the geometric mean rate, and the harmonic mean rate are all considered. The user priority weights are setto ( α , α ) = ( , ) , and the Pareto boundary is obtained by the exhaustive search method in [22]. Onecan see from this figure that, for all rate utility functions, the proposed SCA algorithm first approachesthe Pareto boundary and then converges to the corresponding optimal rate tuple along the boundary. InFigure 2(b), we display similar results with an asymmetric user priority, i.e., ( α , α ) = ( , ) . It can beobserved that the SCA algorithm still converges to the optimal rate tuples in a similar fashion. March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 20 σ (dB) A v e r age W e i gh t ed S u m R a t e ( b i t s / s e c / H z ) SCA Algorithm, η =0.2SCA Algorithm, η =1MRT, η =0.2MRT, η =1 (a) σ (dB) A v e r age W e i gh t ed S u m R a t e ( b i t s / s e c / H z ) SCA Algorithm, η =0.2SCA Algorithm, η =1MRT, η =0.2MRT, η =1ZF (b) Fig. 3: Simulation results of average achievable sum rate versus /σ ; (a) K = N t = 4 , and fullrank { Q ki } , (b) K = 4 , N t = 8 and rank( Q ki ) = 2 for all k, i . The priority weights are set to ( α , α , α , α ) = ( , , , ) . The results are obtained by averaging over 500 realizations of { Q ki } . Example 2:
To further demonstrate the effectiveness of the proposed SCA algorithm, we evaluate theperformance of the SCA algorithm for the case of K = N t = 4 in this example. (Since under this setting,the exhaustive search method in [22] is too complex to implement, and, to the best of our knowledge,there is no existing method for comparison, we can only compare the proposed SCA algorithm withthe heuristic MRT and ZF schemes.) Figure 3(a) shows the simulation results of the average achievablesum rate versus /σ . From this figure, one can observe that the proposed SCA algorithm yields bettersum rate performance than the MRT scheme, especially when /σ > dB. For /σ ≤ dB, thetwo methods exhibit comparable performance. In Figure 3(b), we have shown the simulation results for K = 4 , N t = 8 and rank( Q ki ) = 2 for all k, i . Under this setting, the ZF scheme is feasible and itsaverage sum rate performance is also shown in Figure 3(b). It can be observed from this figure that theZF scheme outperforms the MRT scheme for high /σ or when the cross-link interference is strong( η = 1 ). Nevertheless, as can be seen from Figure 3(b), the proposed SCA algorithm still outperformsboth the MRT and the ZF schemes.Figure 4 demonstrates the simulation results for the weighted geometric mean rate and the weightedharmonic mean rate, for K = N t = 4 and for an asymmetric weighting ( α , α , α , α ) = ( , , , ) .Performance comparison results similar to those in Figure 3 can also be observed in this figure. Inaddition, it is interesting to note from Figure 4 that, in contrast to the sum rate performance as shown inFigure 3, the weighted geometric mean rates and weighted harmonic mean rates achieved by the proposed March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 21 σ (dB) A v e r age W e i gh t ed G eo m e t r i c M ean R a t e ( b i t s / s e c / H z ) SCA Algorithm, η =0.2SCA Algorithm, η =1MRT, η =0.2MRT, η =1 (a) σ (dB) A v e r age W e i gh t ed H a r m on i c M ean R a t e ( b i t s / s e c / H z ) SCA Algorithm, η =0.2SCA Algorithm, η =1MRT, η =0.2MRT, η =1 (b) Fig. 4: Simulation results of the proposed SCA algorithm (Algorithm 1), for K = N t = 4 and ( α , α , α , α ) = ( , , , ) ; (a) weighted geometric mean rate versus /σ , (b) weighted harmonicmean rate versus /σ . Each of the results is obtained by averaging over 500 realizations of { Q ki } .SCA algorithm in Figure 4(a) and Figure 4(b) saturate for high /σ . These phenomena might resultfrom the fact that user fairness plays a more prominent role in the geometric mean rate and the harmonicmean rate; and thereby in the interference dominated region (i.e., when /σ or η is large), the geometricmean rate and the harmonic mean rate cannot increase as fast as the weighted sum rate. Example 3:
In this example, we examine the performance of the proposed distributed SCA algorithm(Algorithm 2). Figure 5(a) shows the convergence behaviors (the evolution of sum rate at each round)of the distributed SCA algorithm for N t = 8 , K = 4 , , and for N t = 12 , K = 6 , where /σ = 10 dB, η = 0 . . Each curve in Figure 5(a) is obtained by averaging over 500 sets of randomly generated { Q ki } .It can be observed from Figure 5(a) that the sum rate performance of the distributed SCA algorithm isalmost the same as its centralized counterpart for N t = 8 , K = 4 ; whereas there is a gap between the sumrates achieved by the centralized and distributed SCA algorithms for N t = 8 , K = 6 . One explanationfor this gap is that, when the system is nearly fully loaded (i.e., when K is close to N t ), the distributedSCA algorithm, which updates only the variables associated with one transmitter at a time, is more likelyto get stuck at a stationary point that is not as good as that achieved by the centralized SCA algorithmwhich optimizes all the variables in each iteration. As also shown in Fig. 5(a), when we increase N t to 12, the decentralized algorithm again converges to the centralized solution. Figure 5(b) shows that,for N t = 8 , K = 4 the distributed SCA algorithm yields performance similar to that achieved by itscentralized counterpart for almost all of the 30 tested problem instances within 10 round-robin iterations. March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 22 A v e r age S u m R a t e ( b i t s / s e c / H z ) K=4, N t =8 (centralized solution=1.571)K=6, N t =8 (centralized solution=1.573)K=6, N t =12 (centralized solution=1.612) (a) Index of Problem Instance S u m R a t e ( b i t s / s e c / H z ) Centralized7 iterations10 iterations (b)
Fig. 5: Performance of Algorithm 2, for /σ = 10 dB and η = 0 . ; (a) convergence curves versus roundnumber for N t = 8 , K = 4 , , and for N t = 12 , K = 6 , averaged over 500 sets of randomly generated { Q ki } , (b) comparison with Algorithm 1 for N t = 8 , K = 4 over 30 sets of randomly generated { Q ki } .VI. C ONCLUSIONS
In this paper, we have presented two efficient approximation algorithms for solving the rate outageconstrained coordinated beamforming design problem in (2). In view of the fact that the original designproblem involves complicated nonconvex constraints, we first presented an efficient SCA algorithm(Algorithm 1) based on SDR and first-order approximation techniques. We have shown that the proposedSCA algorithm, which involves solving convex problem (20) iteratively, can yield a stationary point ofthe outage constrained beamforming design problem, provided that problem (20) can yield a rank-onebeamforming solution. We further presented a distributed SCA algorithm (Algorithm 2) that can yieldapproximate beamforming solutions of problem (2) in a distributed, round-robin fashion, using only localCDI and a small amount of messages exchanged among the transmitters. The distributed SCA algorithmwas also shown to provide a stationary point of (2) provided that problem (36) can yield a rank-onebeamforming solution. Finally, our simulation results demonstrated that the proposed SCA algorithmyields near-optimal performance for K = 2 , and significantly outperforms the heuristic MRT and ZFschemes. Furthermore, the distributed SCA algorithm was also shown to exhibit performance comparableto its centralized counterpart within 10 rounds of round-robin iterations for most of the problem instances. March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 23 A PPENDIX AP ROOF OF C LAIM ˆ R i [ n ] = log (1 + e ¯ y i [ n − ) + e ¯ y i [ n − (ˆ y i [ n ] − ¯ y i [ n − · (1 + e ¯ y i [ n − ) ≤ log (1 + e ˆ y i [ n ] ); (A.1) similarly, from (20c), we have e ¯ x ik [ n ] = Tr( ˆ W i [ n ] Q ik ) = e ¯ x ik [ n − (ˆ x ik [ n ] − ¯ x ik [ n −
1] + 1) ≤ e ˆ x ik [ n ] , (A.2)for all k ∈ K ci , i = 1 , . . . , K . We also note from (20d) and (19) that ˆ x ii [ n ] = ¯ x ii [ n ] for all i, n . Onthe other hand, by (18), the definition of ¯ x ik [ n ] , ¯ y i [ n ] in (19), and the fact that (20b), (20f) hold withequality at the optimum, we can obtain ρ i exp( σ i e ¯ y i [ n ] − ¯ x ii [ n ] ) Y k = i (1 + e − ¯ x ii [ n ]+¯ x ki [ n ]+¯ y i [ n ] )= ρ i exp( σ i e ˆ y i [ n ] − ˆ x ii [ n ] ) Y k = i (1 + e − ˆ x ii [ n ]+ˆ x ki [ n ]+ˆ y i [ n ] ) . (A.3)Combining the above observations, i.e., (A.2), (A.3) and ˆ x ii [ n ] = ¯ x ii [ n ] , and by the monotonicity of theexponential function, we obtain that ¯ y i [ n ] ≥ ˆ y i [ n ] , which implies ˆ R i [ n ] ≤ e ˆ y i [ n ] ) ≤ e ¯ y i [ n ] ) = ˜ R i [ n ] ∀ i, n. (A.4) Suppose that e ˆ x ki [ n ] − e ¯ x ki [ n ] does not converge to zero for some i and k ∈ K ci , then there exists an ǫ > such that, for all N ≥ , e ˆ x ki [ n ] > e ¯ x ki [ n ] + ǫ for some n ≥ N . From (A.1) to (A.4), we must have e ¯ y i [ n ] > e ˆ y i [ n ] + ǫ ′ and thus ˜ R i [ n ] > ˆ R i [ n ] + ǫ ′′ , where ǫ ′ , ǫ ′′ > , which, together with (22), implies thatthe utility U ( ˜ R [ n ] , . . . , ˜ R k [ n ]) diverges as n goes to infinity however. Therefore, we must have lim n →∞ ( e ˆ x ik [ n ] − e ¯ x ik [ n ] ) = 0 ∀ i, k, (A.5) lim n →∞ ( ˜ R i [ n ] − ˆ R i [ n ]) = 0 ∀ i. (A.6) Now we use (A.5) to prove (31). It follows from (A.2) and (A.5) that lim n →∞ ( e ˆ x ik [ n ] − e ¯ x ik [ n − (ˆ x ik [ n ] − ¯ x ik [ n −
1] + 1)) = 0 (A.7) for all i and k ∈ K ci . Consider the 2nd-order Taylor series expansion [34] of e ˆ x ik [ n ] at ¯ x ik [ n − , i.e., e ˆ x ik [ n ] = e ¯ x ik [ n − (ˆ x ik [ n ] − ¯ x ik [ n −
1] + 1) + e θ [ n ]ˆ x ik [ n ]+(1 − θ [ n ])¯ x ik [ n − (ˆ x ik [ n ] − ¯ x ik [ n − , where ≤ θ [ n ] ≤ for all n ≥ . Substituting it into (A.7) gives rise to lim n →∞ e θ [ n ]ˆ x ik [ n ]+(1 − θ [ n ])¯ x ik [ n − (ˆ x ik [ n ] − ¯ x ik [ n − = 0 . Since both ¯ x ik [ n ] and ˆ x ik [ n ] are bounded by Claim 2, we conclude that (31) is true. March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 24
To show (32), we note from (A.1), (A.4) and (A.6) that lim n →∞ (cid:18) ln(1 + exp(ˆ y i [ n ])) − ln(1 + exp(¯ y i [ n − − exp(¯ y i [ n − y i [ n − y i [ n ] − ¯ y i [ n − (cid:19) = 0 . (A.8) Analogously, by considering the 2nd-order Taylor series expansion of ln(1 + e ˆ y i [ n ] ) at ¯ y i [ n − , i.e., ln(1 + e ˆ y i [ n ] ) = ln(1 + e ¯ y i [ n − ) + exp(¯ y i [ n − y i [ n − y i [ n ] − ¯ y i [ n − θ [ n ]ˆ y i [ n ] + (1 − θ [ n ])¯ y i [ n − θ [ n ]ˆ y i [ n ] + (1 − θ [ n ])¯ y i [ n − (ˆ y i [ n ] − ¯ y i [ n − , where ≤ θ [ n ] ≤ for all n ≥ , and substituting it into (A.8), we obtain lim n →∞ exp( θ [ n ]ˆ y i [ n ] + (1 − θ [ n ])¯ y i [ n − y i [ n ] − ¯ y i [ n − (1 + exp( θ [ n ]ˆ y i [ n ] + (1 − θ [ n ])¯ y i [ n − = 0 . Again, since ¯ y i [ n ] and ˆ y i [ n ] are bounded by Claim 2, we obtain (32). (cid:4) A PPENDIX BP ROOF OF T HEOREM ¯ z k [ n, i −
1] = e ¯ y k [ n,i − − ¯ x kk [ n − u k ( i − ] for all k = 1 , . . . , K . Then it can be shown that ¯ u [ n − , i ] , (cid:16) ˆ W i [ n − , { ˜ R k [ n, i − } k , { ¯ x ik [ n − } k , { ¯ y k [ n, i − } k , { ¯ z k [ n, i − } k (cid:17) , is a feasible point of (36). Hence, U ( ˆ R [ n, i ] , . . . , ˆ R K [ n, i ]) ≥ U ( ˜ R [ n, i − , . . . , ˜ R K [ n, i − forall i = 1 , . . . , K. In addition, analogous to (15), we have ˜ R j [ n, i ] ≥ ˆ R j [ n, i ] for i, j, n , and thus U ( ˜ R [ n, i ] , . . . , ˜ R K [ n, i ]) ≥ U ( ˜ R [ n, i − , . . . , ˜ R K [ n, i − , i = 1 , . . . , K, which implies that thesequence { U ( ˜ R [1 , , . . . , ˜ R K [1 , , . . . , U ( ˜ R [1 , K ] , . . . , ˜ R K [1 , K ]) , U ( ˜ R [2 , , . . . , ˜ R K [2 , , . . . } isnondecreasing. Since it is also bounded, U ( ˜ R [ n, i ] , . . . , ˜ R K [ n, i ]) , i = 1 , . . . , K , converge as n → ∞ .Now let us look at the KKT conditions of problem (36). Recall the definitions of ¯Ψ ki ( · ) and ¯Φ j ( · ) in(23) and (24) and their inner approximation properties in (25) to (30). Let Θ [ i ] i ( x ii , y i , z i , { ¯ x ki [ n − u ki ] } k = i ) , ρ i e σ i z i Y k = i (cid:16) e − x ii +¯ x ki [ n − u ki ]+ y i (cid:17) − , (A.9) Θ [ i ] j ( x ij , y j , z j , { ¯ x kj [ n − u ki ] } k = i ) , ρ j e σ j z j (cid:16) e − ¯ x jj [ n − u ji ]+ x ij + y j (cid:17) × Y k = j,k = i (cid:16) e − ¯ x jj [ n − u ji ]+¯ x kj [ n − u ki ]+ y j (cid:17) − , j ∈ K ci . (A.10)Moreover, let ˆ u [ n, i ] , ( ˆ W i [ n ] , { ˆ R k [ n, i ] } , { ˆ x ik [ n ] } k , { ˆ y k [ n, i ] } , { ˆ z k [ n, i ] } ) be the optimal solution of (36), and let λ [ n, i ] , ( λ b i [ n, i ] , { λ b k [ n, i ] } k = i , λ d [ n, i ] , { λ e k [ n, i ] } k = i , { λ f k [ n, i ] } k ,λ g i [ n, i ] , { λ g k [ n, i ] } k = i , λ P [ n, i ] , { λ δk [ n, i ] } k ) (cid:23) , March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 25 where λ b i [ n, i ] , { λ b k [ n, i ] } k = i , λ d [ n, i ] , { λ e k [ n, i ] } k = i , { λ f k [ n, i ] } k , λ g i [ n, i ] , { λ g k [ n, i ] } k = i denote the dualvariables associated with constraints in (36b) to (36g), and λ P [ n, i ] , λ δk [ n, i ] , denote the dual variablesassociated with constraint Tr( W i ) ≤ P i and Tr( W i Q ik ) ≥ δ , respectively. Let L [ i ] ( ˆ u [ n, i ] , λ [ n, i ]) bethe Lagrangian function. We can write the KKT conditions of (36) as follows: ∂ L [ i ] ( ˆ u [ n, i ] , λ [ n, i ]) ∂ W i = λ P [ n, i ] I N t − ( λ d [ n, i ] + λ δi [ n, i ]) Q ii + X k = i λ e k [ n, i ] ∂ ¯Ψ ik ( ˆ W i [ n ] , ˆ x ik [ n ] | ¯ x ik [ n − ∂ W i − λ δk [ n, i ] Q ik ! (cid:23) , (A.11a) ∂ L [ i ] ( ˆ u [ n, i ] , λ [ n, i ]) ∂R j = − ∂U ( ˆ R [ n, i ] , . . . , ˆ R K [ n, i ]) ∂R j + λ fj [ n, i ] ∂ ¯Φ j ( ˆ R j [ n, i ] , ˆ y j [ n, i ] | ¯ y j [ n, i − ∂R j ≥ ∀ j, (A.11b) ∂ L [ i ] ( ˆ u [ n, i ] , λ [ n, i ]) ∂x ii = λ b i [ n, i ] ∂ Θ [ i ] i (ˆ x ii [ n ] , ˆ y i [ n, i ] , ˆ z i [ n, i ] , { ¯ x ki [ n − u ki ] } k = i ) ∂x ii + λ d [ n, i ] e ˆ x ii [ n ] − λ g i [ n, i ] e ˆ y i [ n,i ] − ˆ x ii [ n,i ] = 0 , (A.11c) ∂ L [ i ] ( ˆ u [ n, i ] , λ [ n, i ]) ∂x ij = λ b j [ n, i ] ∂ Θ [ i ] j (ˆ x ij [ n ] , ˆ y j [ n, i ] , ˆ z j [ n, i ] , { ¯ x kj [ n − u ki ] } k = i ) ∂x ij + λ e j [ n, i ] ∂ ¯Ψ ij ( ˆ W i [ n ] , ˆ x ij [ n ] | ¯ x ij [ n − ∂x ij = 0 ∀ j ∈ K ci , (A.11d) ∂ L [ i ] ( ˆ u [ n, i ] , λ [ n, i ]) ∂y i = λ b i [ n, i ] ∂ Θ [ i ] i (ˆ x ii [ n ] , ˆ y i [ n, i ] , ˆ z i [ n, i ] , { ¯ x ki [ n − u ki ] } k = i ) ∂y i + λ f i [ n, i ] ∂ ¯Φ i ( ˆ R i [ n, i ] , ˆ y i [ n, i ] | ¯ y i [ n, i − ∂y i + λ g i [ n, i ] e ˆ y i [ n,i ] − ˆ x ii [ n ] = 0 , (A.11e) ∂ L [ i ] ( ˆ u [ n, i ] , λ [ n, i ]) ∂y j = λ b j [ n, i ] ∂ Θ [ i ] j (ˆ x ij [ n ] , ˆ y j [ n, i ] , ˆ z j [ n, i ] , { ¯ x kj [ n − u ki ] } k = i ) ∂y j + λ f j [ n, i ] ∂ ¯Φ j ( ˆ R j [ n, i ] , ˆ y j [ n, i ] | ¯ y j [ n, i − ∂y j + λ g j [ n, i ] e ˆ y j [ n,i ] − ¯ x jj [ n − u ji ] = 0 ∀ j ∈ K ci , (A.11f) ∂ L ( i ) ( u [ n, i ] , λ [ n, i ]) ∂z i = λ b i [ n, i ] ∂ Θ [ i ] i (ˆ x ii [ n ] , ˆ y i [ n, i ] , ˆ z i [ n, i ] , { ¯ x ki [ n − u ki ] } k = i ) ∂z i − λ g [ n, i ] = 0 , (A.11g) ∂ L [ i ] ( ˆ u [ n, i ] , λ [ n, i ]) ∂z j = λ b j [ n, i ] ∂ Θ [ i ] j (ˆ x ij [ n ] , ˆ y j [ n, i ] , ˆ z j [ n, i ] , { ¯ x kj [ n − u ki ] } k = i ) ∂z j − λ g j [ n, i ] = 0 , j ∈ K ci , (A.11h) and λ P [ n, i ] · (Tr( ˆ W i [ n ]) − P i ) = 0 , ∂ L [ i ] ( ˆ u [ n, i ] , λ [ n, i ]) ∂R j ˆ R j [ n, i ] = 0 ∀ j, (A.12a) ∂ L [ i ] ( ˆ u [ n, i ] , λ [ n, i ]) ∂ W i · ˆ W i [ n ] = , λ δj [ n, i ] · ( δ − Tr( ˆ W i [ n ] Q ik )) = 0 ∀ j. (A.12b) March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 26
Note that we have omitted the complementary slackness conditions for constraints (36b)-(36g) since theyare trivially satisfied at ˆ u [ n, i ] .To show the desired results, we also need the following two claims: Claim 3
It holds true that lim n →∞ | ˆ x ik [ n ] − ¯ x ik [ n − | = 0 ∀ i, k, (A.13a) lim n →∞ | ˆ x ik [ n ] − ¯ x ik [ n ] | = 0 ∀ i, k, (A.13b) lim n →∞ | ˆ y k [ n, − ¯ y k [ n − , K ] | = 0 , lim n →∞ | ˆ y k [ n, i ] − ¯ y k [ n, i − | = 0 ∀ i, k, (A.13c) lim n →∞ | ˆ y k [ n, i ] − ¯ y k [ n, i ] | = 0 ∀ i, k, (A.13d) lim n →∞ | ˜ R k [ n, − ˜ R k [ n − , K ] | = 0 , lim n →∞ | ˜ R k [ n, i ] − ˜ R k [ n, i − | = 0 ∀ i, k, (A.13e) lim n →∞ | ˆ R k [ n, i ] − ˜ R k [ n, i ] | = 0 ∀ i, k. (A.13f) Claim 4
For each i , ˆ u [ n, i ] generated by Algorithm 2 is bounded for all n . The proof of Claim 3 is presented in Appendix C. Similar to Claim 1, (A.13a) to (A.13d) imply that the re-strictive approximations in (36e) and (36f) are asymptotically tight as n → ∞ . Since problem (36) satisfiesthe Slater’s condition, the dual variable vector λ [ n, i ] is bounded [33]. Moreover, ˆ u [ n, i ] is also bounded byClaim 4. Now let us consider the primal-dual solution pair ( ˆ u [ n, i ] , λ [ n, i ]) for all i = 1 , . . . , K . Since theyare all bounded, there exists a subsequence { n , . . . , n ℓ , . . . } ⊆ { , . . . , n, . . . } and limit points ˆ u ⋆ [ i ] , ( ˆ W ⋆i , { ˆ R ⋆k [ i ] } , { ˆ x ⋆ik } k , { ˆ y ⋆k [ i ] } , { ˆ z ⋆k [ i ] } ) and λ ⋆ [ i ] , ( λ b ⋆i [ i ] , { λ b ⋆k [ i ] } k = i , λ d ⋆ [ i ] , { λ e ⋆k [ i ] } k = i , { λ f ⋆k [ i ] } k , λ g ⋆i [ i ] , { λ g ⋆k [ i ] } k = i , λ P ⋆ [ i ] , { λ δ⋆k [ i ] } k ) (cid:23) for all i , such that lim ℓ →∞ ˆ u [ n ℓ , i ] = ˆ u ⋆ [ i ] , lim ℓ →∞ λ [ n ℓ , i ] = λ ⋆ [ i ] (A.14)for i = 1 , . . . , K . By (A.13e) and (A.13f), we see that both ˆ R k [ n ℓ , i ] and ˜ R k [ n ℓ , i ] converge to the samelimit point, and they are the same for all i , i.e., ˆ R ⋆k [1] = ˆ R ⋆k [2] = · · · = ˆ R ⋆k [ K ] , ˜ R ⋆k , k = 1 , . . . , K. (A.15)Analogously, by (A.13a) to (A.13d), we have that ˆ y ⋆k [1] = ˆ y ⋆k [2] = · · · = ˆ y ⋆k [ K ] , ˆ y ⋆k , k = 1 , . . . , K, (A.16) ˆ z ⋆k [1] = ˆ z ⋆k [2] = · · · = ˆ z ⋆k [ K ] , ˆ z ⋆k , k = 1 , . . . , K. (A.17) March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 27
Then, it follows from the inner approximation properties in (25) to (30), (A.13a), (A.13c), and (A.14)to (A.17) that the KKT conditions in (A.11) and (A.12) converge along the subsequence { n , . . . , n ℓ , . . . } to ∂ L [ i ] ( ˆ u ⋆ [ i ] , λ ⋆ [ i ]) ∂ W i = λ P ⋆ [ i ] I N t − ( λ d ⋆ [ i ] + λ δ⋆i [ i ]) Q ii + X k = i λ e ⋆k [ i ] ∂ Ψ ik ( ˆ W ⋆i , ˆ x ⋆ik ) ∂ W i − λ δ⋆k [ i ] Q ik ! (cid:23) , (A.18a) ∂ L [ i ] ( ˆ u ⋆ [ i ] , λ ⋆ [ i ]) ∂R j = − ∂U ( ˜ R ⋆ , . . . , ˜ R ⋆K ) ∂R j + λ f ⋆ j [ i ] ∂ Φ j ( ˜ R ⋆j , ˆ y ⋆j ) ∂R j ≥ ∀ j, (A.18b) ∂ L [ i ] ( ˆ u ⋆ [ i ] , λ ⋆ [ i ]) ∂x ii = λ b ⋆i [ i ] ∂ Θ [ i ] i (ˆ x ⋆ii , ˆ y ⋆i , ˆ z ⋆i , { ˆ x ⋆ki } k = i ) ∂x ii + λ d ⋆ [ i ] e ˆ x ⋆ii − λ g ⋆i [ i ] e ˆ y ⋆i − ˆ x ⋆ii = 0 , (A.18c) ∂ L [ i ] ( ˆ u ⋆ [ i ] , λ ⋆ [ i ]) ∂x ij = λ b ⋆j [ i ] ∂ Θ [ i ] j (ˆ x ⋆ij , ˆ y ⋆j , ˆ z ⋆j , { ˆ x ⋆kj } k = i ) ∂x ij + λ e ⋆j [ i ] ∂ Ψ ij ( ˆ W ⋆i , ˆ x ⋆ij ) ∂x ij = 0 ∀ j ∈ K ci , (A.18d) ∂ L [ i ] ( ˆ u ⋆ [ i ] , λ ⋆ [ i ]) ∂y j = λ b ⋆j [ i ] ∂ Θ [ i ] j (ˆ x ⋆ij , ˆ y ⋆j , ˆ z ⋆j , { ˆ x ⋆kj } k = i ) ∂y j + λ f ⋆j [ i ] ∂ Φ j ( ˜ R ⋆j , ˆ y ⋆j ) ∂y j + λ g ⋆j [ i ] e ˆ y ⋆j − ˆ x ⋆jj = 0 ∀ j, (A.18e) ∂ L [ i ] ( ˆ u ⋆ [ i ] , λ ⋆ [ i ]) ∂z j = λ b ⋆j [ i ] ∂ Θ [ i ] j (ˆ x ⋆ij , ˆ y ⋆j , ˆ z ⋆j , { ˆ x ⋆kj } k = i ) ∂z j − λ g ⋆j [ i ] = 0 ∀ j, (A.18f) and λ i ⋆ [ i ] · (Tr( ˆ W ⋆i ) − P i ) = 0 , ∂ L [ i ] ( ˆ u ⋆ [ i ] , λ ⋆ [ i ]) ∂R j ˜ R ⋆j = 0 ∀ j, (A.19a) ∂ L [ i ] ( ˆ u ⋆ [ i ] , λ ⋆ [ i ]) ∂ W i · ˆ W ⋆i = , λ δ⋆j [ i ] · ( δ − Tr( ˆ W ⋆i Q ik )) = 0 ∀ j. (A.19b) It can be observed from (39) that, for ρ i < , ˜ R j [ n, i ] is strictly greater than zero for all i, j, n ; therefore, ˜ R ⋆j > for all j , which indicates that ∂ L [ i ] ( ˆ u ⋆ , ˆ λ ⋆ [ i ]) ∂R j = 0 for all i, j by (A.19a). Substituting this into(A.18b) for all i = 1 , . . . , K , gives rise to λ f ⋆j [1] = · · · = λ f ⋆j [ K ] = ∂U ( ˜ R ⋆ , . . . , ˜ R ⋆K ) ∂R j ∂ Φ j ( ˜ R ⋆j , ˆ y ⋆j ) ∂R j ! − , λ f ⋆j , j = 1 , . . . , K. (A.20)In addition, one can verify that ∂ Θ [1] j (ˆ x ⋆ j , ˆ y ⋆j , ˆ z ⋆j , { ˆ x ⋆kj } k =1 ) ∂y j = · · · = ∂ Θ [ K ] j (ˆ x ⋆Kj , ˆ y ⋆j , ˆ z ⋆j , { ˆ x ⋆kj } k = K ) ∂y j ,∂ Θ [1] j (ˆ x ⋆ j , ˆ y ⋆j , ˆ z ⋆j , { ˆ x ⋆kj } k =1 ) ∂z j = · · · = ∂ Θ [ K ] j (ˆ x ⋆Kj , ˆ y ⋆j , ˆ z ⋆j , { ˆ x ⋆kj } k = K ) ∂z j , which, together with (A.18e) (A.18f) and (A.20), lead to λ b ⋆j [1] = · · · = λ b ⋆j [ K ] , λ b ⋆j , λ g ⋆j [1] = · · · = λ g ⋆j [ K ] , λ g ⋆j ∀ j. (A.21)Finally, by (A.18), (A.19), (A.20) and (A.21), we conclude that ( { ˆ W ⋆i } , { ˜ R ⋆k } , { ˆ x ⋆ik } k , { ˆ y ⋆k } , { ˆ z ⋆k } ) and ( { λ b ⋆i } , { λ d ⋆ [ i ] } , {{ λ e ⋆k [ i ] } k = i } i , { λ f ⋆k } , { λ g ⋆i } , { λ P ⋆ [ i ] } , { λ δ⋆k [ i ] } ) satisfy the KKT conditions of problem(8). The proof is completed. (cid:4) March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 28 A PPENDIX CP ROOF OF C LAIM e ¯ x ik [ n ] ≤ e ˆ x ik [ n ] , ˆ R k [ n, i ] ≤ e ˆ y k [ n,i ] ) (A.22)for all k ∈ K ci , i = 1 , . . . , K . Also by (36c), (36g) and (39), we have ρ j exp( σ j e ˆ y j [ n,i ] − ¯ x jj [ n − u ji ] ) (cid:16) e − ¯ x jj [ n − u ji ]+ˆ x ij [ n ]+ˆ y j [ n,i ] (cid:17) Y k = jk = i (cid:16) e − ¯ x jj [ n − u ji ]+¯ x kj [ n − u ki ]+ˆ y j [ n,i ] (cid:17) = ρ j exp( σ j e ¯ y i [ n,i ] − ¯ x ii [ n − u ji ] ) (cid:16) e − ¯ x jj [ n − u ji ]+¯ x ij [ n ]+¯ y j [ n,i ] (cid:17) Y k = jk = i (cid:16) e − ¯ x jj [ n − u ji ]+¯ x kj [ n − u ki ]+¯ y j [ n,i ] (cid:17) , for all j ∈ K ci . Using the above equation and (A.22) and the monotonicity of exponential function, weobtain ˆ y j [ n, i ] ≤ ¯ y j [ n, i ] . Thus, ˆ R j [ n, i ] ≤ e ˆ y j [ n,i ] ) ≤ e ¯ y j [ n,i ] ) = ˜ R j [ n, i ] ∀ j ∈ K ci , (A.23)Similarly, by (36b), (36g) and (39), we have ˆ R i [ n, i ] ≤ e ˆ y i [ n,i ] ) = 1ln 2 ln(1 + e ¯ y i [ n,i ] ) = ˜ R i [ n, i ] . (A.24)Using the same arguments as in obtaining (A.3) to (A.6) in Appendix A, we can show that (A.13f),(A.13b), (A.13d), (A.13c) and lim n →∞ | ˆ x ik [ n ] − ¯ x ik [ n − | = 0 ∀ i, k ∈ K ci , (A.25)which is (A.13a) for k = i , are true. What remains is to prove (A.13e) and lim n →∞ | ˆ x ii [ n ] − ¯ x ii [ n − | = 0 ∀ i .It follows from (A.13c), (A.13d) and the triangle inequality that lim n →∞ | ¯ y k [ n, − ¯ y k [ n − , K ] | = 0 , lim n →∞ | ¯ y k [ n, i ] − ¯ y k [ n, i − | = 0 ∀ i, k, which, by the definition in (38), is equivalent to (A.13e). By considering (39) for transmitter i − , andthe fact that (36b) holds with equality at the optimal point for transmitter i , we can obtain ρ i exp( σ i e ¯ y i [ n,i − − ¯ x ii [ n − ) Y k = i (1 + e − ¯ x ii [ n − x ki [ n − u ki ]+¯ y i [ n,i − )= ρ i exp( σ i e ˆ y i [ n,i ] − ˆ x ii [ n ] ) Y k = i (1 + e − ˆ x ii [ n ]+¯ x ki [ n − u ki ]+ˆ y i [ n,i ] ) ∀ i. March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 29
Since both { ˆ y i [ n, i ] } ∞ n =1 and { ¯ y i [ n, i − } ∞ n =1 are bounded, and by (A.13c), we obtain from the aboveequation that lim n →∞ | ˆ x ii [ n ] − ¯ x ii [ n − | = 0 ∀ i. Thus the proof of Claim 3 has been completed. (cid:4) R EFERENCES [1] W.-C. Li, T.-H. Chang, C. Lin, and C.-Y. Chi, “A convex approximation approach to weighted sum rate maximization ofmultiuser MISO interference channel under outage constraints,” in
Proc. IEEE ICASSP , Progue, Czech, May 22-27, 2011,pp. 3368–3371.[2] H. Zhang, N. B. Mehta, A. F. Molisch, J. Zhang, and H. Dai, “Asynchronous interference mitigation in cooperative basestation systems,”
IEEE Trans. Wireless Commun. , vol. 7, pp. 155–165, Jan. 2008.[3] D. Gesbert, S. Hanly, H. Huang, S. S. Shitz, O. Simeone, and W. Yu, “Multi-cell MIMO cooperative networks: A newlook at interference,”
IEEE J. Sel. Areas Commun. , vol. 28, pp. 1380–1408, Dec. 2010.[4] E. Bj¨ornson, N. J. Jald´en, M. Bengtsson, and B. Ottersten, “Optimality properties, distributed strategies, and measurement-based evaluation of coordinated multicell OFDMA transmission,”
IEEE Trans. Signal Process. , vol. 59, pp. 6086–6101,Dec. 2011.[5] H. Dahrouj and W. Yu, “Coordinated beamforming for the multicell multi-antenna wireless system,”
IEEE Trans. WirelessCommun. , vol. 9, pp. 1748–1759, May 2010.[6] L. Ventruino, N. Prasad, and X.-D. Wang, “Coordinated linear beamforming in downlink multi-cell wireless networks,”
IEEE Trans. Wireless Commun. , vol. 9, pp. 1451–1461, Apr. 2010.[7] ——, “Coordinated scheduling and power allocation in downlink multicell OFDMA networks,”
IEEE Trans. VehicularTech. , vol. 58, pp. 2835–2848, July 2009.[8] A. B. Carleial, “Interference channels,”
IEEE Trans. Inf. Theory , vol. 24, pp. 60–70, Jan. 1978.[9] X. Shang and B. Chen, “Achievable rate region for downlink beamforming in the presence of interference,” in
Proc.Asilomar Conference on Signals, Systems, and Computers , Pacific Grove, CA, Nov. 4-7, 2007, pp. 1684–1688.[10] V. S. Annapureddy and V. V. Veeravalli, “Sum capacity of MIMO interference channels in the low interference regime,”
IEEE Trans. Inf. Theory , vol. 57, pp. 2565–2581, May 2011.[11] E. G. Larsson, E. A. Jorswieck, J. Lindblom, and R. Mochaourab, “Game theory and the flat-fading Gaussian interferencechannel,”
IEEE Signal Process. Mag. , vol. 26, pp. 18–27, Sep. 2009.[12] E. A. Jorswieck, E. G. Larsson, and D. Danev, “Complete characterization of the Pareto boundary for the MISO interferencechannel,”
IEEE Trans. Signal Process. , vol. 56, pp. 5292–5296, July 2008.[13] X. Shang, B. Chen, and H. V. Poor, “Multiuser MISO interference channels with single-user detection: Optimality ofbeamforming and the achievable rate region,”
IEEE Trans. Inf. Theory , vol. 57, pp. 4255–4273, July 2011.[14] R. Mochaourab and E. A. Jorswieck, “Optimal beamforming in interference networks with perfect local channelinformation,”
IEEE Trans. Signal Process. , vol. 59, pp. 1128–1141, Mar. 2011.[15] Y.-F. Liu and Z.-Q. Luo, “Coordinated beamforming for MISO interference channel: Complexity analysis and efficientalgorithms,”
IEEE Trans. Signal Process. , vol. 59, pp. 1142–1157, Mar. 2011.[16] R. Zakhour and D. Gesbert, “Coordination on the MISO interference channel using the virtual SINR framework,” in
Proc.Int. ITG Workshop Smart Antennas , Berlin, Germany, Feb. 16-18, 2009.
March 3, 2018 DRAFTCCEPTED BY IEEE TRANSACTIONS ON SIGNAL PROCESSING, NOV. 2012 30 [17] R. Zhang and S. Cui, “Cooperative interference management with MISO beamforming,”
IEEE Trans. Signal Process. ,vol. 58, pp. 5450–5458, Oct. 2010.[18] E. G. Larsson and E. A. Jorswieck, “Competition versus cooperation on the MISO interference channel,”
IEEE J. Sel.Areas Commun. , vol. 26, pp. 1059–1069, Sep. 2008.[19] D. A. Schmidt, C. Shi, R. A. Berry, M. L. Honig, and W. Utschick, “Distributed resource allocation schemes: Pricingalgorithms for power control and beamformer design in interference networks,”
IEEE Signal Process. Mag. , vol. 26, pp.53–63, Sep. 2009.[20] E. Bj¨ornson, R. Zakhour, D. Gesbert, and B. Ottersten, “Cooperative multicell precoding: Rate region characterization anddistributed strategies with instantaneous and statistical CSI,”
IEEE Trans. Signal Process. , vol. 58, pp. 4298–4310, Aug.2010.[21] E. Karipidis, A. Gr¨undinger, J. Lindblom, and E. G. Larsson, “Pareto-optimal beamforming for the MISO interferencechannel with partial CSI,” in
Proc. 3rd IEEE International Workshop on Computational Advances in Multi-Sensor AdaptiveProcessing , Aruba, Dutch Antilles, Dec. 13-16, 2009, pp. 5–8.[22] J. Lindblom, E. Karipidis, and E. G. Larsson, “Outage rate regions for the MISO IFC,” in
Proc. Asilomar Conference onSignals, Systems and Computers , Pacific Grove, CA, Nov. 1-4, 2009, pp. 1120–1124.[23] ——, “Outage rate regions for the MISO interference channel: Definitions and interpretations,” http://arxiv.org/abs/1106.5615v1 .[24] S. Kandukuri and S. Boyd, “Optimal power control in interference-limited fading wireless channels with outage-probabilityspecifications,”
IEEE Trans. Wireless Commun. , vol. 1, pp. 46–55, Jan. 2002.[25] S. Ghosh, B. D. Rao, and J. R. Zeidler, “Outage-efficient strategies for multiuser MIMO networks with channel distributioninformation,”
IEEE Trans. Signal Process. , vol. 58, pp. 6312–6324, Dec. 2010.[26] Z.-Q. Luo, W.-K. Ma, A. M.-C. So, Y. Ye, and S. Zhang, “Semidefinite relaxation of quadratic optimization problems,”
IEEE Signal Process. Mag. , vol. 27, pp. 20 –34, May 2010.[27] J. Mo and J. Walrand, “Fair end-to-end window-based congestion control,”
IEEE/ACM Trans. Networking , vol. 8, pp.556–567, Oct. 2000.[28] T. Bonald and L. Massouli´e, “Impact of fairness on internet performance,” in
Proc. ACM SIGMETRICS , Cambridge, MA,June 16-20, 2001, pp. 82–91.[29] Z.-Q. Luo and T.-H. Chang, “SDP relaxation of homogeneous quadratic optimization: Approximation bounds andapplications,” in
Convex Optimization in Signal Processing and Communications , D. P. Palomar and Y. C. Eldar, Eds.Cambridge, U.K.: Cambridge Univ. Press, 2010, ch. 4.[30] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 1.21,” http://cvxr.com/cvx,Apr. 2011.[31] Y. Huang and D. P. Palomar, “Rank-constrained separable semidefinite programming with applications to optimalbeamforming,”
IEEE Trans. Signal Process. , vol. 58, pp. 664–678, Feb. 2010.[32] B. R. Marks and G. P. Wright, “A general inner approximation algorithm for nonconvex mathematical programs,”
OperationsResearch , vol. 26, pp. 681–683, 1978.[33] D. P. Bertsekas, A. Nedi ´ c, and A. E. Ozdaglar, Convex Analysis and Optimization . Cambridge, MA: Athena Scientific,2003.[34] R. A. Gamboa and B. Middleton, “Taylor’s formula with remainder,” in
Proc. 3rd International Workshop of the ACL2Theorem Prover and Its Applications , Grenoble, France, April 8-9, 2002., Grenoble, France, April 8-9, 2002.