Expectation Propagation as Turbo Equalizer in ISI Channels

Depósito de Investigación de la Universidad de Sevilla https://idus.us.es/ This is an Accepted Manuscript of an article published by Cambridge University Press Journal of Fluid Mechanics 897 (2020), available at: https://doi.org/10.1017/jfm.2020.373 Copyright 2020. The Autors. En idUS Licencia Creative Commons CC BY-NC-ND 1 Expectation Propagation as Turbo Equalizer in ISI Channels Irene Santos, Juan Jos´ e Murillo-Fuentes, Rafael Boloix-Tortosa, Eva Arias-de-Reyna, and Pablo M. Olmos Abstract—In probabilistic equalization of channels with intersymbol interference, the BCJR algorithm and its approximations become intractable for high order modulations, even for moderate channel dispersions. In this paper we introduce a novel soft equalizer to approximate the symbol a posteriori probabilities (APP) where the expectation propagation (EP) algorithm is used to provide an accurate estimation. This new soft equalizer is presented as a block solution, denoted as block-EP (BEP), where the structure of the matrices involved is exploited to reduce the complexity order to OpLN 2q, i.e., linear in the length of the channel, L, and quadratic in the frame length, N. The solution is presented in complex-valued formulation within a turbo equalization scheme. This algorithm can be cast as a linear minimummean-squared-error (LMMSE) turbo equalization with double feedback architecture where constellations being discrete is a restriction exploited by the EP that provides a first refinement of the APP. In the experiments included, the BEP exhibits a robust performance, regardless of the channel response, with gains in the range 1.5-5dB compared to the LMMSE equalization. Index Terms—Expectation propagation (EP), BCJR, complexvalued, turbo equalization, ISI. I. INTRODUCTION SOFT or probabilistic channel equalization [1] is a technique to mitigate the interference between symbols (ISI) provoked by the dispersive nature of the channel [2], [3]. It provides the posterior probabilities of the estimated transmitted symbols given the observation, from which nowadays decoders highly benefit [4]. These two tasks, equalization and decoding, were initially considered separately, but the performance was remarkably improved by joining them into a turbo equalization scheme [5]–[7]. In turbo equalization, an equalizer and a decoder exchange information in terms of log-likelihood ratios (LLRs). After one or more iterations, the channel decoder generates the LLRs which are delivered back to the equalizer as updated a priori information. The optimal BCJR algorithm [8] computes the a posteriori probability (APP) for each transmitted symbol providing maximum a posteriori (MAP) probabilistic decisions. It works I. Santos, J.J. Murillo-Fuentes, R. Boloix-Tortosa and E. Ariasde-Reyna are with the Dept. Teor´ ıa de la Se˜ nal y Comunicaciones, Escuela T´ ecnica Superior de Ingenier´ ıa, Universidad de Sevilla, Camino de los Descubrimiento s/n, 41092 Sevilla, Spain. E-mail: {irenesantos,murillo,rboloix,earias}@us.es. P. M. Olmos is with the Dept. Teor´ ıa de la Se˜ nal y Comunicaciones, Universidad Carlos III de Madrid, Avda. de la Universidad 30, 28911, Legan´ es (Madrid), Spain. E-mail: [email protected]. This work was partially funded by Spanish government (Ministerio de Econom´ ıa y Competitividad TEC2016-78434-C3-{2-3}-R and Juan de la Cierva Grant No. IJCI-2014-19150), by the European Union (FEDER and Marie Curie Initial Training Network “Machine Learning for Personalized Medicine” MLPM2012, Grant No. 316861), and by Comunidad de Madrid in Spain (project ’CASI-CAM-CM’, id. S2013/ICE-2845). on a trellis representation, assuming perfect knowledge of the channel impulse response (CIR) and a channel with finite memory [9]. The BCJR complexity is proportional to the number of trellis branches, ML, increasing with the number of taps of the channel, L, and the size of the constellation used, M. The BCJR memory requirements per step also grow with the number of states. Therefore, for a few taps and a multilevel constellation the complexity becomes intractable. To reduce the complexity of the BCJR some suboptimal algorithms, based on performing a simplified trellis search with only Mestates, have been proposed in the literature. They can be divided into two different families. The first one consists in reducing the effective length of the CIR, such as the reduced-state BCJR (RS-BCJR) algorithm [10], which is based on the reduced-state sequence detection (RSSD) [11]– [13]. The key idea is to cancel the final channel taps, by truncating the memory of the channel, to reduce the number of states. On the other hand, other algorithms only keep the states with highest APP, i.e., unlike the previous algorithms, they perform a reduced search on the original full trellis, instead of a full search on a reduced-state trellis. This is the case of the M-BCJR algorithm [14]. Some approaches try to join both families to improve the results. We mention the M*-BCJR in [15], that outperforms both RS-BCJR and M-BCJR algorithms. These approaches have some important limitations. A first issue is that they are usually designed for some types of channels [10], [14], [15]. Secondly, they are unable to merge the paths determined by forward and backward trellises, since these paths generally do not match in both procedures. To overcome these problems, a variation of the M-BCJR algorithm is proposed in [16], where the authors use a different active state selection criterion. In [17] the channels are equalized differently according to their minimum, maximum or mixed phase nature [18]. Depending on the channel realization, they resort to the forward recursion or forward-trellis, backward-trellis or an optimized mixture of them. This approach is referred to as nonzero (NZ) completion. Thirdly, all approximated methods above ignore some paths in the trellis, which means that the explored paths tend to be considered much more reliable than they really are. To reduce this overestimation, the NZ with output saturation (NZ-OS) is proposed in [17]. The above issues can be mitigated with channel shortening approaches, such as [19], where the authors perform a full search over an optimized and reduced trellis. However, it is important to remark that the performance of this approach and approximated BCJR solutions degrades if the number of survivor paths, Me, does not grow accordingly with the total number of states. Hence they are computationally 2 unfeasible for large trellises and equalizers of the MMSE type are preferred [20]. Message passing approaches have been also investigated, see [21]–[25] and references therein. In [21], an equalizer based on the belief propagation (BP) algorithm is developed to reduce the inference complexity in sparse channels. However, the complexity of the method still grows exponentially with the size of the modulation and the number of nonzero channel interferers. In [24] the graphical model of the system is rewritten to end with a graph with loops equal or larger than 6, and the loopy BP is then applied. Its output is an approximation reported to provide good results in ISI channel equalization for BPSK modulations, although the method can be applied to other modulations. Its complexity scales linearly with the frame length and the channel memory, but quadratically with the constellation order. A different approach is proposed in [23], where a successive interference canceler applied to equalization is developed by considering the interference plus noise as Gaussian distributed. In [22] the authors develop an approach that introduces expectation propagation (EP) approximate inference [26]–[28] to incorporate into the BP algorithm the information from the BPSK symbol estimates coming from the channel decoder. In [25] the authors develop a different EP implementation for the equalizer in [22]. These two works develop particular instances of EP to project the BP messages into the right distribution (Gaussian/discrete), allowing feasible updates of the BP messages. The main difference between the EP method proposed by [22] and [25] is the procedure to tackle numerical instabilities during the EP updates related to variance parameters taking extremely small values or even negative ones. While in [22] the authors introduce a damping approach, in [25] negative variance parameters are replaced by their absolute values. At this point it is important to remark that [22], [23], [25] consider BPSK transmissions and the EP is just used to better approximate messages from the decoder to the equalizer. Soft linear equalization, such as the linear minimum-mean- squared-error (LMMSE) [29], is a suboptimal but low-cost alternative. Its complexity is dominated by the inversion of a matrix of size Nwhen all observations are processed in a block or batch approach. When no turbo scheme is used the complexity can be reduced to OpNlog Nqby exploiting the circulant nature of the involved matrix and the fast Fourier transform (FFT). However, after the feedback in the turbo equalization the matrix is no longer circulant and we can not use the FFT, yielding a cubic complexity in N[30]. In [31], some approximations are introduced to lower the complexity. Some windowed versions have also been developed in order to reduce this complexity. Specifically, in [32] a sliding window LMMSE algorithm is proposed and in [33] these results are improved by replacing the sliding window with an extending window. In [5], [32], [34] the authors also propose some approximated windowed solutions to further reduce the complexity. In this work we focus on a novel block solution to improve the equalization algorithms above when dealing with multidimensional constellations, with linear complexity in the constellation size and for any channel realization. We present an EP-based algorithm that approximates the joint posterior probability of the transmitted bits in a centralized manner, i.e., we do not use EP to project the BP messages into a different distribution, as [22], [25] would do. EP has been already successfully applied to multiple-input multipleoutput (MIMO) detection [35]–[37], low-density parity-check (LDPC) channel decoding [38], [39], tracking of flat-fading channels [40] and equalization of BPSK transmissions with message passing BP approaches [22], [25]. We exploit the key idea in [36], where compared to [35] the EP is used to better approximate the full posterior rather than to improve message passing algorithms at some points. Preliminary results for equalization of real-valued systems were discussed in [41]. In this paper we propose a turbo equalization scheme where we implement EP to obtain a complex-valued Gaussian approximation to the probability of the transmitted symbols conditioned to the received signal. We denote this approach as turbo BEP equalizer (T-BEP). We discuss the interpretation of this T-BEP as a double turbo LMMSE equalizer, where the discrete nature of the transmitted symbols is used as a first feedback. To deal with negative variance we avoid updates whenever EP provides negative variance. The computational complexity per step of the T-BEP is dominated by the inversion of an N-dimensional banded covariance matrix, which has exactly the same structure as in the turbo LMMSE scheme [30]. We exploit the structure of the matrices to reduce it to OpLN2q, i.e., quadratic complexity in Nand independent of M. In addition, we address the estimation of the mutual information between the detected and transmitted symbols to explain the obtained gain. Finally, we compare it to the most relevant approximated BCJR approaches discussed above. As a result, at low dimensional scenarios we achieve a performance close to the optimal BCJR solution, regardless of the channel realization. For multilevel constellations the BEP outperforms the turbo LMMSE equalizer (T-LMMSE). This performance is further improved with the T-BEP. Gains in the 2-5 dB range are reported for 16 and 64-QAM constellations. The paper is organized as follows. We first describe in Section II the structure and model of the communication system at hand and review the formulation for the block LMMSE turbo equalizer. In Section III the EP algorithm is introduced. Section IV is devoted to describe the novel proposed T-BEP equalizer and develop its formulation for complex numbers. We also study the convergence of the algorithm to propose values for the parameters of the EP equalizer and exploit the structure of the matrices involved to propose efficient computations. In Section V, we include several experiments to show the good performance of the T- BEP. We end with some conclusions. The following specific notation is used throughout the paper. If uis a vector, ui:jrefers to a column vector with the entries of vector uindexed by the set ti, i ´1, i ´2, ..., ju. We use u˚to denote the complex conjugate of u. To denote a normal distribution of a random proper complex vector uwith mean vector µand covariance matrix Σwe use CNpu:µ,Σq. 3 Channel Coding śModulator Channel Equalizer Demap ś´1Channel Decoder Map ś mibtctuk wk ykpEpuk|yqLEpct|yqLEpbt|yqˆmi LDpbtq LDpctqpDpukq Fig. 1: System model and turbo equalization. II. SYSTEM MODEL AND LMMSE SOLUTION A. System model In Fig. 1 we represent the discrete-time communication system with turbo equalization. It can be divided into four parts: transmitter, channel model, equalizer and turbo equalization. 1) Transmitter: A block of Kmessage bits, m“ rm1, ..., mKsJ, is encoded with a rate R“K{Vinto the codeword b“ rb1, ..., bVsJand permuted with an interleaver to c“ rc1, ..., cVsJ. An M-ary modulation is considered to obtain N“rV{log2Mssymbols, in u. Then, the block frames u“ ru1, ..., uNsJ“Rpuq`jIpuqare transmitted over the channel, where each component uk“Rpukq`jIpukq P A. Hereafter, Adenotes the set of symbols of the constellation of order |A| “ M. The mean transmitted symbol energy and energy per bit are denoted by Esand Eb, respectively. 2) Channel model: The channel is completely specified by the CIR, i.e., h“ rh1, ..., hLsJ, where Lis the length of the CIR, and the noise variance σ2 w, that we assume known at the receiver. The received signal y“ ry1, ..., yN`L´1sJP CN`L´1is given by » — — — — — – y1 . . . yN`L´1 fi ffi ffi ffi ffi ffi fl “ » — — — — — — — – h10 . . .... hL ...h1 .... . . 0hL fi ffi ffi ffi ffi ffi ffi ffi fl » — — — — — – u1 . . . uN fi ffi ffi ffi ffi ffi fl `» — — — — — – w1 . . . wN`L´1 fi ffi ffi ffi ffi ffi fl(1) or y“Hu `w(2) in matrix form, where yk“ L ÿ l“1 hluk´l`1`wk“hJuk:k´L`1`wk,(3) uk“0for kă1and kąN, and w„CN`w:0, σ2 wI˘. In this case, it is circular complex AWGN due to its zero mean [42]. 3) Equalization: Given the model above, the posterior probability of the transmitted symbol vector uyields ppu|yq “ ppy|uqppuq ppyq9CN`y:Hu, σ2 wI˘ N ź k“1 IukPA,(4) where IukPAis the indicator function that takes value one if ukPAand zero otherwise. The optimal equalizer to estimate ppuk|yq,k“1, . . . , N, is the BCJR algorithm, unaffordable for multilevel constellations when the channel memory, L, grows. 4) Turbo Equalization: The key point of turbo equalization is the iterative exchange of information between the equalizer and the decoder for the same set of received symbols [5], [7], [34]. Given the extrinsic log-likelihood ratios (LLR) from the equalizer to the decoder, LEpbt|yq, the latter computes (after one or more iterations) an estimation of the information bits, ˆm, and an extrinsic LLR on the coded bits LDpbtq “ log ppbt“0|LEpb|yqq ppbt“1|LEpb|yqq´LEpbt|yq.(5) These LLRs are mapped again and delivered back to the equalizer as updated a priori probability, with a slight abuse of notation we denote it by pDpuq. This process is repeated iteratively for a given maximum number of iterations, T, or until convergence. B. Turbo LMMSE Equalization Algorithm We include next the description of the turbo equalization with LMMSE, that we denote by T-LMMSE, since we will use it as benchmark and its structure is a good starting point to develop our solution. Given the CIR, the posterior approximation provided by the LMMSE equalizer is obtained by replacing the discrete uniform prior ppuqin (4) by a product of independent Gaussian distributions with mean Eruks P C and variance Vruks P R`, qMMSEpuq “ CN`y:Hu, σ2 wI˘ N ź k“1 CNpuk:Eruks,Vruksq. (6) This distribution is a proper complex Gaussian [42] qMMSEpuq “ CNpu:µMMSE,ΣMMSEq(7) where ΣMMSE “`σ´2 wHHH`Σ´1 u˘´1,(8) µMMSE “ΣMMSE `σ´2 wHHy`Σ´1 uµu˘(9) and µu“ rEru1s, ..., EruNssJ,Σu“ diag pVru1s, ..., VruNsq. The complexity of this solution is dominated by the matrix inversion in (8). For the first iteration, we set Eruks “ 0and Vruks “ Es, and the matrix yields a circulant matrix whose inverse can be computed with complexity of order OpNlog Nqby means of the FFT. Then, the symbol probability of each entry is computed by independently deciding on each component. In turbo equalization, the extrinsic information is passed to the channel decoder. The extrinsic information for symbol kis computed from (7), assuming equally probable symbols 4 for the a priori information of the k-th one. By using this turbo scheme in the receiver, the equalizer is fed back with the statistics Eruksand Vruks, which are obtained from the updated a priori probability pDpukqin Fig. 1 as Epukq “ ÿ uPA u¨pDpuk“uq,(10) Vpukq “ ÿ uPApu´Eruksq˚pu´Eruksq¨pDpuk“uq.(11) At this point, it is important to remark that since the variances are no longer equal, the matrix to invert in (8) is not circulant. Therefore, the complexity of its inversion is of order OpN2Lq, as discussed later in this paper. III. EXPECTATION PROPAGATION Expectation propagation [26]–[28], [43] is a technique in Bayesian machine learning to approximate an intractable probability distribution, in which inference is unfeasible, by exponential family distributions. Suppose we are given some statistical distribution with hidden variables xand observables Dthat factors as ppx|Dq 9 fpxq I ź i“1 tipxq,(12) where fpxqbelongs to an exponential family Fwith sufficient statistics Φpxq, and tipxqare nonnegative factors that do not belong to F, making direct inference over (12) not possible. EP provides a feasible approximation to ppx|Dqby an exponential distribution qpxqfrom Fwhich factorizes as qpxq 9 fpxq I ź i“1 ˜ tipxq,(13) where factors ˜ tipxq P Fare optimized to achieve an accurate global approximation qpxq Ð ppx|Dq, which optimally satisfies EqpxqrΦpxqs “ Eppx|DqrΦpxqs. This is known as the moment matching solution. A feasible algorithm to approximate this solution is the sequential EP algorithm [26], [27], which optimizes each factor ˜ tipxqin turns independently in the context of all of the remaining factors. A sketch of the EP algorithm is given in Algorithm 1 where qr`spxqis the approximation to qpxqin (13) at iteration `. IV. BLOCK-EP TURBO EQUALIZER The EP is endowed with a great flexibility, given by the model in (13). In order to improve the accuracy of the EP solution, it is important to retain as much structure as possible from the true distribution and separate it from the latent (unknown) factors, tipxq[26]. Bearing this in mind, we develop an EP approximation to (4), namely the posterior distribution of the transmitted symbols given the channel outcome y. A. The BEP equalizer The following Gaussian exponential family will be considered to find a suitable approximation to (4): qpuq 9 CN`y:Hu, σ2 wI˘ N ź k“1 exp pu˚ kγk`γ˚ kuk´Λku˚ kukq, (15) Algorithm 1 The EP algorithm Initialize approximating factors ˜ tipxqand then qpxqin (13). repeat for i“1, ..., Ido 1) Compute the distribution ˜pipxq 9 tipxqqr`szipxq “ tipxqqr`spxq{˜ tr`s ipxq(14) and its moments, where qr`szipxqis the so called cavity function. 2) Compute the refined factor ˜ tr``1s ipxqby setting the moments of the distribution ˜ tr``1s ipxqqr`szipxqequal to the moments of ˜pipxq. end for until convergence (or stopped criterion) where the product of indicator functions is replaced by a product of univariate proper complex Gaussians, each parameterized by a pγk,Λkqpair, k“1, . . . , N. For any value γkPCand ΛkPR`,qpuqis also proper complex Gaussian CNpu:µ,Σqwith Σ“R´1“ pσ´2 wHHH`diagpΛqq´1,(16) µ“Σpσ´2 wHHy`γq,(17) exhibiting a similar structure to the LMMSE in (7). Based on these definitions in mind, we adapt the EP algorithm in Algorithm 1 to our setting. We denote the resulting algorithm block-EP (BEP) soft equalizer, and it can be found in Algorithm 2. A few remarks: ‚Step 1) is the “LMMSE” step of the algorithm, as it requires computing Σand µin (16) and (17) for the current configuration of pΛ,γq. Note that this step is equivalent to the LMMSE method in (10) and (11). ‚Steps 2.1)-2.3) can be done in parallel for k“ 1, . . . , N. They can be seen as a refinement to the estimate computed at Step 1) by enforcing a discrete distribution. Given the factorization in (15), marginal qr`spukq,k“1, . . . , N is proportional to exp pu˚ kγk`γ˚ kuk´Λku˚ kukq. This term is canceled out in ˆpr`spukqand the“true” discrete factor is introduced. Parameters canceled are recomputed by moment matching in (22)-(23). Eqn. (22) and (23) are proposed following the guidelines in [36, Eq. 35-36]. The parameter update in (22) may return a negative value for some k’s. For those k’s, we keep the values from the previous iteration. We introduce a smoothing parameter βP r0,1sand a small constant . To avoid numerical instabilities, constant  is the minimum allowed variance at each iteration, i.e., σ2r`s pk“maxp, σ2r`s pkq. A block diagram of the BEP equalizer is included in Fig. 2 (a). The LMMSE block corresponds to Step 1). The grey block represents the refinement of the current marginals qr`spukq through projection over discrete alphabet, whose output is fed back to the LMMSE block. 5 LMMSE Demap ykppuk|yqLpbt|yq Sloops ukPAChannel Decoder ˆmi Block-EP equalizer LMMSE Compute pEpuk|yq + Demap Map ykppuk|yqLEpbt|yq pDpukq Sloops Tloops ukPAChannel Decoder ˆmi LDpbtq Inner Turbo Equalizer (BEP) Outer Turbo Equalizer (a) (b) Fig. 2: In (a), BEP equalizer block diagram. In (b), Turbo BEP block diagram. B. The Turbo BEP equalizer The BEP equalizer can be further improved by using a turbo scheme. The resulting proposed equalizer, denoted as T-BEP, Algorithm 2 Block-EP equalizer (BEP) Input: pγr1s k,Λr1s kqinitialization. for `“1, ..., S do 1) Calculate the moments of qr`spuqin (15) for the current values of γkÐγr`s kand ΛkÐΛr`s k. for k“1, ..., N do 2.1) Compute the k-th marginal of the distribution qr`spuq, denoted as qr`spukq “ CN´uk:µr`s k, σ2r`s k¯, and the distribution qr`szkpukq “ qr`spukq exp ´u˚ kγr`s k`γr`s k ˚uk´Λr`s ku˚ kuk¯∼ ∼CN´uk:zr`s k, v2r`s k¯,(18) namely cavity marginal function, where v2r`s k“σ2r`s k 1´σ2r`s kΛr`s k , zr`s k“v2r`s k˜µr`s k σ2r`s k´γr`s k¸. 2.2) Obtain the distribution ˆpr`spukq 9 qr`szkpukqIukPA and estimate its mean µr`s pkand variance σ2r`s pk. 2.3) Set the mean and variance of the unnormalized Gaussian distribution qr`szkpukqexp ´u˚ kγr``1s k,new `γr``1s k,new ˚uk´Λr``1s k,newu˚ kuk¯ (19) equal to µr`s pkand σ2r`s pk. To this end, compute: Λr``1s k,new “´σ´2r`s pk´v´2r`s k¯,(20) γr``1s k,new “´µr`s pkσ´2r`s pk´zr`s kv´2r`s k¯.(21) 2.4) Update the values as Λr``1s k“βΛr``1s k,new `p1´βqΛr`s k,(22) γr``1s k“βγr``1s k,new `p1´βqγr`s k.(23) end for end for With the values γrS`1s,ΛrS`1sobtained after EP algorithm, calculate the final distribution qpuqin (15). Algorithm 3 Block-EP turbo equalizer (T-BEP) 0) Initialize pγr1s k,Λr1s kq“p0, E´1 sq. for t“1, . . . , T do 1) With the current initialization to pγr1s k,Λr1s kq, run BEP equalizer in Algorithm 2. 2) Compute an estimate to the extrinsic LLRs, LEpbt|yq, by feeding pEpuk|yq “ qrS`1szkpukq(24) to the demapper. Feed LEpbt|yqto the channel decoder. 3) From the channel decoder per-bit soft output, recompute a probability distribution for each symbol pDpukq from the decoder and compute its mean Eruksand variance Vruksgiven by (10) and (11). 4) Re-initialize pγr1s k,Λr1s kqto pEruksVruks´1,Vruks´1q. end for can be easily described as the T-LMMSE in Subsection II-B by just replacing (7) with the result of the BEP, Algorithm 2, in (15). A detailed implementation of the T-BEP is included in Algorithm 3. Also, a block diagram is included in Fig. 2(b). Note that the T-BEP can actually be seen as a turbo equalizer with two loops. First we run an inner turbo scheme, i.e. BEP. After Siterations of the BEP iterative procedure, the extrinsic LLR is given to the decoder in an outer loop, which is repeated Ttimes. We propose to approximate the extrinsic probabilities by the cavity functions at the end of the EP algorithm, see (18). In this second stage the restrictions from the channel coding are exploited. The output of the channel decoder is used to initialize the BEP iterative procedure, whose outputs are then fed forward to the channel decoder. This is a major difference with respect to the previous scheme used in [5], [32], [34], where the estimates were refined only using the output of the decoder, as shown in Fig. 1, and the proposed inner loop was not present. C. Efficient implementation All pγr`s k,Λr`s kqpairs for k“1, ..., N can be updated in parallel. The most involved step is the computation of an N-dimensional inverse matrix in (16) for each `-iteration, whose complexity is dominated by its size, i.e., OpN3q. Once this inverse is computed, the parallel update of all pairs pγr`s k,Λr`s kq Ð pγr``1s k,Λr``1s kqby means of step 2 and 3 in 6 Algorithm 3 has a smaller computational complexity, linear in NM. To reduce the complexity of the matrix inversion, we propose to exploit the banded structure of the channel matrix along with the short length of the channel compared to N. The matrix Rin (16) is a symmetric, positive-definite and banded matrix with bandwidth 2L´1. We can decompose Rusing the band Cholesky factorization [44] such that R“GGJ,(25) where Gis a lower triangular banded matrix with bandwidth L that can be computed with NL2operations. Then, the inverse of the covariance matrix can be rewritten as Σ“R´1“G´1JG´1.(26) We invert matrix Gby Gauss-Jordan elimination. For every diagonal element, say Gk,k, we divide row kof less than N non-null elements by Gk,k and cancel all the L´1lower elements of its column. Repeated for the whole diagonal yields a complexity of OpN2Lq. D. Convergence Although the convergence is not guaranteed, we concluded empirically that in about S“10 iterations the distribution qpuqconstructed in (15) tipically reaches a stationary value. Experimental results show that controlling numerical instabilities in the parameter updates in the turbo case (Tą0) is simply done by setting “10´9and β“0.1. This solution is robust regardless the constellation order, SNR or channel realization. In the equalization case (T“0), we have found that results can be improved if is first kept constant to a relatively high constant (0.5) to then reduce it. More precisely, we have used “2´maxp`´5,1q. For equalizing 64- QAM constellations, the best performance has been reported by setting fixed to “0.9. We have selected the previous values after extensive experimentation as a trade off between convergence speed and accuracy. We emphasize again that the latter heuristics to control stability in the equalization case are not needed if the turbo scheme is used, as the feedback loop naturally stabilizes the BEP output. In Fig. 3, we include a representative example of the convergence by depicting the evolution of some components of the mean vector µin (17) and covariance matrix Σin (16) along different values of S in the low-Eb{N0regime (specifically, Eb{N0“3dB) for a given observation, y, with a 16-PAM modulation, L“4 and T“0. As shown in Fig. 3, approximately after S“10 iterations the EP equalizer reaches a stationary value for the mean and variance, equal to the value provided by the optimal BCJR approach. E. Performance analysis We compare probabilistic equalizers using the mutual information between the transmitted symbol ukand detected symbol ˆuk, distributed according to the estimation of the posterior distribution of ukgiven y: Ipuk,ˆukq “ ÿ ukÿ ˆuk ppuk,ˆukqlog2 ppˆuk|ukq ppˆukq,(27) 0 5 10 15 20 25 30 ´3 ´2 ´1 EP iteration ℓ (a) Mean BEP BCJR 0 5 10 15 20 25 30 7 8 9 10 EP iteration ℓ (b) Variance BEP BCJR Fig. 3: Evolution of the mean (a) and variance (b) of 3 randomly chosen entries of the approximate posterior as EP iterates with 16-PAM, L“4and Eb{N0“3dB. 0 5 10 15 20 25 1 2 3 4 5 6 Eb{N0(dB) Ipuk,ˆukq BEP LMMSE 10 11 12 13 3 3.5 4 4.5 5 Fig. 4: Mutual information for BEP and LMMSE for 64-QAM and random channels with L“7. where note that ppuk,ˆukqis the joint probability distribution of ukand ˆukafter marginalizing the channel output y, the channel impulse response (in the case we consider it random) and the rest of symbols in the sequence u. We resort to Monte Carlo estimates to generate samples from ppuk,ˆukqand then evaluate Ipuk,ˆukqfrom these samples. We assume that the channel taps are Gaussian distributed, but are perfectly known at the receiver. First, we collect NPZ`samples from the joint distribution of u,h,yand ˆukusing standard sampling techniques in directed graphical models [43]. The key aspect is to note that, during the sampling procedure, ppˆuk|yqis dependent on the considered probabilistic detection method. After Nsamples of possible puk,ˆukqpairs have been collected, we use these samples to estimate ppˆukq,ppˆuk|ukq and, finally, Ipuk,ˆukqin (27). The higher Ipuk,ˆukqis for each probabilistic equalizer, the closer we perform to channel capacity. In Fig. 4 we depict the mutual information in (27) 7 computed for N“106samples per Eb{N0point and averaged over all the symbols in the transmitted sequence, for the BEP and LMMSE algorithms, considering a 64-QAM constellation and channels of L“7complex-valued taps. Observe that, for very small Eb{N0, the noise is so large that both methods achieve a small mutual information. On the other hand, for very large Eb{N0values, both methods converge to 6 bits, i.e., the number of bits transmitted per QAM symbol. Any probabilistic method will eventually saturate to this value. However, the key aspect is to be able to design a probabilistic equalizer able to improve the mutual information for intermediate Eb{N0values. Observe that before the saturation, the BEP achieves a gain w.r.t. LMMSE of around 1.5dB. F. Computational Complexity A detailed comparison of the complexity for the T-BEP and the T-LMMSE is included in Table I. We also include the computational complexity of the T-BCJR and its approximated approaches with turbo. From the computational point of view, as Mand/or Lgrow, the BCJR and approximated approaches are unaffordable. Algorithm Complexity T-BEP S′(LN2+NM) + α+S′T(LN2+α) T-LMMSE Nlog N+α+T(LN2+α) T-BCJR (T+ 1)(NML+α) Approx. T-BCJR (T+ 1)(NMeM+α) TABLE I: Complexity comparison between algorithms, where S1“ S`1and αis the complexity of the LDPC decoder. V. EXPERIMENTAL RESULTS In this section, we illustrate the good performance of the BEP equalizer for different scenarios. Each channel tap is i.i.d. complex circular Gaussian distributed with zero mean and variance equal to 1{L. The channel response is normalized. We average the BER over 1000 random frames per channel realization. We limit to 5 the absolute value of LLRs given to the decoder in order to avoid very confident probabilities which negatively affect its estimations. In the following, when mentioning approximated BCJR algorithms we refer to the M- BCJR [14], M*-BCJR [15], RS-BCJR [10], NZ and NZ-OS [17] solutions. We denote by NZ the approach consisting in running FT, BT or DT-NZ in [17] depending on the phase of the channel. If output saturation is used we refer to them as NZ-OS. We use a (3,6)-regular LDPC of rate 1/2, for a maximum of 200 iterations using the belief propagation as decoder [45], [46]. The codes were generated using the progressive edge-growth algorithm [47]. The maximum number of iterations of the LDPC decoder is set to 100 in every iteration of the turbo equalizer. A. BEP soft equalizer In Fig. 5 and Fig. 6 we show a comparison for some typical channels found in the literature using codewords of V“1024 bits. In Fig. 5 we simulate the following scenario in [17]: 6.5 7 7.5 8 10´3 10´2 10´1 Eb{N0(dB) BER Fig. 5: BER for LMMSE (Ź), BEP (˝), BCJR (˛), M-BCJR (˝), M*-BCJR (˚), RS-BCJR (‹), NZ (Ÿ) and NZ-OS (‚) equalizers for BPSK and the minimum phase channel h“1 ?140 r7 6 5 4 3 2 1sJ. 6 6.5 7 7.5 8 10´3 10´2 10´1 Eb{N0(dB) BER Fig. 6: BER for LMMSE (Ź), BEP (˝), BCJR (˛), M-BCJR (˝), M*-BCJR (˚), RS-BCJR (‹), NZ (Ÿ) and NZ-OS (‚) equalizers for BPSK and the maximum phase channel h“1 ?140 r1 2 3 4 5 6 7sJ. minimum phase channel h“ r7654321sJ{?140, BPSK symbols and Me“4states out of 64. In Fig. 6 we simulate the maximum phase channel h“ r1 2 3 4 5 6 7sJ{?140 and BPSK symbols with Me“8states for the approximated solutions as in [17]. To study the performance for other channels, in Fig. 7 and Fig. 8 we include the averaged BER over 100 random channels with L“5real-valued taps. We fix to Me“8the number of states for the approximated solutions. In Fig. 7, we consider BPSK modulation, hence the BCJR has 16 states per step, while in Fig. 8 a 4-PAM is used, increasing the number of states to 256. In Fig. 5 and Fig. 6 we can observe that approximated methods M-BCJR, M*-BCJR, RS-BCJR are quite sensitive to the channel realization. This is not only caused because they are based on just a forward or a backward strategy, but also because their parameters need to be tuned according to the particular channel. For this reason, the approximations M- BCJR, M*-BCJR and RS-BCJR fail with maximum phase 8 5 6 7 8 9 10 11 12 10´4 10´3 10´2 10´1 Eb{N0(dB) BER Fig. 7: BER for LMMSE (Ź), BEP (˝), BCJR (˛), M-BCJR (˝), M*-BCJR (˚), RS-BCJR (‹), NZ (Ÿ) and NZ-OS (‚) equalizers for BPSK and 100 random channels with L“5. 10 12 14 16 18 20 10´5 10´4 10´3 10´2 10´1 Eb{N0(dB) BER Fig. 8: BER for LMMSE (Ź), BEP (˝), BCJR (˛), M-BCJR (˝), M*-BCJR (˚), RS-BCJR (‹), NZ (Ÿ) and NZ-OS (‚) equalizers for 4-PAM and 100 random channels with L“5. channels, as observed in Fig. 6. For averaged channels, see Fig. 7, NZ and NZ-OS exhibit a good performance at a low computational complexity due to the low state reduction ratio (ML{Me). However, these solutions fail if Medoes not grow accordingly, as can be observed in Fig. 8. This means that if the number of states, Me, is kept as a fraction of the total number, ML, they may exhibit a good performance. But if a large number of states is unaffordable, a fraction of it becomes also intractable. The LMMSE in these experiments does not fail as the approximated approaches do, but its performance degrades significantly compared to the BER of the BCJR in Fig. 8. Finally, the BEP exhibits a quite robust behavior, closer to the BCJR performance. Note here that the number of states is low, and that the BCJR can be used with optimal results at a low complexity. For multilevel constellations and a channel with a few taps, the BCJR or their approximations are no longer computationally affordable. In the case of large dimensions, filter-based equalizers of the MMSE type are preferred [20]. In Fig. 9, 8 10 12 14 16 18 20 22 24 10´4 10´3 10´2 10´1 Eb{N0(dB) BER LMMSE BEP CS Fig. 9: BER for LMMSE (Ź), BEP (˝) and CS (‹) equalizers for 16-QAM (solid lines) and 64-QAM (dashed lines) and 100 random channels with L“7. we depict the BER curves for BEP, LMMSE equalization and the algorithm in [19] that we denote as channel shortening (CS), considering codewords of V“4096 bits, 100 random channels of L“7complex-valued taps and two different modulations: 16-QAM (solid lines) and 64-QAM (dashed lines). The CS algorithm has been simulated with a reduced memory ν“2(ν“1) for the 16-QAM (64-QAM) case. We observe in this experiment that the BEP exhibits significant improvements with respect to the LMMSE and CS solution for both constellations. B. Turbo BEP To improve the estimates the turbo-equalization can be used. We simulate scenarios with large dimensions, where the only computationally feasible algorithms are the BEP and LMMSE solutions. In Fig. 10 we simulate a 16-QAM constellation, averaging over 100 random channels with a large memory, L“20 complex-valued taps and codewords of V“1024 bits. We represent the first 3 iterations of the turbo scheme, since we found no further improvement for Tě3. In Fig. 11 and Fig. 12, we depict the BER curves for turbo BEP and LMMSE equalization, averaged over 100 random channels of L“7 complex-valued taps, 64-QAM modulation and codewords of V“1024 and V“4096 bits, respectively, for T“3. The results are remarkable. The BEP with no feedback from the decoder is 1.5-4.5dB away from the LMMSE equalizer for BER=10´4, depending on the channel length, L, and the constellation size, M. It even outperforms the turbo LMMSE equalizer. The BEP estimation can be further improved using the turbo scheme, in about 1or 1.5dB depending on the codeword length. VI. CONCLUSIONS When the number of states involved in the BCJR equalization or its approximations is high, their computational complexity is unaffordable. Approximated strategies are interesting whenever the number of states is not too high. Therefore,