An Efficient Global Algorithm for Single-Group Multicast Beamforming

TL;DR

This paper introduces the first tailored global branch-and-bound algorithm for single-group multicast beamforming, guaranteeing optimal solutions and outperforming existing methods in efficiency and scalability.

Contribution

It presents a novel, guaranteed global optimization algorithm specifically designed for the NP-hard multicast beamforming problem, setting a new benchmark.

Findings

The algorithm guarantees finding the global solution.

It outperforms the Baron solver in computational efficiency.

Performance of existing algorithms declines with larger problem sizes.

Abstract

Consider the single-group multicast beamforming problem, where multiple users receive the same data stream simultaneously from a single transmitter. The problem is NP-hard and all existing algorithms for the problem either find suboptimal approximate or local stationary solutions. In this paper, we propose an efficient branch-and-bound algorithm for the problem that is guaranteed to find its global solution. To the best of our knowledge, our proposed algorithm is the first tailored global algorithm for the single-group multicast beamforming problem. Simulation results show that our proposed algorithm is computationally efficient (albeit its theoretical worst-case iteration complexity is exponential with respect to the number of receivers) and it significantly outperforms a state-of-the-art general-purpose global optimization solver called Baron. Our proposed algorithm provides an important benchmark for performance evaluation of existing algorithms for the same problem. By using it as the benchmark, we show that two state-of-the-art algorithms, semidefinite relaxation algorithm and successive linear approximation algorithm, work well when the problem dimension (i.e., the number of antennas at the transmitter and the number of receivers) is small but their performance deteriorates quickly as the problem dimension increases.

Tables

Table 1. TABLE I: Comparison of our proposed ACR-BB algorithm and Baron

Instance	ACR-BB			Baron
ID	Value	# of Iter.	Time	Value	# of Iter.	Time
01	1.2344	90	0.7	1.2344	6297	327.9
02	1.1412	40	0.2	1.1412	2249	147.4
03	1.3541	68	0.4	1.3541	2917	94.5
04	1.4251	49	0.3	1.4251	2847	261.5
05	1.4136	61	0.4	1.4136	1825	49.8
06	0.8540	40	0.3	0.8540	1303	31.9
07	0.9777	48	0.3	0.9777	457	19.0
08	5.2469	1	0.0	5.2469	4305	58.4
09	11.9555	16	0.1	11.9555	5358	87.1
10	1.3258	62	0.4	1.3258	653	26.5
Average	2.6928	47.5	0.3	2.6928	2821.1	110.4

Table 2. TABLE II: Average and worst-case results of number of iterations and CPU time (in seconds) under different setups

Setting	Average Performance		Worst-Case Performance
$(N,M)$	# of Iter.	Time	# of Iter.	Time
(2,8)	47.9	0.3	103	0.7
(2,16)	72.7	0.5	147	1.0
(2,24)	100.6	0.6	201	1.2
(2,32)	136.0	0.9	225	1.4
(2,40)	151.3	1.0	298	2.0
(2,48)	171.7	1.2	278	1.8
(2,56)	199.9	1.4	343	2.3
(2,64)	215.4	1.5	341	2.4
(4,8)	330.1	2.8	886	8.1
(4,16)	1721.2	16.6	3763	37.8
(4,24)	4492.1	46.4	11627	125.8
(4,32)	8579.6	98.2	17442	216.0
(4,40)	16108.7	184.9	34675	419.9
(6,8)	1728.6	21.5	6388	83.4
(6,16)	13233.1	199.9	36804	593.3
(6,24)	63691.6	1159.7	175041	3679.3
(8,8)	2090.0	30.1	6301	95.1
(8,16)	64168.9	1396.9	235309	5786.3

Table 3. TABLE III: Average objective values and relative gaps obtained by SLA and SDR

Setting	Average Objective Value			Average Relative Gap
$(N,M)$	ACR-BB	SLA	SDR	SLA	SDR
(2,8)	1.629	1.629	1.642	0.02%	0.79%
(2,16)	2.800	2.803	2.827	0.10%	0.96%
(2,24)	3.389	3.398	3.463	0.26%	2.20%
(2,32)	4.457	4.489	4.642	0.73%	4.16%
(2,40)	5.720	5.731	5.871	0.20%	2.65%
(2,48)	5.679	5.739	5.892	1.06%	3.77%
(2,56)	5.461	5.491	5.758	0.55%	5.44%
(2,64)	5.914	5.967	6.291	0.89%	6.37%
(4,8)	0.514	0.525	0.577	2.01%	12.17%
(4,16)	0.837	0.902	1.121	7.83%	33.97%
(4,24)	1.132	1.256	1.710	10.95%	51.03%
(4,32)	1.328	1.587	2.045	19.48%	54.04%
(4,40)	1.525	1.770	2.587	16.08%	69.64%
(6,8)	0.334	0.341	0.393	2.30%	17.75%
(6,16)	0.531	0.578	0.799	8.90%	50.39%
(6,24)	0.667	0.779	1.113	16.77%	66.99%
(8,8)	0.246	0.253	0.278	2.95%	12.93%
(8,16)	0.376	0.420	0.618	11.68%	64.11%

Table 4. TABLE IV: Worst-case relative gaps obtained by SLA and SDR.

$(N,M)$	SLA	SDR
(2,8)	1%	6%
(2,16)	6%	10%
(2,24)	4%	18%
(2,32)	13%	30%
(2,40)	8%	22%
(2,48)	10%	28%
(2,56)	12%	25%
(2,64)	14%	25%
(4,8)	34%	42%
(4,16)	39%	82%
(4,24)	65%	100%
(4,32)	70%	103%
(4,40)	74%	104%
(6,8)	22%	92%
(6,16)	62%	126%
(6,24)	73%	128%
(8,8)	23%	46%
(8,16)	51%	124%

Table 5. TABLE V: Probability of obtaining the global solution by SLA and SDR

$(N,M)$	SLA	SDR
(2,8)	98%	74%
(2,16)	96%	68%
(2,24)	92%	48%
(2,32)	88%	28%
(2,40)	92%	30%
(2,48)	82%	22%
(2,56)	86%	24%
(2,64)	88%	14%
(4,8)	80%	38%
(4,16)	50%	6%
(4,24)	42%	0%
(4,32)	28%	0%
(4,40)	14%	0%
(6,8)	74%	40%
(6,16)	44%	0%
(6,24)	18%	0%
(8,8)	72%	28%
(8,16)	42%	2%

Equations

$$y_k(t)=s(t){\mathbf{h}}_k^H\mathbf{w}+n_k(t),$$

$${\text{SNR}}_k=\frac{|{\mathbf{h}}_k^H\mathbf{w}{|}^2}{{\sigma }_k^2},k=1,2,\dots ,M.$$

$$\begin{array}[]{cl}\displaystyle\min_{\mathbf{w}}&\|\mathbf{w}\|^{2}\\[10.0pt] \mbox{s.t.}&\displaystyle\frac{|\mathbf{h}_{k}^{H}\mathbf{w}|^{2}}{\sigma_{k}^{2}}\geq\gamma_{k},~{}k=1,2,...,M,\end{array}$$

$$\frac{|{\mathbf{h}}_k^H\mathbf{w}{|}^2}{{\sigma }_k^2}\ge {\gamma }_k⟺|{\tilde{\mathbf{h}}}_k^H\mathbf{w}{|}^2\ge 1.$$

$$\begin{array}[]{cl}\displaystyle\min_{\mathbf{w}}&\|\mathbf{w}\|^{2}\\[10.0pt] \mbox{s.t.}&{\displaystyle|\mathbf{h}_{k}^{H}\mathbf{w}|^{2}\geq 1},~{}k=1,2,\ldots,M.\end{array}$$

$$\begin{array}[]{cl}\displaystyle\max_{\mathbf{w}}&\displaystyle\min_{k=1,2,\ldots,M}\left\{|\mathbf{h}_{k}^{H}\mathbf{w}|^{2}\right\}\\[10.0pt] \mbox{s.t.}&\displaystyle\|\mathbf{w}\|^{2}\leq 1.\end{array}$$

$$\text{Trace}({\mathbf{H}}_k\mathbf{W})\ge 1,\mathbf{W}⪰0,\text{and}\text{Rank}(\mathbf{W})=1,$$

$$\begin{array}[]{cl}\displaystyle\min_{\mathbf{W}}&\textrm{Trace}(\mathbf{W})\\[5.0pt] \mbox{s.t.}&\displaystyle\textrm{Trace}(\mathbf{H}_{k}\mathbf{W})\geq 1,~{}k=1,\ldots,M,\\[5.0pt] &\mathbf{W}\succeq\mathbf{0}.\end{array}$$

$${\mathbf{v}}_k:={\left[\mathrm{Re}({\mathbf{h}}_k^H\mathbf{w}),\mathrm{Im}({\mathbf{h}}_k^H\mathbf{w})\right]}^T,k=1,2,\dots ,M.$$

$$\Vert {\mathbf{v}}_k^n{\Vert }^2+2({\mathbf{v}}_k^n{)}^T({\mathbf{v}}_k-{\mathbf{v}}_k^n)\ge 1,$$

$$\begin{array}[]{cl}\displaystyle\min_{\mathbf{w},\,\mathbf{v}}&\|\mathbf{w}\|^{2}\\[5.0pt] \mbox{s.t.}&\displaystyle\|\mathbf{v}_{k}^{n}\|^{2}+2(\mathbf{v}_{k}^{n})^{T}(\mathbf{v}_{k}-\mathbf{v}_{k}^{n})\geq 1,~{}k=1,...,M,\\[5.0pt] &\displaystyle\mathbf{v}_{k}=\left[\mathcal{R}(\mathbf{h}_{k}^{H}\mathbf{w}),\mathcal{I}(\mathbf{h}_{k}^{H}\mathbf{w})\right]^{T},~{}k=1,...,M,\end{array}$$

$$\Vert {\mathbf{v}}_k{\Vert }^2\ge \Vert {\mathbf{v}}_k^n{\Vert }^2+2({\mathbf{v}}_k^n{)}^T({\mathbf{v}}_k-{\mathbf{v}}_k^n),\forall {\mathbf{v}}_k.$$

$$\mathbf{c}={\left[c_1,c_2,\dots ,c_{M-1}\right]}^T$$

$$\begin{array}[]{cl}\displaystyle\min_{\mathbf{w},\,\mathbf{c}}{}&\|\mathbf{w}\|^{2}\\[5.0pt] \mbox{s.t.}&c_{k}=\mathbf{h}_{k}^{H}\mathbf{w},~{}k=1,...,M-1,\\[5.0pt] &|c_{k}|\geq 1,~{}k=1,...,M-1,\\[5.0pt] &\mathbf{h}_{M}^{H}\mathbf{w}\geq 1.\end{array}$$

$$\{c_k\in \mathbb{C}||c_k|\ge 1\}.$$

$$\mathcal{D}_{[l_{k},u_{k}]}=\left\{(x_{k},y_{k})\,\left|\,\begin{array}[]{@{}lll}c_{k}=x_{k}+y_{k}\textbf{i},~{}|c_{k}|\geq 1\\ \arg{(c_{k})}\in[l_{k},u_{k}]\end{array}\right.\right\}$$

$${\mathcal{D}}_{[0,\pi /2]}=\left\{(x,y)|x\ge 0,y\ge 0,\sqrt{x^2+y^2}\ge 1\right\}$$

$$\text{Conv}({\mathcal{D}}_{[0,\pi /2]})=\{(x,y)|x\ge 0,y\ge 0,x+y\ge 1\}$$

$$\textrm{Conv}(\mathcal{D}_{[l_{k},u_{k}]})=\left\{(x,y)\,\left|\,\begin{array}[]{@{}lll}\sin(l_{k})x-\cos(l_{k})y\leq 0\\ \sin(u_{k})x-\cos(u_{k})y\geq 0\\ a_{k}x+b_{k}y\geq a_{k}^{2}+b_{k}^{2}\end{array}\right.\right\}.$$

$$a_k=\frac{\cos (l_k)+\cos (u_k)}{2}\text{and}{b}_k=\frac{\sin (l_k)+\sin (u_k)}{2}.$$

$${\mathcal{F}}_{[l,u]}:=\left\{x+y\text{i}|(x,y)\in \text{Conv}({\mathcal{D}}_{[l,u]})\right\}$$

$$\left\{x+y\text{i}|(x,y)\in {\mathcal{D}}_{[l,u]}\right\}.$$

$$\begin{array}[]{cl}\displaystyle\min_{\mathbf{w},\,\mathbf{c}}&\|\mathbf{w}\|^{2}\\[5.0pt] \mbox{s.t.}&c_{k}=\mathbf{h}_{k}^{H}\mathbf{w},~{}k=1,\dots,M-1,\\[5.0pt] &c_{k}\in\mathcal{F}_{[l_{k},u_{k}]},~{}k=1,\ldots,M-1,\\[5.0pt] &\mathbf{h}_{M}^{H}\mathbf{w}\geq 1.\end{array}$$

$$\{x+y\text{i}|(x,y)\in {\mathcal{D}}_{[l_k,u_k]}\}\supseteq {\mathcal{F}}_{[l_k,u_k]}$$

$$\underset{c_k\in {\mathcal{F}}_{[l_k,u_k]}}{\min }|c_k|=\cos \left(\frac{u_k-l_k}{2}\right).$$

$$\left(\frac{\cos (l_k)+\cos (u_k)}{2},\frac{\sin (l_k)+\sin (u_k)}{2}\right)$$

$$\begin{array}[]{cl}\displaystyle\min_{\mathbf{w}}&\|\mathbf{w}\|^{2}\\[5.0pt] \mbox{s.t.}&\mathbf{h}_{M}^{H}\mathbf{w}\geq 1.\end{array}$$

$${\mathbf{w}}^0=\frac{{\mathbf{h}}_M}{\Vert {\mathbf{h}}_M{\Vert }^2},L^0=\frac{1}{\Vert {\mathbf{h}}_M{\Vert }^2},$$

$$c_k^0=\frac{{\mathbf{h}}_k^H{\mathbf{h}}_M}{\Vert {\mathbf{h}}_M{\Vert }^2},k=1,2,\dots ,M-1.$$

$$(U^t-L^t)/L^t\le ϵ,$$

Full text

I Introduction

Physical layer multicasting via transmit beamforming has been recognized as a powerful technique for efficient audio and video streaming in multi-user multi-antenna wireless communication networks. For instance, multicast beamforming is a part of the Evolved Multimedia Broadcast Multicast Service (eMBMS) in the Long-Term Evolution (LTE) standard, commercially known as LTE Broadcast [1]. Multicast beamforming exploits channel state information at the transmitter and utilizes multiple transmit antennas to broadcast common information to a preselected group of users.

One scenario of particular interest in this paper is single-group multicast beamforming, where all users receive the same data stream from the transmitter and the data rate is determined by the minimum received signal-to-noise-ratio (SNR). The earliest mathematical formulation of the single-group multicast beamforming problem is to maximize the average SNR subject to the total transmission power constraint [2]. This formulation is simple, and can be solved efficiently. However, the solution that maximizes the average SNR does not consider the SNR of each individual user, so that the minimum received SNR may be significantly lower than the average SNR. Hence, the solution obtained by maximizing the average SNR does not always achieve a satisfactory common data rate. To overcome this drawback, reference [3] proposed two new problem formulations, the quality of service (QoS) constrained problem formulation and the max-min fairness problem formulation, where the former one minimizes the total transmission power subject to SNR constraints of all receivers and the latter one maximizes the minimum SNR among all users subject to the total transmission power constraint. These two new formulations can guarantee QoS of each user and are shown to be equivalent from an optimization point of view [3]. Unfortunately, the two problems are NP-hard in general [4] and thus it is impossible to solve them to global optimality in polynomial time (unless P=NP) [5]. This is sharply different from the downlink unicast beamforming problem [6, 7, 8] and the uplink/downlink coordinated beamforming problem [9, 10, 11], which can be equivalently reformulated as a convex problem and can be efficiently solved to global optimality such as by using the uplink-downlink duality theory.

Various algorithms have been proposed to solve the single-group multicast beamforming problem; see [4, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27] and references therein. Based on our knowledge, these algorithms are either convex relaxation based algorithms or convex approximation based algorithms. Due to the NP-hardness of the problem, none of them are guaranteed to find the global solution of the problem (except for very special problem instances [18, 28]). Moreover, since there might exist large gaps between the original problem and its convex relaxations/approximations, the quality of the returned solutions by the aforementioned algorithms based on these convex relaxations/approximations might be poor.

The goal of this paper is to develop a computationally efficient global algorithm for solving the single-group multicast beamforming problem that is guaranteed to find the global solution of the problem.

I-A Related Works

One of the state-of-the-art algorithms for the single-group multicast beamforming problem is the semidefinite relaxation (SDR) algorithm [4]. The main observation behind the SDR algorithm is that the beamforming problem can be equivalently formulated as a rank-one constrained semidefinite program (SDP). The SDR algorithm drops the rank-one constraint, solves the SDR, and then applies a Gaussion randomization strategy to generate a rank-one approximate solution of the original problem based on the obtained solution of the SDR. The SDR algorithm is capable of finding high quality approximate solutions when the number of antennas at the transmitter and the number of users are small. However, the performance of the SDR algorithm deteriorates quickly as the number of users increases. As justified in [29], the provable (worst-case) approximation accuracy of the solution returned by the SDR algorithm degrades linearly with the number of users.

One of the best algorithms for the problem of interest is the successive linear approximation (SLA) algorithm [12]. At each iteration, the SLA algorithm approximates the original nonconvex problem by using its first-order Talor series expansion at the current iterate and then solves the resulting convex quadratic program to generate the next iterate. It has been shown in [12] that the sequence of points generated by the SLA algorithm converges to a Karush-Kuhn-Tucker (KKT) point of the original beamforming problem. Moreover, numerical simulation results in [12] have shown that the SLA algorithm performs better than the SDR algorithm (in the sense that the SLA algorithm finds a better solution with a lower objective value than the SDR algorithm).

Recently, an alternating maximization (AM) algorithm for the single-group multicast beamforming problem was proposed in [13]. Recall that the beamforming problem is equivalent to an SDP with a rank-one constraint [4]. The main contribution of [13] is that it reformulated the problem as a nonconvex problem without the rank-one constraint, which is naturally amenable to AM. Numerical simulation results in [13] show that the AM algorithm performs better than the SDR algorithm. However, further simulation results in [14] show that the AM algorithm performs worse and has higher complexity than the SLA algorithm.

We also mention some interesting related works, which extend the previously mentioned works [4, 12, 13] on the single-group multicast beamforming problem from different aspects. One extension is to develop more efficient algorithms for solving the problem when computational efficiency is a big issue, i.e., when the channel changes fast over time. Along this direction, low-complexity algorithms [15, 16, 17, 18, 19] as well as adaptive (online) algorithms (which learn the channel correlation matrices) [20, 21, 22, 23] have been proposed. These algorithms generally are very fast, but their performance is worse than the ones of the SDR algorithm [4] and the SLA algorithm [12]. Recently, a new adaptive multiplicative update (MU) algorithm and a hybrid MU-SLA algorithm based on it have been proposed in [14]. The MU-SLA algorithm runs the SLA algorithm for only one iteration with the returned point by the MU algorithm as the initial point.

Another extension is to develop new physical-layer transmit strategies instead of using the beamforming transmit strategy. Along this direction, a beamformed Alamouti scheme has been independently proposed in [24, 25, 26], which can be seen as a rank-two generalization of the previous (rank-one) SDR beamforming framework. It has been shown in [27] that the worst-case approximation accuracy of the beamformed Alamouti scheme degrades only at a rate of the square root of the number of the users. This improves over the beamforming strategy, where the approximation accuracy degrades at a rate of the number of users [29].

I-B Our Contributions

This is the first paper that proposes a tailored efficient global algorithm for solving the single-group multicast beamforming problem, which is in sharp contrast to all existing works that focus on the design of approximation algorithms or local optimization algorithms. Our proposed algorithm is based on the branch-and-bound strategy combined with a new argument cut technique. The argument cuts (in Definition 1) are used to design effective convex relaxations of nonconvex constraints in the single-group multicast beamforming problem and therefore play an important role in our proposed branch-and-bound algorithm for solving the problem.

Since the single-group multicast beamforming problem is NP-hard, there does not exist a polynomial time algorithm which can solve it to global optimality (unless P=NP) [5]. Therefore, our proposed algorithm has an exponential worst-case iteration complexity with respect to the number of users (see Theorem 1). However, our simulation results show that our proposed algorithm is highly efficient and it significantly outperforms the state-of-the-art general-purpose global optimization solver Baron [30, 31]. The high efficiency of our proposed algorithm is mainly due to the new argument cuts. More specifically, our proposed algorithm can globally solve small-scale problem instances with $2$ antennas at the transmitter and $8$ users within $0.3$ seconds on average while Baron needs $110.4$ seconds; our proposed algorithm can solve median-scale problem instances with $4$ antennas at the transmitter and $8$ users within $2.8$ seconds on average while Baron fails to solve the same problem instances within $10$ minutes.

Even though in some scenarios (i.e., both the number of antennas at the transmitter and the number of users are large) it takes our proposed algorithm a (relatively) large number of iterations (and thus long time) to find the global solution, our proposed algorithm is still attractive from the following two aspects. First, when the channel changes slowly with time and computational efficiency is not an issue, our proposed algorithm is capable of finding the global solution and therefore provides (potentially much) better performance compared to the existing algorithms. Second, when the channel changes fast with time and computational efficiency becomes an issue in this case, our proposed algorithm can serve as a benchmark to evaluate the performance of the existing heuristic/local optimization algorithms. Our simulation results show that the relative gaps between the objective values at the solutions returned by the SDR algorithm and the SLA algorithm and the optimal objective value (achieved by our proposed algorithm) exceed 50% for some problem instances with a (relatively) large number of antennas at the transmitter and a (relatively) large number of receivers.

I-C Organization and Notations

The organization of this paper is as follows. In Section II, we introduce the QoS constrained formulation of the single-group multicast beamforming problem, and review two state-of-the-art algorithms for solving the problem. In Section III, we propose a global branch-and-bound algorithm for solving the problem111The MATLAB codes of the proposed algorithm are available at [https://www.dropbox.com/s/safrgm97emdgyl9/ARC-BB.rar?dl=0].. Simulation results are presented in Section IV to illustrate the efficiency of our proposed algorithm and the paper is concluded in Section V.

We adopt the following notations in this paper. We use lowercase boldface and uppercase boldface letters to denote (column) vectors and matrices, respectively. We use $\mathbb{C}$ ($\mathbb{R}$) to denote the complex (real) domain, $\mathbb{C}^{n}$ ($\mathbb{R}^{n}$) to denote the set of the $n$-dimensional complex (real) column vectors, and $\mathbb{C}^{m\times n}$ to denote the set of $m\times n$ complex matrices. For a given complex number $c,$ the notations $\mathrm{Re}(c),~{}\mathrm{Im}(c),~{}{|c|},$ and $\arg{(c_{i})}$ stand for its real part, its imaginary part, its modulus, and its argument, respectively. For a given (complex) vector $\mathbf{x}$, $\|\mathbf{x}\|$ denotes its Euclidean norm, $\mathbf{x}^{H}$ denotes its Hermitian, and $\mathbf{x}^{T}$ denotes its transpose. The same notations $\mathbf{A}^{H}$ and $\mathbf{A}^{T}$ also apply to the matrix $\mathbf{A}.$ For a given complex Hermitian matrix $\mathbf{A}$, $\mathbf{A}\succeq\mathbf{0}$ means $\mathbf{A}$ is positive semidefinite, $\textrm{Trace}(\mathbf{A})$ denotes the trace of the matrix $\mathbf{A}$, and $\textrm{Rank}(\mathbf{A})$ denotes the rank of the matrix $\mathbf{A}$. For two given Hermitian matrices $\mathbf{A}$ and $\mathbf{B}$, $\mathbf{A}\succeq\mathbf{B}$ means $\mathbf{A}-\mathbf{B}\succeq\mathbf{0}$. Finally, we use $\mathbf{i}$ to denote the imaginary unit which satisfies the equation $\mathbf{i}^{2}=-1$ and use $\mathbf{I}_{N}$ to denote the $N\times N$ identity matrix.

II Problem Formulation and Review

In this section, we first introduce the QoS constrained formulation of the single-group multicast beamforming problem and then review two state-of-the-art algorithms for solving the problem.

II-A Problem Formulation

Consider the single-group multicast beamforming problem for the multi-user multi-input single-output (MISO) downlink channel, where the base station (transmitter) is equipped with $N$ antennas, and broadcasts a common data stream to $M$ users with a single antenna. Let $\mathbf{h}_{k}\in\mathbb{C}^{N}$ denote the channel vector between the base station and the $k$-th receiver and let $\mathbf{w}\in\mathbb{C}^{N}$ denote the beamforming vector used by the base station. Assume that $s(t)$ is the broadcasting data stream. Then the transmitted signal by the base station is given by $s(t)\mathbf{w}$ and the received signal at the $k$-th receiver is given by

[TABLE]

where $n_{k}(t)$ is the additive white Gaussian noise (AWGN) with variance $\sigma_{k}^{2}$. Assume that the transmitted signal $s(t)$ has a unit power. Then, the SNR of the $k$-th user can be written as

[TABLE]

This paper is interested in minimizing the total transmission power at the base station while satisfying the SNR constraints of all users. Mathematically, the problem can be formulated as follows:

[TABLE]

where $\gamma_{k}$ is the desired transmission SNR target of user $k.$ Let $\tilde{\mathbf{h}}_{k}=\mathbf{h}_{k}/\sqrt{\gamma_{k}\sigma_{k}^{2}}$, then

[TABLE]

For ease of notation, we drop the ~ and study the following single-group multicast beamforming problem in this paper:

[TABLE]

A closely related problem is to maximize the minimum SNR among all users subject to the total transmission power constraint:

[TABLE]

It has been shown in [3] that problems (P) and (2) are equivalent from an optimization perspective of view, i.e., the global solution of problem (P) can be obtained by appropriately scaling the global solution of problem (2), and vice versa. Therefore, we focus on problem (P) in this paper.

II-B Review of Two State-of-the-Art Algorithms

In this subsection, we briefly review two state-of-the-art algorithms for solving problem (P). We shall compare our proposed algorithm with these two algorithms later in Section IV.

One of the state-of-the-art algorithm for solving problem (P) is the SDR algorithm [4]. The SDR algorithm is based on the SDR technique, which has been widely used to solve optimization problems arising from signal processing and wireless communications (see [28] and references therein). The main observation in the SDR algorithm is that the constraint $|\mathbf{h}_{k}^{H}\mathbf{w}|^{2}\geq 1$ can be equivalently rewritten as

[TABLE]

where $\mathbf{H}_{k}=\mathbf{h}_{k}\mathbf{h}_{k}^{H}\in\mathbb{C}^{N\times N}$ and $\mathbf{W}=\mathbf{w}\mathbf{w}^{H}\in\mathbb{C}^{N\times N}.$ By dropping the rank-one constraint, problem (P) is relaxed to

[TABLE]

The above SDR problem (3) can be solved in polynomial time by using the interior-point algorithm [32]. If the optimal solution $\mathbf{W}^{*}\succeq\mathbf{0}$ of the SDR problem (3) is of rank one, i.e., $\mathbf{W}^{*}$ admits the decomposition $\mathbf{W}^{*}={\mathbf{w}^{*}}(\mathbf{w}^{*})^{H},$ then ${\mathbf{w}^{*}}$ is a global solution of problem (P). However, the solution $\mathbf{W}^{*}$ of the SDR problem (3) is not always of rank one and thus the global solution of problem (P) might not be obtained. In this case, the SDR algorithm employs the Gaussian randomization techniques to randomly generate approximate solutions based on $\mathbf{W}^{*},$ then scales the approximate solutions to satisfy all SNR constraints, and finally picks the one that has the smallest norm as the final solution. The SDR algorithm can find high quality approximate solutions when the problem dimension, especially the number of users, is small. The worst-case approximation accuracy of the SDR algorithm was shown in [4, 29].

Another state-of-the-art algorithm for solving problem (P) is the SLA algorithm [12]. The basic idea of the SLA algorithm is to approximate the nonconvex constraint $|\mathbf{h}_{k}^{H}\mathbf{w}|^{2}\geq 1$ by a linear constraint. Specifically, the SLA algorithm introduces the auxiliary variables

[TABLE]

Given the current point $\left\{\mathbf{v}_{k}^{n}\in\mathbb{R}^{2}\right\}$ (at the $n$-th iteration), the SLA algorithm first approximates the nonconvex constraint $|\mathbf{h}_{k}^{H}\mathbf{w}|^{2}\geq 1$ by the linear constraint

[TABLE]

and then solves the following approximation problem to obtain the next iterate $\left\{\mathbf{v}_{k}^{n+1}\right\}:$

[TABLE]

where $\mathbf{v}$ is a collection of $\left\{\mathbf{v}_{k}\right\}_{k=1}^{M}.$ The convergence of the SLA algorithm to a KKT point has been established in [12]. It is worthwhile remarking that the subproblem (4) in the SLA algorithm is a linearly constrained convex quadratic program and thus can be solved efficiently to global optimality [32]. Moreover, for any given $\mathbf{v}_{k}^{n},$ there holds

[TABLE]

Therefore, the feasible region of subproblem (4) is a subset of that of the original problem (P) and the SLA algorithm is a convex (inner) approximation algorithm. This differs from the SDR algorithm which is a convex relaxation algorithm. One potential drawback of the SLA algorithm is that its performance depends on the choice of the initial point $\left\{\mathbf{v}_{k}^{0}\right\}.$ To overcome this, [12] proposed to randomly generate many points, scale them, and pick the best one as the initial point.

III Proposed Global Branch-and-Bound Algorithm

In this section, we propose a global optimization algorithm for solving problem (P). Our proposed algorithm is based on the branch-and-bound scheme, which is a general framework for designing global optimization algorithms [33].

A typical branch-and-bound algorithm (for the minimization problem) is generally based on an enumeration procedure, which partitions the feasible region to smaller subregions and constructs sub-problems over the partitioned subregions recursively. In the enumeration procedure, a lower bound for each subproblem is estimated by solving a relaxation problem. Meanwhile, an upper bound is obtained from the best known feasible solution generated by the enumeration procedure or by some other local optimization/heuristic algorithms. A subproblem with a lower bound being larger than the obtained upper bound is called as an inactive subproblem, which does not contain the global solution of the original problem in its feasible region, and thus will not be further enumerated. The procedure terminates until all active subproblems have been enumerated, and then an optimal solution within a given relative error tolerance can be obtained.

The efficiency of a branch-and-bound algorithm considerably relies on the quality of the lower bound as well as the upper bound. The quality of the lower bound depends on the tightness of the convex relaxation and the one of the upper bound depends on the local optimization or heuristic algorithms which are employed to generate the feasible solutions (to the original problem). With better lower and upper bounds, more inactive subproblems can be detected and more unnecessary enumerations can be avoided.

In the remaining part of this section, we first propose a convex quadratic programming relaxation in Section III-A, which provides efficient lower bounds in our branch-and-bound algorithm. Then, we present our proposed branch-and-bound algorithm for solving problem (P) in Section III-B. Finally, we show that our proposed branch-and-bound algorithm indeed can find the global solution of problem (P) within any given positive relative error tolerance and analyze its worst-case iteration complexity in Section III-C.

III-A New Argument Cut based Relaxation

As is well known, the SDR (3) is the tightest convex relaxation of problem (P) and its optimal value provides high quality lower bounds on that of problem (P). However, the decision variable in the SDR is lifted to an $N\times N$ matrix, whose dimension is much larger than the dimension $N$ of the original variable $\mathbf{w}$. In comparison, the subproblem (4) is a convex linearly constrained quadratic program with $N$ variables, which can be solved much more efficiently than the SDR (3). However, the subproblem (4) is a convex inner approximation of problem (P) but not a convex relaxation, which makes its objective value not suitable for serving as a lower bound. In the next, we propose a new convex relaxation, which achieves much higher computational efficiency than the SDR and provides valid lower bounds with satisfactory tightness.

Without loss of generality, we first change the constraint $|\mathbf{h}_{M}^{H}\mathbf{w}|^{2}\geq 1$ (which is equivalent to $|\mathbf{h}_{M}^{H}\mathbf{w}|\geq 1$) to $\mathbf{h}_{M}^{H}\mathbf{w}\geq 1$ in problem (P). This is because that, for any $\mathbf{w}$ satisfying $|\mathbf{h}_{M}^{H}\mathbf{w}|\geq 1,$ we can always find an appropriate $\theta\in\mathbb{R}$ such that $\exp({\textbf{i}\theta})\mathbf{h}_{M}^{H}\mathbf{w}\geq 1.$ Second, we introduce a new variable

[TABLE]

with $c_{k}=\mathbf{h}_{k}^{H}\mathbf{w}$, and represent the constraint $|\mathbf{h}_{k}^{H}\mathbf{w}|^{2}\geq 1$ by $|c_{k}|^{2}\geq 1$, which is equivalent to $|c_{k}|\geq 1$. Then, problem (P) is transformed to the following problem (P’):

[TABLE]

In problem (P’), $|c_{k}|\geq 1~{}(k=1,...,M-1)$ are the only nonconvex constraints. To develop the branch-and-bound algorithm, we need to relax these nonconvex constraints to convex ones. It is worthwhile remarking here that we treat the last SNR constraint different from the others in problem (P’). The purpose of doing this is to reduce the number of nonconvex constraints in it and this will make the corresponding relaxation problem more tight and thus improve the efficiency of our proposed branch-and-bound algorithm.

Let us consider the set

[TABLE]

It is simple to see that the convex envelope of set (5) is the whole complex set $\mathbb{C}$. Obviously, it is loose, if we directly relax set (5) to its convex envelope, i.e., relax the constraint $|c_{k}|\geq 1$ to $c_{k}\in\mathbb{C}$.

Next, we develop a tighter convex relaxation for (5). To do so, we introduce $x_{k}=\mathrm{Re}(c_{k})$ and $y_{k}=\mathrm{Im}(c_{k})$, and assume that the argument of $c_{k}$ satisfies $\arg{(c_{k})}\in[l_{k},u_{k}]$. Let

[TABLE]

and let $\textrm{Conv}(\mathcal{D}_{[l_{k},u_{k}]})$ be the convex envelope of the set $\mathcal{D}_{[l_{k},u_{k}]}.$

Let us first give an illustration on how $\mathcal{D}_{[l_{k},u_{k}]}$ and $\textrm{Conv}(\mathcal{D}_{[l_{k},u_{k}]})$ look like by using Fig. 1 where $[l_{k},u_{k}]=[0,\pi/2].$ In Fig. 1,

[TABLE]

is the light blue region outside the arc AB, which is obviously a nonconvex set; and

[TABLE]

is the light blue region, which is convex and is determined by three lines $\overline{\text{OA}}$, $\overline{\text{OB}}$, and $\overline{\text{AB}}$.

The following proposition characterizes $\textrm{Conv}(\mathcal{D}_{[l_{k},u_{k}]})$ in the general case if $u_{k}-l_{k}\leq\pi.$

Proposition 1**.**

Suppose $l_{k}$ and $u_{k}$ in (6) satisfying $u_{k}-l_{k}\leq\pi$. Then

[TABLE]

*where *

[TABLE]

Proof.

It is simple to see that the two extreme points of the set $\textrm{Conv}(\mathcal{D}_{[l_{k},u_{k}]})$ are $(\cos(l_{k}),\sin(l_{k}))$ and $(\cos(u_{k}),\sin(u_{k}))$. With a little abuse of notation, we denote these two points by B and A (which reduce to the points B and A in Fig. 1 where $[l_{k},u_{k}]=[0,\pi/2]$). To prove the proposition, it is sufficient to determine the three lines $\overline{\text{OA}}$, $\overline{\text{OB}}$, and $\overline{\text{AB}}.$ The equations of these three lines are given by $\sin(l_{k})x-\cos(l_{k})y=0,$ $\sin(u_{k})x-\cos(u_{k})y=0,$ and $a_{k}x+b_{k}y=a_{k}^{2}+b^{2}_{k},$ respectively, where $a_{k}$ and $b_{k}$ are given in (8). The proof is completed. ∎

The linear inequalities in (7) are introduced due to the argument constraints. Hence, we name these linear inequalities the argument cuts in this paper.

Definition 1** (Argument Cuts).**

The linear inequalities in $\textrm{Conv}(\mathcal{D}_{[l_{k},u_{k}]})$ are called the argument cuts.

To the best of our knowledge, the argument cuts have not been used in the literatures. With the help of the argument cuts, we are able to develop efficient convex relaxations for problem (P’). More specifically, it follows from Proposition 1 that

[TABLE]

is the convex envelope of the nonconvex set

[TABLE]

Assume $\arg(c_{k})\in[l_{k},u_{k}]$ for all $k=1,2,\ldots,M-1$ in problem (P’). Then we obtain the following convex argument cut based relaxation (ACR) of problem (P’):

[TABLE]

We consider the following two cases: whether or not the width of the interval $[l_{k},u_{k}],$ i.e., $d_{k}:=u_{k}-l_{k},$ is greater than $\pi.$

Case I: $d_{k}>\pi.$ In this case, $c_{k}\in\mathcal{F}_{[l_{k},u_{k}]}$ is equivalent to $c_{k}\in\mathbb{C}$ and thus the corresponding argument cuts do not take effect. Therefore, we can drop the constraint $c_{k}\in\mathcal{F}_{[l_{k},u_{k}]}$ in problem (9).

• -

Case II: $d_{k}\leq\pi.$ It follows from Proposition 1 that $c_{k}\in\mathcal{F}_{[l_{k},u_{k}]}$ can be represented by at most222Notice that when $d_{k}=\pi$, the first two linear constraints in (7) are equivalent to each other, and the third linear constraint becomes redundant. three (real) linear inequality constraints. It is simple to see that the argument cuts are effective in this case.

To sum up, we know that ACR problem (9) is a linearly constrained convex quadratic program with $2N+2M-2$ (real) variables, at most $2M-1$ (real) linear equality constraints333Recall that the constraint $\mathbf{h}_{M}^{H}\mathbf{w}\geq 1$ corresponds a (real) linear inequality constraint and a (real) linear equality constraint., and at most $3M-2$ (real) linear inequality constraints. Therefore, ACR problem (9) can be solved efficiently and globally by using the interior-point algorithm within ${\cal{O}}\left(N^{3}M^{3.5}\right)$ operations [32, Page 423].

We conclude this subsection with a tightness measure of the ACR. Let us focus on the case $d_{k}\leq\pi.$ Note that the inclusion

[TABLE]

generally is not true, which implies that there exist some point $c_{k}\in\mathcal{F}_{[l_{k},u_{k}]}$ such that $|c_{k}|<1.$ Therefore, the smallest norm of the points in $\mathcal{F}_{[l_{k},u_{k}]}$ can be used to measure the tightness of the ACR. The following theorem shows that the smallest norm of the points in $\mathcal{F}_{[l_{k},u_{k}]}$ can be computed in a closed form.

Proposition 2**.**

Given any interval $[l_{k},u_{k}]$ with $u_{k}-l_{k}\leq\pi$. We have

[TABLE]

Proof.

The point that has the smallest norm in a convex set is the projection of the origin onto the corresponding set. We first consider the special case where $[l_{k},u_{k}]=[0,\pi/2],$ as shown in Fig. 1. In this case, the projection of the origin $O$ onto the set $\mathcal{F}_{[0,\pi/2]}$ is the point $C,$ which satisfies $\overline{\text{OC}}\bot\overline{\text{AB}}.$ Notice that the coordinates of $A$ and $B$ in Fig. 1 are $(0,1)$ and (1,0), respectively. It is simple to see that the coordinate of the point $C$ is $(0.5,0.5).$ In the general case of $\mathcal{F}_{[l_{k},u_{k}]},$ the coordinate of the point $c_{k}$ is

[TABLE]

and its norm is $\cos\left(\frac{u_{k}-l_{k}}{2}\right).$ The proof is completed. ∎

We can see from Proposition 2 that: the smaller the width of the interval, the tighter the ACR; as the width of the interval goes to zero, the set $\mathcal{F}_{[l_{k},u_{k}]}$ becomes $\{x+y\textbf{i}\,|\,(x,y)\in\mathcal{D}_{[l_{k},u_{k}]}\}$ and the ACR becomes tight. Hence, an effective approach to tightening the ACR is to reduce the width of the corresponding interval.

III-B Proposed Branch-and-Bound Algorithm

In this subsection, we propose a branch-and-bound algorithm for globally solving problem (P’). The basic idea of the proposed algorithm is to relax the original problem (with appropriate argument constraints) to ACR (9) and gradually tighten the relaxation by reducing the width of the associated intervals.

For ease of presentation, we introduce the following notations. Let $\mathcal{A}=\prod_{k=1}^{M-1}[l_{k},u_{k}]$ and let ACR($\mathcal{A}$) denote the ACR problem defined over the set $\mathcal{A};$ let $\mathcal{P}$ denote the constructed problem list and let $\{\mathcal{A},\mathbf{c},L\}$ denote a problem instance from the list $\mathcal{P},$ where $L$ is the optimal value of ACR($\mathcal{A}$) and $\mathbf{c}$ is its optimal solution; let superscript $t$ denote the iteration number; let ${U^{t}}$ denote the upper bound at the $t$-th iteration; and let $\mathbf{w}^{*}$ denote the best known feasible solution and let $U^{*}$ denote the objective value of problem (P’) at $\mathbf{w}^{*}$.

We are now ready to present the main steps of the proposed branch-and-bound algorithm.

Initialization. We initialize all intervals $[l^{0}_{k},u^{0}_{k}]$ for all $k=1,2,\ldots,M-1$ to be $[0,2\pi],$ i.e., set $\mathcal{A}^{0}=[0,2\pi]^{M-1}.$ In this case, problem (9) reduces to

[TABLE]

Its optimal solution $(\mathbf{w}^{0},\mathbf{c}^{0})$ and its optimal value $L^{0}$ are

[TABLE]

[TABLE]

Termination. Let $\{\mathcal{A}^{t},\mathbf{c}^{t},L^{t}\}$ be the problem instance that has the least lower bound in the problem list $\mathcal{P}.$ If

[TABLE]

where $\epsilon$ is the preselected relative error tolerance, we terminate the algorithm; otherwise we branch some interval of the above problem instance according to some rule. We can see from (11) that, both lower and upper bounds are important to avoid unnecessary branches and enumerations and good lower and upper bounds can significantly improve the efficiency of our proposed algorithm. Below, we shall introduce our branch rule as well as lower and upper bounds one by one.

Branch. Again, let $\{\mathcal{A}^{t},\mathbf{c}^{t},L^{t}\}$ be the problem instance that has the least lower bound in the problem list $\mathcal{P}$ and suppose that the stopping criterion (11) is not satisfied. In this case, we first select the interval that leads to the largest gap to be branched to smaller sub-intervals, i.e., the interval with index

[TABLE]

then we partition $\mathcal{A}^{t}$ into two sets (denoted as $\mathcal{A}_{l}^{t}$ and $\mathcal{A}_{r}^{t}$), with its $k^{*}$-th interval partitioned into two equal intervals and all the others unchanged. It follows from Proposition 2 that the ACR problems defined over the newly obtained two sets, i.e., the two children problems, are tighter than the one defined over the original set $\mathcal{A}^{t}.$ Once $\mathcal{A}^{t}$ has been branched into two sets, the problem instance defined over it will be deleted from the problem list $\mathcal{P}$ and the two children problems will be added into $\mathcal{P}$ if their optimal objective values are less than or equal to the current upper bound.

Lower Bound. Obviously, for any problem instance $\{\mathcal{A},\mathbf{c},L\},$ $L$ is a lower bound of the optimal value of the original problem (P’) defined over $\mathcal{A}.$ Therefore, the smallest lower bound among all bounds is a lower bound of the optimal value of the original problem. This statement will be formally summarized as Lemma 1 and proved in Section III-C.

Upper Bound. An upper bound of the original problem (P’) can be obtained by appropriately scaling the solution of any ACR problem instance. More specifically, let $\{\mathbf{w}^{t},\mathbf{c}^{t}\}$ be the solution of problem ACR($\mathcal{A}^{t}$). Then,

[TABLE]

is feasible to the original problem and $\|\hat{\mathbf{w}}^{t}\|^{2}$ is an upper bound of the original problem. In our proposed algorithm, the upper bound ${U^{t}}$ is the best objective values at all of the known feasible solutions at the $t$-th iteration.

By judiciously combining the above main steps, we can obtain our proposed branch-and-bound algorithm for solving problem (P) (equivalent to problem (P’)). The pseudo-code of our proposed algorithm can be found below. We will call the algorithm ACR-BB for short from now on. To make the ACR-BB algorithm more clear, an illustration on how it works is given in Appendix A.

We emphasize again that both lower and upper bounds play important roles in the efficiency of the proposed ACR-BB algorithm, because good lower and upper bounds can effectively detect inactive subproblems and avoid unnecessary branches and enumerations. Lines 9 and 11 of the proposed ACR-BB algorithm ensure that it will never branch and enumerate subproblems with the lower bound $L$ being greater than or equal to the upper bound ${U^{t}}.$ The efficiency of the proposed ACR-BB algorithm will be shown in Section IV.

III-C Global Convergence and Worst-Case Iteration Complexity

In this subsection, we present some theoretical results of our proposed ACR-BB algorithm. We first define the $\epsilon$-optimal solution of problem (P).

Definition 2** ($\epsilon$-Optimal Solution).**

Given any $\epsilon>0,$ a feasible point $\mathbf{w}$ is called an $\epsilon$-optimal solution of problem (P) if it satisfies

[TABLE]

where $\nu^{*}$ is the optimal value of problem (P).

The following Lemma 1 shows that the sequence $\left\{L^{t}\right\}$ generated by the ACR-BB algorithm is a lower bound of the optimal value of the original problem.

Lemma 1**.**

For any given instance of problem (P), let $\nu^{*}$ be its optimal value. Then, we have

[TABLE]

Proof.

The optimal value of problem (9) is always positive because zero is not a feasible solution to it. Hence, $L^{t}>0$ for all $t\geq 1.$ Next, we show $L^{t}\leq\nu^{*}$ for all $t\geq 1.$ At the beginning of the $t$-th iteration of the ACR-BB algorithm, the set $\mathcal{A}^{0}$ has been partitioned into $t$ small subsets, and the global solution of the given problem instance must lie in one of them (we denote the corresponding subset as $\mathcal{A}^{*}$). Then, the optimal value of problem ACR($\mathcal{A}^{*}$) must be less than or equal to $\nu^{*}.$ Moreover, since $L^{t}$ is the smallest lower bound of all subproblems in $\mathcal{P}$ at the $t$-th iteration, it follows that $L^{t}$ is less than or equal to the optimal value of problem ACR($\mathcal{A}^{*}$). Therefore, we get $L^{t}\leq\nu^{*}$ for all $t\geq 1.$ ∎

From lemma 1, we immediately get

[TABLE]

This further shows that, if the ACR-BB algorithm terminates, i.e., condition (11) is satisfied, the returned solution $\mathbf{w}^{*}$ by the algorithm is an $\epsilon$-optimal solution of problem (P).

The following Lemma 2 shows that the ACR-BB algorithm will terminate.

Lemma 2**.**

For any given instance of problem (P), let $\{\mathcal{A}^{t},\mathbf{c}^{t},L^{t}\}$ be the subproblem chosen in Line 9 and let $k^{*}$ be the index chosen in Line 14. If

[TABLE]

where

[TABLE]

*then condition (11) holds true and the ACR-BB algorithm will terminate in Line 12. *

Proof.

It follows from (14) and Proposition 2 that

[TABLE]

By the above inequality and the definition of $k^{*}$ (see Line 14 of the ACR-BB algorithm), we obtain

[TABLE]

Let $(\mathbf{w}^{t},\mathbf{c}^{t})$ be the solution of problem ACR($\mathcal{A}^{t}$). Then, it follows from (12) and (16) that the scaled feasible solution $\hat{\mathbf{w}}^{t}=\mathrm{Scale}(\mathbf{w}^{t},\mathbf{c}^{t})$ satisfies

[TABLE]

Moreover, since ${U^{t}}$ is the objective value at the best known feasible solution at the $t$-th iteration, we get

[TABLE]

Now, we can use the fact $L^{t}=\|\mathbf{w}^{t}\|^{2}$ and (17) to obtain

[TABLE]

The proof is completed. ∎

Based on Lemmas 1 and 2, we obtain the main result of this subsection.

Theorem 1** (Iteration Complexity).**

For any given $\epsilon>0$ and any given instance of problem (P) with $M$ users, the ACR-BB algorithm will return an $\epsilon$-optimal solution of the given instance within at most

[TABLE]

iterations, where $\delta$ is defined in (15).

Proof.

It follows from Lemmas 1 and 2 that, if (14) is satisfied, the ACR algorithm will terminate and return an $\epsilon$-optimal solution of the given problem instance. To show the theorem, it remains to show that the algorithm will terminate within $T$ iterations, where $T$ is defined in (18). Next, we show this based on the contradiction principle.

Suppose that the algorithm does not terminate within $T$ iterations. This fact, together with Lemma 2, implies that the interval that is chosen to be partitioned at the $t$-th iteration must satisfy $u^{t}_{k^{*}}-l^{t}_{k^{*}}>2\delta$ for all $t=1,2,\ldots,T.$ Then, after the partition, the width of the two sub-intervals $[l_{k^{*}}^{t},z_{k^{*}}^{t}]$ and $[z_{k^{*}}^{t},u_{k^{*}}^{t}]$ is greater than $\delta.$ Based on this, we can conclude that, for each subset $\mathcal{A}=\prod_{k=1}^{M-1}[l_{k},u_{k}]$ partitioned from the original set $\mathcal{A}^{0}$, there holds $u_{k}-l_{k}>\delta$ for all $k=1,...,M-1.$ Hence, the volume of each subset $\mathcal{A}$ is not less than $\delta^{M-1}$ and the total volume of all $T$ subsets is not less than $T\delta^{M-1}.$ Obviously, the volume of $\mathcal{A}^{0}$ is $(2\pi)^{M-1}.$ By the choice of $T,$ we get $T\delta^{M-1}>(2\pi)^{M-1},$ which further implies that the total volume of all $T$ subsets is greater than the one of the original set $\mathcal{A}^{0}.$ This is a contradiction. Hence, the algorithm will terminate within at most $T$ iterations. ∎

Some remarks on Theorem 1 are in order. First, from Theorem 1, we can immediately obtain the following convergence result, which shows that both the sequences of the upper bounds and the lower bounds generated by the ARC-BB algorithm with $\epsilon=0$ converge to the optimal objective value of problem (P). Therefore, the ACR-BB algorithm with $\epsilon=0$ indeed is a global algorithm and is able to find the global solution of problem (P).

Corollary 1** (Global Optimality).**

For any given instance of problem (P), let $\nu^{*}$ be its optimal value and let $\left\{U^{t}\right\}$ and $\left\{L^{t}\right\}$ be the iterates generated by the ARC-BB algorithm with $\epsilon=0$. Then, we have $U^{t}\rightarrow\nu^{*}$ and $L^{t}\rightarrow\nu^{*}.$

Proof.

Theorem 1 shows that, for any given $\epsilon>0,$ there exists an integer $T\geq 1$ (given in (18)) such that $(U^{t}-L^{t})/L^{t}\leq\epsilon$ for all $t\geq T.$ This statement is equivalent to ${(U^{t}-L^{t})}/{L^{t}}\rightarrow 0.$ Since the sequence $\left\{L^{t}\right\}$ is uniformly bounded away from zero for any given instance of problem (P), we further obtain ${U^{t}-L^{t}}\rightarrow 0.$ This, together with Lemma 1 and the definition of $U^{t}$, immediately implies the desired results. The proof is completed. ∎

In practice, we need to preselect a positive relative error tolerance $\epsilon$ in our proposed ACR-BB algorithm, as most of iterative optimization algorithms [34, Chapter 9]. Theorem 1 shows that the total number of iterations for our proposed ACR-BB algorithm to return an $\epsilon$-optimal solution of any instance of problem (P) is exponential with respect to the number of receivers $M.$ The iteration complexity of our proposed algorithm seems high at first sight. However, as will be shown in Section IV, its practical iteration complexity is actually significantly less than the worst-case bound in (18). It is also worthwhile remarking that there is no polynomial time algorithm which can globally solve the problem (unless P=NP) [5], because the problem is NP-hard.

IV Simulation Results

In this section, we present some numerical simulation results of our proposed ACR-BB algorithm for solving problem (P). More specifically, we first present some simulation results to show the convergence behaviors of the proposed algorithm in Section IV-A. Then, we show the efficiency of the proposed algorithm by comparing it with a state-of-the-art general-purpose global optimization solver called Baron [30, 31]. Last but not least, we use the proposed algorithm as the benchmark to evaluate the performance of two state-of-the-art algorithms [4, 12] for solving the same problem.

In all of our simulations, the channel vectors $\mathbf{h}_{k}$ are randomly generated according to the distribution $\mathcal{C}\mathcal{N}(\mathbf{0},\mathbf{I}_{N})$ as done in [4, 12, 27]. We performed all numerical experiments on a PC with a 3.40-GHz Intel Core i7-2600 processor with access to 4 GB of RAM. We implemented our proposed ACR-BB algorithm in Matlab 7.10. We use the built-in function “quadprog” in Matlab to solve all linearly constrained convex quadratic programs, and use Sedumi [35] to solve all SDPs.

IV-A Convergence Behaviors of the ACR-BB Algorithm

In this subsection, we generate a problem instance with $(N,M)=(4,40),$ apply our proposed ACR-BB algorithm to solve it, and study the convergence behaviors of our proposed algorithm, i.e., the convergence behaviors of the lower bounds $\left\{L^{t}\right\}$ and the upper bounds $\left\{U^{t}\right\}.$ In this simulation, we set $\epsilon=1$e$-7$ in our proposed algorithm such that the algorithm can find the “real” optimal solution when it terminates. We denote the objective value at the returned solution by $\bar{\nu}$ and define

[TABLE]

In the above, $E_{1}^{t}$ and $E_{2}^{t}$ are the relative errors of the lower bound $L^{t}$ and the upper bound $U^{t}$, respectively; $E_{3}^{t}$ is an easily computable upper bound of $E_{1}^{t}+E_{2}^{t}$. In this paper, we will call all of them relative errors. The above three relative errors versus the number of iterations are illustrated as Fig. 2.

It can be seen from Fig. 2 that $E_{2}$ converges much faster than $E_{1}$ and $E_{3}$ at the beginning stage of our proposed algorithm. For instance, $E_{2}$ becomes smaller than $1$e$-1$ at the $6702$-th iteration, while $E_{1}$ and $E_{3}$ become smaller than $1$e$-1$ at the $9653$-th and $12221$-th iteration, respectively. The above results imply that an $0.1$-optimal solution has been found at the $6702$-th iteration but it is verified (to be an $0.1$-optimal solution) at the $12221$-th iteration. These results show that:

The upper bounds $\left\{U^{t}\right\}$ converge faster than the lower bounds $\left\{L^{t}\right\}$ at the beginning stage of our proposed algorithm; and more importantly,

• -

A low-accuracy solution can be more easier found than verified by our proposed algorithm.

We can also see from Fig. 2 that it takes our proposed ACR-BB algorithm $3339$ iterations to reduce $E_{3}$ from $1$e$-1$ to $1$e$-3$ and $96$ iterations to reduce it from $1$e$-3$ to $1$e$-6,$ respectively. This shows that the performance of our proposed algorithm in terms of the number of iterations is not sensitive to the choice of the relative error tolerance $\epsilon$ (when it is smaller than $10^{-3}$).

In the remaining part of this section, we set $\epsilon=5$e$-3.$

IV-B Efficiency of the ACR-BB Algorithm

In this subsection, we study the efficiency of our proposed ACR-BB algorithm for solving problem (P).

We first compare our proposed algorithm with Baron444To the best of our knowledge, our proposed ACR-BB algorithm is the first tailored global algorithm for problem (P) and there is no existing global algorithms specially designed for the problem that we can compare our proposed algorithm with. [30, 31], which is a state-of-the-art general-purpose global optimization solver and has been widely applied to solve problems arising from various applications [36, 37]. Although both of our proposed algorithm and Baron lie in the branch-and-bound framework, the difference between them is remarkable. First, the branch strategy in the two algorithms is different, i.e., our proposed algorithm branches in the (argument) ranges of $\arg(c_{k})$ but Baron transforms problem (P) into a real (nonconvex) quadratic program by representing the real and imaginary parts of each complex variable with two independent real variables and branches in the ranges of these real variables. Second, convex relaxation used in the two algorithms is also different, i.e., our proposed algorithm is based on ACR (9), which is a linearly constrained quadratic program, but Baron is based on linear programming relaxation [30, 31]. Since the lower bound based on ACR (9) is generally better than the one based on linear programming relaxation, our proposed algorithm is much more efficient than Baron, as shown below.

We apply our proposed ACR-BB algorithm and Baron to solve 10 randomly generated problem instances with $(N,M)=(2,8).$ The comparison results are summarized in Table I, where the first column shows the IDs of the corresponding problem instances, the second column shows the objective values and the number of iterations and CPU time results obtained by our proposed algorithm, the third column shows the objective values and the number of iterations and CPU time results obtained by Baron, and the last row shows the average results over the $10$ instances. We can observe from Table I that our proposed ACR-BB algorithm performs $366$ times faster than Baron (on average) to find the same solutions. In fact, we have tried to use Baron to solve problem instances with larger $N$ and/or $M.$ Unfortunately, we found that Baron fails to solve all problem instances with $N\geq 4$ and $M\geq 8$ within $10$ minutes. These observations demonstrate that our specially designed ACR-BB algorithm achieves significantly higher efficiency on globally solving problem (P) than Baron.

Next, we present more numerical results on applying our proposed ACR-BB algorithm to solve problem (P) with larger $N$ and/or $M$ without comparing it with other general-purpose global optimization solvers. All the results to be shown from now on are obtained by averaging over or choosing from $50$ randomly generated problem instances for each pair $(N,M)$. Table II reports the number of iterations and CPU time results, where the first column shows the setup of the problem, the second column shows the average results (over the $50$ problem instances), and the third column shows the worst-case result (among the $50$ problem instances). Figs. 3 and 4 plot the average CPU time and the average number of iterations versus the number of users with $N=4.$

We can observe from the above simulation results (Table II) that our proposed ACR-BB algorithm is capable of finding the global solution (with the relative error being not greater than $5$e$-3$) within several seconds when $N$ is small and $M$ is not too large. In particular, the proposed algorithm can solve all generated problem instances with $N=2$ and $M$ ranging from $8$ to $64$ within $2.4$ seconds. When $N=4,$ the proposed algorithm can solve problem instances with $M\leq 16$ within $20$ seconds on average and problem instances with $M$ ranging from $24$ to $40$ within $185$ seconds on average. The proposed algorithm can solve all $50$ problem instances with $N=4$ and $M$ ranging from $8$ to $40$ within $7$ minutes. Moreover, we can observe from Fig. 4 that the number of iterations of our proposed algorithm is extremely smaller than the worst-case bound in (18). These results demonstrate that our proposed ACR-BB algorithm is highly efficient for problem (P) when $N$ is small and $M$ is not too large. The high efficiency of our proposed ACR-BB algorithm is due to the newly developed ACR (9), which effectively avoids unnecessary branches that do not contain global optimal solutions.

When $N$ and $M$ (especially $N$) are large, it takes our proposed algorithm relatively long time (and a relatively large number of iterations) to find the global solution. For instance, our proposed algorithm needs $1396.9$ seconds on average to solve problem instances with $(N,M)=(8,16).$ Even though the efficiency of our proposed algorithm is not very high when $N$ and $M$ (especially $N$) are large, our proposed algorithm is still useful in some scenarios. For instance, in the scenario where the channel is stationary/constant and computational efficiency is not a big issue, our proposed algorithm can find the global solution and thus can provide the best performance in terms of transmitting the minimum total power. Another key application of our proposed ACR-BB algorithm is that it provides an important benchmark for performance evaluation of other heuristic/local optimization algorithms for the same problem, as will be done in the next subsection.

IV-C Performance Evaluation of Existing Heuristic/Local Optimization Algorithms

In this subsection, we evaluate the performance of two state-of-the-art algorithms for problem (P), the SDR algorithm [4] and the SLA algorithm [12], by using our proposed ACR-BB algorithm as the benchmark. Our evaluation metrics are the objective values and the relative gaps. Note that the objective function in our problem is the total transmission power and hence the smaller the objective value returned by an algorithm, the better the algorithm. The relative gap of an algorithm is defined as $(\nu_{A}-\bar{\nu})/\bar{\nu}$, where $\nu_{A}$ is the objective value returned by the algorithm and $\bar{\nu}$ is the objective value returned by our proposed algorithm.

In our simulations, the SDR algorithm and the SLA algorithm are implemented by following the descriptions in [4] and [12], respectively. In particular, $1000$ points are randomly generated based on the solution of the SDR (3) to obtain the final approximate solution in [4] and also $1000$ points are randomly generated to obtain a good initial point for the SLA algorithm in [12]. Recall that we have generated ${900}$ instances of problem (P) in Table II and obtained their optimal solution (with the relative errors being not greater than $5$e$-3$) by using our proposed ACR-BB algorithm. We apply the SDR algorithm and the SLA algorithm to solve these instances and evaluate their performance.

Table III reports the average objective values and average relative gaps of the SDR algorithm and the SLA algorithm. We can see from Table III that the average relative gaps of both algorithms are small when $N$ and $M$ (especially $N$) are small and the average relative gaps of the SLA algorithm is smaller than the ones of the SDR algorithm. This shows that, the quality of the returned solutions by both algorithms is high when $N$ and $M$ (especially $N$) are small and the quality of the returned solution by the SLA algorithm is generally better than the one of the returned solution by the SDR algorithm. However, as $N$ and/or $M$ increase, the average relative gaps of both algorithms quickly become large. For instance, when $(N,M)=(4,40)$, the relative gap of the SDR algorithm is $69.64\%$; when $(N,M)=(4,32)$, the relative gap of the SLA algorithm is $19.48\%$. These results show that the quality of the solutions obtained by the two algorithms are not good when $N$ and/or $M$ are large and the performance of both algorithms degrades quickly as $N$ and/or $M$ increase. This makes sense, because both of the algorithms are not global algorithms and they are more likely to get stuck at a local solution when $N$ and/or $M$ are large.

To study the robustness of the two algorithms, we list the worst-case relative gaps of the two algorithms in Table IV. We can see from Table IV that the worst-case relative gap of the SLA algorithm exceeds $70\%$ in two setups $(4,40)$ and $(6,24);$ and the worst-case relative gap of the SDR algorithm exceeds $100\%$ in five setups $(4,32),~{}(4,40),~{}(6,16),~{}(6,24),$ and $(8,16).$ These results clearly demonstrate that both of the algorithms are not robust (especially when $N$ and/or $M$ are relatively large) because their performance for some problem instance is bad and the SLA algorithm is more robust than the SDR algorithm.

We also test the probability that the two algorithms can find the global solution. We call that an algorithm finds the global solution of a problem instance if the relative gap is less than or equal to $1$e$-3.$ The probability that an algorithm finds the global solution is defined as the ratio of the number of instances of which it finds the global solution and the total number of instances of which it is applied to solve. The latter is $50$ in our case. The probability that the two algorithms find the global solution is reported in Table V. As can be seen from the table, the probability that the SLA algorithm and the SDR algorithm can find the global solution is (relatively) large when both $N$ and $M$ are small and the probability that the SLA algorithm finds the global solution is generally much larger than the one that the SDR algorithm finds the global solution. However, the probability that the SLA algorithm and the SDR algorithm find the global solution decreases quickly as $N$ and/or $M$ increase. These results are intuitive and agree well with our previous simulation results.

In summary, we can make the following conclusions on the performance of the two state-of-the-art algorithms for problem (P), the SLA algorithm and the SDR algorithm, by using our proposed ACR-BB algorithm as benchmark:

both of the algorithms perform well when $N$ and $M$ are small;

• -

the performance of the SLA algorithm is generally better than the one of the SDR algorithm in terms of (average and worst-case) relative gaps, robustness, and probability of finding the global solutions; and

• -

the performance of both of the algorithms degrades quickly as $N$ and/or $M$ increase.

It is worthwhile remarking that both of the SLA algorithm and the SDR algorithm are more efficient than our proposed ACR-BB algorithm. For instance, to solve problem (P) with $N=8$ and $M=16,$ both of the algorithms need less than 0.25 seconds (on average) but it takes our proposed algorithm 1396.9 seconds. However, our proposed ACR-BB algorithm is a global algorithm which is guaranteed to find the global solution of the problem within any given relative error tolerance whereas the SLA algorithm is a local optimization algorithm and the SDR algorithm is an approximation algorithm and both of them cannot guarantee global optimality of their returned solutions.

V Conclusion

In this paper, we proposed the ACR-BB algorithm for solving the single-group multicast beamforming problem. The proposed algorithm is guaranteed to find the global solution of the problem within any given relative error tolerance. To the best of our knowledge, the proposed ACR-BB algorithm is the first specially designed global algorithm for solving the single-group multicast beamforming problem. The proposed algorithm is based on the branch-and-bound strategy as well as a newly developed argument cut technique. Simulation results show that the proposed algorithm significantly outperforms a state-of-the-art general-purpose solver Baron. The high efficiency of the proposed algorithm is mainly due to the use of the argument cut based relaxation, which can potentially be used in other related algorithms/problems. An important role that the proposed algorithm can play is that it can be used as a benchmark for performance evaluation of other heuristic/local optimization algorithms for the same problem.

Acknowledgment

The authors would like to thank three anonymous reviewers for their constructive comments, which significantly improved the quality and presentation of the paper.

Appendix A An illustration of the ACR-BB algorithm

To make the proposed ACR-BB algorithm clear, an illustration of applying it to solve the following problem instance with $(N,M)=(2,3)$ is given:

[TABLE]

[TABLE]

[TABLE]

We set $\epsilon=0.1.$

At the [math]-th iteration, we initialize $\mathcal{A}^{0}=\left[0,2\pi\right]^{2}.$ The ACR-BB algorithm solves problem ACR($\mathcal{A}^{0}$) and obtains its optimal solution

[TABLE]

[TABLE]

[TABLE]

Then, $\mathbf{w}^{0}$ is scaled to obtain the feasible point

[TABLE]

Now,

[TABLE]

and

[TABLE]

At the $1$-th iteration, we have

[TABLE]

Since

[TABLE]

the ACR-BB algorithm starts executing Line 14. Recall $|\mathbf{c}^{1}|=\left[1.8330,0.7270\right],$ where $|\cdot|$ denotes the component-wise absolute value operator. We have

[TABLE]

[TABLE]

The ACR-BB algorithm solves problems ACR($\mathcal{A}_{l}^{1}$) and ACR($\mathcal{A}_{r}^{1}$) and obtains their optimal solutions

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

Then, $\mathbf{w}_{l}^{1}$ and $\mathbf{w}_{r}^{1}$ are scaled to obtain two feasible points

[TABLE]

[TABLE]

respectively. Now, since $\|\hat{\mathbf{w}}_{l}^{1}\|^{2}=1.1015>{U^{1}}$ and $\|\hat{\mathbf{w}}_{r}^{1}\|^{2}={0.8573}={U^{1}},$ ${U^{1}}$ and $\mathbf{w}^{*}$ are not updated; since $L_{l}^{1}<{U^{1}}$ and $L_{r}^{1}<{U^{1}},$ we have

[TABLE]

At the $2$-th iteration, we have ${U^{2}=U^{1}};$ since $L_{r}^{1}<L_{l}^{1},$ we have

[TABLE]

Since

[TABLE]

the ACR-BB algorithm starts executing Line 14. Recall $|\mathbf{c}^{2}|=\left[1.8330,0.7270\right].$ We have

[TABLE]

[TABLE]

The ACR-BB algorithm solves problems ACR($\mathcal{A}_{l}^{2}$) and ACR($\mathcal{A}_{r}^{2}$) and obtains their optimal solutions

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

Then $\mathbf{w}_{l}^{2}$ and $\mathbf{w}_{r}^{2}$ are scaled to obtain two feasible points

[TABLE]

[TABLE]

respectively. Now, since $\|\hat{\mathbf{w}}_{l}^{2}\|^{2}=0.8023<{U^{2}}$ and $\|\hat{\mathbf{w}}_{r}^{2}\|^{2}=0.7811<{U^{2}},$ $U^{2}$ is updated to $0.7811$ and $\mathbf{w}^{*}$ is updated to $\hat{\mathbf{w}}_{r}^{2}$; since $L_{l}^{2}<{U^{2}}$ and $L_{r}^{2}<{U^{2}},$ we have

[TABLE]

At the $3$-th iteration, we have ${U^{3}=U^{2}};$ since $L_{l}^{2}<L_{l}^{1}<L_{r}^{2},$ we have

[TABLE]

Since

[TABLE]

the ACR-BB algorithm starts executing Line 14. Because $|\mathbf{c}^{3}|=\left[1.8500,0.7517\right],$ we have

[TABLE]

[TABLE]

The ACR-BB algorithm solves problems ACR($\mathcal{A}_{l}^{3}$) and ACR($\mathcal{A}_{r}^{3}$) and obtains their optimal solutions

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

Then $\mathbf{w}_{l}^{3}$ and $\mathbf{w}_{r}^{3}$ are scaled to obtain two feasible points

[TABLE]

[TABLE]

respectively. Now, since $\|\hat{\mathbf{w}}_{l}^{3}\|^{2}=0.5453<{U^{3}}$ and $\|\hat{\mathbf{w}}_{r}^{3}\|^{2}=0.5072<{U^{3}},$ ${U^{3}}$ is updated to $0.5072$ and $\mathbf{w}^{*}$ is updated to $\hat{\mathbf{w}}_{r}^{3};$ since $L_{l}^{3}<{U^{3}}$ and ${L_{r}^{3}=U^{3}},$ we have

[TABLE]

At the $4$-th iteration, we have $U^{4}=U^{3};$ since $L_{l}^{3}<L_{l}^{1}<L_{r}^{3}<L_{r}^{2},$ we have

[TABLE]

Since

[TABLE]

the ACR-BB algorithm terminates and return

[TABLE]