# Structure-Preserving Transformers for Sequences of SPD Matrices

###### Abstract

In recent years, Transformer-based auto-attention mechanisms have been
successfully applied to the analysis of a variety of context-reliant data
types, from texts to images and beyond, including data from non-Euclidean
geometries. In this paper, we present such a mechanism, designed to classify
sequences of Symmetric Positive Definite matrices while preserving their
Riemannian geometry throughout the analysis. We apply our method to automatic
sleep staging on timeseries of EEG-derived covariance matrices from a standard
dataset, obtaining high levels of stage-wise performance.

Index Terms—  Transformers, SPD Matrices, Structure-Preserving,
Electroencephalography, Sleep Staging

##  1 Introduction

When analyzing the relationship between feature vectors or concurrent signals,
correlation and covariance matrices are a useful tool. Such matrices are at
least Positive Semi-Definite, and often fully Symmetric Positive Definite
(SPD). The set of $n\times n$ SPD matrices ($SPD(n)$) is a non-Euclidean,
Riemannian (i.e. metric) manifold, and the regular Euclidean operations of
most Neural Network (NN)-based models seldom preserve that geometric
structure, introducing deformations such as the “swelling effect” [1].
Structure-preserving NN-based approaches have been introduced [2, 3], deriving
their layers from one of two geodesic-defining metrics on $SPD(n)$. Affine
invariant metrics offer the best properties, but present computational
challenges (e.g. no closed-form formula for averaging) [4]. LogEuclidean
metrics are less isotropic, but still prevent swelling while being easier to
compute [1]. With $A,B\in SPD(n)$, we chose this LogEuclidean distance:

$\delta_{LE}(A,B)=\lVert log_{mat}(A)-log_{mat}(B)\rVert_{2}$ (1)

This metric relies on the matrix logarithm $log_{mat}(\cdot)$, bijectively and
isometrically mapping $SPD(n)$ onto $Sym(n)$, the vector space of $n\times n$
symmetric matrices (with $exp_{mat}(\cdot)$ being its inverse). Here, $\lVert
X\rVert_{2},X\in Sym(n)$ is the $\mathcal{L}_{2}$ norm applied to the upper
triangular of $X$. LogEuclidean operations are thus the Riemannian equivalent
to Euclidean operations on $Sym(n)$.

QLinearLinearKMatmulSoftmaxLinear CombinationMatmulV$(N\times L\times
d)$$(N\times L\times d)$$(N\times L\times L)$$h$$\begin{array}[]{c}([N\times
h]\times L\\\ \newline \times[d/h])\end{array}$ $([N\times h]\times L\times
L)$$(N\times L\times d)$

Fig. 1: The SP-MHA architecture. In parentheses are tensor dimensions at every
step, with $N$ the batch size. Fig. 2: SPDTransNet global architecture, with
$t=3$ feature tokens per epoch.

In this paper, we present a structure-preserving self-attention mechanism
applicable to sequences of SPD matrices, derived from the aforementioned
LogEuclidean metric. We embed said mechanism into a Transformer-based
architecture, and apply it to a biomedical classification problem.
Transformer-based technology has exploded in popularity ever since its
introduction in [5], with self-attention mechanisms being applied to very
different problems. With regards to Riemannian geometry, innovations seem
centered around the computation and application of attention maps,
specifically. For instance, Konstantinidis et al. [6] combine the standard
attention maps with Grassmann and SPD manifold-valued maps, to enrich their
computer vision model’s descriptive capabilities. By contrast, both He et al.
[7] and Li et al. [8] developed architectures to analyze 2D-manifold-valued
data in 3D space, the former achieving rotational equivariance with respect to
surfaces on the manifold and the latter developing two geodesic distances
applicable to point clouds, and building attention maps from these distances.
More generally, Kratsios et al. [9] provide a mathematical framework to apply
attention mechanisms on a variety of constrained sets, including manifolds.
While the latter approaches share our interest in preserving geometric
information, little to no focus is given to a Transformer’s other components.
As far as we are aware, ours is the only approach to apply structure-
preserving Transformers to SPD manifold-valued data.

##  2 SPD Structure-Preserving Attention

Let $B_{m}=\\{e_{i,j}\\}_{0<i\leq j}\subset\mathbb{R}^{m\times m}$ be the the
canonical basis of $Sym(m)$, with $(e_{i,j})_{i,j}=(e_{i,j})_{j,i}=1$, and all
other coefficients at 0. Let the triangular number $d(m)=\frac{m(m+1)}{2}$ be
the dimension of $Sym(m)$. Any matrix $M$ of $Sym(m)$ can be written in the
basis $B_{m}$ as a vector (a.k.a. “token”) of coordinates in
$\mathbb{R}^{d(m)}$. Therefore, any linear combination of these tokens would
equate to the same linear combination in $Sym(m)$, and thus to a LogEuclidean
weighted sum in $SPD(m)$, preserving its manifold structure.

###  2.1 Structure-Preserving Multihead Attention (SP-MHA)

In the original Linear Multihead Attention (L-MHA) component of Transformers
[5], the input tokens in the Q, K and V tensors are processed in parallel in
$h$ attention heads, then recombined through concatenation. There is no
guarantee that any underlying SPD structure in our tokens would survive this
concatenation. Echoing the similar concerns, Li et al. [8] decided to forego
having multiple heads. We chose instead to redefine the bloc, keeping the
parallel attention maps computation without sacrificing our data’s structure.

Let $d(m)$ be the dimension of input tokens. As seen in Figure 1, our SP-MHA
bloc does the following:

$MHA_{SP}(Q,K,V)=C\left(sm\left(\dfrac{\mathcal{L}_{Q}(Q)\cdot\mathcal{L}_{K}(K)^{T}}{\sqrt{d(m)/h}}\right)\right)\cdot
V$ (2)

with $\mathcal{L}_{Q}(\cdot)$ and $\mathcal{L}_{K}(\cdot)$ banks of $h$ linear
maps from $\mathbb{R}^{d(m)}$ to $\mathbb{R}^{\frac{d(m)}{h}}$, $sm(\cdot)$
the softmax function, and $C(\cdot)$ the weighted linear combination of the
$h$ post-softmax attention maps. Although said attention maps are identical to
their L-MHA counterparts, the only operation applied to V is the final matrix
multiplication, i.e. linear combinations of V’s tokens weighted by the
combined attention map, which do not threaten our tokens’ vector space
geometry.

###  2.2 Triangular linear maps

Let $Sym(n)$ and $Sym(m)$ have the canonical bases $B_{n}$ and $B_{m}$,
respectively. Let $\mathcal{L}_{n,m}(\cdot)$ be a linear map from $Sym(n)$ to
$Sym(m)$, represented by the matrix $W$ in $\mathbb{R}^{d(m)\times d(n)}$ with
respect to the bases (implemented in code through a fully connected NN layer
between tokenized matrices). We shall refer to such a map as a “triangular”
linear map.

Let $A^{*},B^{*}$ be in $SPD(n)$, mapped to $A,B\in Sym(n)$ through
$log_{mat}(\cdot)$. As $\mathcal{L}_{n,m}(\cdot)$ is a continuous linear map:

$\lVert\mathcal{L}_{n,m}(A)-\mathcal{L}_{n,m}(B)\rVert_{2}\leq\lVert
W\rVert_{*}\cdot\lVert A-B\rVert_{2}$ (3)
$\delta_{LE}(prj_{n,m}(A^{*}),prj_{n,m}(B^{*}))\leq\lVert
W\rVert_{*}\cdot\delta_{LE}(A^{*},B^{*})$ (4)

with $\lVert\cdot\rVert_{*}$ the matrix norm induced by the norm
$\lVert\cdot\rVert_{2}$, and
$prj_{n,m}(\cdot)=exp_{mat}\circ\mathcal{L}_{n,m}\circ log_{mat}(\cdot)$
mapping $SPD(n)$ onto $SPD(m)$. By definition of $\delta_{LE}$ (Equation 1),
Equations 3 and 4 are strictly identical. Hence, applying
$\mathcal{L}_{n,m}(\cdot)$ on our tokens is equivalent to applying
$prj_{n,m}(\cdot)$ on matrices in $SPD(n)$. The output tokens exhibit the
Riemannian structure of $SPD(m)$, and relations of proximity are preserved.
Therefore, so is the overall structure of our data.

Note that while other SPD-to-SPD NN-based mappings have been proposed [10, 2],
they rely on full-rank weight tensors, whereas $prj_{n,m}(\cdot)$ does not
require special constraints.

|  Model  | MF1 | Macro Acc. |  N1 F1  | Valid. metric |  Token dim. $d(m)$ |  # Feat. Tokens $t$
---|---|---|---|---|---|---|---
1 |  DeepSleepNet [11] |  78.14 $\pm$ 4.12 |  80.05 $\pm$ 3.47 |  53.52 $\pm$ 8.24  | N/A | N/A | N/A
2 |  IITNet [12] |  78.48 $\pm$ 3.15 |  81.88 $\pm$ 2.89 |  56.01 $\pm$ 6.54  | N/A | N/A | N/A
3 |  GraphSleepNet [13] |  75.58 $\pm$ 3.75 |  79.75 $\pm$ 3.41 |  50.80 $\pm$ 8.06  | N/A | N/A | N/A
4 |  Dequidt et al. [14] |  81.04 $\pm$ 3.26 |  82.59 $\pm$ 3.45 |  58.42 $\pm$ 6.09  | N/A | N/A | N/A
5 |  Seraphim et al. [15] |  79.78 $\pm$ 4.56 |  81.76 $\pm$ 4.61 |  58.43 $\pm$ 6.41  | MF1 | Concatenation | 1
6 |  SPDTransNet, $L=13$ |  81.06 $\pm$ 3.49 |  84.87 $\pm$ 2.47 |  60.39 $\pm$ 6.77  | MF1 |  351 ($m=26$) | 7
7 |  SPDTransNet, $L=21$ |  81.24 $\pm$ 3.29 |  84.40 $\pm$ 2.61 |  60.50 $\pm$ 6.18  | MF1 |  351 ($m=26$) | 10
8 |  SPDTransNet, $L=29$ |  80.83 $\pm$ 3.40 |  84.29 $\pm$ 2.65 |  60.35 $\pm$ 6.01  | N1 F1 |  351 ($m=26$) | 5
Table 1: Results obtained from both our model and the re-trained literature.
Best results are in bold.

##  3 Application to EEG Sleep Staging

The study of sleep most often requires the analysis of electrophysiological -
including electroencephalographic (EEG) - signals, subdivided into fixed-
length windows (“epochs”) and manually labeled with the appropriate sleep
stages, inferred from properties of the signal in and around each epoch [16].

As seen in a recent survey by Phan et al. [17], state-of-the-art automatic
sleep staging models typically use two-step architectures - given a sequence
of epochs, epoch-wise features are extracted before being compared at the
sequence-wise level, utilizing this contextual information to improve
classification. Since epochs often contain markers indicative of multiple
stages, two-step architectures tend to subdivide them further, extracting
features from subwindows using convolutional NNs [12] and/or recurrent NNs
[18, 19] \- the latter utilizing RNNs for both steps. Multiple authors have
adapted this context-inclusive approach to Transformer-based architectures
[20, 21, 15], with auto-attention mechanisms at both the intra- and inter-
epoch levels, taking advantage of the high performance they offer when applied
to sequence-based data.

###  3.1 Improving stage-wise classification

According to the aforementioned survey [17], current sleep staging models have
attained a sufficient performance level to replace manual staging in some
contexts. However, we have found that class-wise performance was often
lacking, particularly with regards to the N1 sleep stage [16], universally
difficult to classify (Section 4). Most EEG datasets are heavily imbalanced,
with the N1 stage often underrepresented (Section 3.3) - models optimized for
high overall accuracy may thus sacrifice N1 classification if it improves
global performance. To account for this, recent approaches [14, 15] elected to
primarily evaluate their performance through the macro-averaged F1 (MF1)
score, a class-wise balanced metric widely used in the literature. They also
rebalance their training sets through oversampling, so that all stages within
have the same number of classification targets. While the survey states that a
sequence-to-sequence classification scheme (classifying each epoch in the
input sequence) might lead to better performance, having multilabel inputs is
nonsensical for this rebalancing - hence their use of a sequence-to-epoch
scheme (classifying one epoch per sequence).

Seraphim et al. [15] hypothesized that an analysis through functional
connectivity - the activation correlations between different brain regions
[22] \- enhances stage-wise performance. Such an analysis was first done by
Jia et al. [13], using epoch-wise graph learning to estimate said connectivity
and sequence-wise spatio-temporal graph NNs to compare them. By contrast,
Seraphim et al. estimate it through covariance matrices. Their two-step model
uses standard Transformer encoders at each step, reminiscent of [21]. Each
input epoch is described as a multichannel timeseries of SPD matrices, which
are then tokenized bijectively. However, their approach does not guarantee the
preservation of their data’s SPD structure, as they operate a channel-wise
concatenation of their tokens, in addition to the concatenations found within
their encoders (Section 2.1). Hence, we propose a Transformer-based model
capable of analyzing EEG-derived functional connectivity through SPD matrices
without sacrificing the SPD structure of our data throughout the analysis.

###  3.2 The SPDTransNet model

As can be seen in Figure 2, our SPDTransNet model takes as input a sequence of
$L$ epochs, composed of a central epoch to classify and surrounding epochs to
provide context. Given $\ell$ the context size, we have $L=2\cdot\ell+1$.

Each EEG signal is decomposed into $C$ channels, divided into epochs, and
further subdivided into $S$ subwindows per epoch. After preprocessing (Section
3.3), each epoch is described by $S\times C$ matrices in $SPD(n)$. Each matrix
is mapped onto $Sym(n)$ logarithmically (Section 1), tokenized (Section 2),
and linearly mapped onto $Sym(m)$ (with $m>n$, as we have found that larger
tokens improve performance). The $S\times C$ grid of tokens is then arranged
into a 1D sequence, with the $S$ tokens in the channel 1 followed by the $S$
tokens in channel 2, etc.

At the intra-epoch level, a first positional encoding is applied to the
tokens, which pass through the first Transformer encoder. The $S\times C$
output tokens are then uniformly divided into $t$ groups, with each group
averaged into a single token. The $L$ sets of $t$ tokens are then regrouped at
the inter-epoch level, and passed through another positional encoding and
Transformer encoder pair. Finally, the $t$ tokens corresponding to the central
epoch (of index $\ell+1$ in Figure 2) go through two FC blocs111Fully
connected layers followed by ReLU activation and dropout layer., and are
mapped onto $\hat{y}_{\ell+1}\in\mathbb{R}^{c}$ by a final classification
linear map, with $c$ the number of classes.

We ensure structure preservation by using the SP-MHA bloc in all Transformer
encoders, and choosing all linear maps within said encoders’ Feed-Forward (FF)
components [5] and the aforementioned FC blocs to be triangular (Section 2.2).
The ReLU and dropout layers in the FF and FC blocs do not cause issue, as
setting a values within a token to 0 won’t remove the corresponding matrix
from $Sym(m)$. Same for the positional encodings, average poolings and in-
encoder layer normalizations, which all qualify as linear combinations.

As such, our model preserves the SPD structure of its input up to the final
classification map.

###  3.3 Dataset and preprocessing

We utilize the MASS SS3 dataset [23] due to its large number of available EEG
electrode-derived signals and widespread use in the literature. It is composed
of 62 full-night recordings of healthy subjects, segmented into 30s epochs.
Due to its nature, it is unbalanced, with the largest and smallest of its
$c=5$ classes (stages N2 and N1) composed of 50.24% and 8.16% of the dataset,
respectively. As do [14] and [15], we selected the 8 electrodes F3, F4, C3,
C4, T3, T4, O1 and O2.

To estimate functional connectivity from those signals, we apply the same
preprocessing pipeline as [15]222More details at
github.com/MathieuSeraphim/SPDTransNet., computing $S\times C=30\times 7$
covariance matrices in $SPD(8)$, with $S$ the sequence length and $C$ the
number of frequency-based channels. We then augment our matrices with signal-
derived information before whitening them, leading to more uniformly
distributed matrices in $SPD(9)$ (i.e. $n=9$). Said whitening requires the
computation of average covariance matrices per recording and channel, which
was done in [15] by computing the covariances over the entire recording.
Instead, we average all relevant matrices using the standard affine invariant
metric [4], improving performance.

##  4 Experiments & Results

To maximize class-wise performance, we operate a hyperparameter research per
configuration, followed by a 31-fold cros-validation. As do [14, 15] (Section
3.1), we rebalance all training sets and maximize the MF1 score. To explore
the importance of the context length $\ell$ (Section 3.2) within our model, we
ran hyperparameter researches with $\ell$ = 6, 10 or 14 (i.e. $L$ = 13, 21 or
29), with hyperparameter research configuration unchanged between them.

Our hyperparameter researches use the Optuna tool [24], with 5 simultaneous
runs and 50 total runs per configuration. Hyperparameters include the token
size $d(m)$, set by the first linear map (Section 3.2) and chosen in {351,
378} (i.e. $m$ $\in$ {26, 27})333Since $\frac{d(m)}{h}$ must be an integer,
potential values for those are limited.; the $h$ parameter of each Transformer
encoder, in {3, 9}; and the number of epoch feature tokens $t$ (Section 3.2),
chosen among {1, 3, 5, 7, 10} - with in particular $t=1$ akin to describing
each epoch with a single token, and $t=7$ corresponding to one token being
preserved per channel. We train all folds on the hyperparameters giving the
best validation MF1, as well as those with the best F1 score for the N1 stage.
Out of those two sets, the results from the set yielding the best average test
MF1 is presented in lines 6 to 8 of Table 1, with the corresponding
hyperparameter set, $d(m)$ and $t$ in the final three columns.

We compare ourselves to five models: DeepSleepNet [11], often used as a
benchmark, with a pre-trained epoch-wise global feature map submodel followed
by a sequence-to-sequence RNN; IITNet [12], the source of our 31 folds,
extracting multiple features per epoch through CNNs and comparing them through
sequence-wise RNNs; GraphSleepNet [13], expliciting epoch-wise functional
connectivity through graph learning; Dequidt et al. [14], utilizing a single-
step pretrained visual CNN, who both maximize MF1 performance and rebalance
training sets; and Seraphim et al. [15], with a similar approach to ours
lacking in structural preservation (Section 3.1). These models were re-trained
using our methodology - except for oversampling in DeepSleepNet’s sequence-to-
sequence submodel - though we use only their published hyperparameters.
Finally, as test sets vary between models due to sequence-based recording-wise
border effects, we trim test set borders to enforce uniformity. All these
changes cause the results we obtained to differ somewhat from those initially
published. Our results, averaged over all folds, are displayed in lines 1 to 5
of Table 1.

As shown in Table 1, we obtain the best MF1 and N1 F1 scores for $L=21$,
whereas the best macro-averaged accuracy is obtained for $L=13$. For all
values of $L$, we outperform the state-of-the-art on the considered metrics
(except for the MF1 score for $L=29$). Moreover, all three configurations have
around a two-point lead in both macro accuracy and N1 F1 score. While our
model favors the smaller token size of $d(m)=351$ for all values of $L$, it
seems that having a large number of tokens to describe each epoch (at least
$t=5$) is necessary for best performance. Overall, $L=21$ seems to be a good
compromise to capture enough contextual information without burdening our
model with irrelevant data.

##  5 Conclusion

We presented SP-MHA, a novel, structure-preserving Multihead Attention bloc,
and integrated it into our SPDTransNet model, designed to analyze SPD matrix
sequences. We proved said model’s capabilities through automatic EEG sleep
staging, obtaining a high level of per-stage performance relative to the
literature. Beyond this two-step analysis, SPDTransNet can be easily adapted
to a variety of problems, for instance by using only a single encoder step
and/or implementing a sequence-to-sequence classification scheme.

## References

  * [1] V. Arsigny et al.,  “Log-euclidean metrics for fast and simple calculus on diffusion tensors,”  Magnetic Resonance in Medicine, vol. 56, no. 2, pp. 411–421, 2006.
  * [2] Z. Huang and L. Van Gool,  “A riemannian network for spd matrix learning,”  in Association for the Advancement of Artificial Intelligence (AAAI), 2017.
  * [3] R. Chakraborty et al.,  “Manifoldnet: A deep neural network for manifold-valued data with applications,”  IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 44, no. 2, pp. 799–810, 2022.
  * [4] X. Pennec et al.,  “A riemannian framework for tensor computing,”  International Journal of Computer Vision, vol. 66, no. 1, pp. 41–66, Jan 2006.
  * [5] A. Vaswani et al.,  “Attention is all you need,”  in Advances in Neural Information Processing Systems, I. Guyon et al., Eds. 2017, vol. 30, Curran Associates, Inc.
  * [6] D. Konstantinidis et al.,  “Multi-manifold attention for vision transformers,”  ArXiv, vol. abs/2207.08569, 2022.
  * [7] L. He et al.,  “Gauge equivariant transformer,”  in Advances in Neural Information Processing Systems, M. Ranzato et al., Eds. 2021, vol. 34, pp. 27331–27343, Curran Associates, Inc.
  * [8] Z. Li et al.,  “Geodesic self-attention for 3d point clouds,”  in Advances in Neural Information Processing Systems, A. H. Oh et al., Eds., 2022.
  * [9] A. Kratsios et al.,  “Universal approximation under constraints is possible with transformers,”  in International Conference on Learning Representations, 2022.
  * [10] M. T. Harandi et al.,  “From manifold to manifold: Geometry-aware dimensionality reduction for spd matrices,”  in Computer Vision – ECCV 2014, D. Fleet et al., Eds., Cham, 2014, pp. 17–32, Springer International Publishing.
  * [11] A. Supratak et al.,  “Deepsleepnet: a model for automatic sleep stage scoring based on raw single-channel eeg,”  IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 25, no. 11, pp. 1998–2008, Nov 2017.
  * [12] H. Seo et al.,  “Intra- and inter-epoch temporal context network (iitnet) using sub-epoch features for automatic sleep scoring on raw single-channel eeg,”  Biomedical Signal Processing and Control, vol. 61, pp. 102037, 2020.
  * [13] Z. Jia et al.,  “Graphsleepnet: Adaptive spatial-temporal graph convolutional networks for sleep stage classification.,”  in IJCAI, 2020, pp. 1324–1330.
  * [14] P. Dequidt et al.,  “Automatic sleep stage classification on eeg signals using time-frequency representation,”  in Artificial Intelligence in Medicine, J. M. Juarez et al., Eds., Cham, 2023, pp. 250–259, Springer Nature Switzerland.
  * [15] M. Seraphim et al.,  “Temporal sequences of eeg covariance matrices for automated sleep stage scoring with attention mechanisms,”  in Computer Analysis of Images and Patterns, N. Tsapatsoulis et al., Eds., Cham, 2023, pp. 67–76, Springer Nature Switzerland.
  * [16] R. B. Berry et al.,  “Aasm scoring manual updates for 2017 (version 2.4),”  Journal of Clinical Sleep Medicine, vol. 13, no. 05, pp. 665–666, 2017.
  * [17] H. Phan and K. Mikkelsen,  “Automatic sleep staging of eeg signals: recent development, challenges, and future directions,”  Physiological Measurement, vol. 43, no. 4, pp. 04TR01, apr 2022.
  * [18] H. Phan et al.,  “Xsleepnet: Multi-view sequential model for automatic sleep staging,”  IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 44, no. 9, pp. 5903–5915, 2022.
  * [19] H. Phan et al.,  “L-seqsleepnet: Whole-cycle long sequence modelling for automatic sleep staging,”  IEEE Journal of Biomedical and Health Informatics, pp. 1–10, 2023.
  * [20] T. Zhu et al.,  “Convolution-and Attention-Based Neural Network for Automated Sleep Stage Classification,”  Int J Environ Res Public Health, vol. 17, no. 11, Jun 2020.
  * [21] H. Phan et al.,  “Sleeptransformer: Automatic sleep staging with interpretability and uncertainty quantification,”  IEEE Transactions on Biomedical Engineering, vol. 69, no. 8, pp. 2456–2467, 2022.
  * [22] S. Eickhoff and V. Müller,  “Functional connectivity,”  in Brain Mapping, A. W. Toga, Ed., pp. 187–201. Academic Press, Waltham, 2015.
  * [23] C. O’reilly et al.,  “Montreal archive of sleep studies: an open-access resource for instrument benchmarking and exploratory research,”  Journal of sleep research, vol. 23, no. 6, pp. 628–635, 2014.
  * [24] T. Akiba et al.,  “Optuna: A next-generation hyperparameter optimization framework,”  in Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining, 2019, pp. 2623–2631.