Message Passing Based Two-Timescale Bayesian Learning for Joint Channel and Memory Hardware Impairments Tracking

# Message Passing Based Two-Timescale Bayesian Learning for Joint Channel and Memory Hardware Impairments Tracking
Source: [https://arxiv.org/html/2607.01660](https://arxiv.org/html/2607.01660)
Wei Xu,Graduate Student Member, IEEE,and An Liu,Senior Member, IEEEThis work was supported in part by the National Key Research and Development Program of China under Grant 2025ZD1301900; in part by National Key Laboratory of Millimeter\-Wave and Terahertz Remote Sensing and the Zhejiang Provincial Key Laboratory of Information Processing, Communication and Networking \(IPCAN\), Hangzhou, China; in part by the Zhejiang Provincial Key Laboratory of Multi\-Modal Communication Networks and Intelligent Information Processing, Hangzhou, China\. \(Corresponding authors: An Liu\.\)Wei Xu and An Liu are with the College of Information Science and Electronic Engineering, Zhejiang University, Hangzhou 310027, China \(e\-mails: 12231077@zju\.edu\.cn and anliu@zju\.edu\.cn\)\.

###### Abstract

Hardware impairments in massive multiple\-input multiple\-output \(MIMO\) receivers introduce inter\-symbol memory and inter\-element coupling, severely degrading channel estimation\. This paper employs a residual recurrent gated unit \(RGRU\) to model the intra\-slot memory of the hardware impairments and proposes a message\-passing\-based two\-timescale Bayesian deep learning \(MP\-TTBDL\) framework for joint channel and impairment tracking\. Owing to small\-scale fading, the wireless channel varies rapidly across slots, whereas hardware impairments drift slowly due to hardware aging and environmental variations\. To capture these distinct physical timescales, a fast\-varying Markov prior and a slow\-varying Gaussian Markov prior are assigned to the sparse channel and the network parameters, respectively\. Based on a multi\-slot factor graph formulation, a message\-passing algorithm is developed\. Specifically, the inter\-slot messages admit closed\-form updates, while the intra\-slot factor graph, due to its complex recurrent structure, is partitioned into a channel tracking module and an impairments calibration module\. The channel tracking module performs sparse channel estimation via turbo orthogonal approximate message passing \(Turbo\-OAMP\), and the impairments calibration module updates the impairment parameters via a specially designed deep approximate message passing \(DAMP\) procedure, with the two modules iteratively exchanging extrinsic information through expectation propagation \(EP\) until convergence\. Simulation results show that the proposed framework robustly achieves lower channel estimation error than conventional compensators followed by channel estimation across different online impairment scenarios and signal\-to\-noise ratio \(SNR\) conditions\.

## IIntroduction

Massive MIMO has become a cornerstone of modern wireless communications, with base stations \(BSs\) deploying large antenna arrays to exploit spatial diversity and multiplexing gains\. To realize these gains, accurate channel estimation and beamforming are indispensable\. However, as the number of antennas grows, RF front\-end hardware impairments including LNA nonlinearity\[[20](https://arxiv.org/html/2607.01660#bib.bib83),[30](https://arxiv.org/html/2607.01660#bib.bib84)\], crosstalk between RF chains\[[4](https://arxiv.org/html/2607.01660#bib.bib7),[3](https://arxiv.org/html/2607.01660#bib.bib6)\], and IQ imbalance\[[27](https://arxiv.org/html/2607.01660#bib.bib80),[1](https://arxiv.org/html/2607.01660#bib.bib81)\]become increasingly severe, introducing unknown phase and amplitude errors, distorting pilot signals, and severely degrading channel estimation performance\. These impairments have gradually become a critical bottleneck in system design\.

To mitigate these impairments, conventional methods typically train an offline compensator to learn the inverse mapping of the hardware impairments, and then apply the learned model with fixed parameters during online reception\[[36](https://arxiv.org/html/2607.01660#bib.bib75),[15](https://arxiv.org/html/2607.01660#bib.bib18),[24](https://arxiv.org/html/2607.01660#bib.bib54),[31](https://arxiv.org/html/2607.01660#bib.bib65)\]\. Since practical RF front\-end hardware exhibits both nonlinearity and memory effects across symbols\[[2](https://arxiv.org/html/2607.01660#bib.bib4),[22](https://arxiv.org/html/2607.01660#bib.bib52)\], the Volterra series has been widely used to characterize these effects and their inverses\[[21](https://arxiv.org/html/2607.01660#bib.bib50),[13](https://arxiv.org/html/2607.01660#bib.bib82),[9](https://arxiv.org/html/2607.01660#bib.bib8)\]\. However, the number of Volterra parameters grows combinatorially with the nonlinearity order, memory depth, and the number of RF chains, making it prohibitive for large\-scale systems\. As a practical alternative, simplified models such as the generalized memory polynomial \(GMP\)\[[21](https://arxiv.org/html/2607.01660#bib.bib50)\]are often adopted for compensation\. Yet, the exact inverse mapping may not be accurately representable by polynomials, and GMP\-based compensators can therefore suffer from significant model mismatch\. Recently, recurrent neural networks such as long short\-term memory \(LSTM\) and GRU\[[5](https://arxiv.org/html/2607.01660#bib.bib14),[16](https://arxiv.org/html/2607.01660#bib.bib79),[11](https://arxiv.org/html/2607.01660#bib.bib33)\]have been employed to compensate for the combined impairments, achieving superior performance over GMP in many digital predistortion tasks thanks to their universal approximation capability\.

Existing compensation approaches suffer from two fundamental issues\. First, in the presence of inter\-symbol memory and inter\-chain crosstalk, each observation depends on multiple past transmitted symbols, and learning the inverse mapping that recovers the original signals from a mixed and overlapping set of contributions becomes an ill\-posed problem, leading to poor compensation and limited generalization\. In contrast, the forward hardware\-impairment model possesses a natural causal and recursive structure; it describes how hardware distorts original signals symbol by symbol, a structure that is far easier to approximate with a recurrent network\. With an accurate forward model, signal recovery can then be cast as a Bayesian inference problem together with priors, yielding more reliable results\. Second, conventional methods typically train the compensator offline, and the resulting static model may fail to track the actual hardware behavior during online operation due to temperature variations\[[26](https://arxiv.org/html/2607.01660#bib.bib61)\], component aging, and other environmental factors\[[37](https://arxiv.org/html/2607.01660#bib.bib76)\]\. Online adaptation of the impairment\-aware receiver can mitigate such drift by learning from live data\[[14](https://arxiv.org/html/2607.01660#bib.bib38)\], but most existing pipelines perform point estimation without an adaptive mechanism to regulate the update step size, and typically require buffering additional calibration samples and executing a separate retraining stage after coarse detection, which introduces substantial memory and computational overhead beyond the main receiver processing\.

Channel estimation error is critical to wireless communication system performance\. In recent years, Bayesian inference has been widely applied to various channel estimation scenarios and has achieved superior estimation performance\. For example, subspace\-constrained variational Bayesian inference \(SC\-VBI\) has been proposed for structured compressive sensing with a dynamic grid by replacing high\-dimensional matrix inversions with low\-dimensional subspace\-constrained updates\[[38](https://arxiv.org/html/2607.01660#bib.bib86)\], an alternating maximum a posteriori \(MAP\) framework jointly detects visibility regions and estimates channels in near\-field spatially non\-stationary extra\-large MIMO \(XL\-MIMO\) systems\[[35](https://arxiv.org/html/2607.01660#bib.bib87)\], and a Turbo\-CS\-based sparse\-Markov scheme enables robust multi\-user uplink channel tracking in 5G New Radio \(NR\) under hopping sounding reference signal \(SRS\) patterns and system imperfections\[[29](https://arxiv.org/html/2607.01660#bib.bib88)\]\. However, few existing works consider channel estimation under non\-ideal transceivers with hardware impairments\.

Message\-passing\-based Bayesian inference has also been widely applied to symbol detection, where embedding neural impairment models in a factor graph yields significant performance gains when the channel is assumed known\[[7](https://arxiv.org/html/2607.01660#bib.bib16),[6](https://arxiv.org/html/2607.01660#bib.bib17)\], validating the benefit of principled inference over a compensator\. Nevertheless, such approaches typically neglect the memory effects of hardware impairments, do not account for time\-varying impairments during online operation, and do not consider multipath channel estimation in practical wireless communication scenarios, leaving the additional gains from jointly inferring the sparse multipath channel and impairment parameters in a unified factor graph unexploited\. Therefore, it is imperative to establish a joint Bayesian estimation framework that accounts for both the wireless channel and hardware impairments, whereas directly applying conventional compensation and then estimating the channel from the compensated signal inevitably leads to significant degradation in channel estimation performance\.

Building upon our prior work that considered the joint tracking problem under memoryless impairments\[[34](https://arxiv.org/html/2607.01660#bib.bib72)\], this paper further investigates the uplink channel estimation problem for a massive MIMO receiver with memory hardware impairments\. We use an RGRU to model the forward impairment mapping, and integrate sparse channel estimation with RGRU\-based impairment learning in a unified factor graph\. Then we propose the MP\-TTBDL framework that jointly tracks the fast\-varying channel and the slow\-varying impairments\.

The main contributions of this paper are summarized as follows:

- •A Unified Bayesian Framework for Joint Channel and Memory Impairment Estimation:We employ a residual RGRU with probit activation functions to model the memory hardware impairments at the receiver\. By embedding this RGRU as the observation likelihood in a unified factor graph together with the sparse channel model, we formulate the problem as a closed\-loop joint Bayesian inference problem in which the channel and impairment parameters are estimated jointly from distorted pilots\. By exploiting the natural causal and recursive structure of the forward impairment process, this framework avoids the ill\-posed inverse mapping learned by conventional compensators and removes the need for a separately trained inverse compensator\.
- •Two\-Timescale Markov Transition Models for Joint Tracking:The wireless channel varies rapidly across slots due to small\-scale fading, whereas hardware impairments drift slowly owing to temperature variations and component aging\. To reflect these distinct physical timescales, we design the transition probabilities of the sparse angular\-domain channel as a fast\-varying Markov chain, and those of the RGRU impairment parameters as a slow\-varying Gaussian Markov chain\. This two\-timescale Markov transition model design aligns the statistical model with the practical transmission system and allows the receiver to continuously track both channel variations and hardware aging in a principled manner\.
- •The Proposed MP\-TTBDL Algorithm:Based on the proposed Bayesian tracking framework, we develop the MP\-TTBDL algorithm that combines Turbo\-OAMP for sparse channel tracking with a novel DAMP scheme tailored to the RGRU architecture\. Given the complex factor graph structure of the RGRU, we design a dedicated DAMP procedure for message updating that specifies a coordinated message update schedule and closed\-form approximation methods\. Within each slot, the two modules iteratively exchange extrinsic information via EP, and the updated messages are propagated across slots under the two\-timescale Markov priors, enabling continuous online joint tracking\.
- •Simulation Validation:Extensive simulations are conducted to validate the proposed contributions\. On synthetically generated impairment data, the dedicated DAMP converges faster and achieves higher fitting accuracy than conventional optimizers for learning the RGRU parameters\. When the online impairments match the offline pretraining configuration, the proposed joint Bayesian framework yields significantly lower channel estimation normalized mean squared error \(NMSE\) than compensation\-based baselines across a wide signal\-to\-noise ratio \(SNR\) range, validating the advantage of the joint Bayesian framework\. Under slowly time\-varying impairments, the MP\-TTBDL algorithm maintains robust channel tracking performance, confirming the effectiveness of the two\-timescale Markov prior and the MP\-TTBDL algorithm\.

The remainder of this paper is organized as follows\. Section II presents the system model, the RGRU\-based impairment surrogate, and the two\-timescale Markov priors\. Section III details the MP\-TTBDL framework, including the factor graph design and the inference algorithm\. Section IV provides simulation results, and Section V concludes the paper\.

Notations:𝟎N\\boldsymbol\{0\}\_\{N\}refers to an all\-zero column vector of dimensionNN\.𝑰N\\boldsymbol\{I\}\_\{N\}refers to anN×NN\\times Nidentity matrix\.\(⋅\)¯\\overline\{\\left\(\\cdot\\right\)\},\(⋅\)T\\left\(\\cdot\\right\)^\{\\mathrm\{T\}\},\(⋅\)H\\left\(\\cdot\\right\)^\{\\mathrm\{H\}\}, and\(⋅\)−1\\left\(\\cdot\\right\)^\{\\mathrm\{\-1\}\}denote conjugate, transpose, conjugate transpose, and inverse, respectively\.ℜ⁡\(⋅\)\\Re\(\\cdot\)denotes the real part of a complex argument\.ℝ\\mathbb\{R\}andℂ\\mathbb\{C\}refer to the sets of real and complex numbers\.\|a\|\|a\|and‖𝒂‖2\\\|\\boldsymbol\{a\}\\\|\_\{2\}denote the amplitude of scalaraaand the 2\-norm of vector𝒂\\boldsymbol\{a\}, respectively\.ana\_\{n\}andAm​nA\_\{mn\}refer to thenn\-th entry of vector𝒂\\boldsymbol\{a\}and the\(m,n\)\(m,n\)\-th entry of matrix𝑨\\boldsymbol\{A\}, respectively\.\[⋅\]i,j\\left\[\\cdot\\right\]\_\{i,j\}denotes the\(i,j\)\(i,j\)\-th entry of a matrix\. The Hadamard \(element\-wise\) product is denoted by⊙\\odot, and the Kronecker product is denoted by⊗\\otimes\.𝒪​\(⋅\)\\mathcal\{O\}\(\\cdot\)denotes the order of complexity\.o​\(⋅\)o\(\\cdot\)denotes higher\-order infinitesimals\.ℛ​\(⋅\)=\[ℜ⁡\(⋅\)ℑ⁡\(⋅\)\]\\mathcal\{R\}\(\\cdot\)=\\begin\{bmatrix\}\\Re\(\\cdot\)\\\\ \\Im\(\\cdot\)\\end\{bmatrix\}facilitates computation in the real domain\.

## IISystem Model

### II\-AIdeal Uplink Pilot Signal Model

In a 5G wideband communication system, a massive MIMO BS equipped withNNantennas serves multiple single\-antenna users\. Without loss of generality, we focus on the uplink channel estimation problem for a single user under flat fading\. The proposed framework admits straightforward extensions to frequency\-selective channels and multi\-user scenarios\. We denote byPPthe number of pilot symbols per slot and byp∈\{1,…,P\}p\\in\\\{1,\\ldots,P\\\}the pilot symbol index\. Then, the received signal𝒚p​\(t\)∈ℂN\\boldsymbol\{y\}\_\{p\}\(t\)\\in\\mathbb\{C\}^\{N\}at symbolppof slotttis expressed as:

𝒚p​\(t\)=𝒉​\(t\)​rp​\(t\)\+𝒏p​\(t\),p=1,…,P,\\boldsymbol\{y\}\_\{p\}\\left\(t\\right\)=\\boldsymbol\{h\}\\left\(t\\right\)r\_\{p\}\\left\(t\\right\)\+\\boldsymbol\{n\}\_\{p\}\\left\(t\\right\),\\quad p=1,\\ldots,P,\(1\)whererp​\(t\)∈ℂr\_\{p\}\(t\)\\in\\mathbb\{C\}is the uplink pilot symbol at symbolpp,𝒉​\(t\)∈ℂN\\boldsymbol\{h\}\(t\)\\in\\mathbb\{C\}^\{N\}is the channel vector assumed constant within slottt, and𝒏p​\(t\)∈ℂN\\boldsymbol\{n\}\_\{p\}\(t\)\\in\\mathbb\{C\}^\{N\}is the additive white Gaussian noise \(AWGN\)\.

Due to sparse multipath propagation, we express the channel𝒉​\(t\)\\boldsymbol\{h\}\(t\)as a linear combination of off\-grid angular\-domain basis vectors\[[17](https://arxiv.org/html/2607.01660#bib.bib40)\]to exploit angular\-domain sparsity for improved channel estimation accuracy:

𝒉=𝑨​\(ϑ\)​𝒙,\\boldsymbol\{h\}=\\boldsymbol\{A\}\\left\(\\boldsymbol\{\\vartheta\}\\right\)\\boldsymbol\{x\},whereϑ∈ℝN\\boldsymbol\{\\vartheta\}\\in\\mathbb\{R\}^\{N\}is a dynamic grid withNNgrid points in angular domain,𝒙∈ℂN\\boldsymbol\{x\}\\in\\mathbb\{C\}^\{N\}is the sparse angular\-domain channel vector, and columns in𝑨​\(ϑ\)\\boldsymbol\{A\}\\left\(\\boldsymbol\{\\vartheta\}\\right\)are steering vectors

𝑨​\(ϑ\)\\displaystyle\\boldsymbol\{A\}\\left\(\\boldsymbol\{\\vartheta\}\\right\)=\[𝒂​\(ϑ1\),…,𝒂​\(ϑN\)\],\\displaystyle=\\left\[\\boldsymbol\{a\}\\left\(\\vartheta\_\{1\}\\right\),\\ldots,\\boldsymbol\{a\}\\left\(\\vartheta\_\{N\}\\right\)\\right\],𝒂​\(ϑ\)\\displaystyle\\boldsymbol\{a\}\\left\(\\vartheta\\right\)≜\[1,ej​π​sin⁡\(ϑ\),…,ej​\(N−1\)​π​sin⁡\(ϑ\)\]T\.\\displaystyle\\triangleq\\left\[1,e^\{j\\pi\\sin\\left\(\\vartheta\\right\)\},\\ldots,e^\{j\\left\(N\-1\\right\)\\pi\\sin\\left\(\\vartheta\\right\)\}\\right\]^\{T\}\.For clarity, we adopt the half\-wavelength spaced uniform linear array \(ULA\) steering\-vector model above\. Extension to other array geometries, such as uniform planar arrays \(UPA\), is straightforward by replacing the steering\-vector dictionary𝑨​\(ϑ\)\\boldsymbol\{A\}\(\\boldsymbol\{\\vartheta\}\)with the corresponding array response, while the overall Bayesian inference and message\-passing framework remains unchanged\. The dynamic angle gridϑ\\boldsymbol\{\\vartheta\}is often initialized using a uniform grid and then updated based on gradient method\.

### II\-BUplink Signal Model with Memory Hardware Impairments

![Refer to caption](https://arxiv.org/html/2607.01660v1/x1.png)Figure 1:Block diagram of a practical massive MIMO uplink receiver with hardware impairments; RF\-chain crosstalk, LNA nonlinearity, and IQ imbalance are shown as representative examples\.In a practical massive MIMO receiver, the signal captured by the antenna array is inevitably distorted by the non\-ideal RF front\-end\. As illustrated in Fig\.[1](https://arxiv.org/html/2607.01660#S2.F1), the hardware impairments can take various forms\. For instance: \(i\) electromagnetic coupling among closely spaced RF chains introduces linear and nonlinear crosstalk; \(ii\) the LNA exhibits nonlinear amplification with memory, where the output depends on current and several past inputs; and \(iii\) IQ imbalance arising from mismatched local oscillator branches results in phase and amplitude distortion of the complex baseband signal\[[24](https://arxiv.org/html/2607.01660#bib.bib54),[36](https://arxiv.org/html/2607.01660#bib.bib75),[28](https://arxiv.org/html/2607.01660#bib.bib63)\]\. Together, these effects impose a nonlinear mapping with inter\-symbol memory and inter\-element coupling on the received pilot, making hardware distortion a severe bottleneck for channel estimation as the antenna count grows\.

In prior digital predistortion works, such memory input\-output behavior is typically modeled by memory polynomials\[[21](https://arxiv.org/html/2607.01660#bib.bib50)\]or recurrent neural networks\[[5](https://arxiv.org/html/2607.01660#bib.bib14),[16](https://arxiv.org/html/2607.01660#bib.bib79),[11](https://arxiv.org/html/2607.01660#bib.bib33)\], where GRU\-type networks in particular achieve robustly low fitting error\. Besides, neural networks are more general in modeling nonlinear mappings and can be used to characterize more complicated hardware impairments in practice\. In such tasks, a single\-layer RGRU already provides sufficient modeling accuracy, as demonstrated in\[[32](https://arxiv.org/html/2607.01660#bib.bib66)\]; deeper RGRUs, due to the large number of extra parameters, tend to overfit the training data and suffer performance degradation during online deployment\. Therefore, in this paper, we only consider a single\-layer RGRU\. For each slottt, let𝝅p​\(t\)∈ℝNπ\\boldsymbol\{\\pi\}\_\{p\}\(t\)\\in\\mathbb\{R\}^\{N\_\{\\pi\}\}denote the hidden state of the RGRU after processing thepp\-th pilot symbol\. Given𝒚p​\(t\)\\boldsymbol\{y\}\_\{p\}\(t\)and𝝅p​\(t\)\\boldsymbol\{\\pi\}\_\{p\}\(t\), the RGRU outputs an estimate ofℛ​\(𝒚~p​\(t\)\)\\mathcal\{R\}\(\\tilde\{\\boldsymbol\{y\}\}\_\{p\}\(t\)\)with AWGN:

ℛ​\(𝒚~p​\(t\)\)=ℛ​\(𝒚p​\(t\)\)\+𝑾y~​\(t\)​𝝅p​\(t\)\+𝒃y~​\(t\)\+𝒏~p,\\mathcal\{R\}\\left\(\\tilde\{\\boldsymbol\{y\}\}\_\{p\}\\left\(t\\right\)\\right\)=\\mathcal\{R\}\\left\(\\boldsymbol\{y\}\_\{p\}\\left\(t\\right\)\\right\)\+\\boldsymbol\{W\}\_\{\\tilde\{y\}\}\\left\(t\\right\)\\boldsymbol\{\\pi\}\_\{p\}\\left\(t\\right\)\+\\boldsymbol\{b\}\_\{\\tilde\{y\}\}\\left\(t\\right\)\+\\tilde\{\\boldsymbol\{n\}\}\_\{p\},\(2\)where we employℛ​\(⋅\)\\mathcal\{R\}\(\\cdot\)to facilitate computing in the real domain, and denote the weights/bias in the output layer as𝑾y~​\(t\)\\boldsymbol\{W\}\_\{\\tilde\{y\}\}\(t\)/𝒃y~​\(t\)\\boldsymbol\{b\}\_\{\\tilde\{y\}\}\(t\), respectively\.𝒏~p∼∏n=12​N𝒩​\(nn;0,σ~2\)\\tilde\{\\boldsymbol\{n\}\}\_\{p\}\\sim\\prod\_\{n=1\}^\{2N\}\\mathcal\{N\}\(n\_\{n\};0,\\tilde\{\\sigma\}^\{2\}\)represents the fitting error\. The hidden state𝝅p​\(t\)\\boldsymbol\{\\pi\}\_\{p\}\(t\)encodes the symbol\-level impairment memory accumulated along the pilot chain within slottt\. It is recurrently updated by a single\-layer GRU:𝒛p​\(t\)\\boldsymbol\{z\}\_\{p\}\(t\)and𝒄p​\(t\)\\boldsymbol\{c\}\_\{p\}\(t\)are the update and reset gates,𝝅ˇp​\(t\)\\check\{\\boldsymbol\{\\pi\}\}\_\{p\}\(t\)filters the past memory, and𝝅~p​\(t\)\\boldsymbol\{\\tilde\{\\pi\}\}\_\{p\}\(t\)extracts a new memory component from the current real\-domain pilotℛ​\(𝒚p​\(t\)\)\\mathcal\{R\}\(\\boldsymbol\{y\}\_\{p\}\(t\)\):

𝒛p​\(t\)=𝑾z​\(t\)​\[𝝅p−1​\(t\);ℛ​\(𝒚p​\(t\)\)\]\+𝒃z​\(t\),\\displaystyle\\boldsymbol\{z\}\_\{p\}\\left\(t\\right\)=\\boldsymbol\{W\}\_\{z\}\\left\(t\\right\)\[\\boldsymbol\{\\pi\}\_\{p\-1\}\\left\(t\\right\);\\mathcal\{R\}\\left\(\\boldsymbol\{y\}\_\{p\}\\left\(t\\right\)\\right\)\]\+\\boldsymbol\{b\}\_\{z\}\\left\(t\\right\),\(3\)𝒄p​\(t\)=𝑾c​\(t\)​\[𝝅p−1​\(t\);ℛ​\(𝒚p​\(t\)\)\]\+𝒃c​\(t\),\\displaystyle\\boldsymbol\{c\}\_\{p\}\\left\(t\\right\)=\\boldsymbol\{W\}\_\{c\}\\left\(t\\right\)\[\\boldsymbol\{\\pi\}\_\{p\-1\}\\left\(t\\right\);\\mathcal\{R\}\\left\(\\boldsymbol\{y\}\_\{p\}\\left\(t\\right\)\\right\)\]\+\\boldsymbol\{b\}\_\{c\}\\left\(t\\right\),\(4\)𝝅ˇp​\(t\)=𝒬​\(𝒄p​\(t\)\)⊙𝝅p−1​\(t\),\\displaystyle\\check\{\\boldsymbol\{\\pi\}\}\_\{p\}\\left\(t\\right\)=\\mathcal\{Q\}\\left\(\\boldsymbol\{c\}\_\{p\}\\left\(t\\right\)\\right\)\\odot\\boldsymbol\{\\pi\}\_\{p\-1\}\\left\(t\\right\),\(5\)𝝅~p​\(t\)=𝑾π~​\(t\)​\[𝝅ˇp​\(t\);ℛ​\(𝒚p​\(t\)\)\]\+𝒃π~​\(t\),\\displaystyle\\boldsymbol\{\\tilde\{\\pi\}\}\_\{p\}\\left\(t\\right\)=\\boldsymbol\{W\}\_\{\\tilde\{\\pi\}\}\\left\(t\\right\)\[\\check\{\\boldsymbol\{\\pi\}\}\_\{p\}\\left\(t\\right\);\\mathcal\{R\}\\left\(\\boldsymbol\{y\}\_\{p\}\\left\(t\\right\)\\right\)\]\+\\boldsymbol\{b\}\_\{\\tilde\{\\pi\}\}\\left\(t\\right\),\(6\)𝝅p​\(t\)=\(1−𝒬​\(𝒛p​\(t\)\)\)⊙𝝅p−1​\(t\)\\displaystyle\\boldsymbol\{\\pi\}\_\{p\}\\left\(t\\right\)=\\left\(1\-\\mathcal\{Q\}\\left\(\\boldsymbol\{z\}\_\{p\}\\left\(t\\right\)\\right\)\\right\)\\odot\\boldsymbol\{\\pi\}\_\{p\-1\}\\left\(t\\right\)\+𝒬​\(𝒛p​\(t\)\)⊙\(2​𝒬​\(𝝅~p​\(t\)\)−𝟏\),\\displaystyle\+\\mathcal\{Q\}\\left\(\\boldsymbol\{z\}\_\{p\}\\left\(t\\right\)\\right\)\\odot\\left\(2\\mathcal\{Q\}\\left\(\\boldsymbol\{\\tilde\{\\pi\}\}\_\{p\}\\left\(t\\right\)\\right\)\-\\boldsymbol\{1\}\\right\),\(7\)where𝒬​\(⋅\)\\mathcal\{Q\}\(\\cdot\)is an activation function applied element\-wise, which is often chosen as the sigmoid function in traditional GRU networks\. In this paper, we set𝒬\\mathcal\{Q\}to the probit function, i\.e\.,𝒬​\(x\)=∫−∞x𝒩​\(t;0,1\)​dt\\mathcal\{Q\}\(x\)=\\int\_\{\-\\infty\}^\{x\}\\mathcal\{N\}\(t;0,1\)\\,\\mathrm\{d\}t, to facilitate message passing\. The RGRU parameters are summarized in a set𝝎​\(t\)=\{𝑾y~​\(t\),𝑾z​\(t\),𝑾c​\(t\),𝑾π~​\(t\),𝒃y~​\(t\),𝒃z​\(t\),𝒃c​\(t\),𝒃π~​\(t\)\}\\boldsymbol\{\\omega\}\(t\)=\\left\\\{\\boldsymbol\{W\}\_\{\\tilde\{y\}\}\(t\),\\boldsymbol\{W\}\_\{z\}\(t\),\\boldsymbol\{W\}\_\{c\}\(t\),\\boldsymbol\{W\}\_\{\\tilde\{\\pi\}\}\(t\),\\boldsymbol\{b\}\_\{\\tilde\{y\}\}\(t\),\\boldsymbol\{b\}\_\{z\}\(t\),\\boldsymbol\{b\}\_\{c\}\(t\),\\boldsymbol\{b\}\_\{\\tilde\{\\pi\}\}\(t\)\\right\\\}which are estimated via approximate message passing\.

### II\-CTwo\-Timescale Markov Priors

We adopt a random\-process framework over time slots to jointly capture the angular\-domain channel𝒙\\boldsymbol\{x\}and impairment\-related parameters𝝎\\boldsymbol\{\\omega\}\. Owing to small\-scale fading, the channel is treated as approximately stationary within each slot, yet it varies rapidly across successive slots\[[39](https://arxiv.org/html/2607.01660#bib.bib77)\]\. The parameters of the memory hardware impairment model, on the other hand, are assumed to evolve slowly throughout the channel tracking process\. Consequently, the channel dynamics are modeled by a fast\-varying Markov chain across slots, while the RGRU weights follow a slowly\-varying Gaussian Markov chain\.

Specifically, for the channel𝒙​\(t\)\\boldsymbol\{x\}\\left\(t\\right\), we define two hidden variables𝒔¯​\(t\)\\bar\{\\boldsymbol\{s\}\}\\left\(t\\right\)and𝒙¯​\(t\)\\bar\{\\boldsymbol\{x\}\}\\left\(t\\right\), where𝒔¯​\(t\)\\bar\{\\boldsymbol\{s\}\}\\left\(t\\right\)denotes the support vector to indicate activity, and𝒙¯​\(t\)\\bar\{\\boldsymbol\{x\}\}\\left\(t\\right\)represents𝒙​\(t\)\\boldsymbol\{x\}\\left\(t\\right\)when active\[[17](https://arxiv.org/html/2607.01660#bib.bib40)\]\. Therefore,𝒙​\(t\)\\boldsymbol\{x\}\\left\(t\\right\)can be written as \([8](https://arxiv.org/html/2607.01660#S2.E8)\):

𝒙​\(t\)=𝒔¯​\(t\)⊙𝒙¯​\(t\),\\boldsymbol\{x\}\\left\(t\\right\)=\\bar\{\\boldsymbol\{s\}\}\\left\(t\\right\)\\odot\\bar\{\\boldsymbol\{x\}\}\\left\(t\\right\),\(8\)where𝒔¯​\(t\)\\bar\{\\boldsymbol\{s\}\}\\left\(t\\right\)represents the binary support vector, which describes the dynamic sparsity of the channel support pattern, and the complex\-valued vector𝒙¯​\(t\)\\bar\{\\boldsymbol\{x\}\}\\left\(t\\right\)represents the hidden value\. We assume independent dynamic priors for𝒔¯​\(t\)\\bar\{\\boldsymbol\{s\}\}\\left\(t\\right\)and𝒙¯​\(t\)\\bar\{\\boldsymbol\{x\}\}\\left\(t\\right\), i\.e\.:

p​\(𝒔¯​\(t\)\|𝒔¯​\(t−1\)\)=∏n=1Np​\(s¯n​\(t\)\|s¯n​\(t−1\)\),p\\left\(\\bar\{\\boldsymbol\{s\}\}\\left\(t\\right\)\|\\bar\{\\boldsymbol\{s\}\}\\left\(t\-1\\right\)\\right\)=\\prod\_\{n=1\}^\{N\}p\\left\(\\bar\{s\}\_\{n\}\\left\(t\\right\)\|\\bar\{s\}\_\{n\}\\left\(t\-1\\right\)\\right\),p​\(𝒙¯​\(t\)\|𝒙¯​\(t−1\)\)=∏n=1Np​\(x¯n​\(t\)\|x¯n​\(t−1\)\),p\\left\(\\bar\{\\boldsymbol\{x\}\}\\left\(t\\right\)\|\\bar\{\\boldsymbol\{x\}\}\\left\(t\-1\\right\)\\right\)=\\prod\_\{n=1\}^\{N\}p\\left\(\\bar\{x\}\_\{n\}\\left\(t\\right\)\|\\bar\{x\}\_\{n\}\\left\(t\-1\\right\)\\right\),where

p​\(s¯n​\(t\)\|s¯n​\(t−1\)\)\\displaystyle p\\left\(\\bar\{s\}\_\{n\}\\left\(t\\right\)\|\\bar\{s\}\_\{n\}\\left\(t\-1\\right\)\\right\)\(9\)=\{ρ01​δ​\(s¯n​\(t\)−1\)\+ρ00​δ​\(s¯n​\(t\)\)s¯n​\(t−1\)=0ρ11​δ​\(s¯n​\(t\)−1\)\+ρ10​δ​\(s¯n​\(t\)\)s¯n​\(t−1\)=1,\\displaystyle=,p​\(x¯n​\(t\)\|x¯n​\(t−1\)\)\\displaystyle p\\left\(\\bar\{x\}\_\{n\}\\left\(t\\right\)\|\\bar\{x\}\_\{n\}\\left\(t\-1\\right\)\\right\)\(10\)=𝒞​𝒩​\(x¯n​\(t\);\(1−α\)​x¯n​\(t−1\)\+α​ξ,α2​κ\),\\displaystyle=\\mathcal\{CN\}\\left\(\\bar\{x\}\_\{n\}\\left\(t\\right\);\\left\(1\-\\alpha\\right\)\\bar\{x\}\_\{n\}\\left\(t\-1\\right\)\+\\alpha\\xi,\\alpha^\{2\}\\kappa\\right\),withρ00=1−ρ01,ρ10=1−ρ11\\rho\_\{00\}=1\-\\rho\_\{01\},\\rho\_\{10\}=1\-\\rho\_\{11\}, and we defined¯s,nt=p​\(s¯n​\(t\)\|s¯n​\(t−1\)\),d¯x,nt=p​\(x¯n​\(t\)\|x¯n​\(t−1\)\)\\bar\{d\}\_\{s,n\}^\{t\}=p\\left\(\\bar\{s\}\_\{n\}\\left\(t\\right\)\|\\bar\{s\}\_\{n\}\\left\(t\-1\\right\)\\right\),\\bar\{d\}\_\{x,n\}^\{t\}=p\\left\(\\bar\{x\}\_\{n\}\\left\(t\\right\)\|\\bar\{x\}\_\{n\}\\left\(t\-1\\right\)\\right\)in the following parts\. Note thatρ01/ρ11\\rho\_\{01\}/\\rho\_\{11\}represents the probability that the support turns active from inactive/active,α\\alpharepresents the temporal correlation between time slots,ξ\\xirepresents the mean of the angular\-domain channel process, andκ\\kapparepresents the variance of Gaussian perturbation\. Largerρ01\\rho\_\{01\},ρ10\\rho\_\{10\},α\\alpha, andκ\\kappamean that the channel varies faster over time\.

For the RGRU hidden states, we abbreviate the updating rules in \([3](https://arxiv.org/html/2607.01660#S2.E3)\)–\([7](https://arxiv.org/html/2607.01660#S2.E7)\) at symbolppto𝝅p​\(t\)=GRU​\(𝝅p−1​\(t\),ℛ​\(𝒚p​\(t\)\);𝝎​\(t\)\)\\boldsymbol\{\\pi\}\_\{p\}\(t\)=\\mathrm\{GRU\}\(\\boldsymbol\{\\pi\}\_\{p\-1\}\(t\),\\mathcal\{R\}\(\\boldsymbol\{y\}\_\{p\}\(t\)\);\\boldsymbol\{\\omega\}\(t\)\)in order to construct the top\-level factor graph\. To track slow impairment drift, we model each elementωi​\(t\)∈𝝎​\(t\)\\omega\_\{i\}\(t\)\\in\\boldsymbol\{\\omega\}\(t\)with a first\-order Gaussian Markov transition,

p​\(ωi​\(t\)\|ωi​\(t−1\)\)=𝒩​\(ωi​\(t\);ωi​\(t−1\),σωi2\),p\\left\(\\omega\_\{i\}\\left\(t\\right\)\|\\omega\_\{i\}\\left\(t\-1\\right\)\\right\)=\\mathcal\{N\}\\left\(\\omega\_\{i\}\\left\(t\\right\);\\omega\_\{i\}\\left\(t\-1\\right\),\\sigma\_\{\\omega\_\{i\}\}^\{2\}\\right\),\(11\)whereσωi2\\sigma\_\{\\omega\_\{i\}\}^\{2\}controls the per\-slot drift variance ofωi\\omega\_\{i\}, withσωi2\\sigma\_\{\\omega\_\{i\}\}^\{2\}chosen much smaller than the channel transition variance so that𝝎​\(t\)\\boldsymbol\{\\omega\}\(t\)evolves on a slower timescale than𝒔¯​\(t\)\\bar\{\\boldsymbol\{s\}\}\(t\)and𝒙¯​\(t\)\\bar\{\\boldsymbol\{x\}\}\(t\)\.

## IIIProposed Algorithm

In this section, we detail the MP\-TTBDL framework to jointly track the channel and the impairments from online pilots, based on the models in Section II\. MP\-TTBDL consists of two components: inter\-slot forward message passing, which leverages the two\-timescale Markov priors to propagate temporal messages across slots, and per\-slot iterative message passing, which performs joint Bayesian inference of the channel and impairment parameters\. The two components are elaborated in the following subsections\.

### III\-ATop\-Level Factor Graph and Message Passing Overview

To establish a Bayesian framework for joint channel estimation and impairment tracking, we formulate the probabilistic model overTTslots\. The objective is to estimate the marginal posteriors of channel variables\{𝒔¯​\(t\),𝒙¯​\(t\),𝒙​\(t\)\}\\\{\\bar\{\\boldsymbol\{s\}\}\(t\),\\bar\{\\boldsymbol\{x\}\}\(t\),\\boldsymbol\{x\}\(t\)\\\}and impairment parameters\{𝝎​\(t\)\}\\\{\\boldsymbol\{\\omega\}\(t\)\\\}from distorted observations\{𝒚~p​\(t\)\}p=1P\\\{\\tilde\{\\boldsymbol\{y\}\}\_\{p\}\(t\)\\\}\_\{p=1\}^\{P\}fort=1,…,Tt=1,\\ldots,T\. Combining the RGRU\-based observation model in \([2](https://arxiv.org/html/2607.01660#S2.E2)\)–\([7](https://arxiv.org/html/2607.01660#S2.E7)\), the two\-timescale Markov prior in \([9](https://arxiv.org/html/2607.01660#S2.E9)\)–\([11](https://arxiv.org/html/2607.01660#S2.E11)\), and the channel observation model in \([1](https://arxiv.org/html/2607.01660#S2.E1)\), the joint distribution of all random variables is written as:

p​\(\{𝒔¯​\(t\),𝒙¯​\(t\),𝒙​\(t\),\{𝒚p​\(t\),𝝅p​\(t\),𝒚~p​\(t\)\}p=1P,𝝎​\(t\)\}t=1,…,T\)\\displaystyle p\\Bigl\(\\bigl\\\{\\bar\{\\boldsymbol\{s\}\}\(t\),\\bar\{\\boldsymbol\{x\}\}\(t\),\\boldsymbol\{x\}\(t\),\\\{\\boldsymbol\{y\}\_\{p\}\(t\),\\boldsymbol\{\\pi\}\_\{p\}\(t\),\\tilde\{\\boldsymbol\{y\}\}\_\{p\}\(t\)\\\}\_\{p=1\}^\{P\},\\boldsymbol\{\\omega\}\(t\)\\bigr\\\}\_\{t=1,\\ldots,T\}\\Bigr\)∝p​\(𝒔¯​\(1\)\)​p​\(𝒙¯​\(1\)\)​p​\(𝝎​\(1\)\)\\displaystyle\\propto p\\left\(\\bar\{\\boldsymbol\{s\}\}\(1\)\\right\)p\\left\(\\bar\{\\boldsymbol\{x\}\}\(1\)\\right\)p\\left\(\\boldsymbol\{\\omega\}\(1\)\\right\)×∏t=2Tp\(𝒔¯\(t\)\|𝒔¯\(t−1\)\)p\(𝒙¯\(t\)\|𝒙¯\(t−1\)\)p\(𝝎\(t\)\|𝝎\(t−1\)\)\\displaystyle\\times\\prod\_\{t=2\}^\{T\}p\\left\(\\bar\{\\boldsymbol\{s\}\}\(t\)\|\\bar\{\\boldsymbol\{s\}\}\(t\-1\)\\right\)p\\left\(\\bar\{\\boldsymbol\{x\}\}\(t\)\|\\bar\{\\boldsymbol\{x\}\}\(t\-1\)\\right\)p\\left\(\\boldsymbol\{\\omega\}\(t\)\|\\boldsymbol\{\\omega\}\(t\-1\)\\right\)×∏t=1Tδ\(𝒙\(t\)−𝒔¯\(t\)⊙𝒙¯\(t\)\)∏p=1Pp\(𝒚p\(t\)\|𝒙\(t\)\)\\displaystyle\\times\\prod\_\{t=1\}^\{T\}\\delta\\left\(\\boldsymbol\{x\}\(t\)\-\\bar\{\\boldsymbol\{s\}\}\(t\)\\odot\\bar\{\\boldsymbol\{x\}\}\(t\)\\right\)\\prod\_\{p=1\}^\{P\}p\\left\(\\boldsymbol\{y\}\_\{p\}\(t\)\|\\boldsymbol\{x\}\(t\)\\right\)×∏t=1T∏p=1Pδ\(𝝅p\(t\)−GRU\(𝝅p−1\(t\),ℛ\(𝒚p\(t\)\);𝝎\(t\)\)\)\\displaystyle\\times\\prod\_\{t=1\}^\{T\}\\prod\_\{p=1\}^\{P\}\\delta\\bigl\(\\boldsymbol\{\\pi\}\_\{p\}\(t\)\-\\mathrm\{GRU\}\(\\boldsymbol\{\\pi\}\_\{p\-1\}\(t\),\\mathcal\{R\}\(\\boldsymbol\{y\}\_\{p\}\(t\)\);\\boldsymbol\{\\omega\}\(t\)\)\\bigr\)×∏t=1Tp\(𝝅0\(t\)\)∏p=1Pp\(𝒚~p\(t\)\|𝝅p\(t\),𝒚p\(t\),𝝎\(t\)\)\.\\displaystyle\\times\\prod\_\{t=1\}^\{T\}p\\left\(\\boldsymbol\{\\pi\}\_\{0\}\(t\)\\right\)\\prod\_\{p=1\}^\{P\}p\\left\(\\tilde\{\\boldsymbol\{y\}\}\_\{p\}\(t\)\|\\boldsymbol\{\\pi\}\_\{p\}\(t\),\\boldsymbol\{y\}\_\{p\}\(t\),\\boldsymbol\{\\omega\}\(t\)\\right\)\.\(12\)wherep​\(𝒚p​\(t\)\|𝒙​\(t\)\)p\(\\boldsymbol\{y\}\_\{p\}\(t\)\|\\boldsymbol\{x\}\(t\)\)andp​\(𝒚~p​\(t\)\|𝝅p​\(t\),𝒚p​\(t\),𝝎​\(t\)\)p\(\\tilde\{\\boldsymbol\{y\}\}\_\{p\}\(t\)\|\\boldsymbol\{\\pi\}\_\{p\}\(t\),\\boldsymbol\{y\}\_\{p\}\(t\),\\boldsymbol\{\\omega\}\(t\)\)come from \([1](https://arxiv.org/html/2607.01660#S2.E1)\) and \([2](https://arxiv.org/html/2607.01660#S2.E2)\), respectively:

p​\(𝒚p​\(t\)\|𝒙​\(t\)\)=𝒞​𝒩​\(𝒚p​\(t\);𝑨​\(ϑ\)​𝒙​\(t\),σ2​IN\)\.p\\left\(\\boldsymbol\{y\}\_\{p\}\(t\)\|\\boldsymbol\{x\}\(t\)\\right\)=\\mathcal\{CN\}\\bigl\(\\boldsymbol\{y\}\_\{p\}\(t\);\\boldsymbol\{A\}\(\\boldsymbol\{\\vartheta\}\)\\boldsymbol\{x\}\(t\),\\sigma^\{2\}I\_\{N\}\\bigr\)\.\(13\)p​\(𝒚~p​\(t\)\|𝝅p​\(t\),𝒚p​\(t\),𝝎​\(t\)\)\\displaystyle p\\left\(\\tilde\{\\boldsymbol\{y\}\}\_\{p\}\(t\)\|\\boldsymbol\{\\pi\}\_\{p\}\(t\),\\boldsymbol\{y\}\_\{p\}\(t\),\\boldsymbol\{\\omega\}\(t\)\\right\)=𝒩​\(ℛ​\(𝒚~p​\(t\)\);ℛ​\(𝒚p​\(t\)\)\+𝑾y~​\(t\)​𝝅p​\(t\)\+𝒃y~​\(t\),σ~2​𝑰2​N\)\.\\displaystyle=\\mathcal\{N\}\\bigl\(\\mathcal\{R\}\\left\(\\tilde\{\\boldsymbol\{y\}\}\_\{p\}\(t\)\\right\);\\mathcal\{R\}\\left\(\\boldsymbol\{y\}\_\{p\}\(t\)\\right\)\+\\boldsymbol\{W\}\_\{\\tilde\{y\}\}\(t\)\\boldsymbol\{\\pi\}\_\{p\}\(t\)\+\\boldsymbol\{b\}\_\{\\tilde\{y\}\}\(t\),\\tilde\{\\sigma\}^\{2\}\\boldsymbol\{I\}\_\{2N\}\\bigr\)\.\(14\)In addition,p​\(𝝅0​\(t\)\)=𝒩​\(𝝅0​\(t\);𝟎Nπ,σπ02​𝑰Nπ\)p\(\\boldsymbol\{\\pi\}\_\{0\}\(t\)\)=\\mathcal\{N\}\(\\boldsymbol\{\\pi\}\_\{0\}\(t\);\\boldsymbol\{0\}\_\{N\_\{\\pi\}\},\\sigma\_\{\\pi\_\{0\}\}^\{2\}\\boldsymbol\{I\}\_\{N\_\{\\pi\}\}\)assigns a zero\-mean prior with small variance to the initial RGRU hidden state at the first pilot of each slot\.

To compute the marginal posterior distributions, we construct the multi\-slot factor graph shown in Fig\.[2](https://arxiv.org/html/2607.01660#S3.F2), with the factor nodes summarized in Table[I](https://arxiv.org/html/2607.01660#S3.T1), and develop a corresponding approximate message passing algorithm\. For notational convenience, we define the stacked observation vector𝐲=\[𝒚1​\(t\)T,…,𝒚P​\(t\)T\]T∈ℂN​P\\mathbf\{y\}=\[\\boldsymbol\{y\}\_\{1\}\(t\)^\{\\mathrm\{T\}\},\\ldots,\\boldsymbol\{y\}\_\{P\}\(t\)^\{\\mathrm\{T\}\}\]^\{\\mathrm\{T\}\}\\in\\mathbb\{C\}^\{NP\}and the stacked noise vector𝐧\\mathbf\{n\}similarly, under which the per\-slot observation model can be rewritten in the compact form𝐲=\(𝒓⊗𝑨​\(ϑ\)\)​𝒙\+𝐧\\mathbf\{y\}=\(\\boldsymbol\{r\}\\otimes\\boldsymbol\{A\}\(\\boldsymbol\{\\vartheta\}\)\)\\,\\boldsymbol\{x\}\+\\mathbf\{n\}, where𝒓\\boldsymbol\{r\}collects thePPpilot symbols\. At each slottt, inter\-slot messages from slott−1t\-1are first combined with the two\-timescale Markov transition factors to form equivalent priors on𝒔¯​\(t\)\\bar\{\\boldsymbol\{s\}\}\(t\),𝒙¯​\(t\)\\bar\{\\boldsymbol\{x\}\}\(t\), and𝝎​\(t\)\\boldsymbol\{\\omega\}\(t\)\. Within the slot, the factor graph is partitioned into a channel tracking \(CT\) module and an impairments calibration \(IC\) module, which exchange extrinsic information overτmax\\tau\_\{\\max\}turbo iterations: the CT module uses Turbo\-OAMP to estimate the sparse channel𝒙​\(t\)\\boldsymbol\{x\}\(t\), while the IC module uses DAMP to update the RGRU parameters𝝎​\(t\)\\boldsymbol\{\\omega\}\(t\)\. After convergence, the updated posteriors are forwarded to slott\+1t\+1via inter\-slot message passing, and the process repeats\. Meanwhile, based on the channel estimation results from the CT module, an EM\-based surrogate function for the angular\-domain grid can be constructed, and the grid points are dynamically updated via gradient descent as detailed in\[[17](https://arxiv.org/html/2607.01660#bib.bib40)\]\. We shall useΔa→b\\Delta\_\{a\\rightarrow b\}to denote the message from nodeaato nodebb\.

![Refer to caption](https://arxiv.org/html/2607.01660v1/x2.png)Figure 2:Top\-level factor graph for \([12](https://arxiv.org/html/2607.01660#S3.E12)\) with factor definitions in Table[I](https://arxiv.org/html/2607.01660#S3.T1)\.TABLE I:Illustration of factor nodes in Fig\.[2](https://arxiv.org/html/2607.01660#S3.F2)\.FactorDistributionFunctional formdstd\_\{s\}^\{t\}p​\(𝒔¯​\(t\)\|𝒔¯​\(t−1\)\)p\(\\bar\{\\boldsymbol\{s\}\}\(t\)\|\\bar\{\\boldsymbol\{s\}\}\(t\-1\)\)\([9](https://arxiv.org/html/2607.01660#S2.E9)\)dxtd\_\{x\}^\{t\}p​\(𝒙¯​\(t\)\|𝒙¯​\(t−1\)\)p\(\\bar\{\\boldsymbol\{x\}\}\(t\)\|\\bar\{\\boldsymbol\{x\}\}\(t\-1\)\)\([10](https://arxiv.org/html/2607.01660#S2.E10)\)dtd^\{t\}p​\(𝒙​\(t\)\|𝒔¯​\(t\),𝒙¯​\(t\)\)p\(\\boldsymbol\{x\}\(t\)\|\\bar\{\\boldsymbol\{s\}\}\(t\),\\bar\{\\boldsymbol\{x\}\}\(t\)\)δ​\(𝒙​\(t\)−𝒔¯​\(t\)⊙𝒙¯​\(t\)\)\\delta\(\\boldsymbol\{x\}\(t\)\-\\bar\{\\boldsymbol\{s\}\}\(t\)\\odot\\bar\{\\boldsymbol\{x\}\}\(t\)\)ftf^\{t\}p​\(𝐲​\(t\)\|𝒙​\(t\)\)p\(\\mathbf\{y\}\(t\)\|\\boldsymbol\{x\}\(t\)\)\([13](https://arxiv.org/html/2607.01660#S3.E13)\)ΔC​T→I​C​\(t\)\\Delta\_\{CT\\rightarrow IC\}\(t\)Extrinsic message from CT to IC on𝐲​\(t\)\\mathbf\{y\}\(t\)\([21](https://arxiv.org/html/2607.01660#S3.E21)\)dωtd\_\{\\omega\}^\{t\}p​\(𝝎​\(t\)\|𝝎​\(t−1\)\)p\(\\boldsymbol\{\\omega\}\(t\)\|\\boldsymbol\{\\omega\}\(t\-1\)\)\([11](https://arxiv.org/html/2607.01660#S2.E11)\)ΔI​C→C​T​\(t\)\\Delta\_\{IC\\rightarrow CT\}\(t\)Extrinsic message from IC to CT on𝐲​\(t\)\\mathbf\{y\}\(t\)\([25](https://arxiv.org/html/2607.01660#S3.E25)\)G​R​UtGRU^\{t\}p​\(𝐲^​\(t\)\|𝐲​\(t\),𝝎​\(t\)\)p\(\\hat\{\\mathbf\{y\}\}\(t\)\|\\mathbf\{y\}\(t\),\\boldsymbol\{\\omega\}\(t\)\)δ​\(𝐲^​\(t\)−GRU​\(𝐲​\(t\);𝝎​\(t\)\)\)\\delta\(\\hat\{\\mathbf\{y\}\}\(t\)\-\\mathrm\{GRU\}\(\\mathbf\{y\}\(t\);\\boldsymbol\{\\omega\}\(t\)\)\)

gtg^\{t\}p​\(𝐲~​\(t\)\|𝐲^​\(t\),𝐲​\(t\),𝝎​\(t\)\)p\(\\tilde\{\\mathbf\{y\}\}\(t\)\|\\hat\{\\mathbf\{y\}\}\(t\),\\mathbf\{y\}\(t\),\\boldsymbol\{\\omega\}\(t\)\)

\([14](https://arxiv.org/html/2607.01660#S3.E14)\)
### III\-BForward Message Passing Across Slots

Here we only consider forward message passing across time slots\. Based on sum\-product rule, the messages are computed as follows:

Δ𝒔¯​\(t\)→dst\+1=Δdst→𝒔¯​\(t\)⋅Δdt→𝒔¯​\(t\),\\Delta\_\{\\bar\{\\boldsymbol\{s\}\}\\left\(t\\right\)\\rightarrow d\_\{s\}^\{t\+1\}\}=\\Delta\_\{d\_\{s\}^\{t\}\\rightarrow\\bar\{\\boldsymbol\{s\}\}\\left\(t\\right\)\}\\cdot\\Delta\_\{d^\{t\}\\rightarrow\\bar\{\\boldsymbol\{s\}\}\\left\(t\\right\)\},\(15\)Δ𝒙¯​\(t\)→dxt\+1=Δdxt→𝒙¯​\(t\)⋅Δdt→𝒙¯​\(t\),\\Delta\_\{\\bar\{\\boldsymbol\{x\}\}\\left\(t\\right\)\\rightarrow d\_\{x\}^\{t\+1\}\}=\\Delta\_\{d\_\{x\}^\{t\}\\rightarrow\\bar\{\\boldsymbol\{x\}\}\\left\(t\\right\)\}\\cdot\\Delta\_\{d^\{t\}\\rightarrow\\bar\{\\boldsymbol\{x\}\}\\left\(t\\right\)\},\(16\)Δ𝝎​\(t\)→dωt\+1=Δdωt→𝝎​\(t\)⋅ΔG​R​Ut→𝝎​\(t\),\\Delta\_\{\\boldsymbol\{\\omega\}\(t\)\\rightarrow d\_\{\\omega\}^\{t\+1\}\}=\\Delta\_\{d\_\{\\omega\}^\{t\}\\rightarrow\\boldsymbol\{\\omega\}\(t\)\}\\cdot\\Delta\_\{GRU^\{t\}\\rightarrow\\boldsymbol\{\\omega\}\(t\)\},\(17\)
For messages within thett\-th time slot, we first denote equivalent input prior factorsdˇst,dˇxt\\check\{d\}\_\{s\}^\{t\},\\check\{d\}\_\{x\}^\{t\}anddˇωt\\check\{d\}\_\{\\omega\}^\{t\}that respectively represent the prior information of𝒔¯​\(t\)\\bar\{\\boldsymbol\{s\}\}\(t\),𝒙¯​\(t\)\\bar\{\\boldsymbol\{x\}\}\(t\)and𝝎​\(t\)\\boldsymbol\{\\omega\}\(t\)extracted from all previous observations:

dˇst=Δdst→𝒔¯​\(t\)=∫𝒔¯​\(t−1\)p​\(𝒔¯​\(t\)\|𝒔¯​\(t−1\)\)​Δ𝒔¯​\(t−1\)→dst,\\check\{d\}\_\{s\}^\{t\}=\\Delta\_\{d\_\{s\}^\{t\}\\rightarrow\\bar\{\\boldsymbol\{s\}\}\\left\(t\\right\)\}=\\int\_\{\\bar\{\\boldsymbol\{s\}\}\\left\(t\-1\\right\)\}p\\left\(\\bar\{\\boldsymbol\{s\}\}\\left\(t\\right\)\|\\bar\{\\boldsymbol\{s\}\}\\left\(t\-1\\right\)\\right\)\\Delta\_\{\\bar\{\\boldsymbol\{s\}\}\\left\(t\-1\\right\)\\rightarrow d\_\{s\}^\{t\}\},\(18\)dˇxt=Δdxt→𝒙¯​\(t\)=∫𝒙¯​\(t−1\)p​\(𝒙¯​\(t\)\|𝒙¯​\(t−1\)\)​Δ𝒙¯​\(t−1\)→dxt,\\check\{d\}\_\{x\}^\{t\}=\\Delta\_\{d\_\{x\}^\{t\}\\rightarrow\\bar\{\\boldsymbol\{x\}\}\\left\(t\\right\)\}=\\int\_\{\\bar\{\\boldsymbol\{x\}\}\\left\(t\-1\\right\)\}p\\left\(\\bar\{\\boldsymbol\{x\}\}\\left\(t\\right\)\|\\bar\{\\boldsymbol\{x\}\}\\left\(t\-1\\right\)\\right\)\\Delta\_\{\\bar\{\\boldsymbol\{x\}\}\\left\(t\-1\\right\)\\rightarrow d\_\{x\}^\{t\}\},\(19\)dˇωt=Δdωt→𝝎​\(t\)=∫𝝎​\(t−1\)p​\(𝝎​\(t\)\|𝝎​\(t−1\)\)​Δ𝝎​\(t−1\)→dωt\.\\check\{d\}\_\{\\omega\}^\{t\}=\\Delta\_\{d\_\{\\omega\}^\{t\}\\rightarrow\\boldsymbol\{\\omega\}\(t\)\}=\\int\_\{\\boldsymbol\{\\omega\}\(t\-1\)\}p\\left\(\\boldsymbol\{\\omega\}\(t\)\|\\boldsymbol\{\\omega\}\(t\-1\)\\right\)\\Delta\_\{\\boldsymbol\{\\omega\}\(t\-1\)\\rightarrow d\_\{\\omega\}^\{t\}\}\.\(20\)Based on the Markov priors in \([9](https://arxiv.org/html/2607.01660#S2.E9)\)–\([11](https://arxiv.org/html/2607.01660#S2.E11)\), \([18](https://arxiv.org/html/2607.01660#S3.E18)\)–\([20](https://arxiv.org/html/2607.01660#S3.E20)\) can be evaluated in closed form\.

### III\-CMessage Passing Within Each Slot

At each time slottt, the receiver performs joint inference by iteratively exchanging messages between the CT and IC modules forτmax\\tau\_\{\\max\}turbo iterations\. In each turbo iteration, the CT module updates messages via Turbo\-OAMP, while the IC module updates messages via DAMP\. Then, the two modules exchange extrinsic messages\. In particular,ΔI​C→C​T\\Delta\_\{IC\\rightarrow CT\}aggregates the IC backward Gaussian messages on\{𝒚p​\(t\)\}p=1P\\\{\\boldsymbol\{y\}\_\{p\}\(t\)\\\}\_\{p=1\}^\{P\}and is treated as equivalent AWGN side information in the CT module, whereasΔC​T→I​C\\Delta\_\{CT\\rightarrow IC\}is formed via EP and is treated as equivalent AWGN side information in the IC module\. After convergence, the messages passed to time slott\+1t\+1are updated based on the sum\-product rule as \([15](https://arxiv.org/html/2607.01660#S3.E15)\)–\([17](https://arxiv.org/html/2607.01660#S3.E17)\)\. The message passing procedures with respect to the CT and IC modules are elaborated in the following; the slot indexttis dropped for brevity, as no ambiguity arises\.

#### III\-C1CT Module \(Turbo\-OAMP\)

The CT module estimates the sparse angular\-domain channel𝒙​\(t\)\\boldsymbol\{x\}\(t\)from the stacked observations𝐲\\mathbf\{y\}within slottt\. We apply Turbo\-OAMP\[[17](https://arxiv.org/html/2607.01660#bib.bib40)\]to perform this estimation efficiently\. Specifically, the CT module is further divided into sub\-module A and sub\-module B, which extract information from the observations and the sparse prior factors, respectively\. Firstly, sub\-module B treats the extrinsic message w\.r\.t\.𝒙\\boldsymbol\{x\}from sub\-module A as equivalent AWGN observations, computes the posterior distribution of𝒙\\boldsymbol\{x\}based on standard sum\-product rule, and then updates the extrinsic message to sub\-module A via EP\. Next, sub\-module A treats the extrinsic message w\.r\.t\.𝒙\\boldsymbol\{x\}from sub\-module B as Gaussian prior, computes the posterior for𝒙\\boldsymbol\{x\}via linear minimum mean square error \(LMMSE\), and then updates the extrinsic message to sub\-module B via EP\. The two sub\-modules exchange extrinsic messages to fuse the observations and sparse prior information until convergence\. The detailed updating rules for messages within the CT module are the same as Turbo\-OAMP\. Please refer to\[[18](https://arxiv.org/html/2607.01660#bib.bib39),[17](https://arxiv.org/html/2607.01660#bib.bib40)\]for the details\.

Here we only give the updating rule for the top\-level extrinsic messageΔC​T→I​C\\Delta\_\{CT\\rightarrow IC\}from the CT module to the IC module\. Plugging in the messageΔ𝒙→f=𝒞​𝒩​\(𝒙;𝝁𝒙→f,Σ𝒙→f\)\\Delta\_\{\\boldsymbol\{x\}\\rightarrow f\}=\\mathcal\{CN\}\(\\boldsymbol\{x\};\\boldsymbol\{\\mu\}\_\{\\boldsymbol\{x\}\\rightarrow f\},\\varSigma\_\{\\boldsymbol\{x\}\\rightarrow f\}\)computed via Turbo\-OAMP, the mean and variance of messageΔf→𝒚p\\Delta\_\{f\\rightarrow\\boldsymbol\{y\}\_\{p\}\}for each symbolppare given by:

𝝁f→𝒚p=𝑨​\(ϑ\)​𝝁𝒙→f,Σf→𝒚p=𝑨​\(ϑ\)​Σ𝒙→f​𝑨H​\(ϑ\)\+σ2​IN\.\\boldsymbol\{\\mu\}\_\{f\\rightarrow\\boldsymbol\{y\}\_\{p\}\}=\\boldsymbol\{A\}\(\\boldsymbol\{\\vartheta\}\)\\boldsymbol\{\\mu\}\_\{\\boldsymbol\{x\}\\rightarrow f\},\\quad\\varSigma\_\{f\\rightarrow\\boldsymbol\{y\}\_\{p\}\}=\\boldsymbol\{A\}\(\\boldsymbol\{\\vartheta\}\)\\varSigma\_\{\\boldsymbol\{x\}\\rightarrow f\}\\boldsymbol\{A\}^\{H\}\(\\boldsymbol\{\\vartheta\}\)\+\\sigma^\{2\}I\_\{N\}\.The covariance matrix is not diagonal, and thus the extrinsic message is further approximated via EP:

ΔC​T→I​C≈𝒫​\{ΔI​C→C​T​∏p=1PΔf→𝒚p\}ΔI​C→C​T\.\\Delta\_\{CT\\rightarrow IC\}\\approx\\frac\{\\mathcal\{P\}\\left\\\{\\Delta\_\{IC\\rightarrow CT\}\\prod\_\{p=1\}^\{P\}\\Delta\_\{f\\rightarrow\\boldsymbol\{y\}\_\{p\}\}\\right\\\}\}\{\\Delta\_\{IC\\rightarrow CT\}\}\.\(21\)where𝒫​\{⋅\}\\mathcal\{P\}\\left\\\{\\cdot\\right\\\}approximates the distribution as a product of independent Gaussian via moment matching\.

#### III\-C2IC Module \(DAMP for RGRU\)

Under the turbo framework, the IC module within each slot is separated for processing the messages associated with the RGRU\-based impairment model, where the memory hardware impairments evolve over thePPpilot symbols\. Based on the observation model given in \([2](https://arxiv.org/html/2607.01660#S2.E2)\)–\([7](https://arxiv.org/html/2607.01660#S2.E7)\), we take the processing of thepp\-th symbol as an illustrative example and construct the detailed factor graph of the IC module, as depicted in Fig\.[3](https://arxiv.org/html/2607.01660#S3.F3), with the definitions of the factor nodes provided in Table[II](https://arxiv.org/html/2607.01660#S3.T2)\. Owing to the recurrent state transitions and gating nonlinearities inherent in the RGRU architecture, the resulting factor graph contains numerous loops, making exact marginal inference computationally intractable\. While the DAMP framework in\[[33](https://arxiv.org/html/2607.01660#bib.bib71)\]provides a powerful tool for Bayesian learning in feedforward deep neural networks \(DNNs\), it cannot be directly applied to our recurrent architecture\. Extending DAMP to the RGRU presents two fundamental design challenges: \(i\) The recurrent state transition creates a non\-trivial, loopy factor graph that requires a carefully coordinated message schedule, unlike the straightforward layer\-wise schedule in feedforward networks\. \(ii\) Both the probit activation functions𝒬​\(⋅\)\\mathcal\{Q\}\(\\cdot\)and the element\-wise products between gates and states introduce non\-Gaussian messages, requiring tailored approximations to keep the message passing tractable\. To address these challenges, we develop a custom DAMP procedure tailored to the RGRU\. Our solution establishes a principled intra\-symbol forward message pass \(FMP\) and backward message pass \(BMP\) schedule, as summarized in Table[III](https://arxiv.org/html/2607.01660#S3.T3), that systematically propagates messages along thePPpilot symbols while respecting the complex loopy structure of the recurrent graph\. Furthermore, for the nonlinear submodules, we derive closed\-form variational Bayesian inference \(VBI\) updates in Appendix[A](https://arxiv.org/html/2607.01660#A1)to approximate the otherwise intractable messages\. This extension of the DAMP framework to a recurrent, gated architecture is a key methodological contribution that enables efficient online Bayesian learning of the IC module\. After completing message passing, the extrinsic messageΔI​C→C​T\\Delta\_\{IC\\to CT\}is approximated as the product of the backward Gaussian messages over all symbols, which is fed to the CT module as equivalent observations\. The entire IC module naturally decomposes into four bilinear submodules and two nonlinear submodules, each with its own equivalent input priors and observation model\. In the following, we name each submodule by its output and detail the message passing rules\.

![Refer to caption](https://arxiv.org/html/2607.01660v1/x3.png)Figure 3:Factor graph of the RGRU\-based IC module at symbolpp\.TABLE II:Illustration of factor nodes in Fig\.[3](https://arxiv.org/html/2607.01660#S3.F3)\.FactorDistributionFunctional formΔC​T→I​C\\Delta\_\{CT\\rightarrow IC\}Extrinsic message from CT to IC on𝒚p\\boldsymbol\{y\}\_\{p\}\([21](https://arxiv.org/html/2607.01660#S3.E21)\)dˇπp−1\\check\{d\}\_\{\\pi\_\{p\-1\}\}Equivalent prior for𝝅p−1\\boldsymbol\{\\pi\}\_\{p\-1\}\([24](https://arxiv.org/html/2607.01660#S3.E24)\)dˇωz\\check\{d\}\_\{\\omega\_\{z\}\}Equivalent prior for𝝎z\\boldsymbol\{\\omega\}\_\{z\}\([20](https://arxiv.org/html/2607.01660#S3.E20)\)dˇωc\\check\{d\}\_\{\\omega\_\{c\}\}Equivalent prior for𝝎c\\boldsymbol\{\\omega\}\_\{c\}dˇωπ\\check\{d\}\_\{\\omega\_\{\\pi\}\}Equivalent prior for𝝎π\\boldsymbol\{\\omega\}\_\{\\pi\}dˇωy~\\check\{d\}\_\{\\omega\_\{\\tilde\{y\}\}\}Equivalent prior for𝝎y~\\boldsymbol\{\\omega\}\_\{\\tilde\{y\}\}gzpg\_\{z\_\{p\}\}p​\(𝒛p\|𝝅p−1,𝒚p,𝝎z\)p\(\\boldsymbol\{z\}\_\{p\}\|\\boldsymbol\{\\pi\}\_\{p\-1\},\\boldsymbol\{y\}\_\{p\},\\boldsymbol\{\\omega\}\_\{z\}\)δ\\deltafunction from observation models \([3](https://arxiv.org/html/2607.01660#S2.E3)\)–\([7](https://arxiv.org/html/2607.01660#S2.E7)\) and \([14](https://arxiv.org/html/2607.01660#S3.E14)\), respectively\.gcpg\_\{c\_\{p\}\}p​\(𝒄p\|𝝅p−1,𝒚p,𝝎c\)p\(\\boldsymbol\{c\}\_\{p\}\|\\boldsymbol\{\\pi\}\_\{p\-1\},\\boldsymbol\{y\}\_\{p\},\\boldsymbol\{\\omega\}\_\{c\}\)gπˇpg\_\{\\check\{\\pi\}\_\{p\}\}p​\(𝝅ˇp\|𝒄p,𝝅p−1\)p\(\\check\{\\boldsymbol\{\\pi\}\}\_\{p\}\|\\boldsymbol\{c\}\_\{p\},\\boldsymbol\{\\pi\}\_\{p\-1\}\)gπ~pg\_\{\\tilde\{\\pi\}\_\{p\}\}p​\(𝝅~p\|𝝅ˇp,𝒚p,𝝎π\)p\(\\tilde\{\\boldsymbol\{\\pi\}\}\_\{p\}\|\\check\{\\boldsymbol\{\\pi\}\}\_\{p\},\\boldsymbol\{y\}\_\{p\},\\boldsymbol\{\\omega\}\_\{\\pi\}\)gπpg\_\{\\pi\_\{p\}\}p​\(𝝅p\|𝒛p,𝝅p−1,𝝅~p\)p\(\\boldsymbol\{\\pi\}\_\{p\}\|\\boldsymbol\{z\}\_\{p\},\\boldsymbol\{\\pi\}\_\{p\-1\},\\tilde\{\\boldsymbol\{\\pi\}\}\_\{p\}\)gy~pg\_\{\\tilde\{y\}\_\{p\}\}p​\(ℛ​\(𝒚~p\)\|𝝅p,𝒚p,𝝎y~\)p\(\\mathcal\{R\}\(\\tilde\{\\boldsymbol\{y\}\}\_\{p\}\)\|\\boldsymbol\{\\pi\}\_\{p\},\\boldsymbol\{y\}\_\{p\},\\boldsymbol\{\\omega\}\_\{\\tilde\{y\}\}\)Four bilinear submodules𝒛\\boldsymbol\{z\},𝒄\\boldsymbol\{c\},𝝅~\\tilde\{\\boldsymbol\{\\pi\}\}, and𝒚~\\tilde\{\\boldsymbol\{y\}\}, share the same equivalent observation model:

𝒚^eq=𝑾Y​𝒖eq\+𝒃Y\+𝒏eq,\\hat\{\\boldsymbol\{y\}\}\_\{\\mathrm\{eq\}\}=\\boldsymbol\{W\}\_\{Y\}\\boldsymbol\{u\}\_\{\\mathrm\{eq\}\}\+\\boldsymbol\{b\}\_\{Y\}\+\\boldsymbol\{n\}\_\{\\mathrm\{eq\}\},\(22\)where𝑾Y\\boldsymbol\{W\}\_\{Y\},𝒖eq\\boldsymbol\{u\}\_\{\\mathrm\{eq\}\}, and𝒃Y\\boldsymbol\{b\}\_\{Y\}are jointly estimated, the entries of𝒏eq\\boldsymbol\{n\}\_\{\\mathrm\{eq\}\}are i\.i\.d\. AWGN\. At symbolppof slottt, their corresponding equivalent inputs are\[𝝅p−1;ℛ​\(𝒚p\)\]\[\\boldsymbol\{\\pi\}\_\{p\-1\};\\mathcal\{R\}\(\\boldsymbol\{y\}\_\{p\}\)\],\[𝝅p−1;ℛ​\(𝒚p\)\]\[\\boldsymbol\{\\pi\}\_\{p\-1\};\\mathcal\{R\}\(\\boldsymbol\{y\}\_\{p\}\)\],\[𝝅ˇp;ℛ​\(𝒚p\)\]\[\\check\{\\boldsymbol\{\\pi\}\}\_\{p\};\\mathcal\{R\}\(\\boldsymbol\{y\}\_\{p\}\)\], and𝝅p\\boldsymbol\{\\pi\}\_\{p\}, respectively\. Since each bilinear submodule is identical to the linear layer model in the DAMP framework of\[[33](https://arxiv.org/html/2607.01660#bib.bib71)\]\. The forward and backward message updates are similar to Bayesian generalized approximate message passing \(BiGAMP\) in\[[23](https://arxiv.org/html/2607.01660#bib.bib68)\]\(refer to\[[33](https://arxiv.org/html/2607.01660#bib.bib71)\]for details\)\.

Two submodules,𝝅ˇ\\check\{\\boldsymbol\{\\pi\}\}and𝝅\\boldsymbol\{\\pi\}, deal with observation models

𝒚^eq=f​\(𝒖eq\)\+𝒏eq,\\hat\{\\boldsymbol\{y\}\}\_\{\\mathrm\{eq\}\}=f\(\\boldsymbol\{u\}\_\{\\mathrm\{eq\}\}\)\+\\boldsymbol\{n\}\_\{\\mathrm\{eq\}\},\(23\)wheref​\(⋅\)f\(\\cdot\)denotes a nonlinear transformation and𝒖eq\\boldsymbol\{u\}\_\{\\mathrm\{eq\}\}collects the submodule input variables\. Specifically, in submodule𝝅ˇ\\check\{\\boldsymbol\{\\pi\}\},ffis defined by \([5](https://arxiv.org/html/2607.01660#S2.E5)\),𝒖eq=\{𝒄p,𝝅p−1\}\\boldsymbol\{u\}\_\{\\mathrm\{eq\}\}=\\\{\\boldsymbol\{c\}\_\{p\},\\boldsymbol\{\\pi\}\_\{p\-1\}\\\}, and the equivalent observation refers to the Gaussian backward messageΔ𝝅ˇ→g𝝅ˇt\\Delta\_\{\\check\{\\boldsymbol\{\\pi\}\}\\rightarrow g\_\{\\check\{\\boldsymbol\{\\pi\}\}\}^\{t\}\}\. In submodule𝝅\\boldsymbol\{\\pi\},ffis defined by \([7](https://arxiv.org/html/2607.01660#S2.E7)\),𝒖eq=\{𝒛p,𝝅p−1,𝝅~p\}\\boldsymbol\{u\}\_\{\\mathrm\{eq\}\}=\\\{\\boldsymbol\{z\}\_\{p\},\\boldsymbol\{\\pi\}\_\{p\-1\},\\tilde\{\\boldsymbol\{\\pi\}\}\_\{p\}\\\}, and the equivalent observation refers to the Gaussian backward messageΔ𝝅→g𝝅t\\Delta\_\{\\boldsymbol\{\\pi\}\\rightarrow g\_\{\\boldsymbol\{\\pi\}\}^\{t\}\}\. Because the observation models involve nonlinear mappings through𝒬​\(⋅\)\\mathcal\{Q\}\(\\cdot\), some message computations are intractable; we therefore introduce mean\-field VBI\. We take submodule𝝅\\boldsymbol\{\\pi\}as an example to describe the computation as detailed in Appendix[A](https://arxiv.org/html/2607.01660#A1)\.

The detailed updating rules for the submodules are omitted for brevity\. Each intra\-symbol DAMP iteration consists of an FMP and a BMP: the FMP computes equivalent priors for the submodules sequentially in the order listed in Table[III](https://arxiv.org/html/2607.01660#S3.T3), while the BMP computes equivalent observations in the reverse order\. Messages among submodules are combined via the sum\-product rule, and non\-Gaussian extrinsic messages are approximated by EP\.

After the above message passing for symbolppcompletes, the messages on the hidden state𝝅p\\boldsymbol\{\\pi\}\_\{p\}are aggregated to form the equivalent prior for the message computation of the next symbol, i\.e\.,

dˇπp=Δgπp→𝝅p​Δgy~p→𝝅p,\\check\{d\}\_\{\\pi\_\{p\}\}=\\Delta\_\{g\_\{\\pi\_\{p\}\}\\to\\boldsymbol\{\\pi\}\_\{p\}\}\\;\\Delta\_\{g\_\{\\tilde\{y\}\_\{p\}\}\\to\\boldsymbol\{\\pi\}\_\{p\}\},\(24\)where the right\-hand side collects the messages from the state transition factorgπpg\_\{\\pi\_\{p\}\}and the output factorgy~pg\_\{\\tilde\{y\}\_\{p\}\}that point to𝝅p\\boldsymbol\{\\pi\}\_\{p\}\. The entire forward message\-passing procedure is executed sequentially forp=1,…,Pp=1,\\ldots,P, propagating the accumulated impairment memory along the pilot chain within slottt\.

TABLE III:Equivalent priors and observations for IC submodules\.After message passing within the IC module, the output extrinsic message on\{𝒚p\}p=1P\\\{\\boldsymbol\{y\}\_\{p\}\\\}\_\{p=1\}^\{P\}is updated by aggregating backward messages based on the sum\-product rule:

ΔI​C→C​T=∏p=1PΔg𝒛,pt→𝒚p​Δg𝒄,pt→𝒚p×Δg𝝅~,pt→𝒚p​Δg𝒚~,pt→𝒚p,\\begin\{split\}\\Delta\_\{IC\\to CT\}&=\\prod\_\{p=1\}^\{P\}\\Delta\_\{g\_\{\\boldsymbol\{z\},p\}^\{t\}\\rightarrow\\boldsymbol\{y\}\_\{p\}\}\\,\\Delta\_\{g\_\{\\boldsymbol\{c\},p\}^\{t\}\\rightarrow\\boldsymbol\{y\}\_\{p\}\}\\\\ &\\times\\Delta\_\{g\_\{\\tilde\{\\boldsymbol\{\\pi\}\},p\}^\{t\}\\rightarrow\\boldsymbol\{y\}\_\{p\}\}\\,\\Delta\_\{g\_\{\\tilde\{\\boldsymbol\{y\}\},p\}^\{t\}\\rightarrow\\boldsymbol\{y\}\_\{p\}\},\\end\{split\}\(25\)which is sent to the CT module as the estimated observation without impairments\.

The MP\-TTBDL algorithm is summarized in Alg\.[1](https://arxiv.org/html/2607.01660#alg1)\.

Algorithm 1Online joint channel and impairments tracking algorithm1:for

t=1,2,…t=1,2,\\ldotsdo

2:Compute input prior messages

dˇst,dˇxt,dˇωt\\check\{d\}\_\{s\}^\{t\},\\check\{d\}\_\{x\}^\{t\},\\check\{d\}\_\{\\omega\}^\{t\}\. \(Initialization if

t=1t=1\)

3:Set

ΔI​C→C​T\\Delta\_\{IC\\rightarrow CT\}to a non\-informative \(flat\) Gaussian extrinsic message for the CT observation model\.

4:for

τ=1​…​τmax\\tau=1\\ldots\\tau\_\{\\max\}do

5:

∙\\bulletCT Module:

6:Update messages within CT Module via Turbo\-OAMP with input prior

dˇst,dˇxt\\check\{d\}\_\{s\}^\{t\},\\check\{d\}\_\{x\}^\{t\}and extrinsic message

ΔI​C→C​T\\Delta\_\{IC\\rightarrow CT\}\.

7:Update output extrinsic message

ΔC​T→I​C\\Delta\_\{CT\\rightarrow IC\}using \([21](https://arxiv.org/html/2607.01660#S3.E21)\)\.

8:

∙\\bulletIC Module:

9:Update messages within IC Module via DAMP with input prior

dˇωt\\check\{d\}\_\{\\omega\}^\{t\}and extrinsic message

ΔC​T→I​C\\Delta\_\{CT\\rightarrow IC\}\.

10:Update output extrinsic message

ΔI​C→C​T\\Delta\_\{IC\\rightarrow CT\}using \([25](https://arxiv.org/html/2607.01660#S3.E25)\)\.

11:endfor

12:Update the output messages to the next time slot as \([15](https://arxiv.org/html/2607.01660#S3.E15)\)–\([17](https://arxiv.org/html/2607.01660#S3.E17)\)\.

13:Update angular domain grid points based on gradient descent\.

14:endfor

### III\-DComplexity Analysis

To reduce model complexity, we exploit the localized crosstalk structure in practical massive MIMO receivers\. Strong crosstalk is typically confined to adjacent RF chains\[[36](https://arxiv.org/html/2607.01660#bib.bib75),[28](https://arxiv.org/html/2607.01660#bib.bib63)\]\. Consequently, the distorted sample on thenn\-th RF chain depends only on𝒚𝒮n\\boldsymbol\{y\}\_\{\\mathcal\{S\}\_\{n\}\}, where𝒮n\\mathcal\{S\}\_\{n\}is the index set of the RF chains that have strong crosstalk with thenn\-th chain\. The RGRU can therefore be decomposed intoNNparallel sub\-networksGRUn\\mathrm\{GRU\}\_\{n\}, each with inputℛ​\(𝒚𝒮n\)\\mathcal\{R\}\(\\boldsymbol\{y\}\_\{\\mathcal\{S\}\_\{n\}\}\)and output inℝ2\\mathbb\{R\}^\{2\}, as conceptually illustrated in Fig\.[4](https://arxiv.org/html/2607.01660#S3.F4), where the weight matrices and the hidden state𝝅p\\boldsymbol\{\\pi\}\_\{p\}are omitted for simplicity\. Equivalently, constraining𝑾z​\(t\)\\boldsymbol\{W\}\_\{z\}\(t\),𝑾c​\(t\)\\boldsymbol\{W\}\_\{c\}\(t\),𝑾π~​\(t\)\\boldsymbol\{W\}\_\{\\tilde\{\\pi\}\}\(t\), and𝑾y~​\(t\)\\boldsymbol\{W\}\_\{\\tilde\{y\}\}\(t\)to be block diagonal reduces the model parameters from𝒪​\(N2​Nsub2\)\\mathcal\{O\}\(N^\{2\}N\_\{\\mathrm\{sub\}\}^\{2\}\)to𝒪​\(N​Nsub2\)\\mathcal\{O\}\(NN\_\{\\mathrm\{sub\}\}^\{2\}\)without loss of fitting performance\. In this parallel arrangement, each sub\-networkGRUn\\mathrm\{GRU\}\_\{n\}maintains its own hidden state; letNsubN\_\{\\mathrm\{sub\}\}denote its sub\-network hidden dimension\. The collection of these local states constitutes the overall memory of the impairment model\.

![Refer to caption](https://arxiv.org/html/2607.01660v1/x4.png)Figure 4:Factor graph of the parallel sub\-network structure, where weights and hidden states are omitted for simplicity\.We compare the model sizes of the proposed RGRU with the two conventional compensators evaluated in Sec\.[IV](https://arxiv.org/html/2607.01660#S4): a GMP compensator and a GRU compensator that adopts the same parallel per\-chain sub\-network architecture as the RGRU\. From \([2](https://arxiv.org/html/2607.01660#S2.E2)\) and \([3](https://arxiv.org/html/2607.01660#S2.E3)\)–\([6](https://arxiv.org/html/2607.01660#S2.E6)\), the proposed RGRU satisfies\|𝝎\|=∑n=1N\(3​Nsub2\+6​Nsub​\|𝒮n\|\+5​Nsub\+2\)∼𝒪​\(N​Nsub2\)\|\\boldsymbol\{\\omega\}\|=\\sum\_\{n=1\}^\{N\}\\bigl\(3N\_\{\\mathrm\{sub\}\}^\{2\}\+6N\_\{\\mathrm\{sub\}\}\|\\mathcal\{S\}\_\{n\}\|\+5N\_\{\\mathrm\{sub\}\}\+2\\bigr\)\\sim\\mathcal\{O\}\(NN\_\{\\mathrm\{sub\}\}^\{2\}\)\. The parallel GRU compensator reuses the same per\-chain sub\-network architecture and therefore has the same\|𝝎\|\|\\boldsymbol\{\\omega\}\|\. For the GMP compensator, each chain builds GMP bases with memory tapsm=0,…,Mm=0,\\ldots,Mand odd nonlinearity orders\{1,3,…,K\}\\\{1,3,\\ldots,K\\\}, givingNb,n=\(M\+1\)​⌈K/2⌉​\(2​\|𝒮n\|−1\)N\_\{b,n\}=\(M\+1\)\\lceil K/2\\rceil\\,\(2\|\\mathcal\{S\}\_\{n\}\|\-1\)basis functions per chain; with widely\-linear fitting, the real coefficient count is\|𝜷\|=4​∑n=1NNb,n∼𝒪​\(N​K​M\)\|\\boldsymbol\{\\beta\}\|=4\\sum\_\{n=1\}^\{N\}N\_\{b,n\}\\sim\\mathcal\{O\}\(NKM\)\.

We next compare the per\-slot computational complexity of Alg\. 1 with the GMP and GRU compensators\. All schemes share the same Turbo\-OAMP channel estimator, which requiresτmax\\tau\_\{\\max\}internal iterations between sub\-modules A and B and contributes𝒪​\(τmax​P​N2\)\\mathcal\{O\}\(\\tau\_\{\\max\}PN^\{2\}\)operations per slot\. Under the parallel per\-chain architecture in Fig\.[4](https://arxiv.org/html/2607.01660#S3.F4), the proposed MP\-TTBDL additionally executes the IC module \(DAMP\) in each of theτmax\\tau\_\{\\max\}CT–IC turbo iterations, introducing an extra𝒪​\(τmax​P​N​Nsub2\+τmax​P​N​Nsub\)\\mathcal\{O\}\(\\tau\_\{\\max\}PNN\_\{\\mathrm\{sub\}\}^\{2\}\+\\tau\_\{\\max\}PNN\_\{\\mathrm\{sub\}\}\)per slot, plus𝒪​\(N​Nsub2\)\\mathcal\{O\}\(NN\_\{\\mathrm\{sub\}\}^\{2\}\)for the cross\-slot message updates, yielding𝒞TTBDL=𝒪​\(τmax​P​N2\+τmax​P​N​Nsub2\+τmax​P​N​Nsub\+N​Nsub2\)\\mathcal\{C\}\_\{\\mathrm\{TTBDL\}\}=\\mathcal\{O\}\(\\tau\_\{\\max\}PN^\{2\}\+\\tau\_\{\\max\}PNN\_\{\\mathrm\{sub\}\}^\{2\}\+\\tau\_\{\\max\}PNN\_\{\\mathrm\{sub\}\}\+NN\_\{\\mathrm\{sub\}\}^\{2\}\)\. For the baselines, the GMP compensator applies a fixed polynomial transform costing𝒪​\(P​N​K​pmd\)\\mathcal\{O\}\(PNKp\_\{\\mathrm\{md\}\}\)and the GRU compensator uses a feedforward pass of complexity𝒪​\(P​N​Nsub2\+P​N​Nsub\)\\mathcal\{O\}\(PNN\_\{\\mathrm\{sub\}\}^\{2\}\+PNN\_\{\\mathrm\{sub\}\}\), each followed by the same Turbo\-OAMP stage, giving𝒞GMP=𝒪​\(P​N​K​pmd\+τmax​P​N2\)\\mathcal\{C\}\_\{\\mathrm\{GMP\}\}=\\mathcal\{O\}\(PNKp\_\{\\mathrm\{md\}\}\+\\tau\_\{\\max\}PN^\{2\}\)and𝒞GRU=𝒪​\(P​N​Nsub2\+P​N​Nsub\+τmax​P​N2\)\\mathcal\{C\}\_\{\\mathrm\{GRU\}\}=\\mathcal\{O\}\(PNN\_\{\\mathrm\{sub\}\}^\{2\}\+PNN\_\{\\mathrm\{sub\}\}\+\\tau\_\{\\max\}PN^\{2\}\)\. AsNNgrows large whilePP,τmax\\tau\_\{\\max\}, andNsubN\_\{\\mathrm\{sub\}\}remain fixed, the per\-slot complexity of all three schemes asymptotes to𝒪​\(τmax​P​N2\)\\mathcal\{O\}\(\\tau\_\{\\max\}PN^\{2\}\)\.

## IVPerformance Evaluation

In this section, we present simulation results to evaluate the performance of the proposed algorithm\. We first conduct offline pretraining to demonstrate the superiority of our message\-passing\-based neural network training algorithm over conventional methods\. We then consider the case where the online impairments exactly match the offline pretraining configuration, and verify the effectiveness of the proposed joint Bayesian deep learning framework\. Finally, we evaluate the scenario where the online impairments drift slowly over time, confirming both the effectiveness of the two\-timescale Markov transition model and the robustness of the MP\-TTBDL algorithm design\.

### IV\-AOffline Pretraining

The receiver is equipped withN=64N=64antennas in the following simulations\. To generate representative training and test data, we synthesize the distorted observations by cascading a linear crosstalk model, a nonlinear amplifier with memory, and a phase\-only IQ imbalance model, all of which are widely adopted in the literature for hardware impairment simulation\. Crosstalk is assumed only between adjacent RF chains\[[36](https://arxiv.org/html/2607.01660#bib.bib75),[28](https://arxiv.org/html/2607.01660#bib.bib63)\], with linear coupling coefficientϵ=0\.1778\\epsilon=0\.1778\(−15\-15dB\)\. The IQ imbalance is modeled as a phase mismatchϕ=0\.1778\\phi=0\.1778\(−15\-15dB\)\. The amplifier nonlinearity follows a GMP with memory depthpmd=4p\_\{\\mathrm\{md\}\}=4and nonlinearity orderK=3K=3, whose coefficients are taken from a measured\-data example in MATLAB\[[19](https://arxiv.org/html/2607.01660#bib.bib89)\]\. These specific parametric models are used only to generate evaluation data and are not assumed by the proposed method, which can handle more general hardware impairments\.

We generate10,00010\{,\}000independent sequences of 10 symbols, randomly split into5,0005\{,\}000training and5,0005\{,\}000test samples\. For each symbol, the entries of the ideal received vector𝒚p∈ℂN\\boldsymbol\{y\}\_\{p\}\\in\\mathbb\{C\}^\{N\}are independently drawn from𝒞​𝒩​\(0,1\)\\mathcal\{CN\}\(0,1\)\.

To demonstrate the advantages of our expectation\-maximization turbo deep approximate message passing \(EM\-TDAMP\) algorithm\[[33](https://arxiv.org/html/2607.01660#bib.bib71)\]in neural network training, we conducted simulations with two sub\-network hidden dimensions,Nsub=32N\_\{\\mathrm\{sub\}\}=32andNsub=64N\_\{\\mathrm\{sub\}\}=64, for each sub\-network shown in Fig\.[4](https://arxiv.org/html/2607.01660#S3.F4)\. We used a batch size of100100and trained for55epochs\. For the Adam baseline\[[12](https://arxiv.org/html/2607.01660#bib.bib37)\], the initial learning rate was set to0\.010\.01\. The fitting performance is presented in Fig\.[5](https://arxiv.org/html/2607.01660#S4.F5), which clearly shows that our algorithm achieves faster convergence and superior fitting accuracy than Adam under both configurations\.

![Refer to caption](https://arxiv.org/html/2607.01660v1/x5.png)

\(a\)Nsub=64N\_\{\\mathrm\{sub\}\}=64

![Refer to caption](https://arxiv.org/html/2607.01660v1/x6.png)

\(b\)Nsub=32N\_\{\\mathrm\{sub\}\}=32

Figure 5:Fitting performance comparison of EM\-TDAMP and Adam forNsub=64N\_\{\\mathrm\{sub\}\}=64andNsub=32N\_\{\\mathrm\{sub\}\}=32\.Our method performs Bayesian inference on the top\-level factor graph during online tracking and requires no separately trained compensator\. For fair comparison, we train two conventional compensators on the same pretraining dataset: a GMP compensator \(nonlinearity orderK=3K=3, memory depthpmd=4p\_\{\\mathrm\{md\}\}=4\) with15,16815\{,\}168real coefficients obtained by least squares, and a GRU compensator that adopts the same parallel per\-chain sub\-network architecture withNsub=32N\_\{\\mathrm\{sub\}\}=32, giving a total of243,456243\{,\}456real parameters identical to the proposed RGRU\.

### IV\-BOnline Joint Tracking

For online joint tracking, we used the same RGRU ensemble as in offline pretraining, withNsub=32N\_\{\\mathrm\{sub\}\}=32per sub\-network\. The time\-varying flat\-fading uplink channel𝒉​\(t\)\\boldsymbol\{h\}\(t\)was generated by the QuaDRiGa ray\-tracing simulator\[[10](https://arxiv.org/html/2607.01660#bib.bib78)\]under the 3GPP TR 38\.901 urban macrocell non\-line\-of\-sight \(NLOS\) scenario, with carrier frequency66GHz, two propagation clusters, and per\-cluster azimuth/elevation angular spreads of1∘1^\{\\circ\}\. The user speed was set to11m/s, withP=4P=4pilot symbols per slot and a slot duration of11ms\.

The proposed method is compared against the following benchmark schemes:

- •MP\-BDL:joint inference with a static𝝎\\boldsymbol\{\\omega\}prior; used to validate the necessity of the slow Markov prior for tracking impairment drift\.
- •MP\-BDL\-Frozen:CT–IC message passing with frozen𝝎\\boldsymbol\{\\omega\}; used to confirm the benefit of online impairment adaptation\.
- •w/o compensation:Turbo\-OAMP is applied directly to the distorted observations𝒚~\\tilde\{\\boldsymbol\{y\}\}to quantify the degradation caused by hardware impairments\.
- •GMP for compensation:an offline\-trained GMP compensator is followed by Turbo\-OAMP to demonstrate the advantage over a conventional polynomial compensator\.
- •GRU for compensation:an offline\-trained parallel GRU compensator of the same model size is followed by Turbo\-OAMP to assess the gain of the proposed joint Bayesian framework over a decoupled approach that performs compensation and channel estimation separately\.

We first consider the scenario where the online impairments are identical to those used in offline pretraining\. Since no impairment drift occurs, only MP\-BDL\-Frozen and the conventional compensators are compared, so as to validate the intra\-slot joint Bayesian inference framework\. The channel estimation NMSE is averaged over100100consecutive slots withτmax=5\\tau\_\{\\max\}=5turbo iterations between the CT and IC modules\. Fig\.[6](https://arxiv.org/html/2607.01660#S4.F6)reports the channel NMSE versus SNR for all schemes, along with an upper bound obtained from impairment\-free observations and a lower bound corresponding to the uncompensated case\. Hardware impairments are seen to severely degrade channel estimation accuracy\. Although conventional compensators mitigate this degradation to some extent, a considerable gap to the impairment\-free upper bound persists\. By contrast, the proposed MP\-BDL\-Frozen achieves the lowest channel NMSE among all compared schemes and closely approaches the upper bound across the entire SNR range, with the advantage being especially pronounced at high SNR\.

To reveal the source of this gain, Fig\.[7](https://arxiv.org/html/2607.01660#S4.F7)plots the NMSE of the recovered signal𝐲\\mathbf\{y\}versus the turbo iteration indexτ\\tau, where SNR=20=20dB is taken as a representative case\. Under the joint Bayesian inference framework, the IC module reconstructs𝐲\\mathbf\{y\}by exploiting the prior information supplied by the CT module\. As the turbo iterations proceed, the two modules progressively refine their extrinsic messages, yielding an estimate of𝐲\\mathbf\{y\}that is more accurate than that produced by the baseline compensators and, in turn, a lower channel estimation error\. This gain becomes more prominent at high SNR, where the prior variance from the CT module is smaller and the prior information is therefore more precise, further amplifying the benefit of joint inference\. It is observed that the NMSE converges withinτmax=3\\tau\_\{\\max\}=3iterations, and this value is consequently adopted as the default in all subsequent simulations\.

![Refer to caption](https://arxiv.org/html/2607.01660v1/x7.png)Figure 6:Channel estimation NMSE versus SNR when the online impairments coincide with the offline pretraining\.![Refer to caption](https://arxiv.org/html/2607.01660v1/x8.png)Figure 7:NMSE of the estimatedy\\mathrm\{y\}versusτ\\tauat SNR=20=20dB when the online impairments coincide with the offline pretraining\.The performance under slowly time\-varying impairments is evaluated next, so as to validate the effectiveness of the proposed two\-timescale Markov transition model MP\-TTBDL algorithm\. Let𝜽​\(t\)\\boldsymbol\{\\theta\}\(t\)collect all time\-varying impairment parameters at slottt, including the coefficients for LNA, the IQ phase mismatches, and the crosstalk coupling coefficients\. Att=0t=0, the crosstalk and IQ mismatch amplitudes are set toϵn​\(0\)=ϕn​\(0\)=0\.2239\\epsilon\_\{n\}\(0\)=\\phi\_\{n\}\(0\)=0\.2239\(−13\-13dB\)\. The LNA coefficients are initialized by scaling the fixed coefficients used in the matched\-online experiment by1\.021\.02,1\.051\.05, and1\.11\.1for the linear, third\-order, and fifth\-order taps, respectively, yielding a degraded amplifier state\. Fort=1,…,600t=1,\\ldots,600, the parameters evolve according to𝜽​\(t\)=ρ​𝜽​\(t−1\)\+\(1−ρ\)​𝜽⋆\+η​\|𝜽​\(t−1\)\|2⊙𝝃​\(t\)\\boldsymbol\{\\theta\}\(t\)=\\rho\\,\\boldsymbol\{\\theta\}\(t\-1\)\+\(1\-\\rho\)\\,\\boldsymbol\{\\theta\}^\{\\star\}\+\\eta\\,\|\\boldsymbol\{\\theta\}\(t\-1\)\|^\{2\}\\odot\\boldsymbol\{\\xi\}\(t\), withρ=0\.995\\rho=0\.995,η=0\.01\\eta=0\.01, and𝝃​\(t\)∼𝒩​\(𝟎,𝑰\)\\boldsymbol\{\\xi\}\(t\)\\sim\\mathcal\{N\}\(\\boldsymbol\{0\},\\boldsymbol\{I\}\)element\-wise independent across slots\. The target𝜽⋆\\boldsymbol\{\\theta\}^\{\\star\}drives the IQ and crosstalk amplitudes from−13\-13dB toward−12\-12dB, and drives the third\- and fifth\-order LNA taps further upward by approximately2%2\\%and5%5\\%relative to their initial degraded values, while the linear LNA taps remain largely unchanged\. Over the600600slots, the parameters therefore evolve gradually from the degraded initial condition toward the more severe impairment state specified by𝜽⋆\\boldsymbol\{\\theta\}^\{\\star\}, with the small per\-slot step size governed byρ\\rhoandη\\etacapturing the slow drift characteristic of practical hardware aging and environmental variations\.

To evaluate the tracking capability under time\-varying impairments, the channel estimation NMSE and the NMSE of the recovered signal𝒚\\boldsymbol\{y\}are first examined at SNR=20=20dB over600600slots\. Fig\.[8](https://arxiv.org/html/2607.01660#S4.F8)presents the channel estimation NMSE versustt\. The MP\-BDL\-family algorithms \(MP\-BDL, MP\-BDL\-Frozen, and the proposed MP\-TTBDL\) consistently outperform the conventional schemes that apply compensation followed by Bayesian channel estimation\. Within the MP\-BDL family, MP\-TTBDL achieves an increasingly larger gain over both MP\-BDL and MP\-BDL\-Frozen as tracking proceeds, validating the effectiveness of the two\-timescale Markov prior in refining the impairment parameter estimates over time\. Fig\.[9](https://arxiv.org/html/2607.01660#S4.F9)reports the NMSE of the recovered signal𝒚\\boldsymbol\{y\}versustt\. It is observed that MP\-TTBDL produces a more accurate estimate of𝒚\\boldsymbol\{y\}than the competing schemes, confirming that the channel estimation gain observed in Fig\.[8](https://arxiv.org/html/2607.01660#S4.F8)originates from the improved recovery of the input signal𝒚\\boldsymbol\{y\}, which in turn provides higher\-quality input to the CT module\. The forward RGRU fitting NMSE, defined as‖𝐲^−𝐲~‖F2/‖𝐲~‖F2\\\|\\hat\{\\mathbf\{y\}\}\-\\tilde\{\\mathbf\{y\}\}\\\|\_\{\\mathrm\{F\}\}^\{2\}/\\\|\\tilde\{\\mathbf\{y\}\}\\\|\_\{\\mathrm\{F\}\}^\{2\}, is further examined in Fig\.[10](https://arxiv.org/html/2607.01660#S4.F10)to isolate the contribution of the online RGRU parameter update\. This metric quantifies how accurately the RGRU forward model reproduces the impaired observations\. As online pilot data accumulate, MP\-TTBDL progressively refines the RGRU parameters via inter\-slot message passing, and the fitting NMSE exhibits an overall decreasing trend over time\. MP\-BDL, which retains the inter\-slot message passing but lacks the slow Markov prior to model the gradual parameter drift, demonstrates a limited tracking capability and thus falls short of the proposed MP\-TTBDL\. The baselines without online adaptation, i\.e\., MP\-BDL\-Frozen and other baseline schemes, are unable to adjust the impairment model based on live observations and show no improvement over time\. These results confirm that the two\-timescale message passing mechanism effectively exploits live observations to enhance the fidelity of the impairment surrogate, which drives the channel estimation gain observed in Fig\.[8](https://arxiv.org/html/2607.01660#S4.F8)\.

Fig\.[11](https://arxiv.org/html/2607.01660#S4.F11)reports the channel estimation NMSE averaged over the last100100slots versus SNR\. MP\-TTBDL consistently achieves the lowest NMSE among all compared schemes, confirming its robustness\. The performance gain is particularly pronounced at high SNR, since observations at higher SNR provide more informative messages for the online RGRU parameter update and thus amplify the advantage of the two\-timescale tracking framework\.

![Refer to caption](https://arxiv.org/html/2607.01660v1/x9.png)Figure 8:Channel tracking NMSE versusttunder time\-varying impairments at SNR=20=20dB, where each data point is averaged over100100slots \(0\.10\.1s\) and the total tracking duration is600600slots \(0\.60\.6s\)\.![Refer to caption](https://arxiv.org/html/2607.01660v1/x10.png)Figure 9:NMSE of estimated𝐲\\mathbf\{y\}versustt\.![Refer to caption](https://arxiv.org/html/2607.01660v1/x11.png)Figure 10:Forward RGRU fitting NMSE versustt\.![Refer to caption](https://arxiv.org/html/2607.01660v1/x12.png)Figure 11:Channel tracking NMSE versus SNR under time\-varying impairments \(averaged over the last100100slots\)\.

## VConclusion

This paper proposed an MP\-TTBDL framework for joint channel and hardware impairment tracking in massive MIMO receivers\. An RGRU was adopted to capture the intra\-slot memory of the impairments\. Distinct Markov priors were placed on the fast\-varying sparse channel and the slow\-varying network parameters, allowing temporal messages to propagate across slots at two different timescales and thereby track both channel dynamics and hardware aging\. To handle the complex intra\-slot structure, the factor graph within each slot was split into a CT module and an IC module\. Turbo\-OAMP handled the sparse channel estimation, while a proposed DAMP procedure updated the impairment parameters, and the two modules refined their estimates by exchanging extrinsic information through EP\. Specifically, a mean\-field VBI scheme was employed to approximate the non\-Gaussian messages produced by the nonlinear observation models\. Simulation results demonstrated that the proposed framework achieves consistently lower channel estimation error than the baselines that apply compensation followed by Bayesian channel estimation across both static and slowly time\-varying impairment conditions\. This advantage is particularly pronounced at high SNR and accumulates over time during online tracking, while the per\-slot computational complexity stays within a constant factor of these baselines\.

Several directions deserve further investigation\. First, extending the framework to multi\-user scenarios, where multi\-user pilot design may be jointly optimized to further improve the overall channel estimation performance, represents a natural generalization of this work\. Second, more advanced recurrent architectures, e\.g\., Transformers, may be explored to address more complex hardware impairments expected in future transceivers\. In such scenarios, message passing algorithms for the corresponding factor graphs should be developed\. Third, since the per\-slot computational cost of the proposed framework is dominated by the online tracking of slow\-varying network parameters, low\-complexity implementations are desirable for practical deployment to achieve a favorable tradeoff between tracking performance and algorithm complexity\.

## Appendix AVBI Updates for Submodule𝝅\\boldsymbol\{\\pi\}

This appendix derives the closed\-form mean\-field VBI updates for the nonlinear submodule𝝅\\boldsymbol\{\\pi\}in Sec\.[23](https://arxiv.org/html/2607.01660#S3.E23)\. The update formulas for𝝅ˇ\\check\{\\boldsymbol\{\\pi\}\}are structurally similar and can be obtained by adapting the same procedure to its specific observation model \([5](https://arxiv.org/html/2607.01660#S2.E5)\)\. Since the observation model is element\-wise, we consider the derivation in the scalar case\. We denote byy^eq\\hat\{y\}\_\{\\mathrm\{eq\}\}the equivalent observation with noise varianceσeq2\\sigma\_\{\\mathrm\{eq\}\}^\{2\}\.

##### Forward Updates \(FMP\)

During FMP, the posterior distribution ofπp\\pi\_\{p\}is approximated by a Gaussian distribution\. Since Eq\. \([7](https://arxiv.org/html/2607.01660#S2.E7)\) is a deterministic mapping, the corresponding mean and variance are expressed based on the variational posteriorsq​\(zp\)q\(z\_\{p\}\),q​\(π~p\)q\(\\tilde\{\\pi\}\_\{p\}\), andq​\(πp−1\)q\(\\pi\_\{p\-1\}\):

μπp\\displaystyle\\mu\_\{\\pi\_\{p\}\}=\(1−𝔼​\[𝒬​\(zp\)\]\)​μπp−1\+𝔼​\[𝒬​\(zp\)\]​\(2​𝔼​\[𝒬​\(π~p\)\]−1\),\\displaystyle=\\bigl\(1\-\\mathbb\{E\}\[\\mathcal\{Q\}\(z\_\{p\}\)\]\\bigr\)\\mu\_\{\\pi\_\{p\-1\}\}\+\\mathbb\{E\}\[\\mathcal\{Q\}\(z\_\{p\}\)\]\\bigl\(2\\mathbb\{E\}\[\\mathcal\{Q\}\(\\tilde\{\\pi\}\_\{p\}\)\]\-1\\bigr\),vπp\\displaystyle v\_\{\\pi\_\{p\}\}=\(1−𝔼​\[𝒬​\(zp\)\]\)2​vπp−1\+μπp−12​Var​\[𝒬​\(zp\)\]\\displaystyle=\\bigl\(1\-\\mathbb\{E\}\[\\mathcal\{Q\}\(z\_\{p\}\)\]\\bigr\)^\{2\}v\_\{\\pi\_\{p\-1\}\}\+\\mu\_\{\\pi\_\{p\-1\}\}^\{2\}\\mathrm\{Var\}\[\\mathcal\{Q\}\(z\_\{p\}\)\]\+Var​\[𝒬​\(zp\)\]​vπp−1\+\(𝔼​\[𝒬​\(zp\)\]\)2⋅4​V​a​r​\[𝒬​\(π~p\)\]\\displaystyle\\quad\+\\mathrm\{Var\}\[\\mathcal\{Q\}\(z\_\{p\}\)\]v\_\{\\pi\_\{p\-1\}\}\+\\bigl\(\\mathbb\{E\}\[\\mathcal\{Q\}\(z\_\{p\}\)\]\\bigr\)^\{2\}\\cdot 4\\mathrm\{Var\}\[\\mathcal\{Q\}\(\\tilde\{\\pi\}\_\{p\}\)\]\+\(2​𝔼​\[𝒬​\(π~p\)\]−1\)2​Var​\[𝒬​\(zp\)\]\\displaystyle\\quad\+\\bigl\(2\\mathbb\{E\}\[\\mathcal\{Q\}\(\\tilde\{\\pi\}\_\{p\}\)\]\-1\\bigr\)^\{2\}\\mathrm\{Var\}\[\\mathcal\{Q\}\(z\_\{p\}\)\]\+4​V​a​r​\[𝒬​\(zp\)\]​Var​\[𝒬​\(π~p\)\],\\displaystyle\\quad\+4\\mathrm\{Var\}\[\\mathcal\{Q\}\(z\_\{p\}\)\]\\mathrm\{Var\}\[\\mathcal\{Q\}\(\\tilde\{\\pi\}\_\{p\}\)\],where for any Gaussian variableXX, the expectation𝔼​\[𝒬​\(X\)\]\\mathbb\{E\}\[\\mathcal\{Q\}\(X\)\]and varianceVar​\[𝒬​\(X\)\]\\mathrm\{Var\}\[\\mathcal\{Q\}\(X\)\]can be expressed in terms of the bivariate normal cumulative distribution function, and the specific expressions can be found in\[[25](https://arxiv.org/html/2607.01660#bib.bib69),[8](https://arxiv.org/html/2607.01660#bib.bib70)\]\. The forward message is𝒩​\(μπp,vπp\)\\mathcal\{N\}\(\\mu\_\{\\pi\_\{p\}\},v\_\{\\pi\_\{p\}\}\), and is computed once after the BMP procedure below converges\.

##### Backward Updates \(BMP\)

Under the mean\-field assumption, the variational distribution factorizes asq​\(zp,π~p,πp−1\)=q​\(zp\)​q​\(π~p\)​q​\(πp−1\)q\(z\_\{p\},\\tilde\{\\pi\}\_\{p\},\\pi\_\{p\-1\}\)=q\(z\_\{p\}\)q\(\\tilde\{\\pi\}\_\{p\}\)q\(\\pi\_\{p\-1\}\), with each factor constrained to a Gaussian form:q​\(zp\)=𝒩​\(zp;μzp,vzp\)q\(z\_\{p\}\)=\\mathcal\{N\}\(z\_\{p\};\\mu\_\{z\_\{p\}\},v\_\{z\_\{p\}\}\),q​\(π~p\)=𝒩​\(π~p;μπ~p,vπ~p\)q\(\\tilde\{\\pi\}\_\{p\}\)=\\mathcal\{N\}\(\\tilde\{\\pi\}\_\{p\};\\mu\_\{\\tilde\{\\pi\}\_\{p\}\},v\_\{\\tilde\{\\pi\}\_\{p\}\}\), andq​\(πp−1\)=𝒩​\(πp−1;μπp−1,vπp−1\)q\(\\pi\_\{p\-1\}\)=\\mathcal\{N\}\(\\pi\_\{p\-1\};\\mu\_\{\\pi\_\{p\-1\}\},v\_\{\\pi\_\{p\-1\}\}\), whereμzp\\mu\_\{z\_\{p\}\},vzpv\_\{z\_\{p\}\},μπ~p\\mu\_\{\\tilde\{\\pi\}\_\{p\}\},vπ~pv\_\{\\tilde\{\\pi\}\_\{p\}\},μπp−1\\mu\_\{\\pi\_\{p\-1\}\}, andvπp−1v\_\{\\pi\_\{p\-1\}\}denote the updated posterior parameters\. The equivalent priors are written as𝒩​\(⋅;μzp→gπt,vzp→gπt\)\\mathcal\{N\}\(\\cdot;\\mu\_\{z\_\{p\}\\rightarrow g\_\{\\pi\}^\{t\}\},v\_\{z\_\{p\}\\rightarrow g\_\{\\pi\}^\{t\}\}\),𝒩​\(⋅;μπ~p→gπt,vπ~p→gπt\)\\mathcal\{N\}\(\\cdot;\\mu\_\{\\tilde\{\\pi\}\_\{p\}\\rightarrow g\_\{\\pi\}^\{t\}\},v\_\{\\tilde\{\\pi\}\_\{p\}\\rightarrow g\_\{\\pi\}^\{t\}\}\), and𝒩​\(⋅;μπp−1→gπt,vπp−1→gπt\)\\mathcal\{N\}\(\\cdot;\\mu\_\{\\pi\_\{p\-1\}\\rightarrow g\_\{\\pi\}^\{t\}\},v\_\{\\pi\_\{p\-1\}\\rightarrow g\_\{\\pi\}^\{t\}\}\), respectively\. During BMP, the KL divergence between this variational distribution and the true posterior under \([7](https://arxiv.org/html/2607.01660#S2.E7)\) is minimized with respect to the variational parameters\. For a given variable, e\.g\.,πp−1\\pi\_\{p\-1\}, the optimal factor satisfies

log⁡q​\(πp−1\)=𝔼q​\(zp\)​q​\(π~p\)​\[log⁡p​\(y^eq,zp,π~p,πp−1\)\]\+const,\\log q\(\\pi\_\{p\-1\}\)=\\mathbb\{E\}\_\{q\(z\_\{p\}\)q\(\\tilde\{\\pi\}\_\{p\}\)\}\\bigl\[\\log p\(\\hat\{y\}\_\{\\mathrm\{eq\}\},z\_\{p\},\\tilde\{\\pi\}\_\{p\},\\pi\_\{p\-1\}\)\\bigr\]\+\\text\{const\},where

p​\(y^eq,zp,π~p,πp−1\)\\displaystyle p\(\\hat\{y\}\_\{\\mathrm\{eq\}\},z\_\{p\},\\tilde\{\\pi\}\_\{p\},\\pi\_\{p\-1\}\)=𝒩​\(y^eq;\(1−𝒬​\(zp\)\)​πp−1\+𝒬​\(zp\)​\(2​𝒬​\(π~p\)−1\),σeq2\)\\displaystyle=\\mathcal\{N\}\\bigl\(\\hat\{y\}\_\{\\mathrm\{eq\}\};\\,\(1\-\\mathcal\{Q\}\(z\_\{p\}\)\)\\pi\_\{p\-1\}\+\\mathcal\{Q\}\(z\_\{p\}\)\(2\\mathcal\{Q\}\(\\tilde\{\\pi\}\_\{p\}\)\-1\),\\sigma\_\{\\mathrm\{eq\}\}^\{2\}\\bigr\)×Δzp→gπt​Δπ~p→gπt​Δπp−1→gπt,\\displaystyle\\times\\Delta\_\{z\_\{p\}\\rightarrow g\_\{\\pi\}^\{t\}\}\\,\\Delta\_\{\\tilde\{\\pi\}\_\{p\}\\rightarrow g\_\{\\pi\}^\{t\}\}\\,\\Delta\_\{\\pi\_\{p\-1\}\\rightarrow g\_\{\\pi\}^\{t\}\},\(26\)Based on this expression, the update ofπp−1\\pi\_\{p\-1\}yields a Gaussian distribution whose mean and variance are obtained in closed form, whereas those ofzpz\_\{p\}andπ~p\\tilde\{\\pi\}\_\{p\}require linearization of𝒬​\(⋅\)\\mathcal\{Q\}\(\\cdot\)around the current variational means\. The three factors are updated alternately until convergence, and the specific updating rules are provided below\.

- •Update forπp−1\\pi\_\{p\-1\}: vπp−1\\displaystyle v\_\{\\pi\_\{p\-1\}\}=\(1vπp−1→gπt\+1−2​𝔼​\[𝒬​\(zp\)\]\+𝔼​\[𝒬​\(zp\)2\]σeq2\)−1,\\displaystyle=\\left\(\\frac\{1\}\{v\_\{\\pi\_\{p\-1\}\\rightarrow g\_\{\\pi\}^\{t\}\}\}\+\\frac\{1\-2\\mathbb\{E\}\[\\mathcal\{Q\}\(z\_\{p\}\)\]\+\\mathbb\{E\}\[\\mathcal\{Q\}\(z\_\{p\}\)^\{2\}\]\}\{\\sigma\_\{\\mathrm\{eq\}\}^\{2\}\}\\right\)^\{\-1\},μπp−1\\displaystyle\\mu\_\{\\pi\_\{p\-1\}\}=vπp−1\(μπp−1→gπtvπp−1→gπt\+1σeq2\(y^eq\(1−𝔼\[𝒬\(zp\)\]\)\\displaystyle=v\_\{\\pi\_\{p\-1\}\}\\Biggl\(\\frac\{\\mu\_\{\\pi\_\{p\-1\}\\rightarrow g\_\{\\pi\}^\{t\}\}\}\{v\_\{\\pi\_\{p\-1\}\\rightarrow g\_\{\\pi\}^\{t\}\}\}\+\\frac\{1\}\{\\sigma\_\{\\mathrm\{eq\}\}^\{2\}\}\\bigl\(\\hat\{y\}\_\{\\mathrm\{eq\}\}\\bigl\(1\-\\mathbb\{E\}\[\\mathcal\{Q\}\(z\_\{p\}\)\]\\bigr\)−\(𝔼\[𝒬\(zp\)\]−𝔼\[𝒬\(zp\)2\]\)\(2𝔼\[𝒬\(π~p\)\]−1\)\)\)\.\\displaystyle\\quad\-\\bigl\(\\mathbb\{E\}\[\\mathcal\{Q\}\(z\_\{p\}\)\]\-\\mathbb\{E\}\[\\mathcal\{Q\}\(z\_\{p\}\)^\{2\}\]\\bigr\)\\bigl\(2\\mathbb\{E\}\[\\mathcal\{Q\}\(\\tilde\{\\pi\}\_\{p\}\)\]\-1\\bigr\)\\bigr\)\\Biggr\)\.
- •Update forzpz\_\{p\}:Since the variational meanμzp\\mu\_\{z\_\{p\}\}varies slowly across successive iterations, a first\-order Taylor expansion aroundμzp\\mu\_\{z\_\{p\}\}is sufficiently accurate, i\.e\., 𝒬​\(zp\)≈𝒬​\(μzp\)\+ϕ​\(μzp\)​\(zp−μzp\),\\mathcal\{Q\}\(z\_\{p\}\)\\approx\\mathcal\{Q\}\(\\mu\_\{z\_\{p\}\}\)\+\\phi\(\\mu\_\{z\_\{p\}\}\)\(z\_\{p\}\-\\mu\_\{z\_\{p\}\}\),whereϕ​\(⋅\)\\phi\(\\cdot\)denotes the standard normal probability density function\. With this approximation, the posterior mean and variance ofzpz\_\{p\}are given by vzp=\(1vzp→gπt\+ϕ​\(μzp\)2σeq2​𝔼​\[\(2​𝒬​\(π~p\)−1−πp−1\)2\]\)−1\\displaystyle v\_\{z\_\{p\}\}=\\Bigl\(\\frac\{1\}\{v\_\{z\_\{p\}\\\!\\rightarrow\\\!g\_\{\\pi\}^\{t\}\}\}\+\\frac\{\\phi\(\\mu\_\{z\_\{p\}\}\)^\{2\}\}\{\\sigma\_\{\\mathrm\{eq\}\}^\{2\}\}\\,\\mathbb\{E\}\\bigl\[\(2\\mathcal\{Q\}\(\\tilde\{\\pi\}\_\{p\}\)\-1\-\\pi\_\{p\-1\}\)^\{2\}\\bigr\]\\Bigr\)^\{\-1\}μzp=vzp​μzp→gπtvzp→gπt\+vzp​ϕ​\(μzp\)σeq2​\(𝒬​\(μzp\)−ϕ​\(μzp\)​μzp\)\\displaystyle\\mu\_\{z\_\{p\}\}=v\_\{z\_\{p\}\}\\,\\frac\{\\mu\_\{z\_\{p\}\\\!\\rightarrow\\\!g\_\{\\pi\}^\{t\}\}\}\{v\_\{z\_\{p\}\\\!\\rightarrow\\\!g\_\{\\pi\}^\{t\}\}\}\+v\_\{z\_\{p\}\}\\,\\frac\{\\phi\(\\mu\_\{z\_\{p\}\}\)\}\{\\sigma\_\{\\mathrm\{eq\}\}^\{2\}\}\\bigl\(\\mathcal\{Q\}\(\\mu\_\{z\_\{p\}\}\)\-\\phi\(\\mu\_\{z\_\{p\}\}\)\\mu\_\{z\_\{p\}\}\\bigr\)×𝔼​\[\(2​𝒬​\(π~p\)−1−πp−1\)2\]\\displaystyle\\times\\mathbb\{E\}\\bigl\[\(2\\mathcal\{Q\}\(\\tilde\{\\pi\}\_\{p\}\)\-1\-\\pi\_\{p\-1\}\)^\{2\}\\bigr\]\+vzp​ϕ​\(μzp\)σeq2​\(y^eq−𝔼​\[πp−1\]\)×\(2​𝔼​\[𝒬​\(π~p\)\]−1\)\\displaystyle\+v\_\{z\_\{p\}\}\\,\\frac\{\\phi\(\\mu\_\{z\_\{p\}\}\)\}\{\\sigma\_\{\\mathrm\{eq\}\}^\{2\}\}\\bigl\(\\hat\{y\}\_\{\\mathrm\{eq\}\}\-\\mathbb\{E\}\[\\pi\_\{p\-1\}\]\\bigr\)\\times\\bigl\(2\\mathbb\{E\}\[\\mathcal\{Q\}\(\\tilde\{\\pi\}\_\{p\}\)\]\-1\\bigr\)
- •Update forπ~p\\tilde\{\\pi\}\_\{p\}:Similarly, linearizing2​𝒬​\(π~p\)−12\\mathcal\{Q\}\(\\tilde\{\\pi\}\_\{p\}\)\-1aroundμπ~\\mu\_\{\\tilde\{\\pi\}\}gives vπ~p\\displaystyle v\_\{\\tilde\{\\pi\}\_\{p\}\}=\(1vπ~p→gπt\+𝔼​\[𝒬​\(zp\)2\]​\(2​ϕ​\(μπ~p\)\)2σeq2\)−1,\\displaystyle=\\left\(\\frac\{1\}\{v\_\{\\tilde\{\\pi\}\_\{p\}\\rightarrow g\_\{\\pi\}^\{t\}\}\}\+\\frac\{\\mathbb\{E\}\[\\mathcal\{Q\}\(z\_\{p\}\)^\{2\}\]\\,\(2\\phi\(\\mu\_\{\\tilde\{\\pi\}\_\{p\}\}\)\)^\{2\}\}\{\\sigma\_\{\\mathrm\{eq\}\}^\{2\}\}\\right\)^\{\-1\},μπ~p\\displaystyle\\mu\_\{\\tilde\{\\pi\}\_\{p\}\}=vπ~p\(μπ~p→gπtvπ~p→gπt\+2​ϕ​\(μπ~p\)σeq2\(y^eq−𝔼\[𝒬\(zp\)\]𝔼\[πp−1\]\\displaystyle=v\_\{\\tilde\{\\pi\}\_\{p\}\}\\Biggl\(\\frac\{\\mu\_\{\\tilde\{\\pi\}\_\{p\}\\rightarrow g\_\{\\pi\}^\{t\}\}\}\{v\_\{\\tilde\{\\pi\}\_\{p\}\\rightarrow g\_\{\\pi\}^\{t\}\}\}\+\\frac\{2\\phi\(\\mu\_\{\\tilde\{\\pi\}\_\{p\}\}\)\}\{\\sigma\_\{\\mathrm\{eq\}\}^\{2\}\}\\Bigl\(\\hat\{y\}\_\{\\mathrm\{eq\}\}\-\\mathbb\{E\}\[\\mathcal\{Q\}\(z\_\{p\}\)\]\\mathbb\{E\}\[\\pi\_\{p\-1\}\]−\(2𝒬\(μπ~p\)−1\)𝔼\[𝒬\(zp\)\]\)\)\.\\displaystyle\\quad\-\\bigl\(2\\mathcal\{Q\}\(\\mu\_\{\\tilde\{\\pi\}\_\{p\}\}\)\-1\\bigr\)\\mathbb\{E\}\[\\mathcal\{Q\}\(z\_\{p\}\)\]\\Bigr\)\\Biggr\)\.

The three updates are performed alternately until convergence\. In our simulations, this process typically requires only two iterations\.

## References

- \[1\]\(2022\-04\)Analysis of receiver IQ\-imbalance in IM\-OFDMA uplink systems\.IEEE Commun\. Lett\.26\(4\),pp\. 917–921\.External Links:[Document](https://dx.doi.org/10.1109/LCOMM.2022.3143856)Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p1.1)\.
- \[2\]S\. Amin, P\. N\. Landin, P\. Händel, and D\. Rönnow\(2014\-04\)Behavioral modeling and linearization of crosstalk and memory effects in RF MIMO transmitters\.IEEE Trans\. Microw\. Theory Tech\.62\(4\),pp\. 810–823\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p2.1)\.
- \[3\]S\. A\. Bassam, M\. Helaoui, S\. Boumaiza, and F\. M\. Ghannouchi\(2008\-06\)Experimental study of the effects of RF front\-end imperfection on MIMO transmitter performance\.In2008 IEEE MTT\- Int\. Microw\. Symp\. Dig\.,pp\. 1187–1190\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p1.1)\.
- \[4\]S\. A\. Bassam, M\. Helaoui, and F\. M\. Ghannouchi\(2009\-05\)Crossover digital predistorter for the compensation of crosstalk and nonlinearity in MIMO transmitters\.IEEE Trans\. Microw\. Theory Techn\.57\(5\),pp\. 1119–1128\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p1.1)\.
- \[5\]A\. Fawzy, S\. Sun, T\. J\. Lim, Y\. Yu, C\. Yu, and Y\. Guo\(2023\-06\)Iterative learning control based digital pre\-distortion for mitigating impairments in MIMO wireless transmitters\.IEEE Trans\. Veh\. Technol\.72\(6\),pp\. 6933–6947\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p2.1),[§II\-B](https://arxiv.org/html/2607.01660#S2.SS2.p2.6)\.
- \[6\]D\. Gao, Q\. Guo, M\. Jin, Z\. Yuan, G\. Liao, W\. Li, and Y\. Wu\(2026\)Bayesian joint nonlinear system model learning, sensing and signal detection in ISAC with hardware imperfections\.IEEE Trans\. Wirel\. Commun\.25,pp\. 5685–5700\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p5.1)\.
- \[7\]D\. Gao, Q\. Guo, G\. Liao, Y\. C\. Eldar, Y\. Li, Y\. Yu, and B\. Vucetic\(2024\-01\)Signal detection in MIMO systems with hardware imperfections: message passing on neural networks\.IEEE Trans\. Wirel\. Commun\.23\(1\),pp\. 820–834\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p5.1)\.
- \[8\]A\. Genz\(2004\)Numerical computation of rectangular bivariate and trivariate normal and t probabilities\.Statistics and Computing14\(3\),pp\. 251–260\.External Links:[Document](https://dx.doi.org/10.1023/B%3ASTCO.0000035304.20635.31)Cited by:[Appendix A](https://arxiv.org/html/2607.01660#A1.SS0.SSS0.Px1.p1.8)\.
- \[9\]M\. Goffioul, P\. Wambacq, G\. Vandersteen, and S\. Donnay\(2002\-03\)Analysis of nonlinearities in RF front\-end architectures using a modified Volterra series approach\.InProc\. Design, Automation and Test in Europe Conf\. Exhib\. \(DATE\),pp\. 352–356\.External Links:[Document](https://dx.doi.org/10.1109/DATE.2002.998297)Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p2.1)\.
- \[10\]S\. Jaeckel, L\. Raschkowski, M\. Rürup, L\. Thiele, K\. Borner, K\. Zwikiera, S\. Fries, J\. Jemai, C\. Priebe, D\. Hampicke,et al\.\(2019\)QuaDRiGa: a 3\-d multi\-layer channel model for multi\-user massive MIMO\.IEEE Access7,pp\. 168171–168185\.Cited by:[§IV\-B](https://arxiv.org/html/2607.01660#S4.SS2.p1.7)\.
- \[11\]P\. Jaraut, M\. Rawat, and F\. M\. Ghannouchi\(2018\-11\)Composite neural network digital predistortion model for joint mitigation of crosstalk, I/Q imbalance, nonlinearity in MIMO transmitters\.IEEE Trans\. Microw\. Theory Tech\.66\(11\),pp\. 5011–5020\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p2.1),[§II\-B](https://arxiv.org/html/2607.01660#S2.SS2.p2.6)\.
- \[12\]D\. P\. Kingma and J\. Ba\(2015\)Adam: a method for stochastic optimization\.In3rd Int\. Conf\. Learn\. Represent\. \(ICLR\),Cited by:[§IV\-A](https://arxiv.org/html/2607.01660#S4.SS1.p3.5)\.
- \[13\]M\. Kong, J\. Shi, B\. Sang, J\. Ding, K\. Wang, W\. Li, F\. Wang, C\. Liu, Y\. Wang, Y\. Wei, L\. Zhao, W\. Zhou, X\. Xin, B\. Liu, B\. Ye, W\. Chen, and J\. Yu\(2022\)800\-Gb/s/carrier WDM coherent transmission over 2000 km based on truncated PS\-64QAM utilizing MIMO Volterra equalizer\.IEEE J\. Lightw\. Technol\.40\(9\),pp\. 2789–2797\.External Links:[Document](https://dx.doi.org/10.1109/JLT.2022.3148336)Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p2.1)\.
- \[14\]J\. Kwon, S\. Jeon, Y\. Jeon, and H\. V\. Poor\(2024\-06\)MIMO detection under hardware impairments: learning with noisy labels\.IEEE Trans\. Wirel\. Commun\.23\(6\),pp\. 6015–6029\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p3.1)\.
- \[15\]J\. Li, C\. You, C\. Zhou, Y\. Zeng, and Z\. Feng\(2026\)Near\-field target localization: effect of hardware impairments\.IEEE Trans\. Wirel\. Commun\.\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p2.1)\.
- \[16\]Z\. X\. Li and R\. Lv\(2024\)High\-relative\-bandwidth multiband LNA distortion elimination with neural network\.IEEE Access12,pp\. 81398–81409\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p2.1),[§II\-B](https://arxiv.org/html/2607.01660#S2.SS2.p2.6)\.
- \[17\]L\. Lian, A\. Liu, and V\. K\. N\. Lau\(2019\-04\)Exploiting dynamic sparsity for downlink FDD\-massive MIMO channel tracking\.IEEE Trans\. Signal Process\.67\(8\),pp\. 2007–2021\.Cited by:[§II\-A](https://arxiv.org/html/2607.01660#S2.SS1.p2.1),[§II\-C](https://arxiv.org/html/2607.01660#S2.SS3.p2.7),[§III\-A](https://arxiv.org/html/2607.01660#S3.SS1.p2.17),[§III\-C1](https://arxiv.org/html/2607.01660#S3.SS3.SSS1.p1.7)\.
- \[18\]A\. Liu, L\. Lian, V\. K\. N\. Lau, and X\. Yuan\(2018\-12\)Downlink channel estimation in multiuser massive MIMO with hidden markovian sparsity\.IEEE Trans\. Signal Process\.66\(24\),pp\. 6436–6451\.External Links:[Document](https://dx.doi.org/10.1109/TSP.2018.2862420)Cited by:[§III\-C1](https://arxiv.org/html/2607.01660#S3.SS3.SSS1.p1.7)\.
- \[19\]MathWorks\(2024\)Generalized memory polynomial power amplifier model fitting\.Note:MATLAB Example\[Online\]\. Available:[https://www\.mathworks\.com/help/simrf/ug/generalize\-memory\-polynomial\-power\-amplifier\-model\-fitting\.html](https://www.mathworks.com/help/simrf/ug/generalize-memory-polynomial-power-amplifier-model-fitting.html)Cited by:[§IV\-A](https://arxiv.org/html/2607.01660#S4.SS1.p1.7)\.
- \[20\]M\. H\. Moghaddam, S\. Rezaei Aghdam, N\. Mazzali, and T\. Eriksson\(2022\-02\)Statistical modeling and analysis of power amplifier nonlinearities in communication systems\.IEEE Trans\. Commun\.70\(2\),pp\. 822–835\.External Links:[Document](https://dx.doi.org/10.1109/TCOMM.2021.3139630)Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p1.1)\.
- \[21\]D\.R\. Morgan, Z\. Ma, J\. Kim, M\.G\. Zierdt, and J\. Pastalan\(2006\-10\)A generalized memory polynomial model for digital predistortion of RF power amplifiers\.IEEE Trans\. Signal Process\.54\(10\),pp\. 3852–3860\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p2.1),[§II\-B](https://arxiv.org/html/2607.01660#S2.SS2.p2.6)\.
- \[22\]B\. O’Brien, J\. Dooley, A\. Zhu, and T\. J\. Brazil\(2006\-09\)Estimation of memory length for RF power amplifier behavioral models\.In2006 Eur\. Microw\. Conf\.,pp\. 680–682\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p2.1)\.
- \[23\]J\. T\. Parker, P\. Schniter, and V\. Cevher\(2014\-11\)Bilinear generalized approximate message passing—part I: derivation\.IEEE Trans\. Signal Process\.62\(22\),pp\. 5839–5853\.External Links:[Document](https://dx.doi.org/10.1109/TSP.2014.2357776)Cited by:[§III\-C2](https://arxiv.org/html/2607.01660#S3.SS3.SSS2.p2.14)\.
- \[24\]J\. Qi and S\. Aissa\(2011\-03\)Analysis and compensation for the joint effects of HPA nonlinearity, I/Q imbalance and crosstalk in MIMO beamforming systems\.In2011 IEEE Wirel\. Commun\. Netw\. Conf\.,Cancun, Mexico,pp\. 1562–1567\.External Links:ISBN 978\-1\-61284\-255\-4Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p2.1),[§II\-B](https://arxiv.org/html/2607.01660#S2.SS2.p1.1)\.
- \[25\]C\. E\. Rasmussen and C\. K\. I\. Williams\(2006\)Gaussian processes for machine learning\.MIT Press,Cambridge, MA\.Cited by:[Appendix A](https://arxiv.org/html/2607.01660#A1.SS0.SSS0.Px1.p1.8)\.
- \[26\]P\. Taghikhani, K\. Buisman, M\. Versleijen, J\. Perez\-Cisneros, and C\. Fager\(2022\-02\)Temperature\-dependent characterization of power amplifiers using an efficient electrothermal analysis technique\.IEEE Trans\. Microwave Theory Techn\.70\(2\),pp\. 1349–1360\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p3.1)\.
- \[27\]A\. Tarighat, R\. Bagheri, and A\. H\. Sayed\(2005\-08\)Compensation schemes and performance analysis of IQ imbalances in OFDM receivers\.IEEE Trans\. Signal Process\.53\(8\),pp\. 3257–3268\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p1.1)\.
- \[28\]T\. Uthayakumar, A\. Abdelhafiz, X\. Wang, and M\. Jian\(2023\-05\)Scalable single\-input behavioral modeling architecture for MIMO systems with crosstalk\.IEEE Trans\. Circuits Syst\. Regul\. Pap\.70\(5\),pp\. 2179–2188\.Cited by:[§II\-B](https://arxiv.org/html/2607.01660#S2.SS2.p1.1),[§III\-D](https://arxiv.org/html/2607.01660#S3.SS4.p1.17),[§IV\-A](https://arxiv.org/html/2607.01660#S4.SS1.p1.7)\.
- \[29\]Y\. Wan, G\. Liu, A\. Liu, and M\. Zhao\(2024\-06\)Robust multi\-user channel tracking scheme for 5G new radio\.IEEE Trans\. Wirel\. Commun\.23\(6\),pp\. 5878–5894\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p4.1)\.
- \[30\]H\. Wang, M\. Eleraky, B\. Abdelaziz, B\. Lin, E\. Liu, Y\. Liu, M\. Ghorbanpoor, C\. Chu, A\. Ruffino, J\. Xu,et al\.Power amplifiers performance survey 2000–present\.Note:\[Online\]\. Available:[https://ideas\.ethz\.ch/surveys/pa\-survey\.html](https://ideas.ethz.ch/surveys/pa-survey.html)Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p1.1)\.
- \[31\]C\. Wu, Y\. Lin, C\. Wu, and C\. Lee\(2023\)Deep learning\-based end\-to\-end design for OFDM systems with hardware impairments\.IEEE Open J\. Commun\. Soc\.4,pp\. 2468–2482\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p2.1)\.
- \[32\]Y\. Wu, G\. D\. Singh, M\. Beikmirza, L\. C\. N\. de Vreede, M\. Alavi, and C\. Gao\(2024\-05\)OpenDPD: an open\-source end\-to\-end learning & benchmarking framework for wideband power amplifier modeling and digital pre\-distortion\.In2024 IEEE Int\. Symp\. Circuits Syst\. ISCAS,pp\. 1–5\.External Links:2401\.08318Cited by:[§II\-B](https://arxiv.org/html/2607.01660#S2.SS2.p2.6)\.
- \[33\]W\. Xu, A\. Liu, Y\. Zhang, and V\. Lau\(2024\)Bayesian deep learning via expectation maximization and turbo deep approximate message passing\.IEEE Trans\. Signal Process\.72,pp\. 3865–3878\.Cited by:[§III\-C2](https://arxiv.org/html/2607.01660#S3.SS3.SSS2.p1.5),[§III\-C2](https://arxiv.org/html/2607.01660#S3.SS3.SSS2.p2.14),[§IV\-A](https://arxiv.org/html/2607.01660#S4.SS1.p3.5)\.
- \[34\]W\. Xu and A\. Liu\(2025\)Message passing based two\-timescale bayesian deep learning for joint channel and hardware impairments estimation\.In2025 IEEE 101st Veh\. Technol\. Conf\. \(VTC2025\-Spring\),Oslo, Norway\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p6.1)\.
- \[35\]W\. Xu, A\. Liu, M\. Zhao, and G\. Caire\(2024\)Joint visibility region detection and channel estimation for XL\-MIMO systems via alternating MAP\.IEEE Trans\. Signal Process\.72,pp\. 4827–4842\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p4.1)\.
- \[36\]L\. Zhao, H\. Qian, W\. Feng, M\. Li, and K\. Kang\(2023\-10\)A low\-complexity digital predistorter for large\-scale MIMO systems with crosstalk\.IEEE Trans\. Veh\. Technol\.72\(10\),pp\. 13203–13213\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p2.1),[§II\-B](https://arxiv.org/html/2607.01660#S2.SS2.p1.1),[§III\-D](https://arxiv.org/html/2607.01660#S3.SS4.p1.17),[§IV\-A](https://arxiv.org/html/2607.01660#S4.SS1.p1.7)\.
- \[37\]X\. Zheng and V\. Lau\(2021\-12\)Online DNN\-based channel estimator for massive MIMO systems with nonlinear distortion\.In2021 IEEE Glob\. Commun\. Conf\. GLOBECOM,Madrid, Spain,pp\. 1–6\.External Links:ISBN 978\-1\-7281\-8104\-2Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p3.1)\.
- \[38\]Y\. Zhou, A\. Liu, and W\. Xu\(2025\)Subspace constrained variational Bayesian inference for structured compressive sensing with a dynamic grid\.IEEE Trans\. Signal Process\.73,pp\. 781–794\.Cited by:[§I](https://arxiv.org/html/2607.01660#S1.p4.1)\.
- \[39\]J\. Ziniel and P\. Schniter\(2013\-11\)Dynamic compressive sensing of time\-varying signals via approximate message passing\.IEEE Trans\. Signal Process\.61\(21\),pp\. 5270–5284\.Cited by:[§II\-C](https://arxiv.org/html/2607.01660#S2.SS3.p1.2)\.