# A novel dual-stream time-frequency contrastive pretext tasks framework for
sleep stage classification

Sergio Kazatzidis, Siamak Mehrkanoon1<EMAIL_ADDRESS><EMAIL_ADDRESS>*Corresponding author. Department of Information and
Computing Sciences, Utrecht University, Utrecht, Netherlands

###### Abstract

Self-supervised learning addresses the challenge encountered by many
supervised methods, i.e. the requirement of large amounts of annotated data.
This challenge is particularly pronounced in fields such as the
electroencephalography (EEG) research domain. Self-supervised learning
operates instead by utilizing pseudo-labels, which are generated by pretext
tasks, to obtain a rich and meaningful data representation. In this study, we
aim at introducing a dual-stream pretext task architecture that operates both
in the time and frequency domains. In particular, we have examined the
incorporation of the novel Frequency Similarity (FS) pretext task into two
existing pretext tasks, Relative Positioning (RP) and Temporal Shuffling (TS).
We assess the accuracy of these models using the Physionet Challenge 2018
(PC18) dataset in the context of the downstream task sleep stage
classification. The inclusion of FS resulted in a notable improvement in
downstream task accuracy, with a 1.28 percent improvement on RP and a 2.02
percent improvement on TS. Furthermore, when visualizing the learned
embeddings using Uniform Manifold Approximation and Projection (UMAP),
distinct clusters emerge, indicating that the learned representations carry
meaningful information.

###### Index Terms:

Self-supervised learning, Sleep Staging, Pretext tasks

##  I Introduction

Large annotated data sets are generally necessary to train deep neural
networks. However, obtaining these data sets can be time-consuming and costly
[1]. Self-supervised learning has recently garnered considerable attention due
to its capacity to circumvent the need for labeled data, addressing this
fundamental challenge. This paradigm utilizes self-generated pseudo-labels
during the pre-training phase. These pseudo-labels are generated by pretext
tasks in order to obtain a meaningful representation of the data. The choice
of the appropriate pretext tasks is important in capturing the general
features of the data. Furthermore, these learned features could be utilized in
downstream tasks such as classification and segmentation [2]. Self-supervised
learning has been successfully applied in diverse fields such as Natural
Language Processing [3, 4, 5], and Computer Vision [6, 7, 8]. However, the
utilization of Self-supervised learning in the field of EEG processing is
still in its early phases. This is notable given the resource-intensive nature
of obtaining annotations in EEG research [9]. For instance, annotating a
24-hour EEG sleep signal segment requires roughly five hours of concentrated
effort by an expert [10].

The significance of sleep quality for overall health highlights the need for
thorough investigation, including the classification of sleep stages, which
can aid in detecting sleep disorders like sleep apnea. Deep learning networks,
particularly Convolutional Neural Networks (CNNs), offer a viable approach for
this classification task [11] and various other tasks [12, 13]. However, this
method still requires a substantial volume of labeled data. The authors in
[14] devised three pretext tasks, Relative Positioning (RP), Temporal
Shuffling (TS), and Contrastive Predictive Coding (CPC), to facilitate the
learning of meaningful representations of the data without the requirement for
large volumes of labeled data. In this paper, we specifically focus on RP and
TS. RP is a pretext task where the model determines whether a second segment
of EEG-data is temporally close or distant from the first segment, which is
referred to as the anchor segment, operating on the assumption that nearby EEG
data is more similar, while more distant data is less similar. In TS the model
determines whether three windows of EEG-data are ordered or shuffled in time.
In this paper we present an improvement of the aforementioned pretext tasks,
by incorporating the novel Frequency Similarity (FS) task. To this end, we
propose a dual-stream pretext tasks architecture that concatenates the learned
features of RP or TS with those of FS, introducing a frequency component to
the features to enhance the data representation. In the FS pretext task the
model evaluates whether a window exhibits greater similarity in the frequency
domain to the anchor window than another window.

This paper is organized as follows. Section II provides a brief overview of
related work, while Section III covers preliminaries. Section IV introduces
the FS pretext task and the novel architecture. Section V showcases the
obtained results, which are subsequently discussed. The conclusion is drawn in
Section VI.

##  II Related works

###  II-A Self-supervised learning

The major challenge of supervised machine learning techniques is their heavy
reliance on the amount and the quality of the labeled data. Data annotation
can be a laborious and expensive endeavor, motivating the exploration of self-
supervised learning as a viable alternative. In self-supervised learning one
uses the unannotated data to learn the underlying structure and to obtain a
rich representation of the data. The adoption of the right pretext tasks is
often pivotal in learning meaningful data representations [2]. Self-supervised
learning has been used mainly in the field of Computer Vision [6, 7, 8] and
Natural Language Processing [3, 4, 5], but in other fields such as Audio
Processing [15, 16], Video Processing [17, 18, 19] and EEG Processing [20, 21]
self-supervised learning is rapidly gaining popularity. Self-supervised
learning approaches can be broadly categorized into generative and contrastive
methods [2]. Generative methods aim to reconstruct the original input while
learning a meaningful representation [22]. In contrast, contrastive methods
adopt a discriminative approach, seeking to group similar samples closely
while pushing diverse samples apart [2]. An example of the generative method
can be observed in Bert [3], where certain tokens within a sentence are masked
and predicted by the model. An early instance of the contrastive method in
Computer Vision involves utilizing a jigsaw puzzle to learn the data’s
structure by reassembling shuffled patches back into the original image [7].

In this paper, we will focus on contrastive learning, a method that clusters
similar samples together while increasing the distance between dissimilar data
pairs. Positive samples, signifying similarity, are created through various
data augmentation techniques such as random cropping, resizing, color
distortions, and rotations [2]. These techniques are applied to the original
data to generate positive data pairs. Any other data is considered negative
and dissimilar with respect to the original data. A successful framework in
contrastive learning that employs these data augmentations is SimCLR [6],
which uses two different data augmentation techniques on the original image
and maximizes their agreement in order to find the underlying structure of the
original image. Another prominent contrastive framework is Momentum Contrast
(MoCo) [8], which introduces momentum updates to stabilize the rate of
representation change, enhancing representation consistency.

###  II-B EEG-signals

This central focus of this paper revolves around EEG data, with the primary
objective of uncovering its underlying structures and patterns for sleep stage
classification. In EEG research, data annotation is a laborious and costly
process, primarily carried out by domain experts. Moreover, the noise in the
data can complicate the visual interpretation. Self-supervised learning offers
a promising solution to these challenges by utilizing pseudo-labels generated
through pretext tasks [14]. Various self-supervised learning approaches have
been employed to unveil the underlying patterns in EEG data. For instance, in
studies like [23] and [24] multiple data augmentations are applied to EEG
data, including techniques such as time warping, Gaussian noise, additive
noise, amplitude scale, average filter and more. These augmentations are used
to create positive pairs from the same data. The aim is to bring their
representations close together, while increasing the gap between these
positive pairs and other negative pairs, which consist of the original data
and any other dissimilar data. Another approach, as described in [25],
involves devising a pretext task for the network to solve. In this case, a
task called “Temporal Rearrange” is introduced, where a window of EEG data is
divided into three blocks and shuffled (or not). The model’s task is to
predict whether shuffling occurred, thereby learning meaningful
representations from the EEG data. In some cases, multiple pretext tasks are
learned together, a strategy known as multi-task learning. For instance, in a
recent study [21], the authors simultaneously trained three distinct pretext
tasks: one focusing on spatial aspects, another centered around frequency-
related features, and the third utilizing contrastive learning. In contrast,
here we focus on learning features separately through the time-based and
frequency-based pretext tasks and subsequently concatenating them. In
particular, we utilize Relative Positioning (RP) and Temporal Shuffling (TS)
proposed in [14] for the time-based pretext tasks. We will then combine the
learned features of RP or TS with features learned by a novel frequency-based
task that will be introduced in Section 4. After learning and extracting the
valuable and necessary information from these pretext tasks, the knowledge is
applied to a downstream task. An example of a downstream task where the
information from the pretext tasks can be utilized is sleep stage
classification. In this task the model must classify the EEG data into one of
the five stages of sleep [9], which is crucial for identifying sleep
disorders. Currently, sleep stage scoring in EEG data is predominantly
performed manually by experts, a time-consuming and error-prone process [11].
The proposed self-supervised learning methods in this paper hold the potential
to automate and improve this process.

##  III Preliminaries

In this section the two utilized pretext tasks, i.e. Relative Positioning and
Temporal Shuffling [14], will be briefly explained.

###  III-A Relative Positioning

The Relative Positioning task (RP) requires a pair of two windows. The
underlying assumption of this pretext task is that neighboring windows
correspond to similar sleep stages, while distant windows represent different
sleep stages. The first window, referred to as the anchor window, is denoted
by $x_{t}$. The second window, denoted by $x_{t^{\prime}}$, is selected based
on two key parameters $\tau_{pos}$ and $\tau_{neg}$. The parameter
$\tau_{pos}$ indicates the time frame around the anchor window that is
considered close by. Positive pairs are composed of the anchor window and
another window within this positive context. The parameter $\tau_{neg}$
defines the negative context. This parameter indicates the minimum time frame
from the anchor window that is deemed sufficiently distant. Windows falling
inside this negative context are eligible to be utilized as negative samples,
forming negative pairs with the anchor window. The index $t$ associated with
each window indicates their starting time, i.e. the window $x_{t^{\prime}}$
has starting time $t^{\prime}$ . The labeling of the $i$-th pair is then
determined as follows [14]:

$y_{i}=\begin{cases}0,\ |t_{i}-t^{\prime}_{i}|>\tau_{neg}\\\ or\\\ 1,\
|t_{i}-t^{\prime}_{i}|\leq\tau_{pos}.\end{cases}$ (1)

In the study of RP, here in this paper a window of 30 seconds is utilized. In
addition, the parameters $\tau_{pos}$ and $\tau_{neg}$ are set to 1 minute and
15 minutes respectively.

(a)

(b)

(c)

Figure 1: (a) Architectural framework for pretext task learning. (b) Dual-
stream time-frequency framework with RP and FS for the downstream task. (c)
Dual-stream time-frequency framework with TS and FS for the downstream task.

To learn from this task two operations are defined: $h_{\theta}$ and $g_{RP}$.
$h_{\theta}$ is a feature extractor (embedder), in our case StagerNet [26],
which maps a window to a representation in the feature space. $g_{RP}$ is a
contrastive module that combines feature embeddings using element-wise
absolute differences. An overview of the RP process is shown in Fig. 1(a)
within the blue dashed box. Ultimately, a linear model with a dropout rate of
50% is employed to predict whether a given pair is positive or negative.

###  III-B Temporal Shuffling

The Temporal Shuffling (TS) pretext task requires three windows, denoted as
$x_{t}$, $x_{t^{\prime}}$, $x_{t^{\prime\prime}}$. Initially, we randomly
select the first window, denoted as the anchor window $x_{t}$. Subsequently,
another window $x_{t^{\prime\prime}}$ is chosen to be in close temporal
proximity to the anchor window, guided by the parameter $\tau_{pos}$. The
third window, $x_{t^{\prime}}$, can be positioned either in between windows
$x_{t}$ and $x_{t^{\prime\prime}}$ or in the negative context based on
$\tau_{neg}$. The parameter $\tau_{neg}$ indicates the minimum time frame from
the anchor window that is deemed sufficiently distant. Following this, we
proceed to create window triplets that can take one of two forms: they are
either arranged in temporal order ($t<t^{\prime}<t^{\prime\prime}$) or
shuffled ($t<t^{\prime\prime}<t^{\prime}$ or $t^{\prime}<t<t^{\prime\prime}$),
with $t$ representing the starting time of window $x_{t}$. The labeling of the
$i$-th triplet is then determined as follows [14]:

$y_{i}=\begin{cases}0,\ t^{\prime}_{i}\notin\ [t_{i},t^{\prime\prime}_{i}]\ \
\\\ or\\\ 1,\ t^{\prime}_{i}\in\ ]t_{i},t^{\prime\prime}_{i}[\;.\ \
\end{cases}$ (2)

In the study of TS, here in this paper a window of 30 seconds is utilized. In
addition, the parameters $\tau_{pos}$ and $\tau_{neg}$ are set to 30 minutes
and 120 minutes respectively. The TS learning process closely mirrors that of
Relative Positioning, and is visually represented within the red dashed box in
Fig. 1(a). After obtaining the embeddings $h_{\theta}(x_{t})$,
$h_{\theta}(x_{t^{\prime}})$ and $h_{\theta}(x_{t^{\prime\prime}})$, the
contrastive module $g_{TS}$ concatenates the absolute differences,
$|h_{\theta}(x_{t})-h_{\theta}(x_{t^{\prime}})|$ and
$|h_{\theta}(x_{t^{\prime}})-h_{\theta}(x_{t^{\prime\prime}})|$, as follows
[14]:

$\begin{IEEEeqnarraybox}[]{[}][c]{rCl}g_{TS}(h_{\theta}(x_{t}),h_{\theta}(x_{t^{\prime}}),h_{\theta}(x_{t^{\prime\prime}}))&=&(|h_{\theta}(x_{t})-h_{\theta}(x_{t^{\prime}})|,\\\
&&|h_{\theta}(x_{t^{\prime}})-h_{\theta}(x_{t^{\prime\prime}})|).\end{IEEEeqnarraybox}$
(3)

##  IV Method

###  IV-A Frequency Similarity

We propose a novel pretext task, referred to as Frequency Similarity (FS),
that operates in the frequency domain. Three distinct non-overlapping windows
($x_{t}$,$x_{t^{\prime}}$,$x_{t^{\prime\prime}}$) are randomly sampled from
the EEG recording, with $x_{t}$ serving as the anchor window. The only
requirement is that these windows must be distinct from one another.
Subsequently, the windows $x_{t}$, $x_{t^{\prime}}$ and $x_{t^{\prime\prime}}$
are transformed from the time domain into the frequency domain utilizing the
Welch method [27], resulting in $x_{f}$, $x_{f^{\prime}}$ and
$x_{f^{\prime\prime}}$. The Welch method is a well-established PSE technique
that estimates the distribution of power across different frequency components
within an EEG signal. This transformation allows the pretext task to operate
in the frequency domain, enabling the assessment of relationships between
these frequency representations. The Hellinger Distance (HD) is chosen as the
metric to assess the distance between the transformed windows $x_{f}$ and
$x_{f^{\prime}}$, as well as between $x_{f}$ and $x_{f^{\prime\prime}}$. This
choice is motivated by a study in [28], which demonstrated the effectiveness
of the Hellinger Distance for analyzing EEG data in the frequency domain. The
assignment of the label $y_{i}$ is based on whether the distance from
$x_{f^{\prime\prime}}$ to $x_{f}$ is less than or greater than the distance
from $x_{f^{\prime}}$ to $x_{f}$.

The labeling of the $i$-th triplet then is determined as follows:

$\begin{IEEEeqnarraybox}[]{[}][c]{rCl}y_{i}=\begin{cases}0,\
\frac{1}{N}\sum_{i=1}^{N}\frac{1}{\sqrt{2}}\cdot\lVert
x_{f}^{ch_{i}}-x_{f^{\prime}}^{ch_{i}}\rVert_{2}\\\ \ \ \ \
<\frac{1}{N}\sum_{i=1}^{N}\frac{1}{\sqrt{2}}\cdot\lVert
x_{f}^{ch_{i}}-x_{f^{\prime\prime}}^{ch_{i}}\rVert_{2}\\\ or\\\ 1,\
\frac{1}{N}\sum_{i=1}^{N}\frac{1}{\sqrt{2}}\cdot\lVert
x_{f}^{ch_{i}}-x_{f^{\prime}}^{ch_{i}}\rVert_{2}\\\ \ \ \ \
>\frac{1}{N}\sum_{i=1}^{N}\frac{1}{\sqrt{2}}\cdot\lVert
x_{f}^{ch_{i}}-x_{f^{\prime\prime}}^{ch_{i}}\rVert_{2}\;,\end{cases}\end{IEEEeqnarraybox}$
(4)

where N represents the total number of channels.

The feature extraction process in FS operates similarly to RP and TS. This
process is depicted in Fig. 1(a) within the dashed purple box. The windows are
mapped to a representation in the feature space using $h_{\theta}$, which, in
this case, is the StagerNet [26]. It’s important to note that the StagerNet
receives the three windows $x_{t}$, $x_{t^{\prime}}$ and
$x_{t^{\prime\prime}}$ of the time domain as input, but the pseudo-labels are
determined in the frequency domain. $g_{FS}$ aggregates the window
representations by concatenating the absolute differences between them,
similar to $g_{TS}$ as described in Eq. (3). Subsequently, the output of the
$g_{FS}$ is fed into a linear model with weights $w$ and bias weight $w_{0}$
which is utilized to predict whether the distance from $x_{f^{\prime\prime}}$
to $x_{f}$ is less than or greater than the distance from $x_{f^{\prime}}$ to
$x_{f}$. A dropout rate of 50% is applied to the linear model layer. The
selected loss function is binary cross-entropy, combined with a sigmoid
function. This choice is based on the nature of the output of the linear
model. The mathematical expression of the loss function is as follows:

$\mathcal{L}=-\left[y\cdot\log(\sigma(\hat{y}))+(1-y)\cdot\log(1-\sigma(\hat{y}))\right],$
(5)

where $y$ is the target value, $\hat{y}$ is the prediction of the linear model
and $\sigma$ is the sigmoid function. A visualization of the pretext task
Frequency Similarity is provided in Fig. 2.

Figure 2: Visual explanation of the introduced Frequency Similarity (FS)
pretext task. The EEG inputs are transformed with the Welch method into the
frequency domain. Subsequently, the distances between the windows are
computed, and labels are assigned based on these distances.

###  IV-B The proposed dual-stream architecture

In our approach, each pretext task, i.e. RP, TS and FS, trains its own
separate feature extractor denoted as $h_{\theta}$. Following pre-training,
these feature extractors utilize the original EEG dataset as input and
generate their new feature representation. Subsequently, the features obtained
from RP and TS are then individually concatenated with the features of FS, as
depicted in Fig. 1(b) and Fig. 1(c). This results in twice as many learned
features, encompassing information from both the original time-based pretext
tasks and frequency-based FS. As a result, the concatenated features
encapsulate knowledge of both time and frequency domains of the EEG data.
These combined features are subsequently frozen and utilized in the logistic
regression for the downstream task.

###  IV-C Downstream task

The downstream task addressed in this paper is sleep staging, which has been
previously studied in the literature, primarily using supervised deep learning
methods [29, 30]. However, more recent research papers, such as [14, 23, 24],
have demonstrated the successful application of Self-Supervised Learning (SSL)
for this task. Sleep staging involves classifying epochs of sleep EEG-data
into one of the five sleep stages: W (awake), N1, N2, N3 (different sleep
levels), and R (REM-sleep). The epochs considered here have a duration of 30
seconds. Automating this task is crucial as it currently relies heavily on
experts’ manual work [31].

###  IV-D Stagernet embedder

The chosen architecture for the embedder, denoted as $h_{\theta}$, is the
StagerNet [26]. StagerNet is a convolutional neural network (CNN) designed to
process EEG windows as input. It begins by applying a spatial convolution to
the input, followed by a permutation and a temporal convolution. Subsequently,
it employs a max-pooling layer, followed by another temporal convolution,
another max-pooling layer. Finally, it flattens the data and applies dropout,
resulting in a feature vector with 100 dimensions. To ensure a fair comparison
with the embedder in [14], we employed the same version of StagerNet as the
embedder. This modified version includes twice the number of convolutional
channels (16) compared to the original StagerNet [26] and also includes batch
normalization. The StagerNet employs an Adam Optimizer [32] with parameters
$\beta_{1}=0.9$ and $\beta_{2}=0.999$, and a learning rate $\alpha=$ 5 ×
$10^{-4}$. The training process was set to run for a maximum of 70 epochs or
until the validation loss stopped decreasing for a minimum of 10 consecutive
epochs.

###  IV-E Data

The models in this study are trained on a subset of the Physionet Challenge
2018 EEG [33] dataset, specifically the 994 participants from the training set
for whom sleep annotations were available. To create the training, validation,
and testing sets, a random split of 60% for training, 20% for validation, and
20% for testing is performed, resulting in 595, 199, and 199 recordings in
each respective set. For the pretext tasks RP, TS, and FS, 2,000 pairs or
triplets of windows were randomly sampled from each recording. This EEG
dataset comprises six EEG channels (F3-M2, F4-M1,C3-M2, C4-M1, O1-M2 and
O2-M1) and employs a sampling frequency of 200 Hz. The stages are categorized
as W (wakefulness), N1, N2, N3 (non-REM sleep stages) and R (REM sleep).
Following the lines of [14], first the EEG-waves are filtered applying a 30 Hz
Finite Impulse Response lowpass filter with a Hamming window, to filter out
higher frequencies, which are not relevant for sleep staging. In addition, as
in [14], only the channels F3-M2 and F4-M1 are retained to reduce data
dimensionality. Subsequently, the data is downsampled to 100 Hz. Non-
overlapping 30-second windows are subsequently extracted and then normalized
per channel to have zero-mean and unit standard deviation. The proposed code
is available at GitHub111https://github.com/Serg99io/DSTF-sleepstaging.

##  V Results and discussion

###  V-A Sleep stage classification with different pretext tasks

In this section, we evaluate the performance of the discussed pretext tasks,
Relative Positioning (RP), Temporal Shuffling (TS), and Frequency Similarity
(FS), within the context of the downstream task of sleep stage classification.
This evaluation is conducted utilizing all available labeled examples,
amounting to a total of 521,943 labeled examples. During the initial phase of
pre-training, we train distinct embedders for each of the pretext tasks RP, TS
and FS. Subsequently, the training, validation and test sets are provided as
input to these pre-trained embedders to obtain their respective
representations. Each distinct embedder yields a feature vector with 100
dimensions, capturing the essential information found in an EEG window. The
features obtained with FS are then concatenated with those of RP and TS
separately, resulting in a feature vector with a dimensionality of 200. For
sleep stage classification, a linear logistic regression model with
L2-regularization is employed. We have selected weighted average precision,
weighted average recall and balanced accuracy as our primary evaluation
metrics. The metrics are formally expressed as follows:

$\displaystyle Weighted\ avg\
precision=\frac{1}{N}\sum_{i=1}^{C}n_{c}\frac{TP_{i}}{TP_{i}+FP_{i}}\ ,$ (6)
$\displaystyle Weighted\ avg\
recall=\frac{1}{N}\sum_{i=1}^{C}n_{c}\frac{TP_{i}}{TP_{i}+FN_{i}}\ ,$ (7)
$\displaystyle Balanced\
accuracy=\frac{1}{C}\sum_{i=1}^{C}\frac{TP_{i}}{TP_{i}+FN_{i}}\ \ ,$ (8)

where $N$ represents the total number of instances, $C$ represents the number
of classes, $n_{c}$ represents the number of instances per class, $TP_{i}$,
$FP_{i}$ and $FN_{i}$ represent the number of true positives, false positives
and false negatives for class $i$ respectively. The obtained results are
summarized in Table LABEL:accuracy.

TABLE I:  Balanced accuracy and additional metrics for sleep stage
classification on the Physionet Challenge 2018 EEG dataset utilizing all
labeled examples.

| Pretext
---
tasks
| Balanced
---
accuracy
| Weighted avg
---
precision
| Weighted avg
---
recall
RP [14] | 0.6942 | 0.69 | 0.67
TS [14] | 0.6829 | 0.68 | 0.66
FS | 0.6365 | 0.65 | 0.62
RP + FS | 0.7070 | 0.70 | 0.69
TS + FS | 0.7031 | 0.70 | 0.68

One can observe that incorporating FS with either RP or TS leads to noticeable
improvements in accuracy compared to using RP or TS individually.

In addition, we also examine the performance of the models by varying the
number of labeled examples per class. In particular, we evaluate the models on
1, 10, 100, 1000, 10000, 50000, or all labeled examples per class and average
the results over five iterations. The obtained results are depicted in Fig. 3.

(a)

(b)

Figure 3: (a) Test balanced accuracy of sleep stage classification with RP and
FS features at different levels of labeled examples. (b) Test balanced
accuracy of sleep stage classification with TS and FS features at different
levels of labeled examples.

The results show that the incorporation of FS consistently improves
performance across different levels of labeled examples, except when utilizing
one labeled example per class. This discrepancy may be attributed to increased
variability in the selection of labeled examples with a single instance, a
phenomenon also observed in a previous study [14]. The improvement observed
across different levels of labeled examples can be attributed to FS offering
valuable additional and distinct information associated with the frequency
domain. Fig. 4 displays the accuracy achieved for each individual sleep stage
category. Notably, when considering the integration of FS, the results show
either equal or improved accuracy in comparison to the original RP and TS
pretext tasks. This finding underlines the significant value added by FS as an
enhancement over RP and TS. It’s important to emphasize that while FS alone
may not reach the same performance levels as RP and TS, it provides additional
insights from the frequency domain, ultimately leading to an improvement in
overall performance.

Figure 4: Test accuracy of each of the different pretext tasks and their
combinations over the five sleep stages, W (wakefulness), N1, N2, N3 (non-REM
sleep stages) and R (REM sleep).

Fig. 5 illustrates the confusion matrices of the pretext tasks. It can be
observed that when incorporating FS into RP or TS, the accuracy remains equal
or increases. This improvement primarily results from a slight enhancement in
the classification of each sleep stage. Instead of solely improving the
classification of one class, it contributes to a modest improvement across all
classes. Furthermore, an interesting observation from Fig. 5(a) is that the
majority of misclassifications of the wakefulness (W) stage occur when it is
predicted as N1. This could be due to the close temporal proximity between the
N1 and W sleep stages.

(a)

(b)

(c)

(d)

Figure 5: (a) The confusion matrix of RP. (b) The confusion matrix of TS. (c)
The confusion matrix of RP and FS. (d) The confusion matrix of TS and FS.

###  V-B U-MAP visualizations of learned features

To gain deeper insights into the results and the features acquired through
these pretext tasks, we utilize UMAP (Uniform Manifold Approximation and
Projection) [34], a dimensionality reduction technique designed to preserve
the global data structure. UMAP visualizations provide insight into the
meaningful representations acquired through self-supervised learning. In Fig.
6, the UMAP visualizations of the features of RP and the features of RP and FS
are presented. These visualizations clearly reveal the formation of distinct
clusters in both the features of RP and of RP and FS, indicating that the
learned features have meaningful representations of the original EEG data.
However, some overlap between clusters suggests that additional information
may be required for more precise sleep stage predictions, especially for
distinguishing between N1 and N2 stages. The features of RP and FS offer an
improvement over solely the features of RP as it enhances the separation
between stages W (wakefulness) and N3. Notably, in Fig. 4, stage N3 exhibits
the highest accuracy and a similar pattern is observed in Fig. 6 where the
clearest clusters are predominantly associated with this stage.

(a)

(b)

Figure 6: The UMAP visualizations of the five sleep stages W (wakefulness),
N1, N2, N3 (non-REM sleep stages) and R (REM sleep). (a) Relative Positioning
(RP). (b) Relative Positioning + Frequency Similarity (RP + FS).

(a)

(b)

Figure 7: The UMAP visualizations of the five sleep stages W (wakefulness),
N1, N2, N3 (non-REM sleep stages) and R (REM sleep). (a) Temporal Shuffling
(TS). (b) Temporal Shuffling + Frequency Similarity (TS + FS).

UMAP visualizations of TS and FS are depicted in Fig. 7. Various clusters are
visible, and the incorporation of Frequency Similarity features appears to
improve the distinction between stages N3 and R. This observation is further
supported by Fig. 5(b) and Fig. 5(d), where the incorporation of FS results in
reduced misclassifications between N3 and R stages.

##  VI Conclusion

In conclusion, in this paper we propose an improvement to the Relative
Positioning (RP) and Temporal Shuffling (TS) pretext tasks by introducing the
Frequency Similarity (FS) task, which integrates frequency information into
time-based features. This augmentation leads to notable performance
improvements, demonstrating the potential of self-supervised learning in sleep
staging tasks. Furthermore, UMAP visualizations reveal clusters within the
obtained EEG representations, confirming their meaningfulness. The proposed
approach facilitates the identification of sleep disorders like sleep apnea
and narcolepsy while reducing the need for costly and time-consuming manual
annotations.

## References

  * [1] L. Jing and Y. Tian, “Self-supervised visual feature learning with deep neural networks: A survey,” IEEE transactions on pattern analysis and machine intelligence, vol. 43, no. 11, pp. 4037–4058, 2020.
  * [2] A. Jaiswal, A. R. Babu, M. Z. Zadeh, D. Banerjee, and F. Makedon, “A survey on contrastive self-supervised learning,” Technologies, vol. 9, no. 1, p. 2, 2020.
  * [3] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep bidirectional transformers for language understanding,” arXiv preprint arXiv:1810.04805, 2018.
  * [4] T. Mikolov, E. Grave, P. Bojanowski, C. Puhrsch, and A. Joulin, “Advances in pre-training distributed word representations,” arXiv preprint arXiv:1712.09405, 2017.
  * [5] L. Floridi and M. Chiriatti, “Gpt-3: Its nature, scope, limits, and consequences,” Minds and Machines, vol. 30, pp. 681–694, 2020.
  * [6] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A simple framework for contrastive learning of visual representations,” in International conference on machine learning, pp. 1597–1607, PMLR, 2020.
  * [7] M. Noroozi and P. Favaro, “Unsupervised learning of visual representations by solving jigsaw puzzles,” in European conference on computer vision, pp. 69–84, Springer, 2016.
  * [8] K. He, H. Fan, Y. Wu, S. Xie, and R. Girshick, “Momentum contrast for unsupervised visual representation learning,” in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 9729–9738, 2020.
  * [9] X. Jiang, J. Zhao, B. Du, and Z. Yuan, “Self-supervised contrastive learning for eeg-based sleep staging,” in 2021 International Joint Conference on Neural Networks (IJCNN), pp. 1–8, IEEE, 2021.
  * [10] R. K. Malhotra and A. Y. Avidan, “Sleep stages and scoring technique,” Atlas of sleep medicine, pp. 77–99, 2013.
  * [11] H. W. Loh, C. P. Ooi, J. Vicnesh, S. L. Oh, O. Faust, A. Gertych, and U. R. Acharya, “Automated detection of sleep stages using deep learning techniques: A systematic review of the last decade (2010–2020),” Applied Sciences, vol. 10, no. 24, p. 8963, 2020.
  * [12] L. Alzubaidi, J. Zhang, A. J. Humaidi, A. Al-Dujaili, Y. Duan, O. Al-Shamma, J. Santamaría, M. A. Fadhel, M. Al-Amidie, and L. Farhan, “Review of deep learning: Concepts, cnn architectures, challenges, applications, future directions,” Journal of big Data, vol. 8, pp. 1–74, 2021.
  * [13] S. Mehrkanoon, “Deep shared representation learning for weather elements forecasting,” Knowledge-Based Systems, vol. 179, pp. 120–128, 2019.
  * [14] H. Banville, O. Chehab, A. Hyvärinen, D.-A. Engemann, and A. Gramfort, “Uncovering the structure of clinical eeg signals with self-supervised learning,” Journal of Neural Engineering, vol. 18, no. 4, p. 046020, 2021.
  * [15] H. Nguyen, F. Bougares, N. Tomashenko, Y. Esteve, and L. Besacier, “Investigating self-supervised pre-training for end-to-end speech translation,” in Interspeech 2020, 2020.
  * [16] S. Liu, A. Mallol-Ragolta, E. Parada-Cabaleiro, K. Qian, X. Jing, A. Kathan, B. Hu, and B. W. Schuller, “Audio self-supervised learning: A survey,” Patterns, vol. 3, no. 12, 2022.
  * [17] I. Misra, C. L. Zitnick, and M. Hebert, “Shuffle and learn: unsupervised learning using temporal order verification,” in Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11–14, 2016, Proceedings, Part I 14, pp. 527–544, Springer, 2016.
  * [18] R. Qian, T. Meng, B. Gong, M.-H. Yang, H. Wang, S. Belongie, and Y. Cui, “Spatiotemporal contrastive video representation learning,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 6964–6974, 2021.
  * [19] M. C. Schiappa, Y. S. Rawat, and M. Shah, “Self-supervised learning for videos: A survey,” ACM Computing Surveys, vol. 55, no. 13s, pp. 1–37, 2023.
  * [20] Z. Zhang, S.-h. Zhong, and Y. Liu, “Ganser: A self-supervised data augmentation framework for eeg-based emotion recognition,” IEEE Transactions on Affective Computing, 2022.
  * [21] Y. Li, J. Chen, F. Li, B. Fu, H. Wu, Y. Ji, Y. Zhou, Y. Niu, G. Shi, and W. Zheng, “Gmss: Graph-based multi-task self-supervised learning for eeg emotion recognition,” IEEE Transactions on Affective Computing, 2022.
  * [22] X. Liu, F. Zhang, Z. Hou, L. Mian, Z. Wang, J. Zhang, and J. Tang, “Self-supervised learning: Generative or contrastive,” IEEE transactions on knowledge and data engineering, vol. 35, no. 1, pp. 857–876, 2021.
  * [23] X. Jiang, J. Zhao, B. Du, and Z. Yuan, “Self-supervised contrastive learning for eeg-based sleep staging,” in 2021 International Joint Conference on Neural Networks (IJCNN), pp. 1–8, IEEE, 2021.
  * [24] M. N. Mohsenvand, M. R. Izadi, and P. Maes, “Contrastive representation learning for electroencephalogram classification,” in Machine Learning for Health, pp. 238–253, PMLR, 2020.
  * [25] Y. Ou, S. Sun, H. Gan, R. Zhou, and Z. Yang, “An improved self-supervised learning for eeg classification,” Math. Biosci. Eng, vol. 19, pp. 6907–6922, 2022.
  * [26] S. Chambon, M. N. Galtier, P. J. Arnal, G. Wainrib, and A. Gramfort, “A deep learning architecture for temporal sleep stage classification using multivariate and multimodal time series,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 26, no. 4, pp. 758–769, 2018.
  * [27] P. K. Rahi, R. Mehra, et al., “Analysis of power spectrum estimation using welch method for various window techniques,” International Journal of Emerging Technologies and Engineering, vol. 2, no. 6, pp. 106–109, 2014.
  * [28] G. Chen, G. Lu, Z. Xie, W. Shang, et al., “Anomaly detection in eeg signals: a case study on similarity measure,” Computational intelligence and neuroscience, vol. 2020, 2020.
  * [29] M. Abou Jaoude, H. Sun, K. R. Pellerin, M. Pavlova, R. A. Sarkis, S. S. Cash, M. B. Westover, and A. D. Lam, “Expert-level automated sleep staging of long-term scalp electroencephalography recordings using deep learning,” Sleep, vol. 43, no. 11, p. zsaa112, 2020.
  * [30] H. Phan, K. Mikkelsen, O. Y. Chén, P. Koch, A. Mertins, and M. De Vos, “Sleeptransformer: Automatic sleep staging with interpretability and uncertainty quantification,” IEEE Transactions on Biomedical Engineering, vol. 69, no. 8, pp. 2456–2467, 2022.
  * [31] H. Phan and K. Mikkelsen, “Automatic sleep staging of eeg signals: recent development, challenges, and future directions,” Physiological Measurement, vol. 43, no. 4, p. 04TR01, 2022.
  * [32] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
  * [33] M. M. Ghassemi, B. E. Moody, L.-W. H. Lehman, C. Song, Q. Li, H. Sun, R. G. Mark, M. B. Westover, and G. D. Clifford, “You snooze, you win: the physionet/computing in cardiology challenge 2018,” in 2018 Computing in Cardiology Conference (CinC), vol. 45, pp. 1–4, IEEE, 2018.
  * [34] L. McInnes, J. Healy, and J. Melville, “Umap: Uniform manifold approximation and projection for dimension reduction,” arXiv preprint arXiv:1802.03426, 2018.


# Some monoids of Pisot matrices

Artur Avila and Vincent Delecroix

###### Abstract

A matrix norm gives an upper bound on the spectral radius of a matrix.
Knowledge on the location of the dominant eigenvector also leads to upper
bound of the second eigenvalue. We show how this technique can be used to
prove that certain semi-group of matrices arising from continued fractions
have a Pisot spectrum: namely for all matrices in this semi-group all
eigenvalues except the dominant one is smaller than one in absolute value.

##  1 Introduction

A _dominant eigenvalue_ of a real square matrix is an eigenvalue of maximum
modulus. We call a square matrix _Pisot_ if it has non-negative integer
entries, its dominant eigenvalue is simple and all eigenvalues different from
the dominant one have absolute values less than one. We prove that several
monoids of non-negative matrices enjoy the property of all being Pisot.

Our first family of matrices is related to the so called _fully subtractive_
(multidimensional) continued fraction algorithm. For an integer $d\geq 2$ we
define for each $k=1,\ldots,d$ the matrix $A_{FS,d}^{(k)}$ by

$(A_{FS,d}^{(k)})_{ij}=\left\\{\begin{array}[]{ll}1&\text{if $j=k$ or
$i=j$},\\\ 0&\text{otherwise}\end{array}\right.$

For $d=3$ this boils down to the three matrices

$A_{FS,3}^{(1)}=\begin{pmatrix}1&0&0\\\ 1&1&0\\\ 1&0&1\end{pmatrix},\qquad
A_{FS,3}^{(2)}=\begin{pmatrix}1&1&0\\\ 0&1&0\\\ 0&1&1\end{pmatrix},\qquad
A_{FS,3}^{(3)}=\begin{pmatrix}1&0&1\\\ 0&1&1\\\ 0&0&1\end{pmatrix}.$

All non-degenerate products of the matrices $A_{FS,d}^{(k)}$ satisfy the Pisot
property.

######  Theorem 1.

Let $A=A_{FS,d}^{(i_{1})}A_{FS,d}^{(i_{2})}\ldots A_{FS,d}^{(i_{n})}$ be a
product of the fully subtractive matrices in dimension $d$. Then the matrix
$A$ is primitive if and only if all letters $\\{1,\ldots,d\\}$ appear in the
sequence $(i_{1},i_{2},\ldots,i_{n})$. Moreover, if the matrix $A$ is
primitive then it is Pisot.

Recall that a non-negative square matrix $A$ is primitive if there exists a
positive integer $n$ so that $A^{n}$ has all its entries positive. The case
$d=3$ of Theorem 1 was proved in [ArIt01]. The authors used an induction on
characteristic polynomials and our approach is radically different.

The same result holds for another set of $3\times 3$ matrices related to the
Brun multidimensional continued fractions. Let

$A_{Br}^{(1)}=\begin{pmatrix}1&1&0\\\ 0&1&0\\\ 0&0&1\end{pmatrix},\qquad
A_{Br}^{(2)}=\begin{pmatrix}1&1&0\\\ 1&0&0\\\ 0&0&1\end{pmatrix},\qquad
A_{Br}^{(3)}=\begin{pmatrix}1&0&1\\\ 1&0&0\\\ 0&1&0\end{pmatrix}.$

######  Theorem 2.

Let $A=A_{Br}^{(i_{1})}A_{Br}^{(i_{2})}\ldots A_{Br}^{(i_{n})}$ be a product
of the $3\times 3$ Brun matrices. Then, $B$ is primitive if and only if the
matrix $B^{(3)}$ appears in the product. Moreover, if $B$ is primitive then it
is Pisot.

This result was already known since the work of Brun [Br57].

The proofs of Theorem 1 and 2 are elementary and uses the following
inequality. Given a non-negative primitive $d\times d$ matrix $A$ and its
Perron-Frobenius eigenvector $v\in\mathbb{R}^{d}_{+}$ the absolute value
$\lambda_{2}$ of its second largest eigenvalue satisfies

$\lambda_{2}\leq\sup_{x\in v^{\perp}\backslash\\{0\\}}\frac{\|Ax\|}{\|x\|}.$

where $\|.\|$ is any norm on $\mathbb{R}^{d}$. In our proof, the information
we have on the localization of the Perron-Frobenius eigenvector comes from the
dynamical systems induced by the matrices; in other words the fully
subtractive and Brun continued fraction algorithms.

From a diophantine approximation point of view, the Pisot property is
particularly interesting because it provides the so called exponential
convergence of the continued fraction expansion for almost every vectors (see
[La93]). We show that the above results naturally extends to this situation in
Section 6.

Beyond continued fractions, Pisot matrices are of special interest in
substitutive dynamical systems. More precisely, replacing matrices with so
called substitutions, the associated dynamical systems admit $d-1$ eigenvalues
and in many cases the dynamical system can be proved to have purely discrete
spectrum (see [Py02] Chapter 7). This was the main motivation for the study of
the fully subtractive matrices in [ArIt01].

##  2 Fully subtractive and Brun continued fractions

Let us consider a finite or countable set $\mathcal{A}$ that we call alphabet
and for each $i\in\mathcal{A}$ a matrix
$A^{(i)}\in\operatorname{SL}(d,\mathbb{Z})$. We already saw two examples of
this with the fully subtractive algorithm where
$\mathcal{A}_{FS,d}=\\{1,2,\ldots,d\\}$ and for Brun algorithm where
$\mathcal{A}_{Br}=\\{1,2,3\\}$.

To the data $(\mathcal{A},(A^{(i)})_{i\in\mathcal{A}})$ we associate the set
of infinite words $\Delta=\mathcal{A}^{\mathbb{N}}$, the shift map
$T:\Delta\rightarrow\Delta$ and a cocycle

$\forall x\in\Delta,\forall n\geq 0,\quad
A_{n}(x)=A^{(x_{0})}A^{(x_{1})}\ldots A^{(x_{n-1})}.$

The maps $A_{n}:\Delta\rightarrow\operatorname{SL}(d,\mathbb{Z})$ satisfy the
so called _cocycle property_ : $A_{m+n}(x)=A_{m}(x)A_{n}(T^{m}x)$.

######  Definition 3.

Let $(A^{(i)})_{i\in\mathcal{A}}$ be a set of matrices in
$\operatorname{SL}(d,\mathbb{Z})$ where $\mathcal{A}$ is a finite or countable
alphabet. We say that a set $D\subset\mathbb{P}(\mathbb{R}^{d})$ is _adapted_
to these matrices if it is non-empty, it is the closure of its interior and
for all $i\in\mathcal{A}$ we have $A^{(i)}D\subset D$.

For example $PP(\mathbb{R}^{d})$ is always adapted. But we will be interested
in the somewhat smallest adapted set in order to localize the dominant
eigenvector.

Given $(A^{(i)})_{i\in\mathcal{A}}$ and $D\subset\mathbb{P}(\mathbb{R}^{d})$
adapted, we define $D^{(i)}=A^{(i)}D$ and more generally for a finite word
$w=i_{0}i_{1}\ldots i_{n-1}$ we define $D^{(w)}=A^{(i_{0})}A^{(i_{1})}\ldots
A^{(i_{n-1})}D$. Note that for any $w$ the set $D^{(w)}$ is _not_ empty. Given
an infinite word $x=x_{0}x_{1}\ldots\in\Delta$ we also set
$D_{n}(x)=A_{n}(x)D$ and $D_{\infty}(x)=\bigcap_{n\geq 0}D_{n}(x)$.

Let $A\in\operatorname{SL}(d,\mathbb{Z})$ and let $\lambda$ be its spectral
radius. Consider its Jordan decomposition over $\mathbb{C}$ and the Jordan
blocks associated with an eigenvalue of modulus $\lambda$ and being of maximal
dimension. To each of these maximal Jordan block is associated exactly one
eigenvector $v_{i}$. The _dominant eigenspace_ of $A$ is
$\mathbb{R}^{d}\cap(\mathbb{C}v_{1}\oplus\mathbb{C}v_{2}\oplus\ldots\mathbb{C}v_{r})$.
We have the following elementary result.

######  Lemma 4.

Let $(A^{(i)})_{i\in\mathcal{A}}$ be a finite or countable set of matrices in
$\operatorname{SL}(d,\mathbb{Z})$ and let $D$ be adapted. Let
$x=x_{0}x_{1}\ldots\in\Delta$ be an infinite word over $\mathcal{A}$. Then for
any $n$, the set $D_{n}(x)$ contains a basis of the dominant eigenspace of
$A_{n}(x)$.

We omit the proof that only uses the fact that the maximum growth of
$\|A^{n}u\|$ is $\lambda^{n}n^{k}$ where $k$ is the maximal dimension of a
Jordan block associated with an eigenvalue of maximal modulus of $A$.

Let

$D_{FS,d}=\\{(x_{1},\ldots,x_{d})\in\mathbb{P}(\mathbb{R}_{+}^{d}):\forall
i,j,k\quad x_{i}<x_{j}+x_{k}\\}$

and

$D_{Br}=\\{(x,y,z)\in\mathbb{P}(\mathbb{R}_{+}^{3}):\ x>y>z\\}.$

Then it is easily seen that $D_{FS,d}$ is adapted for the fully subtractive
matrices in dimension $d$ and $D_{Br}$ is adapted for the Brun matrices. In
figures 1 and 2 one can see the projective picture of the domains $D^{(1)}$,
$D^{(2)}$ and $D^{(3)}$. Note that in these cases, the domains $D^{(i)}$ are
disjoint but that it is not a requirement in our definition. Moreover, one can
see that in the Brun case the $D^{(i)}$ form a partition while it is not the
case for the fully subtractive.

$\displaystyle A_{FS,3}^{(1)}=\begin{pmatrix}1&0&0\\\ 1&1&0\\\
1&0&1\end{pmatrix}\qquad A_{FS,3}^{(2)}=\begin{pmatrix}1&1&0\\\ 0&1&0\\\
0&1&1\end{pmatrix}$ $\displaystyle A_{FS,3}^{(3)}=\begin{pmatrix}1&0&1\\\
0&1&1\\\ 0&0&1\end{pmatrix}.$

Figure 1: Fully subtractive partition of the domains with $d=3$.

$\displaystyle A_{Br}^{(1)}=\begin{pmatrix}1&1&0\\\ 0&1&0\\\
0&0&1\end{pmatrix}\qquad A_{Br}^{(2)}=\begin{pmatrix}1&1&0\\\ 1&0&0\\\
0&0&1\end{pmatrix}$ $\displaystyle A_{Br}^{(3)}=\begin{pmatrix}1&0&1\\\
1&0&0\\\ 0&1&0\end{pmatrix}.$

Figure 2: The matrices and domains for the Brun algorithm.

If the $D^{(i)}$ are disjoint one can define a continued fraction algorithm as
follows. One defines a partial map $f:D\dashrightarrow D$ by setting
$fx=(A^{(i)})^{-1}x$ on $D^{(i)}$.

One can compute that for Brun one has

$f_{Br}(x)=\operatorname{sort}(x-y,y,z)$

where $\operatorname{sort}:\mathbb{P}(\mathbb{R}_{+}^{3})\rightarrow D$ is the
map which permutes the coordinates in order to sort them. While for the fully
subtractive one has

$f_{FS,d}(x)=(x_{1}-x_{i},x_{2}-x_{i},\ldots,x_{i-1}-x_{i},x_{i},x_{i+1}-x_{i},\ldots,x_{d}-x_{i})\qquad\text{if
$x_{i}=\min(x_{1},\ldots,x_{d})$}.$

##  3 Strategy

The proofs of Theorems 1 and 2 follow a general strategy that we describe now.
We let $\|.\|$ be the $L^{\infty}$ norm on $\mathbb{R}^{d}$ and the associated
operator norm on matrices. That is for a vector $v$ and a matrix $A$

$\|v\|=\max(|v_{1}|,\ldots,|v_{d}|)\quad\text{and}\quad\|A\|=\max_{i=1,\ldots,d}\sum_{j=1,\ldots,d}|a_{ij}|.$

In this section the norm used on $\mathbb{R}^{d}$ has no importance. But it
turns out that, to apply the results to continued fraction algorithms, the
most convenient one was always the $L^{\infty}$ norm.

To a non-zero vector $v$ in $\mathbb{R}^{d}$, we associate its dual hyperplane
$H_{v}=\\{z\in\mathbb{R}^{d};\langle v,z\rangle=0\\}$. Given a non zero vector
$v$ in $\mathbb{R}^{d}$ we define the following semi-norm on $d\times d$
matrices

$\|B\|_{v}=\sup_{z\in
H_{v}\backslash\\{0\\}}\frac{\|Bz\|}{\|z\|}=\max_{\stackrel{{\scriptstyle\|z\|\leq
1}}{{z\in H_{v}}}}\|Bz\|.$

More generally, if $\Lambda\subset\mathbb{R}^{d}$ is a cone, we define

$\|B\|_{\Lambda}=\sup_{v\in\mathbb{P}(\Lambda)}\|B\|_{v}.$

Let $(A^{(i)})_{i\in\mathcal{A}}$ be a finite or countable set of matrices as
in Section 2. Let also $\Delta=\mathcal{A}^{\mathbb{N}}$ and
$D\subset\mathbb{P}(\mathbb{R}^{d})$ be adapted. Recall that
$D_{n}(x)=A_{n}(x)D=D^{(x_{0}x_{1}\ldots x_{n-1})}$ (in particular,
$D_{0}(x)=D$ and $D_{1}(x)=D^{(x_{0})}$) and that
$D_{\infty}(x)=\bigcap_{n\geq 0}D_{n}(x)$.

######  Lemma 5.

Let $(A^{(i)})_{i\in\mathcal{A}}$ be a finite or countable set of matrices in
$\operatorname{SL}(d,\mathbb{Z})$ and let $D\subset\mathbb{P}(\mathbb{R}^{d})$
be adapted. Let $B^{(i)}$ (respectively $B_{n}(x)$) denote the transposed of
$A^{(i)}$ (resp. $A_{n}(x)$). If for all $i\in\mathcal{A}$ we have

$\left\|B^{(i)}\right\|_{D^{(i)}}\leq 1.$

Then for any point $x=x_{0}x_{1}\ldots\in\Delta$ we have

$\|B_{n}(x)\|_{D_{n}(x)}\leq 1.$

In particular, if $x=(x_{0}x_{1}\ldots x_{p-1})^{\infty}$ is periodic the
matrix $A_{p}(x)=A^{(x_{0})}A^{(x_{1})}\ldots A^{(x_{p-1})}$ has at most one
eigenvalue greater than one in absolute value.

###### Proof.

Let $B$ be the transposed cocycle. The hypothesis is just the case $n=1$.
Assume that this inequality holds for $n$. By definition
$A_{n+1}(x)=A_{1}(x)A_{n}(Tx)$ and $D_{n+1}(x)=A_{1}(x)D_{n}(Tx)$. Hence $v\in
D_{n}(Tx)$ if and only if $A_{1}(x)v\in D_{n+1}(x)$. Let us choose $v\in
D_{n}(Tx)$, then

$\|B_{n+1}(x)\|_{A_{1}(x)v}=\|B_{n}(Tx)B_{1}(x)\|_{A_{1}(x)v}\leq\|B_{n}(Tx)\|_{v}\cdot\|B_{1}(x)\|_{A_{1}(x)v}.$

Now, let $x=(x_{0}x_{1}\ldots x_{p-1})^{\infty}$ be a periodic point. Because
$D_{p-1}(x)\supset D_{\infty}(x)$, the matrix $B_{p}(x)$ satisifies

$\|B_{p}(x)\|_{D_{\infty}(x)}\leq 1.$

By Lemma 4 a basis of the dominant eigenspace of $B_{p}$ belongs to
$D_{\infty}(x)$. Consequently, the union of the orthogonals of $D_{\infty}(x)$
contains all eigenspaces corresponding to the non-dominant eigenvalues. From
the above inequality, we deduce that the absolute value of all eigenvalues
different from the first one are bounded by $1$. ∎

##  4 Pisot property for Arnoux-Rauzy matrices

In this section we prove Theorem 1. Let $d\geq 2$ be an integer and let
$e_{1}$, $e_{2}$, …$e_{d}$ be the canonical basis of $\mathbb{R}^{d}$. Let
$e=e_{1}+e_{2}+\ldots+e_{d}$ and for $i=1,\ldots,k$ let $f_{i}=e-e_{i}$. The
domain $D_{FS,d}$ is the convex hull of the rays vectors
$\mathbb{R}_{+}f_{i}$.

Let as usual $B^{(w)}=(A^{(w)})^{*}$ and $B_{n}(x)=(A_{n}(x))^{*}$. We claim
that we even have a stronger property than what is required in Lemma 5

$\forall i=1,\ldots,d,\qquad\|B^{(i)}\|_{D}\leq 1.$

Let us prove this claim. Let $v\in D_{\infty}$, then we may write
$v=\mu_{1}f_{1}+\mu_{2}f_{2}+\ldots+\mu_{d}f_{d}$ for some non-negative
numbers $\mu_{i}$ that satisfy $\mu_{1}+\mu_{2}+\ldots+\mu_{d}=1$. We hence
have $v=e-\sum\mu_{i}e_{i}$ and

$H_{v}=\\{z\in\mathbb{R}^{d};\langle
z,v\rangle=0\\}=\\{z\in\mathbb{R}^{d};\sum z_{j}=\sum\mu_{j}z_{j}\\}.$

Given $z\in H_{v}$ we have
$B^{(i)}z=(z_{1},\ldots,z_{i-1},\sum\mu_{z}z_{j},z_{i+1},\ldots,z_{d})$, In
other words, $B^{(i)}$ acts on $H_{v}$ as a stochastic matrix $P(i,v)$ which
is the identity except its $i$-th row which is
$(\mu_{1},\mu_{2},\ldots,\mu_{d})$. In particular, $\|B^{(i)}\|_{v}\leq 1$.

Now for a given finite product $A=A^{(i_{0})}A^{(i_{1})}\ldots A^{(i_{p-1})}$
if one of the letter $\\{1,\ldots,d\\}$ is missing in the sequence
$(i_{0},i_{1},\ldots,i_{p-1})$ then $Ae_{i}=e_{i}$ and so the matrix is not
primitive. On the other hand, if all letters appear it is easy to see that all
entries in $A$ are positive.

Now let $x=(x_{0}x_{1}\ldots x_{p-1})^{\infty}$ be a periodic point that
contains all letters from $\mathcal{A}$. Because of positivity, all the orbits
is contained in the interior of $D$ and the dominant eigenvalue is simple. Let
$v_{0}\in D$ be a dominant eigenvector and let $v_{n}=A_{n}(x)v_{0}$. Because
all $(v_{n})$ belongs to the interior of $D$ the coefficients $\mu_{1}$,
$\mu_{2}$, …, $\mu_{d}$ that appear in the stochastic matrices
$P(x_{i},v_{i})$ are all positive. Now, the product $P(x_{p-1},v_{p-1})\ldots
P(x_{1},v_{1})P(x_{0},v_{0})$ is a stochastic matrix with all its entries
positive. Hence, its second eigenvalue, which is also the second eigenvalue of
$A_{p}(x)$, is less than $1$ in absolute value.

##  5 Pisot property for Brun algorithm (in dimension 3)

We now turn to the proof of Theorem 2. Let $\mathcal{A}=\\{1,2,3\\}$ and
$A^{(1)}$, $A^{(2)}$, $A^{(3)}$ be the matrix of the Brun algorithm. We let
$\Delta=\mathcal{A}^{\mathbb{N}}$ and denote by $A$ and $B$ respectively the
cocycle and the transposed cocycle. We claim that, as in the case of the fully
subtractive, we have the stronger property that

$\forall i\in\\{1,2,3\\},\quad\|B^{(i)}\|_{D^{(i)}}\leq 1.$

We only need to consider the matrix $\displaystyle
B^{(1)}=\begin{pmatrix}1&0&0\\\ 1&1&0\\\ 0&0&1\end{pmatrix}$ since the other
two are obtained by multiplying by a permutation matrix which will not change
the $L^{\infty}$-norm.

Let $v=\mu_{1}(1:0:0)+\mu_{2}(1:1:0)+\mu_{3}(1:1:1)\in D$ for some
$\mu_{1},\mu_{2},\mu_{3}$ such that $\mu_{1}+\mu_{2}+\mu_{3}=1$ and
$H_{v}=\\{z\in\mathbb{R}^{d};z_{1}+z_{2}=\mu_{1}z_{2}-\mu_{3}z_{3}\\}$. Now,
for any $z\in H_{v}$ we have
$B^{(1)}(z_{1},z_{2},z_{3})=(z_{1},z_{1}+z_{2},z_{3})=(z_{1},\mu_{1}z_{2}-\mu_{3}z_{3},z_{3})$.
In other words
$\|B^{(1)}(z_{1},z_{2},z_{3})\|_{1}\leq\|(z_{1},z_{2},z_{3})\|$.

Now, given a product $A=A_{B}r^{(i_{1})}\ldots A_{B}r^{(i_{n})}$ it is easy to
see that if $3$ does not appear in the sequence $(i_{1},i_{2},\ldots,i_{n})$
then $Ae_{3}=e_{3}$ and hence the matrix $A$ can not be irreducible.
Conversly, if $3$ appears then $A^{3}$ is easily seen to be positive.

Now consider the matrix $P$ built from the begining of the proof. As in the
case of the fully subtractive algorithm for a primitive product we got that
the $\mu_{i}$ are all positive. Given a product $A$ where the matrix
$A_{Br}^{(3)}$ appears, the matrix $A^{3}$ is then such that all rows are such
that sum of their absolute values is strictly less than one. In other words
$\|A^{3}\|<1$.

##  6 Lyapunov exponents

Let $(A^{(i)})_{i\in\mathcal{A}}$ be a finite or countable set of matrices.
Let $\Delta$, $T$, $A$, $B$ denote as before the infinite words, the shift map
the cocycle and the transposed cocycle. Let also $D$ be adapated to these
matrices.

The asymptotic of the cocycle (or the transposed cocycle) are studied through
_Lyapunov exponents_. Given a $T$-invariant ergodic probability measure $\mu$
on $\Delta$, we associate the real numbers
$\gamma^{\mu}_{1}\geq\gamma^{\mu}_{2}\geq\ldots\geq\gamma^{\mu}_{d}$ defined
by

$\forall
k\in\\{1,2,\ldots,d\\},\quad\gamma_{1}+\gamma_{2}+\ldots+\gamma_{k}=\lim_{n\to\infty}\int_{\Delta}\frac{\log\|\wedge^{k}A_{n}(x)\|}{n}d\mu(x).$

In order to be well defined we assume that

$\int_{\Delta}\max\left(\log\|A_{1}(x)\|,\log\|A_{1}(x)^{-1}\|\right)d\mu(x)<\infty$
(1)

and we refer to this condition as the _$\log$ -integrability_ of the cocycle.
If the alphabet $\mathcal{A}$ is finite the cocycle is automatically
$\log$-integrable. If $x$ is a periodic point of $T$ and
$\mu=(\delta_{x}+\delta_{Tx}+\ldots+\delta_{T^{n-1}x})/n$ is the sum of Dirac
masses distribued along its orbit, then the associated Lyapunov exponents are
the logarithms of the absolute values of eigenvalues of $A_{n}(x)$ where $n$
is the period of $x$. In that sense, Lyapunov exponents generalize
eigenvalues.

Given a measure $\mu$ for which the cocycle is $\log$-integrable, we say that
$(\Delta,T,A,\mu)$ has _Pisot spectrum_ if the associated Lyapunov exponents
satisfy $\gamma_{1}^{\mu}>0>\gamma_{2}^{\mu}$. This property is related to the
strong convergence of higher dimensional continued fraction algorithm [La93].

Now we restate Lemma 5 in a more dynamical context.

######  Lemma 6.

Let $(A^{(i)})_{i\in\mathcal{A}}$ be a finite or countable set of non-negative
matrices in $\operatorname{SL}(d,\mathbb{Z})$. Let $(\Delta,T,A,B)$ be the
associated full shift with its cocycle and its transposed cocycle. Let also
$D$ be adapated. Assume that

$\forall
i\in\mathcal{A},\quad\left\|\left(A^{(i)}\right)\right\|_{D^{(i)}}\leq 1.$

Let $\mu$ be a $T$-invariant and ergodic measure on $D$ so that

  * •

the cocycle $A_{n}$ is $\log$-integrable,

  * •

there exists a cylinder $[w]$ such that $\mu([w])>0$, $A^{(w)}$ is positive
and $\|A^{(w)}\|_{D^{(w)}}<1$.

Then two first Lyapunov exponents of the cocycle $A_{n}$ for the measure $\mu$
satisfies $\gamma_{1}^{\mu}>0>\gamma_{2}^{\mu}$.

###### Proof.

Let us first prove that $\gamma_{1}>0$.

Now, by definition, for $\mu$-almost every $x$

$\gamma_{1}=\lim_{n\to\infty}\frac{\log\|A_{n}(x)\|}{n}$

Let $m=|w|$ be the length of $w$ and consider the position which are multiple
of $m$. For a $\mu$-generic $x$ we have by Birkhoff theorem that

$\lim_{n\to\infty}\frac{\\#\\{i\leq n:\ T^{i}(x)\in[w]\\}}{n}=\mu([w])$

In other words, given sequence of length $n$ large enough we can find a linear
number of disjoint occurrences of $w$ (up to a sublinear error). Let $k_{n}$
be the number of these occurrences, then necessarily each entry of $A_{n}(x)$
is larger than the corresponding one in $C^{k_{n}}$ where $C$ is the matrix
which contains a $1$ in every position. In particular $\gamma_{1}>0$.

From the existence of $w$ it also follows that for a $\mu$-generic $x$ the
cone $D_{\infty}(x)$ is reduced to a line contained in the interior of
$\mathbb{R}_{+}^{d}$. We can hence define $\mu$-almost everywhere a function
$v:\Delta\rightarrow\mathbb{R}_{+}^{d}$ by $D_{\infty}(x)=\mathbb{R}_{+}v(x)$
and $\|v(x)\|=1$. We then have the following formulas which holds for
$\mu$-almost every $x$

$\gamma_{1}=\lim_{n\to\infty}\frac{-\log\|A_{n}(x)^{-1}v(x)\|}{n}d\mu(x)\quad\text{and}\quad\gamma_{2}=\lim_{n\to\infty}\frac{\log\|B_{n}(x)\|_{v(x)}}{n}d\mu(x).$

It is then easy to derive the estimate for $\gamma_{2}$. The map $x\mapsto
v(x)$ and the dual hyperplanes $H_{v(x)}$ satisfy the following covariance
properties

$D_{\infty}(Tx)=A(x)^{-1}D_{\infty}(x)\quad\text{and}\quad
H_{A^{-1}v}=A^{*}H_{v}.$

Hence as in the proof of Lemma 5, we deduce that for all $v\in D_{\infty}(x)$

$\|B_{m+n}(x)\|_{v(x)}\leq\|B_{m}(x)\|_{v(x)}\|B_{n}(T^{m}x)\|_{v(T^{m}x)}$

In particular, for $\mu$-almost every $x\in[w]$, any $n\geq|w|$ we get that
$\|B_{n}(x)\|_{v(x)}<1$. Let $\delta=\|B^{(w)}\|_{D^{(w)}}<1$. Using the same
argument as in the estimation of $\gamma_{1}$ we get that

$\gamma_{2}\leq\liminf_{n\to\infty}\frac{k_{n}\log\delta}{n}.$

And the above limit is strictly negative. ∎

## References

  * [Ar88] P. Arnoux, Un exemple de semi conjugaisaon entre un échanges d’intervalles et une translation sur le tore, Bull. SMF 116 no 4 (1988), 489–500.
  * [ArIt01] P. Arnoux, S. Ito, Pisot substitutions and Rauzy fractals, Bull. Belg. Math. Soc. Simon Stevin 8 no 2 (2001), 181–207.
  * [ArRa91] P. Arnoux, G. Rauzy, Représentation géométriques des suites de complexité $2n+1$, Bull. SMF 119 no 2 (1991) 199–215.
  * [ArSt] P. Arnoux, Š. Starosta, The Rauzy gasket, Further developments in fractals and related fields, Trends Math. Birkhäuser/Springer, New York (2013).
  * [Br57] V. Brun, Algorithmes euclidiens pour trois et quatre nombres, Congr. Math. Scand. XIII, (1957), 45–64,
  * [Br81] A. Brentjes, Multidimensional continued fraction algorithms, Mathematical Centre Tracts, 145, Mathematisch Centrum, Amsterdam, (1981).
  * [KrMe95] C. Kraaikamp, R. Meester, Ergodic properties of a dynamical system arising from percolation theory, Ergod. Th. and Dynam. Sys., 15 (1995) 653–661.
  * [La93] J.C. Lagarias, The quality of the diophantine approximations found by the Jacobi-Perron algorithm and related algorithms, Mh. Math. 115 (1993) 299–328.
  * [Py02] N. Pytheas Fogg, Substitutions in Dynamics, Arithmetic and Combinatorics, Lecture Notes in Mathematics, (2002).


11institutetext: Department of Biomedical Engineering, King’s College London,
SE1 7EH, UK
11email:
<EMAIL_ADDRESS>WMG, University of Warwick, Coventry, CV4 71AL
22email<EMAIL_ADDRESS>
# V-FCNN: Volumetric Fully Convolution Neural Network For Automatic Atrial
Segmentation

Nicoló Savioli1 Giovanni Montana1,2 Pablo Lamata1

###### Abstract

Atrial Fibrillation (AF) is a common electro-physiological cardiac disorder
that causes changes in the anatomy of the atria. A better characterization of
these changes is desirable for the definition of clinical biomarkers,
furthermore, thus there is a need for its fully automatic segmentation from
clinical images. In this work, we present an architecture based on
3D-convolution kernels, a Volumetric Fully Convolution Neural Network
(V-FCNN), able to segment the entire volume in a one-shot, and consequently
integrate the implicit spatial redundancy present in high-resolution images. A
loss function based on the mixture of both Mean Square Error (MSE) and Dice
Loss (DL) is used, in an attempt to combine the ability to capture the bulk
shape as well as the reduction of local errors products by over-segmentation.
Results demonstrate a reasonable performance in the middle region of the atria
along with the impact of the challenges of capturing the variability of the
pulmonary veins or the identification of the valve plane that separates the
atria to the ventricle. A final dice of $92.5\%$ in $54$ patients ($4752$
atria test slices in total) is shown.

###### Keywords:

Cardiac Imaging, Segmentation, FCNN, Atria, Fibrillation, Shape, Clinical
Biomarkers, Anatomy, Deep Learning.

##  1 Introduction

Atrial Fibrillation (AF) is a common electro-physiological cardiac disorder
with a large prevalence worldwide [1] that causes changes in the anatomy of
the atria. A better characterization of these changes as a consequence of the
substrate that causes and sustains fibrillation is desirable for the
definition of clinical biomarkers. These biomarkers can be directly started
from the image (i.e. the shape of the atria [2], or the fibrotic burden from
late gadolinium enhanced (LGE) magnetic resonance imaging (MRI) [3]) or from
mechanistic simulations of the function (i.e. computation of the risk of
arrhythmia perpetuation [4]).

There is thus a need for a fully automatic segmentation of the atria from
clinical images, especially in LGE studies. The current state of the art is
based in tedious along with error-prone manual procedures, and the main
difficulty is the lack of contrast from atrial tissue including the
surrounding background. Fully automated solutions are desirable to speed up
the process and remove inter- and intra-observer variability. In this
direction, a combination of multi-atlas registration within 3D level-set has
been proposed, reporting a reasonable performance in the main atrial body and
pulmonary vein regions [5]. The large computational burden of this multi-atlas
approach can be alleviated by the use of convolutional neural networks (CNNs),
as it has been illustrated in [6] for the analysis of 2D MRI slices.

The first idea followed in this work is the use of 3D CNNs as the effective
solution for the segmentation of the 3D LGE datasets. The starting point is
the v-net [7], which takes into consideration the spatial redundancy naturally
present on the entire volumetric stack with 3D-kernels, showing good benefits
in different cardiac segmentation problems [8],[9]. This architecture is
modified to reduce the memory burden and speed-up its training. The last idea
explored in this work is the choice of a sensible loss function, the joint
combination of both the mean squared error (MSE) and Dice Loss, which has been
reported to be beneficial as the MSE minimizes global image details while the
Dice Loss reduces local over-segmentation errors [7].

Figure 1: V-FCNN architecture. Input is the (XYZ) 3D MRI volume of size
($127\times 127\times 88$), also passed through the down-sampling path (blue
arrow), represented by a 3D kernels Convolution Neural Network (CNN) able to
progressively reduce the input volume slices. Then, the hidden features, at
the end of it, are restored within 3D up-sampling kernels (red-arrow), ending
in an output being a 3D mask of size ($127\times 127\times 88$). Both down-
sampling and up-sampling paths consist of four 3D-convolutions blocks (blue
boxes) followed by PreLU plus 3D-Batch Normalisation (BN). The number of
feature maps for each convolution layers are $16,32,64,128$ both in down and
up-sampling.

##  2 Atrial Datasets

A population of 100 3D GE-MRIs, data and masks, were made available through
the [2018 Atrial Segmentation Challenge] and used in this work. Images had a
acquisition resolution of 0.625 x 0.625 x 0.625 $mm^{3}$.

##  3 Method

A volumetric fully convolutional network, V-FCNN, is designed with two main
paths (see Fig.1): the volumetric down-sampling path as well as the volumetric
up-sampling path. The volumetric down-sampling path has four 3-D convolutions
blocks, each following by PreLU along with 3-D Batch Normalisation (BN)
layers. This path takes as input the entire (XYZ) volume and progressively
reduces the size of each slice (XY) together with the number of stack spatial
slices (Z). In this phase, the volume is compressed and presents both (XY) and
(Z) reduction. In a complementary fashion, the volumetric up-sampling path
restores the compressed volume to its initial size, with every 3-D up-sampling
convolution blocks being followed by PreLU and 3-D Batch Normalisation (BN).
Each sampling paths of V-FCNN contains four blocks within 16,32,64,128 -3D
kernels respectively (with the size fix to $3\times 3$).

Image down-sampling during a segmentation task presents the problem of feature
map reduction, followed by a strong spatial information loss. This problem has
been addressed in v-net [7], [8] by adding skips layers between down-sampling
and up-sampling layers in order to fuse low-level features within high level
features.

Our work explores an alternative approach to the skip path connections. In
particular, we boost our model to capture fine details in two ways. First, the
use of max-pooling operations is avoided in order to prevent the loss of
spatial resolution (i.e if pools do not overlap well, pooling operation loses
appreciable information where the objects are located in the image). The
second idea is to use a loss layer that combines the $MSE_{loss}$ and
$DICE_{loss}$ metrics as an attempt to recover details in the image. The
rationale is that the $DICE_{loss}$ term searches for local details in the
volume data, as a consequence, the $MSE_{loss}$ can be seen as a
regularization that instead focuses on global features on the MRI volume. This
bimodal loss also prevents the V-FCNN to fall in a local minimum; especially
within small atrial regions. Specifically, the loss is defined as:

$\begin{split}Loss=MSE_{loss}+\lambda\cdot
DICE_{loss}=\frac{1}{Z}\sum_{s=0}^{Z}(y[s]-\hat{y}[s])^{2}\\\
+\lambda\cdot\frac{1}{Z}\sum_{s=0}^{Z}\frac{2\sum_{i}^{N}y[s]_{i}\hat{y}[s]_{i}}{\sum_{i}^{N}y[s]_{i}+\sum_{i}^{N}\hat{y}[s]_{i}},\end{split}$
(1)

where $y[s]$ are the Ground Truth (GT) binary slices, as well as $\hat{y}[s]$
are the correspondent’s prediction masks. $s$ is an index through the spatial
slices in the Z (z-axis). Then, the sums for $i$ in $DICE_{loss}$ runs over
all N pixels of the prediction masks $\hat{y}[s]$ and the $y[s]$ GT masks.
Finally, $\lambda$ controls the amount of $DICE_{loss}$ during the
optimisation training process (set to $1e-3$).

The size of each MRI slice is reduced to $127\times 127$ pixels using down-
sampling bi-cubic interpolation, allowing the network to be faster in training
stage moreover enabling a solution without expensive GPU hardware (i.e. large
images need more GPU global memory $\sim 32$ GB). An up-sampling bi-cubic
interpolation is finally used for restoring the mask to the size of $640\times
640$ pixels.

Figure 2: Visual comparison of the segmentations obtained from V-FCNN (green
line) vs clinical ground truth (red line) in three different test patients
(number 1, 2 or 4). The comparison is made at three different sections of the
atrium: top, middle and bottom. Note how the V-FCNN is able to segment not
only visually simpler slices (middle section) but also more complex cases (top
and bottom sections).

###  3.1 Implementation details

For our experiments, we train the network with a number of epochs of 1000 up
to convergence. The Stochastic Gradient Descent (SDG) was used with a learning
rate of $1e-4$, while the momentum and weight decay are $0.9$, $1e-5$
respectively. Furthermore, to increase the generalization of the network, data
argumentation was used, finding particularly effective the random vertical and
horizontal flip in close combination with plane translation. Input sequences
were equalized in grayscale intensity with CLAHE (Contrast Limited Adaptive
Histogram Equalization) [10], and the noise was minimized with a combination
of High-Pass Filters and Gaussian blurring filters.

Figure 3: Visual comparison between ground truth (red) compared with those
obtained by proposed V-FCNN (green). Note that, the mesh coarse resolution is
related to the low number of triangles used.

##  4 Experiments

The experimentation phase is divided into two phases. In the first phase
(preparation phase) 5 patients are used for validation ($440$ slices in
total), and the rest (90) is used for training. In the second phase
(competition), our algorithm was evaluated in the test-set with 54 patients
($4752$ atria slices in total).

Segmentation accuracy was measured with the Dice Metric (DM) [11], which was
subdivided into 3 regions of the atria: top (including the pulmonary veins),
middle and bottom (including the valve plane that divides the atria and
ventricle). Besides, DM and the surface Hausdorff Distance (HD) [12] are
computed for the entire atrial anatomy.

Particularly, given a point $p$ and a surface $S$, the distance
$\upvarepsilon(p.S)$ is defined as:

$\upvarepsilon(p,S)=min_{p^{{}^{\prime}}\in S}d(p,p^{{}^{\prime}})$ (2)

whereas the $d(\cdot)$ is the Euclidean distance between two points in a
Euclidean space. Then, the HD between mesh surfaces $S_{GroundTruth}$ and
$S_{vfcnn}$ is given as:

$HD(S_{GroundTruth},S_{vfcnn})=max_{p\in
S_{GroundTruth}}\upvarepsilon(p,S_{vfcnn})$ (3)

Proposed V-FCNN achieved in the first phase an average DM of $69.6\pm 16.1$,
$82.1\pm 1.4$ and $78.0\pm 6.0$ in the top, middle and bottom regions,
respectively (see Table 1 and an illustrative example in Fig. 2). The HD
ranged from 0.31 to 0.86 mm. The average DM and HD were $76.58(7.87)$ and
$0.59$ respectively. The DM in the competition phase was $92.5\%$.

Visual inspection of the contours reveals how a loss of accuracy is seen on
the first top slices for two of the cases, creating an artificial flattening
of the shape (see patients 1 and 2 of Fig. 3), and on the top of the atria
with the more variable anatomy of the pulmonary veins.

| TOP | MID | BOTTOM | 3D-MESH
---|---|---|---|---
Patient number | DM (%) | DM (%) | DM (%) | HD
1 | 77.74 (8.46) | 84.62 (2.40) | 76.99 (10.89) | 0.86
2 | 72.26 (9.50) | 81.25 (1.57) | 54.18 (33.70) | 0.66
3 | 74.10 (8.92) | 77.75 (1.36) | 68.27 (14.75) | 0.31
4 | 81.60 (0.74) | 81.54 (1.20) | 72.32 (14.32) | 0.55
5 | 84.48 (2.58) | 85.36 (0.66) | 76.33 (7.05) | 0.58

Table 1: Automatic segmentation results (Dice Metric and Hausdorff Distance)
for all five test patients. Results report the mean and standard deviation,
and are divided into three different atrium sections: top, middle and bottom.
Figure 4: Exemplary result in 3 cases of the competition cases, illustrated
with the 3D reconstruction (red, top) and 3 different sections of atrium
segmentation: top, middle, and bottom.

##  5 Discussion and Conclusions

The exploitation of spatial coherence across different volumetric slices is an
important resource for improving the accuracy of fully automated segmentation.
Proposed V-FCNN achieves a good segmentation performance, mostly in the middle
atrial section. Limitations occur in the top and bottom sections, caused by
the presence of the pulmonary veins and the difficulty to identify where the
atria and ventricle split. In the preparation phase the V-FCNN achieved a DM
between $82.05(3.43)$ (top case) and $69.23(14.92)$ (worst case), along with a
DM of $92.5$ at the competition phase.

A cardiac imaging study has a huge level of correlation in space, as well as
many segmentation algorithms have exploited the spatial coherence within
different deep-learning techniques [7], [8], [9], [13], [14]. Some of them
have used 3D convolutional kernels [7], [8], [9] that have the advantage of
directly capturing the spatial information in each convolutional layer without
adding extra parameters (i.e. the addition of a recurring network). This is
the main reason that has driven us to use them in our optimized V-FCNN.

The limitation of the 3D convolutional kernels is the large memory burden, and
this is the reason why the original detail of the images was reduced (from 640
to 127 pixels). However, this image reduction is required due to the current
hardware limitations available. In fact, while the 3D kernels have a good
ability to extract the spatial correlation of the MRI sequence, on the other
they need a high capacity in GPU global memory size. Thus, an excessive
reduction of the MRI to adapt to the available hardware causes an inevitable
loss in segmentation details (i.e probable motivation of $92.5\%$ dice).
Therefore, we recommend future tests with 3D kernel at more higher MRI
resolution with the help of NVIDIA V100 GPUs [15] that presents $16-32$ GB of
global memory.

Furthermore, the concatenation of two CNNs, the second working at the full
resolution cropping the region of interests containing the atria, is a good
strategy that will minimize the impact of this choice. An alternative is to
work with the information at two levels of detail, which has shown to achieve
a better performance [16].

Proposed V-FCNN simplifies the v-net [7] by removing the skip-paths, in an
attempt to achieve a much quicker training convergence. The constitutive
advantage of the skip-paths is to increase the localization, but this comes at
the cost of a considerable slow down of the speed of the network (i.e GPU
global memory leak) while propagating the gradient, at every iteration,
forward and back through those paths [17].

The choice of the loss function is an important consideration in any CNN, and
the joint minimization of MSE and Dice Loss has been adopted in our solution.
The interpretation is that the MSE looks at global volumetric features, while
the Dice Loss (DL) regularize it trying to fit local details. This choice
allowed the learning to avoid local minima and the corresponding slow
convergence of using only MSE. The optimal weight between them, and the
inclusion of further criteria such as an L1/L2 loss or a statistical distance
to a library of existing cases, should be addressed in future extensions of
this work.

An alternative strategy to 3D kernels is the use of recurrent units [13],[14],
where recurrence is used to capture the redundancy between adjacent slices.
The use of recurrence is a much more memory efficient approach (i.e fewer
parameters to capture recurring partners), also reporting improvements at the
challenging apical slices of the left ventricle [13]. But these solutions are
limited by the vanishing gradient that unfortunately occurs within long
sequences, what can be partially avoided by imposing upper bound constraints
on the backward gradient (ie. gradient clipping) or with a regularisation term
[18]. Future work is thus still needed to maximize the synergies between 3D
kernels and recurrent units within a more higher input MRI resolution.

## Acknowledgements

This work was supported by the Wellcome/EPSRC Centre for Medical Engineering
at King’s College London [g.a. 203148/Z/16/Z]. PL holds a Wellcome Trust
Senior Research Fellowship [g.a. 209450/Z/17/Z].

## References

  * [1] EN Prystowsky, DW Benson, and V Fuster.  Management of patients with atrial fibrillation?: A statement for healthcare professionals from the subcommittee on electrocardiography and electrophysiology.  American Heart Association. Circulation, 93:15, 1996.
  * [2] Marta Varela, Felipe Bisbal, Ernesto Zacur, Antonio Berruezo, Oleg V. Aslanidi, Lluis Mont, and Pablo Lamata.  Novel computational analysis of left atrial anatomy improves prediction of atrial fibrillation recurrence after ablation.  Frontiers in Physiology, 8:68, 2017.
  * [3] RJ Kim, E Wu, A Rafael, and et al.  The use of contrast‐enhanced magnetic resonance imaging to identify reversible myocardial dysfunction.  N Engl J Med, 343:1445–1453, 2000.
  * [4] Patrick M. Boyle, Sohail Zahid, and Natalia A. Trayanova.  Towards personalized computational modelling of the fibrotic substrate for atrial arrhythmia.  EP Europace, 18:136–145, 2016.
  * [5] Qian Tao, Rahil Shahzad, Esra Gucuk Ipek, Floris F. Berendsen, Saman Nazarian, and Rob J. van der Geest.  Fully automated segmentation of left atrium and pulmonary veins in late gadolinium enhanced mri.  Journal of Cardiovascular Magnetic Resonance, 18(1):O84, Jan 2016\.
  * [6] Aliasghar Mortazi, Rashed Karim, Kawal S. Rhode, Jeremy Burt, and Ulas Bagci.  Cardiacnet: Segmentation of left atrium and proximal pulmonary veins from MRI using multi-view CNN.  CoRR, abs/1705.06333, 2017.
  * [7] Fausto Milletari, Nassir Navab, and Seyed-Ahmad Ahmadi.  V-net: Fully convolutional neural networks for volumetric medical image segmentation.  CoRR, abs/1606.04797, 2016.
  * [8] Fabian Isensee, Paul Jaeger, Peter M. Full, Ivo Wolf, Sandy Engelhardt, and Klaus H. Maier-Hein.  Automatic cardiac disease assessment on cine-mri via time-series segmentation and domain specific features.  CoRR, abs/1707.00587, 2017.
  * [9] Hinrich B. Winther, Christian Hundt, Bertil Schmidt, Christoph Czerner, Johann Bauersachs, Frank K. Wacker, and Jens Vogel-Claussen.  $$\nu$$-net: Deep learning for generalized biventricular cardiac mass and function parameters.  CoRR, abs/1706.04397, 2017.
  * [10] Hardeep Kaur and J. Sharmila Rani.  Mri brain image enhancement using histogram equalization techniques.  2016 International Conference on Wireless Communications, Signal Processing and Networking (WiSPNET), pages 770–773, 2016.
  * [11] M.R. Avendi, Arash Kheradvar, and Hamid Jafarkhani.  A combined deep-learning and deformable-model approach to fully automatic segmentation of the left ventricle in cardiac mri.  Medical Image Analysis, 30:108 – 119, 2016.
  * [12] Paolo Cignoni, Claudio Rocchini, and Roberto Scopigno.  Metro: Measuring error on simplified surfaces.  CoRR, 1996.
  * [13] Rudra P. K. Poudel, Pablo Lamata, and Giovanni Montana.  Recurrent fully convolutional neural networks for multi-slice mri cardiac segmentation.  In Maria A. Zuluaga, Kanwal Bhatia, Bernhard Kainz, Mehdi H. Moghari, and Danielle F. Pace, editors, Reconstruction, Segmentation, and Analysis of Medical Images, pages 83–94, Cham, 2017. Springer International Publishing.
  * [14] Jianxu Chen, Lin Yang, Yizhe Zhang, Mark S. Alber, and Danny Z. Chen.  Combining fully convolutional and recurrent neural networks for 3d biomedical image segmentation.  CoRR, abs/1609.01006, 2016.
  * [15] Zhe Jia, Marco Maggioni, Benjamin Staiger, and Daniele Paolo Scarpazza.  Dissecting the NVIDIA volta GPU architecture via microbenchmarking.  CoRR, abs/1804.06826, 2018.
  * [16] Z. Xiong, V. V. Fedorov, X. Fu, E. Cheng, R. Macleod, and J. Zhao.  Fully automatic left atrium segmentation from late gadolinium enhanced magnetic resonance imaging using a dual fully convolutional neural network.  IEEE Transactions on Medical Imaging, pages 1–1, 2018.
  * [17] Adam Paszke, Abhishek Chaurasia, Sangpil Kim, and Eugenio Culurciello.  Enet: A deep neural network architecture for real-time semantic segmentation.  CoRR, abs/1606.02147, 2016.
  * [18] Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio.  Understanding the exploding gradient problem.  CoRR, abs/1211.5063, 2012.