Pilot Assignment for Cell-Free Massive MIMO: A Spectral Clustering Approach

In this letter, we study the pilot assignment scheme in the uplink cell-free massive multiple-input multiple-output system. To accurately depict the potential co-pilot interference (CPI) among the user terminals (UTs), we introduce the normalized contribution value of AP to the received signal strength at the UT. A weighted conflict clustering graph is constructed by using the potential CPI to define the weight of the edge that connects the co-pilot UTs. Then, the pilot assignment optimization is transformed into a Min-Cut problem, which is further converted as a more effective normalized cut to find a partition of the UTs into disjoint clusters for minimizing the overall CPI. We develop a two-stage algorithm via spectral clustering to solve the relaxed problem efficiently. Numerical results are provided to validate the superiority of the proposed algorithm as compared to the existing methods in terms of the uplink average sum-rate and the quality of channel estimation.

in coherence interval, which causes the mutually correlated co-pilot interference (CPI) known as pilot contamination that degrades the sum-rate performance of the CF mMIMO systems.Thereby, it is of great significance to and timely to design pilot assignment (PA) schemes to mitigate the pilot contamination.
Recent progress has been made to explore the PA solutions in CF mMIMO systems.For instance, through discussing a random policy as the worst-case baseline, the authors in [1] presented a greedy PA scheme that iteratively updates the pilot sequence for the UT with minimum rate.However, the greedy policy can only increase the worst UT's rate without achieving the globally optimal performance.A location-based greedy (LBG) PA algorithm was designed in [5] by integrating the geographical circle area into the pilot reuse.However, the main idea behind LBG is to avoid the UTs with similar geographical locations using the same pilot.Since the conflict relations among the UTs are actually related to both locations of the UTs and APs, it is not accurate to describe the CPI among the UTs only through the locations of UTs.The authors in [6], [7] selected the service APs via the method proposed in [8], and quantified the interference relationship between UTs through the overlap of the service AP matrix.However, once the UTs are densely distributed, the primary service APs of each UT are very similar, so it is of little significance to construct interference graph through the service APs.Besides, it is difficult to effectively utilize all pilot sequences in [6], and it is even impossible for the interference aware user-group (IA-UG) scheme in [7] to determine the optimal threshold for the service AP selection in real system conditions.The authors in [4] defined an effective metric to quantify the potential CPI by treating the ratio of the large-scale fading coefficients between the target UT and the interfering UT at the service AP as the representation of interference from the co-pilot UTs.Moreover, a heuristic algorithm was proposed to solve the converted Max k-Cut problem for the intended PA solution.However, due to the diversity of large-scale fading coefficients, the APs that contribute little to the target UT may gain large weights in interference representation, which makes such a metric difficult to accurately depict the inter-UT CPI.
Against this background, we introduce the normalized contribution value of AP to accurately describe the potential CPI among all UTs.We particularly focus on a more practical case that there exist the centers of user activities in a large area and a large number of UTs are distributed near these centers.To our best knowledge, we are the first to study the PA scheme for a CF mMIMO system with tools from spectral clustering.The specific contributions of this letter are three-fold: • We utilize the normalized contribution value of AP to the received signal strength at the UT to quantify the potential CPI among all UTs.problem, by partitioning the UTs with low potential CPI into a cluster to minimize the overall CPI.• We further convert the Min-Cut problem into an effective normalized cut, and employ spectral clustering to develop an efficient two-stage algorithm to solve it.Notations: Upper (lower) case boldface letters are used for matrices (vectors).I L denotes L = L dimensional identity matrix.(•) * , (•) T , (•) H , tr(•), and E{•} denote conjugate, transpose, conjugated-transpose, trace, and expectation operations, respectively.CN (0, X) denotes the circularly symmetric complex Gaussian distribution with covariance matrix X.

II. SYSTEM MODEL
Consider the uplink of a CF mMIMO system, where M APs, each equipped with L antennas, coherently serve N (< M) single-antenna UTs over the same time and frequency resource blocks.All APs and UTs are randomly distributed within a large area.An error-free fronthaul network connects all APs with a central processing unit (CPU) for enabling the coordination of the coherent and joint transmissions.We adopt the more realistic spatially correlated Rayleigh fading channels [9].The correlated Rayleigh channel vector between UT n and AP m can be accordingly modeled as g n,m ∼ CN (0, R n,m ), where R n,m ∈C L×L is the spatial correlation matrix reflecting the effects of macroscopic wave propagation.With the covariance matrix R n,m , the large-scale fading coefficient can be written by β n,m tr(R n,m )/L, which is assumed to be known a priori at AP m.

A. Channel Estimation
Denote by τ c the coherence interval length (in samples), of which the first portion τ p is reserved for channel estimation.We focus on τ p <N , where a pilot sequence can be assigned to multiple UTs.Let ψ n ∈ C τp×1 be the τ p -dimensional pilot sequence assigned to UT n, and define P(n) as the set of UTs assigned with the same pilot as UT n.We represent the pilot set by P τp [ψ 1 , . . ., ψ τp ] ∈ C τp×τp , wherein every pair of pilots is orthogonal, i.e., ψ H ψ =1 if = , and ψ H ψ = 0, otherwise.In the uplink training phase, the received pilot signal y p m ∈C L×τp at AP m can be written as where ρ p is the normalized signal-to-noise ratio (SNR) of pilot sequence, and w p m ∈ C L×τp is the additive Gaussian noise vector of AP m, with all the elements being the independent and identically distributed (i.i.d.) CN (0, σ 2 I L ).Here, σ 2 is the noise power during the training phase.Via the projection of y p m onto ψ H n , the minimum mean-square error (MMSE) estimate of g n,m is specified by where yp n,m = ψ H n y p m .We shall note that the MMSE estimate ĝn,m is distributed as CN (0, G n,m ), where

B. Uplink Data Transmission
Let s n be the data signal transmitted by UT n with a power control coefficient η n ∈ [0, 1], for E{|s n | 2 } = 1.Due to the simultaneous data transmissions in the uplink across all N UTs, the data signal received by AP m is then expressed by where ρ u is the normalized uplink SNR, and w u m is the i.i.d.CN (0, σ 2 I L ) additive Gaussian noise at AP m.
To detect signal s n , with the decoding of maximal ratio combining (MRC), AP m sends ĝ * n,m y u m to the CPU via the fronthaul network.The combined signal received at the CPU for UT n is given as As detailed in [1], the combined received signal r u n can be then decomposed as (5) where DS n and BU n account for the strength of desired signal (DS) and the beamforming uncertainty (BU) for UT n, and CPI n denotes the CPI caused by UT n with the same pilot as UT n, for n = n .Following the worst-case of the uncorrelated Gaussian noise in (5), the uplink signal-to-interference-plusnoise ratio (SINR) for UT n can be modeled by ( 6), shown at the bottom of the next page.

III. SPECTRAL CLUSTERING BASED PILOT ASSIGNMENT A. Problem Formulation
This letter aims at improving the uplink sum-rate of the system by optimizing the PA across all UTs.We note from (2) that the accuracy about the MMSE estimate of channel vector for each UT is closely related to the corresponding CPI as in (5).More explicitly, the CPI for UT n is largely determined by the large-scale fading coefficient of UT n using the same pilot as UT n.Since the large-scale fading coefficients may vary greatly for different UTs, the potential CPI can differ as well for the different UTs.An efficient PA can avoid the same pilot reused by the UTs with the larger potential CPI, thus markedly reducing the estimation error for all UTs.This may result in the larger uplink SINR for the UT in (6), accordingly increasing its achievable rate.Therefore, the uplink sum-rate maximization problem can be thus transformed into an overall CPI minimization problem, which is formulated as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where C {ψ n , n ∈ {1, 2, . . ., N }} is a feasible PA solution.Problem ( 7) is an NP-hard problem [4], meaning that no polynomial algorithm exists to solve it optimally.An intuitive solution is to use the brute-force search.However, the computational cost for this method grows exponentially with system scale N, e.g., τ N p PA results [4], making it inefficient in practical systems.To solve this problem efficiently, we recognize that the assignment of pilots with limited dimensions to multiple UTs can be regarded as the partition of multiple vertices into different clusters from the graph theory perspective.Specifically, the UTs are set as the vertices of graph, and the potential CPI between the co-pilot UTs is set as the weight of the edge connecting them.By help of the graph cut paradigm [10], pilot contamination mitigation can be solved by finding a Min-Cut solution of the graph to approximate the optimal solution of (7).

B. Problem Transformation 1) Potential CPI Representation:
In order to characterize the potential CPI among all UTs, we first need to determine the interference relation between two UTs.The potential pilot contamination as defined in [4] can approximately quantify the potential CPI between two UTs under the user-centric approach.In this way, a UT is served only by the APs that it receives with best average channel conditions.However, in the CF approach of our interest, each AP communicates with all UTs, which results in the case that there exist the APs that contribute a negligible amount to the received signal strength among all the APs serving the UTs.We need a novel paradigm to balance the error caused by the difference in the dominant contributions of the APs to the received signal strength of a specific UT due to the large number of service APs. 1 With this in mind, we incorporate the normalized contribution value of AP to a specific UT.From (6), it follows that the CPI is closely related to the large-scale fading coefficient ratio of βn,m for two UTs served by AP m.Then, the potential CPI between UT n and UT n using the same pilot, for n = n , can be calculated as Remark: This potential CPI measurement is similar to evaluating the channel similarity of UTs.In the most extreme case, the channel coefficients of UT n and UT n are highly similar, and the potential CPI between them is 2, meaning that the potential CPI between them is large and the effect is bidirectional.
2) Conflict Clustering Graph Construction: We focus on the case that there are only small number of UTs assigned with the same pilot.In that case, the interference caused by multiple co-pilot UTs can be represented by the potential CPI.Then, the number of orthogonal pilots can be regarded as the number of clusters.Thereby, the PA can be achieved by dividing the UTs into different clusters, aiming to minimize the overall potential CPI.Before moving on, we provide the formal definition about the conflict clustering graph for tractability.
Definition 1: V} is a collection of subsets of V called the set of clusters, for i , j = 1, 2, . . ., |V| and x = 1, 2, . . ., |C|.In a weighted graph G, each edge {v i , v j } is labeled by a weight function For G, we note that each UT is treated as a vertex in V, and thus represent V {1, 2, . . ., N } as the set of vertices.Two co-pilot UTs causing the CPI form an edge in E, implying an unordered pair of two conflicting vertices in V. We then denote by |C| τ p the number of clusters.Based on the potential CPI between any two co-pilot UTs in (8), the weight of edge {n, n } can be given as ({n, n }) = 1/μ n,n .Note that the weight of an edge reflects the degree of similarity between two co-pilot UTs.The higher degree of similarity for two UTs is, the smaller weight is assigned to the edge, indicating the more serious conflict between them.
Regarding two disjoint clusters C x and C x in G, for x = x , the total weights of the edges that have been removed can be defined by the degree of dissimilarity between them, which is also called the cut [10].With (8), the cut between C x and C x in G can be calculated by Given a weighted graph G with vertex set V with N vertices, nonnegative edge weight ({n, n }), and pilot set P τp , our target for the overall potential CPI minimization in problem ( 7) is transformed to find a partition of V into τ p disjoint clusters C = {C 1 , C 2 , . . ., C τp } with the minimum total weights of the edges in E crossing the different clusters.To this end, the most direct way to obtain τ p disjoint clusters C = {C 1 , C 2 , . . ., C τp } is to solve the Min-Cut problem, which can be written as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

C. Algorithm Design via Spectral Clustering
Spectral clustering [10] is a graph theory based method, which can cut a graph into disjoint clusters with vertices in a cluster having large weights and vertices in different clusters having small weights.This observation motivates us to employ spectral clustering to solve the Min-Cut problem in (10) by partitioning the UTs with low potential CPI to form a cluster.
Note that the minimum cut criterion in problem (10) may cut small sets of isolated vertices in G.To form a reasonably large groups of vertices into the clusters, we adopt the normalized cut [10] as an unbiased measure, which calculates the cut cost as a fraction of the total edge connections to all the vertices in G.As such, the normalized cut Θ(•) to find a partition of where ) stands for the size of cluster C x , measured by the weights of the edges from vertices in C x to all vertices in V of G.
To relax (11) as a solvable form, we define a clustering matrix by Q = [q 1 , . . ., q τp ] N ×τp , with each column as an index vector q x = (q 1,x , . . ., q N ,x ) T whose elements are given by In this case, q T x Lq x = IV.NUMERICAL RESULTS This section verifies the performance of the proposed algorithm by the Monte-Carlo simulations.We configure a CF mMIMO system that all the APs and the UTs are positioned inside a square area of 0.5 km × 0.5 km, and their locations follow a random distribution and a standard bivariate normal distribution, respectively.The square area is wrapped around at the edges to imitate the boundary-free condition.The AP is equipped with L = 2 antennas, and the angular standard deviation in R n,m is 10 • .We use the three-slope model-based path loss in [1] to model the large-scale fading coefficients.The shadow fading coefficient S n,m ∼N (0, 8 2 ) is considered when the distance is larger than 50 m, and the coefficients are correlated as where d n,n is the distance between UT n and UT n , and d m,m is the distance between AP m and AP m .By integrating the channel estimation overhead, the per-UT uplink rate is expressed by , where τ c = 200 samples and bandwidth B = 20 MHz.The max-min power control scheme [1] is used for uplink data transmission.Unless otherwise stated, the values of other simulation parameters are set as ρp =100 mW, ρu =100 mW, and σ 2 =−96 dBm.
For comparison, we provide the performance obtained by the existing PA schemes as benchmarks: random scheme [1], greedy scheme [1], graph coloring scheme [6], LBG scheme [5], and IA-UG scheme with user grouping int.and AP selection threshold δ = 0.5 [7].The complexity of the benchmarks and Algorithm 1 are compared in Table I.The complexity of the proposed algorithm is only slightly higher than that of the IA-UG scheme and the graph coloring scheme Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE I COMPLEXITY COMPARISON OF THE BENCHMARKS
when τ p < N .We also consider a perfect case of CPI n =0 in (5), where all UTs are assigned to the pairwisely orthogonal pilots, called the orthogonal scheme.The quality of channel estimation for the benchmarks is evaluated via the metric of the normalized mean squared error (NMSE), defined by Fig. 1(a) shows the cumulative distribution of R u n under case 1 with M = 50 and case 2 with M = 150.We can see that in both cases, the proposed algorithm can uniformly improve the per-UT uplink rate, and only has 15.63% and 9.92% sum-rate losses, respectively, compared with the orthogonal scheme.The findings can be explained from the representation of potential CPI in (8), which can effectively reflect the CPI among the co-pilot UTs.Besides, the proposed algorithm can select the most suitable UTs to use the non-orthogonal pilots.
Fig. 1(b) plots the uplink average sum-rate of the system w.r.t.τ p .As shown, the average sum-rate first increases markedly and then decreases slowly as τ p grows.To explain, the increase of τ p results in the higher pilot reuse possibility, and however, incurs the larger channel estimation overheads in coherence interval.We also observe that as τ p increases, the performance curve of the proposed algorithm can approach the orthogonal scheme more significantly.This is due to the fact that the representation of potential CPI in (8) is tailored for two UTs, meaning that the smaller the number of UTs that reuse the pilots, the more accurate PA solution can be obtained to approach the global optimal solution of problem (7).
Fig. 1(c) depicts the NMSE performance w.r.t.M. From the results, we can see that the channel NMSE of all different schemes declines slowly with M. The reason is that the deployment of more APs for the UTs presents a decreasing trend of the distances between the UTs and the APs, resulting in the increased signal strengths received at the APs.As expected, the channel NMSE of the proposed algorithm is always lower than that of other PA schemes under different M.This verifies the effectiveness of the proposed algorithm under different AP layout conditions.
V. CONCLUSION An efficient spectral clustering-based PA scheme for the CF mMIMO uplink system was explored in this letter.We quantified the potential CPI among all UTs by help of the normalized contribution value of AP.The weighted conflict clustering graph was designed by using the potential CPI to calculate the edge's weight.We transformed the PA optimization problem into an effective normalized cut via spectral clustering, and proposed a two-stage algorithm to solve it.Validated via simulations, three useful remarks can be drawn from numerical results, providing insights for the PA design: 1) Our proposed algorithm has an uniform improvement of the per-UT uplink rate with low sum-rate loss compared to the orthogonal scheme; 2) Our proposed algorithm can significantly approach the orthogonal scheme in terms of the uplink average sumrate when the number of pilots is large enough; 3) The number of APs is not a main influencing factor on the NMSE performance, which, however, affects the cost of dedicated hardware for the practical system design.
value of AP m w.r.t.UT n.Note that ξ n,m is used to balance the potential CPI calculation error caused by the difference in the dominant contributions of different APs to the signal strength received by the UT.

7 :
2cut(Cx ,Cx ) vol(Cx ) , where L = D − Ξ is a Laplacian matrix representing the eliminated potential CPI between clusters.Here, Ξ=[1/μ n,n ] N ×N is a square matrix reflecting the potential CPI among the UTs, and D is a diagonal matrix to store the degree of vertex in V, where the n-th diagonal element d n is equal to d n = N n =1 ({n, n }).Thus, the Min-Cut problem in (10) can be relaxed into min {C1,C2,...,Cτ p } tr Q T LQ s.t.Q T LQ = I. (13) By relaxing the discreteness condition and substituting S = D 1/2 Q, problem (13) can be further relaxed as min {C1,C2,...,Cτ p } tr S T D −1/2 LD −1/2 S s.t.S T S = I. (14) Problem (14) is shown to be a standard trace minimization problem, which can be solved by finding k eigenvectors in matrix S. We first start with calculating k eigenvectors b 1 , b 2 , . . ., b k corresponding to the smallest k eigenvalues of D −1/2 LD −1/2 .Accordingly, we can form the matrix B [b 1 , . . ., b k ] ∈ R N ×k by stacking the eigenvectors in columns.After normalizing each row of matrix B via f i,j = b i,j / k l=1 b 2 i,l , for b i,j ∈ b j , i = 1, 2, . . ., N and j = 1, 2, . . ., k , we can then construct the matrix F=[f i,j ] N ×k ∈ R N ×k .By letting ζ i ∈ R k be the vector with regard to the i-th row of F, we propose an efficient two-stage algorithm via spectral clustering to derive a partition of C = {C 1 , C 2 , . . ., C τp } in G, which enables an effective solution C * = {ψ n , ∀n} to problem (14).The procedure of the proposed algorithm is elaborated in Algorithm 1. Complexity: The complexity of Algorithm 1 is from two stages.At stage I, the construction of conflict clustering graph is dominated by obtaining matrix Ξ, which entails O(N 2 M ) operations.At stage II, the computationally most expensive Algorithm 1 Proposed Spectral Clustering Based Algorithm Require: N, M, P τp , β n,m , ∀n, m.Ensure: C * = {ψ n , ∀n}.1: Set k = τ p = P τp .Stage I: Conflict Clustering Graph Construction 2: Obtain potential CPI matrix Ξ= 1/μ n,n based on (8).Stage II: Spectral Clustering 3: Obtain diagonal matrix D via d n = N n =1 n, n .4: Obtain Laplacian matrix L via L = D − Ξ. 5: Calculate the k eigenvectors b 1 , b 2 , . . ., b k corresponding to the smallest k eigenvalues of D −1/2 LD −1/2 .6: Normalize B with k eigenvalues via f i,j = Cluster the points ζ i ∈ R k , for i = 1, 2, . . ., N , with the k-means algorithm into C = {C 1 , C 2 , . . ., C k }. part is dominated by calculating the eigenvectors of normalized graph Laplacians D −1/2 LD −1/2 in Step 5 and using the k-means algorithm to perform the eigenvector classification in Step 7. The former requires O(N 3 ) operations [10], and the latter has the complexity of O(ιN τ 2 p ) [11], where ι is the maximum number of iterations needed for k-means algorithm.The complexity of obtaining matrices D in Step 3, L in Step 4, and B in Step 6 would be O(N 2 ), O(1), and O(N τ p ), respectively.Therefore, the total computational complexity of Algorithm 1 is O(N (NM + N 2 + ιτ 2 p )).

Fig. 1 .
Fig. 1.(a) The cumulative distribution of the per-UT uplink rate R u n , with N = 20, τp = 10; (b) The uplink average sum-rate versus the number of pilots τp, with M = 80, N = 40; (c) The NMSE performance versus the number of APs M, with N = 20, τp = 10.