{'fraction_non_alphanumeric': 0.03702663671373556, 'fraction_numerical': 0.008384146341463415, 'mean_word_length': 4.332406417112299, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Geometry-aware stationary subspace analysis

In many real-world applications data exhibits non-stationarity, i.e., its distribution changes over time. One approach to handling non-stationarity is to remove or minimize it before attempting to analyze the data. In the context of brain computer interface (BCI) data analysis this may be done by means of stationary subspace analysis (SSA). The classic SSA method finds a matrix that projects the data onto a stationary subspace by optimizing a cost function based on a matrix divergence. In this work we present an alternative method for SSA based on a symmetrized version of this matrix divergence. We show that this frames the problem in terms of distances between symmetric positive definite (SPD) matrices, suggesting a geometric interpretation of the problem. Stemming from this geometric viewpoint, we introduce and analyze a method which utilizes the geometry of the SPD matrix manifold and the invariance properties of its metrics. Most notably we show that these invariances alleviate the need to whiten the input matrices, a common step in many SSA methods which often introduces errors. We demonstrate the usefulness of our technique in experiments on both synthesized and real-world data.


Introduction
A common assumption in statistical modeling is that the distribution of observed data does not change over time, i.e., that it is stationary. In most cases it is this assumption of stationarity which allows results to be effectively generalized from the sample to the population. When the stationarity assumption is violated, as is often the case in real-world applications such as speech enhancement [11] or neurological data analysis [25], specialized machine learning methods must be developed in order to maintain adequate prediction capabilities.
A relatively well-studied non-stationary setting is covariate-shift [26], in which the input distribution changes but the conditional distribution of the outputs does not. The problem of covariate-shift has received growing attention in recent years, and many theoretical and practical aspects have been addressed (see [28] for an in-depth exploration of this topic).
These works typically do not aim to remove or reduce non-stationarity in the data, but rather they try to cope with its existence. A different approach is to remove, or minimize, any existing non-stationarities before attempting to analyze the collected data. In the context of brain computer interface (BCI) data analysis, two such note-worthy methods are stationary subspace analysis (SSA) [30] and stationary common spacial patterns (sCSP) [32].
Similar in spirit to independent component analysis (ICA) [16], SSA statistically models the data as a mixture of stationary and non-stationary signals. Unlike ICA, however, the signals are not assumed to be independent of each other. The data is first split into (possibly overlapping) time frames called epochs. Then a projection matrix is found by optimizing a cost function based on the divergence between distributions in various epochs.
The other method, by now quite a standard step for classification tasks in BCI systems, is the (supervised) sCSP method. Its goal is to project the data onto a subspace in which the various data classes are more separable. The sCSP method directs this subspace towards a stationary subspace by means of regularization.
Although SSA is essentially an unsupervised method, variations of it exist which are useful for supervised tasks such as classification [25]. These methods attempt to remove non-stationarity while keeping the discriminative inter-class variations intact. In this work we present an alternative method for stationary subspace extraction, focusing for the moment on the unsupervised setting.
Unlike SSA, the inputs to our method are not the raw signals themselves, but rather covariance matrices, computed from the signals (or from features based on the signals) using one of the many existing covariance estimators (e.g., [20]). Covariance matrices have gained increasing attention in recent years, and are now commonly used in many machine learning and signal processing applications such as computer vision applications [29], brain imaging [24] and BCI data analysis [6]. Their rich mathematical structure has been extensively studied [7], and advances in optimization methods on matrix manifolds in recent years have motivated the development of geometric methods for various machine learning tasks such as dictionary learning [9], metric learning [19] and dimensionality reduction [13]. This in turn motivates a covariance-based approach which utilizes the geometric properties of the symmetric positive definite (SPD) matrix manifold for SSA.
2 Geometry-aware stationary subspace analysis As discussed in the introduction, the task of extracting the stationary part from an observed mixture of stationary and non-stationary signals is essential in various applications. In this section we present our approach to this problem. We name it geometry-aware SSA (gaSSA) since we utilize the geometric properties of covariance matrices. To find a stationary subspace, SSA uses a cost function which is based on a matrix divergence. As a first step, we suggest a formulation that uses a symmetrized version of the same matrix divergence. We then show that this symmetrized matrix divergence is a distance between SPD matrices and offer a geometric interpretation to the problem of SSA.
We begin with a formal statement of the problem. To this end we provide a review of the original SSA model and framework [30]. Next, we present the symmetrized matrix divergence and discuss its relation to the Riemannian geometry of SPD matrices. Finally, we end this section by introducing and analyzing a general formulation for our geometry-aware SSA.

Stationary subspace analysis
where A ∈ GL D (R), the general linear group of size D over R, i.e., the set of all D × D invertible matrices with entries in R. The spaces spanned by the column vectors A s and A n are referred to as the s-space and n-space, respectively. The SSA model makes relatively few assumptions on the sand n-sources. First, the s-sources are stationary only in the weak (or wide) sense [17]. That is, their first and second moments are required to be constant in time. For the n-sources, their first two moments may vary between epochs of length T denoted by τ i = [t 0 (i), . . . , t 0 (i) + T ], where t 0 (i) is the start time of the i-th epoch. The sources do not necessarily follow a Gaussian distribution, but non-stationarities are assumed to be visible in the first two moments. Furthermore, this model does not assume that the sources are independent, namely, their covariance matrix is given by where Σ s ∈ R m×m , Σ n ∈ R (D−m)×(D−m) and Σ sn ∈ R m×(D−m) . Note that since Σ s is time independent we have dropped the notation (τ i ).
The goal of SSA is to find a de-mixing transformation A −1 that separates the s-sources from the n-sources. This matrix A −1 is not unique, but rather undetermined up to scaling, sign and linear transformations within the sand n-spaces. So, without loss of generality, the data may be centered and whitened such that the s-sources have a zero mean and a diagonal covariance matrix with unit variance. 1 Put differently, the de-mixing matrix is written as A −1 = BZ where Z = Cov(x) −1/2 is a whitening matrix created by a covariance estimator Cov(·) (in this case it is the empirical estimator) and B ∈ O D = V ∈ R D×D : V V = I , the set of all D × D orthogonal matrices. To find the matrix B, the signals are split into N epochs τ 1 , . . . , τ N of length T . Each epoch is characterized by its empirical mean µ i and covariance Σ i . Then for each epoch, the mean and covariance of the s-sources may be written as where I m D is the D × D identity matrix, truncated to the first m columns. Since the true µ s and Σ s are by definition stationary, the matrix B is that for which µ s i and Σ s i vary the least across all epochs. Owing to the maximum entropy principle, SSA uses the Kullback-Leibler (KL) divergence between Gaussian distributions to compare the epoch distributions up to their second moment. The matrix B is thus found by minimizing the following cost function: where D KL is the KL divergence and N (µ, Σ) denotes a multivariate Gaussian distribution with mean µ and covariance Σ.

Symmetrized matrix divergence
SSA and its variants use in their cost function the KL divergence between Gaussian distributions. In what follows we assume that these distributions have a zero mean.Under this assumption the KL divergence is a Bregman matrix divergence [4] between covariance matrices. The family of Bregman matrix divergences is generally defined as D KL is obtained for Φ(X) = − log det X, so it is often called the log-determinant divergence [18]. Bregman matrix divergences are useful in machine learning and have a number of useful properties [18], such as linearity and convexity in the first argument (and, in the case of the KL divergence, also in the second). However, as can be seen from their definition, they are asymmetric and do not satisfy the triangle equality. In particular, in our case we have that D KL (X || Y ) = D KL (Y || X) for two arbitrary matrices X = Y . Subsequently, symmetrized versions of the Bregman matrix divergence, in particular Jensen-Bregman divergences, have been studied in recent years [23]. In the case of the KL divergence (for zero mean distributions), this gives and is called the Jensen-Bregman log-determinant (JBLD) divergence [10]. The JBLD has many favorable properties (see [27]), primarily that its square root comprises a metric on the SPD matrix manifold. Moreover, in Sra [27] it has been shown that it is a close approximation to the affine invariant Riemannian metric (AIRM) [7] and shares many of its mathematical properties. The practical properties of both metrics, in particular their invariance properties, will be stationary sources non-stationary sources Figure 1: An illustration of our covariance-based approach. A set of input covariance matrices Σ i ∈ S + D are made up of a mixture (purple) of stationary (blue) and non-stationary (red) parts. Due to non-stationarities, in the original space the matrices are spread out. The matrices are mapped to two lower dimensional spaces, namely the s− (stationary) and n− (non-stationary) spaces, via the matrices Q and W , respectively. In the s-space the matrices are now more localized compared to the n-space in which they have a high variance. The spaces spanned by the columns of the matrices are orthogonal to each other, that is Q ⊥ W. discussed in Section 2.5. In the context of SPD matrices, the JBLD is referred to as the symmetric Stein divergence or the (square of the) log-determinant metric [27]. In the sequel we adopt the notation δ s (X, Y ).
Motivated by the above, we formulate SSA using a new cost function based on δ s (cf. Eqs. (3) and (4)): where Q = I m D B andΣ i = ZΣ i Z are the matrices whitened with Z =Σ −1/2 ,

A geometric interpretation
By replacing the cost function with one based on the symmetrized divergence we gain not only the beneficial properties of the symmetric divergence, but also new insight into the problem of SSA. First, note that in Eq. (7) the problem is ultimately framed in terms of distances between SPD matrices. This suggests adopting a geometric perspective, whereby the notion of stationarity is captured by the dispersion of the matrices Σ i . In this view, the assumption that the covariance matrices of stationary signals do not vary much between epochs translates to them having small distances between them. An illustration of this idea is presented in Fig. 1. In this figure, the matrices Σ i are seen as points on the SPD matrix manifold S + D . The goal of our method is to find transformations Q and W that map the matrices onto two separate manifolds of lower dimension -the stationary and non-stationary space, respectively. The matrices in the stationary space will exhibit small variation, while the non-stationarities will be captured in the non-stationary space where the variation of the matrices will be greater. The transformations Q and W may be chosen to be orthogonal to each other, producing well separated sand n-spaces. That is, W ∈ Q ⊥ for Q = span Q (likewise Q ∈ W ⊥ ), where ⊥ denotes the orthogonal complement. Put formally, our objective is to find a rank-m transformation matrix Q ∈ R D×m which maps Σ i ∈ S + D to Σ s i ∈ S + m for m < D such that the logdeterminant distance between the compressed centered matrices Σ s i = Q Σ i Q and their mean, which for the centered matrices is I, is minimized. Note that the space spanned by the columns of Q is of importance, and not the specific columns themselves. So, we may optimize Q over the Grassmann manifold [12], G = span(Q) : Q ∈ R D×m , Q Q = I , the set of all m-dimensional linear subspaces of R D×D .
In practice, for optimization we employ a Riemannian trust-regions method described in Absil et al. [2] and implemented efficiently in Boumal et al. [8].

A generic geometric formulation
Given the strong relation between the log-determinant metric and the AIRM [7], a natural progression is to incorporate the AIRM into the cost function. To understand why it would be beneficial to use the AIRM it is necessary first to briefly discuss the geometry of the SPD matrix manifold.
When equipped with the Frobenius inner product A, B F = tr(A B), the set S + n of SPD matrices of size n × n, belongs to a Euclidean space. In this case, similarity between SPD matrices can be measured simply by using the Euclidean distance derived from the Euclidean norm. This is readily seen in the following example for 2 × 2 SPD matrices. A matrix A ∈ S + 2 can be written as A = [ a c c b ] with ab − c 2 > 0, a > 0 and b > 0. Then matrices in S + 2 can be represented as points in R 3 and the constraints can be plotted as a convex cone whose interior is populated by the SPD matrices (see Fig. 2). In this representation, the Euclidean geometry of symmetric matrices then implies that distances are computed along straight lines.
Despite its simplicity, the Euclidean geometry has several drawbacks and is not always well suited for SPD matrices [13,3]. For example, due to an artifact referred to as the swelling effect [3], for a task as simple as averaging two matrices, it may occur that the determinant of the average is larger than any of the two matrices. Another drawback, illustrated in Fig. 2 and documented by Fletcher et al. [13], is the fact that this geometry forms a non-complete space. Hence, in this Euclidean space interpolation between SPD matrices is possible, but extrapolation may produce indefinite matrices, leading to uninterpretable solutions.
An efficient alternative which addresses these issues is to consider the space of SPD matrices as a curved space, namely a Riemannian manifold. Of the possible Riemannian distances, the AIRM, due to its favorable mathematical properties, is widely used in many applications (see, for example [13,24]). It is defined for any X, Y ∈ S + D as: where log(·) is the matrix logarithm and X 2 F = tr X X is the Frobenius norm.
In the curved space, the geodesics between matrices obtained by the AIRM are computed on curved lines as illustrated in Fig. 2 for the space S + 2 . Symmetric matrices with null and infinite eigenvalues (i.e., those which lie on the boundary of the convex cone, but not in it) are both at an infinite distance from any SPD matrix on the manifold (within the cone).
So, let us also consider a cost function of the same form defined w.r.t. the AIRM. 2 A general expression for our geometry-aware SSA (gaSSA) is then: Definition 2.1 (gaSSA) gaSSA w.r.t. a metric δ is defined as whereΣ i = ZΣ i Z are the matrices whitened with Z =Σ −1/2 andΣ = argmin is the matrix mean w.r.t. δ. We note that this cost function is similar in spirit to the one in Harandi et al. [14] and can be considered an unsupervised version of it. In the next section we will show that the need for matrix whitening can be alleviated. The Euclidean gradient w.r.t. Q of the cost function, used for the optimization, can be found in [15].

Symmetries and invariance properties
We now discuss the symmetries of our optimization problem and the invariance properties of our chosen metrics. These properties will enable us to significantly simplify our problem. For brevity we will state the results in terms of the AIRM, but the same holds true for the log-determinant metric.
Our key observation stems from the fact that δ r and δ s are invariant to congruent transformations of the form X → P H XP for P ∈ GL D (C) and P H = P is the conjugate transpose [7]. However, since our discussion is limited to real matrices, we have for X, Y ∈ S + D and a real-valued invertible matrix P . This is a crucial point since the whitening matrix Z and, more importantly, the mixing matrix A, act on the covariance matrices in this way.
n be a set of SPD matrices of size n×n and let Σ i = AΛ i A for some real-valued invertible matrix A. Denoting the Riemannian mean of Λ byΛ, the Riemannian meanΣ of the set is given byΣ = AΛA .
Using the above we obtain several useful equivalence relations.

Corollary 2.3
The following expressions are equivalent: where Γ i is the covariance matrix of the unmixed sources in the i-th epoch.
We have essentially shown that both the whitening operation and the mixing matrix A do not affect the distance between the covariance matrices of the original unmixed signals. We can then re-write our optimization problem as One may remark that A Q no longer has orthonormal columns and so does not belong to the Grassmann manifold. Indeed this is true. The final transition is due to the observation that the solution to our optimization problem is not unique. Rather, since we are interested in recovering the stationary subspace and not the exact sources themselves, the solution is invariant to any transformation (e.g., subspace scaling and rotation) acting within each of the sand n-spaces separately. Furthermore, we have chosen W ∈ Q ⊥ , and so the s-space is orthogonal to the n-space. Now, choosing orthogonal bases within each of the subspaces we may restrict ourselves to orthogonal mixing matrices A and find a transformation Q which lies in the Grassman manifold.
The final result is quite remarkable. First, it shows that our problem is essentially agnostic to the mixing matrix.Secondly, it eliminates the need to pre-whiten the matrices. This is useful in situations where it is not appropriate to whiten the data in advance, for example, when working with nearly illconditioned matrices. Furthermore, this eliminates any error that may be caused by inaccurate estimation of the matrix mean.
In conclusion, we have two variations of gaSSA given in the first and last terms of Eq. (12). The difference between the two is whether or not the input covariance matrices are whitened. Our analysis shows that whitening does not improve performance, and may in fact lead to a degradation of the results in certain cases. So, we claim that it is in general preferable not to whiten the matrices. In terms of the chosen metric, we do not expect a significant difference when using δ r vs. δ s . In the following section we will present experimental evidence to support these claims.

Experimental results
In this section we present experimental results on synthetic and data taken from a real BCI experiment. We compare the performance of gaSSA to the existing SSA and investigate the effects of matrix whitening and the choice of metric.

Toy data
For our first experiment we generated data following the SSA model as a mixture of stationary and non-stationary sources. To generate non-stationarity in the data we used a slightly modified version of the scheme provided in the SSA toolbox [22] and detailed in its user manual. Here we bring only a brief description: The elements of the mixing matrix A are chosen uniformly from the range [−0.5, 0.5] and its columns are normalized to 1. The distribution of the s-sources is constant over all epochs, namely s s (t) ∼ N (0, Λ s ). In the SSA toolbox, Λ s is taken to be the identity matrix, however, we choose Λ s to be a random matrix of the form Λ s = BΓB for an orthogonal matrix B and diagonal matrix Γ.
The n-sources are correlated with the s-sources, and for the i-th epoch τ i = [t 0 (i), . . . , t 0 (i) + T ] they are given by s The covariance matrices Λ n i are generated for each epoch in the same way as Λ s . So, the covariance matrix of the (unmixed) sources in the i-th epoch may be written as Using the data generated by the scheme above we compared the performance of our method to that of SSA. We used the AIRM and log-determinant metric both with and without matrix whitening.
As a performance measure we used the distance between the estimated nspace A n and the true n-space A n . This is owing to the fact that, as discussed in von Bünau et al. [31], the n-space and s-sources are identifiable, while the s-space and n-sources are not. To illustrate this, note that to be stationary, the s-sources must consist strictly of stationary sources, while the n-sources will remain non-stationary even if they include a mixture of stationary signals. The distance between sub-spaces is computed using δ G , the metric on the Grassmann manifold [1]. Shortly, this metric is based on the principal angles between the two spaces.
We generated 50 epochs of length T = 250 for several values of D and m. For each pair (D, m) we conducted the experiment 25 times. At each iteration the optimization procedure was restarted 5 times with different initial guesses and the transformation matrix which obtained the lowest cost was selected.
Due to length restrictions we describe only the results for D = 19 and m = 12. Other choices of parameters exhibited the same behavior. The n-space errors were identical for δ r and δ s (0.0067 ± 0.0001). They were consistently lower than SSA (0.0115 ± 0.0004). For both metrics the scheme without matrix whitening performed slightly better than those including whitening (0.0063 ± 0.0001 vs. 0.0067 ± 0.0001). These results are consistent with the analysis of Section 2.5.

Brain-computer interface
Next, we applied our method to data taken from the BCI competition IV dataset II. This dataset contains motor imagery (MI) EEG signals affected by eye movement artifacts. It was collected in a multi-class setting, with the subjects performing more than 2 different MI tasks. However, as in Lotte and Guan [21], we evaluate our algorithms on two-class problems by selecting only signals of left-and right-hand MI trials.
We applied the same pre-processing as described in Lotte and Guan [21]. EEG signals were band-pass filtered in 8 − 30 Hz, using a 5 th order Butterworth filter. For each trial, we extracted features from the time segment located from 0.5s to 2.5s after the cue instructing the subject to perform MI.
The data was initially divided into two parts: a training data set and a test data set. Similarly to von Bünau et al. [31] the first 20% of the test trials were set aside for adaptation. The aim of the adaptation part is to mitigate any non-stationarities between the test and the training session. We then learned the s-space in an unsupervised manner over the training and adaptation part. As before, our method was reinitialized 5 times and the transformation attaining the lowest cost was chosen.
The performance was measured by means of the classification rate on the test set. We used the following naive classifier, referred to as minimum distance to the mean (MDM) in Barachant et al. [5]: Using the labels of the training set, we compute the mean (in the s-space) for each of the two classes. Then, we classify the compressed covariance matrices in the test set according to their distance to the class means; each matrix is assigned the class to which it is closer. The original data is comprised of 22 signals. Since the true number of stationary signals is unknown, we repeated the experiment for several values in the range m ∈ [10,18]. Due to lack of space, we bring the results only for the values m = 10 and m = 14, the former being an example of challenging task and the latter an example of an intermediate one. The results for the nine subjects in the dataset are summarized in Table 1. The full results can be found in the supplementary material.
The results show that our method outperforms SSA for most subjects. As predicted, the methods without pre-whitening generally performed better than those which included whitening of the covariance matrices. In this more complex setting, more accurate estimation of the mixing matrix does not guarantee better classification. In terms of the metric, we see that δ r and δ s perform roughly the same.