Exploring Uplink Achievable Rate for HPO MIMO Through Quasi-Monte Carlo and Variance Reduction Techniques

The power consumption at the receiver side will be dramatically increased in the millimetre-wave and massive multiple-input-multiple-output (MIMO) communication systems due to the wide bandwidth and a large number of antennas adopted. A half phase-only MIMO (HPO MIMO) scheme, in which the base station (BS) acquires $\pi $ -periodic phase measurements of the complex envelop signals was proposed very recently to overcome the above problem. Due to the non-linear nature of HPO MIMO, the valuation of the achievable rate is very challenging. The purpose of the paper is to provide an efficient method for calculating the achievable rate of the HPO MIMO system. By the mutual information theory, we transform the achievable rate into a sum of two high-dimensional integrations. However, calculating those integrations suffers from the enormous computational burden when using the traditional Monte-Carlo method. In order to increase efficiency, a new method by combining quasi-Monte Carlo with a variance reduction technique is proposed. Besides, we derive the probability density function (PDF) of the HPO MIMO system and analyze the uplink achievable rate of the HPO MIMO scheme. Numerical results show that our proposed method is efficient for calculating the achievable rate of the HPO MIMO system. With the proposed method we confirm that HPO MIMO is a promising technology in future low-power communication scenarios.

the above problems. According to our survey [15]- [17], non-linear receivership (each RF link has the 4.135 Watt) have less power consumption and cost than low-resolution ADCs (each RF link has the 9.62 Watt). This is because non-linear receivers require only one RF link, while lowresolution ADCs take advantage of low power consumption. 1-bit ADC also requires two RF links (traditional receivers require Inphase (I) and Quadrature (Q) signal components channels to modulate signals separately).
Recently, half phase-only MIMO (HPO MIMO) scheme as one of the above nonlinear MIMO systems was first proposed in [12]. This scheme is designed by using lowprecision low-noise-amplifiers (LNAs), π -phase detector, and only one high-resolution ADC on each RF chain at BS side. Because cumbersome superheterodyne circuits and inphase-quadrature (IQ) demodulation have been removed, the circuit power consumption and fabrication cost of the HPO MIMO system could be much lower than that of the conventional MIMO system. Besides, in existing non-linear MIMO systems, due to the fact that the observed information is missing half dimension information, the problem of 'observation ambiguity' is encountered in channel estimation part. However, existing references show that HPO MIMO suffers less ambiguity than envelop modulation systems (also named amplitude or magnitude MIMO) [12], [18]. The ambiguity of HPO-MIMO systems has been solved [13], but the ambiguity of magnitude MIMO systems have not been solved. Therefore, it is more meaningful to study HPO MIMO system. More details about the HPO MIMO can be found in [12].
Equipping π -phase detectors on multiple antennas seems to be a promising technique for the future low power consumption communication systems [19]. Channel estimation (CE) problem in HPO MIMO system has been discussed in [13], [20]. Multiuser detection (MUD) problem in HPO MIMO system has been discussed in [12], [13], [20]. Besides, the synchronization problem in which joint time and frequency in HPO MIMO have been discussed in [13], which exploits the repeatability of a pilot preamble. Despite this growing interest in the HPO MIMO system, to explore the low-power and low-cost system from the theory aspect has received little attention to date. This paper is to fill in this gap.
The purpose of the paper is to provide an efficient method for calculating the achievable rate of the HPO MIMO system. Due to the non-linear nature of HPO MIMO, it has no closed-form solution, then we can not calculate it in the same way as traditional MIMO systems. Therefore, we seek the achievable rate of HPO MIMO from the perspective of mutual information theory. However, when we use mutual information theory, high-dimensional integration is inevitable, especially when the number of N R , N T is relatively large. Therefore, numerical analysis methods play a critical role here. Unfortunately, Monte Carlo (MC) method has outrageous low computational efficiency characteristics, especially in the HPO MIMO system. Therefore, it is essential to find an efficient numerical analysis method to calculate highdimensional integrals. Based on number theory and abstract algebra, we propose a quasi-Monte Carlo (QMC) method to calculate the achievable rate of the HPO MIMO system. Antithetic variables technology as one of variance reduction techniques is used to further increase the efficiency of calculation.
The main contributions of this paper are three-fold. Firstly, the paper proposes an efficient method to calculate the uplink achievable rate for HPO MIMO. In addition, the proposed algorithm can also be used to calculate other problems with high-dimensional integrals. Secondly, this paper derives the PDF of the HPO MIMO system. Thirdly, the paper compares the performance of HPO MIMO and conventional MIMO.
The rest of this paper is organized as follows. Section II presents the details of the proposed RF chain structure and the model of HPO MIMO. Section III derives the uplink achievable rate for HPO MIMO. Section IV develops a new method for calculating achievable rate of HPO MIMO. Section V reports some simulation results which reveal the efficiency of our proposed method and performance of HPO MIMO. Section VI concludes.

II. SYSTEM MODEL
A single-cell HPO MIMO system with N T single-antenna terminals and a BS with N R antennas in Fig. 1 is considered [13], where each antenna is equipped with a HPO-detector (containing π -phase detector and combines with only one high revolution ADC) at the BS side. For the uplink, we have three assumptions: 1. all N T users transmit independent data symbols to the BS simultaneously through the stationary memory-less flat fading channel H ; 2. the sampled signal constitutes a sufficient statistic data at VOLUME 8, 2020 the receiver; 3. in this paper, we assume H is long enough to be estimated, and changes to an independent realization come in the next block. A discrete-time baseband is considered in our model. The signal received at BS after obtaining π -phase operator g HPO (·) is formulated as y = g HPO {(Hs + w)} = g HPO {(z + w)} = g HPO (r), (1) where the function g HPO (·) means the operator of obtaining π -phase of the received signal y as shown in equation (2). We denote Hs by z, and denote Hs + w by r, and s ∼ CN (0, σ 2 s I N T ) is the transmitting signals by each user. It obeys the complex Gaussian distribution with 0 mean and σ 2 s variance. I N T represents the identity matrix with N T dimension. In our simulation, we fix ε (|s| 2 ) = 1. w ∼ CN (0, σ 2 w I N R ) is the N R * 1 additive white Gaussian noise vector with 0 mean and σ 2 w variance. H represents the N R * N T channel matrix which remains constant for coherence period. Therefore, r ∼ CN (0, σ 2 s HH H +σ 2 w I N R ). Due to the character of our HPO-detector, formula (1) can be extended to where r i is the i-th element of the vetctor r, i = 1, 2, . . . , N R . The inverse tangent function atan(·) in (2) places the angle in the correct quadrant. The function π (·) is used to replace the function atan(·) in the following content. Due to the physical characteristics of the HPO detecter, the output signal falls within the interval of (−π/2, π/2].

III. ACHIEVABLE RATE FOR HPO MIMO
We compute the achievable rates instead of capacity since the capacity of the HPO MIMO system is too cumbersome to find. The key characteristic in our HPO MIMO system is that the received signal y at BS side is missing half the dimensions information after HPO-detector g HPO (·) process. Since we are considering a MIMO system, we need to consider the number of receiving antennas when calculating the achievable rate of the system. The probability of the output y conditioned on z plays a crucial role. Since the noise vector w is white across the receive antennas and the nonlinear function g HPO (·) acts on each antenna separately, the conditional distribution of y given z can be factorized as where N R denotes the number of receiving antennas, y i and z i denote the i-th elements in the real vector y and the complex vector z, respectively. r ∼ CN (z, σ 2 w I N R ) when z is given. To compute the achievable rate for HPO MIMO, we first propose the following result Theorem 1: The PDF of π -phase MIMO y i,HPO = π (r i ) is as follows where i = 1, . . . , N R , Proof: According to [9], the PDF of 2π -phase MIMO f (y i |z i ) has the form of (5). After transforming angle, one can easily obtain the result of Theorem 1.
Since we assumed that the channel H is memoryless, the capacity is given by the maximum mutual information between the channel input s and the detector output y. Furthermore, since H is known to the receiver, the channel output is the pair (y, H). The mutual information can be written as From (1), we obtain that s and H are statistically independent. We can observe that s, z and y constitute a Markov chain s → z → y under a given channel realization H = H. Because s and y are conditionally independent, (6) can be upgraded to where according to the information theory, the mutual information of interest I (z; y|H ) is extended to From (7) to (8), the fact of h(y|z, H) = h(y|z) is used. In order to calculate the mutual information of interest I (z; y|H), the PDFs f (y|H) and f (y|z) are needed. Given f (y|z), f (y|H) can by rewritten as Combining (8) and (9), we have where f (y|z) is computed by (3) and (4). Remark: Deriving the expression of f (y|H), f (z, y) and f (z|H) is nontrivial. Instead of trying to explore the analytical expression of them, we evaluate (8) by the idea of stochastic simulation sample mean algorithm such as Monte Carlo (MC).

MC approximate h(y|H), f (y|z, H) and f (y|H) by
respectively, where z j as the MC samples are generated according to complex Gaussian distribution CN (0, σ 2 s HH H ). Combining (3), (10), (11) and (13). We approximate the achievable rate for HPO MIMO as Finally, averaging (11) -(13) over many different channel response H , we can obtain the mutual information in (6).

IV. ANTITHETIC-QMC METHOD FOR CALCULATING ACHIEVABLE RATE OF HPO MIMO
In this section, we will answer why we need to use QMC sampling to calculate multi-dimensional integrals instead of traditional MC sampling methods. Also we give the details of QMC and antithetic variates method to calculate the mutual information of interest I (z; y|H) in HPO MIMO system. As can be seen from formula (8), to calculate the achievable rate of the HPO MIMO system, we have to make two sampling approximations for the first term and one sampling approximation for the second term. Also, a subtraction operation is made after the above sampling approximations. As we knew, the MC sampling principle is to achieve an infinitely close to multi-dimensional integral by increasing the sampling point density, which means the more numbers of sampling points, the higher the accuracy of the multi-dimensional integral approximation. In our simulation, to obtain a stable achievable rate, the number of sampling points is at least 10,000. Therefore, our simulation requires at least 10,000*10,000 times to obtain smooth convergence. The above summarizes only under the case of the single input single output (SISO) situation. If the situation is extended to a multi-dimensional MIMO scenario, the number of simulation will increase drastically as the growth of antenna array numbers. This will make computational efficiency very low. Therefore, to increase the efficiency of multi-dimensional integration calculation becomes very necessary and essential in the HPO MIMO system.

A. MONTE CARLO
We consider MC method to calculate an integral problem, such as formula (12), (11) and (13) in Section II. We assume f (·) is an integral function on the unit cube C q = [0, 1) q with q-dimension. In order to calculate integral I (f ) as follows [21] we first extract random samples Then we calculate an estimateÎ K of I (f ) by the following formula: The key point of the MC method is that it is necessary to extract an independent, uniformly distributed random sequence in C q [22]. Since a certain algorithm generates such a random sequence, it also can be called the pseudorandom sequence. According to the law of large numbers,Î K converges to I (f ) in probability. Let where Var denotes variance. By the central limit theorem, the integration error produced by the MC satisfies the following inequality with probability approximate 1 − α, where z α/2 denotes the α/2 quantile of Gaussian distribution N (0, 1). Therefore, the convergency rate of MC method is O(1/ √ K ).

B. QUASI-MONTE CARLO
Different from MC, QMC methods draw on number theory and abstract algebra [23]. The basic idea of QMC methods is to replace the pseudo-random sequences in MC with low discrepancy sequences (LDS) or quasi-random sequences which are more uniform than random sequences. To measure the uniformity of a sequence, we introduce the following definition [24], [25] Definition 1: Let B be a collection of all rectangles in where #{x i ∈ B} denotes the number of x i contained in B.
According to the law of iterated logarithms, if a sequence is random, the expectation of the discrepancy of the sequence is bounded by log log K / √ K . In contrast, the discrepancy of many low discrepancy sequences such as Halton sequence is bounded by constant times (log K ) q /K .

1) HALTON' LDS
Any non-negative integer k can be constituted as the following form based on the prime number b: where We define the radical inverse function ϕ b (k) based on b as follows: it's easy to find for any integer k 0, In fact, it is no need to generate the Halton sequence starting from k = 0. That is to say, for d nonnegative integers n 1 , . . . , n d , the above sequence can be taken as We can generate Halton' LDS by the following algorithm

Algorithm 1 An Algorithm for Generating Halton' LDS
Niederreiter [26] states that the discrepancy of Halton sequence satisfies that where c q is a constant which only depends on dimension q. (18) implies that the Halton sequence is significantly more uniform than a random sequence.

2) ERROR BOUND OF QMC
The integration error produced by QMC satisfies the following Koksma-Halwka inequality.
where D K is the discrepancy of the sequence P K . The variation V (f ) in (22) of f (x) in the sense of Hardy and Krause is defined as: If f (x 1 , . . . , x s ) is sufficiently differentiable, for all positive k ≤ s and k integers 1 ≤ i 1 < i 2 < . . . < i k ≤ s, define the quantity as follows: (23) and (24) imply that once f (· ) is given, V (f ) would be a constant. Therefore, the integration error produced by QMC only depends on D K . From (18), the convergency rate of QMC method based on Halton sequence is O (log K ) q−1 K . Since (log K ) q−1 can be absorbed into any power of K , the convergency rate of QMC method can be thought as near O(1/K ). Therefore, QMC methods accelerate convergence from O(1/ √ K ) of MC methods to nearly O(1/K ). Statisticians focus on how to construct the best low discrepancy (quasi-random) point sets [27]- [30]. At the same time, variance reduction techniques [31], [32] are widely studied from another side for improving the efficiency of such sampling methods. In our work, we combine the flexibility of variance reduction techniques with the characteristic of effectiveness and fast convergence of low discrepancy sequences. We not only use the Halton sequences as our low discrepancy point sets, but also combine antithetic variates technique as the variance reduction techniques, and the detail can be found in the following part.

C. ANTITHETIC VARIATES TECHNIQUE
Antithetic variates is one of variance reduction techniques [33] for MC sample [34], [35]. It attempts to reduce variance by introducing negative dependence between pairs of sampling points [36]. If U is uniformly distributed over [0, 1], then 1 − U obeys uniform distribution, too. According to this principle, if we generate U 1 , . . . , U n as the first path, then we still can generate 1 − U 1 , . . . , 1 − U n as the second path without changing the law of this simulation process. Antithetic variates technique reduce the variance based on the following fact where Cov[·, ·] denotes the covariance. In particular, for Gaussian distributions commonly used in telecommunications, antithetic variates can be implemented by pairing the sequences Y 1 , Y 2 , . . . of independently and identically distributed (i.i.d.) CN (0, 1) variables. The key characteristics of the antithetic variates technologies is that for each i, Y i andỸ i have the same distribution, and all the pairs (Y 1 ,Ỹ 1 ), (Y 2 ,Ỹ 2 ), . . . , (Y n ,Ỹ n ) are i.i.d. The antithetic variates estimatorŶ AV is the average of all 2n observations from the common distribution of Y i andỸ i , the variance can be written as 75878 VOLUME 8, 2020 because Y i andỸ i have the same distribution, the variance of Y i andỸ i are the same. The condition for antithetic sampling to reduce variance becomes therefore, (26) implies that the efficiency of the MC method is increased after introducing an antithetic variate. This technique can also be applied in the QMC method.

D. ANTITHETIC-QMC METHOD
We combine antithetic variates technique and QMC method and propose the Antithetic-QMC method to calculate uplink achievable rate for HPO MIMO systems. We develop Antithetic-QMC algorithm as follows:

Algorithm 2 Antithetic-QMC Algorithm for Calculating the Uplink Achievable Rate for HPO MIMO System
Input:

V. NUMERICAL RESULTS
In this section, we present some numerical experiments to confirm the analysis in the previous parts. To verify the correctness of the PDF of HPO MIMO system which we derived, we first compare the histogram of HPO MIMO distribution and the PDF of HPO MIMO. In the following numerical simulations, the input s of transmitter follows by Complex Gaussian distribution s ∼ CN (0, σ 2 s I N T ). The channel matrix entries H are modeled as i.i.d. zero-mean unitvariance circularly symmetric Gaussian random variables. Where we assume σ 2 s = 1 and σ 2 h = 1. The received vector z is perturbed by a zero-mean circularly symmetric Gaussian noise vector w, with autocorrelation function E[w k w l ] = σ 2 w I N R δ(k − l), where k and l are symbol instants, and we We assume all users have the same level of SNR in our experiments. Fig. 2 (a) shows our derived 3D joint PDF f (y i , z i ) of HPO MIMO, imaginary part and real part of z i are used to illustrate f (y i , z i ). In order to describe our derived PDF f (y i , z i ) of HPO MIMO more clearly, we plot the corresponding 2D histogram of PDF y π −phase under the different variance σ 2 w of Gaussian noise in Fig. 2 (b), (c) and (d), respectively. We randomly generate the HPO MIMO distributions y HPO according to the phase-obtaining operation g π {(Hs + w)}, and draw their normalized histograms as shown in the blue columns. At the same time, we plot the HPO MIMO PDF f HPO (y i |z i ) as shown in the red line which we derived from (3). It can be seen that these blue columns and the red line are closely matched under the different noise σ 2 w . Those results confirm the correctness of our derived PDF f HPO (y i |z i ) in (3).
To test the uniformity of Halton' LDS, we implement algorithm 1 and compare distributions of N 1 = 4 (Fig. 3) and N 2 = 16 (Fig. 4) points on a unit square n = 2 given by two different sampling techniques: MC and Halton' LDS. This provides a qualitative picture of the uniformity properties of these sampling techniques. In the first case, the unit square is divided into 4 and 4 2 squares of measure 1/4 and 1/4 2 , respectively. In the second case, the unit square is divided into 16 and 16 2 squares of measure 1/16 and 1/16 2 , respectively. Fig. 3/ Fig. 4 (c), (d) show 2-dimensional distributions of Halton' points, and each of the 4 small squares contains exactly one Halton' point. Random sampling MC (Fig. 3/  Fig. 4 (a), (b)) do not possess either of these properties. Fig. 5 shows distributions of 64 points in two dimensions. From Fig. 5 (c), it is clear that each of the 64 2 subsquares  contains exactly 1 Halton' point ( Fig. 5 (d)). And this phenomenon is existing in all types of MC (Fig. 5 (a)) samplings: clustering and empty subsquares are clearly visible from these plots (Fig. 5 (b)). From Fig. 3 to Fig. 5, in summary it  would appear that Halton' LDS sampling gives a better way of arranging N points in n dimensions than MC method. Furthermore, we extracted 2000 points from the two-dimensional unit square using the MC method and the Halton' LDS, respectively. From Fig. 6, we can see that the points using the Halton' LDS method are more uniformly, whereas the points using the MC sampling method are denser in some places and sparse in others. Fig. 6 shows that Halton' LDS has smaller standard error than MC.
To test the efficiency of the Antithetic-QMC method, we implement algorithm 2 and compare Antithetic-QMC with MC and QMC. MC and QMC for computing achievable rate of HPO MIMO can be implemented by modifying algorithm 2. The MC estimateÎ converges to the true value of I as sample points N approaches to infinity by the law of large numbers. Therefore we calculate the accurate value of I by MC method. A comparison of methods requires a figure of merit. For sampling-based methods, standard error is an appropriate figure of merit [37]. As our figure of merit, we take log 2 of standard error δ, where The smaller log 2δ, the higher efficiency of calculating.  Fig. 7 compares the log 2 of standard error among MC, QMC and Antithetic-QMC for calculating the achievable rate for HPO MIMO systems with simulation trials N = 100, 1000, 5000 and 10000, respectively. The reason for the choice of 10000 as the limit is that in the variance reduction technique, the standard error δ = 0.02 marks the result of the sampling method close to the multi-dimensional integral result. Fig. 7 shows that QMC converges fast, and the standard error is significantly smaller than MC. Antithetic-QMC further increases the convergency rate compared to QMC.
The performance of Antithetic-QMC is slightly better than QMC at a high number of samples. However, these two methods have much better performance than traditional MC methods, especially at higher sampling numbers. This is precise because both of these QMC methods use a more uniform LDS sequence instead of a random sequence in calculating the multi-dimensional integration in the HPO MIMO system. We can conclude that Antithetic-QMC is superior to QMC and MC in calculating uplink achievable rate for HPO MIMO system.
With the proposed method we compare the achievable rate between the HPO MIMO and the conventional MIMO under the SISO and MIMO scenario. Fig. 8 reports the main outcomes. Fig. 8 presents two key phenomena. Firstly, with the increase of the signal-to-noise ratio (SNR) from −10dB to 30dB, achievable rate keeps increasing. Secondly, HPO MIMO always achieves more than half the performance of its corresponding conventional MIMO. Although the performance of HPO MIMO is inferior to that of conventional MIMO, we do prefer π -phase detection. This is because cumbersome demodulation is removed, and no analog combination of IQ signals is required as conventional MIMO, its circuit power and costs are much lower than its corresponding MIMO [11]. This result indirectly indicates that HPO MIMO can achieve the rates of conventional MIMO by increasing the number of antennas at the receiver. The results show that the HPO MIMO system can obtain 50% − 67% of the capacity of traditional MIMO with only little power consumption compare to the conventional MIMO system. This conclusion proves that HPO MIMO is a promising technology in future low-power communication scenarios.

VI. CONCLUSION AND FUTURE
We derive a semi-analytical expression of the achievable rate for HPO MIMO system, in which high-dimensional integrals are approximated numerically by quasi-Monte Carlo and antithetic variates technique. Our proposed method speeds up the computational efficiency of the achievable rate for the HPO MIMO system compared to Monte Carlo and quasi-Monte Carlo methods. The high accuracy can be obtained by increasing sampling points numbers. Therefore, our method can be used as a benchmark for testing efficiency of other analytical approximation or numerical methods to calculate the achievable rate.
With our proposed method, we provide a comparison on the achievable rates of HPO MIMO and conventional MIMO systems. Our comparison results show that although the performance of HPO MIMO is not as good as its corresponding linear detection (which can be achieved more than half), these low-power consumption techniques unexpectedly utilize additional receive antennas to achieve higher spatial multiplexing gain. This result seems promising since it could save more power and cost to use more HPO-receiver RF chains than the conventional ones. He is also the Leader and a coauthor of the comprehensive book Design and Implementation of Information-Centric Networking (Cambridge University Press, 2020). He was involved in many standardization activities organized by ITU-T and ICNRG of IRTF and contributed to the ITU-T Standards ITU-T Y.3071: Data Aware Networking (Information Centric Networking) Requirements and Capabilities and Y.3033-Data Aware Networking-Scenarios and Use Cases. His research interests include smart grids, information-centric networking, the Internet of Things, blockchain, and information security. He has had experience with editorial and conference organizations. Moreover, he has served as a TPC Member for the ITU Kaleidoscope 2020, the IEEE VTC2020-Spring, the IEEE CCNC 2020, the IEEE WCNC 2020, He is currently a Professor and the Head of the Discipline of Network and Cybersecurity, University of Technology Sydney. He is also the Co-Founder and a CTO of Ultimo Digital Technologies Pty Ltd., developing the Internet of Things (IoT) and blockchain. Prior to that, he was a Principal Scientist and the Research Leader with CSIRO, where he led wireless networking research activities. He specializes in system design and modeling and has delivered networking solutions to a number of government agencies and industry customers. He has over 150 research publications and has supervised over 30 Ph.D. students. His research interests include wireless networking, cybersecurity, and blockchain.
Dr. Liu was a recipient of the Australian Engineering Innovation Award and the CSIRO Chairman Medal. He was the Founding Chair of the IEEE NSW VTS Chapter. He has served as the Technical Program Committee Chair and the Organizing Committee Chair in a number of IEEE conferences. VOLUME 8, 2020