Gaussian processes regressors for complex proper signals in digital communications

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ This is an Accepted Manuscript of an article published by 2014 IEEE 8th Sensor Array and Multichannel Signal Processing Workshop (SAM) on 25/06/2014, available at: https://doi.org/10.1109/SAM.2014.6882359 “© 2014 IEEE. Personal use of this material is permitted. Permission from IEEEmust be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other Works” 1 Gaussian Processes Regressors for Complex Proper Signals in Digital Communications Rafael Boloix-Tortosa, F. Javier Pay´ an-Somet and Juan Jos´ e Murillo-Fuentes Dep. Teor´ ıa de la Se˜ nal y Comunicaciones Escuela Superior de Ingeniera. Universidad de Sevilla Camino de los Descubrimientos s/n, 41092, Sevilla, Spain Email: [email protected] Abstract In this paper we develop the complex-valued version of the Gaussian processes for regression (GPR) for proper complex signals. This tool has proved to be useful in the nonlinear detection in digital communications in real valued models. GPRs can be cast as nonlinear MMSE where hyperparameters can be tuned optimizing a marginal likelihood (ML). This feature allows for a flexible kernel that can easily adapt either to a linear or nonlinear solution. We introduce the complex-valued form of the GPR, and develop it for the proper complex case. We also deal with the optimization of the ML. Some experiments included illustrate the good performance of the proposal. I. INTRODUCTION In digital communications the optimal detection is usually nonlinear and unfeasible due to its computational complexity and the need of non-observable data. Hence approximate solutions are considered. The MMSE criteria is widely used, where the linear form is preferred if computational complexity is to be kept to a minimum. Gaussian processes for regression (GPR) have been proposed as nonlinear MMSE solutions that provide a good compromise between performance and the pair complexity plus training data needed [1], [2]. Their major advantage being the optimization of its hyperparameters through the marginal likelihood. These solutions have been developed for real valued models, assuming that an augmented input space with real and complex parts could be used to cope with the complex-valued model [3]. In [4] we may find a good description of the digital communication model deeply revised from the complex-valued point of view. Besides, results on the kernel design for the complex case have been reported recently [5], [6]. Based on these works, we propose a full description of the GPR for proper complex-valued systems. We develop the marginal likelihood needed in the optimization of the hyperparameters, and its derivation to optimize it. Some experimental results have been included to illustrate the good performance of this solution. April 12, 2015 DRAFT 2 II. SYSTEM MODEL AND DETECTION We propose a general model where we transmit Ksymbols every sample time, bt= [bt(1), bt(2), ..., bt(K)]>, and receive Nsignals arranged in the column-vector rt: rt=P·vt+nt,(1) where vt= [b> t,b> t−1. . . b> t−Ms+1]>and ntis an N-dimensional column-vector with additive white Gaussian noise, where Mstakes into account the memory or length of the channel response. Unless needed, subindex twill be omitted hereafter. As an example, in synchronous DS-CDMA the matrix Pin (1) summarizes the effect of the channel, the length- Nspreading codes and the different amplitudes for each of the Kusers, and vtis a KMs-dimensional vector including the transmitted bits (Ms= 2 if 2≤Mc≤K, where Mcis the length of the channel in chip periods) [7], [1]. Many other communication systems can be modeled as (1) aswell. At the receiver, we learn the system model using training data. In wireless communications short training sequences are mandatory to increase the number of information bits. After the training stage, for every new received vector, denoted as r∗, a detector estimates the j-th transmitted symbol as ˆ b∗(j) = f(r∗)for some j, where this index will be omitted hereafter. Linear detectors work well when the inter-symbolic interference (ISI) and the users interference are negligible. In the general case, the solution is highly nonlinear and nonlinear detectors are necessary to get a good estimation. In digital communications, linear MMSE estimators are typically used, due to their simplicity and nearly optimal results in many standard scenarios. The linear, f(r) = r>w, MMSE detector can be expressed as a Wiener filter [8]: wmmse = argmin w Ehb−r>w2i= (Crr )−1Crb,(2) where Crr and Crbare, respectively, the autocorrelation of the inputs and the cross-correlation between the inputs and outputs. This solution can be estimated from training data by replacing Crr and Crbby their sampled versions: wsmmse = (R·R>)−1R·b,(3) where the columns of Rare the received signals at every symbol in the training set. We denote this solution as linear sampled MMSE (SMMSE). However, in the general case, the MMSE is not linear. It is widely known that the MMSE solution for AWGN channels is given by the conditional mean of bgiven r: ˆ b∗=fmmse(r∗) = E[b∗|r∗],(4) which is a nonlinear function of the inputs, unless the inputs are Gaussian distributed, that it is not the case. The computation of the optimal solution is hard to accomplish, given the great number of possible combinations of transmitted symbols due to the number of points in the transmitted constellations, the number of users and the channel memory. In [1], [2] the real valued GPR is proposed as a nonlinear tool to estimate the MMSE response. April 12, 2015 DRAFT 3 However, to our knowledge, there is no formulation available for complex-valued systems. In the following we develop a complete formulation of the GPR and hyperparameters adaptation for the complex proper case, which is the main contribution of this work. III. COMPLEX GAUSSIAN PROCESSES A complex Gaussian process f(x)is a collection of complex random variables, any finite number of which have a joint complex Gaussian distribution. The general pdf of a complex Gaussian random vector z: Ω →CDis [4] p(z) = 1 πndet1/2Rzz exp −1 2(z−µz)HR−1 zz (z−µz),(5) where the augmented vector zis obtained by stacking zon top of its complex conjugate z∗. Similarly, the augmented mean vector µzis obtained by stacking vector µzon top of its complex conjugate µ∗ z.Rzz is the augmented covariance matrix of z Rzz =  Rzz ˜ Rzz ˜ R∗ zz R∗ zz  ,(6) where Rzz =E(z−µz)(z−µz)His the usual covariance matrix and ˜ Rzz =E(z−µz)(z−µz)>is the complementary covariance or pseudo-covariance matrix. It is important to note that both Rzz and ˜ Rzz are required for a complete second-order characterization of the complex Gaussian random vector z∼ N µz,Rzz,˜ Rzz. A complex Gaussian process f(x)is completely specified by its augmented mean function, µ(x), its covariance function, k(xi,xj), and its complementary covariance function, ˜ k(xi,xj), i.e. f(x)∼ GP(µ(x), k(xi,xj),˜ k(xi,xj)). There is an important special case in which the complementary covariance for a complex Gaussian vector vanishes. Definition 3.1: [4] If the complementary covariance matrix vanishes, ˜ Rzz = 0,zis called proper, otherwise it is called improper. The pdf (5) of a complex proper Gaussian random vector z: Ω →CDsimplifies to p(z) = 1 πndet Rzz exp −(z−µz)HR−1 zz (z−µz).(7) We can propose now the following definition: Definition 3.2: If the complementary covariance function of a complex Gaussian process vanishes, the process is called proper, otherwise it is called improper. IV. COMPLEX PROPER GAUSSIAN PROCESSES FOR REGRESSION We consider the model f(x) = φ(x)Hw,(8) where xis the complex input vector of dimension Dand we specify a prior over the parameters w: a zero mean proper complex Gaussian with covariance matrix Σp, i.e., w∼ N (0,Σp,0). This linear regression model gives an example of a zero mean complex Gaussian process, µ(x) = φ(x)HE[w] = 0, with covariance function k(xi,xj) = φ(xi)HEwwHφ(xj) = φ(xi)HΣpφ(xj) = ϕ(xi)Hϕ(xj). The complementary covariance function April 12, 2015 DRAFT 4 vanishes, ˜ k(xi,xj) = φ(xi)HEww>φ∗(xj)=0, i.e., f(x)is a complex proper Gaussian process: f(x)∼ GP(0, k(xi,xj),0). We can specify a covariance function, take ninput vectors xi,i= 1, ..., n aggregated column-wise in the D×n matrix X, and write out the corresponding n×ncovariance matrix K(X,X)elementwise from k(xi,xj). For those inputs we generate a complex proper Gaussian random vector f∼ N (0,K(X,X),0). In the more realistic situation y=f(x) + =φ(x)Hw+, (9) the observed values ydiffer from f(x)by additive complex proper Gaussian noise with zero mean and variance σ2 n, i.e., ∼ N 0, σ2 n,0. Hence, the complex random vector ythat results from the inputs in Xis also proper Gaussian: y∼ N (0,C,0), with C=K(X,X)+σ2 nI. Furthermore, yand fare cross-proper, as ˜ Ryf =Eyf>= Eff>+Ef>= 0. Hence, yand fare jointly proper [4], i.e., the composite complex random vector y>,f>> is Gaussian proper. If we take now a number of test input vectors aggregated column-wise in matrix X?, the joint distribution of the training outputs yand the test outputs f?is   y f? ∼ N  0,  C K(X,X?) K(X?,X)K(X?,X?) ,0 ,(10) and the conditional distribution of f?given yis f?|x?,X,y∼ N ¯ f?,cov(f?),0,(11) where ¯ f?=K(X?,X)C−1y,(12) cov(f?) = K(X?,X?)−K(X?,X)C−1K(X,X?).(13) In the case that there is only one test input vector x?, ¯ f?=k> ?C−1y,(14) cov(f?) = k(x?,x?)−k> ?C−1k?,(15) where k(x?) = k?is the vector of covariances between the test point and the ntraining points. A. The marginal likelihood The marginal likelihood (or evidence) is the integral of the likelihood times the prior p(y|X) = Zp(y|f,X)p(f|X)df,(16) Here, the prior is complex proper Gaussian, f|X∼ N (0,K(X,X),0). Since yand fare jointly proper, the likelihood is the factorized proper Gaussian, y|f∼ N f, σ2 nI,0. Now, using (7) in (16) the log marginal likelihood yields L= log p(y|X) = −yHC−1y−log det C−nlog π, (17) i.e., y|X∼ N (0,C,0). April 12, 2015 DRAFT 5 B. Model selection for complex proper GP regression In order to turn complex Gaussian processes into powerful practical tools it is essential to develop methods that address the model selection problem. It is common to use a hierarchical specification of models, with the parameters, wat the first level and some hyperparameters θiat the second level. In our model, the covariance function k(xi,xj)can be parameterized in terms of the hyperparameters θi. It is a property of the marginal likelihood that it automatically incorporates a trade-off between model fit and model complexity [9]. This is the reason why the marginal likelihood is valuable in solving the model selection problem. The log marginal likelihood in (17) can be explicitly written conditioned on the hyperparameters, L(θ), with C=C(θ) = K(X,X) + σ2 nI, where θis the vector of hyperparameters. We can maximize (17) in order to set those hyperparameters. Notice that the variance of the additive noise, σ2 n, is also treated as another hyperparameter and can be estimated from this optimization. L(θ)in (17) is a real function of a complex Hermitian matrix C(θ)that depends on the hyperparameters. The chain rule can be used to seek the derivative of the log marginal likelihood with respect to C(θ)and then the derivative of C(θ)with respect to the hyperparameters. However, the Cauchy-Riemann conditions dictate that nonconstant real valued functions that are defined on complex domains are not holomorphic, so the complex derivative cannot be used. Instead, we must seek formal or Wirtinger derivatives [5], [10]. The chain rule for the Wirtinger derivative of the function L(θ)with respect to the hyperparameter θiyields [10] ∂L ∂θi = n X k=1 n X l=1 ∂L ∂Clk ∂Clk ∂θi +∂L ∂C∗ lk ∂C∗ lk ∂θi,(18) where Clk is the element (l, k)of matrix C. When finding the Wirtinger derivatives ∂L ∂Clk and ∂L ∂C∗ lk the elements Clk and C∗ lk are treated as independent variables. Therefore, ∂L ∂θi = n X k=1 n X l=1 ∂L ∂Clk ∂Clk ∂θi = Tr "∂L ∂C>∂C ∂θi#.(19) The derivative of the first term of the function Lin (17) with respect to Cyields ∂ ∂C−yHC−1y=−∂ ∂CTr yHC−1y =C−> yyH>C−>.(20) The derivative of the second term of the function Lin (17) with respect to Cyields ∂ ∂C(−log det C) = −C−>.(21) Substitution of (20) and (21) in (19) yield ∂L ∂θi = Tr C−1yyHC−1−C−1∂C ∂θi.(22) April 12, 2015 DRAFT 6 C. Complex covariance functions We recall here the complex kernel function used in the Experiments section of this study, the independent Gaussian kernel [6]. This kernel is defined as κ(xi,xj) = κR(αi, αj) + κR(βi, βj)) +j(κR(αi, βj)−κR(βi, αj)) ,(23) where xl=αl+jβland κRis the well-known real Gaussian kernel. This kernel is a generic extension of κRand inherits the intuitive physical meaning of a measure of similarity of the samples that κRhas. V. EXPERIMENTS By replacing rand bin (1) by xand y, respectively, in (9) we yield the complex GPR model to estimate b. That estimate can be found as ¯ f?in (14) for each new received vector x?=r∗. In these experiments we focus on different scenarios to show the good behavior of this detector. We consider here QPSK modulations. Note that they are proper complex-valued and so are the received signals in (1). In Fig. 1 we have included the BER for the complex GPR with the kernel in (23) for the experiment in Section V.A in [1]. In this experiment we have a DS-CDMA system with K= 2 equal power users and signatures [+1+1−1−1] and [+1 −1−1 + 1]. The channel is given by c(z) = 0.3+0.7z−1+ 0.3z−2. The four received chips after the two matched filters are the inputs, x, i.e., N= 2. The number of training samples used was n= 64. In Fig. 2 the decision boundary is included for the complex GPR with 64 training samples and SNR = 16 dB, along with the points corresponding to the bits of the real part of the inputs for both users. This experiment illustrates a scenario where the GPR is of major interest, since the linear response fails to detect the transmitted symbols. The complex GPR successfully computes the nonlinear solution, even for the short number of training samples used. In Fig. 3 we include a DS-CDMA scenario with K= 10 equal power users and signatures with N= 31 chips and a channel of 10 zero mean complex Gaussian iid taps. We depict the BER along the number of training samples, n, for a SNR = 12 dB. In this experiment, the inputs to the detector are N= 31 received chips for one symbol period. Since the large dimension, the solution is nearly linear and the linear SMMSE exhibits a good performance. The GPR with the nonlinear complex kernel adapts well to this linear solution, even for a short number of samples. In the short length scenario the regularization term learned is high to compensate for the lack of information, and the GPR performs better than the SMMSE. As the number of training samples increases, the nonlinear part of the kernel learns the linear solution and the GPR is able to follow the SMMSE. The slight mismatch between both curves is again due to the regularization term. VI. CONCLUSION In this paper we develop a complete formulation of the GPR and hyperparameters adaptation for the complex proper case. Complex proper GPR can adapt either to linear or nonlinear solutions with short training sequences, therefore they are of interest as detectors in communication applications, as shown in the included experiments. April 12, 2015 DRAFT 7 4 6 8 10 12 14 16 18 20 10−5 10−4 10−3 10−2 10−1 SNR(dB) BER MF SMMSE GPR Fig. 1. BER along the SNR in a DS-CDMA scenario with K= 2 users and signatures [+1 + 1 −1−1] and [+1 −1−1 + 1] for SMMSE (◦), the matched filter () and GPR () detectors for n= 64 training samples. −1 0 1 −1 −0.5 0 0.5 1 <(x(1)) <(x(2)) Fig. 2. Decision boundary for the GPR in a DS-CDMA scenario with K= 2 users and signatures [+1 + 1 −1−1] and [+1 −1−1 + 1] for n= 64 training samples and SNR = 16 dB. ACKNOWLEDGMENT This work was partially funded by Spanish government (Ministerio de Educaci´ on y Ciencia, Consolider-Ingenio 2010 CSD2008-00010, TEC2012-38800-C03-{02}) and European Union (FEDER). April 12, 2015 DRAFT 8 50 100 150 200 250 10−5 10−4 10−3 10−2 10−1 n BER MF SMMSE GPR Fig. 3. BER along n, in a DS-CDMA scenario with K= 10 users and signatures with N= 31 for SMMSE (◦), the matched filter () and GPR () detectors for SNR = 12 dB. REFERENCES [1] J. Murillo-Fuentes and F. P´ erez-Cruz, “Gaussian process regressors for multiuser detection in DS-CDMA systems,” IEEE Transactions on Communications, vol. 57, no. 8, pp. 2339–2347, 2009. [2] F. P´ erez-Cruz, S. Van Vaerenbergh, J. Murillo-Fuentes, M. Lazaro-Gredilla, and I. Santamaria, “Gaussian processes for nonlinear signal processing: An overview of recent advances,” Signal Processing Magazine, IEEE, vol. 30, no. 4, pp. 40–50, July 2013. [3] F. P´ erez-Cruz, J. Murillo-Fuentes, and S. Caro, “Nonlinear channel equalization with Gaussian Processes for Regression,” IEEE Transactions on Signal Processing, vol. 56, no. 10, pp. 5283–5286, Oct. 2008. [4] P. J. Schreier and L. L. Scharf, Statistical Signal Processing of Complex-Valued Data. The Theory of Improper and Noncircular Signals. Cambridge, UK: Cambridge University Press, 2010. [5] P. Bouboulis and S. Theodoridis, “Extension of Wirtinger’s calculus to reproducing kernel Hilbert spaces and the complex kernel LMS,” Signal Processing, IEEE Transactions on, vol. 59, no. 3, pp. 964–978, 2011. [6] F. Tobar, A. Kuh, and D. Mandic, “A novel augmented complex valued kernel LMS,” in Sensor Array and Multichannel Signal Processing Workshop (SAM), 2012 IEEE 7th, 2012, pp. 473–476. [7] S. Chen, A. K. Samingan, and L. Hanzo, “Support vector machine multiuser receiver for DS-CDMA signals in multipath channels,” IEEE Transactions on Neural Network, vol. 12, no. 3, pp. 604–611, December 2001. [8] S. Verd´ u, Multiuser Detection. Cambridge University Press, 1998. [9] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning. Massachusetts Institute of Technology, Cambridge, Massachusetts: MIT Press, 2006. [10] A. Hjorungnes and D. Gesbert, “Complex-valued matrix differentiation: Techniques and key results,” Signal Processing, IEEE Transactions on, vol. 55, no. 6, pp. 2740–2746, June 2007. April 12, 2015 DRAFT