# Embedding Tangent Space Extreme Learning Machine for EEG Decoding in Brain Computer Interface Systems

Authors: Mingwei Zhang; Yao Hou; Rongnian Tang; Youjun Li
Journal: Journal of Control Science and Engineering
 (2021)

---

## Abstract

In motor imagery brain computer interface system, the spatial covariance matrices of EEG signals which carried important discriminative information have been well used to improve the decoding performance of motor imagery. However, the covariance matrices often suffer from the problem of high dimensionality, which leads to a high computational cost and overfitting. These problems directly limit the application ability and work efficiency of the BCI system. To improve these problems and enhance the performance of the BCI system, in this study, we propose a novel semisupervised locality-preserving graph embedding model to learn a low-dimensional embedding. This approach enables a low-dimensional embedding to capture more discriminant information for classification by efficiently incorporating information from testing and training data into a Riemannian graph. Furthermore, we obtain an efficient classification algorithm using an extreme learning machine (ELM) classifier developed on the tangent space of a learned embedding. Experimental results show that our proposed approach achieves higher classification performance than benchmark methods on various datasets, including the BCI Competition IIa dataset and in-house BCI datasets.


---

## Body

## 1. Introduction

Brain computer interface (BCI) systems have been developed to provide a novel and promising alternative method for humans to interact with their environment. Among currently operational BCI systems, many have been based on electroencephalography (EEG) signals, owing to the noninvasive nature of the technique and affordable data acquisition equipment, both of which also facilitate real-time operation [1]. Although many EEG signals such as P300 waves or steady-state visual evoked potentials (SSVEP) could be used for BCI [2], in this study, we focus on motor imagery-based BCI systems, with which trained subjects can voluntarily generate EEG control signals by imagining movements of different parts of their bodies [3]. With the advantage not requiring external stimulation, motor imagery-based BCI has proven effective in practice.The common spatial pattern (CSP) algorithm has been shown to be one of the most effective methods for feature extraction from motor imagery signals [4]. The success of CSP techniques in BCI applications depends heavily on learning a proper spatial filter from a training dataset. Some extensions of CSP were proposed to improve the performance of feature extraction in motor imagery classification. In [5], a novel regularized CSP combined with correlation-based channel selection was proposed to extract effective features of EEG. In [6], a discriminative canonical pattern matching combined with CSP was used to extract the movement-related features. In [7], to capture the shared salient information across multiple related spatial patterns, a multiview learning-based sparse optimization was proposed to jointly extract robust CSP features with L1L2-norm regularization. In [8], a generic model training by weighted linear discriminant analysis was proposed to reduce the correction time of the BCI system. In [9], a fusion method for internal feature selection was proposed to fuse different feature selection rules. It can consume less computational cost and obtain high performance of motor imagery classification. Typically, in such motor imagery classification system, some discriminant information contained in the testing dataset remains unutilized. However, graph-based models are able to utilize such information because they naturally incorporate diverse types of information, such as the relationship between training and testing datasets [10].Graph construction is an essential step in graph-based models. Two conventional methods are commonly used, including thek-nearest neighbor method and the ε-ball-based method [11]:(1)
k-nearest neighbor methods consider data samples xi and xj as neighbors if xi is among the k nearest neighbors of xi or xj is among the k nearest neighbors of xi, where the nearest neighbors are determined as the samples with the smallest Euclidean distance(2)
ε-ball neighborhood methods consider data samples xi and xj as neighbors if xi−xj2<ε, where ⋅2 is the Euclidean norm and ε∈ℝThe above traditional graph construction methods are relatively simple to implement and have been applied effectively in practice. However, if the data samples considered do not lie in a Euclidean space, thek-nearest neighbor and ε-ball-based methods cannot guarantee the connectivity of the whole graph and often result in several separated subgraphs because neighbors defined by Euclidean distance cannot characterize the real geometrical relations among samples well [12]. An obvious example is the covariance feature data used in motor imagery-based BCI. Because covariance features of EEG signals with the form of symmetric positive-definite (SPD) matrices lie on a Riemannian manifold [13], we propose an extension of the graph construction approach based on Riemannian geodesic distance to construct a proper graph for EEG signals, namely, a Riemannian graph. In fact, EEG classification techniques on high-dimensional Riemannian manifolds have recently received increased attention [14]. For example, an efficient classification algorithm called tangent space linear discriminant analysis (TS + LDA), was proposed in [14]. TS + LDA mapped all the data points in a Riemannian manifold into its tangent space, which is known as the best hyperplane for classification tasks [15], and then applied linear discriminant analysis (LDA) for classification in the tangent space.After constructing a Riemannian graph for an EEG signal, many graph-based algorithms can be exploited on the constructed Riemannian graph. Such graph embedding algorithms are among the most successful in graph-based data processing on high-dimensional Riemannian manifolds because the most important information is always contained in a low-dimensional embedding of a high-dimensional Riemannian manifold. Many state-of-the-art semisupervised and unsupervised algorithms have been proposed. Semisupervised learning methods can utilize both training and test data to boost the algorithmic performance. Many related works, such as [16–18], were proposed to classify EEG based on training and test data in the context of semisupervised learning. In [16], a self-training support vector machine (SVM) method was proposed to use the test data for updating the parameter of SVM model. In [17], an iterative semisupervised SVM method was proposed for channel selection and classifier training. In [18], a semisupervised expectation maximization (EM) was incorporated into CSP to extract more effective features. On the contrary, many unsupervised dimensionality reduction methods also are proposed to learn a low-dimensional embedding preserving similarities between vertex pairs on such graphs, including locally linear embedding [19], locality-preserving projection [20], and Laplacian eigenmaps [21].Considering locality preservation, in this study, we propose a bilinear locality-preserving embedding learning algorithm on a constructed Riemannian graph. Because the learned embedding was also an SPD matrix space, inspired by the TS + LDA algorithm, we applied an extreme learning machine (ELM) [22, 23] method on the embedding tangent space for classification. Many related works of graph embedding methods and ELM-based methods have been recently proposed. In [24], a hybrid-graph discriminant learning method was proposed to reveal the complex high-order relationship of the original data and to extract the discriminant features. In [25], a local constrained manifold structure collaborative preserving embedding was proposed to enhance the intraclass compactness and the interclass separability. In [26], to improve the classification performance of motor imagery, a sparse Bayesian ELM algorithm was proposed to automatically control the model complexity and exclude redundant hidden neurons. In [27], a multikernel ELM method integrating two types of kernels was proposed to effectively explore the supplementary information from multiple nonlinear feature spaces.In [28], a clustering-based multitask feature learning method was proposed to improve EEG pattern decoding. Most of these studies have achieved great success in graph embedding and classification. However, the above methods learn embedding or classifier on Euclidean space and ignore that the covariance matrices lie on Riemannian manifold. The contributions of this study are threefold.(1)
A graph model of the covariance features of motor imagery signals based on Riemannian geodesic distance was constructed. The neighborhood relationships among data samples were characterized by Riemannian geodesics because the EEG covariance features endowed with a Riemannian metric formed a Riemannian manifold.(2)
A novel framework of graph embedding is proposed for the constructed Riemannian graph. The calculation of graph embedding is formulated as an eigenvalue problem.(3)
A classification algorithm, called an embedding tangent space extreme learning machine (ETS + ELM), is proposed to learn low-dimensional embeddings. Experimental results show that the proposed algorithm achieved higher classification performance than prior methods used for comparisonThe remainder of this paper is organized as follows. In Section2, some basic concepts concerning SPD matrix spaces and BCI data models are briefly reviewed. In Section 3, we introduce the construction of Riemannian graphs and propose a bilinear locality-preserving graph embedding algorithm and an embedding tangent space ELM for classification. Extensive experimental results are presented in Section 4 to demonstrate the effectiveness of the proposed method. Finally, some conclusions are presented in Section 5.


## 2. Concepts of Riemannian Geometry and BCI Data Models

### 2.1. Riemannian Geometry

In mathematics, Riemannian manifoldℳ is defined as a smooth manifold equipped with a Riemannian metric ,, i.e., a smoothly varying inner product ,p,p∈M, on its tangent space [29]. A smooth manifold is a topological manifold with a smooth structure locally similar to Euclidean space. Every point on this topological manifold has a neighborhood for which there exists a homeomorphism mapping the neighborhood to a Euclidean space. The materials and methods section should contain sufficient detail so that all procedures can be repeated. It may be divided into headed subsections if several methods are described.For many applications, the space of SPD matrices endowed with a Riemannian metric is one of most important Riemannian manifolds [30]. Owing to the properties of SPD matrices, explicit formulae exist for the major concepts of Riemannian geometry. Therefore, these concepts can be applied relatively easily in the context of pattern recognition.


### 2.2. Space of Symmetric Positive-Definite Matrices

Denoting the space of symmetric matrices by(1)SN=P∈ℝN×N,P=PT,and the space of positive-definite matrices by(2)PN=P∈ℝN×N,uTPu>0,∀u∈ℝN,the space of SPD matrices is defined as(3)SPDN=SN∩PN.The space of SPD matrices lies on a differentiable Riemannian manifold [30]. Hence, all the mathematical tools defined in Riemannian geometry can be applied to any SPDN.The Riemannian distance between two matricesP1,P2∈calSPDN is the minimum length of the curves connecting P1 and P2 on a Riemannian manifold [31]. There are many possible mathematical definitions of the Riemannian distance. In this study, we adopt the following formula for Riemannian distance:(4)δRP1,P2=logP1−1P2F=∑i=1Nlog2ηi1/2,where ⋅F is the Frobenius norm of a matrix and ηi is the ith real eigenvalue of P1−1P2.The associated Euclidean tangent space is critical in the analysis of a Riemannian manifold, as illustrated in Figure1. As a preliminary step to define the tangent space, we first, respectively, introduce the logarithmic mapping and exponential mapping operators as(5)LogPPj=Sj=P1/2logP−1/2PjP−1/2P1/2,ExpPSj=Pj=P1/2expP−1/2SjP−1/2P1/2,where LogPPj is a mapping from a point on the manifold to the tangent space at P, while ExpPSj is a mapping from the tangent space at P to the manifold. These two operators are both one-to-one mapping operators between the Riemannian manifold and the tangent space.Figure 1
Illustration of Riemannian manifold and tangent space.Because the tangent spaceTPN=LogPPj,Pj∈SPDN is the space of symmetric matrices, it comprises only N(N + 1)/2 independent elements. A minimal representation of the tangent space TPN [10] is defined as(6)TN=sj=upperP−1/2LogPPjP−1/2∈ℝNN+1/2,where the upper ⋅ operator is used to retain the upper triangular portion of the symmetric matrix and obtain its vectorization.The mean of the SPD matrices plays an important role in classification. Since the neighborhood of Riemannian manifold is local homeomorphic to its tangent space, the trinational classifier can be performed on tangent space to obtain high classification performance [14]. However, large neighborhood will lead to large distortion between neighborhood and its tangent space, particularly when the selection of tangent point is not representative. The mean of SPD matrices is the best representation of all data points. It is defined as the point Pm=SPDN with the minimum sum of squared distances to all SPD matrices in a dataset C:(7)Pm=argminP∈calSPDN∑Pi∈CδR2P,Pi.


### 2.3. BCI Data Model

The recorded EEG signal was a multilag/multichannel signal:(8)Xt=xt,…,xt+L−1∈ℝN×L,where N and L denote the number of channels and sampled points, respectively, and xt=x1t,…,xNtT∈ℝN×1 is the snapshot vector. In motor imagery BCIs, the second-order statistical information (covariance feature) of Xt often provides discriminative information for classification. The covariance feature of Xt is represented by its sample covariance matrix (SCM) [32] as(9)Pt=1L−1XtXTt.


## 2.1. Riemannian Geometry

In mathematics, Riemannian manifoldℳ is defined as a smooth manifold equipped with a Riemannian metric ,, i.e., a smoothly varying inner product ,p,p∈M, on its tangent space [29]. A smooth manifold is a topological manifold with a smooth structure locally similar to Euclidean space. Every point on this topological manifold has a neighborhood for which there exists a homeomorphism mapping the neighborhood to a Euclidean space. The materials and methods section should contain sufficient detail so that all procedures can be repeated. It may be divided into headed subsections if several methods are described.For many applications, the space of SPD matrices endowed with a Riemannian metric is one of most important Riemannian manifolds [30]. Owing to the properties of SPD matrices, explicit formulae exist for the major concepts of Riemannian geometry. Therefore, these concepts can be applied relatively easily in the context of pattern recognition.


## 2.2. Space of Symmetric Positive-Definite Matrices

Denoting the space of symmetric matrices by(1)SN=P∈ℝN×N,P=PT,and the space of positive-definite matrices by(2)PN=P∈ℝN×N,uTPu>0,∀u∈ℝN,the space of SPD matrices is defined as(3)SPDN=SN∩PN.The space of SPD matrices lies on a differentiable Riemannian manifold [30]. Hence, all the mathematical tools defined in Riemannian geometry can be applied to any SPDN.The Riemannian distance between two matricesP1,P2∈calSPDN is the minimum length of the curves connecting P1 and P2 on a Riemannian manifold [31]. There are many possible mathematical definitions of the Riemannian distance. In this study, we adopt the following formula for Riemannian distance:(4)δRP1,P2=logP1−1P2F=∑i=1Nlog2ηi1/2,where ⋅F is the Frobenius norm of a matrix and ηi is the ith real eigenvalue of P1−1P2.The associated Euclidean tangent space is critical in the analysis of a Riemannian manifold, as illustrated in Figure1. As a preliminary step to define the tangent space, we first, respectively, introduce the logarithmic mapping and exponential mapping operators as(5)LogPPj=Sj=P1/2logP−1/2PjP−1/2P1/2,ExpPSj=Pj=P1/2expP−1/2SjP−1/2P1/2,where LogPPj is a mapping from a point on the manifold to the tangent space at P, while ExpPSj is a mapping from the tangent space at P to the manifold. These two operators are both one-to-one mapping operators between the Riemannian manifold and the tangent space.Figure 1
Illustration of Riemannian manifold and tangent space.Because the tangent spaceTPN=LogPPj,Pj∈SPDN is the space of symmetric matrices, it comprises only N(N + 1)/2 independent elements. A minimal representation of the tangent space TPN [10] is defined as(6)TN=sj=upperP−1/2LogPPjP−1/2∈ℝNN+1/2,where the upper ⋅ operator is used to retain the upper triangular portion of the symmetric matrix and obtain its vectorization.The mean of the SPD matrices plays an important role in classification. Since the neighborhood of Riemannian manifold is local homeomorphic to its tangent space, the trinational classifier can be performed on tangent space to obtain high classification performance [14]. However, large neighborhood will lead to large distortion between neighborhood and its tangent space, particularly when the selection of tangent point is not representative. The mean of SPD matrices is the best representation of all data points. It is defined as the point Pm=SPDN with the minimum sum of squared distances to all SPD matrices in a dataset C:(7)Pm=argminP∈calSPDN∑Pi∈CδR2P,Pi.


## 2.3. BCI Data Model

The recorded EEG signal was a multilag/multichannel signal:(8)Xt=xt,…,xt+L−1∈ℝN×L,where N and L denote the number of channels and sampled points, respectively, and xt=x1t,…,xNtT∈ℝN×1 is the snapshot vector. In motor imagery BCIs, the second-order statistical information (covariance feature) of Xt often provides discriminative information for classification. The covariance feature of Xt is represented by its sample covariance matrix (SCM) [32] as(9)Pt=1L−1XtXTt.


## 3. Proposed Bilinear Locality-Preserving Embedding and ELM Classification

### 3.1. Riemannian Graph Construction

In SPD matrix spaces, because the data lie on a Riemannian manifold, it is natural to construct a Riemannian graphGR=V,W via Riemannian geodesic distance, as illustrated in Figure 2. V represents the set of graph vertices, and W represents the graph weight matrix. For a given dataset, the Riemannian graph is constructed as follows.Figure 2
Riemannian graph from Riemannian manifold.First, to construct an adjacency graph, a Riemannian distance-basedk nearest neighbor method is used. The points Pi,Pj∈V are considered neighbors if Pi is among the k nearest neighbors of Pi or Pj is among the k nearest neighbors of Pi, where k is a positive integer and k nearest neighbors are measured by the Riemannian geodesic distance dij=δRPi,Pj.Second, to calculate graph weights, many approaches may be used, such as heat kernels [33], inverse Euclidean distance [34], and local linear reconstruction coefficients [19]. In this study, we use a heat kernel based on the Riemannian geodesic distance to construct graph weights:(10)Wije−dij2/2σ2,ifPiandPjare neighbors,0,otherwise,where Wij is the i,jth element of W and the parameter σ controls the width of the neighborhoods.


### 3.2. Bilinear Locality-Preserving Embedding and ELM on Tangent Space

Because the covariance features studied in this paper lie in SPD matrix spaces, we propose a bilinear mapping withV∈ℝM×N and VT under the constraint VVT=IM, to construct a system learning a low-dimensional embedding. The low-dimensional embedding Pe=VPVT∈ℝM×M was expected to preserve the locality geometry of a high-dimensional Riemannian manifold. Locality preservation attempts to ensure that if Pi and Pj in a high-dimensional manifold are “close,” then VPiVT and VPjVT in low-dimensional embedding are also “close.” A reasonable criterion for locality preservation is to minimize the following objective function:(11)min∑i,j=1CVPiVT−VPjVTF2Wij,where the F-norm ⋅F is used to measure the “closeness” of data points in the embedding.There is no closed-form optimal solution for the optimization problem (11). It should be noted that the objective function (11) is the trade-off model between optimized solution and Riemannian locality preserving. In this study, we propose an iterative strategy to obtain a local optimal solution of (11). In the proposed strategy, the matrix V is first initialized as V0=IM,0, and the new optimal Vt+1 of the t+1th iteration can be obtained by solving the following eigendecomposition problem, given the tth iteration of V and Vt:(12)max−trVt+1LVt+1T,where L=∑i,j=1CPi−PjVtTVtPi−PjTWij. Negating (11), the minimization problem becomes a maximization problem. Because the optimization problem in (12) is a generalized eigenvector problem [35, 36], the optimal transformation matrix Vt+1 is the first M generalized eigenvectors of L corresponding to the M largest eigenvalues. The iteration is repeated until the solution converges, and a local optimal solution of (11) is obtained. Further details are provided in Algorithm 1.Algorithm 1: Embedding tangent space extreme learning machine (ETS + ELM).
Input: given SPD matrix set C with samples Ρi∈ℝN×N,i=1,…,C, dimensions of embedding M, iterative number Nite, stop threshold τ, and ELM parameter ai,bi,Nh, and No;(1)
Construct a Riemannian graphGR over all samples based on Riemannian geodesic distance (Section 3.1);(2)
Initialize:V0=IM,0;(3)
Fort = 1 : 1:(Nite), do;Calculate the matricesL (Section 3.2);Obtain mapping matrixVt+1∈ℝM×N by (12);IfVt−Vt−1F2≤τbreak;end if;end for;(4)
ConstructM×M-dimensional embedding Ρei=VoptΡiVoptT,i=1,…,C, where Vopt is convergent matrix in step 4;(5)
Calculate the Riemannian mean of all data on embedding:Ρcm=argminΡc∑iδR2Ρc,Ρci(6)
Project embedding point onto tangent space on Riemannian mean by (6) TN=sj;(7)
Apply a three-layer ELM in the tangent space, learning the weightβ from (14), using the training data;(8)
The label of the testing data may be obtained byHtβ, where Ht is the mapping matrix of the testing data;Output: the label of the testing data.After learning the optimized bilinear mapping matrixVopt, the proposed approach maps all points in the embedding into the tangent space according to (6). It has been reported that the tangent space located at the Riemannian mean is the best hyperplane for classification, as it is associated with the minimal distance loss [15]. A three-layer ELM network was constructed on a tangent space for classification because the ELM classifier is more flexible to handle complex and multivariate data than LDA and SVM classifier. The mapping function from the jth input tangent vector sj to the ith hidden node was defined as(13)hisj=gsj;θ=11+exp−aiTsj+bi,where θ=ai,bi denotes the parameters of the sigmoid mapping function. Based on the principles of ELMs, we randomly chose parameters ai and bi within the range [−1, 1]. We denote by hsj∈ℝ1×Nh the output vector of the hidden layer with respect to sj and by β∈ℝNh×No the output weights that connect the hidden layer to the output layer, where Nh and No are the number of hidden neurons and output neurons, respectively. Moreover, the output weights β can be learned by solving an optimization problem [22]:(14)minHβ−Y22+λβ22,where H=hs1T,…,hsjT,…,hsCTT and λ parameters are used to balance the loss and L2 is a regularizer. The parameter Y is a matrix of the target output, and the jth row of Y is a binary vector with only one entry (corresponding to the class to which sj belongs) equal to one.The framework of ETS + ELM algorithm is shown in Figure3; the proposed ETS + ELM algorithm includes a Riemannian graph construction for data representation, locality preserving with bilinear mapping for dimensionality reduction, tangent space mapping for feature extraction, and 3-layer extreme learning machine network for classification.Figure 3
The framework of ETS + ELM algorithm.


## 3.1. Riemannian Graph Construction

In SPD matrix spaces, because the data lie on a Riemannian manifold, it is natural to construct a Riemannian graphGR=V,W via Riemannian geodesic distance, as illustrated in Figure 2. V represents the set of graph vertices, and W represents the graph weight matrix. For a given dataset, the Riemannian graph is constructed as follows.Figure 2
Riemannian graph from Riemannian manifold.First, to construct an adjacency graph, a Riemannian distance-basedk nearest neighbor method is used. The points Pi,Pj∈V are considered neighbors if Pi is among the k nearest neighbors of Pi or Pj is among the k nearest neighbors of Pi, where k is a positive integer and k nearest neighbors are measured by the Riemannian geodesic distance dij=δRPi,Pj.Second, to calculate graph weights, many approaches may be used, such as heat kernels [33], inverse Euclidean distance [34], and local linear reconstruction coefficients [19]. In this study, we use a heat kernel based on the Riemannian geodesic distance to construct graph weights:(10)Wije−dij2/2σ2,ifPiandPjare neighbors,0,otherwise,where Wij is the i,jth element of W and the parameter σ controls the width of the neighborhoods.


## 3.2. Bilinear Locality-Preserving Embedding and ELM on Tangent Space

Because the covariance features studied in this paper lie in SPD matrix spaces, we propose a bilinear mapping withV∈ℝM×N and VT under the constraint VVT=IM, to construct a system learning a low-dimensional embedding. The low-dimensional embedding Pe=VPVT∈ℝM×M was expected to preserve the locality geometry of a high-dimensional Riemannian manifold. Locality preservation attempts to ensure that if Pi and Pj in a high-dimensional manifold are “close,” then VPiVT and VPjVT in low-dimensional embedding are also “close.” A reasonable criterion for locality preservation is to minimize the following objective function:(11)min∑i,j=1CVPiVT−VPjVTF2Wij,where the F-norm ⋅F is used to measure the “closeness” of data points in the embedding.There is no closed-form optimal solution for the optimization problem (11). It should be noted that the objective function (11) is the trade-off model between optimized solution and Riemannian locality preserving. In this study, we propose an iterative strategy to obtain a local optimal solution of (11). In the proposed strategy, the matrix V is first initialized as V0=IM,0, and the new optimal Vt+1 of the t+1th iteration can be obtained by solving the following eigendecomposition problem, given the tth iteration of V and Vt:(12)max−trVt+1LVt+1T,where L=∑i,j=1CPi−PjVtTVtPi−PjTWij. Negating (11), the minimization problem becomes a maximization problem. Because the optimization problem in (12) is a generalized eigenvector problem [35, 36], the optimal transformation matrix Vt+1 is the first M generalized eigenvectors of L corresponding to the M largest eigenvalues. The iteration is repeated until the solution converges, and a local optimal solution of (11) is obtained. Further details are provided in Algorithm 1.Algorithm 1: Embedding tangent space extreme learning machine (ETS + ELM).
Input: given SPD matrix set C with samples Ρi∈ℝN×N,i=1,…,C, dimensions of embedding M, iterative number Nite, stop threshold τ, and ELM parameter ai,bi,Nh, and No;(1)
Construct a Riemannian graphGR over all samples based on Riemannian geodesic distance (Section 3.1);(2)
Initialize:V0=IM,0;(3)
Fort = 1 : 1:(Nite), do;Calculate the matricesL (Section 3.2);Obtain mapping matrixVt+1∈ℝM×N by (12);IfVt−Vt−1F2≤τbreak;end if;end for;(4)
ConstructM×M-dimensional embedding Ρei=VoptΡiVoptT,i=1,…,C, where Vopt is convergent matrix in step 4;(5)
Calculate the Riemannian mean of all data on embedding:Ρcm=argminΡc∑iδR2Ρc,Ρci(6)
Project embedding point onto tangent space on Riemannian mean by (6) TN=sj;(7)
Apply a three-layer ELM in the tangent space, learning the weightβ from (14), using the training data;(8)
The label of the testing data may be obtained byHtβ, where Ht is the mapping matrix of the testing data;Output: the label of the testing data.After learning the optimized bilinear mapping matrixVopt, the proposed approach maps all points in the embedding into the tangent space according to (6). It has been reported that the tangent space located at the Riemannian mean is the best hyperplane for classification, as it is associated with the minimal distance loss [15]. A three-layer ELM network was constructed on a tangent space for classification because the ELM classifier is more flexible to handle complex and multivariate data than LDA and SVM classifier. The mapping function from the jth input tangent vector sj to the ith hidden node was defined as(13)hisj=gsj;θ=11+exp−aiTsj+bi,where θ=ai,bi denotes the parameters of the sigmoid mapping function. Based on the principles of ELMs, we randomly chose parameters ai and bi within the range [−1, 1]. We denote by hsj∈ℝ1×Nh the output vector of the hidden layer with respect to sj and by β∈ℝNh×No the output weights that connect the hidden layer to the output layer, where Nh and No are the number of hidden neurons and output neurons, respectively. Moreover, the output weights β can be learned by solving an optimization problem [22]:(14)minHβ−Y22+λβ22,where H=hs1T,…,hsjT,…,hsCTT and λ parameters are used to balance the loss and L2 is a regularizer. The parameter Y is a matrix of the target output, and the jth row of Y is a binary vector with only one entry (corresponding to the class to which sj belongs) equal to one.The framework of ETS + ELM algorithm is shown in Figure3; the proposed ETS + ELM algorithm includes a Riemannian graph construction for data representation, locality preserving with bilinear mapping for dimensionality reduction, tangent space mapping for feature extraction, and 3-layer extreme learning machine network for classification.Figure 3
The framework of ETS + ELM algorithm.


## 4. Experimental Results

### 4.1. Experimental Setup

#### 4.1.1. Data Description

In this study, the BCI Competition 2008 motor imagery data [37] and our in-house experimental data are used for performance evaluation.(1)
The BCI Competition 2008 IIa EEG motor imagery data were recorded from nine subjects who performed four kinds of motor imagery actions (imagined right hand, left hand, foot, and tongue movements) with 22 EEG channels. For each subject, there were 144 trials per class, including training and testing trials. The EEG signals were filtered by an 8−30 Hz bandpass filter.(2)
Our in-house EEG data were recorded from 12 subjects with 64 EEG channels. For each subject, there were 234 trials per class. Our in-house experiments tested right and left hand imagined movements.The sampling rate of the EEG signal was 250 Hz. The bandpass filter was the same as the Competition dataset.


#### 4.1.2. Algorithms Evaluated

We evaluate the proposed algorithms (ETS + ELM) against five competing algorithms.(1)
CSP + LDA: the CSP in [4] and linear discriminant analysis (LDA) in [38] were applied for evaluation(2)
CSP + SVM: the number of CSP spatial filters was set to 8, as suggested in [4], for feature extraction, and a support vector machine (SVM) [39] was used for classification(3)
MDRM: minimum distance to Riemannian mean [14] was used for evaluation(4)
TS + LDA : LDA was used to classify signals on a tangent space [14]


#### 4.1.3. Parameter Setting

In all experiments, we empirically setk=40 and σ=1 during Riemannian graph construction. In ETS + ELM, the number of iterations was set to Nite=8, and the iterative stop threshold to τ=0.01. In fact, most information of the original space was lost during the process, while the number of embedding dimensions was too small. However, all of the methods evaluated suffer from computational cost problems under conditions of excessively high dimensionality. Hence, the dimension M was empirically set to 10. Because the ELM method suffers from an overfitting problem in the case of excessively hidden neurons, a regularizer parameter was used to control the overfitting problem. Therefore, the regularizer parameter λ was empirically set to 0.001, while the number of hidden neurons was 2000.


### 4.2. Experiments for EEG Classification

To evaluate the performance of the proposed algorithms, ETS + ELM, some competing algorithms (CSP + LDA, CSP + SVM, MDRM [14], and TS + LDA [14]) were compared on the EEG datasets from the 2008 BCI Competition and our in-house BCI system.First, we evaluated the convergence performance of the iterative algorithm for optimization problem (12). Because the Kappa coefficient [40] scales linearly from 0 to 1 in the range between random and perfect classification, it was selected as a performance measure [40] for the classification algorithms in this study. In each iteration step of (12), we were able to obtain a learned mapping matrix Vt and a corresponding low-dimensional embedding. The performance of the ETS + ELM classification algorithm under the same parameter settings was evaluated on these embeddings. The average performance of nine subjects on the BCI Competition IIa dataset is shown in Figure 4. It is clear that the mapping matrix Vt learned from the proposed algorithm converged quickly. Normally, it required only two iteration steps.Figure 4
Iteration vs. mean Kappa value of ETS + ELM for BCI competition 2008 IIa (the Kappa value is the mean of nine subjects).We also compared the performance of various methods on classification of the BCI Competition 2008 data. In Table1, the performance of the proposed methods, MDRM, TS + LDA, and the top three Competition winners for dataset IIa of BCI Competition 2008 (1st of BCI08, 2nd of BCI08, and 3rd of BCI08) are presented for comparison. From the results shown in Table 1, it may be observed that the proposed ETS + ELM method achieves a mean performance of 0.577, which was the best performance of all the methods compared.Table 1
Classification results of BCI competition 2008 IIa in kappa values.
MethodsMean KappaSubjectsS01S02S03S04S05S06S07S08S09ETS + ELM0.5770.740.270.790.520.310.310.690.770.79TS + ELM0.5730.730.320.770.510.320.310.670.760.771st of BCI080.5700.680.420.750.480.400.270.770.750.61TS + LDA [14]0.5670.740.380.720.500.260.340.690.710.76MDRM [14]0.5210.750.370.660.530.290.270.560.580.682nd of BCI080.5200.690.340.710.440.160.210.660.730.693rd of BCI080.3100.380.180.480.330.070.140.290.490.44We also tested the performance of the proposed algorithms in terms of classification accuracy on the BCI Competition 2008 IIa dataset with a 30-fold cross-validation procedure. The cross-validation reflects the generalization performance of the proposed method. From the results in Table2, it is clear that, for all subjects, ETS + ELM outperformed the reference methods. This indicates that, by finding a low-dimensional embedding with the property of locality preservation, the classification method ETS + ELM exhibited the best performance.Table 2
Classification accuracy of BCI competition 2008 IIa in 30-fold cross-validation (%).
SubjectMethodsETS + ELMTS + ELMTS + LDA [14]MDRM [14]CSP + LDA [35]CSP + SVM [36]S0181.580.580.577.878.376.3S0260.358.051.344.144.750.7S0388.186.887.576.882.285.1S0464.160.059.354.959.152.9S0555.459.34543.839.748.8S0657.153.255.347.150.149.2S0782.880.782.1728178.1S0885.080.784.875.268.577.4S0987.387.486.176.677.482.2Mean ± std73.5 ± 13.271.8± 13.870.2 ± 17.163.2 ± 15.264.6 ± 16.666.7 ± 15.7Moreover, we evaluated the performance of all algorithms on the dataset recorded from our in-house BCI systems. As shown in Table3, the proposed ETS + ELM showed a significant performance improvement over the prior methods (an approximately 2% improvement in classification accuracy), even when compared with the benchmark algorithm CSP [4]. In conclusion, the experimental results on the competition dataset and in-house dataset demonstrate that the proposed ETS + ELM approached achieved a higher performance than the other methods. In additional, we provide the ANOVA and t-test analysis of Tables 1–3. The repeated measures ANOVA results for Tables 1–3 are p = 0.2 × 10−10, 0.7 × 10−7, and 0.03, respectively. It indicates that all studied methods have statistically significant difference on classification performance. The t-test analysis of Tables 1–3 are shown in Table 4.Table 3
Classification accuracy on in-house BCI dataset (%).
SubjectMethodsETS + ELMTS + ELMTS + LDA [14]MDRM [14]CSP + LDA [45]CSP + SVM [36]A0110010010010098.296.4A0298.493.797.493.598.796.1A0396.689.093.596.787.083.8A0491.085.091.284.289.482.4A0587.180.184.280.778.968.4A0690.591.273.267.882.175.0A0773.271.266.067.867.866.0A0876.371.078.971.970.173.6A0982.275.572.667.579.479.4A1094.497.789.690.993.596.1A11100100100100100100A1272.985.388.785.477.474.1Mean ± std88.4 ± 10.386.6± 10.486.2 ± 11.383.8 ± 12.685.2 ± 10.582.6 ± 11.3Table 4
Thet-test results of Tables 1–3.
Table 1Table 2Table 3ETS + ELM vs. TS + LDA∽∗∽ETS + ELM vs. MDRM∽∗∗∽ETS + ELM vs. CSP + LDA—∗∗∗ETS + ELM vs. CSP + SVM—∗∗∗ETS + ELM vs. 1st winner∽——ETS + ELM vs. 2nd winner∗——ETS + ELM vs. 3rd winner∗∗——Note. ∽, nonsignificant; ∗p<0.05; ∗∗p<0.01; —, nonexistent.We also compared the computational loads of the proposed algorithm with the competing methods. In the BCI competition dataset, the training time of proposed algorithm is 24.7 s, which is far longer than the training time of MDRM (0.12 s), TSLDA (1.17 s), CSPLDA (0.18 s), and CSPSVM (0.19 s). It is because that the proposed algorithm is a semisupervised learning model where the training and test dataset are simultaneously used for learning features, but the competing methods only used training dataset to learn feature. Other reason is that the proposed algorithm takes a lot of time in iterative strategy of optimization (equation (11)), but the competing methods solve the optimization problem by quick eigenvalue decomposition. However, the test time of proposed algorithm is 0.06 s, which is close to the test time of MDRM (0.17 s), TSLDA (0.03 s), CSPLDA (0.04 s), and CSPSVM (0.0.3 s). The computational loads on in-house dataset are similar with those of competition dataset.Finally, we compared our bilinear locality-preserving embedding with other state-of-the-art manifold learning algorithms. We show the 2-dimensional embedding learning from ISOMAP, LLE, and the proposed algorithm in Figure5. The motor imagery EEG data are a spatiotemporal signal, and the correlations of EEG signals in high-dimensional space are expected to be destroyed in 2-dimensional embedding. However, compared with feature distributions extracted by other manifold learning methods, the proposed method can also obtain a moderately good separability of classes in the feature space.Figure 5
2-dimensional embedding learning from ISOMAP, LLE, and proposed algorithm for subject 1 of BCI competition IV IIa. Right hand vs. left hand show in the first row, and foot vs. tongue show in the second row.


## 4.1. Experimental Setup

### 4.1.1. Data Description

In this study, the BCI Competition 2008 motor imagery data [37] and our in-house experimental data are used for performance evaluation.(1)
The BCI Competition 2008 IIa EEG motor imagery data were recorded from nine subjects who performed four kinds of motor imagery actions (imagined right hand, left hand, foot, and tongue movements) with 22 EEG channels. For each subject, there were 144 trials per class, including training and testing trials. The EEG signals were filtered by an 8−30 Hz bandpass filter.(2)
Our in-house EEG data were recorded from 12 subjects with 64 EEG channels. For each subject, there were 234 trials per class. Our in-house experiments tested right and left hand imagined movements.The sampling rate of the EEG signal was 250 Hz. The bandpass filter was the same as the Competition dataset.


### 4.1.2. Algorithms Evaluated

We evaluate the proposed algorithms (ETS + ELM) against five competing algorithms.(1)
CSP + LDA: the CSP in [4] and linear discriminant analysis (LDA) in [38] were applied for evaluation(2)
CSP + SVM: the number of CSP spatial filters was set to 8, as suggested in [4], for feature extraction, and a support vector machine (SVM) [39] was used for classification(3)
MDRM: minimum distance to Riemannian mean [14] was used for evaluation(4)
TS + LDA : LDA was used to classify signals on a tangent space [14]


### 4.1.3. Parameter Setting

In all experiments, we empirically setk=40 and σ=1 during Riemannian graph construction. In ETS + ELM, the number of iterations was set to Nite=8, and the iterative stop threshold to τ=0.01. In fact, most information of the original space was lost during the process, while the number of embedding dimensions was too small. However, all of the methods evaluated suffer from computational cost problems under conditions of excessively high dimensionality. Hence, the dimension M was empirically set to 10. Because the ELM method suffers from an overfitting problem in the case of excessively hidden neurons, a regularizer parameter was used to control the overfitting problem. Therefore, the regularizer parameter λ was empirically set to 0.001, while the number of hidden neurons was 2000.


## 4.1.1. Data Description

In this study, the BCI Competition 2008 motor imagery data [37] and our in-house experimental data are used for performance evaluation.(1)
The BCI Competition 2008 IIa EEG motor imagery data were recorded from nine subjects who performed four kinds of motor imagery actions (imagined right hand, left hand, foot, and tongue movements) with 22 EEG channels. For each subject, there were 144 trials per class, including training and testing trials. The EEG signals were filtered by an 8−30 Hz bandpass filter.(2)
Our in-house EEG data were recorded from 12 subjects with 64 EEG channels. For each subject, there were 234 trials per class. Our in-house experiments tested right and left hand imagined movements.The sampling rate of the EEG signal was 250 Hz. The bandpass filter was the same as the Competition dataset.


## 4.1.2. Algorithms Evaluated

We evaluate the proposed algorithms (ETS + ELM) against five competing algorithms.(1)
CSP + LDA: the CSP in [4] and linear discriminant analysis (LDA) in [38] were applied for evaluation(2)
CSP + SVM: the number of CSP spatial filters was set to 8, as suggested in [4], for feature extraction, and a support vector machine (SVM) [39] was used for classification(3)
MDRM: minimum distance to Riemannian mean [14] was used for evaluation(4)
TS + LDA : LDA was used to classify signals on a tangent space [14]


## 4.1.3. Parameter Setting

In all experiments, we empirically setk=40 and σ=1 during Riemannian graph construction. In ETS + ELM, the number of iterations was set to Nite=8, and the iterative stop threshold to τ=0.01. In fact, most information of the original space was lost during the process, while the number of embedding dimensions was too small. However, all of the methods evaluated suffer from computational cost problems under conditions of excessively high dimensionality. Hence, the dimension M was empirically set to 10. Because the ELM method suffers from an overfitting problem in the case of excessively hidden neurons, a regularizer parameter was used to control the overfitting problem. Therefore, the regularizer parameter λ was empirically set to 0.001, while the number of hidden neurons was 2000.


## 4.2. Experiments for EEG Classification

To evaluate the performance of the proposed algorithms, ETS + ELM, some competing algorithms (CSP + LDA, CSP + SVM, MDRM [14], and TS + LDA [14]) were compared on the EEG datasets from the 2008 BCI Competition and our in-house BCI system.First, we evaluated the convergence performance of the iterative algorithm for optimization problem (12). Because the Kappa coefficient [40] scales linearly from 0 to 1 in the range between random and perfect classification, it was selected as a performance measure [40] for the classification algorithms in this study. In each iteration step of (12), we were able to obtain a learned mapping matrix Vt and a corresponding low-dimensional embedding. The performance of the ETS + ELM classification algorithm under the same parameter settings was evaluated on these embeddings. The average performance of nine subjects on the BCI Competition IIa dataset is shown in Figure 4. It is clear that the mapping matrix Vt learned from the proposed algorithm converged quickly. Normally, it required only two iteration steps.Figure 4
Iteration vs. mean Kappa value of ETS + ELM for BCI competition 2008 IIa (the Kappa value is the mean of nine subjects).We also compared the performance of various methods on classification of the BCI Competition 2008 data. In Table1, the performance of the proposed methods, MDRM, TS + LDA, and the top three Competition winners for dataset IIa of BCI Competition 2008 (1st of BCI08, 2nd of BCI08, and 3rd of BCI08) are presented for comparison. From the results shown in Table 1, it may be observed that the proposed ETS + ELM method achieves a mean performance of 0.577, which was the best performance of all the methods compared.Table 1
Classification results of BCI competition 2008 IIa in kappa values.
MethodsMean KappaSubjectsS01S02S03S04S05S06S07S08S09ETS + ELM0.5770.740.270.790.520.310.310.690.770.79TS + ELM0.5730.730.320.770.510.320.310.670.760.771st of BCI080.5700.680.420.750.480.400.270.770.750.61TS + LDA [14]0.5670.740.380.720.500.260.340.690.710.76MDRM [14]0.5210.750.370.660.530.290.270.560.580.682nd of BCI080.5200.690.340.710.440.160.210.660.730.693rd of BCI080.3100.380.180.480.330.070.140.290.490.44We also tested the performance of the proposed algorithms in terms of classification accuracy on the BCI Competition 2008 IIa dataset with a 30-fold cross-validation procedure. The cross-validation reflects the generalization performance of the proposed method. From the results in Table2, it is clear that, for all subjects, ETS + ELM outperformed the reference methods. This indicates that, by finding a low-dimensional embedding with the property of locality preservation, the classification method ETS + ELM exhibited the best performance.Table 2
Classification accuracy of BCI competition 2008 IIa in 30-fold cross-validation (%).
SubjectMethodsETS + ELMTS + ELMTS + LDA [14]MDRM [14]CSP + LDA [35]CSP + SVM [36]S0181.580.580.577.878.376.3S0260.358.051.344.144.750.7S0388.186.887.576.882.285.1S0464.160.059.354.959.152.9S0555.459.34543.839.748.8S0657.153.255.347.150.149.2S0782.880.782.1728178.1S0885.080.784.875.268.577.4S0987.387.486.176.677.482.2Mean ± std73.5 ± 13.271.8± 13.870.2 ± 17.163.2 ± 15.264.6 ± 16.666.7 ± 15.7Moreover, we evaluated the performance of all algorithms on the dataset recorded from our in-house BCI systems. As shown in Table3, the proposed ETS + ELM showed a significant performance improvement over the prior methods (an approximately 2% improvement in classification accuracy), even when compared with the benchmark algorithm CSP [4]. In conclusion, the experimental results on the competition dataset and in-house dataset demonstrate that the proposed ETS + ELM approached achieved a higher performance than the other methods. In additional, we provide the ANOVA and t-test analysis of Tables 1–3. The repeated measures ANOVA results for Tables 1–3 are p = 0.2 × 10−10, 0.7 × 10−7, and 0.03, respectively. It indicates that all studied methods have statistically significant difference on classification performance. The t-test analysis of Tables 1–3 are shown in Table 4.Table 3
Classification accuracy on in-house BCI dataset (%).
SubjectMethodsETS + ELMTS + ELMTS + LDA [14]MDRM [14]CSP + LDA [45]CSP + SVM [36]A0110010010010098.296.4A0298.493.797.493.598.796.1A0396.689.093.596.787.083.8A0491.085.091.284.289.482.4A0587.180.184.280.778.968.4A0690.591.273.267.882.175.0A0773.271.266.067.867.866.0A0876.371.078.971.970.173.6A0982.275.572.667.579.479.4A1094.497.789.690.993.596.1A11100100100100100100A1272.985.388.785.477.474.1Mean ± std88.4 ± 10.386.6± 10.486.2 ± 11.383.8 ± 12.685.2 ± 10.582.6 ± 11.3Table 4
Thet-test results of Tables 1–3.
Table 1Table 2Table 3ETS + ELM vs. TS + LDA∽∗∽ETS + ELM vs. MDRM∽∗∗∽ETS + ELM vs. CSP + LDA—∗∗∗ETS + ELM vs. CSP + SVM—∗∗∗ETS + ELM vs. 1st winner∽——ETS + ELM vs. 2nd winner∗——ETS + ELM vs. 3rd winner∗∗——Note. ∽, nonsignificant; ∗p<0.05; ∗∗p<0.01; —, nonexistent.We also compared the computational loads of the proposed algorithm with the competing methods. In the BCI competition dataset, the training time of proposed algorithm is 24.7 s, which is far longer than the training time of MDRM (0.12 s), TSLDA (1.17 s), CSPLDA (0.18 s), and CSPSVM (0.19 s). It is because that the proposed algorithm is a semisupervised learning model where the training and test dataset are simultaneously used for learning features, but the competing methods only used training dataset to learn feature. Other reason is that the proposed algorithm takes a lot of time in iterative strategy of optimization (equation (11)), but the competing methods solve the optimization problem by quick eigenvalue decomposition. However, the test time of proposed algorithm is 0.06 s, which is close to the test time of MDRM (0.17 s), TSLDA (0.03 s), CSPLDA (0.04 s), and CSPSVM (0.0.3 s). The computational loads on in-house dataset are similar with those of competition dataset.Finally, we compared our bilinear locality-preserving embedding with other state-of-the-art manifold learning algorithms. We show the 2-dimensional embedding learning from ISOMAP, LLE, and the proposed algorithm in Figure5. The motor imagery EEG data are a spatiotemporal signal, and the correlations of EEG signals in high-dimensional space are expected to be destroyed in 2-dimensional embedding. However, compared with feature distributions extracted by other manifold learning methods, the proposed method can also obtain a moderately good separability of classes in the feature space.Figure 5
2-dimensional embedding learning from ISOMAP, LLE, and proposed algorithm for subject 1 of BCI competition IV IIa. Right hand vs. left hand show in the first row, and foot vs. tongue show in the second row.


## 5. Conclusions

Inspired by the fact that the covariance feature of EEG signals in the form of SPD matrices lies on a Riemannian manifold, in this study, we have presented an approach to construct a Riemannian graph for motor imagery classification. A bilinear locality preservation method has been proposed to find low-dimensional embeddings carrying discriminant information. A classification algorithm, called ETS + ELM, was further proposed on a tangent space of the learned embedding. Experimental results on BCI competition and our in-house datasets show that the proposed bilinear locality-preserving embedding was able to capture the useful locality and discriminant information of the original manifold, and the proposed classification algorithms have higher performance than the other competing algorithms compared. The proposed methods can also be applied to many other pattern-recognition problems with input data in the form of SPD matrices, such as pedestrian detection problem [15] or hyperspectral image processing [41]. In future, we will further study weight calculation method to construct preferable Riemannian graph and study the interpretation of the feature learned by Riemannian-based methods.



---

*Source: 9959195-2021-09-17.xml*