A review of classification algorithms for EEG-based brain–computer interfaces: a 10 year update

The introduction of Riemannian geometry in the field of BCI has challenged some of the conventions adopted in the classic classification approaches; instead of estimating spatial filters and/or select features, the idea of a Riemannian geometry classifier (RGC) is to map the data directly onto a geometrical space equipped with a suitable metric. In such a space, data can be easily manipulated for several purposes, such as averaging, smoothing, interpolating, extrapolating and classifying. For example, in the case of EEG data, mapping entails computing some form of covariance matrix of the data. The principle of this mapping is based on the assumption that the power and the spatial distribution of EEG sources can be considered fixed for a given mental state and such information can be coded by a covariance matrix. Riemannian geometry studies smooth curved spaces that can be locally and linearly approximated. The curved space is named a manifold and its linear approximation at each point is the tangent space. In a Riemannian manifold the tangent space is equipped with an inner product (metric) smoothly varying from point to point. This results in a non-Euclidean notion of distance between any two points (e.g. each point may be a trial) and a consequent notion of centre of mass of any number of points (figure 2). Therefore, instead of using the Euclidean distance, called the extrinsic distance, an intrinsic distance is used, which is adapted to the geometry of the manifold, and thus to the manner in which the data have been mapped [47, 232].

Zoom In Zoom Out Reset image size

Amongst the most common matrix manifolds used for BCI applications, we encountered the manifold of Hermitian or symmetric positive definite (SPD) matrices [19] when dealing with covariance matrices estimated from EEG trials, and the Stiefel and Grassmann manifolds [62] when dealing with subspaces or orthogonal matrices. Several machine learning problems can be readily extended to those manifolds by taking advantage of their geometrical constraints (i.e. learning on manifold). Furthermore, optimization problems can be formulated specifically on such spaces, which is leading to several new optimization methods and to the solution of new problems [2]. Although related, manifold learning, which consists of empirically attempting to locate the non-linear subspace in which a dataset is defined, is different in concept and will not be covered in this paper. To illustrate these notions, consider the case of SPD matrices. The square of the intrinsic distance between two SPD matrices $\newcommand{\C}{{\bf C}} \C_1$ and $\newcommand{\C}{{\bf C}} \C_2$ has a closed-form expression given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\C}{{\bf C}} \displaystyle \delta^2 \left ( \C_1,\C_2 \right)=\sum_n \log^2\lambda_n\left ( \C_1^{-1}\C_2 \right), \label{eq:Rdist} \nonumber \end{align} \tag{ 1 }$$

where $\newcommand{\M}{{\bf M}} \lambda_n(\M)$ denotes the nth eigenvalue of matrix $\newcommand{\M}{{\bf M}} \M$. For $\newcommand{\C}{{\bf C}} \C_1$ and $\newcommand{\C}{{\bf C}} \C_2$ SPDs, this distance is non-negative, symmetric and is equal to zero if and only if $\newcommand{\C}{{\bf C}} \C_1=\C_2$. Interestingly, when $\newcommand{\C}{{\bf C}} \C_1$ and $\newcommand{\C}{{\bf C}} \C_2$ are the means of two classes, the eigenvectors of matrix $\newcommand{\C}{{\bf C}} (\C_1^{-1}\C_2)$ are used to define CSP filters, while its eigenvalues are used for computing their Riemannian distance [47]. Using the distance in equation (1), the centre of mass $\newcommand{\G}{{\bf G}} \G$ of a set $\newcommand{\C}{{\bf C}} \left \{\C_1, ..., \C_K \right \}$ of K SPD matrices (figure 3), also called the geometric mean, is the unique solution to the following optimization problem

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\G}{{\bf G}} \newcommand{\C}{{\bf C}} \displaystyle \arg\!\min\limits_{\G} \sum_k \delta^2 ( \C_k,\G). \label{eq:Rmean} \nonumber \end{align} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

As discussed thoroughly in [47], this definition is analogous to the definition of the arithmetic mean $\newcommand{\C}{{\bf C}} 1/K\sum_k \C_k$, which is the solution of the optimization problem (2) when the Euclidean distance is used instead of the Riemannian one. In contrast to the arithmetic mean, the geometric mean does not have a closed-form solution. A fast and robust iterative algorithm for computing the geometric mean has been presented in [48]. The simplest RGC methods allow immediate classification of trials (mapped via some form of covariance matrix) by simple nearest neighbour methods, using exclusively the notion of Riemannian distance (equation (1)), and possibly with the notion of geometric mean (2). For instance, the Riemannian minimum distance to mean (RMDM) classifier [15, 13] computes a geometric mean for each class using training data and then assigns an unlabelled trial to the class corresponding to the closest mean (figure 3). Another class of RGCs consists of methods projecting the data points to a tangent space followed by a classification, thereafter using standard classifiers such as LDA, SVM, logistic regression, etc [13, 14]. These methods take advantage of both the Riemannian geometry and the possibility of executing complex decision functions using dedicated classifiers. An alternative approach is to project the data in the tangent space, filter the data there (for example by LDA), and map the data back onto the manifold to finally carry out the RMDM.

As described above, Riemannian classifiers either operate directly on the manifold (e.g. the RMDM) or by the projection of the data in the tangent space. Simple RGCs on the manifold have been shown to be competitive as compared to previous state-of-the-art classifiers used in BCI as long as the number of electrodes is not very large, providing better robustness to noise and better generalization capabilities, both on healthy users [13, 46, 100] and clinical populations [158]. RGCs based on tangent space projection clearly outperformed the other state-of-the-art methods in terms of accuracy [13, 14], as demonstrated by the first place they have been awarded in five recent international BCI predictive modelling data competitions, as reported in [47]. For a comprehensive review of the Riemannian approaches in BCI, the reader can refer to [47, 232]. The various approaches using Riemannian geometry classifiers for EEG-based BCIs are summarized in table 3.

As highlighted in [232], the processing procedures of Riemannian approaches such as RMDM is simpler and involves fewer stages than more classic approaches. Also, Riemannian classifiers apply equally well to all BCI paradigms (e.g. BCIs based on mental imagery, ERPs and SSVEP); only the manner in which data points are mapped in the SPD manifold differs (see [47] for details). Furthermore, in contrast to most classification methods, the RMDM approach is parameter-free, that is, it does not require any parameter tuning, for example by cross-validation. Hence, Riemannian geometry provides new tools for building simple, more robust and accurate prediction models.

Several reasons have been proposed to advocate the use of the Riemannian geometry. Due to its logarithmic nature the Riemannian distance is robust to extreme values, that is, noise. Also, the intrinsic Riemannian distance for SPD matrices is invariant both to matrix inversion and to any linear invertible transformation of the data, e.g. any mixing applied to the EEG sources does not change the distances among the observed covariance matrices. These properties in part explain why Riemannian classification methods provide a good generalization capability [224, 238], which enabled researchers to set up calibration-free adaptive ERP-BCIs using simple subject-to-subject and session-to-session transfer learning strategies [6].

Interestingly, as illustrated in [94], it is possible to not only interpolate along geodesics (figure 2) on the SPD manifolds, but also to extrapolate (e.g. forecast) without leaving the manifold and respecting the geometrical constraints. For example, in [99] interpolation has been used for data augmentation by generating artificial covariance matrices along geodesics but extrapolation could also have been used. Often, the Riemannian interpolation is more relevant than its Euclidean counterpart as it does not suffer from the so-called swelling effect [232]. This effect describes the fact that a Euclidean interpolation between two SPD matrices does not involve the determinant of the matrix as it should (i.e. the determinant of the Euclidean interpolation can exceed the determinant of the interpolated matrices). In the spirit of [231], the determinant of a covariance matrix can be considered as the volume of the polytope described by the column of the matrix. Thus, a distance that is immune to the swelling effect will respect the shape of the polytope along geodesics.

As equation (1) indicates, computing the Riemannian distance between two SPD matrices involves adding squared logarithms, which may cause numerical problems; the smallest eigenvalues of matrix $\newcommand{\C}{{\bf C}} (\C_1^{-1}\C_2)$ tend towards zero as the number of electrodes increases and/or the window size for estimating $\newcommand{\C}{{\bf C}} \C_1$ and $\newcommand{\C}{{\bf C}} \C_2$ decreases, making the logarithm operation ill-conditioned and numerically unstable. Further, note that the larger the dimensions, the more the distance is prone to noise. Moreover, Riemannian approaches usually have high computational complexities (e.g. growing cubically with the number of electrodes for computing both the geometric mean and the Riemannian distance). For these reasons, when the number of electrodes is large with respect to the window size, it is advocated to reduce the dimensions of the input matrices. Classical unsupervised methods such as PCA or supervised methods such as CSP can be used for this purpose. Recently, Riemannian-inspired dimensionality reduction methods have been investigated as well [94, 95, 189].

Interestingly, some approaches have tried to bridge the gap between Riemannian approaches and more classical paradigms by incorporating some Riemannian geometry in approaches such as CSP [12, 233]. CSP was the previous golden standard and is based on a different paradigm than Riemannian geometry. Taking the best of those two paradigms is expected to gain better robustness while compressing the information.

As mentioned previously, the classification pipeline in BCI typically involves spatial filtering of the EEG signals followed by classification of the filtered data. This results in the independent optimization of several sets of parameters, namely for the spatial filters and for the final classifier. For instance, the typical linear classifier decision function for an oscillatory activity BCI would be the following:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\w}{{\bf w}} \newcommand{\X}{{\bf X}} \displaystyle f(\X,\w,{\bf S}) = \sum_i w_i {\rm log}({\rm var}({\bf s}_i^T \X)) + w_0 \label{equ:stdDecisionClassif} \nonumber \end{align} \tag{ 3 }$$

where $\newcommand{\X}{{\bf X}} \X$ is the EEG signals matrix, $\newcommand{\w}{{\bf w}} \w = [w_0, w_1, \ldots, w_N]$ is the linear classifier weight vector, and ${\bf S} = [{\bf s}_1, {\bf s}_2, \ldots, {\bf s}_N]$ is a matrix of spatial filter ${\bf s}_i$. Optimizing w and s_i separately may thus lead to suboptimal solutions, as the spatial filters do not consider the objective function of the classifier. Therefore, in addition to RGC, several authors have shown that it is possible to formulate this dual optimization problem as a single one, where the parameters of the spatial filters and the linear classifier are optimized simultaneously, with the potential to obtain improved performance. The key principle of these approaches is to learn classifiers (either linear vector classifiers or matrix classifiers) that directly use covariance matrices as input, or their vectorised version. We briefly present these approaches below.

In [214], the EEG data were represented as an augmented covariance matrix $\newcommand{\A}{{\bf A}} \A$, containing as block diagonal terms both the first order term $\newcommand{\X}{{\bf X}} \X$, i.e. the signal time course, and as second order terms the covariance matrices of EEG trials band-pass filtered in various frequency bands. The learned classifier is thus a matrix of weights $\newcommand{\W}{{\bf W}} \W$ (rather than a vector), with the decision function $\newcommand{\W}{{\bf W}} \newcommand{\A}{{\bf A}} f(\A, \W) = \langle \A, \W\rangle + b$. Due to the large dimensionality of the augmented covariance matrix, a matrix regularization term is necessary with such classifiers, e.g. to obtain sparse temporal or spatial weights. Note that this approach can be applied to both ERP and oscillatory-based BCI, as the first order terms capture the temporal variation, and the covariance matrices capture the EEG signals' band power variations.

Following similar ideas in parallel, [65] represented this learning problem in tensor space by constructing tensors of frequency-band specific covariance matrices, which can then be classified using a linear classifier as well, provided appropriate regularization is used.

Finally, [190] demonstrated that equation (3) can be rewritten as follows, if we drop the log-transform:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\w}{{\bf w}} \displaystyle f(\Sigma, \w_{\Sigma}) = {\rm vec}({\boldsymbol{\Sigma}})^T \w_{\Sigma} + w_0 \nonumber \end{align} \tag{ 4 }$$

with $\newcommand{\w}{{\bf w}} \w_{\Sigma} = \sum_i w_i {\rm vec}({\bf s}_i {\bf s}_i^T)$, $\newcommand{\X}{{\bf X}} {\boldsymbol{\Sigma}} = \X^T \X$ being the EEG covariance matrix, and $\newcommand{\M}{{\bf M}} {\rm vec}(\M)$ being the vectorisation of matrix $\newcommand{\M}{{\bf M}} \M$. Thus, equation (3) can be optimized directly in the space of vectorised covariance matrices by optimizing the weights $\newcommand{\w}{{\bf w}} \w_{\Sigma}$. Here as well, owing to the usually large dimensionality of ${\rm vec}(\sigma)$, appropriate regularization is necessary, and [190] explored different approaches to do so.

These different approaches all demonstrated higher performance than the basic CSP+LDA methods on motor imagery data sets [65, 190, 214]. This suggests that such formulations can be worthy alternatives to the standard CSP+LDA pipelines.

By simultaneously optimizing spatial filters and classifiers, such formulations usually achieve better solutions than the independent optimization of individual sets of components. Their main advantage is thus increased classification performance. This formulation nonetheless comes at the expense of a larger number of classifier weights due to the high increase in dimensionality of the input features (covariance matrix with $(N_c*(N_c+1))/2$ unique values versus N_c values when using only the channels' band power). Appropriate regularization is thus necessary. It remains to be evaluated how such methods perform for various amounts of training data, as they are bound to suffer more severely from the curse of dimensionality than simpler methods with fewer parameters. These methods have also not been used online to date. From a computational complexity point of view, such methods are more demanding than traditionnal methods given their increased number of parameters, as mentioned above. They also generally require heavy regularization, which can make their calibration longer. However, their decision functions being linear, they should be easily applicable in online scenarios. However, it remains to be seen whether they can be calibrated quickly enough for online use, and what their performance will be for online data.

Tensors (i.e. multi-way arrays) provide a natural representation for EEG data, and higher order tensor decompositions and factorizations are emerging as promising (but not yet very well established and not yet fully explored) tools for analysis of EEG data; particularly for feature extraction, clustering and classification tasks in BCI [38–40, 42, 43].

The concept of tensorization refers to the generation of higher-order structured tensors (multiway arrays) from lower-order data formats, especially time series EEG data represented as vectors or organized as matrices. This is an essential step prior to tensor (multiway) feature extraction and classification [41, 42, 182].

The order of a tensor is the number of modes, also known as ways or dimensions (e.g. for EEG BCI data: space (channels), time, frequency, subjects, trials, groups, conditions, wavelets, dictionaries). In the simplest scenario, multichannel EEG signals can be represented as a 3rd-order tensor that has three physical modes: space (channel)  ×  time  ×  frequency. In other words, S channels of EEG which are recorded over T time samples, can produce S matrices of $F \times T$ dimensional time-frequency spectrograms stacked together into an $F \times T \times S$ dimensional third-order tensor. For multiple trials and multiple subjects, the EEG data sets can be naturally represented by higher-order tensors: e.g. for a 5th-order tensor: space  ×  time  ×  frequency  ×  trial  ×  subject.

It should be noted that almost all basic vector- and matrix-based machine learning algorithms for feature extraction and classification have been or can be extended or generalized to tensors. For example, the SVM for classification has been naturally generalized to the tensor support machine (TSM), Kernel TSM and higher rank TSM. Furthermore, the standard LDA method has been generalized to tensor Fisher discriminant analysis (TFDA) and/or higher order discriminant analysis (HODA) [41, 183]. Moreover tensor representations of BCI data are often very useful in mitigating the small sample size problem in discriminative subspace selection, because the information about the structure of data is often inherent in tensors and is a natural constraint which helps reduce the number of unknown feature parameters in the description of a learning model. In other words, when the number of EEG training measurements is limited, tensor-based learning machines are expected often to perform better than the corresponding vector- or matrix-based learning machines, as vector representations are associated with problems such as loss of information for structured data and over-fitting for high-dimensional data.

To ensure that the reduced data sets contain maximum information about input EEG data, we may apply constrained tensor decomposition methods. For example, this could be achieved on the basis of orthogonal or non-negative tensor (multi-array) decompositions, or higher order (multilinear) discriminant analysis (HODA), whereby input data are considered as tensors instead of more conventional vector or matrix representations. In fact, tensor decomposition models, especially PARAFAC (also called CP decomposition), TUCKER, hierarchical tucker (HT) and tensor train (TT) are alternative sophisticated tools for feature extraction problems by capturing multi-linear and multi-aspect structures in large-scale higher-order data-sets [39, 183]. Using this type of approach, we first decompose multi-way data using TUCKER or CP decompositions, usually by imposing specific constraints (smoothness, sparseness, non-negativity), in order to retrieve basis factors and significant features from factor (component) matrices. For example, wavelets/dictionaries allow us to represent the data often in a more efficient way, i.e. a sparse manner with different sparsity profiles [43, 183].

Moreover, in order to increase performance of BCI classification, we can apply two or more time-frequency representations or the same frequency transform but with two or more different parameter settings. Different frequency transforms (or different mother wavelets) allow us to obtain different sparse tensor representations with various sparsity profiles and some complimentary information. For multichannel EEG signals we can generate a block of at least two tensors, which can be concatenated as a single data tensor: space  ×  time  ×  frequency  ×  trial [43, 182, 183].

The key problem in tensor representation is the choice of a suitable time-frequency representation (TFR) or frequency transform and the selection of optimal, or close to optimal, corresponding transformation parameters. By exploiting various TFRs, possibly with suitably selected different parameter settings for the same data, we may potentially improve the classification accuracy of BCI due to additional (partially redundant) information. Such approaches have been implemented e.g. for motor imagery (MI) BCI by employing different complex Morlet (Gabor) wavelets for EEG data sets with 62 channels [183]. For such data sets, the authors selected different complex Morlet wavelets with two different bandwidth frequency parameters fb  =  1 Hz and fb  =  6 Hz for the same centre frequency fc  =  1 Hz. For each mother wavelet the authors constructed a 4th-order tensor: 62-channels  ×  23-frequency bins  ×  50-time frames  ×  120-trials for both training and test EEG data. The block of training tensor data can be concatenated as the 5th-order tensor: 62-channels  ×  23-frequency bins  ×  50-time frames  ×  2-wavelets  ×  120-trials.

The HODA algorithm was used to estimate discriminant bases. The four most significant features were selected to classify the data, and led to an improved accuracy higher than 95%. Thus, it appears that by applying tensor decomposition for suitably constructed data tensors, considerable performance improvement in comparison to the standard approaches can be achieved for both motor-imagery BCI [183, 223] and P300 paradigms [175].

In this approach, transformation of data with a dictionary aims to un-correlate the raw data and express them in a sparse domain. Different dictionaries (transformations) contribute to obtaining different sparse representations with various sparsity profiles. Moreover, augmentation of dimensionality to create samples with additional modes improved the performance.

To summarize, tensor decompositions with nonnegative, orthonormal or discriminant bases improved the classification accuracy for the BCI dataset by almost 10%. A comparison of all methods mentioned is provided in table 4.

From the time-frequency analysis perspective, tensor decompositions are very attractive, even for a single channel, because they simultaneously take into account temporal and spectral information and variability and/or consistency of time frequency representations (TFRs) for trials and/or subjects. Furthermore, they provide links among various latent (hidden) variables (e.g. temporal, spectral and spatial components) often with physical or physiological meanings and interpretations [40, 183].

Furthermore, standard canonical correlation analysis (CCA) was generalized to tensor CCA and multiset CCA and was successfully applied to the classification of SSVEP for BCI [242, 243, 245, 246]. Tensor canonical correlation analysis (TCCA) and its modification multiset canonical correlation analysis (MsetCCA) have been one of the most efficient methods for frequency recognition in SSVEP-BCIs. The MsetCCA method learns multiple linear transforms that implement joint spatial filtering to maximize the overall correlation amongst canonical variates, and hence extracts SSVEP common features from multiple sets of EEG data recorded at the same stimulus frequency. The optimized reference signals are formed by combination of the common features and are completely based on training data. Extensive experimental study with EEG data demonstrated that the tensor and MsetCCA method improve the recognition accuracy of SSVEP frequency in comparison with the standard CCA method and other existing methods, especially for a small number of channels and a short time window length. The superior results indicate that the tensor MsetCCA method is a very promising candidate for frequency recognition in SSVEP-based BCIs [243].

In summary, the recent advances in BCI technologies have generated massive amounts of brain data exhibiting high dimensionality, multiple modality (e.g. physical modes such as frequency or time, multiple brain imaging techniques or conditions), and multiple couplings as functional connectivity data. By virtue of their multi-way nature, tensors provide powerful and promising tools for BCI analysis and fusion of massive data combined with a mathematical backbone for the discovery of underlying hidden complex (space-time-frequency) data structures [42, 183].

Another of their advantages is that, using tensorization and low-rank tensor decomposition, they can efficiently compress large multidimensional data into low-order factor matrices and/or core tensors which usually represent reduced features. Tensor methods can also analyze linked (coupled) blocks of trials represented as large-scale matrices into the form of tensors in order to separate common/correlated from independent/uncorrelated components in the observed raw EEG data.

Finally, it is worth mentioning that tensor decompositions are emerging techniques not only for feature extraction/selection and BCI classification, but also for pattern recognition, multiway clustering, sparse representation, data fusion, dimensionality reduction, coding, and multilinear blind brain source separation (MBSS). They can potentially provide convenient multi-channel and multi-subject space-time-frequency sparse representations, artefact rejection, feature extraction, multi-way clustering and coherence tracking [39, 40].

On the cons side, the complexity of tensor methods is usually much higher than standard matrix and vector machine learning methods. Moreover, since tensor methods are just emerging as potential tools for feature extraction and classification, existing algorithms are not always mature and are still not fully optimized. Thus, some efforts are still needed to optimize and test them for real-life large scale data sets.

A restricted Boltzmann machine (RBM) is a Markov random field (MRF) [120] associated with a bipartite undirected graph. It is composed of two sets of units: m visible ones $V = (V_1, \cdots, V_m)$ and n hidden ones $H = (H_1, \cdots, H_n)$. The visible units are used for representing observable data whereas the hidden ones capture some dependencies between observed variables. For the usual type of RBM such as those discussed in this paper, units are considered as random variables that take binary values $({\bf v}, {\bf h})$ and $\newcommand{\W}{{\bf W}} \W$ is a matrix whose entries w_i,j are the weights associated with the connection between unit v_i and h_j. The joint probability of a given configuration $({\bf v}, {\bf h})$ can be modelled according to the probability $p({\bf v}, {\bf h})= \frac{1}{Z} {\rm e}^{-E({\bf v}, {\bf h})}$ with the energy function $E({\bf v}, {\bf h})$ being

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\W}{{\bf W}} \displaystyle E({\bf v},{\bf h}) = {\bf a}^\top {\bf v} - b^\top {\bf h} - {\bf v}^\top \W {\bf h} \nonumber \end{align*}$$

where ${\bf a}$ and ${\bf b}$ are bias weight vectors. Note Z is a normalizing factor in order that $p({\bf v}, {\bf h})$ sums to one for all possible configurations. Owing to the undirected bipartite graph property, hidden (respective to the visible) variables are independent given the visible (hidden) ones leading to:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle p({\bf v} | {\bf h}) = \Pi_{i=1}^m p(v_i | {\bf h}) \quad \quad p({\bf h} | {\bf v}) = \Pi_{j=1}^n p(h_j | {\bf v}) \nonumber \end{align*}$$

and marginal distributions over the visible variables can be easily obtained as [69]:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle p({\bf v})= \frac{1}{Z} \sum_h {\rm e}^{-E({\bf v},{\bf h})}. \nonumber \end{align*}$$

Hence, by optimizing all model parameters $\newcommand{\W}{{\bf W}} (\W, {\bf b}, {\bf a})$, it is possible to model the probability distribution of the observable variables. Other properties of RBMS as well as connections of RBMs with stochastic neural networks are detailed in [69, 90]

To learn the probability distribution of the input data, RBMs are usually trained according to a procedure denoted as contrastive divergence learning [89]. This learning procedure is based on a gradient ascent of the log-likelihood of the training data. The derivative of the log-likelihood of an input v can be easily derived [69] and the mean of this derivative over the training set leads to the rule:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\W}{{\bf W}} \displaystyle \sum_v \frac{\partial L(\W|{\bf v})}{\partial w_{i,j}} \propto \langle v_i h_j \rangle_{\rm data} - \langle v_i h_j \rangle_{\rm model} \nonumber \end{align*}$$

with the two brackets respectively denoting expectation over $p({\bf h}|{\bf v})q({\bf v})$ and over the model ($p({\bf v}, {\bf h})$) with q being the empirical distribution of the inputs. While the first term of this gradient is tractable, the second one has exponential complexity. Contrastive divergence aims at approximating this gradient using a Gibbs chain procedure that computes the binary state of ${\bf h}$ using $p({\bf h}|{\bf v})$ and then obtaining an estimation of ${\bf v}$ using $p({\bf v}|{\bf h})$ [89]. There exist other methods for approximating the gradient of RBMs log-likelihood that may lead to better solutions as well as methods for learning with continuous variables [17, 213].

The above procedure allows one to learn a generative model of the inputs using a simple layer of RBMs. A deep learning strategy can be obtained by stacking several RBMs with the hidden units of one layer used as inputs of the subsequent layers. Each layer is usually trained in a greedy fashion [90] and fine-tuning can be performed depending on the final objective of the model. A RBM is illustrated in figure 5 (left).

Zoom In Zoom Out Reset image size

A convolutional neural network (ConvNet or CNN) is a feedforward neural network (a network in which information flows uni-directionally from the input to the hidden layers to the output) which has at least one convolutional layer [71, 117, 117]. Such a convolutional layer maps its input to an output through a convolution operator. Suppose that the input is a 1D signal {x_n} with N samples, its convolution through a 1D filter {h_m} of size M is given by:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle y(n) = \sum_{i=0}^{M-1} h_i x_{n-i} \quad \forall n = 0, \cdots, N -1. \nonumber \end{align*}$$

This equation can be extended to higher dimensions by augmenting the number of summations in accordance with the dimensions. Several filters can also be independently used in convolution operations leading to an increased number of channels in the output. This convolutional layer is usually followed by nonlinearities [75] and possibly by a pooling layer that aggregate the local information of the output into a single value, typically through an average or a max operator [25]. Standard ConvNet architectures usually stack several of these layers (convolution  +  non-linearity (+ pooling)) followed by other layers, typically fully connected, that act as a classification layer. Note however that some architectures use all convolutional layers as classification layers. Given some architectures, the parameters of the models are the weights of all the filters used for convolution and the weights of the fully connected layers.

ConvNets are usually trained in a supervised fashion by solving an empirical risk minimization problem of the form:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\w}{{\bf w}} \newcommand{\x}{{\bf x}} \displaystyle \hat \w = \arg\min_w \frac{1}{\ell}\sum_i L(y_i, f_\w(\x_i)) + \Omega(\w) \nonumber \end{align*}$$

where $\newcommand{\x}{{\bf x}} \newcommand{\e}{{\rm e}} \{\x_i, y_i\}_{i=1}^\ell$ are the training data, $\newcommand{\w}{{\bf w}} f_\w$ is the prediction function related to the ConvNet, $L(\cdot, \cdot)$ is a loss function that measures any discrepancy between the true labels of $\newcommand{\x}{{\bf x}} \x_i$ and $\newcommand{\x}{{\bf x}} \newcommand{\w}{{\bf w}} f_\w(\x)_i$, and Ω is a regularization function for the parameters of the ConvNet. Owing to the specific form of the global loss (average loss over the individual samples), stochastic gradient descent and its variants are the most popular means for optimizing deep ConvNets. Furthermore, the feedforward architecture of $\newcommand{\w}{{\bf w}} f_\w(\cdot)$ allows the computation of the gradient at any given layer using the chain rule. This can be performed efficiently using the back-propagation algorithm [193].

In several domain applications, ConvNets have been very successful because they are able to learn the most relevant features for the task at hand. However, their performances strongly depend on their architectures and their learning hyper-parameters. A ConvNet is illustrated in figure 5 (right).

In order to classify more than two mental tasks, two main approaches can be used to obtain a multiclass classification function [215]. The first approach consists in directly estimating the class using multiclass techniques such as decision trees, multilayer perceptrons, naive Bayes classifiers or k-nearest neighbours. The second approach consists of decomposing the problem into several binary classification problems [5]. This decomposition can be accomplished in different ways using (i) one-against-one pairwise classifiers [20, 84], (ii) one-against-the-rest (or one-against-all) classifiers [20, 84], (iii) hierarchical classifiers similar to a binary decision tree and (iv) multi-label classifiers [154, 215]. In the latter case, a distinct subset of L labels (or properties) is associated to each class [58]. The predicted class is identified according to the closest distance between the predicted labels and each subset of labels defining a class.

The number of commands provided by motor imagery-based BCIs depends on the number of mental imagery states that the system is able to detect. This, in turn, is limited by the number of body parts that users can imagine moving in a manner that generate clear and distinct EEG patterns. Multi-label approaches can thus prove useful for detecting combined motor imagery tasks, i.e. imagination of two or more body parts at the same time [125, 192, 226], with each body part corresponding to a single label (indicating whether that body part was used). Indeed, in comparison with the standard approach, this approach has the advantage of considerably increasing the number of different mental states while using the same number of body parts: 2^P compared to P, where P is the number of body parts. Thus, EEG patterns during simple and combined motor imagery tasks were investigated to confirm the separability of seven different classes of motor imagery for BCI [126, 226, 249]. For the purpose of achieving continuous 3D control, both hands motor imagery was adopted to complement the set of instructions in a simple limb motor imagery based-BCI to go up (and rest to go down) [114, 192]. The up/down control signal was the inverted addition of left and right autoregressive spectral amplitudes calculated for each of the electrodes and 3 Hz-frequency bins. Another method converted circular ordinal regression to a multi-label classification approach to control a simulated wheelchair, using data set IIIa of the third BCI competition, with as motor tasks imagination of left hand, right hand, foot and tongue movements [57]. Multiclass and multi-label approaches have been compared to discriminate height commands from the combination of three motor imagery tasks (left hand, right hand and feet) to control a robotic arm [125]. A first method used a single classifier applied to the concatenated features related to each activity source (C3, Cz, C4), with one source for each limb involved. A second approach consisted of a hierarchical tree of three binary classifiers to infer the final decision. The third approach was a combination of the first two approaches. All methods used the CSP algorithm for feature extraction and Linear discriminant analysis (LDA) for classification. All methods were validated and compared to the classical one-versus-one (OVO) and one-versus-rest (OVR) methods. Results obtained with the hierarchical method were similar to the ones obtained with the OVO and OVR approaches. The performance obtained with the first approach (single classifier) and the last (combined hierarchical classifier) were the best for all subjects. The various multi-label approaches explored are mentioned in table 7.

Multiclass and multi-label approaches therefore aim to recognize more than two commands. In both cases, the resulting increase in the number of recognized classes potentially provides the user with a greater number of commands to interact more quickly with the system, without the need for a drop-down menu, for example. The multi-label approach can make learning shorter and less tiring, as it requires learning only a small number of labels. The many possible combinations of these labels leads to a large number of classes and therefore to more commands. In addition, the multi-label approach allows redundancy in the labels describing a class, which can lead to better class separation. Usually, the number of labels to produce is less than the number of classes. Finally, as compared to standard methods, multiclass and multilabel approaches usually have a lower computational complexity since they can share parameters, e.g. using multilayer perceptron, or class descriptors (especially if no redundancy is introduced).

However, there might be a lack of relationship between the meaning of a label and the corresponding mental command, e.g. two hand imagery to go up. This may generate a greater mental workload and therefore fatigue. It is therefore necessary to choose carefully the mapping between mental commands and corresponding labels. Finally, classification errors remain of course possible. In particular, the set of estimated labels may sometimes not correspond to any class, and several classes may be at equal distances, thus causing class confusion.

As previously discussed, most EEG-based BCIs are currently optimized for each subject. Indeed, this has been shown to lead in general to substantially higher classification performances than subject-independent classifiers. Typical BCI systems can be optimized by using only a few training data, typically 20–100 trials per class, as subjects cannot be asked to produce the same mental commands thousands of times before being provided with a functional BCI. Moreover, collecting such training data takes time, which is inconvenient for the subjects, and an ideal BCI would thus require a calibration time as short as possible. This calls for classifiers that can be calibrated using as little training data as possible. In the following we present those classifiers that were shown to be effective for this purpose. They rely on using statistical estimators dedicated to small sample size or on dividing the input features between several classifiers to reduce the dimensionality, thus reducing the amount of training data needed by each classifier.

The three main classifiers that have been shown to be effective with little training data, and thus effective for EEG-based BCI design, are the shrinkage LDA classifier [22, 137, 142], random forest [3, 54] and Riemannian classifiers [47, 232].

The shrinkage LDA (sLDA), is a standard LDA classifier in which the class-related covariance matrices used in its optimization were regularized using shrinkage [22]. Indeed, covariance matrices estimated from little data tend to have larger extreme eigenvalues than the real data distribution, leading to poor covariance estimates. This can be solved by shrinking covariance matrices ${\boldsymbol{\Sigma}}$ as $\newcommand{\I}{{\bf I}} \hat{{\boldsymbol{\Sigma}}} = {\boldsymbol{\Sigma}} - \lambda \I$, with $\newcommand{\I}{{\bf I}} \I$ the identity matrix, and λ the regularization parameter. Interestingly enough, there are analytical solutions to automatically determine the best λ value (see [118]. The resulting sLDA classifier has been shown to be superior to the standard LDA classifier for BCI, both for ERP-based BCI [22] and for oscillatory activity BCI [137]. It has also been shown that such classifier can be calibrated with much fewer data than LDA to achieve the same performance [137, 142]. For instance, for mental imagery BCI, an sLDA has been shown to obtain similar performance with ten training trials per class than a standard LDA with 30 training trials per class, effectively reducing the calibration time three-fold [137].

Random forest (RF) classifiers are ensembles of several decision tree classifiers [26]. The idea behind this classifier is to randomly select a subset of the available features, and to train a decision tree classifier on them, then to repeat the process with many random feature subsets to generate many decision trees, hence the name random forest. The final decision is taken by combining the outputs of all decision trees. Because each tree only uses a subset of the features, it is less sensitive to the curse-of-dimensionality, and thus requires fewer training data to be effective. Outside of BCI research, among various classifiers and across various classification problems and domains, random forest algorithms were actually often found to be among the most accurate classifiers, including problems with small training data sets [26, 36]. RFs were used successfully even online both for ERP-based BCI [3] and for motor imagery BCI [54]. They outperformed designs based on LDA classifiers for motor imagery BCI [54].

Riemannian classifiers have been discussed in section 4.2.1. Typically, a simple Riemannian classifier such as the RMDM requires less training data as compared to optimal filtering approaches such as the CSP for motor imagery [46] and xDAWN for P300 [15]. This is due to the robustness of the Riemannian distance, which the geometric mean inherits directly, as discussed in [47]. Even more robust mean estimations can be obtained computing Riemannian medians or trimmed Riemannian means. Shrinkage and other regularization strategies can also be applied to a Riemannian framework to improve the estimation of covariance matrices when a small number of data points is considered [100]. These methods are summarized in table 8.

sLDA, RF and the RMDM are simple classifiers that are easy to use in practice and provide good results in general, including online. We thus recommend their use. sLDA and RMDM do not have any hyper-parameters, which makes them very convenient to use. sLDA have been shown to be superior to LDA both for ERP and oscillatory activity-based BCI across a number of data sets [22, 137]. There is thus no reason to use classic LDA; instead sLDA should be preferred. RMDM performs as well as CSP+LDA for oscillatory activity-based BCI [13, 46], as well as xDAWN+LDA but better than a step-wise LDA on time samples for ERP-based BCI [15, 46] and better than CCA for SSVEP [100]. Note that because LDA is a linear classifier, it may be suboptimal in the hypothetic future case when vast amounts of training data will be available. RF on the other hand is a non-linear classifier that can be effective both with small and large training sets [36]. RMDM is also non-linear and performs well with small as well as large training sets [46]. In terms of computational complexity, while RF can be more demanding than RMDM or sLDA since it uses many classifiers, all of them are fairly simple and fast methods, and all have been used online successfully on standard computers.