Information Theoretic Approaches for Motor-Imagery BCI Systems: Review and Experimental Comparison

entropy Review Information Theoretic Approaches for Motor-Imagery BCI Systems: Review and Experimental Comparison Rubén Martín-Clemente ID 1,*, Javier Olias 1ID , Deepa Beeta Thiyam 1,2, Andrzej Cichocki 3,4,5 ID and Sergio Cruces 1,*ID 1Departamento de Teoría de la Señal y Comunicaciones, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Seville, Spain; [email protected] (J.O.); [email protected] (D.B.T.) 2Department of Sensor and Biomedical Technology, School of Electronics Engineering, VIT University, Vellore, Tamil Nadu 632014, India 3 Skolkovo Institute of Science and Technology (Skoltech), Moscow 143026, Russia; [email protected] or [email protected] 4Laboratory for Advanced Brain Signal Processing, Brain Science Institute, RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan 5Systems Research Institute, Polish Academy of Sciences, 01-447 Warsaw, Poland *Correspondence: [email protected] (R.M.-C.); ser[email protected] (S.C.); Tel.: +34-954-487-475 Received: 31 October 2017; Accepted: 19 December 2017; Published: 2 January 2018 Abstract: Brain computer interfaces (BCIs) have been attracting a great interest in recent years. The common spatial patterns (CSP) technique is a well-established approach to the spatial filtering of the electroencephalogram (EEG) data in BCI applications. Even though CSP was originally proposed from a heuristic viewpoint, it can be also built on very strong foundations using information theory. This paper reviews the relationship between CSP and several information-theoretic approaches, including the Kullback–Leibler divergence, the Beta divergence and the Alpha-Beta log-det (AB-LD)divergence. We also revise other approaches based on the idea of selecting those features that are maximally informative about the class labels. The performance of all the methods will be also compared via experiments. Keywords: common spatial patterns; generalized divergences; brain computer interfaces 1. Introduction The electroencephalogram (EEG) is a record over time of the differences of potential that exist between different locations on the surface of the head [ 1 , 2 ]. It originates from the summation of the synchronous electrical activity of millions of neurons distributed within the cortex. In recent years, there has been a growing interest in using the EEG as a new communication channel between humans and computers. Brain-computer interfaces (BCIs) are computer-based systems that enable us to control a device with the mind, without any muscular intervention [ 3 – 6 ]. This technology, though not yet mature, has a number of therapeutic applications, such as the control of wheelchairs by persons with severe disabilities, but also finds use in fields as diverse as gaming, art or access control. There are several possible approaches for designing a BCI [ 1 , 4 , 7 ]. Among them, motor imagery (MI)-based BCI systems seem to be the most promising option [ 6 , 8 – 10 ]. In MI-based BCI systems, the subject is asked to imagine the movement of different parts of his or her body, such as the hands or the feet. The imagined actions are then translated into different device commands (e.g., when the subject imagines the motion of the left hand, the wheelchair is instructed to turn to the left). What makes this possible is that the spatial distribution of the EEG differs between different imagined movements. More precisely, since each brain hemisphere mainly controls the opposite side of the body, Entropy 2018,20, 7; doi:10.3390/e20010007 www.mdpi.com/journal/entropy Entropy 2018,20, 7 2 of 29 the imagination of right and left limb movements produces a change of power over the contralateral left and right brain motor areas. These fluctuations, which are due to a pair of phenomena known as event-related desynchronization (ERD) or power decrease and event-related synchronization (ERS) or power increase [ 11 , 12 ], can be detected and converted into numerical features. By repeating the imagined actions several times, a classifier can be trained to determine which kind of motion the subject is imagining (see [ 13 ] for a review). In practice, three classes of MI are used in BCIs, namely the movements of the hands, the feet and the tongue. Left hand movement imagery is more prominent in the vicinity of the electrode C4 (see Figure 1), while right hand imagined actions are detected around electrode C3 [ 14 ]. The imagery of feet movements appears in the electrode Cz and its surrounding area; nevertheless, it is not usually possible to distinguish between left foot or right foot motor imagery because the corresponding activation areas are too close in the cortex [ 11 , 14 ]. Finally, imagery of tongue movements can be detected on the primary motor cortex and the premotor cortex [ 15 ]. One of the inherent difficulties of designing a BCI is that the EEG features are highly non-stationary and vary over sessions. To cope with this problem, the background state of the subject (i.e., his or her motivation, fatigue, etcetera) and the context of the experiment can be both modeled as latent variables, whose parameters can be estimated using the expectation-maximization (EM) algorithm [ 16 , 17 ]. Overall, current BCI approaches achieve success rates of over 90%, although much depends on the person from whom the EEG data are recorded [14]. Cz T4C4 C3 T3 Pz Fz T6 O2 T5 F7 F8 O1 Fp1 Fp2 F4F3 P3 P4 Figure 1. Electrode locations of the international 10–20 system for EEG recording. The letters “F”, “T”, “C”, “P” and “O” stand for frontal, temporal, central, parietal and occipital lobes, respectively. Even numbers correspond to electrodes placed on the right hemisphere, whereas odd numbers refer to those on the left hemisphere. The “z” refers to electrodes placed in the midline. The common spatial patterns (CSP) method [ 4 , 18 – 22 ] is a method of dimensionality reduction that is widely used in BCI systems as a preprocessing step. Basically, assuming two classes of MI-EEG signals (e.g., left hand and right hand MI tasks), CSP projects the EEG signals onto a low-dimensional subspace, which captures the variability of one of the classes while, at the same time, trying to minimize the variance in the other class. The goal is to enhance the ability of the BCI to discriminate between the different MI tasks, and it has been shown that CSP is able to reduce the dimension of the data significantly without decreasing the classification rate. It is noteworthy that CSP admits an interesting probabilistic interpretation. Under the assumption of Gaussian distributed data, CSP is equivalent to maximizing the symmetric Kullback–Leibler (KL) divergence between the probability distributions of the two classes after the projection onto the low dimensional space [ 23 , 24 ]. As a generalization of this idea in the context of BCI, it is interesting to investigate the dimensionality reduction ability of other different divergence-based criteria, which is drawing a lot of interest among the computational neuroscience community. Entropy 2018,20, 7 3 of 29 The present manuscript is a review of the state of the art of information theoretic approaches for motor imagery BCI systems. The article is written as a guideline for researchers and developers both in the fields of information theory and BCI, and the goal is to simplify and organize the ideas. We will present a number of approaches based on Kullback–Leibler divergence, Beta divergence (which is a generalization of Kullback–Leibler’s) and Alpha-Beta log-det (AB-LD) divergence (which include as special cases Stein’s loss, the S-divergence or the Riemannian metric), as well as their relation to CSP. We will also review a technique based on the idea of selecting those features that are maximally informative about the class labels. Complementarily, for the purpose of comparison, several non-information theoretic variants of CSP and their different regularization schemes are revised in the paper. The performance of all approaches will be evaluated and compared through simulations using both real and synthetic datasets. The paper is organized as follows: The CSP algorithm is introduced in Section 2. Section 3 introduces the main characteristics of the Kullback–Leibler divergence, the Beta divergence and the Alpha-Beta log-det divergence, respectively, as well as their application to the problem of designing MI-BCI systems and the algorithms used to optimize them. Section 4reviews an information-theoretic feature extraction framework. Section 5presents, as has been said before, several extensions of CSP not based on information-theoretic principles. Finally, Section 6presents the results of some experiments in which the performances of the above criteria are tested, in terms of their accuracy, computational burden and robustness against errors. EEG Measurement and Preprocessing For measuring the EEG, several different standardized electrode placement configurations exist. The most common among them is the International 10–20 system, which uses a set of electrodes placed at locations defined relative to certain anatomical landmarks (see Figure 1). The ground reference electrode is usually positioned at the ears or at the mastoid. To obtain a reference-free system, it is common practice to calculate the average of all the electrode potentials and subtract it from the measurements [1,2]. The EEG is usually contaminated by several types of noise and artifacts. Eye blinks, for example, elicit a large potential difference between the cornea and the retina that can be several orders of magnitude greater than the EEG. In the rest of the paper, it is assumed that the signals have already been pre-processed to remove noise and interferences. To this end, several techniques [ 25 ], such as autoregressive modeling [ 26 ], the more complex independent component analysis (ICA) [ 27 ], or the signal space projection (SSP) method [ 28 ], have shown good or excellent results (see also [ 29 ] and the references therein). Signal preprocessing includes also the division of the EEG into several frequency bands that are separately analyzed [ 30 , 31 ]. The “mu” band (8–15 Hz) and the “beta” band (16–31 Hz) are particularly useful in BCIs, as they originate from the sensorimotor cortex, i.e., the area that controls voluntary movements [2]. 2. The Common Spatial Pattern Criterion In this section, we present the common spatial patterns (CSP) method [ 4 , 18 – 22 , 32 , 33 ]. Consider a two-class classification problem, where the EEG signals belong to exactly one of two classes or conditions (e.g., left-/right-hand movement imagination). To fix notation, let Xi,k∈RD×T be the matrix that contains the EEG data of class i∈ { 1, 2 } in the k -th trial or experiment, where D is the number of channels and T the number of samples in a trial. The corresponding sample covariance estimator is defined by: Σi,k=1 T−1Xi,kX> i,k, (1) Entropy 2018,20, 7 4 of 29 where (·)> denotes “transpose”. Here, the EEG signals are assumed to have zero-mean, which is fulfilled as they are band-pass filtered (see the previous section). If L trials per class are performed, the spatial covariance matrix for class i is usually calculated by averaging the trial covariance matrices as: Σi=1 L L ∑ k=1 Σi,k(2) In practice, these covariance matrices are often normalized in power with the help of the following transformation: Σi←Σi/ tr (Σi), (3) where tr(·)denotes the trace operator. After the BCI training phase, in which matrices Σ1 and Σ2 are estimated using training data, suppose that a new, not previously observed, data matrix X∈RD×T of imagined action is captured. The problem that arises is to develop a rule to allocate these new data to one class or the other. A useful approach is to define a weight vector w∈RD (also known as a ‘spatial filter’) and allocate X to one class if the variance of w>X exceeds a certain predefined threshold and to the other if not; this relates to the fact that event-related desynchronizations and event-related synchronizations, i.e., the phenomena underlying the MI responses, are associated with power decreases/increases of the ongoing EEG activity [12]. Of course, not just any spatial filter is of value. To enhance the discrimination of the MI tasks, CSP proposes using spatial filters that maximize the variance of the band-pass filtered EEG signals in one class while, simultaneously, minimizing it for the other class. Mathematically, CSP aims at maximizing an objective function based on the the following Rayleigh quotient: J(w) = w>Σ1w w>Σ2w=σ2 1 σ2 2 , (4) where σ2 iis the variance of the i-th projected class and Σiis the covariance matrix of the i-th class. It is a straightforward derivation to obtain that the spatial filters that hierarchically maximize (4) can be computed by solving the generalized eigenvalue problem: Σ1w=λΣ2w. (5) Each eigenvector wigives a different solution. Observe that: w> iΣ1wi=λiw> iΣ2wi→λi=w> iΣ1wi w> iΣ2wi =J(wi), where λi is the generalized eigenvalue corresponding to wi . Therefore, the larger (or smaller) the eigenvalue, the larger the ratio between the variances of the two classes and the better the discrimination accuracy of the filter. The latter readily suggests selecting the spatial filters among the principal and the minor eigenvectors (i.e., the eigenvectors associated with the largest and smallest eigenvalues, respectively). Let: WCSP = [w1, . . . , wd]∈RD×d(6) be the matrix that collects these d≤D top (i.e., most discriminating) spatial filters. Given a data matrix X∈RD×Tof observations, all of the same class, the outputs of the spatial filters are defined as: yi=w> iX,i=1, . . . , d,(d≤D)(7) Entropy 2018,20, 7 5 of 29 which can be gathered at the d×T output matrix Y=W> CSPX . Denoting by Σ the sample covariance matrix of X , it follows that the covariance matrix of the outputs is given by W> CSPΣWCSP , while the variance of the output of the i -th spatial filter is equal to w> iΣwi . Finally, not the sample variances, but the log transformed sample variances of the outputs, i.e., Fi=log(w> iΣwi),i=1, . . . , d, (8) are used as features for the classification of the imagined movements. Observe that, as long as d<D , the dimensionality of the data is reduced. CSP admits an interesting neurological interpretation. First note that the scalp EEG electrodes measure the addition of numerous sources of neural activity, which are spread over large areas of the neocortical surface, and this does not always allow a reliable localization of the cortical generators of the electrical potentials. It has been suggested that CSP linearly combines the EEG signals so that the sources of interest are enhanced while the others are suppressed [34]. Another interpretation of (4) may be as follows: the basic theory of principal component analysis (PCA) states that maximizing w>Σiw finds the direction vector that best fits, in the least-squares sense, the data of class i in the D -dimensional space. Similarly, minimizing this ratio obtains the opposite effect. Thus, we can interpret that CSP seeks directions that fit well with the data in one class, but are not representative of the data in the other class. By projecting the EEG data onto them, a significant reduction of the variance of one of the classes, while preserving the information content of the other, can be thus obtained. An interesting generative model perspective has been proposed in [ 35 , 36 ]. Here, the above data matrices are assumed to be generated by a latent variable model: Xi(:, k) = A Yi(k) + Ni(k), where we have used the notation Xi(: , k)∈RD for the k -th column of the data matrix Xi , i.e., it is the observation vector at time k for class i , i= 1, 2; A∈RD×s is a mixing matrix, the same for both classes; Yi(k)∼ N(0 , Γi) is an s -dimensional column vector of latent variables ( s has to be estimated from the data) and Ni(k)∼ N(0 , ∆i) is a D -dimensional vector of noise, independent of the data. Here, the covariance matrices Γi and ∆i are assumed to be diagonal matrices, implying that the latent factors are also independent of each other. Under this model, the columns of matrix A can be regarded as the “spatial patterns” that explain how the EEG data are formed at each electrode location, where the latent variables represent the degree to which each “spatial pattern” appears in the data. Under the assumptions that the noise is negligible and matrix A is square, it is noteworthy that the CSP spatial filters are precisely the columns of the matrix A−> [35]. 3. Divergence-Based Criteria CSP produces quite good results in general, but also suffers from various shortcomings: e.g., it is sensitive to artifacts [ 37 , 38 ] and its performance is degraded for non-stationary data [ 39 ]. For these reasons, CSP is still an active line of research, and a number of variants have been proposed in the literature. In particular, in this paper, we are interested in reviewing CSP-variants based on an information-theoretic framework. There is a common assumption in the literature that the classes can be modeled by multivariate Gaussian distributions with zero-means and different covariance matrices. This assumption is based on the principle of maximum entropy, not in actual measures of EEG data. By projecting the data onto the principal generalized eigenvectors, CSP transforms them onto a lower dimensional space where the variance of Class 1 is maximized, while the variance for Class 2 is minimized. Conversely, the projection onto the minor generalized eigenvectors has the opposite effect. Since a zero-mean univariate normal variable is completely determined by its variance, we can understand the ratio (4) as a measure of how much the distributions of the projected classes differ from each other (the larger Entropy 2018,20, 7 6 of 29 the ratio between the variances, the more different the distributions). By accepting this viewpoint, it is interesting to investigate the ability of other measures of dissimilarity between statistical distributions, rather than the ratio of the corresponding variances, to help in discriminating between the classes. In fact, the most interesting features for classification often belong to those subspaces where there is a large dissimilarity between the conditional densities of the considered classes, which is another justification for proposing a divergence maximization framework in the context of MI-BCI. In the following sections, we review the main information-theoretic-based approaches. 3.1. Criterion Based on the Symmetric Kullback–Leibler Divergence Divergences are functions that measure the dissimilarity or separation between two statistical distributions. Given two univariate Gaussian densities N1( 0, σ1) and N2( 0, σ2) , their Kullback–Leibler divergence (the KL divergence between two distributions f1 and f2 is defined as DivKL (f1kf2)=R∞ −∞f1(x)log f1(x) f2(x)dx) is easily found to be: DivKL (N1(0, σ1)kN2(0, σ2))=1 2 log σ2 2 σ2 1 +σ2 1 σ2 2−1!. If the densities have interchangeable roles, it is reasonable to consider the use of a symmetrized measure like the one provided by the symmetrized Kullback–Leibler (sKL) divergence. This is defined simply as: sDivKL (N1kN2)=DivKL (N1kN2)+DivKL (N2kN1) =1 2 σ2 1 σ2 2 +σ2 2 σ2 1!−1. (9) The resemblance to the CSP criterion (4) is quite obvious, as was already noted, e.g., in [ 24 ]. In particular, note that, since z+1 z increases when z goes to either infinity or zero, (9) is maximized by either maximizing or minimizing the ratio of the variances σ1and σ2. The generalization to multivariate data is straightforward. Let Y=W>X , where X is the observed data matrix and W= [w1 , . . . , wd]∈RD×d denotes an arbitrary matrix of spatial filters with 1 ≤d≤D . Under the assumption that the EEG data are conditionally Gaussian distributed for each class ck∈ { 1, 2 } , i.e., X|ck∼ N(0 , Σk) , the spatially-filtered data are also from a normal distribution, i.e., Y|ck∼ N(0,¯ Σk), where: ¯ Σk=W>ΣkW∈Rd×d, k= 1, 2. The KL divergence between two d -dimensional multivariate Gaussian densities f1=N1(0 , ¯ Σ1) and f2=N2(0,¯ Σ2), that is, DivKL (f1kf2)=ZN1(0,¯ Σ1)log N1(0,¯ Σ1) N2(0,¯ Σ2)dy, can be shown to be (after some algebra): DivKL (f1kf2)=1 2log |¯ Σ2| |¯ Σ1|−d+trace(¯ Σ2)−1(¯ Σ1), (10) where |·| stands for “determinant”. The symmetrized Kullback–Leibler (sKL) divergence between the probability distributions of the two classes is now defined as: Entropy 2018,20, 7 7 of 29 DivsKL (f1kf2)=DivKL (f1kf2)+DivKL (f2kf1) =1 2trace (W>Σ1W)−1(W>Σ2W) +(W>Σ2W)−1(W>Σ1W)−d, where we show again the explicit dependency on W. We can naturally extend this formula to define the equivalent sKL matrix divergence: DsKL(W>Σ1WkW>Σ2W) = 1 2trace (W>Σ1W)−1(W>Σ2W) +(W>Σ2W)−1(W>Σ1W)−d. (11) It has been shown in [ 23 ] that the subspace of the filters that maximize the sKL matrix divergence, WsKL =arg maxWDsKL(W>Σ1WkW>Σ2W), (12) coincides with the subspace of those that maximize the CSP criterion, in the sense that the columns of WsKL and WCSP span the same subspace: span(WsKL) = span(WCSP), (13) that is, every column of WsKL is a combination of the top spatial filters of WCSP and vice versa. In practice, WsKL is first used to project the data onto a lower dimensional subspace, and then, WCSP is determined by applying CSP to the projected data. Some advantage can be gained, compared to using CSP only, if in the first step the optimization of the sKL matrix divergence is also combined with some suitable regularization scheme. For example, to fight against issues caused by the non-stationarity of the EEG data, it has been proposed to maximize the regularized objective function [23]: LsKL(W) = (1−φ)DsKL(W>Σ1WkW>Σ2W)−φ∆(W), (14) where 0 ≤φ<1 and: ∆(W) = 1 2L 2 ∑ i=1 L ∑ k=1 DivKL N(0, W>Σi,kW)kN(0, W>ΣiW)(15) is a regularization term, where we have assumed that L trials per class have been performed and Σc,k is the covariance matrix in the k -th trial of class c∈ { 1, 2 } . This proposed regularization term enforces the transformed data in all the trials to have the same statistical distribution. Other ideas have been proposed in [ 23 ], and a related approach can be found in [ 40 ]. Observe also that (15) is defined on the basis of the KL divergence, not on its symmetrized version. The KL divergence is calculated by a formula similar to (10), giving: DivKL N(0, W>Σi,kW)kN(0, W>ΣiW)=1 2log |W>ΣiW| |W>Σi,kW|−d+trace(W>ΣiW)−1(W>Σi,kW).(16) The inverse of Σi,k does not appear in (16) , which makes sense if this matrix is ill-conditioned due to insufficient sample size. For this reason, the KL divergence is preferred to its symmetric counterpart. In addition, the logarithm in (16) downweights the effect of |W>Σi,kW|−1 in case Σi,k is nearly singular. Entropy 2018,20, 7 8 of 29 3.2. Criterion Based on the Beta Divergence The beta divergence, which is a generalization of the Kullback–Leibler’s, seems to be an obvious alternative measure of discrepancy between Gaussians. Given two zero-mean multivariate probability density functions f1(y)and f2(y), the beta divergence is defined for β>0 as: Divβ(f1(y)kf2(y)) = 1 βZfβ 1(y)−fβ 2(y)f1(y)dy−1 β+1Zfβ+1 1(y)−fβ+1 2(y)dy. As limβ→0 fβ 1−fβ 2 β=log f1 f2 , it can be shown that the beta divergence converges to the KL divergence for β→0. Let f1=N( 0, ¯ Σ1) and f2=N( 0, ¯ Σ2) , with ¯ Σi=W>ΣiW∈Rd×d , i= 1, 2, be the zero-mean Gaussian distributions of the spatially-filtered data. In this case, the symmetric beta divergence between them yields the following closed form formula [41]: Dsβ(W>Σ1WkW>Σ2W) = γ |¯ Σ1|−β/2 +|¯ Σ2|−β/2 −(β+1)d/2 |¯ Σ2|1−β 2 |β¯ Σ1+¯ Σ2|1/2 +|¯ Σ1|1−β 2 |β¯ Σ2+¯ Σ1|1/2 !!,(17) where γ=1 βq1 (2π)βd(β+1)d . Observe that Dsβ is somewhat protected against possible large increases in the elements of Σ1 or Σ2 caused by outliers or estimation errors. For example, if Σi (resp. ¯ Σi ) grows, i∈ { 1, 2 } , then the contribution of all the terms containing Σi (resp. ¯ Σi ) in (17) tends to vanish. Compared with the previous case, if Σ1(for example) increases, then the term: trace (W>Σ2W)−1(W>Σ1W) may dominate (11). With the necessary changes of divergences being made, the regularizing framework previously defined by Equations (14) and (15) can be easily adapted to the present case [ 23 ]. It has been argued in [ 23 ] that small values of β penalize abrupt changes in the covariance matrices caused by single extreme events, such as artifacts, whereas a large β is more suitable to penalize the gradual changes over the dataset from trial to trial. Alternatively, supposing that L trials per class are performed, it has been also proposed in [ 41 ] to use as the objective function the sum of trial-wise divergences: ¯ Dsβ(W) = L ∑ i=1 Dsβ(W>Σ1,iWkW>Σ2,iW), where Σ1,iand Σ2,iare the covariance matrices in the i-th trial of Class 1 and Class 2, respectively. 3.3. Criterion Based on the Alpha-Beta Log-Det Divergence Given the covariance matrices of each class, Σ1 and Σ2 , an extension of the Kullback–Leibler symmetric matrix divergence given in Equation (11) is the Alpha-Beta log-det (AB-LD) divergence, defined as [42,43]: D(α,β) LD (Σ1kΣ2) = 1 αβ log  α(Σ−1 2 2Σ1Σ−1 2 2)β+β(Σ−1 2 2Σ1Σ−1 2 2)−α α+β+ (18) for α6=0, β6=0, α+β6=0, Entropy 2018,20, 7 9 of 29 where: |x|+=(x x ≥0, 0, x<0, denotes the non-negative truncation operator. For the singular cases, the definition becomes: D(α,β) LD (Σ1kΣ2) =                                1 α2tr (Σ 1 2 2Σ−1 1Σ 1 2 2)α−I−αlog |Σ 1 2 2Σ−1 1Σ 1 2 2|for α6=0, β=0, 1 β2tr (Σ−1 2 2Σ1Σ−1 2 2)β−I−βlog |Σ−1 2 2Σ1Σ−1 2 2|for α=0, β6=0, 1 α2log  (Σ−1 2 2Σ1Σ−1 2 2)α(I+log(Σ−1 2 2Σ1Σ−1 2 2)−α)+ for α=−β, 1 2||log(Σ 1 2 2Σ−1 1Σ 1 2 2)||2 Ffor α,β=0. (19) It can be easily checked that D(α,β) LD (Σ1kΣ2) = 0 iff Σ1=Σ2 . The interest in the AB-LD divergence is motivated by the fact that, as can be observed in Figure 2, it generalizes several existing log-det matrix divergences, such as the Stein’s loss (the Kullback–Leibler matrix divergence), the S-divergence, the Alpha and Beta log-det families of divergences and the geodesic distance between covariance matrices (the squared Riemannian metric), among others [43]. 1 1 α β DKL (Σ1|| Σ2) DR(Σ1|| Σ2) KL matrix div. Squared Riemannian metric Dual KL matrix div. DKL (Σ2|| Σ1) Power log-det div. DLD ( α , α )(Σ1|| Σ2) Beta log-det div. S-divergence 4× Alpha log-det div. DLD (1, β )(Σ1|| Σ2) DLD ( α ,1− α )(Σ1|| Σ2) Figure 2. Illustration of the Alpha-Beta log-det divergence (AB-LD) divergence D(α,β) LD (Σ1kΣ2) in the (α , β) -plane. Note that the position of each divergence is specified by the value of the hyperparameters (α , β) . This parameterization smoothly connects several positive definite matrix divergences, such as the squared Riemannian metric ( α= 0, β= 0), the KL matrix divergence or Stein’s loss ( α= 1, β= 0), the dual KL matrix divergence (α=0, β=1) and the S-divergence (α=1 2,β=1 2), among others. Entropy 2018,20, 7 16 of 29 The second regularization approach consists of including a penalty term in the CSP objective function [64]. The regularized CSP objective functions can be represented as: ˜ J1(w) = w>Σ1w w>Σ2w+α∆(w)(47) ˜ J2(w) = w>Σ2w w>Σ1w+α∆(w)(48) where α is the regularization parameter. The regularized Tikhonov-CSP approach (RTCSP) penalizes the solutions with large weights by using a penalty term ∆(w)of the form: ∆(w) = kwk. The filters w can computed by solving an eigenvalue problem similar to that of the standard CSP algorithm. Specifically, the stationary points of ˜ J1(w)verify [64]: (Σ2+αI)−1Σ1w=λw. Similarly, the stationary points of ˜ J2(w) are the eigenvectors of matrix (Σ1+αI)−1Σ2 . Observe that it is necessary to optimize both objective functions, as the stationary points of any of them alone maximize the variance of one class, but do not minimize the variance of the other class. Finally, all the previous approaches admit the following generalization: in traditional CSP, the EEG data is usually band-pass pre-filtered using one single filter between 8 and 30 Hz, which is a range that covers the so-called “alpha”, “beta” and “mu” EEG bands. An straightforward extension, known as the filter bank CSP (FBCSP) technique, was proposed in [ 30 ], where the input MI-EEG signals are bandpass filtered between different bands of frequency ((4–8 Hz), (8–12 Hz), . . . , (36–40 Hz)) and the CSP algorithm, or any of its variants, is applied to each band for the computation of the spatial filters. The results of all analyses are then combined to form the final response (see Figure 4). Similar approaches have been proposed in [ 10 , 65 , 66 ]. An extension to the multiclass problem can be found in [ 67 ]. Since the optimal frequency bands can vary from subject to subject, several alternative approaches have been proposed that combine the time-frequency characteristics of the EEG data [ 68 , 69 ] for improving the classification accuracy and reducing the number of electrodes [70]. Input EEG signal LDA Class Output 15-30 Hz 7-15 Hz 4-7 Hz 30-40 Hz Frequency filtering CSP Algorithm CSP Algorithm CSP Algorithm CSP Algorithm Spatial filtering Figure 4. Architecture of filter bank CSP. LDA is shorthand for Linear Discriminant Analysis. Entropy 2018,20, 7 17 of 29 6. Experimental Results Initially, we will test the algorithms using real datasets obtained from the BCI competition III (dataset 3a) and BCI competition IV (dataset 2a), which are publicly available at [ 71 ]. On the one hand, the dataset 3a from BCI competition III consists of EEG data acquired from three subjects (k3b, k6b and l1b) at a sampling frequency of 250 Hz using a 60-channel EEG system. In each trial, an arrow to the left, right, up or down was shown on a display for a few seconds, and in response to the stimulus, the subject was asked to respectively perform left hand, right hand, tongue and foot MI movements. The dataset consists of 90 trials per class for Subject k3b and 60 trials per class for Subjects k6b and l1b. On the other hand, the dataset 2a from BCI competition IV was acquired by using 22 channels from nine subjects (A01–A09) while also performing left hand, right hand, tongue and foot MI movements following a similar procedure. The signals were also sampled at 250 Hz and were recorded in two sessions on different days, each of them with 72 trials per each class. For a total of four possible motor-imagery (MI) movements, (4 2)= 6 different combinations of pairs of MI movements (i.e., left hand-right hand, left hand-foot, left hand-tongue, right hand-foot, right hand-tongue, foot-tongue) can be formed. The experiments below consider all possible combinations: since 12 users are available and for nine of them we have recordings performed on two different days, this makes a total of 3 × 6 + 9 × 6 × 2 = 126 different experiments. We repeated eight times each of these 126 possible experiments, and results were averaged. For each repetition, 60 trials were selected at random from each MI movement, which were split into 40 trials for training and 20 trials for testing. Additionally, in the case of the BCI competition IV, we averaged over the two sessions conducted for each user to avoid biasing the statistical tests. As a result, 3 × 6 + 9 × 6 × 2 / 2 = 72 averaged performance measures are finally available for each algorithm. The data have been initially bandpass filtered between the cut-off frequencies of 8–30 Hz, except before using the FBCSP method, which as we explained in Section 5, considers four bands for covering the frequency range between 4 and 40 Hz. The information of the classes in each trial is summarized by their respective covariance matrices. These matrices are estimated, normalized by their trace and used as input to the algorithms that carry out the calculation of the spatial filters prior to the MI classification, which is performed by using linear discriminant analysis (LDA). The only parameter of the CSP algorithm is the number of spatial filters that one would like to consider. Although, this number d is usually fixed a priori for each dataset, it is advantageous to estimate automatically the best number of spatial filters for each user by using the combination of cross-validation and hypothesis testing proposed in [ 72 ]. Figure 5a illustrates this fact. The figure represents the scatter plot of the accuracies, expressed as a percentage, that have been respectively obtained by the CSP algorithm for a fixed value of d= 8 (x-axis) and for the estimation of the best value of d (y-axis). These estimated accuracies have been obtained by averaging eight test samples, as explained above. The accuracies obtained for different individuals or for different pairs of conditions can be reasonably considered approximately independent and nearly Gaussian. Under this hypothesis, a one-sided paired t-test of statistical significance can be used to compare the results obtained by both alternatives. Let δf(m) = fy(m)−fx(m) be the paired differences of accuracy ((y-axis value) vs. (x-axis value)) for m=1, . . . , M, where M=72 is the number of samples. Then, the averaged difference is: ∆f=1 M M ∑ m=1 δf(m)(49) and the unbiased estimate of its variance is: s2=s2 ∆f M, (50) Entropy 2018,20, 7 18 of 29 where s2 ∆f=1 M−1∑M m=1(δf(m)−∆f)2 . Under the null hypothesis ( H 0) that the expected performance values coincide, i.e., E[fy(m)] = E[fx(m)], the t-statistic: T−STAT =∆f s∆f/√M. (51) follows a Student’s tdistribution with M− 1 degrees of freedom. Thus, the probability that the null hypothesis can generate a t-statistic larger than T−STAT gives the p-value of the right-sided test: P−VAL =Prob(t>T−STAT|H0). (52) The more positive is T−STAT , the smaller is the P−VAL , and the probability of observing a t-statistics larger than T−STAT decreases under the null hypothesis. When the p-value falls below the 0.05 threshold of significance, the hypothesis of not having a performance improvement when using the alternative procedure can be rejected, because this would correspond to a quite improbable situation. On the contrary, if the p-value of the right-sided test is above 0.05, the null hypothesis cannot be rejected. In this particular case, the p-value of the test in Figure 5a is below 0.05; therefore, one can reject the hypothesis that the automatic estimation of d does not improve the results over the method that a priori selects d=8 filters. We briefly name and describe below some of the implementations that optimize the already mentioned criteria for dimensionality reduction in MI-BCIs. Because of the substantially higher computational complexity of most of the alternatives to CSP (see Table 1), it is not practical to develop a specific automatic estimation procedure of the number of spatial filters for each of them. For this reason, we will consider in their implementations the same number of spatial filters that was automatically estimated for CSP. • CSP (see Section 2) and ITFE (see Section 4): apart from the number of spatial filters, these two methods do not have hyper-parameters to tune. Their respective algorithms have been implemented according to the specifications given in [18,49]. • RTCSP (see Section 5): RTCSP has a regularization parameter, which has been selected by five-fold cross-validation in { 0, 0.1, 0.2, . . . , 1 } . The MATLAB implementation of this algorithm has been obtained from [73]. • FBCSP (see Section 5): In this case, we have used a variation of the algorithm in [ 30 ]. The selected frequency bands correspond to the brainwaves theta (4–7 Hz), alpha (8–15 Hz), beta (16–31 Hz) and low gamma (32–40 Hz), where five-fold cross-validation has been used to select the best combination of these frequency bands. We extract d features from each band, where d is selected using the method in [72]. • DivCSP (see Sections 3.2 and 3.4). The values of β and φ (the regularization parameter) have been selected by five-fold cross-validation, β∈[ 0, 1 ] , φ∈[ 0, 0.5 ] . This divergence includes the KL divergence as a particular case when β= 0. MATLAB code of the algorithm has been downloaded from [ 74 ] and used without any modification. Optimization has been performed using the so-called subspace method (see Section 3.4). • Sub-LD (sub-space log-det): this algorithm, which also belongs to the class of the subspace methods, is based on the criterion in [ 42 ] to maximize the Alpha-Beta log-det divergence (see Sections 3.3 and 3.4). In this paper, the implementation of the algorithm is based on the BFGS method on the Stiefel manifold of semi-orthogonal matrices and takes as the initialization point the solution obtained by the CSP algorithm. The regularization parameter η has been chosen by five-fold cross-validation in the range of values (− 0.2, 0.2 ) , which are not far from zero. The negative values of η favor the expansion of the clusters, while the positive values favor their contraction. For η close to zero, the solution of this criterion should not be far from that of CSP, Entropy 2018,20, 7 19 of 29 which improves the convergence time of the algorithm and reduces the impact of the values of α , β in the results, so both parameters have been fixed to 0.5. 55 60 65 70 75 80 85 90 95 100 CSP + FIXED d=8 55 60 65 70 75 80 85 90 95 100 CSP + ESTIMATED d P-VAL: 0.01 T-STAT: 2.22 (a) 2 4 6 8 10 12 NUMBER OF SPATIAL FILTERS 0 0.05 0.1 0.15 0.2 0.25 0.3 SAMPLE PROBABILITY (b) Figure 5. Illustration of the advantages in performance of using an automatic cross-validation method to estimate the best even number of features d with respect to using an a priori fixed value of d . The automatic method relies on the technique proposed in [ 72 ], which was implemented here using one-sided t-tests of significance instead of the original two-sided tests. ( a ) Scatter plot comparison of the accuracies (in percentage) obtained by the CSP algorithm for fixed d= 8 (x-axis) and for the automatic estimation of d(y-axis); (b) histogram of the estimated best even number of features d. Table 1shows the typical execution time of a single run of each algorithm, programmed in MATLAB language, in a PC with Intel I7-6700 CPU @ 3.4-GHz processor and 16 GB of RAM. The algorithms that use cross-validation for selecting the hyper-parameters need more iterations, hence the run time has to be multiplied by the number of the hyper-parameters combinations that are evaluated. Entropy 2018,20, 7 20 of 29 Table 1. Computational burden of the considered algorithms, which are sorted in increasing value of their respective execution times without using cross-validation. FBCSP, filter bank CSP; ITFE, information theoretic feature extraction. Algorithm Time (s) CSP 0.0017 FBCSP 0.0050 ITFE 0.3070 Sub-LD 1.0538 DivCSP 4.6696 Figure 6represent the boxplot of the accuracy of the algorithms, considering together all the combinations of the motor imagery movements from all subjects in datasets III 3a and IV 2a. The p-values and t-statistics shown below the box-plots of Figure 6are above the 5% threshold of significance, revealing that, in this experiment, one cannot reject the null hypotheses. It follows that the expected accuracies of the alternative algorithms are not significantly higher than the expected accuracies obtained with CSP. Supporting this conclusion, Figure 7represents the specific boxplots that corresponds to MI movements involving the right hand. Additionally, we have tested, in the case “left hand versus right hand”, whether the improvement obtained by using the alternative algorithms is significant or not. The accuracy in the classification and the corresponding p-values of the tests are shown in Figure 8. The results reveal that, in general and except in a few isolated cases, the null hypothesis that the other methods do not significantly improve performance over CSP cannot be discarded. CSP Sub-LD RTCSP ITFE FBCSP DivCSP ALGORITHM 30 40 50 60 70 80 90 100 ACCURACY (%) REFERENCE P-VAL: 0.3 T-STAT: 0.51 P-VAL: 0.98 T-STAT: -2.05 P-VAL: 1 T-STAT: -3.6 P-VAL: 1 T-STAT: -4.9 P-VAL: 1 T-STAT: -5.9 Figure 6. Comparison of the expected accuracy percentages obtained by each of the considered algorithms. The figure shows box-plot illustrations where the median is shown in red line, while the 25% and 75% percentiles are respectively at the bottom and top of each box. Larger positive values T−STAT  0 and smaller P−VAL  1 / 2 would correspond with greater expected improvements over CSP. However, none of the p-values, which are shown below their respective box-plots, is able to attain the 5% threshold level of significance ( P−VAL < 0.05), so the possible improvements cannot be claimed to be statistically significant with respect to those obtained by CSP. Entropy 2018,20, 7 21 of 29 CSP Sub-LD RTCSP FBCSP ITFE DivCSP ALGORITHM 50 60 70 80 90 100 ACCURACY (%) (a) CSP Sub-LD RTCSP FBCSP ITFE DivCSP ALGORITHM 50 60 70 80 90 100 ACCURACY (%) (b) CSP Sub-LD RTCSP FBCSP ITFE DivCSP ALGORITHM 50 60 70 80 90 100 ACCURACY (%) (c) Figure 7. Performance of the algorithms for different motor imagery combinations involving the right hand. ( a ) Right-hand versus left-hand motor imagery classification; ( b ) right-hand versus feet motor imagery classification; (c) right-hand versus tongue motor imagery classification. The results of Figure 6were obtained by choosing through cross-validation the best possible values for the different parameters of the algorithms. Figures 9and 10 show how many times each value of the parameters has been selected after cross-validation. They also show the number of times that CSP outperformed the corresponding algorithm, the number of times that the algorithm outperformed CSP or the cases in which both of them were equivalent. Without limiting the foregoing, it must be also remarked that the alternative algorithms perform better than CSP for some subjects and MI movements. Entropy 2018,20, 7 22 of 29 k3b k6b l1b A01 A02 A03 A04 A05 A06 A07 A08 A09 USER 50 55 60 65 70 75 80 85 90 95 100 ACCURACY CSP Sub-LD RTCSP ITFE FBCSP DivCSP (a) k3b k6b l1b A01 A02 A03 A04 A05 A06 A07 A08 Al09 USER 10-2 10-1 100 P-VALUE Sub-LD RTCSP ITFE FBCSP DivCSP Threshold (b) Figure 8. Accuracy percentages and p-values for the testing of an improvement in performance over CSP when the right hand versus left hand movement imagination are discriminated. The results reveal that, in general and except in a few isolated cases, the null hypothesis that the other methods do not significantly improve the performance over CSP cannot be discarded. ( a ) Average accuracy obtained by the algorithms for each subject; ( b )p-values of the t-tests that compare whether the performance of the alternative algorithms is significantly better than the one obtained by CSP. The horizontal dashed line represents the threshold level of significance of 5%. Entropy 2018,20, 7 23 of 29 Figure 9. Histogram of the values of the regularization parameter in the Sub-LD algorithm that have been chosen by cross-validation. (a) (b) Figure 10. Histogram of the hyper-parameters of the DivCSP algorithm selected by cross-validation. (a) Case with β∈[0, 0.5]and φ=0; (b) case with β=0.5 and φ∈[0, 0.5]. 6.1. Results on Artificially Perturbed Data In order to study the performance of the algorithms under artificial perturbations of the datasets we have conducted two experiments. The first one consists of introducing random label changes in the real datasets, while the second one defines sample EEG covariance matrices for each condition and artificially introduces outlier covariance matrices in the training procedure to quantify the resulting deterioration in performance. Exchanging labels of the training set at random is one of the most harmful perturbations that one can consider in a real experiment. It models the failure of the subjects to imagine the correct target MI movements due to fatigue or lack of concentration. For this experiment, we selected a subject who has a relatively good performance in absence of perturbations. Figure 11 presents the progressive degradation of the accuracy of the algorithms as the percentage of mismatched labels increases. Entropy 2018,20, 7 24 of 29 0 5 10 15 20 25 30 PERCENTAGE OF MISMATCHED LABELS 65 70 75 80 85 90 95 100 ACCURACY (%) CSP Sub-LD RTCSP FBCSP ITFE Figure 11. Comparison of the accuracy percentages obtained by each of the considered algorithms with respect to the percentage of mismatched labels in the training set. This experiment illustrates deterioration of the performance of the algorithms with respect to the increase of the percentage of randomly switched labels of the motor imagery movements. In the second experiment, we have created artificial EEG data and consider the effect of adding random outliers. The artificial data were generated starting from two auxiliary covariance matrices Ck , k= 1, 2 for the construction of the conditional covariance matrices of each class. These covariances were generated randomly by drawing two random Gaussian matrices A(k) with i.i.d. elements a(i) ij ∼ N( 0, 1 ) and forming the covariance matrices with Ck=A(k)(A(k))> , k= 1, 2. In order to control the difficulty of the classification problem, we introduce a dissimilitude parameter δ∈[ 0, 1 ] that interpolates between the two auxiliary covariance matrices as follows: Σ1=C1/2 1(C−1/2 1C2C−1/2 1)(1−δ) 2C1/2 1(53) Σ2=C1/2 2(C−1/2 2C1C−1/2 2)(1−δ) 2C1/2 2(54) In this way, when δ= 0, the two interpolated covariance matrices coincide Σ1=Σ2 , and it is impossible to distinguish between them. On the contrary, when δ= 1, we obtain the original randomly generated matrices Σ1=C1 and Σ2=C2 . The matrices Σk are used as the expected covariance matrix of the observations for class k , while the sample covariance matrices for each trial are generated from a Wishart distribution with scale matrix 1 TΣk and T degrees of freedom (where T denotes the trial length). The outlier matrices have been generated following a similar scheme, though interpolation is not used and the resulting covariances are scaled by a factor of five. In our simulations with artificial data, we have set the dissimilitude parameter to δ= 0.1. The results obtained for artificial data and with different percentages of outlier covariance matrices in the training set are shown in Figure 12. One can observe how the performance progressively deteriorates with the number of outliers, similarly for all the methods, although at a smaller rate than in the case having the same percentage of mismatched labels. The parameters of the algorithms have been selected by cross-validation. Entropy 2018,20, 7 25 of 29 0 5 10 15 20 25 30 PERCENTAGE OF OUTLIERS 65 70 75 80 85 90 95 100 ACCURACY (%) CSP Sub-LD RTCSP ITFE Figure 12. Accuracy percentages versus the percentage of training trials with outliers in a synthetic classification experiment. 7. Conclusions In this paper, we have reviewed several information theoretic approaches for motor-imagery BCI systems. In particular, we have focused on those based on the Kullback–Leibler divergence, Beta divergence, Alpha-Beta log-det divergence and information theoretic feature extraction, exploring the existing links with common spatial patterns, which is a widely-used technique for spatial filtering in BCI applications. The performance of all these methods has been evaluated through experimental simulations using real and synthetic data. In general, the results obtained for real data from BCI competitions reveal a similar performance for all the considered criteria in terms of their percentages of accuracy. However, CSP clearly outperforms the other methods when comparing the required computational burdens. In the case of synthetic data with outliers, a comparison of the divergence-based methods with small regularization parameters reveals that they can slightly increase the frequency of obtaining a better performance, although the average accuracy results are still similar to those obtained with CSP. Therefore, although these divergence-based methods are not yet a practical alternative to CSP, this line of research is in its infancy, and divergence-based methods can have an underlying potential for improvements in performance that remains to be explored. Acknowledgments: Part of this work was supported by the Spanish Government under MICINN project TEC2014-53103-P. Andrzej Cichocki was partially supported by the MES Russian Federation grant 14.756.31.0001. We also thank the anonymous reviewers for their careful reading of our manuscript and their many insightful comments and suggestions. Author Contributions: Rubén Martín-Clemente and Sergio Cruces collaborated in writing the paper and coordinating the study. Andrzej Cichocki critically revised the manuscript by providing inspiring comments. Javier Olias and Deepa Beeta Thiyam conducted the experimental work. All authors have read and approved the final manuscript. Conflicts of Interest: The authors declare no conflict of interest.