{'fraction_non_alphanumeric': 0.046527407551010515, 'fraction_numerical': 0.0156787670236999, 'mean_word_length': 4.289284798769546, 'pattern_counts': {'":': 0, '<': 0, '<?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': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Asymptotic Performance Analysis of a K-Hop Amplify-and-Forward Relay MIMO Channel

The present paper studies the asymptotic performance of multi-hop amplify-and-forward relay multiple-antenna communication channels. Each multi-antenna terminal in the network amplifies the received signal, sent by a source, and retransmits it upstream towards a destination. Achievable ergodic rates under both jointly optimal detection and decoding and practical separate decoding schemes for arbitrary signaling schemes, along with the average bit error rate for various receiver structures are derived in the regime where the number of antennas at each terminal grows large without a bound. To overcome the difficulty of averaging over channel realizations we apply large-system analysis based on the replica method from statistical physics. The validity of the large-system analysis is further verified through Monte Carlo simulations of realistic finite-sized systems.


I. INTRODUCTION
M ULTIPLE-input multiple-output (MIMO) relaying has been proved a promising technology that enables reliable communication with increased coverage and data rates [1], [2]. A wireless sensor network serves as a practically relevant example, where multi-hop (and cluster-based) relaying helps to overcome the low-power-budget constraints [3]. Thus, understanding of the fundamental limits of multi-hop relay MIMO channels has been regarded as an important milestone in the research efforts within the field of complex cooperative networks [4], [5]. In spite of the considerable efforts, however, the capacity of a multi-hop MIMO relay channel, in its most general formulation, remains an open problem.
The relays in a cooperative network may realize different cooperative strategies. A regenerative strategy (e.g., decodeand-forward or compress-and-forward [6], [7]) involves decoding or quantization of the received signal, re-encoding the underlying message and subsequent retransmission upstream. A non-regenerative strategy (e.g., amplify-and-forward The research leading to these results has received funding from the Swedish Research Council (VR). The material of this paper was presented in parts at the IEEE International Symposia on Information Theory in 2013, Istanbul, Turkey, July 2013 and Honolulu, U.S.A., July 2014. The work was done while M. A. Girnyk  (AF) [4]) involves direct amplification of the received noisy signal with subsequent retransmission towards the destination. Non-regenerative relaying, being simple in implementation and independent of the modulation schemes at the source terminal, is of particular interest.
The present paper is devoted to a general K-hop AF relay MIMO channel in the presence of fast fading, whose achievable rates are characterized by the end-to-end ergodic mutual information (MI). Efficient evaluation of the latter is problematic due to computation of the expectation of the MI over the channel realizations, as well as over non-Gaussian input signals in the case the latter are used.
To overcome these difficulties, several asymptotic approaches, based on the large-system assumption, have been recently proposed. In [8], various asymptotic limits are considered; namely, where the number of antennas grows large either at the source, the destination or the relay terminal, while the number of antennas at the other terminals stays fixed. In other works, the authors apply techniques from random matrix theory to obtain an explicit approximation of the ergodic MI. For instance, in [9], such an explicit expression for the achievable rate of a two-hop AF relay MIMO channel with Gaussian inputs is obtained. Further, assuming that the number of antennas at each terminal within the relay MIMO network grows very large, the authors of [10] analyze a multi-hop relay channel with noise only at the destination terminal. Meanwhile, in [11] the capacity of a general K-hop AF relay MIMO channel with Gaussian inputs and noise at every hop is expressed in terms of the limiting eigenvalue distribution of a product of MIMO channel matrices, obtained via a set of recursive equations. In [12,Sec. 3.5], a recursive expression for ergodic MI between the input and output of a K-hop AF relay MIMO channel is derived by means of deterministic equivalents introduced in [13]. Another recent paper, [14], considers a similar setup while evaluating the average throughput of a multi-hop network under the assumption of full channel state information (CSI) and an SVD-based linear precoder applied at the source terminal.
Multi-hop AF networks have also been considered in terms of reliability. The latter is usually evaluated in terms of the average bit error rate (BER), which, too, yields no closedform expression for the AF relay MIMO scenario with fading. Thus, in [15], the instantaneous BER of BPSK transmission is derived in closed-form and the average performance is then obtained through subsequent Monte Carlo estimation.
Most of the existing work on closed-form average BER is, however, limited either to one-antenna terminals [16], [17], or to bounds on the BER performance [18]. The authors of [16] derive a closed-form expression for average BER of a dual-hop relaying system using various finite-alphabet constellations via the harmonic mean. An exact closed-form expression for the average BER of single-antenna AF relay system in terms of Lauricella multivariate hypergeometric functions is provided in [17]. In turn, in [18], an upper-bound on the average BER of a MIMO AF relay channel is derived using the arithmetic-geometric mean inequality. Furthermore, the largesystem BER performance of the linear minimum mean-square error (MMSE) detector in the dual-hop AF relay scenario is analyzed in [19]. Nonetheless, the problem in its general formulation, accounting for other conventional detectors and arbitrary number of hops, K, remains open. The above being said, the present paper provides a framework for efficient performance analysis of a K-hop AF relay MIMO communication channel under general conditions. Namely, • We provide an explicit expression for the end-to-end ergodic MI for a multi-hop channel with arbitrary signaling at the source terminal in the regime, where the numbers of antennas at each terminal grow without bounds at constant ratios. This allows the evaluation of the system spectral efficiency under optimal joint detection and decoding (JDD). For Gaussian signals, our results degenerate to those reported in [12,Sec. 3.5], [20]. For non-Gaussian signals our results partially reduce to those in [19], [21]. • In addition, we show that the decoupling principle, reported by Guo and Verdú [22] in the context of a DS-CDMA system, holds also for the multi-hop AF relay MIMO setting. Namely, in the large-system regime, an AF relay MIMO channel, where joint spatial detection at the receiver of the destination terminal is followed by separate decoding (SD), decouples into a bank of scalar Gaussian per-stream channels. This allows the characterization of the system performance in the corresponding practically motivated scenario. The obtained results degenerate to those presented in [19], [23], [24]. • Furthermore, the framework allows us to determine the average uncoded BER of the system for various types of detection schemes. The corresponding results degenerate to those provided in [19], [23], [24], [25]. The results also align with those presented in [22], [26], [27] for CDMA and MIMO-CDMA systems. Our analysis is based on the replica-symmetric (RS) ansatz of the replica method from statistical physics. The method was invented in early 50's by Kac [28], and it provides a powerful framework for efficient analysis of macroscopic quantities of large many-body systems (e.g., spin glasses [29], [30]). Albeit not yet rigorously justified, the method proves efficient in the cases where all other methods fail (e.g., traveling salesman problem [31]). It was introduced to the field of communications by Tanaka [26], which inspired lots of subsequent research efforts [22], [32], [25], [27], [33]. The method provides a powerful framework for efficient performance analysis of average performance of channels described by a linear vector model (such as CDMA or MIMO systems 1 ) under general conditions. The method has also gained popularity in other fields of engineering, such as compressed sensing [36], watermarking [37] and machine learning [38].
The remainder of the paper is organized as follows. In the following section, we describe the system model and formulate the problem of interest. Next, in Section III, we present the main results of the paper, viz., an expression for the Helmholtz free energy, along with the decoupling principle, for an equivalent statistical-physics system. In Section IV, we then provide an application of the main result to the AF relay MIMO channel of interest, which enables the evaluation of the performance in terms of the achievable ergodic rate and average BER for various detection schemes. Then, Section V presents the results of numerical simulations alongside relevant discussion. Finally, in Section VI, we draw conclusions. The derivations of the claims are postponed to the appendices.
Notation: Throughout the paper we use upper case boldfaced letters to denote matrices, e.g., X, with elements denoted by [X] i,j , lower case bold-faced letters to denote column vectors, e.g., x, with elements x i , and lower case light-faced letters to denote scalar variables, e.g., x. Superscripts (·) T and (·) H denote the transpose and Hermitian adjoint operators, respectively. Meanwhile, tr{X} and det(X) denote the trace and the determinant of matrix X. Also, I M , 0 M and 1 M denote the identity matrix, the all-zeros vector and the all-ones vector of length M . Operator E{·} denotes the expectation, δ(·) denotes the Dirac delta function, while ⊗ represents the Kronecker product. Operators Re{·} and Im{·} provide the real and imaginary parts of the argument, respectively. With ρ being the signal-to-noise ratio (SNR), we denote the ergodic MI as I(ρ), the differential entropy of the received signal as h s (ρ) and the conditional differential entropy as h n (ρ). The respective quantities with bars on top denote the corresponding asymptotic approximation based on the RS ansatz. Finally, I η (ρ) denotes the MI associated with an equivalent decoupled scalar channel with inverse noise variance η.

II. SYSTEM MODEL
Consider a K-hop channel, consisting of K − 1 multiantenna relay terminals that assist a multi-antenna source to communicate with a multi-antenna destination. The corresponding setup is depicted in Fig. 1, and it corresponds to that of [11], [12,Sec. 3.5]. There is no direct link between the source and the destination, and the terminals operate under a time-division multiple-access (TDMA) protocol, so that a single transmitter-receive pair is active at a given time. 2 Furthermore, we assume that each relay receives only the signals from the preceding hop. Namely, a symbol sent by the source has to traverse all the K hops before it reaches the destination, where K is a fixed finite number. The source, the destination and the kth relay terminal are equipped with M 0 , M K and M k antennas, respectively. AF relaying is employed at each relay, so that the relay simply amplifies and retransmits the received signal upstream without decoding it.
In our flat-fading model, the received signal at the kth terminal is given by where y k−1 and y k are the input and output of the kth channel, respectively. Moreover, ρ k is the SNR at terminal k and β k−1 is the normalization constant chosen so that the long-term transmit power constraint at terminal k − 1 is satisfied, that is, For later convenience, let us define also a set including all the hops K {1, . . . , K} and group the corresponding SNR values into a vector ρ [ρ 1 , . . . , ρ K ] T . The channel matrix between terminals k − 1 and k, H k , is assumed to have circularly symmetric complex Gaussian (CSCG) ran- The end-to-end input-output relation of the K-hop channel can be written as follows: where x y 0 is the input to the relay channel, assumed for the sake of simplicity to have i.i.d. zero-mean unit-variance components distributed according to some probability density 2 The assumption is motivated by the following reasoning. If all nodes that are not transmitting were to listen for a transmission at all times, backward links in the network would appear. This would lead to an additional problem of scheduling of the transmissions which is outside of the scope of the present paper. function p(x), while y y K , and the corresponding matrices are defined as 3 In this paper, we assume that the destination has full channel state information (CSI), 4 whereas the source and the relays are not aware of the channel states. Hence, the longterm maximum achievable rate is given by the average MI between the input and output of the channel. Define the set of channel matrices H {H 1 , . . . , H K } and the set of noise realizations N {n 1 , . . . , n K−1 }. The end-to-end (normalized) achievable ergodic rate of the channel for a given input distribution p(x) is where the differential entropy terms are given by 5 with the conditional distribution of the channel being Although (5) represents the achievable rate of a relay MIMO channel, it assumes optimal joint detection and decoding (JDD), which in practice may be prohibitively complex. A more plausible alternative involves joint spatial detection, followed by a bank of single-user decoders and is referred to as separate decoding (SD) hereafter. In this case, the detector estimates the symbol vector based on the generalized posterior mean estimator (GPME) [22], given by where subscript q reflects the fact that the receiver uses some postulated channel law q(y|x , H) and distribution q(x ) for postulated inputs x . The GPME is, in general, suboptimal if the latter two do not match the conditional density p(y|x, H) and the prior distribution p(x) of the actual channel (3), respectively. However, as was shown in [22], most of the practically relevant detectors can be regarded as a GPME, optimal for the postulated channel law q(y|x , H). It was further shown that in order to capture suboptimality of the aforementioned detectors it suffices to postulate a channel with 3 The matrix product notation is defined as following [39], where a similar multi-hop channel was analyzed in terms of the eigenvalue distribution. 4 CSI acquisition at the destination can, in principle, be done via conventional pilot-based methods [40]. Moreover, specific algorithms for pilot design in AF relay MIMO systems, which allow for estimation of the per-hop channels directly at the destination, have recently been proposed [41], [42]. The case where the receiver is ignorant of the intermediate channels leads to mismatched decoding [43], [44]. For an application of covariance mismatched decoding in the context of large-scale MIMO systems, see for example, [45]. 5 To account for the TDMA protocol, a factor of 1/K is applied to (5). a mismatch only in the noise variance σ 2 and prior distribution q(x ), i.e., where x is a postulated channel input vector, and n k ∼ CN (0 M k , σ 2 I M k ), k ∈ K are postulated CSCG noise vectors with variance σ 2 . That is, the conditional density of the postulated channel (9) is given by (10) The corresponding achievable communication rate of the SD scheme is given by An individually optimal maximum a posteriori probability (MAP) detector performs an exhaustive search over all possibly transmitted vectors [46], yielding the estimatê Unfortunately, the implementation of such a detector is computationally prohibitive in practice and hence from the practical view-point reduced-complexity alternatives are preferable. The three conventional linear detection schemes considered in this paper are [47]: • Matched filter (MF), which maximizes the SNR at the output of the detector, disregarding the interference between the streams. • Zero forcing (ZF) filter, which removes the interference, while at the same time enhancing the noise. • Linear MMSE (LMMSE) filter, which minimizes the mean-square error (MSE) without constraints on interference, being an optimal linear detector. In general, average performance of these detectors over fading channels exhibits no closed-form expression. The analysis is further complicated by the fact that in the AF relay MIMO setting the resulting channel matrix in (3) becomes a product of the channel matrices of the K hops, yielding non-i.i.d. entries. Note here that in [39] the asymptotic eigenvalue distribution of such a product matrix has been derived, while in [10] the ergodic MI has been characterized for a similar large-matrix multi-hop relay channel. However, in both respective models, the noise is added only at the last hop of the relay channel. In contrast, this paper considers a model where the noise is present at each hop, and therefore the two aforementioned results are not applicable. In the next section, we present the developed framework describing the large-system behavior of a K-hop AF relay MIMO communication system.

III. MAIN RESULTS
Since direct computation of the ergodic MI in (5) is prohibitive due to the necessity of averaging of the logarithmic terms in (6), some simplifying assumptions have to be invoked to make the problem tractable. In the present paper, we investigate the performance of the system in the large-system limit (LSL), meaning that for each intermediate hop k the number of transmit and receive antennas tend to infinity at some constant ratio, viz., ∀i, j M j = α i,j M i → ∞, where α i,j are finite constants. Moreover, we assume that self-averaging property holds, i.e., the randomness of the channel state vanishes in the LSL. This property, being a challenging problem per se, is yet to be proved in the AF relay MIMO context. 6 Hence, the assumption that the property holds is adopted following the existing replica calculus literature [22], [25], [27]. The above assumptions allow us to state the following two results.

A. Free Energy
The main problem in evaluation of (5) is due to the expectation operators over the channel states H in the differential entropy terms (6a) and (6b). The evaluation becomes particularly difficult in the case of non-Gaussian priors. One way to approach this problem is to use the replica method, which allows us to formulate the upcoming Claim 1.
As a preliminary step, the differential entropy in (6a) is rewritten as follows where the partition function of the corresponding many-body system reads with the argument of the exponent being referred to as the Hamiltonian of the system. Then, the normalized differential entropy (6a) has a meaning of the Helmholtz free energy of the system, which can be obtained as [22] Note that in the above expression, the expectation has been moved inside the logarithm and it therefore remains to evaluate the uth moment of Z(y, H) for u ∈ R. Considering the LSL, the free energy (15) is evaluated in the claim below. 7 Claim 1. In the LSL, the free energy (15) is given by (16) at the top of the page, where the parameters are obtained from the solution of the following system of equations Here E z,x |x − x | 2 and E z,x |x − x | 2 denote the MSE and the posterior variance, respectively, associated with the two fixed scalar Gaussian channels given below where w, w ∼ CN (0, 1). In (17h) and (17g), x denotes the MMSE estimate of (18b). In the case of multiple solutions, only that solution minimizing (16) is valid.
Proof: The derivation of the claim is in Appendix A. According to (13), the above result allows the characterization of the differential entropy (6a) of the source under JDD by setting ξ k = η k and ν k = ε k for all k ∈ K. Meanwhile, the conditional entropy term (6b) is found in a similar way by setting η 1 = ε 1 = 0, as will be shown later in Section IV. We also note that a similar result was derived in [48] for a system model that is equivalent to a two-hop channel without noise at the relay. The variables in (17) describe implicitly similar coupled virtual channels to those explicitly obtained in [48].

B. Decoupling
Claim 1 enables the characterization of the performance of an AF relay MIMO system under JDD. To characterize the performance of the SD scheme we formulate the following result.
Claim 2. Let x m , x m and x m q denote the mth entries of x, x and x q . In the LSL, the joint distribution of (x m , x m , x m q ) of channels (1) and (9) converges to the joint distribution (x, x , x q ), associated with channels (18a) and (18b).
Proof: The derivation is presented in Appendix B.
The above claim extends the decoupling principle, reported in [22], to the K-hop AF relay MIMO scenario. It reveals that in the LSL for each data stream m the joint distribution of the input x m , postulated input x m and output of the GPME x m q , associated with the original and postulated AF relay MIMO channels (1) and (9), converges to the joint distribution of the same set of quantities related to the single-user Gaussian scalar channels (18a) and (18b). Namely, the channel with the GPME receiver in the LSL decouples into a bank of scalar Gaussian channels fully characterizing its performance.

A. Achievable Rates
In this section, we apply the results of the asymptotic analysis and evaluate the performance of a K-hop MIMO AF relay channel (1). Firstly, we evaluate the achievable rate under the assumption of full CSI at the receiver. For this, we evaluate the two differential entropy terms, (6a) and (6b), with help of Claim 1.
1) Joint Detection and Decoding: If the JDD scheme is employed at the receiver, the postulated input distribution is set exactly the same as the actual distribution, i.e., q(x) = p(x) and σ = 1. Consequently, we have the following result. Corollary 1. In the LSL, the normalized differential entropy (6a) is given bȳ whereĨ η1 (ρ 1 ) I (z; x) is the MI between the input x and the output z of the fixed scalar channel given by with w ∼ CN (0, 1). The parameters η k , ε k , ∀k ∈ K satisfy the set of fixed-point equations shown in (17) with ξ k = η k and ν k = ε k for all k ∈ K, where the term ε 1 = β 0 ρ 1 E z,x |x − x | 2 reflects the MMSE of the fixed scalar channel (20). Furthermore, the normalized conditional entropy (6b) in the LSL may be computed directly by (19) and (17) by setting η 1 = ε 1 = 0 beforehand. In the case of multiple solutions, the one minimizing the differential entropy of interest is valid.
Proof: The result follows from Claim 1 for the special case of q(x) = p(x) and σ = 1. Due to the RS assumption we set η k = ξ k and ε k = ν k for all k ∈ K. 8 Then, in the conditional entropy term (6b), the inner expectation is taken only over N . Hence one has to proceed in exactly the same way as before, but without replication of x and x , as done in (32) in Appendix A. When the free energy is derived, one sets η k = ξ k and ε k = ν k for all k ∈ K and η 1 = ε 1 = 0 to obtain the asymptotic approximation for (6b).
The expectation in (17h) is taken over the joint distribution p(z, x; η 1 ), and hence ε 1 can be seen as the MMSE of the fixed scalar channel (18a). The entropy termh s (ρ), given in (19), represents the amount of information contributed by the transmitted signal, x, and by the noise components, n k , added at each hop. Note that both the MMSE and MI terms are relatively easy to compute since they are associated with fixed scalar channels (18a) and (18b). Meanwhile, the differential entropyh n (ρ) represents the amount of information discarded at the destination terminal due to noise removal, and hence does not contain terms related to the signal vector x.
2) Joint Detection and Separate Decoding: Based on the decoupling result of Claim 2, along with the fact that for the SD scheme statistical properties of the MIMO channel (1) are completely characterized by the conditional joint distribution p(x m , x m , x m q |H), we have the following result.
Corollary 2. In the LSL, the single-user achievable rate under SD converges tō Knowing both the differential entropy terms in Claim 1 and MI (21) above, one can directly compute the achievable rate as a function of SNRs ρ for both the JDD and SD schemes. For instance, the following examples present the expressions for the MMSE and MI for three particularly relevant signal constellations.
Example 1 (Gaussian inputs). When p(x) is the standard complex Gaussian density, the MMSE term (17h) is given by and the MI between the input and output of (18a) reads Example 2 (QPSK inputs). For the QPSK constellation we have p( and by the I-MMSE relation [50], the MI between the output and the input of (18a) is evaluated as 8 It is noteworthy that p(x) = q(x) and σ = 1 do not necessarily imply η k = ξ k and ε k = ν k . The latter is based on the assumption that RS holds for the JDD scheme, which has been rigorously justified for the CDMA systems with BPSK inputs in [49]. In the unlikely case where the assumption does not hold, there may exist other solutions for which η k = ξ k and ε k = ν k .
Finally, comparing the results obtained for JDD and SD, we are able to quantify the loss due to separation of decoding.

Corollary 3.
In the LSL, the information loss due to separation of detection and decoding is given bȳ where η k and ε k correspond to the entropy term (6a), while η k and ε k correspond to (6b); both obtained according to Claim 1.

B. Bit Error Rate
Since by Claim 2 the performance of the given system with SD is fully described by the equivalent scalar channels (18), we can evaluate the average BER of uncoded transmission. The following result summarizes the finding for a QPSK constellation.
Corollary 4. In the LSL, the average BER of the AF relay MIMO system operating with QPSK inputs is given as where η 1 is a solution to the fixed-point equation system (17).
Example 3 (Individually optimal detector). To realize MAP detection, the receiver needs to postulate the actual distribution, i.e., q(x) = p(x) and σ = 1. This leads to ξ k = η k and ν k = ε k for k ∈ K. Moreover, ε 1 in (17h) depends on the input constellation and has to be evaluated numerically. For instance, for QPSK signals, the parameter is given by (24). In case of multiple solutions, the one minimizing the free energy (16) should be chosen.
Example 4 (Linear detectors). For a linear detector, the postulated input distribution q(x) is set to be standard Gaussian, which yields the MSE of (18a) and posterior variance of (18b) Then, by tuning the parameter σ, we can obtain the three popular linear detectors: • When σ → ∞, the output of GPME converges to the output of the MF. • When σ → 0, the output of GPME tends to that of the ZF detector. • When σ → 1, the output of GPME provides the output of the LMMSE detector. Finally, to obtain an approximation for the average BER, one has to plug the corresponding value of σ into (17) and find the unique solution for the system of fixed-point equations. Then, the resulting η 1 is plugged into (25) and the average BER is acquired.

A. Achievable Rates
To support the main result, we simulate a network with K = 3 hops. To satisfy the power constraint (2), we set the normalization coefficients β k , according to

1) Joint detection and decoding:
For the first simulations setup let us fix M 0 = M 1 = M 2 = M 3 = 8, so that each terminal has eight antennas, and let SNRs be ρ 2 = ρ 3 = 20 dB, while varying the SNR of the first hop ρ 1 . In Fig. 2, we plot the ergodic rates achievable with JDD at the receiver for two types of channel inputs: Gaussian and QPSK signals. The rate loss of 1/K here is due to the TDMA protocol. Solid lines reflect the asymptotic results, while markers represent the numerical averaging via Monte Carlo simulations over at least 1000 channel realizations. 9 We note that for the case of Gaussian inputs the approximation matches the simulations perfectly. Moreover, quite expectedly, it reproduces the result obtained in [12,Sec. 3.5.3]. Meanwhile, for the case of QPSK inputs there is a slight gap in the middle-and high-SNR region, which decreases with increasing the actual number of antennas at each terminal.
To illustrate the latter argument, we plot in Fig. 3 the simulated values of the ergodic rates achievable with QPSK inputs and a JDD receiver vs.  Fig. 2. The star markers 9 Discrete constellations usually require a much higher number of Monte Carlo iterations to obtain stable results. For instance 2.5 · 10 10 iterations have been simulated to produce each point in Fig. 2 and Fig. 3  at 1/M → 0 represent the asymptotic results obtained using Corollary 1 in the LSL and quadratic curves are fitted to the simulated data using non-linear least-squares regression. From the extrapolation one can see that in both cases the simulated ergodic rate approaches the RS solution as M → ∞.
Next, to illustrate the usefulness of the result we incorporate pathloss γ k = d −α k into the SNRs ρ k (set to 10 dB for all k ∈ K) where d k is the distance between terminals k − 1 and k, and α = 4 is the pathloss exponent. In Fig. 4, we plot the achievable rate of JDD as a function of the distance, d, between the source and destination terminals. The three curves correspond to K = 1, 2, 3. The relays are added in such a way that all the terminals of the network are equidistant. Notably, for different values of d, different numbers of hops provide higher achievable ergodic rate (depicted by a dashed line). Thus, one could, in principle, use this information to select the most suitable number of relays or their positions. This selection, however, falls outside the scope of the present paper.
2) Separate decoding: In Fig. 5, for the same antenna setup we plot achievable rates of both JDD and SD schemes as functions of SNR for different signaling schemes (Gaussian, QPSK and 8-PSK). We observe that at low SNR the performance curves of discrete constellations tend to follow the respective curves (JDD and SD), related to Gaussian signaling schemes. Then at certain SNR values, the performance curves, corresponding to the discrete inputs, suddenly switch to the entropy-limited regime. Such a sudden change in performance indicates the occurrence of a phase transition at a certain SNR threshold. A physical analogy to such a behavior is freezing water or the hysteresis of a ferro-magnetic material [51]. From a practical point of view, this can be explained by the sparsity of a discrete constellation, which helps to identify symbols perfectly once the SNR is sufficiently high. Similar behavior was also observed in the context of MIMO [25] and CDMA systems [26], as well as iterative turbo coding [52]. Notably, for the SD scheme with antenna ratios α 0,1 , α 1,2 > 1, Gaussian signaling is no longer optimal, and it is outperformed by discrete constellations, as previously observed in [22], [53].
To quantify the performance of linear detectors, Fig. 6 plots the achievable ergodic rate as a function of SNR for the network with three hops. To enable illustration of the performance of the ZF detector, we set the numbers of antennas unequal. For instance, we set M 0 = 4, M 1 = 6, M 2 = 8 and M 3 = 12. Meanwhile, the SNRs of the hops are set as ρ 1 = ρ, ρ 2 = 0.7ρ and ρ 3 = 0.5ρ. Fig. 6(a) depicts the performance of various detection schemes under Gaussian channel inputs, while Fig. 6(b) illustrates the behavior under QPSK inputs. From the figure one can see that, quite expectedly, the JDD scheme outperforms the other depicted detectors. This is because optimal JDD is essentially the best one can do. In addition, from Fig. 6(a) one notes that for Gaussian signaling performance of the MAP detector matches to that of the LMMSE, which highlights the optimality of the latter in the case of Gaussian signals.
Another observation is that the MF detector demonstrates near-optimal performance at low SNR and becomes increasingly inefficient as the SNR increases. In turn, the ZF detector demonstrates significant degradation in performance at low SNR, whilst becoming progressively efficient with increasing SNR. The same holds for the case of QPSK inputs, as depicted in Fig. 6(b). However, in contrast to Gaussian signaling, where the loss due to separation of decoding (when comparing optimal JDD and MAP detection schemes) grows with SNR, for QPSK constellation this loss vanishes in the high-SNR region due to saturation of the MI in this setup. Moreover, while for the JDD scheme Gaussian signaling always outperforms QPSK, for the MAP detector scheme this is not always the case. Accurate comparison of the two figures reveals that for the small mid-SNR region (roughly between 5 and 10 dB) the QPSK constellation actually performs better than Gaussian signaling under SD. Also, as expected, linear detection schemes (MF, ZF and LMMSE) perform worse.

B. Bit Error Rate
To verify the obtained approximation for BER (25), we simulate the average BER of a two-hop AF relay MIMO communication system with QPSK signaling. We fix the numbers of antennas at terminals M 0 = 24, M 1 = 30 and M 2 = 36 and relations between SNRs at each hop ρ 1 = ρ 2 = ρ. Fig. 7 plots the results of numerical simulations, averaged over 500 channel realizations. From the figure, we note that the asymptotic result predicts well the behavior of the system at low SNR even for a system with finite (albeit quite large) size. However, as SNR increases the approximation becomes less accurate, and henceforth the diversity order effects, usually visible in that region, are not captured. The reason for this is that at high SNR the analysis becomes much more sensitive to the large-system assumption. It becomes increasingly precise as the system size grows large, while the average BER in such regime tends to a Gaussian Q-function. At the same time, for the sake of comparison Fig. 7 also depicts average BER of the optimal MAP detector. At low SNR the corresponding performance curve follows the upper bound (given by the performance of a linear MMSE detector), whereas at high SNR it tends to the lower bound 10 given by 10 The lower bound is attained with interference-free transmission. That is, since in (3)    where Meanwhile, Fig. 8 depicts the average BER of a MAP detector for a similar AF relay MIMO scenario but with K = 3 hops. Here, the numbers of antennas are set to M 0 = 10, M 1 = 9, M 2 = 8 and M 3 = 7, while SNRs is colored, a whitening filter of form is applied, which yields a signal with covariance E H CC H and white noise, where C is given in (29).   of the hops are still equal, i.e., ρ 1 = ρ 2 = ρ 3 = ρ. The figure shows a waterfall-like behavior for the performance of the MAP detector. Namely, the globally stable (true) solution, which minimizes the free energy, instantly switches between the upper and lower bounds at the transition point. At the same time, there exist metastable solutions, which are local minimizers for the free energy. Since the detection algorithm does not know the true detection results initially, the system may get trapped in a metastable solution. When taking into account the latter, a Z-shaped curve is obtained for the BER of the MAP detector (depicted by a dash-dotted line in the figure). Note that there is a region where the curve exhibits an increase in BER with increasing the SNR. The true system BER, however, never gets to that point and instead manifests a sharp transition when reaching the transition point. This behavior is somewhat similar to that of Fig. 5, also indicating the occurrence of a phase transition.
VI. CONCLUSIONS In this paper, we have developed a framework for the asymptotic performance analysis of a K-hop AF relay MIMO system with arbitrary K and discrete channel inputs. The framework captures the effects of separation between detection and decoding, as well as suboptimality of linear detectors. More precisely, we have evaluated the performance of the system under separate detection in terms of achievable ergodic data rate, as well as average bit error rate, in the limit where antenna arrays grow large without bounds. The main result states that the K-hop AF relay MIMO channel with the GPME detector at the destination terminal decouples into a bank of per-stream scalar channel with a GPME detector front end. Comparing to Monte Carlo simulations, it has been confirmed that the results provide an accurate approximation for a finitesized system at low and moderate SNR. It has been further shown via a numerical example that the number of hops in a multi-hop AF network might have a significant impact on the system performance and might be properly adjusted using the present results. Moreover, the individually optimal detection scheme is shown to experience a phase transition at certain SNR values. The obtained compact expressions may also be useful for the design of coding schemes improving the system performance. The presented results are potentially extendable to more sophisticated channel models of interest (e.g., Kronecker model or Rician fading) allowing for further performance optimization as in [54]. Another practical problem of interest would be to consider the case of covariance mismatched decoding [45], where the receiver does not know the instantaneous realizations of the intermediate channels.

APPENDIX A DERIVATION OF CLAIM 1
Recall the expression of the partition function given in (14) (30) The free energy in (15) can thus be rewritten as Since for real-valued u computing E y,H {Z u (y, H)} is very difficult, we make a non-rigorous assumption of replica continuity. That is, it is postulated that the uth moment of the partition function can be first evaluated for integer u and then generalized to real-valued u assuming analytic continuation 11 . We proceed with the so-called replication, i.e., we introduce u replicas of the postulated channel as follows where x (a) and n (a) k denote the ath replica vectors that are assumed to be i.i.d., while x (0) and n ] T ∈ C M k (u+1) . In (32), the averaging should be performed over all possible channel inputs X, channel gains H and noise realizations N . According to the Fubini theorem [56,Theorem 18.3], provided that the expectation in (32) exists, the multiple integral can be computed via repeated integrals. In other words, averaging over the channel matrices and noise vectors of all hops can be done iteratively hop-by-hop. That is, at each hop, we average out the randomness of the corresponding channel matrix and the noise vector, while keeping the variables related to other hops fixed. To do this, define the following set of vectors v (0) 1 For each k, stack these vectors as V k [v k−1 }, a ∈ {1, . . . , u}, vector V k is a complex Gaussian random vector (vide [26]) with the covariance matrix given by where and E 1 01×u 0u×1 σ 2 Iu ∈ R (u+1)×(u+1) . Starting from hop K, we proceed in a similar way as in [22], [26] and write (32) as 1 M 0 ln E y,H {Z u (y, H)} = 1 M 0 ln e M K G (u) K (Q K ) dµ (u) T (u) 1 (Q 1 ,Q 1 ) = α 0,K ln 1 + u σ 2 + α 0,K ln det (I u+1 + Q K Σ) − 1 M 0 ln E X exp ρ 1 β 0 X H (Q 1 ⊗ I M0 )X + K−1 k=1 α 0,k ln det I u+1 − ρ k+1 β kQk+1 (Q k + E) + K k=1 α 0,k−1 tr{Q kQk } + uα 0,K ln πσ 2 (51) T (u) 1 (Q 1 ,Q 1 ) = (u − 1)α 0,K ln 1 + p K − q K σ 2 + α 0,K ln 1 + p K − q K σ 2 + u σ 2 (1 + r K − m K − m * K + q K ) + α 0,K u ln πσ 2 + K k=1 α 0,k−1 (r k r k + um k m k + um * k m * k + up k p k + u(u − 1)q k q k ) α 0,k ln 1 − β k ρ k+1 um * k+1 m k + um k+1 m * k + (p k+1 + (u − 1)q k+1 ) σ 2 + p k + (u − 1)q k +r k+1 (1 + r k ) where we have omitted all the vanishing terms, and the probability measure of Q K is given by which yields the rate function I (u)  Combining all together, we get at hop K K (Q K ) + α 0,K ln π(u + 1) + α 0,K ln det(I u+1 + Q K Σ), (44) where Σ I u+1 − 1 u+1 1 u+1 1 T u+1 .
We now need to evaluate the second term in (44) following the same procedure as above. Namely, we rewrite it, omitting the vanishing terms, as where, using the Gaussian integral, we can obtain G (u) K−1 (Q K−1 ) = α 0,K−1 ln det I u+1 − ρ K β K−1QK Q K−1 + E .
Now, proceeding with the same steps as before, we arrive at 1 M 0 ln M (u) K−1 (Q K−1 ) + α 0,K−1 tr{Q KQK } + α 0,K−2 tr{Q K−1QK−1 } + α 0,K−1 ln det I u+1 − ρ K β K−1QK Q K−1 + E + α 0,K ln π(u + 1) + α 0,K ln det(I u+1 + Q K Σ). (48) Proceeding with the same procedure for all the hops, we can show that for k ∈ {2, . . . , K} the corresponding log-MGF term can be written as where T (u) 1 (Q 1 ,Q 1 ) reads as (51) at the top of the next page. Finding the fixed point of (50) is a complicated task and may not be realizable directly. Hence, a simplifying replicasymmetry (RS) assumption is made in order to proceed. Namely, all k matrices Q k andQ k are postulated to have the following structure The assumption has been widely accepted in the literature following the reasoning that the physics of the whole system should not depend on the artificially introduced replica indices. 12 With the RS assumption, (51) is simplified via the Gaussian linearization (based on the Hubbard-Stratonovich transform [60], [61]), as carried out in [22], to (53) on the top of the page, where η 1