Distributed Multi-Stream Beamforming in MIMO Multi-Relay Interference Networks

TL;DR

This paper proposes a distributed beamforming optimization method for multi-stream MIMO relay interference networks, achieving power minimization with convergence guarantees and improved efficiency over existing solutions.

Contribution

It introduces a novel ADMM-based distributed algorithm for joint transmit and relay beamforming optimization in interference networks.

Findings

ADMM algorithm converges to the optimal solution.

Proposed method outperforms existing solutions in convergence rate and message exchange.

Joint transmit and relay beamforming further reduces power consumption.

Abstract

In this paper, multi-stream transmission in interference networks aided by multiple amplify-and-forward (AF) relays in the presence of direct links is considered. The objective is to minimize the sum power of transmitters and relays by beamforming optimization under the stream signal-to-interference-plus-noise-ratio (SINR) constraints. For transmit beamforming optimization, the problem is a well-known non-convex quadratically constrained quadratic program (QCQP) that is NP-hard to solve. After semi-definite relaxation (SDR), the problem can be optimally solved via alternating direction method of multipliers (ADMM) algorithm for distributed implementation. Analytical and extensive numerical analyses demonstrate that the proposed ADMM solution converges to the optimal centralized solution. The convergence rate, computational complexity, and message exchange load of the proposed algorithm outperforms the existing solutions. Furthermore, by SINR approximation at the relay side, distributed joint transmit and relay beamforming optimization is also proposed that further improves the total power saving at the cost of increased complexity.

Figures

Tables

Table 1. Table I : LIST OF CHANNEL NOTATIONS.

${𝐉}_{ki}$	Channel from the $i^{th}$ transmitter to the $k^{th}$ receiver
${𝐆}_{kr}$	Channel from the $r^{th}$ relay to the $k^{th}$ receiver
${𝐇}_{ri}^{\prime \prime }$	Channel from the $i^{th}$ transmitter to the $r^{th}$ relay
${𝐇}_{k.i}^{\prime }$	Effective channel from the $i^{th}$ transmitter to the $k^{th}$ receiver
	through all relays
${𝐇}_{ki}$	Aggregate effective channel from the $i^{th}$ transmitter to the
	$k^{th}$ receiver, i.e., stacked matrix of ${𝐉}_{ki}$ and ${𝐇}_{k.i}^{\prime }$

Table 2. Table II : LIST OF SOME SIMPLIFICATIONS.

${𝐑}_{rk}$	${𝐇}_{rk}^{\prime \prime H}{𝐅}_r^H{𝐅}_r{𝐇}_{rk}^{\prime \prime }$
${𝐑}_{.k}^{\prime }$	${\sum }_{r=1}^R{𝐇}_{rk}^{\prime \prime H}{𝐅}_r^H{𝐅}_r{𝐇}_{rk}^{\prime \prime }$
${𝐑}_{.k}$	${𝐈}_{M_k}+{𝐑}_{.k}^{\prime }$
${𝐑}_{k,l}^i$	${𝐇}_{ki}^H{\overline{𝐯}}_{k,l}{\overline{𝐯}}_{k,l}^H{𝐇}_{ki}$

Table 3. Table III : AVERAGE CPU TIMES OF THE PROPOSED AND DISTRIBUTED ALGORITHMS FOR THE K ​ 3 ​ - M ​ 6 ​ - N ​ 8 ​ - R ​ 14 K 3 - M 6 - N 8 - R 14 \text{K}3\text{-}\text{M}6\text{-}\text{N}8\text{-}\text{R}14 NETWORK.

Algorithm	Proposed	ADMM-BG	ADAL
Average CPU time (sec.)	1.48	1.98	2.13

Table 4. Table IV : PERFORMANCES OF DISTRIBUTED TRANSMIT VS. DISTRIBUTED JOINT TRANSMIT AND RELAY BEAMFORMING FILTER DESIGNS FOR THE K ​ 3 ​ - M ​ 10 ​ - N ​ 10 ​ - R ​ 3 K 3 - M 10 - N 10 - R 3 \text{K}3\text{-}\text{M}10\text{-}\text{N}10\text{-}\text{R}3 NETWORK AT 21 ​ - ​ 21 21 - 21 21\text{-}21 dB.

	Dist.	Dist. joint
	$\left\{{𝐔}_k\right\}$	$\{{𝐔}_k,{𝐅}_r\}$
Average sum-SINR (dBm)	28.32	28.32
Average total power (dBm)	46.05	45.23
Average iteration number	16	23
​​​​Total complexity per iteration​​​​	​ $𝒪\left(90\right)$	​$𝒪\left(90+450=540\right)$
Average total complexity	​​​​​ $𝒪\left(16\cdot 90=1440\right)$ ​​​	​​​​​ $𝒪\left(23\cdot 540=12420\right)$​​​

Equations

$$p_{i}=E(||\mathbf{x}_{i}||^{2})=\text{tr}\big{(}\mathbf{U}_{i}\mathbf{U}_{i}^{H}\big{)}\leq p_{i}^{\text{max}},$$

$${\mathbf{y}}_r=\sum _{i=1}^K{\mathbf{H}}_{ri}^{\prime \prime }{\mathbf{x}}_i+{\mathbf{n}}_r\text{ and }{\mathbf{y}}_k(1)=\sum _{i=1}^K{\mathbf{J}}_{ki}{\mathbf{x}}_i+{\mathbf{n}}_k(1),$$

$$p_{r}=\text{tr}\Big{(}\!\!\sum_{i=1}^{K}\!\!\mathbf{U}_{i}^{H}\mathbf{H}_{ri}^{\prime\prime H}\mathbf{F}_{r}^{H}\mathbf{F}_{r}\mathbf{H}^{\prime\prime}_{ri}\mathbf{U}_{i}\Big{)}+\sigma_{r}^{2}\text{tr}\Big{(}\mathbf{F}_{r}\mathbf{F}_{r}^{H}\Big{)}.$$

$${\mathbf{y}}_k(2)=$$

$$+\sum _{r=1}^R{\mathbf{G}}_{kr}{\mathbf{F}}_r\sum _{\begin{array}{c}j=1,j\ne k\end{array}}^K{\mathbf{H}}_{rj}^{\prime \prime }{\mathbf{x}}_j+\sum _{r=1}^R{\mathbf{G}}_{kr}{\mathbf{F}}_r{\mathbf{n}}_r+{\mathbf{n}}_k(2),$$

$${\mathbf{H}}_{k.i}^{\prime }\triangleq \sum _{r=1}^R{\mathbf{G}}_{kr}{\mathbf{F}}_r{\mathbf{H}}_{ri}^{\prime \prime }.$$

$${\mathbf{y}}_k(2)={\mathbf{H}}_{k.k}^{\prime }{\mathbf{x}}_k+\sum _{\begin{array}{c}j=1,j\ne k\end{array}}^K{\mathbf{H}}_{k.j}^{\prime }{\mathbf{x}}_j+\sum _{r=1}^R{\mathbf{G}}_{kr}{\mathbf{F}}_r{\mathbf{n}}_r+{\mathbf{n}}_k(2).$$

$$\mathbf{H}_{ki}\triangleq\left[\begin{array}[]{c}\!\!\mathbf{J}_{ki}\!\!\\ \!\!\mathbf{H}^{\prime}_{k.i}\!\!\\ \end{array}\right]\!\text{ and }\mathbf{n}_{k}\triangleq\!\left[\begin{array}[]{c}\!\!\mathbf{n}_{k}(1)\!\!\\ \!\!\sum_{r=1}^{R}\mathbf{G}_{kr}\mathbf{F}_{r}\mathbf{n}_{r}+\mathbf{n}_{k}(2)\!\!\\ \end{array}\right].$$

$$\mathbf{y}_{k}\triangleq\left[\begin{array}[]{c}\mathbf{y}_{k}(1)\\ \mathbf{y}_{k}(2)\\ \end{array}\right]=\mathbf{H}_{kk}\mathbf{U}_{k}\mathbf{s}_{k}+\sum_{\begin{subarray}{c}j=1,j\neq k\end{subarray}}^{K}\mathbf{H}_{kj}\mathbf{U}_{j}\mathbf{s}_{j}+\mathbf{n}_{k}.$$

$${\text{SINR}}_{k,l}=\frac{{\tilde{\zeta }}_{k,l}^{k,l}}{{\sum }_{(j,m)\ne (k,l)}{\tilde{\zeta }}_{k,l}^{j,m}+{\sigma }_{n_{k,l}}^2},$$

$$\mathbf{R}_{n_{k}}\!\!\triangleq\!\!\left(\!\begin{array}[]{cc}\!\!\sigma_{k}^{2}(1)\mathbf{I}_{M_{k}}\!\!&\!\!\!\bm{0}\!\!\!\\ \!\!\bm{0}\!\!&\!\!\!\sum_{r=1}^{R}\!\sigma_{r}^{2}\mathbf{G}_{kr}\mathbf{F}_{r}\mathbf{F}_{r}^{H}\mathbf{G}_{kr}^{H}\!+\!\sigma_{k}^{2}(2)\mathbf{I}_{M_{k}}\!\!\\ \end{array}\!\right)$$

$$\sum _{k=1}^K\sum _{l=1}^{d_k}{\mathbf{u}}_{k,l}^H{\mathbf{R}}_{.k}{\mathbf{u}}_{k,l},$$

$${\mathbf{R}}_{rk}\triangleq {\mathbf{H}}_{rk}^{\prime \prime H}{\mathbf{F}}_r^H{\mathbf{F}}_r{\mathbf{H}}_{rk}^{\prime \prime }.$$

$${\text{SINR}}_{k,l}=\frac{{\mathbf{u}}_{k,l}^H{\mathbf{R}}_{k,l}^k{\mathbf{u}}_{k,l}}{{\sum }_{(j,m)\ne (k,l)}{\mathbf{u}}_{j,m}^H{\mathbf{R}}_{k,l}^j{\mathbf{u}}_{j,m}+{\sigma }_{n_{k,l}}^2},$$

$$\underset{\{{\mathbf{u}}_{k,l}{\}}_{\forall k\in \mathcal{K},\forall l\in {\mathcal{L}}_k}}{\min }$$

$${\text{SINR}}_{k,l}\ge {\gamma }_{k,l},\text{leavevmode}\forall k\in \mathcal{K},\text{leavevmode}\forall l\in {\mathcal{L}}_k$$

$$\underset{\{{\mathbf{X}}_{k,l}{\}}_{\forall k,l}}{\min }$$

$$\frac{1}{{\gamma }_{k,l}}\text{tr}({\mathbf{X}}_{k,l}{\mathbf{R}}_{k,l}^k)-\sum _{(j,m)\ne (k,l)}\text{tr}({\mathbf{X}}_{j,m}{\mathbf{R}}_{k,l}^j)\ge {\sigma }_{n_{k,l}}^2,$$

$$\forall k\in \mathcal{K},\forall l\in {\mathcal{L}}_k$$

$${\mathbf{X}}_{k,l}\in {\mathbb{S}}_{+}^{M_k},\text{leavevmode}\forall k\in \mathcal{K},\text{leavevmode}\forall l\in {\mathcal{L}}_k$$

$$\text{rank}({\mathbf{X}}_{k,l})=1,\text{leavevmode}\forall k\in \mathcal{K},\text{leavevmode}\forall l\in {\mathcal{L}}_k.$$

$$\underset{\mathbf{x}\in {\mathbb{C}}^m,\mathbf{y}\in {\mathbb{C}}^n}{\min }$$

$$\mathbf{Ax}+\mathbf{By}=\mathbf{c}$$

$$\mathbf{x}\in {\mathcal{S}}_1,\mathbf{y}\in {\mathcal{S}}_2,$$

$${\mathfrak{L}}_{\rho }(\mathbf{x},\mathbf{y},\mathbf{z})=$$

$$+\frac{\rho }{2}\Vert \mathbf{Ax}+\mathbf{By}-\mathbf{c}{\Vert }_2^2,$$

$${\mathbf{x}}^{s+1}=$$

$${\mathbf{y}}^{s+1}=$$

$${\mathbf{z}}^{s+1}=$$

$${\zeta }_{k,l}\triangleq$$

$${\zeta }_{k,l}^{‵}\triangleq$$

$$\underset{\{{\mathbf{X}}_{k,l},{\zeta }_{k,l},{\zeta }_{k,l}^{‵}{\}}_{\forall k,l}}{\min }$$

$${\zeta }_{k,l}+{\zeta }_{k,l}^{‵}=0,\forall k\in \mathcal{K},\forall l\in {\mathcal{L}}_k$$

$${\zeta }_{k,l}=\frac{1}{{\gamma }_{k,l}}\text{tr}({\mathbf{X}}_{k,l}{\mathbf{R}}_{k,l}^k)-{\sigma }_{n_{k,l}}^2,\forall k\in \mathcal{K},\forall l\in {\mathcal{L}}_k$$

$${\zeta }_{k,l}^{‵}=-\sum _{(j,m)\ne (k,l)}\text{tr}({\mathbf{X}}_{j,m}{\mathbf{R}}_{k,l}^j),\forall k\in \mathcal{K},\forall l\in {\mathcal{L}}_k$$

$$\text{The constraint }\text{eqref}op:0727_{1}915_{S}PSDConstraint.$$

$$\displaystyle\mathfrak{L}_{\rho}\big{(}\{\mathbf{X}_{k,l},\zeta_{k,l},\zeta_{k,l}^{\backprime},\lambda_{k,l}\}_{\forall k,l}\big{)}=\!\sum_{k=1}^{K}\sum_{l=1}^{d_{k}}\!\Big{(}\text{tr}(\mathbf{X}_{k,l}\mathbf{R}_{.k})$$

$$\displaystyle\hskip 71.13188pt+\lambda_{k,l}(\zeta_{k,l}\!+\!\zeta_{k,l}^{\backprime})\!+\!\frac{\rho}{2}(\zeta_{k,l}\!+\!\zeta_{k,l}^{\backprime})^{2}\Big{)},$$

$$\mathbf{x}=$$

Full text

I Introduction

The advent of future wireless networks bearing new components in large numbers including eNodeBs, cloud servers, relays, smart grids, massive multiple-input multiple-output (MIMO), and big data nodes, and the recent advances both in software and hardware architectures surging the applicability of parallel and distributed computations [1] have brought innovative solutions and paradigm shifts in recent years[2, Chap. 10][3, 4, 5]. However, the distributed solutions based on conventional dual decomposition and other methods lack the effectiveness on the numerical stability, fast convergence rates [6], and scalability to high dimensional problems [7] compared to alternating direction method of multipliers (ADMM) which combines the strengths of dual decomposition and augmented Lagrangian methods [8].

ADMM studies in relay [9, 10, 11] and point-to-point [12, 13, 14, 15, 16, 7] networks are limited. In [9], an improved version of ADMM is proposed for power minimization under SINR constraints in multi-cluster relay networks with single antenna nodes. In [10], max-min SINR optimization is studied for decode-and-forward (DF) relay networks with single antenna transmitters and receivers under the constraints of total transmitter power, total relay power, and total number of relays. In [11], a distributed transmit power control algorithm via ADMM is proposed for a single relay aided network. ADMM in point-to-point networks is slightly more investigated in the literature. In [7], cloud radio access networks (C-RAN) with single antenna receivers are optimized for power minimization under SINR constraints. In [12], power minimization with the worst-case SINR constraints due to the channel state information (CSI) errors in multi-cell coordinated multiple-input single-output (MISO) networks is solved via semi-definite relaxation (SDR). In [13], power minimization problem with rate constraints in wireless sensors networks is studied. In [14], beamforming design for power minimization under SINR constraints is proposed in MISO downlink systems by solving second order cone problems (SOCP). In [15], max-min flow rate optimization problem is considered for software defined radio access networks (SD-RAN) with single antenna nodes under wired and wireless link constraints. In [16], a distributed power control algorithm is proposed for the utility maximization without SINR constraints in interference networks with single antenna nodes.

The three prominent features of future high performance wireless networks are multi-stream transmissions, multi-antenna nodes, and the line-of-sights between transmitters and receivers. The existence of direct links between transmitters and receivers is an immediate outcome of deploying a large number of intermediate nodes in the network to close the distances between transmitters and receivers. In C-RANs, SD-RANs, and wireless relay networks, the wireless intermediate nodes are called radio access units [7], base stations [15], and relays, respectively. Wireless relay networks can be regarded as the wireless communications parts of C-RANs and SD-RANs, which embody wired communications parts in their architectures as well. Although the three mentioned features are critical in practical systems, many studies on relay networks are based on communications settings with lesser features due to the difficulties that arise from the coexistence of all three features[17]. Among many objective functions [18, Chap. 8], the mean square error (MSE) minimization problems are more tractable in relay networks with the three features. Nevertheless, the studies on MSE are still limited[19], and to the best of our knowledge, the MSE problem in a relay network with all the three features is studied only in [20].

The limited attributes of aforementioned researches on relay networks can diminish their applications in future practical systems. In this paper, we consider amplify-and-forward (AF) multi-relay interference networks with all the three mentioned attributes. The problem of interest in this network setting is to minimize the total power consumption of transmitters and relays with guaranteed quality of service (QoS) by distributed optimization of the transmit beamforming filters. The QoS metric chosen in this paper is stream SINR, which is interconnected to the data rate and bit error rate (BER) expressions. The problem is a member of non-convex quadratically constrained quadratic programming (QCQP) problems [21, 22], which is not directly amenable to distributed optimization due to the intricate stream SINR constraints. In particular, the SINR constraints are not in a linear form and also are not decoupled over streams for parallel and distributive implementation. We fit the problem into the ADMM framework, a potent tool for distributed optimization.

When designing a solution for the problem, the feasibility of problem, i.e., the SINR constraints (targets) must be jointly supported, must be assured in the first stage before solving the main problem in the second stage. There are two approaches to assure the success of the first stage, namely, deriving the feasibility conditions [23] and relaxing the initial conditions, e.g., searching for feasible SINR targets[24, 4] and reducing the number of users [25]. The derivation of the feasibility condition is challenging even in simpler networks. In [26], an approximate condition is derived for a multi-relay network with a single transmitter and a receiver. The feasibility search, on the other hand, can be as costly as the solving the main problem in terms of the number of iterations and computational complexity per iteration.

In this paper, we propose to use random initializations of beamforming vectors to automatically determine feasible SINR targets with a high probability. All cases where the convergence is slow, fluctuant, and infeasible, i.e., the SINR targets are infeasible, make up less than $2\%$ of the simulations presented in this paper. The elimination of these mentioned cases is beneficial in two important applications: (1) Accurate and extensive cross-analyses of crucial network parameters, and benchmarking with competitive schemes over these varying parameters can be executed in short times. (2) By weighting the auto assigned SINR targets with scalar variables, the feasibility search problem can be reduced to a simpler linear search problem as demonstrated in Section VII-F. The cross-analyses are accurate since no approximations are needed in contrast to the approximately derived feasibility conditions [26].

Relay beamforming design in the existence of direct links has been a long standing open problem due to the challenge in the expression of SINR metric in terms of relay filters as detailed in Section VI. In this paper, an SINR reformulation is proposed that gives good approximation when the direct channels and the effective channels between the transmitters and receivers are independent. The proposed distributed joint transmit and relay beamforming assures improved total power saving than that of only the distributed transmit beamforming optimization. However, since each relay serves all streams in the network, the complexity substantially increases.

The main contributions of this paper are summarized as follows:

• •

A generic network model with multiple multi-antenna nodes at all sides, i.e., the transmitter, relay, and receiver sides, to carry out multi-stream transmissions in the presence of direct links is transformed into a compact matrix system model.

• •

The proposed distributed ADMM algorithm achieves the optimal centralized solutions in the given generic network model. In addition, in terms of convergence rate, computational complexity, and message exchange load performance metrics, the proposed solution surpasses other optimal distributed algorithms.

• •

By eliminating the feasibility of problem stage automatically, an extensive evaluation of the effects of crucial network parameters on the system performance metrics are attained that reveals new insights in this paper.

• •

Analytical and numerical results demonstrate the lower computational complexity and message exchange load, higher convergence rate, i.e., lesser number of iterations, and finally the convergence of proposed distributed algorithm.

• •

SINR approximation at the relay side is proposed to implement distributed joint transmit and relay beamforming optimization that further improves the total power saving at the cost of increased complexity.

The closest to our contribution in this paper is given in [11], where ADMM is also applied. The major difference between [11] and this paper is that a scalar power variable in a single relay network and a beamforming vector in a multi-relay aided network are optimized, respectively. Clearly, the extension in this paper is nontrivial. Further differences between [11] and this paper are discussed in later sections.

ADMM framework is applicable to many problems under mild conditions. The methodology is different than the distributed algorithms that are designed for particular problems. In [27], a distributed algorithm is proposed for interference alignment in signal space. To achieve the alignment in a distributed manner, each user minimizes the interference covariance matrix. Such distributed solutions that are designed for particular areas are in contrary to the distributed ADMM solutions that have broad application areas. In fact, ADMM solution for interference alignment is already exploited in [28].

The rest of the paper is organized as follows. The multi-stream transmission capable multi-relay interference network is introduced in Section II. In Section III, the transmit beamforming design problem for power minimization under stream SINR constraints is formulated. The distributed ADMM solution and benchmark distributed solutions are presented in Section IV. The attributes of distributed solutions including convergence, computational complexity, and message exchange load are studied in Section V. Distributed joint transmit and relay beamforming filter optimization via SINR approximation at the relay side is provided in Section VI. The numerical results and discussions are presented in Section VII, and finally, the paper is concluded with the summary of main results in Section VIII.

II System Model

Consider a $K$-user two-hop MIMO multi-relay interference network aided by $R$ relays. The $i^{th}$ source (transmitter) and the $i^{th}$ destination (receiver) each has $M_{i}$ antennas while the $r^{th}$ relay has $N_{r}$ antennas as shown in Fig. 1. Each transmitter communicates with its corresponding receiver with the aid of all relays. Without loss of generality, $i^{th}$ transmitter and receiver pair can be called $i^{th}$ user. We assume that all relay nodes work in half-duplex mode. Thus the communication between the users is completed in two time slots and there are non-negligible direct links between all transmitters and receivers. In the first time slot, the $i^{th}$ transmitter transmits the ${M_{i}\times 1}$ signal vector $\mathbf{x}_{i}=\mathbf{U}_{i}\mathbf{s}_{i}$, where $\mathbf{U}_{i}\in\mathbb{C}^{M_{i}\times d_{i}}$ and $\mathbf{s}_{i}\in\mathbb{C}^{d_{i}}$ are the transmit beamforming matrix and symbol vector with E\big{(}\mathbf{s}_{i}\mathbf{s}_{i}^{H}\big{)}=\mathbf{I}_{d_{i}} and {E\big{(}\mathbf{s}_{i}\mathbf{s}_{j}^{H}\big{)}=\bm{0}} for $j\neq i$, respectively. Here, $d_{i}$ is the number of streams of the $i^{th}$ user, i.e., the number of independent data streams to be transmitted between the $i^{th}$ transmitter and receiver pair. We assume $d_{i}\leq M_{i}$ for sufficient degrees of freedom in signal detection. The transmitted signal from the user has a power constraint

[TABLE]

where $p_{i}^{\text{max}}$ is the maximum power of the $i^{th}$ transmitter. The received signal at the $r^{th}$ relay and $k^{th}$ receiver in the first time slot is given by

[TABLE]

respectively, where $\mathbf{H}^{\prime\prime}_{ri}$ is the channel from the $i^{th}$ transmitter to the $r^{th}$ relay, and $\mathbf{J}_{ki}$ is the channel between the $i^{th}$ transmitter and the $k^{th}$ receiver. $\mathbf{n}_{r}$ and $\mathbf{n}_{k}(1)$ are the complex additive white Gaussian noise (AWGN) at the $r^{th}$ relay and at the $k^{th}$ receiver in the first time slot with zero mean, and with the covariances {E\big{(}\mathbf{n}_{r}\mathbf{n}_{r}^{H}\big{)}=\sigma_{r}^{2}\mathbf{I}_{N_{r}}} and {E\big{(}\mathbf{n}_{k}(1)\mathbf{n}_{k}^{H}(1)\big{)}=\sigma_{k}^{2}(1)\mathbf{I}_{M_{k}}}, respectively. In the second time slot, the received signal $\mathbf{y}_{r}$ is precoded at the $r^{th}$ relay by the ${N_{r}\times N_{r}}$ relay filter $\mathbf{F}_{r}$, ${\mathbf{x}_{r}=\mathbf{F}_{r}\mathbf{y}_{r}}$. The $r^{th}$ relay transmit power is {p_{r}=E(||\mathbf{x}_{r}||^{2})=\text{tr}\big{(}\mathbf{F}_{r}\mathbf{R}_{y_{r}}\mathbf{F}_{r}^{H}\big{)}\leq p_{r}^{\text{max}}}, where {\mathbf{R}_{y_{r}}=E\big{(}\mathbf{y}_{r}\mathbf{y}_{r}^{H}\big{)}} and $p_{r}^{\text{max}}$ is the covariance matrix of the received signal at and the maximum power of the $r^{th}$ relay, respectively. The relay transmit power can be explicitly written as

[TABLE]

The received signal at the $k^{th}$ receiver in the second time slot is given by

[TABLE]

where $\mathbf{G}_{kr}$ is the channel between the $r^{th}$ relay and the $k^{th}$ receiver.

Define the effective channel from transmitter $i$ to receiver $k$ through all relays as

[TABLE]

Then, $\mathbf{y}_{k}(2)$ in (II) is rewritten as

[TABLE]

To obtain the SINR expression, the received signal at a receiver is written in terms of the desired signal, interference signal, and noise summands. Thus, define the aggregate channel matrix and noise vector as

[TABLE]

The list of channel notations used in the paper is given in Table I for convenience. Then, the aggregate received signal at receiver $k$ can be written as

[TABLE]

For the sake of linear decoding complexity, the intra- and inter-user stream interferences are treated as noise in this paper. The receive filter is $\bar{\mathbf{V}}_{k}\in\mathbb{C}^{2M_{k}\times d_{k}}$ since the receive filters for the time slots $1$ and $2$, $\mathbf{V}_{k}(1)\in\mathbb{C}^{M_{k}}$ and $\mathbf{V}_{k}(2)\in\mathbb{C}^{M_{k}}$, respectively, are stacked in this matrix, i.e., \bar{\mathbf{V}}_{k}\triangleq\left[\begin{array}[]{c}\mathbf{V}_{k}(1)\\ \mathbf{V}_{k}(2)\\ \end{array}\right]. After applying the receive filter $\bar{\mathbf{V}}_{k}$ to the received signal $\mathbf{y}_{k}$, the SINR of the $l^{th}$ stream of the $k^{th}$ user is obtained as

[TABLE]

where \tilde{\zeta}_{k,l}^{i,n}\triangleq\big{|}\bar{\mathbf{v}}_{k,l}^{H}\mathbf{H}_{ki}\mathbf{u}_{i,n}\big{|}^{2} and $\bar{\mathbf{v}}_{k,l}\triangleq[\mathbf{v}_{k,l}^{T}(1)\,\mathbf{v}_{k,l}^{T}(2)]^{T}$ is the $2M_{k}\times 1$ receive beamforming vector for the $l^{th}$ stream of the $k^{th}$ user, i.e., $l^{th}$ column of the receiver filter $\bar{\mathbf{V}}_{k}$ (similarly $\mathbf{v}_{k,l}(i)$, ${i=1,2}$ are the $l^{th}$ columns of the receiver filters $\mathbf{V}_{k}(i)$, ${i=1,2}$, respectively), $\sigma_{n_{k,l}}^{2}=\bar{\mathbf{v}}_{k,l}^{H}\mathbf{R}_{n_{k}}\bar{\mathbf{v}}_{k,l}$ is the power of the aggregate noise after the receive beamforming

[TABLE]

is the covariance matrix of the aggregate noise, and the notation $(j,m)\neq(k,l)$ denotes that $j\neq k$ and/or $m\neq l$.

III Problem Formulation

In this paper, the total power of transmitters and relays under SINR per stream constraints is minimized via the distributed ADMM algorithm [6] by optimizing all stream beamforming filters $\mathbf{u}_{k,l},\,\forall k,l$ in parallel at each stream, aka processor, in the network. From (1) and (3), the total power, i.e., the sum of transmitter and relay powers, can be rewritten in terms of stream filters $\mathbf{u}_{k,l}$ ($l^{th}$ column of the transmitter filter $\mathbf{U}_{k}$) as

[TABLE]

where ${\mathbf{R}_{.k}\triangleq\mathbf{I}_{M_{k}}+\mathbf{R}_{.k}^{\prime},\mathbf{R}_{.k}^{\prime}\triangleq\sum_{r=1}^{R}\mathbf{R}_{rk}\text{, and}}$

[TABLE]

Similarly, SINR of a stream can be rewritten as

[TABLE]

where $\mathbf{R}_{k,l}^{i}\triangleq\mathbf{H}_{ki}^{H}\bar{\mathbf{v}}_{k,l}\bar{\mathbf{v}}_{k,l}^{H}\mathbf{H}_{ki}$. The superscript index in $\mathbf{R}_{k,l}^{i}$ indicates the $i^{th}$ transmitter, since $\mathbf{H}_{ki}$ is the aggregate channel between the $i^{th}$ transmitter and $k^{th}$ receiver. The simplifications of $\mathbf{R}$ matrices are listed in Table II for convenience.

The problem formulation is given as

P1

[TABLE]

where ${\mathcal{K}\triangleq\{1,2,\ldots,K\}}$ and ${\mathcal{L}_{k}\triangleq\{1,2,\ldots,d_{k}\}}$ is the set of users and set of streams of the $k^{th}$ user, respectively.

IV Multi-Stream Beamforming Under Stream SINR Constraints

PIII is a non-convex quadratically constrained quadratic programming (QCQP) problem [21, 22]. To obtain the $\mathbf{u}_{k,l}$ filters, PIII can be equivalently rewritten by using the fact that $\mathbf{u}_{k,l}^{H}\mathbf{R}_{k,l}^{i}\mathbf{u}_{k,l}=\text{tr}(\mathbf{X}_{k,l}\mathbf{R}_{k,l}^{i})$, where $\mathbf{X}_{k,l}\triangleq\mathbf{u}_{k,l}\mathbf{u}_{k,l}^{H}$, and by rewriting the SINR constraint as a summation inequality rather than a division inequality. PIII can be rewritten as

P2

[TABLE]

The transmit beamforming covariance matrix constraint (13c) imposes the convex constraint that $\mathbf{X}_{k,l}$ matrix belongs to the cone of symmetric and positive semi-definite matrices of dimension $M_{k}$ (denoted by $\mathbb{S}_{+}^{M_{k}}$). Note that, since the covariance matrices in (13b), i.e., $\mathbf{R}_{k,l}^{k}$ and $\mathbf{R}_{k,l}^{j}$, are Hermitian, $\text{tr}(\mathbf{X}_{k,l}\mathbf{R}_{k,l}^{i})\in\mathbb{R}$, which means that the SINR inequality constraint (13b) is well defined. However, PII is still non-convex due to the last constraint (13d), hence SDR can be applied, i.e., the last constraint can be relaxed [21, 22]. Nonetheless, the resulting SDR of PII still cannot be solved distributively. To obtain $\mathbf{X}_{i,n}$ at the stream ${(i,n)}$ via a parallel and distributed approach, both the objective function and constraints in PII must be separable with respect to each stream. The objective function of PII is separable; however, the reformulated SINR constraints (13b) are coupled. Moreover, since the SINR constraints are not linear, ADMM is not directly applicable to PII.

IV-A Proposed ADMM Algorithm

Before presenting the proposed ADMM algorithm, we briefly review the main steps of ADMM. ADMM can solve the following convex problem

P($\mathbf{x},\mathbf{y}$)

[TABLE]

where $\mathbf{A}\in\mathbb{C}^{k\times m}$, $\mathbf{B}\in\mathbb{C}^{k\times n}$, $\mathbf{c}\in\mathbb{C}^{k}$, the functions $f$ and $g$ are convex, $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are nonempty convex sets. Then the augmented Lagrangian for P($\mathbf{x},\mathbf{y}$) is given as

[TABLE]

where $\mathbf{z}\in\mathbb{C}^{k}$ is the Lagrange multiplier of the constraint (14b), re(.) is the real part operator, and $\rho$ is again the Lagrangian dual update step size. P($\mathbf{x},\mathbf{y}$) is solved by ADMM via three steps at each iteration $s$ as follows

[TABLE]

In order to transform the reformulated SINR constraints to a linear form in PII, we initially introduce an auxiliary variable for each summand in the SINR constraint (13b) as

[TABLE]

Then, since the SINR constraints are active at the optimal point, i.e., the SINR constraints (13b) must hold with equality [29], the resulting SDR of PII can be equivalently rewritten as

P3

[TABLE]

Now, the SINR constraint (13b) is in a linear form via the constraints (18b), (18c), and (18d), and the coupling constraint (18b) is a simple linear constraint that is viable for the ADMM algorithm. The partial augmented Lagrangian for PIV-A can be written as

[TABLE]

where $\lambda_{k,l}$ is the Lagrange multiplier of the constraint (18b), and $\rho\in\mathbb{R}_{+}$ is a positive constant parameter for adjusting the convergence speed, i.e., Lagrangian dual update step size.

P($\mathbf{x},\mathbf{y}$) and PIV-A have the following correspondences

[TABLE]

where $x_{kl}\triangleq f_{x}\left(\text{tr}(\mathbf{X}_{k,l}\mathbf{R}_{.k})\right)$ and $f_{x}$ is a mapping function between $x_{kl}$ and $\mathbf{X}_{k,l}$.

As mentioned earlier, in [11], total power minimization under SINR constraints problem is solved via a distributed power control algorithm in a simpler network. In contrast, in this work, the problem is solved via beamforming vectors in a generic network. Hence, the problem in this work is more challenging and also the results are more effective than [11], i.e., higher sum-SINRs can be achieved with lesser power consumptions. Another distinction between [11] and this paper is the difference between the augmented Lagrangian functions, where both of the proposed solutions are fundamentally based on. As seen in [11, P2], the auxiliary definitions [11, Eqs. (3c), (3d)] are augmented in the Lagrangian function [11, Eq. (4)] to obtain closed-form solutions, whereas in this work, the summation of auxiliary terms (18b) are augmented as seen in (IV-A). In [11], power control algorithm is proposed for a simpler relay network architecture. In other words, transmit power scalar variables $p_{k,l}$ are optimized in [11] as opposed to the transmit beamforming filters $\mathbf{u}_{k,l}$ in this work. The summation of auxiliary terms in [11, Eq. (3b)] cannot be augmented in the Lagrangian function. Otherwise, the $p_{k,l}$ variable disappears in the Lagrangian differentiation process, hence closed-form solutions cannot be obtained. On the other hand, the auxiliary definitions (18c) and (18d) cannot be augmented due to the trace operator in this work. Therefore, in this work, (18b) is augmented in the Lagrangian function and the covariance matrix solutions are obtained via CVX. Resorting to CVX for the solutions of covariance matrices is a common practice in the literature [12, 9]. After obtaining the transmit covariance matrices $\mathbf{X}_{k,l}$, obtaining the transmit vectors $\mathbf{u}_{k,l}$ is a well-known process [21], i.e., if the obtained covariance matrix is not rank-one, then additional rank-one approximate solution methods can be applied including the Gaussian randomization method. However, as observed by the numerical results in Section VII, the optimal $\mathbf{X}_{k,l}$ matrices are always rank-one. This indicates that the proposed distributed optimal solution to PIV-A serves as a global optimal solution to PIII.

Remark 1*.*

PIV-A is separable with respect to the streams, thus it requires solving $B$ problems of P{\big{(}\mathbf{X}_{k,l},\zeta_{k,l},\zeta_{k,l}^{\backprime}\big{)}}, where ${B\triangleq\sum_{i=1}^{K}d_{i}}$ is the total number of streams in the network. However, in order to apply the two-block ADMM algorithm which is better understood than $x$-block ADMM algorithms in terms convergence [30], where ${x\!>\!2}$, the variables {\big{\{}\mathbf{X}_{k,l},\zeta_{k,l},\zeta_{k,l}^{\backprime}\big{\}}_{\forall k\in\mathcal{K},\,\forall l\in\mathcal{L}_{k}}} can be divided into two groups $\{\mathbf{X}_{k,l}\}$ and {\big{\{}\zeta_{k,l},\zeta_{k,l}^{\backprime}\big{\}}}. Hence, PIV-A is separated into two simpler parts: P($\mathbf{X}_{k,l}$) and P{\big{(}\zeta_{k,l},\zeta_{k,l}^{\backprime}\big{)}}. Therefore, ADMM can distributively solve PIV-A by solving $2B$ simpler subproblems in parallel.

The Main Steps of ADMM

ADMM consists of sequential updates of primal variables ${\mathbf{X}_{k,l}}$, $\zeta_{k,l}$, and $\zeta^{\backprime}_{k,l}$, and the dual variables ${\lambda_{k,l},\mu_{k,l},\text{and }\mu_{k,l}^{\backprime}}$ [1] as presented below.

Update of $\mathbf{X}_{k,l}$

The smaller problem of $\mathbf{X}_{k,l}$ is

P($\mathbf{X}_{k,l}$)

[TABLE]

where the superscript ′ denotes the constraints to be considered without the symbol $\forall$, e.g., (18c)′ denotes the constraint only for the stream $(k,l)$, not for all streams ${\forall k\in\mathcal{K}}$, $\forall l\!\in\!\mathcal{L}_{k}$ as seen in (18c). The above linear problem can be solved locally at each processor $(k,l)$ in parallel.

Updates of $\zeta_{k,l}\text{ and }\zeta_{k,l}^{\backprime}$

The smaller problem of $\zeta_{k,l}\text{ and }\zeta_{k,l}^{\backprime}$ is

P($\zeta_{k,l},\zeta_{k,l}^{\backprime}$)

[TABLE]

Then according to the Karush-Kuhn-Tucker (KKT) condition, the updates of auxiliary variables are given as

[TABLE]

where $\mu_{k,l}^{s}$ and $\mu_{k,l}^{\backprime s}$ are the Lagrangian multipliers of the constraints (18c) and (18d), respectively.

Updates of $\lambda_{k,l},\mu_{k,l},\text{and }\mu_{k,l}^{\backprime}$

The updates of dual variables ${\lambda_{k,l},\mu_{k,l},\text{and }\mu_{k,l}^{\backprime}}$ are given as

[TABLE]

where ${\rho_{c}\in\mathbb{R}_{+}}$ is the conventional Lagrangian dual update step size.

IV-B Pseudocode

The pseudocode of the proposed ADMM algorithm is given in Algorithm 1, where Matlab scripting language is used. As seen in step 8 of the algorithm, when the absolute deviation from the SINR target of each stream is within $\Delta_{k,l}^{\text{max}}$ accuracy, the algorithm is terminated with a success flag.

IV-C Feasibility of Problem

PIV-A can be infeasible for the given SINR targets $\gamma_{k,l}$. Therefore, the feasibility of PIV-A must be assured in the first stage before solving PIV-A in the second stage. There are basically two techniques in the literature: 1) the feasibility conditions are derived [23], and 2) the initial conditions are relaxed, e.g., the SINR targets are tested and the feasible targets are searched [24, 4]. Deriving the feasibility condition of PIV-A is challenging even for a simpler network and a simpler problem [26]. In [26], power minimization under the worst stream SINR condition is studied for a multi-antenna relay network with a single transmitter and a receiver, and the direct links are neglected. Furthermore, in [26], the feasibility condition is derived based on the signal-to-interference ratio (SIR) instead of SINR. After assuming the targets are feasible based on the approximate feasibility condition derived from the SIR metric, the beamforming vectors are solved based on SINR. The approximate feasibility condition is more accurate in the high SNR regime, where noise can be neglected. On the other hand, testing the SINR targets and searching for the feasible SINR targets [24, 4] are as costly as solving problem PIV-A, i.e., both require high iteration numbers and high computational complexities per iteration, which severely impedes the cross-analysis due to the long simulation durations.

In this work, we adopt a new technique. We randomly initialize the transmit and relay beamforming vectors with full powers, i.e., ${p_{k}^{\text{max}}}$, ${\forall k\in\mathcal{K}}$ and $p_{r}^{\text{max}}$, ${\forall r\in\mathcal{R}}$ to determine feasible SINR targets with a high probability. The infeasible cases that make up a small portion of tests along with slow and fluctuant converging cases are filtered through step 11 and the detection window as detailed earlier. Hence, by testing randomly and automatically generated feasible SINR targets for any permutation of the network parameters, the cross-analysis is executed systematically and extensively within short simulation durations. Since cross-analysis at this comprehensive level has not been performed in the literature yet, many new insights are revealed in this paper. For all simulations, the receive filters are also randomly initialized but normalized to unity since including the receiver power as another parameter in the already large set of cross-varying network parameters substantially complicates the cross-analysis in Section VII.

In fact, the transmit and relay power constraints

[TABLE]

respectively, where ${p_{r}^{\prime\text{max}}\triangleq p_{r}^{\text{max}}-\sigma_{r}^{2}\text{tr}(\mathbf{F}_{r}\mathbf{F}_{r}^{H})}$, can be incorporated into P$(\mathbf{X}_{k,l})$. The Matlab script for CVX solution is given in Algorithm 2. However, instead of incorporating the power constraints into CVX as shown in Algorithm 2, checking the power constraints at step 11 of Algorithm 1 has the following advantages. Firstly, the CVX algorithm runs faster without the incorporated power constraints, particularly for networks with high number of relays and relay antennas. Secondly, if the SINR targets are infeasible, the algorithm starts fluctuating, i.e., different SINRs that do not meet SINR targets are achieved over the iterations while the power constraints are still assured by the constraints incorporated in Algorithm 2. Therefore, it is time-consuming to detect whether the fluctuation is due to infeasibility or it is a rare case where the algorithm still converges after the fluctuation. These two cases significantly hinder the simulation durations and executing extensive cross-analysis in Section VII becomes impractical. As seen in Section VII, in the worst case, our proposed distributed algorithm requires $35$ iterations on average. Therefore, instead of plugging the power constraints (IV-C) into the CVX optimization, letting the algorithm converge to the SINR targets and then checking the power constraints by step 11 of Algorithm 1 can swiftly determine the feasibility of the targets.

IV-D Benchmark Distributed Algorithms

IV-D1 ADMM with Bounded Guarantee (ADMM-BG)

The conventional dual decomposition method [31] for a distributed solution of PIV-A is not applicable since the dual of PIV-A can be unbounded [12]. This can be demonstrated by considering to solve the dual problem

P($\mu_{k,l},\mathbf{X}_{k,l}$)

[TABLE]

where

[TABLE]

instead of solving the subproblem P($\mathbf{X}_{k,l}$) in step a) of our proposed algorithm presented earlier via CVX. The inner optimization problem of the above dual problem can be unbounded below given the dual variable $\mu_{k,l}$ so that $-\eta_{k,l}^{\prime}\text{tr}(\mathbf{X}_{k,l}\mathbf{R}_{k,l}^{k})\rightarrow-\infty$, i.e., the solution can go to minus infinity. This problem can be avoided by the extra quadratic penalty term added in the ADMM scheme

[TABLE]

As seen in [12, Eq. (23)], authors enforce the quadratic penalty term also for the first term $\text{tr}(\mathbf{X}_{k,l}\mathbf{R}_{.k})$ in (IV-D1) by defining an auxiliary variable

[TABLE]

and then introducing the slack variables $t_{k,l}\geq 0,\forall k\in\mathcal{K},\forall l\in\mathcal{L}_{k}$ [12, Eq. (29d)]

[TABLE]

Note that the first term $\text{tr}(\mathbf{X}_{k,l}\mathbf{R}_{.k})$ in (IV-D1) is the objective term in PIV-A, similar to [12, Eq. (17a)]. Hence, the final local partial augmented Lagrangian function is given as

[TABLE]

Following the same principle from [12] as explained above, the modified local partial augmented Lagrangian function of (IV-A) is

[TABLE]

where $\lambda_{k,l}$ is the notation used in (IV-A) that corresponds to $\eta_{k,l}^{\prime}$ used in this subsection. Hence step a) of our proposed distributed algorithm is modified as follows

P($p_{k,l},\mathbf{X}_{k,l}$)

[TABLE]

which can be again solved via CVX. Accordingly, the algorithm 1 needs small modifications. In summary, in ADMM-BG, P($p_{k,l},\mathbf{X}_{k,l}$) is solved instead of P($\mathbf{X}_{k,l}$) in step a) of our proposed algorithm to guarantee bounded solutions.

Due to the added superfluous variables and constraints in (28) and (29), the modified approach, ADMM-BG algorithm, has a slower convergence rate, i.e., more iterations are needed, and also has a higher computational complexity per iteration than the original approach in step a), i.e., the proposed algorithm in this paper. As explained in more detail in Section V, as the problem size increases, i.e., more constraints in this case, the complexity of problem solution increases since CVX depends on the number of constraints. The iteration number also increases as a direct consequence of a larger problem size. Moreover, it is observed that ADMM-BG scheme is significantly more unstable, i.e., the rare frequencies of slow-, fluctuant-, and non-convergent cases significantly increase. Finally, more constraints directly necessitates more variables, i.e., more antennas at the transmitter side. Due to the constraint (28), the minimum number of transmitter antennas is more for ADMM-BG.

IV-D2 Accelerated Distributed Augmented Lagrangians (ADAL)

In [9], an improved version of ADMM is proposed, which is coined as ADAL since it converges faster than the conventional ADMM. ADAL introduces an auxiliary variable for each interference term inside the summation [9, Eq. (11), 3rd constraint] in contrast to the summation of interference terms (18d) as proposed in this work. Hence, instead of two constraints (18c) and (18d) per stream as proposed in this work, ADAL has 1 (desired signal) + ${B-1}$ (interference signals)${=B}$ constraints per stream. More constraints increase 1) the iteration number (the problem size increases, similar to more number of users and antennas), 2) the computational complexity per iteration (the complexity of CVX is dependent on the number of constraints), and 3) the minimum number of antennas (the transmit covariance is involved in more constraints, thus more antennas-variables are needed). Following the approach in [9], the resulting SDR of PII at the optimal point can be equivalently rewritten as

P4

[TABLE]

where ${\bm{\zeta}_{k,l}\triangleq[\zeta_{1,1}^{k,l},\ldots,\zeta_{i,n},\ldots,\zeta_{K,d_{K}}^{k,l}]^{T}}$ is the ${B\times 1}$ vector that stacks all auxiliaries $\zeta_{i,n}^{k,l}$, ${\forall i\in\mathcal{K},\forall n\in\mathcal{L}_{i}}$. The ADAL-Direct algorithm proposed in [9] is shortly referred as ADAL in our work.

The partial augmented Lagrangian for PIV-A can be written as

[TABLE]

where $\bm{\lambda}\triangleq[\lambda_{1,1},\ldots,\lambda_{K,d_{K}}]^{T}$ is the vector of Lagrange multipliers of the constraint (18b), $\rho\in\mathbb{R}_{+}$ is the Lagrangian dual update step size, and $\|.\|_{2}$ is the $L_{2}$ vector norm operator. ADAL proposes to distribute the problem over streams as follows, i.e., the local partial augmented Lagrangian is given as

[TABLE]

where ${\vec{\bm{\zeta}}_{j,m}}$ is the auxiliary variable transmitted from stream ${(j,m),\forall j\neq k,\,\forall m\neq l}$ to all streams ${(i,n)}$, ${\forall i\in\mathcal{K}\setminus\{j\}}$, ${\forall n\in\mathcal{L}_{i}\setminus\{m\}}$. In other words, the ${\vec{\bm{\zeta}}_{j,m}}$ variable is evaluated at stream ${(j,m)}$ and is transmitted to all streams ${(i,n),\forall i\neq j,\,\forall n\neq m}$, thus it is a constant vector for stream ${(k,l)}$.

The Main Steps of ADAL

Updates of $\bm{\zeta}_{k,l}$ and $\mathbf{X}_{k,l}$

The smaller problem of $\bm{\zeta}_{k,l}$ and $\mathbf{X}_{k,l}$ is

P($\bm{\zeta}_{k,l},\mathbf{X}_{k,l}$)

[TABLE]

where (33d)′ this time denotes the constraint (33d) without $\forall j,m$; that is, for a particular stream $(k,l)$. For instance, assume $K=3$ and $d=2$, then for stream $(k,l)=(1,2)$, the set of $\zeta$ variables is $Z=\{\zeta_{1,1}^{1,2},\zeta_{1,2}^{1,2},\ldots,\zeta_{3,1}^{1,2},\zeta_{3,2}^{1,2}\}$, where ${|Z|=B}$. The above linear problem can be solved locally at each processor in parallel via CVX in Matlab.

Update of $\vec{\bm{\zeta}}_{k,l}$

The $\bm{\zeta}_{k,l}^{s}$ output of P($\bm{\zeta}_{k,l},\mathbf{X}_{k,l}$) is further updated as follows

[TABLE]

where $\tau$ is a step size. $\vec{\bm{\zeta}}_{k,l}^{s+1}$ is to be transmitted to the other nodes as mentioned earlier.

Update of $\bm{\lambda}$

As the final step, the dual variables are updated as follows

[TABLE]

Note that the above update is also local to a processing node $(k,l)$ since $\vec{\bm{\zeta}}_{j,m}^{s+1},\forall(j,m)\neq(k,l)$ updates are provided from other nodes than $(k,l)$.

The ${B\!-\!1}$ constraints in (33d)′ of P($\bm{\zeta}_{k,l},\mathbf{X}_{k,l}$) are overwhelming as opposed to our proposed ADMM algorithm, where the summation of these terms is utilized. Hence, the disadvantages of ADAL follow as mentioned earlier.

Remark 2*.*

$\bm{\zeta}_{k,l}$ and $\mathbf{X}_{k,l}$ variables in P($\bm{\zeta}_{k,l},\mathbf{X}_{k,l}$) are jointly solved via CVX. the closed-form solution of $\bm{\zeta}_{k,l}$ can be obtained from the KKT conditions, then only the $\mathbf{X}_{k,l}$ variables are solved via CVX to reduce the overall computational complexity.

Consider whether the joint solution of $\bm{\zeta}_{k,l}$ and $\mathbf{X}_{k,l}$ variables via CVX or closed-form solution of $\bm{\zeta}_{k,l}$ and CVX solution of $\mathbf{X}_{k,l}$ option should be chosen. As mentioned earlier, there are $B$ equality constraints per stream due to the (33c)′ and (33d)′ constraints. With $\zeta_{k,l}^{i,n}$ (joint solution) each constraint is a bivariate equation of $\zeta_{k,l}^{i,n}$ (scalar) and $\mathbf{X}_{k,l}$ (matrix), whereas without $\zeta_{k,l}^{i,n}$ (no joint solution), each equality is a univariate equation of $\mathbf{X}_{k,l}$. Thus, in the univariate case, $B$ is the minimum number of antennas required at a node, and in the bivariate case, the requirement is lower than $B$ by the help of additional variables $\zeta_{k,l}^{i,n}$. As shown earlier, our proposed method has significantly lower number of equality constraints, thus both the iteration number (can be determined from numerical results since the computation of iteration number is not possible) and computational complexity (can be approximately computed in terms of limiting behavior) per iteration are significantly low. Therefore, we choose joint solution of $\bm{\zeta}_{k,l}$ and $\mathbf{X}_{k,l}$ via CVX for ADAL so that the iteration number comparison with our proposed distributed algorithm and the number of antennas in Section VII is more fair and convenient to compare.

From the above discussion, it can be also concluded that centralized algorithms can be less restrictive than the distributed algorithms in terms of minimum number of antenna requirements in general. The SDR of the QCQP problem PII can be solved via CVX [32] in a centralized manner. Hence, for the centralized algorithm, there are $Kd$ constraints and $KM$ variables. On the other hand, for ADAL, there are ${Kd+1\approx Kd}$ constraints and $M$ variables per processor. Therefore, the number of variables for the distributed algorithm is approximately shrunk $1/K$ times compared to the centralized solution. However, a centralized solution is less attractive when the dimension of the problem is large as mentioned earlier.

Similarly, it can be concluded that the power control algorithms are also more restrictive than the beamforming optimization problems. In particular, there is only 1 power variable $p_{k,l}$ per stream, on the other hand, there are $M_{k}$ variables per stream in the transmit beamforming vector $\mathbf{u}_{k,l}\in\mathbb{C}^{M_{k}}$. Since the number of variables in the beamforming vectors that directly determine the SINR targets is significantly larger than the number of power variables, automatic generation of feasible SINR targets by random initializations of beamforming vectors is not possible for power control algorithms.

IV-E Competitive SINR Targets

PIII, where the objective function is minimizing the total power and the QoS guarantee is meeting the SINR targets, is one of the most important optimization problems for resource allocation in wireless networks [12, 14, 9, 7, 11]. This problem can be considered as the dual interest of network utility maximization problems [33], [18, Chap. 8]. Depending on the system design goals, the trade-off between the consumed total power and the achieved sum-SINR in the network can be controlled by the given SINR targets.

The SINR targets can be determined at the physical or higher layers of the network. One of the approaches in the physical layer is to search for feasible SINR targets between the lower and upper bounds set by beamforming filters. For the lower bound, random filter initializations are good choices as supported by the numerical results in this paper. For the upper bound, SINR maximizing (max-SINR) filter initializations can be good choices as shown in [27]. In Section VII-F, we present the numerical results showing that by searching for higher feasible SINR targets between the lower and upper bounds, the system can achieve higher SINRs at the cost of increased power consumptions. In Appendix A, closed-form solutions of max-SINR filters for distributed implementation are provided.

V Attributes of the Algorithms

V-A Convergence

The global convergence of two-block ADMM is proven based on mild conditions in [30]. The optimal $\mathbf{X}_{k,l}$ solution for each subproblem of PIV-A is a global optimal solution for the complete SDR problem in PIV-A since at each iteration, each subproblem of PIV-A is convex, and the iterations are pursued until convergence. We analytically and numerically show in this section and in Section VII, respectively, that the proposed distributed algorithm is convergent.

The convergence of ADMM to the global optimum of P($\mathbf{x},\mathbf{y}$) under mild conditions is given by the following lemma.

Lemma 1*.*

[1, Proposition 4.2] Assuming that the optimal solution set of P($\mathbf{x},\mathbf{y}$) is nonempty, and $\mathbf{A}^{T}\mathbf{A}$ and $\mathbf{B}^{T}\mathbf{B}$ are invertible, the sequence of solutions $\{\mathbf{x}^{s+1},\mathbf{y}^{s+1},\mathbf{z}^{s+1}\}$ obtained from (16a), (16b), and (16c), respectively, is bounded and every limit point of $\{\mathbf{x}^{s+1},\mathbf{y}^{s+1}\}$ is an optimal solution of P($\mathbf{x},\mathbf{y}$).

Proposition 1*.*

The iterates $\{\mathbf{X}_{k,l}^{s+1},\zeta_{k,l}^{s+1},\zeta_{k,l}^{\backprime\,s+1}\}_{\forall k,l}$ and $\{\lambda_{k,l}^{s+1},\mu_{k,l}^{s+1},\mu_{k,l}^{\backprime s+1}\}_{\forall k,l}$ in Algorithm 1 converge to the optimal solution of primal and dual problems as $s\rightarrow\infty$. When the algorithm converges, the optimal $\{\mathbf{X}_{k,l}\}_{\forall k,l}$ solutions obtained in step 5 of Algorithm 1 are the global optimal solutions of PIV-A.

Proof.

Since ${g(\mathbf{y})=0}$ and ${\mathbf{A}=\bm{0}}$, the Lagrangian function (IV-A) is modified correspondingly, and also, the ${\mathbf{A}^{T}\mathbf{A}}$ condition in Lemma 1 is dropped. Based on Lemma 1 and the correspondences in (IV-A), the proposition on the convergence of the proposed distributed algorithm follows. ∎

V-B Computational Complexity

The most computationally intensive part of the proposed ADMM algorithm is step 5 of Algorithm 1, where the general-rank positive semi-definite matrices $\{\mathbf{X}_{k,l}\}$ are obtained explicitly via CVX without closed-form solutions. Obtaining closed-form beamforming solutions is challenging even in simpler network architectures [12, 9]. After obtaining the optimal general-rank positive semi-definite matrices $\{\mathbf{X}_{k,l}\}$ via the proposed ADMM algorithm, rank-one solutions $\mathbf{u}_{k,l}$, where $\mathbf{u}_{k,l}\mathbf{u}_{k,l}^{H}=\mathbf{X}_{k,l}$, can be efficiently derived [21]. Although there are no closed-form solutions, semi-definite programming (SDP) has still many real-time applications [22, 34].

Remark 3*.*

Similar to the subproblem of proposed distributed algorithm P($\mathbf{X}_{k,l}$) in (IV-A), note that the SDP solutions, e.g., via CVX, are also needed for the subproblems of benchmark schemes ADMM-BG P($p_{k,l},\mathbf{X}_{k,l}$) in (IV-D1) and ADAL P($\bm{\zeta}_{k,l},\mathbf{X}_{k,l}$) in (IV-D).

The complexity of SDP solution for a generic SDR of a QCQP problem [21, Eq. (5)] is given as $\mathcal{O}(\max\{M,c\}^{4}M^{1/2}\log(1/\epsilon))$, where $M$ is the antenna number, $c$ is the number of constraints, and ${\epsilon>0}$ is the solution accuracy. Note that $c=1,2,\text{and }B$ for the proposed algorithm, ADMM-BG, and ADAL, respectively. Hence, for the former two algorithms, ${M\geq c}$, while for the last algorithm, ${M\geq c}$ or ${M<c}$. For the sake of simplicity, we evaluate the complexity of algorithms based on $\mathcal{O}(Mc)$; that is, the limiting behavior of simple multiplication of the number of variables and constraints. The complexity benchmark of algorithms based on this simplified formula also matches with CPU time based simulation benchmarks as shown in the end of next section.

For the proposed distributed algorithm, $P(\mathbf{X}_{k,l})$ has $M$ variables and $1$ constraint, P($\zeta_{k,l},\zeta_{k,l}^{\backprime}$) has $2$ variables that have closed-form solutions (23), and finally $3$ more variables in (24), thus the complexity of proposed distributed algorithm is ${\mathcal{O}(M+5)}$ per processor. For ADMM-BG, P($p_{k,l},\mathbf{X}_{k,l}$) has $M$ (from $\mathbf{X}_{k,l}$) $+$ $1$ (from $p_{k,l}$) variables and $2$ constraints. Due to the similar steps to our proposed algorithm, ADMM-BG also has ${2+3=5}$ variables, and additional $1$ variable $t_{k,l}$. Hence, the complexity of ADMM-BG is ${\mathcal{O}(2(M+1)+6)}$. Finally, for ADAL, P($\bm{\zeta}_{k,l},\mathbf{X}_{k,l}$) has $M+B$ variables and $B$ constraints, and $B$ variables in each of (37) and (38). Thus, the complexity of ADAL is ${\mathcal{O}(B(M+B)+2B)}$. In summary, the proposed algorithm has the lowest complexity, followed by ADMM-BG, and then ADAL. The numerical CPU benchmarks of the proposed and existing algorithms are presented in Section VII-D.

V-C Message Exchange Load

For the update of (24c) at each processor, a single value is needed, i.e., the summation term, which can be delivered by a collector node. This means that each processor needs to send out ${B\!-\!1}$ scalar values, one for each other processor, to the collector node. In total, ${B(B\!-\!1)}$ and $B$ numbers of scalar values need to be exchanged from the processors to the collector node and from the collector node to the processors, respectively. As a result, in total, $B^{2}$ scalar values need to be exchanged in the network. Since the update (24c) is also needed for ADMM-BG, the message exchange loads of ADMM-BG and the proposed algorithm are same. For ADAL, the updates (38) are achieved by sending out (37) from each stream to the collector node. Therefore, for both message exchange directions between the collector node and processors, exchange of $B^{2}$ scalar values is needed, making a sum of $2B^{2}$ exchange of scalar values in the network. In summary, ADAL has the highest message exchange load, followed by ADMM-BG and the proposed algorithm.

Note that a deployed central node or any node among the existing nodes of relay network can serve as a collector node. Moreover, parallel updates of the variables in the ADMM algorithm are robust to delays and errors [6, 28].

VI Distributed Joint Transmit and Relay Beamforming Filter Optimization

Relay filter design has been a long standing open problem when the direct links exist. To pinpoint the hurdle that direct links cause in relay filter design, SINR (8) is reformulated as an explicit function of relay filters next.

VI-A SINR Reformulation

Consider the signal power from stream $n$ of transmitter $i$ to stream $l$ of receiver $k$

[TABLE]

where

[TABLE]

${\mathbf{v}_{k,l}(1)}$ and ${\mathbf{v}_{k,l}(2)}$ is the $M_{k}\times 1$ receive beamforming vector for the $l^{th}$ stream of the $k^{th}$ user in the first and second time slot, respectively.

The first and second summand in the last equality of (43) is independent of relay filters and can be rewritten in terms of Hermitian matrices. ${\mathbf{A}^{H}+\mathbf{A}}$ in the third summand is a Hermitian matrix, thus the third summand can be rewritten as ${\text{re}(2\mathbf{u}_{i,n}^{H}\mathbf{A}\mathbf{u}_{i,n})}$, where

[TABLE]

Hence, the signal power from stream $n$ of transmitter $i$ to stream $l$ of receiver $k$ via relay $r$ is given as

[TABLE]

where

[TABLE]

is a non-Hermitian matrix. Due to the asymmetry that the direct links bring through the third summands of signal powers in (43), relay filter design in the existence of direct link is a non-trivial problem. The third summand is particularly not negligible when the direct channels $\mathbf{J}_{ki}$ and the effective channels $\mathbf{H}^{\prime}_{k.i}$ are correlated. However, when these channels are independent, the affect of this term is less prominent. By omitting this term, the approximate SINR of the $l^{th}$ stream of the $k^{th}$ user is obtained as

[TABLE]

where

[TABLE]

$\sigma_{n_{k,l}}^{2}$ is the sum of noises in the first and second time slots that can be obtained from (9).

Next, (47) can be rewritten as an explicit function of relay filter

[TABLE]

where ${\mathbf{f}_{r}\triangleq\text{vec}(\mathbf{F}_{r})},$ \bar{\mathbf{f}}_{r}\triangleq\left[\begin{array}[]{c}\mathbf{f}_{r}\\ t\\ \end{array}\right], ${|t|^{2}=1}$, and ${\mathbf{Y}_{r}\triangleq\bar{\mathbf{f}}_{r}\bar{\mathbf{f}}_{r}^{H}}$. Further derivation details of (48) and the definition of $\tilde{\mathbf{C}}_{klrin}$ are given in Appendix B.

VI-B Total Power Reformulation

Thus far, the exact SINR (11) is approximated by (47), and the approximate SINR is rewritten as a function of relay filter in homogenous quadratic form (48). Next, the total power is rewritten as a function of relay filters.

Noting the relay power (3) and using the equality $\text{tr}(\mathbf{A}\mathbf{B}\mathbf{A}^{H})=\text{vec}(\mathbf{A})^{H}(\mathbf{B}\otimes\mathbf{I})\text{vec}(\mathbf{A})$, the total power is obtained as

[TABLE]

where

[TABLE]

Furthermore, (49) can be rewritten as

[TABLE]

where

[TABLE]

VI-C Problem Formulation

Using (48) and (50), the relay beamforming filter design problem is given as

P5

[TABLE]

where \bm{\Theta}\triangleq\left[\begin{array}[]{cc}\bm{0}&\bm{0}\\ \bm{0}&1\\ \end{array}\right]. The partial augmented Lagrangian is given as

[TABLE]

where $\dot{\lambda}_{r}^{k,l}\in\mathbb{R}$ is the Lagrange multiplier of the constraint (51b), and $\dot{\rho}\in\mathbb{R}_{+}$ is the Lagrangian dual update step size. From hereafter, the distributed solution follows the similar steps of the distributed transmit beamforming filter design proposed in Section IV-A.

VI-D Pseudocode

The algorithm for the distributed joint solution is obtained by adding two more steps after step 5 of Algorithm 1. Similar to the transmitter side optimization, the first step after step 5 of Algorithm 1 is obtaining $\{\mathbf{Y}_{r}\}$. The second step is obtaining the auxiliary variables $\{\dot{\zeta}_{r}^{k,l},\dot{\zeta}_{r}^{k,l\,\backprime}\}$, and the Lagrangian multipliers $\{\dot{\mu}_{k,l}\}$ and $\{\dot{\mu}_{k,l}^{k,l\,\backprime}\}$ for the constraints (51c) and (51d), respectively.

The discussions of distributed transmit beamforming filter design hold for the distributed relay beamforming filter design as well. For instance, CVX optimization and closed-form solutions are utilized for $\{\mathbf{Y}_{r}\}$ and $\{\dot{\zeta}_{r}^{k,l},\dot{\zeta}_{r}^{k,l\,\backprime},\dot{\mu}_{k,l},\dot{\mu}_{k,l}^{k,l\,\backprime}\}$, respectively, and the optimal $\mathbf{Y}_{r}$ matrices are observed to be always rank-one. However, in contrast to the transmit side optimization as discussed at the end of Section IV-C, the relay power constraint (51e) is better to be incorporated into the CVX optimization, and the transmit power constraint (IV-C) at the transmitter side as well, since the distributed joint algorithm is significantly more costly in computations due to the distributed relay filter optimization as shown in Section VII-G.

VII Numerical Results

In this section, the performance of the proposed scheme is compared with the centralized, ADMM-BG, and ADAL solutions by cross-varying the network parameters. All distributed algorithms in this section are optimal, thus they achieve the same SINR targets and also have the same total power consumptions that are achieved by the centralized solution. The network parameters to be cross-varied are the number of transmit ($M$) and relay ($N$) antennas, the number of relays ($R$), and the transmit ($\text{SNR}_{\text{T}}$) and relay ($\text{SNR}_{\text{R}}$) SNRs.

VII-A Simulation Settings

In our simulations, the distances between transmitter and relay sides, and relay and receiver sides are same and fixed to 1 km. The shadowing coefficient follows log-normal distribution with zero mean and standard deviation equal to $8$. The power of receive filter $\bar{\mathbf{v}}_{k,l}$ is normalized to 1 by equal power allocation over the receive filters $\mathbf{v}_{k,l}(i),\,i=1,2$, i.e., the receive filter power for each time slot is normalized to $\sqrt{0.5}$. The transmitters and relays have the same transmit power constraints among themselves, i.e., $p_{k}^{\text{max}}=p_{\text{T}}^{\text{max}}$, $\forall k\in\mathcal{K}$ and $p_{r}^{\text{max}}=p_{\text{R}}^{\text{max}}$, $\forall r\in\mathcal{R}$, where ${\mathcal{R}\triangleq\{1,2,\ldots,R\}}$ is the set of relays, but $p_{\text{T}}^{\text{max}}$ and $p_{\text{R}}^{\text{max}}$ are not necessarily equal. Hence, the SNR at the transmitter and relay side can be defined as

[TABLE]

respectively, and $\sigma_{k}^{2}=\sigma^{2},\,\forall k\in\mathcal{K}$ is assumed for simplicity. Without loss of generality, the streams of a user are assumed to have the same SINR target, which is equal to the average SINR of the user, i.e., the average of user’s stream SINRs. For both the proposed algorithm and ADMM-BG, the step sizes $\rho$ and $\rho_{c}$ are empirically tuned to $1.2$ and $0.5$, respectively. For ADAL, the step sizes $\rho$, $\rho_{c}$, and $\tau$ are tuned to $9$, $0.5$, and $0.3$, respectively. The step sizes are kept constant for the simulations. If $|\text{SINR}_{k,l}-\gamma_{k,l}|\leq\Delta_{k,l}^{\text{max}}=10^{-4}$ is met by the stream, the stream is regarded as achieving the SINR target.

Finally, for the benchmarks of the proposed algorithm and centralized solution, $100$ Monte Carlo channels are tested. Due to the slow convergence rates of ADMM-BG and ADAL, $20$ Monte Carlo channels are tested for their benchmarks with the proposed algorithm.

VII-B Some Notes

For all simulations, random initializations of filters and channels are pre-saved in files to feed the same inputs to the algorithms. Thus, reliable conclusions with less number of channel tests can be reached. Especially for large network sizes, this approach significantly reduces the simulation durations. Therefore, the networks with many varying parameters are tested in reasonable time durations by the use of pre-saved inputs.

The numerical results in this section are presented for ${K\!=\!3}$ networks since cross-varying $K$ as well significantly increases the data sets to be analyzed. In all tests, each user is assumed to have $d_{k}=2$ streams. Thus, a total of ${B=6}$ streams are transmitted in all network configurations. Inline with the objective function of PIII, the total power curves are plotted in this section. On the other hand, although QoS assurance of PIII imposes an SINR target for each of the streams, a single curve, sum-SINR, is plotted instead of plotting $6$ individual stream SINR curves for the sake of clarity of the figures.

For all simulations of the proposed distributed algorithm presented in this section excluding the Section VII-F, the percentage of all cases including slow-, fluctuant-, and non-convergent, i.e., infeasible SINR targets, makes up in total $1.93\%$. We note that in Section VI-F, we propose the linear search method to determine higher feasible SINR targets than those are set by random filter initializations. The next iteration of linear search method is proceeded depending on whether the result of the previous iteration is feasible or infeasible.

The conclusions drawn in this section by cross-varying these parameters are not straightforward, e.g., increasing $R$ can assist in achieving the SINR targets while also increasing the power consumption. Moreover, as discussed later, the ratios of these parameters also significantly affect the outputs, e.g., the ratios of $\text{SNR}_{\text{T}}$ and $\text{SNR}_{\text{R}}$ have different effects on the outputs.

VII-C Proposed ADMM vs. Centralized Solution (Figs. 2 – 4)

VII-C1 Performances over Network Sizes and SNRs (Fig. 2)

In Figs. 2(a), 2(b), and 2(c), the numerical results of average total power, sum-SINR, and iteration number vs. $M$ and $R$ are plotted, respectively.

As seen in Fig. 2(a), the proposed algorithm achieves the optimal centralized algorithm solutions, i.e., they have the same power consumptions. As mentioned earlier, all algorithms achieve the given SINR targets as seen in Fig. 2(b), and the optimality of an algorithm is determined whether it has the same power consumption with the centralized solution or not. In all cases, the numerical sum-SINR results in Section VII equal to the sum of the given feasible SINR targets of each stream, which are determined by random initializations of beamforming vectors.

As seen in Figs. 2(a) and 2(b), the network $\text{K}3\text{-}\text{M}4\text{-}\text{N}8\text{-}\text{R}9$ has a higher sum-SINR and a lower total power consumption than $\text{K}3\text{-}\text{M}3\text{-}\text{N}8\text{-}\text{R}10$, respectively. Since the optimization problem, power minimization while achieving the SINR targets, is solved at the transmitter side, providing more resources, i.e., more antennas, to the transmitters returns better results, i.e., higher SINRs are achieved with lower power consumptions. As expected, the sum-SINRs of the networks $\text{K}3\text{-}\text{M}10\text{-}\text{N}8\text{-}\text{R}3$, $\text{K}3\text{-}\text{M}10\text{-}\text{N}8\text{-}\text{R}4$, and $\text{K}3\text{-}\text{M}15\text{-}\text{N}8\text{-}\text{R}10$ gradually increase as the total numbers of relays in the networks increase. However, their total power consumptions remain similar because although $R$ is increased, $M$ is chosen sufficiently large for these networks.

As seen in Fig. 2(c), the iteration number increases when the SNR increases. When the transmitter and relay SNRs are $12$ (${\text{SNR}_{\text{T}}\!=\text{SNR}_{\text{R}}\!=\!12}$), the iteration numbers are similar around 20$\pm$1. Thus we focus on the ${\text{SNR}_{\text{T}}\!=\text{SNR}_{\text{R}}\!=\!21}$ results. In our problem setting, where only transmit beamforming vectors are optimized, increasing $M$ and decreasing $R$, decreases the iteration number. Lower $R$ yields lower SINR, in other words, the range of numbers in consideration is lower, thus the iteration number is lower. Although higher $M$ increases SINR, the iteration number decreases because, again, providing more resources, i.e., more antennas, to the transmitters returns better results, i.e., less number of iterations. The algorithm converges faster because the degree of freedom is richer due to more resources, variables, available for the optimization of transmit beamforming vectors.

VII-C2 Performances over a Wider Range of SNRs for a Particular Network Size (Fig. 3)

In Figs. 3(a), 3(b), and 3(c), the numerical results of average total power, sum-SINR and iteration number vs. transmit and relay SNRs for the $\text{K}3\text{-}\text{M}10\text{-}\text{N}8\text{-}\text{R}3$ network are presented, respectively. Again, the proposed algorithm achieves the optimal centralized solutions as seen in Fig. 3(a). In general, increasing SNR increases the power consumption and also sum-SINR. However, consider the blue trend lines in Figs. 3(a) and 3(b). When $\text{SNR}_{\text{T}}$ is sufficiently large, e.g., $\text{SNR}_{\text{T}}\!=\!42$, increasing $\text{SNR}_{\text{R}}$ can have a small effect on sum-SINR as seen in the $42\text{-}12$, $42\text{-}21$, and $42\text{-}42$ dB results of Fig. 3(b). Next, consider the green trend lines in Figs. 3(a) and 3(b). When $\text{SNR}_{\text{R}}$ is sufficiently large, increasing $\text{SNR}_{\text{T}}$ is advantageous as seen in $12\text{-}42$, $21\text{-}42$, and $42\text{-}42$ dB results of Figs. 3(a) and 3(b). $\text{SNR}_{\text{R}}=42$ is sufficiently large to assist in achieving higher SINR targets, and providing more resources to the transmitter side returns better results, i.e., higher SINRs are achieved with lower power consumptions. On the other hand, $\text{SNR}_{\text{R}}=12\text{ and }21$ values are not sufficiently large. Thus, both the total power consumption and sum-SINR are increased by the increasing $\text{SNR}_{\text{T}}$, e.g., consider the $12\text{-}21$, $21\text{-}21$, and $42\text{-}21$ dB results.

In Figs. 3(a) and 3(b), the numerical results for another benchmark that assumes the direct links between transmitters and receivers do not exist are presented. When $\text{SNR}_{\text{T}}$ is high, more power is consumed while less sum-SINR is achieved as seen in the $42\text{-}12$ and $42\text{-}21$ dB results. However, when $\text{SNR}_{\text{R}}$ is high, the effect of nonexisting direct links vanishes as seen in the $21\text{-}42$ and $42\text{-}42$ dB results.

As seen in Fig. 3(c), in general, the iteration number increases as SNR increases, by increasing both $\text{SNR}_{\text{T}}$ and $\text{SNR}_{\text{R}}$, e.g., $12\text{-}12$, $21\text{-}21$, and $42\text{-}42$ dB, and by only increasing $\text{SNR}_{\text{T}}$, e.g., $12\text{-}12$, $21\text{-}12$, and $42\text{-}12$ dB. However, when only $\text{SNR}_{\text{R}}$ is increased and reaches up to $42$ dB, the iteration number decreases, i.e., compare $12\text{-}12$ and $12\text{-}21$ with $12\text{-}42$, and compare $21\text{-}12$ and $21\text{-}21$ with $21\text{-}42$ dB results. As seen in Fig. 3(a), there are notable peaks in power consumptions at the $12\text{-}42$ and $21\text{-}42$ dB points. The disproportionately high $\text{SNR}_{\text{R}}$ value with respect to $\text{SNR}_{\text{T}}$ helps in rapidly achieving the SINR targets, which cannot be controlled by the transmitter side, before the transmit beamforming optimization can further reduce the total power consumption. When the marginal ${\text{SNR}_{\text{T}}=42}$ results are compared, as explained earlier, similar sum-SINR results are achieved for these dB points while total power consumptions are reduced from left to right as displayed in Fig. 3(a). As mentioned earlier, smaller range of numbers of interest, i.e., smaller power consumption values, decreases the iteration number. Thus the iteration number decreases gradually in the last three dB points of Fig. 3(c).

The total power savings can be obtained from the figures, i.e., total power saving (dB)${=(Kp_{\text{T}}^{\text{max}}+Rp_{\text{R}}^{\text{max}})}$ (dB)$-$(y-axis value) (dB). As seen in the results, e.g., Figs. 2(a) and 3(a), the power savings are large as explained next. Due to the random initializations of beamforming vectors, the initial power budgets are likely to be highly redundant to achieve the SINR targets determined by these beamforming vectors. On the other hand, if max-SINR filters are used for the initializations instead of random initializations to determine the SINR targets, the SINR targets are likely to be infeasible.

VII-C3 Supporting Results for Disproportionate SNRs over Network Sizes (Fig. 4)

In Fig. 4, the early convergence of the algorithm when $\text{SNR}_{\text{R}}$ is disproportionately high is shown again over different network sizes. Since the algorithm achieves the SINR targets rapidly due to the high $\text{SNR}_{\text{R}}$, the transmitter side lacks the opportunity to further reduce the power consumption similar to the case observed in Fig. 3.

VII-C4 Scrutinized Performances over Multiple and Single Channel Realizations (Figs. 5 and 6)

In Fig. 5, total power consumptions vs. channel realizations are presented to demonstrate that the proposed algorithm achieves the optimal centralized solutions at each channel realization. In Figs. 6(a) and 6(b), the total power consumptions and total absolute target SINR deviations, $|\Delta|\triangleq|\sum_{k,l}^{K,d_{k}}\Delta_{k,l}|$ vs. iteration numbers are presented for a randomly selected channel realization.

VII-D Proposed ADMM vs. Existing Distributed Algorithms (Table III, Figs. 7 and 8)

The average CPU times for the aforementioned distributed algorithms in the $\text{K}3\text{-}\text{M}6\text{-}\text{N}8\text{-}\text{R}14$ network at ${\text{SNR}_{\text{T}}\!=\!\text{SNR}_{\text{R}}\!=\!21}$ dB are shown in Table III, based on a desktop computer with 64-bit operating system, Intel i7, CPU 3.40 GHz, and 16 GB RAM. A similar comparative trend is exhibited for other network configurations in our extensive experiments. The numerical results in Table III match with the analytical results obtained in Section V-B.

As seen in Fig. 7, the iteration number of the proposed algorithm is significantly lower than ADMM-BG. The iteration numbers of the two algorithms are more distinct when disproportionate resources are allocated at the transmitter side, i.e., higher $M$ and $\text{SNR}_{\text{T}}$. As expected, both algorithms achieve the SINR targets and the optimal centralized solutions, which are not plotted to avoid duplicate results.

In Fig. 8, the iteration numbers of the proposed and ADAL algorithms are compared, respectively. Each network configuration is tested at three different SNRs, particularly, from left to right at $12\text{-}12$, $21\text{-}21$, and $42\text{-}42$ dBs. For the simplicity of figures, SNRs are not noted. As seen in Fig. 8, due to the significantly higher number of constraints in ADAL than the proposed algorithm, the iteration number of our proposed algorithm is always significantly lower than ADAL.

VII-E BER Performance of Proposed ADMM (Fig. 9)

The important application areas of the SINR apparatus are deriving the closed-form BER and outage expressions [35, 36]. Therefore, as mentioned earlier, stream SINR is interconnected with the BER metric. Contrary to the rate results, the BER results cannot be obtained from the numerical SINR results presented in this section. In general, BER and outage performances improve as SINR improves, e.g., SINR outage is the probability of SINR descending below a preset SINR target. To illustrate this general trend, the sum-SINR and the corresponding BER improvement percentages of the system are demonstrated in Fig. 9. For instance, at $\text{SNR}_{\text{R}}=21$ dB, BER is improved, i.e., decreased, by $1.59\%$ by increasing the $\text{SNR}_{\text{T}}$ from $6$ to $12$ dB, respectively. At $\text{SNR}_{\text{R}}=12$ dB, $0.37\%$ BER improvement is observed by increasing the $\text{SNR}_{\text{T}}$ from $3$ to $6$ dB.

The results in Fig. 9 are obtained over $300$ Monte Carlo channels and uncoded QPSK modulation is used. For the earlier numerical results, the extraction of beamforming vectors $\mathbf{u}_{k,l}$ from $\mathbf{X}_{k,l}$ is not needed. However, for BER simulations, $\mathbf{u}_{k,l}$ needs to and can be extracted by rank-one decomposition [21] since the ranks of $\mathbf{X}_{k,l}$ matrices are observed to be always $1$ [37, 38]. Based on this observation, the proposed SDR solution is also the global optimal solution to PIII.

Contrarily, the maximum transmit powers $\text{SNR}_{\text{T}}$ and $\text{SNR}_{\text{R}}$ can be fixed and higher SINRs can still be achieved by searching for competitive SINR targets as demonstrated in the next section.

VII-F Searching for Competitive SINR Targets (Fig. 10)

In this section, random filter and max-SINR filter initializations are considered to determine the lower and upper bounds of SINR targets, respectively. In Fig. 10, the leftmost bars indicate the results of the proposed distributed algorithm with random filter initializations, which is applied in all simulations of earlier sections. The rightmost bars present the results of distributed max-SINR algorithm. The bars on the left and right sides of the middle indicate the results of the proposed distributed algorithm with $1$ iteration and $2$ iterations of linear SINR target search, respectively. As the iteration number is increased, results closer to the SINRs achieved by the distributed max-SINR algorithm can be obtained. In summary, by fixing the transmit power constraints, but by searching for higher feasible SINR targets, higher SINRs can be achieved at the cost of increased power consumptions. This is an important design trade-off in wireless networks. The results in Fig. 10 are obtained over $20$ Monte Carlo channels.

VII-G Distributed Joint Transmit and Relay Beamforming Filter Design (Table IV)

Distributed joint transmit and relay beamforming filter design improves the total power saving to achieve the same sum-SINR at the cost of increased complexity. Along the similar lines of the complexity analysis for transmitter side in Section V-B, the additional complexity due to the added distributed relay beamforming design is obtained as {\mathcal{O}\big{(}B(2N+5)\big{)}} per relay, aka processor. Hence, as seen in Table IV for the $\text{K}3\text{-}\text{M}10\text{-}\text{N}10\text{-}\text{R}3$ network, the additional total complexity for all relays per iteration is given by {\mathcal{O}\big{(}RB(2N+5)\big{)}} or {\mathcal{O}\big{(}3\!\cdot\!6\!\cdot\!(2\!\cdot\!10\!+\!5)\!=\!450\big{)}}, where ${B=6}$ since each user transmits $2$ streams as noted in the beginning of the section. On the other hand, as obtained in Section V-B, the total complexity of distributed transmitter optimization for all streams per iteration is ${\mathcal{O}(B(M\!+\!5)\!=\!6\!\cdot\!(10\!+\!5)\!=\!90)}$. Each relay serves all streams to achieve the SINR targets. Thus, the distributed optimization at the relay side has higher number of constraints that substantially increases the total complexity. For the special case ${M\!=\!N\!\gg\!5}$, the ratio of the total complexities of the relay and transmitter sides is obtained as

[TABLE]

Due to the joint optimization, the average iteration number also increases as seen in Table IV. Therefore, the average total complexities of the distributed transmit and the distributed joint transmit and relay beamforming design algorithms are given as ${\mathcal{O}(16\cdot 90=1440)}$ and ${\mathcal{O}(23\cdot 540=12420)}$, respectively. The first three rows of results in Table IV are obtained via numerical results averaged over $20$ Monte Carlo channels. The total power saving increases as more resources are allocated to the relay side, i.e., ${R>K}$, ${N>M}$, and ${\text{SNR}_{\text{R}}>\text{SNR}_{\text{T}}}$.

The conclusions drawn from the numerical results of distributed transmit beamforming optimization are also valid for distributed joint transmit and relay beamforming optimization. For instance, recall the total power and sum-SINR results of increasing $\text{SNR}_{\text{T}}$ ${12\text{-}42}$, ${21\text{-}42}$, ${42\text{-}42}$ dB and increasing $\text{SNR}_{\text{R}}$ ${12\text{-}12}$, ${12\text{-}21}$, ${12\text{-}42}$ dB points in Fig. 3. For the former case, higher sum-SINRs are achieved with lower power consumptions. On the other hand, for the latter case, increasing $\text{SNR}_{\text{R}}$ simply increases both the total power consumption and the sum-SINR. Clearly, for distributed transmit beamforming optimization, providing more resources to the transmitter side returns better results. On the other hand, for distributed joint transmit and relay beamforming optimization, providing more resources to the dominant side returns better results. For instance, when the relay side is dominant ${R\gg K}$, clearly, increasing $\text{SNR}_{\text{R}}$ returns better results than increasing $\text{SNR}_{\text{T}}$.

VIII Conclusion

A distributed ADMM algorithm is proposed to design transmit beamforming matrices for a generic wireless relay network that has been hardly studied in the literature due to the challenges raised by the coexistence of multi-stream transmissions, multiple multi-antenna nodes, and the presence of direct links. The traits of the proposed algorithm are low complexity, iteration number, and message exchange loads. Due to the challenge that the direct links bring, an approximate SINR formulation at the relay side is proposed to design distributed joint transmit and relay beamforming filters that further improve the total power saving at the cost of increased complexity.

Appendix A Closed-Form Solutions of MAX-SINR Filters

Both downlink and uplink iterations of the distributed max-SINR algorithm are based on the covariance matrix of interference signals plus noise, which is second-order statistics information. The covariance matrix of downlink direction is similar to the one given for interference channels, i.e., there are no relay nodes, in [27]. In downlink direction, based on the covariance matrix at stream $(k,l)$

[TABLE]

the receive filter is obtained as

[TABLE]

For the uplink iteration, we initially need to obtain the reciprocal channel $\overleftarrow{\mathbf{H}}_{k.i}^{\prime}$ of $\mathbf{H}_{k.i}^{\prime}$ in (5). The effective channel in the uplink direction from transmitter $i$ to receiver $k$ through all relays is given as

[TABLE]

where $\overleftarrow{\mathbf{G}}_{ri}\triangleq{\mathbf{G}}_{ir}^{H}$ and $\overleftarrow{\mathbf{H}}_{kr}^{\prime\prime}\triangleq\mathbf{H}_{rk}^{H}$ are the ${N_{r}\times M_{i}}$ and ${M_{k}\times N_{r}}$ channel matrices between transmitter $i$ and relay $r$, and relay $r$ and receiver $k$, respectively. The covariance matrix of uplink direction can be obtained by interpreting the time domain, i.e., the two time slots due to the relay network communication, as the space domain. Thus, the covariance matrix at stream $(k,l)$ is given as

[TABLE]

where $\overleftarrow{\mathbf{J}}_{kj}\triangleq\mathbf{J}_{jk}^{H}$ is the $M_{k}\times M_{j}$ channel matrix for the direct link between transmitter $j$ and receiver $k$, and

[TABLE]

is the noise covariance matrix. Then,

[TABLE]

and finally, the receive filter is obtained as

[TABLE]

As detailed in [27], after obtaining the receive filter (56) in the downlink direction, the uplink direction is iterated by setting the transmit filter as ${\overleftarrow{\mathbf{u}}_{k,l}=\mathbf{v}_{k,l}}$. After obtaining the receive filter (61) in the uplink direction, the downlink direction is re-iterated by setting the transmit filter as ${\mathbf{u}_{k,l}=\overleftarrow{\mathbf{v}}_{k,l}}$. The uplink and downlink iterations are continued until the convergence.

Appendix B Derivation Details of SINR Approximation (48)

The term ${\zeta_{k,l}^{i,n}(2)}$ in (47) can be rewritten as

[TABLE]

where

[TABLE]

$\otimes$ and $\text{vec}(.)$ denotes the Kronecker product and denotes staking the columns of a matrix in a column vector, respectively. Therefore,

[TABLE]

where

[TABLE]

Now consider the following term in (9)

[TABLE]

where

[TABLE]

$\mathbf{e}_{n}$ is an $N\times 1$ zero vector except 1 at the $n^{th}$ row. Therefore,

[TABLE]

where ${t_{kls}\triangleq\mathbf{t}_{kls}^{H}\mathbf{t}_{kls}}$. Hence, the noise power in the second time slot is

[TABLE]

where

[TABLE]

Hence, the approximate SINR as a function of relay filter $\mathbf{f}_{r}$ is given as

[TABLE]

where

[TABLE]

Furthermore,

[TABLE]

where

[TABLE]

Note that ${\widehat{\text{SINR}}_{k,l}(\mathbf{f}_{r})}$ is in an inhomogeneous quadratic form. To rewrite it in homogeneous quadratic form, let

[TABLE]

Thus, (48) is obtained.