{'fraction_non_alphanumeric': 0.04569325373670272, 'fraction_numerical': 0.012612774361506351, 'mean_word_length': 4.067091198544462, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accelerated Projected Gradient Method for the Optimization of Cell-Free Massive MIMO Downlink

We consider the downlink of a cell-free massive multiple-input multiple-output (MIMO) system where large number of access points (APs) simultaneously serve a group of users. Two fundamental problems are of interest, namely (i) to maximize the total spectral efficiency (SE), and (ii) to maximize the minimum SE of all users. As the considered problems are non-convex, existing solutions rely on successive convex approximation to find a sub-optimal solution. The known methods use off-the-shelf convex solvers, which basically implement an interior-point algorithm, to solve the derived convex problems. The main issue of such methods is that their complexity does not scale favorably with the problem size, limiting previous studies to cell-free massive MIMO of moderate scales. Thus the potential of cell-free massive MIMO has not been fully understood. To address this issue, we propose an accelerated projected gradient method to solve the considered problems. Particularly, the proposed solution is found in closed-form expressions and only requires the first order information of the objective, rather than the Hessian matrix as in known solutions, and thus is much more memory efficient. Numerical results demonstrate that our proposed solution achieves far less run-time, compared to other second-order methods.


I. INTRODUCTION
Cell-free massive multiple-input multiple-output (MIMO) was introduced in [1] to overcome the inter-cell interference which is the main inherent limitation of cellular-based networks. In cell-free massive MIMO, many access points (APs) distributed over the whole network serve many users in the same time-frequency resource. There are no cells, and hence, no boundary effects. Unlike colocated massive MIMO, each AP is equipped with just a few antennas. But when the number of APs is very large, cell-free massive MIMO is still able to exploit the favorable propagation and channel hardening properties, like the colocated massive MIMO. In particular, in the downlink, each AP uses its local channel estimates acquired during the uplink training phase to perform simple beamforming techniques. With this way, there is no need for exchanging the instantaneous channel state information (CSI) among the APs or the central prcoessing unit (CPU) [2].
The research on cell-free massive MIMO is still in its infancy and thus deserves more extensive studies. We discuss here some of the noticeable and related studies in the literature. In [1], Ngo et al. considered the problem of minimum rate maximization to provide uniformly good services to all users. The problem was then solved using a bisection search and a sequence of linear feasibility problems. In [3], Nguyen et al. adopted zeroforcing precoding and studied the energy efficiency maximization (EEmax) problem. In this work, an iterative method based on successive convex approximation (SCA) was derived. In [4], both the max-min fairness and sum-rate maximization problems were considered and solved by SCA. The SCA-based method was also used in [2] and [5] to solve the EEmax and max-min fairness power controls with different cell-free massive MIMO setups, respectively.
A common feature of all the above mentioned pioneer studies on cell-free massive MIMO is the use of a second-order interiorpoint method which requires the computation of the Hessian matrix of the objective, and thus their computational complexity and memory requirement makes them infeasible in practical power control of large-scale cell-free massive MIMO. It only allows us to characterize the performance of cell-free massive MIMO for a relatively small area. For example, the work of [1] was able to consider an area of 1 km×1 km, consisting of 100 APs serving 40 users. Numbers with the same order of magnitude were also observed in the above mentioned papers. The performance of these scenarios fractionally represents the full potential of cell-free massive MIMO.
To fully understand the performance limits of cell-free massive MIMO, we need to devise more scalable resource allocation methods. To this end, we propose in this paper a first order method to maximize the total spectral efficiency and the minimum rate of the downlink where the conjugate beamforming is adopted. A similar method has been used in [6] to solve the EEmax problem. In proposed method, we customize the accelerated proximal gradient (APG) method presented in [7], which mainly requires the computation of the gradient of the objective. The proposed method is provably convergent and is numerically shown to achieve much lower run-time compared to an existing second-order method given in [2].
Notations: Bold lower and upper case letters represent vectors and matrices. CN (0, a) denotes a complex Gaussian random variable with zero mean and variance a. X T and X † stand for the transpose and Hermitian of X, respectively. x i is the i-th entry of vector x; [X] i,j is the (i, j)-th entry of X. ∇f (x) represents the gradient of f (x) and ∂ ∂xi f (x) is the partial gradient with respect to x i . [x] + denotes the projector onto the positive orthant. || · || represents the Euclidean norm; | · | is the absolute value of the argument.

A. System Model
We consider the downlink of a cell-free massive MIMO system model as in [2]. In particular, there are M APs serving K single-antenna users in time division duplex (TDD) mode. Each AP is equipped with N antennas. All the APs and the users are assumed to be distributed in a large area. As TDD operation is adopted, APs first estimate the channels using pilot sequences from the uplink (commonly known as uplink training) and then apply a beamforming technique to transmit signals to all users in the downlink, or use a matched filter technique to detect signals in the uplink. Since this work focuses on the downlink transmission, we neglect the uplink phase. Let us denote by T c and T p the length of the coherence time and the uplink training phase in samples. The uplink training and downlink payload transmission phases are summarized as follows. The details can be found in [2].

1) Uplink Training and Downlink Payload Transmission:
We assume the channel is reciprocal, i.e., the channel gains on the uplink and on the downlink are the same. Consequently, APs can estimate the downlink channel based on the pilot sequences sent by all users on the uplink. Let T p ψ k ∈ C Tp×1 , where ||ψ k || 2 = 1, be the pilot sequence transmitted from the k-th user, k = 1, . . . , K. As in [2], we model the channel between the m-th AP and the k-th user as where β mk represents the large-scale fading and h mk ∈ C N ×1 comprises of small-scale fading coefficients between the N antennas of the m-th AP and the k-th user. We further assume that the entries of h mk follow i.i.d. CN (0, 1).
Let r mk = R up,m ψ k , where R up,m is the received signal at the m-th AP. Given r mk , the minimum mean-square error (MMSE) estimate of g mk is calculated as [2] g mk = ζ p T p β mk where ζ p is the normalized transmit signal-to-noise ratio (SNR) of each pilot symbol. The mean square of any element ofĝ mk is given by For downlink payload data transmission, the APs use the channel estimates (1) and conjugate beamforming to form separate radio beams to the K users. Denote the symbol to be sent to the k-th user by c k and the power control coefficient between the m-th AP and the k-th user by η mk . Then the transmitted signal at the m-th AP is where ζ d is the maximum downlink transmit power at each AP normalized to the noise power. Note that the total power at each AP is The received signal at the k-th user is written as where w k ∼ CN (0, 1) is the white Gaussian noise.
2) Spectral Efficiency: By using the use-and-then-forget capacity bounding technique, we can obtain the following spectral efficiency of the k-th user [2] where and

B. Problem Formulation
To formulate the considered problem and to facilitate the development of the proposed algorithm, we define µ m ∈ R K + to be the vector of all power control coefficients associated with the m-th AP as µ m to include the power control coefficients of all APs. To express the spectral efficiency in (6) as a function of µ, we denote byμ k = [µ 1k ; µ 2k ; . . . ; µ Mk ] the vector of power control coefficients associated with user k. Thus we can write ν T ik η i asν ikμk , wherē Similarly, we can write D ikηi asD iμi , whereD i is a diagonal matrix with the m-th diagonal entry equal to √ β mi . Now the spectral efficiency of the k-th user (in nat/s/Hz) can be expressed as where γ k (µ) is the SINR of the k-th user given by The total spectral efficiency of the system is defined as In this paper, we consider a total power constraint at each AP which is given by E{||x m || 2 } ≤ ζ d or ||µ m || 2 ≤ 1 N , m = 1, 2, . . . , M , which follows from (4). For the problem formulation purpose, we define the following set S = µ|µ ≥ 0; ||µ m || 2 ≤ 1/N, m = 1, 2, . . . , M . In this paper, we consider the following two common power control optimization problems, namely • The problem of maximizing the total spectral efficiency (SEmax) • The problem of maximizing the minimum rate (MRmax) among all users (also known as max-min fairness maximization) For the above problems, SCA has gradually become a standard mathematical tool [1], [2]. In all known solutions for the considered problems or related ones, interior point methods (through the use of off-the-shelf convex solvers) are invoked to solve these convex problems [2], [4], [5], which do not scale with the problem size. In the next section, we propose methods that can tackle this scalability problem.

III. PROPOSED SOLUTIONS
In this section, we present solutions to (P 1 ) and (P 2 ), using the accelerated proximal gradient (APG) methods introduced in [7], which concerns the following problem where f (µ) is differentiable (but possibly nonconvex) and g(x) can be both nonconvex and nonsmooth. If we let g(µ) be the indicator function of S, defined as then (13) is actually equivalent to (P 1 ). Basically the objective of (P 2 ) is nonsmooth, but we can still solve (P 2 ) by applying a proper smoothing technique. We also note that when g(µ) is the indicator function of S, the proximal operator of g(µ) becomes the Euclidean projection onto S. In the following, we customize the APG methods to solve the considered problems.
A. Proposed Solution to (P 1 ) Since f (µ) for (P 1 ) is differentiable, the proposed algorithm for solving (P 1 ) follows closely the monotone APG method in [7], which is outlined in Algorithm 1. In Algorithm 1, α > 0 is called the step size which should be sufficiently small to guarantee its convergence. Also, the notation P S (u) denotes the projection onto S, i.e., P S (u) = arg min ||x − u|| | x ∈ S . From a given operating point, we move along the direction of the gradient of f (µ) with the step size α, and then project the resulting point onto the feasible set. In particular, y n in Step 4 is an extrapolated point which is used for convergence acceleration. However, y n can be a bad extrapolation unlike APG methods for the convex case, and thus Step 7 is there to fix this issue. We now give the details for the two main operations of Algorithm 1, namely: the projection onto the feasible set S and the gradient of f (µ).
Algorithm 1 General Description of Proposed Algorithm for Solving (P 1 ) and (P 2 ) 1: Input: µ 0 >= 0, t0 = t1 = 1, α > 0 2: µ 1 = z 1 = µ 0 3: for n = 1, 2, . . . do 4: tn+1 = 0.5 √ 4t 2 n + 1 + 1 9: end for 10: Output: µ * 1) Projection onto S: We show that the projection in Steps 5 and 6 in Algorithm 1 can be done in parallel and by closedform expressions. Recall that for a given x ∈ R MK×1 , P S (x) is the solution to the following problem It is easy to see that the above problem can be decomposed into sub-problems at each AP m as The above problem is in fact the projection onto the intersection of a ball and the positive orthant. Interestingly, the analytical solution to this problem can be found by applying [8, Theorem 7.1], which produces 2) Gradient of f (µ) for (P 1 ): To implement Algorithm 1, we also need to compute ∇ µ f (µ), which is found as where The gradient of SE k (µ) with respect toμ i , i = 1, 2, . . . , K, is found as Now we recall the following equality ∇||Ax|| 2 = 2A T Ax for any symmetric matrix A, and thus ∇μ i b k (µ) and ∇μ i c k (µ) are respectively given by (22)

B. Proposed Solution to (P 2 )
We recall that for (P 2 ) the objective function is which is non-differentiable. Thus a straightforward application of the APG method is impossible. To overcome this issue, we adopt a smoothing technique. In particular, f (µ) is approximated by the following log-sum-exp function given by [9] f τ (µ) = − 1 τ where τ > 0 is the positive smoothness parameter. To obtain (24), we have used the fact that min In other words, f τ (µ) is a differentiable approximation of f (µ) with a numerical accuracy of log K τ . Thus, with a sufficiently high τ , we can find an approximate solution to (P 2 ) by running Algorithm 1 with f (µ) being replaced by f τ (µ) in (24). In this regard, the gradient of f τ (µ) is easily found

A. Complexity Analysis
We now provide the complexity analysis of the proposed algorithm for one iteration using the big-O notation. It is clear that the complexity of Algorithm 1 is dominated by the computation of three quantities: the objective, the gradient, and the projection. It is easy to see that KM multiplications are required to compute SE k (µ). Therefore, the complexity of finding Similarly, we can find that the complexity of ∇SE(µ) which is O(K 2 M ). The projection of µ onto S is given in (17), which requires the computation of the l 2 -norm of K × 1 vector x m at each AP, and thus the complexity of the projection is O(KM ). In summary, the per-iteration complexity of the proposed algorithm for solving (P 1 ) is O(K 2 M ). We can also conclude that the per-iteration complexity for solving (P 2 ) is also O(K 2 M ).

B. Convergence Analysis
We now discuss the convergence result of Algorithm 1 for solving (P 1 ) and (P 2 ). For (P 1 ), from (20), (21), (22), we can check that ∇f (µ) is Lipschitz continuous, or equivalently f (µ) has Lipschitz continuous gradient. That is, there exists a constant L > 0 such that For the considered problems, it is possible to find a Lipschitz constant for ∇f (µ) but we skip the details here for the sake of brevity. Let us simply denote by L the Lipschitz constant of ∇f (µ). The convergence of Algorithm 1 is stated in the following lemma.
Lemma 1. Let the step size α in Algorithm 1 satisfy α < 1 L . Then the iterates {µ n } generated by Algorithm 1 are bounded. Let µ * be any accumulation point of {µ n }, then µ * is a critical point of (P 1 ).
Proof: Please see Appendix A. The same convergence result applies to Algorithm 1 for solving (P 2 ) with the differential approximation in (24).

C. Improved Convergence with Line Search
In practice, we do not need to compute a Lipschitz constant of ∇f (µ) for two reasons. First, the best Lipschitz constant of ∇f (µ) (i.e. the smallest L such that (26) holds) is hard to find. Second, the conditions α < 1 L is sufficient but not necessary for Algorithm 1 to converge. Thus, we can allow α to take on larger values to speed up the convergence of Algorithm 1 by means of a linear search procedure. Specifically, let L 0 > 0, δ > 0 and ρ ∈ (0, 1). Then the step size α in Step 5 of Algorithm 1 is set to α = L o ρ mn , where m n is the first nonnegative smallest integer such that The above line search follows the Armijo rule. The same line search procedure can be used for Step 6 of Algorithm 1. We refer to these modifications as Algorithm 1 with line search.

V. NUMERICAL RESULTS
In this section, we evaluate the performance of the proposed method in terms of computational complexity and achieved spectral efficiency. The users and the APs are uniformly dropped over a D × D km 2 . The large-scale fading coefficient between the m-th AP and the k-th user is generated as where PL mk and z mk represent the path loss and log-normal shadowing with mean zero and standard deviation σ sh , respectively. In this paper we adopt the three-slope path loss model as in [2] in which PL mk (in dB) is given by where L is a constant dependent on carrier frequency, d mk (in km) is the distance between the m-th AP and the k-th user, and d 0 and d 1 (both in km) are reference distances.
Similar to [2], we choose L = 140.7 dB, d 0 = 0.01 km and d 1 = 0.05 km. We consider a system having a bandwidth of B = 20 MHz, the noise power density is N 0 = −174 (dBm/Hz), and a noise figure of 9 dB. The length of the coherence time and the uplink training phase are set to T p = 20, T c = 200, respectively. If not otherwise mentioned, we set the power transmit power for downlink data transmission and uplink training phase (before normalization) as ζ d = 1 W and ζ p = 0.2 W. Single antenna AP is considered in all numerical experiments. These parameters are taken from [2].
For comparison purpose, we modify the SCA-based method presented in [2], which is dedicated to the EEmax problem. However, it can be easily modified to deal with the problems considered in this paper. In this regard we use convex conic solver MOSEK [10] through the modeling tool YALMIP [11].
In the first numerical experiment, we compare the convergence rate of the proposed method with the SCA-based method. Figures 1(a) and 1(b) show the convergence of the proposed method and the SCA-based method for the total spectral efficiency and the min-rate maximization problem, respectively. We can see that the proposed and the SCA-based methods achieve the same performance but the SCA-based method requires fewer iterations. However, the main advantage of our proposed method over the SCA-based method is that each iteration of the proposed method is very memory efficient, and hence, performs very fast. As a result, the run-time of the proposed method is far less than that of the SCA-based method as shown in Table I. In Table I, we report the actual run-time of both methods to solve the SEmax problem. Here, we execute our codes on a 64-bit Windows operating system with 16 GB RAM and Intel CORE i7, 3.7 GHz. Both iterative methods are terminated when the difference of the objective for the last 5 iterations is less than 10 −3 . We next take advantage of the proposed method to explore the spectral efficiency performance of cell-free massive MIMO for a large metropolitan area. In particular, we investigate the performance for two cases D = 1 and D = 10. To obtain a fair comparison we keep the AP density, defined as the number of APs per square kilometer, the same for both cases. Note that the AP density of 1000 means 10 000 APs for the case of D = 10, which has not been studied in the literature previously.
To appreciate the proposed method for this large-scale scenario, we compare it with the SCA method and the equal power allocation (EPA) method where the power control coefficient η mk is given by, η mk = ( k i=1 ν mi ) −1 . The results in Figure 2 are interesting. First, increasing the AP density improves the sum spectral efficiency of the system. Second, for the same AP density, a larger area provides a better sum spectral efficiency. The reason is that for a larger area, the users that are served by the APs become far apart each other. As a result, the inter-user inference becomes weaker, leading to an improved sum spectral efficiency. On the other hand, the EPA method yields smaller spectral efficiency as the coverage area is larger because more power should be to spent to the users having small path loss. The SCA method produces the same spectral efficiency as the proposed APG method but it is unable to run for D = 10 on the system specifications mentioned above. Thus, the proposed scheme outclasses both SCA-based and EPA methods in terms of improved sum spectral efficiency and larger coverage area.  Fig. 3. CDF of per-user spectral efficiency for (P 1 ) and (P 2 ). The parameters other than the mentioned are kept the same. expected, the spread of CDF for the SEmax power control is the larger than the max-min fairness power control. On the other hand, the order is reversed in terms of fairness as the CDF curve of the max-min problem has the largest slope which shows that the spectral efficiencies of all the users are very close to each other. The per-user spectral efficiencies in the large-scale scenario are lesser than the ones in the small-scale scenario. This is because when the number users increases, the inter-user interference increases which yields to the reduced performance.

VI. CONCLUSION
We have considered the downlink of cell-free massive MIMO and aimed to maximize the total and minimum spectral efficiencies, subject to a sum power constraint at each AP. Conjugate beamforming has been adopted, resulting in a power control problem for which an accelerated project gradient method has been proposed. The proposed solution only requires the first order information of the objective and, in particular, can be founded by closed-form expressions. We have numerically shown that the proposed method can achieve the same spectral efficiency as a known SCA-based method but with much lesser run-time. For the first time, we have evaluated the SE performance of cell-free massive MIMO for an area of 10 km×10 km, consisting of up to 10 000 APs. In this case, the achieved sum spectral efficiency can be up to 200 (bit/s/Hz).
Applying (28) yields It is easy to see that f (v n+1 ) ≥ f µ n if α < 1 L . From Step Similarly, if f (z n+1 ) < f (v n+1 ), then From (31), (32), and (33) we have Since the feasible set of the considered problems is compact convex, the iterates {v n } and {µ n } are both bounded and thus, {µ n } has accumulation points. We can prove that each accumulation point of Algorithm 1 is indeed a critical point of (P 1 ), following the same arguments as those in [7] but skip the details due to the space limitation.
ACKNOWLEDGMENT This publication has emanated from research supported in part by a grant from Science Foundation Ireland under grant number 17/CDA/4786.