Hamiltonian Monte Carlo for Vector Perturbation Precoding in MU-MIMO via Continuous Relaxation

Hamiltonian Monte Carlo for Vector Perturbation Precoding in MU-MIMO via Continuous Relaxation Junichiro Hagiwara∗ , Toshihiko Nishimura† , Yasutaka Ogawa† , and Takeo Ohgane† ∗ School of Social Informatics, Mukogawa Women’s University, Nishinomiya, Japan

Email: hagiwara junichiro [email protected] † Faculty of Information Science and Technology, Hokkaido University, Sapporo, Japan

arXiv:2609.23011v1 [cs.IT] 19 Sep 2026

Email: {nishim, ogawa, ohgane}@ist.hokudai.ac.jp

Abstract—Multi-user multiple-input multiple-output (MUMIMO) is a key technology that improves wireless capacity through multiple antennas. In MU-MIMO downlink precoding, vector perturbation (VP) is a representative nonlinear method that achieves high performance. However, its search for the integer perturbation vector reduces to a closest vector problem, whose complexity grows rapidly as the number of users increases. We propose a method that relaxes the discrete structure of the integer perturbation into a continuous mixture of t-distributions, enabling efficient search via gradient-based Hamiltonian Monte Carlo (HMC). Complexity analysis and numerical experiments demonstrate the effectiveness of the proposed method. Its search complexity scales as O(N 2 ) in the number of users N . At a symbol error rate of 10−3 , it performs within 2.4 dB of a hypersphere approximation benchmark, which approximates the performance limit of VP. This paper reframes the VP perturbation search as a probabilistic inference problem, providing a general formulation for handling high-dimensional discrete search in a continuous space. Index Terms—continuous relaxation, Hamiltonian Monte Carlo (HMC), multi-user multiple-input multiple-output (MUMIMO), precoding, vector perturbation (VP).

I. I NTRODUCTION HE growing use of wireless communications, exemplified by 5G, calls for further improvements in spectral efficiency. Multi-user multiple-input multiple-output (MUMIMO) [1] is a promising multi-antenna technology that improves cell capacity without heavily relying on the implementation capabilities of user terminals. In the MUMIMO downlink, each user terminal cannot directly observe the channel states of other users. Transmit precoding [2], which suppresses inter-user interference in advance at the base station, is therefore essential. Precoding methods are broadly classified into linear and nonlinear methods. Typical linear methods include zero forcing (ZF) and minimum mean squared error (MMSE) [3], while typical nonlinear methods include Tomlinson–Harashima precoding (THP) [4], [5] and vector perturbation (VP) [6], [7]. In general, there is a tradeoff: linear methods have low complexity but limited performance, whereas nonlinear methods achieve high performance at the cost of high complexity. Although current deployments rely mainly on linear methods [8], [9], developing nonlinear methods that achieve high performance with practical complexity is important for further capacity gains. The theoretical performance limit of nonlinear precoding is given by dirty

T

paper coding (DPC) [10], and VP is recognized as a practical approximation of it. In VP, an integer perturbation vector is added to the transmit symbols, and its effect is removed by a modulo operation at the receiver. However, an exhaustive search for the integer perturbation vector requires a search cost that grows exponentially with the number of users, becoming infeasible in large-scale systems such as massive MIMO [11]. This perturbation search problem reduces to a closest vector problem, and various methods have been studied to solve it efficiently [2]. For example, sphere search [7] achieves good performance by restricting the search candidates to lattice points within a hypersphere of a given radius, but it still incurs high computational complexity. More recently, several methods have been proposed, including sequential perturbation updates based on a correlation metric [12], approximate message passing (AMP) applied to ZF and ZF with successive interference cancellation (SIC) [13], and quantum annealing that assumes dedicated hardware [14], [15]. These methods reportedly achieve good performance with high computational efficiency. However, to the best of our knowledge, most existing methods treat the perturbation search essentially as a discrete optimization problem, and the direction of exploiting powerful algorithms developed for continuous spaces has not been sufficiently explored. In this paper, we reframe the VP perturbation search as a probabilistic inference problem and develop a formulation for handling high-dimensional discrete search in a continuous space. Specifically, we propose a method that relaxes the discrete structure of each integer perturbation element into a continuous mixture of t-distributions, enabling efficient search via gradient-based Hamiltonian Monte Carlo (HMC) [16]. This method extends our preliminary study [17]. In particular, we introduce lattice reduction (LR) [18] as preprocessing. LR improves the orthogonality of the lattice basis and thereby the conditioning of the perturbation space explored by HMC. In related work, we have applied a similar probabilistic inference framework to MIMO signal detection at the receiver [19]. In contrast, this paper addresses a different problem, namely the VP perturbation search at the transmitter. It provides a VPspecific formulation, including the construction of a likelihood based on the perturbation search criterion and the introduction

of LR preprocessing. Returning to our proposed method for VP, complexity analysis shows that its search scales as O(N 2 ) in the number of users N . Numerical experiments show that, at a symbol error rate (SER) of 10−3 , it performs within 2.4 dB of a hypersphere approximation benchmark [20]. This suggests that near-limit VP performance can be achieved in polynomial complexity on general-purpose computing platforms. These results are expected to contribute to research and development toward 6G, as well as to broader insights into the continuous relaxation of combinatorial optimization. The main contribution of this paper lies not in a new HMC algorithm, but in formulating the VP integer perturbation search as a differentiable posterior that is tractable with standard HMC. The rest of this paper is organized as follows. Section II formulates the problem, Section III describes the proposed method, Section IV presents numerical results and discussion, and Section V concludes the paper. Notation: 0 and I denote the zero vector and the identity matrix, respectively. AH and AT denote the Hermitian transpose and the transpose of a matrix A, respectively. N (µ, Σ) denotes the normal distribution with mean vector µ and covariance matrix Σ. tν (µ, σ) denotes the t-distribution with location parameter µ, scale parameter σ, and degrees of freedom ν. ∥·∥ denotes the Euclidean norm. Q(·) denotes the nearest-integer quantizer. II. P ROBLEM FORMULATION A. System Model We consider a MU-MIMO downlink system in which a base station with M transmit antennas serves N singleantenna user terminals. The information symbol vector u ∈ CN has elements drawn independently and uniformly from the specified quadrature amplitude modulation (QAM) constellation. Let ũ ∈ CM the precoded transmit symbol vector, H ∈ CN ×M the channel matrix, y ∈ CN the received symbol vector, and w ∈ CN the noise vector. Their relationship is given by y = H ũ + w,

(1)

where w follows the complex normal distribution CN (0, 2 σw I). The total transmit power of the precoded vector is normalized to a given constant Pt , i.e., ∥ũ∥2 = Pt . In this paper, we assume that H is known at the transmitter. B. VP [6], [7] Fig. 1 shows the processing diagram of VP. The precoder takes u as input and produces ũ, while the detector takes y as input and produces û, the estimate of u. The precoder adds an integer perturbation vector to the information vector, and the detector removes the effect of the perturbation by a modulo operation. The perturbation vector is chosen as a whole to improve the detection performance. For symbols adversely affected by H, the perturbation increases their transmit power to mitigate this effect. Conversely, for symbols not adversely affected by H, no perturbation is added. Since the average

Precoder

Detector

Perturbation vector search

Transmitter

Channel

Receiver

Fig. 1. Processing diagram of VP.

transmit power is normalized to a constant, this effectively redistributes the transmit power to compensate for the adverse effect of H. ũ is given by q Pt (2) ũ = ∥G(u+τ ℓ)∥2 G(u + τ ℓ), where ℓ is the N -dimensional Gaussian integer perturbation vector, G = H H (HH H )−1 is the pseudo-inverse of H, and τ = 2cmax + ∆. Here, cmax is the absolute value of the constellation symbol(s) with largest magnitude, and ∆ is the spacing between constellation points. Similarly, û is given by q ℓ)∥2 û = modτ (u + τ ℓ + ∥G(u+τ w), (3) Pt where modτ (·) denotes the modulo operation with modulus τ . In this paper, for simplicity, we assume that ∥G(u + τ ℓ)∥2 is known at the receiver. From (3), the coefficient of w causes noise enhancement, and a smaller value leads to more accurate recovery of u. Therefore, the search criterion for ℓ at the precoder is given by argminℓ ∥G(u + τ ℓ)∥2 .

(4)

2

Hereafter, ∥G(u + τ ℓ)∥ is referred to as the VP cost. Applying lattice reduction [18] as preprocessing to improve the conditioning of G, (4) becomes argminℓ ∥GT T −1 (u + τ ℓ)∥2 = argminℓLR ∥GLR (uLR + τ ℓLR )∥2 ,

(5)

where GLR = GT , uLR = T −1 u, and ℓLR = T −1 ℓ, and T is an N × N unimodular matrix, i.e., det(T ) = ±1. Lattice reduction transforms the basis of G into that of GLR , whose basis vectors are shorter and closer to orthogonal. Therefore, searching for ℓLR in (5) and then transforming it back to ℓ improves the search performance, compared with directly searching for ℓ in (4). We adopt this approach in our method. To obtain T , we use the well-known LLL algorithm [21]. III. P ROPOSED M ETHOD Minimizing the VP cost in (5) is a high-dimensional combinatorial optimization problem, whose solution is not easily obtained. To address this, we take an efficient approach that finds the solution through a probabilistic interpretation. Finding ℓLR for a given uLR is probabilistically equivalent to inferring the posterior distribution p(ℓLR | uLR ). From Bayes’ theorem, p(ℓLR | uLR ) ∝ p(uLR | ℓLR )p(ℓLR ). We will examine each component in turn. In this section, to simplify

Algorithm 1 HMC sampling 1: Initialize ℓLR at random 2: for k = 1, . . . , LHMC do 3: Draw r from N (0, I) 4: Numerically solve Hamilton’s equations to obtain r ′ and ℓ′LR 5: Update ℓLR ← ℓ′LR with probability min[1, exp(HHMC (r, ℓLR ) − HHMC (r ′ , ℓ′LR ))] 6: Regard the updated ℓLR as a sample from the posterior distribution p(ℓLR | uLR ) 7: end for

(a) Discrete categorical distribution

(b) Continuous mixture of t-distributions Fig. 2. Three-peaked categorical distribution and mixture of t-distributions.

the representation of the perturbation elements, we use a realvalued representation equivalent to the complex one. Here, the symbol notation remains unchanged, but the dimensions N and M in the complex representation are regarded as 2N and 2M , respectively.

c

1 2 2 ∥GLR (uLR + τ ℓLR )∥ σlik 2N 1 X X −(ν+1)/2 + log an,i ,

log p(ℓLR | uLR ) = −

n=1

an,i = 1 +

A. Likelihood: p(uLR | ℓLR ) Since the VP cost has the form of a squared Euclidean norm, it can be probabilistically interpreted as a likelihood following a normal density, as follows: 2 −1 p(uLR | ℓLR ) = N (uLR ; −τ ℓLR , σlik (GT ). LR GLR )

the HMC described later operates on differences of the logposterior distribution, we define the log-posterior distribution, retaining only the terms that depend on ℓLR , as follows:

(6)

Although σlik equals 1 by the definition of the VP cost, it affects the numerical search efficiency of the posterior distribution. We therefore retain it as a tunable parameter and set an appropriate value in advance according to the problem. B. Prior Distribution: p(ℓLR ) It is known from [13] that most effective ℓLR concentrate around the reference point Q(−uLR /τ ). Considering the search efficiency, we restrict the deviation from the reference point to {−1, 0, 1} in this paper. The prior distribution therefore follows a three-peaked categorical distribution (Fig. 2 (a)). We relax it into a three-peaked continuous mixture of t-distributions (Fig. 2 (b)), as follows: Q2N p(ℓLR ) = n=1 p(ℓLR,n ), P1 p(ℓLR,n ) = 31 i=−1 tν (ℓLR,n ; Q(−uLR,n /τ ) + i, σpri ) . (7) ν and σpri are tunable parameters, whose appropriate values must be determined in advance according to the problem. We use a t-distribution rather than the typical normal distribution as the component distribution to improve the numerical search efficiency of the posterior distribution. Since the t-distribution has a sharper mode and heavier tails than the normal distribution, it enables both intensive search around a mode and active movement across modes [19]. C. Posterior Distribution: p(ℓLR | uLR ) Since both the likelihood and the prior distribution are continuous, the posterior distribution is also continuous. Since

(8)

i=−1

(ℓLR,n − (Q(−uLR,n /τ ) + i))2 . 2 νσpri

Because the mixture prior is nonconjugate to the likelihood, an analytical point estimate of ℓLR is not available. We therefore explore the posterior numerically. In this case, powerful algorithms developed for continuous spaces can be applied. In this paper, we apply efficient HMC [16], a type of Markov chain Monte Carlo (MCMC), to approximate the posterior distribution with a fixed number of samples. The obtained samples are only candidates, and the perturbation search requires a point estimate. Therefore, after quantizing them to the nearest integers, we finally select, as the final estimate, the one that minimizes the VP cost. Algorithm 1 outlines the HMC sampling procedure. Here, LHMC is the number of HMC iterations, and r is an artificially introduced auxiliary momentum variable that mitigates the risk of the ℓLR search falling into a local optimum. Letting t denote a fictitious time, Hamilton’s equations are dℓLR /dt = r and dr/dt = ∂ log p(ℓLR | uLR )/∂ℓLR . This gradient information of the log-posterior distribution improves the efficiency of the ℓLR search. The Hamiltonian HHMC (r, ℓLR ) is defined as the sum of the kinetic energy 12 ∥r∥2 and the potential energy − log p(ℓLR | uLR ). The HMC procedure itself, including the numerical integration and step-size/pathlength settings, follows the standard implementation based on the Stan language [22]. Note that our contribution lies not in the HMC algorithm but in the differentiable posterior for the VP perturbation vector. IV. N UMERICAL RESULTS AND DISCUSSION The purpose of this paper is to demonstrate the effectiveness of the proposed method based on the continuous relaxation approach. To this end, we evaluate the SER through computer simulations and analyze the computational complexity.

TABLE II PARAMETER S ETTINGS FOR THE P ROPOSED M ETHOD Parameter Likelihood scale σlik Prior scale σpri Degrees of freedom ν Number of parallel HMC chains Number of HMC iterations LHMC

Setting√ 0.72/ N 2.3/N 1.8 5 400

A. Common Assumptions In our simulations, the settings in Table I are held constant across all SER evaluations. Owing to space constraints, we report representative results for 16QAM; in preliminary experiments, we confirmed similar results with the same settings for other modulation orders such as 4QAM and 64QAM. B. Parameters of the Proposed Method (Table II) The tunable parameters of the proposed method were determined by minimizing the SER in a preliminary validation. When comparable SERs were obtained, we adopted the setting with lower complexity. Since the VP search criterion does not depend on noise, the settings in Table II were used commonly across all signal-to-noise ratio (SNR) points. The likelihood and prior scale parameters σlik and σpri control how sharply the gradient of the log-posterior distribution is reflected in the HMC search, and thus strongly affect the search efficiency. In this paper, we confirmed through preliminary validation that their optimal values exhibit a loglog linear trend with respect to the number of users N , and used the empirical formulas shown in Table II. The degrees of freedom of the mixture of t-distributions were searched over the grid {0.5, 0.6, . . . , 4.9, 5.0} and set to 1.8. To reduce the dependence on initial values and improve the stability of the search, we run multiple HMC chains in parallel. The number of parallel chains was searched over the grid {1, 2, . . . , 9, 10, 20, 50} and set to 5. Finally, the number of HMC iterations LHMC was set to 400, considering the balance between search performance and computational complexity. We verify its validity in the following. A larger LHMC improves the search performance but increases the computational complexity, so their balance must be considered. Fig. 3 shows the evolution of the minimum log VP cost achieved up to each HMC iteration. The values in the figure are the medians over 600 trials, obtained under the conditions in Tables I and II (except for Number of SER evaluation trials and LHMC ). The saturation level increases

log VP cost

Item Setting Number of SER evaluation trials 10,000 for N ∈ {8, 16}; 5,000 for N = 32 Transmit antennas and users M = N , with N ∈ {8, 16, 32} Modulation 16QAM Total transmit power Pt = 1 Fading model Uncorrelated quasi-static Rayleigh fading Channel coding None

0.5 1.0 1.5 2.0 2.5

TABLE I C OMMON SIMULATION SETTINGS

N=32 N=16 N= 8

0

200

400

600

800

1000

HMC Iterations Fig. 3. Minimum log VP cost achieved up to each HMC iteration (median over 600 trials).

with N because the VP cost is an unnormalized quantity that includes the dimensional scale. This does not affect the relative saturation assessment of the improvement, where the value at 1000 iterations is regarded as the convergence reference. In Fig. 3, the additional improvement beyond 400 iterations was 0.3%, 1.3%, and 3.3% of the total improvement up to 1000 iterations for N = 8, 16, 32, respectively (reflecting the greater search difficulty at larger N ). We therefore judged that the update of the minimum has largely saturated at 400 iterations, and stopped at 400 iterations as a sufficient value considering the computational complexity. For completeness, convergence of the Markov chains was verified using standard diagnostics (e.g., the Gelman–Rubin statistic [23]); details are omitted as they are not central to this paper. C. Comparison Methods As a representative basic method, we begin with linear ZF, which does not apply VP. As representative statistical VP methods, based on [13], we consider four methods: LRaided ZF, LR-aided ZF-SIC, LR-aided ZF with AMP, and LRaided ZF-SIC with AMP, hereafter referred to as LR-ZF, LRZF-SIC, LR-ZF-AMP, and LR-ZF-SIC-AMP, respectively. For AMP, we adopt the better-performing setting reported in [13]. As a performance benchmark for VP, we consider the characteristic based on the hypersphere approximation (HA) [20], which replaces the residual region of the closest vector problem solved by exhaustive-search VP with a hypersphere of the same volume. This is hereafter referred to as HA. Although HA is not a strict lower bound, it is known to accurately approximate the performance limit of VP [20], and we use it as a reference for the achievable performance in this paper. Specifically, HA is evaluated by a large-system asymptotic expression with a fixed ratio of transmit to receive antennas, and is therefore shown as a constant SNR shift from the SISO AWGN (single-input single-output additive white Gaussian noise) characteristic. Since the proposed method combines HMC with LR, it is hereafter referred to as HMC w/ LR. For comparison, we also consider two variants without LR, obtained by setting T = I: one with the σ settings in Table II, and one with σlik = σpri = 1, unlike in Table II. These are hereafter referred to as HMC w/o LR and HMC w/o LR (σ = 1), respectively. Note that quantum annealing

(a) M = N = 8

(b) M = N = 16

(c) M = N = 32

Fig. 4. Average SER versus average received SNR.

relies on dedicated hardware and different assumptions from this paper, so it is excluded from the comparison. D. SER Performance (Fig. 4) The performance of linear ZF degrades as N increases. This is because, as N increases, small channel eigenvalues are more likely to appear, which in turn increases the risk of noise enhancement [24]. Among the methods proposed in [13], the performance improves in the order LR-ZF, LR-ZF-AMP, LR-ZF-SIC, and LR-ZF-SIC-AMP. For all methods, the performance degrades as N increases, and AMP becomes less effective. In particular, at N = 32, the performances of LR-ZF-SIC and LR-ZFSIC-AMP almost overlap. This is presumably because small channel eigenvalues are more likely to appear as N increases. Specifically, the Gaussian assumption on the AMP residual is increasingly violated. For the proposed method, we first examine HMC w/o LR (σ = 1). The degradation becomes more severe as N increases, and the performance remains around that of LRZF. This is because, as N increases, the channel hardening effect [3] reduces the relative channel-dependent fluctuation in the VP cost, so that a narrower search around the prior modes becomes sufficient. As a result, the optimal values of σlik and σpri are considered to decrease as N increases, as shown in Table II, whereas fixing both to 1 degrades the performance. Next, we examine HMC w/o LR, with σlik and σpri set to their optimal values. The performance improves substantially over HMC w/o LR (σ = 1), but a floor tends to remain in the highSNR region. At high SNR, the AWGN amplitude becomes smaller, so the VP cost has a relatively larger effect on noise enhancement. As a result, the search performance for the perturbation vector becomes more apparent. Next, examining HMC w/ LR, we find that, except for a marginal reversal in the low-to-mid SNR region at N = 32, it largely improves the performance of HMC w/o LR and also eliminates the floor in the high-SNR region. This indicates that applying LR reduces the difficulty of the search and enables better perturbation vectors to be found. Furthermore, comparing HMC w/ LR with LR-ZF-SIC-AMP, we find that HMC w/

LR outperforms it at all SNRs. Specifically, the gain at SER = 10−3 is 0.7, 1.2, and 1.0 dB for N = 8, 16, 32, respectively. We further compare HMC w/ LR with HA. The degradation from HA at SER = 10−3 is limited to 2.0, 1.5, and 2.4 dB for N = 8, 16, 32, respectively. Thus, HMC w/ LR shows good results, although its performance improves from N = 8 to 16 and degrades at N = 32. This may be attributed to the competing effects of channel hardening and small eigenvalues, combined with the fixed LHMC . In particular, we preliminarily confirmed that the degradation at N = 32 is partly related to the limited LHMC , as suggested by the weaker saturation in Fig. 3; the full clarification of this N dependence is left for future work. In summary, we confirmed that, on general-purpose computing platforms, the proposed continuous relaxation approach outperforms conventional AMP-based methods and achieves an SER close to the performance limit of VP. E. Computational Complexity We compare the number of multiplications in the iterative search between the proposed method and LR-ZF-SIC-AMP, excluding one-time preprocessing common to both. In the proposed method, the dominant cost is the matrix-vector product (GH LR GLR )ℓLR in the evaluation of the gradient of the logposterior distribution, which requires N 2 operations, since (GH LR GLR ) can be computed once in advance. This gradient evaluation is performed L (≥ 1) times when numerically solving Hamilton’s equations in line 4 of Algorithm 1. Since L varies adaptively under the default settings of Stan, we use a representative value for the estimation. In this paper, since the number of parallel HMC chains is 5 and LHMC = 400, the complexity is 5 · 400 · LN 2 . In contrast, the dominant cost in LR-ZF-SIC-AMP is the matrix-vector product in the AMP processing, which requires M N operations; since M = N in our numerical experiments, this amounts to N 2 operations. This operation appears at three places in the AMP processing loop [13, Alg. 1]: lines 6&7, 8, and 10. Based on [13], the number of AMP loops is set to 20, so the complexity is 20·3·N 2 . Comparing the two methods, both are of polynomial order O(N 2 ), but the proposed method is computationally

Fig. 5. LR-ZF-SIC-AMP with matched complexity versus HMC w/ LR (M = N = 8).

heavier than LR-ZF-SIC-AMP. Here, we examine the performance of LR-ZF-SIC-AMP when the number of AMP loops is increased so that the complexity is approximately matched. For a conservative comparison with LR-ZF-SICAMP, we consider M = N = 8, where the AMP processing is most effective. Since the median of L in the simulation of Fig. 4(a) was 15, we increase the number of AMP loops up to 5 · 400 · 15/3 = 10000. This performance comparison is shown in Fig. 5. The performance of LR-ZF-SIC-AMP is almost saturated and improves only marginally, so the proposed method still outperforms it. This suggests that the performance gain of the proposed method is not merely due to the allocation of computational resources, but stems from the inherent effectiveness of the method. That is, the proposed method extends the achievable SER performance. V. C ONCLUSION This paper has shown that the proposed continuous relaxation approach can achieve performance close to the performance limit of VP in the MU-MIMO downlink, on generalpurpose computing platforms and with polynomial complexity. This result indicates a promising direction toward the practical use of high-performance nonlinear precoding for 6G wireless communications. In this paper, the parameters affecting the performance of the proposed method were optimized empirically, and clarifying their mathematical background is left for future work. This paper also recasts the integer perturbation search as a probabilistic problem in a continuous space, achieving good performance even in high-dimensional systems. We have applied a similar probabilistic inference framework to MIMO signal detection at the receiver [19], demonstrating that this framework has the versatility to be applied to both transmitter- and receiver-side problems. These results are expected to contribute to broadening the knowledge of continuous relaxation for combinatorial optimization. R EFERENCES [1] D. Gesbert, M. Kountouris, R. W. Heath, Jr., C.-B. Chae, and T. Sälzer, “Shifting the MIMO paradigm,” IEEE Signal Process. Mag., vol. 24, no. 5, pp. 36–46, Sep. 2007.

[2] M. A. Albreem, A. H. Al Habbash, A. M. Abu-Hudrouss, and S. S. Ikki, “Overview of precoding techniques for massive MIMO,” IEEE Access, vol. 9, pp. 60 764–60 801, 2021. [3] A. Chockalingam and B. Rajan, Large MIMO Systems. Cambridge, UK: Cambridge University Press, 2014. [4] M. Tomlinson, “New automatic equaliser employing modulo arithmetic,” Electron. Lett., vol. 7, no. 5-6, pp. 138–139, Mar. 1971. [5] H. Harashima and H. Miyakawa, “Matched-transmission technique for channels with intersymbol interference,” IEEE Trans. Commun., vol. 20, no. 4, pp. 774–780, Aug. 1972. [6] C. B. Peel, B. M. Hochwald, and A. L. Swindlehurst, “A vectorperturbation technique for near-capacity multiantenna multiuser communication — part I: channel inversion and regularization,” IEEE Trans. Commun., vol. 53, no. 1, pp. 195–202, Jan. 2005. [7] B. M. Hochwald, C. B. Peel, and A. L. Swindlehurst, “A vectorperturbation technique for near-capacity multiantenna multiuser communication — part II: perturbation,” IEEE Trans. Commun., vol. 53, no. 3, pp. 537–544, Mar. 2005. [8] 3GPP, “NR; Physical layer procedures for data,” 3rd Generation Partnership Project (3GPP), Technical Specification TS 38.214, Mar. 2024, version 17.9.0, Release 17. [Online]. Available: https: //www.3gpp.org/ftp/Specs/archive/38 series/38.214/ [9] Signals Research Group, “5G Multi-User MIMO — it isn’t just for downlink anymore,” Samsung Networks White Paper, Feb. 2024, [Online]. Available: https://www.samsung.com/global/business/networ ks/insights/white-papers/5g-multi-user-mimo-it-isnt-just-for-downlin k-anymore/. [10] M. H. M. Costa, “Writing on dirty paper (corresp.),” IEEE Trans. Inf. Theory, vol. 29, no. 3, pp. 439–441, May 1983. [11] L. Lu, G. Y. Li, A. L. Swindlehurst, A. Ashikhmin, and R. Zhang, “An overview of massive MIMO: Benefits and challenges,” IEEE J. Sel. Top. Signal Process., vol. 8, no. 5, pp. 742–758, Oct. 2014. [12] A. Li, Y. Ma, and R. Tafazolli, “Correlation matching pursuit for MIMO vector perturbation,” in Proc. IEEE Global Commun. Conf. (GLOBECOM), Dec. 2023, pp. 6664–6669. [13] S. Lyu and C. Ling, “Hybrid vector perturbation precoding: The blessing of approximate message passing,” IEEE Trans. Signal Process., vol. 67, no. 1, pp. 178–193, Jan. 2019. [14] S. Kasi, A. K. Singh, D. Venturelli, and K. Jamieson, “Quantum annealing for large MIMO downlink vector perturbation precoding,” in Proc. IEEE Int. Conf. Commun. (ICC), Jun. 2021, pp. 1–6. [15] S. Winter, Y. Zhang, G. Zheng, and L. Hanzo, “A lattice-reduction aided vector perturbation precoder relying on quantum annealing,” IEEE Wireless Commun. Lett., vol. 13, no. 5, pp. 1225–1229, May 2024. [16] S. Duane, A. Kennedy, B. J. Pendleton, and D. Roweth, “Hybrid Monte Carlo,” Phys. Lett. B, vol. 195, no. 2, pp. 216–222, Sep. 1987. [17] S. Kakuta, J. Hagiwara, T. Nishimura, T. Sato, Y. Ogawa, and T. Ohgane, “A study on vector perturbation method using MCMC in massive MIMO,” IEICE Tech. Rep., vol. 123, no. 76, pp. 13–18, Jun. 2023, Paper RCS2023-30 (in Japanese without peer review). [18] D. Wübben, D. Seethaler, J. Jaldén, and G. Matz, “Lattice reduction: A survey with applications in wireless communications,” IEEE Signal Process. Mag., vol. 28, no. 3, pp. 70–91, May 2011. [19] J. Hagiwara, K. Matsumura, H. Asumi, Y. Kasuga, T. Nishimura, T. Sato, Y. Ogawa, and T. Ohgane, “Near-optimal stochastic MIMO signal detection with a mixture of t-distributions prior,” in Proc. IEEE Global Commun. Conf. (GLOBECOM), Dec. 2023, pp. 6627–6633. [20] D. J. Ryan, I. B. Collings, I. V. L. Clarkson, and R. W. Heath, Jr., “Performance of vector perturbation multiuser MIMO systems with limited feedback,” IEEE Trans. Commun., vol. 57, no. 9, pp. 2633– 2644, Sep. 2009. [21] A. K. Lenstra, H. W. Lenstra, Jr., and L. Lovász, “Factoring polynomials with rational coefficients,” Mathematische Annalen, vol. 261, no. 4, pp. 515–534, Dec. 1982. [22] Stan Development Team, Stan Reference Manual, Stan Development Team, 2025, version 2.37.0. [Online]. Available: https://mc-stan.org [23] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin, Bayesian Data Analysis, 3rd ed. Boca Raton, FL, USA: CRC Press, 2013. [24] D. Tse and P. Viswanath, Fundamentals of Wireless Communication. Cambridge, U.K.: Cambridge Univ. Press, 2005.