Multi-Group Multicast Beamforming by Superiorized Projections onto Convex Sets

In this paper, we propose an iterative algorithm to address the nonconvex multi-group multicast beamforming problem with quality-of-service constraints and per-antenna power constraints. We formulate a convex relaxation of the problem as a semidefinite program in a real Hilbert space, which allows us to approximate a point in the feasible set by iteratively applying a bounded perturbation resilient fixed-point mapping. Inspired by the superiorization methodology, we use this mapping as a basic algorithm, and we add in each iteration a small perturbation with the intent to reduce the objective value and the distance to nonconvex rank-constraint sets. We prove that the sequence of perturbations is bounded, so the algorithm is guaranteed to converge to a feasible point of the relaxed semidefinite program. Simulations show that the proposed approach outperforms existing algorithms in terms of both computation time and approximation gap in many cases.


I. INTRODUCTION
M ANY applications in wireless networks involve multicast communication, which can be defined as the transmission of identical information to multiple receivers.One example is connected driving, where applications such as platooning can benefit from transmitting the same status or control information to a group of vehicles [1].Another example is the transmission of audio signals for live events, where each spectator can select from a variety of audio streams.Both use cases can benefit considerably from physical layer precoders that ensure a given quality-of-service (QoS) level for the requested stream at each receiver while reusing the same time and frequency resources for all receivers.
Physical layer multicasting schemes have been extensively investigated in the last two decades.The authors of [2] show that the performance of multicast transmission can be greatly improved by exploiting channel state information (CSI) at the transmitter.They consider two beamforming problems for single-group multicast beamforming, the max-min-fair (MMF) multicast beamforming problem and the QoS constrained multicast beamforming problem.While the MMF formulation aims at maximizing the lowest signal-to-noise ratio (SNR) among a group of users subject to a unit power constraint on the beamforming vector, the objective of the QoS constrained formulation is to minimize the transmit power subject to SNR constraints for the individual users.Moreover, the authors of [2] show that the solutions to both problems are equivalent up to a scaling factor.
The more general case with multiple cochannel multicast groups is considered in [3].Unlike the single-group case, the QoS constrained and MMF versions of the multi-group multicast beamforming problem are different in the sense that a solution to one version cannot generally be obtained by scaling a solution to the other.However, algorithms for the QoS constrained formulation can be straightforwardly extended to approximate the MMF version, by performing a bisection search over the target signal-to-interference plus noise ratio (SINR) values.In this paper, we will therefore restrict our attention to the QoS constained formulation.
The QoS-constrained multi-group multicast beamforming problem is a well-studied nonconvex quadratically constrained quadratic programming (QCQP) problem, for which various algorithmic approximations have been proposed.Existing approaches such as semidefinite relaxation with Gaussian randomization and successive convex approximation (SCA) algorithms -also known as convex-concave-procedures (CCP) -involve solving a sequence of convex subproblems.Solutions to these subproblems can be approximated either using offthe-shelf interior-point methods or using first-order algorithms such as the alternating direction method of multipliers (ADMM).While the use of interior-point methods typically results in a high computational complexity, the ADMM can require a large number of iterations to achieve a certain accuracy.Regardless of the algorithm used to approximate each subproblem, the CCP results in nested approximation loops.Terminating the inner iteration after a finite number of steps can hinder the feasibiltiy of estimates, which is required to ensure that the CCP converges.By contrast, if we assume the singular value decomposition of a matrix to be computable, 1 the algorithm proposed in this paper is free of nested optimization loops.
In cases where constrained minimization becomes too costly, the superiorization methodology (see, e.g., [5], [6]) constitutes a promising alternative.Whereas the goal of constrained minimization is to find a feasible point (i.e., a point satisfying all constraints) for which the objective value is minimal, superiorization typically builds upon a simple fixed-point algorithm that produces a sequence of points which provably converges to a feasible point.This fixed-point algorithm serves as the so-called basic algorithm, which is then modified by adding small perturbations in each iteration with the intent to find a feasible point with reduced (not necessarily minimal) objective value.By showing that the basic algorithm is bounded perturbation resilient, its convergence guarantee towards a feasible point can be extended to this modified algorithm called a superiorized version of the basic algorithm.
In this paper, we consider the QoS-constrained multi-group multicast beamforming problem in [3] with optional perantenna power constraints as introduced in [7].We propose an algorithmic approximation based on superiorization of a bounded perturbation resilient fixed point mapping.To do so, we formulate the problem in a product Hilbert space composed of subspaces of Hermitian matrices.This allows us to approximate a feasible point of the relaxed problem with the well-known projections onto convex sets (POCS) algorithm [8], which iteratively applies a fixed-point mapping comprised of the (relaxed) projections onto each constraint set.We show that this operator is bounded perturbation resilient, which allows us to add small perturbations in each iteration with the intent to reduce the objective value and the distance to the nonconvex rank-one constraints.Simulations show that, compared to existing methods, the proposed approach can provide better approximations at a lower computational cost in many cases.

A. Preliminaries and Notation
Unless specified otherwise, lowercase letters denote scalars, lowercase letters in bold typeface denote vectors, uppercase letters in bold typeface denote matrices, and letters in calligraphic font denote sets.The sets of nonnegative integers, nonnegative real numbers, real numbers, and complex numbers are denoted by N, R + , R, and C, respectively.The real part, imaginary part, and complex conjugate of a complex number x ∈ C are denoted by Re{x}, Im{x}, and x * , respectively.The nonnegative part of a real number x ∈ R is denoted by (x) + := max(x, 0).
We denote by Id the identity operator and by I N the N ×Nidentity matrix.The all-zero vector or matrix is denoted by 0 and the ith Cartesian unit vector is denoted by e i , where the dimension of the space will be clear from the context.The Euclidean norm of a real or complex column vector x is denoted by x 2 =

√
x H x. The ith singular value of a matrix A ∈ C N ×N is denoted by σ i (A), where the singular values are ordered such that σ 1 (A) ≥ • • • ≥ σ N (A).For square matrices A we define diag(A) to be the column vector composed of the diagonal of A, and for row or column vectors a we define diag(a) to be a square diagonal matrix having a as its diagonal.We write A 0 for positive semidefinite (PSD) matrices A.
The distance between two points x, y ∈ H in a real Hilbert space (H, •, • ) is d(x, y) = x − y , where • is the norm induced by the inner product •, • .The distance between a point x ∈ H and a nonempty set C ⊂ H is defined as d(x, C) = inf y∈C x − y .Following [9], we define the projection of a point x ∈ H onto a nonempty subset C ⊂ H as the set and denote by P C : H → C an arbitrary but fixed selection of Π C , i.e., (∀x ∈ H) P C (x) ∈ Π C (x).If C is nonempty, closed, and convex, the set Π C (x) is a singleton for all x ∈ H, so Π C has a unique selection P C , which itself is called a projector.For closed nonconvex sets C = ∅ in finite-dimensional Hilbert spaces, Π C (x) is nonempty for all x ∈ H, although it is not generally a singleton.Nevertheless, we will refer to the selection P C as the projector, as the distinction from the setvalued operator Π C will always be clear.
A fixed point of a mapping T : For the following statements, let (H, •, • ) be a real Hilbert space with induced norm • .Fact 2. Let T : H → H be a nonexpansive mapping with Fix(T ) = ∅.Then for any initial point x 0 ∈ H and α ∈ (0, 1), the sequence (x n ) n∈N ⊂ H generated by converges weakly2 to an unspecified point in Fix(T ).This fact is a special case of [13,Proposition 17.10b].

II. PROBLEM STATEMENT
In Section II-A, we define the system model and state the multi-group multicast beamforming problem with QoS-and per-antenna-power-constraints, and we reformulate it in terms of a nonconvex semidefinite program (SDP).A well-known approach to approximating solutions to such problems resorts to solving a convex relaxation: First, the original problem is relaxed and solved using, e.g., interior point methods.Subsequently, randomization techniques are applied to obtain candidate solutions to the original problem [3], [14].However, in real-time applications, the complexity of interior point solvers becomes prohibitive as it grows very fast with the system size (i.e., the number of users and the number of antennas).
Therefore, in Section II-B, we formulate the problem in a real product Hilbert space composed of complex (Hermitian) matrices.This formulation makes the problem accessible by a variety of first-order algorithms with low complexity and provable convergence properties.

A. System Model and Original Problem
Following the system model in [3], we consider the downlink in a network with a transmitter equipped with N antenna elements, each of them represented by an element of the set N := {1, . . ., N }.Each user k ∈ K := {1, . . ., K} is equipped with a single receive antenna.The users are grouped into M disjoint multicast groups G m ⊆ K indexed by m ∈ M := {1, . . ., M }, such that M m=1 G m = K.Each member of a multicast group G m is intended to receive the same information-bearing symbol x m ∈ C. The receive signal for the kth user can be written as , where w m ∈ C N is the beamforming vector for the mth multicast group, h k ∈ C N is the instantaneous channel to user k, and n k ∈ C -drawn independently from the distribution CN (0, σ 2 k ) -is the noise sample at the receiver.Consequently, the transmit power for group G m is proportional to w m 2 2 .In this paper, we consider the multi-group multicast beamforming problem with QoS-constraints [3], which has the objective to minimize the total transmit power subject to constraints on the QoS expressed in terms of SINR requirements.We use the following problem formulation from [7], with an individual power-constraint for each transmit antenna: The objective function in (1a) corresponds to the total transmit power.The inequalities in (1b) constitute the SINR-constraints, where γ k is the SINR required by user k.The inequalities in (1c) correspond to the per-antenna power constraints, where e i ∈ R N is the ith Cartesian unit vector.The problem in (1) is a nonconvex QCQP, which is known to be NP-hard [2].A well-known strategy for approximating solutions to such problems is the semidefinite relaxation technique [3], [14].By this technique, we obtain a convex relaxation of the original problem by reformulating it as a nonconvex semidefinite program and by dropping the nonconvex rank constraints.More precisely, using the trace identity tr(AB) = tr(BA) for matrices A, B of compatible dimensions, we can write w m This formulation is equivalent to (1) 2) can be obtained by simply dropping the rank-constraints in (2e).The approach in [2], [3] solves this relaxed problem and, subsequently, generates candidate approximations for Problem (2) (and hence (1)) using randomization techniques.A solution to the relaxed problem is typically found using general-purpose interior point solvers, which results in high computational cost for large-scale problems.In the multi-group setting [3], each randomization step involves solving an additional power control problem, which further increases the computational burden.

B. Problem Formulation in a Real Hilbert Space
The objective of this section is to show that Problem (2) can be formulated in a real Hilbert space, which enables us to approach the problem by means of efficient projectionbased methods.To this end, we consider the real vector space V := C N ×N of complex N × N -matrices.More precisely, we define vector addition in the usual way, and we restrict scalar multiplication to real scalars a ∈ R, where each coefficient of a vector X ∈ V is multiplied by a to obtain the vector aX ∈ V.In this way, V is a real vector space, i.e., a vector space over the field R.
If we equip the space V with a real inner product3 which induces the standard Frobenius norm we obtain a real Hilbert space (V, •, • ).
In the remainder of this paper, we restrict our attention to the subspace H := {X ∈ V | X = X H } of Hermitian matrices.Following the notation in [8], we define a product space H M as the M -fold Cartesian product of H.In this vector space, the sum of two vectors X = (X 1 , . . ., X M ) and Y = (Y 1 , . . ., Y M ) ∈ H M is given by X + Y := (X 1 + Y 1 , . . ., X M + Y M ) and scalar multiplication is restricted to real scalars a ∈ R, where a (X 1 , . . ., X M ) := (aX 1 , . . ., aX M ).We equip the space H M with the inner product which induces the norm where • is also a real Hilbert space.
In order to pose Problem (2) in this Hilbert space, we express the objective function in (2a) and the constraints in (2b)-(2e) in terms of a convex function and closed sets in H M , •, • as shown below: 1) The objective function in (2a) can be written as the following inner product: where J = (I N , . . ., I N ).This follows from ( 3), ( 4), and the fact that to the closed half-space where .
Here, we introduced indices {g k } k∈K that assign to each receiver k ∈ K the multicast group G m to which it belongs (i.e., In order to verify that the set Q k in (6) indeed represents the SINR constraint for user k in (2b), we rearrange4 Using the definition of the inner product in (3), and the fact that which corresponds to the kth SINR constraint in (2b).
3) The per-antenna power constraints in (2c) are expressed by the closed convex set where This follows immediately from ( 3) and ( 4). 4) The PSD constraints in (2d) correspond to the closed convex cone C + given by 5) The rank constraints in (2e) can be represented by the nonconvex set Consequently, we can pose Problem (2) as minimize The problems in ( 2) and ( 9) are equivalent in the sense that {X m ∈ V} m∈M solves Problem (2) if and only if (X 1 , . . ., X M ) ∈ H M solves Problem (9).The advantage of the formulation in ( 9) is that it enables us to (i) streamline notation, (ii) express the updates of the algorithm proposed later in Section III in terms of well-known projections, and (iii) simplify proofs by using results in operator theory in Hilbert spaces, as we show in the following.
It is worth noting that all constraint sets described above are closed, so a projection onto each of the sets exists for any point X ∈ H M .This property is crucial to derive projection-based algorithms, such as the proposed algorithm.In particular, note that we cannot replace the inequality in (2e) with an equality, as commonly done in the literature.The reason is that, with an equality, the corresponding set is not closed, as shown in Remarks 1 and 2, and the practical implication is that the projection may not exist everywhere.Specifically, this happens whenever X = (X 1 , . . ., X M ) satisfies X m = 0 for some m ∈ M, which would leave the update rule at such points undefined in projection-based methods.This is illustrated for the case ).Since a sequence of zeros can only converge to zero, the singular value decomposition The above shows that R contains all its limit points, so it is closed.
is not a closed set, since for all X ∈ R and α ∈ (0, 1), the sequence , Z is any of the closest points of the set R to the zero vector 0.
Therefore, for any α ∈ (0, 1), αZ ∈ R and d(0, αZ) < d(0, Z), i.e., αZ ∈ R is closer to the zero vector than Z ∈ R , thus contradicting our assumption that Z is one of the closest points in R to the vector 0.

III. ALGORITHMIC SOLUTION
The main difficulty in solving (9) is the presence of the nonconvex rank constraint.A well-known technique for approximating rank-constrained semidefinite programs using convex optimization methods is the semidefinite relaxation approach [2], [3], [14].This approach first solves (9) without the rank constraint, and then it applies heuristics to obtain rank-one approximations based on the solution to this relaxed problem.Similarly, we can obtain a convex relaxation minimize of Problem ( 9) by dropping the nonconvex constraint set R. In principle, we could solve this relaxed problem using first-order techniques for constrained convex minimization.For instance, we could apply a projected (sub-)gradient method (see, e.g., [15, Section 3.2.3]),which interleaves (sub-)gradient steps for the objective function with projections onto the feasible set of Problem (10).However, computing the projection onto the intersection of all constraint sets in Problem (10) typically requires an inner optimization loop because no simple expression for this projection is known.As it was shown in [16], superiorization can significantly reduce the computation time compared to the projected gradient method in some applications if the projection onto the feasible set is difficult to compute.
The superiorization methodology typically relies on an iterative process that solves a convex feasibility problem (i.e., that produces a sequence of points converging to a point within the intersection of all constraint sets) by repeatedly applying a computationally simple mapping.This iterative algorithm is called the Basic Algorithm.Based on this Basic Algorithm, the superiorization methodology automatically produces a Superiorized Version of the Basic Algorithm, by adding bounded perturbations to the iterates of the Basic Algorithm in every iteration.
n∈N be a bounded sequence in a real Hilbert space and let The perturbations are typically generated based on subgradient steps for a given objective function, in a way that ensures the sequence of perturbations to be bounded.By showing that the Basic Algorithm is bounded perturbation resilient (i.e., that the resulting sequence is guaranteed to converge to a feasible point, even when bounded perturbations are added in each iteration), one can ensure that the sequence produced by the Superiorized Version of the Basic Algorithm also converges to a feasible point.In contrast to constrained minimization, superiorization does not guarantee that the objective value of the resulting approximation is minimal.However, the limit point of the superiorized algorithm typically has a lower objective value than the limit point of the unperturbed Basic Algorithm [6].
To apply the superiorization methodology to Problem (9), we proceed as follows.In Section III-A, we propose a Basic Algorithm by defining a mapping T : H M → H M .Given any point X (0) ∈ H M , this mapping generates a sequence of points converging to a feasible point of Problem (10) by In Section III-B, we define a sequence β (n) Y (n) n∈N of bounded perturbations, with the intent to reduce slightly (i) the objective value of Problem ( 9) and (ii) the distance to the nonconvex rank constraint R in every iteration.As we show in Proposition 2 below, the proposed perturbations can achieve both goals simultaneously.The sequence of perturbations yields a Superiorized Version of the Basic Algorithm given by In Section III-C, we prove that the algorithm in (12) converges to a feasible point of Problem (10) by showing that the mapping T is bounded perturbation resilient, and that n∈N is a sequence of bounded perturbations.The relation between the proposed method and the superiorization methodology is discussed in detail in Section III-D.Finally, the proposed algorithm is summarized in Section III-E.

A. Feasibility-Seeking Basic Algorithm
A feasible point for the relaxed SDP in (10) can be found by solving the convex feasibility problem According to Fact 2 and Definition 2, given any X (0) ∈ H M , the iteration in (11) generates a sequence of points converging to a point in C if T is α-averaged nonexpansive with Fix(T ) = C .A particular case of such a mapping, which is used in the well-known projections onto convex sets (POCS) algorithm [8], is given by (see also Fact 1) where for a nonempty closed convex set C ∈ H M , T µ C = Id + µ(P C − Id) denotes the relaxed projector onto C with relaxation parameter µ ∈ (0, 2).The formal expressions for the projections of X ∈ H M onto each of the sets in (13) are given below.
1) The SINR constraint sets Q k ∈ H M are half-spaces, the projections onto which are given by [11,Example 29.20] 2) The per-antenna power constraint set P is an intersection of the N half-spaces defined by the normal vectors D i in (7) for i ∈ N .Since these vectors are mutually orthogonal, i.e., (∀i ∈ N )(∀j ∈ N \{i}) D i , D j = 0, the projection onto P can be written in closed form as This follows from [8, Thm 4.3-1] and Halperin's Theorem (see [17], [18,Thm. 4.2]).
3) The set C + is the intersection of PSD cones in orthogonal subspaces of H M .The projection of X ∈ H M onto C + is therefore given component-wise by P C+ (X) = P H+ (X 1 ), . . ., P H+ (X M ) , where, n∈N of vectors X (n) ∈ H M produced by the update rule in ( 14) is guaranteed to converge to a solution of the feasibility problem in (13) for any X (0) ∈ H M , if a solution exists (i.e., if C = ∅).Note that this is the case if the relaxed semidefinite program in (10) is feasible.Alternatively, we can derive this convergence guarantee immediately from Remark 4 and Fact 2.

B. Proposed Perturbations
In the following, we devise perturbations that steer the iterates of the fixed point algorithm in (12) towards a solution to the nonconvex problem in (2) and (9).To do so, we introduce a mapping that reduces the objective value and a mapping that reduces the distance to rank constraint sets.Then we define the proposed perturbations based on the composition of these two mappings As proven in Proposition 2 below, the resulting perturbations can achieve both goals simultaneously.
1) Power Reduction by Bounded Perturbations: In the literature on superiorization, the perturbations are typically defined based on subgradient steps of the objective function (see, e.g., [6]).For the linear objective function in (10), this would result in perturbations of the form −α (I N , . . ., I N ) for some α > 0. These perturbations are problematic for the problem considered here because we are interested in solutions comprised of positive semidefinite rank-one matrices, and adding these perturbations to an iterate X = (X 1 , . . ., X M ) may result in indefinite full-rank component matrices X m − αI N .To avoid this problem, we introduce the function f 1 : 5 For the case of real symmetric matrices, see, e.g., [19,Lemma 2.1].The result in [19] is based on [20, Corollary 7.4.9.3], which assumes complex Hermitian matrices.The generalization of [19, Lemma 2.1] to complex Hermitian matrices is straightforward.
where • * is the nuclear norm.Since C ⊂ C + by ( 13), we have where λ i (X m ) and σ i (X m ) denote the ith eigenvalue and singular value of the mth component matrix of X, respectively.Hence we can write Therefore, by (5), minimizing f 1 over C is equivalent to minimizing the linear objective function in (9) (or ( 10)) over C , in the sense that the solution sets to both formulations are the same.As we will show below, this surrogate objective function gives rise to power-reducing perturbations, which are guaranteed not to increase the rank of their arguments' component matrices (see Remark 3).The power-reducing perturbations are designed according to two criteria.Firstly, they should decrease the value of the surrogate function f 1 .Secondly, they should not be too large in order to avoid slowing down convergence of the Basic Algorithm.For a given point X ∈ H M we derive a perturbation Y τ satisfying these two criteria by solving the problem (17) Here, |||Y||| 2 acts as a regularization on the perturbations' magnitude, and the parameter τ ≥ 0 balances the two design criteria.The next proposition shows that Y τ can be easily computed.
Proposition 1.The unique solution to (17) is given by where D τ : H → H is the singular value shrinkage operator [21] and Proof: Denote the perturbed point for a given choice of τ by Z τ := X+Y τ .By substituting Y = Z−X in (17), we can identify this point as Z τ = prox τ f1 (X), where the proximal mapping is given by prox τ f1 (X) ∈ arg min Note that the function is separable over m.Consequently, we can compute the proximal mapping in (20) by solving According to [21, Thm.2.1], the unique solution to (21) is given by Z τ | m = D τ (X m ). 6Substituting Y τ = Z τ − X yields (18), which is the desired result.

By defining (∀X ∈ H
we can express the power-reducing perturbation for a point X ∈ H M as Y = T α P (X) − X, where the mapping T α P := prox ασmax(X)f1 is given component-wise by (∀m ∈ M) Note that T 0 P (X) = X, and (∀α ≥ 1) T α P (X) = 0. Therefore, the magnitude of the power-reducing perturbations can be controlled by choosing the parameter α ∈ [0, 1].Moreover, in contrast to performing subgradient steps for the original cost function in (9), applying the perturbations in (23) cannot increase the rank: ).This follows immediately from (19).

2) Incorporating the Rank Constraints by Bounded
Perturbations: Next, we define perturbations that steer the iterate towards the rank constraint set R in (8).While objective functions used for superiorization are usually convex, the function f 2 : i.e., the distance to the set R, constitutes a nonconvex superiorization objective, so our approach does not follow exactly the superiorization methodology in [6] (but we can still prove convergence).
As the perturbations may steer the iterates away from the feasible set, their magnitude should not be unnecessarily large.Therefore, we choose the rank-reducing perturbations as P R (X) − X, where P R (X) ∈ Π R (X) denotes a (generalized) projection of a given point X ∈ H M onto the closed nonconvex set R. Since R is a closed set, the set-valued projection Π R (X) is nonempty for all X ∈ H M .A projection onto R can be computed by truncating all but the largest singular value of each component matrix to zero.We formally state this fact below.
3) Combining Power-and Rank Perturbations: Since both T α P in (23) and P R in (25) operate on the singular values of the component matrices, their composition is given by (∀m ∈ M) where, (∀m ∈ M) U m = [u m1 , . . ., u mN ] and V m = [v m1 , . . ., v mN ].Moreover, it is easy to verify that (∀X ∈ H M )(∀α ≥ 0), T α P P R (X) = P R T α P (X).We will now use the composition of T α P and P R to define a mapping Y α : Finally, we define the sequence n∈N of perturbations in (12) by where α (n) n∈N is a sequence in [0, 1] and β (n) n∈N is a summable sequence in [0, 1].The following proposition shows that the perturbations in (27) can simultaneously reduce the objective value and the distance to the rank constraint set. 1) The perturbations cannot increase the distance to the set C + , i.e., (∀X 3) If α > 0 and X ∈ C + , then the perturbations decrease the objective value of Problem (9), i.e., J, X+ λY α (X) < J, X whenever J, X > 0. 4) If λ > 0, the perturbations decrease the distance to the rank constraint set R. More precisely, ∀X ∈ H M f 2 (X + λY α (X)) < f 2 (X) whenever f 2 (X) > 0. Proof: 1) This is an immediate consequence of (26).
2) It follows from (19) and (23) . Moreover, by (25) we have that (∀λ 3) This result follows from 1) and 2), since With the perturbations defined in (27), the iteration in (12) yields the update rule of the proposed algorithm, where

C. Convergence of the Proposed Algorithm
We will now examine the convergence of the proposed algorithm in (28).For this purpose, let β (n)  n∈N be a summable sequence in [0, 1], let α (n)  n∈N be a sequence of nonnegative numbers, and denote by Y (n)  n∈N the sequence of perturbations according to (27).Then the sequence X (n)  n∈N produced by the algorithm in (28) converges to a feasible point of Problem (10) for all X (0) ∈ H M .To show this, we prove the following facts.
1) The mapping T in ( 14) is bounded perturbation resilient.
2) The sequence Y (n)  n∈N is bounded, such that n∈N is a sequence of bounded perturbations.
Consequently, the bounded perturbation resilience of T follows directly from [12,Thm. 3.1].We summarize this fact in the following Lemma.Lemma 1. [12] The algorithm in (12) is guaranteed to converge to a point in the solution set C of the feasibility problem in (13) n∈N is a sequence of bounded perturbations.
Proof: The authors of [12] have proved the bounded perturbation resilience of α-averaged nonexpansive mappings with nonempty fix-point set in a real Hilbert space.Consequently, this lemma follows from Remark 4 and [12, Thm.3.1].
2) Boundedness of the Perturbations: It remains to show that the sequence Y (n)  n∈N is bounded for all sequences α (n)  n∈N of nonnegative numbers and To this end, we note that (∀n ∈ N) Y (n) ≤ X (n) for any sequence α (n)  n∈N of nonnegative numbers: Lemma 2. The mapping Y α in (26) satisfies denote the singular value decomposition of the mth component matrix of X.According to (26), the mth component matrix of which concludes the proof.
The following known result, which is a special case of [11,Lemma 5.31], will be used in Lemma 3 to prove that the proposed perturbations are bounded: n∈N is a summable sequence in [0, 1] and that (∀n ∈ N) α (n) ≥ 0. Then the sequence of perturbations where (a) follows from the nonexpansivity of T , and (b) is a consequence of the triangle inequality.By Lemma 2, the perturbations defined in (27) satisfy (∀n ∈ N) Y (n) ≤ X (n) .Consequently, applying the triangle inequality again yields we can deduce from Fact 4 that the sequence a (n)  n∈N converges.This implies that there exists r ∈ R such that (∀n ∈ N) X (n) − Z ≤ r.
Consequently, we have where (a) follows from Lemma 2, (b) follows from the triangle inequality, and (c) follows from Fact 4.
Combining Lemmas 1 and 3 shows that the proposed algorithm converges to a feasible point of the relaxed semidefinite program in (10).This is summarized in the following proposition.
Proposition 3. The sequence produced by the algorithm in (12) with perturbations given by ( 27) is guaranteed to converge to a feasible point of Problem (10) n∈N is a summable sequence in [0, 1] and α (n)  n∈N is a sequence in R + .Proof: Follows immediately from Lemma 1 and Lemma 3.

D. Relation to the Superiorization Methodology
The authors of [6] define superiorization as follows: 'The superiorization methodology works by taking an iterative algorithm, investigating its perturbation resilience, and then, using proactively such permitted perturbations, forcing the perturbed algorithm to do something useful in addition to what it is originally designed to do.'Although our proposed algorithm matches this informal definition, there are some slight differences to the formal definition in [6], where the perturbations are required to be nonascending vectors for a convex superiorization objective function.
Definition 4 (Nonascending Vectors [6]).Given a function φ : R J → R and a point y ∈ R J , a vector d ∈ R J is said to be nonascending for φ at y iff d ≤ 1 and there is a δ > 0 such that for all λ ∈ [0, δ] we have φ(y + λd) ≤ φ(y).
In our case, the goal of superiorization is two-fold, in the sense that it is expressed by two separate functions f 1 : 15) is convex, the function f 2 in (24) (i.e., the distance to nonconvex rank constraint set R in ( 8)) is a nonconvex function.Moreover, we use perturbations that are not restricted to a unit ball, and therefore they are not necessarily nonascending vectors.However, as we have shown in Proposition 2, the proposed perturbations simultaneously reduce the values of f 1 and f 2 .Keeping these slight distinctions in mind, we will refer to the proposed algorithm in (12) as Superiorized Projections onto Convex Sets.

E. Summary of the Proposed Algorithm
The proposed multi-group multicast beamforming algorithm is summarized in Algorithm 1.It is defined by the relaxation parameters µ 1 , . . ., µ K+2 of the operator T in ( 14), a scalar a ∈ (0, 1) controlling the decay of the power-reducing perturbations, a scalar b ∈ (0, 1) controlling the decay of the sequence of perturbation scaling factors, i.e., (∀n ∈ N) α (n) = a n and β (n) = b n .The stopping criterion is based on a tolerance value > 0, and a maximum number n max of iterations.
The arguments of the algorithm are the indices g 1 , . . ., g K assigning a multicast group to each user, the channel vectors h 1 , . . ., h K ∈ C N , SINR requirements γ 1 , . . ., γ K , and noise powers σ 1 , . . ., σ K of all users as well as the per-antenna power constraints p 1 , . . ., p N .At each step, the algorithm computes a perturbation according to (26) and applies the feasibility seeking operator T in (14).It terminates when the relative variation of the estimate falls within the tolerance , or when the maximum number n max of iterations is reached.Finally, the beamforming vectors w = {w m } m∈M are computed by extracting the strongest principal component where Eq. ( 26) 7: Eq. ( 14) break 10: Eq. ( 30) IV. NUMERICAL RESULTS In this section, we compare Algorithm 1 (S-POCS) to several other methods from the literature.We choose identical noise levels and target SINRs for all users, i.e., (∀k ∈ K) σ k = σ and γ k = γ.For each problem instance, we generate K i.i.d.Rayleigh-fading channels In the first simulation, we drop the per-antenna power constraints, i.e., we set (∀i ∈ N ) p i = ∞, and we consider the following algorithms: • The proposed method summarized in Algorithm 1 (S-POCS) • Semidefinite relaxation with Gaussian randomization [3] (SDR-GauRan) • The successive convex approximation algorithm from [23], [24] (FPP-SCA) • The ADMM-based convex-concave procedure from [7] (CCP-ADMM) The S-POCS algorithm is as described in Algorithm 1, with parameters a = 0.95, b = 0.999, = 10 −6 , n max = 10 5 .For the QoS-constraint sets, we use relaxation parameters (∀k ∈ K) µ k = 1.9, and for the per-antenna power constraint set P and the PSD constraint C + , we use unrelaxed projections, i.e., µ K+2 = µ K+1 = 1.We initialize the S-POCS algorithm with X (0) = 0.The convex optimization problems in the SDR-GauRan and FPP-SCA algorithms are solved with the interior point solver SDPT3 [25].The parameters of the CCP-ADMM algorithm are as specified in [7].Achieving a fair comparison between these methods is difficult because the structure of the respective algorithms is quite different.
The SDR-GauRan algorithm begins by solving the relaxed problem in (10), and, subsequently, generates random candidate beamforming vectors using the RandA method [2], [3].In the multi-group setting, where M > 1, an additional convex optimization problem (multigroup multicast power control (MMPC), [3]) needs to be solved for each candidate vector.If no feasible MMPC problem is found during the RandA procedure, we define the output of the SDR-GauRan algorithm to be {ψ(X m )} m∈M , where X ∈ H M is a solution to the relaxed SDP in (10).
The FPP-SCA algorithm from [23] works by solving a sequence of convex subproblems.By introducing slack variables, the feasibility of each subproblem is ensured.This obviates the need for a feasible initialization point, which is typically required to ensure convergence of CCP/SCA algorithms.
The CCP-ADMM algorithm uses an ADMM algorithm to find a feasible starting point for the CCP.Subsequently, a similar ADMM algorithm is used to approximate each subproblem of the CCP.Because the ADMM is a first-order method, the performance of CCP-ADMM is heavily dependent on the stopping criteria of the inner ADMM algorithm.
By contrast, the S-POCS algorithm does not require an initialization phase, and it works by iteratively applying a sequence of operators, which can be computed in a fixed number of steps.Therefore, we compare the performance based on computation time.Although we exclude the time required for evaluating the performance, we note that the computation time required by each of the methods severely depends on the particular implementation.
The authors of [7] assess the performance of the considered algorithms based by comparing the transmit power achieved by the resulting beamformers.However, none of the methods considered here can guarantee feasibility of the beamforming vectors, when the algorithms are terminated after a finite number of iterations.Furthermore, in the multi-group case, it may not be possible to scale an arbitrary candidate beamformer w = {w m ∈ C N } m∈M such that it satisfies all constraints in Problem (1).In principle, we could evaluate the performance by observing both the objective value (i.e., the transmit power of the beamformers) and a measure of constraints violation such as the normalized proximity function used in [26].However, defining this measure of constraints violation is not straightforward, as the considered methods approach the problem in different spaces.Moreover, we are interested in expressing the quality of a beamforming vector by a single value to simplify the presentation.Therefore, we will compare the performance based on the minimal SINR achieved by the beamformer ρ(w) • w with The scaled vector ρ(w) • w satisfies all power constraints, and its total power is bounded by the optimal objective value P SDR of the relaxed SDP in (10).More compactly, given a candidate beamformer w = {w m ∈ C N } m∈M for Problem (1), we assess its performance based on the function7 Since P SDR is a lower bound on the objective value of the original problem in (1), it holds (∀{w m ∈ C N } m∈M ) that SINR min ρ (w) ≤ γ, where equality can only be achieved, if the relaxed problem in (10) has a solution composed of rankone matrices.

A. Performance vs. Computation Time
We will now examine how the performance metric in (31) evolves over time for beamforming vectors produced by the respective algorithms.Figure 1 shows the performance comparison for an exemplary scenario with N = 20 antennas, and K = 20 users split evenly into M = 2 groups, where σ = 1, γ = 1, and (∀i ∈ N ) p i = ∞.It can be seen that the S-POCS algorithm quickly converges to a point achieving an SINR close to the specified target value γ.The discontinuities in the SINR curve for the CCP-ADMM algorithm are due to the inner-and outer optimization loops.For the SDR-GauRan algorithm, the SINR increases whenever the randomization produces a beamformer with better performance than the previous one.The SINR of the FPP-SCA algorithm improves continuously, albeit more slowly than the S-POCS and CCP-ADMM algorithms.Next, we evaluate the performance over 100 randomly generated problems.Since the SINR does not increase monotonically for all of the methods considered, we assume that each algorithm can keep track of the best beamformer produced so far.In this way, the oscillations in the SINR metric for the CCP-ADMM algorithm do not have a negative impact on its average performance.
Figure 2 shows the performance of the beamforming vectors computed with the respective algorithms over time for a system with N = 20 transmit antennas, and K = 20 users split evenly into M = 2 multicast groups.The shaded regions correspond to the 100%, 75%, 50%, and 25% quantiles over all randomly generated problems.More precisely, the margins of the shaded regions correspond to the 1st, 13th, 26th, 38th, 63rd, 75th, 88th, and 100th out of 100 sorted y-axis values.For each algorithm, the median is represented by a bold line.The S-POCS algorithm achieves the highest median SINR, while requiring the lowest computation time among all methods considered.Moreover, it can be seen that the variation around this median value is less severe compared to the remaining approaches.Put differently, the time required for reaching a certain SINR varies much less severely for the S-POCS algorithm than for the remaining methods.This can be of particular interest in delay sensitive applications, where a beamforming vector for a given channel realization must be computed within a fixed time period.

B. Varying number of antennas
In this subsection, we investigate the impact of the transmit antenna array size N on the performance of the respective beamforming algorithms.To do so, we generate 100 random problem instances for each array size N with K = 20 users split evenly in to M = 2 multicast groups.We choose unit target SINR and unit noise power for all users, and unit perantenna power constraints, i.e., γ = 1, σ = 1 and (∀i ∈ N) p i = 1.For the SDR-GauRan algorithm, we generate 200 candidate beamforming vectors for each problem instance.We use the CCP-ADMM algorithm with parameters as specified in [7].Since the inner ADMM iteration converges slowly for some problem instances, we set the maximal number of steps of the ADMM to j max = 300.For the outer CCP loop, we use the stopping criteria specified in [7], i.e., we stop the algorithm once the relative decrease of the objective value is below 10 −3 or t max = 30 outer iterations are exceeded.For the FPP-SCA algorithm, we use a fixed number of 30 successive convex approximation steps.Figure 3 shows the performance metric in (31) for different numbers N of transmit antennas, averaged over 100 random problem instances each.For all N , S-POCS achieves highest value for SINR min ρ (•), followed by the FPP-SCA, CCP-ADMM, and SDR-GauRan algorithms.For N ≥ 80, the S-POCS algorithm achieves an SINR of SINR min ρ (w S-POCS ) ≥ −0.05 dB.By contrast, the remaining methods do not exceed SINR min ρ (w FPP-SCA ) = −0.12dB, SINR min ρ (w CCP-ADMM ) ≥ −0.15 dB ,SINR min ρ (w SDR-GauRan ) ≥ −1.18 dB, respectively.The corresponding average computation times are shown in Figure 4.The S-POCS algorithm requires 0.26 %-2.38 % of the computation time required by SDR-GauRan, 0.95 %-11.64 % of the computation time required by FPP-SCA, and 6.49 %-233.6 % of the computation time required by CCP-ADMM.For N ≥ 80, the computation time of S-POCS exceeds that of CCP-ADMM.

C. Varying number of users
In the following simulation, we fix an array size of N = 50 antenna elements, and we evaluate the performance of each method for K ∈ {4, 8, 16, 32, 48, 64} users split evenly into M = 4 multicast groups.Figure 5 shows the performance metric in (31) averaged over 100 random problem instances for each K.As before, we choose γ = 1, σ = 1, and (∀i ∈ N ) p i = 1.
While all algorithms achieve close to optimal performance for small numbers of users, the SINR in (31) decreases considerably faster for SDR-GauRan than for the remaining methods.For all values of K, S-POCS achieves the highest value for SINR min ρ (•) among all methods.The corresponding average computation times are shown in Figure 6.S-POCS requires 1.76 %-6.12 % of the computation time required by SDR-GauRan, 3.75 %-5.41 % of the computation time required by FPP-SCA, and 20.18 %-1626 % of the computation time required by CCP-ADMM.While the CCP-ADMM takes only a fraction of the time required by S-POCS for small K, it slows down considerably as K increases.For moderate and large numbers of users, S-POCS outperforms the remaining methods in terms of both approximation gap and computation time.

D. Varying Target SINR
In the following simulation, we evaluate the impact of the target SINR on the respective algorithms in a system with N = 30 antenna elements, K = 20 users split evenly into M = 2 multicast groups, and unit noise power σ = 1.Since the target SINR has a strong impact on the transmit power, we set (∀i ∈ N ) p i = ∞, to avoid generating infeasible instances of Problem (1). Figure 7 shows the performance metric in (31) achieved by each method for the respective target SINR.Except for the SDR-GauRan algorithm, which exhibits a gap of about 2 dB to the target SINR, all methods achieve close to optimal performance for each target SINR.Figure 8 shows the computation time required by each algorithm for varying target SINR γ.The average computation time of FPP-SCA is almost constant.For SDR-GauRan and CCP-ADMM, the computation decreases slightly with an increasing target SINR.While the proposed S-POCS algorithm converges quickly for low target SINR levels, its computation time exceeds that of the CCP-ADMM for target SINRs above 8 dB.This indicates that the best choice of first-order algorithms for multicast beamforming depends on the regime in which the system is operated.V. CONCLUSION In this paper, we proposed an algorithm for multi-group multicast beamforming with per-antenna power constraints.We showed that the sequence produced by this algorithm is guaranteed to converge to a feasible point of the relaxed semidefinite program, while the perturbations added in each iteration reduce the objective value and the distance to the nonconvex rank constraints.Numerical comparisons show that the proposed method outperforms state-of-the-art algorithms in terms of both approximation gap and computation time in many cases.Its advantage over existing algorithms is particularly pronounced in the low target SINR regime as well as for large numbers of receivers.This makes the proposed method particularly relevant for low-energy or massive access applications.
In comparison to other techniques, the computation time of the proposed method varies less severely across different problem instances of the same dimension.In communication systems, which are typically subject to strict latency constraints, the iteration can be terminated after a fixed number of steps without suffering severe performance loss.Moreover, the simple structure of the proposed method allows for a straightforward implementation in real-world systems.
The applicability of the proposed algorithm is not restricted to the multicast beamforming problem considered here.A slight modification of the rank-constraint naturally leads to an algorithm for the general rank multicast beamforming problem considered in [27].Future research could apply superiorized projections onto convex sets to other nonconvex QCQP problems such as MIMO detection or sensor network localization [14].VI.APPENDIX Remark 5.The function •, • defined in (3) is a real inner product.

Fig. 1 .
Fig. 1.SINR min ρ (w (t) ) over time in a system with N = 20 antennas and K = 20 users users split evenly into M = 2 multicast groups.

Fig. 2 .
Fig. 2. SINR minρ (w(t) ) over time in a system with N = 20 antennas and K = 20 users split evenly into M = 2 multicast groups.The shaded regions include the outcomes for 100%, 75%, 50%, and 25% out of 100 problem instances, respectively, and the bold line represents the median.

Fig. 4 .
Fig. 4. Computation time for K = 20 users split evenly into M = 2 groups for varying antenna array sizes N .

Fig. 5 .
Fig. 5. SINR min ρ (w) for a system with N = 50 transmit antennas and a varying number of users split evenly into M = 4 multicast groups.

Fig. 6 .
Fig. 6.Computation time for a system with N = 50 transmit antennas and a varying number of users split evenly into M = 4 multicast groups.

Fig. 8 .
Fig. 8. Computation time for a system with N = 30 transmit antennas and K = 20 users split evenly into M = 2 multicast groups.
and replacing the expression w m w H m by a positive semidefinite rank-one matrix X m ∈ C N ×N for all m ∈ M, we obtain the nonconvex semidefinite program