{'fraction_non_alphanumeric': 0.05730868928209174, 'fraction_numerical': 0.019065329276080994, 'mean_word_length': 3.966787946636813, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accelerated and Deep Expectation Maximization for One-Bit MIMO-OFDM Detection

In this paper we study the expectation maximization (EM) technique for one-bit MIMO-OFDM detection (OMOD). Arising from the recent interest in massive MIMO with one-bit analog-to-digital converters, OMOD is a massive-scale problem. EM is an iterative method that can exploit the OFDM structure to process the problem in a per-iteration efficient fashion. In this study we analyze the convergence rate of EM for a class of approximate maximum-likelihood OMOD formulations, or, in a broader sense, a class of problems involving regression from quantized data. We show how the SNR and channel conditions can have an impact on the convergence rate. We do so by making a connection between the EM and the proximal gradient methods in the context of OMOD. This connection also gives us insight to build new accelerated and/or inexact EM schemes. The accelerated scheme has faster convergence in theory, and the inexact scheme provides us with the flexibility to implement EM more efficiently, with convergence guarantee. Furthermore we develop a deep EM algorithm, wherein we take the structure of our inexact EM algorithm and apply deep unfolding to train an efficient structured deep net. Simulation results show that our accelerated exact/inexact EM algorithms run much faster than their standard EM counterparts, and that the deep EM algorithm gives promising detection and runtime performances.


Introduction
In massive multiple-input multiple-output (MIMO), we have recently seen enormous interest in techniques that are related to low-resolution or coarsely quantized (CQ) signals from a fundamental signal processing viewpoint. This trend is driven by our interest in replacing high-resolution analogto-digital converters (ADCs) and digital-to-analog converters (DACs) with lower resolution ones, so that we can keep the hardware cost and energy consumption of the radio-frequency (RF) front ends at a lower level. In the conventional scenario of high-resolution ADCs/DACs, quantization errors are insignificant and are typically treated as part of the background noise. In the low-resolution scenario, recent research has suggested that we should exploit the low-resolution or quantization structure to better harness the problem. For example, for MIMO downlink precoding with one-bit DACs, a popular direction is to directly design the one-bit MIMO transmit signals by discrete optimization; see, e.g., [1][2][3][4] and the references therein. In addition to exploiting the low-resolution structure, there have been endeavors that apply Sigma-Delta modulation-a widely-used technique in temporal DACs/ADCs-to spatially shape the quantization errors such that users in a designated angular sector will experience much less quantization effects [5,6].
Another meaningful problem in massive MIMO with low-resolution ADCs/DACs is MIMO detection. The problem typically appears in the uplink, and the goal is to detect signals transmitted by different users over the same RF band. While numerous MIMO detection methods were already developed in the conventional high-resolution ADC scenario, they do not assume the presence of quantization and may not be directly applied to the low-resolution ADC scenario. We have indeed seen research on applying linear receivers to the low-resolution scenario [7][8][9], but recent research appears to favor exploiting the quantization structure. The key idea that underlies the majority of low-resolution or CQ MIMO detection methods is that we consider maximum likelihood (ML) or related inference approaches, and the likelihood function is derived from quantized signals (quantization is treated as a function, and we do not model quantization errors as noise). This rationale gave rise to a variety of CQ MIMO detection methods; e.g., expectation maximization (EM) [10], convex and non-convex proximal gradient (PG) methods [11][12][13], deep learning [13][14][15], and posterior inference via approximate message passing [16] and variational Bayes [17].

The Problem of Interest: CQ MIMO-OFDM Detection
In MIMO detection it is common to assume frequency-flat channels. Yet, in the current mobile systems, we consider frequency-selective channels. This is not a problem in the high-resolution scenario, because the current systems usually use the orthogonal frequency division multiplexing (OFDM) scheme, and OFDM breaks the frequency-selective channels into many parallel frequencyflat channels. This means that the MIMO detection methods for frequency-flat channels can be directly applied. Unfortunately, the same does not hold for the CQ case. The presence of quantization destroys the structure required to turn the problem into parallel frequency-flat channels. Subsequently, when we consider MIMO-OFDM detection with low-resolution ADCs, we have to engage with the full form of the problem. The difficulty arising is that the problem is a massive-scale problem; the problem size scales with the OFDM size, which is of the order of hundreds. While in principle we can directly apply the CQ MIMO detection methods in the frequency-flat channel case to the MIMO-OFDM case, it is computationally infeasible to do so due to the massive problem size. It is necessary to exploit the OFDM structure to efficiently realize MIMO detection methods.
There are far fewer studies for CQ MIMO-OFDM detection [10,[17][18][19]. In [18], the PG method for the box ML relaxation formulation was developed. The PG method is able to exploit the OFDM structure to have per-iteration efficient implementation. In [10,19], the EM method for the Gaussian maximum a posteriori formulation was built. The EM method is iterative, and at each iteration we solve a problem-the so-called M-step problem-that takes the same form as MIMO-OFDM without quantization. Hence we can exploit the OFDM structure to have per-iteration efficient implementation, and this makes the EM method appealing in OMOD. The same benefit was also seen in the application of the variational Bayes method (which is related to EM) to CQ MIMO-OFDM detection [17].

This Work
In this paper we study the EM technique for CQ MIMO-OFDM detection. In particular, we consider MIMO-OFDM detection with one-bit ADCs, or simply one-bit MIMO-OFDM detection (OMOD). We are interested in EM convergence rate analysis for a class of approximate ML OMOD formulations. This aspect has not been seriously studied in CQ MIMO detection. Our analysis draws insights from the analysis techniques in PG methods [20][21][22][23]. We will reveal a connection between the PG and EM methods, which is unobvious at first sight, and from that connection we establish EM convergence results by taking key ideas from the PG analysis techniques. We will show that (a) EM converges at a rate of at least 1/k, where k denotes the iteration number; that (b) EM should converge faster than PG, in theory. Also, we show how the SNR and channel conditions can have an impact on the convergence rate. In particular, our analysis indicates that the convergence rate is slower as the SNR is higher; this applies to both PG and EM. This high-SNR slow convergence phenomenon was noticed by numerical results in previous studies [13,18], and, to the best of our knowledge, this is the first time we see a theoretical result for such phenomenon.
Our EM convergence study also leads to new OMOD algorithms. In first-order optimization, it is well-known that the PG method can be accelerated to a convergence rate of at least 1/k 2 if we employ extrapolation [20,21]; we assume convex problems. Using the connection between the PG and EM methods, we apply the acceleration idea to EM. The result is an accelerated EM scheme for OMOD, which has the 1/k 2 convergence rate guarantee and runs faster than the standard EM as will be shown in our numerical results. In addition we consider an inexact EM scheme, where we want to solve the M-step problems not too precisely to save complexity. By taking insight from the study of inexact PG [22], we provide theoretical conditions under which the inexact EM scheme will have convergence guarantees. The inexact EM scheme can be used together with acceleration, and our numerical results will demonstrate that the accelerated inexact EM algorithm for OMOD has significant merits in runtime savings.
Lately we have seen much interest in applying deep learning to MIMO detection [13,15,[24][25][26][27]. The most successful approach at present is deep unfolding [28]. The principle is to see a modelbased algorithm as a deep net, take the structure, and use data-driven learning to try to learn a better deep net. We apply deep unfolding to our inexact EM algorithm for OMOD; that is, simply speaking, we take the structure of how efficient EM can handle MIMO-OFDM, and we use that to build an efficient highly-structured deep net. Our numerical results will show that the deep inexact EM algorithm is competitive in both detection and runtime performances, compared to the previous OMOD algorithms.
It is worth noting that the convergence analysis and the accelerated and inexact EM schemes in this study are, in a broad sense, also applicable to a variety of formulations arising from regression from CQ data, such as CQ MIMO channel estimation [29,30] and CQ compressive sensing [31].

Related Works
We should discuss related works. Our EM convergence analysis plays a large part in this study, and one may question if the existing convergence analyses [32][33][34][35] provide the same or better results. The existing studies mostly consider a general framework, covering a variety of optimization schemes and a rich class of problems. Our study, on the other hand, exploits the specific structure of our OMOD problem; e.g., the previously mentioned result of high-SNR slow convergence is unique to our problem. In terms of the analysis techniques used, our analysis is considered closest to Mairal's study [34], which considers a general majorization minimization (MM) framework and can cover EM. Mairal's study takes insight from the PG analysis techniques through some assumptions. We also draw insight from the PG analysis techniques, but we do so by best using our problem's specific structure. As we will discuss in details, our result arguably provides better characterization of the convergence behaviors.
One key result in this study is accelerated EM for OMOD. Indeed, Mairal's study [34] has described the same kind of acceleration for MM. To the best of our knowledge, no application was provided in Mairal's study and in subsequent studies. Our application is probably the first, showing the efficiency of such accelerated EM or MM. We should also discuss the application of deep unfolding. While deep unfolding has been widely considered in MIMO detection, it has not been used in CQ MIMO-OFDM detection.
The preliminary version of this study was presented in a conference [36]. It examines inexact EM as a practical algorithm for OMOD. It neither provides convergence analysis nor develops the accelerated and deep EM algorithms.

Organization and Notations
This paper is organized as follows. Section 2, 3 and 4 reviews the EM basics, the OMOD problem and the state of the arts of OMOD, respectively. In Section 5, we perform EM convergence analysis and establish the accelerated and inexact EM methods. Section 6 develops an accelerated inexact EM algorithms for OMOD. Section 7 build a deep inexact EM algorithm via deep unfolding. Section 8 provides the numerical results, and we conclude in Section 9.
Some basic notations are described as follows. The sets of real, non-negative, non-positive and complex numbers are denoted by R, R + , R − and C, respectively; we denote j = √ −1; boldface lowercase letters, e.g., x, represent column vectors; the ith element of a vector x is denoted by x i ; boldface capital letters, e.g., X, represent matrices; the superscripts ⊤ , H , −1 , and † denote the transpose, Hermitian transpose, inverse, and pseudoinverse, respectively; ℜ(x) and ℑ(x) denote the real and imaginary parts of x, respectively; x = (x 1 , . . . , x m ) means that x = [ x ⊤ 1 , . . . , x ⊤ m ] ⊤ ; · denotes the Euclidean norm; 0 and I denote an all-zero vector and an identity matrix, respectively; Diag(x) denotes a diagonal matrix whose (i, i)th element is x i ; x ∼ D means that x is a random vector following a probability distribution D; N (µ, Σ) denotes a multivariate Gaussian distribution with mean µ and covariance Σ; N (x; µ, Σ) denotes the distribution function corresponding to N (µ, Σ).

EM for Linear Regression from CQ Data
In this section, we review the expectation maximization (EM) method [37] for the problem of linear regression from truncated data [38]. Consider the following data model where y ∈ {±1} m is the observed data; sgn : R m → {±1} m is the sign function; A ∈ R m×n is known; θ ∈ R n is a model parameter; v ∼ N (0, σ 2 I) is noise; σ 2 is known. Note that x is unobserved. The problem is to estimate θ from y, and we consider maximum-likelihood (ML) estimationθ ∈ arg max θ∈R n L(θ) := log p(y|θ), where p(y|θ) is the distribution of y given θ; the same notation will be used to denote the distributions of other random variables. It can be shown that p(y|θ) = m i=1 p(y i |θ), where a i denotes the ith row of A; X (y i ) = R + if y = 1, and X (y i ) = R − if y = −1. This leads to where Φ(z) = z −∞ e −t 2 /2 dt/ √ 2π is the standard cumulative Gaussian distribution. It is known that L is concave [11,39].
The ML estimation problem (2) has no closed-form solution in general, and EM is a popular way to solve the ML problem. Intuitively, the idea is to leverage on the fact that estimating θ is easy if x is observable. To describe the EM method, let .
Here, E x [f (x)|y, θ ′ ] = X f (x)p(x|y, θ ′ )dx denotes a conditional expectation of x; X denotes the support of p(x|y, θ ′ ). It is known that G(θ|θ ′ ) ≤ L(θ) for all θ, θ ′ , and G(θ ′ |θ ′ ) = L(θ ′ ) for all θ ′ ; i.e., G(θ|θ ′ ) is a lower bound of L(θ). EM considers a successive lower-bound approximation An important question is whether Q(θ|θ ′ ) is tractable. By noting p(x|θ) = N (x; Aθ, σ 2 I), one can verify that where "const." does not depend on θ. It can be shown that where φ(t) = e −t 2 /2 / √ 2π; see, e.g., [29], for details. To prove (6), one needs to show that p(x i |y i , θ) is a truncated Gaussian distribution, namely, p(x i |y i , θ) = p(x i |θ)/p(y i |θ), with X (y i ) being the support; that E x i [x i |y, θ ′ ] is the mean of p(x i |y i , θ ′ ); that a truncated Gaussian mean has an explicit expression [29]. The EM method can hence be written as where x k = E x [x|y, θ k ] is computed by (6). The result is an iterative linear regression: we estimate x by the conditional mean x k , inferred from y and the previous estimate θ k ; and then we estimate θ by the linear regression in (7), which is easy to solve. Intuitively, this alternating process of estimating x and θ makes sense. It is worth noting that the EM method for linear regression from truncated data dates back to as early as 1980 in statistics, with application to blood lead data [38]. It also arose in signal processing, with applications to MIMO channel estimation and detection from coarsely quantized signals [10,19,29,30].

One-Bit MIMO-OFDM Detection Background
In this section, we review the background of one-bit MIMO-OFDM detection.

MIMO-OFDM Model
We start with the signal model. Consider a wireless uplink scenario where a number of users simultaneously send signals to a multi-antenna base station (BS). The transmitted signals undergo frequency selective channel effects, upon reception at the BS; and the user side employs the OFDM transmission scheme. The received signals at the BS can be modeled as where r m ∈ C W is a block of discrete-time signals received by the BS's mth antenna; W is the block length, or the OFDM size; M is the number of antennas at the BS; F ∈ C W ×W is the unitary discrete Fourier transform (DFT) matrix; s n ∈ C W is a block of symbols sent by user n; N is the number of users; H m,n ∈ C W ×W is a circulant matrix, constructed from a discrete-time impulse response h m,n = (h 0,m,n , . . . , h W ′ −1,m,n , 0, . . . , 0) ∈ C W of the channel from user n to the BS's mth antenna, with W ′ ≪ W ; ν m is circular complex Gaussian noise with mean 0 and covariance σ 2 C I, and with ν 1 , . . . , ν M being independent. The symbols s w,n 's are drawn from a QAM constellation set where D is a positive integer. The above model is standard, and we refer the reader to the literature (e.g., [40]) for details. For conciseness, we let r = (r 1 , . . . , r M ) and express (8) as where s = (s 1 , . . . , s N ), ν = (ν 1 , . . . , ν M ), The above model is based on an assumption that the signal acquisitions at the BS have negligible quantization errors; to put it another way, the BS is equipped with high-resolution ADCs. In this study we consider the one-bit ADC case wherein q m = sgn(ℜ(r m )) + j · sgn(ℑ(r m )), m = 1, . . . , M, are the acquired signals. As mentioned earlier, this one-bit case is motivated by the massive MIMO scenario where we want to use cheap ADCs to reduce the hardware cost and energy consumption of the RF front ends.

Subcarrier Decoupling in Unquantized MIMO-OFDM
In this subsection we review the MIMO-OFDM subcarrier decoupling concept, which we will need later. To describe it, consider the ML detection of s from the unquantized signal r s ∈ arg min s∈S NW r − Hs 2 .
At first look, this is a massive-scale problem; the OFDM size W is often of the order of hundreds. But it is widely known that (12) can be decoupled into many smaller problems. Since H m,n is circulant, we have the eigendecomposition whereh m,n = F h m,n . Letr m = F r m . One can show that whereř w = (r 1,w , . . . ,r M,w ) ∈ C M ,Ȟ w = [h w,m,n ] m,n ∈ C M ×N andš w = (s w,1 , . . . , s w,N ) ∈ C N are the received signal, MIMO channel response and symbol vector at subcarrier w, respectively. We can hence rewrite (12) as We see that each problem in (14) has a much smaller size, and can be dealt with independently.

One-Bit MIMO-OFDM Detection (OMOD)
Let us consider one-bit MIMO-OFDM detection, or OMOD for short. From (9) and (11), we can write where q = (q 1 , . . . , q M ), v ∼ N (0, σ 2 I), σ 2 = σ 2 C /2. The above model takes the same form as the one in Section 2. Assume that H and σ 2 are known. Following the formulation in Section 2, the ML detection for OMOD is formulated as where b i = y i a i /σ; a i is the ith row of A; or, There is an important aspect to note. The subcarrier decoupling in Section 3.2, which allows us to decouple the massive-scale unquantized MIMO-OFDM problem into many smaller problems, does not work for OMOD. This is because f in (17) does not possess the structure required to apply subcarrier decoupling. This means that we will have to confront the massive-scale nature of MIMO-OFDM in the one-bit case.

Problem Statement, and State of the Arts
In this section we provide the main problem statement and review the existing methods for OMOD.

Problem Formulation
Our study will revolve around the following formulation for some penalty function h : The above formulation is an approximation of the ML OMOD detection in (17). Consider the following two examples.

2) Box [18]
Consider where This corresponds to a relaxation of the ML OMOD detector (17) by replacing the constellation set D with an interval [−U, U ].
The main issue is how we can handle problem (19), which has a massive scale, in an efficient fashion. The next two subsections describe two existing methods to deal with (19).

Proximal Gradient
One way to handle problem (19) is the proximal gradient (PG) method [23]. The PG method for (19) is given by where η k > 0 is a step size; ∇f (θ) is the gradient of f at θ, and it can be shown that is the proximal operator associated with h. We typically consider cases for which prox h (x) is easy to compute. For example, for the box We want to examine the complexity of the PG method. The most expensive part lies in computing the gradients ∇f (θ k ). At first sight, computing ∇f (θ) in (24) seems expensive since B is very large. But one can use the OFDM structure to reduce the complexity. It was shown in [18] that ∇f (θ) can be computed with a complexity of where C denotes the complexity to numerically compute Φ, given a precision 1 . For the reader's convenience, in the supplemental material we describe the details with the computation of ∇f (θ). Simply speaking, to compute ∇f (θ) efficiently, we need to run fast Fourier transforms (FFTs) or inverse FFTs (IFFTs) for a total of 2M times. To sum up, the per-iteration complexity of the PG method is (26). The PG method for OMOD was studied in [18], wherein the efficient FFT/IFFT-based method for computing the gradients was proposed, and the box relaxation was considered.

Expectation Maximization
Another way to handle problem (19) is the EM method. Following the EM review in Section 2, let be the EM surrogate function of f in problem (19), where r ′ = E r [r|q, s ′ ] is the conditional mean estimate of the unquantized received signal r. The EM method is given by where, for convenience, we denote h(s) = h(ℜ(s)) + h(ℑ(s)). The merit of the EM method is that it allows us to apply subcarrier decoupling. To describe this, let By the subcarrier decoupling concept in Section 3.2, we can decompose (28) intǒ for w = 1, . . . , W . The basic description of the EM method is complete. Intuitively, the EM method seems more promising than the PG method in the previous subsection. The EM method has no parameter, such as step size, to select; it inherits the merits of subcarrier decoupling, which seems to suggest better efficiency; it performs more optimization in one iteration, specifically, (30), which seems to suggest faster convergence.
The story is actually more complex, as one tries to unravel. We should not forget the fact that the conditional means r k m 's need to be computed. Examining the equations closely (see (6), (8), (9), and (15)), one will find that r k m is given by for all m, n, where ⊙ denotes the element-wise product. It can be verified that the complexity of (29) and (31) is particularly, we need FFTs and IFFTs for a total of 2M times. The complexity in (32) is the per-iteration complexity of EM, without counting the complexity required to solve (30). We see that (32) is not much better than the per-iteration complexity (26) for the PG method. When the complexity of solving (30) factors in, the per-iteration complexity of the EM method is likely to be higher than that of PG. The next question then goes to whether EM leads to faster convergence. We will try to answer this question analytically in Section 5, and empirically in Section 8. The EM method for OMOD was studied in [10], wherein the efficient subcarrier-decoupling implementation was proposed, and GMAP was considered. For GMAP, the problems in (30) have closed-form solutions, and this advantage was exploited.

Convergence Analysis and New Insights
In this section, we analyze the theoretical convergence of the PG and EM methods in the context of OMOD. As will be shown, the analysis leads to new insights and methods.

Problem Setup and Preliminary Concepts
Let us write down our main problem, problem (19), again: where proper, lower semicontinuous and convex. We assume that problem (33) has an optimal solution, and we will use the notation θ ⋆ to denote an optimal solution to problem (33). We also recall some key formulas. The gradient ∇f (θ) can be written as The EM surrogate function of f can be written as Some concepts and notations should be introduced. Let ϕ : Consider a minimization problem min θ∈R n ϕ(θ). A pointθ ∈ R n is said to be a critical point of the problem if 0 ∈ ∂ϕ(θ), where ∂ϕ(θ) is the limiting subdifferential of ϕ at θ; see [41] and the references therein. A pointθ is said to be an ε-critical point of the problem if dist(0, ∂ϕ(θ)) ≤ ε, where dist(y, X ) = min x∈X y − x denotes the distance between y and X . If ϕ is convex, then any critical point of the problem is an optimal solution to the problem. The notations σ max (X) and σ min (X) denote the largest and smallest singular values of X, respectively; ·, · denotes the inner product.
Due to space limitation, all the new theoretical results will have their proofs relegated to the supplemental materials.

Proximal Gradient
Consider the PG method for problem (33): where we assume a constant step size η > 0. The following result is well known.
Fact 1 (see, e.g., [23]) Consider problem (33) and the accompanied problem setup, and consider the PG method (38). Suppose that f has L f -Lipschitz continuous gradient on dom F . If the step size is chosen as η = 1/L f , then we have Fact 1 suggests that the PG iterate θ k approaches the optimal solution, in terms of the objective value, at a rate of at least 1/k. Our question is whether f possesses the desired Lipschitz continuous gradient property. This appears to have not been answered before, and we show that the answer is yes.
Consequently, the PG method (38) with step size η = 1/σ max (B) 2 has the following convergence result for problem (33): It is worth noting that, from (18), we can show that We can use (40) to select the step size η; note that the prior work [18] selects the step size manually.

Expectation Maximization
Next, consider the EM method for problem (33): Let us first show our EM convergence result.
Proposition 2 Consider problem (33) and the accompanied problem setup, and consider the EM method (41). We have We will describe the insight behind Proposition 2 in the next next subsection. Eq. (42) shows that the EM method has a convergence rate of at least 1/k-same as PG. It is interesting to compare the convergence results of PG and EM. Since the EM convergence result (42) is seen to be at least no worse than the PG result in (39). In fact, the intuitive impression is that the EM method should converge faster than the PG method, particularly when the gap in (43) is large.
One may wonder if, by applying the existing convergence analysis results [32][33][34][35] to our problem, we can obtain the same or better claim than (42). This aspect is complex, and we provide a detailed account of it in the supplemental material. In summary, the best result we can find is for some constant R > 0; it comes from Mairal's analysis [34]. Our result in (42) appears to be better than (44); e.g., it is unconvincing to use (44) to argue that EM should converge faster than PG.

A Slow Convergence Phenomenon
Before we show the insight behind Proposition 2, we want to examine another aspect. The PG and EM theoretical bounds in (39) and (42) are seen to be loose if B has a large magnitude. From (18), the magnitude of B is seen to increase with the SNR; see also (40). This gives an impression that, under high SNRs, the PG and EM methods could converge slowly. Further analysis shows that this is true for our interested choices of h.
Proposition 3 Consider problem (33) and the accompanied problem setup. Consider either the PG method (38) with η = 1/σ max (B) 2 , or the EM method (41). Suppose that θ k is ε-optimal in the sense that for some optimal solution θ ⋆ to problem (33) and for ε > 0.
(a) Suppose that h(θ) = Á X (θ) for some non-empty convex closed set X . Then k must satisfy for the PG method; and for the EM method, where σ + min (B) denotes the smallest positive singular value of B, and P = B ⊤ (B ⊤ B) † B is an orthogonal projection matrix, which performs projection onto the subspace spanned by b 1 , . . . , b m .
for the PG method; and for the EM method.
Proposition 3 shows that the number of iterations required by the PG or EM method must increase with the magnitude of B, or the SNR, assuming that the other variables are held fixed.

Insight Behind The Convergence Results
In this subsection we reveal the insight behind the EM convergence result in Proposition 2. The story begins with how PG convergence results were established in the literature. The idea is to see the PG method as a majorization-minimization (MM) method. Specifically, consider the MM method The PG method (38) with step size η = 1/L f can be regarded as an MM, with Many PG convergence proofs make strong use of the structure of (46). As an important discovery, we found that the EM surrogate function g in (36) possesses the following structure We see that (47) resembles (46). This has tremendous insights. It suggests that we may use the PG convergence proof concepts in the literature to show convergence for EM. In fact, the proof of Proposition 2 is conceptually the same as that of Fact 1.
We also want to show the insight behind Proposition 3. As an example, consider h(θ) = 0 and the PG method (38) with step size η = 1/σ max (B) 2 . From (34) one can show that where the second inequality is due to 0 < ψ(z) < 1 (see (35)). Applying (48) to the PG method (38) gives We see that, if the magnitude of B is large, then every θ k changes little, and we have slow convergence. This is the idea that underlies the results in Proposition 3.

Accelerated EM
The PG-EM resemblance revealed in the previous subsection opens new possibilities. In first-order convex optimization, it is well known that the convergence of the PG method can be accelerated by applying extrapolation [20,21]. We can accelerate the EM method in the same manner. Consider the following modified EM method for k ≥ 0, where ϑ 0 = θ 0 ; {α k } k≥1 is an extrapolation coefficient sequence. We employ the FISTA sequence where t 0 = 1, This modified EM method is the direct consequence of the EM-PG resemblance and accelerated PG concepts. In particular, if we replace g(θ|ϑ k ) in (50) with u(θ|ϑ k ) in (46), we return to the widely-known accelerated PG method. Same as the accelerated PG method, the modified EM method (50)-(53) has an accelerated convergence rate of at least 1/k 2 .
Proposition 4 Consider problem (33) and the accompanied problem setup, and consider the modified EM method (50)-(53). We have

Inexact EM
The EM method (41) requires us to solve the problems in (41) exactly. When the exact solutions are computationally non-negligible to compute, we may want to accept economical, but inexact, solutions. This motivates the study of inexact EM for k ≥ 0, where θ k+1 is an approximate solution, produced by a solver that may not give high solution precision due to its limited complexity. The question is how robust the EM method is under solution errors. Intriguingly, the same kind of questions has been addressed in the context of inexact PG [22]. We take insight from the aforementioned context to answer the question at hand. We characterize the inexactness of θ k+1 by dist(0, ∂G(θ k+1 |θ k )) ≤ ε k+1 , for some ε k+1 ≥ 0; ε k+1 = 0 means that θ k+1 is an optimal solution to the problem in (55). Also, we assume the following descent condition Eq. (57) makes sense in the following way: If the inexact solver for the problem in (55) is warmstarted by taking θ k as the initialization, the resulting approximate solution θ k+1 should yield an objective value G(θ k+1 |θ k ) better than the initial objective value G(θ k |θ k ). We have the result below.
Proposition 5 Consider problem (33) and the accompanied problem setup. Consider the inexact EM method (55), with (56)-(57) being satisfied. Suppose that B has full column rank. For k ≥ 1, we have Proposition 5 reveals two aspects. First, if the solution accuracy ε k improves with k in such a way that ∞ i=1 ε i < ∞, i.e., ε k decreases at a rate of k −p for p > 1, or faster, then the inexact EM method has a convergence rate of at least 1/k. Second, the inexact EM method may be more sensitive to solution errors if the smallest singular value σ min (B) of B is smaller, i.e., B is more ill-conditioned. In massive MIMO, we should expect that B is reasonably well-conditioned. We can also consider the accelerated version of inexact EM. (33) and the accompanied problem setup. Consider the accelerated version of the inexact EM method (55), wherein θ k is replaced by ϑ k in (51)-(53). Suppose that (56) holds, and that B has full column rank. For k ≥ 1, we have

Proposition 6 Consider problem
Note that the above result does not require the descent condition (57). Proposition 6 looks similar to its non-accelerated counterpart in Proposition 5. The notable difference is that, to achieve a convergence rate of at least 1/k 2 in Proposition 6, we need the solution accuracies ε k 's to satisfy ∞ i=1 i ε i < ∞; i.e., ε k decreases at a rate of k −p , p > 2, or faster. This means that the accelerated inexact EM method may require higher solution accuracies than the non-accelerated counterpart.

The Case of Non-Convex h
Our preceding analyses assume convex h. We can also handle non-convex h. Consider the same inexact EM scenario in (55). The new difficulty is that the problems in (55) are non-convex. We assume that the inexact solver for the problems in (55) can only guarantee finding an ε-critical point, rather than an optimal or ε-optimal solution. The ε-criticality of θ k+1 is characterized by (56). Also, we assume that the descent condition (57) holds; our justification is the same as before. Furthermore, we modify G(θ|θ k ) in (55) as for some τ > 0. We have the following result.
Eq. (59) employs an ellipsoidal extension of the ε-critical point definition to pin down the convergence. As an easier way to understand, consider the following expression which can be shown to be an implication of (59): Eq. (60) suggests that the EM method leads to a δ-critical

Accelerated EM Algorithms for OMOD
In this section we use the new theoretical insights gained in the previous section to build efficient algorithms for OMOD. Our attention is paid to the implementation of the accelerated inexact EM (AIEM) method.

An AIEM Algorithm for a General Convex h
Consider the OMOD problem (19) with a general convex h. By modifying EM OMOD method in Section 4.3 as the AIEM method in Section 5.7, we have the following AIEM OMOD method: Let s 0 be the initialization, and let s 0 ex = s 0 . We perform, for k = 0, 1, 2, . . ., Here {š k+1 w } W w=1 satisfies a solution quality condition for some given ε k+1 > 0, where e k+1 ∈ R W is given by 4. Extrapolation for acceleration: where {α k } k≥1 is given by (52)-(53).
Note that if we solve the problems in (61) exactly, the above method becomes the accelerated (exact) EM; that if we set α k = 0 for all k, the above method becomes the regular non-accelerated EM; and that (62) is identical to the solution quality requirement in (56). Proposition 6 suggests that we can choose ε k = C k −p for some C > 0 and p > 2.
We have not specified the inexact solver for the M-step problems in (61). The AIEM method provides us with the freedom to choose the solver, and in this study we employ the accelerated PG (APG) method. Fixing k and w, the APG method for the problem in (61) is given by for j ≥ 0, where η w = 1/σ max (Ȟ w ) 2 is the step size, and ∇ϕ k w (s) =Ȟ H wȞw s −Ȟ H wř k w . We stop the APG loop when {s j+1 w } W w=1 meets the solution quality requirement (62). Assembling the above components together, we obtain the AIEM algorithm in Algorithm 1.
Let us also consider the GMAP case h(θ) = λ θ 2 /2. The problems in (61) can actually be solved in closed form [10]:š k+1 w Hence we can replace the APG solver (Line 10-17 of Algorithm 1) with (66), and the result is an accelerated EM.

A Deep Inexact EM Algorithm for OMOD
In this section we apply deep unfolding [24,28] to the inexact EM algorithm to build a data-driven detector for OMOD. The idea of deep unfolding is to alter an existing iterative algorithm by datadriven learning, thereby hoping to learn a better algorithm. Specifically, we see each iteration of compute ζ k m for all m via set r k m = z k m + ζ k m for all m if a stopping rule is met, exit and output s k+1 21: end for the predecessor algorithm as a network layer, untie some of the parameters or alter some of the structures, and learn the untied parameters or new structures from data. We tried various ways to deeply unfold Algorithm 1. The successful one, by our empirical tests, is shown in Algorithm 2; we will call the algorithm DIEM. The untied parameters are marked red, and some notable changes marked blue.
The most notable alternation of DIEM is Line 9 of Algorithm 2. We replace the proximal operator prox ηwσ 2 h (Line 12 of AIEM) with a nonlinear activation function Ω γ . We use a multilevel sigmoid function; its real-valued scalar version is where γ > 0, ̺ : R → [−1, 1] is a 0-centered sigmoid function: we have Ω γ (s) = [ℜ(Ω γ (s i ))+j ℑ(Ω γ (s i )] i for the complex-valued vector case. The sigmoid functions compute ζ k m for all m via set r k m = z k m + ζ k m for all m setř k w = (r k 1,w , . . . ,r k M,w ) for all w 8: set s k+1 such thatš k+1 w = s k+1 w for all w 12: end for were used in deep MIMO detection [25] as the activation function. The intuitive idea is that Ω γ approaches the decision function of the constellation set as γ → ∞, and we want to use it to encourage the symbol estimates to lie closer to the constellation points.
The other alterations of DIEM are as follows: we abandon acceleration in the EM loop; we consider one-step APG with the inexact solver; we untie the step size η k and extrapolation coefficient α k of the APG in a layer-varying fashion; we add a parameter β k to partially alter the noise standard deviation.
The training of DIEM is standard. Let p = {s, H, σ, q} be a problem instance, let π = {α, β, γ, η} be the set of network parameters, and let o π (p) be the network output. We generate a large number of problem instances p 1 , . . . , p T according to the signal model. Then we learn π from p 1 , . . . , p T using a deep learning software; the objective is where s t is the symbol vector of the problem instance p t .
We want to provide an interpretation for DIEM. Consider a case where we choose a nonconvex penalty h with the OMOD (19) to better approximate the difficult constellation constraints. As studied in Section 5.8, the corresponding nonconvex OMOD problem (19) can be handled by inexact EM. Following Section 6, we can show that the nonconvex inexact EM can be implemented by nearly the same procedure as the non-accelerated version of Algorithm 1. From this view, DIEM can be regarded as the deep unfolding of the nonconvex inexact EM-the constellation-promoting h is learnt bottom-up. It is worth noting that a similar viewpoint was used to develop deep learning methods for inverse problems in imaging; see [42] and the references therein.

Numerical Results
In this section numerical results are provided to illustrate the detection and runtime performance of the new algorithms.

Simulation Settings
We use Monte Carlo simulations to test the performance of the algorithms. The signals are generated by the model in Section 3.1. The channel is generated by a multipath model; see, e.g., [43] and the reference therein. To describe the model, leth l,n = (h l,1,n , . . . , h l,M,n ), l = 0, . . . , W ′ − 1, be the lth coefficient of the time-domain channel impulse response from user k to all the BS antennas. Assuming that the BS antennas are arranged as a uniform linear array, we model eachh l,n as h l,n = J j=1 α j,l,n a(θ j,l,n ), where a(θ) = (1, e −j 2πd λ sin(θ) , . . . , e −j(M −1) 2πd λ sin(θ) ) is the response of the array to a signal coming from an angle θ ∈ (π/2, π/2); d is the inter-antenna spacing; λ is the wavelength; J is the number of paths; α j,l,n and θ j,l,n are the complex channel gain and angle of a path, respectively. In our simulations, we set W ′ = 16, J = 4 and d = λ/2; each α j,l,n is randomly generated, following a circular complex Gaussian distribution with mean zero and variance 1/J; each θ j,l,n is randomly generated, following a uniform distribution on (π/2, π/2). The number of trials of our simulations is 500. In addition, the SNR is defined as . Our simulations consider the following algorithms.
1. GMAP-EM: the standard EM method for GMAP OMOD, proposed previously in [10] and reviewed in Section 4.3; 2. GMAP-AEM: the accelerated EM method for GMAP OMOD, built in Section 6; 3. BOX-PG: the PG method for box OMOD, proposed previously in [18] and reviewed in Sec- We tested all the algorithms by MATLAB 8.5 on a desktop with Intel i7-7700 processor and 16GB RAM memory. The parameter settings with the PG and EM algorithms (the non-deep ones) are as follows. The stopping condition is either s k+1 − s k / s k ≤ 5 × 10 −4 or k > 1000. For BOX-AIEM, we set the M-step solution accuracies as ε k = 2N W k −2.1 . For BOX-EM, we set ε k = 2N W × 10 −4 . For BOX-PG, we choose the step size η as the reciprocal of (40); see Proposition 1. We set the noise power parameter σ 2 as σ = σ actual + σ 0 , where σ 2 actual is the actual noise power; σ 0 > 0 is some given constant. Doing so is to avoid slow convergence. As shown in Proposition 3, the PG and EM methods can converge slowly for small σ 2 . In our simulations, we set σ 0 = 3. In the supplemental material, we provide additional numerical results to demonstrate the impacts of σ 0 on the performance.
The settings with DIEM are as follows. The training is implemented on Pytorch 1.2 platform for Python 3.5.2. The training optimizer is stochastic gradient descent, with step size 10 −3 . The SNR range in our training is the SNR range in our simulations, with an additional margin of 2 to 5dB. The number of layers is K = 20. We will release the source codes of DIEM at a later stage.

Comparison of EM and Accelerated Exact/Inexact EM
In this subsection we focus on comparison of the standard EM method and our modified EM methods. Fig. 1 shows the bit error rate (BER) performance of the GMAP-and box-based algorithms. GMAP-AEM is seen to yield nearly the same performance as GMAP-EM. This is expected since both solve the same problem, and the problem is convex. Similarly, we see the same behaviors with the box-based algorithms. Fig. 2 shows the average numbers of EM or PG iterations used by the algorithms. As can be seen, the accelerated EM methods use much less iterations than the PG and standard EM methods. This demonstrates the benefit of accelerated EM.
From Fig. 2, we also observe a phenomenon, namely, that the numbers of iterations of the PG and EM algorithms rise with the SNR. This agrees with the theoretical results in Proposition 3. Fig. 3 shows the average runtimes. It is more interesting to examine the box case. While we see in Fig. 2(b) that BOX-EM uses less numbers of iterations than BOX-PG, we see in Fig. 3(b) that BOX-EM runs slower than BOX-PG. This is because the M-step in EM requires non-negligible computations. Moreover, BOX-AIEM is seen to be faster than BOX-EM (and also BOX-PG); the use of inexact M-step is a key factor for its runtime advantage.

Performance Comparison for DIEM
In this subsection we demonstrate the performance of DIEM. Figs. 4-5 show the BER performance of DIEM under various system settings. We benchmark DIEM against ZF, GMAP-AEM and BOX-AIEM. DIEM is seen to exhibit remarkably better BER performance than the other algorithms. Tables 1-2 show the runtimes of the various algorithms. As can be seen, DIEM is much faster than the other algorithms. The reason for the efficiency of DIEM is that it uses a fixed number of layers, specifically, 20; and each layer has simple operations, specifically, one-step APG-like operations.

Conclusion
To conclude, the EM convergence rate was analyzed for a class of approximate ML formulations for OMOD. The analysis revealed the possibility of slow convergence under high SNRs. It also led to accelerated and inexact EM schemes. The implementation of these schemes for OMOD was    studied, and the use of deep unfolding in EM for OMOD was also considered. The resulting EM algorithms were numerically shown to have competitive performance in the OMOD application. 1 Computation of the Gradient in (24) We describe how the gradient ∇f (θ) in (24) can be efficiently computed. Recall θ = (ℜ(s), ℑ(s)), and let f (s) = f (θ) for convenience. Let ∇ sn f (s) = ∇ ℜ(sn) f (s) + j · ∇ ℑ(sn) f (s) denote the gradient of f w.r.t. s n , where ∇ ℜ(sn) f (s) and ∇ ℑ(sn) f (s) denote the gradient w.r.t. the real and imaginary parts of s n , respectively. Working through the details (see (24), (8), (9), and (15)), one can show that for all m, n, where ψ is defined in (35). Eqs. (69a) and (69d) are more expensive to compute, at first sight, but one can utilize the OFDM structure to compute them efficiently. Using (13), we can rewrite (69a) and (69d) as respectively. We can use IFFT and FFT to implement (70)-(71). In summary, the gradient ∇f (θ) in (24) can be computed using 2M FFTs/IFFTs; we also need numerical computations of Φ for 2M W times (see (69b)-(69c)).

Proof of Proposition 1
We will need the following result.
Proof of Lemma 1: Moreover, ψ(z) is closely related to Mill's ratio, and it was shown in the context therein that dψ(z)/dz ≤ 1 for all z [44]. The desired result follows.
From the expression of ∇f (θ) in (34), we have where we apply Lemma 1 to obtain the second inequality. The proof is complete.

Proof of Propositions 2 and 4-7
The proofs of Propositions 2, 4, 5, 6 and 7 follow the same principles as their PG counterparts [22,23], although a careful verification remains necessary. The following lemma, which considers the structure of our problem, is important.

Proof of Proposition 2
We start with the basic EM method. Applying Lemma 2 to the EM method (41), withθ = θ ⋆ and ε k = 0, we have for k ≥ 0. It follows that, for k ≥ 1, Also, by the basic EM property F (θ k+1 ) ≤ F (θ k ), we have

Proof of Proposition 5
Next, we consider the inexact EM scenario in (55)-(57). For convenience, let w k = k−1 i=0 (F (θ i+1 )− F (θ ⋆ )). Following the same proof in the preceding subsection, but with ε k = 0 in general, one can show that, for k ≥ 1, Note that (81) is due to the descent condition F (θ k+1 ) ≤ F (θ k ), enabled by the condition G(θ k+1 |θ k ) ≤ G(θ k |θ k ) in (57). The new challenge lies in the last term of (80). To deal with it, assume that B has full column rank. We have where we have used the result Bx ≥ σ min (B) x . Since w k ≥ 0, Eqs. (80) and (82) imply that where we denote σ min = σ min (B) for convenience. Consider the following lemma.
where λ i ≥ 0 for all i. Then, for k ≥ 1, we have Applying Lemma 3 to (83) gives for k ≥ 1, where Now we use (82) and (84) to bound the last term of (80). From (80), we have where the second inequality is obtained by applying (82) and (84); and the third inequality is due to A k ≥ A i for i ≤ k. Applying (86) and (85) to (81) leads to the final result in Proposition 5.

Proof of Proposition 4
To describe the proof for the accelerated EM method, we first write down some basic results [21]. The sequence {t k } k≥0 in (53) satisfies (a) t 2 k − t k = t 2 k+1 ; (b) t k ≥ (k + 2)/2; and (c) t k ≤ k + 1. Also, the following result is important.
where the first inequality is due to f (θ) ≤ g(θ|θ ′ ), and the second inequality is due to the descent condition (57). Summing the above equation over i = 0, . . . , k − 1 yields and further, Applying k−1 i=0 u i ≥ k min i=0,...,k−1 u i to the left-hand side of the above equation leads to the result in (92).
Second, we show that where σ min = σ min (B). From the ε-criticality condition (56), the structure of G(θ|θ ′ ) in (58), and the expressions of g(θ|θ ′ ) in (77), it can be shown that for some v k+1 ∈ ∂h(θ k+1 ) and e k+1 ≤ ε k+1 . This, together with R = (B ⊤ B) −1 , give rise to where we have used the fact that x R ≤ σ max (R) 1/2 x and σ max (R) 1/2 = 1/σ min (B). It can be verified that Let u = ψ(Bθ k ) − ψ(Bθ k+1 ) for convenience. We have where the first inequality is due to the fact that B(B ⊤ B) −1 B ⊤ is an orthogonal projection matrix, and the second inequality is due to Lemma 1. Putting the above components together, we have (93). Third, we use (92) and (93) to show the desired result. From (93), we have which is due to (x + y) 2 ≤ 2x 2 + 2y 2 . By appropriately setting the δ i 's in (92), we have min i=0,...,k−1 Combining the above two equations gives min i=0,...,k−1 which leads to the final result of Proposition 7.

A Discussion on the Applications of the Existing Convergence Analysis Results
In this section we examine and discuss the applications of the existing convergence analysis results to the EM method, with an emphasis on our problem. The EM convergence was first studied by Wu [32] in the statistics field. EM is a special case of the MM method, and there is a rich collection of literature on the convergence of the MM and related first-order methods in the optimization field; see, e.g., [33][34][35] for studies closer to the context of EM. In these MM studies we often consider a general framework that covers a variety of optimization schemes and a rich class of problems. In that regard, we should point out that we focus on a specific problem arising from MIMO detection or regression from quantized data, and exploit opportunities unique to our problem. Let us be specific. The convergence analyses by Wu [32] and Razaviyayn et al. [33] show that the EM or MM method can converge to a stationary point. They consider the non-convex case with fairly general assumptions, but the convergence rate is not considered. The MM convergence rate in the convex case was studied by Mairal [34] and Hong et al. [35]. If we apply the analysis framework by Hong et al. to the EM method (41) for our problem (33), we have the following result: for some constant C ≥ 0, where κ(B) = σ max (B)/σ min (B) is the condition number of B. Also, the application of Mairal's framework leads to for some constant R > 0. We will give the details of (94)-(95) shortly. The result in (94) is weak compared to our EM result in (42); in fact, it is even weaker than the PG result in (39). Also, as discussed in the main manuscript, the result in (95) is arguably not as good as our EM result in (42).
The framework by Hong et al. [35] gives the following result where The question is what are the values of γ and β for our problem. It is seen from (47) that γ = σ min (B) 2 . Also, it can be shown from (77) and Proposition 1 that β = 2σ max (B) 2 . Applying the above results to (96) leads to (94).
Mairal's framework [34] gives the following result To show the value of L r for our problem, we first note that f is convex and has σ max (B) 2 -Lipschitz continuous gradient (cf. Proposition 1 for the latter). From the expression of g in (47), one can easily verify that g(θ|θ ′ ) is convex w.r.t. θ and has σ max (B) 2 -Lipschitz continuous gradient w.r.t. θ, given any θ ′ . By Lemma B.9 in the supplementary material of [34], r(θ|θ ′ ) has σ max (B) 2 -Lipschitz continuous gradient w.r.t. θ. We therefore have (95). and subsequently where the second inequality is due to (99), and the last equation is obtained by putting α = 1 + ηλ into the equation. Next, we apply (100) to (97a) to obtain Assume C < 1. Then (101) can be rewritten as Applying log(z) ≤ z − 1 for z > 0 to both sides of (102) gives By setting η = 1/σ 2 1 and putting the expression of C in (101) into (103), we obtain the desired result What remains to show is that C < 1 is true for any λ > 0. Since ∇F (θ ⋆ ) = 0, and ∇F (θ) = ∇f (θ) + λθ, we have the inequality is due to (99). Consequently, and it is seen that the C in (101) satisfies C < 1 for any ε > 0. The proof is complete.
The proof of Lemma 6 will be shown later. Applying the above results to the equivalent form of EM in (106), we obtain where we have also used σ max (BW † B ⊤ ) = 1, which can be verified from (105), and ψ(Bθ) ≤ √ m (see (99)). Applying (109) to (108b), and using Bx ≥ σ r V ⊤ 1 x , we obtain ε ≥ V ⊤ Finally, it should be noted that σ r is the smallest positive singular value of B; and that V ⊤ 1 x = B ⊤ (BB ⊤ ) † Bx , which can be verified from (105). Applying the above results to (110) leads to the desired result.

Proof of Lemma 6
The proof is similar to that of Lemma 5; see, e.g., [23]. By the optimality condition of the problem in (107), θ 1 and θ 2 must satisfy Since h is convex, we have the following property g 1 − g 2 , θ 1 − θ 2 ≥ 0, for any g 1 ∈ ∂h(θ 1 ), g 2 ∈ ∂h(θ 2 ); see, e.g., [39]. Applying (114) to the above equation, we obtain where the Cauchy-Schwarz inequality is used to obtain the third equation. The proof is complete.

Impact of Noise Power Loading
As a supplemental numerical result, we examine the impact of noise power loading in (68). For brevity we consider GMAP-EM only; the other algorithms exhibit similar behaviors, as we observed by numerical study. Fig. 6 shows the BER performance for various values of the loading parameter σ 0 . We see performance loss when loading is applied; the loss increases with σ 0 . Table 3 shows the average number of iterations for various values of σ 0 and the SNR. We see that, without noise power loading, convergence is particularly slow for high SNRs; and that such undesirable effects are alleviated when loading is applied.