Massive MIMO Multicast Beamforming Via Accelerated Random Coordinate Descent

TL;DR

This paper introduces a novel approach using accelerated random coordinate descent combined with majorization minimization to efficiently solve the complex nonconvex multicast beamforming problem in massive MIMO systems, significantly reducing computation time.

Contribution

It presents a new method that adapts ARCD for nonconvex, nonsmooth, and nonseparable problems in massive MIMO beamforming, achieving lower computational complexity.

Findings

Reduces execution time by an order of magnitude.

Successfully applies ARCD to a nonconvex problem.

Outperforms state-of-the-art methods in efficiency.

Abstract

One key feature of massive multiple-input multiple-output systems is the large number of antennas and users. As a result, reducing the computational complexity of beamforming design becomes imperative. To this end, the goal of this paper is to achieve a lower complexity order than that of existing beamforming methods, via the parallel accelerated random coordinate descent (ARCD). However, it is known that ARCD is only applicable when the problem is convex, smooth, and separable. In contrast, the beamforming design problem is nonconvex, nonsmooth, and nonseparable. Despite these challenges, this paper shows that it is possible to incorporate ARCD for multicast beamforming by leveraging majorization minimization and strong duality. Numerical results show that the proposed method reduces the execution time by one order of magnitude compared to state-of-the-art methods.

Equations

$$\mathrm{P}:{\min }_{\begin{array}{c}\mathbf{v}\end{array}}||\mathbf{v}|{|}_2^{2}s.t.|{\mathbf{g}}_k^H\mathbf{v}{|}^2\ge \gamma ,\forall k=1,\cdots ,K,$$

$$\mathrm{P}[n]:{\mathbf{v}}^{[n+1]}={\mathrm{argmin}}_{\begin{array}{c}\mathbf{v}\end{array}}$$

$$\displaystyle-|\mathbf{g}^{H}_{k}\mathbf{v}^{[n]}|^{2}\geq\gamma,~{}~{}\forall k\Bigg{\}}.$$

$$\mathrm{D}:{\max }_{\begin{array}{c}\mathbf{q}⪰0\end{array}}-{{\rm Y}}^{[n]}\left(\mathbf{q}\right),$$

$${{\rm Y}}^{[n]}\left(\mathbf{q}\right)=$$

$$-\sum _{k=1}^K{q}_k\left(\gamma +|{\mathbf{g}}_k^H{\mathbf{v}}^{[n]}{|}^2\right).$$

$${\mathbf{v}}^{\ast }=\sum _{k=1}^K{q}_k^{\ast }\cdot {\mathbf{g}}_k{\mathbf{g}}_k^H{\mathbf{v}}^{[n]}.$$

$$\mathcal{L}(\mathbf{v},\mathbf{q})=||\mathbf{v}|{|}_2^2$$

$$+{\sum }_{k=1}^K{q}_k\left\{\gamma -2\mathrm{Re}\left[{\left({\mathbf{v}}^{[n]}\right)}^H{\mathbf{g}}_k{\mathbf{g}}_k^H\mathbf{v}\right]+|{\mathbf{g}}_k^H{\mathbf{v}}^{[n]}{|}^2\right\},$$

$${\max }_{\mathbf{q}⪰0}{\min }_{\mathbf{v}}\mathcal{L}\left(\mathbf{v},\mathbf{q}\right).$$

$${\min }_{\mathbf{v}}\mathcal{L}\left(\mathbf{v},\mathbf{q}\right)$$

$$\mathbf{F}=\left[{\mathbf{g}}_1{\mathbf{g}}_1^H{\mathbf{v}}^{[n]},\cdots ,{\mathbf{g}}_K{\mathbf{g}}_K^H{\mathbf{v}}^{[n]}\right]\in {\mathbb{C}}^{N\times K},$$

$$\mathbf{d}={\left[\gamma +|{\mathbf{g}}_1^H{\mathbf{v}}^{[n]}{|}^2,\cdots ,\gamma +|{\mathbf{g}}_K^H{\mathbf{v}}^{[n]}{|}^2\right]}^T\in {\mathbb{C}}^{K\times 1},$$

$$L_k$$

$$\displaystyle\Big{|}\left[\nabla_{\mathbf{q}}\Upsilon^{[n]}(\mathbf{q}+t\cdot\mathbf{e}_{k})\right]_{k}-\left[\nabla_{\mathbf{q}}\Upsilon^{[n]}(\mathbf{q})\right]_{k}\Big{|}\leq L_{k}|t|,$$

$${\nabla }_{\mathbf{q}}{{\rm Y}}^{[n]}\left(\mathbf{q}\right)=2\mathrm{Re}\left({\mathbf{F}}^H\mathbf{Fq}\right)-\mathbf{d},$$

$${\left[{\mathbf{q}}^{[m+1]}\right]}_l={\left({\left[{\mathbf{q}}^{[m]}\right]}_l-\frac{1}{L_l}{\left[{\nabla }_{\mathbf{q}}{{\rm Y}}^{[n]}({\mathbf{q}}^{[m]})\right]}_l\right)}^{+},$$

$${\left[{\nabla }_{\mathbf{q}}{{\rm Y}}^{[n]}({\mathbf{q}}^{[m]})\right]}_l=$$

$$|{\mathcal{Y}}^{[m]}|=Y,{\mathcal{Y}}^{[m]}\subseteq \{1,\cdots ,K\},$$

$$\mathrm{Pr}(1\in {\mathcal{Y}}^{[m]})=\cdots =\mathrm{Pr}(K\in {\mathcal{Y}}^{[m]}),$$

$$\displaystyle\left[\mathbf{q}^{[m+1]}\right]_{i}=\Bigg{\{}\left[\mathbf{q}^{[m]}\right]_{i}-\frac{1}{Kc^{[m]}L_{i}}$$

$$\displaystyle\times\left[\nabla_{\mathbf{q}}\Upsilon^{[n]}\left(\mathbf{q}^{[m]}+(c^{[m]})^{2}\mathbf{z}^{[m]}\right)\right]_{i}\Bigg{\}}^{+},~{}~{}\forall i\in\mathcal{Y}_{m},$$

$$c^{[m+1]}=$$

$${\left[{\mathbf{z}}^{[m+1]}\right]}_i=$$

$$\times \left({\left[{\mathbf{q}}^{[m+1]}\right]}_i-{\left[{\mathbf{q}}^{[m]}\right]}_i\right),\forall i\in {\mathcal{Y}}_m,$$

Full text

**Index Terms— ** Acceleration, beamforming, large-scale, massive MIMO, random coordinate descent.

1 Introduction

Massive multiple-input multiple-output (MIMO) facilitates the concentration of wireless beams towards target directions [1], and is a promising technology for 5G communication systems and beyond [2, 3]. On the other hand, in many emerging applications such as video streaming [4] and computation offloading [5], a large number of users could be interested in the same data, making massive MIMO multicast beamforming indispensible.

In the context of massive MIMO multicasting, a fundamental criterion for beamforming optimization is to minimize the power consumption subject to quality-of-service (QoS) constraints [6, 7, 8, 9]. However, due to high dimensionality (e.g., the number of antennas $N$ and the number of users $K$ can be in the range of hundreds or more [2, 3]), traditional semidefinite relaxation (SDR) becomes extremely time-consuming, since the complexity of SDR is at least $O(KN^{3.5})$ [6, 7, 10]. To reduce the computational complexity of beamforming optimization in large-scale settings, first-order methods (FOMs), which only involve the computation of gradients [11], are recently proposed. In particular, the alternating direction method of multipliers (ADMM) has been derived in [12, 13] by introducing slack variables and solving the augmented Lagrangian problem. Nonetheless, the per-iteration complexity of ADMM is still $O(KN)$, and a fundamental question is: can we achieve a lower per-iteration complexity for the multicast beamforming optimization problem?

It turns out that this is possible if we only update one coordinate in each iteration [14]. This leads to the coordinate descent method, which involves a per-iteration complexity of $O(K)$. But unfortunately, the naive way of cyclicly updating the coordinates may diverge [15]. Even if it converges, the convergence rate can be slow [14], thus offsetting the benefit brought by the low per-iteration complexity. In fact, this is the reason why coordinate descent method received less attention compared to its full-gradient counterpart.

However, there has been a revival of interest in coordinate descent method recently [16, 17, 18]. In particular, it is proved in [16] that if we randomly update one coordinate in each iteration, the resultant random coordinate descent (RCD) is guaranteed to converge with a rate of $O(1/m)$, where $m$ is the iteration counter. Moreover, by adopting coordinate-wise momentum updates, RCD can be further accelerated to a faster convergence rate of $O\left(1/m^{2}\right)$ [19]. This results in the parallel accelerated random coordinate descent (ARCD) method, which reduces the computation time by orders of magnitude compared to traditional FOMs in extensive applications, e.g., inverse problem [19], $l_{1}$-regularized least squares problem [19], supervised learning [20], etc.

Nonetheless, ARCD can only be used to solve convex and smooth problems with separable constraints [16, 17, 18, 19, 20]. In contrast, the multicast beamforming optimization problem is nonconvex with nonseparable constraints. Therefore, we cannot directly apply ARCD to this application. To this end, leveraging majorization minimization [21], the nonconvex problem is transformed into a sequence of convex problems but with nonseparable constraints. Furthermore, to resolve the coupling among different coordinates, the nonseparable problem is equivalently transformed into its Lagrangian dual counterpart, which is proved to be a coordinate-wise Lipschitz smooth problem with separable constraints, thus allowing ARCD to work on this dual problem. Numerical results are further presented to demonstrate the low complexity nature of the proposed method.

Notation. Italic letter $x$, small bold letter $\mathbf{x}$, and capital bold letter $\mathbf{X}$ represent scalar, vector, and matrix, respectively. The operator $[x]^{+}=\mathrm{max}(x,0)$, $\mathrm{Re}(x)$ takes the real part of $x$, and $|x|$ takes the modulus of $x$. The symbol $[\mathbf{x}]_{k}$ takes the $k^{\mathrm{th}}$ element of vector $\mathbf{x}$, $[\mathbf{A}]_{k,:}$ is a row vector taking the $k^{\mathrm{th}}$ row of matrix $\mathbf{A}$, and $[\mathbf{A}]_{k,k}$ takes the diagonal element at the $k^{\mathrm{th}}$ row and the $k^{\mathrm{th}}$ column of matrix $\mathbf{A}$. Finally, $\mathbb{E}(\cdot)$ represents the expectation of a random variable and $O(\cdot)$ represents the order of arithmetic operations.

2 Problem Formulation and Existing Methods

We consider a massive MIMO system consisting of a base station with $N$ antennas, and $K$ single-antenna users. In particular, the base station transmits a signal $s$ with $\mathbb{E}[|s|^{2}]=1$ to all the users through the beamforming vector $\mathbf{v}\in\mathbb{C}^{N\times 1}$ with power $||\mathbf{v}||^{2}$. Accordingly, the received signal $r_{k}\in\mathbb{C}$ at the user $k$ is $r_{k}=\mathbf{h}^{H}_{k}\mathbf{v}s+n_{k}$, where $\mathbf{h}^{H}_{k}\in\mathbb{C}^{1\times N}$ is the downlink channel vector from the base station to user $k$, and $n_{k}\in\mathbb{C}$ is the Gaussian noise at the $k^{\mathrm{th}}$ user with power $\sigma_{k}^{2}$. Based on the expression of $r_{k}$, the received SNR at user $k$ is $|\mathbf{g}^{H}_{k}\mathbf{v}|^{2}$, where $\mathbf{g}_{k}:=\mathbf{h}_{k}/\sigma_{k}$.

In multicast systems, our aim is to provide guaranteed SNR for all the users, while minimizing the total transmit power at base station: [6]:

[TABLE]

where $\gamma$ is the common SNR target. Problem $\mathrm{P}$ has been proved to be NP-hard in general [6, Claim 1]. To solve $\mathrm{P}$, a traditional way is to apply SDR for convexification [6]. However, since SDR needs to solve a semidefinite programming (SDP) problem with $K$ variables and one semidefinite constraint of dimension $N\times N$, SDR requires a complexity of $O\left(\sqrt{N}(K^{3}+K^{2}N^{2}+KN^{3})\right)$ [22], which is too demanding when $N$ or $K$ is large.

To reduce the computational complexity of beamforming optimization, the majorization minimization (MM) framework [23, 21, 13, 24] can be adopted to transform $\mathrm{P}$ into a sequence of surrogate problems, Then, an iterative algorithm can be obtained with the following update at the $n^{\mathrm{th}}$ iteration:

[TABLE]

It has been proved in [23] that the sequence $\{\mathbf{v}^{[0]},\mathbf{v}^{[1]},\cdots\}$ converges to a Kruash-Kuhn-Tucker solution to $\mathrm{P}$.

However, even capitalizing on the MM framework, $\mathrm{P}[n]$ is still large-scale, and the interior point method (IPM) adopted in [23] would lead to time-consuming computations since the complexity of IPM is $O\left(\sqrt{K}(N^{3}+2NK)\right)$ [22]. To overcome this challenge, ADMM has been proposed for solving $\mathrm{P}[n]$ [12, 13]. Such a method reformulates $\mathrm{P}[n]$ into an augmented Lagrangian problem and then uses hybrid gradient method to solve it. Therefore, its per-iteration complexity is only $O(KN)$.

3 ARCD: Accelerated Random Coordinate Descent

While the per-iteration complexity of ADMM is lower than that of SDR, a natural question is: can we achieve a lower complexity than that of ADMM for solving $\mathrm{P}[n]$? This section will show that it is possible under the framework of ARCD. To begin with, the following property is established.

Property 1**.**

Strong duality holds for $\mathrm{P}[n]$.

Proof.

To prove the strong duality of $\mathrm{P}[n]$, it suffices to show that $\mathrm{P}[n]$ is convex and satisfies Slater’s condition [25]. Since the objective $||\mathbf{v}||_{2}^{2}$ of $\mathrm{P}[n]$ is convex quadratic and the constraints of $\mathrm{P}[n]$ are linear, $\mathrm{P}[n]$ is convex. On the other hand, showing that $\mathrm{P}[n]$ satisfies Slater’s condition is equivalent to finding a feasible point of $\mathrm{P}[n]$ satisfying all the constraints with strict inequality. To this end, consider the solution $\mathbf{v}=(1+\chi)\cdot\mathbf{v}^{[n]}$ with $\chi>0$, and it can be shown that $2\mathrm{Re}\left[\left(\mathbf{v}^{[n]}\right)^{H}\mathbf{g}_{k}\mathbf{g}^{H}_{k}\mathbf{v}^{\prime}\right]-|\mathbf{g}^{H}_{k}\mathbf{v}^{[n]}|^{2}>\gamma$ for all $k$. This completes the proof. ∎

Based on the result of Property 1, the dual problem of $\mathrm{P}[n]$ must have the same optimal value as $\mathrm{P}[n]$ [25]. Therefore, we propose to transform $\mathrm{P}[n]$ into its Lagrangian dual domain, which gives the following proposition.

Proposition 1**.**

The dual problem of $\mathrm{P}[n]$ is

[TABLE]

where $\mathbf{q}=[q_{1},\cdots,q_{K}]^{T}\in\mathbb{R}^{K\times 1}$ and

[TABLE]

Moreover, denoting the optimal solution of $\mathbf{q}$ to $\mathrm{D}$ as $\mathbf{q}^{*}$, the optimal $\mathbf{v}^{*}$ of $\mathrm{P}[n]$ is

[TABLE]

Proof.

The Lagrangian of $\mathrm{P}[n]$ is

[TABLE]

With the Lagrangian $\mathcal{L}$, the dual problem of $\mathrm{P}[n]$ is given by [25]

[TABLE]

To compute $\mathop{\mathrm{min}}_{\mathbf{v}}~{}\mathcal{L}\left(\mathbf{v},\mathbf{q}\right)$, we set $\partial\mathcal{L}/\partial~{}\mathrm{conj}(\mathbf{v})=0$, and obtain $\mathbf{v}^{*}=\sum_{k=1}^{K}q_{k}\cdot\mathbf{g}_{k}\mathbf{g}^{H}_{k}\mathbf{v}^{[n]}$, which gives (5). Finally, putting $\mathbf{v}^{*}$ in (5) into (6), the objective function of (6) is

[TABLE]

and the problem (6) is equivalently written as $\mathrm{D}$. ∎

Based on Proposition 1 and by defining

[TABLE]

the function $\Upsilon^{[n]}\left(\mathbf{q}\right)$ in $\mathrm{D}$ can be re-written as $\Upsilon^{[n]}\left(\mathbf{q}\right)=||\mathbf{F}\mathbf{q}||^{2}-\mathbf{d}^{T}\mathbf{q}$, which is quadratic. Therefore, FOMs (e.g., gradient descent, accelerated gradient projection, etc.) can be adopted by computing $\nabla\Upsilon^{[n]}\left(\mathbf{q}\right)$ in each iteration, which require a complexity of $O(KN)$. However, in the following, a stronger property of $\Upsilon^{[n]}(\mathbf{q})$ will be established, thus allowing more efficient updates with per-iteration complexity smaller than $O(KN)$.

Property 2**.**

The function $\Upsilon^{[n]}(\mathbf{q})$ in $\mathrm{D}$ is coordinate-wise $L_{k}$-smooth, where

[TABLE]

Proof.

According to [16, 17], $\Upsilon^{[n]}(\mathbf{q})$ is coordinate-wise $L_{k}$-smooth if and only if

[TABLE]

where $t\in\mathbb{R}$ and $\mathbf{e}_{k}$ represents the unit vector with the $k^{\mathrm{th}}$ element being 1 and others being zero. By computing

[TABLE]

and putting $\nabla_{\mathbf{q}}\Upsilon^{[n]}(\mathbf{q}+t\cdot\mathbf{e}_{k})$ and $\nabla_{\mathbf{q}}\Upsilon^{[n]}(\mathbf{q})$ into (10), the left hand side of (10) becomes $\left[2\mathrm{Re}\left(\mathbf{F}^{H}\mathbf{F}\right)\right]_{k,k}|t|$ and the property is immediately proved. ∎

Based on Property 2 and since the constraint $\mathbf{q}\succeq\mathbf{0}$ is separable (it can be re-written as $q_{k}\geq 0$ for all $k$), problem $\mathrm{D}$ can be solved by the coordinate descent method. In particular, at iteration $m$, we first choose a coordinate, say $l$, and then update

[TABLE]

where the $l^{\mathrm{th}}$ coordinate-wise gradient

[TABLE]

It can be seen from (12) and (13) that the complexity of updating (12) is only $O(K)$.

The next question is how to choose the index $l$. A straightforward idea is to choose $l$ in a cyclic manner, i.e., $l=1,2,\cdots,K,1,2,\cdots,K$. But this may lead to slow convergence of the algorithm [14]. A better way is to choose $l$ as a random integer number ranging from $1$ to $K$ with equal probability, i.e., $\mathrm{Pr}(l=1)=\cdots=\mathrm{Pr}(l=K)=1/K$. With such a random schedule and the coordinate update in (12), the sequence $\{\mathbf{q}^{[0]},\mathbf{q}^{[1]},\cdots\}$ converges to the optimal $\mathbf{q}^{*}$ to $\mathrm{D}$ with a convergence rate of $O(1/m)$ [16, 17]. This is the so-called RCD method.

However, the RCD method can still be improved in the following two ways. First, the update of (12) is sequential, which means that the next iteration must wait before the current iteration is completed. If a block of coordinates is updated in parallel, the running time of RCD can be further reduced [17, 18]. Second, there is a gap between the convergence rate $O(1/m)$ of RCD and the best known convergence rate $O(1/m^{2})$ [26] for solving smooth problems. This indicates that we can consider adding momentums to accelerate the convergence of RCD [19].

Based on the above observations, the following ARCD is adopted for solving $\mathrm{D}$. More specifically, at the $m^{\mathrm{th}}$ iteration, instead of generating a random number $l$, we generate a random set $\mathcal{Y}^{[m]}$ with

[TABLE]

and update all the coordinates in $\mathcal{Y}^{[m]}$ as

[TABLE]

where the momentum $(c^{[m]})^{2}\mathbf{z}^{[m]}$ is added, with $\mathbf{z}^{[m]}$ being the direction at the $m^{\mathrm{th}}$ iteration and $c^{[m]}$ being the step-size to control the importance of $\mathbf{z}^{[m]}$.

It can be seen from (15) that different coordinates in $\mathcal{Y}^{[m]}$ are updated in parallel. On the other hand, by choosing

[TABLE]

with the initial $c^{[0]}=Y/K$ and $\mathbf{z}^{[0]}=\mathbf{0}$, $\mathbf{q}^{[m+1]}$ computed using (15) is guaranteed to converge to the optimal solution to $\mathrm{D}$ with a convergence rate $O(1/m^{2})$ [19, Theorem 3]. Since $\mathrm{D}$ is equivalent to $\mathrm{P}[n]$ and $\mathrm{P}[n]$ represents a surrogate problem for $\mathrm{P}$, problem $\mathrm{P}$ can be solved via MM and ARCD in the Lagrangian dual domain. The entire procedure is summarized in Algorithm 1.

4 Simulation Results

This section presents simulation results to verify the performance of the proposed ARCD. In particular, each random channel is generated according to $\mathcal{CN}(\mathbf{0},\varrho\mathbf{I})$ [24], where the path-loss is $\varrho=-90~{}\mathrm{dB}$. It is assumed that the noise power $\sigma^{2}_{1}=\cdots=\sigma^{2}_{K}=-80~{}\mathrm{dBm}$, which includes thermal noise and receiver noise [27]. Each point in the figures is obtained by averaging over $100$ simulation runs, with independent channels between consecutive runs. All problem instances are solved by Matlab R2015b on a desktop with Intel Core i5-4570 CPU at 3.2 GHz and 8GB RAM. For comparison, we also simulate the MM-IPM111The MM-IPM [23] is implemented using the Matlab software CVX Mosek [25]. [23], the ADMM 222The ADMM is implemented by introducing slack variables $\{\mathbf{w}_{k}=\mathbf{v}\}_{k=1}^{K}$ [12] and adding a consensus penalty $a/2\sum_{k=1}^{K}||\mathbf{v}-\mathbf{w}_{k}+\bm{\eta}_{k}||_{2}^{2}$, with $a=2/\sqrt{N}$ [13] and the dual variables being $\{\bm{\eta}_{k}\}$. [12, 13], and the asymptotic method333The asymptotic solution is implemented by assuming $N\rightarrow+\infty$ [1], and this method is used as initialization for the other simulated methods. [1].

To evaluate the solution quality and the running time of Algorithm 1, we simulate the case of $N=200$ with $K\in\{50,100,200,500\}$. Notice that the case of large $K$ is very important for future crowd sensing applications. For the methods based on MM, the number of MM iterations is $20$ [23, 12, 13]. Furthermore, the ADMM for solving $\mathrm{P}[n]$ stops when the change of objective functions between consecutive iterations is smaller than $\tau=10^{-5}$ [12]. If ADMM fails to reach the above condition within $2000$ iterations, it stops and outputs the result at iteration $2000$ [12, 13]. It can be observed from Fig. 1a that Algorithm 1 significantly outperforms the asymptotic solution [1], and slightly outperforms the ADMM [12, 13]. In fact, Algorithm 1 achieves the same power consumption as the MM-IPM [23]. However, as illustrated in Fig. 1b, Algorithm 1 only requires less than $2$ seconds to finish for all the simulated value of $K$. Compared to MM-IPM and ADMM, Algorithm 1 saves at least $90\%$ (one order of magnitude) of the computation times.

5 Conclusions

This paper studied the massive MIMO multicast beamforming optimization problem. With majorization minimization and strong duality, the primal problem was transformed into a coordinate-wise Lipschitz smooth problem with separable constraints. By further adopting ARCD, lower complexity than that of existing algorithms was achieved. Simulation results showed that the proposed method reduces the execution time by one order of magnitude compared to existing methods while guaranteeing the same performance.