HAL Id: hal-01394253
https://hal.inria.fr/hal-01394253
Submitted on 9 Nov 2016
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diusion de documents
scientiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Riemannian approaches in Brain-Computer Interfaces: a
review
Florian Yger, Maxime Berar, Fabien Lotte
To cite this version:
Florian Yger, Maxime Berar, Fabien Lotte. Riemannian approaches in Brain-Computer Interfaces:
a review. IEEE Transactions on Neural Systems and Rehabilitation Engineering, IEEE Institute of
Electrical and Electronics Engineers, 2017. �hal-01394253�
SPECIAL ISSUE ON BMI/BCI SYSTEMS IN IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING 1
Riemannian approaches in Brain-Computer
Interfaces: a review
Florian Yger, Maxime Berar and Fabien Lotte
Abstract—Although promising from numerous applications,
current Brain-Computer Interfaces (BCIs) still suffer from a
number of limitations. In particular, they are sensitive to noise,
outliers and the non-stationarity of ElectroEncephaloGraphic
(EEG) signals, they require long calibration times and are not
reliable. Thus, new approaches and tools, notably at the EEG
signal processing and classification level, are necessary to address
these limitations. Riemannian approaches, spearheaded by the
use of covariance matrices, are such a very promising tool slowly
adopted by a growing number of researchers. This article, after a
quick introduction to Riemannian geometry and a presentation
of the BCI-relevant manifolds, reviews how these approaches
have been used for EEG-based BCI, in particular for feature
representation and learning, classifier design and calibration time
reduction. Finally, relevant challenges and promising research
directions for EEG signal classification in BCIs are identified,
such as feature tracking on manifold or multi-task learning.
Index Terms—Riemannian geometry, Brain-Computer Inter-
face (BCI), covariance matrices, subspaces, source extraction,
Electroencephalography (EEG), classification
I. INTRODUCTION
B
RAIN-COMPUTER INTERFACES (BCIs) enable their
users to interact with computers via brain activity only,
this activity being typically measured by ElectroEncephaloG-
raphy (EEG) [1]. For instance, a BCI can enable a user to move
a cursor leftwards or rightwards on a computer screen, by
imagining left or right hand movements respectively [2]. BCIs
have proven very promising, e.g., to provide communication
to severely paralyzed users [3], as a new control device for
gaming [4] or to design adaptive human-computer interfaces
that can react to the user’s mental state [5], to name a few
[6]. However, most of these applications are prototypes and
current BCI are still scarcely used outside laboratories.
The main reason that prevents EEG-based BCIs from being
widely used is their low robustness and reliability [1][7].
Indeed, current BCIs too often recognize erroneous commands
from the user, which results in rather low accuracy and
information transfer rate, even in laboratory conditions [7][8].
Moreover, EEG-based BCIs are very sensitive to noise, e.g.,
user motions [9], as well as to the non-stationarity of EEG
signals [10]. Indeed, a BCI calibrated in a given context is
very likely to have much lower performances when used in
F. Yger is with the laboratory LAMSADE, Universit
´
e Paris-Dauphine, Place
du Mar
´
echal de Lattre de Tassigny, 75775 Paris Cedex 16, France, e-mail:
florian.yger@dauphine.fr.
M. Berar is with the laboratory LITIS, Universit
´
e de Rouen, Av-
enue de l’Universit
´
e, 76800 Saint Etienne du Rouvray, France e-mail:
maxime.berar@univ-rouen.fr.
F. Lotte is with Potioc / LaBRI, Inria Bordeaux Sud-Ouest, 200 avenue de
la vieille tour, 33405, Talence Cedex, France e-mail: fabien.lotte@inria.fr.
another context, see, e.g., [11][12]. Finally, in addition to
this reliability issue, BCIs also suffer from long calibration
times. This is due to the need to collect numerous training
EEG examples from each target user, to calibrate the BCI
specifically for this user, to maximize performances [13].
Therefore, for BCIs to be usable in practice, they must be
robust across contexts, time and users, and with calibration
times as short as possible. These challenges can be addressed
at multiple levels, e.g., at the neuroscience level, by identifying
new neurophysiological markers that are more reliable than the
ones currently used, at the human level, by training users to
gain accurate and stable control over the EEG patterns they
produce [14] or at the signal processing level, by building
features and classifiers that are robust to context changes,
and that can be calibrated with as little data as possible.
Regarding EEG signal processing, some recent results sug-
gest that a new family of approaches is very promising to
address the multiple challenges mentioned above: Riemannian
approaches [15][16][17]. These approaches enable the direct
manipulation of EEG signal covariance matrices and sub-
spaces, with an appropriate and dedicated geometry, the Rie-
mannian geometry. They have recently shown their superiority
to other classical EEG signal processing approaches based on
feature vector classification, by being the winning methods
on a couple of recent brain signal classification competi-
tions, notably the ”‘DecMEG2014”’ (https://www.kaggle.com/
c/decoding-the-human-brain) and the ”‘BCI challenge 2015”’
(http://neuro.embs.org/2015/bci-challenge/) competitions.
Since then, these approaches have witnessed an increased
enthusiasm from the research community, and have been used
to explore new feature representations, to learn features and
design robust classifiers as well as to generate artificial EEG
data [15][18][19][17][20]. Thus, time is now ripe to review
what these methods are, what they can already contribute to
BCI design, and how they should be further explored in the
future. This is what this review paper proposes.
This paper is organized as follows: Section II first presents
the standard design of a BCI, to illustrate how it employs
covariance matrices. Then, Section III proposes a brief tuto-
rial on Riemannian geometry, that makes possible the direct
manipulation of such covariance matrices. Then Section IV
reviews how such approaches have been used for EEG-based
BCI, in particular for subspace methods (e.g., spatial filtering,
see Section IV-A) and for covariance methods (to represent
and classify EEG signals, to learn metrics or reduce calibration
time, see Section IV-B). Finally, Section V identifies relevant
challenges and promising research directions for EEG signal
classification based on Riemannian approaches.
SPECIAL ISSUE ON BMI/BCI SYSTEMS IN IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING 2
II. STANDARD BCI DESIGN
A prominent type of BCIs is oscillatory activity-based
BCIs, that exploit amplitude changes in EEG oscillations.
They notably include motor imagery-based BCIs. Such a BCI
is typically designed around the Common Spatial Patterns
(CSP) algorithm to optimize spatial filters and the Linear
Discriminant Analysis (LDA) classifier. CSP aims at learning
spatial filters such that the variance of spatially filtered signals
is maximized for one class and minimized for the other class
[21]. Formally, optimizing CSP spatial filters w (w being a
column vector) consists in extremizing the following function:
J
CSP
(w) =
w
T
C
1
w
w
T
C
2
w
(1)
where C
j
is the average spatial covariance matrix
1
of the band-
pass filtered EEG signals from class j, and (.)
>
denotes the
transposition. Typically, these spatial covariance matrices are
obtained by computing the spatial covariance matrix S
j
i
from
each trial Z
j
i
from class j, and then averaging them:
C
j
=
1
N
j
N
j
X
i
S
j
i
=
1
N
j
N
j
X
i
Z
j
i
Z
j
i
>
(2)
with N
j
the number of trials in class j and Z
j
i
∈ R
N
c
×N
s
is
the j
th
EEG trial from class i, with N
s
the number of samples
in a trial, and N
c
the number of channels. Note that EEG
signals are usually band-pass filtered and thus have a zero
mean. Equation 1 being a generalized Rayleigh quotient, it
can be extremized by Generalized Eigen Value Decomposition
(GEVD) of the average covariance matrices C
1
and C
2
[21].
The spatial filters which maximize/minimize J
CSP
(w) are
the eigenvectors corresponding to the largest and smallest
eigenvalues of this GEVD, respectively. It is common to select
3 pairs of CSP filters w
i
, corresponding to the 3 largest and
smallest eigenvalues [21]. Once the filters w
i
are obtained,
a feature f
i
is computed as f
i
= log(w
T
i
Sw
i
), where S is
the current trial covariance matrix. The LDA classifier uses a
linear hyperplane to separate feature vectors from two classes
[22]. The intercept b and normal vector a of this hyperplane
are computed as follows:
a = C
−1
(µ
1
− µ
2
) and b = −
1
2
(µ
1
+ µ
2
)
T
a (3)
with µ
1
and µ
2
being the mean feature vectors for each class
and C the covariance matrix of both classes. It is interesting
to note that both CSP and LDA require estimating covariance
matrices. Although other types of BCI such as Event-Related
Potentials (ERP)-based BCI may not manipulate covariance
matrices as explicitly, we will see later that they can also be
represented using covariance matrices (see Section IV-B1). As
mentioned earlier, Riemannian geometry provides tool to di-
rectly manipulate those covariance matrices, without the need
for spatial filters. The next section provides an introductory
tutorial on this Riemannian geometry and associated tools.
1
throughout this manuscript, average covariance matrices estimated over
several trials will be denoted as C, whereas covariance matrices estimated
from a single trial will be denoted as S
Fig. 1. On the differential manifold M, the tangent space at X
0
is the set
of the velocities ˙γ(0) of the curves γ(t) passing through X
0
at t = 0.
III. CONCEPTS OF RIEMANNIAN GEOMETRY
To phrase it in an almost over-simplistic way, Riemannian
geometry is the branch of mathematics that studies smoothly
curved spaces that locally behave like Euclidean spaces.
Although it may look like an exotic mathematical tool, the
concept of Riemannian manifold is more common than one
can expect as the Earth is an example of a Riemannian
manifold and we will use this as a pedagogical example in
what follows. As this section will only give basic concepts
about Riemannian geometry, the reader can refer to [23], [24]
for an in depth exploration of the topic.
A. Riemannian manifolds
Riemannian manifolds are defined as the result of imbri-
cated mathematical structures in the manner of Russian dolls.
First, we need to define topological manifolds, which are
spaces in which every point has a neighborhood homeomor-
phic to R
n
. To simplify, a topological manifold can be seen
as a space that locally looks flat.
Endowed with a differential structure -also called atlas-
(i.e. a set of bijections called charts between a collection of
subsets of the topological manifold and a set of open subsets of
R
n
), the topological manifold becomes a differential manifold.
Within differential manifolds, there exist smooth manifolds
which are differential manifolds for which the transitions
between maps are smooth. To simplify again, this step aims at
giving rules for locally translating a point on the manifold to
its linear approximation. Those rules are local but on a smooth
manifold, the rules slightly change from one point to another.
On every point of a smooth differential manifold, the notion
of tangent space can be defined as the velocity of the curves
passing the point as illustrated in Fig. 1. A Riemannian
manifold is then a real smooth manifold equipped with an
inner product on the tangent space at each point.
For any Riemannian manifold, there exists a pair of map-
pings transporting points from the manifold to any given
tangent space and vice versa. More precisely, the Exponential
mapping transports a tangent vector (i.e. a point in a tangent
space) to the manifold and the Logarithmic mapping is locally
defined to transport a point in the neighborhood of a point to
the tangent space defined at this point. As a consequence, Rie-
mannian manifolds can be locally approximated by Euclidean
SPECIAL ISSUE ON BMI/BCI SYSTEMS IN IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING 3
Fig. 2. The Earth can be seen as a Riemannian manifold as it can be locally
linearly approximated. On earth and for Riemannian manifolds in general, a
map is only valid locally and implies deformations when used further.
spaces via their tangent spaces, however deformations occur
for points mapped far from where the tangent point is defined.
As depicted in Fig. 2, in our example of the Earth, the
topological manifold is the sphere that is locally homeomor-
phic to R
2
. Defining a differential structure on the Earth boils
down to defining geographic maps and to gather them in an
atlas. Once a scalar product is defined on each map, distances
can be computed locally and then extended to a global notion
of distance computed along curves on the manifold, called
geodesics. On an Earth map, the deformations of the shapes
of countries far from the center of the map results from the
approximation made by choosing a tangent plan defined at the
center of the map.
For example, as described in [25], on the set of Symmetric
Positive Definite (SPD) matrices P
n
-that will be presented
next-, the tangent space at a point X, T
X
P
n
, is the space
of symmetric matrices. Then, one choice (for making a Rie-
mannian manifold out of this space) is to equip every tangent
space with the following metric:
∀A, B ∈ T
X
P
n
hA, Bi
X
= Tr
X
−1
AX
−1
B
.
From this choice of metric follows the definition of a distance
in Eq. 8. Then, any symmetric matrix A belonging to T
X
P
n
,
the tangent space at X, can be mapped on P
n
(with the
reciprocal operation) as
S
A
= exp
G
(A) = X
1
2
exp
X
−
1
2
AX
−
1
2
X
1
2
, (4)
A = log
X
(S
A
) = X
1
2
log
X
−
1
2
S
A
X
−
1
2
X
1
2
(5)
with log(.) and exp(.) the matrix logarithm and exponential.
Given a differential manifold, a common way to create a
Riemannian manifold is to embed the tangent space with the
usual scalar product of the ambient space. Doing so creates a
Riemannian submanifold. However, for some cases (like the
space of SPD matrices), several scalar products are available,
leading as we will see, to different Riemannian geometries.
B. Bestiary of manifolds for BCI
When dealing with EEG signals, one has to manipulate
matrices under certain types of constraints. Most of the
0
0.5
1
1.5
0
0.2
0.4
0.6
0.8
1
1.2
−1.5
−1
−0.5
0
0.5
1
1.5
b
Euclidean and Riemannian geometry of 2x2 covariance matrices
a
c
Fig. 3. Comparison between Euclidean (straight dashed lines) and Riemannian
(curved solid lines) distances measured between 2 × 2 SPD matrices.
classical constraints used in BCI applications are smooth.
Hence, the space under constraints can be interpreted
as a smoothly curved space. Equipped with the proper
mathematical structure, those spaces can be handled as
Riemannian manifold. In terms of applications, Stiefel and
Grassman manifolds are well-suited to model subspaces
methods as in [26], [27], [28], [20] and Symmetric Positive
Definite (SPD) matrices naturally models covariance
matrices [15], [19], [16].
1) Manifolds of subspaces: In BCI applications, the
most prominent constraints are orthogonality constraints
(that occur in every eigenproblem for example). More
formally, the space of orthonormal matrices is defined as
the Stiefel manifold St(n, p) = {X ∈ R
n×p
|X
>
X = I
n×p
}.
At every point of this space, we have the tangent space
T
X
St(n, p) = {V ∈ R
n×p
|X
>
V + V
>
X = 0} equipped
with the usual Euclidean scalar product of R
n×p
.
The Stiefel manifold comprises several special
cases -like the orthogonal group O
p
= St(p, p)-
and variants -like the generalized Stiefel manifold
St
G
(n, p) = {X ∈ R
n×p
|X
>
GX = I
n×p
} [26]. Many
pre-processing, blind source separation or classification
methods currently used in practice (Principal Component
Analysis (PCA), Canonical Correlation Analysis (CCA), ...)
share a manifold structure of this kind. Moreover, equipped
with a group structure over the space O
p
, the Stiefel manifold
becomes a quotient space named the Grassman manifold
Gr(n, p). This manifold is particularly recommended when
a cost function is invariant over O
p
in order to reduce the
search space. For a complete overview about those manifolds,
the interested reader can refer to [29] and [24].
2) Manifold of covariances matrices: The space of SPD
matrices, noted P
n
= {X ∈ R
n×n
|X = X
>
, X 0} and
composed of symmetric matrices of strictly positive eigenval-
ues, can be successfully applied for manipulating covariance
matrices from EEG signals. At every point of the space of
SPD matrices, we have the tangent space T
X
P
n
homogeneous
to the space of n × n symmetric matrices. Depending on
the choice of scalar product to equip the tangent spaces,
one Euclidean and two different Riemannian geometries are
SPECIAL ISSUE ON BMI/BCI SYSTEMS IN IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING 4
0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
t
1.5
1.6
1.7
1.8
1.9
2.0
2.1
2.2
determinant
Comparison of Euclidean and Riemannian
interpolation along geodesics
Euclidean geometry
Riemannian geometry
Fig. 4. Illustration of the swelling effect : the evolution of the determinant
of matrices interpolationg two other matrices in the Euclidean sense (in blue)
and in the Riemannian sense (in red). In the Euclidean interpolation, the
determinant of the interpolation is bigger than the determinant of the points.
In the Riemannian interpolation, the determinant is constant.
possible and subsequently two different Riemannian distances
(and one Euclidean) can be defined:
• the Euclidean distance
2
δ
e
(X
A
, X
B
) = kX
A
− X
B
k
F
=
p
hX
A
− X
B
, X
A
− X
B
i
F
,
(6)
• the LogEuclidean distance [30], [19]
δ
l
(X
A
, X
B
) = klog (X
A
) − log (X
B
)k
F
, (7)
• the Affine Invariant Riemannian Metric (AIRM) dis-
tance [25]
δ
r
(X
A
, X
B
) = klog(X
−
1
2
A
X
B
X
−
1
2
A
)k
F
, (8)
where the log(.) corresponds to the matrix logarithm.
The difference between the Euclidean and Riemannian
geometries is illustrated in Fig. 3, where 2 × 2 SPD matrices
are represented as points in R
3
. The positivity constraint is
a cone, inside which SPD matrices lie strictly. The minimal
path between the two points for Euclidean distance is a straight
line, whereas the computed AIRM distances draw curves.
Despite its apparent simplicity, the Euclidean geometry
has several drawbacks and is not always well suited for
SPD matrices [31], [30], [32], which motivates the use of
Riemannian geometries. For example, the Euclidean geometry
induces some artifacts like the so-called swelling effect [30].
As illustrated in Fig. 4, this effect is observed in task as
simple as averaging two matrices, where the determinant of
the average is larger than any of the two matrices. Hence,
such an effect can be particularly harmful for data analysis
as it adds spurious variation to the data. Another drawback,
illustrated in Fig. 3 and already noticed in [31], is the fact that
Euclidean geometry for SPD matrices forms a non-complete
space. Hence, in this Euclidean space interpolation between
matrices is possible (but affected by the swelling effet), but
extrapolation may produce uninterpretable indefinite matrices.
2
Note that, in order to make it a Euclidean space, P
n
is equipped with the
Frobenius inner product hX
A
, X
B
i
F
= Tr
X
>
A
X
B
Fig. 5. Gradient descent on a Riemannian manifold : at the point X
0
, the
Euclidean gradient is projected on the tangent space at X
0
and then mapped
back on the manifold M.
As a consequence, some methods can be transformed into
their Riemannian counter-part by substituting the Euclidean
distance by a Riemannian distance. For example, in a metric
space -like a Riemannian manifold-, the Fr
´
echet mean extends
the concept of mean
3
. It can be defined as the element
¯
E of
the space minimizing the sum of its squared distances
4
δ to
the points E
i
of the dataset:
¯
X = argmin
X
N
X
i=1
δ
2
(X
i
, X). (9)
Contrary to the Euclidean mean, the Fr
´
echet mean usually does
not have any closed-form solution and it must be obtained by
solving the optimization problem in Eq 9. Such an averaging
method has attracted a lot of attention from the optimization
community and for example, for SPD matrices, several algo-
rithms have been proposed [34], [35], [36].
C. Optimization on Riemannian manifolds
Optimization on matrix manifolds is by now a mature field
with most of the classical optimization algorithms having a
(Riemannian) geometric counterpart [24]. In spirit, it consists
in optimizing a function over a curved space instead of
optimizing it under a smooth constraint in a big Euclidean
space. In Riemannian setting, descent directions are no longer
straight lines, but curves on the manifold and at every step,
an admissible solution is created. To account for this, the
algorithms must be slightly modified. Fig. 5 sums up the
general recipe for the implementation of Riemannian gradient
descent for a function f:
1) At each iteration, at point X
t
, transform the Euclidian
gradient D
X
t
f into a Riemannian gradient ∇
X
t
f.
2) Perform a line search along geodesics at X
t
in the
direction H = ∇
X
t
f. To do so the Riemannian gradient
∇
X
t
f is mapped to the geodesic using the exponential
map. For computational reasons, an approximation of
the exponential map, called a retraction, is often used
instead [24].
The reader interested by the details about optimization can
refer to [24] and to [29] for the special case of Stiefel and
Grassman manifolds. Most of the examples given in [24] have
been implemented efficiently in a MATLAB toolbox [37].
3
In the litterature, the Fr
´
echet mean is also sometimes called the Karcher
mean or the geometric mean.
4
Note that using a simple distance instead of the squared distance would
give a Fr
´
echet median [33].