A) Precoded MU-MIMO OFDM system

We consider uplink communications of IoT environments, which is modeled as a precoded MU-MIMO OFDM system. Figure 1 shows the system model, where the number of transmit terminals is $N$, the number of receiving antennas at the base station is $M$, and the number of subcarriers is $Q$. Given that the number of transmit terminal is typically large in IoT environments, we focus on the overloaded scenario and assume $M<N$ in this paper. The symbol alphabet and the frequency domain transmitted OFDM symbol vector from the $n$-th transmit IoT terminal are denoted by $\mathcal {S}$ and $\boldsymbol {s}_{n}$, respectively. Here, taking IoT environment-specific feature into consideration, we assume only $N_{\text {act}}$ IoT terminals out of $N$ terminals are active meaning that only $N_{\text {act}}$ terminals transmit OFDM signal blocks. Non-active $N-N_{\text {act}}$ terminals actually keep silent, but we can regard they transmit all zero signal block $\textbf{0}_{Q}$. We thus have $\boldsymbol {s}_{n} \in \mathcal {S}^{Q}$ when the $n$-th terminal is active, and otherwise $\boldsymbol {s}_{n}=\textbf{0}_{Q}$. When we use the cyclic prefix with the length greater than or equal to the channel order, the received signal vector after the removal of the cyclic prefix is given by

(1)\begin{align} \hspace{-4mm} \begin{bmatrix} \boldsymbol {r}_{1}^{\text{f}, \text{OFDM}} \\ \vdots \\ \boldsymbol {r}_{M}^{\text{f}, \text{OFDM}} \end{bmatrix} &= \begin{bmatrix} \boldsymbol {\Lambda}^{(1,1)}\boldsymbol {P} & \hspace{-2mm} \dotsb & \hspace{-2mm} \boldsymbol {\Lambda}^{(1,N)}\boldsymbol {P} \\ \vdots & \hspace{-2mm} \ddots & \hspace{-2mm} \vdots \\ \boldsymbol {\Lambda}^{(M,1)}\boldsymbol {P} & \hspace{-2mm} \dotsb & \hspace{-2mm} \boldsymbol {\Lambda}^{(M,N)}\boldsymbol {P} \end{bmatrix} \begin{bmatrix} \boldsymbol {s}_{1} \\ \vdots \\ \boldsymbol {s}_{N} \end{bmatrix} + \begin{bmatrix} \boldsymbol {v}_{1}^{\text{f}} \\ \vdots \\ \boldsymbol {v}_{M}^{\text{f}} \end{bmatrix}, \end{align}

where $\boldsymbol {r}_{m}^{\text {f}, \text {OFDM}} \in \mathbb {C}^{Q}$ is the frequency domain received OFDM signal block at the $m$-th receiving antenna [Reference Hayashi, Nakai, Hayakawa and Ha10]. The diagonal matrix $\boldsymbol {\Lambda }^{(m,n)}={\rm diag} ({\lambda _{1}^{(m,n)},\dotsc ,\lambda _{Q}^{(m,n)}}) \in \mathbb {C}^{Q \times Q}$ is composed of the channel frequency responses with the order of $L-1$ between the $n$-th IoT terminal and the $m$-th receiving antenna. The diagonal elements can be written as

(2)\begin{align} \begin{bmatrix} \lambda_{1}^{(m,n)} \\ \vdots \\ \lambda_{Q}^{(m,n)} \end{bmatrix} &= \sqrt{Q}\boldsymbol {D} \begin{bmatrix} h_{1}^{(m,n)} \\ \vdots \\ h_{L}^{(m,n)} \\ \textbf{0}_{Q-L} \end{bmatrix}, \end{align}

where $\boldsymbol {D} \in \mathbb {C}^{Q \times Q}$ is a $Q$-point unitary DFT matrix defined as

(3)\begin{align} \hspace{-4mm} \boldsymbol {D} &= \frac{1}{\sqrt{Q}} \begin{bmatrix} 1 & \hspace{-2mm} 1 & \hspace{-2mm} \dotsb & \hspace{-2mm} 1 \\ 1 & \hspace{-2mm} e^{-j\frac{2\pi \times 1 \times 1}{Q}} & \hspace{-2mm} \dotsb & \hspace{-2mm} e^{-j\frac{2\pi \times 1 \times (Q-1)}{Q}} \\ \vdots & \hspace{-2mm} \vdots & \hspace{-2mm} & \hspace{-2mm} \vdots \\ 1 & \hspace{-2mm} e^{-j\frac{2\pi \times (Q-1) \times 1}{Q}} & \hspace{-2mm} \dotsb & \hspace{-2mm} e^{-j\frac{2\pi \times (Q-1) \times (Q-1)}{Q}} \\ \end{bmatrix} \end{align}

and $h_{1}^{(m,n)},\dotsc ,h_{L}^{(m,n)}$ denotes the impulse response of the frequency-selective channel between the $n$-th IoT terminal and $m$-th receiving antenna. $\boldsymbol {P} \in \mathbb {C}^{Q \times Q}$ is a precoding matrix required to achieve good detection performance by the convex optimization-based detection scheme. Note that, in [Reference Hayashi, Nakai, Hayakawa and Ha10], we have numerically confirmed that we do not have to use different precoding matrices among IoT terminals, and a common Hadamard precoding matrix leads to good SER performance. $\boldsymbol {v}_{m}^{\text {f}} \in \mathbb {C}^{Q}$ is the frequency domain additive white noise vector at the $m$-th receiving antenna with mean $\textbf{0}_{Q}$ and covariance matrix $\sigma _{\text {v}}^{2}\boldsymbol {I}_{Q}$.

Fig. 1. Uplink MU-MIMO OFDM system for IoT environment.

B) Non-precoded MU-MIMO SC-CP system

Here, we show a non-precoded MU-MIMO SC-CP signal model. Assuming the length of cyclic prefix is greater than or equal to the channel order $L-1$, the time domain received signal block at the $m$-th receiving antenna of the base station is written as

(4)\begin{align} \boldsymbol {r}_{m}^{\text{t}, \text{SC-CP}} &= \sum_{n=1}^{N} \boldsymbol {D}^{\mathsf{H}} \boldsymbol {\Lambda}^{(m,n)} \boldsymbol {D} \boldsymbol {s}_{n} + \boldsymbol {v}_{m}^{{\rm t}}, \end{align}

where $\boldsymbol {v}_{m}^{{\rm t}} \in \mathbb {C}^{Q}$ is the time domain additive white noise vector at the $m$-th receiving antenna having mean $\textbf{0}_{Q}$ and covariance matrix $\sigma _{\text {v}}^{2} \boldsymbol {I}_{Q}$ [Reference Hayashi, Nakai-Kasai and Hayakawa12, Reference Wang and Giannakis20]. By stacking from $\boldsymbol {r}_{1}^{\text {t}, \text {SC-CP}}$ to $\boldsymbol {r}_{M}^{\text {t}, \text {SC-CP}}$ in (4), and multiplying a unitary matrix of

(5)\begin{align} \begin{bmatrix} \boldsymbol {D} & \textbf{0} & \dotsb & \textbf{0} \\ \textbf{0} & \boldsymbol {D} & & \vdots \\ \vdots & & \ddots & \textbf{0} \\ \textbf{0} & \dotsb & \textbf{0} & \boldsymbol {D} \end{bmatrix} \in \mathbb{C}^{QM \times QM} \end{align}

from the left of both sides, we have the frequency domain received non-precoded SC-CP IoT signal vector at the base station as

(6)\begin{align} \begin{bmatrix} \boldsymbol {r}_{1}^{\text{f}, \text{SC-CP}} \\ \vdots \\ \boldsymbol {r}_{M}^{\text{f}, \text{SC-CP}} \end{bmatrix} &= \begin{bmatrix} \boldsymbol {D} & \textbf{0} & \dotsb & \textbf{0} \\ \textbf{0} & \boldsymbol {D} & & \vdots \\ \vdots & & \ddots & \textbf{0} \\ \textbf{0} & \dotsb & \textbf{0} & \boldsymbol {D} \end{bmatrix} \begin{bmatrix} \boldsymbol {r}_{1}^{\text{t}, \text{SC-CP}} \\ \vdots \\ \boldsymbol {r}_{M}^{\text{t}, \text{SC-CP}} \end{bmatrix} \end{align}

(7)\begin{align} &= \begin{bmatrix} \boldsymbol {\Lambda}^{(1,1)} \boldsymbol {D} & \hspace{-2mm} \dotsb & \hspace{-2mm} \boldsymbol {\Lambda}^{(1,N)} \boldsymbol {D} \\ \vdots & \hspace{-2mm} & \hspace{-2mm} \vdots \\ \boldsymbol {\Lambda}^{(M,1)} \boldsymbol {D} & \hspace{-2mm} \dotsb & \hspace{-2mm} \boldsymbol {\Lambda}^{(M,N)} \boldsymbol {D} \end{bmatrix} \begin{bmatrix} \boldsymbol {s}_{1} \\ \vdots \\ \boldsymbol {s}_{N} \end{bmatrix} + \begin{bmatrix} \boldsymbol {v}_{1}^{\text{f}} \\ \vdots \\ \boldsymbol {v}_{M}^{\text{f}} \end{bmatrix}, \end{align}

where $\boldsymbol {r}_{m}^{\text {f}, \text {SC-CP}} \in \mathbb {C}^{Q}$ is the frequency domain received SC-CP signal vector at the $m$-th base station antenna and $\boldsymbol {v}_{m}^{\text {f}}=\boldsymbol {D}\boldsymbol {v}_{m}^{\text {t}}$ ($m=1,\dotsc , M$). It should be noted here that this received signal model can be regarded as a special case of (1), where the precoding matrix $\boldsymbol {P}$ is set to be $\boldsymbol {D}$. Thus, if DFT matrix $\boldsymbol {D}$ is appropriate for the precoding matrix of overloaded MU-MIMO OFDM system with the convex optimization-based signal detection, then the choice of non-precoded SC-CP signaling is extremely suited for IoT environments because this approach requires neither the IDFT operation nor the precoding operation at the IoT node (transmitter side).

A) Discreteness and group sparsity of transmitted symbol vector

One of the main ideas of the proposed method is based on the fact that the transmitted symbol has the discreteness, i.e. the element of $\boldsymbol {s}$ is in the finite set $\mathcal {S} \cup \{ 0 \}$. For the reconstruction of such discrete-valued vector, several methods have been proposed [Reference Tan, Rasmussen and Lim21–Reference Nagahara24]. These methods reconstruct the discrete-valued vector in the real domain. For the reconstruction of the complex discrete-valued vector, SCSR optimization [Reference Hayakawa and Hayashi13] has been proposed by extending the approach in [Reference Nagahara24]. For details of the related methods, see Section III-F.

When only a few transmit terminals are active, the transmitted symbol vector $\boldsymbol {s}$ has not only the discreteness but also the group sparsity, i.e. $\boldsymbol {s}_{n}=\textbf{0}_{Q}$ holds for most transmit terminals, because the transmitted symbol vector of non-active terminals can be regarded as $\textbf{0}_{Q}$. Note that not only is $\boldsymbol {s}$ sparse, but multiple transmitted symbols from the non-active IoT terminal become zero simultaneously when OFDM or SC-CP signaling is employed. Such group sparsity has been utilized as prior knowledge in the literature of compressed sensing and sparse regression [Reference Yuan and Lin25–Reference Stojnic, Parvaresh and Hassibi27] as well as wireless communications [Reference Messai, Amis, Guilloud and Aïssa-El-Bey28, Reference Messai, Aïssa-El-Bey, Amis and Guilloud29]. Thus, we can expect a certain improvement of the detection performance by using the group sparsity of the transmitted signal vector $\boldsymbol {s}$.

B) Proposed DGS optimization

The proposed method utilizes both the discreteness and the group sparsity of $\boldsymbol {s}$ discussed in the previous subsection. The proposed DGS optimization problem is given by

(12)\begin{align} &\mathop{\mathrm{minimize}}\limits_{\boldsymbol {x} \in \mathbb{C}^{QN}} \left\{{ \sum_{\ell=1}^{S} q_{\ell} g_{\ell} \left({\boldsymbol {x} - c_{\ell}\textbf{1}}\right) + \alpha \sum_{n=1}^{N} \left\lVert{\boldsymbol {x}_{n}}\right\rVert_{2} }\right\} \nonumber\\ &\ \mathrm{subject\ to\ } \left\lVert{ \boldsymbol {r}-\boldsymbol {A}\boldsymbol {x} }\right\rVert_{2} \le \varepsilon, \end{align}

where $c_{1},\dotsc ,c_{S}$ denotes all elements in $\{0\} \cup \mathcal {S}$. For example, when we use quadrature phase shift keying (QPSK) and set $\mathcal {S}=\left \{{1+j, -1+j, -1-j, 1-j}\right \}$, we have $S=5$ and $c_{1}=0, c_{2}=1+j, c_{3}=-1+j, c_{4}=-1-j, c_{5}=1-j$. The vector $\boldsymbol {x}_{n} \in \mathbb {C}^{Q}$ is the $n$-th subvector of $\boldsymbol {x}=\left [{\boldsymbol {x}_{1}^{\mathsf {T}}\ \dotsb \ \boldsymbol {x}_{N}^{\mathsf {T}}}\right ]^{\mathsf {T}} \in \mathbb {C}^{QN}$. $q_{\ell }$ ($>0$), $\alpha$ ($>0$), and $\varepsilon$ ($>0$) are parameters ($\ell =1,\dotsc ,S$), where $\sum _{\ell =1}^{S}q_{\ell }=1$. The function $g_{\ell }(\cdot )$ is a sparse regularizer for a sparse vector in the complex-valued domain. In this paper, we consider two convex regularizers given by

(13)\begin{align} h_{1}(\boldsymbol {u}) &= \lVert \boldsymbol {u} \rVert_{1} \end{align}

(14)\begin{align} &= \sum_{i=1}^{QN} \sqrt{\mathrm{Re}\{{u_{i}}\}^{2} + \mathrm{Im}\{{u_{i}}\}^{2}}, \end{align}

(15)\begin{align} h_{2}(\boldsymbol {u}) &= \left\lVert \mathrm{Re}\{{\boldsymbol {u}}\} \right\rVert_{1} + \left\lVert \mathrm{Im}\{{\boldsymbol {u}}\} \right\rVert_{1} \end{align}

(16)\begin{align} &= \sum_{i=1}^{QN} \left( \left\lvert \mathrm{Re}\{{u_{i}}\} \right\rvert + \left\lvert \mathrm{Im}\{{u_{i}}\} \right\rvert \right) \end{align}

in accordance with [Reference Hayakawa and Hayashi13]. The first one $h_{1}(\cdot )$ is based on the modulus of the complex number, whereas the second one $h_{2}(\cdot )$ considers the real part and the imaginary part separately. We can also use some non-convex sparse regularizers as in [Reference Hayakawa and Hayashi30], though the global optimizer cannot be necessarily obtained.

The objective function of the DGS optimization (12) is the weighted sum of two regularizers. The first term $\sum _{\ell =1}^{S} q_{\ell } g_{\ell } \left ({\boldsymbol {x} - c_{\ell }\textbf{1}}\right )$ in (12) has been proposed in [Reference Hayakawa and Hayashi13] as a regularizer for the discrete-valued vector. On the basis of the idea of compressed sensing [Reference Donoho8, Reference Hayashi, Nagahara and Tanaka9], the sparse regularizer $g_{\ell }(\cdot )$ is applied by using the fact that $\boldsymbol {s}-c_{\ell }\textbf{1}$ has some zero elements. A good choice of the sparse regularizer depends on the symbol alphabet $\mathcal {S}$. For example, when we use QPSK and set $c_{1}=0, c_{2}=1+j, c_{3}=-1+j, c_{4}=-1-j, c_{5}=1-j$, it is reasonable to define $g_{1}(\cdot )=h_{1}(\cdot )$ and $g_{\ell }(\cdot )=h_{2}(\cdot )$ ($\ell =2,\dotsc ,5$). For more detailed discussion, see [Reference Hayakawa and Hayashi13]. In addition to the first regularizer, we also use the regularizer $\alpha \sum _{n=1}^{N} \left \lVert {\boldsymbol {x}_{n}}\right \rVert _{2}$ to promote the group sparsity of $\boldsymbol {s}$. Note that we use the $\ell _{2}$ norm $\left \lVert {\boldsymbol {x}_{n}}\right \rVert _{2}=\sqrt {\sum _{i=1}^{Q} \left \lvert {x_{n,i}}\right \rvert ^{2}}$ of the complex-valued vector $\boldsymbol {x}_{n} = \left [{x_{n,1}\ \dotsb \ x_{n,Q}}\right ]^{\mathsf {T}} \in \mathbb {C}^{Q}$, whereas the $\ell _{2}$ norm of the real valued vector is commonly used in the literature [Reference Yuan and Lin25–Reference Stojnic, Parvaresh and Hassibi27].

It should be noted that the conventional method [Reference Hayashi, Nakai, Hayakawa and Ha10] considers an optimization problem in the real domain, whereas the DGS optimization problem (12) is the optimization in the complex domain $\mathbb {C}^{QN}$. Such approach enables us to directly utilize the discrete nature of the transmitted symbols in the complex domain. Moreover, the DGS optimization can utilize the group sparsity of the transmitted symbol vector, which has not been used in the conventional methods.

C) Proposed algorithm for DGS optimization based on ADMM

We here propose an algorithm for the DGS optimization. Since the objective function of the DGS optimization problem (12) is not differentiable, simple gradient-based methods cannot be applied. However, we can derive an optimization algorithm for the DGS optimization (12) on the basis of ADMM [Reference Gabay and Mercier15–Reference Li, Wang and Wang19] by appropriate reformulation as shown below. We firstly rewrite the DGS optimization problem (12) with new variables $\boldsymbol {z}_{1}, \dotsc , \boldsymbol {z}_{S}, \boldsymbol {z}_{\text {GS}} \in \mathbb {C}^{QN}$ and $\boldsymbol {z}_{\text {B}} \in \mathbb {C}^{QM}$ as

(17)\begin{align} \mathop{\mathrm{minimize}}\limits_{\substack{\boldsymbol {x}, \boldsymbol {z}_{1}, \dotsc, \boldsymbol {z}_{S}, \boldsymbol {z}_{\text{GS}} \in \mathbb{C}^{QN}\\ \boldsymbol {z}_{\text{B}} \in \mathbb{C}^{QM}}}&\quad \left\{ \sum_{\ell=1}^{S} q_{\ell} g_{\ell}\left({\boldsymbol {z}_{\ell} - c_{\ell}\textbf{1}}\right)\right.\notag\\ &\qquad +\left. \alpha \sum_{n=1}^{N} \left\lVert{\boldsymbol {z}_{\text{GS},n}}\right\rVert_{2} + \chi_{\mathcal{B}} \left({\boldsymbol {z}_{\text{B}}}\right) \right\} \nonumber\\ \mathrm{subject\ to}&\quad \boldsymbol {x} = \boldsymbol {z}_{1} = \dotsb = \boldsymbol {z}_{L} = \boldsymbol {z}_{\text{GS}}, \nonumber\\ &\quad \boldsymbol {A} \boldsymbol {x} = \boldsymbol {z}_{\text{B}}. \end{align}

In (17), we express the constraint $\left \lVert { \boldsymbol {r}-\boldsymbol {A}\boldsymbol {x} }\right \rVert _{2} \le \varepsilon$ in (12) by $\chi _{\mathcal {B}} \left ({\boldsymbol {z}_{\text {B}}}\right ) = \chi _{\mathcal {B}} \left ({\boldsymbol {A} \boldsymbol {x}}\right )$ in the objective function, where we define $\mathcal {B} := \left \{{ \boldsymbol {u} \in \mathbb {C}^{QM} \,\left |\, \left \lVert {\boldsymbol {r}-\boldsymbol {u}}\right \rVert _{2} \le \varepsilon \right . }\right \}$ and the corresponding indicator function

(18)\begin{align} \chi_{\mathcal{B}}(\boldsymbol {z}_{\text{B}}) &= \begin{cases} 0 & (\boldsymbol {z}_{\text{B}} \in \mathcal{B}) \\ \infty & (\boldsymbol {z}_{\text{B}} \notin \mathcal{B}) \end{cases}. \end{align}

Since the objective function in (17) becomes infinity when $\boldsymbol {x}$ does not satisfy the constraint $\left \lVert { \boldsymbol {r}-\boldsymbol {A}\boldsymbol {x} }\right \rVert _{2} \le \varepsilon$, the optimization problem (17) is equivalent to the original optimization problem (12). The vector $\boldsymbol {z}_{\text {GS},n} \in \mathbb {C}^{Q}$ is the $n$-th subvector of $\boldsymbol {z}_{\text {GS}}=[{\boldsymbol {z}_{\text {GS},1}^{\mathsf {T}}\ \dotsb \ \boldsymbol {z}_{\text {GS},N}^{\mathsf {T}}}]^{\mathsf {T}} \in \mathbb {C}^{QN}$. To further rewrite the optimization problem (17), we denote the objective function as

(19)\begin{align} \hspace{-2mm} g_{\text{GS}}(\boldsymbol {z}_{\text{GS}}) &= \sum_{n=1}^{N} \left\lVert{\boldsymbol {z}_{\text{GS},n}}\right\rVert_{2}, \end{align}

(20)\begin{align} \hspace{-2mm} f(\boldsymbol {z}) &= \sum_{\ell=1}^{S} q_{\ell} g_{\ell}\left({\boldsymbol {z}_{\ell} - c_{\ell}\textbf{1}}\right) + \alpha g_{\text{GS}}(\boldsymbol {z}_{\text{GS}}) + \chi_{\mathcal{B}} \left({\boldsymbol {z}_{\text{B}}}\right). \end{align}

The constraint in (17) can also be summarized as $\boldsymbol {\Phi } \boldsymbol {x} = \boldsymbol {z}$, where we define

(21)\begin{align} \boldsymbol {z} &= \left[{\boldsymbol {z}_{1}^{\mathsf{T}}\ \dotsb\ \boldsymbol {z}_{S}^{\mathsf{T}}\ \boldsymbol {z}_{\text{GS}}^{\mathsf{T}}\ \boldsymbol {z}_{\text{B}}^{\mathsf{T}}}\right]^{\mathsf{T}} \nonumber\\ &\quad \in \mathbb{C}^{(S+1)QN+QM}\end{align}

(22)\begin{align} \boldsymbol {\Phi} &= \left[{\boldsymbol {I}_{QN}\ \dotsb\ \boldsymbol {I}_{QN}\ \boldsymbol {I}_{QN}\ \boldsymbol {A}^{\mathsf{T}}}\right]^{\mathsf{T}} \nonumber\\ &\quad \in \mathbb{C}^{((S+1)QN+QM) \times QN}. \end{align}

By using the above expression, the DGS optimization problem (17) can be finally represented as

(23)\begin{align} \mathop{\mathrm{minimize}}\limits_{\substack{\boldsymbol {x} \in \mathbb{C}^{QN}\\ \boldsymbol {z} \in \mathbb{C}^{(S+1)QN+QM}}} &\quad f(\boldsymbol {z}) \nonumber\\ \mathrm{subject\ to} & \quad \boldsymbol {\Phi} \boldsymbol {x} = \boldsymbol {z}. \end{align}

The optimization problem (23) is a standard form for ADMM, and hence we can derive the optimization algorithm based on ADMM.

The update equations of ADMM for (23) are given by

(24)\begin{align} \boldsymbol {x}^{k+1} &= \mathop{{\mathrm{arg\ min}}}\limits_{\boldsymbol {x} \in \mathbb{C}^{QN}} \left\{{\rho \left\lVert{\boldsymbol {\Phi}\boldsymbol {x}-\boldsymbol {z}^{k}+\boldsymbol {w}^{k}}\right\rVert_{2}^{2}}\right\}, \end{align}

(25)\begin{align} \boldsymbol {y}^{k+1} &= \boldsymbol {\Phi}\boldsymbol {x}^{k+1}+\boldsymbol {w}^{k}, \end{align}

(26)\begin{align} \boldsymbol {z}^{k+1} &= \mathop{\mathrm{arg\ min}}\limits_{\boldsymbol {z} \in \mathbb{C}^{(S+1)QN+QM}} \left\{{f(\boldsymbol {z}) + \rho \left\lVert{\boldsymbol {y}^{k+1}-\boldsymbol {z}}\right\rVert_{2}^{2}}\right\}, \end{align}

(27)\begin{align} \boldsymbol {w}^{k+1} &= \boldsymbol {y}^{k+1} - \boldsymbol {z}^{k+1}, \end{align}

where $\rho$ ($>0$) is a parameter and $\boldsymbol {w}^{k} \in \mathbb {C}^{(S+1)QN+QM}$. In the algorithm, $\boldsymbol {x}^{k}$ is the estimate for the transmitted signal vector $\boldsymbol {s}$ at the $k$-th iteration. The Wirtinger derivative [Reference Brandwood31] of $\left \lVert {\boldsymbol {\Phi }\boldsymbol {x}-\boldsymbol {z}^{k}+\boldsymbol {w}^{k}}\right \rVert _{2}^{2}$ in (24) with respect to $\boldsymbol {x}$ is given by

(28)\begin{align} &\frac{\partial}{\partial \boldsymbol {x}^{\mathsf{H}}} \left\{{\left\lVert{\boldsymbol {\Phi}\boldsymbol {x}-\boldsymbol {z}^{k}+\boldsymbol {w}^{k}}\right\rVert_{2}^{2}}\right\} \nonumber\\ &\quad = \left({(S+1)\boldsymbol {I}_{QN}+\boldsymbol {A}^{\mathsf{H}}\boldsymbol {A}}\right) \boldsymbol {x} \nonumber\\ & \qquad - \left(\sum_{\ell=1}^{S} ({\boldsymbol {z}_{\ell}^{k} - \boldsymbol {w}_{\ell}^{k}}) + ({\boldsymbol {z}_{\text{GS}}^{k}-\boldsymbol {w}_{\text{GS}}^{k}}) + \boldsymbol {A}^{\mathsf{H}} ({\boldsymbol {z}_{\text{B}}^{k}-\boldsymbol {w}_{\text{B}}^{k}})\right). \end{align}

Here, in the same manner as $\boldsymbol {z}$ in (21), we define $\boldsymbol {w}_{1}^{k}, \dotsc , \boldsymbol {w}_{S}^{k}, \boldsymbol {w}_{\text {GS}}^{k} \in \mathbb {C}^{QN}$ and $\boldsymbol {w}_{\text {B}}^{k} \in \mathbb {C}^{QM}$ as the subvectors of the vector

(29)\begin{align} \boldsymbol {w}^{k} &= \left[{{\boldsymbol {w}_{1}^{k}}^{\mathsf{T}}\ \dotsb\ {\boldsymbol {w}_{S}^{k}}^{\mathsf{T}}\ {\boldsymbol {w}_{\text{GS}}^{k}}^{\mathsf{T}}\ {\boldsymbol {w}_{\text{B}}^{k}}^{\mathsf{T}}}\right]^{\mathsf{T}} \nonumber\\ &\in \mathbb{C}^{(S+1)QN+QM}. \end{align}

From (28), $({\partial }/{\partial \boldsymbol {x}^{\mathsf {H}})} \{\left \lVert {\boldsymbol {\Phi }\boldsymbol {x}-\boldsymbol {z}^{k}+\boldsymbol {w}^{k}}\right \rVert _{2}^{2}\} = \textbf{0}$ gives the explicit formula for the update of $\boldsymbol {x}^{k}$ in (24) as

(30)\begin{align} \boldsymbol {x}^{k+1} &= \left({(S+1) \boldsymbol {I}_{QN} + \boldsymbol {A}^{\mathsf{H}}\boldsymbol {A}}\right)^{-1} \Biggl( \sum_{\ell=1}^{S} \left({\boldsymbol {z}_{\ell}^{k} - \boldsymbol {w}_{\ell}^{k}}\right) \nonumber\\ & \quad + \left({\boldsymbol {z}_{\text{GS}}^{k}-\boldsymbol {w}_{\text{GS}}^{k}}\right) + \boldsymbol {A}^{\mathsf{H}} \left({\boldsymbol {z}_{\text{B}}^{k}-\boldsymbol {w}_{\text{B}}^{k}}\right) \Biggr). \end{align}

The update of $\boldsymbol {y}^{k}$ in (25) can be expressed separately as

(31)\begin{align} \boldsymbol {y}^{k+1} &= \begin{bmatrix} \boldsymbol {y}_{1}^{k} \\ \vdots \\ \boldsymbol {y}_{S}^{k} \\ \boldsymbol {y}_{\text{GS}}^{k} \\ \boldsymbol {y}_{\text{B}}^{k} \end{bmatrix} = \begin{bmatrix} \boldsymbol {x} + \boldsymbol {w}_{1}^{k} \\ \vdots \\ \boldsymbol {x} + \boldsymbol {w}_{S}^{k} \\ \boldsymbol {x} + \boldsymbol {w}_{\text{GS}}^{k} \\ \boldsymbol {A} \boldsymbol {x} + \boldsymbol {w}_{\text{B}}^{k} \end{bmatrix}. \end{align}

Next, we consider the update of $\boldsymbol {z}^{k}$ in (26), which can be denoted as

(32)\begin{align} \boldsymbol {z}^{k+1} &= \mathrm{prox}_{({1}/{2\rho}) f} \left({\boldsymbol {y}^{k+1}}\right) \end{align}

(33)\begin{align} &= \begin{bmatrix} c_{1}\textbf{1} + \mathrm{prox}_{({q_{1}}/{2\rho}) g_{1}} \left({\boldsymbol {y}_{1}^{k+1}-c_{1}\textbf{1}}\right) \\ \vdots \\ c_{S}\textbf{1} + \mathrm{prox}_{({q_{S}}/{2\rho}) g_{S}} \left({\boldsymbol {y}_{L}^{k+1}-c_{L}\textbf{1}}\right) \\ \mathrm{prox}_{({\alpha}/{2\rho}) g_{\text{GS}}} \left({\boldsymbol {y}_{\text{GS}}^{k+1}}\right) \\ \mathrm{prox}_{({1}/{2\rho}) \chi_{\mathcal{B}}} \left({\boldsymbol {y}_{\text{B}}^{k+1}}\right) \end{bmatrix} \end{align}

because the function $f(\cdot )$ is separable as in (20), where we have also used the property of proximity operator for translation [Reference Combettes and Pesquet17]. As discussed in [Reference Hayakawa and Hayashi13], the proximity operator of $\gamma g_{\ell }(\cdot )$ can be computed by using

(34)\begin{align} \left[{\mathrm{prox}_{\gamma h_{1}} \left({\boldsymbol {u}}\right)}\right]_{i} &= \begin{cases} \left({\left\lvert{u_{i}}\right\rvert - \gamma}\right) \dfrac{u_{i}}{\left\lvert{u_{i}}\right\rvert} & (\left\lvert{u_{i}}\right\rvert \ge \gamma) \\ 0 & (\left\lvert{u_{i}}\right\rvert < \gamma)\end{cases},\end{align}

(35)\begin{align} \left[{\mathrm{prox}_{\gamma h_{2}} \left({\boldsymbol {u}}\right)}\right]_{i} &= \mathrm{sign}\left({\mathrm{Re}\{{u_{i}}\}}\right) \max\left({\left\lvert{\mathrm{Re}\{{u_{i}}\}}\right\rvert - \gamma, 0}\right) \nonumber\\ & \quad + j \cdot \mathrm{sign}\left({\mathrm{Im}\{{u_{i}}\}}\right) \max\left({\left\lvert{\mathrm{Im}\{{u_{i}}\}}\right\rvert - \gamma, 0}\right), \end{align}

where $\gamma >0$ and $\boldsymbol {u} = [{u_{1}\ \dotsb \ u_{QN}}]^{\mathsf {T}} \in \mathbb {C}^{QN}$. We denote the $i$-th element of the vector by $[\cdot ]_{i}$. The proximity operator of $\gamma g_{\text {GS}}(\cdot )$ is given by

(36)\begin{align} \hspace{-4mm} \left[{\mathrm{prox}_{\gamma g_{\text{GS}}} (\boldsymbol {u})}\right]_{n} &= \begin{cases} \left({\left\lVert{\boldsymbol {u}_{n}}\right\rVert_{2}-\gamma}\right) \dfrac{\boldsymbol {u}_{n}}{\left\lVert{\boldsymbol {u}_{n}}\right\rVert_{2}} & \hspace{-2mm} (\left\lVert{\boldsymbol {u}_{n}}\right\rVert_{2} \ge \gamma) \\ \textbf{0}_{Q} & \hspace{-2mm} (\left\lVert{\boldsymbol {u}_{n}}\right\rVert_{2} < \gamma) \end{cases}, \end{align}

where $\boldsymbol {u}_{n} \in \mathbb {C}^{Q}$ and $[{{\rm prox}_{\gamma g_{\text {GS}}} (\boldsymbol {u})}]_{n} \in \mathbb {C}^{Q}$ are the $n$-th subvectors of $\boldsymbol {u}=\left [{\boldsymbol {u}_{1}^{\mathsf {T}}\ \dotsb \ \boldsymbol {u}_{N}^{\mathsf {T}}}\right ]^{\mathsf {T}} \in \mathbb {C}^{ QN}$ and ${\rm prox}_{\gamma g_{\text {GS}}} (\boldsymbol {u})$, respectively. The proximity operator of $\gamma \chi _{\mathcal {B}}(\cdot )$ can be written as the projection onto $\mathcal {B}$, i.e.

(37)\begin{align} \mathrm{prox}_{\gamma \chi_{\mathcal{B}}} (\boldsymbol {u}) &= \begin{cases} \boldsymbol {r} + \varepsilon \dfrac{\boldsymbol {u}-\boldsymbol {r}}{\left\lVert{\boldsymbol {u}-\boldsymbol {r}}\right\rVert_{2}} & (\left\lVert{\boldsymbol {u} - \boldsymbol {r}}\right\rVert_{2} \ge \varepsilon) \\ \boldsymbol {u} & (\left\lVert{\boldsymbol {u} - \boldsymbol {r}}\right\rVert_{2} < \varepsilon) \end{cases}. \end{align}

We summarize the proposed algorithm for the DGS optimization (23) in Algorithm 1.

Algorithm 1 Proposed Algorithm for DGS Optimization (23)

D) IW-DGS

By using the weighting technique in [Reference Hayakawa and Hayashi5, Reference Hayakawa and Hayashi13, Reference Candès, Wakin and Boyd32], we here extend the DGS optimization into the W-DGS optimization so that we can use prior information. Moreover, we propose an iterative approach named IW-DGS, where we iteratively solve the W-DGS optimization with the updated parameters.

Assuming that the sparse regularizer $g_{\ell }(\cdot )$ is separable as $h_{1}(\cdot )$ or $h_{2}(\cdot )$, we extend the DGS optimization (12) to the W-DGS optimization given by

(38)\begin{align} &\mathop{\mathrm{minimize}}\limits_{\boldsymbol {x} \in \mathbb{C}^{QN}} \left\{{ \sum_{\ell=1}^{S} \sum_{i=1}^{QN} q_{i,\ell} g_{\ell} \left({x_{i} - c_{\ell}}\right) + \sum_{n=1}^{N} \alpha_{n} \left\lVert{\boldsymbol {x}_{n}}\right\rVert_{2} }\right\} \nonumber\\ &\ \mathrm{subject\ to\ } \left\lVert{\boldsymbol {r}-\boldsymbol {A}\boldsymbol {x}}\right\rVert_{2} \le \varepsilon, \end{align}

where $q_{i,\ell },\alpha _{n}$ ($>0$) are the parameters. When $q_{i,\ell }=q_{\ell }$ ($i=1,\dotsc ,QN$) and $\alpha _{n}=\alpha$ ($n=1,\dotsc ,N$), the W-DGS optimization (38) is equivalent to the DGS optimization (12). In the W-DGS optimization (38), we can use the different coefficients $q_{i,\ell }$ and $\alpha _{n}$ for each transmitted symbol $s_{i}$ and transmitted subvector $\boldsymbol {s}_{n}$, respectively. We can thus tune these parameters if we have some prior information about $\boldsymbol {s}$, e.g. a tentative estimate of $\boldsymbol {s}$.

We provide the ADMM-based algorithm for the W-DGS optimization (38). Letting $\tilde {g}_{\text {GS},n}(\boldsymbol {z}_{\text {GS}})=\lVert {\boldsymbol {z}_{\text {GS},n}}\rVert _{2}$ and $\tilde {f}(\boldsymbol {z})=\sum _{\ell =1}^{S} \sum _{i=1}^{QN} q_{i,\ell } g_{\ell } ({z_{i,\ell } - c_{\ell }}) + \sum _{n=1}^{N} \alpha _{n} \tilde {g}_{\text {GS},n}(\boldsymbol {z}_{\text {GS}}) + \chi _{\mathcal {B}} ({\boldsymbol {z}_{\text {B}}})$, we can rewrite the W-DGS optimization problem (38) as

(39)\begin{align} \mathop{\mathrm{minimize}}\limits_{\substack{\boldsymbol {x} \in \mathbb{C}^{QN}\\ \boldsymbol {z} \in \mathbb{C}^{(S+1)QN+QM}}} &\quad \tilde{f}(\boldsymbol {z}) \nonumber\\ \mathrm{subject\ to} &\quad \boldsymbol {\Phi}\boldsymbol {x}=\boldsymbol {z}, \end{align}

where $\boldsymbol {z}=\left [{\boldsymbol {z}_{1}^{\mathsf {T}}\ \dotsb \ \boldsymbol {z}_{S}^{\mathsf {T}}\ \boldsymbol {z}_{\text {GS}}^{\mathsf {T}}\ \boldsymbol {z}_{\text {B}}^{\mathsf {T}}}\right ]^{\mathsf {T}} \in \mathbb {C}^{(S+1)QN+QM}$, $\boldsymbol {z}_{\ell }=\break \left [{z_{1,\ell }\ \dotsb \ z_{QN,\ell }}\right ]^{\mathsf {T}} \in \mathbb {C}^{QN}$ ($\ell =1,\dotsc ,S$), and $\boldsymbol {z}_{\text {GS}}=[\boldsymbol {z}_{\text {GS},1}^{\mathsf {T}}\ \dotsb \break \boldsymbol {z}_{\text {GS},N}^{\mathsf {T}}]^{\mathsf {T}} \in \mathbb {C}^{QN}$. The ADMM-based algorithm for (39) can be obtained by replacing the update of $\boldsymbol {z}_{\ell }^{k}$ and $\boldsymbol {z}_{\text {GS},n}^{k}$ in Algorithm 1 with

(40)\begin{align} \boldsymbol {z}_{\ell}^{k+1} &= \begin{bmatrix} c_{\ell} + \mathrm{prox}_{\frac{q_{1,\ell}}{2\rho} g_{\ell}} \left({y_{1,\ell}^{k+1}-c_{\ell}}\right) \\ \vdots \\ c_{\ell} + \mathrm{prox}_{\frac{q_{CN,\ell}}{2\rho} g_{\ell}} \left({y_{CN,\ell}^{k+1}-c_{\ell}}\right) \end{bmatrix}, \end{align}

(41)\begin{align} \boldsymbol {z}_{\text{GS},n}^{k+1} &= \mathrm{prox}_{\frac{\alpha_{n}}{2\rho}\tilde{g}_{\text{GS},n}} \left({\boldsymbol {y}_{\text{GS},n}^{k+1}}\right) \end{align}

(42)\begin{align} &= \begin{cases} \left({\left\lVert{\boldsymbol {y}_{\text{GS},n}^{k+1}}\right\rVert_{2}-\frac{\alpha_{n}}{2\rho}}\right) \frac{\boldsymbol {y}_{\text{GS},n}^{k+1}}{\left\lVert{\boldsymbol {y}_{\text{GS},n}^{k+1}}\right\rVert_{2}} & \left({\left\lVert{\boldsymbol {y}_{\text{GS},n}^{k+1}}\right\rVert_{2}\ge\frac{\alpha_{n}}{2\rho}}\right) \\ \textbf{0}_{Q} & \left({\left\lVert{\boldsymbol {y}_{\text{GS},n}^{k+1}}\right\rVert_{2}<\frac{\alpha_{n}}{2\rho}}\right) \end{cases}, \end{align}

respectively, where the vector $\boldsymbol {y}_{\text {GS},n}^{k+1} \in \mathbb {C}^{Q}$ is the $n$-th subvector of $\boldsymbol {y}_{\text {GS}}^{k+1}=[{{\boldsymbol {y}_{\text {GS},1}^{k+1}}^{\mathsf {T}}\ \dotsb \ {\boldsymbol {y}_{\text {GS},N}^{k+1}}]^{\mathsf {T}}}^{\mathsf {T}}$.

We propose an iterative algorithm named IW-DGS in Algorithm 2, where we iteratively compute the solution of the W-DGS optimization with the update of the parameters $q_{i,\ell }$ and $\alpha _{n}$.

Algorithm 2 IW-DGS

In Fig. 2, we show the illustration of the proposed IW-DGS. At each outer iteration $t$, we update the parameters by using the tentative estimate at the previous iteration $t-1$. In this paper, we consider the update given by

(43)\begin{align} q_{i,\ell} &= \frac{\left\lvert{\hat{s}_{i}^{\text{pre}}-c_{\ell}}\right\rvert^{-1}}{\sum_{\ell^{\prime}=1}^{S} \left\lvert{\hat{s}_{i}^{\text{pre}}-c_{\ell^{\prime}}}\right\rvert^{-1}}, \end{align}

(44)\begin{align} \alpha_{n} &= \frac{\left\lVert{\hat{\boldsymbol {s}}_{n}^{\text{pre}}}\right\rVert_{2}^{-1}}{\sum_{n^{\prime}=1}^{N}\left\lVert{\hat{\boldsymbol {s}}_{n^{\prime}}^{\text{pre}}}\right\rVert_{2}^{-1}} N\alpha, \end{align}

where we define $\hat {\boldsymbol {s}}^{\text {pre}}=\left [{(\hat {\boldsymbol {s}}_{1}^{\text {pre}})^{\mathsf {T}}\ \dotsb \ (\hat {\boldsymbol {s}}_{N}^{\text {pre}})^{\mathsf {T}}}\right ]^{\mathsf {T}} =[\hat {s}_{1}^{\text {pre}}\ \dotsb \break \hat {s}_{QN}^{\text {pre}}]^{\mathsf {T}} \in \mathbb {C}^{QN}$ as the estimate of $\boldsymbol {s}$ at the previous iteration. The denominators in (43) and (44) play the role of normalization to satisfy $\sum _{\ell =1}^{S} q_{i,\ell } = 1$ for any $i=1,\dotsc ,QN$ and $\sum _{n=1}^{N} \alpha _{n}=N\alpha$, respectively. By using (43) and (44), we can utilize the tentative estimate $\hat {\boldsymbol {s}}^{\text {pre}}$ as the prior information. For example, the parameter $q_{i,\ell }$ becomes large when the corresponding tentative estimate $\hat {s}_{i}^{\text {pre}}$ is close to $c_{\ell }$. In this case, the new estimate of $s_{i}$ also tends to be close to $c_{\ell }$. Similarly, the parameter $\alpha _{n}$ becomes large when the norm of corresponding tentative estimate $\hat {\boldsymbol {s}}_{n}^{\text {pre}}$ is small. In this case, the new estimate of $\boldsymbol {s}_{n}$ also tends to become almost $\textbf{0}_{Q}$.

Fig. 2. Illustration of IW-DGS.

E) Computational complexity

We here discuss the computational complexity of IW-DGS in Algorithm 2. The order of the complexity is dominated by the computation of the inverse matrix $({(S+1) \boldsymbol {I}_{QN} + \boldsymbol {A}^{\mathsf {H}}\boldsymbol {A}})^{-1}$, which requires $\mathcal {O}(Q^{3}N^{3})$ complexity in a direct calculation [Reference Boyd and Vandenberghe33, Ch. 11]. However, by using the structure of $\boldsymbol {A}$, the order of the complexity can be reduced.

We firstly consider the computational complexity for $({(S+1) \boldsymbol {I}_{QN} + \boldsymbol {A}^{\mathsf {H}}\boldsymbol {A}})^{-1}$ in the case of MU-MIMO OFDM systems (11) without precoding, i.e. $\boldsymbol {P}=\boldsymbol {I}_{Q}$. In this case, $\boldsymbol {\Theta } := (S+1) \boldsymbol {I}_{QN} + \boldsymbol {A}^{\mathsf {H}}\boldsymbol {A} \in \mathbb {C}^{QN \times QN}$ is composed of $N^{2}$ diagonal matrices as

(45)\begin{align} \boldsymbol {\Theta} &= \begin{bmatrix} \boldsymbol {\Theta}^{(1,1)} & \dotsb & \boldsymbol {\Theta}^{(1,N)} \\ \vdots & \ddots & \vdots \\ \boldsymbol {\Theta}^{(N,1)} & \dotsb & \boldsymbol {\Theta}^{(N,N)} \end{bmatrix}, \end{align}

where $\boldsymbol {\Theta }^{(n_{1},n_{2})} = {\rm diag} ({\theta _{1}^{(n_{1},n_{2})}, \dotsc , \theta _{Q}^{(n_{1},n_{2})}}) = \delta _{n_{1}n_{2}}(S+1)\break \boldsymbol {I}_{Q} + \sum _{m=1}^{M} ({\boldsymbol {\Lambda }^{(m,n_{1})}})^{\mathsf {H}} \boldsymbol {\Lambda }^{(m,n_{2})} \in \mathbb {C}^{Q \times Q}$ and $\delta _{n_{1}n_{2}}$ denotes the Dirac delta. Since $\boldsymbol {\Theta }^{(n_{1},n_{2})}$ is a diagonal matrix as well as $\boldsymbol {\Lambda }^{(m,n_{1})}$ and $\boldsymbol {\Lambda }^{(m,n_{2})}$, each $\boldsymbol {\Theta }^{(n_{1},n_{2})}$ can be computed with $\mathcal {O}(QM)$. We define

(46)\begin{align} \tilde{\boldsymbol {\Theta}}_{c} &= \begin{bmatrix} \theta_{c}^{(1,1)} & \dotsb & \theta_{c}^{(1,N)} \\ \vdots & \ddots & \vdots \\ \theta_{c}^{(N,1)} & \dotsb & \theta_{c}^{(N,N)} \end{bmatrix} \in \mathbb{C}^{N \times N} \end{align}

and its inverse matrix

(47)\begin{align} \tilde{\boldsymbol {\Omega}}_{c} &= \begin{bmatrix} \omega_{c}^{(1,1)} & \dotsb & \omega_{c}^{(1,N)} \\ \vdots & \ddots & \vdots \\ \omega_{c}^{(N,1)} & \dotsb & \omega_{c}^{(N,N)} \end{bmatrix} \end{align}

(48)\begin{align} &:= \tilde{\boldsymbol {\Theta}}_{c}^{-1} \in \mathbb{C}^{N \times N}, \end{align}

which correspond to the $c$-th subcarrier. The inverse matrix of $\boldsymbol {\Theta }$ in (45) can be written as

(49)\begin{align} \boldsymbol {\Theta}^{-1} &= \begin{bmatrix} \boldsymbol {\Omega}^{(1,1)} & \dotsb & \boldsymbol {\Omega}^{(1,N)} \\ \vdots & \ddots & \vdots \\ \boldsymbol {\Omega}^{(N,1)} & \dotsb & \boldsymbol {\Omega}^{(N,N)} \end{bmatrix}, \end{align}

where $\boldsymbol {\Omega }^{(n_{1},n_{2})} = {\rm diag} ({\omega _{1}^{(n_{1},n_{2})}, \dotsc , \omega _{Q}^{(n_{1},n_{2})}}) \in \mathbb {C}^{Q \times Q}$. Since each $\tilde {\boldsymbol {\Omega }}_{c} \in \mathbb {C}^{N \times N}$ ($c=1,\dotsc ,Q$) in (48) can be obtained with $\mathcal {O}(N^{3})$, the overall complexity for the inverse matrix $({(S+1) \boldsymbol {I}_{QN} + \boldsymbol {A}^{\mathsf {H}}\boldsymbol {A}})^{-1}$ becomes $\mathcal {O}(QN^{3})$, which is much lower than $\mathcal {O}(Q^{3}N^{3})$ in the direct calculation. It should be noted that we can compute each $\tilde {\boldsymbol {\Omega }}_{c}$ in parallel. Moreover, the computation of the inverse matrix is required only once in the algorithm. Since the matrix-vector multiplication related to $\boldsymbol {A}$ requires $\mathcal {O}(Q^{2}MN)$, the overall complexity of IW-DGS with respect to the system size is $\mathcal {O}(Q^{2}MN+QN^{3})$ in this case. Note that, in practice, we can compute the matrix-vector multiplication more efficiently by using the sparsity of $\boldsymbol {A}$, i.e. the number of non-zero elements in $\boldsymbol {A}$ is only $QMN$, whereas the number of all elements is $Q^{2}MN$.

Even when we use some precoding, we can also reduce the order of the complexity as long as the precoding matrix $\boldsymbol {P}$ is orthogonal, i.e. $\boldsymbol {P}\boldsymbol {P}^{\mathsf {H}}=\boldsymbol {I}_{Q}$. From the Sherman–Morrison–Woodbury formula [Reference Hager34], we have

(50)\begin{align} ({(S+1) \boldsymbol {I}_{QN} + \boldsymbol {A}^{\mathsf{H}}\boldsymbol {A}})^{-1} &= \frac{1}{S+1}\boldsymbol {I}_{QN} - \frac{1}{(S+1)^{2}} \boldsymbol {A}^{\mathsf{H}}\nonumber\\ & \quad \times\left({\boldsymbol {I}_{QM}+\frac{1}{S+1}\boldsymbol {A}\boldsymbol {A}^{\mathsf{H}}}\right)^{-1}\boldsymbol {A}. \end{align}

By using the orthogonality of $\boldsymbol {P}$, the matrix $\boldsymbol {\Xi } := \boldsymbol {I}_{QM}+\frac {1}{S+1}\boldsymbol {A}\boldsymbol {A}^{\mathsf {H}} \in \mathbb {C}^{QM \times QM}$ in (50) can be written with $M^{2}$ diagonal matrices as

(51)\begin{align} \boldsymbol {\Xi} &= \begin{bmatrix} \boldsymbol {\Xi}^{(1,1)} & \dotsb & \boldsymbol {\Xi}^{(1,M)} \\ \vdots & \ddots & \vdots \\ \boldsymbol {\Xi}^{(M,1)} & \dotsb & \boldsymbol {\Xi}^{(M,M)} \end{bmatrix}, \end{align}

where $\boldsymbol {\Xi }^{(m_{1},m_{2})}$ is given by $\boldsymbol {\Xi }^{(m_{1},m_{2})}=\delta _{m_{1}m_{2}}\boldsymbol {I}_{Q} + ({1}/({S+1}))\break \sum _{m=1}^{M} \boldsymbol {\Lambda }^{(m,n_{2})}\boldsymbol {P}({\boldsymbol {\Lambda }^{(m,n_{1})}\boldsymbol {P}})^{\mathsf {H}} =\delta _{m_{1}m_{2}}\boldsymbol {I}_{Q} + ({1}/({S+1}))\sum _{m=1}^{M}\break \boldsymbol {\Lambda }^{(m,n_{2})}({\boldsymbol {\Lambda }^{(m,n_{1})}})^{\mathsf {H}} \in \mathbb {C}^{Q \times Q}$ ($m_{1},m_{2}=1,\dotsc ,M$). Since each $\boldsymbol {\Xi }^{(m_{1},m_{2})}$ is a diagonal matrix, we can calculate the inverse matrix $\boldsymbol {\Xi }^{-1}$ with $\mathcal {O}(QM^{3})$ in the same manner for $\boldsymbol {\Theta }$ in (45). Once we have $\boldsymbol {\Xi }^{-1}=({\boldsymbol {I}_{QM}+({1}/({S+1}))\boldsymbol {A}\boldsymbol {A}^{\mathsf {H}}})^{-1}$, we can update $\boldsymbol {x}^{k}$ only with some matrix-vector multiplications by using (50). Given that the matrix-vector multiplication requires $\mathcal {O}(Q^{2}MN)$, the overall complexity is given by $\mathcal {O}(Q^{2}MN+QM^{3})$ in this case. One drawback of such approaches is that we require more number of matrix-vector multiplications compared to the case where we obtain the inverse matrix $({(S+1) \boldsymbol {I}_{QN} + \boldsymbol {A}^{\mathsf {H}}\boldsymbol {A}})^{-1}$ in advance. However, when the complexity for $({(S+1) \boldsymbol {I}_{QN} + \boldsymbol {A}^{\mathsf {H}}\boldsymbol {A}})^{-1}$ is prohibitive, the approach with the Sherman–Morrison–Woodbury formula (50) would be a promising candidate.

F) Related work

We here discuss some related optimization problems proposed for the discrete-valued vector reconstruction. Table 2 summarizes the conventional optimization problems and the proposed DGS optimization.

Table 2. Convex optimization-based methods for discrete-valued vector reconstruction

The use of the discreteness in the optimization problem has been proposed as the SOAV optimization [Reference Aïssa-El-Bey, Pastor, Sbaï and Fadlallah23, Reference Nagahara24], where the sum of $\ell _{1}$ norms is used as the regularizer for the discrete-valued vector. The idea of the SOAV optimization has been extended to the discrete-valued vector reconstruction in the complex-valued domain [Reference Hayakawa and Hayashi13]. Moreover, the extension to non-convex optimization has been proposed for both real- and complex-valued cases [Reference Hayakawa and Hayashi30]. The non-convex optimization in the real-valued domain is called sum of sparse regularizers (SSR) optimization. Such non-convex optimization-based approaches can achieve better performance than the convex ones in some cases. However, the convergence of the algorithm is not guaranteed for the non-convex optimization, and the performance is sensitive to the choice of the parameters.

The works in [Reference Hayakawa and Hayashi13, Reference Aïssa-El-Bey, Pastor, Sbaï and Fadlallah23, Reference Nagahara24, Reference Hayakawa and Hayashi30] do not focus on the signal detection in MU-MIMO systems, and hence the structure of the measurement matrix is different from the channel matrix considered in this paper. In [Reference Hayakawa and Hayashi30], for example, the ideal i.i.d. Gaussian matrix is used. To evaluate the performance in the MU-MIMO signal detection, the convex SOAV optimization and SCSR optimization have been applied to the scenario in [Reference Hayashi, Nakai, Hayakawa and Ha10, Reference Hayashi, Nakai-Kasai and Hayakawa12], respectively. The results show that the use of the discreteness by the convex optimization problems is effective also in the MU-MIMO signal detection.

The main advantage of the proposed DGS optimization against the above methods is to utilize the group sparsity of the unknown vector, which is preferable for the signal detection in IoT environments considered in this paper. Moreover, unlike [Reference Hayakawa and Hayashi30], the convergence of the proposed algorithm based on ADMM is guaranteed because the DGS optimization is convex.