Low-Complexity Robust MISO Downlink Precoder Design With Per-Antenna Power Constraints

Low-Complexity Robust MISO Downlink Precoder Design With Per-Antenna Power Constraints

Mostafa Medra  Timothy N. Davidson
This work was supported in part by the Natural Sciences and Engineering Research Council of Canada under grant RGPIN-2015-06631.A preliminary version of a portion of this work appears in Proc. 2015 IEEE Int. Wkshp Signal Process. Adv. Wireless Commun.The authors are with the Department of Electrical and Computer Engineering, McMaster University, Hamilton, Ontario, Canada. (E-mail: medramm@mcmaster.ca; davidson@mcmaster.ca).
Abstract

This paper considers the design of the beamformers for a multiple-input single-output (MISO) downlink system that seeks to mitigate the impact of the imperfections in the channel state information (CSI) that is available at the base station (BS). The goal of the design is to minimize the outage probability of specified signal-to-interference-and-noise ratio (SINR) targets, while satisfying per-antenna power constraints (PAPCs), and to do so at a low computational cost. Based on insights from the offset maximization technique for robust beamforming, and observations regarding the structure of the optimality conditions, low-complexity iterative algorithms that involve the evaluation of closed-form expressions are developed. To further reduce the computational cost, algorithms are developed for per-antenna power-constrained variants of the zero-forcing (ZF) and maximum ratio transmission (MRT) beamforming directions. In the MRT case, our low-complexity version for systems with a large number of antennas may be of independent interest. The proposed algorithms are extended to systems with both PAPCs and a total power constraint. Simulation results show that the proposed robust designs can provide substantial gains in the outage probability while satisfying the PAPCs.

Broadcast channel, downlink beamforming, robust precoding, outage, per-antenna power constraints, massive MIMO, zero-forcing, maximum ratio transmission.

I Introduction

The spatial multiplexing capabilities of base stations (BSs) with multiple antennas offer the potential for substantial gains in the quality of service (QoS) that can be offered to users in a downlink system; e.g., [1]. In particular, linear beamforming schemes have been developed to simultaneously serve multiple users at their requested signal-to-interference-and-noise ratio (SINR) targets [2, 3, 4, 5, 6]. However, the performance of those beamforming schemes can be quite sensitive to the accuracy of the channel state information (CSI) that is available at the BS. Since that information is typically obtained through estimation on the uplink (in time division duplexing, TDD, systems) or through estimation on the downlink and quantized feedback (in frequency division duplexing, FDD, systems), the CSI at the BS is inherently uncertain. That observation has spawned the development of a variety of design strategies that incorporate different models for the uncertainty into the design. One strategy is to require the requested SINR to be met for all channels that are within a specified distance of the BS’s model for the channel [7, 8, 9, 10, 11, 12]. However, in many scenarios that is a rather conservative approximation of the outage that occurs in practice. Furthermore, although this strategy, or a mild approximation thereof, often results in a convex optimization problem for finding the beamformers, the computational cost of solving those problems can be quite significant. Fortunately, different approaches to approximating the outage probability can yield alternative design strategies that provide excellent performance in practice, even when the uncertainties in the CSI are quite substantial, and do so in a computationally inexpensive way. One such strategy is the offset maximization algorithm [13], in which the beamformers are designed to maximize a carefully structured offset on the performance specification (see Section II-B).

The above-mentioned design strategies seek to jointly design the beamforming directions and the power allocated to each direction. However, significant reductions in the computational cost can be obtained by computing the beamforming directions using a (computationally cheap) conventional technique and then developing a robust power loading algorithm. The beamforming directions in this approach are typically chosen to be either the maximum ratio transmission (MRT) [14] or zero-forcing (ZF) directions [15]. For the case of additive Gaussian uncertainties in the BS’s CSI, single-integral expressions for the outage probability can be obtained [16] and an effective algorithm for finding the power loading that minimizes the power required to meet the specified outage constraint has been developed [17]. However, that algorithm is rather computationally expensive. In [18], insights from bounds on the cumulative distributive function were used to develop a new robust power loading technique that provides performance close to that of the optimal algorithm in [17], but has significantly lower computational cost.

The existing literature on robust downlink beamforming has tended to focus on designs that impose a constraint on the total power transmitted by the BS. In practice, each antenna will typically be driven by its own power amplifier, and hence the design ought to include constraints on the power transmitted from each antenna, as well as the total power. In the case of perfect CSI, a number of downlink beamforming algorithms that incorporate per-antenna power constraints (PAPCs) have been developed [19, 20, 21, 22, 23]. For robust beamforming designs that can be formulated as convex problems (e.g., [9, 8, 12]) and are solved using generic solvers, incorporating these additional constraints is quite straightforward. However, doing so increases the computational cost of what are, in comparison to the perfect CSI case, already quite expensive algorithms. The goal of this paper is to develop robust beamforming designs that incorporate PAPCs and have reasonable computational costs. Our technique is based on insights developed from the offset maximization approach to robust beamforming [13], a closely related power loading technique [18], and observations regarding the structure of the optimality conditions for the design problem. These observations enable us to develop a low-complexity dual update optimization startegy related to that in [23] that involves the evaluation of a sequence of closed-form expressions. After extending that algorithm to systems that have both PAPCs and a total power constraint, we make the observation that a large fraction of the computational cost arises from the design of the beamforming directions. To reduce that cost, we develop PAPCed variants of the ZF and MRT directions, and show how these can be incorporated into our design approach. Furthermore, we develop a low-complexity version of our PAPCed MRT beamforming algorithm for “massive MIMO” systems with a large number of antennas. As scaling techniques for large MRT beamformers have been recently proposed [24], that algorithm may be of independent interest.

Ii System model and design approach

We consider a narrowband multiple-input single-output (MISO) downlink in which an -antenna BS sends independent messages to single-antenna users. The transmitted signal at a given signalling instant is constructed using linear beamforming as , where is the power-normalized data symbol for user , and is the associated beamformer. In some settings we will refer to as the direction of the beamformer, and as the power allocated to that direction. That enables us to write

The received signal at user can be written as

(1)

where denotes the channel between the BS and receiver , and represents the additive zero-mean circular complex Gaussian noise at that user.

In the problems that we will consider, each user specifies the SINR that it will require in order to support the service that it desires. This constraint takes the form

(2)

where is the noise variance at receiver , and is the required SINR. We will find it convenient to rewrite that constraint as

where

(3)

If we denote the signal transmitted from antenna by , then the power constraint on the BS as a whole can be written as , where we have used the assumptions that the messages are independent and that the symbols are normalized. If we let denote the maximum power that can be transmitted from antenna , the PAPC can be written as , where denotes the (,)th entry of the given matrix.

In order for a BS to be able to evaluate whether a candidate set of beamformers satisfies the SINR constraints in (2), the BS must know each channel vector ; e.g.,[2]. However, typically the BS will only have an estimate of each channel, denoted . To incorporate the uncertainty in that channel estimate into the design, we will postulate a conditional distribution, , and convert the deterministic QoS constraint into the chance constraint , where is the required outage probability. In this paper, we will model the uncertainty additively; i.e.,

(4)

with having zero-mean and being independent of the channel and data. Our results will focus on the case where is a zero-mean circular Gaussian random variable of covariance . Among a number of scenarios, that model is appropriate in certain TDD systems in which channels are estimated during the uplink training phase.

Ii-a Design approach

With the uncertainty modeled as described above, one approach to the design of the downlink beamformers is to seek to minimize the probability of outage of the SINR targets, subject to a total power constraint and PAPCs; i.e.,

(5a)
s.t. (5b)
(5c)
(5d)

This problem is hard to solve even without the PAPCs. However, in the case that the PAPCs are omitted, the offset maximization algorithm [13] is a low-complexity algorithm that has been shown to provide good performance. The goal of this paper is to use insights from the development of the offset maximization approach to develop an effective low-complexity algorithm for the PAPCed case. One observation that we will use is that the performance of the offset maximization approach can be improved by applying the robust power loading algorithm in [18] to the beamforming directions generated by the offset maximization. Doing so reveals that robust beamformers can be obtained with a computational cost that is similar to that of beamformer design in the perfect CSI case. (Many existing approaches to robust beamforming are much more expensive than the perfect CSI case; e.g., [11, 12].) However, like the perfect CSI case, it is the computation of the directions that dominate the computational cost. Therefore, we also propose to apply the principles that underlie the power loading in [18] to beamforming directions that can be computed more efficiently, such as PAPCed variants, derived herein, of the classical ZF and MRT directions; see Sections IV, and V. In the latter case, a further approximation that is suitable for scenarios with a large number of antennas at the BS substantially reduces the computational cost, and has almost the same outage performance.

To lay the groundwork for the development of the proposed beamforming schemes, in the following subsections we briefly review the offset maximization approach to beamformer design under a total power constraint [13], and the low-complexity robust power loading technique for systems with a total power constraint that was developed in [18].

Ii-B Offset maximization beamforming directions

The offset maximization beamformers [13] can be found by solving the following problem:

(6a)
s.t. (6b)
(6c)

It is implicit in (6c) that this algorithm tries to find the largest noise-plus-interference power each user can endure, under the total power constraint. In [13] an efficient method to solve (6) was developed by considering the following problem, in which, for now, it is assumed that the optimal value for (6), , is known:

(7a)
s.t. (7b)

It can be shown [13] that the optimal value of the problem in (7) is , and that any set of beamformers that optimize (7) are also optimal for (6). Also, at optimality, all the constraints are satisfied with equality.

The advantage of the connection between problems (6) and (7) is that a highly efficient algorithm for the problem in (7) with (i.e., the perfect CSI case) was developed in [3]; see also [6]. That algorithm can be extended to jointly find the optimal beamformers and the optimal offset, , for the problem in (6). In particular, if we let denote the Lagrange multiplier for the SINR constraint in (7b), then from the KKT conditions of (7) we can find the offset maximization directions by solving the eigen problem

(8)

where the Lagrange multipliers must satisfy the fixed-point relation

(9)

Since (8) can be solved using a power method, the complexity of finding the directions is dominated by the matrix inversion in (9), which requires operations. Having found those directions, the offset maximization power loading and the optimal offset can be found by solving the linear equations that arise when the constraints in (6b) and (6c) hold with equality.

Ii-C Robust power loading

The offset maximization algorithm described above uses the same offset to increase the robustness of each user to channel uncertainty. The goal of robust power loading approach in [18] is to provide a computationally-efficient way to adapt the offset to the characteristics of each user’s channel. For an arbitrary set of beamforming directions , the generic power loading problem can be stated as

(10a)
s.t. (10b)
(10c)

The derivation of the algorithm developed in [18] for producing good solutions to (10) begins by observing the under the additive uncertainty model in (4), the probability that is equal to the probability that

(11)

If we assume that the norms of the errors are small, as they will need to be for reliable operation [25], then we can approximate the quadratic term by a Gaussian random variable of the same mean and variance. In that case, the distribution of becomes Gaussian. (Recall that we are focusing on the case where is Gaussian, with zero mean and of covariance ; cf. (4).) Under that approximation, if we design the power loading so that the mean, , of is a significant multiple of its standard deviation, , then that user will achieve a low outage probability. Indeed, we can choose a value for that multiple so that the target outage probability is guaranteed to be satisfied; see, e.g., [7]. We also note that the optimal solution of (10) has equal values for . If that were not the case, the user(s) with higher outage probability could be allocated more power and the other user(s) less, which would reduce the objective value, and thus contradict the assumed optimality. Therefore, it is natural to choose the same multiple, , for each user in the approximation of the outage constraint in (10c). The resulting approximation of the problem in (10) can be written as [18]

(12a)
s.t. (12b)
(12c)

From the definition of in (11) and the channel uncertainty model in (4), it can be shown that

(13)

which is linear in the design variables . (Recall from (3) that ) Similarly, we have that

(14)

The structure of the problem in (12) is such that the constraints hold with equality at optimality [18]. Since is not a linear function in , that results in a set of non-linear equations for the power loading. The following iterative linearization technique has been shown in [18] to be an effective way to obtain good solutions to (12):

  1. Initialize each .

  2. Find and by solving the set of linear equations that arise from equality in (12b) and (12c) for the current values of , where is defined in (13).

  3. Update each using (14).

  4. Return to (2) until a convergence criterion is satisfied.

We note that the matrix that relates to and in step 2 is constant, and, accordingly, we need only invert this matrix once [18]. In practice, this algorithm converges quickly with a high probability [18]. In [26], the performance of this algorithm was shown to provide very similar performance to the optimal power loading in [17], and at a cost that is dominated by the operations that result from the initial matrix inversion.

Iii Offset maximization designs with PAPCs

To simplify the development of the proposed robust beamforming technique, we will first consider the addition of PAPCs to the offset maximization problem in (6). We will then modify the resulting algorithm using insights from the above robust power loading algorithm.

When we add the PAPCs to the offset maximization problem in (6), the design problem becomes

(15a)
s.t. (15b)
(15c)
(15d)

Although the formulation in (15) is not convex, it can be transformed in a straightforward way into a second order cone program, using the technique that was used for the case of perfect CSI; cf. [4, 23]. While that formulation can be solved using a generic interior point method (e.g., [27]), such generic methods do not exploit the structure of the problem, and the development of tailored algorithms that do exploit the structure offers the potential for improved computational efficiency.

In the following subsections, we will first develop a low-complexity algorithm for the case where we have PAPCs only, with no total power constraint. Then we will tackle the general problem with both types of power constraints. The development will use insights from algorithms developed for the perfect CSI case [23] and insights from the robust power loading algorithm described in Section II-C.

Iii-a Dominant PAPCs

If , the total power constraint can never be active and the problem in (15) can be rewritten as

(16a)
s.t. (16b)
(16c)

Motivated by the way that a customized algorithm for (7) was adapted [13] to solve the problem in (6), we consider the following problem in which, for now, is presumed to be known,

(17a)
s.t. (17b)
(17c)

In the context of (17), the constant term in the objective is superfluous, but it will simplify the interpretation of the Lagrangian. Using arguments analogous to those in [13, 23], it can be shown that any set of beamformers that is optimal for (17) is also optimal for (16), and the optimal value of in (17) is one.

Now, let denote the dual variable of the th condition in (17b) and denote the dual variable of the th condition in (17c). Let us also define the diagonal matrix , such that . These definitions enable us to write the Lagrangian of the problem in (17) as

(18)

Using the notion of complementary slackness, since the optimal value of is one, at optimality we have that . Also, at optimality we have , with lying in the null space of this matrix. This can be simplified to show that and

(19)

where denotes the Moore-Penrose pseudo-inverse, should be in the same direction. Further simplifications show that the dual variable in (19) should satisfy the fixed point equation

(20)

From (20) we observe that if we were given the optimal , we could find the optimal values for using (20) and then the optimal directions by normalizing the obtained using (19). After doing so, we could complete the solution of (17) by finding the optimal values for . That can be done by solving the set of linear equations that arise from the fact that at optimality (17c) holds with equality. (If this were not the case for condition in (17c), then the amplitude of could be decreased which would allow a smaller value of while satisfy all the other constraints.) To adapt that approach to solve (16), in the final step we must simultaneously solve for and . To do so we observe that enters linearly into (17c), and hence all we need is one more linearly independent equation. To obtain that equation we observe that if , then the th component of (17b) holds with equality. By summing over all the active constraints in (17b) we obtain the following equation

(21)

In the case that all the are positive — a case that happens quite often — the equation in (21) simplifies to

To complete the algorithm, we need to develop a technique to determine the optimal . One strategy for doing so is to apply the projected subgradient technique developed in [23]. That involves applying the update equation , where proj() denotes the projection of a matrix on the space of diagonal positive semidefinite matrices that satisfy and, consistent with the syntax used in Matlab, when diag() operates on a matrix it produces a vector containing the diagonal elements and when it operates on a vector it produces a diagonal matrix with the elements of the vector on the diagonal. The initialization parameters used in [23] were chosen to be and the step size chosen to be . Although this strategy converges, it can be quite slow [23]. In this paper, we will refine the approach in two ways. First, in Appendix A we develop a computationally cheap quasi-closed-form expression for the projection of in a 2-norm sense. Second, based on insights from [28] we will choose a step size of the form , for some positive scalar . In addition, in Section III-C we will identify a prediction step that can be used in the first iteration to accelerate the algorithm. One simple termination strategy is to stop the algorithm when , where is the maximum allowable violation of the power constraint for the th antenna. Following the above development, the algorithm can be summarized as shown in Algorithm 1.

1:Initialize the diagonal matrix such that each element is non-negative and . Set .
2:while  for any  do
3:     Find using (20).
4:     Solve for the directions by normalizing the obtained using (19).
5:     Find the power loading and by solving the set of linear equations arising from (17c) holding with equality and (21).
6:     Update using the results in Appendix A.
7:     Increment .
8:end while
Algorithm 1 Offset maximization with PAPCs

Having developed an efficient algorithm for the offset maximization problem with PAPCs, we now seek to incorporate the principles of the robust power loading discussed in Section II-C. To do so, we note that in the offset maximization design, the directions are independent of the offset term in (16c); cf. (19) and (20). That suggests that we could simply modify the power loading step. Indeed, once the directions have been obtained in step 4 of Algorithm 1, we can replace the power loading in step 5 by the and that solve (12). Those values can be found using the algorithm in Section II-C; see [18]. Incorporating that robust power loading algorithm into the framework of Algorithm 1 results in Algorithm 2.

1:Initialize the diagonal matrix such that each element is non-negative and . Set .
2:while  for any  do
3:     Find using (20).
4:     Solve for the directions by normalizing obtained using (19).
5:     Find and by solving and (21) using the method provided in Section II-C.
6:     Update using the results in Appendix A.
7:     Increment .
8:end while
Algorithm 2 PAPCed offset maximization with robust power loading

Iii-B Total and PAPCed algorithm

Using the principles outlined in Section II-B and the previous section, we can develop an algorithm for solving the general problem in (15), which has PAPCs and a total power constraint. In this section, we will focus on the case when is sufficiently smaller than to ensure that the total power constraint is active. (Otherwise, the problem can be solved by the techniques in the previous section.) Similar to the previous section, we will obtain the beamforming directions by normalizing the beamformers resulting from the following problem

(22a)
s.t. (22b)
(22c)

and then we will refine the power loading using the method described in Section II-C. As in the previous development, the Lagrangian of (22) plays a key role. It can be written as

(23)

Using the KKT conditions, for a given value for we can compute the corresponding directions and then the robust power loading in Section II-C. Furthermore, the subgradient used in the previous section remains a subgradient in this case. However, the structure of the KKT conditions is simpler in this case, which results in a more straightforward projection for the matrix. Indeed, since the only constraint on in this case is that it is non-negative, the update equation for can be written as

(24)

where the maximum operator is defined element-wise, and is the vector whose th element is . Therefore, we can construct an algorithm that has a similar structure to that in Algorithm 2. Having said that, in the case of PAPCs only there is a strong likelehood that the PAPCs will be active at optimality, and hence it makes sense to initialize the algorithm with a positive definite matrix . In the general case, the PAPCs are less likely to be active at optimality, and hence we will initialize the algorithm with The resulting algorithm is provided in Algorithm 3.

1:Initialize . Set .
2:while  for any  do
3:     Find using the fixed point equations
4:
5:     Solve for the directions , where
6:.
7:     Find and by solving and using the method provided in Section II-C.
8:     Update using (24).
9:     Increment .
10:end while
Algorithm 3 Generalized offset maximization

Iii-C Algorithm acceleration

As will be apparent in the simulations in Section VI, the modified update in Appendix A and the improved step size selection result in a substantial reduction of the number of iterations required over the number required using the choices made in [23]. Furthermore, we have observed that and the corresponding matrix at the termination of the algorithm are typically closely related. If that relationship can be determined with reasonable accuracy, this observation suggests that a predictive step could be used to further reduce the number of iterations. As an example of what can be done, in Section VI we illustrate how replacing with a simple affine prediction, , of the terminating matrix results in substantial reduction in the number of iterations.

Iv Conventional ZF beamforming with per-antenna power constraints

Even though the computational cost of each iteration of the PAPCed offset maximization beamforming algorithms in the previous section is dominated by terms that are only , when the BS has a large number of antennas the resulting computational load can still be substantial. The dominating components arise from determining the beamforming directions, and the fact that these directions are updated at each iteration. That suggests that we may be able to develop lower cost algorithms for systems with a large number of antennas if we could find a way to simplify the computation of the beamforming directions. In this section we will do that by developing variants of the nominal ZF directions, and we will integrate them with the robust power loading technique while ensuring that the required PAPCs are satisfied. In the following section we will develop analogous techniques based on variants of the MRT directions. For the ZF case, the beamforming directions are obtained using techniques developed in [29], but in the MRT case, the design of the beamforming directions appears to be new.

To develop PAPCed variants of the conventional ZF and MRT beamformers, we observe that in contrast to QoS-based designs, in which the SINR is controlled directly (e.g., (15d)), the conventional ZF and MRT designs focus on the desired signal power and interference components of the SINR separately. In particular, given that the SINR for user is , if we were to maximize the minimum nominal received signal power subject to a total power constraint (i.e., subject to ) we would obtain beamformers that are a particular power loading of the nominal MRT directions. If we were to add the nominal ZF constraints on the interference into that problem (i.e., ), then we would obtain beamformers that are a particular power loading of the ZF directions [29]. Due to the structure of the total power constraint, in many simple beamforming problems the optimization of the beamforming directions decouples from the power loading. That is indeed the case for our formulation for MRT and ZF beamforming directions. As an example, if we were to maximize the minimum value of , which is the power of the signal transmitted in the direction of user , rather than the power received by that user, we would obtain a set of beamformers in the MRT or ZF directions, but with a different power loading.

When the total power constraint is replaced by PAPCs, the optimization of the beamforming directions becomes coupled with the power loading and hence the choice of the metric to optimize changes both the power loading and the directions. While our approach will work for either metric, and indeed for several others, we will focus on the second metric . The rationale for this choice is that while the received signal power is suitable for the ZF problem in the perfect CSI case, where the ZF constraints will eliminate the interference [29], it can be quite sensitive to the interference incurred due to channel estimation errors. (This is illustrated in our simulation results in Section VI.) Accordingly, we define the normalized channel directions and we formulate the following generic problem to obtain PAPCed versions of the conventional beamformers

(25a)
s.t. (25b)
(25c)
(25d)

The value of determines whether the problem is of the ZF type, the MRT type, or a variant thereof. When is negligible compared to the noise power, the formulation describes a ZF-based approach, and when is of the order of the noise power this represents a regularized ZF-based approach; cf. [15]. When is sufficiently large, the constraints in (25d) become inactive, and accordingly the formulation describes an MRT-based approach.

One strategy for solving (25) is to employ a semidefinite relaxation [30]. As in related beamforming methods based on semidefinite relaxation (e.g., [2]), that approach involves the solution of a convex optimization problem for a set of matrices and a post-processing step that extracts good beamformers from these matrices. However, the computational cost of solving the convex optimization problem is even higher than that of the offset maximization algorithm, and that is only the cost of determining the beamforming directions. Accordingly, in the following sections we will present low-cost algorithms for robust beamforming with PAPCed variants of the ZF and MRT beamforming directions.

Iv-a ZF beamforming with PAPCs only

When , the problem in (25) involves finding the beamforming vectors that remove the interference at the receivers under the nominal channel conditions and satisfy the PAPCs. The essence of this problem was addressed in [29] using a re-parametrization technique. In particular, let us define the matrix as the matrix whose th column is and the matrix . The th column of , denoted is a zero-forcing direction for the th user with a unit signal gain; i.e., . If we let denote a matrix whose columns form a basis for the null space of , then the set of all ZF directions for the th user is given by the th column of , for an arbitrary scaling matrix . Accordingly, the solution to the problem in (25) takes the form

(26)

where is the th column of matrix [29]. Note that the constraints in (25c) and (25d) (with ) are automatically satisfied by designing the precoding vectors in the form in (26). The conditions that remain to be met are the PAPCs, and that can be done by adjusting the scaling matrix . In [29], this problem was formulated as a convex quadratically-constrained program that can be efficiently solved

(27a)
s.t. (27b)

where is the th column of the identity matrix. To complete the design, we choose the largest value for such that , which means that the beamformers of the form in (26) satisfy the remaining constraints; i.e., those in (25b).

If we let the th entry of the diagonal matrix denote the dual variable of the th PAPC in (27b), then the KKT conditions of the dual problem of (27) show that the scaling matrix should satisfy . Although such a relation does not allow for a closed-form solution, as we do not know , it does allow for the integration of the robust power loading method in [18], as an alternative to giving all the users the same nominal signal strength . Furthermore, the explicit relation between and allows us to use the sub-gradient algorithm for and to calculate accordingly. The proposed algorithm is summarized as Algorithm 4.

1:Find and . Initialize . Set .
2:while  for any  do
3:     Compute
4:     Find the beamformers directions by normalizing .
5:     Find and by solving and (21) using the method provided in Section II-C.
6:     Update using the results in Appendix A.
7:     Increment .
8:end while
Algorithm 4 ZF with PAPCs and robust power loading

From a computational respective, the key steps in the initialization of this algorithm are the finding of the ZF directions and the null space of , which requires operations. Each iteration of the algorithm involves the iterative solution of the linear equations in step 5, which, as explained in Section II-C, requires operations, and the matrix operations required to update in step 3, which require operations. When the number of antennas is close to the number of users , the dimensions of the matrix are small, which means that in that case the computational cost of this algorithm is dominated by finding the ZF directions and the null space in the initialization step.

Iv-B Generalized ZF beamforming

The extension of the ZF design with PAPCs to accommodate a total power constraint is straightforward, and follows the same steps that were used in the generalized offset maximization problem; see Section III-B. The generalized ZF problem can be formulated by adding the total power constraint to the constraints in (25). Then we consider the equivalent power minimization problem, assuming, for now, that the optimal is known

(28a)
s.t. (28b)

Consistent with our previous analysis, we will let denote the diagonal matrix with the dual variables of the PAPCs on its diagonal. From the KKT conditions we can then show that Furthermore, as in the previous algorithm we replace the uniform power loading, , with the robust power loading from [18]. The resulting modified version of Algorithm 4 is stated in Algorithm 5. As is apparent from Algorithm 5, the order of its computational cost is the same as that of Algorithm 4.

1:Find and . Initialize . Set .
2:while  for any  do
3:     Compute
4:     Find the beamformers directions by normalizing .
5:     Find and by solving and using the method provided in Section II-C.
6:     Update using (24).
7:     Increment .
8:end while
Algorithm 5 Generalized ZF

V Conventional MRT with per-antenna power constraints

As we have seen in the previous section, our approach to imposing PAPCs on the class of ZF beamformers can result in an algorithm of lower computational cost than that of offset maximization with PAPCs. However, any advantage is dependent on the size of the null space of the channel matrix. In settings with a large number of antennas and a small number of users, such as those arise in massive MIMO, the size of the null space can be quite large. In this section, we will show how the complexity can be further reduced by using an MRT-based approach rather than the ZF-based approach.

V-a MRT with PAPCs

In the MRT case, the interference conditions are omitted from the problem in (25), and the problem of finding nominal MRT-based beamformers that satisfy PAPCs can be written as

(29a)
s.t. (29b)
(29c)

Following a similar analysis to those performed earlier, if we let denote the dual variable for the th PAPC, define the diagonal matrix such that , and define to be the dual variable for the th condition in (29c), then the Lagrangian of the problem in (29) can be written as:

(30)

Accordingly, we can state the KKT conditions in a simplified form as , , and The last condition can be re-written as which means that and have the same direction. We note that at optimality is equal to the th element of , scaled by then divided by the th element of . This equation does not allow any optimal to be zero except if the channel vector contains a zero, which, under most reasonable channel models, is a “zero-probability” event. Since each is positive, the constraints in (29b) are all active, and accordingly .

Similar to the analysis of the previous problems, if we know , then we can find the beamforming directions using and subsequently solve for the power loading using the linear equations that arise when (29c) holds with equality. The value of in (29c) can be calculated using the KKT equation . If the equations in (29c) were not satisfied with equality at optimality, we could rescale the beamforming vectors to get a larger value of , which would contradict the assumed optimality. This observation is similar to the observation in the offset maximization section that enabled the use of the subgradient algorithm to find . Accordingly, we can suggest the iterative algorithm in Algorithm 6.

1:Initialize . Set .
2:while  for any  do
3:     Solve for the directions using .
4:     Find the beamformer magnitudes and using the linear equations that arise when the constraints in (29c) are satisfied with equality and .
5:     Update using Appendix A.
6:     Increment .
7:end while
Algorithm 6 Nominal MRT with PAPCs

Algorithm 6 provides an iterative way to find the values of , and, accordingly, the optimal precoding vectors. Its complexity per iteration is no more than linear in . However, we will now develop a closed-form expression that approximates the optimal solution of Algorithm 6 when the PAPCs are the same; i.e., . This closed-form removes the need for any iterations, which allows for an algorithm that is suitable for massive MIMO settings. To develop the approximation, we first note that the PAPCs and the MRT constraints hold with equality at optimality. That means that at optimality