EEG Signal Processing in MI-BCI Applications With Improved Covariance Matrix Estimators

1534-4320 (c) 2018 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2019.2905894, IEEE Transactions on Neural Systems and Rehabilitation Engineering IEEE TRANSACTIONS ON ..., VOL. X, NO. X, AUGUST 2018 1 EEG Signal Processing in MI-BCI Applications with Improved Covariance Matrix Estimators Javier Olias, Rub´ en Mart´ ın-Clemente, MaAuxiliadora Sarmiento-Vega and Sergio Cruces Abstract—In brain-computer interfaces the typical models of the EEG observations usually lead to a poor estimation of the trial covariance matrices, given the high non-stationarity of the EEG sources. We propose the application of two techniques that significantly improve the accuracy of these estimations and can be combined with a wide range of motor imagery BCI methods. The first one scales the observations in such a way that implicitly normalizes the common temporal strength of the source activities. When the scaling applies independently to the trials of the observations the procedure justifies and improves the classical preprocessing for the EEG data. Additionally, when the scaling is instantaneous and independent for each sample, the procedure particularizes to Tyler’s method in statistics for obtaining a distribution-free estimate of scattering. In this case, the proposal provides an original interpretation of this existing method as a technique that pursuits an implicit instantaneous power-normalization of the underlying source processes. The second technique applies to the classifier and improves its performance through a convenient regularization of the features covariance matrix. Experimental tests reveal that a combination of the proposed techniques with state-of-the-art algorithms for motor-imagery classification provides a significant improvement in the classification results. Index Terms—Common spatial pattern, brain-computer interfaces, motor-imagery classification, covariance matrix estimation. I. INTRODUCTION Brain-computer interfaces (BCI) have a great potential for enabling the communication between machine and humans by means of the analysis of the electroencephalographic activity. Nowadays, almost all the Motor Imagery BCI (MI-BCI) systems summarize most of the relevant information about the measurements in two kinds of covariance matrices: the covariance matrices of the filtered observations (employed for dimensionality reduction) and the covariance matrices of the features (which are required for classification). In the dimensionality reduction stage one tries to select those subspaces of the observations that retain most of the discriminative power, for instance, using the technique of Common Spatial Patterns (CSP) [1]. After that, the features are usually chosen as a non-linear transformation of the band-power statistics of the projected observations onto the previously selected subspaces [2]. The covariance matrices of these features (together with their class-conditional expectations) play a relevant role in the classification stage of MI-BCI [3]. Although CSP was All the authors are with the Department of Teor´ ıa de la Se˜ nal y Comunicaciones, Universidad de Sevilla, Camino de los Descubrimientos s/n, Seville 41092, Spain. E-mails: {folias, ruben, sarmiento, sergio}@us.es This work has been led by the corresponding author: Dr. Cruces. This work was supported in part by the Spanish Government under MINECO grant TEC2017-82807-P. only suitable for two-class classification problems, some later alternatives have been also proposed for multi-class settings (see, for instance, [4], [5]). There are several sources of difficulty in the processing of EEG signals. Among them, we may cite: the inevitable presence of noise and interference at the sensors, the low spatial resolution of the BCI headsets [1], the possible presence of outliers in the measurements [6], the difficulty in gathering sufficient data trials for training [7], the need to determine the suitable number of features in those method that apply to dimensionality reduction [8], and the non-stationarity of the EEG signals [9]–[11]. The non-stationary can happen at different levels. The classical inter-subject and inter-session variabilities have been frequently addressed in the literature [11]. In this work, we will shift our attention to the less studied variabilities that happen between trials, and also within samples of the same trial. The signals generated by the brain are non-stationary in power at the trial and sample levels. We will show later that this power variability hinders the correct estimation of the covariance matrices of the trials, which are the most used statistics in the existing MI-BCI implementations. Our experimental results avail the hypothesis that the correction of this EEG signal variability leads to improved covariance matrix estimates, which allow transversal improvements in accuracy for the tested classification algorithms. The main contributions of the article are the following: •We show that the standard power normalization of the observations, which is widely used in the preprocessing of the EEG data for MI-BCI, is useful but suboptimal. •We propose the power-normalization of the effective EEG source activities. This normalization has no hyperparameters and, in general, improves the quality of the covariance matrix estimates during training and testing. •The shrinkage of the feature covariance matrices in MIBCI was shown to be beneficial when the number of training trials is small [12]. We propose the application of an alternative shrinkage estimate (gLDA) that is based on the Gaussianity of the features [13]. Our experimental results confirm that the proposed powernormalization and gLDA implementation lead to a transversal improvement in the performance of the existing MI-BCI algorithms. In addition, the proposal seems to be much less sensitive with respect to the number of features employed in the dimensionality reduction stage. The article is organized as follows. Section II introduces the basic model of the EEG measurements and section III discusses some classical and state-of-the-art approaches for 1534-4320 (c) 2018 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2019.2905894, IEEE Transactions on Neural Systems and Rehabilitation Engineering 2 IEEE TRANSACTIONS ON ..., VOL. X, NO. X, AUGUST 2018 MI-BCI. Section IV presents an overcomplete model of the observations and defines the effective sources of the mixture. Section V describes the proposal for the normalization in power of the EEG sources and also analyzes its links with the standard preprocessing of the observations. This method is extended in section VI to the case of the instantaneous powernormalization of the effective sources. Section VII presents some variations in the implementations of the classifier using shrinkage estimates of the feature covariance matrices. The experimental results are provided and discussed in section VIII, while section IX is devoted to the conclusions. II. BASIC MODEL OF THE EEG OBSERVATIONS The EEG headset is based on an array of sensors that measures the electromagnetic activity on the scalp. At time t, the variations of the activities of the sensors are measured with respect to a given referential system (see EEG referencing in [1]) and passband filtered to retain the 8Hz-32Hz band. After that, they are centered at the origin by subtracting the estimated mean of each trial and collected into the observation vector x(t) = [x1(t), . . . , xNx(t)]T∈RNx. The physiological nature of the problem allows one to model the ith-element of the observation vector as a superposition of contributions from: some desired latent EEG source activities sj(t), j = 1, . . . , Ns, and some filtered additive interference or noise component which we denote by nj(t), j = 1, . . . , Nx. We will not assume any specific value for Nswhich, depending on the experiment, could be greater or lower than Nx. The contribution of ith-source sj(t) to the jth-observation xj(t)is modeled as aijsj(t), where the factor aij refers to the attenuation of the almost instantaneous propagation of the source activity to the sensor position. In vector form, the filtered observations are known to follow the linear instantaneous mixing model [14] x(t) = As(t) + n(t),(1) where A= [aij]ij ∈RNx×Nsrefers to the mixing matrix. In those cases where we would like to make explicit the trial to which the observations belong to, we will use the notation xτ(t)that refers to the vector of observations of the trial τ at time t. The global power of the non-stationary process of filtered observations is defined by Px≡ hE[kxτ(t)k2]it,τ =1 NτT T X t=1 Nτ X τ=1 E[kxτ(t)k2].(2) A column-wise concatenation of the observed vector samples from a trial results in the matrix model of the observations Xτ=ASτ+Nτ,(3) where Xτ,Nτ∈RNx×Tand Sτ∈RNs×T. In the following, the class of a trial τwill be denoted by c(τ)∈ {c1, . . . , cK}. III. THE COMMON SPATIAL PATTERNS AND OTHER SUCCESSFUL APPROACHES FOR MI-BCI The Common Spatial Patterns (CSP) is a method designed for the case of having two classes (K= 2) [15]. Let the classconditional covariance matrices of the classes be Σx|c1and Σx|c2. The CSP algorithm (see Table I ) tries to reduce the dimensionality of the observations by finding a p-dimensional subspace for which the two classes are maximally separated in a certain divergence sense [6], [16]. This goal is achieved by setting the p < Nxspatial filters w1, . . . , wp(for the sake of simplicity pis assumed to be even) equal to the p/2principal and p/2minor eigenvectors of the following generalized eigenvalue problem Σx|c1w=Σx|c2wλ . (4) The selected eigenvectors are grouped in the matrix of spatial filters W= [w1, . . . , wp], which is used to perform the dimensionality reduction of the observations Yτ=WTXτ∈Rp×T.(5) There are several possible extension of CSP to multi-class (K > 2) scenarios. Some are based on the joint approximated diagonalization of the covariances matrices of the observations for each class [14] WTΣx|ckW=Dckk= 1, . . . , K, (6) where Dckrefers to an approximately diagonal matrix. The one proposed in [4] combines this approximated diagonalization with a method to choose the most relevant filters based on an Information Theoretic Feature Extraction criterion (ITFE). Although the dimensionality reduction stage (implemented by CSP and ITFE) sometimes is omitted, in general, as we will see in the simulations, it is a recommended procedure for datasets with moderate or relatively large number of sensors. After the dimensionality reduction, some basic linear classification results can be obtained using Fisher’s Linear Discriminant Analysis (LDA). However, other state-of-the-art proposals are nowadays preferable. This is the case of sLDA [7] a shrinkage variant of LDA and also of the classifiers that exploit the Riemmanian geometry of the manifold of symmetric and positive definite (SPD) matrices. Among these classifiers, we can mention the Riemannian Minimum Distance to Mean (RMDM) [5], which is based on the minimization of the Riemmanian distance between the sample covariance matrices of the test trials and the Riemmanian mean of the classes. Other improved classification methods are obtained by using as features the projection of the sample covariance matrices onto the tangent space (of the Riemmanian SPD manifold) at the referential Riemmanian mean of the set of covariance matrices. In this way, an LDA classifier applied to tangent space (TS) features give rise to a TSLDA implementation. Similarly, the logistic regression (LR) classification of TS features leads to a TSLR implementation [17]. The interested reader in Riemmanian approaches for Brain-Computer Interfaces can find in [17] and [18] respective tutorial reviews on this topic. A. Classical estimation of the class covariance matrices As the observations have been already centered, the EEG spatial covariance matrix of trial τis given by C(0) Xτ=1 TXτXT τ.(7) 1534-4320 (c) 2018 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2019.2905894, IEEE Transactions on Neural Systems and Rehabilitation Engineering OL´ IAS et al.: EEG PROCESSING WITH IMPROVED COVARIANCE MATRIX ESTIMATORS 3 The notation C(i) Xτis adopted in this paper in order to allow the possibility to refine this estimate through additional iterations. Then, since the trials may have unequal power, the standard CSP implementation [1] normalizes the EEG covariance matrices as C(1) Xτ= C(0) Xτ Tr{C(0) Xτ}/Nx ≡XτXT τ 1 NxTr{XτXT τ}.(8) One may note that this definition only differs from the classical normalization in the following irrelevant 1 Nxscaling term, which is mainly adopted here for notational convenience. Finally, the class-conditional covariance matrices are usually estimated by means of the arithmetic mean of the trials ˆ Σ(1) x|ck=1 NckX τ:c(τ)=ck C(1) Xτk= 1, . . . , K. (9) In the following sections, we propose an alternative normalization for the training and test covariance matrices. It has no additional hyperparameters and, in general, outperforms the standard one considered in (8). In particular, we will show that the standard normalization can be regarded as a first approximation to the proposed approach. At this point, it is worth to comment other estimators for Σx|ckwhich have been suggested for MI-BCI applications according to various strategies. The adaptation with respect to differences between the training and testing distributions of the data has been considered in [19], which suggests the weighting of the samples according to their estimated importance. In [20], class covariance matrices estimators of minimum βdivergence for a Wishart model have been proposed to ensure the robustness with respect to data outliers. The solution, which is based on an iteratively weighting of the trial covariance matrices of each class, uses cross-validation (CV) for the determination of the hyper-parameter of the divergence. The use of CV is also required in [21], which proposed several regularized covariance matrix estimates with the aim to avoid overfitting. One regularized estimate, which has the remarkable advantage of avoiding CV, was proposed by Ledoit and Wolf in [22]. IV. OVERCOMPLETE MODEL OF THE OBSERVATIONS AND EFFECTIVE COMPONENT OF THE SOURCES The linear mixing model of equation (1) provides an overcomplete representation of the observations. This can be seen by integrating the noise/interference contribution into an extended sources vector s0 τ(t)to obtain xτ(t)=(A I)sτ(t) nτ(t)=A0s0 τ(t).(10) Moreover, there is an inherent linear indeterminacy between the sources and the columns of the mixing matrix. In this sense, note that, for any arbitrary invertible matrix M∈ R(Ns+Nx)×(Ns+Nx), the model satisfies xτ(t) = A0s0 τ(t)=(A0M−1) (Ms0 τ(t)).(11) We avoid this indeterminacy by assuming, from here on, that the global covariance matrix of the source signal process is equal to the identity matrix. As we initially considered the centering of the observations, this matrix is then given by Σ0 s=hE[s0 τ(t)(s0 τ(t))T]it,τ =I,and the global covariance matrix of the observations is Σx=hE[xτ(t)(xτ(t))T]it,τ =A0A0T.(12) The fact that the resulting mixing matrix A0∈ RNx×(Ns+Nx)is wide and of rank Nx, implies that not all the components of the extended vector of sources s0(t)will contribute to the observations. Only the component of the sources that is aligned with the range space of the rows of A0will have an effective contribution, while the orthogonal component to this subspace will be discarded. To see this, consider the orthogonal decomposition of the extended sources s0 τ(t) = ΠA0Ts0 τ(t) + Π⊥ A0Ts0 τ(t),(13) where the proyection matrix onto the rows of the extended mixing matrix is ΠA0T=A0T(A0A0T)−1A0and the orthogonal projection matrix is given by Π⊥ A0T=I−ΠA0T. Since these projection matrices satisfy A0ΠA0T=A0and A0Π⊥ A0T=0, it is easily observed that xτ(t) = A0s0 τ(t) = A0˜ sτ(t)(14) where ˜ sτ(t) = ΠA0Ts0 τ(t)represents the effective sources, i.e., the component of the extended sources with a non-negligible influence in the value of the observations xτ(t). Moreover, it is straightforward to check that the global covariance matrix of the effective sources coincides with the following projection matrix Σ˜ s= ΠA0T, which is unitary similar to the identity matrix of dimension Nx. Hence, the global power of the effective sources is P˜ s= Tr{Σ˜ s}= Tr{INx}=Nx.(15) V. NORMALIZATION OF THE POWER OF THE SOURCES Let’s define the power of the effective sources for trial τas P˜ Sτ= Tr{C˜ Sτ}=1 TTr{˜ Sτ˜ ST τ}.(16) When ˜ Sτfor τ= 1, . . . , Nτhave dissimilar powers, their contribution to the class-conditional covariance matrices is not homogeneous. In this situation, a fraction of the trials may dominate the estimation, implying a higher variance in the estimates. The covariance normalization by the power of the observations P(0) Xτ= Tr{C(0) Xτ}/Nxin (8) only partially alleviates the previous effect, since, due to the equivalence P(0) Xτ= Tr{AC(0) ˜ Sτ AT}, it depends on the interaction between the mixing matrix and the trial covariance of the sources. Instead, we propose to equalize the power of the effective sources in each trial in such a way that they all coincide with the global power of the process, which was defined in (15). Although we don’t have direct access to ˜ Sτ, we explain in the sequel a method that allows us to iteratively equalize its power, contributing in this way to obtain more reliable estimates of the covariance matrices. 1534-4320 (c) 2018 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2019.2905894, IEEE Transactions on Neural Systems and Rehabilitation Engineering 4 IEEE TRANSACTIONS ON ..., VOL. X, NO. X, AUGUST 2018 TABLE I PSEUDO-CODE OF CSP+LDA ALGORITHM FOR MI-BCI. PREPROCESSING FOR TRAINING & TESTING Freq. filtering & centering of the data ∀τ. Compute C(0) Xτand C(1) Xτ,∀τ, using (7)-(8). Determine ˆ Σ(1) x|c1and ˆ Σ(1) x|c2with (9). METHOD FOR DIM. REDUCTION (STANDARD-CSP) % Obtain the spatial filters solving (4) [V,D] = eig(ˆ Σ(1) x|c1,ˆ Σ(1) x|c2) % Sort the Nxsolutions [∼, ind] = sort(diag(D)),V=V(:, ind), % Select the extreme peigenvectors W= [V(:,1 : p/2),V(:, Nx−p/2) + 1 : Nx] % Spatial filtering CYτ=WTC(0) Xτ W, τ = 1,...,Nτ. % Transf. for obtaining normal-like features fτ= log ( diag(CYτ)/sum(diag(CYτ)) ) (F1) LEARNING THE LDA (BINARY) CLASSIFIER % Using the training pairs (c(τ),fτ) µk=1 Nck Pτ:c(τ)=ck fτfor k= 1,2. ˆ Σf|ck=1 Nck Pτ:c(τ)=ck(fτ−µk)(fτ−µk)T, k = 1,2. p(ck) = Nck/(Nc1+Nc2), k = 1,2. ˆ Σf=p(c1)ˆ Σf|c1+p(c2)ˆ Σf|c2 α=ˆ Σ−1 f(µ1−µ2) β=1 2(µ1+µ2)−log p(c1)−log p(c2) (µ1−µ2)Tˆ Σ−1 f(µ1−µ2)(µ1−µ2) CLASSIFICATION OF TEST DATA % Implementation using LDA with equal cov. CYτ=WTC(0) Xτ W, τ ∈Set of test trials. Evaluate fτusing the same formula as in (F1) ˆc(τ) = ckwhere k = 1.5 + 0.5 sign(αT(fτ−β)). Consider the notation for the inner product between two symmetric positive definite matrices of dimension Nx, hCXτ,Σ−1 xi=1 Nx Tr{CXτΣ−1 x}.(17) Lemma 1 (Power of the effective sources): The power of the effective sources for each trial τis given by the scaled inner product between the covariance matrix of the trial and the inverse of the global covariance matrix of the observations P˜ Sτ=NxhCXτ,Σ−1 xi.(18) This lemma, which provides an exact formula for the evaluation of the power of the effective sources, is proved in Appendix A. However, the determination of Σx= hE[xτ(t)(xτ(t))T]it,τ involves an expectation operation and, as a consequence, is not feasible. Instead, we can estimate it from the available samples at a given iteration i−1. Under a Gaussian mixture model for the observations, a natural estimate of Σx(built from a combination of maximum likelihood estimates) is given by the arithmetic mean of the covariances of the trials ˆ Σ(i−1) x=hC(i−1) Xτiτ=1 Nτ Nτ X τ=1 C(i−1) Xτ.(19) After substituting ˆ Σ(i−1) xfor Σxin (18), the estimated power of the effective sources at iteration i−1is ˆ P(i−1) ˜ Sτ . The ratio between the power of the effective sources at iteration (i−1) and their global power (obtained in (15)) is given by (ˆσ(i−1) τ)2≡1 Nx ˆ P(i−1) ˜ Sτ=hC(0) Xτ,(ˆ Σ(i−1) x)−1i.(20) In order to equalize the power across trials at the ith iteration, we should normalize the observations as X(i) τ=Xτ/ˆσ(i−1) τ,(21) since this scaling replaces the estimated power ˆ P(i−1) ˜ Sτ of the effective sources in each trial by the global average power P˜ s=Nx. After that, the scaled observations X(i) τlead to normalized estimates of the trial covariance matrices C(i) Xτ=1 TX(i) τ(X(i) τ)T= (ˆσ(i−1) τ)−2C(0) Xτ∀τ(22) and to an improved estimate of the global covariance matrix ˆ Σ(i) x=hC(i) Xτiτ=1 Nτ Nτ X τ=1 C(i) Xτ.(23) This new estimate can still help in improving the normalization of the sources, so the estimation procedure can continue in a recursive manner until the relative variation in the estimate of the global covariance matrix falls below a tolerance threshold . For instance, by continuing with the iteration until the following condition is met: kˆ Σ(i) x−ˆ Σ(i−1) xkF/kˆ Σ(i) xkF< . After the convergence of the iteration, the following average covariance matrices of each class are used as inputs to the method of dimensionality reduction ˆ Σ(i) x|ck=1 NckX τ:c(τ)=ck C(i) Xτk= 1, . . . , K. (24) A. Expliciting the link with the standard preprocessing of the EEG observations At iteration i= 0, before having access to the observed data, we may consider an initial isotropic estimate for the covariance matrix of the observations ˆ Σ(0) x=I. Hence, the estimates of the covariance matrices in (22) are, for i= 1, equal to C(1) Xτ= (ˆσ(0) τ)−2C(0) Xτ=XτXT τ 1 NxTr{XτXT τ}∀τ, (25) which exactly coincide with those provided by the standard normalization of the trials in (8). Next, the global covariance matrix ˆ Σ(1) xis evaluated using (23) and used, in another iteration (i= 2), to improve the normalization of the observations in each trial. Then, the new trial covariance matrices are C(2) Xτ= (ˆσ(1) τ)−2C(0) Xτ= C(0) Xτ hC(0) Xτ,(ˆ Σ(1) x)−1i ∀τ(26) 1534-4320 (c) 2018 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2019.2905894, IEEE Transactions on Neural Systems and Rehabilitation Engineering OL´ IAS et al.: EEG PROCESSING WITH IMPROVED COVARIANCE MATRIX ESTIMATORS 5 and again the new global covariance matrix ˆ Σ(2) xis evaluated. One can continue with the iterations of the procedure until it convergences. In the section of simulations, we will later illustrate with a controlled experiment (see Figure 1) the improvement of the estimates of the trial covariance matrices with respect to the number of iterations. Although we have previously suggested the initialization of the iteration with ˆ Σ(0) x=Ifor revealing the link between the proposal and the classical preprocessing of CSP, in practice, it is better to choose as initial estimate the sample covariance matrix of the trials ˆ Σ(0) x=<C(0) Xτ>τ. This latter estimate is more informative than the identity matrix, which contributes to a faster convergence of the iteration. VI. INSTANTANEOUS POWER NORMALIZATION LEADS TO AN EXISTING ESTIMATOR OF SCATTER Until now, in order to illustrate the links of the proposed power-normalization with the preprocessing used in classical CSP, we have only addressed the equalization of the power across trials. However, the technique is easily extended for equalizing the power of the sources over temporal juxtaposed (or overlapped) windows of arbitrary length. For signals like the EEG sources, which are highly non-stationary, one can improve the estimates of the covariance matrices by equalizing the power across samples, i.e., considering windows of one sample length. Let us consider the instantaneous correlation matrix estimate of the observations at the trial τand time t C(0) xτ(t)=xτ(t)(xτ(t))T,(27) which is based on a single sample. In similarity with (20), given ˆ Σxat iteration (i−1), we obtain the power ratio for each trial and time sample (ˆσ(i−1) τ,t )2≡ hC(0) xτ(t),(ˆ Σ(i−1) x)−1i(28) =1 Nx Tr{xτ(t)(xτ(t))T(ˆ Σ(i−1) x)−1}(29) =1 Nx Tr{(xτ(t))T(ˆ Σ(i−1) x)−1xτ(t)}.(30) Its evaluation with (30), is recommended in the instantaneous case because of the computational advantages over (29). The instantaneous power-normalization of ˜ sτ(t)is simply obtained by scaling the observations x(i) τ(t) = xτ(t)/ˆσ(i−1) τ,t ∀τ, t. (31) Therefore, the covariance matrices estimates of each trial C(i) Xτ=1 T T X t=1 x(i) τ(t)(x(i) τ(t))T(32) are, in general, more reliable and contribute, using (23), to an improved estimation of the averaged covariance matrix ˆ Σ(i) x. The whole iteration over the set of training trials is summarized in the top part of Table II. The procedure for the estimation of the covariance matrix of the test trials, which is shown in the second part of Table II, is coherent with the updates performed in the last iteration for the training trials. TABLE II PSEUDOCODE OF THE INSTANTANEOUS POWER-NORMALIZATION,WHICH PARTICULARIZES TO A VERSION OF TYLER’S METHOD IN STATISTICS FOR OBTAINING A ROBUST ESTIMATOR OF SCATTER. PREPROCESSING FOR TRAINING TRIALS Freq. filtering & centering of the data ∀τ. Set i=0, ˆ Σ(0) x=<C(0) Xτ>τ Repeat i=i+1 (ˆσ(i−1) τ,t )2=1 NxTr{(xτ(t))T(ˆ Σ(i−1) x)−1xτ(t)} ∀τ, t C(i) Xτ=1 TPT t=1(ˆσ(i−1) τ,t )−2xτ(t)(xτ(t))T∀τ ˆ Σ(i) x=1 NτPNτ τ=1 C(i) Xτ Until kˆ Σ(i) x−ˆ Σ(i−1) xkF/kˆ Σ(i) xkF<  Return C(i) Xτ∀τand ˆ Σ(i) x|ck=1 Nck Pτ:c(τ)=ck C(i) Xτk= 1, . . . , K. PREPROCESSING FOR A TESTING TRIAL τ Freq. filtering & centering of the trial. Given the last used estimate ˆ Σ(i−1) xin training... Evaluate (ˆσ(i−1) τ,t )2=1 NxTr{(xτ(t))T(ˆ Σ(i−1) x)−1xτ(t)} ∀t Return C(i) Xτ=1 TPT t=1(ˆσ(i−1) τ,t )−2xτ(t)(xτ(t))T∀τ ∗Note: after this preprocessing the evaluation of the features no longer needs normalization, i.e., fτ= log ( diag(CYτ) ) (F2) The combination of the instantaneous power-normalization with CSP will be referred, hereinafter, as nCSP. However, since the proposed normalization aims to recover the stationarity in power of the effective sources vector, it is unnecessary to apply any additional normalization of the features. Hence, the recommended evaluation of the features is simply given by formula (F2) of Table II, i.e., fτ= log ( diag(CYτ) ), which replaces all the instances of formula (F1) in Table I. In Appendix B, we discuss the link between the instantaneous power-normalization iteration and the method proposed by Tyler in [23] for obtaining a distribution-free estimator of scatter within the class of elliptically distributed data. Both techniques use complementary arguments that arrive at a similar final result. However, Tyler’s method assumes that the observations are drawn from an elliptical distribution, a hypothesis that may be true for a single trial (Nτ= 1) and a unique class (K= 1). For multiple classes, the previous hypothesis can no longer be true, whereas the proposal based on the power-normalization of the effective sources still provides an admissible statistical interpretation for the iteration. VII. GAUSSIAN SHRINKAGE LDA Under the hypotheses of p-dimensional Gaussian features for each class, with means µk, k = 1,2, and homoscedastic covariance matrices Σf|c1=Σf|c2=Σf, the LDA classifier considered in the Table I implements the maximum a posteriori (MAP) Bayesian classification [24]. However, when the number of feature vectors fτfor training is not sufficiently large with respect to their dimension p, this method can be prone to 1534-4320 (c) 2018 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2019.2905894, IEEE Transactions on Neural Systems and Rehabilitation Engineering 6 IEEE TRANSACTIONS ON ..., VOL. X, NO. X, AUGUST 2018 overfitting. Moreover, the implementation of the classifier uses the within-class precision matrix of the features ˆ Σ−1 f, which in this situation may be poorly conditioned. To address this problem we should resort to some form of regularization of the averaged within-class covariance ˆ Σf=p(c1)ˆ Σf|c1+p(c2)ˆ Σf|c2.(33) Regularized Discriminant Analysis [25] considered the projection of the sample covariance estimate (in our case, the ˆ Σfdefined previously) onto the identity matrix to obtain hˆ Σf,IiI≡vI, and then estimate the true covariance matrix Σfwith the convex combination ˜ Σf= (1 −ρ)ˆ Σf+ρ(vI).(34) The shrinkage of the sample covariance matrix towards the projection can improve the matrix conditioning and provide a closer estimate to the true covariance matrix for a carefully chosen parameter ρ. The problem consists in finding the optimal value for ρ. Ledoit and Wolf in [22] studied how to automatically determine it by approximately minimizing the minimum quadratic error between the unknown covariance matrix Σfand its shrunken estimation ˜ Σf min ρE  Σf−˜ Σf   2 Fs.t. ˜ Σf= (1 −ρ)ˆ Σf+ρ(vI).(35) The estimator of ρobtained by Ledoit and Wolf is given by ˆρLW =min(PNτ τ=1k(fτ−µkτ)(fτ−µkτ)T−ˆ Σfk2 F N2 τkˆ Σf−(Tr{ˆ Σf}/p)Ik2 F ,1),(36) where µkτrefers to the mean of the class to which the feature fτbelongs. This choice for the estimate guarantees an asymptotically optimal combination of the sample covariance matrix and the identity matrix, is asymptotically consistent and makes no assumption over the data distribution. In the context of MI-BCI, Lotte considered in [7] the Shrunken LDA (sLDA) classification. This method replaces the sample covariance matrix of the features ˆ Σfin Linear Discriminant Analysis with the Ledoit and Wolf regularized covariance matrix ˜ Σffor ρ= ˆρLW. sLDA obtained significant accuracy improvements over standard LDA so its use was highly recommended [3]. Note, however, that the LDA classifier assumes conditional Gaussian classes and, under this assumption, the Ledoit and Wolf regularization technique usually does not provide the best possible mean-square error for finite samples. Chen et al. recognized in [13] that ˆρLW uses statistics of the features up to order four, while under Gaussian hypothesis the mean and covariance condense all the relevant information. They developed an Oracle Approximate Shrinkage (OAS) procedure for small samples that exploits the Gaussian hypothesis. The estimator of ρprovided by the OAS method is ˆρOAS = min        1−2 pTr( ˆ Σ2 f) + Tr2(ˆ Σf) Nτ+1−2 pTr( ˆ Σ2 f)−Tr2(ˆ Σf) p,1       ,(37) and it was shown to attain a better mean square error in simulations than ˆρLW and other alternatives. In what follows we denote by gLDA the implementation of the LDA classifier in combination with the Oracle Approximate Shrinkage estimator of the feature covariance matrix, which is obtained from equation (34) with ρ= ˆρOAS. The added “g” refers to the Gaussian hypothesis of the centered features. According to our experiments in MI-BCI, the classification with gLDA provides relevant gains in accuracy with respect to both the standard LDA and sLDA techniques. VIII. EXPERIMENTAL RESULTS In this section, we will try to corroborate through illustrative simulations the good performance of the proposed covariance estimators that form part of the proposals nCSP and gLDA. The first simulation reveals the expected improvement in the estimation of the covariances with a set of synthetic data since, for its evaluation, the true underlying covariance matrices of the classes have to be known. The remaining simulations consider real datasets from the BCI competitions and test the possible combination of the proposals with state-of-theart techniques. A. Testing the improvement in the estimation of Σx|ck In this experiment, we design a synthetic simulation for corroborating the improvement that can be obtained with the proposed estimation method for the class covariance means. The centroids for the right-hand and left-hand classes have been set equal to the estimated class covariances of user A01 from the dataset IV-2a [29]. We used 25 training trials per class, each with 22 sensors and a length of 500 samples. The samples xτ(t)of each trial τwere drawn from a multidimensional Gaussian density N(0,˜ Cτ), where ˜ Cτwas generated from a local perturbation of the conditional-class mean Σx|ck of the trial. The details of the procedure for the generation of local random covariance matrices in the neighborhood of its centroids are described in Appendix C. The proposed power-normalization technique does not help to guess the absolute scales of the underlying covariance matrix centroids, because these scales are subordinated to the objective of equalizing the power of the effective sources. Fortunately, it is well known that they are irrelevant in the evaluation of the common principal directions. Hence, a good measure of similarity between the true and estimated covariance centroids should be invariant with respect to the scaling of the compared arguments. A natural measure of dissimilarity between covariance matrices with arbitrary scaling is the scaleinvariant version of the Riemmanian distance DR(ˆ Σx|ck,Σx|ck) = min αδR(αˆ Σx|ck,Σx|ck)(38) = min α Nx X i=1 log2λi α! 1 2 (39) = Nx X i=1 log2λi e1 NxPNx j=1 log λj! 1 2 (40) where δR(·,·)denotes the standard Riemmanian distance and λi, i = 1, . . . , Nx, refers to the eigenvalues of ˆ Σ−1 x|ckΣx|ck. 1534-4320 (c) 2018 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2019.2905894, IEEE Transactions on Neural Systems and Rehabilitation Engineering OL´ IAS et al.: EEG PROCESSING WITH IMPROVED COVARIANCE MATRIX ESTIMATORS 7 Classic State-of-the-art: CSP+... Proposed: nCSP+... Dataset User CSP+LDA sLDA RMDM TSLR gLDA (p-value) RMDM (p-value) TSLR (p-value) III-3a k3b 94.31 95.09 94.90 96.13 95.38 (1.7e-02) 94.97 (1.2e-01) 96.31 (2.8e-02) k6b 75.77 77.77 77.23 78.94 78.78 (7.1e-04) 78.15 (7.7e-05) 79.69 (1.0e-03) l1b 86.73 88.18 88.46 89.25 89.48 (1.6e-07) 89.90 (1.4e-13) 90.29 (2.7e-08) −−− mean 85.60 87.01 86.87 88.11 87.88 (4.0e-09) 87.68 (6.9e-13) 88.76 (4.3e-09) III-4a aa 67.68 68.56 63.12 72.37 68.0 (8.6e-01) 63.62 (1.3e-01) 70.31 (1.0e+00) al 96.81 96.5 96.25 96.56 96.18 (9.3e-01) 95.75 (1.0e+00) 96.31 (9.8e-01) av 62.12 62.25 58.31 66.25 63.06 (8.7e-02) 59.75 (4.7e-03) 67.87 (5.1e-03) aw 85.12 83.37 80.62 87.62 82.93 (8.7e-01) 80.75 (1.9e-01) 87.75 (1.9e-01) ay 87.68 90.87 87.62 91.0 89.31 (1.0e+00) 87.68 (2.0e-01) 91.62 (3.7e-02) −−− mean 79.88 80.31 77.18 82.76 79.89 (9.4e-01) 77.51 (5.8e-02) 82.77 (2.3e-01) IV-2a A01 86.97 88.18 88.15 89.02 89.0 (1.3e-06) 88.74 (1.0e-05) 89.23 (3.2e-02) A02 73.15 74.63 76.20 77.65 75.20 (7.3e-03) 74.78 (1.0e+00) 76.15 (1.0e+00) A03 87.34 88.53 88.30 89.75 89.78 (1.1e-13) 90.08 (3.3e-28) 90.60 (6.1e-11) A04 68.77 69.99 70.63 71.36 70.95 (3.0e-05) 70.93 (4.7e-02) 71.38 (2.3e-01) A05 56.89 58.61 58.78 59.27 60.07 (4.9e-07) 60.19 (5.5e-08) 59.82 (1.1e-02) A06 61.41 62.36 62.41 63.32 63.12 (2.6e-03) 62.81 (2.5e-02) 63.26 (8.1e-01) A07 88.75 89.92 90.44 91.09 91.20 (5.4e-17) 91.39 (5.5e-09) 91.70 (9.0e-06) A08 86.40 87.65 86.33 88.32 88.73 (2.0e-10) 88.72 (8.2e-48) 89.18 (2.1e-09) A09 82.70 83.95 82.64 84.25 84.67 (6.7e-05) 84.59 (7.5e-27) 85.26 (2.5e-09) −−− mean 76.93 78.20 78.21 79.34 79.19 (8.8e-39) 79.14 (6.0e-41) 79.62 (6.4e-06) TABLE III EXPECTED USER ACCURACY FOR THE binary MI CLASSIFICATION PROBLEM IN EACH OF THE CONSIDERED DATASETS. THE BEST PERFORMANCES ARE MARKED IN BOLD. ONE CAN OBSERVE THAT IN THE MAJORITY OF THE CASES THE IMPROVEMENTS OBTAINED WHEN COMBINING THE STATE-OF-THE-ART METHODS WITH THE PROPOSED TECHNIQUES CAN BE REGARDED AS STATISTICALLY SIGNIFICANT (p-value < 5e-02). 0 1 2 3 4 5 6 3 4 5 6 Iterations Scale-inv. Riemannian distance Block-based normalization Instantaneous normalization Fig. 1. Variations of the scale-invariant Riemannian distance (between the reference Σx|ckand estimated ˆ Σx|ckcovariance matrices) with respect to the number of iterations of the proposed power-normalization procedures. The solid lines represent average distances while the bars represent the 25% and 75% percentiles. Iteration 0refers to the absence of normalization, iteration 1 coincides with the standard trace-based normalization used in CSP, while the remaining iterations are instances of the proposed normalization. Figure 1 illustrates the improvement in the estimation of the class-conditional covariance matrix means for the block (Section V) and instantaneous (Section VI) power-normalization procedures, when both share the initialization ˆ Σ(0) x=I. The x-axis represents the iteration iat which the covariance matrix estimate ˆ Σ(i) x|ckis evaluated, whereas the y-axis represents the average across classes of the scale-invariant Riemannian distances between Σx|ckand ˆ Σ(i) x|ck. The simulation results confirm the expected improvement of these normalizations with respect to the classical one, which corresponds with the result obtained for iteration 1. Being, in this case, the power-instantaneous equalization method slightly more precise than the block based-implementation. B. Experiments using the BCI competitions datasets This subsection is devoted to the experimental comparison of the proposals on real BCI datasets. The normalization scheme for the estimation of the class-conditional covariance matrices, proposed in Section VI, can be combined with a variety of MI-BCI techniques to improve their accuracy. In particular, we compare the differences of performance, between classical CSP and nCSP (our proposal), when they are used in combination with the following classifiers: LDA, its shrinkage variants sLDA and gLDA, RMDM and TSLR. The Python code for the RMDM and tangent space (TS) implementations can be downloaded from [5]. For TSLR, the Logistic Regression (LR) classifier was implemented according to version 0.19.2 of [26] with its default parameters. The experiments in this section have been carried out using three datasets from BCI competitions. Dataset 3a from BCI competition III [27], which contains 60 EEG channels, three users and four classes of motor imagery movements (MIM); dataset 4a from BCI competition III [28] with 108 EEG channels, four users and two classes of MIM; finally, dataset 2a from BCI competition IV [29] has 22 EEG channels, nine users with two sessions per user and four classes of MIM. Each experiment consists of 40 Monte-Carlo simulations where the whole set of available trials for each session, user and pair of movements is randomly split into testing and training groups. After that, the averaged performance over the test trials is reported. The simulations report the average classification accuracy over all the possible confrontations of pairs of classes (K= 2) for each user. By default, the number 1534-4320 (c) 2018 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2019.2905894, IEEE Transactions on Neural Systems and Rehabilitation Engineering 8 IEEE TRANSACTIONS ON ..., VOL. X, NO. X, AUGUST 2018 of training and testing trials is set to 40 and the number of spatial filters is set 8, except for those cases where a range of these values is specified. In Table III, we show the accuracy results for each subject in each of the three datasets. We also report the mid-p values of one-sided McNemar’s tests of hypotheses [30] for paired data that allows to check whether the proposals have significant advantages in accuracy with respect to their respective stateof-the-art approaches. One can observe in Table III that for two of the datasets the proposal nCSP leads to significant improvements in expected accuracy with respect to classical CSP, whereas, its performance remains equivalent for the dataset III-4a. We also compare the algorithms when the number of training trials varies from 4to 80, while the number of testing trials remains equal to the default value of 40. For this purpose, we have employed the dataset IV-2a. Figures 2(a) and 2(b) represent the improvement of nCSP+gLDA, nCSP+RMDM and nCSP+TSLR with respect to their respective baselines: CSP+sLDA, CSP+RMDM and CSP+TSLR. In both figures, the best performance over the whole range of training trials is obtained for the proposed normalization. Figure 2(b) reveals that the use of nCSP instead of CSP progressively increases the improvement with the number of training trials. In Figure 2(a), the combination of nCSP with gLDA sustains the improvement across the number of training trials. A disaggregated analysis reveals that gLDA improves greatly over sLDA for a small number of training trials. In the last experiment, we analyze the sensitivity of the methods with respect to the chosen number of spatial filters p for dimensionality reduction. Figure 3 enables us to compare the accuracy of the proposals nCSP+gLDA and nCSP+TSLR with respect the existing approaches, for the datasets IV-2a and III-3a. These figures reveal that the standard method CSP+LDA (orange dashed-line) is quite sensitive to the choice p. Its performance attains a maximum at a relatively small value of pand greatly decreases as this number increases. This finding supports the necessity of employing automatic selection techniques to determine the right number of spatial filters for each user [8]. Although, use of Ledoit and Wolf covariance shrinkage estimates (green dashed-line) partially alleviates the previous drawback, the accuracy for the proposed nCSP+gLDA (green solid-line) is more robust with respect to a misspecification of the optimum number of spatial filters. In our simulations, the best performance was obtained for the nCSP procedure in combination with the Tangent Space Logistic Regression (TSLR) classifier (blue continuous-line). This method has outperformed CSP+TSLR (blue dashed-line) in expected user accuracy over all the range of the number of spatial filters and training trials. Similar results have been obtained for multiclass scenarios. We refer the interested reader to the supplementary material [31] that accompanies this manuscript and includes an illustrative Python demo. IX. CONCLUSION In this work, we have studied the problem of obtaining improved covariance matrix estimators for the processing of the MI-BCI signals. We have proposed the application of two techniques that improve the accuracy of these estimations. To counter the inter and intra-trial non-stationarity that hinders the correct estimation of the trial covariance matrices, we propose a power normalization of the EEG source activities. When this is implemented across trials, it improves the classical normalization of the observations used for the EEG trials. Furthermore, the instantaneous power-normalization of the sample source vector seems to enable superior classification results. In this latter case, the proposal extends Tyler’s method (for obtaining an estimate of scatter) to the context of heterogeneous trial observations. The second technique refers to a convenient regularization of the feature covariance matrix of the classifiers. Both proposals are transversal, in the sense that they can be easily combined with the existing MI-BCI algorithms to boost their performance. Experimental tests on several BCI competition datasets reveal that a combination of the proposed techniques with state-of-the-art algorithms for motor-imagery classification provides a significant improvement in the classification results. APPENDIX A. Proof of the formula for the power of the effective sources We start by noting that there is a one to one correspondence between Xτand ˜ Sτ, which is given by ˜ Sτ=ΠA0T˜ Sτ=A0T(A0A0T)−1Xτ.(41) Recalling the invariance of the trace of the product of compatible matrices with respect to cyclic permutations in the matrix positions, i.e., Tr{˜ Sτ˜ ST τ}= Tr{˜ ST τ˜ Sτ}, and using (41) to substitute the value of ˜ Sτin (16), we obtain P˜ Sτ=1 TTr{XT τ(A0A0T)−1Xτ}.(42) As we have seen in equation (12), A0A0Tcoincides with the global average covariance matrix of the observations Σx, hence, we can write without any approximations that P˜ Sτ=1 TTr{XT τΣ−1 xXτ}= Tr{Σ−1 xC(0) Xτ}.(43) B. Equivalence with Tyler’s method for estimation of scatter The algorithmic solution provided by the proposed instantaneous power-normalization technique may be regarded as a variation of Tyler’s method used in statistics for obtaining a robust m-estimator of scatter [23]. As it will be shown, for a single trial (Nτ= 1) and a single class (K= 1), both techniques use complementary arguments to arrive by different paths to a similar final result. To trace back the equivalence, we review the problem considered by Maronna in [32], where he studied how to obtain robust affine-invariant estimates of mean and scatter from a set {x(1), . . . , x(T)}of multivariate i.i.d. samples, drawn from an elliptical distribution. Let the density of x(t)for a given scatter matrix Cxbe p(x(t);Cx)=κ|Cx|−1/2 φ((x(t)−µx)TC−1 x(x(t)−µx)) (44) where φ(·)is an integrable and non-negative function with domain R+and κis the normalization constant. For simplicity, 1534-4320 (c) 2018 IEEE. Translations and content mining are permitted for academic research only. Personal use is also permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TNSRE.2019.2905894, IEEE Transactions on Neural Systems and Rehabilitation Engineering OL´ IAS et al.: EEG PROCESSING WITH IMPROVED COVARIANCE MATRIX ESTIMATORS 9 4 8 12 20 30 40 50 60 70 80 72 74 76 78 80 82 Number of training trials % Accuracy nCSP+gLDA CSP+sLDA CSP + LDA (a) Improvement of nCSP+gLDA over the baselines. 4 8 12 20 30 40 50 60 70 80 72 74 76 78 80 82 Number of training trials % Accuracy nCSP+TSLR CSP+TSLR nCSP+RMDM CSP+RMDM (b) Improvement of nCSP+TSLR and nCSP-RMDM over the baselines. Fig. 2. This experiment shows the accuracy of the binary classification methods with respect to the number of training trials for dataset IV-2a. The arrows in Subfigures (a) and (b) represent the improvement in performance of nCSP+gLDA, nCSP+RMDM and nCSP+TSLR with respect to their baselines. 2 4 6 8 10 12 14 16 18 20 22 76 77 78 79 80 Number of spatial filters (Dataset IV-2a) % Accuracy CSP+TSLR nCSP+TSLR CSP+sLDA nCSP+gLDA CSP + LDA 2 4 6 8 12 16 20 24 28 32 36 40 44 82 83 84 85 86 87 88 89 Number of spatial filters (Dataset III-3a) % Accuracy CSP+TSLR nCSP+TSLR CSP+sLDA nCSP+gLDA CSP + LDA Fig. 3. Variations in performance of the MI-BCI binary classification methods with respect to p, the number of spatial filters. The results confirm the advantages of using nCSP in combination with the state-of-the-art classifiers to improve the expected user accuracy. we assume in this exposition that the mean µxis known (or can be reasonably estimated from the data) and focus on the steps for the estimation of Cx. The normalized log-likelihood of the observations is hlog p(x(t); Cx)it= log κ−1 2log |Cx|+hlog φ(αt)it where αt= (x(t)−µx)TC−1 x(x(t)−µx). In [32] the maximization of the log-likelihood leads to an M-estimator of scatter ˆ Cxthat satisfies the estimating equation ˆ Cx−u(ˆαt)(x(t)−µx)(x(t)−µx)Tt=0(45) where u(αt) = −2dlog φ(αt) dαt. Although there is no close-form solution to this equation because of the coupling between ˆαtand CX, there is a general set of conditions that guarantees its uniqueness (see [32]). Years later, Tyler considered in [23] the same problem. He studied the properties of the specific weighting function u(αt) = Nx/αtand showed that this choice gives the “most robust estimator of the scatter matrix of an elliptical distribution in the sense of minimizing the maximum asymptotic variance”. He also proposed to iteratively solve the estimation equation through a fixed point iteration. In our particular case, xτ(t)≡x(t)−µxand Tyler’s iteration for the estimation of the trial covariance matrices