{'fraction_non_alphanumeric': 0.04536627923338868, 'fraction_numerical': 0.010947569864632876, 'mean_word_length': 4.169333721268548, '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': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Energy Efficiency Maximization in Large-Scale Cell-Free Massive MIMO: A Projected Gradient Approach

This paper considers the fundamental power allocation problem in cell-free massive mutiple-input and multiple-output (MIMO) systems which aims at maximizing the total energy efficiency (EE) under a sum power constraint at each access point (AP) and a quality-of-service (QoS) constraint at each user. Existing solutions for this optimization problem are based on solving a sequence of second-order cone programs (SOCPs), whose computational complexity scales dramatically with the network size. Therefore, they are not implementable for practical large-scale cell-free massive MIMO systems. To tackle this issue, we propose an iterative power control algorithm based on the frame work of an accelerated projected gradient (APG) method. In particular, each iteration of the proposed method is done by simple closed-form expressions, where a penalty method is applied to bring constraints into the objective in the form of penalty functions. Finally, the convergence of the proposed algorithm is analytically proved and numerically compared to the known solution based on SOCP. Simulations results demonstrate that our proposed power control algorithm can achieve the same EE as the existing SOCPs-based method, but more importantly, its run time is much lower (one to two orders of magnitude reduction in run time, compared to the SOCPs-based approaches).

[2]- [7]. Basically, cell-free massive MIMO uses a massive number of distributed access points (APs) together with simple linear processing to coherently serve many users using the same time and frequency resources. Each APs can be equipped with several antennas. It relies on the favorable propagation and channel hardening property of massive MIMO technology [8], and the macro-diversity of network MIMO technique [9]. Therefore, it can provide universally good service to all users in the network regardless their locations.
In cell-free massive MIMO, because the APs and users are distributed over a large area, power controls are very important to control the near-far effect, and hence, can significantly improve the system performance as well as to save the radiated powers from the APs in the downlink and the users in the uplink. Thus, many research works on power allocations in cell-free massive MIMO have been studied [2], [3], [10]- [13]. In [2], [3], power control coefficients at the APs and users were optimally chosen to maximize the minimum spectral efficiency of all users. In [10], the max-min power power control under limited backhaul was investigated. The downlink transmission power optimization of cell-free massive MIMO with spatially correlated Rayleigh fading channels for noncoherent joint and coherent joint transmission was exploited in [11]. A joint downlink transmit powers and the number of active APs optimization was proposed and solved in [12]. In [13], the total energy efficiency maximization taken into account the hardware and backhaul power consumption was proposed and exploited. In the context of cell-free massive MIMO, the solutions to power control problems in most, if not all, of previous work are based on successive convex approximation principle, which approximates a non-convex problem by a sequence of convex second order cone programs (SCOPs) [14]. As a result, these methods have very high computational complexity, as they rely on interior point methods (through the use of off-the-shelf convex solvers) to solve these convex problems. Therefore, they are not implementable for large-scale cell-free massive MIMO with many APs and users (e.g. in stadium or shopping malls where we may have thousands of APs and active users).
To deal with the large-scale problem in cell-free massive MIMO, [15] proposed a scalable framework, which uses AP selection to cope with computational complexity and backhaul requirements. However, the power allocation is quite simple and heuristic, and thus, it may underestimate the capacity of the system. Another approach to deal with the large-scale problem is presented in [16]- [18]. The main idea of this approach is to decompose a large optimization problem into smaller optimization subproblems to reduce computational complexity. However, this method cannot be applied to solve optimization problems in which the variables are coupled, such as those with quality of service (QoS) constraints.
In this paper, we consider the energy efficiency maximization problem with QoS constraints with an emphasis on large-scale settings. For such scenarios, the number of power control coefficients (i.e., the product of the number of APs and the number of users) which need to be optimized can be extremely large. Thus, our goal is to propose a novel power control algorithm for the energy efficiency maximization problem, which has much lower computational complexity compared to the traditional SOCP-based method. It is apparent from the above discussions that an efficient numerical method for this particular problem is still demanding. To this end, we combine the penalty method and the accelerated projected gradient (APG) method. More specifically, the penalty method is used to handle the QoS constraints in our problem, resulting in more tractable subproblems. Note that the penalty method is widely used to deal with constrained optimization [19]. In principle, the penalty method penalizes a set of constraints by proper terms and adds the penalty terms into the objective, creating the so-called penalized objective. In this way, an optimal problem with sophisticated constraints can be converted into a regularized optimization problem with simple constraints for which efficient solutions are easier to derive. By increasing the penalty parameter, the solutions to these regularized problems converge to a solution of the original problem. In this paper, to solve the regularized problems obtained from the penalty method, we then apply the APG method which is a variant of the accelerated proximal gradient method proposed in [20] for nonconvex programming. As shall be numerically shown in Section V, compared to the sequential SOCP-based method, the proposed method achieves the same total EE but with much lower run time and computational complexity since it is entirely based on first order oracle (i.e the value of the objective and its gradient). Thus, the proposed method can be readily modified to tackle the high complexity of other resource allocation problems for large-scale cell-free massive MIMO. The main contributions of this paper are as follows.
• We provide the mathematical background of the APG method with a detailed proof as an alternative solution for sequential SOCP-based method to deal with many resource allocation problems in large-scale cell-free massive MIMO.
• In our proposed APG method, no external optimization solver is needed as the projection, which is the main operation of the proposed method, is done by closed form expression.
As the result, it is much faster to output a solution, compared to the known sequential SOCP-based method.
• We customize the presented APG method to solve the total EE maximization problem in cell-free massive MIMO.
• For our specific problem, we first transform the problem of total EE maximization, subject to transmit power constraints at APs and the individual quality-of-service (QoS) constraints at each user, into a form amenable to the application of the APG method.
• We then combined the penalty method and the APG method to achieve a good and lowcomplexity power control algorithm. In particular, the QoS constraints are penalized by a proper smooth penalty term which in controlled by a penalty parameter. The penalty term is then added to the original objective, giving rise to the penalized problem. The APG is applied to solve the penalized problem, whereby each iteration admits closed-form expressions. The computational complexity of the proposed algorithm is provided.
• We provide numerical results to show that the proposed algorithm can achieve the same performance as an SCOPs-based method but with much reduced run time.
• We also verify that our proposed algorithm converges to a feasible solution regardless the choice of starting point.
The rest of this paper is organized as follows. Section II provide the preliminaries of the APG method. Next, Section III recalls total EE optimization problem in cell-free massive MIMO.
Then, Section IV proposes to use APG method for total EE optimization problem. Section V evaluates the system performance by using numerical results. Finally, the conclusion is drawn in Section VI.
Notation: Standard notations are used in this paper. The superscripts (·) T and (·) H stand for the transpose and the Hermitian, respectively. Notation [x] + denotes the projector onto the positive orthant. Notation ∇, and ∂ denote the gradient, and sub-gradient, respectively. Notation domf denotes a domain of function f , and ||.|| denotes the l 2 -norm. Notation ≡ is used to define an equivalent quantity. Finally, we use ⊙, ⊗, and ., . to denote the Hadamard, the Kronecker products and an inner product, respectively.

II. MATHEMATICAL PRELIMINARIES: ACCELERATED PROJECTED GRADIENT METHOD
In this section, we first provide the general framework of an accelerated proximal gradient method and then present a variant of the accelerated proximal gradient method, which is termed the accelerated projected gradient (APG) method, to deal with the EE maximization problem to be considered in Section III. First, we recall some definitions. A function f is said to be proper the function is said to be bounded from below by x 0 . If f (x) ≥ x 0 , ∀x ∈ X, then the function is said to be bounded from below by x 0 . A function f is lower semicontinuous at point We now present a general mathematical framework, called the accelerated proximal gradient method for non-convex problems presented in [20], which concerns the following optimization problem where f (x) is L-Lipschitz continuous gradient, g(x) is proper and lower semicontinuous, and T (x) is coercive. Then the accelerated proximal gradient method for solving (1), consists of the following iterations: where α x , α y are step sizes, and prox αg is the proximal operator defined as In general, accelerated proximal gradient method in (2) is designed to cope with the unconstrained optimization problem (1). However, most of resource allocation problems in cell-free massive MIMO are constrained optimization problems. Therefore, in this paper, we present a special case of the accelerated proximal gradient method, which is called the APG method, to deal with those problems. Specifically, we consider the following optimization problem where C is the feasible set of the considered problem, which is often defined by a set of constraints. Again, assume that f (x) is a proper function with Lipschitz continuous gradient, and bounded from below. We remark that g(x) in (1) is not necessarily smooth. Thus, to apply the iterations in (2) to solve (4), we can let g(x) in (1) be the indicator function of the feasible set C. In this way, the proximal operator in (3) reduces to the Euclidean projection onto C [20].
As a result, the APG method for solving (4), consists of the following iterations: where P C (x) denotes the Euclidean projection of x onto C, which is defined as

III. TOTAL ENERGY EFFICIENCY OPTIMIZATION IN CELL-FREE MASSIVE MIMO
In this section, we first briefly introduce the system model of cell-free massive MIMO and formulate the total EE optimization problem with conjugate beamforming at the APs, taking into account arbitrary pilot sequence assignments, and imperfect channel estimation. Then, by using penalty functions (PFs), we reformulate the total EE optimization problem into the form so that the APG method can be applied.

1) Spectral Eficiency:
We consider a cell-free massive MIMO downlink, which is shown in between AP m and user k is modeled as where β mk is the large-scale fading, and h mk is the small-scale fading, whose elements are i.i.d. CN(0, 1) RVs. The downlink transmission needs two phases: uplink training and downlink payload data transmission phases. In the training phase, all users send their pilot sequences, √ τ p ϕ k ∈ C τp×1 , ∀k, where ϕ k 2 = 1, to all APs in the system. Then, pilot signal received at AP m is where ρ p is the normalized transmit signal-to-noise ratio (SNR) of each pilot symbol, W p,m is the noise matrix whose elements are i.i.d. CN(0, 1) RVs. After that, each AP uses its received pilot signals from all K users to estimate its local channels using the minimum mean-square error (MMSE) technique [21]. The channel estimate of g mk iŝ In the downlink data transmission phase, APs use conjugate beamforming technique and the channels estimated in the training phase to precode the symbols intended for all users. The vector of transmitted signals from AP m is where q k is the symbol intended for user k which satisfies E |q k | 2 = 1, ∀k, ρ d is the normalized transmit power at each AP, and η mk is power coefficient between the AP m and the user k, where Then the signal received at user k is given by By applying the same technique as the one in [13], the downlink achievable spectral efficiency (SE) of the k-th user can be expressed as where τ c is the length of each coherence interval, τ p is the length of the training phase for each co-

2) Power Consumption Model:
In this paper, the total power consumption is modeled as [13] where P m is the power consumption at the m-th AP, and P bh,m is the power consumed by the backhaul link connecting the CPU and the m-th AP. Specifically, P m is modeled as where 0 ≤ α m ≤ 1 is the power amplifier efficiency, N 0 is the noise power, and P tc,m is the internal power required to run the circuit components at each antenna of the m-th AP. Next, P bh,m is modeled as where P 0,m is a fixed power consumption of each backhaul, P bt,m is the traffic-dependent power, and B is the system bandwidth.

3) Total Energy Efficiency :
While spectral efficiency has been a common performance measure for wireless communication design, in this paper we aim to maximize the total energy efficiency of the system, which is defined as how many bits can be transmitted by one Joule.
Specifically, the total EE (bit/Joule) can be calculated as

B. Optimization Problem Formulation
Our problem is to maximize the total EE (19) by allocating the power coefficients {η mk }, under a sum power constraint at each AP and a QoS constraint, i.e., SE constraint S ok at each user. The optimization problem is stated as In this paper (P) is assumed to be feasible. An equivalent form of problem (P) can be rewritten as (P 1 ) : whereP fix M m=1 (NP tc,m +P 0,m ). Note that the objective function of the problem (P 1 ) is nonconvex. The common method to tackle such such a nonconvex problem is to iteratively approximate a nonconvex function by a series of convex functions under the framework of successive convex approximation. In fact, this is the method presented in [13], in which the EE maximization problem is solved by a sequence of SOCPs. However, as the complexity dramatically increases when the system scales up (i.e. the numbers of APs and users increase), such method cannot provide a solution for large-scale optimization problems in cell-free massive MIMO systems with thousands of APs and users. In the next section, we will propose a new algorithm based on the APG method to solve problem (P 1 ). Our proposed algorithm has very low complexity, and hence, and can efficiently deal with the systems with many APs and users.

IV. PROPOSED APG ALGORITHM FOR THE TOTAL ENERGY EFFICIENCY OPTIMIZATION
In this section, we apply the APG method to efficiently solve the total energy efficiency optimization problem in Section III. We first reformulate the optimization problem (P 1 ) by change of variables such that the gradient of the function and the resulting projection can be computed more efficiently. We also apply a penalty method to convert the total EE maximization problem into the form which is amenable to applying the APG method. Note that our algorithm is done over large-scale fading time scale and is performed at the CPU. The details of these steps are described in the following.

A. Problem Reformulation
In order to apply the APG method for solving problem (P 1 ), we need to reformulate (P 1 ) into the form of (4) and to make sure that the gradient of the objective is Lipschitz continuous. To this end, we first introduce a new variable θ mk = √ η mk γ mk and define new notations as follows: is the vector of all power control coefficients associated with AP m.
where e k ∈ R K denotes the k-th unit vector, i.e., the vector such that e k = 1 and e j = 0, ∀j = k.
Then, (P 1 ) can be rewritten as where and As shall be seen shortly, the projection onto C can be done by closed form expressions. Thus, the main obstacle in deriving an efficient algorithm for solving (P 2 ) is the QoS constraints.
To overcome this issue we propose to combine the penalty method and the APG method as described in the next subsection.

B. Proposed Algorithm
The overall structure of the proposed method is as follows: • The penalty method is invoked to bring the QoS constraints into the objective by some form of a loss function through a penalty parameter, leading to the penalized problem.
• The APG method is then applied to solve the penalized problem. This process is repeated until a stopping criterion is achieved.

1) Penalty Method:
The constraint u k (θ θ θ) ≥ S ok can be written as where a k 2 So k τc . Then, for each QoS constraint we introduce the following quadratic loss function [19] where Note that g k (θ θ θ) is convex and Ψ k (θ θ θ) is smooth. Then, for a given penalty coefficient ξ, the penalized objective function of (P 2 ), denoted by f ξ (θ θ θ), is given by We remark that the above regularized objective is formed in the context of maximization. Also note that the value of the penalty coefficient ξ should be selected appropriately. If this parameter is large, the feasibility is guaranteed but the resulting optimization problem is numerically illconditioned. On the other hand, if it is too small, it may produce a suboptimal solution or even converges to an infeasible point, i.e. the constraints are violated.
In practice, to avoid the above issues, we can first solve the penalized optimization problem for a small value of ξ and check if the stopping criterion is met. If not, we can increase ξ by ρ > 1 times and repeat this process until the stopping criterion is met. In this iterative process, it is critical to use the solution of the previous iteration as the starting point of the next. In essence, the key to the penalty method is to solve the following regularized optimization problem for a given ξ which has the same form as (4). We are now in a position to apply the APG method to solve (31) which is detailed next section.

2) APG Method:
We first show that f ξ (θ θ θ) is a proper function with Lipschitz continuous gradient and bounded from above, and thus, the APG method is applicable to solve (31). Towards this end, It is easy to see that the function f (θ θ θ) is proper and bounded from above 1 , which is shown by the following inequalities The above inequality holds since the total EE is bounded from above as its numerator, i.e. total SE, is limited by inter-user interference and total power consumption. The Lipschitz continuity of the gradient of f ξ (θ θ θ) is stated in Proposition 1.
Lipschitz contant L f given by (84) in Appendix B.
Proof. See Appendix B.
we note that we have modified (5), accounting for the maximization context, where we move along the gradient to increase the objective of the current point. Note also that, for a practical purpose we stop the APG procedure when the relative increase in the objective during the last 10 iterations is less than ς.
It is clear that the key operations in the implementation of Algorithm 1 are the computation of the gradient ∇f ξ (θ θ θ) and the projections in (32a) and (32b). In particular, these two operations can be done in closed-form as shown in Proposition 2 and Proposition 3, respectively. Proposition 2. ∇f ξ (θ θ θ) can be calculated as where and where Proof: See Appendix C.
Proof: See Appendix D.

C. Proposed Algorithm with Line Search
In Algorithm 1 is guaranteed to converge for any fixed step sizes smaller than L f . However, it is possible that L f given in (84) is significantly larger than the best Lipschitz constant of the gradient of f ξ (θ θ θ) which is practically difficult to find. In order to find a larger step size, and thus faster convergence, we can carry out a line search to tune the step size in (32a) and (32b). In this paper, we can perform a line search as described in Algorithm 2, inspired from [20], which works as follows. In each iteration, the backtracking line search starts with a large Input: ν < 1, δ > 0 / * step size for (32a) * / z (n+1) = P C (y (n) + α y ∇f ξ (y (n) )), repeat / * step size for (32b) * / v (n+1) = P C (θ θ θ (n) + α θ ∇f ξ (θ θ θ (n) )), Algorithm 2: Backtracking line search for finding a step size for (32a) and (32b) step size, and then decrease it until a better feasible solution is found. As we can see from Algorithm 2, the algorithm will always terminate with a better point, in the sense of maximizing the objective function f ξ (θ θ θ). Note that the backtracking line search in Algorithm 2 follows the Barzilai-Borwein (BB) rule [22]. As ∇f ξ (θ θ θ) is Lipschitz continuous with a Lipschitz constant L f given in (84), the line search procedure is guaranteed to terminate after finite steps.

D. Convergence Analysis of Proposed Method
The convergence of Algorithm 1 is guaranteed that of the APG method and the penalty method. Specifically, for a given ξ, similar to Theorem 1 of [20], we can show that the objective sequence {f ξ (θ (n) )} is monotonically increasing. Also, the sequence {θ (n) } is bounded and thus has accumulation points. Each accumulation point is also a stationary solution to (31).
Furthermore, following the arguments in [19,Chap. 10] we can show that the iterate sequence {θ ξm } converges (in the subsequence sense) to a feasible point of (P 2 ) when ξ m → ∞. Thus the obtained solution of Algorithm 1 is also a stationary point of (P 2 ). The proof of these claims is given in Appendix E. We note however that since Algorithm 1 will terminate for some finite ξ m when a pre-determined error tolerance is met, it can only produce an approximate stationary solution of (P 2 ).

E. Computational Complexity Analysis
It is obvious that the complexity of Algorithm 1 in each iteration is dominated by that of (32a), and (32b) and the computation of the objective. Here we use the big-O notation to analyse the complexity of (32a) and (32b). From Proposition 2, it is easy to see that the complexity to calculate the gradient of ∇f ξ (θ θ θ) is O(KM 2 ). The projection P C (u) requires the complexity of O(KM) which is obvious from Proposition 3. Similarly, the complexity of computing f ξ (θ θ θ) is O(KM 2 ). As a result, the overall complexity of Algorithm 1 is O(I P I AP G KM 2 ) where I P and I AP G are the number of iterations of the outer loop (i.e. the penalty method) and the inner loop (i.e. the APG method) in Algorithm 1, respectively. Note that the line search procedure contributes negligible complexity since the gradient can be reused and the projection requires much less complexity. In [13], a SCA method based on solving a sequence of SOCPs was presented. We remark that the complexity of solving an SOCP in each SCA iteration is O( [23]. Thus the complexity of the SCA method in [13] is O( the number of SCA iterations. It is apparent that the computational complexity of our proposed method is much lower than the SCA method in [13].This point is numerically demonstrated in the next section.

V. NUMERICAL RESULTS
In this section, numerical results will be provided to evaluate as well as show the benefits of our proposed algorithm.

A. System Setup
We consider cell-free massive MIMO systems, where locations of M APs and K users randomly uniformly generated within an area of 1 × 1 km 2 . The wrapped around technique is used. The large-scale fading coefficient is modeled as: where z mk is the log-normal shadowing with the standard derivation σ sh = 8 dB, and PL mk is the three-slope-based path loss, which is modeled (in dB) as

B. Convergence of Proposed Algorithm
In the first experiment, we show the performance of the proposed APG method in comparison with the sequential SOCPs-based method in [13], with D = 1 km, τ c = 200, and τ p = K. As can be seen in Fig. 2, the proposed method achieves the same EE performance as the SOCPs- Next, we need to verify that the proposed APG algorithm will not violate any PFs, or equivalently, the total loss, K k Ψ k (θ θ θ), will converge to 0, regardless the starting point. The numerical results are shown in Fig. 4 using two different scenarios with the number of APs, M = 100 and M = 400, respectively. In this figure, blue curves and orange curves represent total PFs (total loss) and total EE, respectively. In the both scenarios, total PFs starts in infeasible domain, and gradually converges to 0 when the algorithm terminate, as the result of increasing the penalty parameter ξ.      The result shows the trade-off between the spectral efficiency and energy efficiency. While the total SE proportionally increases when extending the system scale (because the number of users increases), the total EE shows the opposite trend (because the energy consumption increases).

D. System setup based on number of users
Depending on practical applications and requirements, suitable numbers of APs and users can be designed.

VI. CONCLUSION
We applied a APG method to deal with a large-scale EE optimization problem in cellfree massive MIMO systems, taking into account of power contraints at the APs and QoS contraints at each user. To compare with the sequential SOCPs-based approaches, our proposed method achieves the same performance, while its run time is much faster, i.e., one to two orders of magnitude reduction. Therefore, this method has a great potential to deal with largescale optimization problems in cell-free massive MIMO, and hence, can be applied to practical scenarios. In addition, our optimization problem is done over large-scale fading time scale. It means the optimal power control coefficients can be updated only once for each large-scale fading realization. Since the large-scale fading coefficient changes very slowly with time, the system can fully operate in a bursty communication scenario with some random activation of the users.

A. Useful properties and Lemmas
In this section, we provide some useful properties and lemmas related to Lipschitz continuity that shall be used to analyse the Lipschitz continuity of the gradient of the objective in (31).
Some of the following properties and lemmas are in fact extension of the results in [24, Section 1.5] for scalar-valued functions. We also provide proofs to these for the sake of completeness.

1) Linear
Combinations: Let f 1 : R n → R n and f 2 : R n → R n be Lipschitz continuous with a Lipschitz constants L 1 and L 2 , respectively. Then f 1 ± f 2 is Lipschitz continuous with a Lipschitz constant L 1 + L 2 . Proof: where (a) is the triangle inequality. Proof: 3) The Composition of Functions: Let f 1 : R → R and f 2 : R n → R be Lipschitz continuous with Lipschitz constants L 1 and L 2 , respectively. Then the composite function f 1 •f 2 is Lipschitz continuous with a Lipschitz constant L 1 L 2 . Proof:

4) The Quotient of Functions:
Let f 1 : R n → R n and f 2 : R n → R be Lipschitz continuous with Lipschitz constants L 1 and L 2 , respectively on a bounded set S, such that ||f 1 (x)|| ≤ M for all x ∈ S, and further assume that there is a constant c > 0 such that |f 2 (x)| ≥ c for all x ∈ S. Then f 1 /f 2 is Lipschitz continuous with a Lipschitz constant (ML 1 + L 2 /c 3 ).
Proof: First, since f 2 is Lipschitz continuous with a Lipschitz constants L 2 , and |f 2 (x)| ≥ c, we have Thus 1/f 2 is Lipschitz continuous with a Lipschitz constant L 2 c 2 . Next, applying the product property in A2 for f 1 and 1/f 2 , we have Lemma 1. Let f : R n → R, and assume that ||x + y|| ≤ c, ∀x, y ∈ S. Then f (x) = (a T x) 2 is Lipschitz continuous with a Lipschitz constant c||a|| 2 . Proof: where (b1) and (b2) base on the Cauchy-Schwarz inequality.
As ∇f ξ (θ θ θ) is computed from u k (θ θ θ), v(θ θ θ), g k (θ θ θ), ∇u k (θ θ θ), ∇v(θ θ θ), and ∇g k (θ θ θ) we now need to find the Lipschitz constants of these terms, and then apply properties in Appendix A to conclude the Lipschitz constant of ∇f ξ (θ θ θ). To this end the following results are in order 1) ∇u k (θ θ θ) is Lipschitz continuous, and ||∇u k (θ θ θ)|| is bounded from above: Recall that where and In order to show that ∇u k (θ) is Lipschitz continuous, we will prove that the first term and the second term of the right hand side in (45) are Lipschitz continuous, respectively. First, note that by applying Lemma 1 and Lemma 2, it is easy to see that n k (θ θ θ), d k (θ θ θ) and n k (θ θ θ) + d k (θ θ θ) are Lipschitz continuous and their Lipschitz constants are and respectively. Next, we have ∇n k (θ θ θ) = 2ρ d N 2 A T kγ γ γ kkγ γ γ T kk A k θ θ θ C n k θ θ θ is Lipschitz continuous, as where λ max (C n k )is the largest eigenvalue of C n k . Moreover, and ln2 (n k (θ θ θ) + d k (θ θ θ)) ≥ ln 2.
Then, from (52), (53), (54), and (55), by applying the quotient of functions in Section A4, we have that ln2(n k (θ θ θ)+d k (θ θ θ)) is Lipschitz continuous with a Lipschitz constant as follow By applying the similar method to ∇n k (θ θ θ), we can prove that ∇d k (θ θ θ) is Lipschitz continuous with a Lipschitz constant written as and where (d1) is based on the fact that |n k (θ θ θ)| is bounded by the power constraints. Then, from (50), (57), (58), and (59), by applying the product of functions in Section A2, it is easy to see that n k (θ θ θ)∇d k (θ θ θ) is Lipschitz continuous with a Lipschitz constant given by and Similar to (59), we have for some d max , and thus Then, from (51), (52), (62), (63) and by applying the product of functions in Section A2, we can see that (n k (θ θ θ) + d k (θ θ θ)) d k (θ θ θ) is Lipschitz continuous with a Lipschitz constant given by Next, we have We now can conclude, from (60), (61), (64), and (65) and by applying the quotient of functions in the Section A4, that n k (θ θ θ)∇d k (θ θ θ) ln2(n k (θ θ θ)+d k (θ θ θ))d k (θ θ θ) is Lipschitz continuous with the following Lipschitz constant Finally, from (56), and (66), by applying the linear combination of the functions in Section A1, it follows that ∇u k (θ θ θ) is Lipschitz continuous with a Lipschitz constant given by and ||∇u k (θ θ θ)|| is bounded as 2) u k (θ θ θ) is bounded and Lipschitz continuous: First, we have Next, to prove u k (θ θ θ) is Lipschitz continuous, we first consider the function h 1 (x) = log 2 (x) over the domain 1 ≤ x ≤ x max . Note that h 1 (x) is continuously differentiable and thus we have By the mean value theorem, there exists some ξ between y and z, y, z ∈ (1, x max ), such that In other words, h 1 (x) = log 2 (x) is Lipschitz continuous with a constant log 2 e. Next we consider the function h 2 (x) = 1 + n k (θ θ θ) d k (θ θ θ) , for which we have where (d3) is based on the fact that |d k (θ θ θ)| is bounded by the power constraints. Then, from (50), (51), (59), and (72), the quotient property in Section A4 implies that h 2 (x) is Lipschitz continuous with a Lipschitz constant found as Finally, from (71), and (73) by applying composition property in Section A3, we can prove that u k (θ θ θ) is Lipschitz continuous with the following Lipschitz constant 3) v(θ θ θ) and ∇v(θ θ θ) are bounded and Lipschitz continuous: It is easy to see that v(θ θ θ) is Next, let us rewrite v(θ θ θ) as Then, following the same steps in Section B1, we can show that v(θ θ θ) is Lipschitz continuous with a Lipschitz constant L v expressed as Next, we have and thus ||∇v(θ θ θ)|| is bounded by and ∇v(θ θ θ) is Lipschitz continuous with the following Lipschitz constant 4) g k (θ θ θ) and ∇g k (θ θ θ) are bounded and Lipschitz continuous: By following the same method in the Section B1 and using the fact that d k (θ θ θ) ≥ 1, we can show that g k (θ θ θ) is Lipschitz continuous with L g k = a k L d k + ||γ γ γ T kk A k || and bounded as ∇g k (θ θ θ) is also Lipschitz continuous with L ∇g k = a k L ∇d k and bound as which is due to the fact that the second term in (24) is always greater than or equal to zero.
Then, following the same method to find the Lipschitz constant of ∇u k (θ) in Section B1, we can prove that ∇f ξ (θ θ θ) is Lipschitz continuous with a Lipschitz constant given by where ; and L 3 = 2ξ K k=1 (ζ g k L g k + ζ ∇g k L ∇g k ).

C. Proof of Proposition 2
Recall that the PF is then which can be rewritten as Using the quotient rule we can write ∇f (θ θ θ) as where and To find the gradient of ∇g k (θ θ θ) we recall the following equalities where B k ∈ R M ×M + is a diagonal matrix whose m-th element is [B k ] m = β mk . Then, by applying the chain rule, we can easily compute the gradient of ∇g k (θ θ θ) as shown in (37). To find the gradient of ∇u k (θ θ θ), we first apply the chain rule together with the quotient rule, we have where and The substitution of (96) and (97) into (93) yields (38).

D. Proof of Proposition 3
We now show that the projection onto C admits an analytical solution and is parallelizable.
Recall that P C (u) is explicitly written as Note that the objective in (98)

E. Convergence Proof of Algorithm 1
The convergence proof is divided into two parts. In the first part of the proof we show that, for a given ξ m , the APG iterations converge to a stationary solution to the penalized problem (31).
In the second part we show that {θ θ θ ξ } converges to a feasible point of (P 2 ). Thus the convergent point of Algorithm 1 is indeed a stationary point of (P 2 ).
In the second part of the proof, we now assert that any accumulation point is a stationary solution of (31). As f ξ (θ θ θ (n) ) is non-decreasing, the objective at all the accumulation points is the same which is denoted by f * ξ . Then, from (103) we have Summing over n = 1, 2, · · · , ∞, we have Since α x < 1 L , we can conclude that v (n+1) → θ θ θ (n) as n → ∞.
This inequality simply means that θ θ θ * is a stationary solution to (31) which completes the first part of the proof. Now we show that {θ θ θ ξm } indeed converges to a feasible point of (P 2 ). Note that for small ξ m , θ θ θ ξm may not be feasible to (P 2 ), and that the following inequalities always hold for lager where f * is the optimal objective of (P 2 ). In the above, the first equality is due to the negativity of the penalty term and the second inequality is true for two reasons. First, the problem max θ∈C f ξm (θ) becomes a convex problem for large ξ m since f ξm (θ) becomes concave. Second, the APG method can find the optimal solution since the problem is now convex. Thus the second inequality in (111) holds because the optimal objective is no less than the objective at any feasible solution.
Thus for these k we would have Thus it is easy to see that f ξm (θ ξm ) → −∞ when ξ m → ∞, which contradicts (111) and thus completes the proof.