Channel Estimation By Transmitting Pilots From Reconfigurable Intelligent Surface

The emerging reconfigurable intelligent surface (RIS) [1], which is also known as intelligent reflecting surface (IRS) [2], has attracted upsurging interests from both academia and industry, and is envisioned as a promising technology for the next-generation wireless communication systems [3]. The RIS can boost the network communication performance from various aspects while consuming rather low energy. Potentials of RIS have been extensively explored recently and many exciting applications can be found in [1]-[3] and the references therein.

Besides its versatility, one most critical aspect of RIS technology is the channel state information (CSI) acquisition. CSI is essential to configure RIS’ phase-shifting, which, however, is practically difficult to obtain for RIS due to its incapability of transmitting, receiving or processing signals. Therefore, a multitude of recent works are dedicated to exploring channel estimation (CE) techniques suitable for RIS, e.g., the references in [4]. Among them, most available methods belong to the so-called cascaded channel (CscdChn) estimation scheme, e.g., [5]-[13]. That is, the effective channel cascading the two-hop links by way of RIS is measured as a whole during the training procedure. For instance, the work [5] proposes to turn on one single RIS element at a time and estimate each CscdChn one after another. In [6], the authors design the pilot sequence according to the minimum variance unbiased estimation rule. The authors of [7] decompose the CscdChn into rank-one subchannels and design pilot signals via adjusting RIS phase-shifting to minimize the mean square error (MSE) of the estimate of subchannels. The paper [8] proposes both linear minimum mean square error (LMMSE) based and deep learning based pilot design methods to improve the CE precision of the CscdChn and shows that discrete Fourier transform (DFT) pilots can achieve nearly optimal MSE. The authors of [9] wisely exploit the fact that the link between the base station (BS) and RIS are the same for all users’ CscdChns to reduce the overall training overhead of multi-user system. The work [10] studies the RIS training method for orthogonal frequency division multiplexing (OFDM) system via optimizing the training phase-shifts and demonstrates that the DFT sequence can achieve the optimal performance. The authors of [11] propose a hierarchical training scheme via grouping the RIS elements and performing CE progressively using discrete phase shifters. Besides, a line of research, e.g., [12] and [13], explores to leverage the hidden sparsity in channel for high frequency/large antenna array systems to decrease training overhead.

Although the CscdChn training method predominates the existing literature, there still exist a number of works exploring active methods at the RIS side for channel training [14]-[18]. For example, the work [14] proposes a hybrid RIS structure via embedding sensors in element array to measure channels. The authors of [15] install one single receive (RX) RF-chain in RIS and utilize it to develop a compressive sensing (CS) based channel recovery algorithm via exploiting channel’s sparsity. In [16], active sensing elements are installed on the uniform planar element array of the RIS and arranged in “L” shape to sense the channels. The latest works [17] and [18] consider similar hybrid RIS architecture as proposed in [14] and utilize tensor decomposition and atomic norm methods to recover the sparse channels, respectively. Besides, some works investigate extraction of one-hop CSI of BS-RIS or RIS-user through the reflected pilot signals. For instance, the authors of [19] and [20] propose factor decomposition and matrix completion methods to achieve this goal. The work [21] exploits matrix calibration theory and channel sparsity to recover separate RIS channels.

Along with the deepening research on RIS technology, some critical insights have been obtained very recently. Especially, as pointed out by [22], one predominant defect of RIS is the severe fading loss of the concatenated channels. This feature is named as double-fading effect in [22] and derives from the multiplicative nature of the CscdChn’s attenuation. As analyzed in [22], the attenuation loss experienced by the two-hop channels going through RIS is generally several orders of magnitude larger than that of the direct channel. This effect has also been substantiated by numerical experiments and real-field tests reported in the latest works [23] and [24]. At the same time, according to Cramér-Rao lower bound (CRLB) analysis presented in [1] and [25], the MSE of CE is inversely proportional to the signal-to-noise-ratio (SNR). Therefore, the double-fading loss indeed severely weakens the pilot signals’ receiving SNR and hence bottlenecks the CE performances of the classical CscdChn training schemes, e.g., [5]-[13], and the one-hop CSI extraction methods relying on the reflected pilot signals, e.g., [19]-[21].

Based on the above inspections, to improve the CE precision, one natural thought is to breakthrough the double-fading curse. In fact, a number of works have made efforts toward this end. For instance, the authors of [23] and [24] have proposed a novel active-RIS architecture by introducing amplifiers into RIS elements to magnify the impinging signals. At the same time, the works [14]-[18] introduce a small number of RX RF-chains at the RIS side to perform signal receiving locally.

Enlightened by all the above works, this paper proposes a novel RIS channel training scheme. Specifically, we equip the RIS with one single TX RF-chain [15], [26], as shown in Fig. 1. During the channel training period, the TX RF-chain is connected to all the phase shifters [26], which enables the RIS to broadcast pilot signals to the BS and all users. This novel scheme cherishes the following advantages:

• •

  This scheme can overcome the double-fading curse. In fact, only “one-hop” channels are estimated.

• •

  Only one TX RF-chain is required, which can largely maintain the hardware efficiency of RIS device.

• •

  The BS and all users conduct CE simultaneously. This makes the RIS channel training overhead independent of the number of users.

Motivated by the above observations, this paper proposes a novel RIS-transmitting (RIS-TX) channel training scheme as specified above. We conduct a comprehensive research on this new scheme. Specifically, the contributions of this paper are specified as follows:

• •

  Firstly, we put forward a brand new RIS-TX CE scheme. Via incorporating one TX RF-chain into RIS device, we propose that RIS broadcasts “pilot signals”, which are indeed implemented by RIS elements’ predefined phase shifting, to the BS and all mobile users. This novel CE scheme’s pilot overhead is lower than the conventional CscdChn methods, e.g., [4]-[11], [19], [20], and can effectively overcome double-fading effect. To the best of authors’ knowledge, similar RIS-TX scheme has never been considered in the existing literature, e.g., [4]-[26].

• •

  Next, we study the associated pilot sequence design problem. Specifically, we aim at minimizing the estimation MSE of the effective RIS channels of all users via designing the pilot sequences. This task turns out to be a nonconvex quartic optimization problem and is much more challenging than those in the existing relevant literature, e.g., [6]-[11]. To tackle this difficult problem, we successfully develop a gradient descent (GD)-based solution, which exhibits fast convergence.

• •

  Besides, we also develop an alternative pilot design solution by exploiting the cutting-the-edge penalty duality decomposition (PDD) framework [27]. This PDD-based solution has all its block coordinates updated either in closed form or by solving a structured linear equation and consequently runs highly efficiently.

• •

  Furthermore, we also theoretically analyze the performance of our RIS-TX training scheme. CRLB analysis has uncovered one important fact that the number of RIS elements exerts opposite impact on the RIS-TX scheme and the prevailing CscdChn counterpart. The former yields superior estimation accuracy when RIS has moderate size (e.g. $\leq$ 10000 elements) while the latter is more beneficial for large-scale RIS. Asymptotic analysis reflects that the two competing CE schemes exhibit distinct scaling laws in power and RIS size. These insights provide valuable guidelines for selecting between the RIS-TX and CscdChn scheme in various application scenarios, which have never been discussed in the existing literature, e.g., [4]-[24].

• •

  Last but not least, extensive numerical results demonstrate the effectiveness of our proposed algorithms and the benefit of the novel RIS-TX training scheme. Very interestingly, the pilot sequences optimized by our GD/PDD-based algorithms outperform the DFT pilot sequence, which has been reported to be optimal for the CscdChn scheme by different literature [6]-[8], [10], [11].

As illustrated in Fig. 1, we consider a narrowband time division duplex (TDD) cellular system which consists of a BS equipped with $M_{\mathrm{B}}$ antennas, $K$ single-antenna users and a RIS containing $M_{\mathrm{R}}$ reflecting elements. The sets of users and RIS elements are denoted by $\mathcal{K}\triangleq\{1,\cdots,K\}$ and $\mathcal{M}\triangleq\{1,\cdots,M_{\mathrm{R}}\}$, respectively. Moreover, we denote $\mathbf{G}\in\mathbb{C}^{M_{\mathrm{B}}\times M_{\mathrm{R}}}$ and $\mathbf{h}_{k}\in\mathbb{C}^{M_{\mathrm{R}}},\;k\in\mathcal{K}$, as the channel between the BS and the RIS and the one between the RIS and the $k$th user, respectively. Note that our setting does not consider the direct channels connecting the BS and users for simplicity, which can be easily acquired by following the standard CE procedure for conventional cellular system without RIS [28].

As previously discussed, one RF-chain is built in the RIS device to empower it with transmission capability, as shown in Fig. 1. During the channel training period, every RIS element is switched to connecting the single RF-chain and adjusts the signal’s phase before emitting it over the air. Specifically, denote the RF-chain signal as a complex scalar $\sqrt{\alpha}e^{j\theta_{0}}$ with $\sqrt{\alpha}$ and $\theta_{0}$ representing its amplitude and phase, respectively. Therefore, the emitted signal from the RIS can be expressed as $\bm{\psi}\triangleq[\sqrt{\alpha}e^{j(\theta_{0}+\theta_{1})},\cdots,\sqrt{\alpha}e^{j(\theta_{0}+\theta_{M_{\mathrm{R}}})}]^{\mathrm{T}}$ with $\theta_{i},\;i\in\mathcal{M}$, being the phase shift yielded by the $i$th element. Since the pilot signal is predefined, we can just assume $\theta_{0}=0$ without loss of generality.

Based on the above, we elaborate the proposed RIS-TX CSI acquisition procedure. The training period is composed of $T$ consecutive time slots and is assumed to be shorter than the channel’s coherence time. Within each time slot, the RIS broadcasts one pilot symbol, which is phase-modulated by the RIS elements, to the BS and all users. Specifically, denoting $\mathcal{T}\triangleq\{1,\cdots,T\}$, we can express the transmitted signal in the $t$th time slot as [7]-[11]

where $\bm{\phi}^{(t)}\triangleq[e^{j\theta_{1}^{(t)}},\cdots,e^{j\theta_{M_{\mathrm{R}}}^{(t)}}]^{\mathrm{T}}$.

After collecting all $T$ pilot symbols from the RIS simultaneously, the BS and all users conduct signal processing and obtain the estimate of their own channel connecting the RIS. This procedure will be detailed in the sequel.

In the $t$th time interval, the BS receives the signal from the RIS which is given as

where $\mathbf{n}_{\mathrm{B}}^{(t)},\;t\in\mathcal{T}$, represents the thermal noise at the BS and $\mathbf{n}_{\mathrm{B}}^{(t)}\sim\mathcal{CN}(\bm{0},\bm{\Sigma}_{\mathrm{B}})$ with the positive definite matrix $\bm{\Sigma}_{\mathrm{B}}$ denoting the covariance matrix of the noise which is assumed to be known [8], [9], [29].^1††1Note the noise statistics generally vary slowly in a large timescale [30], which hence can be obtained at a low expense of time overhead. It is also reasonable to assume that the noise covariance matrix is invariant during the coherence time, which indicates that different $\mathbf{n}_{\mathrm{B}}^{(t)}$’s have one same covariance matrix, i.e., $\bm{\Sigma}_{\mathrm{B}}$. Besides, we assume that $\mathbf{n}_{\mathrm{B}}^{(t)}$’s with different $t$ are independent for simplicity [6]-[11]. To fully determine the channel $\mathbf{G}$, the length of the pilot sequence $T$ should satisfy $T\geq M_{\mathrm{R}}$. To reduce training’s overhead, we just assume that $T=M_{\mathrm{R}}$. By collecting all $T$ observations and stacking them in a column-by-column manner, the training data obtained by the BS can be represented as

where $\mathbf{Y}_{\mathrm{B}}\triangleq[\mathbf{y}_{\mathrm{B}}^{(1)},\cdots,\mathbf{y}_{\mathrm{B}}^{(T)}]$, $\bm{\Psi}\triangleq[\bm{\psi}^{(1)},\cdots,\bm{\psi}^{(T)}]$, $\mathbf{N}_{\mathrm{B}}\triangleq[\mathbf{n}_{\mathrm{B}}^{(1)},\cdots,\mathbf{n}_{\mathrm{B}}^{(T)}]$. Via vectorizing $\mathbf{Y}_{\mathrm{B}}$ and utilizing the formula $\mathsf{vec}(\mathbf{ABC})=(\mathbf{C}^{\mathrm{T}}\otimes\mathbf{A})\mathsf{vec}(\mathbf{B})$, we obtain

where $\mathbf{g}\triangleq\mathsf{vec}(\mathbf{G})$, $\mathbf{n}_{\mathrm{B}}\triangleq\mathsf{vec}(\mathbf{N}_{\mathrm{B}})$, $\mathsf{A}_{\mathrm{G}}(\bm{\Psi})\triangleq\bm{\Psi}^{\mathrm{T}}\otimes\mathbf{I}_{M_{\mathrm{B}}}$ and $\mathbf{I}_{M_{\mathrm{B}}}$ is $M_{\mathrm{B}}\times M_{\mathrm{B}}$ identity matrix.

Similarly, the compact form of the signals observed by the $k$th user during the whole training period can be expressed as

where $\mathbf{n}_{\mathrm{U},k}\triangleq[n_{\mathrm{U},k}^{(1)},\cdots,n_{\mathrm{U},k}^{(T)}]^{\mathrm{T}}$ and $n_{\mathrm{U},k}^{(t)},\;k\in\mathcal{K},\;t\in\mathcal{T}$, indicates the thermal noise at the $k$th user in the $t$th interval with $n_{\mathrm{U},k}^{(t)}\sim\mathcal{CN}(0,\sigma_{\mathrm{U},k}^{2})$. For simplicity, it is assumed that each $n_{\mathrm{U},k}^{(t)}$ is independent to others for different users and time slots.

In this paper, we assume that the channels $\mathbf{g}$ and $\{\mathbf{h}_{k}\}_{k=1}^{K}$ have zero mean and known covariance matrices. Specifically, we denote $\mathbf{R}_{\mathrm{G}}$ and $\mathbf{R}_{\mathrm{h},k},\;k\in\mathcal{K}$, as covariance matrices of $\mathbf{g}$ and $\mathbf{h}_{k},\;k\in\mathcal{K}$, respectively. In reality, channel statistics generally vary slowly and can be easily obtained [29]. For instance, the covariance matrix can be empirically estimated at low overhead [31], [32]. Besides, the zero mean assumption can also be readily satisfied via subtracting their nonzero mean values. Besides, it is reasonable to assume that the channels associated with the BS and users are mutually uncorrelated. With the linear receivers $\mathbf{W}_{\mathrm{G}}$ at the BS and $\mathbf{W}_{\mathrm{h},k},\;k\in\mathcal{K}$, at the $k$th user, the BS/every user obtains the linear estimate of their own channels given as follows

whose covariance matrices are given as follows after some manipulations

where $\bm{\Sigma}_{\mathrm{U},k}\triangleq\mathsf{Diag}(\sigma_{\mathrm{U},k}^{2},\cdots,\sigma_{\mathrm{U},k}^{2}),\;k\in\mathcal{K}$.

After performing CE, mobile users feedback their local estimates to the BS. The effective CscdChn by way of the RIS for each user is defined by $\mathbf{H}_{\mathrm{c},k}\triangleq\mathbf{G}\mathsf{Diag}(\mathbf{h}_{k}),\;k\in\mathcal{K}$, and its estimate is given as $\hat{\mathbf{H}}_{\mathrm{c},k}\triangleq\hat{\mathbf{G}}\mathsf{Diag}(\hat{\mathbf{h}}_{k}),\;k\in\mathcal{K}$. Based on $\{\hat{\mathbf{H}}_{\mathrm{c},k}\}_{k=1}^{K}$, the BS will determine the scheduling and beamforming for all users. Therefore, it is critical to minimize the estimation error between the estimate $\{\hat{\mathbf{H}}_{\mathrm{c},k}\}_{k=1}^{K}$ and their true values. To evaluate the estimation performance, the MSE matrix of $\hat{\mathbf{H}}_{\mathrm{c},k},\;k\in\mathcal{K}$, is defined as $\mathbf{C}_{\mathrm{H}_{\mathrm{c}},k}\triangleq\mathbb{E}\{\mathsf{vec}(\mathbf{H}_{\mathrm{c},k}-\hat{\mathbf{H}}_{\mathrm{c},k})\mathsf{vec}^{\mathrm{H}}(\mathbf{H}_{\mathrm{c},k}-\hat{\mathbf{H}}_{\mathrm{c},k})\},\;k\in\mathcal{K}$, and its value is determined by the following theorem, which is derived in Appendix A.

Hence, the estimation MSE of the effective channel $\hat{\mathbf{H}}_{\mathrm{c},k},\;k\in\mathcal{K}$, is given as $\mathsf{Tr}\{\mathbf{C}_{\mathrm{H}_{\mathrm{c},k}}\},\;k\in\mathcal{K}$, which is actually a function of the pilot sequence $\bm{\Psi}$. In the subsequent exposition, we define $\bm{\Psi}=\bm{\Phi}\mathbf{P}$ with $\bm{\Phi}\triangleq[\bm{\phi}^{(1)},\cdots,\bm{\phi}^{(T)}]^{\mathrm{T}}$ and $\mathbf{P}\triangleq\mathsf{Diag}(\sqrt{\alpha_{1}},\cdots,\sqrt{\alpha_{T}})$. Besides, we denote $\bm{\mathcal{W}}_{\mathrm{h}}\triangleq\{\mathbf{W}_{\mathrm{h},k}\}_{k=1}^{K}$.

Based on the above exposition, we aim at minimizing the overall MSE of all users’ effective channels $\{\hat{\mathbf{H}}_{\mathrm{c},k}\}_{k=1}^{K}$ via designing the pilot sequence and the linear receivers $(\bm{\Phi},\mathbf{P},\mathbf{W}_{\mathrm{G}},\bm{\mathcal{W}}_{\mathrm{h}})$. This task is formulated as the following problem

where $P_{\mathrm{max}}$ in (11a) represents the overall transmission power budget for training. The problem $(\mathcal{P}1)$ is highly challenging due to its nonconvex objective. Especially, it is a nonconvex quartic function in $(\bm{\Phi},\mathbf{P})$. In the following, we propose two algorithms to deal with $(\mathcal{P}1)$.

In the following, we resolve $(\mathcal{P}1)$ in a block coordinate descent (BCD) [33] manner.

1) Optimize the phase shift matrix $\bm{\Phi}$

Fixing other variables, the optimization w.r.t. $\bm{\Phi}$ is given by

which is non-convex due to the equality constraint (12a) and challenging to solve. Hence, we turn to investigate another optimization equivalent to $(\mathcal{P}2)$ which can be expressed as

where $\bm{\Theta}\triangleq\angle(\bm{\Phi})$ and $\angle(\bm{\Phi})$ takes the value of phases of $\bm{\Phi}$’s entries in an element-wise manner. By noting that $e^{j(\theta+2n\pi)}=e^{j\theta},\;n\in\mathbb{Z}$, problem $(\mathcal{P}3)$ can be viewed as an unconstrained optimization problem w.r.t. $\bm{\Theta}$. Therefore, we adopt GD algorithm to tackle $(\mathcal{P}3)$. Denote $\mathsf{g}(\bm{\Theta})$ as the objective of $(\mathcal{P}3)$. Then, the gradient descent of $\bm{\Theta}$ is conducted as

where $\eta_{\bm{\Theta}}^{(n)}$ is the step size for updating $\bm{\Theta}$. Generally, constant or diminishing step size rule can yield satisfactory convergence [34]. We relegate the derivation of $\nabla_{\bm{\Theta}}\mathsf{g}(\bm{\Theta})$ to Appendix B. Lastly, $\bm{\Phi}^{(n+1)}$ is given by $\bm{\Phi}^{(n+1)}=e^{j\bm{\Theta}^{(n+1)}}$.

2) Optimize the power allocation $\mathbf{P}$

We proceed to study the optimization of power allocation $\mathbf{P}$ when other variables are given. Notice again $\mathbf{P}\triangleq\mathsf{Diag}(\sqrt{\alpha_{1}},\cdots,\sqrt{\alpha_{T}})$. Therefore, by defining $\mathbf{p}\triangleq\mathsf{diag}(\mathbf{P})$, the sub-problem to optimize power allocation can be expressed as

The objective of $(\mathcal{P}4)$ is non-convex and we still adopt GD-like method to attack it. Note that plain GD algorithm only applies to unconstrained optimization problem. For the constrained problem $(\mathcal{P}4)$, we adopt the gradient projection (GP) method [34]. GP method is an iterative procedure. In each iteration, it performs two steps in order: i) gradient descent and ii) projection onto the feasible domain. Denoting the objective of $(\mathcal{P}4)$ as $\mathsf{h}(\mathbf{p})$, we elaborate these two steps in the following:

• •

  Gradient Descent: The variable $\mathbf{p}$ moves in the direction of negative derivative, i.e.

  $$\mathbf{p}^{(n+1)}:=\mathbf{p}^{(n)}-\eta_{\mathrm{p}}^{(n)}\nabla_{\mathbf{p}}\mathsf{h}(\mathbf{p})|_{\mathbf{p}=\mathbf{p}^{(n)}},$$

  (16)

  where $\eta_{\mathrm{p}}^{(n)}$ is the step size for updating power allocation, which can be set by following the step size rules presented in [34]. The calculation of $\nabla_{\mathbf{p}}\mathsf{h}(\mathbf{p})$ is shown in Appendix C.

• •

  Projection: After performing the update of $\mathbf{p}$ as shown above, there exists a risk that $\mathbf{p}^{(n+1)}$ moves out of the feasible domain identified by (15a) and (15b). Hence, we need project $\mathbf{p}^{(n+1)}$ onto the convex feasible domain of $(\mathcal{P}4)$ [34], which means to solve the following optimization

  	$\displaystyle(\mathcal{P}5):\min_{\mathbf{x}\in\mathbb{R}^{T}}\;$	$\displaystyle\frac{1}{2}\|\mathbf{x}-\mathbf{p}^{(n+1)}\|_{2}^{2}$		(17)
  	$\displaystyle\mathrm{s.t.}\;$	$\displaystyle\mathbf{x}^{\mathrm{T}}\mathbf{x}\leq\frac{P_{\mathrm{max}}}{M_{\mathrm{R}}},$		(17a)
  		$\displaystyle\mathbf{x}\geq\bm{0}.$		(17b)

  The optimal solution to $(\mathcal{P}5)$ is given by the following Theorem 2, whose proof is left to Appendix D.

  Theorem 2.

  If $\|[\mathbf{p}^{(n+1)}]_{+}\|_{2}^{2}\leq\frac{P_{\mathrm{max}}}{M_{\mathrm{R}}}$, we have

  $$\mathbf{x}^{\star}=[\mathbf{p}^{(n+1)}]_{+},$$

  (18)

  otherwise,

  $$\mathbf{x}^{\star}=\sqrt{\frac{P_{\mathrm{max}}}{M_{\mathrm{R}}}}\frac{[\mathbf{p}^{(n+1)}]_{+}}{\|[\mathbf{p}^{(n+1)}]_{+}\|_{2}},$$

  (19)

  where $[\mathbf{p}^{(n+1)}]_{+}\triangleq\max\{\mathbf{p}^{(n+1)},\bm{0}\}$ in an element-wise manner.

3) Optimize the linear receiver at the BS $\mathbf{W}_{\mathrm{G}}$

When other variables are given, the subproblem w.r.t. $\mathbf{W}_{\mathrm{G}}$ is indeed an unconstrained convex quadratic optimization problem. Therefore, by checking the first-order optimality condition of $\mathbf{W}_{\mathrm{G}}$, the newly obtained $\mathbf{W}_{\mathrm{G}}$ can be expressed in the closed-form as follows

where

with $\mathsf{Ddiag}(\mathbf{X})$ constructing a diagonal matrix with its diagonal elements being those of the square matrix $\mathbf{X}$. The derivation of (20) is detailed in Appendix E.

4) Optimize the linear receiver at every user $\bm{\mathcal{W}}_{\mathrm{h}}$

Following the similar procedure of updating $\mathbf{W}_{\mathrm{G}}$, the solution to $\mathbf{W}_{\mathrm{h},k},\;k\in\mathcal{K}$, can also be obtained in a closed-form as follows (details are omitted for brevity)

where

$\mathbf{D}_{\mathrm{h},1}$ and $\mathbf{D}_{\mathrm{h},2}$ are $M_{\mathrm{R}}$-dimensional diagonal matrices with $[\mathbf{D}_{\mathrm{h},1}]_{m,m}=\sum_{n=1}^{M_{\mathrm{B}}}[\mathbf{W}_{\mathrm{G}}^{\mathrm{H}}\mathsf{A}_{\mathrm{G}}(\bm{\Psi})\mathbf{R}_{\mathrm{G}}]_{(m-1)M_{\mathrm{B}}+n,(m-1)M_{\mathrm{B}}+n},\;m\in\mathcal{M}$, and $[\mathbf{D}_{\mathrm{h},2}]_{m,m}=\sum_{n=1}^{M_{\mathrm{B}}}[\mathbf{C}_{\hat{\mathrm{G}}}]_{(m-1)M_{\mathrm{B}}+n,(m-1)M_{\mathrm{B}}+n},\;m\in\mathcal{M}$, respectively.

The overall algorithm of solving $(\mathcal{P}1)$ based on GD is summarized in Algorithm 1.

One potential defect of the previously discussed GD-based solution is its onerous gradient evaluation procedures, as reflected in Appendix B $\&$ C. In this subsection, we adopt the PDD methodology [27] to propose an alternative solution that updates all block coordinates highly efficiently. To this end, we first transform the difficult objective in (11) into a more tractable form amiable to block coordinate update via introducing auxiliary variables. Specifically, we introduce multiple copies of $\bm{\Phi}\mathbf{P}$ by defining

and rewrite the original complicated objective into simple forms of the newly introduced copies as follows

Therefore, the problem $(\mathcal{P}1)$ can be equivalently rewritten as

where the newly introduced notations above are defined as

Following the PDD method [27], to tackle $(\mathcal{P}6)$, we turn to solve its augmented Lagrangian problem given as follows

The PDD method is a two-layer iterative procedure [27], with its inner layer alternatively optimizing the different blocks of variables comprising $\bm{\mathcal{V}}$ and its outer layer selectively updating the dual variables $\bm{\Lambda}_{i}$, $i=1,2,3,4$ or the penalty coefficient $\rho$. The two layers’ updating procedure will be specified in the sequel.

For the inner layer, we adopt the block successive upper-bound minimization (BSUM) to update various blocks [27], [35], with each block’s update being elaborated as follows.

1) Optimize the phase shift matrix $\bm{\Phi}$

When other variables are given, the subproblem w.r.t. $\bm{\Phi}$ can be expressed as

Note that the quadratic term $\|\bm{\Phi}\mathbf{P}\|_{\mathrm{F}}^{2}$ reduces to a constant $T\sum_{t=1}^{T}\alpha_{t}$ due to (35a), which is independent of $\bm{\Phi}$. Hence, utilizing the above fact, the problem $(\mathcal{P}8)$ boils down to maximizing a linear objective given as follows

where the parameter $\mathbf{Q}$ is defined as

and the $(i,j)$th element of the matrix $\tilde{\mathbf{C}}_{\bm{\Phi}}$ is given by $\sum_{m=1}^{M_{\mathrm{B}}}[(\mathbf{P}\!\otimes\!\mathbf{I}_{M_{\mathrm{B}}})(\mathbf{U}_{2}^{\mathrm{T}}\!+\!\rho\bm{\Lambda}_{2}^{\mathrm{T}})]_{(i\!-\!1)M_{\mathrm{B}}\!+\!m,(j\!-\!1)M_{\mathrm{B}}\!+\!m},\;i,j\in\mathcal{M}$.

The problem $(\mathcal{P}9)$, although nonconvex, achieves optimum when the elements’ phases of $\bm{\Phi}$ align with those of $\mathbf{Q}^{\mathrm{H}}$. Consequently, the optimal solution of $(\mathcal{P}8)$ is given analytically by

2) Optimize the power allocation $\mathbf{P}$

After some manipulations, the update of $\mathbf{P}$ with the remaining variables fixed can be written as

where $\mathbf{p}\!\triangleq\!\mathsf{diag}({\mathbf{P}})$, $\mathbf{q}^{\mathrm{H}}\triangleq\big{(}\mathsf{diag}\big{(}\sum_{k=1}^{K}(\mathbf{U}_{1,k}^{\mathrm{T}}+\rho\bm{\Lambda}_{1,k}^{\mathrm{T}})\bm{\Phi}\big{)}\big{)}^{\mathrm{T}}+\mathbf{c}_{2}^{\mathrm{T}}$ and $[\mathbf{c}_{2}]_{m}=\sum_{n=1}^{M_{\mathrm{B}}}[(\mathbf{U}_{2}^{\mathrm{T}}+\rho\bm{\Lambda}_{2}^{\mathrm{T}})(\bm{\Phi}\otimes\mathbf{I}_{M_{\mathrm{B}}})]_{(m-1)M_{\mathrm{B}}+n,(m-1)M_{\mathrm{B}}+n},\;m\in\mathcal{M}$.

Following the similar procedure as solving $(\mathcal{P}5)$, the optimal solution to $(\mathcal{P}10)$ is given by

• •

  if $\big{\|}\frac{[\mathsf{Re}\{\mathbf{q}\}]_{+}}{KM_{\mathrm{R}}+M_{\mathrm{B}}M_{\mathrm{R}}}\big{\|}_{2}^{2}\leq\frac{P_{\mathrm{max}}}{M_{\mathrm{R}}}$, we have

  $$\mathbf{p}^{\star}=\frac{[\mathsf{Re}\{\mathbf{q}\}]_{+}}{KM_{\mathrm{R}}+M_{\mathrm{B}}M_{\mathrm{R}}},$$

  (40)

• •

  otherwise,

  $$\mathbf{p}^{\star}=\sqrt{\frac{P_{\mathrm{max}}}{M_{\mathrm{R}}}}\frac{[\mathsf{Re}\{\mathbf{q}\}]_{+}}{\|[\mathsf{Re}\{\mathbf{q}\}]_{+}\|_{2}},$$

  (41)

which can be proven by exploiting Appendix D.

3) Optimize the auxiliary variable $\bm{\mathcal{U}}_{1}$

Obviously, the optimization w.r.t. $\bm{\mathcal{U}}_{1}$ is an unconstrained convex quadratic minimization problem. Hence, by setting the first-order derivative of (33) w.r.t. $\mathbf{U}_{1,k},\;k\in\mathcal{K}$, to zero, we attain the following equation

where

The linear system in (42) is the well-known Sylvester equation and can be readily solved via the seminal Hessenburg-Schur algorithm, which has the complexity of $\mathcal{O}(M_{\mathrm{R}}^{3})$ [36].^2††2Note in MATLAB, the standard build-in function $\mathsf{sylvester}()$ can be invoked to solve (42) immediately.

4) Optimize the auxiliary variable $\mathbf{U}_{2}$

Following the similar arguments as in updating $\bm{\mathcal{U}}_{1}$, the optimization of $\mathbf{U}_{2}$ also reduces to solving a Sylvester equation given as follows

where

5) Optimize the auxiliary variable $\bm{\mathcal{U}}_{3}$

Next, we proceed to study the update of $\bm{\mathcal{U}}_{3}$, which is an unconstrained convex quadratic problem, whose closed-form solution can be obtained via checking the first-order optimality condition, shown as

where $\tilde{\mathbf{C}}_{\mathrm{d},1}$ and $\tilde{\mathbf{C}}_{\mathrm{d},2}$ are $M_{\mathrm{R}}$-dimensional diagonal matrices whose $i$th diagonal elements, $i\in\mathcal{M}$, equal to $\sum_{j=(i-1)M_{\mathrm{B}}+1}^{iM_{\mathrm{B}}}[\mathbf{U}_{4}^{\mathrm{H}}\mathbf{R}_{\mathrm{G}}^{\frac{1}{2}}]_{jj}$ and $\sum_{j=(i-1)M_{\mathrm{B}}+1}^{iM_{\mathrm{B}}}[\mathbf{U}_{4}^{\mathrm{H}}\mathbf{U}_{4}+\mathbf{W}_{\mathrm{G}}^{\mathrm{H}}(\bm{\Sigma}_{\mathrm{B}}\otimes\mathbf{I}_{M_{\mathrm{R}}})\mathbf{W}_{\mathrm{G}}]_{jj}$, respectively.

6) Optimize the auxiliary variable $\mathbf{U}_{4}$

Similar to the optimization of $\bm{\mathcal{U}}_{3}$, the solution to $\mathbf{U}_{4}$ can be expressed in the closed-form as

where

7) Optimize the linear receiver at every user $\bm{\mathcal{W}}_{\mathrm{h}}$

The update of $\bm{\mathcal{W}}_{\mathrm{h}}$ is an unconstrained convex quadratic optimization problem, whose optimal solution is equivalent to solving the following Sylvester equation

where

8) Optimize the linear receiver at the BS $\mathbf{W}_{\mathrm{G}}$

The update of $\mathbf{W}_{\mathrm{G}}$ is similar to that of $\bm{\mathcal{W}}_{\mathrm{h}}$ and can be conducted by solving the following Sylvester equation (details are skipped to avoid repetition)

For the outer layer, if all the equalities (24)-(27) are approximately satisfied, we will update the Lagrangian multipliers in a gradient ascent manner, i.e., $\bm{\Lambda}_{i}^{(t+1)}:=\bm{\Lambda}_{i}^{(t)}+\rho^{-1}(\mathbf{U}_{i}-\mathbf{U}_{i}^{\mathrm{R}})$, where $\mathbf{U}_{i}^{\mathrm{R}}$ represents the right hand side of the $i$th equality. Otherwise, the penalty coefficient $\rho$ should be decreased to force the equalities to be approached in the subsequent iterations.

The proposed PDD-based algorithm of handling $(\mathcal{P}1)$ is stated in Algorithm 2.