Electroencephalograph (EEG) signal processing techniques for motor imagery Brain Computer interface systems
Abstract
Brain-Computer Interface (BCI) system provides a channel for the brain to control external devices using electrical activities of the brain without using the peripheral nervous system. These BCI systems are being used in various medical applications, for example controlling a wheelchair and neuroprosthesis devices for the disabled, thereby assisting them in activities of daily living. People suffering from Amyotrophic Lateral Sclerosis (ALS), Multiple Sclerosis and completely locked in are unable to perform any body movements because of the damage of the peripheral nervous system, but their cognitive function is still intact. BCIs operate external devices by acquiring brain signals and converting them to control commands to operate external devices. Motor-imagery (MI) based BCI systems, in particular, are based on the sensory-motor rhythms which are generated by the imagination of body limbs. These signals can be decoded as control commands in BCI application. Electroencephalogram (EEG) is commonly used for BCI applications because it is non-invasive. The main challenges of decoding the EEG signal are because it is non-stationary and has a low spatial resolution. The common spatial pattern algorithm is considered to be the most effective technique for discrimination of spatial filter but is easily affected by the presence of outliers. Therefore, a robust algorithm is required for extraction of discriminative features from the motor imagery EEG signals. This thesis mainly aims in developing robust spatial filtering criteria which are effective for classification of MI movements. We have proposed two approaches for the robust classification of MI movements. The first approach is for the classification of multiclass MI movements based on the thinICA (Independent Component Analysis) and mCSP (multiclass Common Spatial Pattern Filter) method. The observed results indicate that these approaches can be a step towards the development of robust feature extraction for MI-based BCI system. The main contribution of the thesis is the second criterion, which is based on Alpha- Beta logarithmic-determinant divergence for the classification of two class MI movements. A detailed study has been done by obtaining a link between the AB log det divergence and CSP criterion. We propose a scaling parameter to enable a similar way for selecting the respective filters like the CSP algorithm. Additionally, the optimization of the gradient of AB log-det divergence for this application was also performed. The Sub-ABLD (Subspace Alpha-Beta Log-Det divergence) algorithm is proposed for the discrimination of two class MI movements. The robustness of this algorithm is tested with both the simulated and real data from BCI competition dataset. Finally, the resulting performances of the proposed algorithms have been favorably compared with other existing algorithms.
Full text
Electroencephalograph (EEG) Signal Processing Techniques for Motor Imagery Brain Computer Interface Systems Submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy by THIYAM DEEPA BEETA June, 2018
DECLARATION I here by declare that the thesis entitled “EEG Signal Processing Techniques for Motor Imagery Brain Computer Interface Systems” submitted by me, for the award of the degree of Doctor of Philosophy to VIT University, is a record of bonafide work carried out by me under the supervision of Dr. E. R. Rajkumar, Associate Professor, SENSE, VIT University, Vellore (currently working as Senior Architect, Robert Bosch Engineering and Business Solutions Private Limited, Bangalore) and Dr. Sergio A. Cruces Alv ´ arez, Associate Professor, Signal Processing Group, University of Seville, Spain. I further declare that the work reported in this thesis has not been submitted and will not be submitted, either in part or in full, for the award of any other degree or diploma in these institutes or any other institute or university. Place: Vellore Date: Signature of the Candidate 11/06/2018
ABSTRACT Brain-Computer Interface (BCI) system provides a channel for the brain to control external devices using electrical activities of the brain without using the peripheral nervous system. These BCI systems are being used in various medical applications, for example controlling a wheelchair and neuroprosthesis devices for the disabled, thereby assisting them in activities of daily living. People suffering from Amyotrophic Lateral Sclerosis (ALS), Multiple Sclerosis and completely locked in are unable to perform any body movements because of the damage of the peripheral nervous system, but their cognitive function is still intact. BCIs operate external devices by acquiring brain signals and converting them to control commands to operate external devices. Motor-imagery (MI) based BCI systems, in particular, are based on the sensory motor rhythms which are generated by the imagination of body limbs. These signals can be decoded as control commands in BCI application. Electroencephalogram (EEG) is commonly used for BCI applications because it is non-invasive. The main challenges of decoding the EEG signal are because it is non-stationary and has low spatial resolution. The common spatial pattern algorithm is considered to be the most effective technique for discrimination of spatial filter but is easily affected by the presence of outliers. Therefore, a robust algorithm is required for extraction of discriminative features from the motor imagery EEG signals. This thesis mainly aims in developing robust spatial filtering criteria which are effective for classification of MI movements. We have proposed two approaches for the robust classification of MI movements. The first approach is for the classification of multiclass MI movements based on the thinICA (Independent Component Analysis) and mCSP (multiclass Common Spatial Pattern Filter) method. The observed results indicate that these approaches can be a step towards the development of a robust feature extraction for MI based BCI system. The main contribution of the thesis is the second criterion, which is based on Alpha- Beta logarithmic-determinant divergence for classification of two class MI movements. A detailed study has been done by obtaining a link between the AB log det divergence and CSP criterion. We propose a scaling parameter to enable similar way for se- i
lecting the respective filters like the CSP algorithm. Additionally, the optimization of the gradient of AB log-det divergence for this application was also performed. The Sub-ABLD (Subspace Alpha-Beta Log-Det divergence) algorithm is proposed for the discrimination of two class MI movements. The robustness of this algorithm is tested with both the simulated and real data from BCI competition dataset. Finally, the resulting performances of the proposed algorithms have been favourably compared with other existing algorithms. ii
ACKNOWLEDGEMENT I am very grateful to my supervisor Dr. E. R. Rajkumar , Associate Professor, SENSE, VIT University (currently working as Senior Architect, Robert Bosch Engineering and Business Solutions Private Limited, Bangalore) for being a constant support throughout my PhD and motivating me to take up challenges. I am also thankful to him for giving me many wonderful opportunities. I owe my sincere gratitude to my co-supervisor Dr. Sergio A. Cruces , Associate Professor, Signal Processing Group, University of Seville for introducing me to this interesting field of signal processing. I am very grateful to him for his efforts to teach me right from scratch till the very end and also for being there to answer all my questions. This work would not have been possible without their support. I would like to thank the Chancellor of VIT University, Vellore, the Rector of University of Seville, Spain and the management of both the universities for allowing me to carry out this co-directed PhD and also for providing me infrastructural facilities and resources needed for my research work. I am also thankful to the Vice chancellor , Dean Academic Research , VIT University, Vellore for all the administrative help and support. I take this opportunity to express my thanks to VIT University , Erasmus Mundus European Commission and the Spanish Government for providing funds for my research work. I would also like to acknowlegde my doctoral committee members, Dr. Venkatesh B. , Professor, Department of Engineering Design, IIT madras, Chennai, Dr. Sriraam N. , Head of Deparment, Department of Medical Electronics, M. S. Ramaiah Institute of Technology, Bangalore and Dr. Arulmozhivarman P. , Dean, School of Electrical Engineering, VIT University, Vellore for their kind and timely suggestions to improve the quality of my research work. I am indebted to my friends Rucha , Poulami , Lalitha and Cassandra for always lending a helping hand in times of my need. I would also like to extend my thanks to my other friends, colleagues and room-mates who have always been there for me, especially during my hard times. My deepest gratitude to my colleagues Juan Antonio , Irene Santos , Sunny Dayal , Maria Jose , Irene Fondon , Pablo Anguilera , Uxi at the University of Seville, iii
for making my stay wonderful and for always being there for me. My stay in Seville would have been difficult without their support. A special thanks goes to Javier for helping me with the experimental studies. My heartfelt gratitude goes to my family and to Amarjit , for their constant encouragement and moral support along with patience and understanding throughout my work. Lastly, I thank God and my heavenly Father for supporting and blessing me to complete this degree. Place: Vellore Date: THIYAM DEEPA BEETA iv 11/06/2018
TABLE OF CONTENTS ABSTRACT .................................... i ACKNOWLEDGEMENT ............................. iii LIST OF FIGURES ................................ ix LIST OF TABLES ................................. xiii LIST OF TERMS AND ABBREVIATIONS ................... xiv LIST OF NOTATIONS .............................. xvi 1 Introduction 2 1.1 Overview .................................. 2 1.1.1 Objective of the Thesis . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.2 Structure of the Thesis . . . . . . . . . . . . . . . . . . . . . . . 3 2 Brain Computer Interface Background 6 2.1 StructureofTheBrain............................ 6 2.2 MotorControl................................ 7 2.3 EEGSignalAcquisition........................... 8 2.3.1 EEGElectrodes........................... 9 2.3.2 EEGRhythm ............................ 10 2.3.3 EEGArtifacts............................ 12 2.3.4 Event Related Synchronization and Desynchronization (ERS/ERD) 13 2.4 MotorImageryBCI............................. 13 2.4.1 Signal Processing Techniques . . . . . . . . . . . . . . . . . . . . 15 2.4.2 Classification ............................ 16 2.5 Conclusions................................. 18 3 Spatial Filtering Methods 20 3.1 Common Spatial Pattern Algorithm . . . . . . . . . . . . . . . . . . . . . 20 v
3.2 Divergence Based CSP Approaches . . . . . . . . . . . . . . . . . . . . . 21 3.3 The information theoretic feature extraction framework . . . . . . . . . . 23 3.4 ProbabilisticCSP .............................. 24 3.5 OtherCSPVariants ............................. 26 3.6 Conclusions................................. 30 4 Study of Other CSP Based Approaches 32 4.1 Simplification of the CSP Objective Function . . . . . . . . . . . . . . . 32 4.2 DiscriminativeCSP............................. 33 4.3 ICACorrectionstoCSP........................... 33 4.4 Experimental dataset and study . . . . . . . . . . . . . . . . . . . . . . . 34 4.5 Results.................................... 35 4.6 Conclusions................................. 35 5 Optimization of Thin Independent Component Analysis-Common Spatial Pattern (ThinICA-CSP) Algorithm 37 5.1 Blind Source Separation Techniques and Its Background . . . . . . . . . 37 5.1.1 Principal Component Analysis (PCA) . . . . . . . . . . . . . . . 38 5.1.2 Independent Component Analysis (ICA) . . . . . . . . . . . . . . 40 5.2 RelatedWork ................................ 43 5.2.1 MulticlassCSP ........................... 44 5.3 Implementation of the Discrimination of the MI-EEG Signals . . . . . . . 44 5.3.1 Experimental Dataset . . . . . . . . . . . . . . . . . . . . . . . . 45 5.3.2 Preprocessing ............................ 45 5.3.3 The Thin ICA-CSP Criterion and its Implementation . . . . . . . 45 5.3.4 Feature Extraction . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.3.5 Classification ............................ 49 5.4 Results.................................... 50 5.5 Conclusions................................. 51 6 Divergence maximization and its Relation with CSP 53 6.1 Background ................................. 54 6.2 Notation and Model of the Measurements . . . . . . . . . . . . . . . . . 55 6.3 The Common Spatial Patterns Algorithm . . . . . . . . . . . . . . . . . . 56 vi
LIST OF TABLES 2.1 Brainrhythms............................... 11 4.1 Performance comparison between CSP and discriminative CSP . . . . . 35 4.2 Performance comparison of ICA corrections to CSP for different values of ⌘=[0,0.5,1] .............................. 35 5.1 Comparative accuracy using LDA and SVM classifier with mCSP and ThinICA-CSP............................... 50 xiii
LIST OF TERMS AND ABBREVIATIONS AB Log-Det Alpha-Beta Logarithmic-Determinant divergence . . . . . . . . . . . 3 ALS Amyotrophic Lateral Sclerosis . . . . . . . . . . . . . . . . . . . . . . . . . 2 BCI Brain-Computer Interfacing . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 BSE BlindSourceExtraction............................ 38 BSS BlindSourceSeparation............................ 15 cCSP Composite Common Spatial Pattern . . . . . . . . . . . . . . . . . . . . . 26 CSP CommonSpatialPattern.............................3 CSSP Common Spatio-Spectral Pattern . . . . . . . . . . . . . . . . . . . . . . 30 CSSSP Common Sparse Spectral-Spatial Pattern . . . . . . . . . . . . . . . . . 30 CV Cross-Validation ................................ 89 ECG Electrocardiograph .............................. 12 EEG Electroencephalograph.............................2 EMG Electromyograph............................... 12 EOG Electrooculograph .............................. 12 FastICA Fast Independent Component Analysis . . . . . . . . . . . . . . . . . . 40 FBCSP FilterBankCSP .............................. 30 i.i.d. independent and identically distributed . . . . . . . . . . . . . . . . . . . . 55 ICA Independent Component Analysis . . . . . . . . . . . . . . . . . . . . . . . 13 InfoMax Information Maximization . . . . . . . . . . . . . . . . . . . . . . . . 40 JADE Joint Approximation Diagonalization of Eigenmatrices . . . . . . . . . . 42 LDA Linear Discriminant Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 17 xiv
MAP-CSP Maximum-a-Posteriori CSP . . . . . . . . . . . . . . . . . . . . . . 25 MEMS Microelectromechanical systems . . . . . . . . . . . . . . . . . . . . . 10 MI MotorImagery ..................................2 MKL Multiple Kernel Learning . . . . . . . . . . . . . . . . . . . . . . . . . . 82 MRP Movement Related Potential . . . . . . . . . . . . . . . . . . . . . . . . . 13 PCA Principal Component Analysis . . . . . . . . . . . . . . . . . . . . . . . . 13 PMA PremotorCortexArea .............................7 PMC PrimaryMotorCortex .............................7 sCSP stationaryCSP ................................ 28 SMA SupplementaryMotorArea ..........................7 SMR Sensorimotorrhythm............................. 14 SOBI Second Order Blind Identification . . . . . . . . . . . . . . . . . . . . . . 40 SPD Symmetric and Positive Definite . . . . . . . . . . . . . . . . . . . . . . . 53 SSEP Steady State Evoked Potentials . . . . . . . . . . . . . . . . . . . . . . . 14 Sub-ABLD Subspace-Alpha Beta Log-Det Algorithm . . . . . . . . . . . . . . . 3 SVM SupportVectorMachine ........................... 17 ThinICA Thin Independent Component Analysis . . . . . . . . . . . . . . . . . . 3 VB-CSP Variational Bayesian CSP . . . . . . . . . . . . . . . . . . . . . . . . . 25 xv
LIST OF NOTATIONS J(w)Criterion ................................... 15 Cov1Average covariance matrix of class 1.................... 15 XEEGobservationmatrix ............................. 15 DBh Bhattacharyyadistance ............................ 23 DGammadivergence ............................... 23 H(C)Entropy of variable C............................ 23 AMixingmatrix .................................. 16 BUnmixingmatrix................................. 16 IIdentitymatrix................................... 22 PPenaltyterm ................................... 22 SSourcematrix................................... 16 TWhitening transform matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 UUnitary matrix that represents eigenvectors . . . . . . . . . . . . . . . . . . . 39 YEstimatedEEGsignal .............................. 16 Eigenvaluesmatrix................................ 39 Cov2Average covariance matrix of class 2.................... 15 ROrthogonalmatrix ................................ 22 WMatrixofspatialfilters ............................. 21 wSpatialfilter ................................... 15 ⌘Regularizingparameter.............................. 26 Scalingfactor................................... 68 Eigenvalue .................................... 21 ⌧Timedelay .................................... 42 ˜ CovcEstimated covariance matrix of class c.................... 26 xvi
cMotorimageryclass................................ 21 cum(.)Cumulants.................................. 42 kurt(.)Kurtosis................................... 24 log(.)Logarithmicfunction............................. 21 pTotal number of filters selected . . . . . . . . . . . . . . . . . . . . . . . . . . 21 p(c)Probability of class c.............................. 55 var(.)Variancefunction .............................. 49 DKL(i, j)KL divergence between two distributions . . . . . . . . . . . . . . . . 26 DsKL(.||.)Symmetric KL divergence between two class covariance matrices . . . 21 Div(.||.)Beta divergence between two distributions . . . . . . . . . . . . . . . 22 DivsKL(.||.)Symmetric KL divergence between two distributions . . . . . . . . 21 I(.;.)Mutual information between two variables . . . . . . . . . . . . . . . . . 23 Contrast function of ThinICA-CSP algorithm . . . . . . . . . . . . . . . . . . 47 ⇥Contrastfunction................................ 47 ¯ xMean of x..................................... 16 GcGeneric matrix of class c............................ 26 RxCorrelation matrix of vector x......................... 41 ⇤Eigenvalues matrix of the pencil matrix . . . . . . . . . . . . . . . . . . . . . 61 ⌦Semi-orthogonalmatrix ............................. 56 sn(t)Non-stationarysource............................. 27 ss(t)Stationarysource ............................... 27 x(t)Observationvector............................... 16 ˆ USemi-orthogonalmatrix ............................. 46 d·eOperator to round the value to the nearest higher integer . . . . . . . . . . . 55 RSetofrealnumbers................................ 21 N(0,I)Normal distribution with zero mean . . . . . . . . . . . . . . . . . . . . 27 N(ˆµs i,ˆ Covs i)Distribution of the stationary sources . . . . . . . . . . . . . . . . 27 xvii
p(yi|c1)Conditional probability distribution for class 1.............. 21 p(yi|c2)Conditional probability distribution for class 2.............. 21 b·cOperator to round the value to the nearest lower integer . . . . . . . . . . . 55 rˆ U[q]Gradient function of ThinICA-CSP . . . . . . . . . . . . . . . . . . . . 47 |.|+Non negative truncation operator . . . . . . . . . . . . . . . . . . . . . . . . 28 ||.||1l1norm..................................... 28 ||.||plpnorm..................................... 28 ||diag(.)||2Sum of squares of diagonal elements . . . . . . . . . . . . . . . . . . 42 xviii
1
CHAPTER 1 Introduction 1.1 Overview The Brain-Computer Interfacing (BCI) system aims at building the bridge between the brain and the computer. The brain produces electrical and magnetic signals while performing different functions. These signals can be recorded using different techniques. Electroencephalograph (EEG) is the most commonly used method for measuring the electrical activity of the brain in BCI applications because of its non-invasive characteristics. BCI system enables a disabled person to operate other assistive devices like a wheelchair or robotic arm by using brain signals. The Motor Imagery (MI) based BCI systems are considered as the most preferable BCI systems. They use MI EEG signals as control commands for external devices without using the peripheral nervous system. Such systems can potentially serve as a communication aid for the people suffering from Amyotrophic Lateral Sclerosis (ALS), Multiple Sclerosis and completely locked-in. The main difference and advantage of BCIs over other assistive devices are the non requirement of any form of muscle control. Unfortunately, the performance accuracy of current BCIs is still very low restricting to use them out of the laboratories. The main reason for the hinder of performance is due to the non-stationary nature of the EEG signals. As a result of this, the signal properties not only changes from person to person but also from trial to trial which gives more challenges in data analysis. In addition to this difficulty, the presence of artifacts such as eye movements, muscle activities and improper placement of electrodes added more challenges in the EEG signal processing. Furthermore, the performance of the BCI system also decreases with the increased in the number of motor imagery movements. Although, the artifacts and non-stationarity cannot be fully removed, a robust signal processing algorithm can be used for better signal analysis and high classification accuracies. Several approaches have been proposed for designing a robust signal processing unit but still, the gap is large to highlight the BCI system for real time applications. Therefore, this motivates the necessity of developing a more robust algorithm for MI movements classification. 2
1.1.1 Objective of the Thesis This thesis focuses on the study of different spatial filtering methods used for classification of MI movements and developing robust algorithms based on it for the application of EEG based MI-BCI system. The first objective is to understand and analyze the various filtering methods in this area. Another objective can be divided into two sections. The first section is to propose a new criterion for the classification of multiclass motor imagery movements. The second section, which is the main contribution of this thesis, will formulate a new criterion based on the Alpha-Beta Logarithmic-Determinant divergence (AB Log-Det) for discrimination of two class MI movements. The list of contributions is the following: • A novel criterion for classification of multiple MI movements is proposed, it is based on the extension of the Thin Independent Component Analysis (ThinICA) method for blind source separation. • ThinICA-Common Spatial Pattern (CSP) algorithm is proposed combining the CSP and the extension of the ThinICA method for multiple MI movements classification. •The relation between CSP and AB Log-Det divergence is determined. • The scaling factor , which provides the equal solution between the CSP and AB Log-Det divergence is obtained. •A novel regularized criterion is proposed based on AB Log-Det divergence. • The optimization of AB Log-Det divergence with proper gradient function is presented • The Subspace-Alpha Beta Log-Det Algorithm (Sub-ABLD) algorithm is proposed to address the robust features extraction problem in MI-BCI systems. 1.1.2 Structure of the Thesis This thesis consists of nine chapters and four appendices. This chapter presents the introduction, motivation, problem statement, objectives and structure of the thesis. In Chapter 2, the basic anatomical and physiological details of the human brain, different types of brain rhythms and the various types of EEG artifacts are introduced. The BCI and MI-based BCI as well as the commonly used filtering and classification algorithms for the discrimination of motor imagery EEG signals are also discussed. 3
In Chapter 3, the existing effective spatial filtering approaches are reviewed in details. The popular CSP algorithm, variants of CSP and other different spatial filtering approaches for discrimination of motor imagery movements are described. In Chapter 4, the simplification of CSP objective function as Rayleigh quotient is shown. Moreover, different CSP based approaches such as discriminative CSP and ICA corrections to CSP based on the existing techniques are presented. In Chapter 5, the extension of existing Thin-ICA criterion is presented. The maximization of the proposed criterion is done for the classification of multiple class MI movements. The Thin-ICA CSP algorithm combines the proposed criterion with the solution of the multiclass CSP algorithm. The performance results using BCI competition dataset are also presented. In Chapter 6, the field of AB Log-Det divergence is introduced for BCI applications. The optimization of this divergence is performed and its interpretation of CSP is obtained. The scaling parameter that provides the equivalent solution between the AB Log-Det divergence and CSP is presented. The gradient of AB Log-Det divergence is computed and validated. The robustness of this divergence based on ↵and is shown. In Chapter 7, the criterion based on AB Log-Det divergence is proposed for addressing the problem for discrimination of two class MI movements. Sub-ABLD algorithm is proposed by optimizing this criterion. This proposed criterion considers both the within class and between class divergence. The algorithm is evaluated using both artificial and BCI competition datasets. In Chapter 8, the simulations of the both the proposed algorithms i.e. ThinICA-CSP and Sub-ABLD algorithm are compared with the other baseline algorithms. The study of the performance of Sub-ABLD algorithm in different scenarios are also presented and the performance results are discussed. In Chapter 9 gives the final conclusion of this thesis and some future research ideas are presented. 4
Table 2.1 Brain rhythms Name Band (Hz) Characteristics Location Waveform Delta () 0.14 Very low frequency waves. For adults observed at the time of deep sleep. Common in infants and children observed during wakefulness frontally in adults, posteriorly in children Theta (✓) 48 Present in adults during drowsiness or an idling state and also normally observed in young children various locations that are not involved in any apparent task Alpha (↵) 813 Caused by closing eyes, relaxation and attenuate as one becomes involved in some mental task posterior regions, occipital and temporal cortex Mu (µ) 8 - 13 Effected by actual movement, motor imagery or stimulation sensorimotor cortex Beta ()1330 Correlated with active thinking, focus, stress, and an alert state frontal regions, somatosensory cortex Gamma () >30 Present during the highly attentive states of consciousness and perception which involves higher mental activity various locations 11
Beta rhythm is generated when planning to execute a movement. Therefore, these two rhythms are most important for MI-based BCI system. 2.3.3 EEG Artifacts Signals which are not generated by the brain but are visible in the EEG recording are termed as EEG artifacts. These artifacts interfere with the analysis of the EEG and hence, the knowledge of their source and nature is essential for their removal (Urig ¨ uen and Garcia-Zapirain, 2015). Typically, the EEG signals are 10 to 100 mV in amplitude. The main types of EEG artifacts are physiological and non-physiological. Physiological artifacts originate from the other regions of the subject’s own body and are often closer in magnitude to the EEG signals than non-physiological artifacts. 2.3.3.1 Physiological Artifacts: The different types of physiological artifacts present in the EEG signals are listed below: Electromyograph (EMG) : EMG is the most common type of artifacts that is picked up during EEG acquisition, due to movement or muscle twitching. However, since its duration is smaller and morphology distinct, it is easily separable except in cases such as Parkinsons or Huntingtons diseases where it is not so easily identifiable due to the high frequency of occurrence. Electrooculograph (EOG) (Urig ¨ uen and Garcia-Zapirain, 2015): Eye movements are generally visible as a part of EEG. They occur due to blinking or the axial movement of the eyeballs. Yet, these movements can be beneficial when dealing with sleep EEG where eye movements can be corresponding to the nature of the EEG waves. Tongue movements : Similar with the eyes, the tongue also can interfere with the EEG signals due to chewing, biting and other tongue movements, especially those involving the tip of the tongue. These artifacts can be removed by recording these movements separately and then eliminating them from the recorded signals. Electrocardiograph (ECG) : ECG artifacts are observed in EEG signals when there are electric changes in the heart, which occurs either because the subject is obese, has heart disease or when the inter-electrode distances are high. Pulse : Pulse interferes with EEG signals when the electrodes are placed directly on the surface of a pulsating vessel on the skull. 2.3.3.2 Non-Physiological Artifacts: In addition to the physiological artifacts, non-physiological artifacts also affect the EEG signals. The common non-physiological artifacts present in EEG signals are described below. 12
Power-line interference : The interference caused by the frequency of the electric supply can be eliminated as the frequency is specific and known beforehand. This changes the baseline of the EEG signal, hence it needs to be removed. Artifacts due to electrodes : These may occur due to the sudden movement or popping of an electrode. This is easily identified on the EEG signal as a sudden vertical transient. The improper contact of the electrodes with scalp can also lead to the distortion of the EEG signals Environmental conditions : Movements of other people in the surroundings may be the cause of artifacts in the EEG signals. Electrical interference from other devices may produce artifacts. Respirators and other equipment can also alter the EEG recordings. Some of the artifacts can be easily removed by filtering the EEG signals with appropriate filters whereas other artifacts are removed by using other signal decomposition techniques such as Independent Component Analysis (ICA), Principal Component Analysis (PCA) etc. 2.3.4 Event Related Synchronization and Desynchronization (ERS/ERD) During movement, two kinds of signals are generated. The low-frequency signal generated in preparation of movement is termed as Movement Related Potential (MRP). The high frequency signals are termed as Event-Related Synchronization/ De-synchronization (ERS/ERD) (Pfurtscheller et al., 2006). These signals are not phase locked, unlike MRP. However, the neural generators for both signals are expected to be different. ERD is a decrease in the oscillation frequency while ERS is an increase in the oscillation frequency in the ongoing pattern of the EEG, in response to an induced event. Such EEG signals are processed in the frequency domain and their effects are expressed as relative increment or decrement in power. Furthermore, the effects of ERS and ERD are varied in different frequency ranges. For example, motor events lead to ERD in the low-frequency bands while the same task may lead to ERS in the higher frequency bands. Fig. 2.4 depicts the occurrence of ERS and ERD within a regular EEG signal. The corresponding increase and decrease in relative power can be observed. 2.4 Motor Imagery BCI BCI can be divided into various categories based on the functionality and the types of the trigger signals. Traditional BCI system can be divided into dependent and independent as well as synchronous and asynchronous. The dependent BCI requires a certain level of motor control from the subjects to operate the BCI system whereas independent BCI does not need any motor control. Another type is synchronous and asynchronous. Synchronous BCI allows the user to operate only during a certain time on the contrary asynchronous BCI can be operated at any time. Furthermore, the BCI can 13
Fig. 2.4 ERS and ERD within a typical EEG signal Source:http://www.bbci.de also be divided based on the types of the input trigger brain signals such as Steady State Evoked Potentials (SSEP), P300, Sensorimotor rhythm (SMR), slow cortical potentials and motor imagery. The SSEP and P300 are generated subconsciously when the subject receives the external stimulus. The sensorimotor rhythm and the motor imagery rhythms are generated when the subject performs the actual or imagery movement of the body limbs. MI-based BCI uses the motor imagery signal which is generated when the users imagine the movement of the body limbs. The steps involved in the MI-based BCI system are summarized in Fig. 2.5. The MI-based BCI is of the synchronous type, where the subject can perform the operation only at a fixed time window. The calibration process starts with the execution of the particular motor imagery movement based on the cue presented in the fixed time window. The subject is asked to repeat this process for several trials. A n number of epoch are extracted and aligned at the start of the cue of each trial from the recorded EEG signals. This extracted epoch of signals is used for further preprocessing. The next step is to select a particular time interval of the acquired signal that represents the ERD/ERS effect. After selecting the particular time interval of the signal, the next step is spectral filtering. The main aim of this is to concentrate on the SMR modulation which is significant in mu and beta rhythm. The filtered signals are 14
used for the computation of the spatial filter. Motor Imagery Movements Preprocessing Feature extraction Classification Control Commands Feedback Fig. 2.5 MI based BCI system Several spatial filtering approaches have been proposed. The commonly used approaches are presented below. 2.4.1 Signal Processing Techniques •Common Spatial Pattern (CSP) (Ramoser et al., 2000; Blankertz et al., 2008): The CSP algorithm computes the spatial filters by maximizing the variance of one class and at the same time minimizing the variance of the other. The obtained spatial filters are used to discriminate the two MI movements. The solution of the CSP is obtained by solving the eigenvalue decomposition. The CSP criterion J(w)can be represented as min w\max wJ(w)=wTCov1w wTCov2w,(2.1) where Cov1 and Cov2 are the covariance matrices of class 1 and class 2 and w is the spatial filter. The equal numbers of the largest and the smallest eigenvalues are selected and the corresponding eigenvectors represent the set of discriminative spatial filters. •Independent Component Analysis (ICA) (Hyv ¨ arinen et al., 2004): ICA is one of the most popular Blind Source Separation (BSS) techniques. Let us consider the EEG observation X which can be denoted as a linear combination of the 15
independent sources Sand the mixing matrix A X=AS.(2.2) The estimated EEG signal Yis given by Y=BX,(2.3) where B is a unmixing matrix and B=A1 . The main aim of ICA is to obtain the unmixing matrix using the statistical information of the observation. The observed EEG signals can be considered as mixing source signals from different regions of the brain. ICA can be used for separating the individual source from the mixing model as well as it can be used for removing the artifacts such as eye movements from the acquired EEG signals. Therefore, ICA keeps only the relevant information which enhances the signal-to-noise ratio of the signal. •Principal Component Analysis (PCA) (Hotelling, 1933): It is used for source extraction as well as dimensionality reduction. The main objective of PCA is to perform a linear transformation of the observation into a set of new components which are known as principal components with less dimension. The constrain of the linear transformation is that the first principal component has to have the largest variance followed by the remaining components. The same is followed for the second principal component. This transformation gives the principal components that are uncorrelated with each other. By performing PCA, the input data are projected into a space of eigenvectors. The eigenvectors are computed using the covariance matrix of the input signal x(t) . The covariance matrix Cov is given by Cov = T X i=1 (xi¯ x)(xi¯ x)T,(2.4) where ¯ xis the mean which is obtained by ¯ x=1 T T X i=1 xi.(2.5) The ranking of the eigenvectors is done based on the eigenvalues. The eigenvectors with the highest eigenvalues are taken as the first principal component. 2.4.2 Classification The next step is to translate the extracted features into control commands. This can be done using classification algorithms. The classifier being used in the field of BCI are 16
linear classifier (Pfurtscheller, 1999), non-linear classifier (Rezaei et al., 2006), neural networks (Hiraiwa et al., 1990), nearest neighbor classifier (Blankertz et al., 2002) and the combination of classifiers (Pfurtscheller et al., 1993). Among these classifiers, the linear classifier is commonly used for MI discrimination. Two popular linear classifiers for the motor imagery signals are the Linear Discriminant Analysis (LDA) and Support Vector Machine (SVM). •Linear Discriminant Analysis (LDA) :The main aim of LDA (Duda and Hart, 1973) is to project the multidimensional data into a reduced dimensional subspace with higher class separability. LDA approach mainly considers the data for each class as a model of probability density functions. The class of the input data is determined by the larger value of probability density function from the others. LDA assumes that all the classes have a normal distribution and have the same covariance matrix. Let us consider there are c classes and x=[x1,...,x n]T be the samples to be classified, where n represents the no. of samples. The mean, ¯x and the global covariance matrix Cov , can be represented as: ¯xc=1 nk c X i=1 xi(2.6) Covc=1 nc c X i=1 (xt¯xi)(xt¯xi)T.(2.7) Then, the classification of data point xis done by g(x)=argmax cxtCov1 c¯xc1 2¯xT cCovc1¯xc(2.8) which decision boundary is a linear function. The class of x is determined by the objective function given in Eqn. 2.8. The LDA is mainly used for binary classification, but it can also be used for multiclass problem. •Support Vector Machine (SVM) : The SVM (Cortes and Vapnik, 1995) approach has wide applications in the fields of machine learning and pattern recognition. It has the ability to deal with high dimensional and non-linear data. The SVM with kernel can generate non-linear decision boundaries which makes it suitable for discriminating non-linearly separable data. In this method, the data was mapped into a high-dimensional space where the data is spread in such a way that a linear hyper-plane can be fitted. The decision function for kernel-based SVM can be defined by: g(x)=sgn(CovT i=1↵icik(x, xi)+b)(2.9) 17
where x=[1,...,t] is the set of training samples, c represents the class labels, ↵i0 is a Lagrangian multiplier which is a solution of the quadratic optimization problem, k(x, xi) represents the kernel and b is the bias. The selection of kernel and setting of the hyperparameters value are important steps in designing SVM classifier. For the experiment in this study, Gaussian RBF is selected which can be defined as k(x, xi)=e||xxi||2(2.10) where =1 22>0 , controls the width of the Gaussian function, ||xxi|| is the norm of x. Moreover, SVM is insensitive to overtraining which makes it suitable for various applications. 2.5 Conclusions This chapter presents the related background of the brain computer interface system. In the first section, the anatomy and physiology of the human brain involved in motor control were discussed. In the following section, EEG acquisition, EEG rhythms, EEG artifacts and ERS/ERD have been explained. The final section dealt with the steps involved in MI-based BCI. The commonly used signal processing techniques and classifiers used for motor imagery signals classification are also presented in this chapter. 18
19
CHAPTER 3 Spatial Filtering Methods Signal processing techniques mainly aim to de-noise the noisy observed signals and enhance the signal to noise ratio in order to extract the relevant/important information from the signals. Generally, EEG signals are very noisy and easily affected by movement of eyes and muscles. However, it is difficult to remove these artifacts without losing relevant information. Furthermore, it is necessary to filter unrelated brain activity and retain only signals of interest. Therefore preprocessing of the signal is necessary before extracting the required features. The preprocessing of motor imagery EEG signals for BCI applications is mainly done by bandpass filtering the signal together with other filtering techniques. The spatial filtering technique is considered to be quite effective for discrimination of motor imagery EEG signals. The motivation of this chapter is to study the existing spatial filtering and other filtering techniques for this application. The CSP algorithm, regularized CSP algorithms and other variants of CSP are discussed in this chapter. 3.1 Common Spatial Pattern Algorithm The CSP algorithm was first presented as a feature classification algorithm in (Fukunaga and KoonTz, 1970). Initially, it was used for detection of abnormalities in clinical EEG (Koles, 1991) and also used for the classification between normal and abnormal EEG (Koles et al., 1994). Later, it was used for the discrimination of two class movements from the single trial EEG (M ¨ uller-Gerking et al., 1999). A similar study had been performed for discrimination of motor imagery movements from multichannel EEG data. Let us consider the EEG signals x(t)=[x1(t),x 2(t),...,x n(t)]T recorded using n channels during the left and right hand MI movements. The main objective of CSP algorithm is to compute spatial filters that discriminate the two MI movements by maximizing the ratio of the variance of signals between the two classes. The CSP objective function can be considered as a Rayleigh quotient maximization problem as: min w\max wJ(w)=wTCov1w wTCov2w(3.1) 20
of the mixing matrix A and the sources which consists of the stationary ss(t) and non stationary sn(t)components x(t)=As(t)=[AsAn]"ss sn#.(3.33) The goal of SSA is to obtain the estimate separation matrix B=RT that separates the stationary sources, here T is the whitening matrix and R is an orthogonal matrix. The optimization is performed by minimizing the objective function in each step which is given by JB(R)= T X i=1 DKL[N(ˆµs i,ˆ Covs i)||N(0,I)],(3.34) where ˆµs i and ˆ Covs i is the mean and covariance of the i -th epoch, N(ˆµs i,ˆ Covs i) is the distribution of the stationary sources in each epoch and N(0,I) represents the normal distribution. The extracted stationary part of the signal is used for the computation of CSP filter (Von B ¨ unau et al., 2010). The extension of this approach known as groupSSA was proposed by Samek et al. (Samek et al., 2011) considering a group of trials and computing the stationary components from each group, which is given by JB(R)= M X i=1 Ni X j=1 DKL[N(ˆµs ij,ˆ Covs ij)||N(¯µs j,¯ Covs j)] (3.35) where M is the number of groups, Ni is the number of epochs in group i , N(ˆµs ij,ˆ Covs ij) is the distribution of epoch j in group i and N(¯µs j,¯ Covs j) is the average distribution. The same group further extended this approach by including a discriminative term in the groupSSA objective function (Samek, M ¨ uller, Kawanabe and Vidaurre, 2012). This method not only considers the stationary components but also considers the discriminative information. The discriminative term is JB(R)=DKL[N(¯µs 1,¯ Covs 1)||N(¯µs 2,¯ Covs 2)].(3.36) This term is subtracted from the objective function given in Eqn. (3.35). After obtaining the stationary and discriminative components, the standard CSP is computed similarly with (Von B¨ unau et al., 2010). Another approach of regularizing the CSP algorithm is by including a penalty term in the standard CSP objective function. Unlike the method present in (Grosse-Wentrup et al., 2009), the authors of (Lotte and Guan, 2010b) proposed an approach to consider spatial information without any priori information. The penalty term P(w)is given by 27
P(w)=wT(DG)w,(3.37) where, Gij =exp(1 2||vivj||2 r2)and Dii =X j Gij vi and vj are the vectors containing the co-ordinates of the ith and jth electrodes, r is the closest distance between the two electrodes. The standard CSP computes the spatial filter based on the variance using l2 norm which makes the performance of CSP affected with outliers. As EEG signals consist of nonstationarities, the utilization of l2 norm leads to magnify the effect of noise which hinders the performance. Therefore to reduce the effect of outliers, the authors of (Wang et al., 2012) uses l1 norm ( ||.||1 ) for the computation of CSP filter which is referred as CSP-l1 ˜ J1(w)=||wTX||1 ||wTY||1 =PTx i=1|wTxi|1 PTy j=1|wTyj|1 .(3.38) Later, a more generalized optimization function was proposed in (Park and Chung, 2013) using lpnorm (||.||p) instead of using l1or l2norm, which is given by ˜ Jp(w)=||wTX||p ||wTY||p =[PTx i=1|wTxi|p](1/p) [PTy j=1|wTyj|p](1/p).(3.39) Another method of attending the nonstationarities nature of EEG signal is to use a penalty term, computed using the additional measurement like EOG or EMG signal, in the denominator of the CSP objective function (Blankertz et al., 2007). The stationary CSP (sCSP) maximizes the difference between class variance and at the same time keeps it stable within trials. The sCSP was first proposed in (Wojcikiewicz et al., 2011b,a) and later extended in (Samek, Vidaurre, M ¨ uller and Kawanabe, 2012). The penalty term is computed using a chunk of trials P(k) sCSP =|(Covk cCovc)|+,(3.40) where |.|+ is the non-negative truncation operator that flips the negative eigenvalues. The sCSP method can only extract the stationary feature but it is unable to address the problem of over-fitting. Therefore, the second extension mainly aims at tackling the problem of over-fitting by modifying the objective function as: ˜ J(w)= wTCov1w wT(Cov1+Cov2)w+⌘1PsCSP (w)+⌘2PTRCSP(w),(3.41) where PTRCSP is the Tikhonov regularizer similar with that used in (Lotte and Guan, 28
2011). Since the sCSP involve flipping of the signs, for some cases it doesn’t lead to an optimal solution. Therefore, the authors of (Arvaneh et al., 2011a, 2013) proposed a well-defined optimization function min wi p/2 X i=1 wiCov2wT i+ p X i=p/2+1 wiCov1wT i,(3.42) subject to: wi(Cov1+Cov2)wT i=1 i=1,2,...,p wi(Cov1+Cov2)wT j=1 i, j =1,2,...,p i6=j This can be reformulated in the presence of the penalty term as min wi (1 ⌘) p/2 X i=1 wiCov2wT i+ p X i=p/2+1 wiCov1wT i+⌘P(w),(3.43) where, P(w)=1 2 2 X c=1 1 T T X i=1 DKL(N(0,wCovi cwT)||N(0,wCovcwT)) subject to: wi(Cov1+Cov2)wT i=1 i=1,2,...,p wi(Cov1+Cov2)wT j=1 i, j =1,2,...,p i6=j where the penalty term for each class is the KL divergence between the trial covariance and the average covariance matrices. The regularization using l1 norm for estimation of covariance matrix proposed in (Wang et al., 2012). It mainly aims at removing the outliers of large deviation, but noises of small deviation are also present in the EEG signals. The CSP- l1 approach fails to consider the effect of noise. Therefore, CSP- l1 had been extended in (Wang and Li, 2016) by modifying the objective function of the CSP-l1approach as J(w)= 1 TxPTx i=1 Pt l=1|wTxi l| 1 TyPTy j=1 Pt l=1|wTyj l|+ TzPTz k=1 Pt1 l=1|wTek l|,ek l=zk lzk l+1,(3.44) where zk l are the trials used for noise modelling and could be obtained using all the trials of both the classes. xl and yl are the EEG signals of the two class. T denotes the number of trials and tis the number of sample points during a trial segment. 29
Another approach for regularizing CSP is by finding sparse spatial filters. The authors of (Arvaneh et al., 2011a) obtained the sparse solution by modifying the normalization term to l1/l2norm min wi (1 ⌘) p/2 X i=1 wiCov2wT i+ p X i=p/2+1 wiCov1wT i!+⌘ p X i=1 ||wi||1 ||wi||2 ,(3.45) where, P(w)=1 2 2 X c=1 1 T T X i=1 D(N(0,wCovi cwT)||N(0,wCovcwT)) subject to: wi(Cov1+Cov2)wT i=1 i=1,2,...,p wi(Cov1+Cov2)wT j=1 i, j =1,2,...,p i6=j where the parameter ⌘ controls sparseness of the data and accuracy. The solution is obtained by optimizing ⌘ . For selecting the channels only two sparse spatial filters are obtained for both the classes. The methods presented in (Yong et al., 2008b) and (Arvaneh et al., 2010) are limited to computation of a single spatial filter whereas the method in (Farquhar et al., 2006) can obtain more than one spatial filter but the correlation between the spatial filters are not considered. The previous approaches mainly concentrate in the selection of channels, whereas the authors of (Arvaneh et al., 2011b) computed the sparse filter without filtering the channel and also obtained uncorrelated filters. The other approaches such as Common Spatio-Spectral Pattern (CSSP) (Lemm et al., 2005), Common Sparse Spectral-Spatial Pattern (CSSSP)(Dornhege et al., 2006), Filter Bank CSP (FBCSP) (Ang et al., 2008) that computes both spectral and spatial filters are also proposed in the literature, which is not presented in this thesis. 3.6 Conclusions In this chapter, the different spatial filtering techniques used for discrimination of motor imagery movements are described. CSP is considered to be the most effective technique for this application. Besides, CSP is easily affected by the presence of the outliers. Hence, several regularized and alternative techniques are proposed which are discussed in this chapter. Moreover, the limitations of existing techniques are also mentioned. 30
31
CHAPTER 4 Study of Other CSP Based Approaches In the previous chapter, the CSP, its variant and other spatial filtering approaches were explained in details. However, the performance of the existing MI based BCI system is still low to use in real time applications. Therefore, different CSP based approaches which include the simplification of CSP solution as a solution of the Rayleigh quotient, the formulation of a discriminative CSP objective function and ICA combined algorithm with CSP are studied and presented in this chapter. Section 4.1 simplifies the CSP objective function as Rayleigh quotient. Section 4.2 presents the discriminative CSP. The ICA corrections to CSP is presented in section 4.3. 4.1 Simplification of the CSP Objective Function The standard CSP computes the spatial filter by maximizing the ratio of variance between the two classes. As discussed before, the objective function of the CSP algorithm is represented as the Rayleigh quotient J(w)=wTCov1w wTCov2w.(4.1) This can also be written as (4.2) wTCov1w wTCov2w=wTCov 1 2 2Cov1 2 2Cov1Cov1 2 2Cov 1 2 2w wTCov 1 2 2Cov 1 2 2w =wTCov 1 2 2MCov 1 2 2w ||Cov 1 2 2w||2 =˜ wTM˜ w where, M=Cov1 2 2Cov1Cov1 2 2(4.3) ˜ w=Cov 1 2 2w ||Cov 1 2 2w|| (4.4) The maximization and minimization of Eqn. 4.2gives the CSP solution. 32
4.2 Discriminative CSP The CSP objective function can be reformulated as discriminative CSP which considers the difference between the covariance of the two classes in the numerator. The objective function of the discriminative CSP is given as min \max wT(Cov1Cov2)w wT(Cov1+Cov2)w.(4.5) The formulated objective function maximizes the difference of the class covariances from the sum of the class covariance which provides more robustness in the presence of outliers. This solution of the discriminative CSP is equivalent to the solution of the standard CSP which is shown below. Theorem 4.2.1. The solution of discriminative CSP is equivalent to the solution of standard CSP max wT(Cov1Cov2)w wT(Cov1+Cov2)w=max¯ wTCov1 2 T(Cov1Cov2)Cov1 2 T¯ w,(4.6) where ¯ wT=wTCov 1 2 T Proof.The equivalence of the solution of CSP and discriminant CSP is shown here (4.7) wT[Cov1(CovTCov1)]w wTCovTw=wT(2Cov1CovT)w wTCovTw =¯ wT(2Cov1 2 TCov1Cov1 2 TI)¯ w and CovT=(Cov1+Cov2). The solution of the maximization problem is the same and equal to the generalized eigenvectors of Cov1and Cov2. 4.3 ICA Corrections to CSP As mentioned in the previous chapter, CSP and ICA are commonly used spatial filtering techniques for discrimination of motor imagery movements. We have tried to incorporate ICA in CSP objective function. It is known that the computation of CSP is easily affected by the presence of the outliers in the data. Therefore, to obtain a more robust objective 33
function, the ICA techniques are used in combination with the CSP approaches. Let X be EEG signals which can be represented by the given equation X=AS.(4.8) The covariance of Xwith zero mean is given by Covx=E[XXT]=AE[SST]AT.(4.9) Applying whitening transformation Tbefore the computation of CSP filter TCovxTT=I,where T=Cov1 2 x(4.10) and the whitened signal is given by Z=TAS =UTS.(4.11) We can further rewrite the above equation as TTATCov 1 2 2TTCov1TCov 1 2 2AT =I.(4.12) The CSP and ICA objective function can be formulated as (4.13) J(w)=wT((1 ⌘)Cov1+⌘CovP)w wTCovTw =wT((1 ⌘)Cov1+⌘CovP)w wT(Cov1+Cov2)w, where, the penalty term CovP can be the reference covariance computed from the orthogonal subspace to the artifact signals and J(w)=(CSP if ⌘=0, ICA if ⌘=1,. (4.14) The solution to this problem is given by Y=WTX(4.15) where, W=[w1,w 2,...,w p]is the discriminative spatial filters. 4.4 Experimental dataset and study This discriminative CSP and ICA corrections to CSP approaches have been tested using the publicly available BCI competition IV dataset 2a (Brunner et al., 2008). The dataset 34
consists of left hand, right hand, foot and tongue motor imagery movements acquired using 22 electrodes from nine healthy subjects. We have considered only the left and right hand motor imagery movements for this study. The MI EEG signals are band-pass filtered between 730Hz and 10- fold cross validation was applied to obtain the training set and the testing set. In the first experiment, we have done a performance comparison between the discriminative CSP and the standard CSP. For the second approach, we have used SOBI algorithm to find the related independent components. In this approach, we have introduced a parameter ⌘ to select between the CSP and ICA algorithm. Therefore, when ⌘=0 then the algorithm is considered as a standard CSP and when ⌘=1 then its ICA. If ⌘=0.5then we considered the combination of ICA with CSP algorithm. 4.5 Results The performance comparison between the discriminative CSP and standard CSP is shown in Table 4.1. From the result, we can observed that both the algorithms give the same performance. The performance results for ICA corrections to CSP approaches for different value of ⌘ are shown in Table 4.2. From the observation, we conclude that the standard CSP algorithm (i.e. ⌘=0) gives the maximum performance. Table 4.1 Performance comparison between CSP and discriminative CSP Methods Average Accuracy(%) CSP 81 discriminative CSP 81 Table 4.2 Performance comparison of ICA corrections to CSP for different values of ⌘=[0,0.5,1] ⌘Average Accuracy(%) ⌘=0, i.e. CSP 81 ⌘=0.580 ⌘=1i.e. ICA 67 4.6 Conclusions This chapter presents different approaches based on CSP. Although the proposed algorithms failed to improve the accuracy, it was interesting to study these. We hope that these approaches can further be developed and modified for obtaining a robust method. Performance results of these studies will not be presented in this thesis. 35
36
find the weight vectors busing the objective function J(b)=E(bTx)43||w||4+P(||w||2).(5.24) The constraint ||w||=1 is taken into consideration. P is a penalty term and several forms of penalty term were proposed in (Hyv¨ arinen and Oja, 1996). 5.2 Related Work As mentioned earlier, the popular algorithms to attend this issue are the CSP algorithm and various ICA techniques. CSP algorithm is a supervised process whereas ICA is an unsupervised technique where the class labels are unknown. CSP algorithm (Blankertz et al., 2008), which is considered to be the most effective method for discrimination of two class MI movements, is further extended to discriminate four MI movements in (Dornhege et al., 2004). Multiclass CSP considers the multiclass problem as a combination of a binary problem. ICA is commonly used for removing the artifacts from the EEG signals (Jung et al., 2000) and also used for discriminating the mu rhythms generated from the sensory motor cortex (Makeig et al., 2004). The authors of (Naeem et al., 2006) and (Brunner et al., 2007) evaluated different ICA algorithms like Infomax, FastICA, SOBI for separation of four-class MI movements. Both the studies conclude that the Infomax performs better than FastICA and SOBI but the performance of Infomax is comparable with the multiclass CSP. Another author proposed a new feature extraction method based on the Infomax algorithm using patterns from the independent components (Zhou et al., 2014). There are various other different approaches proposed for the classification of multiclass EEG signals. Similar with the multiclass CSP, the multiclass FBCSP (Chin et al., 2009) is an extension of the FBCSP method that uses different frequency bands for discrimination of MI movements. Several other approaches have been proposed for the classification of multiclass MI signals. The other approach for multiclass CSP extension is done by using simultaneous JAD techniques (Dornhege et al., 2004), which were further extended by incorporating information theory for extracting the features (Grosse-Wentrup and Buss, 2008). This approach showed the relation between JAD and BSS. It was observed that the JAD based multiclass approach outperforms the one versus rest CSP approach. The connections between JAD and BSS was demonstrated by forming a new formulation based on maximum likelihood framework (Gouy-Pailler et al., 2010). More recently, a group (Barachant et al., 2012) has proposed the separation of multiclass MI signals by exploiting the Riemannian geometry of the manifold of covariance matrices and using the tangent space for the classification of the features. Other authors proposed a new criterion for multiclass problem based on maximizing 43
the harmonic mean of the KL divergence between the class covariance matrices (Wang, 2012). The authors of (Nguyen et al., 2012) applied stationary CSP approach together with JAD to solve the multiclass problem. Other approach includes utilizing Bayesian learning method for multitask classification (Zhang et al., 2013). The authors of (Xu et al., 2011; Llera et al., 2014) and (Nicolas-Alonso et al., 2015) used the adaptive processing in discriminating multiclass problem. Another different approaches in solving the multiclass classification is by performing Fourier transformation (Townsend et al., 2006; Ge et al., 2014). This chapter provides the overview of the BSS techniques. A new criterion for the classification of four class MI movements was proposed based on BSS. 5.2.1 Multiclass CSP The standard CSP method derives the spatial filter by maximizing the variance of one class and simultaneously minimizing the variance of the other class. It is obtained by solving the Rayleigh quotient using the generalized eigenvalue problem J(w)=wTCov1w wTCov2w,(5.25) where Cov1 and Cov2 denotes the covariance matrices of class 1 and class 2. The resulting eigenvalues are sorted in a descending order and the largest and smallest p no. of features are selected. The corresponding eigenvectors of the selected eigenvalues represent the spatial filters. This method can discriminate only two motor imagery classes. Therefore, the same concept had been used to extend for discrimination of multiclass movements (Dornhege et al., 2004). In multiclass CSP, the multiclass problem is considered as a binary problem and spatial filters are obtained between one class upon the rest of the class. The eigenvalues were sorted in descending order and the p number of smallest eigenvalues were selected for each class. 5.3 Implementation of the Discrimination of the MI-EEG Signals The experimental design of the MI based BCI system consists of the following stages: (i) First, EEG signals acquired during the execution of mental tasks were used from BCI competition dataset. (ii) These signals were filtered using a bandpass filter (between particular cut off frequencies) and a small time segment of EEG signals was selected. (iii)The filtered EEG signals were used to determine the spatial filters of the corresponding MI tasks. (iv) The EEG signals were filtered using the obtained spatial filters. (v) The logarithmic variances of the spatially filtered signals were used as features for training the classifiers e.g. LDA and SVM. The overall experimental design of this study is describe briefly in the following sections. 44
5.3.1 Experimental Dataset In this study, dataset 2a from BCI competition IV (Brunner et al., 2008) was used. The data was acquired from nine subjects while performing the MI movement of four different task i.e. left hand, right hand, feet and tongue. The dataset consists of two sessions and each session consists of 288 trials. The EEG signals were acquired using 22 electrodes. The acquired signals were sampled at a sampling frequency of 250 Hz and filtered with a pass band filter of cut off frequencies 0.5 Hz and 100 Hz. A notch filter of 50 Hz was used to remove the line noise. 5.3.2 Preprocessing The mu rhythm represents the ERS and ERD of the sensory motor cortex during the MI movements. In order to select the appropriate information, the acquired EEG signals were filtered using a fifth order Butterworth filter with a cut off frequency between 8-30 Hz for each subject. From the filtered signal, the time segment of 0.5-2.5s after the cue has been selected for each trial. For later convenience, the selected signals were concatenated in a classwise order. 5.3.3 The Thin ICA-CSP Criterion and its Implementation An extension of the ThinICA algorithm, which was proposed in (Cruces, Cichocki and De Lathauwer, 2004) is been proposed in this section. The second and higher order statistics of the observations are combined to obtain a contrast function that we adapt and specialize for EEG processing. The acquired EEG signals with n recording channels at time t are denoted as x(t)2 Rn. Their linear model in term of the latent sources is given by x(t)=As(t)+n(t)(5.26) where A2Rn⇥m is an unknown mixing matrix, s(t)2Rm is the vector of independent sources with covariance matrix equal to the identity Covs=I , while n(t)2Rm denotes additive Gaussian noise of zero mean and covariance matrix Covn=2 nI . The prewhitening of the observations is a useful preprocessing that simplifies the ICA problem. The covariance of the observations is Covx=AAT+2 nI.(5.27) Therefore, pre-multiplying the observations with the whitening matrix T=(Covx2 nI)1 2 45
retains and whitens the principal signal subspace of the observed signals. We denote the prewhitened observations with z(t)=Tx(t)(5.28) =Us(t)+Tn(t)2Rn,(5.29) where the orthogonal matrix U=TA denotes the residual mixing matrix after prewhitenning. The simultaneous estimation of pn source signals (plus a non-separable noise component) can be obtained by multiplying the whitened observations by the transpose of the semi-orthogonal matrix ˆ U=[ ˆ u1,...,ˆ up]2Rn⇥p . Thus, the outputs or estimated sources can be represented as y(t)= ˆ UTz(t).(5.30) Classical ICA algorithms, like (Cruces, Cichocki and Amari, 2004), recover p independent components by maximizing the contrast function based on the sum of the squared kurtosis of the outputs max p X i=1 |Cum(yi(t),yi(t),yi(t),yi(t))|2,(5.31) over the Stiefel manifold, i.e., subject to the constraint. The function is given by: ˆ UTˆ U=I . Note that, for later convenience, the kurtosis of yi(t) is denoted with the fourth-order cumulant Cum(yi(t),yi(t),yi(t),yi(t)) . When the sources are nonstationary, it is convenient to exploit this diversity by splitting the data into K blocks with center at tk,k =1,...,K . The marginal contrast functions can be evaluated at each split of the data and we can simultaneously maximize the accumulated sum of the marginal contrasts. max K X k=1 p X i=1 |Cum(yi(tk),yi(tk),yi(tk),yi(tk))|2 subject to ˆ UTˆ U=I.(5.32) However, it is well known that, generally, the variance of the estimates increases with the order of the statistics. Therefore, for short data records, the variance is high and the previous approach is not very precise. A more reliable method consists in allowing the combination of several low-order statistics of the outputs to build the contrast that estimates the demixing matrix. The third order cumulants are usually not considered because they vanish for symmetric distributions, so it seems reasonable to consider only the combination of second and fourth-order statistics of the outputs. Additionally, when 46
the sources are correlated this additional diversity can be also exploited by combining the second-order cross-cumulants at different delays ⌧2T . All these considerations lead to the proposal of the maximization of the contrast function ⇥ (given below) in the Stiefel manifold (5.33) ⇥(ˆ U)= 4 K X k=1 p X i=1 |Cum(yi(tk),...,y i(tk))|2 +2 K X k=1 p X i=1 X ⌧2T|Cum(yi(tk+⌧),y i(tk))|2. Here, 2 and 4 are weighting proportional to the precision of the respective estimates of the second and higher order statistics. However, the drawback of this contrast function is that it is highly non-linear in ˆ U, and therefore, it is difficult to optimize. Fortunately, the previous contrast function can be generalized by decoupling the extraction candidates that appear at the arguments of the cumulants to obtain a new contrast function which is quadratic in each decoupled extraction matrix. The idea is to use four different estimates of the sources y[q] i(tk)=(ˆ U[q])Tz(tk),q=1,...,4 , to build following the ThinICA-CSP contrast function : (ˆ U[1],..., ˆ U[4])= (5.34) 4 K X k=1 p X i=1 X ⌧2TCum(y[1] i(tk),...,y[4] i(tk)) 2 +2 3 K X k=1 p X i=1 X ⌧2TCum⇣y[1] i(tk+⌧),y[2] i(tk)⌘ 2 +2 3 K X k=1 p X i=1 X ⌧2TCum⇣y[1] i(tk+⌧),y[3] i(tk)⌘ 2 +2 3 K X k=1 p X i=1 X ⌧2TCum⇣y[2] i(tk+⌧),y[3] i(tk)⌘ 2 which is quadratic with respect to each ˆ U[q] and should be independently maximized for each extraction candidate q=1,...,4 under the constraint (ˆ U[q])Tˆ U[q]=I . At its maximum value, all the estimates y[1] i(t)=... =y[4] i(t) will agree and recover p independent sources. We omit here the details of the proof of the contrast function nature of Eqn. 5.34 but the required steps are similar to those presented in (Cruces and Cichocki, 2003). The ThinICA-CSP contrast function in Eqn. 5.34 is sequentially maximized with respect to each extraction candidate ˆ U[q],q=1,...,4 . For this purpose its gradient rˆ U[q] is evaluated and later, with the help of its singular value decomposition, this 47
gradient is orthogonally projected onto the Stiefel manifold to obtain the new value of the extraction matrix ˆ U[q]=VLVT Rwhere [VL,⇤,VR]=svd(rˆ U[q]).(5.35) Following a similar approach to the one presented in (Cruces, Cichocki and De Lathauwer, 2004), it can be shown at each of these iterations we guarantee the monotonous ascent in the contrast function. Although this ICA method can easily recover a subset of the independent components from the observed signals, in our case, it is critical to select only those components which are related to the MI movements. We can easily favor their selection by initializing the unmixing matrix with the solution provided by the multiclass CSP algorithm. At the convergence of Eqn. 5.35 the projection matrix represents the spatial filters for four different MI movements. The parameters of the ThinICA-CSP method can be tuned depending on the dataset, in our simulations, we have set them to: p=8 (eight independent components had been extracted), K=3 (three splits in the data has been considered), T={1,...,20} (the set of delays), while the weighting terms have been set to 4=0.025 and 2=14 . The main steps of the Sub-ABLD iteration are summarized in Algorithm 1. 48
Algorithm 1 ThinICA-CSP algorithm 1: function ThinICA-CSP{X(t)},WmCSP ,p,K,T 2: Compute the whitening transform matrix T,T=(Covx2 nI)1 2X. 3: Whiten the input data, Z=TX. 4: Initialize the semi-orthogonal matrix ˆ Uby computing the SVD of WmCSP T1 5: Compute the estimated output, Yusing Eqn. 5.30 6: Initialize the iteration counter: i=0 7: repeat 8: Compute the second order statistic 9: Compute the fourth order statistic 10: Compute the objective function using Eqn. 5.34 11: Compute the gradient, rˆ U[q] 12: The projection is done by thin SVD factorization of the gradient rˆ U[q] 13: The matrix ˆ U[q]is obtained by using Eqn. 5.35 14: Compute the estimated output 15: until convergence. 16: return WT=ˆ UTT 17: end function 5.3.4 Feature Extraction In this study, 2 filters have been selected for each class which gives a total of p=8 spatial filters for four classes.The training and testing signals were filtered using the obtained spatial filters. The log transformation of the variance ( var(.) ) of the spatially filtered signals gives the required features. Featurei=log(var(ˆ uT iz)).(5.36) 5.3.5 Classification The extracted training features were used for training the classifiers. In order to perform a comparative analysis using the two classifiers, LDA and SVM were used for this study. In this study, the M-SVM 2 method from MSVM package (Lauer and Guermeur, 2011) is used for the discrimination of multiclass MI movements. The M-SVM 2 is the extended 49
version of 2-norm SVM for the multiclass approach. The detailed explanation can be found in (Guermeur and Monfrini, 2011). SVM is also applied in multiclass problem like LDA (Schl¨ ogl et al., 2005). 5.4 Results LDA and SVM are commonly use classifiers in the discrimination of MI movements for BCI applications. Therefore, in order to perform a comparative study of these two classifiers using BCI competition IV dataset 2a, we have tested on mCSP algorithm and the proposed ThinICA-CSP algorithm. The performance observed for all the combinations are shown in Table 5.1 and Fig. 5.3. From the results, it is observed that LDA performs better than SVM for both mCSP and ThinICA-CSP algorithm. Therefore, LDA classifier is chosen for all the remaining studies. Table 5.1 Comparative accuracy using LDA and SVM classifier with mCSP and ThinICA-CSP Methods Classifiers Accuracy(%) mCSP LDA 62 SVM 60 ThinICA-CSP LDA 64 SVM 63 Fig. 5.3 Comparative analysis of LDA and SVM classifier using mCSP and ThinICA- CSP 50
(a) (b) (c) (d) ThinICA_CSP ThinICA_CSP ThinICA_CSP ThinICA_CSP Fig. 5.4 Illustration of the four class (left hand, right hand, foot and tongue) MI movements processed with ThinICA-CSP algorithm, the corresponding filtered signals and spatial patterns for (a) left hand,(b) right hand, (c)foot and (d) tongue for subject A1. Moreover, from the figure, it can also observe that ThinICA-CSP gives the highest performance of 64% than the multiclass CSP. The processed EEG signals for subject A1 at each stage are shown in Fig. 5.4. The first row shows the input EEG signals for left hand, right hand, foot and tongue motor imagery signals which is processed using the proposed ThinICA-CSP algorithm. The third row represents the filtered signals and the fourth row shows the spatial patterns for four different MI movements. 5.5 Conclusions ThinICA-CSP algorithm for discrimination of multiclass MI movements is proposed. The comparative study using LDA and SVM classifiers shows that LDA performs better than SVM for this scenario. Therefore, further studies are done using LDA. The experimented result shows that the proposed algorithm performs better than multiclass CSP algorithm. From the above observation, it can be concluded that utilization of second and higher order statistics and initialization of unmixing matrix with multiclass CSP filter matrix improves the classification performance of multiclass movements. Hence, this contribution provides additional evidence in favor of the use of ICA techniques for the robust classification of multiclass movements. 51
52
6.4 The Divergence Optimization Interpretation of CSP Under the appropriate selection policy for the number of spatial filters for each class, the solution obtained by the CSP algorithm admits an interpretation in terms of the optimization divergence measures (here denoted by Div(·k·) ) between the Gaussian pdfs outputs for each class wi= arg min dim{W}=ni+1 max w2WDiv(p(yi|c1)kp(yi|c2)) ,(6.22) except for a probable permutation in the ordering of some of the spatial filters. The problem can be formulated using the following optimization problem wi= arg min dim{W}=ni+1 max w2WD(Cov(yi|c1)kCov(yi|c2)) (6.23) where D(·k·) refers to a divergence between the covariances of the conditional densities of the outputs. As a consequence of the assumption of zero mean Gaussian densities, the covariances are the only necessary statistics that summarize all the relevant information of the conditional data. In particular, the solution of the CSP algorithm was linked in (Fukunaga and KoonTz, 1970; Samek et al., 2014; Wang, 2012; Samek, Blythe, M ¨ uller and Kawanabe, 2013) with the optimization of the symmetric Kullback–Leibler divergence (sKL) DivsKL(p(yi|c1)kp(yi|c2)) = Zp(yi|c1)logp(yi|c1) p(yi|c2)dyi+Zp(yi|c2)logp(yi|c2) p(yi|c1)dyi, (6.24) =Z(p(yi|c1)p(yi|c2)) log p(yi|c1) p(yi|c2)dyi.(6.25) This divergence measures can be simplified to the symmetric Kullback–Leibler (sKL) divergence between the class conditional covariances DivsKL(p(yi|c1)kp(yi|c2)) = 1 2 Cov(yi|c1) Cov(yi|c2)+1 2 Cov(yi|c2) Cov(yi|c1)1(6.26) ⌘DsKL(Cov(yi|c1)kCov(yi|c2)) .(6.27) In this work, an extension of the existing KL to the criterion of the AB Log-Det divergence ( D(↵,) AB (.k.) ) between the class-conditional covariances defined as (Cichocki 59
et al., 2015) is proposed D(↵,) AB (Cov(yi|c1)kCov(yi|c2)) = 1 ↵ log ↵⇣Cov(yi|c1) Cov(yi|c2)⌘+⇣Cov(yi|c2) Cov(yi|c1)⌘↵ ↵++ (6.28) for ↵6=0,6=0,↵+6=0, where |x|+=(xx0, 0,x<0,(6.29) denotes the non-negative truncation operator. When the arguments covariances are scalars and ↵,>0 , the AB Log-Det divergence can also be rewritten as the logarithmic ratio between the weighted arithmetic mean of the scaled covariances ( Cov↵+(yi|c1) , Cov↵+(yi|c2)) and their weighted geometric mean, i.e., D(↵,) AB (Cov(yi|c1)kCov(yi|c2)) = 1 ↵ log ✓↵ ↵+Cov↵+(yi|c2)+ ↵+Cov↵+(yi|c1)◆ (Cov↵+(yi|c2)) ↵ ↵+(Cov↵+(yi|c1)) ↵+ . (6.30) Additionally if ↵+=1 , the AB Log-det divergence between covariances is proportional to the Alpha–Gamma divergence D(↵,) AB (.k.) (Cichocki, 2010) between the conditional densities D(↵,) AB (Cov(yi|c1)kCov(yi|c2)) ⌘2Div(,↵) AG (p(yi|c1)kp(yi|c2)) =2 ↵ log ✓Z⌦ p(yi|c1)dyi◆✓Zp(yi|c2)dyi◆↵ Zp(yi|c1)p↵(yi|c2)dyi (6.31) for ↵>0,>0,↵+=1. In Section 6.5.3, it is proven that, under certain conditions, the simple optimization of an AB Log-Det divergence also leads to the solution of the CSP algorithm. Although, the potential of these divergences does not rely on their plain optimization but instead rely on their optimization in the presence of some regularization terms that help to specify the desired solutions. Recently, several divergence criteria have been proposed for the extraction of the spatial dimensions with maximum discriminative power. Among these, the multiclass approach based on the maximization of the harmonic mean of Kullback–Leibler divergences (Wang, 2012) and the regularization framework based on the beta divergences 60
(Samek et al., 2014; Samek, Blythe, M ¨ uller and Kawanabe, 2013) are the most noteworthy methods. Another approach based on Bhattacharyya distance and Gamma divergence has also been proposed for classification of motor imagery movements (Brandl et al., 2015). Our work is motivated by the success of these methods in improving the classification accuracy and the robustness against the outliers. The distinctive property of the AB Log-Det divergence is that it smoothly connects (through its hyperparameters) a quite broad family of Log-Det divergences for SPD matrices, covering several relevant classical cases like: the KL divergence, the dual KL divergence, the Beta Log-Det family, the Alpha Log-Det family, the Power Log-Det family, as well as the Affine Invariant Riemannian divergence. 6.5 The Definition of the AB Log-Det Divergence Henceforth, we will work on the multidimensional observation vectors x=[x1,...,x n]T 2Rn . In order to simplify the notation, the covariance matrices of the two classes are renamed as follows P⌘Cov(x|c1),Q⌘Cov(x|c2).(6.32) The AB Log-Det divergence is a directed divergence that evaluates the dissimilarity between two multidimensional covariance matrices. It was defined in (Cichocki et al., 2015) as D(↵,) AB (PkQ)= 1 ↵ log ↵(Q1 2PQ1 2)+(Q1 2PQ1 2)↵ ↵++ (6.33) for ↵6=0,6=0,↵+6=0, while, for the singular cases, its definition is given by D(↵,) AB (PkQ)= 8 > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > : 1 ↵2htr ⇣(Q1 2P1Q1 2)↵I⌘↵log|Q1 2P1Q1 2|ifor ↵6=0,=0, 1 2htr ⇣(Q1 2PQ1 2)I⌘log|Q1 2PQ1 2|ifor ↵=0,6=0, 1 ↵2log (Q1 2PQ1 2)↵(I+ log(Q1 2PQ1 2)↵)+for ↵=, 1 2||log(Q1 2P1Q1 2)||2 Ffor ↵,=0. (6.34) The divergence depends only on the eigenvalues ⇤=diag(1,···, n) of the SPD matrix Q1/2PQ1/2 , which also coincide with the eigenvalues of the matrix Q1P , although 61
their eigenspaces differ. Given the eigenvalue decomposition Q1 2PQ1 2=V1⇤VT 1,(6.35) where V1 is the orthogonal matrix of eigenvectors, and ⇤= diag{1, 2,..., n} is the diagonal matrix with positive eigenvalues i>0,i=1,2,...,n . One of the properties of the AB Log-Det divergence is that it is invariant under a common change of basis on its matrix arguments, i.e., an invertible congruence transformation. Since, with the help of this specific transformation, we have P!(VT 1Q1 2)P(VT 1Q1 2)T=⇤,Q!(VT 1Q1 2)Q(VT 1Q1 2)T=I,(6.36) it can be inferred that the divergence is separable (over the generalized eigenvalues of the matrix pencil (P,Q) ) in a sum of marginal divergences that measure how far are each of the generalized eigenvalues from the unity, i.e., D(↵,) AB (PkQ)=D(↵,) AB (⇤kIn)= n X i=1 D(↵,) AB (ik1).(6.37) Hence, D(↵,) AB (PkQ)= 1 ↵ n X i=1 log ↵ i+↵ i ↵++ ,↵,,↵+6=0.(6.38) Similarly, for the singular cases, the divergence is D(↵,) AB (PkQ)= 8 > > > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > > > : 1 ↵2"n X i=1 ↵ ilog(↵ i)n#for ↵6=0,=0 1 2"n X i=1 ⇣ ilog( i)⌘n#for ↵=0,6=0 1 ↵2"n X i=1 log ↵ i 1+log↵ i+#for ↵=6=0 1 2 n X i=1 log2(i)for ↵,=0. (6.39) This divergence compares two symmetric positive definite matrices and returns its dissimilarity, i.e., a positive value when they are non-coincident and D(↵,) AB (PkQ)=0 iff P=Q . As it can be observed in Fig. 6.1 the AB Log-Det divergence generalizes 62
several existing log-det matrix divergences, like: the Steins loss, 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 (see Table 1 in (Cichocki et al., 2015) for a comprehensive list). Fig. 6.1 This illustration shows the AB Log-Det divergence D(↵,) AB (PkQ) positioned in a plane as a function of their real pair of hyperparameters (↵,) . It is clear from the figure, that the parameterization smoothly connects several relevant positive definite matrix divergences, like: 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. 6.5.1 A Tight Upper-Bound for the AB Log-Det Divergences The divergence D(↵,) AB (PkQ) depends on the generalized eigenvalues 1,···, n of the matrix pencil (P,Q) which, for convenience, are assumed to have a simple spectrum (the eigenvalues are unique or non-coincident) and can be sorted in descending order 1> 2>···> n>0.(6.40) In practice, the assumption is plausible because the real symmetric matrices with unique eigenvalues are known to form an open dense set in the space of all the real symmetric matrices (Tao, 2012). Although the space of the observations is high-dimensional, most of the discriminative information between the two conditions is confined into a low-dimensional subspace. Thus, the spatial filter matrix W2Rn⇥p is used to reduce the dimensionality of the samples from n to p with the linear compression transformation y=WTx2Rp . It is shown in (Cichocki et al., 2015) that, after applying this compression to the arguments 63
of the divergence, the resulting output covariance matrices WTPW and WTQW are more similar than in the original space, as shown in the below equation D(↵,) AB (WTPW kWTQW)= p X i=1 D(↵,) AB (µik1) n X i=1 D(↵,) AB (ik1) = D(↵,) AB (PkQ), (6.41) where µ1··· µp>0 are the generalized eigenvalues of the matrix pencil (WTPW,WTQW) . However, this upper bound is loose for the case of interest (dimensionality reduction), i.e., when p<n . In Appendix A, the possible way to tighten the previous upper-bound with the following new proposal is shown D(↵,) AB (WTPW kWTQW) p X i=1 D(↵,) AB (⇡ik1),(6.42) where ⇡ defines the permutation of the indices 1,...,n that sorts the divergence of the eigenvalues from the unity in descending order D(↵,) AB (⇡1k1) D(↵,) AB (⇡2k1) ···D(↵,) AB (⇡nk1).(6.43) Moreover, the equality with the upper-bound is only obtained for those extraction matrices W that lie within the span of the p generalized eigenvectors of the matrix pencil (P,Q) which are associated with the eigenvalues ⇡1,..., ⇡p that maximize the divergence from unity in Eqn. 6.43. 6.5.2 Relationship between the Generalized Eigenvalues and Eigenvectors of the Matrix Pencils (P,Q)and (p(c1)P,Cov(x)) We have seen in the previous section that the tight upper-bound of the divergence is attained by a subset of the generalized eigenvectors of the matrix pencil (P,Q) , whereas, the CSP solution in Eqn. 6.21 depends on a subset of the generalized eigenvectors of another matrix pencil (p(c1)P,Cov(x)) . In this section, the close relationship between both eigendecompositions is addressed. For this purpose, ⇤ is denoted as the matrix of eigenvalues of Q1P and ⇤(c1) as the matrix of eigenvalues of (Cov(x))1p(c1)P . Then, we write (p(c2)Q)1(p(c1)P)=[(Cov(x))1(p(c2)Q)]1[(Cov(x))1(p(c1)P)],(6.44) and use the decomposition of Cov(x) in Eqn. 6.3 to substitute p(c2)Q=Cov(x) p(c1)Pin the previous equation. In this way, we obtain (p(c2)Q)1(p(c1)P)=[In(Cov(x))1(p(c1)P)]1[(Cov(x))1(p(c1)P)].(6.45) 64
The matrix of eigenvectors V of Q1P diagonalizes both sides of the previous equation ⇤p(c1) p(c2)=V1[p(c1) p(c2)Q1P]V(6.46) =(V1[In(Cov(x))1(p(c1)P)]1V)(V1[(Cov(x))1(p(c1)P)]V)(6.47) =[InV1(Cov(x))1(p(c1)PV)]1(V1[(Cov(x))1(p(c1)P)]V)(6.48) =(In⇤(c1))1⇤(c1).(6.49) Hence, the explicit relationship between the two sets of eigenvalues is obtained ip(c1) p(c2)=(c1) i 1(c1) i⌘g((c1) i),i=1,...,n, (6.50) where g((c1) i) , as can be seen in Fig. 6.2, is a strictly monotonous ascending function over the domain of (c1) i2(0,1) . Moreover, the Equations (6.46)–(6.49) imply that the matrix V of generalized eigenvectors of the matrix pencil (P,Q) exactly coincides with the matrix V(c1) of generalized eigenvectors of the other matrix pencil (p(c1)P,Cov(x)) . 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10−2 10−1 100 101 102 λ( c1 ) g(λ( c1) ) Fig. 6.2 Illustration of the strictly monotonous ascending transformation g(·) that, through Eqn. 6.50, maps eigenvalues of the matrix pencil (p(c1)P,Cov(x)) into the eigenvalues of the matrix pencil (P,Q) , in a case where the sample probabilities of the classes are uniform p(c1)=p(c2) . Note that the eigenvalues of the first pencil are bounded in the interval (0,1) , while the domain of the eigenvalues of the second pencil is (0,1). 6.5.3 Linking the Optimization of the Divergence and the CSP Solution There is a link between the solutions of the CSP method and the solutions obtained with the optimization of the symmetric KL divergence between the class conditional 65
covariances, which was studied in previous works (Fukunaga and KoonTz, 1970; Samek et al., 2014; Wang, 2012). This subsection shows that under the appropriate filter selection criteria the link also extends to the optimization of other divergences, like the AB Log-Det family of divergences. We have previously assumed that generalized eigenvalues are ordered and can be regarded as non-equal. Therefore, they can be clustered in the following three sets of principal, inner and minor eigenvalues of the matrix pencil (P,Q): 1>···> k | {z } kprincipal eigenvalues > k+1 >···> n(pk) | {z } inner eigenvalues > n(pk)+1 >···> n | {z } (pk)minor eigenvalues .(6.51) The following sequence of optimizations induces an alternative ordering of the generalized eigenvalues D(↵,) AB (⇡ik1) = min dim{W}=ni+1 max w2WD(↵,) AB (wT iPwikwT iQwi),i=1,...,n, (6.52) according to a permutation ⇡ that sorts their marginal divergences from 1 in descending order D(↵,) AB (⇡1k1) ···D(↵,) AB (⇡pk1) D(↵,) AB (⇡p+1 k1) ···D(↵,) AB (⇡nk1). (6.53) For building the matrix of spatial filters WDiv ⌘[w1,w2.....wp] , one possible selection policy is to retain only the p most discriminative spatial filters for the considered divergence optimization problem, i.e., those that solve Eqn. 6.52 for i=1,...,p . The filters consist in p eigenvectors ( v⇡i with i=1,...,p ) of the matrix pencil (P,Q) that are arranged according to the permutation ⇡ . From the one-to-one relationship that exists between the generalized eigenvalues and eigenvectors of the matrix pencils (P,Q) and (p(c1)P,Cov(x)) (see the previous subsection) the solution takes the following form WDiv =[v⇡1,...,v⇡p]=[v(c1) ⇡1,...,v(c1) ⇡p].(6.54) This result tells us that the optimization of different divergences (in absence of other regularizing terms) only differs in the selection criteria for the spatial filters, which eventually determine the chosen subindices ⇡1,...,⇡ p. Now, the question of whether these spatial filters that solve the sequence of minimax divergence optimization problems min dim{W}=ni+1 max w2WD(↵,) AB (wT iPwikwT iQwi),i=1,...,p, (6.55) essentially coincide (up to a possible permutation in the order of the spatial filters) with 66
the spatial filters of the CSP solution in Eqn. 6.56 WCSP =[v(c1) 1,...,v(c1) k | {z } kprincipal eigenvectors ,v(c1) n(pk)+1,...,v(c1) n | {z } pkminor eigenvectors ],(6.56) has a simple answer. The straightforward comparison between Eqn. 6.54 and Eqn. 6.56 reveals that both solutions should essentially coincide when the subindices ⇡1,...,⇡ p are a permutation of the integers 1,...,k,n(pk)+1,...,n . Thus, the link between both techniques happens whenever CSP method adopts the filter selection policy of the divergence criterion in Eqn. 6.53. However, many of the CSP implementations find satisfactory to choose the number of spatial filters for each class a priori, respectively as k and pk (we will refer to this case as the original CSP filter selection policy), where k is close to p/2 in order to approximately balance the number of spatial filters for each class (Blankertz et al., 2007; Ramoser et al., 2000). In general, the use of a divergence based selection policy does not ensure a balanced representation of the spatial filters for each of the classes. For instance, consider the synthetic but illustrative situation for n=100 , where we wish to select p=8 spatial filters. If the generalized eigenvalues of the matrix pencil (P,Q) are shifted towards to zero, for instance, equal to {10,0.99,0.98,...,0.03,0.02,0.01} . In most cases, the solution WDiv will select as its columns: only k=1 principal eigenvectors and pk=7 minor eigenvectors, an unbalanced choice. In view of this potential limitation, an interesting question is whether it would be possible to modify the AB Log-Det divergence criterion so as to enforce that its solution essentially coincides with the one obtained by the CSP method with its original filter selection policy. We will show in the following that this requires only a suitable scaling 2R+ in one of the arguments of the divergence. Without loss of generality, we assume scaling in the second argument of the divergence. As it is shown in the Appendix B, there is a permutation ⇡0 of the indices of the spatial filters 1,...,p that links the CSP solution in Eqn. 6.21 with the optimization of the divergence w⇡0 i= arg min dim{W}=ni+1 max w2WD(↵,) AB (wT iPwikwT iQwi),i=1,...,p, (6.57) for any given 2(inf , sup)(6.58) with inf ⌘K(k+1 , n(pk)+1)sup ⌘K(k, n(pk))(6.59) 67
where the function K(a, b)= 8 > > > > > > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > > > > > > : ✓(ab)/ (a↵b↵)/(↵)◆1 ↵+ for ↵,,↵+6=0 ✓log(a/b)) (a↵b↵)/(↵)◆1 ↵ for ↵6=0,=0 ✓(ab)/ log(a/b)) ◆1 for ↵=0,6=0 exp ✓a↵log(eb↵)b↵log(ea↵) ↵(a↵b↵)◆for ↵=6=0 pab for ↵==0 (6.60) determines the value of the constant =K(a, b)2R that equalizes the value of the AB Log-Det divergences between any arbitrary a, b 2R constants (in the first argument) and (in the second argument), i.e., D(↵,) AB (ak)=D(↵,) AB (bk).(6.61) Note that the only role of the scaling factor is to adjust the reference value in one of the arguments of the divergence to ensure the exact balance in the number of spatial filters that are specialized in each class. As it is shown in the Appendix B, this scaling factor prevents that the minimax solution for i=1,...,p , could be attained by some eigenvectors associated with elements of the inner set of eigenvalues in Eqn. 6.51, so the chosen subset of eigenvectors have to essentially coincide with the principal and minor eigenvectors that form the CSP solution in Eqn. 6.56. In practice, a value of which is close to unity and meets the required bounds can be obtained from the truncated choice ?=8 > < > : inf +"for inf 1 1for 12(inf, sup) sup "for sup 1 (6.62) for an arbitrary small value of the constant ✏⌧sup inf . 68
@ @WD(↵,) AB (Cov(Y|c1)kCov(Y|c2)) = 2[PW QW(WTQW)1(WTPW)] (WTQW)1 2Z(WTQW)1 2(6.100) = +2(PW +QW)(WT(↵P+Q)W)1 2 ↵PW[(WT↵PW)+(WT↵PW) (WTQW)1(WT↵PW)]12 ↵QW[(WTQW) +(WTQW)(WT↵PW)1(WTQW)]1. (6.101) With the help of the particular form of the Woodbury identity for the matrix inverse [A+AB1A]1=A1(A+B)1(6.102) we simplify the terms within the brackets. Finally, the fact that ↵+=1 is used to confirm the proportionality with the gradient of the AG divergence given in Eqn. 6.95, @ @WD(↵,) AB (Cov(Y|c1)kCov(Y|c2)) = +2(PW +QW)(WT(↵P+Q)W)1 2 ↵PW[(WT↵PW)1(WT(↵P+Q)W)1] 2 ↵QW[(WTQW)1(WT(↵P+Q)W)1] (6.103) = +2(PW +QW)(WT(↵P+Q)W)1 2 PW[(WTPW)1↵(WT(↵P+Q)W)1] 2 ↵QW[(WTQW)1(WT(↵P+Q)W)1] (6.104) = +2((1 + ↵ )PW +(1+ ↵)QW)(WT(↵P+Q) W)12 PW(WTPW)12 ↵QW(WTQW)1 (6.105) = +2( 1 PW +1 ↵QW)(WT(↵P+Q)W)1 2 PW(WTPW)12 ↵QW(WTQW)1 (6.106) =2 @ @WDiv(↵,) AG (p(x|c2)kp(x|c1)) .(6.107) 75
6.7 Robustness of the AB Log-Det Divergence in Terms of ↵and The squared Riemann metric is known to be the natural distance in the manifold of SPD matrices, as it measures the squared length of the geodesic path between the arguments of the divergence (Harandi et al., 2017). However, in the real data there are usually several model contaminations (mismatches), including outliers or artifacts, that could make other robust divergences preferable. In this section, we study how the hyperparameters ↵ and can influence robustness of the AB Log-Det divergence with respect to the behavior of the squared Riemann metric, which is used as a reference. For convenience, the AB Log-Det divergence is denoted as a function of the spatial filter matrix Wby f(↵,)(W)⌘D(↵,) AB (WTPWkWTQW),(6.108) and its gradient expression given by Eqn. 6.68 is considered. The spatial filters that maximize this divergence should satisfy the following estimating equations @f(↵,)(W) @W= p X i=1 @µi @W (↵,)(µi)=0,(6.109) where µi,i=1,...,p , are the eigenvalues of matrix M , which was defined in Eqn. 6.65, and (↵,)(µi)=@f(↵,)(W) @µi ,i=1,...,p, (6.110) may be regarded as influence functions for each pair (↵,) that account for the penalty variation in the divergence with respect to µi . The complementary term to (↵,)(µi) in Eqn. 6.109, i.e., @µi @W , is a matrix of partial derivatives of the generalized eigenvalues µi with respect to the elements of the spatial filters W and, therefore, it is independent of the considered divergence. It is easy to observe that, in the particular case of ↵==0 , the expression in Eqn. 6.109 represents the estimating equation for the squared Riemann metric @f(0,0)(W) @W= p X i=1 @µi @W (0,0)(µi)=0.(6.111) In order to study the relative robusness to outliers, one can rewrite the estimating equation for a chosen pair of hyperparameters (↵,) in terms of the influence function for the squared Riemannian metric as @f(↵,)(W) @W= p X i=1 ✓@µi @W (0,0)(µi)◆w(↵,)(µi)=0,(6.112) 76
10−210−1100101102 0 2 4 6 8 10 µ D ( 0 , 0 ) ( µ || 1 ) 10−210−1100101102 10−1 100 101 µ w( 0 , 0 ) ( µ ) (a) 10−210−1100101102 0 2 4 6 8 10 µ D ( 1 , 1 ) ( µ || 1 ) 10−210−1100101102 10−1 100 101 µ w( 1 , 1 ) ( µ ) (b) Fig. 6.3 Illustration of the behavior of the AB Log-Det divergence D(↵,) AB (µ, 1) , and of its associated weight function w↵,(µ) , versus µ for different values of ↵= . Note that µ is shown in log-scale. (a) Squared Riemannian metric for ↵==0 (upper plot) and its weight function (lower plot); (b) Power Log-Det divergence for ↵==1 (upper plot) and its weight function (lower plot). 77
10−210−1100101102 0 2 4 6 8 10 µ D ( 1 , 0 ) ( µ || 1 ) 10−210−1100101102 10−1 100 101 µ w( 1 , 0 ) ( µ ) (a) 10−210−1100101102 0 2 4 6 8 10 µ D ( 0 , 1 ) ( µ || 1 ) 10−210−1100101102 10−1 100 101 µ w( 0 , 1 ) ( µ ) (b) Fig. 6.4 Illustration of the behavior of the AB Log-Det divergence D(↵,) AB (µ, 1) , and of its associated weight function w↵,(µ) , versus µ for different values of ↵6= . Note that µ is shown in log-scale. (a) Kullback–Leibler (KL) positive definite matrix divergence for ↵=1,=0, and its weight function (lower plot); (b) Dual KL positive definite matrix div. for ↵=0,=1,and its weight function (lower plot). where the scalar term w(↵,)(µ)= (↵,)(µ) (0,0)(µ)(6.113) 78
acts as a weight function that controls, for a given pair (↵,) , the magnitude of the effect in the estimation equation of departures of µifrom unity. The presence of outliers in the real data, typically results in eigenvalues µi that are too far from unity. However, depending on the problem, the higher prevalence of outliers may be stronger only for the greatest eigenvalues, or for the smallest eigenvalues, or simultaneously for the greatest and smaller eigenvalues. Those hyperparameters (↵,) that are able to down-weight the contribution of the outliers, are considered more robust. Therefore, the shape of the weight functions w(↵,)(µi) is useful to study the relative immunity of the AB Log-Det divergence to outliers. Fig. 6.3a shows the squared Riemannian metric ( ↵==0 ) and its weight function, which is flat since this divergence is taken as reference. Fig. 6.3b presents a similar plot for the Power Log-det divergence with ↵==1 . In this case, the bell shape of the weight function is an indicator of the robustness with respect to the presence of outliers in the greatest and smallest eigenvalues, since they will be down-weighted in the estimating Eqn. 6.112. Similar plots can be done by increasing the magnitude of ↵= , which progressively enhances the robustness. When ↵6= the divergence is asymmetric. Fig. 6.4 respectively present the Kullback–Leibler divergence for SPD matrices ( ↵=1,=0 ) and its dual version ( ↵=0,=1 ), together with their associated weight functions. These plots illustrate the asymmetric cases in situations where ↵+>0 and reveal that, when ↵ , the AB Log-Det divergences tend to be more robust against outliers in the large eigenvalues while, for ↵⌧ , the robustness tends to be with respect to the outliers in small eigenvalues. 6.8 Conclusions The key properties of AB Log-Det divergences have been summarized. We have reexamined the relation between the Common Spatial Pattern criterion with a predefined number of spatial filters for each class and its interpretation as an AB Log-Det divergence optimization problem, to show that a scaling factor in one of the arguments is necessary for the equivalence of the solutions. 79
80
CHAPTER 7 Optimization of Alpha-Beta Log Det Divergence Algorithm It is mentioned in Chapter 6 that the optimization of AB Log-Det divergence is nontrivial. This motivates us to use these divergences with the illustrative application of dimensionality reduction in BCI and explain how to perform their optimization. The EEG data has a typical high-dimensionality, a low signal to noise ratio and may have artifacts/outliers. The dimensionality reduction is then a necessary processing of the EEG signals for extracting those subspaces where the features have highest discriminative power. In this chapter, the Sub-ABLD algorithm is proposed based on the theoretical background presented in Chapter 6 for the discrimination of two class motor imagery movements. The chapter starts with the related reviews of the spatial filtering for MI movements. The proposed criterion and algorithm are presented in Section 7.2 and Section 7.3 describes the experimental study. The results obtained are presented in Section 7.4. 7.1 Review of Some Related Techniques for the Spatial Filtering of Motor Imagery Movements It is well known that the performance of CSP is easily affected by the presence of artifacts. To overcome this drawback several CSP variants algorithm have been proposed which has been discussed in Chapter 3. In this section, the related regularized variants of CSP that have been proposed to improve the classification performance is reviewed. The regularization approaches of CSP are mainly done either in the estimation of the covariance matrices or by modifying the CSP objective function. Most of them combine the estimation of the covariance matrices for each class with the regularization of the CSP objective function using penalty terms. Some of the approaches include the previous information (Lotte and Guan, 2010b), other subject data (Kang et al., 2009; Lotte and Guan, 2010a) and previous session data (Lu et al., 2009) for estimating the class covariance matrix. Another approach used M-estimators to compute the robust class covariance matrices (Xinyi Yong and Birch, 2008) and yet another approach obtained the covariance matrices by finding the minimum squared error (Kawanabe and Vidaurre, 2009). The authors of (Samek, Binder and M ¨ uller, 2013) 81
applied Multiple Kernel Learning (MKL) to combine the information from different subjects. It has been shown in (Lotte and Guan, 2011) that the regularization of the objective function is more useful than regularizing the estimated covariance matrix. Several approaches have been proposed by regularizing the objective function. The authors of (Blankertz et al., 2007) have additionally incorporated the EOG signals for reducing the ocular artifacts. Other authors have tried to ensure robustness by selecting only the important channels and produce sparse spatial filters (Arvaneh et al., 2011b; Farquhar et al., 2006; Yong et al., 2008b). Another approach is to robustify the system by obtaining only the stationary features. A robustify maximin CSP method was proposed that used a set of covariance matrices instead of an individual covariance matrix without using any other user data or data from the previous sessions (Kawanabe et al., 2009, 2014). In order to avoid the presence of the outlier, the CSP objective function has been formulated using lp -norm in (Wang et al., 2012; Park and Chung, 2013). The Stationary Subspace Analysis (SSA) algorithm was proposed to obtain the stationary subspaces of the time series EEG signals by considering only the stationary components of the signals. The limitation of this method is the detection of dissimilarity of the different class as a nonstationary feature (Von B ¨ unau et al., 2009). The group wise SSA (gwSSA) algorithm aims at obtaining the non-stationarities by dividing the dataset into different groups and calculating the minimum KL divergence between estimated source distribution of each trial in a group and the average distribution of the corresponding group. This algorithm not only allows the combining of the multisubject data but also the multiclass data (Samek et al., 2011). But, the gwSSA algorithm cannot find the discriminative information between the classes. The same group proposed a new approach for extracting the discriminative information, by subtracting the inter class divergences from the gwSSA objective function (Samek, M ¨ uller, Kawanabe and Vidaurre, 2012). To overcome the limitation of the SSA algorithm, two-step approaches have been proposed where the initial extraction of the stationary sources was done using the SSA method and later, the CSP was used for the computation of the spatial filters(Von B ¨ unau et al., 2010). Another approach to extract the stationary features is to reduce the nonstationarities between the two sessions. The supervised and unsupervised methods for adaptation of the data space have been proposed using KL divergence between the intersession data (Arvaneh et al., 2013). Recently, the authors of (Wu et al., 2015) presented MAP-CSP algorithm by deriving the probabilistic model of CSP to resolve the issue of overfitting of the baseline CSP algorithm. One of the limitations of the CSP algorithm is that it is mainly suitable only for the discrimination of two classes, while, in general, for an efficient BCI system more than two motor imagery movements are required. In order to formulate it for the multiclass system, the authors of (M ¨ uller-Gerking et al., 1999; Allwein et al., 2001) have 82
reduced the multiclass problem to a binary problem. The authors of (Dornhege et al., 2004) proposed two approaches for the multiclass problem; firstly to find the spatial filters for one class with respect to all the other classes and secondly, by simultaneous diagonalization methods. Other approaches, like (Grosse-Wentrup and Buss, 2008), proposed to solve the multiclass problems by combining information theoretic criteria with joint diagonalization methods. Several other methods have been proposed for the multiclass paradigm using independent component analysis (Naeem et al., 2006) and Riemannian geometry to obtain the spatial filters (Barachant et al., 2012). The authors of (Zhang et al., 2013) derived a relation between Bayes classification error and Rayleigh quotient and used this approach to solve the multiclass problem. In spite of all these different approaches, the performance of MI-based BCI systems is degraded due to the presence of non-stationarities and outliers, which is a challenge for the BCI systems in a real application. Hence, a robust feature extraction algorithm is needed to increase the overall performance of the system. 7.2 Proposed Criterion and Algorithm for Spatial Filtering For the presentation of the proposed criterion some additional notation needs to be defined. Let ˜ x(j)(t)|c denote the output of the passband filtering of the raw observations at time t and for the j th trial of class c2{c1,c 2} . The power of the trials of a given class cis normalized by the operation x(j)(t)= ˜ x(j)(t) qtr{Cov(˜ x(j)|c)} ,(7.1) where Cov(x(j)|c)= 1 L L X t=1 x(j)(t)¯ x(j)x(j)(t)¯ x(j)Twith ¯ x(j)=1 L L X t=1 x(j)(t) (7.2) denotes the sample covariance matrix the jth trial x(j) of class c , and L is the size in samples of each trial. In order to simplify the notation, the covariance matrices of the two classes are renamed as Pj⌘Cov(x(j)|c1)and Qj⌘Cov(x(j)|c2),(7.3) and their averaged versions (the centroids of each class) are denoted as P⌘<Pj>=1 N1 N1 X j=1 Pjand Q⌘<Qj>=1 N2 N2 X j=1 Qj.(7.4) 83
The classification of imagery movements involves extracting the relevant features of the observations and the classification of the observed patterns in the feature space. In the considered application, the data is high-dimensional but only a few features are sufficient to capture the discriminative information about the intended movements. Thus, the extraction of the relevant features involves a dimensionality reduction step for the observations from Rn to Rp where p⌧n . This step is implemented through the spatial filtering, i.e., by projecting the n -dimensional observations onto a p -dimensional subspace which should allow a good discrimination of the cluster centroids and, at the same time, guarantee a compact representation of the clusters. As mentioned earlier, the CSP solution will be obtained by a minimax optimization of the divergence between the projected and scaled centroids of the classes, i.e., D(↵,) AB (wT iPwikwT iQwi) . However, since this solution completely ignores the within-class dispersion of the samples, it is quite sensitive to artifact and outlier in the training dataset. In similarity with the divergence framework presented in (Samek et al., 2014) and with some variants of Fisher LDA (Barber, 2012, pag. 366), one can regularize the previous problem by controlling the dispersion of the trials of each class around their centroids and also by exploiting the degrees of freedom in the selection of the hyperparameters of the divergences. Then, a robust criterion based on the AB Log-Det divergence takes the following form F(W)=D(↵,) AB (WTPWkWTQW)⌘(p(c1)R1+p(c2)R2),(7.5) where the penalties associated to the within-class dispersion involve the averaged divergences R1=1 N1 N1 X j=1 D(↵,) AB (WTPjWkWTPW),(7.6) R2=1 N2 N2 X j=1 D(↵,) AB (WTQjWkWTQW),(7.7) and the parameter ⌘2R+ controls the balance between the maximization of the betweenclass scatter and the minimization of the within-class scatter. Note that in Eqn. 7.7 the fact that the AB Log-Det divergence is invariant under the common scaling of its arguments is used, to simplify D(↵,) AB (WTQjWkWTQW)=D(↵,) AB (WTQjWkWTQW). The optimization of the criterion in Eqn. 7.5 can be performed simultaneously, for all the spatial filters, with the use of subspace techniques (Samek et al., 2014). In the next section, a subspace optimization algorithm based on AB Log-Det divergences is presented. 84
In order to carry out a fair performance comparison, a total of six features (i.e., p=6 ) have been selected for both the algorithms. The performance comparison between both the algorithms is presented in the following subsections. 7.4.1 Observations for Simulated Data To study the performance of the proposed algorithm in the presence of outliers, the experiment was done by increasing the percentage of outlier trials in the training set for both the classes. 0 5 10 15 20 25 30 Percentage of outliers trials 90 91 92 93 94 95 96 97 98 99 100 Accuracy (%) CSP Sub-ABLD Fig. 7.1 Performance comparison of the proposed algorithm Sub-ABLD ( ⌘=1,↵= =1.5) with CSP versus the percentage of outlier trial The performances of the above algorithms with the increasing percentage of outliers in the training set are presented in Fig. 7.1. It can be observed that CSP performs worse in the presence of the outliers than the proposed Sub-ABLD algorithm. 7.4.2 Observations for BCI Competition Datasets In this section, the proposed algorithm is tested using three BCI competition datasets. For each dataset, the performances of the proposed algorithm for the different values of (↵,) and ⌘ were observed. From the observation, the maximum performance of the Sub-ABLD algorithm for the particular (↵,) and ⌘ values was selected. The selected performance is compared with the performances of other existing algorithms. Further analysis is done by using a box plot comparison for all the algorithms. The box plot analysis shows the distribution of the performances. In a box plot representation, the line inside the box represents the median performance. The upper and lower hinge of the 91
CSP Sub−ABLD 50 55 60 65 70 75 80 85 90 95 100 84.5% 89% Methods Accuracy (%) (a) CSP Sub−ABLD 50 55 60 65 70 75 80 85 90 95 100 Accuracy (%) (b) Fig. 7.2 (a) Performance comparison of the proposed algorithm Sub-ABLD ( ⌘=2,↵= =1.5 ) with CSP using BCI competition III dataset 3a and (b) its corresponding boxplot. 92
CSP Sub−ABLD 50 55 60 65 70 75 80 85 90 95 100 88.1% 88.9% Methods Accuracy (%) (a) CSP Sub−ABLD 50 55 60 65 70 75 80 85 90 95 100 Accuracy (%) (b) Fig. 7.3 (a) Performance comparison of the proposed algorithm Sub-ABLD ( ⌘= 0.5,↵==2 ) with CSP using BCI competition datasets III dataset 4a and (b) its corresponding boxplot. 93
CSP Sub−ABLD 50 55 60 65 70 75 80 85 90 95 100 81% 80.8% Methods Accuracy (%) (a) CSP Sub−ABLD 50 55 60 65 70 75 80 85 90 95 100 Accuracy (%) (b) Fig. 7.4 (a) Performance comparison of the proposed algorithm Sub-ABLD ( ⌘= 0.25,↵==1.25 ) with CSP using BCI competition datasets IV dataset 2a and (b) its corresponding boxplot. 94
Fig. 7.5 Results of the Sub-ABLD algorithm for the subject k6 from BCI competition III dataset 3a. This figure illustrates the changes in the average classification performance with respect to the variation of the parameters ↵ and . Relatively good performance results are obtained close to the diagonal and for moderately large values of the parameters. box denote the 75th and 25th percentile of the overall performance distributions. The whiskers are symbolized by the two lines outside the box. The upper and lower whisker represents the maximum and minimum performance observed. For BCI competition III dataset 3a, the Fig. 7.2a shows the comparison of the highest average performance of the Sub-ABLD algorithm with the average performances of other existing algorithms. From the figure, it is observed that the Sub-ABLD algorithm outperforms the CSP algorithm with an average performance accuracy of 89% for this dataset. The box plot comparison is shown in Fig. 7.2b. Although the median performance is slightly higher for CSP, the 25th percentile performance is much smaller than the one of the Sub-ABLD algorithm. As it will be seen later, is a consequence that with the Sub-ABLD algorithm the most difficult subjects have attained a significant improvement in their classification performance. Fig. 7.3 shows the observed average performances using BCI competition III dataset 4a. For this dataset, the proposed algorithm Sub-ABLD performs slightly above than the average performance of CSP, which is 88.1% . From the box plot of the results, it can be observed that the 25th percentiles for both the algorithms are also quite close. Similar results have been obtained for the BCI competition IV dataset 2a, which is shown in Fig. 7.4. Again the algorithms performance of Sub-ABLD is same as CSP, which average performance is 81% . In the box plot, we can observe that the quartiles of both the algorithms coincide. 95
To analyze the effect of performance for different divergences, the parameters (↵,) is varied for a single subject (Subject k6 from BCI competition III dataset 3a, which is one of the subjects with a worst performance for the experiment) and obtained the corresponding performance. The values of (↵,) are varied to cover the interval [0,2] ⇥[0,2] with a mesh of 0.1 spacings. The observed performance is shown in Fig. 7.5. This figure reveals a tendency to improve the classification accuracy of the worst user for values of ↵ and that are close to the diagonal and large enough so they can effectively down-weight the contribution in the estimating equations coming from the largest and smallest generalized eigenvalues. 7.5 Conclusions The Sub-ABLD algorithm has been proposed by optimizing the proposed criterion based on AB Log-Det divergence for discrimination of two class imagery movements. The spatial filters are computed by considering both. This algorithm was tested with synthetic and real datasets and compared with the standard CSP algorithm. The simulations have confirmed the possibility to tune up the hyperparameters of the divergence so as to improve the robustness of the obtained solutions without deteriorating the expected accuracy. 96
97
CHAPTER 8 Simulations This thesis presents two algorithms: ThinICA-CSP algorithm which aims at discriminating four class MI movements and Sub-ABLD for discriminating two class MI movements. The comparison of the performances of the proposed algorithms with the baseline method was already presented in chapter 5 and chapter 7. In this chapter, the performance of the proposed algorithms is compared with existing algorithms. All the coding and execution are performed by using MATLAB 2013a version. Section 8.1 presents the simulation study and comparison of the performance obtained using ThinICA-CSP algorithm. An experimental study of Sub-ABLD is discussed in Section 8.2. 8.1 Simulations using ThinICA-CSP algorithm for discrimination of four class motor imagery movements The comparison of the performance of ThinICA-CSP with the other algorithms like mCSP and JADE for discrimination of four class movements are presented in this section. mCSP is an extension of the standard CSP by regarding the multiclass problem as binary problems for computation of spatial filters (Dornhege et al., 2004). JADE algorithm performs a joint approximate diagonalization of the trial covariance matrices of the classes (Grosse-Wentrup and Buss, 2008). 8.1.1 Experimental Set-up The dataset that was used in Chapter 5 is being considered here. The experimental study is also similar to that explained in section 5.3. For this study, the JADE algorithm is included for a comparative study of performance. The filtered signals are used for the computation of spatial filters using the multiclass CSP, JADE method and the proposed ThinICA-CSP algorithm. Similar to the process presented in chapter 5, two spatial filters have been selected for each class which gives a total of p=8 spatial filters for four classes. The training and testing signals were filtered using the obtained spatial filters. The log transformation of the variance of the spatially filtered signals gives the required features. The extracted training features were used for training LDA classifier. The overall experimental study is shown in Fig 8.1 98
EEG signals (BCI Competition IV dataset 2a) Training data Testing data Obtaining training features (2 from each class i.e. p=8) Obtaining testing features Classification (LDA) Training phase Testing phase Learning classifier Multiclass Spatial filtering techniques (mCSP, JADE, ThinICA-CSP) Band pass filter Concatenate data Fig. 8.1 Experimental study for performance comparison of mCSP, JADE and ThinICA- CSP 8.1.2 Performance Results The performance results of ThinICA-CSP based on the above experimental set up are presented here. Fig. 8.3 shows one spatial pattern of the selected features during MI movements of the left hand, right hand, foot and tongue for subject A1 . The first, second and third column in Fig. 8.3 represents, respectively the spatial patterns obtained using ThinICA-CSP, JADE and multiclass CSP. In the first column of the figure, it can be observed that the ERD activity occurs in the right motor cortex during the left hand MI movement. Similarly, the right hand MI movement results in the ERD activity in the left motor cortex, foot motor imagery shows both ERS and ERD activities in the mid-central region and tongue MI movement denotes more ERS activity on the parietal region. Thus, the above observations agree with the findings in (Pfurtscheller et al., 2006). But, the patterns obtained using JADE and multiclass CSP are not much in agreement with the above findings. This shows that the ThinICA-CSP selects more relevant features than the above mentioned methods. For further analysis, the performance accuracies of the above methods are compared using the same dataset. The performance results are shown in Fig. 8.4. From the figure, it is observed that the ThinICA-CSP gives the highest performance of 64% than the multiclass CSP and JADE algorithms. 99
Nasion Inion POz P1 P2 Fz CP2CP1 Pz CPz C2 C4 C6 CP4 CzC1C3C5 CP3 FC3 FC1 FCz FC2 FC4 Fig. 8.2 Legends of the electrode locations used for acquisition of EEG from BCI competition IV dataset 2a along with nasion and inion. Fig. 8.3 Spatial pattern obtained during MI movements of (a) left hand, (b) right hand, (c) foot and (d) tongue for subject A1 from BCI competition IV dataset 2a using ThinICA- CSP, JADE and multiclass CSP The box plot comparison of all the three methods is shown in Fig. 8.5. The median performance is represented by the middle line inside the box. The 75th percentile and 100
CSP JADE MAPCSP divCSP−WS Sub−ABLD 50 55 60 65 70 75 80 85 90 95 100 88.1% 89.5% 69.8% 88.6% 88.9% Methods Accuracy (%) (a) CSP JADE MAPCSP divCSP−WS Sub−ABLD 50 55 60 65 70 75 80 85 90 95 100 Methods Accuracy (%) (b) Fig. 8.9 (a) Performance comparison of the proposed algorithm Sub-ABLD ( ⌘= 0.5,↵==2 ) with CSP, JADE, MAPCSP and divCSP-WS ( =0.33, 0 ⇤=0.5 ) using BCI competition datasets III dataset 4a and (b) its corresponding boxplot. 107
CSP JADE MAPCSP divCSP−WS Sub−ABLD 50 55 60 65 70 75 80 85 90 95 100 81% 80.5% 75.3% 80.7% 80.8% Methods Accuracy (%) (a) CSP JADE MAPCSP divCSP−WS Sub−ABLD 50 55 60 65 70 75 80 85 90 95 100 Methods Accuracy (%) (b) Fig. 8.10 (a) Performance comparison of the proposed algorithm Sub-ABLD ( ⌘= 0.25,↵==1.25 ) with CSP, JADE, MAPCSP and divCSP-WS ( =0.2, 0 ⇤=0 ) using BCI competition datasets IV dataset 2a and (b) its corresponding boxplot. 8.2.2 To Study the Performance of the Proposed Algorithm in Different Scenarios The following experiments are performed on CSP and the proposed algorithm. (1) Computation of average using Arithmetic and Geometric mean: Arithmetic mean is commonly used for averaging the distributions. However, considering the covariance 108
manifold which is not exactly a plane surface, averaging using geometric mean seems to be more meaningful and accurate. (2) Including outliers in the training trials: Oultier plays a major role in the computation of the efficiency of a system. Therefore, to study the performance of the algorithms, artificial outliers with large variance are included in the EEG dataset. The outlier matrices are generated randomly by drawing a Gaussian matrice Aout with i.i.d. elements and obtained Covout =Aout(Aout)T. Later, the generated Covout is scaled by a factor of 1. (3) Imbalancing the number of training trials for two classes: Generally, equal number of training trials for both the classes are considered for the computation of spatial filters. Therefore, to observe the changes in the performance by imbalancing the number of the training trials, unequal number of training trials for class 2 is considered. (4) Varying the number of training trials: The number of trials used for training the system influences both the performance as well as the processing time. Hence, to select the minimum numbers of trials without compromising the accuracy and the processing time, the experiment is performed by varying the k-folds value. 8.2.2.1 Experimental Set-up These experiments are done using BCI competition III dataset 3a and competition IV dataset 2a.The preprocessing of the data is done similar to the steps described in Section 7.3.2. In the first experiment, both the Arithmetic and Geometric mean are used for computing the average covariance matrices and the performance obtained is compared. The second experiment is done by including the artificial outlier trials in BCI competition dataset and observing the changes in performance by increasing the number of outlier trials. The third experiment is done by varying the number of trials of the two classes used for the computation of spatial filters. The last experiment is done by varying the number of training trials. This is done by varying the k values from [k=5,. . . ,30]. 8.2.2.2 Performance Results All the performance results obtained using BCI competition III dataset 3a and BCI competition IV dataset 2a are presented. The performance results with different methods for averaging the class covariance matrices are shown in Fig. 8.11. From the figure, it is observed that for BCI competition III dataset 3a, the Arithmetic mean performs better than the Geometric mean, whereas both the methods perform almost similar for BCI competition IV dataset 2a. Fig. 8.12 represents the performance of the algorithms in the presence of the outliers. For this experiment, the performance is observed by increasing the number of outlier trials in both the class covariance matrices. The performance degrades with the increased in the number of outlier trials for both the algorithms. 109
Sub-ABLD CSP 0 10 20 30 40 50 60 70 80 90 100 Accuracy (%) Arithmetic Mean Geometric mean (a) Sub-ABLD CSP 0 10 20 30 40 50 60 70 80 90 100 Accuracy (%) Arithmetic Mean Geometric mean (b) Fig. 8.11 (a)Performance comparison of Sub-ABLD algorithm ( ⌘=2,↵==1.5 ) and CSP with Arithmetic mean and Geometric mean using BCI competition datasets III dataset 3a and (b) Performance comparison of Sub-ABLD algorithm ( ⌘=0.25,↵= =1.25 ) and CSP with Arithmetic mean and Geometric mean using BCI competition datasets IV dataset 2a. 110
10 15 20 25 30 35 40 Number of outliers trial 50 55 60 65 70 75 80 85 90 Accuracy (%) Sub-ABLD CSP (a) 10 20 30 40 50 60 70 80 90 No. of outliers trials 50 55 60 65 70 75 80 85 90 Accuracy (%) SubABLD CSP (b) Fig. 8.12 (a) Performance comparison of Sub-ABLD algorithm ( ⌘=2,↵==1.5 ) and CSP by increasing the number of outlier trials using BCI competition datasets III dataset 3a and (b)Performance comparison of Sub-ABLD algorithm ( ⌘=0.25,↵= =1.25 ) and CSP by increasing the number of outlier trials using BCI competition datasets IV dataset 2a. 111
10 15 20 25 30 35 40 45 50 Number of trials included for class 2 50 55 60 65 70 75 80 85 90 Accuracy (%) SubABLD CSP (a) 10 20 30 40 50 60 70 80 90 Number of trials included for class 2 50 55 60 65 70 75 80 85 90 Accuracy (%) SubABLD CSP (b) Fig. 8.13 (a) Performance comparison of Sub-ABLD algorithm ( ⌘=2,↵==1.5 ) and CSP by imbalancing the number of training trials for two class using BCI competition datasets III dataset 3a and (b) Performance comparison Sub-ABLD algorithm ( ⌘= 0.25,↵==1.25 ) and CSP by imbalancing the number of training trials for two class using BCI competition datasets IV dataset 2a. 112
5 10 15 20 25 30 Number of k-folds 50 55 60 65 70 75 80 85 90 Accuracy (%) Sub-ABLD CSP (a) 5 10 15 20 25 30 Number of k folds 40 45 50 55 60 65 70 75 80 85 90 Accuracy (%) SubABLD CSP (b) Fig. 8.14 (a) Performance comparison of Sub-ABLD algorithm ( ⌘=2,↵==1.5 ) and CSP by varying the number of training trials using BCI competition datasets III dataset 3a and (b) Performance comparison Sub-ABLD algorithm ( ⌘=0.25,↵== 1.25 ) and CSP by varying the number of training trials using BCI competition datasets IV dataset 2a. The next experiment is done by imbalancing the number of trials used in the training set. The results obtained during this experiment are shown in Fig. 8.13. It shows that the performance is poor if the number of trials used for computing the spatial filters from 113
the two classes is unequal. Hence, we observed the better performance by increasing the number of trials equally towards class 1 for both the algorithms. This indicates that the algorithms perform better by selecting the same number of trials for both the class. Finally, the performance results presented in Fig. 8.14 are obtained by varying the number of training trials. The observed figures show that the overall performance when k = 10 and k=30 is the same, but the time taken for simulation is more for k=30. Hence, 10 folds cv is appropriate for these two datasets. 8.3 Conclusions In this chapter, different studies have been performed to analyze the robustness of the proposed algorithms. The performance and pattern of ThinICA-CSP is compared with mCSP and JADE. It shows that ThinICA-CSP outperforms the other two algorithms. The Sub-ABLD algorithm was tested with both artificial and real data and the performance comparison with the other existing algorithms like CSP, JADE, MAP-CSP and divCSP- WS was also presented. Sub-ABLD outperforms the other algorithm in the artificial and BCI competition III dataset 3a whereas the performance is same with CSP and JADE for the other two datasets. Furthermore, we have also studied the performance of the algorithm in various scenarios. 114
115
CHAPTER 9 Conclusions This thesis presents two algorithms for the classification of motor imagery movements in BCI applications. At the start, the anatomy and physiology of the human brain, EEG signals, different brain rhythms and artifacts present in the EEG signals were introduced. The MI-based BCI system and steps involved in it such as filtering and classification techniques were also discussed. Later, various existing spatial filtering algorithms and the challenges are studied. Based on this study, two spatial filtering algorithms were proposed. The first algorithm is based on the extension of ThinICA and mCSP to address the problem for discrimination of four class motor imagery movements. In this approach, the contrast function is obtained by combining the second order and higher-order statistics. The existing ThinICA algorithm extracts the independent components by considering only the higher order statistic, but for the non-stationary sources, it is more convenient to analyze the data by splitting it into K blocks. Thus, the marginal contrast function was evaluated for each split and simultaneously maximize the accumulated sum of the marginal contrasts. Furthermore, the combination of the several lower-order statistics of the outputs with delays was also incorporated in the contrast function to estimate the demixing matrix. Although ICA can recover the subset of independent components, it is critical to select only the MI related components. Therefore, it can be done by initializing the unmixing matrix with the solution provided by mCSP. The comparative study using LDA and SVM classifiers was also performed. It is observed that LDA performs better than SVM. The overall performance result shows that the proposed algorithm performs better than mCSP and JADE. This indicates that the utilization of second and higher order statistics and initialization of mixing matrix with the mCSP solutions improves the classification performance. Hence, this provides an additional evidence in favor of ICA techniques for designing a robust classification algorithm. The second algorithm is based on AB Log-Det divergence, which is the main contribution of this thesis. The AB Log-Det divergence optimization problem has been interpreted in terms of CSP criterion. The parameter that provides equivalent solutions between the CSP and AB Log-Det divergence optimization was obtained. The gradient for AB Log-Det divergence was derived. The robustness of the criterion with respect to 116
Cruces, S., Cichocki, A. and De Lathauwer, L. (2004), Thin QR and SVD factorizations for simultaneous blind signal extraction, in ‘2004 12th European Signal Processing Conference,’, IEEE, pp. 217–220. Dornhege, G. (2007), Toward Brain-Computer Interfacing, MIT press. Dornhege, G., Blankertz, B., Curio, G. and M ¨ uller, K. R. (2004), ‘Boosting bit rates in noninvasive EEG single-trial classifications by feature combination and multiclass paradigms’, IEEE Transactions on Biomedical Engineering 51(6), 993–1002. Dornhege, G., Blankertz, B., Krauledat, M., Losch, F., Curio, G. and Muller, K. R. (2006), ‘Combined optimization of spatial and temporal filters for improving brain-computer interfacing’, IEEE Transactions on Biomedical Engineering 53(11), 2274–2281. Duda, R. O. and Hart, P. E. (1973), Pattern Classification and Scene Analysis, Vol. 3, Wiley New York. Edelman, A., Arias, T. A. and Smith, S. T. (1998), ‘The geometry of algorithms with orthogonality constraints’, SIAM Journal on Matrix Analysis and Applications 20 (2), 303– 353. Farquhar, J., Hill, N., Lal, T. N. and Sch ¨ olkopf, B. (2006), Regularised CSP for sensor selection in BCI, in ‘3rd International BCI workshop, Austria’, pp. 1–2. Feder, M. and Merhav, N. (1994), ‘Relations between entropy and error probability’, IEEE Transactions on Information Theory 40(1), 259–266. Fukunaga, K. and KoonTz, W. L. G. (1970), ‘Application of the Karhunen-Loeve expansion to feature selection and ordering’, IEEE Transactions on Computers C- 19(4), 440–447. Gargiulo, G., Calvo, R. A., Bifulco, P., Cesarelli, M., Jin, C., Mohamed, A. and van Schaik, A. (2010), ‘A new eeg recording system for passive dry electrodes’, Clinical Neurophysiology 121(5), 686–693. Ge, S., Wang, R. and Yu, D. (2014), ‘Classification of four-class motor imagery employing single-channel electroencephalography’, PloS one 9(6), e98019. Gouy-Pailler, C., Congedo, M., Brunner, C., Jutten, C. and Pfurtscheller, G. (2010), ‘Nonstationary brain source separation for multiclass motor imagery’, IEEE Transactions on Biomedical Engineering 57(2), 469–478. Griss, P., Tolvanen-Laakso, H. K., Merilainen, P. and Stemme, G. (2002), ‘Characterization of micromachined spiked biopotential electrodes’, IEEE Transactions on Biomedical Engineering 49(6), 597–604. 123
Grosse-Wentrup, M. and Buss, M. (2008), ‘Multiclass common spatial patterns and information theoretic feature extraction’, IEEE Transactions on Biomedical Engineering 55(8), 1991–2000. Grosse-Wentrup, M., Liefhold, C., Gramann, K. and Buss, M. (2009), ‘Beamforming in noninvasive brain–computer interfaces’, IEEE Transactions on Biomedical Engineering 56(4), 1209–1219. Grozea, C., Voinescu, C. D. and Fazli, S. (2011), ‘Bristle-sensorslow-cost flexible passive dry eeg electrodes for neurofeedback and bci applications’, Journal of neural engineering 8(2), 025008. Guermeur, Y. and Monfrini, E. (2011), ‘A quadratic loss multi-class SVM for which a radius–margin bound applies’, Informatica 22(1), 73–96. Harandi, M., Salzmann, M. and Hartley, R. (2017), ‘Dimensionality reduction on SPD manifolds: The emergence of geometry-aware methods’, IEEE Transactions on Pattern Analysis and Machine Intelligence . Harland, C., Clark, T. and Prance, R. (2002), ‘Remote detection of human electroencephalograms using ultrahigh input impedance electric potential sensors’, Applied Physics Letters 81(17), 3284–3286. Hiraiwa, A., Shimohara, K. and Tokunaga, Y. (1990), ‘EEG topography recognition by neural networks’, IEEE Engineering in Medicine and Biology magazine 9(3), 39–42. Horev, I., Yger, F. and Sugiyama, M. (2016), Geometry-aware principal component analysis for symmetric positive definite matrices, in ‘Asian Conference on Machine Learning’, pp. 1–16. Hotelling, H. (1933), ‘Analysis of a complex of statistical variables into principal components.’, Journal of Educational Psychology 24(6), 417. Huang, Z., Wang, R., Shan, S., Li, X. and Chen, X. (2015), Log-Euclidean metric learning on symmetric positive definite manifold with application to image set classification, in ‘International Conference on Machine Learning’, pp. 720–729. Hyv ¨ arinen, A., Karhunen, J. and Oja, E. (2004), Independent component analysis, Vol. 46, John Wiley and Sons. Hyv ¨ arinen, A. and Oja, E. (1996), A neuron that learns to separate one independent component from linear mixtures, in ‘Proceedings IEEE International Conference on Neural Networks. Washington, DC’, pp. 62–67. 124
Jasper, H. (1958), ‘Report of the committee on methods of clinical examination in electroencephalography’, Electroencephalography and Clinical Neurophysiology 10 , 370– 375. Jones, M. C. and Sibson, R. (1987), ‘What is projection pursuit?’, Journal of the Royal Statistical Society. Series A (General) pp. 1–37. Jung, T. P., Makeig, S., Westerfield, M., Townsend, J., Courchesne, E. and Sejnowski, T. J. (2000), ‘Removal of eye activity artifacts from visual event-related potentials in normal and clinical subjects’, Clinical Neurophysiology 111(10), 1745–1758. Kang, H., Nam, Y. and Choi, S. (2009), ‘Composite common spatial pattern for subject- to-subject transfer’, Signal Processing Letters, IEEE 16(8), 683–686. Kawanabe, M., Samek, W., M ¨ uller, K. R. and Vidaurre, C. (2014), ‘Robust common spatial filters with a maxmin approach’, Neural Computation 26(2), 349–376. Kawanabe, M. and Vidaurre, C. (2009), Improving BCI performance by modified common spatial patterns with robustly averaged covariance matrices, in ‘World Congress on Medical Physics and Biomedical Engineering, September 7-12, 2009, Munich, Germany’, Springer, pp. 279–282. Kawanabe, M., Vidaurre, C., Scholler, S. and M ¨ uller, K. R. (2009), Robust common spatial filters with a maxmin approach, in ‘2009 Annual International Conference of the IEEE Engineering in Medicine and Biology Society’, IEEE, pp. 2470–2473. Koles, Z. J. (1991), ‘The quantitative extraction and topographic mapping of the abnormal components in the clinical EEG’, Electroencephalography and Clinical Neurophysiology 79(6), 440–447. Koles, Z., Lind, J. and Flor-Henry, P. (1994), ‘Spatial patterns in the background EEG underlying mental disease in man’, Electroencephalography and Clinical Neurophysiology 91(5), 319–328. Krauledat, M., Grzeska, K., Sagebaum, M., Blankertz, B., Vidaurre, C., M ¨ uller, K.-R. and Schr ¨ oder, M. (2009), Playing pinball with non-invasive BCI, in ‘Advances in Neural Information Processing Systems’, pp. 1641–1648. Lauer, F. and Guermeur, Y. (2011), ‘MSVMpack: a multi-class support vector machine package’, Journal of Machine Learning Research 12(Jul), 2293–2296. Lemm, S., Blankertz, B., Curio, G. and Muller, K. R. (2005), ‘Spatio-spectral filters for improving the classification of single trial EEG’, IEEE Transactions on Biomedical Engineering 52(9), 1541–1548. 125
Li, R. (2013), ‘Rayleigh quotient based optimization methods for eigenvalue problems’, In Summary of Lectures Delivered at Gene Golub SIAM Summer School 2013, Fudan University: Shanghai, China pp. 1–27. Liao, L.-D., Wang, I.-J., Chen, S.-F., Chang, J.-Y. and Lin, C.-T. (2011), ‘Design, fabrication and experimental validation of a novel dry-contact sensor for measuring electroencephalography signals without skin preparation’, Sensors 11(6), 5819–5834. Llera, A., G ´ omez, V. and Kappen, H. J. (2014), ‘Adaptive multiclass classification for brain computer interfaces’, Neural Computation 26(6), 1108–1127. Lotte, F. and Guan, C. (2010a), Learning from other subjects helps reducing braincomputer interface calibration time, in ‘2010 IEEE International Conference on Acoustics Speech and Signal Processing (ICASSP)’, IEEE, pp. 614–617. Lotte, F. and Guan, C. (2010b), Spatially regularized common spatial patterns for EEG classification, in ‘2010 20th International Conference on Pattern Recognition (ICPR)’, IEEE, pp. 3712–3715. Lotte, F. and Guan, C. (2011), ‘Regularizing common spatial patterns to improve BCI designs: unified theory and new algorithms’, IEEE Transactions on Biomedical Engineering 58(2), 355–362. Lu, H., Eng, H. L., Guan, C., Plataniotis, K. N. and Venetsanopoulos, A. N. (2010), ‘Regularized common spatial pattern with aggregation for EEG classification in smallsample setting’, IEEE Transactions on Biomedical Engineering 57(12), 2936–2946. Lu, H., Plataniotis, K. N. and Venetsanopoulos, A. N. (2009), Regularized common spatial patterns with generic learning for EEG signal classification, in ‘2009 Annual International Conference of the IEEE Engineering in Medicine and Biology Society’, IEEE, pp. 6599–6602. Machine Learning in Neural Engineering (2017(last update))). Makeig, S., Debener, S., Onton, J. and Delorme, A. (2004), ‘Mining event-related brain dynamics’, Trends in Cognitive Sciences 8(5), 204–210. M ¨ uller-Gerking, J., Pfurtscheller, G. and Flyvbjerg, H. (1999), ‘Designing optimal spatial filters for single-trial EEG classification in a movement task’, Clinical Neurophysiology 110(5), 787–798. Naeem, M., Brunner, C., Leeb, R., Graimann, B. and Pfurtscheller, G. (2006), ‘Seperability of four-class motor imagery data using independent components analysis’, Journal of Neural Engineering 3(3), 208. 126
Nguyen, T. H., Park, S. M., Ko, K. E. and Sim, K. B. (2012), Multi-class stationary CSP for optimal feature separation of brain source in BCI system, in ‘2012 12th International Conference on Control, Automation and Systems (ICCAS)’, IEEE, pp. 1035–1039. Nicolas-Alonso, L. F., Corralejo, R., Gomez-Pilar, J., ´ Alvarez, D. and Hornero, R. (2015), ‘Adaptive stacked generalization for multiclass motor imagery-based brain computer interfaces’, IEEE Transactions on Neural Systems and Rehabilitation Engineering 23(4), 702–712. Nishimori, Y. (1999), Learning algorithm for independent component analysis by geodesic flows on orthogonal group, in ‘International Joint Conference on Neural Networks, 1999. IJCNN’99’, Vol. 2, IEEE, pp. 933–938. Oehler, M., Neumann, P., Becker, M., Curio, G. and Schilling, M. (2008), Extraction of ssvep signals of a capacitive eeg helmet for human machine interface, in ‘Engineering in Medicine and Biology Society, 2008. EMBS 2008. 30th Annual International Conference of the IEEE’, IEEE, pp. 4495–4498. Park, J. and Chung, W. (2013), Common spatial patterns based on generalized norms, in ‘2013 International Winter Workshop on Brain-Computer Interface (BCI)’, IEEE, pp. 39–42. Pfurtscheller, G. (1999), ‘EEG event-related desynchronization (ERD) and event-related sybchronization (ERS)’, Electroencephalography: Basic principals, clinical applications and related fields pp. 958–967. Pfurtscheller, G., Brunner, C., Schl ¨ ogl, A. and Da Silva, F. L. (2006), ‘Mu rhythm (de) synchronization and EEG single-trial classification of different motor imagery tasks’, NeuroImage 31(1), 153–159. Pfurtscheller, G. and Da Silva, F. L. (1999), ‘Event-related EEG/MEG synchronization and desynchronization: basic principles’, Clinical Neurophysiology 110 (11), 1842–1857. Pfurtscheller, G., Flotzinger, D. and Kalcher, J. (1993), ‘Brain-computer interfacea new communication device for handicapped persons’, Journal of Microcomputer Applications 16(3), 293–299. Quang, M. H. (2016), ‘Infinite-dimensional Log-Determinant divergences II: Alpha-Beta divergences’, arXiv preprint arXiv:1610.08087 . Ramoser, H., M ¨ uller-Gerking, J. and Pfurtscheller, G. (2000), ‘Optimal spatial filtering of single trial EEG during imagined hand movement’, IEEE Transactions on Rehabilitation Engineering 8(4), 441–446. 127
Rezaei, S., Tavakolian, K., Nasrabadi, A. M. and Setarehdan, S. K. (2006), ‘Different classification techniques considering brain computer interface applications’, Journal of Neural Engineering 3(2), 139. Ruffini, G., Dunne, S., Fuentemilla, L., Grau, C., Farres, E., Marco-Pallar ´ es, J., Watts, P. and Silva, S. (2008), ‘First human trials of a dry electrophysiology sensor using a carbon nanotube array interface’, Sensors and Actuators A: Physical 144(2), 275–279. Salakhutdinov, R. and Mnih, A. (2008), Bayesian probabilistic matrix factorization using markov chain monte carlo, in ‘Proceedings of the 25th international conference on Machine learning’, ACM, pp. 880–887. Salvo, P., Raedt, R., Carrette, E., Schaubroeck, D., Vanfleteren, J. and Cardon, L. (2012), ‘A 3d printed dry electrode for ecg/eeg recording’, Sensors and Actuators A: Physical 174, 96–102. Samek, W., Binder, A. and M ¨ uller, K. R. (2013), Multiple kernel learning for braincomputer interfacing, in ‘2013 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)’, IEEE, pp. 7048–7051. Samek, W., Blythe, D., M ¨ uller, K. R. and Kawanabe, M. (2013), Robust spatial filtering with beta divergence, in ‘Advances in Neural Information Processing Systems’, pp. 1007– 1015. Samek, W., Kawanabe, M. and M ¨ uller, K. R. (2014), ‘Divergence-based framework for common spatial patterns algorithms’, IEEE Reviews in Biomedical Engineering 7, 50–72. Samek, W., Kawanabe, M. and Vidaurre, C. (2011), Group-wise stationary subspace analysis–a novel method for studying non-stationarities, in ‘Proceedings International Brain–Computer Interfacing Conference’, pp. 16–20. Samek, W. and M ¨ uller, K. R. (2015), Tackling noise, artifacts and nonstationarity in BCI with robust divergences, in ‘2015 23rd European Signal Processing Conference (EUSIPCO)’, IEEE, pp. 2741–2745. Samek, W., M ¨ uller, K.-R., Kawanabe, M. and Vidaurre, C. (2012), Brain-computer interfacing in discriminative and stationary subspaces, in ‘2012 Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)’, IEEE, pp. 2873–2876. Samek, W., Vidaurre, C., M ¨ uller, K. R. and Kawanabe, M. (2012), ‘Stationary common spatial patterns for brain–computer interfacing’, Journal of Neural Engineering 9(2), 026013. 128
Schl ¨ ogl, A., Lee, F., Bischof, H. and Pfurtscheller, G. (2005), ‘Characterization of four-class motor imagery EEG data for the BCI-competition 2005’, Journal of Neural Engineering 2(4), L14. Spilker, B., Kamiya, J., Callaway, E. and Yeager, C. L. (1969), ‘Visual evoked responses in subjects trained to control alpha rhythms’, Psychophysiology 5(6), 683–695. Sullivan, T. J., Deiss, S. R. and Cauwenberghs, G. (2007), A low-noise, non-contact eeg/ecg sensor, in ‘Biomedical Circuits and Systems Conference, 2007. BIOCAS 2007. IEEE’, IEEE, pp. 154–157. Tangermann, M., M ¨ uller, K. E., clinical neurophysiologyRobert, Aertsen, A., Birbaumer, N., Braun, C., Brunner, C., Leeb, R., Mehring, C., Miller, K. J., Mueller-Putz, G. et al. (2012), ‘Review of the BCI competition IV’, Frontiers in Neuroscience 6, 55. Tao, T. (2012), Topics in Random Matrix Theory, Vol. 132, American Mathematical Soc. The Divergence Methods Web Site (2013). Townsend, G., Graimann, B. and Pfurtscheller, G. (2006), ‘A comparison of common spatial patterns with complex band power features in a four-class BCI experiment’, IEEE Transactions on Biomedical Engineering 53(4), 642–651. Urig ¨ uen, J. A. and Garcia-Zapirain, B. (2015), ‘EEG artifact removal-state-of-the-art and guidelines’, Journal of Neural Engineering 12(3), 031001. Van Erp, J., Lotte, F. and Tangermann, M. (2012), ‘Brain-computer interfaces: beyond medical applications’, Computer 45(4), 26–34. Vidal, J. J. (1973), ‘Toward direct brain-computer communication’, Annual review of Biophysics and Bioengineering 2(1), 157–180. Von B ¨ unau, P., Meinecke, F. C., Kir ´ aly, F. C. and M ¨ uller, K. R. (2009), ‘Finding stationary subspaces in multivariate time series’, Physical Review Letters 103(21), 214101. Von B ¨ unau, P., Meinecke, F. C., Scholler, S. and M ¨ uller, K. R. (2010), Finding stationary brain sources in EEG data, in ‘2010 Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)’, IEEE, pp. 2810–2813. Wang, H. (2012), ‘Harmonic mean of Kullback–Leibler divergences for optimizing multi-class EEG spatio-temporal filters’, Neural Processing Letters 36(2), 161–171. Wang, H. and Li, X. (2016), ‘Regularized filters for L1-norm-based common spatial patterns’, IEEE Transactions on Neural Systems and Rehabilitation Engineering 24(2), 201–211. 129
Wang, H., Tang, Q. and Zheng, W. (2012), ‘L1-norm-based common spatial patterns’, IEEE Transactions on Biomedical Engineering 59(3), 653–662. Wojcikiewicz, W., Vidaurre, C. and Kawanabe, M. (2011a), Improving classification performance of BCIs by using stationary common spatial patterns and unsupervised bias adaptation, in ‘International Conference on Hybrid Artificial Intelligence Systems’, Springer, pp. 34–41. Wojcikiewicz, W., Vidaurre, C. and Kawanabe, M. (2011b), Stationary common spatial patterns: towards robust classification of non-stationary EEG signals, in ‘2011 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)’, IEEE, pp. 577–580. Wolpaw, J. and Wolpaw, E. W. (2012), Brain-computer Interfaces: Principles and Practice, Oxford University Press. Wu, W., Chen, Z., Gao, X., Li, Y., Brown, E. N. and Gao, S. (2015), ‘Probabilistic common spatial patterns for multichannel EEG analysis’, IEEE Transactions on Pattern Analysis and Machine Intelligence 37(3), 639–653. Xinyi Yong, R. K. W. and Birch, G. E. (2008), Robust common spatial patterns for EEG signal preprocessing, in ‘IEEE EMBS 30th Annual International Conference’, IEEE, pp. 2087– 2090. Xu, P., Yang, P., Lei, X. and Yao, D. (2011), ‘An enhanced probabilistic LDA for multi-class brain computer interface’, PloS one 6(1), e14634. Yong, X., Ward, R. K. and Birch, G. E. (2008a), Robust common spatial patterns for EEG signal preprocessing, in ‘30th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 2008. EMBS 2008.’, IEEE, pp. 2087–2090. Yong, X., Ward, R. K. and Birch, G. E. (2008b), Sparse spatial filter optimization for EEG channel reduction in brain-computer interface, in ‘IEEE International Conference on Acoustics, Speech and Signal Processing, 2008. ICASSP 2008’, IEEE, pp. 417–420. Zhang, H., Yang, H. and Guan, C. (2013), ‘Bayesian learning for spatial filtering in an EEG-based brain–computer interface’, IEEE Transactions on Neural Networks and Learning Systems 24(7), 1049–1060. Zhou, B., Wu, X., Zhang, L., Lv, Z. and Guo, X. (2014), ‘Robust spatial filters on three-class motor imagery EEG data using independent component analysis’, Journal of Biosciences and Medicines 2(02), 43. 130
LIST OF PUBLICATIONS 1. Deepa Beeta, Th., Sergio, C., Javier, O., Andrzej C.,2017. Optimization of Aplha- Beta Log-Det Divergences and their Applications in the Spatial Filtering of Two Class Motor Imagery Movements. Entropy ,19(3): 89, pp. 1-40. 2. Deepa Beeta, Th., Rajkumar E.R., 2016. A Comparative Performance Analysis for Classification of Multiclass Motor Imagery Movements. International Journal of Control Theory and Applications, 9(36), pp. 443-450. 3. Deepa Beeta, Th., Rajkumar E.R.,2016. Common Spatial Pattern Algorithm Based Signal Processing Techniques for Classification of Motor Imagery Movements: A Mini Review. International Journal of Control Theory and Applications (Accepted). 4. Deepa Beeta, Th., Sergio, C., Rajkumar E.R., 2016. ThinICA-CSP algorithm for discrimination of multiclass motor imagery movements. In Region 10 Conference (TENCON), IEEE,pp. 2483-2486. 5. Pablo, A.,Deepa Beeta, Th., Auxiliadora, S.,Irene, F., 2014. Cancelacin de interferencias en EEG mediante el Anlisis de Componentes Acotadas. In Congreso Anual de la Sociedad Espaola de Ingeniera Biomdica (CASEIB), pp. 1-4. 6. Deepa Beeta, Th., Rajkumar E.R., 2014. Development of Cost Effective Backscattered Optical Imaging System for Detection of Abnormality in Tissues. In IEEE Engineering in Medicine and Biology Society (EMB), IEEE. (Poster). 131
Appendices 132