Unsupervised Common Spatial Patterns
Abstract
The common spatial pattern (CSP) method is a dimensionality reduction technique widely used in brain-computer interface (BCI) systems. In the two-class CSP problem, training data are linearly projected onto direc tions maximizing or minimizing the variance ratio between the two classes. The present contribution proves that kurto sis maximization performs CSP in an unsupervised manner, i.e., with no need for labeled data, when the classes follow Gaussian or elliptically symmetric distributions. Numerical analyses on synthetic and real data validate these findings in various experimental conditions, and demonstrate the interest of the proposed unsupervised approach.
Full text
IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 27, NO. 10, OCTOBER 2019 2135 Unsupervised Common Spatial Patterns Rubén Martín-Clemente , Member, IEEE , Javier Olias , Sergio Cruces , Senior Member, IEEE , and Vicente Zarzoso , Senior Member, IEEE Abstract —The common spatial pattern (CSP) method is a dimensionality reduction technique widely used in brain-computer interface (BCI) systems. In the two-class CSP problem, training data are linearly projected onto directions maximizing or minimizing the variance ratio between the two classes. The present contribution proves that kurtosis maximization performs CSP in an unsupervised manner, i.e., with no need for labeled data, when the classes follow Gaussian or elliptically symmetric distributions. Numerical analyses on synthetic and real data validate these findings in various experimental conditions, and demonstrate the interest of the proposed unsupervised approach. Index Terms —Common spatial patterns, brain computer interfaces, kurtosis. I. INTRODUCTION COMMON spatial patterns (CSP) is a dimension reduction technique widely used in brain-computer interface (BCI) systems [1]–[4]. Typically, electroencephalogram (EEG) samples acquired under two different experimental conditions provide a multivariate data set with two classes. CSP linearly projects the data onto directions where the variance of the projected data points is significantly higher for one class than for the other [5]–[7]. The projected data variances can then be used as features for classification. CSP is a supervised technique, whose performance relies heavily on the availability of correctly labeled data. The present contribution proves that CSP can be also performed in an unsupervised fashion by maximizing the kurtosis (normalized fourth-order moment) of the projected data. Unsupervised operation spares the need for training labels and is thus immune to erroneous labelling. Apart from its theoretical interest, this result is useful, for instance, in applications where the training labels are not available or may be uncertain. A mathematical proof is derived for data drawn from a mixture of Gaussian densities, and then Manuscript received May 8, 2019; revised July 28, 2019; accepted August 5, 2019. Date of publication September 19, 2019; date of current version October 8, 2019. This work was supported in part by the Spanish Ministry of Economy and Competitiveness under Project TEC2017-82807-P and in part by the European Regional Development Fund (ERFD) of the European Commission. The work of R. Martín-Clemente was supported by the I3S Laboratory through the IFA Scholarship from Université Côte d’Azur in July 2018. Most of this work was carried out while V. Zarzoso was a member of the Institut Universitaire de France . (Corresponding author: Rubén Martín-Clemente.) R. Martín-Clemente, J. Olias, and S. Cruces are with the Department of Signal Theory and Communications, University of Seville, 41092 Seville, Spain (e-mail: [email protected]; [email protected]; [email protected]). V. Zarzoso is with Université Côte d’Azur, CNRS, I3S Laboratory, 06903 Sophia Antipolis, France (e-mail: [email protected]). Digital Object Identifier 10.1109/TNSRE.2019.2936411 generalized to elliptically symmetric distributions. Our experimental evaluation on synthetic and real data corroborates the theoretical findings. Unsupervised techniques are not unknown in EEG processing: e.g., [8] shows that it is possible to perform unsupervised workload classification using EEG spectral features. We can expect that they will become increasingly common in the near future. The paper is organized as follows: Section II reviews the mathematical formulation of the common spatial patterns method. Section III, the core of our contribution, establishes the link between the kurtosis and the CSP criterion under the assumption of a Gaussian mixture model for the data. Section IV extends this result to elliptically distributed classes. Illustrative examples supporting the theoretical derivations are presented and discussed in Section V. The concluding remarks of Section VI bring the paper to an end. For the sake of clarity, proofs of the theoretical results have been deferred to the Appendices. II. COMMON SPATIAL PATTERNS Consider that we are given a set of observations of a random variable Xin Rp, in which each observation belongs to one of two classes C1and C2. CSP is usually applied to problems where the class means are null, andthisassumptionismade in the sequel. A one-dimensional projection of the point cloud can be represented by Y aX,wherea Rp. Denoting by k Cov X Ck the covariance matrix of the data in class Ck, with 1 2, it holds that the variance of class Ck,afterthe projection, equals σ2 k aka,k 1,2.(1) The idea behind CSP is to maximize σ2 1while minimizing σ2 2 or vice versa [9]. To this end, the objective function is defined as the power ratio R a σ2 1 σ2 2 .(2) Note that σ2 k,k 1,2, depend on athrough relation (1). Also, ratio (2) is scale invariant, i.e., R a R ca ,forall c R 0 , and therefore only the direction of the projection is significant but not the overall scaling. To find the optimal a corresponding to the extremum (maximum or minimum) of CSP criterion (2) we set its gradient ∇R a to zero, readily yielding 1a R a 2a .(3) It follows that a is a generalized eigenvector of 1 and 2[10]. Solving this problem we get the eigenvector This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/
2136 IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 27, NO. 10, OCTOBER 2019 corresponding to the maximum (respectively, minimum) eigenvalue, which maximizes (resp. minimizes) the CSP objective function. When more features are needed for the posterior classification stage, a common practice is to project the data onto other eigenvectors from both ends of the eigenvalue spectrum [11]. In the generalized eigenvalue (GEVD) problem (3), aidenote the eigenvectors and their corresponding eigenvalues are assumed to be sorted in decreasing order: R ai R aj , for i j,i 1,2,...,p. It is important to remark that a training set of correctly classified observations is required to estimate matrices k, k 1,2. It is in this sense that CSP can be considered a supervised technique. The remaining of this paper presents a fully data-driven procedure that does not require labeled samples, thus performing CSP in an unsupervised manner. III. KURTOSIS AS A BLIND CSP CRITERION We first focus on the important case where Xis distributed as a mixture of two Gaussian densities [11]: X π1N 0,1 π2N 0,2 (4) where πkstands for the prior probability of class Ck,k 1,2, with π1 π2 1, and 0is a p-dimensional vector of zeros. The distribution of the one-dimensional projection Y aX is also a mixture of Gaussians, that is, Y π1N 0,σ2 1 π2N 0,σ2 2 ,whereσ2 1and σ2 2are defined as in eqn. (1). From the properties of the Gaussian distribution, it follows that E Y2 π1σ2 1 π2σ2 2,E Y4 3π1σ4 1 3π2σ4 2 where E denotes the mathematical expectation operator. Now, the kurtosis is a statistic defined as the normalized fourth-order moment [12] κY a E Y4 E Y2 2.(5a) Under data model (4), the kurtosis can be expressed as: κY a 3π1σ4 1 3π2σ4 2 π1σ2 1 π2σ2 2 2.(5b) Some preliminary manipulations show that ∇κ a 0if and only if we have R a 1 ∇R a 0, meaning that the critical points of R a are also critical points of κ a . Indeed, a thorough analysis detailed in Appendix A leads to the following result: Theorem 1: Under the working assumptions of CSP recalled in Sec. II, the local maximizers of the kurtosis (5b) maximize or minimize the CSP criterion (2). In particular, the maximizers of kurtosis are a1if R a1 1andapif R ap 1. In other words, CSP can be performed blindy, that is, without the need for labeled samples, by projecting the observed data points onto the direction that maximizes the kurtosis of the projections. Consequently, the new method will be referred to as kurtosis-based unsupervised CSP (k-uCSP). Under different working assumptions, namely that the data are a combination of statistically independent variables, the optimization of the kurtosis gives rise to independent component Fig. 1. The kurtosis versus the CSP objective function for π1 π2 1/2 with 0.5 R a p R a R a 1 4. analysis (ICA) [13], [14]. Furthermore, in clustering problems, maximizing the kurtosis can produce the same results as Fisher linear discriminant analysis [15]. We see that the kurtosis is a widely-used tool in signal processing. It is the true model the data follow that determines the outcome achieved. It should be remarked that the expression of kurtosis given in eqn. (5b) is only exploited to prove the equivalence between kurtosis optimization and common spatial pattern analysis in our theoretical derivations, but we use eqn. (5a) in the actual algorithm to find the patterns. Therefore, the proposed technique is fully unsupervised. Finally, Theorem 1 also admits an intuitive interpretation. Dividing both the numerator and denominator of (5b) by σ4 2, the right-hand part of this formula can be expressed in terms of R a , κY a 3π1R2 a 3π2 π1R a π2 2,(6) thus revealing that there exists a relationship between the supervised and unsupervised criteria. Fig. 1 shows an example plot of κYagainst Rwhen π1 π2and Ris in the range 0.5,4 . We make the observation that the maximizers of the kurtosis can lie only in the strictly increasing/decreasing range of (6). Consequently, (6) is also increasing/decreasing in some neighborhood of them and, therefore, invertible. Then, a decrease in the value of the kurtosis in that neighborhood results in a decrease/increase in the value of R.Asitis intuitive, this implies that the local maximizers of the kurtosis maximize or minimize the CSP criterion. In Figure 1, additionally, it is not hard to show that as (6) is decreasing for all vectors sufficiently near ap, then (6) has a maximum at ap. Virtually the same argument shows that there is another maximum at a1. IV. EXTENSION TO ELLIPTIC DISTRIBUTIONS The above result can be extended to the case where the data follow a mixture of zero-mean elliptically symmetric distributions [16]. These distributions generalize the class
MARTÍN-CLEMENTE et al. : UNSUPERVISED CSP 2137 of multivariate Gaussians by allowing for both heavier-than- Gaussian and lighter-than-Gaussian distribution tails. Examples include the tand Laplace distributions. Before continuing, a zero-mean p-dimensional random variable Zis elliptically distributed if its characteristic function ϕZ a E ejaZ , j 1, a Rp, is of the form ϕZ a φ 1 2aa for some nonnegative-definite matrix . Function φ is called characteristic generator of the distribution. For example, for Gaussian variables φ α exp α and for multivariate Laplacian variables φ α 1 1 α . The characteristic function of the univariate random variable W aZis given by ϕW t E ejtW E ejtaZ ϕZ ta φ 1 2t2aa , t R. It follows that, if the variance of Wexists, then σ2 E W2 2 t2ϕW t t 0 φ 0 aa(7) E W4 4 t4ϕW t t 0 3γσ4(8) wherewehavedefinedγ φ ¾ 0 φ ½ 0 2, which is a strictly positive number. Notations φ and φ represent, respectively, the first- and second-order derivative of φ. Now, let us replace the Gaussian distributions in (4) with zero-mean elliptically symmetric distributions having characteristic generators of the same type φ but defined by matrices 1and 2, i.e., φ a1a 2 and φ a2a 2 , respectively. Then, the moments of Ybecome E Y2 π1σ2 1 π2σ2 2,E Y4 3π1γσ4 1 3π2γσ4 2 with σ2 k,k 1,2, given by (7), and their ratio turns out to be a positively scaled version of (5b): κY a γ3π1σ4 1 3π2σ4 2 π1σ2 1 π2σ2 2 2.(9) It readily follows that Theorem 1 also holds for data classes with elliptic distributions defined by the same type of characteristic generator φ . V. EXPERIMENTAL ASSESSMENT A number of experiments are performed to validate the theoretical study of the unsupervised CSP criterion developed in this paper and to test its performance in a variety of experimental conditions. These include synthetically generated as well as real EEG data. To perform the numerical optimization of the kurtosis statistic (5a), we employ the algorithm presented in [24]. We point out that this algorithm does not require any class labels as inputs. A free MATLAB implementation of the algorithm is provided in [25]. A. Simulated Data To illustrate Theorem 1, let us first consider a mixture of equiprobable classes in a two-dimensional space, i.e., p 2. We consider three cases that only differ in the choice of matrix 2. Case 1 assumes that 1 diag 2,1 and 2 diag 0.5,2 ; in case 2, 2is replaced by diag 0.2,0.8 ; in case 3, finally, we set 2 diag 2.5,10 . In all cases, the generalized eigenvector of 1and 2that maximizes the CSP Fig. 2. Criteria R a (dotted) and κ Y a (solid) as a function of angle θ defining the projection direction a θ , θ . Red curves are calculated from Gaussian data, while black curves correspond to Laplacian classes. Dotted curves ( R ) overlap as the pair 1,2 is the same for both distributions. As predicted by Theorem 1, the local maximizers of kurtosis (solid lines) either maximize or minimize the CSP criterion (dotted lines) depending on the generalized eigenvalues R a 1 and R a 2 of 1and 2. Statistics are estimated from T 1000 random samples, with 500 i.i.d. samples per class, and each curve is the average of 100 independent experiments. target function (2) is a1 1,0 , while the corresponding minimizer is a2 0,1 . Figure 2 plots the CSP criterion R[eqn. (2)] and the kurtosis criterion κY[computed from the ratio of expectations in eqn. (5b)] in dotted and solid lines, respectively, for
2138 IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 27, NO. 10, OCTOBER 2019 a cos θ ,sin θ . Red and black solid lines correspond to Gaussian and Laplacian data, respectively, where the latter were generated as explained in [17]. The results are just as predicted by Theorem 1: the maxima of κYare always maxima or minima of the CSP criterion R, with different patterns of correspondence depending on the generalized eigenvalues of 1and 2. In case 1 (Figure 2a), the maxima and minima of Rtransform into maxima of κY,asR a1 4 1 R a2 0.5. In case 2 (Figure 2b), κYhas the same maxima and minima as R, because both eigenvalues are greater than one: R a1 10, R a2 1.25. Finally, in case 3 (Figure 2c), the minima of Rare transformed into maxima of κYand vice versa, since both eigenvalues are lower than one: R a1 0.8, R a2 0.1. Additionally, the point R 1, reached in case 1, defines a global minimum of the kurtosis. As an additional validation experiment, let us also test the case where samples from the same class are correlated, resulting in non diagonal covariance matrices. The following matrices are selected at random: 1 3.8152 3.4131 3.4131 3.3104 ,2 2.8465 0.5267 0.5267 1.2446 . For these covariances, the maximizer a1of Ris the unit vector that makes an angle of θ1 2π 3 radians with the positive x-axis. Similarly, the minimizer a2is at an angle θ2 π 4radians.Figure 3 (top) plots Rand κYin a format similar to that of Figure 2. An alternative representation is given in Figure 3 (bottom). The color-coded circles indicate the kurtosis of the entire projected data (outer circle) and the variance ratio of the projected classes (inner circle) as a function of angle θ. For the reader’s convenience, we complete the figure by drawing a pair of dashed lines in the direction of vectors a1and a2. Additionally, we draw in different colours the scatterplots of the classes and the probability density functions of the data points after being projected onto a1and a2, which illustrates well the disparity between the variances resulting from the projection. Furthermore, it is known that whitening the data reduces the initial generalized eigenvalue problem (3) to a standard eigenvalue problem for Hermitian matrices (see Appendix B), for which the eigenvectors are orthogonal. To show this property, we whiten the previous data before calculating Rand κY.Asa result, Fig. 3, obtained before whitening, transforms into Fig. 4 after whitening. As expected, Fig. 4 (bottom) shows that the directions that maximize the kurtosis become perpendicular. Yet the equivalence between the CSP and kurtosis criteria established by Theorem 1 still holds in this case. B. Real Data 1) Supervised vs. Unsupervised CSP of Brain Data: In typical BCI implementations, subjects are instructed to imagine movements of, e.g., their right or left hand, one after the other, while their p-channel EEG is being recorded. Then, CSP is applied to the data, previously bandpass-filtered in the band of interest, and the power evolution of the resulting time series allows the BCI system to discriminate between the two classes of motor imagery [18]–[20]. Fig. 3. Supervised vs. unsupervised CSP when samples from the same class are correlated, with R a 1 5.4 1 R a 2 0.06. (Top) Criteria R a (dotted) and κ Y a (solid) as a function of angle θdefining the projection direction a θ , θ .AsinFig. 2, red curves are calculated from random Gaussian samples, and black curves from Laplacian data. (Bottom) Alternative representation asa scatterplot of the two correlated Gaussiandistributions. Samples from one class areshown as red dots, and from the other as blue dots. The two color-coded circles around the scatterplots indicate thevalue of R a (inner circle) and κ Y a (outer circle). Straight dashed lines mark the directions representing the CSP and kurtosis projections, which coincide in this experiment. The probability density functions of the projected classes are also shown in red and blue lines. The purpose of this experiment is to illustrate the performance of the k-uCSP approach on real EEG data. The proposed kurtosis-based approach is tested using datasets from the BCI competition IV (dataset 2a) [21]. These contain EEG acquired on two different days from nine subjects with a p 22-channel EEG system at a sampling rate of 250 Hz. As pre-processing, we bandpass filter the EEG to 8 30 Hz. This is usual in BCI and ensures that the data are zeromean. Electrooculogram (EOG) channels are available and ocular artifacts are removed using Signal Space Projections (SSP) [22], [23].
MARTÍN-CLEMENTE et al. : UNSUPERVISED CSP 2139 Fig. 4. Results for the data in Fig. 3 after whitening. Apart from the whitening operation, details about the plots are as in the caption of Fig. 3. In each trial, an arrow pointing either to left, right, down or up is shown on a display for a short time, and the subject is required to respectively imagine left hand, right hand, both feet or tongue movements in response. The imagined action lasts about three seconds but only the final two are kept to avoid initial transient effects. Thus, each trial data matrix consists of 22 rows of channels and 2 (seconds) 250 (samples/seconds) 500 columns of time samples. Data matrices are also normalized as suggested in [26], which generalizes the procedure in [27]. The test statistic proposed in [28] shows that about 2 3 of the data follow an elliptical distribution with a standard level of significance of α 0.01. Finally, 72 trials of each movement are performed per subject and day. Experiments are conducted on groups of trials. Each group comprises 144 trials, half from one of the four classes (e.g., left hand) and half from another (e.g., right hand), all from the same subject and recorded the same day. The corresponding 144 trial data matrices are concatenated along the temporal Fig. 5. One-dimensional projection of EEG data with maximal kurtosis (blueline). Several projected trials from the first user are shown. Labels ‘1’ and ‘2’ correspond to data classes (‘1’ left hand, ‘2’ right hand). Visible to the naked eye, the power of one of the projected classes (namely, class ‘1’) is significantly higher, as would be expected from the application of CSP. For comparison, the projection of the same trials onto the corresponding supervised CSP filter are also shown (bottom black curves). dimension, producing one large observation matrix with dimension 22 72000. In total, there are 9 (subjects) 6 (possible combinations of four classes) 2 (days) 108 groups. The projections of the p-channel EEG onto the p-directions maximizing the kurtosis are computed for each group of trials by the algorithm in [24], [25]. Specifically, after applying a whitening transformation, we maximize p-times the kurtosis, one after the other, under the constraint that the direction obtained in the nth maximization is orthogonal to the previously calculated directions (see Appendix B for details). For illustration purposes, Figure 5 shows several projected trials with maximum kurtosis, calculated from ‘left-hand’ and ‘right-hand’ EEG data recorded from subject 1 during the second session. In this particular experiment, the power of the left-hand motor imagery projections (denoted ‘1’) is clearly higher than that of the other hand (denoted ‘2’), which facilitates class discrimination. A scatter plot of log-variances for all projected trials is also shown in Figure 6,usingthe supervised projection on the x-axis and the unsupervised kurtosis-based projection on the y-axis. As a more formal test of the relationship between both criteria (supervised and unsupervised), we also project the EEG data onto the p 22 generalized eigenvectors of the matrices 1and 2. For each projection direction, we calculate the ratio of variances of the projected classes (σ2 1 σ2 2). These ratios are then arranged in decreasing order along a p-dimensional vector bCSP. A distinct vector is generated for each group of trials, as the projection directions differ from one another. Next, we repeat the experiment with the difference that this second time the data are projected onto the pdirections maximizing the kurtosis. The new variance ratios are stored in a vector bKurt. A strong similarity is found between the vectors bCSP and bKurt calculated from the same data. The Pearson correlation coefficient between the components of these two vectors, averaged for all movements, is 0.9692 for subject 6 (best), and 0.9484 for subject 9 (worst).
2140 IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 27, NO. 10, OCTOBER 2019 Fig. 6. Scatter plot of log-variances for all ‘left-hand’ vs ‘right-hand’ trials recorded the second day from subject 1, using the supervised projection which maximizes R on the horizontal axis and the corresponding unsupervised projection with maximal kurtosis on the vertical axis. Legend: blue circles left-hand, red crosses right-hand. The correlation coefficient averaged over all groups of trials is 0.9586 with standard deviation 0.0257. Both approaches therefore provide very similar variance ratios, though one method is supervised while the other is not. 2) Unsupervised Discrimination: The k-uCSP approach may be combined with some unsupervised classification technique in order to design a full data-driven BCI system. Some preliminary illustrative experiments are presented in the remaining of this Section. Let Y1,...,Ypbe the unsupervised projections that result from projecting a single trial data matrix. These matrices are 22 500, so that each Yiconsists of 500 values. As in [27], their variances are used to calculate the following p 22 features Fi log var Yi p n 1var Yn ,i 1,...,p,(10) which are arranged in a p-vector. For illustrative purposes only, Figure 7 draws the scatter-plot of 80 log-variances Ficomputed from ‘left vs right-hand’ trials from subject 1. Superimposed, we draw the contour plot of the best fitting two-component Gaussian mixture model (GMM) of the point cloud [29]. By assigning a point to the Gaussian distribution it most probably belongs to, we partition the space in two regions with the hope that these regions accurately reflect the original classes. The performance of the classification can be evaluated by using the true labels as ground truth. For example, in Figure 7, when the subject thinks of a ‘left-hand’ movement, 27 times out of 40 (about 2 times out of 3) it is assigned to Region 1. Similarly, ‘right-hand’ movements are 36 times of 40 (9 times out of 10) assigned to Region 2. As both kind of movements are equiprobable, the average probability of success or accurate classification can be Fig. 7. Scatter plot of the k-uCSP features F 1and F 22,calculated by (10) from ‘left vs right-hand’ trials from subject 1 (circles left-hand, crosses right-hand). The features are sorted according a variance criterion: var F 1 ... var F p , where the variance is calculated across all the trials within a group. Superimposed we show the contour plot of the two-component best fitting Gaussian mixture model of the points. calculated as P 27 40 1 2 36 40 1 2 0.7875. Our last experiments are conducted on each group of 144 trials. We run a 10 fold cross-validation. To this end, the trials are divided into training (130 trials) and testing (14 trials). In the training phase, we calculate the unsupervised projection directions from the 130 trials in the training set. Then, we fit a two-component p-dimensional GMM to the resulting 130 log-variance feature vectors. Next, in the testing phase, the 14 trials of the testing set are projected onto the previously calculated directions, and the GMM is used for classifying the corresponding 14 p-vectors of log variance features. This is repeated 10 times, with different training and testing sets, and the results are averaged out. We observe that the performance largely depends on the type of motor imagery. For example, for subject 1, best results are obtained for the classification of ‘right-hand’ vs ‘tongue’ movements: we get 96.48 % of accuracy in session 1 and 97.19 % in session 2, 96.83 % in average. However, ‘feet’ and ‘tongue’ movements are hardly distinguishable one from another: only the 56.19 % (session 1), 61.81 % (session 2) and 59.0 % (average) of the corresponding trials are correctly classified. The following averaged classification accuracies are obtained for the combination of imagined movements with the best performance: Subject 1: 96.83 % (‘right-hand’ vs ‘tongue’), Subject 2: 65.57 % (‘feet’ vs ‘tongue’), Subject 3: 80.45 % (‘right-hand’ vs ‘tongue’), Subject 4: 72.52 % (‘left-hand’ vs ‘feet’), Subject 5: 62.88 % (‘left-hand’ vs ‘tongue’), Subject 6: 64.33 % (‘right-hand’ vs ‘feet’), Subject 7: 79.57 % (‘left-hand’ vs ‘tongue’), Subject 8: 93.07 % (‘left-hand’ vs ‘tongue’), Subject 9: 91.59 % (‘left-hand’ vs ‘tongue’). The average of all these numbers equals 78.31 %, not far from the mean classification accuracy (82.3 %) obtained
MARTÍN-CLEMENTE et al. : UNSUPERVISED CSP 2141 Fig. 8. Box plots of the unsupervised classification accuracy for each user. The boxes extend from the 25th to the 75th percentiles. Lines at either end of the boxes cover the maximum and minimum values. The red line in the middle is the median (50th percentile). Notches represent a confidence interval around the median. The black crosses are the mean values (see also the first column of Table II). TABLE I CLASSIFICATION ACCURACY (IN PERCENTAGE)ASSOCIATED WITH THE COMBINATION OF IMAGINED MOVEMENTS WITH THE BEST PERFORMANCE FOR THE UNSUPERVISED (K-UCSP WITH GMM CLASSIFIER)AND SUPERVISED (CSP AND FISHER LDA) APPROACHES by the unsupervised workload classification approach in [8]. It is a good performance considering its unsupervised nature. For comparison, Table I also gives the results obtained by the supervised CSP approach, where supervised classification is carried out with Fisher linear discriminant analysis (LDA) [30]. The complete results are shown in Figure 8, grouped by subject. For example, the accuracy of subject 1 ranges from 59.19 % to 97.19 %, depending on the motor imagery under consideration. Table II presents the mean value for each subject. For comparison, Table II also shows the accuracy obtained by using the supervised CSP approach and Fisher linear discriminant. Complementarily, Fig. 9 and Table III show the performance of the criterion for each imagined movement, averaged across the subjects. A natural question to ask is why there is a marked difference in some data sets between the performance of the traditional supervised approach and the proposed unsupervised one. To address this question, we have investigated the separation between the classes in the feature space using t-Distributed Stochastic Neighbor Embedding (t-SNE) [31], [32]. This is a popular technique for the visualisation of high-dimensional data. It maps each data point to a location in a low (2 or 3) TABLE II MEAN CLASSIFICATION ACCURACY (IN PERCENTAGE), AVERAGED FOR EACH SUBJECT,COMPARING THE UNSUPERVISED (K-UCSP AND GMM CLASSIFIER)AND SUPERVISED (CSP AND FISHER LDA) APPROACHES Fig. 9. Box plots of the unsupervised classification accuracy for each imaginary movement (‘L-R’: left-hand vs right-hand, ‘L-T’: left-hand vs tongue, ‘R-T’: right-hand vs tongue, ‘L-F’: left-hand vs feet, ‘R-F’: right-hand vs feet, ‘T-F’: tongue vs feet). TABLE III MEAN CLASSIFICATION ACCURACY (IN PERCENTAGE)FOR IMAGINARY MOVEMENTS (‘L-R’: LEFT-HAND VS RIGHT-HAND,‘L-T’:LEFT-HAND VS TONGUE,‘R-T’:RIGHT-HAND VS TONGUE,‘L-F’:LEFT-HAND VS FEET,‘R-F’:RIGHT-HAND VS FEET,‘T-F’:TONGUE VS FEET), AVERAGED ACROSS SUBJECTS,COMPARING UNSUPERVISED AND SUPERVISED APPROACHES dimensional space. Such a mapping preserves the local structure of the data in the sense that if the data points form well-separated clusters in the original high-dimensional space, they will too in the dimension reduced space. Specifically, t-SNE is used to map the vectors F1,...,Fp that contain the log-variance features to a space of two dimensions. Let us see an illustrative example: Figure 10a visualizes the mapping of the ‘tongue’ and ‘left-hand’ k-uCSP log-variance features calculated from user 7 during Session 1. Each point represents one of the feature vectors. For comparison, 10b depicts the corresponding CSP-based log variance features. We see that the classes are separated, but not always well separated. Consequently, the unsupervised classification
2142 IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 27, NO. 10, OCTOBER 2019 Fig. 10. Visualisation of the ‘tongue vs left-hand’ log-variance features for user 7. Each class has a different colour set of spots. algorithm may not be able to differentiate the two groups. In practice, this will lead users to determine on the fly, probably by trial and error, which imagined movements are optimal for themselves (in view of Table III, it may be the combination ‘hand’ vs ‘tongue’ for most people). Before closing the section, it should be remarked that the choice of the GMM classifier based on the features of eqn. (10) is made for illustration purposes only, but is most probably suboptimal. More optimal choices of features and unsupervised classifiers, which are beyond the scope of the present work, are expected to improve unsupervised BCI results. VI. CONCLUSIONS Kurtosis maximization performs CSP in an unsupervised fashion, sparing the need for training data and correct labelling. Further work would also aim at more elaborate BCI systems based on the proposed blind CSP approach. APPENDIX A PROOF OF THEOREM 1 Maximizing (5b) is equivalent to maximizing f a π1σ4 1 π2σ4 2subject to the constraint h a π1σ2 1 π2σ2 2 1 0. Because of that constraint, the objective function f a is defined on a closed interval. This, together with the fact that f a is continuous, ensures that it has both a maximum and a minimum value that it can attain. Let L a,λ f a λh a be the Lagrangian function and let L a,λ be the Hessian matrix of L a,λ with respect to a,havingas i,j -entry L ij 2L ai aj a,λ . Additionally, the tangent plane of h a at a is defined as the set T a v:v∇h a 0 , where ∇represents the gradient with respect to a. We recall the following result [33, Chap. 20]: Theorem 2: Let a be a local maximizer of f a subject to h a 0. Then, there exits λ Rsuch that C1) ∇f a λ ∇h a 0, and C2) for all v T a ,wehavevL a ,λ v 0. If a is a local minimizer, condition C2 becomes vL a ,λ v 0. Next we check conditions C1 and C2 for our particular problem. Condition C1: We begin by computing the gradients ∇f a 4 π1σ2 11a π2σ2 22a ∇h a 2 π11a π22a . Setting the gradient of the Lagrange function to zero, we obtain the equation 2 π1σ2 11a π2σ2 22a λ π11a π22a 0.(11) Premultiplying by awe easily find that λ 2 π1σ4 1 π2σ4 2 π1σ2 1 π2σ2 2 . Substituting into (11), we get σ2 1 σ2 2 σ2 1 2a σ2 2 σ2 2 σ2 1 1a. This equation has two solutions: S1) σ2 21a σ2 12a (12) S2) σ2 1 σ2 2.(13) Observe that S1 corresponds to the CSP solution (3). Condition C2: To ascertain whether these solutions are maximizers or minimizers of the criterion, we need to consider the second-order condition C2. The Hessian matrix of the Lagrangian can also be decomposed as L a,λ F a λH a ,whereFand Hare the Hessian matrices of f a and h a , respectively. We first focus on solution (12), which, after some tedious algebraic manipulations, leads to F a 4 2 i 1 2πiia a i πiσ2 ii H a 2 π11 π22 . In addition, the tangent plane T a is the set of vectors v Rpsuch that v∇h a 0, implying that v π11 π22 a 0 v1a v2a 0 (14) where we have exploited in the last implication the fact that 1a σ2 1 σ2 2 2a , as follows from (12). Let us denote s2 i vivthe variance of class iafter projection onto v,i 1,2. It follows that vL a ,λ v ησ2 2s2 2 σ2 1 σ2 2 R v R a (15)
MARTÍN-CLEMENTE et al. : UNSUPERVISED CSP 2143 where ηdef 4π1π2 π1s2 1 π2s2 2 0 and we have recognized that s2 1 s2 2is, by definition, equal to R v . To check C2, we distinguish the following cases: Case 1: a a1, the dominant eigenvector of the GEVD problem defining CSP solution (3). This vector is also the maximizer of the CSP criterion R defined in eqn. (2), as recalled in Sec. II, and therefore R v R a1 , v T a1 . Now, if R a1 1, then σ2 1 σ2 2, and expression (15) is negative. In other words, a1is a maximizer of the kurtosis (5b). Otherwise, if R a1 1, then σ2 1 σ2 2and expression (15) is positive, so that a1defines a minimum of the kurtosis. Case 2: a ap, the least significant eigenvector of GEVD problem (3). As seen in Sec. II, this vector is also the minimizer of R ,sothatR v R ap , v T ap .Using the same reasoning as in the above case, apis a minimizer of the kurtosis criterion if R ap 1, whereas it defines a maximum if R ap 1. Case 3: a ai,1 i p, any of the remaining eigenvectors of GEVD problem (3). Here we exploit the fact that the generalized eigenvectors enjoy an orthogonality property with respect to covariance matrices 1and 2,i.e., a ikaj 0, for all 1 i j p,k 1,2. From (14), it follows that both a1and apbelong to T ai . Hence, the sign of expression (15) when v a1is different from that when v ap, thus defining a saddle point of the kurtosis criterion (5b). In the light of the above cases, we realize that only the maximizer a1and/or the minimizer apof CSP criterion (2) maximize the kurtosis. More precisely, the maximizer of kurtosis is a1if R a1 1andapif R ap 1. Because, by definition, R a1 R ap (Sec. II), we can never have R a1 1andR ap 1 simultaneously, the only possibility for neither point to maximize kurtosis. To conclude the proof of Theorem 1, it remains to study solution S2 given in eqn. (13). This solution corresponds to the case where the mixture of projected classes simplifies into a single Gaussian distribution, with κY 3, and turns out to be a global minimizer of kurtosis, as follows from next lemma: Lemma 1: κY a 3, a Rp. Proof: The kurtosis (5b) can be expressed as κ a 3 b 2 c 2 b c 2, with b cdef bDc,b σ2 1,σ2 2 ,c 1,1 , D diag π1,π 2 . By the Cauchy-Schwarz inequality, b c b c ,andthenκY a 3, with equality if and only b c, R 0 , implying σ2 1 σ2 2. APPENDIX B COMPUTING SEVERAL PROJECTION DIRECTIONS The whole set of eigenvectors of (12) can be computed as follows. It can be always assumed that E XX I, where Iis the identity matrix, under the assumption that the data has been whitened or sphered in a pre-processing step. Since in general the covariance matrix of the data is given by E XX π11 π22, whitening imposes that π22 I π11. Substituting in (3), we get a usual eigenvalue problem 1a δa,whereδ R a π1R a π2 . Finally, as the eigenvectors of a Hermitian matrix are orthogonal, they can be determined by maximizing the kurtosis under the constraint that the direction obtained in the nth step is orthogonal to the previously computed directions. A possible implementation is presented as pseudo-code in Algorithm 1, where μis the step size and the gradient ∇κY a can be easily obtained from eqn. (5a), i.e., κY a E Y4 E Y2 2,where Y aX,as ∇κY a κY a a 4 E Y2 2 E XY3 E XY E Y4 E Y2 .(16) Observe that this expression can be always evaluated without any knowledge of the data labels by just replacing the expectations with their sample average estimates. For example, note that [24], [25] implement this pseudo-code in Matlab. Furthermore, [24], [25] seek the optimal value of μusing an algebraic line search approach. Algorithm 1 Compute Orthogonal Directions Maximizing the Kurtosis 1: Center and whiten the data, X 1 2 X ¯ X , where ¯ Xis the mean and the covariance matrix of X 2: Choose randomly the initial projection vectors a1,...,ap. 3: repeat 4: for k 1to pdo 5: ak ak μ∇κY ak {Gradient-ascend update rule, see eqn. (16)} 6: ak ak k 1 n 1 a kan an{Gram-Schmidt orthogonalization} 7: ak ak ak 2 {Normalization} 8: end for 9: until convergence ACKNOWLEDGMENT The authors would like to thank the anonymous reviewers for their thoughtful comments and efforts towards improving our manuscript. REFERENCES [1] H. Yuan and B. He, “Brain–computer interfaces using sensorimotor rhythms: Current state and future perspectives,” IEEE Trans. Biomed. Eng., vol. 61, no. 5, pp. 1425–1435, May 2014. [2] F. Lotte, “A tutorial on EEG signal-processing techniques for mentalstate recognition in brain–computer interfaces,” in Guide to Brain- Computer Music Interfacing, E. R. Miranda, Ed. London, U.K.: Springer, 2014. [3] F. Lotte et al., “A review of classification algorithms for EEG-based brain–computer interfaces: A 10 year update,” J. Neural Eng., vol. 15, no. 3, 2018, Art. no. 031005. [4] H. Wang, “Harmonic mean of Kullback–Leibler divergences for optimizing multi-class EEG spatio-temporal filters,” Neural Process. Lett., vol. 36, no. 2, pp. 161–171, Oct. 2012. [5] A. S. Aghaei, M. S. Mahanta, and K. N. Plataniotis, “Separable common spatio-spectral patterns for motor imagery BCI systems,” IEEE Trans. Biomed. Eng., vol. 63, no. 1, pp. 15–29, Jan. 2016. [6] T. Jiang et al., “Characterization and decoding the spatial patterns of hand extension/flexion using high-density ECoG,” IEEE Trans. Neural Syst. Rehabil. Eng., vol. 25, no. 4, pp. 370–379, Apr. 2017.