Student Member, IEEE Meysam Sadeghi

Senior Member, IEEE Luca Sanguinetti
email: luca.sanguinetti@unipi.it

Senior Member, IEEE Romain Couillet
email: romain.couillet@centralesupelec.fr

Senior Member, IEEE Chau Yuen

Reducing the Computational Complexity of Multicasting in Large-Scale Antenna Systems

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

The advent of data-hungry services and applications has significantly increased the amount of data traffic of wireless networks [START_REF]Global mobile data traffic forecast update[END_REF]. A considerable amount of this traffic belongs to the services that are of interest to one or several groups of subscribers such as news headlines, financial data, regular system updates, and video broadcasting [START_REF]Global mobile data traffic forecast update[END_REF], [START_REF] Lecompte | Evolved multimedia broadcast/multicast service (eMBMS) in LTE-advanced: Overview and Rel-11 enhancements[END_REF]. The traditional unicast technology is highly inefficient for these services as it ignores the nature of such a traffic demand [START_REF] Lecompte | Evolved multimedia broadcast/multicast service (eMBMS) in LTE-advanced: Overview and Rel-11 enhancements[END_REF]- [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF].

To address this issue, the multicasting technology has been included in different releases of the third generation partnership project (3GPP) [START_REF] Lecompte | Evolved multimedia broadcast/multicast service (eMBMS) in LTE-advanced: Overview and Rel-11 enhancements[END_REF].

Physical layer multicasting is an efficient multicasting technique designed for wireless networks [START_REF] Sidiropoulos | Transmit beamforming for physical-layer multicasting[END_REF], [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF]. It has been widely studied in the literature either for single or multiple groups [START_REF] Sidiropoulos | Transmit beamforming for physical-layer multicasting[END_REF]- [START_REF] Chen | Weighted max-min fairness oriented load-balancing and clustering for multicast cache-enabled F-RAN[END_REF]. In single-group multicasting, a transmitter exploits its channel state information (CSI) to send out a common stream of data to one group of users, while in multi-group multicasting multiple independent streams of common data are sent to multiple distinct groups of users. In this context, two classes of problems have received particular attention, the so-called quality-ofservice (QoS) problem and the weighted max-min-fairness (MMF) 1 problem. The former aims to minimize the total transmit power while satisfying target signal-to-interference-plus-noise ratios (SINRs) at the active user equipments (UEs). The latter seeks to maximize the minimum weighted SINR among all the UEs in the system, subject to a total transmit power constraint.

A seminal treatment of single-group multicasting for both QoS and MMF problems was first presented in [START_REF] Sidiropoulos | Transmit beamforming for physical-layer multicasting[END_REF]. Therein, it was proved that both QoS and MMF problems are NP-hard and then an approximate solution was presented employing semidefinite relaxation (SDR) technique [START_REF] Luo | Semidefinite relaxation of quadratic optimization problems[END_REF]. This work is then extended to a multi-group single-cell scenario in [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF]. It should be noted that, since in the multi-group case the SINR of every UE depends on the precoding vectors of all other groups, even finding a feasible solution for QoS and MMF problems might be a challenging task [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF], [START_REF] Christopoulos | Multicast multigroup beamforming for per-antenna power constrained large-scale arrays[END_REF]. Therefore, in [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF] the SDR technique is followed by a randomization and a multi-group multicast power control phase. In [START_REF] Christopoulos | Weighted fair multicast multigroup beamforming under per-antenna power constraints[END_REF], the MMF problem is studied under per-antenna power constraint for multi-group single-cell systems. The coordinated physical layer multicasting for single-group multi-cell scenario is investigated in [START_REF] Xiang | Coordinated multicast beamforming in multicell networks[END_REF]. Also, its application to coordinated physical layer multicasting for multi-group multi-cell scenario is studied in [START_REF] Chen | Weighted max-min fairness oriented load-balancing and clustering for multicast cache-enabled F-RAN[END_REF].

The aforementioned works (among many others) are based on the SDR technique, which is characterized by high computational complexity when the system dimensions grow large, especially for large antenna arrays. More precisely, consider a single-cell network wherein a BS with N antennas serves K UEs in G multicasting groups. Then, solving the QoS problem via SDR requires O( √ GN ) iterations of an interior point method with each iteration requiring O(G 3 N 6 + KGN 2 ) arithmetic operations [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF]. The computational cost of finding an approximate solution for the MMF problem is even higher as its solution is achieved by iteratively solving different instances of the QoS problem. Therefore, the SDR-based solutions are not suitable for practical implementation when N , G, or K, grow large, as envisioned in large-scale antenna systems (or more commonly known as massive MIMO systems) wherein N can be of the order of hundreds [START_REF] Marzetta | Noncooperative cellular wireless with unlimited numbers of base station antennas[END_REF]- [START_REF] Arvola | Two-layer precoding for dimensionality reduction in massive MIMO[END_REF].

In the context of massive MIMO multicasting, two possible approaches have been recently

proposed, namely, the asymptotic approach and the successive convex approximation (SCA) approach [START_REF] Christopoulos | Multicast multigroup beamforming for per-antenna power constrained large-scale arrays[END_REF], [START_REF] Xiang | Massive MIMO multicasting in noncooperative cellular networks[END_REF]- [START_REF] Xiang | Massive MIMO multicasting in noncooperative multicell networks[END_REF]. The former exploits the asymptotic orthogonality of the channel vectors, when N grows very large and K is kept fixed, to simplify the SINR expression of each UE and facilitate the design of asymptotically optimal beamforming schemes [START_REF] Xiang | Massive MIMO multicasting in noncooperative cellular networks[END_REF]- [START_REF] Sadeghi | Multi-cell multi-group massive MIMO multicasting: An asymptotic analysis[END_REF]. In particular, [START_REF] Zhou | Joint multicast beamforming and user grouping in massive MIMO systems[END_REF] investigates the MMF problem for the multi-group single-cell multicasting whereas the singlegroup multi-cell case is studied in [START_REF] Xiang | Massive MIMO multicasting in noncooperative cellular networks[END_REF]. The extension to a multi-group multi-cell network is considered in [START_REF] Sadeghi | Multi-cell multi-group massive MIMO multicasting: An asymptotic analysis[END_REF]. The main problem with the asymptotic approach is that an extremely large number of antennas is required to reach the asymptotic orthogonality condition. As a consequence, the performance of the asymptotically optimal precoders is poor when the system does not have an extremely large number of antennas, i.e., N < 1000 [START_REF] Sadeghi | Multi-cell multi-group massive MIMO multicasting: An asymptotic analysis[END_REF].

The SCA approach aims at iteratively solving the non-convex QoS and MMF problems by means of SCA of the original problems around a feasible point [START_REF] Marks | A general inner approximation algorithm for nonconvex mathematical programs[END_REF], [START_REF] Beck | A sequential parametric convex approximation method with applications to nonconvex truss topology design problems[END_REF]. More specifically, the algorithm starts from an initial feasible point, the non-convex constraints are approximated by convex functions around this point, and the resulting convex problem is solved to obtain the next iteration. This procedure is repeated until convergence to a stationary point. In [START_REF] Tran | A conic quadratic programming approach to physical layer multicasting for large-scale antenna arrays[END_REF],

the SCA technique has been applied to reduce the computational complexity of beamforming design in single-group multicasting for large-scale antenna arrays. However, the SCA method is not suitable for multi-group multicasting communications as it requires an initial feasible point, which is hard to compute in these scenarios [START_REF] Christopoulos | Multicast multigroup beamforming for per-antenna power constrained large-scale arrays[END_REF]. To handle this issue, a feasible point pursuit SCA (FPP-SCA) algorithm is proposed in [START_REF] Mehanna | Feasible point pursuit and successive approximation of non-convex QCQPs[END_REF] and applied to multi-group multicasting in [START_REF] Christopoulos | Multicast multigroup beamforming for per-antenna power constrained large-scale arrays[END_REF]. Therein, in order to guarantee the feasibility of the problem, slack variables are added to relax the constraints, and a penalty is used to ensure that slacks are sparingly used. The solution of the resulting optimization problem is then used for another round of approximation and the procedure is repeated until convergence. However, the method itself has two drawbacks.

First, although the solution of the approximated problem is always feasible, it might not be a feasible solution of the original multicasting problem and it is sensitive to the initial point of the algorithm as it is detailed in [START_REF] Mehanna | Feasible point pursuit and successive approximation of non-convex QCQPs[END_REF]. Second, it is still computationally demanding when the number of antennas grows.

In this paper, we address all the aforementioned drawbacks for a multi-group single-cell large-scale antenna system by introducing a two-layer precoding scheme, which is tailored for large-scale antenna systems. Our main contributions are summarized as follows:

1) We present two algorithms for the QoS and MMF problems, that outperform most of the aforementioned solutions while guaranteeing a low computational complexity.

2) We reveal new duality results that allow to solve both QoS and MMF problems simultaneously. This is in sharp contrast with the existing algorithms for which the MMF problem is solved by iteratively solving different instances of the QoS problem.

3) We introduce a heuristic algorithm that significantly improves the computational complexity while only slightly reducing the performance of both QoS and MMF solutions.

The rest of the paper is organized as follows. Section II introduces the system model for a multi-group single-cell large-scale antenna array system and formulates the corresponding QoS and MMF problems. Section III introduces the proposed two-layer precoder, it provides a duality result between transformed versions of the QoS and MMF problems, and then it proposes two algorithms for solving both. Section IV introduces a heuristic solution to further reduce the computational complexity. Section V presents the numerical results whereas conclusions are drawn in Section VI.

Notations:

The following notation is used throughout the paper. Scalars are denoted by lower case letters whereas boldface lower (upper) case letters are used for vectors (matrices). We denote by 0 a matrix of appropriate size where all its elements are zero. The transpose, conjugate transpose, real part, absolute value, and second norm operator are denoted by (•) T , (•) H , Re(•), | • |, and • . The set of all positive real numbers is denoted by R + . A circular symmetric complex Gaussian random vector x is denoted by x ∼ CN (0, C), where 0 and C are its mean and covariance matrix, respectively. The inverse of an invertible function f (.) is shown by f -1 (.).

II. SYSTEM MODEL AND PROBLEM FORMULATION

Consider a single-cell large-scale antenna array system in which a BS equipped with N antennas serves G multicasting groups. Denote by G = {1, . . . , G} the set of indices of all groups and call K j the set of UE indices associated with group j, with cardinality

K j = |K j |
and such that K j ∩ K i = ∅, j = i, i.e., each UE is associated with a single group. Within this setting, we assume that N > K -min j∈G K j , where K = G j=1 K j is the total number of UEs in the network. Since large antenna systems are considered in this work, this technical assumption is naturally in place. A double index notation is used to refer to each UE as e.g., "user k in group j". Under this convention, let g jk ∈ C N be the channel between UE k in group j and the BS and assume that g jk = β jk h jk , where h jk ∼ CN (0 N , I N ) is the small-scale fading channel and β jk accounts for the large-scale channel attenuation (or path loss). We assume that the BS has perfect knowledge of the channel vectors {g jk }.

Denoting by w j ∈ C N the precoding vector associated with group j, the signal y jk received at UE k can be written as:

y jk = g H jk w j s j + G i=1,i =j g H jk w i s i + n jk (1)
where s i ∼ CN (0, 1) is the signal intended to group i, assumed independent across i, and n jk ∼ CN (0, σ 2 jk ) accounts for the additive Gaussian noise. The SINR of UE k in group j, γ jk , can be written as

γ jk = g H jk w j 2 G i=1,i =j g H jk w i 2 + σ 2 jk (2)
and the total average transmit power is G j=1 w j 2 . Under the above assumptions, an instance of the QoS problem can be formulated as follows [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF]:

Q(η) : min {w j } G j=1 w j 2 (3)
s.t. γ jk ≥ η jk ∀j, k (4)
where η jk is the prescribed SINR of UE k in group j and η ∈ C K is the vector collecting all the {η jk }. Accordingly, an instance of the MMF problem can be written as follows [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF]:

F(η, P ) : max {w j } min j min k 1 η jk γ jk (5) s.t. G j=1 w j 2 ≤ P (6)
where P accounts for the power constraint at the BS, and 1 η jk represents the weight of γ jk . As mentioned before, Q(η) and F(η, P ) are NP-hard and the existing algorithms for computing their approximate solutions have either high computational complexity [START_REF] Sidiropoulos | Transmit beamforming for physical-layer multicasting[END_REF], [START_REF] Christopoulos | Weighted fair multicast multigroup beamforming under per-antenna power constraints[END_REF], [START_REF] Christopoulos | Multicast multigroup beamforming for per-antenna power constrained large-scale arrays[END_REF], or poor performance [START_REF] Xiang | Massive MIMO multicasting in noncooperative cellular networks[END_REF], [START_REF] Sadeghi | Multi-cell multi-group massive MIMO multicasting: An asymptotic analysis[END_REF], [START_REF] Xiang | Massive MIMO multicasting in noncooperative multicell networks[END_REF]. A two-layered architecture is proposed next to overcome these drawbacks.

III. THE PROPOSED TWO-LAYER PRECODING SCHEME

In this section, we propose a simple and computationally efficient method to compute approximate solution to the QoS and MMF problems. The method is based on a two-layer precoding scheme: (i) the outer layer restricts the space of valid precoders to those cancelling the inter-group interference, thereby approximating the QoS and MMF problems by simpler (still non-convex) problems, denoted by QoS dec and MMF dec , for which trivial feasible points can be found; (ii) starting from these feasible points, the inner layer is designed to reach a suboptimal solution to the QoS dec and MMF dec problems, which are also feasible solutions of the original QoS and MMF problems. Section III-A presents the outer layer. Section III-B reveals an explicit duality between the QoS dec and MMF dec problems. Section III-C presents the inner layer and the algorithms developed. Section III-D evaluates the complexity of the proposed algorithms.

A. Outer Layer -Removing Multi-group Interference

Denote by G i ∈ C N ×K i the matrix collecting the channel vectors of all the K i UEs in group i.

The complete elimination of the multi-group interference G i=1,i =j g H jk w i s i in (1) is achieved by using the block-diagonalization zero-forcing (BDZF) technique [START_REF] Spencer | Zero-forcing methods for downlink spatial multiplexing in multiuser MIMO channels[END_REF], [START_REF] Choi | A transmit preprocessing technique for multiuser MIMO systems using a decomposition approach[END_REF]. Consider a two-layer precoding vector for group j as follows

w j = F j c j ∀j ∈ G (7)
where c j ∈ C N -τ j with τ j = K -K j is the inner layer, the design of which is discussed later, and F j ∈ C N ×(N -τ j ) is the outer layer. We design F j as an isometric matrix whose columns form a basis for the null space of

G -j = [G 1 , . . . , G j-1 , G j+1 , . . . , G G ] ∈ C N ×τ j , i.e.,
G H -j F j = 0 τ j ×(N -τ j ) . As proposed in [START_REF] Spencer | Zero-forcing methods for downlink spatial multiplexing in multiuser MIMO channels[END_REF], [START_REF] Choi | A transmit preprocessing technique for multiuser MIMO systems using a decomposition approach[END_REF], F j can be obtained through the singular value decomposition (SVD) of G -j . This requires 16(G -1)KN2 + 24N G j=1 (K -K j ) 2 floating point operations (flops) 2 . The same goal can be obtained with lower complexity (linear in the number of BS antennas N ) using the QR-based decomposition approach as shown in [START_REF] Chen | Transmit selection diversity for unitary precoded multiuser spatial multiplexing systems with linear receivers[END_REF]. This produces

G -j = Q j R j = Q j Q j   R j 0   = Q j R j ( 8
)
where

Q j ∈ C N ×(N -τ j )
gives the null space of G -j such that G H -j Q j = 0 τ j ×(N -τ j ) . Therefore we can use Q j , as the outer layer of w j , i.e., F j = Q j . Since the QR decomposition of an m by n matrix can be computed with 8mn 2 -8 3 n 3 flops [START_REF] Arakawa | Computational workloads for commonly used signal processing kernels[END_REF], the total number of flops required to perform the BDZF technique reduces to

8N G j=1 (K -K j ) 2 -8 3 G j=1 (K -K j ) 3 , which increases linearly with N . Plugging F j = Q j into (2) yields γ jk = |g H jk c j | 2 (9)
where

g jk = 1 σ jk (Q j ) H g jk ∈ C N -τ j
denotes the equivalent channel vector of UE k in group j. As ∀j (F j ) H F j = I N -τ j , the proposed outer layer does not change the transmit power, i.e., G j=1 w j 2 = G j=1 c j 2 . Therefore, using the BDZF technique the QoS problem reduces to QoS dec . Note that we define the QoS dec as the transformed version of the QoS problem into G single-group multicasting QoS problems, {Q j (η j )} G j=1 , where the jth problem is given by

Q j (η j ) : min {c j } c j 2 (10)
s.t. g H jk c j 2 ≥ η jk ∀k (11)
where η j ∈ C K j is the vector collecting all the quantities {η jk } in group j. To grasp the relation between the QoS dec problem and the prescribed SINRs, i.e., {η j } G j=1 , we denote an instance of the QoS dec problem by Q(η) = {Q j (η j )} G j=1 . Accordingly, using the BDZF technique, the MMF problem reduces to MMF dec , where we show an instance of it by F(η, P ) and it is given as follows F(η, P ) : max

{c j } min j min k 1 η jk g H jk c j 2 (12) s.t. G j=1 c j 2 ≤ P. ( 13
)
As mentioned in the introduction, finding a feasible point for Q(η) is hard [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF]. The same holds for F(η, P ) since the common approach to solve F(η, P ) relies on iteratively solving Q(η). On the contrary, finding a feasible point for Q(η) = {Q j (η j )} and, thus, for F(η, P ) is a straightforward task, thanks to the single-group nature of each subproblem Q j (η j ). More specifically, for any given group j an initial feasible point can be computed by first choosing an arbitrary beamforming vector c j and then rescaling it so as to meet the most violated SINR constraint with equality. Despite being simpler to solve than Q(η) and F(η, P ), both Q(η)

and F(η, P ) are still NP-hard -as it easily follows observing that for G = 1 they reduce to the single-group problems studied in [START_REF] Sidiropoulos | Transmit beamforming for physical-layer multicasting[END_REF]. The main difficulty in solving Q(η) and F(η, P ) lies in the non-convexity of the SINR constraints. In Section III-C, the SCA technique is used to develop a possible solution capable of overcoming this issue.

Before delving into this, next we further detail the characteristics of Q(η) and F(η, P ) and establish a duality and direct relations between the two problems: the solution of F(η, P ) can be obtained from that of Q(η) (and vice versa). These results will be used in Section III-C and Section IV to compute an approximate solution to F(η, P ) by means of Q(η) without the need of iteratively solving instances of Q(η) as for existing alternatives [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF]- [START_REF] Xiang | Coordinated multicast beamforming in multicell networks[END_REF].

B. On the duality between Q(η) and F(η, P )

Let {c j (η)} and P (η) denote the set of optimal precoding vectors and the optimal objective value of Q(η), respectively. Similarly, let {c • j (η, P )} and t • (η, P ) denote the set of optimal precoding vectors and the optimal objective value of F(η, P ). We then start providing the following result: Lemma 1. For Q(η) and F(η, P ) we have:

c j (αη) = c • j (η, P (αη)) ∀j ∈ G (14)
with α = t • (η, P (αη)). Also, we have that:

c • j (η, P ) = c j (t • (η, P )η) ∀j ∈ G (15)
with P = P (t • (η, P )η).

Proof. The proof proceeds by contradiction. First, notice that, by definition, {c j (αη)} is a feasible solution of F(η, P (αη)) with an objective value equal to α. Now, let us assume there exists a set of precoding vectors {c • j (η, P (αη))} for which t • (η, P (αη)) > α. Clearly, {c • j (η, P (αη))} is also a feasible solution of Q(αη) for which all the SINR constraints are over satisfied. Hence, there exists a constant ν < 1 such that {νc • j (η, P (αη))} meets all the SINR constraints of Q(αη) with equality while providing a smaller objective value than P (αη). This, however, contradicts our assumption and proves that ( 14) is valid with α = t • (η, P (αη)). A similar line of reasoning can be used to prove [START_REF] Luo | Semidefinite relaxation of quadratic optimization problems[END_REF]. By definition the set of precoding vectors {c • j (η, P )} is a feasible solution of Q(t • (η, P )η) with an objective value equal to P . Let us assume there exists {c j (t • (η, P )η)} with P (t • (η, P )η) < P . Then, one can use the remaining power, P -P (t • (η, P )η) to rescale {c j (t • (η, P )η)} and improve F(η, P ). This is in contradiction with our assumption and completes the proof. Also, the following lemma can be simply proved from the definition of Q(η): Lemma 2. For Q(αη) and ∀α ∈ R + we have

P (αη) =αP (η) (16)
and ∀j ∈ G

c j (αη) = √ αc j (η). (17)
We are now ready to state the following explicit duality between Q(η) and F(η, P ):

Theorem 1. Given the set of optimal precoding vectors and the optimal objective value of Q(η),

i.e., {c j (η)} and P (η), the set of optimal precoding vectors and the optimal objective value of F(η, P ), i.e. {c • j (η, P )} and t • (η, P ), are determined as

c • j (η, P ) = P P (η) c j (η) ∀j ∈ G (18) t • (η, P ) = P P (η) (19)
and vice versa as

c j (η) = 1 t • (η, P ) c • j (η, P ) ∀j ∈ G (20)
P (η) = P t • (η, P ) . ( 21
)
Proof. Starting with [START_REF] Hoydis | Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?[END_REF] we have that

P P (η) c j (η) (a) = c j P P (η) η (b) = c • j η, P P P (η) η ( 22
) (c) = c • j η, P P (η) P (η) = c • j (η, P ) (23)
where (a) follows from ( 17), (b) holds because of ( 14) and (c) is obtained using [START_REF] Marzetta | Noncooperative cellular wireless with unlimited numbers of base station antennas[END_REF]. The equality [START_REF] Marzetta | How much training is required for multiuser MIMO?[END_REF] follows from

P (a) = P (t • (η, P )η) (b) = t • (η, P )P (η) (24)
where (a) exploits P = P (t • (η, P )η) (see Lemma 1), and (b) is due to [START_REF] Marzetta | Noncooperative cellular wireless with unlimited numbers of base station antennas[END_REF]. The equality in [START_REF] Arvola | Two-layer precoding for dimensionality reduction in massive MIMO[END_REF] follows from replacing [START_REF] Marzetta | How much training is required for multiuser MIMO?[END_REF] in [START_REF] Hoydis | Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?[END_REF].

Theorem 1 reveals the relation between the optimal precoding vectors and the optimal objective values of Q(η) and F(η, P ). However, as they are NP-hard, any arbitrary algorithm with polynomial complexity can provide an approximate set of precoding vectors, rather than the optimal one. Therefore, it is interesting to establish a relation between the precoding vectors and the objective values of Q(η) and F(η, P ) while they are achieved from any arbitrary sub-optimal algorithm. This relation is given in Propositions 1 and 2.

Proposition 1. Assume {c j,app (η)} is a set of precoding vectors of Q(η) and P app (η) is its associated objective value achieved by any arbitrary algorithm. Then, the set of precoding vectors

{ P P app (η) c j,app (η)} (or { P P app (η) F j c j,app (η)})
is a feasible answer for F(η, P ) (or F(η, P )), and provides an objective value t • app (η, P ) such that t • app (η, P ) ∈ [ P P app (η) , P P (η) ].

Proof. Please refer to the Appendix A.

Proposition 2. Assume {c • j,app (η, P )} is a set of precoding vectors of F(η, P ) and t • app (η, P ) is its associated objective value achieved by any arbitrary algorithm. Then, the set of precod-

ing vectors { 1 √ t • app (η,P ) c • j,app (η, P )} (or { 1 √ t • app (η,P ) F j c • j,app (η, P )}),
is a feasible answer for Q(η) (or Q(η)), and provides an objective value P app (η) such that P app (η) ∈ [ P t • (η,P ) ,

P t • app (η,P ) ].
Proof. Please refer to the Appendix B.

Note that the relation between the QoS and MMF problems was first discovered in [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF], but it was not given in an explicit form. Therefore, the existing works in the literature, as [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF]- [START_REF] Xiang | Coordinated multicast beamforming in multicell networks[END_REF], [START_REF] Christopoulos | Multicast multigroup beamforming for per-antenna power constrained large-scale arrays[END_REF], solve the MMF problem by iteratively solving specific instances of the QoS problem. By virtue of the large number of antennas available in large-scale antenna systems and the BDZF technique, Theorem 1, Proposition 1, and Proposition 2, state that F(η, P ) and Q(η) (also F(η, P ) and Q(η)) can be solved simultaneously. It is also interesting to observe that the upper bound of the objective value of F(η, P ) achieved via Proposition 1 is equal to [START_REF] Marzetta | How much training is required for multiuser MIMO?[END_REF]. Also, the lower bound of the objective value of Q(η) achieved via Proposition 2 is equal to [START_REF] Xiang | Massive MIMO multicasting in noncooperative cellular networks[END_REF].

C. Inner Layer -Successive Convex Approximation

In the sequel, the SCA technique is applied to solve Q(η) and F(η, P ). We begin with Q(η),

and rewrite |ḡ H jk c j | 2 as |ḡ H jk c j | 2 = c H j X jk c j (25)
where X jk = ḡjk ḡH jk is a rank-one positive semi-definite matrix. Thus, for any arbitrary vector z j ∈ C N -τ j we have that (c j -z j ) H X jk (c j -z j ) ≥ 0 from which it follows

c H j X jk c j ≥ 2Re(z H j X jk c j ) -z H j X jk z j . (26)
Now, for any z j the non-convex SINR constraint c H j X jk c j ≥ η jk can be replaced with a tighter convex constraint given by 2Re

(z H j X jk c j ) -z H j X jk z j ≥ η jk . (27)
By replacing [START_REF] Christopoulos | Multicast multigroup beamforming for per-antenna power constrained large-scale arrays[END_REF] with [START_REF] Marks | A general inner approximation algorithm for nonconvex mathematical programs[END_REF], we obtain

Q j (η j , z j ) : min {c j } c j 2 (28)
s.t. 2Re(z H j X jk c j ) -z H j X jk z j ≥ η jk ∀k (29)
which represents a convex approximation of Q j (η j ) for a specific instance of z j . Now, we can introduce Algorithm 1 and its following proposition for the QoS problem.

Algorithm 1 The QoS BDZF-SCA Algorithm 1: Compute F j ∀j ∈ G.

2: for j = 1 to G do 3:

Select an arbitrary z

(1) j

and rescale it such that ∀k z

(1)H j X jk z

(1) j ≥ η jk .

4:

repeat 5:
Solve:

Q j (η j , z (i) j ) : min c (i) j c (i) j 2 (30)
s.t. 2Re(z (i)H j X jk c (i) j ) -z (i)H j X jk z (i) j ≥ η jk ∀k. (31) 6:
Let c

(i) j denote the optimal value obtained from Q j (η j , z

(i) j ), then set z (i+1) j ← c (i) j . 7:
until Convergence 8: end for 9: Compute the precoding vectors w j = F j c j ∀j ∈ G.

Proposition 3. Algorithm 1 converges to a point satisfying the KKT conditions of Q(η), while providing a feasible solution for Q(η).

Proof. Please refer to the Appendix C. Now let us consider F(η, P ) and F(η, P ). A solution to these two problems can be achieved by first applying Algorithm 1 and then using Proposition 1. Besides, we can directly apply the SCA technique to F(η, P ) and find a solution to these two problems, similar to Algorithm 1.

The latter approach is presented in Algorithm 2 and we have the following proposition for Algorithm 2.

Proposition 4. Algorithm 2 converges to a KKT point of F(η, P ), while providing a feasible solution to F(η, P ).

Proof. The proof follows the same lines as the proof of Proposition 3.

Algorithm 2

The MMF BDZF-SCA Algorithm 1: Compute F j ∀j ∈ G.

2: Select an arbitrary set z (1)

:= {z (1) j } G j=1 such that G j=1 z
(1) j 2 ≤ P.

3: repeat 4:

Solve:

F(η, P, z (i) ) : max

{c (i) j } min j min k 1 η jk 2Re(z (i)H j X jk c (i) j ) -z (i)H j X jk z (i) j (32)
s.t. G j=1 c (i) j 2 ≤ P. (33) 5:
Let {c

(i) j } G j=1
denote the optimal value obtained from F(η, P,

z (i) ), then ∀j ∈ G set z (i+1) j ← c (i) j . 6: until Convergence 7:
Generate the precoding vectors of F(η, P ), ∀j ∈ G w j = F j c j .

D. Computational Complexity

The computational load of Algorithms 1 and 2 is now assessed in terms of the number of required flops as follows. Note that both algorithms consist of three steps. The first step computes {F j ; ∀j ∈ G} and requires 8N G j=1 (K -K j ) 2 -8 3 G j=1 (K -K j ) 3 flops using the QR-decomposition [START_REF] Arakawa | Computational workloads for commonly used signal processing kernels[END_REF], [START_REF] Chen | Transmit selection diversity for unitary precoded multiuser spatial multiplexing systems with linear receivers[END_REF]. The second step aims at designing the inner layer precoding vectors {c j ; ∀j ∈ G} -as detailed in lines 2 to 8 (2 to 6) of Algorithm 1 (Algorithm 2). Since Q(η, z) and F(η, P, z) are both convex, they can be solved at each iteration using standard techniques with a worst case complexity of O(N 3 ) [START_REF] Mehanna | Feasible point pursuit and successive approximation of non-convex QCQPs[END_REF]. Therefore, the number of flops required by the second step is O(IN 3 ) with I being the number of iterations required to converge. As it will be observed in Section V, only a few iterations are needed to reach a satisfying solution even for large N . The third step calculates the composite precoding vectors w j = F j c j and requires

IV. A HEURISTIC INNER LAYER OF ORDER O(N )

In the previous section, it was shown that the complexity of the proposed algorithms is of O(N 3 ), which is due to the application of SCA technique to find the inner layer precoding vectors, i.e., {c j } G j=1 ,. Therefore, the inner layer retrieval may still be computationally expensive when N is relatively large. Moreover, it requires optimization packages for solving the convex problems Q(η, z) and F(η, P, z), which may not be available on every hardware platform.

Therefore, in what follows, we present a simple, yet effective, heuristic algorithm for computing the inner layer precoding vectors of Q(η) with a complexity of O(N ). Then, by employing Proposition 1 and the solution obtained for Q(η), we compute an approximate solution for F(η, P ). Therefore, the complexity of simultaneously finding an inner layer precoder for both problems becomes O(N ).

The proposed heuristic algorithm aims at computing the precoding vector c j ∀j ∈ G, in K j sequential steps. The algorithm has two main parts, the ordering part and the successive precoder design part. Assuming that the K j UEs in group j are labeled from 1 to K j , the ordering part will re-label them by a bijective function f j : {1, . . . , K j } → {µ j1 , . . . , µ jK j }, where µ jk = f j (i) means that the UE who was labeled as i is now re-labeled as µ jk and will be served in kth step of the algorithm, k ∈ {1, . . . , K j }. Therefore, the new labels, {µ jk } K j k=1 , will determine the order by which the UEs in group j are served in each step. The successive precoder design part, designs the precoding vector of group j in K j steps such that in kth step the requested SNR of UE µ jk is met with minimum power consumption while the SNR of the previous k -1 ordered UEs, i.e., {µ jt } k-1 t=1 , is not violated. We will detail the successive precoder design and the user ordering in the following two subsections.

A. The Successive Precoder

Assume that ∀j ∈ G the UE ordering is given, i.e., {µ jk } K j k=1 is known. Denote by c (k) j the precoding vector c j at kth step, then it is computed as follows:

c (k) j = c (k-1) j + α (k) j d (k) j k ∈ {1, . . . , K j } (34)
where d

(k) j ∈ C N -τ j is a unit norm vector and α should be designed such that the SNR constraint of µ jk is met with minimum power consumption while the SNR of {µ jt } k-1 t=1 is not violated.

We start by initializing the precoding vector c

(1) j of UE µ j1 such that its own SNR constraint, i.e., g H jµ j1 c

(1) j 2 ≥ η jµ j1 , is met with equality. This yields c

(1) j = √ η jµ j1 g jµ j1 2 g jµ j1 . For k ∈ {2, . . . , K j }, the vectors d (k) j must be chosen such that the previously satisfied k -1 SINR constraints are not violated. This is achieved by selecting d

(k) j orthogonal to {g jµ ji } k-1 i=1 , i.e., g H jµ ji d (k) j = 0 for i = 1, . . . , k -1. To this end, {d (k) j } K j
k=2 are computed using the Gram-Schmidt procedure, which produces d

(k) j = u (k) j u (k) j with u (k) j = g jµ jk - k-1 i=1 u (i)H j g jµ jk u (i) j 2 u (i) j (35)
being the component of g jµ jk orthogonal to the subspace spanned by {u

(i) j } k-1 i=1 . Once the unit norm vectors d (k) j
are computed, we proceed with the design of coefficients {α

(k) j } K j k=2 .
In particular, each α (k) j is chosen such that the power consumption in step k, given by ||c

(k) j || 2 = ||c (k-1) j || 2 + |α (k) j | 2 ,
is minimized while satisfying the kth SINR constraint. More precisely, α (k) j must be computed as the solution of the following problem:

min α (k) j |α (k) j | 2 s.t. g H jµ jk c (k) j 2 ≥ η jµ jk . (36)
As shown in the Appendix D (see also [START_REF] Hunger | Design of single-group multicasting-beamformers[END_REF]), the optimal α (k)

j = |α (k) j | exp(i ∠α (k) j ) is computed as: ∠α (k) j = -∠ρ (k) j (37) |α (k) j | = -ρ (k) j + ρ (k) j 2 -g H jµ jk d (k) j 2 g H jµ jk c (k-1) j 2 -η jµ jk g H jµ jk d (k) j 2 (38) with ρ (k) j = g H jµ jk d (k) j c (k-1) H j g jµ jk .
In the sequel, the above results are used to sort the UEs according to a worst-first policy, which is observed to achieve close-to-optimal performance by means of numerical results in Section V.

B. User Ordering

At this stage, we are only left with the computation of the UE ordering indices {µ jk }. A possible solution is illustrated in [START_REF] Hunger | Design of single-group multicasting-beamformers[END_REF], for the QoS problem in single-group multicasting systems.

More specifically, denote by S the set of indices of the remaining ones, i.e., R

(k-1) j = {1 . . . , K j } \ {f -1 j (µ jt )} k-1 t=1 . Then, in [12] the set S (k) j is computed as S (k) j = S (k-1) j ∪ {µ jk } with µ jk = arg min i∈R (k-1) j |g H ji c (k-1) j | 2 η ji (39)
corresponding to the UE index in R (k-1) j that has the weakest ratio (or also the most violated constraint) between the provided SNR up to step k, given by |g H ji c

(k-1) j

| 2 , and the requested one

given by η ji . The above procedure has the following two drawbacks. Firstly, it needs to calculate at each stage k the quantities |g H ji c

(k-1) j | 2 ∀i ∈ R (k-1) j
, which requires O(K 2 j N ) flops for group j. This is costly if N and K j are large. Secondly, it does not take into account the extra amount of power |α | 2 /η ji takes a very high value. This might happen, for example, because its own channel vector g ji is almost collinear to the channel vectors of the UEs selected in the previous k -1 stages.

According to (39), such a UE will be selected at the very end of the procedure. This, however, would result in a huge power consumption because the Gram-Schmidt procedure will only have a restricted number of degrees of freedom to make c

(K j ) j orthogonal to g ji for i = 1, . . . , K j -1
and at the same time to meet the requested SINR. In other words, the procedure in (39) sorts the UEs according to a best-first criterion such that higher priority is given to the UEs requiring low power to meet their SNR constraints.

Unlike [START_REF] Hunger | Design of single-group multicasting-beamformers[END_REF], we make use of the power increase (38) at each stage k to order the UEs within each group according to a worst-first criterion. As mentioned before, this choice is motivated by the fact that the Gram-Schmidt procedure in [START_REF] Andrews | What will 5G be?[END_REF] progressively reduces the available degrees of freedom as k tends to K j . Therefore, since power consumption is dominated by UEs with the worst conditions (according to some criterion), higher priority should be given to these UEs.

Mathematically, we propose to compute the index µ jk at step k as follows:

µ jk = arg max i∈R (k-1) j |α (k) ji | 2 (40) with |α (k) ji | = -ρ (k) ji + ρ (k) ji 2 -g H ji d (k) j 2 g H ji c (k-1) j 2 -η ji g H ji d (k) j 2 (41)
and ρ

(k) ji = g H ji d (k) j c
(k-1) H j g ji . As seen, µ jk corresponds to the UE index in R (k-1) j for which the incremental power |α

(k)
ji | 2 at stage k takes the maximum value. Note that the computational cost of this operation is still O(K 2 j N ) flops as for [START_REF] Hunger | Design of single-group multicasting-beamformers[END_REF]. To further reduce the computational burden, we propose an alternative approach that exploits the inherent characteristic of large-scale antenna systems. As N is large, each user i ∈ R |, the right-hand-side of (41) reduces to

η ji g ji 2 .
This means that UEs in group j can be ordered by simply sorting the following ratios in a descending order:

η j1 g j1 2 , . . . , η jK j g jK j 2 ∀j ∈ G. ( 42
)
In other words, higher priority should be given to those UEs that have bad channel conditions compared to the target SNRs. In doing so, no greedy strategy is required for UE ordering, thereby reducing the total number of flops to O(K j N ). Based on the above discussion, a heuristic solution is proposed in Algorithm 3 for the inner layer. Numerical results are used in Section V to make comparisons among the above ordering policies in different settings. As it will be seen, the ordering policy of (42) largely outperforms the strategy proposed in [START_REF] Hunger | Design of single-group multicasting-beamformers[END_REF].

C. The Proposed Heuristic Inner Layer Precoder

Collecting the results achieved in Sections IV-A and IV-B, we present the following heuristic algorithm to design the inner layer precoder of Q(η). To emphasis on the simplicity of Algorithm 3 and to enable the reproducibility of our results, its MATLAB code is provided in [START_REF] Sadegi | meysamsadeghi/reducing-the-computational-complexity-of-multicasting-in-large-scale-antenna-systems[END_REF].

Algorithm 3 A heuristic algorithm of the inner layer for solving Q(η)

1: for j = 1 to G do 2:
Sort the UEs in group j in descending order based on { η ji g ji 2 } and label the list as {µ j1 , . . . , µ jK j }, respectively.

for k = 2 to K j do 6: if |g H jµ jk c (k-1) j | 2 < η jµ jk then 7:
Compute α (k) j through (37) and (38).

8:

Update c

(k) j = c (k-1) j + α (k) j d (k) j . 9:
end if

10:
end for 11: end for

The complexity of Algorithm 3 can be evaluated as follows. Observe that evaluating the terms {η ji / g ji 2 } for group j requires 4(N -τ j ) flops whereas sorting a list of size K j needs O(K j log(K j )) flops. Therefore, the flop counts for UE ordering in line 2 is 4K j (Nτ j ) + O(K j log(K j )). The Gram-Schmidt procedure of line 4 can be performed through the QR decomposition, which requires 8(N -τ j )K 2 j -8 3 K 3 j flops [START_REF] Arakawa | Computational workloads for commonly used signal processing kernels[END_REF]. The computation of c | 2 < η jk in line 7 avoids to waste power for those UEs whose requested SNR constraints are already met (more details on this are given in the Appendix). Lines 7 to 10 require O(N -τ j ) flops, and as the condition of line 7 is true at most K j -1 times, the flops required by lines 6 to 11 is O(K j (N -τ j )). Summing all the above terms together, the complexity of Algorithm 3 is found to be O(N ), thereby reducing the complexity of the inner layer precoder by a factor of N 2 .

Note that by jointly employing Proposition 1 and Algorithm 3, the approximated precoding vectors for F(η, P ) can be computed as c • j,BDZF-HEU (η, P ) =

P P BDZF-HEU (η) c j,BDZF-HEU (η) j ∈ G (43)
where {c j,BDZF-HEU (η)} and P BDZF-HEU (η) denotes the precoding vectors and the total power consumption as obtained with Algorithm 3. Therefore the precoding vectors for F(η, P ) are given by {F j c • j,BDZF-HEU (η, P )}.

V. NUMERICAL RESULTS

Monte Carlo simulations are used to assess the performance of the proposed algorithms and to make comparisons with existing alternatives. In particular, we consider the algorithm presented in [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF], which employs the SDR technique followed by a randomization and multicast multigroup power control (MMPC) policy. 3 Comparisons are also made with the asymptotic results of [START_REF] Sadeghi | Multi-cell multi-group massive MIMO multicasting: An asymptotic analysis[END_REF],

the FPP based algorithm presented in [START_REF] Christopoulos | Multicast multigroup beamforming for per-antenna power constrained large-scale arrays[END_REF], and the heuristic algorithms developed in [START_REF] Hunger | Design of single-group multicasting-beamformers[END_REF]. A single cell system with radius of 900 meters is considered with UEs being distributed uniformly and randomly in the cell excluding an inner circular area of radius 100 meters. For each value of N , the average values of power consumption or minimum SINR of the system are obtained from 100 different channel realizations and UE distributions. We assume (if not otherwise specified)

that there are G = 3 multicasting groups, each counting K j = 10 UEs (such that K = 30).

The channel vector g jk between UE k in group j and the BS is modeled as g jk = β jk h jk where h jk ∼ CN (0, I N ) represents the small scale fading and β jk accounts for the large scale attenuation given by β jk = -128.1 -37.6 log 10 d jk dB with d jk being the distance between the UE and the BS expressed in kilometers [START_REF]Physical layer aspects for e-utra[END_REF]. The noise power spectral density is assumed to be -174 dBm/Hz, and the channel bandwidth is 20 MHz [START_REF] Xiang | Massive MIMO multicasting in noncooperative multicell networks[END_REF]. All the simulations are performed on a 64-bit Linux operating system with Intel Xeon processor E5-1680 v3.

Fig. 1 compares the average power consumption of the ordering policies proposed in [START_REF] Hunger | Design of single-group multicasting-beamformers[END_REF] with those given by ( 40) and (42), for G = 1, K = 20 and η = 63, 127, and 255 (which correspond to a target rate for each UE of 6, 7 and 8 bit/s/Hz, respectively). The proposed ordering policies are seen to outperform the ordering of [START_REF] Hunger | Design of single-group multicasting-beamformers[END_REF]. Moreover, the simple ordering policy of (42) has even a slightly better performance than (40). Note that, as the ordering belongs to the heuristic inner layer of the proposed precoder and as the outer layer removes the effect of inter-group interference, the same conclusion holds for G > 1. Based on the above results, the simpler ordering policy presented in (42) will be used in the remainder of this section. SDR of [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF] QoS Alg. of [START_REF] Hunger | Design of single-group multicasting-beamformers[END_REF] BDZF + Alg. 3 (proposed) Alg. 1 (proposed) FPP of [START_REF] Christopoulos | Multicast multigroup beamforming for per-antenna power constrained large-scale arrays[END_REF] Lower Bound lower bound of the problem is also presented as a benchmark [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF]. Observe that, the proposed algorithms outperform the SDR-based solution in [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF] and the heuristic one in [START_REF] Hunger | Design of single-group multicasting-beamformers[END_REF], while they have nearly the same performance as [START_REF] Christopoulos | Multicast multigroup beamforming for per-antenna power constrained large-scale arrays[END_REF]. However, this is achieved at a much lower complexity and computational cost as detailed next. Note that for N ≥ 60 both algorithms are at most 6% away from the lower bound and this gap reduces as N grows large, while for SDR technique this gap is 87% and reduces slowly by adding more antennas. Similar to the results of Fig. 2, the proposed algorithms largely outperform [START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF] and [START_REF] Hunger | Design of single-group multicasting-beamformers[END_REF], while nearly having the same performance as [START_REF] Christopoulos | Multicast multigroup beamforming for per-antenna power constrained large-scale arrays[END_REF]. However, this is achieved for a computational cost that is significantly smaller than other algorithms as detailed next. Observe that Algorithm 2 is within 5% of the upper bound with just N = 50 antennas. Also, Algorithm 3 jointly with BDZF and Proposition 1 achieve the same target with N = 70 antennas.

To assess the computational complexity of the investigated algorithms more intuitively, beside the complexity analysis of Section III-D and Section IV-C, we also present the computation time required to approximately solve Q(η) and F(η, P ) versus N , in Table I For the QoS problem, the requested SINR by each user is assumed to be 255, and for the MMF problem the available power at the BS is considered to be 10 Watt. Unlike Figs. 2 and 3 that provide the average of minimum power consumption or the average of minimum SINR of the system, these figures provide a clear vision on the distribution of these quantities, for the existing and proposed algorithms. It is seen that, as we increase N from 40 to 80, the CDF curves of our proposed algorithms become closer to the optimal bound, and also improve significantly in terms of performance, thanks to the large number of antennas. As an example, for the QoS problem with N = 40, the power consumption is greater than 20 Watt 60% of the times, while with N = 80 our proposed algorithms always meet the requested SINRs with less than 17 Watt.

Also, for the MMF problem none of the algorithms can provide a minimum SINR bigger than 200 with N = 40, while for N = 80 our proposed algorithms can provide a minimum SINR (b) N = 80.

Fig. 5: CDF of minimum SINR (MMF problem) of the system with P = 10.

bigger than 200 in 80% of the times. Fig. 5 also contains the CDF of the asymptotic approach of [START_REF] Sadeghi | Multi-cell multi-group massive MIMO multicasting: An asymptotic analysis[END_REF]. Notice that the asymptotic approach can never provide an SINR which is bigger than 20 (or 100) with N = 40 (or N = 80) antennas. The insufficiency of the asymptotic approach is detailed in [START_REF] Sadeghi | Multi-cell multi-group massive MIMO multicasting: An asymptotic analysis[END_REF].

In Section III, we have elaborated the BDZF-SCA approach and proved the convergence of Algorithms 1 and 2, but we have not specified the number of iterations required by each algorithm to converge. at the kth iteration of either Algorithm 1 or 2 as ε(k), the convergence condition of Table II is

|ε(k+1)-ε(k)| ε(k)
< 10 -3 . Also, Fig. 6 illustrates ε(k) for both algorithms at each iteration index k for different number of antennas N . As it is seen, both algorithms converge in a few iterations for any value of N . So far, we have assumed that the BS has perfect knowledge of the channel vectors {g jk }.

Although some of the analysis can in principle be extended (to some extent) to the scenario with imperfect CSI, we believe that this is out of the scope of this work and thus it is left for the future. To partially fulfill this lack, we now investigate the impact of imperfect CSI on the performance of the proposed algorithms. A time-division-protocol (TDD) is employed such that channel estimation can be performed in the uplink on the basis of UE pilot signals and then used in the downlink. We assume that pilots of length τ p = K are used, with power equal to 1

Watt. The estimates of channel vectors {g jk } are computed at the BS using an MMSE estimator.

This yields [START_REF] Marzetta | Fundamentals of Massive MIMO[END_REF] ĝjk = √ τ p β jk σ 2 + τ p β jk √ τ p g jk + n ∀j, k

where n ∼ CN (0, σ 2 I N ) is the additive Gaussian noise. Fig. 7 reports the performance of the proposed algorithms for MMF when perfect and imperfect CSI (as given in (44)) is available at the BS. We assume that G = 4 and K j = 15 for j = 1, . . . , 4. As expected, imperfect CSI degrades the performance of the proposed algorithms. However, such a performance loss can be compensated by using more antennas at the BS. Quantitively speaking, with imperfect CSI N must be roughly increased by a factor of 15% -30 % compared to the perfect CSI case. The higher N , the larger the factor. For example, to achieve the same performance of the perfect CSI case with N = 90 and 120, then 110 and 150 antennas are respectively needed with imperfect CSI, corresponding to a 22% and 25% increase.

VI. CONCLUSIONS

Multicasting is an efficient technology to transmit distinct common data streams to multiple groups of users. The existing multicasting algorithms are computationally expensive, especially for the large-scale antenna systems. There has been a struggle to present algorithms that can reduce this complexity. However, these algorithms are still too complex to fit in the next generation of wireless systems that enjoy hundreds of antennas, or they have a very poor performance.

In this paper, we designed new algorithms, which are tailored for physical layer multicasting in large-scale antenna systems. These algorithms not only provide good performance, but also this good performance is achieved at a very low computational cost and complexity. Particularly, we used the large number of antennas to cancel the inter-group interference and simplify the QoS and MMF problems. This enabled us to present efficient algorithms for these two problems.

Moreover, we showed that the two problems can be solved simultaneously and at no extra cost, while the baseline method is to solve the MMF problem by iteratively solving the QoS problem.

APPENDIX A PROOF OF PROPOSITION 1

Starting with the power constraint of (13) and replacing c j with P P app (η) c j,app (η) we have which proves the feasibility of the proposed solution. For the achievable objective value of F(η, P ) (or F(η, P )) using { P P app (η) c j,app (η)}, we have

8GN 2 - 8 (

 28 G-1)KN flops. In large-scale antenna array systems, i.e., where N K, the overall complexity of the proposed algorithm is dominated by the second step and it is of O(N 3 ). Taking into account that the complexity of SDR based techniques is greater than O(N 6 )[START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF], a reduction by a factor of N 3 is achieved through Algorithms 1 and 2.

(k- 1

 1 ) j = {µ j1 , . . . , µ j(k-1) } the set of indices of the ordered UEs at step k -1 and call R (k-1) j

j | 2

 2 required at stage k to meet the SNR constraint of the selected UE. To see how this comes about, consider a generic UE i ∈ R (k-1) j such that at stage k the ratio |g H ji c(k-1) j

(k- 1 j | 2

 12 ) j can use the excess degree of freedom, provided by the large number of antennas, to chose d (k) j as collinear as possible to g ji while almost nulling the interference generated to the other UEs, i.e., |g H ji c with ||g ji || 2 and neglecting the term |g H ji c (k-1) j

  using the Gram-Schmidt procedure in[START_REF] Andrews | What will 5G be?[END_REF].

  -τ j + 1) flops. The condition |g H jµ jk c (k-1) j

Fig. 2 Fig. 1 :

 21 Fig.2depicts the average power consumption of the QoS problem versus the number of antennas N at the BS. We assume that η jk = 255 for all UEs (corresponding to 8 bit/s/Hz/UE), and it is chosen in agreement with 5G rate requirements[START_REF] Andrews | What will 5G be?[END_REF], but the conclusions generically hold for all other values of η. The performance of Algorithm 1, and the combination of BDZF and Algorithm 3 is compared to other existing algorithms. As the QoS problem is NP-hard, a

Fig. 2 :

 2 Fig. 2: Average power consumption of the QoS problem with η jk = 255 ∀j, k.

Fig. 3 :

 3 Fig. 3: Average minimum SINR of the MMF problem for P = 10 Watt.

Fig. 3

 3 Fig.3illustrates the average minimum SINR of the MMF problem versus N . The available power at the BS is considered to be 10 Watt. In this figure, the performance achieved by Algorithm 2, and the combination of BDZF, Algorithm 3, and Proposition 1 is compared to other existing algorithms. Also, the upper bound of the problem is depicted as a benchmark.

Figs. 4

 4 Figs.[START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF] and 5 present the cumulative distribution function (CDF) of the precoder power consumption and the minimum SINR of the system for the QoS and MMF problems, respectively.

Fig. 4 :

 4 Fig. 4: CDF of minimum power consumption (QoS problem) of the system.

Fig. 6 :

 6 Fig. 6: The convergence behavior of Algorithms 1 and 2 for different number of antennas N .

Fig. 7 :

 7 Fig. 7: Evaluating the impact of imperfect CSI at the BS for the MMF Problem with G = 4, K j = 15 ∀j and P = 10 Watt. Channel estimation is performed using an MMSE estimator that exploit UL pilot signals of length τ p = K and power equal to 1 Watt.

  . The table presents the average time (in seconds) required to solve the QoS and MMF problems. The second and third columns report the average time required by the SDR and the FPP algorithms to solve an instance of the QoS problem. The fourth and fifth columns present the average required time by the same algorithms to solve an instance of the MMF problem. Note the increase in time from the QoS problem to the MMF problem, as in the SDR and FPP algorithms the MMF is

TABLE I :

 I The average time (in seconds) required to solve QoS, MMF, or both of them. applying the QoS algorithm. The sixth and seventh columns of the table present the average time required to solve both QoS and MMF problems simultaneously using the proposed algorithms. Note that not only we solve both problems at the same time with good performance, but also the required time has reduced significantly. As an example, for N = 90, the combination of BDZF, Algorithm 3, and Proposition 1, solves both problems in less than 5 milliseconds, while the SDR and the FPP algorithms require 2050 or 865 seconds, respectively, just to solve the MMF problem. At the same time, as shown in Figs.2 and 3, the solution provided by joint application of the BDZF, Algorithm 3, and Proposition 1 is nearly as good as the solution achieved from the FPP Algorithm, and significantly outperforms the SDR Algorithm.

		QoS Problem	MMF Problem		QoS and MMF (both problems simultaneously)
	SDR -[4] FPP -[11] SDR -[4] FPP -[11] Alg. 1 + Proposition 1 BDZF + Alg.3 + Proposition 1
	N = 40 55	41	419	356	11.3	2.5 × 10 -3
	N = 50 67	51	579	450	11.6	2.8 × 10 -3
	N = 60 84	61	798	507	11.7	3.1 × 10 -3
	N = 70 110	75	1151	617	11.9	3.5 × 10 -3
	N = 80 146	87	1549	727	12.2	4.0 × 10 -3
	N = 90 182	107	2050	865	12.5	4.5 × 10 -3
	solved by iteratively

TABLE II :

 II Average number of iterations required by Algorithms 1 and 2.

	(G = 2, K = 10)	(G = 3, K = 10)	(G = 3, K = 15)
	Algorithm 2 Algorithm 1 Algorithm 2 Algorithm 1 Algorithm 2 Algorithm 1
	N = 40 15.12	16.56	15.38	16.90	15.95	17.00
	N = 50 15.14	16.61	15.50	16.94	16.12	17.11
	N = 60 15.21	16.74	15.54	16.97	16.23	17.14
	N = 70 15.40	16.88	15.60	17.04	16.46	17.40
	N = 80 15.44	16.94	15.60	17.08	16.52	17.51
	N = 90 15.53	17.10	15.62	17.12	16.68	17.56

  Denote λ := min j min k 1 η jk g H jk c j,app (η) 2 . As {c j,app (η)} is a set of the precoding vectors of P app (η) is an objective value of Q(η) that can be achieved by { 1 √ λ c j,app (η)}, it is bigger than or equal to the optimal objective value of Q(η), i.e.,1 | 2 we should minimize -B + √ B 2 -4AC. Note -B + √ B 2 -4AC always has a negative derivative with respect to B, hence its minimum is achieved for the maximum value of B. Denote Γ jk e j∠θ jk = ḡH jµ jk d ḡjµ jk , we have B = 2 Γ jk Re(e j(∠θ jk +∠α (k) j ) ), the maximum of which achieved if ∠α

	is fixed, to minimize |α	(k)
						(k) j c	(k-1)H j
						(k) j = -∠θ jk and |α	(k)
								2
	t • app (η, P ) = min j	min k	1 η jk	g H jk	P P app (η)	c j,app (η)	=	P P app (η)	min j	min k	1 η jk	g H jk c j,app (η)	2 .
														(46)
	Q(η), λ ≥ 1. Therefore
					t • app (η, P ) =		P λ P app (η) 1	≥	P P app (η)	.	(47)
	As 1 λ P P (η)	.		(48)

λ P app (η) ≥ P (η) , and we have P P app (η) ≤ t • app (η, P ) ≤ j j | is given as in (56).

For brevity, hereafter, we refer to the weighted MMF problem as MMF problem.

The SVD calculation of an m by n matrix requires 16m 2 n + 24mn 2 flops[START_REF] Arakawa | Computational workloads for commonly used signal processing kernels[END_REF].

For the randomization phase, 100 samples are generated using the Gaussian randomization method[START_REF] Karipidis | Quality of service and max-min fair transmit beamforming to multiple cochannel multicast groups[END_REF].

L. Sanguinetti and R. Couillet have been supported by the ERC Starting Grant 305123 MORE. L. Sanguinetti is also funded by the 5GIOTTO project from University of Pisa.

APPENDIX B

PROOF OF PROPOSITION 2 Starting with the SINR constraint of [START_REF] Christopoulos | Multicast multigroup beamforming for per-antenna power constrained large-scale arrays[END_REF] ≥ η jk (49) where in (a) is due to the fact that t • app (η, P ) is the minimum weighted SINR among all UEs. Therefore {

For the objective value we have

As {c • j,app (η, P )} is a set of precoding vectors of F(η, P ), λ ≤ 1. Therefore,

Since 1 λ t • app (η, P ) is an objective value of F(η, P ) achieved by { 1 √ λ c • j,app (η, P )}, it is less than or equal to the optimal objective value of F(η, P ), i.e., 1 λ t • app (η, P ) ≤ t • (η, P ). Therefore we have

where (a) is due to our update rule, z

j which is bigger than η jk due to [START_REF] Xiang | Massive MIMO multicasting in noncooperative multicell networks[END_REF] and [START_REF] Marks | A general inner approximation algorithm for nonconvex mathematical programs[END_REF]. Therefore c

). Hence the objective function of (i + 1)th iteration is less than or equal to the objective function of (i)th iteration. As the objective function is bounded from below, by successively solving the problem we achieve a non-increasing bounded sequence. Therefore the algorithm converges. Due to [START_REF] Xiang | Massive MIMO multicasting in noncooperative multicell networks[END_REF], any internal precoding vector c j that satisfies [START_REF] Marks | A general inner approximation algorithm for nonconvex mathematical programs[END_REF], will also satisfy [START_REF] Christopoulos | Multicast multigroup beamforming for per-antenna power constrained large-scale arrays[END_REF] and as a result, any answer to Q(η, z) is a feasible answer to Q(η) and therefore Q(η). Due to the update rule and the inner approximation in [START_REF] Marks | A general inner approximation algorithm for nonconvex mathematical programs[END_REF], the convergence point satisfies the KKT conditions for Q(η) as detailed in [START_REF] Beck | A sequential parametric convex approximation method with applications to nonconvex truss topology design problems[END_REF].

APPENDIX D

SOLUTION TO [START_REF] Marzetta | Fundamentals of Massive MIMO[END_REF] Hereby we prove the solution of the following problem, i.e. [START_REF] Marzetta | Fundamentals of Massive MIMO[END_REF],

is achieved by α

j ) where:

(k-1) H j g jµ jk . We start with the SNR constraint g H jµ jk c | 2 -η jµ jk as A, B, and C, respectively. The SNR constraint can be represented as

≥ η jµ jk , to minimize the power, no transmission shall be arranged for user µ jk , i.e., α

j = 0 , and the next user shall be served. Otherwise, C < 0. Now we transform (57) to an equality by introducing λ ≥ 0 as follows . Now as A


Marc Dubois
email: marc.f.p.dubois@gmail.com

A note on the normative content of the Atkinson inequality aversion parameter

Keywords:

The aim of the note is to parametrize necessary and sufficient conditions for the Atkinson social welfare functions to satisfy proportional transfer principles. These social welfare functions satisfy the proportional ex post transfer principle if, and only if, the inequality aversion parameter is no smaller than 1. Moreover, these social welfare functions satisfy the proportional transfer principle if, and only if, the inequality aversion parameter is greater than 2. The results are illustrated with French data.

Introduction

The [START_REF] Atkinson | On the Measurement of Inequality[END_REF] inequality indices are defined by a class of social welfare functions (SWFs) exhibiting a downside inequality aversion. Indeed, these SWFs satisfy both the Pigou-Dalton principle of transfer and [START_REF] Kolm | Unequal inequalities II[END_REF] diminishing transfers principle. The former requires that a transfer of a positive amount of utility from one individual to a worseoff one reduces inequality, thus increasing the social welfare. 1 The latter principle requires that "one values more such a transfer between persons with given [utility] difference if these [utilities] are lower than if they are higher" [Kolm (1976, p. 87)]. Such SWFs are well-known to contain a measure of the degree of inequality aversion (ε, hereafter). Kolm (1976, pp. 87-88) points out that the Atkinson SWFs satisfy both principles if, and only if, ε is positive. [START_REF] Roberts | Comparability and Social Choice Theory[END_REF] and [START_REF] Blackorby | Ratio-Scale and Translation-Scale Full Interpersonal Comparability without Domain Restrictions: Admissible Social-Evaluation Functions[END_REF] indicate that the Atkinson SWFs are the only way to make social welfare evaluations whenever the ratio-scale comparability of utility is allowed and inequality aversion is supported. This informational basis asserts that x% of individual i's utility is more or less valuable than y% of individual j's utility. Therefore, it is possible to link these SWFs with some principles involving transfers that are proportional to the benefits of the agents involved in the transfers.

For that purpose, [START_REF] Fleurbaey | Transfer principles and inequality aversion, with an application to optimal growth[END_REF] introduce the proportional transfer principle (PP) and the proportional ex post transfer principle (PEP). The former requires that the social welfare increases as the result of a transfer that reduces the (pre-transfer) utility of a better off by x% and increases that of a worse-off by x%. Its main advantage relies on the desirability of such a transfer even if there is some loss in the transferred benefit. The latter principle follows the same idea as the PP, however it entails a lower loss in the transferred benefit because the benefit and the gift are percentages of post-transfer utility levels in lieu of pre-transfer utility levels.

We aim to parameterize necessary and sufficient conditions for the Atkinson SWFs to satisfy Fleurbaey and Michel's (2001) principles. Such SWFs satisfy the PEP if, and only if, ε is not lower than 1. Moreover, they satisfy the PP if, and only if, ε is higher than 2. For ε > 0, the Atkinson SWFs satisfy the Pigou-Dalton principle of transfer and the diminishing transfers principle. For ε 1, in addition to those principles, the Atkinson SWFs endorse a transfer of a positive amount from an individual to a less well-off one to be desirable even if there is some loss in the transferred benefit. For ε > 2, they endorse a transfer of a positive amount from an individual to a less well-off one to be desirable even if there is a greater loss in the transferred benefit than for ε 1.

This note is structured as follows: Section 2 presents the framework and the results; Section 3 provides an example from the income distribution of Montpellier with data of the National Institute of Statistics and Economic Studies of France (INSEE).

Setup and results

Let R (resp. R ++ and R --) denote the set of real numbers (resp. positive real numbers and negative real numbers). The population N := {1, . . . , n} is assumed to be the same throughout the note. Consider the Atkinson SWF:

(A) W (u) =    n i=1 u 1-ε i 1-ε , if ε > 0 and ε = 1; n i=1 ln (u i ) , if ε = 1.
We interpret

u 1-ε i 1-ε as the ethical value of individual i's utility u i (ln(u i ) being the ethical value if ε = 1). This utility level is an element of the utility distribution u := (u 1 , . . . , u n ) ∈ R n ++ . 2
2 Usually the Atkinson SWFs are introduced as follows:

W A (u) =    1 n n i=1 u 1-ε i 1 1-ε if ε > 0 and ε = 1; 1 n n i=1 ln (u i ) if ε = 1. (1)
Such SWFs are used for making a ranking of a set of utility distributions. For that purpose, consider u, v ∈ R n ++ :

W A (u) > W A (v) ⇐⇒ ln 1 n n i=1 u 1-ε i 1 1-ε > ln 1 n n i=1 v 1-ε i 1 1-ε ⇐⇒ n i=1 u 1-ε i 1 -ε > n i=1 v 1-ε i 1 -ε .
As in [START_REF] Atkinson | On the Measurement of Inequality[END_REF], ε is a parameter displaying some degree of inequality aversion. This parameter is also relevant with the Atkinson inequality index:

(I A ) I(u) = 1 - W (u) µ(u) ,
where µ(u) is the mean utility of the distribution u. The use of I A (or directly of a SWF (A)) implies ratio-scale comparability of utility. Thereby, it is possible to introduce transfers of proportions of utility. In this sense, Fleurbaey and Michel (2001) state the proportional ex post transfer principle.

Definition 2.1. Proportional ex post transfer principle [PEP]. For u, v ∈ R n ++ .
Consider that u is obtained from v by a proportional ex post transfer, i.e.: for some i, j ∈ N and δ > 0,

u i = v i + δu i u j = v j -δu j and for all m ∈ N \ {i, j}, u m = v m . Then, W (u) > W (v).
The PEP is clearly stronger than the Pigou-Dalton principle of transfer. The following proposition examines the constraints that the PEP imposes on the SWFs (A).

Proposition 2.1. For all W defined in (A), the following statements are equivalent: (i) W satisfies the proportional ex post transfer principle; (ii) ε 1.

Proof. Fleurbaey and Michel (2001, pp. 7-8) show that a SWF such as n i m(u i ) satisfies the PEP if, and only if, there is a concave function h : R --→ R such that m(u) = h (ln(u)). Then, h(τ ) = m (e τ ) with τ = ln(u). Let us consider a SWF defined in (A) and ε = 1, so:

h(τ ) = (e τ ) 1-ε 1 -ε = e τ [1-ε] 1 -ε .
The h function is concave if, and only if:

∂ 2 e τ [1-ε] 1-ε ∂τ 2 0 ⇐⇒ ∂ (1-ε)e τ [1-ε] 1-ε ∂τ 0 ⇐⇒ (1-ε)e τ [1-ε] 0 ⇐⇒ 1-ε 0 ⇐⇒ ε 1.
Moreover, for all u i ∈ R ++ , the ethical value of u i is ln (u i ) whenever ε = 1. This case implies that h (τ ) = 0, so h is concave. ++ . Consider that u is obtained from v by a proportional transfer, i.e.: for some i, j ∈ N and δ > 0,

u i = v i (1 + δ) u j = v j (1 -δ) and for all m ∈ N \ {i, j}, u m = v m . Then, W (u) > W (v).
That is why, for instance, [START_REF] Roberts | Comparability and Social Choice Theory[END_REF]

introduces the Atkinson SWF as in (A).

The following proposition examines the constraints that the PP imposes on the SWFs (A). Proposition 2.2. For all W defined in (A), the following statements are equivalent: (i) W satisfies the proportional transfer principle; (ii) ε > 2.

Proof. Fleurbaey et Michel (2001, pp. 5-7) show that a SWF such as n i m(u i ) satisfies the PP if, and only if, there is a function h : R --→ R strictly concave such that m(u

) = h -1 u . Then, h(τ ) = m -1 τ with τ = -1 u i .
Let us consider a SWF defined in (A) and ε = 1, so:

h(τ ) = -1 τ 1-ε 1 -ε ⇐⇒ h(τ ) = (-τ -1 ) 1-ε 1 -ε ⇐⇒ h(τ ) = -τ ε-1 1 -ε .
The h function is strictly concave if, and only if:

∂ 2 -τ ε-1 1-ε ∂τ 2 < 0 ⇐⇒ ∂ -(1-ε)(-τ ) ε-2 1-ε ∂τ < 0 ⇐⇒ ∂τ ε-2 ∂τ < 0 ⇐⇒ (ε -2)τ ε-3 < 0. Given that -1 u i ∈ R --, then τ ∈ R --, we obtain: ε > 2.
According to [START_REF] Fleurbaey | Transfer principles and inequality aversion, with an application to optimal growth[END_REF], the upper bounds on the loss, relative to the donor's gift in the proportional transfer and the proportional ex post transfer are respectively:

1 - v i v j and 1 - u i u j .
From Definitions 2.1 and 2.2, u i and u j are the donor's and the receiver's post-transfer utilities, respectively. Then:

1 - v i v j > 1 - u i u j .
In other words, SWFs that satisfy the PP can support a higher loss in the transferred benefit to reduce inequality than SWFs that purely satisfy the PEP.

An illustration with INSEE data

To illustrate the results shown in the previous section, we have taken the data collected by INSEE-DGFIP-CNAF-CCMSA 3 covering the income deciles (per consumption unit) for Montpellier in 2012. 4 In this section, income is considered as an indicator of household's utility. 5Table 1. Sensitivity of the Atkinson inequality index to ε ε-value 1.0 1.5 2.0 2.5 3.0 3.5 I A of Montpellier 0.09 0.14 0.19 0.23 0.27 0.31 Following [START_REF] Atkinson | On the Measurement of Inequality[END_REF], the inequality indices I A measure income inequality as a social cost. In Montpellier, for ε = 1, we have I A = 0.09, i.e. income inequality is a social cost valued as 9% of the total income. That means that if incomes were equally distributed, the same level of social welfare could be achieved with 91% of the total income. As presented in our results, the upper bound on the loss in the transferred benefit to equalize incomes is increasing over ε. For instance, for ε = 2, ε = 2.5 and ε = 3, if incomes were equally distributed, the same level of social welfare could be achieved with 81%, 77% and 73% of the total income, respectively. -0.020% +0.230% +0.432% +0.574% +0.652% +0.672% Montpellier (PEP4) +0.020% +0.268% +0.469% +0.609% +0.685% +0.703%

In Table 2, "Montpellier (PP4)"indicates an income distribution obtained from the original distribution of Montpellier by a proportional transfer of 4% in which donors are the households in the ninth decile and recipients are those in the second one. For instance, for ε = 3, the social welfare in Montpellier increases by 0.652% as the result of a transfer that reduces incomes in the ninth decile by 4% and increases incomes in the second one by 4%. Moreover, "Montpellier (PEP4)"indicates an income distribution obtained from the original distribution of Montpellier by a proportional ex post transfer of 4% involving the very same households as in "Montpellier (PP4)". For example, for ε = 1, the social welfare in Montpellier increases by 0.02% as the result of a transfer that reduces post-transfer incomes in the ninth decile by 4% and increases post-transfer incomes in the second one by 4%.

For ε 1, all SWFs (A) satisfy the PEP and, by definition, they yield Montpellier (PEP4) as a better distribution than that of Montpellier (original). However, not all SWFs (A) satisfy the PP, and some of them do not yield Montpellier (PP4) as a better distribution than the original one. Indeed, for ε = 1, the social welfare decreases by 0.02% as the result of a transfer that reduces incomes in the ninth decile by 4% and increases incomes in the second one by 4% (see the bold value in Table 2). Following the PEP, the loss in the transferred benefit must not exceed 64.1% of the donor's gift to consider "Montpellier (PP4)"as a better distribution than that of Montpellier (original) (see the bold value in Table 3). In fact, the proportional transfer of 4% involves a loss of 66.9% of the donor's gift. 6 For ε > 2, all SWFs (A) satisfy the PP (and so the PEP); they yield Montpellier (PEP4) and Montpellier (PP4) as better distributions than that of Montpellier (original).

Conclusion

We have provided parametric results derived from [START_REF] Fleurbaey | Transfer principles and inequality aversion, with an application to optimal growth[END_REF]. The Atkinson SWFs with a parameter ε 1 satisfy the proportional ex post transfer principle and such SWFs with ε > 2 satisfy the proportional transfer principle. Studying the example of Montpellier suggests that some Atkinson SWFs with ε 1 yield inequality-reducing reforms as social welfare decreases if such reforms entail higher losses in the transferred benefit than the upper bound determined by the proportional ex post transfer principle.

By analyzing the reforms implemented by a given social planner, it is possible to determine the upper bound on the loss in the transferred benefit he defends. This feature is helpful to estimate his degree of inequality aversion in terms of ε-value.

FleurbaeyDefinition 2 . 2 .

 22 and Michel (2001, p. 4) introduce another proportional transfer which is part of a stronger principle than the PEP. Proportional transfer principle [PP]. For u, v ∈ R n

Table 2 .

 2 Social welfare variations as the result of proportional transfers of 4%

	ε-value	1.0	1.5	2.0	2.5	3.0	3.5
	Montpellier (PP4)

Table 3 .

 3 Losses and upper bounds on the losses (relative to the donor gift)

	Montpellier (PP4) Montpellier (PEP4)

National Institute of Statistics and Economic Studies (INSEE), Public Finances Directorate General (DGFIP), National Family Allowances Office (CNAF), Agricultural Social Insurance Mutual Benefit Fund (CCMSA).

Located on the Mediterranean coast, Montpellier is among the fastest growing cities in France.

It may be worth noting that the following results are purely illustrative; the income deciles are used to make a simplified distribution in which the poorest 10% and the richest 10% of the population are discarded. An income distribution is composed by 8 income levels that are the midpoints of the extreme incomes of 8 groups formed by the deciles. For instance, the group with lowest income in Montpellier embodies 10% of the population with annual income of at least 5,

Euros and less than 10, 068 Euros in 2012. The simplification we make assumes that all households in this group earned 7, 883 euros in 2012 and so on for the other groups.


Via Enrico Fermi 2749, 21027 Ispra (VA), Italy, dejanbrkic0611@gmail.com, ORCID id: 0000-0002-

2502-0601

Analysis of few pipe networks is shown in the discussed paper where the authors conclude that the distribution of flow in a network of pipes with known topology, fixed pipe lengths and diameters, with the known and constant inputs and outputs assigned to nodes (connection between two pipes) for the chosen flow friction model (e.g. Colebrook's) depends among other on the chosen method for calculation of flow and pressure distribution in looped network of pipes. That conclusion is not sustainable. Also friction factors are used as constant which cannot be recommended.

Hydraulic analysis of the network from Figure 4 of the discussed paper (the network with five pipes and two loops) is repeated in this discussion.

Flow Friction

In the discussed paper flow friction factor is given with the constant value of f=0.02 in examples 1 and 2, and f=0.0242 in example 3 (Tables 1 and2 of the discussed paper). It is not acceptable. Flow friction factor is being changed during the calculation between f=0.02012 and f=0.02571. Flow friction is complex variable usually determined by empirical Colebrook's equation (1) or some of its related approximations [START_REF] Colebrook | Turbulent flow in pipes with particular reference to the transition region between the smooth and rough pipe laws[END_REF]Brkić 2011aBrkić , 2012a[START_REF] Brkić | Discussion of 'Exact analytical solutions of the Colebrook-White equation[END_REF]:

            D 3.71 ε f R 2.51 log 2 f 1 10 (1)
It is not obligatory to use this equation. Other equations can be used but f for sure cannot be treated as a constant. As the friction factor f is not a constant and even not always in turbulent zone for which Colebrook's equation ( 1) is only valid, it can be even suggested to use explicit approximation of the Colebrook equation (1a) but which is in addition valid also for laminar, and the transition in between laminar and turbulent zone (Swamee 1993):

                                          R R R f (1a)
The relation between pressure drop "Δp" and flow "Q" through pipe is Δp=r•Q 2 and with the simplification where "r" is with constant value (although it depends on variable "f"), it is similar to Ohm's law for electrical circuits which relates voltage "U" and electrical current "I" through constant thermal resistance "RT" as U=I•RT; note that R in Eqs. 1 and 1a is not thermal resistance but the Reynolds number. With linearization x=Q 2 , system of equations for pipe network in form Δp=r•x can be solved using non iterative methodology for electrical circuits. Knowing "x" for each pipe, the flow "Q" can be calculated easily. Such simplification would not produce accurate results (Liou 1998;Simpson and Elhay 2011;[START_REF] Brkić | Discussion of 'Jacobian matrix for solving water distribution system equations with the Darcy-Weisbach head-loss model' by Angus Simpson and Sylvan Elhay[END_REF][START_REF] Brkić | Discussion of 'Method to cope with zero flows in Newton solvers for water distribution systems[END_REF].

Hydraulic analysis of pipe network with loops -flow distribution through pipes for known inputs and outputs through nodes

The network used for analysis is shown in Figure 1. Pipe length, diameters and node consumption are constants. The first assumed flow pattern through pipes is chosen arbitrarily with only one condition to maintain flow balance in every node (material balance; first Kirchhoff law). The kinematic viscosity of water is prescribed as μ=1.0037·10 -6 m 2 /s, and absolute roughness of pipes is estimated at ε=0.00026m [START_REF] Brkić | Spreadsheet-based pipe networks analysis for teaching and learning purpose[END_REF]. Following the first Kirchhoff's law, system of equations for nodes can be established (2); where Q1-5 are flows through pipes and will be changed during calculation:

node node node node D C B A Q Q Q Q Q Q Q Q Q Q                      0 12 . 0 0 2 . 0 0 4 . 0 0 08 . 0 4 3 2 5 4 5 3 1 2 1 (2)
For two closed paths, following the second Kirchhoff's law, systems of equations for loops can be established (3); where Δp1-5 are pressure drops in pipes and their algebraic sum for each closed path in the network has to be approximately zero.

loop loop II I Q f Q f Q f p p p Q f Q f Q f p p p                                                                0 D L D L D L 8 0 D L D L D L 8 5 5 2 5 5 5 5 4 2 4 4 4 5 3 2 3 3 3 2 5 4 3 5 3 2 3 3 3 5 2 2 2 2 2 5 1 2 1 1 1 2 3 2 1     (3)
First derivative of Δp(Q); pressure drop "Δp" in function of flow "Q" as variable can be calculated as (4):

    Q r Q Q Q Q r Q Q Q f Q p F                                    2 D L 16 D L 8 ' 5 2 5 2      (4)
Original version of Hardy Cross method [START_REF] Cross | Analysis of flow in networks of conduits or conductors[END_REF] in matrix form, accelerated version of Hardy Cross method also known as Newton-Raphson method modified by Epp and Fowler (1971) and node-loop method by Wood and Charles (1972) improved by Wood and Rayes (1981) will be used.

Accelerated Hardy Cross method

Matrix

 

               ' 5 ' 4 ' 3 ' 3 ' 3 ' 3 ' 2 ' 1 ' F F F F F F F F F (5)
Correction of flow "ΔQ" will be calculated in each iteration from [F']x[ΔQ]=[Δp], where [F'] is from (5), and where [Δp] is based on (3). For the presented network it will be (6):

                                       II I Q Q x F F F F F F F F p p p p p p ' 5 ' 4 ' 3 ' 3 ' 3 ' 3 ' 2 ' 1 5 4 3 3 2 1 (6)
ΔQI and ΔQII are correction of flow which will respectfully will receive each pipe in loop I and loop II where pipe 3 common to both loops will receive both correction simultaneously following algebraic rules for that [START_REF] Corfield | Distribution design for increased demand[END_REF][START_REF] Brkić | An improvement of Hardy Cross method applied on looped spatial natural gas distribution networks[END_REF].

Original Hardy Cross method in matrix form

The same procedure as for the accelerated Hardy Cross method will be used with only difference that term ' 3 F  in matrix (5) will be zero. With this term equalized with zero the convergence toward the balanced solution is much slower.

The node-loop method

To avoid calculation of flow "Q" through "ΔQ", Wood and Charles (1972) and Wood and Rayes (1981) introduced the node-loop method. Flow Q can be calculated using [Q]=inv [NL]x[V] where matrix [NL] and [V] are defined as ( 7) and ( 8) respectively as described in [START_REF] Brkić | Iterative methods for looped network pipeline calculation[END_REF][START_REF] Brkić | Spreadsheet-based pipe networks analysis for teaching and learning purpose[END_REF]. To preserve linear independency among rows of [NL], one arbitrarily chosen node has to be omitted (in

this case node D).   node node loop loop loop II I C B A F F F F F F NL                                 ' 5 ' 4 ' 3 ' 3 ' 2 ' 1 0 0 0 0 1 1 0 0 0 1 0 1 0 1 0 0 0 1 1 (7)           loop loop output node input node output node II I C B A Q F Q F Q F p p p Q F Q F Q F p p p V                                                       5 ' 5 4 ' 4 3 ' 3 5 4 3 3 ' 3 2 ' 2 1 ' 1 3 2 1 2 . 0 4 . 0 08 . 0 (8)
In [V], first three rows are with constant values during whole calculation (in every iteration).

Results

Tree methods were used for calculation, the original Hardy Cross, the accelerated Hardy Cross (Newton Raphson) and the node-loop. The final solution is reached after approximately 9 iterations using the accelerated Hardy Cross (Newton Raphson) and the node-loop and after approximately 28 iterations using original Hardy Cross method. The final results are identical and they are listed in Table 1: Table 1. Information about the network, the flows obtained, pressure drops and related head losses All calculations were performed in Microsoft Excel.

Conclusion

Calculation with flow friction factor treated as constant cannot be recommended. Also using the chosen flow friction model (e.g. Colebrook's), final flow distribution does not depend on the chosen method (e.g. Hardy Cross, the node-loop method, etc.) for one defined network with known and constant pipe diameters, lengths and roughness of inner pipe surface and with constant consumptions and inputs assigned to nodes (Gay and Middleton 1971).

Supplementary material:

Complete calculation using the node-loop method extracted from MS Excel is provided.

Notation

The following symbols are used in this discussion:

  Discussion of "Economics and Statistical Evaluations of Using Microsoft Excel Solver in Pipe Network Analysis" by I.A. Oke; A. Ismail; S. Lukman; S.O. Ojo; O.O. Adeosun; and M.O. Nwude, J. Pipeline Syst. Eng. Pract. 06016002; doi. 10.1061/(ASCE)PS.1949-1204.0000240 Dejan Brkić, PhD, Research Scientific Officer; European Commission, DG Joint Research Centre (JRC), Directorate C: Energy, Transport and Climate, Unit C3: Energy Security, Distribution and Markets,

Figure 1 .

 1 Figure 1. Network with five pipes and two loops (adapted from the discussed paper)

  of derivatives [F'] can be established as (5); where on the main diagonal, upper term represent sum of the absolute values of the derivatives from loop I while lower term represents sum of the absolute values of the derivatives from loop II. Terms ' 3 F  are from the pipe 3 because that pipe is common for both loops.

  R -Reynolds number (dimensionless) ε/D -Relative roughness of inner pipe surface (dimensionless) f -Darcy (Moody) flow friction factor (dimensionless) Q -volumetric flow (m 3 /sec) Δp -pressure drop (Pa) ρ -water density (kg/m 3 ) L -length pf pipe (m) D -diameter of pipe (m) π-Ludolph number, π≈3.1415 Disclaimer: Neither the European Commission nor any person acting on behalf of it is responsible for the use which might be made of this publication. Gay, B., and Middleton, P. (1971). "The solution of pipe network problems." Chem. Eng. Sci., 26(1), 109-123. doi:10.1016/0009-2509(71)86084-0 Liou, C. (1998). "Limitations and proper use of the Hazen-Williams equation." J. Hydraul. Eng., 124(9), 951-954. doi. 10.1061/(ASCE)0733-9429(1998)124:9(951) Simpson, A. and Elhay, S. (2011). "Jacobian matrix for solving water distribution system equations with the Darcy-Weisbach head-loss model." J. Hydraul. Eng. 137 (6), 696-700., doi:10.1061/(ASCE)HY.1943-7900.0000341 Swamee, P.K.(1993). "Design of a submarine oil pipeline." J. Transp. Engrg., ASCE, 119(1), 159-170., doi:10.1061/(ASCE)0733-947X(1993)119:1(159) Wood,D., and Charles, C. (1972). "Hydraulic network analysis using linear theory.", J. Hydraul. Div.ASCE, 98(HY7), 1157-1179. Wood, D.J., and Rayes, A.G. (1981). "Reliability of algorithms for pipe network analysis.", J. Hydraul. Div. ASCE, 107(HY10), 1145-1161.

Figure 1 DBFigure

 1 Figure 1 DB

Table 1 .

 1 Information about the network, the flows obtained, pressure drops and related head losses

	Pipe Number	Diameter D (m)	Length L (m)	r	Flow (m 3 •s -1 )	Pressure drop Δp (Pa)	K=r•ρ -1 •g -1	Head loss H (m) H=Δp•ρ -1 •g -1
	1	0.25	1000	16741273.11	0.170123	484523.91	1706.55	49.39
	2	0.20	1500	81276026.10	0.090123	660138.41	8285.01	67.29
	3	0.15	1800	442404818.12	0.050866	1144662.32	45097.33	116.68
	4	0.10	1000	2083940584.61	0.020989	918078.81	212430.23	93.58
	5	0.20	1200	64370456.71	0.179011	2062741.10	6561.71	210.26

Table 1 DB

 1 Click here to download TableTable 1 DB.docx

Dear Editor,

The required file is now attached as "supplemental material". It has been part of submission since the first version of manuscript.

Best regards, Dejan Brkić


Impact du rayonnement au sein de couches limites turbulentes. Etude réalisée à l'aide de simulations numériques directes couplées à une méthode de Monte-Carlo réciproque.


Simulation aux grandes échelles de flamme au cours de la phase d'allumage d'un foyer annulaire équipé de 16 injecteurs swirlés prémélangés. Calculs réalisés dans le cadre du projet PRACE SIMAC.!