Robust MIMO Detection With Imperfect CSI: A Neural Network Solution
In this paper, we investigate the design of statistically robust detectors
for multi-input multi-output (MIMO) systems subject to imperfect channel state
information (CSI). A robust maximum likelihood (ML) detection problem is
formulated by taking into consideration the CSI uncertainties caused by both
the channel estimation error and the channel variation. To address the
challenging discrete optimization problem, we propose an efficient alternating
direction method of multipliers (ADMM)-based algorithm, which only requires
calculating closed-form solutions in each iteration. Furthermore, a robust
detection network RADMMNet is constructed by unfolding the ADMM iterations and
employing both model-driven and data-driven philosophies. Moreover, in order to
relieve the computational burden, a low-complexity ADMM-based robust detector
is developed using the Gaussian approximation, and the corresponding deep
unfolding network LCRADMMNet is further established. On the other hand, we also
provide a novel robust data-aided Kalman filter (RDAKF)-based channel tracking
method, which can effectively refine the CSI accuracy and improve the
performance of the proposed robust detectors. Simulation results validate the
significant performance advantages of the proposed robust detection networks
over the non-robust detectors with different CSI acquisition methods.
Yi Sun, Hong Shen, Wei Xu, Nan Hu, Chunming Zhao
# Robust MIMO Detection With Imperfect CSI: A Neural Network Solution

###### Abstract

In this paper, we investigate the design of statistically robust detectors for multi-input multi-output (MIMO) systems subject to imperfect channel state information (CSI). A robust maximum likelihood (ML) detection problem is formulated by taking into consideration the CSI uncertainties caused by both the channel estimation error and the channel variation. To address the challenging discrete optimization problem, we propose an efficient alternating direction method of multipliers (ADMM)-based algorithm, which only requires calculating closed-form solutions in each iteration. Furthermore, a robust detection network RADMMNet is constructed by unfolding the ADMM iterations and employing both model-driven and data-driven philosophies. Moreover, in order to relieve the computational burden, a low-complexity ADMM-based robust detector is developed using the Gaussian approximation, and the corresponding deep unfolding network LCRADMMNet is further established. On the other hand, we also provide a novel robust data-aided Kalman filter (RDAKF)-based channel tracking method, which can effectively refine the CSI accuracy and improve the performance of the proposed robust detectors. Simulation results validate the significant performance advantages of the proposed robust detection networks over the non-robust detectors with different CSI acquisition methods.

 Imperfect channel state information (CSI), multi-input multi-output (MIMO), robust detector, alternating direction method of multipliers (ADMM), deep unfolding

## I Introduction

Efficient multi-input multi-out (MIMO) detection algorithms are essential to unleashing the full potential of MIMO techniques, which have been extensively studied for decades in literature [2, 3, 4]. As is well known, the maximum likelihood (ML) detectors can achieve the optimal performance while with a huge exponential complexity. Linear detectors, including the zero forcing (ZF) and the linear minimum mean squared error (LMMSE) detectors, enjoy lower complexity but suffer from limited performance. To strike a tradeoff between the performance and the complexity, some sub-optimal detection algorithms can be applied, such as sphere decoding (SD) [5], expectation propagation (EP) [6], and approximate message passing (AMP) [7]. The alternating direction method of multipliers (ADMM) [8] has also been widely used to solve the detection problems from the perspective of convex optimization. In general, the detection performance highly hinges on the quality of the available channel state information (CSI), which is inevitably imperfect due to some practical issues. However, these classical MIMO detectors are developed based on the assumption of perfect CSI and therefore non-robust to the CSI uncertainties. That means, when we apply these "mismatched" detectors by directly regarding the acquired imperfect CSI as if it is perfect, the detection performance can be severely degraded [9]. To handle this, the CSI imperfection should be considered during the design of a MIMO detector to enhance its robustness.

In this paper, we focus on the statistically robust MIMO detection with a random CSI error model. Related works include [10, 11, 12, 13], where the robust detectors using the linear MMSE criterion and the widely-linear MMSE criterion were studied for the MIMO systems without and with in-phase or quadrature-phase imbalance (IQI), respectively. Alternatively, a linear minimum error probability (MEP) detector with a given length of pilots was derived to enable ultra-reliable and low latency communication (URLLC) under imperfect CSI [14]. Furthermore, in [15], the authors investigated a robust orthogonal approximate message passing (OAMP) algorithm as well as its corresponding network by exploiting the statistics of imperfect CSI, which can exhibit stronger robustness against CSI errors than the conventional OAMP algorithm [16]. The optimal robust ML detectors were discussed in [17, 18, 19, 20] by incorporating the statistical distribution of channel estimation errors. Although their performance advantages have been validated by both theoretical analysis and simulation results in [17], the exponential complexity impedes their practicality. On the other hand, in most of these works, the CSI errors were characterized by some relatively simple models, e.g., the independent and identically distributed (i.i.d.) Gaussian model, which, however, is inconsistent with practical systems with spatial correlation. Therefore, a more realistic CSI error model needs to be adopted for the robust MIMO detection design.

Recently, the great advancement of deep learning (DL) has motivated the researches on the DL-aided wireless communications [21, 22, 23, 24, 25]. In particular, DL has been successfully applied in the physical layer communication techniques, such
as channel estimation [26, 27], MIMO precoding [28, 29], and, as our main interest, MIMO detection [30, 31, 32, 33, 34, 35], which have been thoroughly reviewed in [24]. The DL-based designs can be categorized into data-driven and model-driven types. The data-driven DL regards the mapping function from the input to the desired output as a black box and directly learns the function by training a network of a specific structure [36], whose performance highly relies on the training dataset. Consequently, the overfitting problem may occur with an insufficient dataset or a large amount of trainable parameters. To relax the requirement on the training data, the model-driven DL that leverages the expert knowledge has emerged as an alternative [22]. One of the most popular model-driven DL approaches is "deep unfolding", which unfolds an existing iterative algorithm into network layers with trainable parameters [23]. Accordingly, the inherent mechanism of the original algorithm can be maintained and only several parameters need to be learned. The deep unfolding based MIMO detection was first investigated in [30], where a network called DetNet was built by mimicking the projected gradient algorithm for the ML optimization problem. Inspired by this, a variety of MIMO detection networks, such as OAMPNet [16], ADMMNet [31], and LcgNet [32], were proposed based on the idea of deep unfolding. However, these networks are designed based on the availability of perfect CSI at the receiver, which can suffer from performance degradation with mismatched CSI.

Motivated by the advantage of the deep unfolding technique, in this paper, we advocate two model-driven robust MIMO detection networks, called RADMMNet and LCRADMMNet, to combat CSI imperfection. Specifically, inspired by the attractive performances of ADMM-based MIMO detectors under perfect CSI [37, 38], we first apply the ADMM framework to derive the solutions to two robust ML detection problems, based on which RADMMNet and LCRADMMNet are further developed. The performance superiorities of the proposed networks over conventional non-robust MIMO detectors are validated via numerical simulations. The main contributions of this paper include:

* Focusing on a spectrum-efficient frame structure with pilots only used in the first block, we model the corresponding CSI imperfection by including both the channel estimation error and the channel variation, based on which a robust ML detection problem is formulated.
* In order to address the complicated robust ML detection problem, we develop an ADMM-based algorithm that only involves the calculations of closed-form expressions in each ADMM iteration, which, to the best of our knowledge, has not been investigated in previous works. Furthermore, by unfolding the ADMM iterations with non-trivial simplifications, a robust detection network, termed as RADMMNet, is established to learn the layer-wise parameters instead of performing an exhaustive search. The network design is further enhanced by transferring intermediate variables and incorporating trainable convolutional layers, so as to enjoy the advantages of both the model-driven deep unfolding technique and the introduced data-driven structure.
* By exploiting a Gaussian approximation for the CSI error, we obtain a simplified reformulation of the robust ML detection problem. Based on this formulation, an ADMM-based robust detector is proposed to approach the performance of our first detector with lower complexity, where closed-form solutions are also derived in each ADMM iteration. The corresponding deep unfolding network, termed as LCRADMMNet is further constructed in a similar way as RADMMNet.
* A robust data-aided Kalman filter-based channel tracking method is established to mitigate the error propagation caused by channel aging, where the pre-estimated data symbols are utilized to improve the accuracy of the channel estimation, thereby enhancing the robust detection performance.

The rest of this paper is organized as follows. In Section II, we describe the system model and formulate the robust ML detection problem. Section III and Section IV elaborate the design of a robust detector and a low-complexity robust detector, respectively, wherein their corresponding deep unfolding networks are also developed. Section V presents a robust data-aided Kalman filter based channel tracking method. Numerical results are provided in Section VI, followed by the conclusion drawn in Section VII.

## II System Model and Problem Formulation

### _System Model_

Consider an uplink MIMO system, where a base station (BS) equipped with \(M\) receive antennas serves \(K\) single-antenna users. The channel remains constant during a coherent block comprised of \(L\) symbols and varies across blocks. As shown in Fig. 1, we adopt the frame structure that contains \(N+1\) coherence blocks to save the pilot overhead. Specifically, in each frame, \(L_{P}\) pilot symbols followed by \(L_{D}\) data symbols are sent in the first block, while the remaining \(N\) blocks are only used for data transmission. Note that, compared to the conventional frame structure that inserts pilots in each coherence block [39], the pilot overhead reduces from \(\frac{L_{P}}{L}\) to \(\frac{L_{P}}{(N+1)L}\), which is thus much more spectrum-efficient. For a convenient reference, the main symbols used in this paper are summarized in Table I at the top of the next page.

To characterize the channel aging effect, we use the first-order Gauss-Markov process as in [40, 41]. Define \(\mathbf{\Lambda}\triangleq\mathrm{diag}\left(\rho_{1},\rho_{2},\cdots,\rho_{K }\right)\otimes\mathbf{I}_{M}\) and \(\mathbf{\bar{\Lambda}}\triangleq\mathrm{diag}\left(\sqrt{1-\rho_{1}^{2}}, \sqrt{1-\rho_{2}^{2}},\cdots,\sqrt{1-\rho_{K}^{2}}\right)\otimes\mathbf{I}_{M}\), where \(\rho_{k}\) is the channel temporal correlation coefficient of the \(k\)-th user related to the channel aging speed. Then, the channel variation between two neighboring blocks can be represented

\[\mathbf{h}[n]=\mathbf{\Lambda}\mathbf{h}[n-1]+\mathbf{\bar{\Lambda}}\mathbf{w }[n],n=2,3,\cdots,N+1, \tag{1}\]

Fig. 1: Frame structure.
where \(\mathbf{h}[n]\in\mathbb{C}^{MK\times 1}\) denotes the vectorized form of the channel matrix \(\mathbf{H}[n]\in\mathbb{C}^{M\times K}\) at the \(n\)-th block and follows a Gaussian distribution with zero mean and covariance \(\mathbf{C_{h}}\), and \(\mathbf{w}[n]\in\mathbb{C}^{MK\times 1}\) denotes the zero-mean spatially correlated innovation process with covariance \(\mathbf{C_{h}}\). Furthermore, it can be easily inferred that

\[\mathbf{h}[n]=\mathbf{\Lambda}^{n-1}\mathbf{h}[1]+\hat{\mathbf{w}}[n],n=2,3, \cdots,N+1, \tag{2}\]

where \(\mathbf{h}[1]\) represents the channel vector of the first block, and the channel error \(\mathbf{\tilde{w}}[n]\) follows a Gaussian distribution with zero mean and covariance \(\mathbf{C_{\hat{w}}}[n]\), which is given by

\[\mathbf{C_{\hat{w}}}[n]=\sum_{n^{\prime}=1}^{n-1}\mathbf{\Lambda}^{n^{\prime} -1}\bar{\mathbf{\Lambda}}\mathbf{C_{h}}\bar{\mathbf{\Lambda}}\mathbf{\Lambda} ^{n^{\prime}-1},n=2,3,\cdots,N+1. \tag{3}\]

For the first block, letting \(\mathbf{S}_{P}\in\mathbb{C}^{K\times L_{p}}\) denote the pilot sequences, the received pilot signal can be expressed as

\[\mathbf{Y}_{P}=\mathbf{H}[1]\mathbf{S}_{P}+\mathbf{Z}_{P}, \tag{4}\]

where \(\mathbf{Z}_{P}\in\mathbb{C}^{M\times L_{p}}\) is the noise matrix whose elements are independently identically distributed (i.i.d.) Gaussian variables with zero mean and variance \(\sigma^{2}\). Vectorizing (4) yields

\[\mathbf{y}_{P}=\mathbf{Ph}[1]+\mathbf{z}_{P}, \tag{5}\]

where \(\mathbf{y}_{P}=\mathrm{vec}(\mathbf{Y}_{P})\), \(\mathbf{P}=\mathbf{S}_{P}^{T}\otimes\mathbf{I}_{M}\), and \(\mathbf{z}_{P}=\mathrm{vec}(\mathbf{Z}_{P})\). Considering the commonly used orthogonal pilots with the normalized average power of each pilot symbol, i.e., \(\mathbf{S}_{P}\mathbf{S}_{P}^{H}=L_{P}\mathbf{I}_{K}\), it then follows that \(\mathbf{P}^{H}\mathbf{P}=L_{P}\mathbf{I}_{MK}\). From (2), (3), and (5), the linear MMSE (LMMSE) channel estimate of the \(n\)-th block with given \(\mathbf{y}_{P}\) can be obtained via [40]

\[\begin{array}{c}\hat{\mathbf{h}}[n]=\mathbf{\Lambda}^{n-1}\mathbf{C_{h}} \big{(}L_{P}\mathbf{C_{h}}+\sigma^{2}\mathbf{I}_{MK}\big{)}^{-1}\mathbf{P}^{ H}\mathbf{y}_{P},\\ n=1,2,\cdots,N+1.\end{array} \tag{6}\]

Accordingly, we define the CSI uncertainty as \(\mathbf{\Delta}\mathbf{h}[n]=\mathbf{h}[n]-\hat{\mathbf{h}}[n]\), which is Gaussian distributed with zero mean and the covariance written as

\[\mathbf{\Sigma_{h}}[n]=\left\{\begin{array}{c}\sigma^{2}\mathbf{C_{h}}\big{(} L_{P}\mathbf{C_{h}}+\sigma^{2}\mathbf{I}_{MK}\big{)}^{-1},\qquad n=1,\\ \sigma^{2}\mathbf{\Lambda}^{n-1}\mathbf{C_{h}}\big{(}L_{P}\mathbf{C_{h}}+ \sigma^{2}\mathbf{I}_{MK}\big{)}^{-1}\mathbf{\Lambda}^{n-1}+\mathbf{C_{w}}[n ],\\ n=2,3,\cdots,N+1.\end{array}\right. \tag{7}\]

### _Problem Formulation_

We now consider the data transmission. Denote \(\mathbf{x}_{t}[n]\in\mathbb{C}^{K\times 1}\) as the transmitted vector at the \(t\)-th time slot of the \(n\)-th block, each entry of which is drawn from a quadrature amplitude modulation (QAM) constellation set \(\mathcal{X}\). Then, the received signal is given by

\[\begin{array}{c}\mathbf{y}_{t}[n]=\mathbf{H}[n]\mathbf{x}_{t}[n]+\mathbf{z }_{t}[n],\\ t=1,2,\cdots,L,\quad n=1,2,\cdots,N+1,\end{array} \tag{8}\]

where the noise vector \(\mathbf{z}_{t}[n]\sim\mathcal{CN}\left(\mathbf{0},\sigma^{2}\mathbf{I}_{M}\right)\). For simplicity of the notation, the indices \(t\) and \(n\) will be omitted if there is no confusion.

With the imperfect CSI in (6) available, one can recover the transmit signals by directly using the channel estimate as if it is perfect, which corresponds to the mismatched ML criterion:

\[\mathbf{\hat{x}}=\operatorname*{arg\,min}_{\mathbf{x}\in\mathcal{X}^{K}}\Big{\|} \mathbf{y}-\mathbf{\hat{H}}\mathbf{x}\Big{\|}_{2}^{2}, \tag{9}\]

where \(\mathbf{\hat{H}}\) is the matrix form of \(\mathbf{\hat{h}}\). However, this is not the optimal criterion since the CSI uncertainty is neglected. To fill this gap, we first recast (8) by

\[\mathbf{y}=\mathbf{X}\mathbf{h}+\mathbf{z}, \tag{10}\]

where \(\mathbf{X}=\mathbf{x}^{T}\otimes\mathbf{I}_{M}\). Then, by regarding \(\mathbf{h}\) as a conditional Gaussian vector with conditional mean \(\mathbf{\hat{h}}\) and covariance \(\mathbf{\Sigma_{h}}\) and performing some mathematical manipulations as in [18], the robust ML function can be derived as

\[p\left(\mathbf{y}\left|\mathbf{X},\mathbf{\hat{h}}\right.\right)=C\det\left( \mathbf{R}_{\mathbf{X}}^{-1}\right)\exp\left(\mathbf{q}_{\mathbf{X}}^{H} \mathbf{R}_{\mathbf{X}}^{-1}\mathbf{q}\mathbf{x}\right), \tag{11}\]

where \(C\) is a constant, \(\mathbf{q}_{\mathbf{X}}=\frac{1}{\sigma^{2}}\mathbf{X}^{H}\mathbf{y}+\mathbf{ \Sigma}_{\mathbf{h}}^{-1}\mathbf{\hat{h}}\), \(\mathbf{R}_{\mathbf{X}}=\frac{1}{\sigma^{2}}\mathbf{X}^{H}\mathbf{X}+\mathbf{ \Sigma}_{\mathbf{h}}^{-1}\), and \(\mathbf{\Sigma_{h}}\) is given in (7). Furthermore, maximizing (11) yields the subsequent robust ML detection criterion:

\[\mathbf{\hat{x}}=\operatorname*{arg\,min}_{\mathbf{x}\in\mathcal{X}^{K}}\ln\det \left(\mathbf{R}_{\mathbf{X}}\right)-\mathbf{q}_{\mathbf{X}}^{H}\mathbf{R}_{ \mathbf{X}}^{-1}\mathbf{q}_{\mathbf{X}}. \tag{12}\]

Generally, due to the \(K\)-dimensional discrete constraint \(\mathbf{x}\in\mathcal{X}^{K}\), the global optimal solution of the above problem can only be acquired via an exhaustive search, which requires an unacceptable exponential complexity. In the following section, we provide a neural network based solution to the problem with excellent performance and much reduced complexity.
## III Robust ADMM Detector Based Network Design

In this section, we first devise an ADMM-based algorithm to address problem (12). Based on this, we then provide the design of the proposed robust detection network RADMMNet.

### _Robust ADMM Detector_

In order to develop an efficient robust detector, we first reformulate problem (12) into a tractable form. It is known from [38] that the real and the imaginary parts of a \(4^{Q}\)-QAM signal can be decomposed into a weighted sum of \(Q\) binary variables, respectively. Using this representation, problem (12) can be equivalently expressed by

\[\begin{split}&\min_{\mathbf{x},\{\mathbf{v}_{q},\mathbf{v}_{q}\}} \ln\det\left(\mathbf{R}_{\mathbf{X}}\right)-\mathbf{q}_{\mathbf{X}}^{H} \mathbf{R}_{\mathbf{X}}^{-1}\mathbf{q}_{\mathbf{X}}\\ &\text{s.t.}\quad\mathbf{x}-\frac{1}{\alpha_{Q}}\sum\limits_{q=1} ^{Q}2^{q-1}\mathbf{v}_{q}=\mathbf{0},\\ &\quad\mathrm{Re}\{\mathbf{v}_{q}\},\mathrm{Im}\{\mathbf{v}_{q}\} \in\{-1,1\}^{K},q=1,\cdots,Q,\end{split} \tag{13}\]

where \(\alpha_{Q}\) is the power normalization factor of QAM signals and \(\{-1,1\}^{K}\) denotes the binary set of \(K\times 1\) vectors with each entry being \(1\) or \(-1\). Furthermore, we relax the discrete binary constraints by the boxed constraints and then impose a sum of quadratic penalty terms to the objective, yielding the following problem:

\[\begin{split}&\min_{\mathbf{x},\{\mathbf{v}_{q},\mathbf{v}_{q}\}} \ln\det\left(\mathbf{R}_{\mathbf{X}}\right)-\mathbf{q}_{\mathbf{X}}^{H} \mathbf{R}_{\mathbf{X}}^{-1}\mathbf{q}\mathbf{x}-\sum\limits_{q=1}^{Q}\beta_{ q}\left\|\mathbf{v}_{q}\right\|_{2}^{2}\\ &\mathrm{s.t.}\quad\mathbf{x}-\frac{1}{\alpha_{Q}}\sum\limits_{q= 1}^{Q}2^{q-1}\mathbf{v}_{q}=\mathbf{0},\\ &\quad\mathrm{Re}\{\mathbf{v}_{q}\},\mathrm{Im}\{\mathbf{v}_{q}\} \in[-1,1]^{K},q=1,\cdots,Q,\end{split} \tag{14}\]

where \(\beta_{q}>0,q=1,\cdots,Q\) are the penalty parameters, and \([-1,1]^{K}\) denotes the set of \(K\times 1\) vectors with each entry lying in the interval \([-1,1]\). Note that the added penalty term \(-\sum\limits_{q=1}^{Q}\beta_{q}\left\|\mathbf{v}_{q}\right\|_{2}^{2}\) makes the integer solutions more preferable and therefore tightens the boxed constraints [33]. Based on the reformulation, we apply the ADMM algorithm to tackle the problem as follows [8].

To facilitate the application of the ADMM framework, we construct the scaled augmented Lagrangian function of problem (14) as

\[\begin{split}& L_{\mu}\left(\mathbf{x},\{\mathbf{v}_{q}\}_{q=1}^{Q}, \boldsymbol{\lambda}\right)=\ln\det(\mathbf{R}_{\mathbf{X}})-\mathbf{q}_{ \mathbf{X}}^{H}\mathbf{R}_{\mathbf{X}}^{-1}\mathbf{q}\mathbf{x}\\ &-\sum\limits_{q=1}^{Q}\beta_{q}\left\|\mathbf{v}_{q}\right\|_{2} ^{2}+\frac{\mu}{2}\left\|\mathbf{x}-\frac{1}{\alpha_{Q}}\sum\limits_{q=1}^{Q} 2^{q-1}\mathbf{v}_{q}+\boldsymbol{\lambda}\right\|_{2}^{2},\end{split} \tag{15}\]

where \(\boldsymbol{\lambda}\) and \(\mu>0\) denote the scaled dual variables and the corresponding penalty parameter, respectively. Thus, the ADMM procedure can be described as follows:

\[\{\mathbf{v}_{q}^{i}\}_{q=1}^{Q}=\operatorname*{arg\,min}_{\mathbf{v}_{q}\in \mathbb{B}^{K}}L_{\mu}\left(\mathbf{x}^{i-1},\{\mathbf{v}_{q}\}_{q=1}^{Q}, \boldsymbol{\lambda}^{i-1}\right), \tag{16a}\] \[\mathbf{x}^{i}=\operatorname*{arg\,min}_{\mathbf{x}}L_{\mu}\left(\mathbf{x}, \{\mathbf{v}_{q}^{i}\}_{q=1}^{Q},\boldsymbol{\lambda}^{i-1}\right),\] (16b) \[\boldsymbol{\lambda}^{i}=\boldsymbol{\lambda}^{i-1}+\mathbf{x}^{i}-\frac{1}{ \alpha_{Q}}\sum\limits_{q=1}^{Q}2^{q-1}\mathbf{v}_{q}^{i}, \tag{16c}\]

where the superscript \(i\) denotes the index of ADMM iterations and \(\mathbb{B}^{K}\) denotes the set of \(K\times 1\) vectors, both the real and the imaginary parts of whose entries belong to the interval \([-1,1]\).

For subproblem (16a), the difficulty of obtaining the optimal solution lies in the fact that \(\{\mathbf{v}_{q}\}_{q=1}^{Q}\) are coupled with each other. Nevertheless, we find that \(L_{\mu}(\mathbf{x}^{i-1},\{\mathbf{v}_{q}\}_{q=1}^{Q},\boldsymbol{\lambda}^{i -1})\) is convex with respect to (w.r.t.) each \(\mathbf{v}_{q}\) when the condition \(\mu 4^{q-1}-2\alpha_{Q}^{2}\beta_{q}>0\) is satisfied, which inspires us to apply the block coordinate descent (BCD) method [42]. Consequently, by taking the derivative of (15) w.r.t. each \(\mathbf{v}_{q}\) and setting it to zero, a closed-form solution to \(\mathbf{v}_{q}\) can be achieved by

\[\mathbf{v}_{q}^{i}=\mathcal{P}_{\mathbb{B}^{K}}\left\{\frac{2^{q-1}\mu\boldsymbol {\eta}_{q}^{i}}{\mu^{q-1}-2\alpha_{Q}^{2}\beta_{q}}\right\},\quad q=1,2, \cdots,Q, \tag{17}\]

where \(\boldsymbol{\eta}_{q}^{i}=\alpha_{Q}\left(\mathbf{x}^{i-1}+\boldsymbol{\lambda }^{i-1}\right)-\sum\limits_{p<q}2^{p-1}\mathbf{v}_{p}^{i}-\sum\limits_{p>q}2^{ p-1}\mathbf{v}_{p}^{i-1}\) and \(\mathcal{P}_{\mathbb{B}^{K}}\{\cdot\}\) means projecting the real and the imaginary parts of each entry of the input vector onto the interval \([-1,1]\).

The remaining subproblem (16b) is more challenging due to its complicated non-convex objective. To make it tractable, we replace the non-trivial function \(\ln\det\left(\mathbf{R}_{\mathbf{X}}\right)-\mathbf{q}_{\mathbf{X}}^{H} \mathbf{R}_{\mathbf{X}}^{-1}\mathbf{q}_{\mathbf{X}}\) with an appropriate upperbound surrogate function in each ADMM iteration. In this way, the original minimization problem is reduced to minimizing the surrogate function that is easier to handle. Concretely, we first establish an upper bound to the function \(\ln\det\left(\mathbf{R}_{\mathbf{X}}\right)\). Since the function is concave in \(\mathbf{R}_{\mathbf{X}}\), it is majorized by its first-order Taylor expansion as

\[\begin{split}\ln\det\left(\mathbf{R}_{\mathbf{X}}\right)& \leq\ln\det\left(\mathbf{R}_{\mathbf{X}^{i-1}}\right)+\mathrm{tr} \left(\mathbf{R}_{\mathbf{X}^{i-1}}^{-1}\left(\mathbf{R}_{\mathbf{X}}-\mathbf{R }_{\mathbf{X}^{i-1}}\right)\right)\\ &\stackrel{{(a)}}{{=}}\frac{1}{\sigma^{2}}\mathrm{tr} \left(\mathbf{R}_{\mathbf{X}^{i-1}}^{-1}\mathbf{X}^{H}\mathbf{X}\right)+\mathrm{ constant}\\ &\stackrel{{(b)}}{{=}}\mathbf{x}^{H}\mathbf{C}^{i-1} \mathbf{x}+\mathrm{constant},\end{split} \tag{18}\]

where (a) uses the definition of \(\mathbf{R}_{\mathbf{X}}=\frac{1}{\sigma^{2}}\mathbf{X}^{H}\mathbf{X}+ \mathbf{\Sigma}_{\mathbf{h}}^{-1}\), and (b) is acquired based on \(\mathbf{X}=\mathbf{x}^{T}\otimes\mathbf{I}_{M}\) with \(\mathbf{C}^{i-1}\) defined by

\[\mathbf{C}^{i-1}=\frac{1}{\sigma^{2}}\big{(}\mathbf{E}^{H}\left(\mathbf{R}_{ \mathbf{X}^{i-1}}^{-1}\odot\left(\mathbf{1}_{K\times K}\otimes\mathbf{I}_{M} \right)\right)\mathbf{E}\big{)}^{T}, \tag{19}\]

where \(\mathbf{E}=\mathbf{I}_{K}\otimes\mathbf{1}_{M\times 1}\), \(\mathbf{A}\otimes\mathbf{B}\) and \(\mathbf{A}\odot\mathbf{B}\) represent the Kronecker product and the Hadamard product of \(\mathbf{A}\) and \(\mathbf{B}\), respectively, \(\mathbf{I}_{k}\) denotes the identity matrix of size \(k\times k\), and \(\mathbf{1}_{k_{1}\times k_{2}}\) denotes the all-ones vector or matrix of size \(k_{1}\times k_{2}\).

On the other hand, concerning the term \(\mathbf{q}_{\mathbf{X}}^{H}\mathbf{R}_{\mathbf{X}}^{-1}\mathbf{q}_{\mathbf{X}}\) which is jointly convex in \(\mathbf{q}_{\mathbf{X}}\) and \(\mathbf{R}_{\mathbf{X}}\), we obtain its lower bound via the first-order Taylor expansion as

\[\begin{split}&\mathbf{q}_{\mathbf{X}}^{H}\mathbf{R}_{\mathbf{X}}^{-1} \mathbf{q}_{\mathbf{X}}\\ &\geq-\mathrm{tr}\left(\mathbf{R}_{\mathbf{X}^{i-1}\mathbf{q}_{ \mathbf{X}^{i-1}}\mathbf{q}_{\mathbf{X}^{i-1}}^{H}\mathbf{R}_{\mathbf
achieved by using \(\mathbf{X}=\mathbf{x}^{T}\otimes\mathbf{I}_{M}\). In addition, \(\mathbf{A}^{i-1}=\mathbf{y}\mathbf{q}_{\mathbf{X}^{i-1}}^{H}\mathbf{R}_{\mathbf{ X}^{i-1}}^{-1}\), \(\mathbf{B}^{i-1}=\mathbf{R}_{\mathbf{X}^{i-1}}^{-1}\mathbf{q}_{\mathbf{X}^{i-1} }\mathbf{q}_{\mathbf{X}^{i-1}}^{H}\mathbf{R}_{\mathbf{X}^{i-1}}^{-1}\), and \(\mathbf{A}^{i-1}\) and \(\mathbf{F}^{i-1}\) respectively take the forms:

\[\mathbf{d}^{i-1}=\frac{1}{\sigma^{2}}\left(\mathbf{1}_{1\times M}\left( \mathbf{A}^{i-1}\odot(\mathbf{1}_{1\times K}\otimes\mathbf{I}_{M})\right) \mathbf{E}\right)^{T}, \tag{21}\]

and

\[\mathbf{F}^{i-1}=\frac{1}{\sigma^{2}}\big{(}\mathbf{E}^{H}\left(\mathbf{B}^{i- 1}\ \odot(\mathbf{1}_{K\times K}\otimes\mathbf{I}_{M})\right)\mathbf{E}\big{)}^{T}. \tag{22}\]

Combining (18) and (20), problem (16b) boils down to minimizing the surrogate upper bound function as follows:

\[\mathbf{x}^{i}=\operatorname*{arg\,min}_{\mathbf{x}}\mathbf{x}^{ H}\big{(}\mathbf{C}^{i-1}+\mathbf{F}^{i-1}\big{)}\mathbf{x}-2\mathrm{Re}\{ \mathbf{x}^{H}\mathbf{d}^{i-1}\} \tag{23}\] \[+\frac{\mu}{2}\left\|\mathbf{x}+\boldsymbol{\psi}^{i-1}\right\|_ {2}^{2},\]

where \(\boldsymbol{\psi}^{i-1}=-\frac{1}{\alpha_{Q}}\sum\limits_{q=1}^{Q}2^{q-1} \mathbf{v}_{q}^{i}+\boldsymbol{\lambda}^{i-1}\). The objective function of problem (23) is convex quadratic w.r.t. \(\mathbf{x}\) since \(\mathbf{C}^{i-1}\) and \(\mathbf{F}^{i-1}\) are both positive semidefinite and \(\mu>0\). Thus, the update of \(\mathbf{x}\) can be easily obtained by

\[\mathbf{x}^{i}=\left(2\left(\mathbf{C}^{i-1}+\mathbf{F}^{i-1}\right)+\mu \mathbf{I}_{K}\right)^{-1}\left(-\mu\boldsymbol{\psi}^{i-1}+2\mathbf{d}^{i-1} \right). \tag{24}\]

At this step, we are able to address problem (14) by iteratively calculating the closed-form expressions in (17), (24) and (16c) until the convergence is reached, which is also summarized in Algorithm 1. Note that the convergence of ADMM to a general nonconvex problem remains an open issue, which is even harder with inexact optimization of the involved subproblems [8]. In this work, we only obtain inexact solutions to subproblems (16a) and (16b). Nonetheless, our ADMM algorithm can still converge according to the empirical results. In particular, the above ADMM algorithm will be used for developing a neural network based detector via the deep unfolding technique in the next subsection.

```
0:\(\mathbf{R}\) and \(\mathbf{q}\).
1: Initialize: \(\mathbf{s}_{0}\), \(\mathbf{r}_{0}=\mathbf{q}-\mathbf{R}\mathbf{s}_{0}\), \(\mathbf{p}_{0}=\mathbf{r}_{0}\), \(\xi=0\).
2:repeat
3:\(\xi\leftarrow\xi+1\).
4:\(\tau_{\xi}=\mathbf{r}_{\xi-1}^{H}\mathbf{r}_{\xi-1}/\mathbf{p}_{\xi-1}^{H} \mathbf{R}\mathbf{p}_{\xi-1}\).
5:\(\mathbf{s}_{\xi}=\mathbf{s}_{\xi-1}+\tau_{\xi}\mathbf{p}_{\xi-1}\).
6:\(\mathbf{r}_{\xi}=\mathbf{r}_{\xi-1}-\tau_{\xi}\mathbf{R}\mathbf{p}_{\xi-1}\).
7:\(v_{\xi}=\mathbf{r}_{\xi}^{H}\mathbf{r}_{\xi}/\mathbf{r}_{\xi-1}^{H}\mathbf{r} _{\xi-1}\).
8:\(\mathbf{p}_{\xi}=\mathbf{r}_{\xi}+v_{\xi}\mathbf{p}_{\xi-1}\).
9:until convergence.
10:\(\hat{\mathbf{s}}=\mathbf{s}_{\xi}\).
```

**Algorithm 2** CG Algorithm for Calculating \(\mathbf{R}^{-1}\mathbf{q}\)

### _Proposed RADMMNet_

In order to facilitate the proposed network design, we first perform some simplifications for the above ADMM algorithm. Specifically, the \(MK\)-dimensional matrix inversion \(\mathbf{R}_{\mathbf{X}^{i-1}}^{-1}\) involved in the calculation of \(\mathbf{C}^{i-1}\) in (24) requires a high complexity up to \(\mathcal{O}\left(M^{3}K^{3}\right)\), which can lead to a huge computational cost for both offline training and online calculation. To handle this, we further relax the upper bound in (18) by

\[\ln\det\left(\mathbf{R}_{\mathbf{X}}\right) \leq\mathbf{x}^{H}\mathbf{C}^{i-1}\mathbf{x}+\mathrm{constant} \tag{25}\] \[\leq\varepsilon^{i-1}\|\mathbf{x}\|_{2}^{2}+\mathrm{constant},\]

where \(\mathbf{C}^{i-1}\) is given in (19) and \(\varepsilon^{i-1}\) is set to satisfy \(\varepsilon^{i-1}\mathbf{I}_{K}\succeq\mathbf{C}^{i-1}\). We note that we regard \(\varepsilon^{i-1}\) as a trainable parameter in the proposed network, whose value is determined after the network training is completed. In addition, when computing \(\mathbf{x}^{i}\) using (24), we need to calculate \(\mathbf{R}_{\mathbf{X}^{i-1}}^{-1}\mathbf{q}_{\mathbf{X}^{i-1}}\) to obtain \(\mathbf{d}^{i-1}\) and \(\mathbf{F}^{i-1}\) (see (21) and (22)), which will also introduce a very high complexity of \(\mathcal{O}\left(M^{3}K^{3}+M^{2}K^{2}\right)\). To avoid this, we adopt the conjugate gradient descent (CG) method to approximate \(\hat{\mathbf{s}}=\mathbf{R}_{\mathbf{X}^{i-1}}^{-1}\mathbf{q}_{\mathbf{X}^{i-1}}\) by iteratively minimizing the quadratic function \(f(\mathbf{s})=\mathbf{s}^{H}\mathbf{R}_{\mathbf{X}^{i-1}}\mathbf{s}-\mathbf{q} _{\mathbf{X}^{i-1}}^{H}\mathbf{s}\)[43]. The detailed procedure of the CG method is described in Algorithm 2, where the subscript \(\mathbf{X}^{i-1}\) is temporarily omitted for the brevity of the notations. Specifically, line 5, line 6, and line 8 represent the update of the solution \(\mathbf{s}_{\xi}\), the residual \(\mathbf{r}_{\xi}\), and the conjugate direction \(\mathbf{p}_{\xi}\) in the \(\xi\)-th CG iteration,respectively, while line 4 and line 7 give the expressions of the corresponding step sizes \(\tau_{\xi}\) and \(v_{\xi}\). By resorting to the CG method, the complexity of computing \(\mathbf{R}_{\mathbf{X}^{i-1}}^{-1}\mathbf{q}_{\mathbf{X}^{i-1}}\) is reduced to \(\mathcal{O}\left(I_{\mathrm{CG}}M^{2}K^{2}\right)\), where \(I_{\mathrm{CG}}\) denotes the number of CG iterations. Moreover, the iterative nature of the CG method is also friendly to the neural network design1. On the other hand, it is noteworthy that a good initialization of \(\mathbf{s}_{0}\) can improve the convergence of the CG method, which is generally set to \(\mathbf{0}\) if there is no _a prior_ knowledge. Fortunately, we find that the value of \(\mathbf{R}_{\mathbf{X}^{i}}^{-1}\mathbf{q}_{\mathbf{X}^{i}}\) does not change dramatically as the iteration index \(i\) increases. Therefore, when using the CG method to approximate \(\mathbf{R}_{\mathbf{X}^{i-1}}^{-1}\mathbf{q}_{\mathbf{X}^{i-1}}\) in the \(i\)-th ADMM iteration (\(i>1\)), \(\mathbf{s}_{0}\) can be initialized to \(\mathbf{R}_{\mathbf{X}^{i-2}}^{-1}\mathbf{q}_{\mathbf{X}^{i-2}}\) obtained in the \((i-1)\)-th ADMM iteration, which can efficiently reduce the required number of CG iterations down to \(1\). This processing can also be regarded as transferring the intermediate variables to strength the connection between the network layers.

Footnote 1: Although there exist other iterative methods for the same purpose, such as the Newton method and the steepest gradient descent method, they suffer from a high complexity and a low convergence rate, respectively [32]. As for the CG method, it only involves matrix multiplications per iteration with a much lower complexity. Meanwhile, the solution is updated along an independent conjugate direction with an exactly calculated step size in each CG iteration, which means that one will go to the extreme along one direction and there will be no repeated directions in the search procedure, leading to a faster convergence rate. Therefore, the CG method is more preferable.

After performing the above modifications for the proposed ADMM algorithm, we are ready to develop a neural network solution to the robust ML detection problem in (12). We employ the deep unfolding technique to represent the ADMM
iterations by \(I_{\text{max}}\) layers. Note that the involved penalty parameters \(\mu\) and \(\{\beta_{q}\}_{q=1}^{Q}\), as well as \(\varepsilon\) involved in (25) can remarkably affect the performance and convergence of the proposed ADMM algorithm. In the proposed network, we set them to be trainable parameters which can be optimized via offline training, thus effectively avoiding a cumbersome numerical search. Furthermore, we adopt layer-wise parameters, i.e., \(\{\mu^{i},\{\beta_{q}^{i}\}_{q=1}^{Q},\varepsilon^{i}\}_{i=1}^{I_{\text{max}}}\), which can provide more degrees of freedom and potentially accelerate the convergence and improve the performance.

Note that (21) and (22) can be conveniently implemented in Tensorflow by convolutional layers with an identity matrix as the fixed filter. Instead of strictly following the derived expressions, we improve the learning ability of our network by replacing the fixed filter with a trainable one, so that the resulting deep unfolding network can also take advantage of the introduced data-driven structure. It yields:

\[\begin{split}\mathbf{d}^{i-1}&=\frac{1}{\sigma^{2} }\left(\mathbf{1}_{1\times M}\left(\mathbf{A}^{i-1}\odot(\mathbf{1}_{1\times K }\otimes\mathbf{W})\right)\mathbf{E}\right)^{T}\\ &=\mathrm{conv}\left(\frac{1}{\sigma^{2}}\left(\mathbf{A}^{i-1} \right)^{\mathrm{T}};\mathbf{W}\right),\end{split} \tag{26}\]

and

\[\begin{split}\mathbf{F}^{i-1}&=\frac{1}{\sigma^{2 }}\big{(}\mathbf{E}^{H}\left(\mathbf{B}^{i-1}\odot(\mathbf{1}_{K\times K} \otimes\mathbf{W})\right)\mathbf{E}\big{)}^{T}\\ &=\mathrm{conv}\left(\frac{1}{\sigma^{2}}\left(\mathbf{B}^{i-1} \right)^{\mathrm{T}};\mathbf{W}\right),\end{split} \tag{27}\]

where \(\mathbf{E}=\mathbf{I}_{K}\otimes\mathbf{1}_{M\times 1}\), and \(\mathrm{conv}\left(\cdot;\mathbf{W}\right)\) means the convolutional operation with \(\mathbf{W}\in\mathbb{R}^{M\times M}\) being the trainable filter. To guarantee the consistence of the output dimensions, the stride and the depth of the convolutional layers are set to \(M\) and \(1\), respectively. Unlike other layer-wise parameters, we share the same \(\mathbf{W}\) for each network layer to limit the number of trainable parameters, which can effectively alleviate the training difficulty and avoid the overfitting problem.

To summarize, the \(i\)-th layer of the proposed network can be constructed as follows:

\[\mathbf{v}_{q}^{i}=\mathcal{P}_{\mathbb{B}^{K}}\left\{\frac{2^{q-1}\mu^{i} \boldsymbol{\eta}_{q}^{i}}{\mu^{i}4^{q-1}-2\alpha_{Q}^{2}\beta_{q}^{i}}\right\},\quad q=1,2,\cdots,Q, \tag{28a}\] \[\mathbf{x}^{i}=\left(2\mathbf{F}^{i-1}+\left(\varepsilon^{i}+\mu^{i} \right)\mathbf{I}_{K}\right)^{-1}\left(-\mu^{i}\boldsymbol{\psi}^{i-1}+2 \mathbf{d}^{i-1}\right),\] (28b) \[\boldsymbol{\lambda}^{i}=\boldsymbol{\lambda}^{i-1}+\mathbf{x}^{i}-\frac{1}{ \alpha_{Q}}\sum\limits_{q=1}^{Q}2^{q-1}\mathbf{v}_{q}^{i}, \tag{28c}\]

where \(\boldsymbol{\eta}_{q}^{i}=\alpha_{Q}\left(\mathbf{x}^{i-1}+\boldsymbol{\lambda }^{i-1}\right)-\sum\limits_{p<q}2^{p-1}\mathbf{v}_{p}^{i}-\sum\limits_{p>q}2^ {p-1}\mathbf{v}_{p}^{i-1}\), and \(\boldsymbol{\psi}^{i-1}=-\frac{1}{\alpha_{Q}}\sum\limits_{q=1}^{Q}2^{q-1} \mathbf{v}_{q}^{i}+\boldsymbol{\lambda}^{i-1}\). Furthermore, in the expressions of \(\boldsymbol{\psi}^{i-1}\) and (28c), we additionally introduce a trainable parameter \(\varsigma^{i}\) by replacing the original term \(-\frac{1}{\alpha_{Q}}\sum\limits_{q=1}^{Q}2^{q-1}\mathbf{v}_{q}^{i}\) with the term \(-\frac{\varsigma^{i}}{\alpha_{Q}}\sum\limits_{q=1}^{Q}2^{q-1}\mathbf{v}_{q}^{i }-\left(1-\varsigma^{i}\right)\mathbf{x}^{i-1}\). This corresponds to the over-relaxation scheme as suggested by [8, 44], which is expected to improve the convergence of ADMM iterations. Finally, the block diagram of the proposed RADMMNet is illustrated in Fig. 2.

## IV Low-Complexity Robust ADMM Detector Based Network Design

Although we have bypassed the high-dimensional matrix inversion, RADMMNet can still be computationally intensive for the system with a large number of antennas, which imposes huge burdens on both offline and online processes. In this section, we derive a low-complexity robust ADMM detector and further provide the corresponding neural network design.

### _Low-Complexity Robust ADMM Detector_

To begin with, we rewrite (8) as

\[\mathbf{y}=\hat{\mathbf{H}}\mathbf{x}+\boldsymbol{\Delta}\mathbf{H}\mathbf{x}+ \mathbf{z}, \tag{29}\]

where \(\hat{\mathbf{H}}\) is the matrix form of \(\hat{\mathbf{h}}\) and \(\boldsymbol{\Delta}\mathbf{H}=\mathbf{H}-\hat{\mathbf{H}}\) is the corresponding zero-mean Gaussian distributed channel error. The two error covariance matrix of \(\boldsymbol{\Delta}\mathbf{H}\), i.e., \(\boldsymbol{\Sigma}_{\mathbf{H}}\triangleq\mathbb{E}\left\{\left(\boldsymbol{ \Delta}\mathbf{H}\mathbf{)}\!\left(\boldsymbol{\Delta}\mathbf{H}\right)^{H} \right\}\), can be derived from \(\boldsymbol{\Sigma}_{\mathbf{h}}\) in (7) by

\[\boldsymbol{\Sigma}_{\mathbf{H}}=\left(\mathbf{1}_{1\times K}\otimes\mathbf{I}_ {M}\right)\left(\boldsymbol{\Sigma}_{\mathbf{h}}\odot\left(\mathbf{I}_{K} \otimes\mathbf{1}_{M\times M}\right)\right)\left(\mathbf{1}_{K\times 1}\otimes \mathbf{I}_{M}\right). \tag{30}\]

Based on (29), the key idea behind the low-complexity design is to impose a Gaussian approximation on the residual term \(\tilde{\mathbf{r}}=\boldsymbol{\Delta}\mathbf{H}\mathbf{x}+\mathbf{z}\). That is, we assume that \(\tilde{\mathbf{r}}\) follows a Gaussian distribution with zero mean and the covariance \(\mathbf{C}_{\mathbf{r}}\triangleq\mathbb{E}\left\{\left(\boldsymbol{\Delta} \mathbf{H}\mathbf{x}+\mathbf{z}\right)\!\left(\boldsymbol{\Delta}\mathbf{H} \mathbf{x}+\mathbf{z}\right)^{H}\right\}=\boldsymbol{\Sigma}_{\mathbf{H}}+ \sigma^{2}\mathbf{I}_{M}\), where

Fig. 2: Block diagram of RADMMNet.
the expectation is taken w.r.t. \(\mathbf{\Delta}\mathbf{H}\), \(\mathbf{z}\), and \(\mathbf{x}\). Thus, the robust ML criterion can be simplified to

\[\mathbf{\hat{x}}=\operatorname*{arg\,min}_{\mathbf{x}\in\mathcal{K}^{K}}\left( \mathbf{y}-\hat{\mathbf{H}}\mathbf{x}\right)^{H}\mathbf{C}_{\tilde{\mathbf{r}} }^{-1}\left(\mathbf{y}-\hat{\mathbf{H}}\mathbf{x}\right). \tag{31}\]

We also resort to the ADMM framework to address problem (31), where the details are analogous to those in Section III-A. Specifically, following the reformulation of problem (14), we recast problem (31) by

\[\min_{\mathbf{x},\{\mathbf{u}_{q}\},\forall q\}\left(\mathbf{y} -\hat{\mathbf{H}}\mathbf{x}\right)^{H}\mathbf{C}_{\tilde{\mathbf{r}}}^{-1} \left(\mathbf{y}-\hat{\mathbf{H}}\mathbf{x}\right)-\sum\limits_{q=1}^{Q} \kappa_{q}\left\|\mathbf{u}_{q}\right\|_{2}^{2} \tag{32}\] \[\mathrm{s.t.}\quad\mathbf{x}-\frac{1}{\alpha_{Q}}\sum\limits_{q=1 }^{Q}2^{q-1}\mathbf{u}_{q}=\mathbf{0},\]

\[\mathrm{Re}\{\mathbf{u}_{q}\},\mathrm{Im}\{\mathbf{u}_{q}\}\in[-1,1]^{K},q=1, \cdots,Q,\]

where \(\kappa_{q}>0,q=1,\cdots,Q\) are the penalty parameters. Then, the corresponding augmented Lagrangian function can be expressed as

\[L_{\delta}\left(\mathbf{x},\{\mathbf{u}_{q}\}_{q=1}^{Q},\mathbf{ \theta}\right)=\left(\mathbf{y}-\hat{\mathbf{H}}\mathbf{x}\right)^{H}\mathbf{ C}_{\tilde{\mathbf{r}}}^{-1}\left(\mathbf{y}-\hat{\mathbf{H}}\mathbf{x}\right) \tag{33}\] \[-\sum\limits_{q=1}^{Q}\kappa_{q}\left\|\mathbf{u}_{q}\right\|_{2} ^{2}+\frac{\delta}{2}\left\|\mathbf{x}-\frac{1}{\alpha_{Q}}\sum\limits_{q=1}^ {Q}2^{q-1}\mathbf{u}_{q}+\mathbf{\theta}\right\|_{2}^{2},\]

where \(\mathbf{\theta}\) and \(\delta>0\) denote the scaled dual variables and the corresponding penalty parameter, respectively. Therefore, the \(i\)-th ADMM iteration can be given by

\[\{\mathbf{u}_{q}^{i}\}_{q=1}^{Q}=\operatorname*{arg\,min}_{\mathbf{u}_{q}\in \mathbb{B}^{K}}L_{\delta}\left(\mathbf{x}^{i-1},\{\mathbf{u}_{q}\}_{q=1}^{Q}, \mathbf{\theta}^{i-1}\right), \tag{34a}\] \[\mathbf{x}^{i}=\operatorname*{arg\,min}_{\mathbf{x}}L_{\delta}\left(\mathbf{x}, \{\mathbf{u}_{q}^{i}\}_{q=1}^{Q},\mathbf{\theta}^{i-1}\right),\] (34b) \[\mathbf{\theta}^{i}=\mathbf{\theta}^{i-1}+\mathbf{x}^{i}-\frac{1}{\alpha_{Q}}\sum \limits_{q=1}^{Q}2^{q-1}\mathbf{u}_{q}^{i}. \tag{34c}\]

Note that subproblem (34a) is almost the same as subproblem (16a) whose solution can be similarly achieved by

\[\mathbf{u}_{q}^{i}=\mathcal{P}_{\mathbb{B}^{K}}\left\{\frac{2^{q-1}\delta \mathbf{\omega}_{q}^{i}}{\delta 4^{q-1}-2\alpha_{Q}^{2}\kappa_{q}}\right\},\quad q=1,2,\cdots,Q, \tag{35}\]

where \(\mathbf{\omega}_{q}^{i}=\alpha_{Q}\left(\mathbf{x}^{i-1}+\mathbf{\theta}^{i-1}\right) -\sum\limits_{p<q}2^{p-1}\mathbf{u}_{p}^{i}-\sum\limits_{p>q}2^{p-1}\mathbf{u} _{p}^{i-1}\). Moreover, by dropping the terms independent of \(\mathbf{x}\), we equivalently rewrite subproblem (34b) as

\[\mathbf{x}^{i}=\operatorname*{arg\,min}_{\mathbf{x}}\mathbf{x}^{H}\mathbf{\Phi} \mathbf{x}-2\mathrm{Re}\{\mathbf{x}^{H}\mathbf{\gamma}^{i-1}\}, \tag{36}\]

where \(\mathbf{\Phi}=\hat{\mathbf{H}}^{H}\mathbf{C}_{\tilde{\mathbf{r}}}^{-1}\hat{ \mathbf{H}}+\frac{\delta}{2}\mathbf{I}_{K}\) and \(\mathbf{\gamma}^{i-1}=\hat{\mathbf{H}}^{H}\mathbf{C}_{\tilde{\mathbf{r}}}^{-1} \mathbf{y}+\frac{\delta}{2\alpha_{Q}}\sum\limits_{q=1}^{Q}2^{q-1}\mathbf{u}_{q }^{i}-\frac{\delta}{2}\mathbf{\theta}^{i-1}\). Since \(\mathbf{\Phi}\) is positive definite, the optimal solution of \(\mathbf{x}\) in each ADMM update is

\[\mathbf{x}^{i}=\mathbf{\Phi}^{-1}\mathbf{\gamma}^{i-1}. \tag{37}\]

The above procedure is summarized in Algorithm 3. Different from Algorithm 1 for problem (14), the update of \(\mathbf{x}\) in Algorithm 3 is obtained in an exact closed form without the need to construct a surrogate function. The simplification can effectively reduce the computational complexity (see Section VI-C for details).

```
0:\(\mathbf{y}\), \(\hat{\mathbf{H}}\), and \(\mathbf{\Sigma}_{\mathbf{H}}\).
1:Initialize: \(\mathbf{x}^{0}\), \(\{\mathbf{u}_{q}^{0}\}_{q=1}^{Q}\), \(\mathbf{\theta}^{0}\), \(i=0\).
2:repeat
3:\(i\gets i+1\).
4: Sequentially update \(\mathbf{u}_{q}^{i}\) for \(q=1,2,\cdots,Q\) via (35).
5: Update \(\mathbf{x}^{i}\) via (37).
6: Update \(\mathbf{\theta}^{i}\) via (34c).
7:until convergence.
8:\(\hat{\mathbf{x}}=\mathbf{x}^{i}\).
```

**Algorithm 3** Low-Complexity Robust ADMM Detection Algorithm

### _Proposed LCRADMMNet_

Similar to RADMMNet, it is straightforward to build a model-driven network, referred to as LCRADMMNet, by unfolding the iterations of the derived low-complexity robust ADMM detection algorithm. Note that \(\mathbf{\Phi}^{-1}=\left(\hat{\mathbf{H}}^{H}\mathbf{C}_{\tilde{\mathbf{r}}}^{-1} \hat{\mathbf{H}}+\frac{\delta}{2}\mathbf{I}_{K}\right)^{-1}\) in (37) involves an \(M\)-dimensional matrix inversion (\(\mathbf{C}_{\tilde{\mathbf{r}}}^{-1}\)) and a \(K\)-dimensional matrix inversion (\(\mathbf{\Phi}^{-1}\)), which, fortunately, only needs to be calculated once during the ADMM iterations as long as \(\delta\) remains constant. To inherit the complexity advantage, we set \(\delta\) to be trainable and share it for different layers. Besides, a set of layer-wise trainable parameters, denoted by \(\{\{\kappa_{q}^{i}\}_{q=1}^{Q}\}_{i=1}^{I_{\mathbf{\gamma}}}\), are used to take the role of the penalty parameters \(\{\kappa_{q}\}_{q=1}^{Q}\) involved in Algorithm 3.

Therefore, the \(i\)-th layer of the proposed network can be expressed as:

\[\mathbf{u}_{q}^{i}=\mathcal{P}_{\mathbb{B}^{K}}\left\{\frac{2^{q-1}\delta\mathbf{ \omega}_{q}^{i}}{\delta 4^{q-1}-2\alpha_{Q}^{2}\kappa_{q}^{i}}\right\},\quad q=1,2,\cdots,Q, \tag{38a}\] \[\mathbf{x}^{i}=(\mathbf{\Phi})^{-1}\mathbf{\gamma}^{i-1},\] (38b) \[\mathbf{\theta}^{i}=\mathbf{\theta}^{i-1}+\mathbf{x}^{i}-\frac{1}{\alpha_{Q}}\sum \limits_{q=1}^{Q}2^{q-1}\mathbf{u}_{q}^{i}, \tag{38c}\]

where \(\mathbf{\omega}_{q}^{i}=\alpha_{Q}\left(\mathbf{x}^{i-1}+\mathbf{\theta}^{i-1}\right) -\sum\limits_{p<q}2^{p-1}\mathbf{u}_{p}^{i}-\sum\limits_{p>q}2^{p-1}\mathbf{u}_ {p}^{i}\), \(\mathbf{\Phi}=\hat{\mathbf{H}}^{H}\mathbf{C}_{\tilde{\mathbf{r}}}^{-1}\hat{\mathbf{H}} +\frac{\delta}{2}\mathbf{I}_{K}\), and \(\mathbf{\gamma}^{i-1}=\hat{\mathbf{H}}^{H}\mathbf{C}_{\tilde{\mathbf{r}}}^{-1} \mathbf{y}+\frac{\delta}{2\alpha_{Q}}\sum\limits_{q=1}^{Q}2^{q-1}\mathbf{u}_ {q}^{i}-\frac{\delta}{2}\mathbf{\theta}^{i-1}\). Here we can also apply the over-relaxation scheme by substituting \(-\frac{1}{\alpha_{Q}}\sum\limits_{q=1}^{Q}2^{q-1}\mathbf{u}_{q}^{i}\) involved in the expressions of \(\mathbf{\gamma}^{i-1}\) and (38c) with \(-\frac{v^{i}}{\alpha_{Q}}\sum\limits_{q=1}^{Q}2^{q-1}\mathbf{u}_{q}^{i}-\left(1-v^ {i}\right)\mathbf{x}^{i-1}\), where \(v^{i}\) is the trainable relaxation parameter. We illustrate the proposed LCRADMMNet in Fig. 3.

## V RDAKF-Based Channel Tracking

We have developed two robust detection networks to suppress the performance deterioration caused by the imperfect CSI. However, the channel estimates of all
accumulated for the subsequent data blocks without pilots, leading to degraded detection performance. To mitigate the error propagation, we now present a novel robust data-aided Kalman filter (RDAKF)-based channel tracking method, which incorporates not only the channel estimate acquired from the previous block as prior information but also the robust data estimates of the current block for further updating.

The proposed channel tracking method can be described from the perspective of a Kalman filter, which consists of two key stages as detailed in the following.

### _One-step prediction stage_

First, the state equation is utilized for one-step prediction. In this work, the channel variation model (1) can be regarded as the state equation. Then, the channel vector of the \(n\)-th block can be predicted as

\[\hat{\mathbf{h}}[n\left|n-1\right.]=\mathbf{\Lambda}\hat{\mathbf{h}}[n-1 \left|n-1\right.],\quad n=2,3,\cdots,N+1, \tag{39}\]

and the corresponding channel error covariance is given by

\[\mathbf{\Sigma_{h}}[n\left|n-1\right.]=\mathbf{\Lambda}\mathbf{ \Sigma_{h}}[n-1\left|n-1\right.]\mathbf{\Lambda}+\mathbf{\tilde{\Lambda}} \mathbf{C_{h}}\mathbf{\tilde{\Lambda}}, \tag{40}\] \[n=2,3,\cdots,N+1,\]

where \(\hat{\mathbf{h}}[n\left|n-1\right.]\) and \(\mathbf{\Sigma_{h}}[n\left|n-1\right.]\) represent the _a prior_ channel estimate of the \(n\)-th block based on the previous estimate of the \((n-1)\)-th block and the corresponding _a prior_ error covariance, respectively, while \(\hat{\mathbf{h}}[n-1\left|n-1\right.]\) and \(\mathbf{\Sigma_{h}}[n-1\left|n-1\right.]\) represent the _a posterior_ channel estimate of the \((n-1)\)-th block obtained by utilizing the observation of the \((n-1)\)-th block and the corresponding _a posterior_ error covariance, respectively. To guarantee the accuracy of the prediction, \(\hat{\mathbf{h}}[1]\) and \(\mathbf{\Sigma_{h}}[1]\) acquired via (6) and (7) are used for the initialization of \(\hat{\mathbf{h}}[1\left|1\right.]\) and \(\mathbf{\Sigma_{h}}[1\left|1\right.]\).

### _Updating stage_

We now refine the one-step prediction with the Kalman updating stage according to the observation equation in (10). Note that the measurement matrix in (10), i.e., \(\mathbf{X}\), is unknown. To remedy this, we propose to reconstruct a measurement matrix using a robust MMSE detector [10], which yields the data estimates as

\[\hat{\mathbf{x}}_{t}[n\left|n-1\right.]=\hat{\mathbf{H}}^{H}[n\left|n-1 \right|]\Big{(}\hat{\mathbf{H}}[n\left|n-1\right.]\hat{\mathbf{H}}^{H}[n\left| n-1\right.] \tag{41}\] \[+\mathbf{\Sigma_{H}}[n\left|n-1\right.]+\sigma^{2}\mathbf{I}_{M} \Big{)}^{-1}\mathbf{y}_{t}[n],\]

where \(\hat{\mathbf{H}}[n\left|n-1\right.]\) is the matrix form of \(\hat{\mathbf{h}}[n\left|n-1\right.]\) and \(\mathbf{\Sigma_{H}}[n\left|n-1\right.]\) is the row covariance of the channel error \(\mathbf{\Delta H}[n\left|n-1\right.]=\mathbf{H}[n\left|n-1\right.]-\hat{ \mathbf{H}}[n\left|n-1\right.]\), which can be derived from \(\mathbf{\Sigma_{h}}[n\left|n-1\right.]\) by (30). Based on (41), the row covariance of the data estimation error \(\Delta\mathbf{x}_{t}[n\left|n-1\right.]=\mathbf{x}_{t}[n\left|n-1\right.]- \hat{\mathbf{x}}_{t}[n\left|n-1\right.]\) can be expressed by

\[\mathbf{\Sigma_{x_{t}}}[n\left|n-1\right.]=\mathbf{I}_{K}-\hat{ \mathbf{H}}^{H}[n\left|n-1\right.]\Big{(}\hat{\mathbf{H}}[n\left|n-1\right.] \hat{\mathbf{H}}^{H}[n\left|n-1\right.] \tag{42}\] \[+\mathbf{\Sigma_{H}}[n\left|n-1\right.]+\sigma^{2}\mathbf{I}_{M} \Big{)}^{-1}\mathbf{\hat{H}}[n\left|n-1\right.],\]

where \(\mathbf{\Sigma_{x_{t}}}[n\left|n-1\right.]\) is independent of the index \(t\) and thus can be simply denoted by \(\mathbf{\Sigma_{x}}[n\left|n-1\right.]\).

Using the estimate \(\hat{\mathbf{x}}_{t}[n\left|n-1\right.]\), the observation equation in (10) can be rewritten as

\[\mathbf{y}_{t}[n]=\hat{\mathbf{X}}_{t}[n\left|n-1\right.]\mathbf{h}[n]+\tilde{ \mathbf{z}}_{t}[n], \tag{43}\]

where the measurement matrix \(\hat{\mathbf{X}}_{t}[n\left|n-1\right.]=\hat{\mathbf{x}}_{t}^{T}[n\left|n-1 \right.]\otimes\mathbf{I}_{M}\), and the equivalent noise \(\tilde{\mathbf{z}}_{t}[n]=\Delta\mathbf{X}_{t}^{T}[n\left|n-1\right.]\mathbf{h }[n]+\mathbf{z}_{t}[n]=\left(\Delta\mathbf{x}_{t}^{T}[n\left|n-1\right.] \otimes\mathbf{I}_{M}\right)\mathbf{h}[n]+\mathbf{z}_{t}[n]\) is assumed to be zero-mean Gaussian distributed with the covariance given by

\[\mathbf{C_{\tilde{\mathbf{x}}_{t}}}[n] =\mathbb{E}\left\{\tilde{\mathbf{z}}_{t}[n]\tilde{\mathbf{z}}_{t }^{H}[n]\right\} \tag{44}\] \[\overset{(a)}{=}\mathbb{E}\left\{\Delta\mathbf{X}_{t}^{T}[n\left|n -1\right.]\mathbf{C_{h}}\left(\Delta\mathbf{X}_{t}^{T}[n\left|n-1\right.] \right)^{H}\right\}+\sigma^{2}\mathbf{I}_{M}\] \[\triangleq\mathbf{\Xi_{t}}[n\left|n-1\right.]+\sigma^{2}\mathbf{ I}_{M},\]

where (a) is obtained by taking the expectation w.r.t. \(\mathbf{h}[n]\) and \(\mathbf{z}_{t}[n]\) simultaneously. Based on the definition of \(\mathbf{\Sigma_{x}}[n\left|n-1\right.]=\mathbb{E}\left\{\Delta\mathbf{x}_{t}[n \left|n-1\right.]\Delta\mathbf{x}_{t}^{H}[n\left|n-1\right.]\right\}\) and the property of Kronecker product, the \((i,j)\)-th entry of \(\mathbf{\Xi_{t}}[n\left|n-1\right.]\) can be derived as (45) at the top of the next page, where \(\mathbf{0}_{m_{1}\times m_{2}}^{\left\backslash i,j\right.}\) denotes the all-zeros matrix of size \(m_{1}\times m_{2}\) except the \((i,j)\)-th entry being \(1\). Similar to \(\mathbf{\Sigma_{x}}[n\left|n-1\right.]\), the index \(t\) of \(\mathbf{C_{\tilde{\mathbf{z}}_{t}}}[n]\) and \(\mathbf{\Xi_{t}}[n\left|n-1\right.]\) can also be omitted since (44) and (45) are both independent of \(t\).

Fig. 3: Block diagram of LCRADMMNet.
Traditional Kalman updating can be readily performed by stacking the \(L\) observation equations over a block into one \(ML\)-dimensional equation. However, the computational complexity to calculate the Kalman gain is up to \(\mathcal{O}\left(M^{3}L^{3}\right)\) per data block. Note that the variables at different time slots can be considered uncorrelated with each other. Therefore we can update the estimates recursively for \(t=1,2,\cdots,L\) based on the sequential filter method [41] via (46)-(48), where \(\mathbf{K}\) is the Kalman gain matrix and \(\left(\cdot\right)^{(t)}\) stands for the value of the input variable updated at time slot \(t\). \(\mathbf{\hat{h}}^{(0)}[n\left|n\right|]\) and \(\mathbf{\Sigma}_{\mathbf{h}}^{(0)}[n\left|n\right|]\) are set to the one-step prediction \(\mathbf{\hat{h}}[n\left|n-1\right|]\) and \(\mathbf{\Sigma}_{\mathbf{h}}[n\left|n-1\right|]\), respectively. After completing the updating process along the \(L\) time slots over a block, we obtain \(\mathbf{\hat{h}}^{(L)}[n\left|n\right|]\) and \(\mathbf{\Sigma}_{\mathbf{h}}^{(L)}[n\left|n\right|]\), which are then regarded as the _a posterior_ channel estimate and the corresponding error covariance, i.e, \(\mathbf{\hat{h}}[n\left|n\right|]\) and \(\mathbf{\Sigma}_{\mathbf{h}}[n\left|n\right|]\). It can be seen that the high-dimensional Kalman updating operation has been equivalently realized in a sequential manner, resulting in a lower complexity of \(\mathcal{O}\left(M^{3}L\right)\) per data block. The overall procedure of the proposed channel tracking method is summarized in Algorithm 4. Based on this, we acquire the refined channel estimate and the corresponding error covariance, which can be readily utilized by RADMMNet and LCRADMMNet to improve the robust detection performance. The block diagram of the RDAKF-based receiver for the \(n\)-th block is shown in Fig. 4.

```
0:\(\mathbf{Y}_{P}\), \(\mathbf{S}_{P}\), \(\mathbf{A}\), \(\mathbf{\bar{A}}\), \(\mathbf{C}_{\mathbf{h}}\), and \(\mathbf{y}_{t}[n]\) for \(t=1,2,\cdots,L\), \(n=2,\cdots,N+1\).
0:\(\mathbf{\hat{x}}\).
1: Initialize \(\mathbf{\hat{h}}[1\left|1\right|]\) and \(\mathbf{\Sigma}_{\mathbf{h}}[1\left|1\right|\) via (6) and (7), respectively, \(n=1\).
2:for\(n=2,\cdots,N+1\)do
3: Predict \(\mathbf{\hat{h}}[n\left|n-1\right|]\) and \(\mathbf{\Sigma}_{\mathbf{h}}[n\left|n-1\right|]\) via (39) and (40), respectively.
4: Calculate \(\mathbf{\hat{x}}_{t}[n\left|n-1\right|]\) for \(t=1,2,\cdots,L\) via (41).
5: Calculate \(\mathbf{C}_{\mathbf{h}}[n]\) via (44).
6:for\(t=1,2,\cdots,L\)do
7: Update \(\mathbf{K}^{(t)}[n]\), \(\mathbf{\hat{h}}^{(t)}[n\left|n\right|]\) and \(\mathbf{\Sigma}_{\mathbf{h}}^{(t)}[n\left|n\right|]\) via (46)-(48), respectively.
8:endfor
9:\(\mathbf{\hat{h}}[n\left|n\right|]\leftarrow\mathbf{\hat{h}}^{(L)}[n\left|n\right|]\), \(\mathbf{\Sigma}_{\mathbf{h}}[n\left|n\right|]\leftarrow\mathbf{\Sigma}_{ \mathbf{h}}^{(L)}[n\left|n\right|]\).
10:endfor
11:\(\left\{\mathbf{\hat{h}}[n\left|n\right|]\right\}_{n=1}^{N+1}\) and \(\left\{\mathbf{\Sigma}_{\mathbf{h}}[n\left|n\right|]\right\}_{n=1}^{N+1}\).
```

**Algorithm 4** RDAKF-Based Channel Tracking Algorithm

## VI Numerical Results

An uplink multiuser MIMO system is considered, where the number of the receive antennas at the base station and the number of the served users are set to \(M=8\) and \(K=4\), respectively, unless otherwise specified. The spatially correlated channels between the users and the base station are described by the Kronecker model [45], i.e., \(\mathbf{H}=\mathbf{R}_{r}^{1/2}\mathbf{H}_{\mathrm{i.i.d}}\mathbf{R}_{t}^{1/2}\), where \(\mathbf{H}_{\mathrm{i.i.d}}\) is the i.i.d. Rayleigh fading channel matrix, each entry of which follows a Gaussian distribution with zero mean and unit variance, \(\mathbf{R}_{t}\) and \(\mathbf{R}_{r}\) are the exponential correlation matrices of the users and the receiver, which can be characterized by the correlation coefficients \(\rho_{t}\) and \(\rho_{r}\), respectively. In the simulation, we set \(\rho_{t}=0.3\), \(\rho_{r}=0.7\), and the length of a coherence block \(L=10\) with \(L_{P}=K\). Without loss of generality, the channel temporal correlation coefficients of different users are assumed to be the same, i.e., \(\rho_{1}=\rho_{2}=\cdots=\rho_{K}=\rho\), where a layer \(\rho\) means a slower channel variation. The signal-to-noise ratio (SNR) of the system is defined as \(\mathrm{SNR}=\frac{K}{\sigma^{2}}\).

Regarding the network training, both RADMMNet and LCRADMMNet are trained offline in an end-to-end manner with 10000 training samples and 2000 validation samples for each data block under each SNR point. We propose to employ a two-step training strategy, which includes a pre-training step and a fine-tuning step. The Adam optimizer is applied for the parameter learning with the initial learning rates of 0.01 for pre-training and 0.001 for fine-tuning. The batchsize for the two steps are set to 500 and 200, respectively, and the number of epochs are both set to 200. To prevent the overfitting and save the training time, we adopt the early-stopping mechanism, which will be triggered if the loss on the validation dataset does not decrease in successive 5 epochs. In order to take the advantage of the supervised learning, we choose the MSE as

Fig. 4: Block diagram of the proposed RDAKF-based receiver.
the loss function for both networks, which is given by

\[\mathcal{L}=\frac{1}{|\mathcal{S}|}\sum_{(\mathbf{x}^{\prime}\text{ }\text{\tiny{mm}},\mathbf{x})\in\mathcal{S}}\big{\|}\mathbf{x}^{I_{\text{mm}}}- \mathbf{x}\big{\|}_{2}^{2}, \tag{49}\]

where \(\mathbf{x}^{I_{\text{mm}}}\) and \(\mathbf{x}\) denote the network output and the true transmitted vector (also known as the label), respectively, and \(\mathcal{S}\) is the training dataset. By minimizing the loss function, the deep unfolding networks are prompted to learn from the data in a supervised fashion, while maintaining the inherent mechanism of the original iterative algorithm. During the online simulation, we generate random testing samples continuously until 1000 symbol errors are collected so as to obtain more reliable SER results. On the other hand, the values of \(L\), \(N\), and other parameters used in the testing stage, are set to be identical as those in the training stage, which is reasonable since the system parameters such as the frame structure and the number of antennas are generally fixed in a given scenario.

Apart from the proposed networks RADMMNet and LCRADMMNet, we also present the performances of the following baselines for comparisons: 1) the linear MMSE detector that directly uses imperfect CSI (Mismatched MMSE), 2) the ML detector that directly uses imperfect CSI in (9) (Mismatched ML), 3) The robust linear MMSE detector proposed in [10] (Robust MMSE), 4) The ADMM-based detection network proposed in [34] (ADMM-PSNet), 5) the inexact ADMM-based detection network proposed in [31] (ADMMNet), 6) the orthogonal approximate message passing-based deep unfolding detection network proposed in [16] (OAMPNet), 7) the robust orthogonal approximate message passing algorithm (ROAMPP) and its corresponding deep unfolding detection network (ROAMPNet) [15]2, and 8) The robust ML detector in (12) (Robust ML).

Footnote 2: Note that ROAMP and ROAMPNet refer to OAMP and OAMPNet2 in [15], respectively. We rename them to emphasize that they belong to the robust designs and clearly distinguish them from the original algorithms OAMP and OAMPNet [16], which are designed under the perfect CSI assumption.

**Remark 1**: _Note that ADMM-PSNet, ADMMNet, and OAMPNet are all designed under the perfect CSI assumption but trained with imperfect CSI, while the proposed RADMMNet and LCRADMMNet are specially designed for the imperfect CSI case by explicitly incorporating the statistics of the imperfect CSI into the network structure. Consequently, both the proposed networks can exhibit stronger robustness in the presence of CSI errors. As for the robust detection network ROAMPNet, it regards the sum of the CSI error and the receiver noise \(\tilde{\mathbf{r}}=\mathbf{\Delta Hx}+\mathbf{z}\) as an equivalent colored Gaussian noise, which is similar to the design principle of LCRADMMNet. However, LCRADMMNet requires a much lower complexity than ROAMPNet while with a comparable or even better performance. On the other hand, different from the above two robust detection networks that are developed based on a suboptimal robust ML metric using the Gaussian approximation, RADMMNet is established according to the optimal robust ML metric, which exploits not only the model-driven deep unfolding technique but also the introduced data-driven structure. Therefore, it can achieve better performance than LCRADMMNet and ROAMPNet. The advantages of the proposed networks over the existing DL structures will be validated by the following numerical results._

### _Convergence Property_

Different from iterative algorithms where the number of iterations can be adjusted adaptively, the numbers of layers for deep unfolding-based networks needs to be pre-determined. For this purpose, we first investigate the convergence property of the proposed RADMMNet and LCRADMMNet. Fig. 5 plots the symbol error rate (SER) performances of the proposed networks versus the number of layers under different modulation orders, where the corresponding two ADMM-based algorithms, denoted by "RADMM" and "LCRADMM", are also provided for comparison. Note that in the model considered in our work, the channel variation becomes more serious as the data index increases or \(\rho\) gets smaller. Meanwhile, it is intuitive that higher-order modulated signals are more vulnerable to the channel variation, which means that a milder condition is required in the 16QAM case to guarantee an acceptable performance. Therefore, for the QPSK case, we take the fifth data block with SNR = 15 dB and \(\rho=0.99\) for illustration, while for the 16-QAM case, we consider the third data block with SNR = 30 dB and \(\rho=0.995\). It

Fig. 5: SER performances versus the number of layers under different modulation orders.
can be observed that both RADMMNet and LCRADMMNet can considerably outperform the corresponding ADMM-based algorithms with faster convergence rates, i.e., within 10 layers, indicating the effectiveness of the designed network and the advantage of the parameters learned via offline training. Meanwhile, we would like to note that, for the same number of layers/iterations, LCRADMMNet requires nearly the same level of complexity as LCRADMM, while RADMMNet can efficiently reduce the complexity of RADMM due to the introduced simplifications. Consequently, in the following simulation, we set the number of layers of both RADMMNet and LCRADMMNet to 10 and do not present the results of RADMMM and LCRADMM for the simplicity of the curves. We also mention that the number of CG iterations involved in the first layer of RADMMNet, i.e., \(I_{\text{CG}}\), is fixed to 15 empirically, which can reach almost the same SER performance as that with an exact matrix inversion.

### _Performance Evaluation_

In this subsection, we compare the performances of different detectors using the LMMSE channel estimation presented in Section II and the RDAKF-based channel tracking presented in Section IV, respectively. Then, the performances with larger MIMO systems are also provided.

#### Iv-B1 LMMSE channel estimation

We first consider the case that the imperfect CSI is acquired by the LMMSE channel estimation. Fig. 6(a) shows the SER performances of different detection algorithms versus the index of data blocks for SNR = 15 dB and \(\rho=0.99\) under QPSK modulation. Since the channel estimates of all data blocks are obtained using the pilot in the first block, the channel errors will be accumulated due to the time variation of channels. As a result, the performances of all the detection algorithms degrade as the index increases, among which even the mismatched ML detector suffers from an unbearable performance loss. In contrast, the robust detectors can effectively compensate the deterioration caused by the channel aging and achieve satisfying performance, validating the necessity of the robust design. In particular, the proposed RADMMNet yields a relatively low SER and approaches the optimal robust ML detector. LCRADMMNet performs a little worse than RADMMNet while with lower complexity (see Section VII-C). With the significant performance advantage of RADMMNet and LCRADMMNet over the non-robust detectors, more data blocks can be transmitted with only one pilot block, which can improve the spectrum efficiency.

The SER performances versus SNR for the fifth data block are illustrated in Fig. 6(b). As can be seen, the detectors that neglect the imperfection of CSI suffer from high error floors in the high-SNR region. We notice that both ADMM-PSNet and ADMMNet show slightly better performances than the mismatched ML detector. The reason is that they are trained to approach the true transmitted vectors using the imperfect CSI via supervised learning. In this way, the embedded statistical information of CSI errors can be implicitly utilized, which consequently provides the networks partial robustness against imperfect CSI. On the other hand, by fully incorporating the statistical information of the imperfect CSI into the network design, the proposed RADMMNet exhibits a much better performance. In addition, the lower-complexity LCRADMMNet can also considerably outperform other baselines including the

Fig. 6: SER performances with LMMSE channel estimation for the scenario of (a)(b) \(M=8\), \(K=4\), and \(\rho=0.99\) under QPSK modulation, (c)(d) \(M=8\), \(K=4\), and \(\rho=0.995\) under 16QAM modulation, and (e)(f) \(M=8\), \(K=8\), and \(\rho=0.995\) under QPSK modulation.
existed non-robust detectors and the robust ones.

Next we present the SER performances under 16-QAM constellation with \(\rho=0.995\), which corresponds to a lower-mobility scenario. Figs. 6(c) and 6(d) present the SER performances versus the index of data blocks for SNR = 30 dB and those versus SNR for the third data block, respectively. It can be observed that, even with such a slower channel variation, the increasing index will lead to a serious performance degradation for each detection algorithm. Hence, in the 16-QAM case, the pilot block can only be followed by a few data blocks in a frame to guarantee the performance. Nevertheless, the proposed robust detection networks still perform much better than the non-robust detectors. Besides, as shown in Fig. 6(d), compared with the robust MMSE detector, noticable performance gains up to 12 dB and 15 dB in the high-SNR region can be achieved by RADMMNet and LCRADMMNet, respectively. Although the performances of LCRADMMNet and ROAMPNet are close, the complexity of LCRADMMNet is much lower than that of ROAMPNet, which will be analyzed in detail in Section VII-C.

To verify the advantages of our proposed networks in the challenging scenario when the ratio of the number of receive antennas to that of transmit antennas equals \(1\), we conduct simulations by setting \(M=K=8\) and \(\rho=0.995\) under the QPSK modulation. The SER performances versus the index of data blocks for SNR = 20 dB and those versus SNR for the third data block are plotted in Figs. 6(e) and 6(f), respectively, where the mismatched ML detector and the robust ML detector are not incorporated due to their high complexities. Compared with the previous results with \(M=8\) and \(K=4\) in Figs. 6(a) and 6(b), the performance gains of the proposed RADMMNet and LCRADMMNet over other baselines, especially the mismatched MMSE detector and the robust MMSE detector, are more distinct. We also note that, for the high SNR region where the CSI error is the dominant factor of deteriorating the performance, the non-robust detectors and detection networks suffer from an error floor or even an abnormally rising SER, while the proposed robust detection networks can still yield satisfying performances.

#### Vii-C2 Rdmaf-based channel tracking

As an improved alternative of the LMMSE channel estimation, the RDAKF-based channel tracking method is evaluated in terms of the NMSE performances of channel estimates and the SER performances, respectively, where the NMSE is defined as \(\frac{\left\|\mathbf{H}-\mathbf{H}\right\|_{\text{F}}^{2}}{\left\|\mathbf{H} \right\|_{\text{F}}^{2}}\). Figs. 7(a) and 7(b) compare the NMSE performances between the two CSI acquisition methods for different indices of data blocks under different modulation orders. As expected, the RDAKF-based method can effectively track the channel variation with the aid of the pre-estimated data symbols in both the QSPK and 16-QAM cases. For the larger index of data blocks, the estimation performance gap between the RDAKF and the LMMSE-based methods becomes more evident. This is due to the fact that the error propagation caused by the coarse prediction of the LMMSE method gets more severe, which can be suppressed in the updating stage of the RDAKF method. We also notice from Fig. 7(a) that, in the low-SNR region, the inaccurate data estimates may deteriorate

Fig. 7: (a)(b) NMSE performances for different indices of data blocks for the scenario of \(M=8\) and \(K=4\) under different modulation orders, and SER performances with RDAKF channel tracking for the scenario of (c)(d) \(M=8\), \(K=4\), and \(\rho=0.99\) under QPSK modulation and (e)(f) \(M=8\), \(K=4\), and \(\rho=0.995\) under 16QAM modulation.
the accuracy of channel estimates, leading to an even poorer NMSE performance of the RDAKF method than that of the LMMSE-based method.

The SER performances versus the index of data blocks are plotted in Fig. 7(c) with the same simulation settings as in Fig. 6(a). Based on the refined CSI, all the detection algorithms enjoy a performance improvement compared with the results in Fig. 6(a). Therefore, a frame structure containing more data blocks can be supported with the RDAKF-based channel tracking method. However, there is almost no advantage of the robust MMSE-based detector over the mismatched MMSE-based one, indicating the limitation of a linear robust design when the channel estimate is relatively accurate. Besides, the performance gaps between non-robust detection networks and LCRADMMNet are narrowed with the reduced channel error. Fig. 7(d) shows the performances versus SNR for the fifth data block. High error floors can be observed for the non-robust detectors as well as the robust MMSE detector, since they are unable to make good use of the relatively small CSI uncertainty. In contrast, RADMMNet still shows a comparable performance to the robust ML detector by fully exploiting the statistics of CSI errors.

Similar trends can be found in the 16-QAM case as shown in Figs. 7(e) and 7(f), where the simulation settings are the same as in Figs. 6(c) and 6(d), respectively. The proposed RADMMNet can maintain the superiority over other baselines, demonstrating its stronger robustness against imperfect CSI. However, the SER performance gains obtained by utilizing the RDAKF-based channel tracking method are less obvious than those in the QPSK case, albeit with the similar level of improvement in the NMSE performance as illustrated in Fig. 7(b). This is because the detection of 16-QAM signals requires more accurate channel estimates. In addition, the performance gap from the robust ML detector is more evident than in the QPSK case, since 16-QAM signals are more susceptible to CSI errors. Hence, we conclude that the proposed RDAKF-based method is more applicable for QPSK modulation. The results for the more challenging scenario of \(M=K=8\), corresponding to the simulation settings in Figs. 6(e) and 6(f), also exhibit a similar behavior with improved performances benefited from the RDAKF method, which are not presented here due to the space limitation.

#### V-B3 Large-scale MIMO

We now investigate the SER performances for MIMO system with larger sizes. Since a larger array can effectively make up the performance loss caused by the imperfect CSI, it can be inferred that a faster channel variation can be well handled with more receive antennas. To validate this, we test the detection algorithms under higher-speed scenarios by setting smaller \(\rho\)'s for larger-scale MIMO systems, i.e., \(\rho=0.98\) for \(M=32\) and \(K=8\), and \(\rho=0.95\) for \(M=128\) and \(K=16\) in Figs. 8(a) and 8(b), respectively. The results of the mismatched ML detector, the robust ML detector, and RADMMNet are not presented due to their unaffordable complexities. It can be observed that the proposed LCRADMMNet clearly outperforms most benchmark schemes and even the MMSE detector with larger \(\rho\)'s. Additionally, LCRADMMNet can yield a slightly better performance than ROAMPNet while with a much lower complexity, which will be illustrated in the next subsection. On the other hand, the performance gains brought by the RDAKF-based channel tracking are conspicuous, especially for the scenario of \(M=128\) and \(K=16\), thanks to the large antenna array gain that results in a high probability of obtaining correct pre-estimated symbols.

### _Complexity Analysis_

The computational complexities of different detectors are analyzed in Table II. For the proposed RADMMNet, the major computations lie in the CG iterations and the matrix inversion3, which leads to a complexity of \(\mathcal{O}\left(\left(I_{\mathrm{CG}}+I_{\mathrm{max}}-1\right)M^{2}K^{2}+I_{ \mathrm{max}}K^{3}\right)\), while the other proposed network LCRADMMNet only needs to calculate one \(M\)-dimensional matrix inversion, one \(K\)-dimensional matrix inversion, and a few matrix multiplications with a total complexity of \(\mathcal{O}\left(M^{3}+K^{3}+I_{\mathrm{max}}K^{2}\right)\). The robust detection network ROAMPNet, as one of the baselines, requires a higher complexity of \(\mathcal{O}\left(I_{\mathrm{max}}M^{3}\right)\). On the other hand, the

Fig. 8: SER performances versus SNR for the fifth data block under QPSK modulation in large-scale MIMO scenarios.
complexities of the mismatched ML detector and the robust ML detector are \(\mathcal{O}\left(4^{QK}K^{3}\right)\) and \(\mathcal{O}\left(4^{QK}M^{3}K^{3}\right)\), respectively, which can be prohibitively high for a high-modulation order or a large number of antennas. The average CPU time requirements of different detectors under different scenarios are also provided for a more intuitive comparison, which show a consistent trend as our analysis. In particular, the low-complexity design LCRADMMNet requires nearly an order of magnitude less CPU time than RADMMNet. Concerning the SER performances presented before, although the non-robust detection networks may require less complexity than the proposed networks, they suffer from serious performance degradation. Meanwhile, it can be found that RADMMNet can achieve comparable performance to the optimal robust ML detection within only a few layers, while LCRADMMNet can always outperform ROAMPNet with much less complexity. Hence, we conclude that both the proposed networks achieve attractive tradeoffs between the performance and the complexity. Regarding the considered two channel acquisition methods, the complexity of LMMSE channel estimation is \(\mathcal{O}\left(M^{3}K^{3}\right)\), while the proposed RDAKF-based channel tracking method with much better NMSE performance additionally requires a complexity of \(\mathcal{O}\left(M^{3}L\right)\) per data block, due to the matrix inversion involved in the calculation of the Kalman gain.

## VII Conclusion

This paper investigated the statistically robust detector design for MIMO systems by taking into account both channel estimation error and channel variation. We first derived an ADMM-based robust detection algorithm, which admits the calculations of closed-form expressions in each iteration. Then, by deep unfolding the ADMM iterations and introducing some data-driven structures, we advocated a robust detection network RADMMNet with trainable parameters. Furthermore, by adopting a Gaussian approximation for the CSI error, a low-complexity robust MIMO detector was further developed along with the corresponding deep unfolding network LCRADMMNet. In addition, as a complementary way to combat the channel variation and enhance the detection performance, we also presented a Kalman-filter-based channel tracking method by fully exploiting the pre-estimated data symbols. Simulation results confirmed that, the proposed two networks can considerably outperform the non-robust detectors and even approach the optimal robust ML detector with much lower complexities under both the LMMSE channel estimation and the proposed RDAKF-based channel tracking method.
