Multiuser Detection with Sparse Spectrum ... - Semantic Scholar
164
IEEE COMMUNICATIONS LETTERS, VOL. 16, NO. 2, FEBRUARY 2012
Multiuser Detection with Sparse Spectrum Gaussian Process Regression Shaowei Wang, Member, IEEE, and Hualai Gu
Abstract—Multiuser detection in direct-sequence code-division multiple access (DS-CDMA) systems can be implemented by using Gaussian process (GP) for regression in the sense of minimum mean square error criterion. In this Letter we investigate the application of sparse spectrum Gaussian process (SSGP) to the multiuser detection problem. The key point of the SSGP is that the sparsity of spectral representation of Gaussian process leads to an algorithm with much lower complexity than the full GP, while keeping almost the same bit error rate (BER) for DS-CDMA systems. Experimental results validate our proposed SSGP based multiuser detection method. It achieves greater efficiency and good BER performance. Index Terms—Code-division multiple access, machine learning, multiuser detection, and sparse spectrum gaussian process.
I. I NTRODUCTION
M
ULTIUSER detection (MUD) is an important technique to mitigate multiple access interference in directsequence code-division multiple access (DS-CDMA) systems. Though optimum MUD can attain the lowest bit-error-rate (BER) performance [1], its computational complexity increases exponentially with the number of active users, which makes it impractical for applications. Minimum mean square error (MMSE), which is intensively studied in machine learning field, is usually taken as an effective criterion to tackle the problems arising in communications. For this reason, there are many MUD schemes which are based on machine learning methods, such as neural networks [2–4], support vector machines [5], and Gaussian process (GP) for regression [6– 8]. Among these methods, GP for regression is more attractive because it needs only small training samples to attain the same BER performance as the others, which is more promising for wireless systems. The time and memory requirements of the GP are O(𝑁 3 ) and O(𝑁 2 ) [9, 10], respectively, where 𝑁 is the number of training samples. These computational loads, especially the time complexity, are still too high for wireless systems, especially for time sensitive applications, such as the MUD. It is necessary to adopt computationally efficient methods to approximate the GP. Some approximations have been proposed in [11–13], among which sparse spectrum Gaussian process (SSGP) [11] seems appropriate to the MUD problem because it can approximate any desired stationary full GP as discussed in [11]. Manuscript received July 11, 2011. The associate editor coordinating the review of this letter and approving it for publication was H. Shafiee. This work was partially supported by the Jiangsu Science Foundation (No. BK2011051) and the Fundamental Research Funds for the Central Universities (No.1095021029, No.1118021011). The authors are with the School of Electronic Science & Engineering, Nanjing University, Nanjing 210093, China (e-mail: [email protected]). Digital Object Identifier 10.1109/LCOMM.2011.120211.111508
In this Letter we adopt the SSGP to implement the MUD in DS-CDMA systems. We have found that, with appropriate number of basis functions, the SSGP can obtain almost the same BER performance as the full GP for the MUD problem. The remarkable characteristic of the SSGP is that it reduces the time cost to 𝑂(𝑁 𝑀 2 ) and the storage requirement to 𝑂(𝑁 𝑀 ), where 𝑀 is the number of basis functions and generally much smaller than 𝑁 . The complexity is reduced dramatically therefor. II. MMSE M ULTIUSER D ETECTION Consider a synchronous DS-CDMA communication system with binary phase shift keying modulation through an additive white Gaussian noise channel, the number of users is 𝐾, the length of spreading sequences is 𝑀 , and the spreading sequence for the 𝑘th user is 𝒔𝑘 = [𝑠𝑘,1 , . . . , 𝑠𝑘,𝑀 ]⊤ , 𝑠𝑘,𝑖 ∈ {−1, +1} for 𝑖 = 1, 2, ..., 𝑀 . Let 𝑏𝑡 (𝑘) denote the vector of the transmitted bit of the 𝑘th user at time 𝑡, the received signals can be represented in matrix notation as 𝒙𝑡 = 𝑺𝑨𝒃𝑡 + 𝒏𝑡 ,
(1)
where 𝑺 = [𝒔1 , . . . , 𝒔𝐾 ] is an 𝑀 × 𝐾 matrix whose 𝑘th column is the spreading sequence of the 𝑘th user, 𝑨 is a 𝐾 × 𝐾 diagonal matrix whose 𝑘th diagonal element is the amplitude of the 𝑘th user, 𝒃𝑡 = [𝑏𝑡 (1), 𝑏𝑡 (2), ..., 𝑏𝑡 (𝐾)]⊤ and 𝒏𝑡 is a real-valued zero-mean Gaussian random vector with 2 2 covariance 𝔼[𝒏𝑡 𝒏⊤ 𝑡 ] = 𝜎 𝑰, 𝜎 is the power spectral density of Gaussian noise. We estimate the transmitted bit 𝑏𝑡 (𝑘) as follows, ˆ𝑏𝑡 (𝑘) = sign(𝒘 ⊤ 𝒙𝑡 ), 𝑘
(2)
where 𝒘𝑘 is the weight vector for the 𝑘th user and the optimal value of 𝒘 𝑘 in the sense of MMSE is following, −1 2 𝒘 ∗𝑘 = arg min 𝐸[(𝑏𝑡 (𝑘) − 𝒘⊤ 𝑘 𝒙𝑡 ) ] = 𝑹𝒙𝒙 𝑹𝒙𝒃 , 𝒘𝑘
(3)
where 𝑹𝒙𝒙 = 𝐸[𝒙𝑡 𝒙⊤ 𝑡 ] is the autocorrelation of the inputs, and 𝑹𝒙𝒃 = 𝐸[𝒙𝑡 𝑏𝑡 (𝑘)] is the cross-correlation between the inputs and outputs, respectively. Given a set of training data, we can obtain 𝒘∗𝑘 by solving (3). Insert 𝒘∗𝑘 into (2), ˆ𝑏𝑡 (𝑘) can be worked out subsequently. III. S PARSE S PECTRUM G AUSSIAN P ROCESS FOR THE MUD P ROBLEM A. Gaussian Process for Regression Collect a labeled training data set with 𝑁 samples 𝒟 = 𝑑 {𝒙𝑖 , 𝑏𝑖 (𝑘)}𝑁 𝑖=1 , where 𝒙𝑖 ∈ ℝ is the received symbols vector and 𝑏𝑖 (𝑘) is the transmitted symbol of the 𝑘th user. Let 𝒃(𝑘) = [𝑏1 (𝑘), . . . , 𝑏𝑁 (𝑘)]⊤ denote the transmitted symbol
c 2012 IEEE 1089-7798/12$31.00 ⃝
WANG and GU: MULTIUSER DETECTION WITH SPARSE SPECTRUM GAUSSIAN PROCESS REGRESSION
vector, the linear regression model of the GP can be written as follows [14], 𝒃(𝑘) = Φ⊤ 𝒘 + 𝒏, 𝒏 ∼ 𝒩 (0, 𝜎𝑛2 𝑰 𝑁 ),
(4)
where Φ = [𝝓(𝒙1 ), . . . , 𝝓(𝒙𝑁 )] and 𝝓(𝒙𝑖 ) is the vector of fixed basis functions which depend on the input vector 𝒙𝑖 , 𝒘 is the weight parameter vector, 𝜎𝑛2 is the variance of Gaussian noise of the samples, and 𝑰 𝑁 is an identity matrix of size 𝑁 . Denote 𝑿 = [𝒙1 , 𝒙2 , ..., 𝒙𝑁 ]⊤ , the distribution of 𝒃(𝑘) conditioned on 𝑿 is 𝑝(𝒃(𝑘)∣𝑿, 𝒘) = 𝒩 (𝒃(𝑘)∣Φ⊤ 𝒘, 𝜎𝑛2 𝑰 𝑁 ).
(5)
Assume 𝒘 is a zero-mean Gaussian random variable with 2 𝑝(𝒘) = 𝒩 (𝒘∣0, 𝜎𝒘 𝑰), the posterior of 𝒘 can be computed by Bayes’ rule, 𝑝(𝒘∣𝒃(𝑘), 𝑿) =
𝑝(𝒘)𝑝(𝒃(𝑘)∣𝒘, 𝑿) . 𝑝(𝒃(𝑘)∣𝑿)
(6)
165
complexity O(𝑁 2 ). The computation cost is too high for application because the MUD in DS-CDMA systems is time sensitive. It stimulates us to explore more efficient ways to implement the GP based MUD scheme. We propose a low complexity method, named sparse spectrum Gaussian process (SSGP) [11], meet this requirement. B. Sparse spectrum Gaussian process Assume the labeled training samples are drawn from a stationary Gaussian process, the autocorrelation function is equivalent to the stationary covariance function. Furthermore, the power spectrum and the autocorrelation of a stationary random process constitute a Fourier pair [16], ∫ ⊤ 𝑘(𝝉 ) = ℝ𝑑 𝑒2𝜋𝑖𝒇 𝝉 𝐹 (𝒇 )𝑑𝒇 ∫ (13) ⊤ 𝐹 (𝒇 ) = ℝ𝑑 𝑒−2𝜋𝑖𝒇 𝝉 𝑘(𝝉 )𝑑𝝉 ,
where 𝐹 (𝒇 ) is a positive finite measure based on Bochner’s theorem [17]. We can sparsity the power spectral density Since only the likelihood function and the prior probability and obtain a sparse covariance function. Denote 𝐹 (𝒇 ) = depend on the weights, we can write the posterior of 𝒘 as 𝑘(0)𝑝 (𝒇 ), where 𝑝 (𝒇 ) is a probability measure [11], the 𝐹 𝐹 follows [10], kernel function, also the covariance function, can be written 𝑝(𝒘∣𝒃(𝑘), 𝑿) as ∝ 𝑒𝑥𝑝(− 2𝜎12 (𝒃(𝑘) − Φ⊤ 𝒘)⊤ (𝒃(𝑘) − Φ⊤ 𝒘))𝑒𝑥𝑝(− 2𝜎12 𝒘⊤ 𝒘) 𝑘(𝒙 , 𝒙 ) = 𝑘(0) ∫ 𝑒2𝜋𝑖𝒇 ⊤ 𝒙𝑖 (𝑒2𝜋𝑖𝒇 ⊤ 𝒙𝑗 )∗ 𝑝 (𝒇 )𝑑𝒇 𝑛 𝒘 𝑖 𝑗 𝐹 ℝ𝑑 ⊤ ⊤ ¯ ⊤ ( 𝜎12 ΦΦ⊤ + 𝜎12 𝑰)(𝒘 − 𝒘)), ¯ ∝ 𝑒𝑥𝑝(− 12 (𝒘 − 𝒘) = 𝑘(0)𝔼𝑝𝐹 [𝑒2𝜋𝑖𝒇 𝒙𝑖 (𝑒2𝜋𝑖𝒇 𝒙𝑗 )∗ ]. 𝑛 𝒘 (7) (14) ¯ = 𝜎12 ( 𝜎12 ΦΦ⊤ + 𝜎12 𝑰)−1 Φ𝒃(𝑘). So the posterior where 𝒘 We estimate the covariance function by a Monte Carlo 𝑛 𝑛 𝒘 distribution of 𝒘 can be rewritten as method, which takes an average of a few samples correspond−1 ¯ 𝑩 ), 𝑝(𝒘∣𝒃(𝑘), 𝑿) = 𝒩 (𝒘∣𝒘, (8) ing to a finite set of frequencies. The integral in the right-hand side of (14) can be approximated as shown in (15), where 𝒇 𝑟 where 𝑩 = 𝜎12 𝑰 + 𝜎12 ΦΦ⊤ . Then the maximum a posteriori is drawn from 𝑝𝐹 (𝒇 ) independently. And the vector 𝝓(𝒙) 𝒘 𝑛 (MAP) estimate of 𝒘 can be computed as follows, contains 𝑀 = 2𝑢 nonlinear basis functions, as shown in (16). Using (10), (15) and (16), we can calculate the estimation of 𝒘 𝑀𝐴𝑃 = arg max log 𝑝(𝒘∣𝒃(𝑘), 𝑿) 𝒘 𝑏(𝑘) as follows, ¯ ⊤ 𝑩(𝒘 − 𝒘)} ¯ = arg min{ 21 (𝒘 − 𝒘) (9) 𝒘 ˆ𝑏(𝑘) = 𝝓(𝒙)⊤ 𝑨−1 Φ𝒃(𝑘), (17) = arg min{∥ 𝒃(𝑘) − Φ⊤ 𝒘 ∥2 +𝜆 ∥ 𝒘 ∥2 }, 𝒘 ⊤ where Φ = [𝝓(𝒙1 ), . . . , 𝝓(𝒙𝑁 )] and 𝑨 = ΦΦ + 2 . The last term of (9), 𝜆 ∥ 𝒘 ∥2 , can be 𝑢𝜎 2 /𝜎 2 𝑰 . where 𝜆 = 𝜎𝑛2 /𝜎𝒘 𝑛 𝑤 2𝑢 regarded as a regularization term to avoid overfitting [15]. By working out the 𝒘𝑀𝐴𝑃 , 𝑏(𝑘) can be estimated as C. On the Choice of the Hyperparameters ˆ𝑏(𝑘) = 𝝓(𝒙)⊤ Φ(Φ⊤ Φ + 𝜆𝑰)−1 𝒃(𝑘) All 𝒇 𝑖 s in (16), 𝑖 = 1, 2, ..., 𝑢, are hyperparameters, as (10) = 𝒌𝑪 −1 𝒃(𝑘), well as 𝜎𝑛 and 𝜎𝑤 . We have to determine the values of these where 𝑘(𝒙𝑖 , 𝒙𝑗 ) is the kernel function that can be expressed as 𝑘(𝒙𝑖 , 𝒙𝑗 ) = 𝝓(𝒙𝑖 )⊤ 𝝓(𝒙𝑗 ), 𝑪 𝑖𝑗 = 𝑘(𝒙𝑖 , 𝒙𝑗 ) + 𝜆𝛿𝑖𝑗 and 𝒌 = [𝑘(𝒙, 𝒙1 ), . . . , 𝑘(𝒙, 𝒙𝑁 )]. Obviously, it is necessary to know the kernel function 𝑘(⋅) to calculate ˆ𝑏(𝑘). Denote 𝜽 = [𝜃1 , 𝜃2 , 𝜆], where 𝜃1 , 𝜃2 and 𝜆 are hyperparameters, and assume that 𝑘(⋅) has the following form [14], 𝑘(𝒙𝑖 , 𝒙𝑗 ) = 𝜃1 exp(𝜃2 ∥ 𝒙𝑖 − 𝒙𝑗 ∥2 ),
(11)
the optimal solution of the hyperparameters can be computed as follows, 𝜽∗
= arg min 𝑝(𝒃(𝑘)∣𝑿, 𝜽) 𝜽
1 ∣2𝜋𝑪∣
= arg min √ 𝜽
exp(− 12 𝒃(𝑘)⊤ 𝑪 −1 𝒃(𝑘)).
(12)
From (12) we know that the inversion of 𝑪 is needed to work out the 𝜽 ∗ . Since 𝑪 is an 𝑁 × 𝑁 matrix, the complexity of matrix inversion is generally O(𝑁 3 ) with storage
hyperparameters to predict the mean or the variance of each ⊤ point in the test set. Denote 𝜻 = [𝒇 ⊤ 1 , . . . , 𝒇 𝑢 , 𝜎𝑛 , 𝜎𝑤 ], the logarithmic marginal likelihood of the SSGP follows [11], log 𝑝(𝒃(𝑘)∣𝑿, 𝜻) =
[𝒃(𝑘)⊤ 𝒃(𝑘) − 𝒃(𝑘)⊤ Φ⊤ 𝑨−1 Φ𝒃(𝑘)] 1 − log ∣𝑨∣ 2𝜎𝑛2 2 2 𝑢𝜎 𝑁 log 2𝜋𝜎𝑛2 . + 𝑢 log 2𝑛 − (18) 𝜎𝑤 2
−
Calculating the partial derivatives of (18) with respect to 𝜻 and equating them to zeros, we get a set of nonlinear equations whose solutions are the settings of the hyperparameters vector 𝜻. We can solve these equations with numerical optimization techniques to obtain promising settings of 𝜻. It is worth noting that the logarithmic marginal likelihood (18) is nonconvex, the solution of the hyperparameters vector
166
IEEE COMMUNICATIONS LETTERS, VOL. 16, NO. 2, FEBRUARY 2012
𝑢
𝑘(𝒙𝑖 , 𝒙𝑗 ) ≃
⊤ ⊤ ⊤ 𝑘(0) ∑ 2𝜋𝑖𝒇 ⊤ 𝑟 𝒙𝑖 (𝑒2𝜋𝑖𝒇 𝑟 𝒙𝑗 )∗ + (𝑒2𝜋𝑖𝒇 𝑟 𝒙𝑖 )∗ 𝑒2𝜋𝑖𝒇 𝑟 𝒙𝑗 ] [𝑒 2𝑢 𝑟=1
𝑢
=
𝑘(0) ∑ cos(2𝜋𝒇 ⊤ 𝑟 (𝒙𝑖 − 𝒙𝑗 )), 𝑢 𝑟=1
(15)
⊤ ⊤ ⊤ ⊤ 𝝓(𝒙) = [cos(2𝜋𝒇 ⊤ 1 𝒙), sin(2𝜋𝒇 1 𝒙), . . . , cos(2𝜋𝒇 𝑢 𝒙), sin(2𝜋𝒇 𝑢 𝒙)] .
(16)
0
0.075
10
Full GP SSGP
Full GP SSGP MMSE
0.07
BER
BER
0.065 −1
10
0.06
0.055
0.05
−2
0
20
40
60
80
Number of Basis Functions
10
100
Fig. 1. BER as a function of the number of basis functions, 200 training samples and SNR=4dB.
𝜻 obtained by numerical optimization methods is generally a local optimal one, which may lead to poor BER performance in some cases for the DS-CDMA systems. To obtain robust hyperparameters, a multi-start technique is introduced. We initialize 𝜻 to multiple sensible values and optimize each 𝜻 by using numerical optimization method to produce the corresponding local optimal solution. Then we take the best 𝜻 from these local optimal solutions as the hyperparameters to predict 𝑏(𝑘). Over 100 independent experiments, that is, training data set and test data are generated randomly in each experiment, we have found that the hyperparameters determined in this way achieve better BER performance and are more robust, comparing with the method with single start. The computational complexity is O(𝑁 𝑀 2 ) for training the SSGP algorithm. At prediction time, the cost is O(𝑀 ) for the predictive mean and O(𝑀 2 ) for the predictive variance per test point. In brief, the complexity of the SSGP is much lower than O(𝑁 3 ) when 𝑀 ≪ 𝑁 . We will show it is always the case for the MUD problem. Furthermore, the storage complexity of the SSGP is O(𝑁 𝑀 ), which is also lower than O(𝑁 2 ) if 𝑀 ≪ 𝑁. IV. E XPERIMENTAL R ESULTS AND D ISCUSSIONS We conduct a series of experiments to test the performance of the proposed SSGP based MUD algorithm. Consider a DSCDMA system with 8 users, the spreading codes are random binary sequences with length 31, the channel response is following, 𝐻(𝑧) = 0.4 + 0.9𝑧 −1 + 0.4𝑧 −2 ,
(19)
0
50
100
150
200
Number of Training Samples
250
300
Fig. 2. BER as a function of the number of training samples, SNR=4dB and 𝑀 = 20.
and the training and testing data are generated randomly. Firstly, we verify how many basis functions are appropriate for the SSGP to approach the full GP for the MUD problem. Fig.1 shows the BER as a function of the number of basis functions of the SSGP. The SNR is 4 dB and there are 𝑁 = 200 training samples. All users transmit signals with equal power. As seen from Fig.1, the BER varies in a narrow range for the SSGP when the number of basis functions is larger than 20. It is reasonable to set 𝑀 = 20 to approximate the full GP for the SSGP to solve the MUD problem. Fig.2 shows the BER of the GP, the SSGP and the MMSE [15], as a function of the number of training samples. The number of basis functions of the SSGP is 20. We can see that the GP and the SSGP perform better than the MMSE when the number of training samples is relatively small. It means that we can transmit more information bits in each burst of data with the SSGP or the GP based MUD scheme. It can be also seen that the SSGP and the GP have almost the same BER performance. Fig.3 shows the BER as a function of signal to noise ratios (SNRs) for the SSGP, the full GP and the MMSE. Again, it is shown that the SSGP produces a good approximation to the GP. To test near-far resistance performance of the SSGP, we consider a case where all users but the first one have the same bit energy 𝐸𝑖 . The SNR of the first user is 6dB. The number of basis functions of the SSGP is also 𝑀 = 20. Fig.4 shows the curves of BER versus 𝐸1 /𝐸𝑖 (dB). We can see the SSGP and the GP have almost the same performance, both of which perform better than the MMSE.
WANG and GU: MULTIUSER DETECTION WITH SPARSE SPECTRUM GAUSSIAN PROCESS REGRESSION
some intensive topics, such as which kind of basis functions is suitable for a specific problem and how many basis functions are necessary to obtain an acceptable approximation.
0
10
Full GP SSGP MMSE
−1
10
167
ACKNOWLEDGEMENT −2
10
BER
The authors would like to thank the anonymous reviewers, whose invaluable comments helped improve the presentation of this paper substantially. Shaowei Wang is grateful to Prof. Zhi-Hua Zhou for his helpful suggestions for this research.
−3
10
−4
10
R EFERENCES −5
10
Fig. 3.
−4
−2
0
2
4
SNR(dB)
6
8
10
12
BER as a function of SNR, 200 training samples, 𝑀 = 20. −1
10
BER
Full GP SSGP MMSE
−2
10
−3
10
−2
−1.5
−1
−0.5
0
0.5
Amplitude ratio(dB)
1
1.5
2
Fig. 4. Near-far resistance performance, 200 training samples, SNR=6dB and 𝑀 = 20.
V. C ONCLUSION Gaussian process (GP) for regression exhibits good performance with short training samples for the multiuser detection problem, but the complexity is too high. We adopt a sparse spectrum gaussian process algorithm to approximate the GP. The proposed algorithm has much lower complexity than the GP while keeping almost the same BER performance. The method given in this Letter throws some lights on designing more efficient digital communication receivers if considering
[1] S. Verdu, Multiuser Detection. Cambridge University Press, 1998. [2] B. Aazhang, B. P. Paris, and G. C. Orsak, “Neural networks for multiuser detection in code-division multiple-access communications,” IEEE Trans. Commun., vol. 40, no. 7, pp. 1212–1222, July 1992. [3] U. Mitra and H. V. Poor, “Neural network techniques for adaptive multiuser demodulation,” IEEE J. Sel. Areas Commun., vol. 12, no. 9, pp. 1460–1470, Dec. 1994. [4] B. Aazhang, B. P. Paris, and G. C. Orsak, “Neural networks for multiuser detection in code-division multiple-access communications,” IEEE Trans. Commun., vol. 40, no. 7, pp. 1212–1222, July 1992. [5] S. Chen, A. K. Samingan, and L. Hanzo, “Support vector machine multiuser receiver for DS-CDMA signals in multipath channels,” IEEE Trans. Neural Networks, vol. 12, no. 3, pp. 604–611, May 2001. [6] J. Murillo-Fuentes and F. Perez-Cruz, “Gaussian process regressors for multiuser detection in DS-CDMA systems,” IEEE Trans. Commun., vol. 57, no. 8, pp. 2339–2347, Aug. 2009. [7] J. Murillo-Fuentes, S. Caro, and F. Perez-Cruz, “Gaussian process for multiuser detection in CDMA receivers,” in Proc. 2005 Advances in Neural Information Processing Systems. [8] F. Perez-Cruz, J. Murillo-Fuentes, and S. Caro, “Nonlinear channel equalization with Gaussian processes for regression,” IEEE Trans. Signal Process., vol. 56, no. 10, pp. 5283–5286, Oct. 2008. [9] C. K. I. Williams, “Prediction with Gaussian processes: from linear regression to linear prediction and beyond,” Learning and Inference in Graphical Models, M. Jordan, editor. Kluwer Academic Press, 1998. [10] C. E. Rasmussen and C. K. I. Williams, Gaussian Process for Machine Learning. MIT Press, 2006. [11] M. Lazaro-Gredilla, J. Quinonero-Candela, C. Rasmussen, and A. Figueiras-Vidal, “Sparse spectrum Gaussian process regression,” J. Mach. Learn. Res., vol. 11, pp. 1865–1881, June 2010. [12] E. Snelson and Z. Ghahramani, “Sparse Gaussian processes using pseudo-inputs,” in Proc. 2006 Advances in Neural Information Processing Systems. [13] C. Walder, K. I. Kim, and B. Scholkopf, “Sparse multiscale Gaussian process regression,” in Proc. 2008 International Conference on Machine Learning. [14] M. Jordan, Learning in Graphical Models. MIT Press, 1998. [15] L. Rugini, P. Banelli, and S. Cacopardi, “A full-rank regularization technique for MMSE detection in multiuser CDMA systems,” IEEE Commun. Lett., vol. 9, no. 1, pp. 34–36, Jan. 2005. [16] A. Carlson, Communication Systems. McGraw-Hill, 1986. [17] M. Stein, Interpolation of Spatial Data. Springer Verlag, 1999.
IEEE COMMUNICATIONS LETTERS, VOL. 16, NO. 2, FEBRUARY 2012
Multiuser Detection with Sparse Spectrum Gaussian Process Regression Shaowei Wang, Member, IEEE, and Hualai Gu
Abstract—Multiuser detection in direct-sequence code-division multiple access (DS-CDMA) systems can be implemented by using Gaussian process (GP) for regression in the sense of minimum mean square error criterion. In this Letter we investigate the application of sparse spectrum Gaussian process (SSGP) to the multiuser detection problem. The key point of the SSGP is that the sparsity of spectral representation of Gaussian process leads to an algorithm with much lower complexity than the full GP, while keeping almost the same bit error rate (BER) for DS-CDMA systems. Experimental results validate our proposed SSGP based multiuser detection method. It achieves greater efficiency and good BER performance. Index Terms—Code-division multiple access, machine learning, multiuser detection, and sparse spectrum gaussian process.
I. I NTRODUCTION
M
ULTIUSER detection (MUD) is an important technique to mitigate multiple access interference in directsequence code-division multiple access (DS-CDMA) systems. Though optimum MUD can attain the lowest bit-error-rate (BER) performance [1], its computational complexity increases exponentially with the number of active users, which makes it impractical for applications. Minimum mean square error (MMSE), which is intensively studied in machine learning field, is usually taken as an effective criterion to tackle the problems arising in communications. For this reason, there are many MUD schemes which are based on machine learning methods, such as neural networks [2–4], support vector machines [5], and Gaussian process (GP) for regression [6– 8]. Among these methods, GP for regression is more attractive because it needs only small training samples to attain the same BER performance as the others, which is more promising for wireless systems. The time and memory requirements of the GP are O(𝑁 3 ) and O(𝑁 2 ) [9, 10], respectively, where 𝑁 is the number of training samples. These computational loads, especially the time complexity, are still too high for wireless systems, especially for time sensitive applications, such as the MUD. It is necessary to adopt computationally efficient methods to approximate the GP. Some approximations have been proposed in [11–13], among which sparse spectrum Gaussian process (SSGP) [11] seems appropriate to the MUD problem because it can approximate any desired stationary full GP as discussed in [11]. Manuscript received July 11, 2011. The associate editor coordinating the review of this letter and approving it for publication was H. Shafiee. This work was partially supported by the Jiangsu Science Foundation (No. BK2011051) and the Fundamental Research Funds for the Central Universities (No.1095021029, No.1118021011). The authors are with the School of Electronic Science & Engineering, Nanjing University, Nanjing 210093, China (e-mail: [email protected]). Digital Object Identifier 10.1109/LCOMM.2011.120211.111508
In this Letter we adopt the SSGP to implement the MUD in DS-CDMA systems. We have found that, with appropriate number of basis functions, the SSGP can obtain almost the same BER performance as the full GP for the MUD problem. The remarkable characteristic of the SSGP is that it reduces the time cost to 𝑂(𝑁 𝑀 2 ) and the storage requirement to 𝑂(𝑁 𝑀 ), where 𝑀 is the number of basis functions and generally much smaller than 𝑁 . The complexity is reduced dramatically therefor. II. MMSE M ULTIUSER D ETECTION Consider a synchronous DS-CDMA communication system with binary phase shift keying modulation through an additive white Gaussian noise channel, the number of users is 𝐾, the length of spreading sequences is 𝑀 , and the spreading sequence for the 𝑘th user is 𝒔𝑘 = [𝑠𝑘,1 , . . . , 𝑠𝑘,𝑀 ]⊤ , 𝑠𝑘,𝑖 ∈ {−1, +1} for 𝑖 = 1, 2, ..., 𝑀 . Let 𝑏𝑡 (𝑘) denote the vector of the transmitted bit of the 𝑘th user at time 𝑡, the received signals can be represented in matrix notation as 𝒙𝑡 = 𝑺𝑨𝒃𝑡 + 𝒏𝑡 ,
(1)
where 𝑺 = [𝒔1 , . . . , 𝒔𝐾 ] is an 𝑀 × 𝐾 matrix whose 𝑘th column is the spreading sequence of the 𝑘th user, 𝑨 is a 𝐾 × 𝐾 diagonal matrix whose 𝑘th diagonal element is the amplitude of the 𝑘th user, 𝒃𝑡 = [𝑏𝑡 (1), 𝑏𝑡 (2), ..., 𝑏𝑡 (𝐾)]⊤ and 𝒏𝑡 is a real-valued zero-mean Gaussian random vector with 2 2 covariance 𝔼[𝒏𝑡 𝒏⊤ 𝑡 ] = 𝜎 𝑰, 𝜎 is the power spectral density of Gaussian noise. We estimate the transmitted bit 𝑏𝑡 (𝑘) as follows, ˆ𝑏𝑡 (𝑘) = sign(𝒘 ⊤ 𝒙𝑡 ), 𝑘
(2)
where 𝒘𝑘 is the weight vector for the 𝑘th user and the optimal value of 𝒘 𝑘 in the sense of MMSE is following, −1 2 𝒘 ∗𝑘 = arg min 𝐸[(𝑏𝑡 (𝑘) − 𝒘⊤ 𝑘 𝒙𝑡 ) ] = 𝑹𝒙𝒙 𝑹𝒙𝒃 , 𝒘𝑘
(3)
where 𝑹𝒙𝒙 = 𝐸[𝒙𝑡 𝒙⊤ 𝑡 ] is the autocorrelation of the inputs, and 𝑹𝒙𝒃 = 𝐸[𝒙𝑡 𝑏𝑡 (𝑘)] is the cross-correlation between the inputs and outputs, respectively. Given a set of training data, we can obtain 𝒘∗𝑘 by solving (3). Insert 𝒘∗𝑘 into (2), ˆ𝑏𝑡 (𝑘) can be worked out subsequently. III. S PARSE S PECTRUM G AUSSIAN P ROCESS FOR THE MUD P ROBLEM A. Gaussian Process for Regression Collect a labeled training data set with 𝑁 samples 𝒟 = 𝑑 {𝒙𝑖 , 𝑏𝑖 (𝑘)}𝑁 𝑖=1 , where 𝒙𝑖 ∈ ℝ is the received symbols vector and 𝑏𝑖 (𝑘) is the transmitted symbol of the 𝑘th user. Let 𝒃(𝑘) = [𝑏1 (𝑘), . . . , 𝑏𝑁 (𝑘)]⊤ denote the transmitted symbol
c 2012 IEEE 1089-7798/12$31.00 ⃝
WANG and GU: MULTIUSER DETECTION WITH SPARSE SPECTRUM GAUSSIAN PROCESS REGRESSION
vector, the linear regression model of the GP can be written as follows [14], 𝒃(𝑘) = Φ⊤ 𝒘 + 𝒏, 𝒏 ∼ 𝒩 (0, 𝜎𝑛2 𝑰 𝑁 ),
(4)
where Φ = [𝝓(𝒙1 ), . . . , 𝝓(𝒙𝑁 )] and 𝝓(𝒙𝑖 ) is the vector of fixed basis functions which depend on the input vector 𝒙𝑖 , 𝒘 is the weight parameter vector, 𝜎𝑛2 is the variance of Gaussian noise of the samples, and 𝑰 𝑁 is an identity matrix of size 𝑁 . Denote 𝑿 = [𝒙1 , 𝒙2 , ..., 𝒙𝑁 ]⊤ , the distribution of 𝒃(𝑘) conditioned on 𝑿 is 𝑝(𝒃(𝑘)∣𝑿, 𝒘) = 𝒩 (𝒃(𝑘)∣Φ⊤ 𝒘, 𝜎𝑛2 𝑰 𝑁 ).
(5)
Assume 𝒘 is a zero-mean Gaussian random variable with 2 𝑝(𝒘) = 𝒩 (𝒘∣0, 𝜎𝒘 𝑰), the posterior of 𝒘 can be computed by Bayes’ rule, 𝑝(𝒘∣𝒃(𝑘), 𝑿) =
𝑝(𝒘)𝑝(𝒃(𝑘)∣𝒘, 𝑿) . 𝑝(𝒃(𝑘)∣𝑿)
(6)
165
complexity O(𝑁 2 ). The computation cost is too high for application because the MUD in DS-CDMA systems is time sensitive. It stimulates us to explore more efficient ways to implement the GP based MUD scheme. We propose a low complexity method, named sparse spectrum Gaussian process (SSGP) [11], meet this requirement. B. Sparse spectrum Gaussian process Assume the labeled training samples are drawn from a stationary Gaussian process, the autocorrelation function is equivalent to the stationary covariance function. Furthermore, the power spectrum and the autocorrelation of a stationary random process constitute a Fourier pair [16], ∫ ⊤ 𝑘(𝝉 ) = ℝ𝑑 𝑒2𝜋𝑖𝒇 𝝉 𝐹 (𝒇 )𝑑𝒇 ∫ (13) ⊤ 𝐹 (𝒇 ) = ℝ𝑑 𝑒−2𝜋𝑖𝒇 𝝉 𝑘(𝝉 )𝑑𝝉 ,
where 𝐹 (𝒇 ) is a positive finite measure based on Bochner’s theorem [17]. We can sparsity the power spectral density Since only the likelihood function and the prior probability and obtain a sparse covariance function. Denote 𝐹 (𝒇 ) = depend on the weights, we can write the posterior of 𝒘 as 𝑘(0)𝑝 (𝒇 ), where 𝑝 (𝒇 ) is a probability measure [11], the 𝐹 𝐹 follows [10], kernel function, also the covariance function, can be written 𝑝(𝒘∣𝒃(𝑘), 𝑿) as ∝ 𝑒𝑥𝑝(− 2𝜎12 (𝒃(𝑘) − Φ⊤ 𝒘)⊤ (𝒃(𝑘) − Φ⊤ 𝒘))𝑒𝑥𝑝(− 2𝜎12 𝒘⊤ 𝒘) 𝑘(𝒙 , 𝒙 ) = 𝑘(0) ∫ 𝑒2𝜋𝑖𝒇 ⊤ 𝒙𝑖 (𝑒2𝜋𝑖𝒇 ⊤ 𝒙𝑗 )∗ 𝑝 (𝒇 )𝑑𝒇 𝑛 𝒘 𝑖 𝑗 𝐹 ℝ𝑑 ⊤ ⊤ ¯ ⊤ ( 𝜎12 ΦΦ⊤ + 𝜎12 𝑰)(𝒘 − 𝒘)), ¯ ∝ 𝑒𝑥𝑝(− 12 (𝒘 − 𝒘) = 𝑘(0)𝔼𝑝𝐹 [𝑒2𝜋𝑖𝒇 𝒙𝑖 (𝑒2𝜋𝑖𝒇 𝒙𝑗 )∗ ]. 𝑛 𝒘 (7) (14) ¯ = 𝜎12 ( 𝜎12 ΦΦ⊤ + 𝜎12 𝑰)−1 Φ𝒃(𝑘). So the posterior where 𝒘 We estimate the covariance function by a Monte Carlo 𝑛 𝑛 𝒘 distribution of 𝒘 can be rewritten as method, which takes an average of a few samples correspond−1 ¯ 𝑩 ), 𝑝(𝒘∣𝒃(𝑘), 𝑿) = 𝒩 (𝒘∣𝒘, (8) ing to a finite set of frequencies. The integral in the right-hand side of (14) can be approximated as shown in (15), where 𝒇 𝑟 where 𝑩 = 𝜎12 𝑰 + 𝜎12 ΦΦ⊤ . Then the maximum a posteriori is drawn from 𝑝𝐹 (𝒇 ) independently. And the vector 𝝓(𝒙) 𝒘 𝑛 (MAP) estimate of 𝒘 can be computed as follows, contains 𝑀 = 2𝑢 nonlinear basis functions, as shown in (16). Using (10), (15) and (16), we can calculate the estimation of 𝒘 𝑀𝐴𝑃 = arg max log 𝑝(𝒘∣𝒃(𝑘), 𝑿) 𝒘 𝑏(𝑘) as follows, ¯ ⊤ 𝑩(𝒘 − 𝒘)} ¯ = arg min{ 21 (𝒘 − 𝒘) (9) 𝒘 ˆ𝑏(𝑘) = 𝝓(𝒙)⊤ 𝑨−1 Φ𝒃(𝑘), (17) = arg min{∥ 𝒃(𝑘) − Φ⊤ 𝒘 ∥2 +𝜆 ∥ 𝒘 ∥2 }, 𝒘 ⊤ where Φ = [𝝓(𝒙1 ), . . . , 𝝓(𝒙𝑁 )] and 𝑨 = ΦΦ + 2 . The last term of (9), 𝜆 ∥ 𝒘 ∥2 , can be 𝑢𝜎 2 /𝜎 2 𝑰 . where 𝜆 = 𝜎𝑛2 /𝜎𝒘 𝑛 𝑤 2𝑢 regarded as a regularization term to avoid overfitting [15]. By working out the 𝒘𝑀𝐴𝑃 , 𝑏(𝑘) can be estimated as C. On the Choice of the Hyperparameters ˆ𝑏(𝑘) = 𝝓(𝒙)⊤ Φ(Φ⊤ Φ + 𝜆𝑰)−1 𝒃(𝑘) All 𝒇 𝑖 s in (16), 𝑖 = 1, 2, ..., 𝑢, are hyperparameters, as (10) = 𝒌𝑪 −1 𝒃(𝑘), well as 𝜎𝑛 and 𝜎𝑤 . We have to determine the values of these where 𝑘(𝒙𝑖 , 𝒙𝑗 ) is the kernel function that can be expressed as 𝑘(𝒙𝑖 , 𝒙𝑗 ) = 𝝓(𝒙𝑖 )⊤ 𝝓(𝒙𝑗 ), 𝑪 𝑖𝑗 = 𝑘(𝒙𝑖 , 𝒙𝑗 ) + 𝜆𝛿𝑖𝑗 and 𝒌 = [𝑘(𝒙, 𝒙1 ), . . . , 𝑘(𝒙, 𝒙𝑁 )]. Obviously, it is necessary to know the kernel function 𝑘(⋅) to calculate ˆ𝑏(𝑘). Denote 𝜽 = [𝜃1 , 𝜃2 , 𝜆], where 𝜃1 , 𝜃2 and 𝜆 are hyperparameters, and assume that 𝑘(⋅) has the following form [14], 𝑘(𝒙𝑖 , 𝒙𝑗 ) = 𝜃1 exp(𝜃2 ∥ 𝒙𝑖 − 𝒙𝑗 ∥2 ),
(11)
the optimal solution of the hyperparameters can be computed as follows, 𝜽∗
= arg min 𝑝(𝒃(𝑘)∣𝑿, 𝜽) 𝜽
1 ∣2𝜋𝑪∣
= arg min √ 𝜽
exp(− 12 𝒃(𝑘)⊤ 𝑪 −1 𝒃(𝑘)).
(12)
From (12) we know that the inversion of 𝑪 is needed to work out the 𝜽 ∗ . Since 𝑪 is an 𝑁 × 𝑁 matrix, the complexity of matrix inversion is generally O(𝑁 3 ) with storage
hyperparameters to predict the mean or the variance of each ⊤ point in the test set. Denote 𝜻 = [𝒇 ⊤ 1 , . . . , 𝒇 𝑢 , 𝜎𝑛 , 𝜎𝑤 ], the logarithmic marginal likelihood of the SSGP follows [11], log 𝑝(𝒃(𝑘)∣𝑿, 𝜻) =
[𝒃(𝑘)⊤ 𝒃(𝑘) − 𝒃(𝑘)⊤ Φ⊤ 𝑨−1 Φ𝒃(𝑘)] 1 − log ∣𝑨∣ 2𝜎𝑛2 2 2 𝑢𝜎 𝑁 log 2𝜋𝜎𝑛2 . + 𝑢 log 2𝑛 − (18) 𝜎𝑤 2
−
Calculating the partial derivatives of (18) with respect to 𝜻 and equating them to zeros, we get a set of nonlinear equations whose solutions are the settings of the hyperparameters vector 𝜻. We can solve these equations with numerical optimization techniques to obtain promising settings of 𝜻. It is worth noting that the logarithmic marginal likelihood (18) is nonconvex, the solution of the hyperparameters vector
166
IEEE COMMUNICATIONS LETTERS, VOL. 16, NO. 2, FEBRUARY 2012
𝑢
𝑘(𝒙𝑖 , 𝒙𝑗 ) ≃
⊤ ⊤ ⊤ 𝑘(0) ∑ 2𝜋𝑖𝒇 ⊤ 𝑟 𝒙𝑖 (𝑒2𝜋𝑖𝒇 𝑟 𝒙𝑗 )∗ + (𝑒2𝜋𝑖𝒇 𝑟 𝒙𝑖 )∗ 𝑒2𝜋𝑖𝒇 𝑟 𝒙𝑗 ] [𝑒 2𝑢 𝑟=1
𝑢
=
𝑘(0) ∑ cos(2𝜋𝒇 ⊤ 𝑟 (𝒙𝑖 − 𝒙𝑗 )), 𝑢 𝑟=1
(15)
⊤ ⊤ ⊤ ⊤ 𝝓(𝒙) = [cos(2𝜋𝒇 ⊤ 1 𝒙), sin(2𝜋𝒇 1 𝒙), . . . , cos(2𝜋𝒇 𝑢 𝒙), sin(2𝜋𝒇 𝑢 𝒙)] .
(16)
0
0.075
10
Full GP SSGP
Full GP SSGP MMSE
0.07
BER
BER
0.065 −1
10
0.06
0.055
0.05
−2
0
20
40
60
80
Number of Basis Functions
10
100
Fig. 1. BER as a function of the number of basis functions, 200 training samples and SNR=4dB.
𝜻 obtained by numerical optimization methods is generally a local optimal one, which may lead to poor BER performance in some cases for the DS-CDMA systems. To obtain robust hyperparameters, a multi-start technique is introduced. We initialize 𝜻 to multiple sensible values and optimize each 𝜻 by using numerical optimization method to produce the corresponding local optimal solution. Then we take the best 𝜻 from these local optimal solutions as the hyperparameters to predict 𝑏(𝑘). Over 100 independent experiments, that is, training data set and test data are generated randomly in each experiment, we have found that the hyperparameters determined in this way achieve better BER performance and are more robust, comparing with the method with single start. The computational complexity is O(𝑁 𝑀 2 ) for training the SSGP algorithm. At prediction time, the cost is O(𝑀 ) for the predictive mean and O(𝑀 2 ) for the predictive variance per test point. In brief, the complexity of the SSGP is much lower than O(𝑁 3 ) when 𝑀 ≪ 𝑁 . We will show it is always the case for the MUD problem. Furthermore, the storage complexity of the SSGP is O(𝑁 𝑀 ), which is also lower than O(𝑁 2 ) if 𝑀 ≪ 𝑁. IV. E XPERIMENTAL R ESULTS AND D ISCUSSIONS We conduct a series of experiments to test the performance of the proposed SSGP based MUD algorithm. Consider a DSCDMA system with 8 users, the spreading codes are random binary sequences with length 31, the channel response is following, 𝐻(𝑧) = 0.4 + 0.9𝑧 −1 + 0.4𝑧 −2 ,
(19)
0
50
100
150
200
Number of Training Samples
250
300
Fig. 2. BER as a function of the number of training samples, SNR=4dB and 𝑀 = 20.
and the training and testing data are generated randomly. Firstly, we verify how many basis functions are appropriate for the SSGP to approach the full GP for the MUD problem. Fig.1 shows the BER as a function of the number of basis functions of the SSGP. The SNR is 4 dB and there are 𝑁 = 200 training samples. All users transmit signals with equal power. As seen from Fig.1, the BER varies in a narrow range for the SSGP when the number of basis functions is larger than 20. It is reasonable to set 𝑀 = 20 to approximate the full GP for the SSGP to solve the MUD problem. Fig.2 shows the BER of the GP, the SSGP and the MMSE [15], as a function of the number of training samples. The number of basis functions of the SSGP is 20. We can see that the GP and the SSGP perform better than the MMSE when the number of training samples is relatively small. It means that we can transmit more information bits in each burst of data with the SSGP or the GP based MUD scheme. It can be also seen that the SSGP and the GP have almost the same BER performance. Fig.3 shows the BER as a function of signal to noise ratios (SNRs) for the SSGP, the full GP and the MMSE. Again, it is shown that the SSGP produces a good approximation to the GP. To test near-far resistance performance of the SSGP, we consider a case where all users but the first one have the same bit energy 𝐸𝑖 . The SNR of the first user is 6dB. The number of basis functions of the SSGP is also 𝑀 = 20. Fig.4 shows the curves of BER versus 𝐸1 /𝐸𝑖 (dB). We can see the SSGP and the GP have almost the same performance, both of which perform better than the MMSE.
WANG and GU: MULTIUSER DETECTION WITH SPARSE SPECTRUM GAUSSIAN PROCESS REGRESSION
some intensive topics, such as which kind of basis functions is suitable for a specific problem and how many basis functions are necessary to obtain an acceptable approximation.
0
10
Full GP SSGP MMSE
−1
10
167
ACKNOWLEDGEMENT −2
10
BER
The authors would like to thank the anonymous reviewers, whose invaluable comments helped improve the presentation of this paper substantially. Shaowei Wang is grateful to Prof. Zhi-Hua Zhou for his helpful suggestions for this research.
−3
10
−4
10
R EFERENCES −5
10
Fig. 3.
−4
−2
0
2
4
SNR(dB)
6
8
10
12
BER as a function of SNR, 200 training samples, 𝑀 = 20. −1
10
BER
Full GP SSGP MMSE
−2
10
−3
10
−2
−1.5
−1
−0.5
0
0.5
Amplitude ratio(dB)
1
1.5
2
Fig. 4. Near-far resistance performance, 200 training samples, SNR=6dB and 𝑀 = 20.
V. C ONCLUSION Gaussian process (GP) for regression exhibits good performance with short training samples for the multiuser detection problem, but the complexity is too high. We adopt a sparse spectrum gaussian process algorithm to approximate the GP. The proposed algorithm has much lower complexity than the GP while keeping almost the same BER performance. The method given in this Letter throws some lights on designing more efficient digital communication receivers if considering
[1] S. Verdu, Multiuser Detection. Cambridge University Press, 1998. [2] B. Aazhang, B. P. Paris, and G. C. Orsak, “Neural networks for multiuser detection in code-division multiple-access communications,” IEEE Trans. Commun., vol. 40, no. 7, pp. 1212–1222, July 1992. [3] U. Mitra and H. V. Poor, “Neural network techniques for adaptive multiuser demodulation,” IEEE J. Sel. Areas Commun., vol. 12, no. 9, pp. 1460–1470, Dec. 1994. [4] B. Aazhang, B. P. Paris, and G. C. Orsak, “Neural networks for multiuser detection in code-division multiple-access communications,” IEEE Trans. Commun., vol. 40, no. 7, pp. 1212–1222, July 1992. [5] S. Chen, A. K. Samingan, and L. Hanzo, “Support vector machine multiuser receiver for DS-CDMA signals in multipath channels,” IEEE Trans. Neural Networks, vol. 12, no. 3, pp. 604–611, May 2001. [6] J. Murillo-Fuentes and F. Perez-Cruz, “Gaussian process regressors for multiuser detection in DS-CDMA systems,” IEEE Trans. Commun., vol. 57, no. 8, pp. 2339–2347, Aug. 2009. [7] J. Murillo-Fuentes, S. Caro, and F. Perez-Cruz, “Gaussian process for multiuser detection in CDMA receivers,” in Proc. 2005 Advances in Neural Information Processing Systems. [8] F. Perez-Cruz, J. Murillo-Fuentes, and S. Caro, “Nonlinear channel equalization with Gaussian processes for regression,” IEEE Trans. Signal Process., vol. 56, no. 10, pp. 5283–5286, Oct. 2008. [9] C. K. I. Williams, “Prediction with Gaussian processes: from linear regression to linear prediction and beyond,” Learning and Inference in Graphical Models, M. Jordan, editor. Kluwer Academic Press, 1998. [10] C. E. Rasmussen and C. K. I. Williams, Gaussian Process for Machine Learning. MIT Press, 2006. [11] M. Lazaro-Gredilla, J. Quinonero-Candela, C. Rasmussen, and A. Figueiras-Vidal, “Sparse spectrum Gaussian process regression,” J. Mach. Learn. Res., vol. 11, pp. 1865–1881, June 2010. [12] E. Snelson and Z. Ghahramani, “Sparse Gaussian processes using pseudo-inputs,” in Proc. 2006 Advances in Neural Information Processing Systems. [13] C. Walder, K. I. Kim, and B. Scholkopf, “Sparse multiscale Gaussian process regression,” in Proc. 2008 International Conference on Machine Learning. [14] M. Jordan, Learning in Graphical Models. MIT Press, 1998. [15] L. Rugini, P. Banelli, and S. Cacopardi, “A full-rank regularization technique for MMSE detection in multiuser CDMA systems,” IEEE Commun. Lett., vol. 9, no. 1, pp. 34–36, Jan. 2005. [16] A. Carlson, Communication Systems. McGraw-Hill, 1986. [17] M. Stein, Interpolation of Spatial Data. Springer Verlag, 1999.
