{'fraction_non_alphanumeric': 0.0418923474547702, 'fraction_numerical': 0.010203393639734801, 'mean_word_length': 4.588545591559909, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Message Passing-Based Joint User Activity Detection and Channel Estimation for Temporally-Correlated Massive Access

This paper studies the temporally-correlated massive access system where a large number of users communicate with the base station sporadically and continue transmitting data in the following frames in high probability when being active. To exploit both the sparsity and the temporal correlations in the user activities, we formulate the joint user activity detection and channel estimation problem in multiple consecutive frames as a dynamic compressed sensing (DCS) problem. Particularly, the problem is proposed to be solved under Bayesian inference to fully utilize the channel statistics and the activity evolution process. The hybrid generalized approximate message passing (HyGAMP) framework is leveraged to design a HyGAMP-DCS algorithm, which can nearly achieve the Bayesian optimality with efficient computations. Specifically, a GAMP part for channel estimation and an MP part for activity likelihood update are included in the proposed algorithm, then the extrinsic information is exchanged between them for performance enhancement. Moveover, we develop the expectation maximization HyGAMP-DCS (EM-HyGAMP-DCS) algorithm to adaptively learn the hyperparameters during the estimation procedure when the system statistics are unavailable. Particularly, the analytical tool of state evolution is provided to find the appropriate hyperparameter initialization that ensures EM-HyGAMP-DCS to achieve satisfied performance and fast convergence. From the simulation results, it is validated that our proposed algorithm can significantly outperform the existing methods.


I. INTRODUCTION
Massive machine-type communications (mMTC) is one of the three use cases of the fifthgeneration (5G) cellular networks, for Internet of Things (IoT) applications [2]. It aims to provide wireless connectivity to a massive number of IoT devices, whose traffic is typically sporadic and with short packet size [3], [4]. The main technical challenge of mMTC is to design scalable, efficient and low-latency multiple-access schemes for the massive IoT devices. Since the conventional grant-based random access schemes suffer excessive delay and signaling overhead in such scenario, the grant-free random access scheme is considered as the promising solution to both decrease the access delay and reduce the control overhead for coordination [4].
In the grant-free random access protocol, a unique pilot sequence with the usage of identification and channel estimation is pre-allocated to each device. However, due to the massive number of the IoT devices but the limited coherence time interval, the pilot sequences are usually non-orthogonal. When one device is activated, it directly transmits its dedicated pilot followed by data without waiting for the grant from the base station (BS). Meanwhile, the BS can perform joint user activity detection and channel estimation based on the received signals. Given the sporadic traffic generated from IoT devices, the joint user activity detection and channel estimation problem by nature can be cast into a large-scale sparse signal recovery problem, which is promising to be reliably solved by the compressed sensing (CS) techniques [4].
In many practical IoT systems, once a device is activated, it often transmits continuously over multiple consecutive frames. This suggests that the device activities are often temporallycorrelated. To fully exploit the temporally-correlated user activity, we are motivated to detect the active users and estimate their channels in multiple consecutive frames jointly. The joint user activity detection and channel estimation problem can be formulated from the dynamic compressed sensing (DCS) perspective [5], [6], which takes both the sporadic user activity and the temporal correlation of user activities into account. By building the probabilistic and graphical model of the signals, the joint user activity detection and channel estimation can be realized by the standard message passing (MP) algorithm [7]. To simplify the implementation of MP in large-scale systems, we propose to leverage the hybrid generalized approximate message passing (HyGAMP) [8] framework which employs a hybrid of AMP and MP for the graphical model. Since the perfect system statistics are usually unavailable, the expectation maximization (EM) technique is incorporated with the algorithm to adaptively learn the hyperparameters in the estimation procedure. We show that the proposed algorithm can achieve significantly performance improvement for the massive access system with temporally-correlated user activity.

A. Related Works
To accomplish joint user activity detection and channel estimation for the massive access networks, a variety of advanced sparse signal recovery algorithms have been proposed for different wireless communication settings. For the single-cell scenario, the works [9], [10] utilize the standard CS algorithms of orthogonal matching pursuit (OMP) and basis pursuit denoising (BPDN). However, these two kinds of algorithms usually suffer high computational complexity due to the extremely large amount of devices. Then the works [11], [12] propose the computationally efficient AMP algorithms with Bayesian denoiser for user activity detection and characterize their asymptotic performance, which can also outperform the conventional CSbased algorithms. In [13], the message-scheduling generalized AMP algorithm is developed to further reduce the computational complexity without degrading the detection and estimation performance. The AMP-based algorithms are also extended to perform data detection along with the activity detection in [14], [15]. Besides the AMP-based algorithm, covariance matching pursuit [16], [17], dimension reduction-based optimization [18] and deep learning methods [19]- [21] are also investigated for performance enhancement. In particular, the covariance matching pursuit algorithm [16], [17] and dimension reduced-based optimization algorithm [18] are especially designed for the massive MIMO scenario, and they can significantly outperform the conventional CS methods when the length of the pilot sequences is equal. The deep learning methods [19]- [21] propose to build the artificial neural networks and learn the network parameters based on the pre-collected data, which are able to outperform the traditional hand-designed algorithms and enjoy efficient computation simultaneously. Moreover, a sparsity-constrained method is proposed for the non-ideal scenario where different users have unknown frequency offsets in [22]. In the multi-cell systems, the cooperative user activity detection and channel estimation is studied in [23], [24] by employing the AMP-based algorithms. Recently, unsourced random access as a new random access protocol is also studied in [16], [25], where each user shares a common codebook and the BS only needs to detect the transmit codewords.
The aforementioned works all detect the user activity in each frame individually by assuming the user activity is independent for each frame. Due to the low-rate transmission, the information transmission of one active user may occupy multiple consecutive frames. Such temporal correla-tion of user activity can be exploited for performance enhancement. Assuming the user sparsity level and the channel state information (CSI) are available, the work [26] proposes a DCS-based algorithm where the detected active user set in the last slot is used as the initial set to be detected in the current slot. By adaptively exploiting the prior support based on the corresponding support quality information, the work [27] proposes a prior-information-aided adaptive subspace pursuit (PIA-ASP) algorithm, which can always outperform the DCS-based algorithm in [26]. These algorithms are all based on the greedy algorithms including the matrix inversion calculation in each iteration, and thus suffer high computational complexity in the large-scale problem.
Moreover, both works [26], [27] assume that the CSI is perfectly known in prior, which is unrealistic in massive access. By leveraging the system statistics, the authors in [28] propose a sequential AMP (S-AMP) algorithm to sequentially perform user activity detection and channel estimation in each frame. Then the work [29] proposes to extract the side information (SI) from the estimation results in the previous frame to enhance the user activity detection performance in the current frame based on the AMP framework.
Note that these prior works [26]- [29] all only exploit the temporal correlation by using the historical estimation information. The work [30] proposes to further exploit the temporal correlation by additionally utilizing the estimation information in the next frame, which doubles the computational cost. Compared with these works [26]- [30] that perform activity detection in the frame-by-frame way, this paper considers the multi-frame activity detection problem where the block-by-block detection is performed. Moreover, this paper propose to adaptively learn the parameters of the system statistics by combining the EM technique with the proposed algorithm, while the algorithms in [28]- [30] all work in the case with perfect system statistics. To ensure that the proposed HyGAMP-based algorithm incorporated with EM can converge to the satisfied point with fast convergence, we illustrate that the performance analysis tool of state evolution (SE) plays an important role of finding the appropriate hyperparameter initialization.

B. Main Contributions
In this work, we consider the problem of joint user activity detection and channel estimation in multiple consecutive frames for the temporally-correlated massive access systems. The main contributions and distinctions of this work are summarized as follows: 1) This paper proposes to jointly perform user activity detection and channel estimation in multiple consecutive frames, which is formulated as a DCS problem. As such, not only the sparse activity pattern is considered in the problem, but also the statistical relationships of the activities in these consecutive frames are exploited.
2) In contrast to only making use of the historical estimations in [28], [29], this paper proposes a HyGAMP-DCS algorithm to fully exploit the temporally-correlated user activities in the multiple consecutive frames. The HyGAMP-DCS algorithm can combine the computationally efficient GAMP algorithm for channel estimation and the standard MP algorithm for the likelihood update of user activities. In particular, the likelihoods of the user activities in each frame are refined by aggregating both the forward messages from the previous frame and backward messages from the next frame at each iteration.Numerical results show that HyGAMP-DCS can achieve superior performance to the traditional DCS-based algorithms.
3) To make the proposed algorithm applicable for the practical system with imperfect system statistics, this paper proposes to incorporate the EM algorithm in HyGAMP-DCS to adaptively learn the hyperparameters during the estimation procedure. Compared with the existing algorithms in [28]- [30], the proposed EM-based HyGAMP-DCS algorithm (EM-HyGAMP-DCS) can still achieve favorable performance without perfect system statistics.The simulation results also demonstrate that the proposed EM-HyGAMP-DCS algorithm can realize very similar performance to the HyGAMP-DCS algorithm which requires the perfect system statistics. 4) For performance analysis, this paper provides the SE to predict the asymptotic performance of HyGAMP-DCS and EM-HyGAMP-DCS. In particular, the SE is identified as the essential tool to find the appropriate hyperparameter initialization for EM-HyGAMP-DCS, which ensures the algorithm to achieve satisfied performance with fast convergence.Meanwhile, the computational complexity comparison with the state-of-the-art methods is also given to illustrate the computational efficiency of our proposed algorithms.

C. Organizations and Notations
The remaining part of this paper is organized as follows. Section II introduces the system model of the temporally-correlated massive access network. In Section III, we formulate the problem in the probabilistic perspective and introduce the Bayesian inference methods. In Section IV, the HyGAMP-DCS algorithm is proposed, then the EM algorithm is employed for hyperparameter update in Section V. Next, we analyze the performance and the computational complexity in A. Temporal Correlation of the User Activity As shown in Fig. 1, once a user is activated, its data transmission can occupy multiple consecutive frames. As such, there exist temporal correlations of the user activities in neighboring frames.
Similar to [6], [28]- [30], the dynamic user activity over frames is modeled as a stationary first-order Markov chain with two discrete states. The Markov chains on different devices are assumed to independent and identically distributed (i.i.d.). Let λ n,t ∈ {0, 1} indicate whether user n ∈ {1, 2, . . . , N } is active or not in the tth frame. The steady state probability is given by Pr(λ n,t = 1) = p a , where p a denotes the active probability of each user n in frame t. The transition probability matrix is defined as where p 10 = Pr{λ n,t = 0|λ n,t−1 = 1} and other three transition probabilities are defined similarly. By the steady-state assumption, the Markov chain can be completely characterized by two parameters p a and p 10 . The other three transition probabilities can be easily obtained as p 01 = p a p 10 /(1 − p a ), p 00 = 1 − p 01 , and p 11 = 1 − p 10 . Note that 1/p 10 can be considered as the expected number of the consecutive frames occupied by the data transmission of one user, which means that smaller p 10 leads to stronger temporal correlation of the user activity over frames. Although the temporal correlation can be described more precisely by higher-order Markov chains, we consider the first-order Markov chain here for simplicity.

B. Multi-Frame Dynamic Signal Model
In order to exploit the temporally correlated user activities, we consider T consecutive frames for signal transmission with T ≥ 2. . In each frame t ∈ {1, 2, . . . , T }, the transmission is divided into two phases, the pilot phase and data phase. Each user n is pre-assigned a unique pilot sequence a n = [a n,1 , a n,2 , . . . , a n,L ] T ∈ C L×1 with length L. The BS performs user identification and channel estimation based on the detected pilot sequences in the pilot phase and then decodes the data in the data phase. In this work, we concentrate on the joint user activity detection and channel estimation problem in the pilot phase. Due to the limit of orthogonal resources, the pilot length L is assumed to be much smaller than the number of users N , i.e., L N , so that the pilot sequences of different users cannot be mutually orthogonal. The pilot sequences in this work are generated from the independent and identically distributed (i.i.d.) complex Gaussian distributions with zero mean and variance of 1/L, and then normalized to have unit power.
The channel in each frame is assumed to satisfy independent Rayleigh distribution. We denote the effective channel coefficient between user n and the BS in the tth frame as h n,t = √ β n g n,t , where β n = √ P n γ n is the effective large-scale fading component decided by the transmit power P n and the large-scale attenuation γ n , g n,t ∈ C is the small-scale fading coefficient following i.i.d. circularly symmetric complex Gaussian distributions with zero mean and unit variance, i.e., g n,t ∼ CN (0, 1). We adopt the power control strategy in [14], so that the effective large-scale fading component of each user is the same, i.e., β n = β, ∀n. Then the received pilot signals at the BS in continuous T frames can be written as where A = [a 1 , a 2 , . . . , a N ] ∈ C L×N is the pilot matrix of all users; Λ = [λ λ λ 1 , λ λ λ 2 , . . . , λ λ λ T ] ∈ C N ×T is the activity matrix with λ λ λ t = [λ 1,t , λ 2,t , . . . , is the effective channel matrix defined as the Hadmard product of Λ and H with each entry x n,t = λ n,t h n,t ; W = [w 1 , w 2 , . . . , w T ] ∈ C L×T is the additive white complex Gaussian noise (AWGN) matrix where each element has zero mean and variance σ 2 w . Our objective is to jointly detect the active users and estimate their channels by recovering X from the received signal Y given the pilot matrix A. The problem can be regarded as an underdetermined linear inverse problem in the multiple measurement vector (MMV) form. Due to the sparse user activity pattern, each column of the effective channel X, i.e., x t , only has a small portion of non-zero entries, and thus the powerful CS-based algorithms can be readily applied. However, notice that the positions of the non-zeros entries in each x t slowly changes with t due to the temporally-correlated user activity. Such problem with dynamic structure of support is generally referred to as the dynamic compressed sensing (DCS) problem. In this work, we propose to exploit the temporal correlation of user activity in the DCS problem to enhance the performance, which will be detailed in the following sections.

III. PROBLEM FORMULATION
The DCS problem can be solved by optimization-based methods using convex relaxation and the standard convex problem solver, e.g., Dynamic LASSO [5]. It can also be solved by the greedy algorithms based on OMP and SP with better recovery performance by employing the detected active user set in the previous frame as the prior information to assist the detection in the current frame [26], [27]. These two approaches can outperform the traditional CS-based algorithms, but however, are unable to fully utilize the system statistics including the Markov chain-based user activity and the channel distributions. In this section, we propose to build the probabilistic model of the temporal-correlated massive access system, which enables the usage of the powerful Bayesian inference methods.

A. Probabilistic Model
By taking the system statistics into account, we represent all the variables and their relationships from the probability perspective. Since the block fading Rayleigh channel is assumed, the prior probability of the effective channel coefficient of user n in the tth frame, x n,t , can be characterized by an independent Bernoulli-Gaussian distribution as where δ(·) is the Dirac delta function. Note that the more complicated channel model can also be considered and we can approximate the channel distribution with the Gaussian mixture distribution [31]. The conditional probability p(y t |x t ) is obtained from the AWGN channel as By combining the Markov-chain modeled user activities, the probabilistic signal model describing the linear system of (2) can be finally derived as where it is noted that p(λ n,1 |λ n,0 ) = p(λ n,1 ) = (1 − p a )(1 − λ n,1 ) + p a λ n,1 and p(λ n,t |λ n,t−1 ) is given as (61) in the Appendix.

B. Bayesian Inference
Based on the probabilistic model (5), the optimal performance of user activity detection and channel estimation under the minimum mean square error (MMSE) principle can be realized by Bayesian inference. In specific, we first follow the Bayes' rule to obtain the posterior joint probability of X and Λ as where the marginal probability of the received signal Y is derived by Given the posterior marginal probability p(x n,t |Y) = Λ dΛ X \x n,t p(X, Λ|Y)dX \xn,t , the MMSE estimator for x n,t can be expressed as The factor graph representation of the joint posterior distribution p(X, Λ|Y).
The estimator (8) achieves the minimum (Bayesian) MSE defined as where the expectation is computed over the joint distribution p(X, Λ, Y) expressed in (5). Then the activity likelihood p n,t of user n in frame t can be obtained as where p(λ n,t |Y) = X dX Λ \λ n,t p(X, Λ|Y)dΛ \λn,t . By comparing (10) with a predefined threshold, the user activities can be finally decided.
Due to the fact that the number of users in massive access is very large, both the MMSE estimator (8) and the user activity likelihood calculator (10) involve very high-dimensional integrals and thus are hard to be directly obtained. In the next section, we design a computationally efficient HyGAMP-based algorithm to approximately derive (8) and (10) for our problem.

IV. THE HYGAMP-DCS ALGORITHM
Before illustrating the algorithm design, we first give the factor graph representation of the joint posterior probability p(X, Λ Λ Λ|Y) with decomposition (6) as shown in Fig. 2. Intuitively, the factor graph suggests the usage of the MP algorithm [7] to obtain (8) and (10). The standard MP algorithm iteratively updates the messages conveyed between the adjacent nodes and aggregates the messages arrived at the nodes {x n,t } and {λ n,t } to obtain their posterior probabilities, which avoids the extremely high-dimensional integrals in (8) and (10). However, the calculation of messages µ i (n,t)→(l,t) (x n,t ) and υ i (n,t)←(l,t) (x n,t ) in the dense bipartite graph including the nodes {x n,t } and {p(y l,t |x t )} is still intractable since L and N is large in massive access. Fortunately, the HyGAMP framework [8] can be leveraged to significantly simplify the implementation of the standard MP algorithm. Its key idea is to partition the edges in the factor graph into the strong and weak subsets. Then the messages passed on the weak edges are simplified by the AMP-style approximation, while the messages passed on the strong edge are still updated by standard MP principle.
By following the idea of the HyGAMP framework, we propose a HyGAMP-DCS algorithm specific to the large-scale DCS problem. Compared with the conventional HyGAMP algorithm in [8], the MP equations in HyGAMP-DCS are re-derived based on the statistical dependencies of the signals in our considered system. The HyGAMP-DCS algorithm can be divided into two parts of GAMP and MP based on the message updating rule, where these two parts perform channel estimation and activity likelihood update, respectively. By exchanging the extrinsic messages between these two parts, the performance of user activity detection and channel estimation can both be improved. In the following, the details of the message passing equations in these two part is introduced.

A. GAMP Part
Following the HyGAMP framework, the GAMP part regards the edges between the nodes {x n,t } and {p(y l,t |x t )} as the weak edges thanks to their linearizable coupling. Then the messages µ i (n,t)→(l,t) (x n,t ) and υ i (n,t)←(l,t) (x n,t ) are simplified by using GAMP approximations based on the central limit theorem and Taylor expansion.
The basic GAMP estimation is given in line 9-18 of Algorithm 1, which can incorporate arbitrary distribution on the input X and output Y [32]. Compared with the standard MP algorithm, the number of the update variables in GAMP estimation shrinks from O(LN ) to O(L+N ). Specifically, in line 10-13, the GAMP estimation performs the update of the noiseless signal Z = AX. The distribution of each element z l,t is estimated as CN (z l,t ; p l,t (i), τ p l,t (i)) by the messages from the variable nodes {x n,t }, where p l,t (i) is the plug-in estimate of z l,t with variance τ p l,t (i) in the ith iteration.By accounting the AWGN channel between Z and Y, the closed-from expression of MMSE estimation for z l,t and its variance in line 12-13 are available After updating the noiseless signal z n,t , the mean and the variance of each element in the residual signal that removes the estimate X from Y is given in line 14-15. Then line [16][17] give the update of the plug-in estimate r n,t of the true signal x n,t , which is modeled as r n,t (i) = x n,t + n n,t (i) with n n,t (i) ∼ CN (0, τ r n,t (i)).To obtain the MMSE estimator of X, we can give the approximate posterior marginal distribution of x n,t as given as (27) in the MP part.By simplifying (13), the approximate posterior marginal distribution p i (x n,t |Y) can be written in the form of a spike and slap probability as where Thus, we can obtain explicit expression of the posterior mean and variance of x n,m in the line 18 and 19 as is updated by the MP part at each iteration, the MMSE estimator (18) also changes at each iteration. However, the MMSE estimators on x n,t keep constant at each iteration in the AMP-based algorithms proposed in [28], [29] for temporally-correlated massive access.

B. MP Part
By regarding all the edges in the MP part as strong edges, we utilize the standard MP algorithm to update the activity probability of each user with the extrinsic message υ i (n,t)←(n,t) (λ n,t ) from the GAMP part.
During the forward message passing, the backward message passing can also be performed in parallel to convey the useful information from λ n,t+1 to λ n,t for t = 1, . . . , T − 1. Due to the symmetry, the update formulas for ← − q n,t (i) in the backward message ζ i (n,t)←(n,t+1) (λ n,t ) is similar to that of − → q n,t (i), except that { ← − q n,t (i)} are updated sequentially from frame (T − 1) to frame 1. For t ≤ T − 1, the backward message ζ i (n,t)←(n,t+1) (λ n,t ) can be written as where With both the forward message and backward message in the above, the extrinsic message µ i (n,t)→(n,t) (λ n,t ) to refine the channel estimation in the next iteration is subsequently given by µ i (n,t)→(n,t) (λ n,t ) = ξ i (n,t)→(n,t) (λ n,t ) · ζ i (n,t)←(n,t+1) (λ n,t ) for t = 1, . . . , T − 1.Then the update rule of − → p n,t (i) in message µ i (n,t)→(n,t) (λ n,t ) for t = 1, . . . , T − 1 is obtained as − → p n,t (i + 1) = ← − q n,t (i) − → q n,t (i) But for t = T , the update rule is given by µ i (n,t)→(n,t) (λ n,t ) = ξ i (n,t)→(n,t) (λ n,t ) and − → p n,T (i + 1) = − → q n,T (i), since the variable node λ n,T is connected to only one factor node p(λ n,T |λ n,T −1 ) in the MP part. Though S-AMP [28] has similar forward message definition in (23), there is only once message forwarding from frame (t − 1) to frame t after the estimation for frame (t − 1) finishes in S-AMP. Moreover, the backward message is not included in S-AMP, which thus ignores the temporal correlation between λ n,t and λ n,t+1 during the estimation for each frame t.
Compared with S-AMP, HyGAMP-DCS updates both the forward message from λ n,t−1 to λ n,t and the backward message from λ n,t+1 to λ n,t in each iteration. Therefore, it is promising that HyGAMP can provide a more precise estimation of the activity λ n,t , which also contributes to improved channel estimation performance.The whole procedure of the proposed HyGAMP-DCS algorithm is outlined in Algorithm 1.

C. Discussion
In the existing Bayesian inference methods of S-AMP [28] and SI-aided MMV-AMP [29], the user activity and channel estimation is performed in the frame-by-frame way, which thus results in only once unidirectional message propagation from frame (t−1) to frame t. Compared with them, HyGAMP-DCS performs multi-frame user activity detection and channel estimation, where the activities and channels in the continuous T frames are simultaneously estimated. In other words, HyGAMP-DCS can be regarded to perform block-by-block detection in temporally-correlated massive access. In particular, the proposed HyGAMP-DCS algorithm performs bidirectional message propagation between the adjacent frames in each iteration to refine the estimation in each frame, which can further exploit the temporally-correlated user activity in the T consecutive frames for performance improvement.However, this multi-frame detection scheme by HyGAMP-DCS may result in larger detection delay compared with [28], [29]. The average detection delay of HyGAMP-DCS can be easily obtained as T −1 2 (frame), but we can find that significantly performance improvement can be realized with a moderate value of T from the simulation results in Section VII.

V. HYPERPARAMETER LEARNING VIA EM ALGORITHM
The traditional AMP-based algorithms in [28]- [30] and the proposed HyGAMP-DCS algorithm usually require the prior knowledge of the system statistics characterized by the hyperparameters ϑ ϑ ϑ = {p a , β, p 10 , σ 2 w }, which may not be estimated accurately in advance. On the other hand, these statistical parameters in ϑ ϑ ϑ may be tuned based on the machine learning techniques to better match the currently received signals at the BS, which thus improves the performance. As a result, we propose to integrate the expectation maximization algorithm [33] with the proposed HyGAMP-DCS algorithm to simultaneously learns the statistical parameters from the currently received signals during the procedure of the detection and estimation on X and Λ Λ Λ. Compared with the work [31] that concentrates on the CS problem with single measurement vector (SMV), our proposed algorithm solves the DCS-SMV problem and additionally needs to learn the transition probability matrix P of the Markov chain. Moreover, the parameters of the channel statistics are learned from the received signals Y in multiple consecutive frames, which enables us to have a more accurate hyperparameter learning result.
The EM algorithm aims to find the optimal parameter setting that maximizes the likelihood function p(Y, X|ϑ ϑ ϑ), while the optimization problem is usually non-convex and has untractable computation complexity. Instead, the EM algorithm maximizes a lower bound of the likelihood function p(Y, X|ϑ ϑ ϑ) in each iteration, then p(Y, X|ϑ ϑ ϑ) is guaranteed to be increased to a local maximal point. In specific, two steps of the E-step and M-step are alternatively performed until convergence, which can be written as where p i (X) = p(X|Y, ϑ ϑ ϑ i ) is the estimated posterior marginal distribution of X based on the updated statistical parameters ϑ ϑ ϑ i and E p i (x) [·] denotes the expectation over the posterior distribution p i (X). Since the number of users N is very large, the expectation calculation may brings about unaffordable computation cost. Fortunately, the HyGAMP framework enables us to approximate the posterior probability by p i (X) = T t=1 N n=1 p(x n,t |Y, ϑ ϑ ϑ i ) in the large-scale system, which leads to a computationally efficient EM estimation procedure. After obtaining the expectation, we leverage the EM algorithm in the incremental version [34] to simplify the joint optimization problem (29), where the hyperparameters are optimized in the coordinate-wise manner. So that only one hyperparameter ϑ ∈ ϑ ϑ ϑ is updated at a time while others all keep fixed.
The simplified optimization problem can be written as The solution of the problem (30) can be easily obtained in the closed form by setting the derivative of L i (ϑ ϑ ϑ i \ϑ , ϑ) with respect to ϑ to be zero. The hyperparameters in the ith iteration are updated by After obtaining the updated probabilities p a (i+1) and p 10 (i+1), we can simply obtain the updated probability p 01 (i+1) = p a (i)p 10 (i)/(1−p a (i)), p 11 (i+1) = 1−p 10 (i+1), p 00 (i+1) = 1−p 01 (i+1) and finally complete the update of the transition probability matrix in the Markov chain.To make the work self-contained, the detailed derivations of the above equations (31) to (34) where Φ(·) and φ(·) denote the cumulative distribution function and the probability distribution function of the standard normal distribution, respectively. Since the knowledge of SNR may not be available, the appropriate initial value of SNR is in need. The conventional work usually adopt the value SNR 0 = 20dB [31], but we find that this setting can usually lead to a bad local optimal point for the EM-HyGAMP-DCS algorithm in our system setting. Though the suitable value of SNR 0 can be found by exhaustive search method, this approach needs to collected lots of data with perfectly known user activities and user channels, which is inefficient. To tackle this problem, we will firstly introduce the performance analysis tool of SE in the next section and show that the SE can be an effective tool to help find the suitable value SNR 0 for the specific system.Finally, we also outline the whole procedure of the EM-HyGAMP-DCS algorithm in Algorithm 2.

VI. PERFORMANCE AND COMPLEXITY ANALYSIS
In this section, we firstly analyze the performance of the HyGAMP-DCS algorithm by the SE in the asymptotic regime, where L, N → ∞ but L N is fixed. It is noted that the SE can not only be used to decide the setting on the pilot length L and the joint detection frame number T , but also play a important role of finding the appropriate hyperparameter initialization for the EM-HyGAMP-DCS algorithm. Then we illustrate its computational efficiency. Peform line 8-24 in Algorithm 1 with updated σ w (i), p a (i), β(i), and p 10 (i).

A. State Evolution
SE is a common framework to analyze the estimation performance of the AMP-based algorithm in the asymptotic regime. When the condition that the pilot matrix has i.i.d. sub-Gaussian elements is satisfied, the state evolution can accurately track the mean square error (MSE) of estimator X(i) in the large-scale problem. To simplify the SE analysis, we consider the proposed HyGAMP-DCS algorithm to have scalar variance in each frame t. These scalar variances of the algorithm 1 can be rewritten as then τ x n,t (i+1) can also be calculated according to (19) to update the other variables and variances in the next iteration.
Following the assumptions of the SE for the AMP-based algorithms in [32], the asymptotic MSE of the estimator x n,t (i) is identical to E[τ x t (i)] whose expectation is performed on r n,t (i). Specifically, for the random variable the variable r 0 (i) in each iteration i can be modeled as where n r 0 (i) ∼ CN (0, τ r 0 (i)) is the corrupting complex Gaussian noise being independent to x 0 . Under the proposed algorithm with scalar variance, we have where τ x 0 (i) is calculated by (19) and the expectation operates on both x 0 and n r 0 (i − 1). The explicit expression of the expectation is given as Then the asymptotic MSE of X(i) in each iteration i can be given as However, the active probability ρ 0 is updated at the MP part in each iteration and it is hard to give the probability of ρ 0 in closed form. Thus, we resort to the Monte Carlo simulations to obtain ρ n,t = − → p n,t (i) for each user n in frame t at each iteration and then substitute it into (48). Similarly, the performance of the EM-HyGAMP-DCS algorithm can also be predicted by following the above procedure where the statistical parameters ϑ ϑ ϑ in each iteration are replaced by the learned ones ϑ ϑ ϑ(i).
From the simulation trials, we can find that the EM-HyGAMP-DCS algorithm easily converges to a bad local optimal point with an inappropriate setting on SNR 0 . As usual, we can easily collect  lots of signals {y t } received at the BS but the true channels {x 0 t } may be hard to be obtained in advance, which hinders us to find the appropriate initialization based on empirical results. Thus, the SE is possible to help us find an appropriate value of SNR 0 for the initialization of the EM-HyGAMP-DCS algorithm in different systems. Based on the SE, we can try different values of SNR 0 to find the suitable one that assures the EM-HyGAMP-DCS algorithm converges to a satisfied point. Fig. 3 and Fig. 4 show the impact of the initialized value of SNR 0 on the algorithm performance and the SE under the systems with SNR = β  22 The numerical results in Fig. 3 and Fig. 4 show that the performance of EM-HyGAMP-DCS is sensitive to the value of SNR 0 . We can observe that using SNR 0 < SNR can always ensure EM-HyGAMP-DCS algorithm to converge to a satisfied point but may have a slow convergence rate.
For SNR 0 SNR, the EM-HyGAMP-DCS algorithm usually cannot converge to a stationary point and thus achieve a bad performance. Similarly, the SE by using SNR 0 < SNR can also converge to a stationary point. However, if SNR 0 SNR, the SE has fluctuation of TNMSE before it converges in the case of SNR = −10 dB. The SE also have a slower convergence rate when we set SNR 0 SNR. When we have SNR 0 SNR, the SE also cannot converge to a stationary point though it seems to give a smaller TNMSE. For the SNR 0 that is close to the true SNR, the SE is shown to be smooth without TNME fluctuation and have a fast convergence rate. In summary, the SE can be an essential tools for us to select an appropriate value of SNR 0 , which enables EM-HyGAMP-DCS to achieve satisfied performance and fast convergence.
Besides using the SE to find a suitable value of SNR 0 , we can also find the appropriate initialization of the other hyperparameters by the SE. Thus, the SE is not only used to analyze the asymptotic performance of the AMP algorithm under different system settings, but also a powerful tool to find the appropriate hyperparameter initialization for the AMP-based algorithms incorporated with EM algorithm.

B. Computational Complexity Analysis
We compare the computational complexity of the proposed HyGAMP-DCS algorithm with that of the state-of-the-art methods in Table II. The benchmarks include the DCS-OMP algorithm [26], the PIA-ASP algorithm [27], the GAMP algorithm [32], the S-AMP algorithm [28] and DCS-AMP algorithm [6]. The number of the complex multiplications in each iteration is used as the performance metric. Note that DCS-OMP, PIA-ASP and GAMP all contains two levels of iterations as the outer iteration and the inner iteration. Here, the numbers of the complex multiplications in each inner iteration of these three algorithms are given in Table II  In this section, we give the simulation results of our proposed algorithm in the temporallycorrelated massive access system. To show the advantage of our proposed algorithm, the conventional DCS-based algorithm of PIA-ASP [27], the standard GAMP algorithm [32] and the S-AMP algorithm 1 [28] are also evaluated as the benchmarks.
We consider the simulation scenario where there are N = 10 3 users. We define the signalto-noise ratio as SNR = β σ 2 w and the number of consecutive frames for joint detection is set as T = 4. We set the transition probability p 10 = 0.25. The performance metric used for the user activity detection is defined as the time-averaged activity error ratio (TAER), which includes the missed detection ratio and the false alarm ratio and is calculated as where N m (t) is the number of missed detected users in the tth frame and N f (t) is the number of false alarm users in the tth frame. All the numerical results here are obtained by averaging over 10 4 channel realizations.

A. Comparison with the Benchmarks
We evaluate the performance of the proposed HyGAMP-DCS algorithm for user activity detection and channel estimation in Fig. 5 and Fig. 6. One user is detected to be active when the power of its estimated channel is larger than the previously selected threshold , i.e., | x n,t | 2 > , 1 From our simulation trials, we find that the S-AMP algorithm has very similar performance to the SI-aided MMV-AMP algorithm [29]. Thus, we only show the performance of the S-AMP algorithm in this section. increases, which implies that the AMP-based algorithm can acquire larger performance gain than the greedy algorithm in high SNR scenario. By exploiting the temporal correlations, the S-AMP algorithm and the proposed HyGAMP-DCS algorithm can achieve better performance than the GAMP algorithm. Moreover, the HyGAMP-DCS algorithm has further improved performance by additionally exchanging the information between the current frame and the next frame without much computation complexity increase.
The impact of the pilot length L on the detection and estimation performance of the proposed algorithms and the benchmarks is evaluated in Fig. 7 and Fig. 8. We set p a = 0.1 and SNR = −10dB. We observe that the AMP-based algorithms all have larger performance gap to the PIA-ASP algorithm when L increases, which indicates the essential role of the channel statistics for performance improvement. It is also found that the PIA-ASP algorithm can outperform the GAMP algorithm in the user activity detection performance when L is smaller than 190, though the GAMP algorithm always has better channel estimation performance. When L keeps increasing, the GAMP algorithm can outperform the PIA-ASP algorithm and achieve larger performance gain, which means that the temporal correlation provides limited performance gain for the greedy based algorithm. The EM-HyGAMP-DCS algorithm is shown to achieve similar performance to the HyGAMP-DCS algorithm even if the perfect system statistics are unknown.
So that the EM-HyGAMP-DCS algorithm can be more suitable to be applied in the practical T increases, and eventually the performance saturates. As such, we can set a moderate T , which can take a compromise between the performance and the detection delay.  It is observed that strong temporal correlation help all the DCS-based algorithms achieve smaller TAER and TNMSE and their performance is greatly enhanced when p 11 is approach 1. These results demonstrate that temporal correlation plays an important role in the detection and estimation algorithm design. On the other hand, it shows that our proposed algorithm can greatly outperform the benchmarks especially in the scenario with strong temporal correlation, which indicates the advantage of exploiting the temporal correlation from both the neighboring frames.

B. Validation of the Theoretical Analysis
Finally, the simulation result and the theoretical result under the setup of L = 200 and p a = 0.1 are given in Fig. 13 and Fig. 14. We can find that the converged TNMSE of the proposed HyGAMP-DCS algorithm and the EM-HyGAMP-DCS algorithm can be accurately predicted by the SE. Compared with the GAMP algorithm, the proposed algorithms can have faster convergence in the high SNR scenario, which means that the HyGAMP-DCS-based algorithms can achieve performance improvement in both estimation accuracy and computational complexity. The theoretical results also indicate that the EM-HyGAMP-DCS algorithm can achieve almost the same performance as the HyGAMP-DCS algorithm does in the asymptotic regime even if the perfect system statistics are unavailable. Therefore, in the practical system with unperfect system statistics, we can first employ the SE to find the appropriate hyperparameter initialization value and then utilize the EM-HyGAMP-DCS algorithm to nearly achieve the optimal performance.

VIII. CONCLUSIONS
The grant-free temporally-correlated massive access system is studied in this work, where the active users have a large probability to transmit in multiple continous frames. The problem of joint user activity detection and channel estimation in multiple consecutive frames is formulated as a DCS problem with a probabilistic model. We then propose the HyGAMP-DCS algorithm with efficient complexity to exploit both the statistics of channels and the temporally-correlated user activities. Simulation results show that the proposed algorithm can achieve much better performance than the conventional DCS algorithm and the traditional AMP-based algorithms. To acquire the knowledge of the system statistics, the EM algorithm is proposed to be incorporated in the HyGAMP-DCS algorithm by adaptively updating the hyperparameters during the detection and estimation procedure. To ensure that the EM-HyGAMP-DCS algorithm realizes satisfied performance with fast convergence rate, we illustrate that the SE is the powerful tool for finding the appropriate hyperparameter initialization. Additionally, we observe that only a moderate number of continuous frames for joint detection and estimation is enough to approach the optimal performance.
We can observe that all the terms except ln p(y t |x t ) for all t are independent to the noise variance σ 2 w , so that the partial derivative of (51) with respect to σ w is obtained as Finally, we can give the update rule of σ 2 w in (31) by setting the equation (52) equaling to zero.
To tackle the discontinuity of ∂p(xn,t|β,λn,t) ∂β at the point x n,t = 0, we adopt the strategy in [31] where the closed region B {x||x| ≤ , x ∈ C} and its complementaryB = C \ B are defined with the limit → 0. Therefore, we Then the update rule of β can be derived in (32) by setting (55) to be zero.
Finally, we obtain (33) as the update rule for p a by setting (60) equal to zero.
By inserting (63)and (67) into (34), we finally obtain the closed-form expression of the update rule for p 10 .