{'fraction_non_alphanumeric': 0.050366636113657194, 'fraction_numerical': 0.012144821264894592, 'mean_word_length': 4.077012563983248, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mirror Prox Algorithm for Large-Scale Cell-Free Massive MIMO Uplink Power Control

We consider the problem of max-min fairness for uplink cell-free massive multiple-input multiple-output (MIMO) subject to per-user power constraints. The standard framework for solving the considered problem is to separately solve two subproblems: the receiver filter coefficient design and the power control problem. While the former has a closed-form solution, the latter has been solved using either second-order methods of high computational complexity or a first-order method that provides an approximate solution. To deal with these drawbacks of the existing methods, we propose a mirror prox based method for the power control problem by equivalently reformulating it as a convex-concave problem and applying the mirror prox algorithm to find a saddle point. The simulation results establish the optimality of the proposed solution and demonstrate that it is more efficient than the known methods. We also conclude that for large-scale cell-free massive MIMO, joint optimization of linear receive combining and power control provides significantly better user fairness than the power control only scheme in which receiver coefficients are fixed to unity.


I. INTRODUCTION
The recently evolved form of the massive multiple-input multiple-output (MIMO) for the beyond fifth generation (5G) networks is cell-free massive MIMO [1], in which one or several central processing units (CPUs) control a large number of low-cost and low-power access points (APs) to serve a large number of users. The APs and users are distributed in a large coverage area. Owing to the high array gain, multiplexing gain, and macro-diversity gain, cell-free massive MIMO offers many advantages such as high energy efficiency, uniform quality of service, and flexible and cost-effective deployment [2].
To achieve the above benefits, suitable power allocation algorithms need to be employed to control the near/far effect. Since the numbers of APs and users are large, the resulting power control problems consist of many variables. Thus, simple and low-complexity solutions to power control optimization problems for cell-free massive MIMO are of particular interest. In [1], the problem of maximization of the minimum uplink rate (a.k.a. max-min fairness) for singleantenna APs and single-antenna users was solved using a bisection algorithm in combination with linear programming (LP), where maximum-ratio combining (MRC) was considered at each AP for signal detection. In [3], Bashar  subproblems: (i) the receiver coefficient design which is in fact a generalized eigenvalue problem (EVP), and (ii) the power control problem which was solved using geometric programming (GP). A similar alternating optimization (AO)based scheme was used in [4] to solve a mixed quality of service (QoS) problem including the max-min fairness for a set of users and a fixed QoS for the remaining users. Again, the power control subproblem was solved using GP. In [5], Mai et al. employed the zero-forcing (ZF) technique to detect symbols from multi-antenna users where the power control problem was solved using the bisection algorithm. A common feature of all aforementioned studies is the use of LP or (more complex) GP by means of off-the-shelf convex solvers, and thus are suitable for small-scale scenarios. In particular, a low complexity first-order accelerated projected gradient (APG) algorithm was applied to the max-min fairness problem following Nesterov's smoothing technique [6] in both downlink [7], [8] and uplink [9] channels. However, this method can only produce an approximate solution whose accuracy is inversely proportional to the number of users, and hence, is not suitable for very large numbers of users. This clearly calls for novel solutions to the power control problem, which are more computationally efficient than LP or GP methods and more accurate than the APG method.
In this paper, we consider the max-min fairness problem subject to the power constraint at each user. We utilize the MRC technique to detect users' signals using the channel estimates. Similar to the known approaches, we also split the considered problem into two subproblems which are alternately solved. In particular, for the power control subproblem, we propose a mirror prox (MP)-based algorithm. To achieve this, we first reformulate the nonconvex power control problem into an equivalent convex-concave program and then propose an MP method [10] to find a saddle point. The proposed method only requires the first order information of the objective, and thus, is shown to be much faster than the LP-based and GPbased methods. Moreover, it outperforms other existing firstorder methods in terms of achievable rate performance.
Notations: Bold lower and upper case letters represent vectors and matrices. CN (0, K) denotes the multivariate circularly symmetric complex Gaussian random distribution with zero mean and co-variance matrix K. X T and X H stand for the transpose and Hermitian of X, respectively. Notation e i denotes the i-th column of the identity matrix. ∇f (x) represents the gradient of f (x). · denotes the Euclidean or ℓ 2 -norm and | · | is the absolute value of the argument. P C (x) denotes the projection of x onto the convex set C, i.e. P C (x) = arg min z∈C x − z .

II. SYSTEM MODEL AND PROBLEM FORMULATION
Consider a cell-free massive MIMO uplink scenario where L single-antenna users are served by M APs with K antennas per AP. All APs and users are randomly distributed in a coverage area. The APs are connected to a CPU via highcapacity dedicated links. We denote by ζ ml the large-scale fading coefficient between the l-th user and the m-th AP, andh ml ∼ CN (0, I K ) is the vector of the small-scale fading coefficient for all antennas at the m-th AP. The channel vector between the l-th user and the m-th AP is modeled as (1) The length (in symbols) of the uplink training is denoted by τ p . Let η p be the transmit power of each pilot symbol and φ l be the pilot sequence of unit norm transmitted from user l to the APs. Letĥ ml be the minimum mean square error (MMSE) channel estimate of h ml . Thenĥ ml ∼ CN (0, g ml I K ), where [1] To formulate the problem of interest, we define p = [p 1 , . . . , p L ] as the vector of power control coefficients, q = [q 1 ; q 2 ; . . . ; q L ] as the receiver coefficient where q l = [q 1l ; q 2l ; . . . ; q Ml ] associates with the l-th user, andp as the maximum transmit power at each individual user. With this setup, an achievable signal-to-interference plus noise ratio (SINR) is given by [1] (3) In (3),Ḡ li ∈ R M×M + andG l ∈ R M×M + are diagonal matrices with [Ḡ li ] m,m = g ml ζ mi and [G l ] m,m = g ml , respectively, and Note that the achievable rate of the l-th user is given by R l (p, q) = log 2 (1 + γ l (p, q)) .
(5) Similar to [3], we consider the problem of maximization of the minimum rate, which is mathematically stated as maximize First, we note that (P 1 ) is nonconvex, and thus, finding a globally optimal solution is difficult and not practically useful. The standard framework for solving (P 1 ) which is adopted in previous studies such as [1], [3], [4] is based on AO. In this way, two subproblems arise: the receiver coefficient design and the power control problem. Specifically, the receiver coefficient design is obtained by fixing the power allocation p, which admits a closed-form solution given by [3], [4] where After updating the receiver coefficients for all users, we need to find the power coefficients, leading to the following power control problem: maximize To solve (P 2 ), existing methods [1], [3], [4] involve second order methods, i.e., GP or LP, which require very high complexity, and thus, is not suitable for large-scale cell-free massive MIMO. To overcome the complexity issue, a firstorder algorithm was introduced in [9] but it cannot yield a solution of high accuracy for a large number of users. In the next section, we propose a more efficient method for solving (P 2 ) based on the MP framework.
The proposed method is developed based on an equivalent convex-concave reformulation of (P 2 ). To this end, it is easy to see that we can equivalently rewrite (P 2 ) as minimize where We remark that γ −1 l (p) is non-convex but can be converted to a single convex form as for geometric programming by making a change of variables θ i = log p i [3], [4]. In this letter, we propose the change of variables θ i = log ωp i or p i = 1 ω e θi and reformulate (8) as where Θ {θ|θ ≤θ} andθ = log(ωp). Note that f l (θ) is convex and so is (P 3 ). In the context of projected gradient methods, the introduction of ω effectively scales the gradient of f l (θ). We shall numerically demonstrate that a proper value of ω can speed up the convergence of the proposed algorithm as shown in the next section.
Although (P 3 ) is convex, solving it efficiently is still challenging since its objective is nonsmooth due to the max operator. A solution to deal with the nonsmoothness of (P 3 ) is to adopt a smoothing technique as done in [9] but such a method can only produce an approximate solution. To derive an efficient solution, we recall the following equality max where λ = [λ 1 ; λ 2 ; . . . , λ L ] ∈ R L and ∆ λ | 1 T λ = 1; λ ≥ 0 is the standard simplex. Thus, (P 3 ) can be equivalently rewritten as minimize where φ(θ, λ) = L l=1 λ l f l (θ). Note that φ(θ, λ) is convex w.r.t θ for a given λ and concave (in fact linear) w.r.t λ for a given θ.

B. Proposed Mirror Prox Algorithm for Solving (P 4 )
It is now clear that we can find a saddle point of (P 4 ) to solve the power control problem, which indeed motivates the application of the MP algorithm in this letter. To explain the idea of the MP algorithm, let z = (θ, λ) and F(z) = [∇ θ φ(θ, λ), −∇ λ φ(θ, λ)] which is a monotone operator associated with (P 4 ), i.e., F(z)− F(z ′ ), z− z ′ ≥ 0. The idea of the MP algorithm is to find a point z * such that F(z), z − z * ≥ 0, which is indeed a saddle point of (P 4 ). The proposed method is a special case of the MP algorithm with the Euclidean setup. The motivation is that the Euclidean projections onto Θ and ∆ can be done in closed-form. Under the Euclidean setup, the distance between z = (θ, λ) and Note that D(z, z ′ ) is the sum of individual Euclidean distances.
The step-by-step description of the proposed method for solving (P 4 ) is as follows. Let z n = (θ n , λ n ) be the n-th iterate. To obtain the next iterate, we perform the following two proximal mappingŝ z n = arg min z∈Θ×∆ µ n F z n , z + D(z n , z), where µ n > 0 is a step size which needs to be chosen properly to guarantee the convergence. We remark that the two steps above are similar to a classical gradient-type method since the vector field F(z) acts as a descent direction. From the current iterate z n , we first move along the gradient of individual variables to obtain an intermediate pointẑ n = θ n ,λ n as shown in (12a). The main difference and in fact the novel idea of the MP algorithm are as follows. To obtain the next iterate from z n , we do not use the gradients at z n . Instead, z n+1 is obtained using the gradient at the intermediate point z n as given in (12b). Intuition. To gain further insights into the proposed MP algorithm, we note that (12a) implieŝ and, similarly,λ n = P ∆ λ n + µ n ∇ λ φ(θ n , λ n ) .
An important remark is in order. Since φ(θ, λ) is convex in θ and concave in λ, −∇ θ φ(θ, λ) and ∇ λ φ(θ, λ) are the descent and ascent directions, respectively. In this regard, the MP algorithm is similar to the gradient descent and ascent algorithm. That is, the MP algorithm simultaneously minimizes and maximizes φ(θ, λ) to reach a saddle point. However, the gradients at the intermediate point θ n ,λ n are used to obtain θ n+1 , λ n+1 from θ n , λ n .
To obtain a convergent algorithm, the step size in each iteration n needs to satisfy the following condition Note that the above inequality is true if µ n ≤ L −1 , where L is the Lipschitz constant of F(z). In practice, we find µ n by a back tracking line search. In summary, the proposed MP algorithm for solving (P 4 ) is outlined in Algorithm 1. Algorithm 1 MP Algorithm for Solving (P 4 ) 1: Initialization: µ 0 > 0, ρ ∈ (0, 1) 2: θ 1 ← θ initial , λ 1 ← λ initial 3: for n = 1, 2, . . . do 4: µ n = µ n−1 /ρ 5: until δ n ≤ 0 12: end for 13: Computation of gradients. To implement Algorithm 1, we need to calculate ∇ θ φ(θ, λ) and ∇ λ φ(θ, λ) which are given in closed form as where, from (9), we have (17) Remark 1. As mentioned above, the proposed change of variable θ i = log ωp i is equivalent to scaling the gradient of f l (θ), compared to the standard change of variable θ i = log p i . This scaling effect is reflected by the termc l in the above equation. The main motivation for introducing ω is to balance the gradients ∇ θ φ(θ, λ) and ∇ λ φ(θ, λ), which can accelerate the convergence of Algorithm 1 as numerically illustrated in the next section.
Projections onto Θ and ∆. It is easy to see that P Θ (x) admits a closed form solution as follows Also, the projection onto a standard simplex is given by where β is the solution to the equation which can be found by bisection. In summary, to solve (P 1 ) we keep alternately computing (6) and running Algorithm 1 until convergence.

C. Convergence Analysis
We first provide the convergence analysis of Algorithm 1. Let f * be the optimal objective of (P 3 ) andθ n = N n=1 µ n −1 N n=1 µ n θ n be the obtained solution after N iterations. Note thatθ n is the weighted average of the iterates θ n up to iteration N . Then it is shown that f (θ where Ω is a constant that depends on the distance generating function. By the line search procedure in Algorithm 1, we have that δ n ≤ 0 and µ n ≥ L −1 , and thus, f (θ n ) − f * ≤ ΩL/N . That is, Algorithm 1 can achieve a O(1/N )-rate of convergence. The proof follows closely the steps in [11,Proposition 6.1], and thus, is omitted here for brevity.

D. Complexity Analysis and Comparison
We now present the per-iteration complexity analysis of Algorithm 1. It is easy to see that L multiplications are required to compute f l (θ n ). Therefore, the complexity of finding the objective is O(L 2 ). Similarly, we can find that the complexity of ∇f l (θ) is also O(L 2 ). It is clear from (18) that P Θ (x) has complexity of O(L). For P ∆ (λ), by sorting the elements of λ in the ascending order, the complexity of the projection onto a simplex from [12,Theorem 2.2] is O(L) complexity in the worst case. In summary, the per-iteration complexity of the proposed algorithm for solving (P 2 ) is O(L 2 ). If a bisection search is used with LP as done in [1], the worst-case per-iteration complexity is O(L 3.5 ) [13, Eq. (8.1.6)]. The complexity of the GP-based methods in [3], [4] is O(L 4.5 ) [13, Sec. 6.3.1].

IV. NUMERICAL RESULTS
We evaluate the performance of our proposed method in terms of the achievable rate and the run time. We randomly distribute APs and users over a D × D km 2 . Channel coefficients are generated using (1), in which the large-scale fading coefficient between the m-th AP and the l-th user is modeled as ζ ml = PL ml z ml , where PL ml is the corresponding path loss, and z ml represents the log-normal shadowing between the m-th AP and the l-th user with mean zero and standard deviation σ sh , respectively. In this paper, we adopt the threeslope path loss model and model parameters as in [8]. Noise figure is set to 9 dB. We assume pilot sequences to be pairwisely orthogonal to avoid the effect of pilot contamination. The lengths of the coherence interval and the uplink training phase are set to τ c = 200, τ p = 20, respectively. If not otherwise mentioned, we set η p = 0.2 W andp = 0.2 W. The number of antennas at each AP is K = 1.
In Fig. 1, we show the convergence of Algorithm 1. For each channel realization, we run Algorithm 1 for a fixed q and plot the average achievable rate over 100 randomly generated channel realizations. We compare Algorithm 1 with two other iterative baseline schemes: the gradient descent accent (GDA) method [14], [15] and the APG method combined with the smoothing technique [9]. We also benchmark Algorithm 1 with the optimal solution obtained by either the bisection method with LP [1] or the GP in [3], [4].
We can see that Algorithm 1 based on the MP method has better objective than GDA and APG. Also, Algorithm 1 takes lesser number of iterations. Thus, GDA and APG [9] are not suitable for the applications where a fast convergence rate is required. Note that APG presented in [9] is essentially an approximate solution. On the other hand, Algorithm 1 is an exact method and thus can reach the optimal solution at convergence as clearly seen in Fig. 1. Note that we also demonstrate the impact of introducing ω when reformulating (P 2 ) into (P 4 ). As can be observed clearly, a proper value Iteration index Average achievable rate (bps/Hz) optimal solution [1], [3], [4] APG & GDA [8], [13], [14] mirror prox with ω = 1 mirror prox with ω = M of ω can indeed speed up the convergence of Algorithm 1 very significantly. By extensive simulation settings, we find that ω = M yields a good convergence rate for Algorithm 1 overall, which is shown in Fig. 1. Further, we remark that the proposed algorithm converges within 500 iterations. Also, the average run-time per iteration is very small, as shown in Fig.  1.
The main advantage of Algorithm 1 over second-order methods in [1], [3], [4] is that each iteration of Algorithm 1 is very memory efficient and computationally cheap, and hence, can be executed very fast. To demonstrate this, we compare the run-time of Algorithm 1 with the bisection method with LP in [1], both normalized to the run-time of the GP approach in [3], [4] for solving (P 2 ) in Fig. 2. This shows the factor at which a particular method is faster than the GP method with normalized run-time of 1. We run the codes on a 64bit Windows operating system with 16 GB RAM and Intel CORE i7, 3.7 GHz. All iterative methods are terminated when the difference of the objective for the latest two iterations is less than 10 −4 . We observe that solving (P 2 ) using a bisection search with LP is faster than converting it to a single GP, which is then solved by a dedicated GP solver. Most importantly, Algorithm 1 is 40 times faster than the LP method and approximately 100 times faster than the GP approach as shown in Fig. 2.
In the next experiment, we investigate cell-free massive MIMO for large-scale scenarios. Instead of solving a generalized eigenvalue problem to find q * l as suggested in [3], [4] which incurs a complexity of O(M 3 ), we can reduce this complexity to O(M 2 ) using the method in [9,Appendix]. This complexity reduction and the low complexity of Algorithm 1 indeed allow us to investigate the performance of uplink largescale cell-free massive MIMO, which has never been reported previously. In particular, we plot the cumulative distribution function (CDF) of the per-user achievable rate for 100 channel realizations in Fig. 3. Three large-scale scenarios detailed in the caption of Fig. 3   all three scenarios. In particular, we compare the CDF of the per-user rate where both receiver filtering and power control are considered (i.e. alternately computing (6) and running Algorithm 1 until convergence, referred to as AO in the figure) to that where the power control only scheme (i.e. Algorithm 1) is adopted in which AP-weighting coefficients are all set to one. It can be seen from Fig. 3 that the system performance improves when the number of APs increases (even the APdensity is fixed). Let us discuss scenarios (i) and (iii) first. Although the AP density is the same for both scenarios, the user-density (defined as the number of users per km 2 of area) for scenario (iii) is much smaller than that for scenario (i), and hence, users in scenario (iii) suffer less interference than users in scenario (i) do. Thus, users in scenario (iii) have higher rates than those in scenario (i). If we increase the number of users in scenario (iii) from 50 to 100 as in scenario (ii), then the user-density and thereby inter-user interference increase, which decreases the per-user rate accordingly. Furthermore, it is interesting to note that for scenarios (ii) and (iii), considering both power and receiver filtering control can deliver more universally good services to the users than the power control only scheme. However, this comes at the cost of extra complexity imposed by the computation of the receiver filter coefficients. Therefore, it is important to consider both the receiver coefficient design and power control problems for large-scale scenarios.
V. CONCLUSION In this paper, we have considered the max-min fairness problem for uplink cell-free massive MIMO. As in previous studies, we have decomposed the problem into the receiver filtering design problem and the power control problem. In particular, we have proposed a novel power control scheme based on the MP method. The numerical results have demonstrated that the proposed method is superior to other existing methods such as LP-based, GDA-based, and APG-based methods. More specifically, the proposed method achieves the same objective as these conventional schemes but in much lesser time. The proposed scheme has enabled us to investigate the performance of cell-free massive in large-scale scenarios. We have also numerically shown that the proposed solution provides better fairness among the users compared to the power control only scheme where AP-weighting coefficients are all set to one.
ACKNOWLEDGMENT This publication has emanated from research supported by a Grant from Science Foundation Ireland under Grant number 17/CDA/4786.