# Pediatric Otoscopy Video Screening with Shift Contrastive Anomaly Detection

Weiyao Wang, Aniruddha Tamhane, Christine Santos, John R. Rzasa, James H.
Clark, Therese L. Canares and Mathias Unberath This work was funded in part by
the Leon Lowenstein Foundation.W. Wang, A. Tamane and M. Unberath are with the
Johns Hopkins University School of Engineering, Baltimore, MD 21218, USA
(e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>unberath@jhu.edu).J. Rzasa is
with the Robert E. Fischell Institute for Biomedical Devices at the University
of Maryland, College Park 20742, USA (e-mail: rzasa@umd.edu).C. Santos, J.
Clark and T. Canares are with Johns Hopkins University School of Medicine,
Baltimore, MD 21205, USA (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>therese.canares@jhmi.edu).

###### Abstract

Ear related concerns and symptoms represents the leading indication for
seeking pediatric healthcare attention. Despite the high incidence of such
encounters, the diagnostic process of commonly encountered disease of the
middle and external presents significant challenge. Much of this challenge
stems from the lack of cost effective diagnostic testing, which necessitating
the presence or absence of ear pathology to be determined clinically. Research
has however demonstrated considerable variation among clinicians in their
ability to accurately diagnose and consequently manage ear pathology. With
recent advances in computer vision and machine learning, there is an
increasing interest in helping clinicians to accurately diagnose middle and
external ear pathology with computer-aided systems. It has been shown that AI
has the capacity to analyse a single clinical image captured during
examination of the ear canal and eardrum from which it can determine the
likelihood of a pathognomonic pattern for a specific diagnosis being present.
The capture of such an image can however be challenging especially to
inexperienced clinicians. To help mitigate this technical challenge we have
developed and tested a method using video sequences. The videos were collected
using a commercially available otoscope smartphone attachment in an urban,
tertiary-care pediatric emergency department. We present a two stage method
that first, identifies valid frames by detecting and extracting ear drum
patches from the video sequence, and second, performs the proposed shift
contrastive anomaly detection to flag the otoscopy video sequences as normal
or abnormal. Our method achieves an AUROC of 88.0% on the patient-level and
also outperforms the average of a group of 25 clinicians in a comparative
study, which is the largest of such published to date. We conclude that the
presented method achieves a promising first step towards automated analysis of
otoscopy video.

###### Index Terms:

Anomaly detection, otoscope, self-supervised representation learning, deep
learning.

##  I Introduction

Ear related concerns constitutes the leading cause for seeking pediatric
healthcare attention in the USA. Given the current lack of cost-effective
confirmatory testing, accurate diagnosis and subsequent management depend on
visual detection of characteristic findings during otoscope examination.
Despite the frequency of such encounters the medical literature suggest that
the diagnostic accuracy for such pathology is only around 46-56%. In alignment
with the bias towards over diagnosis of ear disease, it is currently estimated
that between 25-50% of all antibiotics prescribed for ear disease are not
indicated [1, 2, 3]. Beyond risking unnecessary medical complications and the
downstream unintended consequence of potential antibiotic resistance, over-
diagnosis of ear disease adds an estimated $59 million in unnecessary
healthcare spending in the US per annum [4]. Computer-aided diagnosis on
otoscopy images [5] has been suggested as a potential tool to improve care of
ear disease. Previous studies have mainly focused on applying machine learning
to classify images of the eardrum. Viscaino et al. compare support vector
machine (SVM), k-nearest neighbor (k-NN) and decision trees on predicting ear
conditions with feature extracted by filter bank, discrete cosine transform
(DCT) and color coherence vector (CCV) [6]. More recently, deep learning has
been applied to classify otoscopy images. Zafer studies the combination of the
fused fine-tuned deep features and SVM model applied to 3 diagnostic class
image classification [7].

Furthermore, Khan et al. [8] compare the performance of a variety of CNN
architectures and find that DenseNet [9] produces the best result in the
classification of 3 diagnostic classes. The authors also use Grad-CAM [10] to
visualize the important regions in the image for the model.

In [11], the authors compare the performance applying 9 kinds of ImageNet [12]
pretrained CNN networks to 6 class otoscopy image classification and then
further improve the performance by ensembling the output of 2 networks. The
prementioned deep learning based methods overall show great performance in the
respective settings, achieving high accuracy ranging from 94% to 99%.

The accuracy reported by previous work, even in multi-class settings, is
remarkable, yet, all previous approaches process high quality still images of
the eardrum. Unfortunately, in a clinical setting the capture of high quality
still images can be challenging in practice, especially given that pediatric
patients are often moving and uncooperative with examination. The use of a
single image also increased the chance of failing to capture a clinically
accurate representation of the anatomy due incomplete view. A method that
analyzes otoscopy video sequences rather than still images could help overcome
aforementioned shortcomings.

The straightforward approach to have the whole video sequence as input is to
train a 3D (2D images plus temporal dimension) CNN mapping from videos to
labels. However, that would require 1) heuristics to harmonize videos to a
preset length, and 2) large amounts of labeled data to avoid over-fitting,
especially considering that the class frequency in any otoscopy video dataset
is likely to be heavily skewed towards normal cases. The present work
addresses these challenges by framing otoscopy interpretation as a video
anomaly detection problem. Under our setting, the model is trained with normal
videos only and will flag videos that diverge from the normal as an anomaly
during testing.

The contribution of this work is two-fold. First, we devise a two stage model
to apply on full video sequences collected in clinics, which differentiates
our approach from all previous methods that only considered high quality still
images. A region detector first extracts eardrum patches from video frames
that are then used to perform anomaly detection in the following stage. This
architecture restricts anomaly detection to semantically useful regions.
Second, we develop the shift contrastive anomaly detection (SCAD) method that
leverages color-jitter based distributional shift-transformation to improve
the separability between normal and abnormal data. Tailored to our
application, our self-supervision task enforces the model to leverage subtle
color features to identify abnormality during testing. Even when the amount of
training videos is relatively small, our approach generates superior binary
(normal / abnormal) otoscopy video screening results compared to both a
baseline method and the average of a group of clinicians. We believe that this
work constitutes a promising first step towards achieving algorithmic decision
support for the diagnosis of pathology of the middle and external ear, and the
development of smart tool which might aid clinicians in the accurate diagnosis
and management of common ear disease.

##  II Related work

The task of anomaly detection [13, 14], also referred to as out-of-
distribution or novelty detection depending on context, is to identify
unfitted data samples. This is a crucial task in many real world applications
such as detecting system malfunction, financial fraud and health issues.
Anomaly detection is a useful but challenging task, because anomalies are hard
to define. In a supervised anomaly detection setting, abnormal samples are
provided to the model in the training stage to give guidance on what qualifies
as an ”anomaly”. However, there could be many reasons a data sample is
abnormal such that collecting a representative amount of abnormal samples is a
hard problem itself. This is especially challenging in medical applications,
where the dataset is very likely to heavily skewed towards normal class.
Unsupervised anomaly detection (UAD), on the other hand, does not need normal
/ abnormal supervision signal during training. Current approaches usually
utilize only the normal data during training, which makes it more appealing to
applications with skewed dataset. Therefore, in the scope of this paper, we
focus only on the UAD setting where the detector can only access to samples
from the normal data distribution during training. Moreover, in many real
world applications, the task is perform UAD on high-dimensional imagery data
(e.g., in our case, detecting anomaly from endoscopic videos), which makes the
problem even more challenging.

A rich set of methods [13] has been proposed to study UAD problem on imagery
data, approaching the problem with four main paradigms: i) Density-based
methods first estimate the normal data density using methods including
Gaussian Mixture and Energy Based Models, and then detect data with low
estimated density as anomaly. DAGMM [15] and DSEBM [16] are methods that
belongs to this category. ii) One-class classifier-based methods fit a
classifier, e.g., Deep SVDD [17] and DROC [18], to separate normal and all
other data and then use it to detect anomalies. iii) Reconstruction-based
techniques learn a reconstruction model, e.g., AnoGAN [19], of normal images
and detect anomalies as samples with high reconstruction error. iv) Self-
supervised-based methods learn a feature extractor using self-supervised
tasks, such as distinguishing whether or not certain transformation has been
applied to the image [20]. Then, during learning features of augmented
versions of the same image should be closer than features of different images
[21]. Upon convergence, anomaly detection is performed on the extracted
features. Notable methods along this line include SVD-RND [22], CutPaste [23],
CSI [24], SSD [25] and MSC [26]. UAD has also been applied to medical imaging
[27] across many domains, including X-ray [28, 29], CT [30, 31], MRI [32, 33,
34] and endoscopy [35] datasets.

Anomalies in medical imaging tend to be more subtle and only reside in a small
region of the whole image, which makes it particularly challenging to apply
UAD in this field. To address this problem, [36] proposes to connect Context-
encoding (CE) and Variational Autoencoders (VAE) to combine both
reconstruction and density based anomaly score. Doing so shows superior
performances in three brain MRI image datasets. To directly target subtle
anomaly in small regions, there are efforts to develop self-supervised tasks
that are tailored to a certain anomaly group. It is desirable to apply
artificial transformation of data that is similar to real anomaly so that the
self-supervised task can force the model to learn features that are
discriminative enough to detect anomaly in testing. As an example, CutPaste
[23] creates discontinuous defects by cutting a region of an image and pasting
it to another location. While this method works well for industrial visual
inspection, the introduced sharp discontinuity is uncommon in medical imagery
making this method less useful for anomaly detection in medical data.
Alternatively, foreign patch interpolation (FPI) [37] creates artificial
defects by interpolating a small local patch with a foreign patch of the
image. A network is trained to estimate the pixel-wise interpolating factor
from the synthesised image. The estimated interpolating factor is then used as
anomaly score during testing. Success is demonstrated using Brain MRI and
abdominal CT image datasets. Poisson image interpolation (PII) [38] extends
this idea by using Poison image editing to blend in the foreign patch instead
of direct linear interpolation, targeting more subtle and continuous
irregularities. Improvement in performance is shown using chest X-ray and
fetal ultrasound datasets. Meanwhile, both FPI and PII utilize self-supervised
tasks targeted at regional shape anomaly which is less applicable to our
problem where subtle color anomaly is more crucial.

##  III Methodology

Figure 1: The screening architecture of our shift contrastive anomaly
detection (SCAD) is composed of three sub-blocks: eardrum detection,
representation learning and anomaly detection. During training, eardrum
detection is learned by supervised training with provided labels.
Representation learning is conducted with self-supervision to obtain an
embedding function $\phi$. During testing, we use the embedding function
$\phi$ to produce frame level anomaly score and video level anomaly score to
perform anomaly detection.

The otoscope video database is skewed towards the non-disease status given the
practice of capturing bilateral ear examinations and occurrence of non-ear
disease presenting with symptoms such as ear pain. Dataset imbalance such as
this can present significant challenge for the development of robust design
boundaries, particularly when relying on small data-sets. To overcome such
limitations, we propose use unsupervised anomaly detection approach to provide
computer-aided screening for otoscopy video. The formal definition of this
problem is as follows:

We are given a large training dataset $\mathcal{D}_{\emph{train}}$ comprising
only normal video sequences and a smaller testing dataset
$\mathcal{D}_{\emph{test}}$ comprising both normal and abnormal video
sequences. Given the training dataset, our objective is to learn an anomaly
score function $A$(${v}$) with video $v$ as input, and detect the abnormal
videos as anomaly during testing. Ideally, the function learns to map normal
videos to small anomaly scores and abnormal videos to large anomaly scores.
During testing, we threshold the score, where $A$($s$) $>\psi$ indicates
anomaly.

In otoscopy video diagnosis, the temporal relation between frames is not as
influential as in action recognition tasks. Whether or not a video is abnormal
is directly determined by the existence of abnormal frames therein. With this
observation, we turn the video anomaly detection problem into a frame anomaly
detection problem with an additional video level aggregation function.
Moreover, only some of the frames in the video are vital to diagnosis, and in
those frames, only the region that visualizes the eardrum. Therefore, it would
be desirable to only conduct anomaly detection on such frames, and the eardrum
region in particular. As shown in Figure 1, our method has two steps during
training: Step 1) Supervised learning of a detection module that extracts
eardrum patches $x$ from given video $v$ using provided labels. Step 2) Self-
supervised representation learning on the labeled region of interest to obtain
an embedding function that ideally projects normal and abnormal samples to
distinct regions of the embedding space. During testing, we compute anomaly
scores using the detected patches and the learned embedding function.

###  III-A Eardrum detection

To detect the frames and regions showing the eardrum, we use a convolutional
neural network and train it in supervised manner using image-label pairs.
Different from the standard detection problem where each image may have
multiple instances, we can only have at most one eardrum in each frame.
Therefore, we do not use a standard object detection architecture, and rather
choose to have a single convolutional neural network to predict the binary
label and corresponding bounding box. We define $\mathcal{L}_{\emph{cls}}$ as
the cross entropy loss of groundtruth label and the classification score, and
define $\mathcal{L}_{\emph{local}}$ as the $L1$ loss of the groundtruth
bounding box and the predicted bounding box. Then, the overall loss is simply
combining the classification loss and localisation loss.

$\mathcal{L}_{\emph{detect}}=\mathcal{L}_{\emph{cls}}+\mathcal{L}_{\emph{local}}$
(1)

Given our problem setup, we train the detection module on the normal training
videos only. Even though the model does not use abnormal videos during
training, we find it generalizes well to abnormal frames during testing as
shown in section IV-B.

###  III-B Shift Contrastive Anomaly Detection (SCAD)

Figure 2: Top: On the mean-shifted representation space,
$\mathcal{L}{\emph{msc}}$ maximizes the angle between negative pairs. Shift-
transformed samples are treated as additional instances in this step. Bottom:
On the angular representation space $\mathcal{L}{\emph{shift\\_angular}}$
increases separability of anomaly by i) increasing the angle between the
embedding of normalized center and the shift-transformed samples; ii)
decreasing the angle between embedding of normalized center and normalized
samples.

Contrastive Learning: Contrastive learning [21, 39] has become the top
performing self-supervised learning method in recent years. In this paradigm,
the first step of the training procedure is to sample a minibatch of size $N$,
and perform augmentation twice to each sample $x_{i}$ to obtain
$(x_{i^{\prime}},x_{i^{\prime\prime}})$ termed as positive pair, producing
$2N$ samples in total. All images are passed through a feature extractor and
then features are typically scaled to the unit sphere by $l_{2}$ normalization
to get representation $\phi$. The contrastive loss to pair
$(x_{i^{\prime}},x_{i^{\prime\prime}})$ is then defined as

$\mathcal{L}_{\emph{con}}(x_{i^{\prime}},x_{i^{\prime\prime}})=-\log{\frac{\exp(\phi({x_{i^{\prime}}})\cdot{\phi(x_{i^{\prime\prime}})}/\tau)}{\sum_{m=1}^{2N}{\mathbbm{1}[i^{\prime}\neq
m]\cdot\exp(\phi({x_{i^{\prime}}})\cdot{\phi(x_{m})}/\tau)}}}\,,$ (2)

where $\tau$ is a temperature hyper-parameter.

The contrastive learning objective pulls $x_{i^{\prime}}$ close to
$x_{i^{\prime\prime}}$ and pushes all other samples away from
$x_{i^{\prime}}$. Intuitively, this will make the learned embedding capture a
latent structure that is meaningful enough to separate samples from one
another and thus be beneficial for downstream tasks. This paradigm shows great
success in self-supervised training for image recognition, but has inherent
problems for anomaly detection. To minimize the loss function, the angles
between positive and negative samples need to be maximized even though both
samples are from the normal class. This results in a scenario where
representations of normal samples span across the whole unit sphere. During
testing, anomaly sample could therefore be projected onto a location that is
close to normal embeddings, which challenges the paradigm of anomaly
detection.

Figure 3: (a) Normal eardrum region examples. (b) Shift-transformation
variants applied to the second normal example of each row. From left to right:
color-jitter-random-cut (CJ-RC), color-jitter-random-region (CJ-RR), color-
jitter-whole-frame (CJ-WF). (c) Abnormal eardrum region examples.

Mean-Shifted Contrastive Loss: To adapt contrastive learning to anomaly
detection, MSC [26] proposes the alternative mean-shifted contrastive loss.
Rather than directly minimizing the contrastive loss in the representation
space, MSC constructs a mean-shifted counterpart, by subtracting the center
$c$ of the whole training set and then normalizing to the unit sphere. For a
given sample $x$, the mean-shifted embedding is defined as

$\theta(x)=\frac{\phi(x)-c}{\left\lVert\phi(x)-c\right\rVert}\,.$ (3)

The mean-shifted loss is then constructed by applying typical contrastive loss
on this mean-shifted embedding space:

$\mathcal{L}_{\emph{msc}}(x^{\prime},x^{\prime\prime})=\mathcal{L}_{\emph{con}}(\theta(x^{\prime}),\theta(x^{\prime\prime}))\,.$
(4)

Because the mean-shifted representation is normalized around $c$, the training
no longer spans the normal samples across the whole unit sphere of $\phi(x)$.
This makes anomaly samples more separable from normal samples in the $\phi(x)$
space.

Distributional shift angular loss: To further improve separability, MSC [26]
uses additional angular center loss

$\mathcal{L}_{\emph{angular}}=-\phi(x)\cdot c\,,$ (5)

to shrink normal samples around the normalized center. The assumption is that
normal data lying in a small region around the center will be more
discriminative. The total loss that MSC uses then becomes
$\mathcal{L}_{\emph{angular}}+\mathcal{L}_{\emph{angular}}$.

This approach works well in datasets composed of regular images (e.g., CIFAR
dataset [40]), where images from different classes are treated as anomaly.
Empirically we find that this strategy does not perform as well in our medical
imaging application. The main reason is that we are having much smaller
semantic variations between normal and abnormal images so that their
respective embeddings are still close to each other even after training using
mean-shifted contrastive and angular center loss.

To address this issue, we develop the distributional shift angular center loss
to further enhance separability between normal and abnormal samples. Our
method constructs additional shift-transformed samples $z$ from $x$ and then
pushes $\phi(z)$ away from the normalized training center $c$. Shift-
transformed samples $z$ are created by applying distributional shift-
transformations on $x$ so that $z$ are no longer from the original data
distribution. The intuition is straightforward: not only should the normal
samples occupy a small region around the center, but also the samples not
drawn from the original distribution cannot be too close to the center either.
Hinge loss is applied to the angle between center $c$ and shift-augmented
samples $z$. The loss function takes the form:

$\mathcal{L}_{\emph{shift\\_angular}}=-\phi(x)\cdot c+max(0,1-\phi(z)\cdot
c)\,.$ (6)

Ideally, these shift-transformed samples $z$ are semantically similar to real
anomalies. In the context of otoscopy images, color is an important feature in
that infection may lead to change of color of the whole or part of the
eardrum. Thus, we employ three variants of color-jitter as the distributional
shift-transformation. They are color-jitter-random-cut (CJ-RC), color-jitter-
random-region (CJ-RR), and color-jitter-whole-frame (CJ-WF). CJ-RC uses a
random rectangle as a mask of transformation. CJ-RR is constructed by
interpolating the color-jittered image and the original image using pixel-wise
weights for each image. We first initialize the pixel-wise weights to be all
zero, and then pick random points to be filled with one, and apply Gaussian
filter to obtain smoothed weights. CJ-WF simply applies color-jitter to the
whole frame. Figure 3 presents illustrations of the three augmentation
strategies.

Final loss: We combine i) the mean-shifted contrastive loss
$\mathcal{L}_{msc}$ and ii) the angular loss with distributional shift
augmentations $\mathcal{L}_{shift\\_angular}$:

$\mathcal{L}_{final}=\mathcal{L}_{msc}+\mathcal{L}_{shift\\_angular}\,.$ (7)

This loss function is tailored to this specific medical imaging application
and enjoys the best of both objectives. The mean-shifted loss makes the
features to be representative of the images and the angular loss encourages
normal and abnormal instances to be distant to each other in the feature
space. Therefore, the combination of these two losses achieves superior
performances, which is demonstrated through ablation in section IV-C.

Frame level anomaly score: In order to classify a sample as normal or
abnormal, we use a simple criterion based on k-Nearest-Neighbors (kNN). kNN
predicts the label of a data point by ensembling the k closest labeled samples
according to some distance measure. In this work, we use the cosine distance
between the features of the target image $x$ and those of all training images.
The anomaly score is then given by:

$a(x)=\Sigma_{\phi(y)\in N_{k}(x)}(1-\phi(x)\cdot\phi(y))\,,$ (8)

where $N_{k}(x)$ denotes the $k$ nearest features to $\phi(x)$ in the training
feature set $\\{\phi(z)\\}_{z\in\mathcal{X}_{train}}$.

Video level anomaly score: We consider a relatively simple method to aggregate
the frame level anomaly score to construct the video level anomaly score. The
aggregation function is simply taking the average of the frame level anomaly
score across all detected eardrum patches. The video level anomaly score is
defined as

$A(s)=\overline{\\{a(x)\\}}\,,$ (9)

where $\\{a(x)\\}$ represents the set of frame level anomaly scores for all
frames in the given video $s$.

In summary, during training, we learn a frame level embedding function $\phi$
with all individual normal frames in the training set. During testing, we
compute the kNN-based anomaly score using the embedding feature from $\phi$
for all detected frames. We then aggregate the predicted anomaly score for
each detected frame in testing video to produce a video level anomaly score
$A$ through averaging. Such a score is expected to be small for normal videos
and large for abnormal videos. Any frame that is substantially different from
normal frames in the training set (e.g., a red, bulging eardrum) would lead to
a large distance to training examples in embedding space and then produce
large image and video level anomaly score. The threshold can be selected to
meet desired clinical requirement (e.g., a certain specificity value) and in
practice is selected using the validation dataset consisting of both normal
and abnormal videos.

(a) (b)

Figure 4: Our detection module successfully extracts eardrum regions in video
sequences. Green bounding box represents groundtruth eardrum region; orange
bounding box represents detected eardrum region. Top Row: Normal video
example. Bottom Row: Abnormal video example.

##  IV Experimental results

###  IV-A Data preparation

We collected a total of 100 otoscopy videos from pediatric patients that were
seen for various conditions in an urban, tertiary-care pediatric emergency
department that treats over 34,000 patients, aged between 0 and 21 years of
age per year. Patients were recruited as a convenience sample based on
willingness to participate in the study. At this site patient ages range from
0-22 years. Study protocols were approved by the local institutional review
board. The length for each video ranges from 5 to 40 seconds. Twenty seven to
thirty frames per second are extracted from the video. Two otolaryngologists
(ears, nose, throat specialists) annotated every video and consensus agreement
was used as ground truth. While the expert annotations are more granular, due
to the small size of the dataset, we decided to restrict analysis to two
classes and each video is assigned into either the normal or abnormal
category. Additional to binary categorical label on video level, we also
annotate frame level binary labels and bounding boxes for each appearance of
the eardrum, where the eardrum occupies at least 20% of the image size in
width and height. Of the 100 videos, 80 are labeled normal and 20 are labeled
abnormal. We assign 60 normal videos to the training set, 10 to the validation
set and 10 to the testing set. For anomaly videos, we assign 10 to the
validation set and 10 to the testing set. The testing set is used for
evaluation of the algorithms and comparison with human clinicians. We do not
have overly apparent abnormal examples (e.g., the whole eardrum is filled with
blood) during testing which makes it a challenging comparative study with
human clinicians.

Video capture of ear exam was performed using Cellscope Oto attached to a
first-generation iPhone SE. Maintaining relative uniformity of video settings
such as color are critical to the accuracy of the AI algorithm. Clinical
deployment of such technology therefore presents significant challenge given
the ad hoc selection of otoscopes used in the clinical setting and perhaps
more importantly the overwhelming reliance on non-digital otoscopes in
clinical practice. As a result we opted to use Cellscope Oto as an affordable
solution that can potentially be deployed at scale. The Cellscope Oto device
is designed to attache to an ear speculum and provides an additional light
source but relies on the iPhone camera to capture the exam. We do note that
the image and video quality captured on this device is not as high as reported
in previous works [6, 41], especially with respect to as sharpness and
clarity. Unlike Cellscope oto however more advanced high quality otoscopes
with digital video capacity (as used in [6, 41]) are often cost prohibitive to
wider adoption.

###  IV-B Eardrum detection

We first train a model for eardrum detection to extract eardrum patches. We
use Resnet-101 as backbone network and choose the stochastic gradient descent
(SGD) optimizer with $1e-3$ learning rate, $0.9$ momentum and $1e-4$ weight
decay. The model is trained using all frames from the normal training set
videos for 5000 epochs with a batch size of 128. The input images are resized
to 256 by 256 and randomly cropped to 224 by 224, followed by a random
selection of data augmentation of color-jitter, cutout, rotation and shearing.

Results: Since under our setup there is at most one instance of a single
class, we evaluate the performance of our model by accuracy of eardrum
detection under different intersection over union (IoU) [42] thresholds. A
correct detection is defined as a frame that produces a correct frame-level
classification and, if the given frame shows the eardrum, also produces a
bounding box with IoU score larger than a certain threshold in that given
frame. As shown in Table I, our model produces high accuracy across varied IoU
thresholds. A visualization of eardrum detection in video sequences can be
found in Figure 4. Note that even though the model is supervised only on
instances from training normal videos, it generalizes well to both normal and
abnormal videos in the testing set. This indicates that the network extracts
features that are more consistent to the image semantic and agnostics to the
details on eardrum.

| Accu (IoU$>$.50) | Accu (IoU$>$.75) | Accu (IoU$>$.90)
---|---|---|---
Normal | $77.3$ | $69.7$ | $37.8$
Abnormal | $82.1$ | $79.3$ | $34.4$
Overall | $80.8$ | $73.1$ | $35.7$
TABLE I: Detection performance evaluated by accuracy % under different IoU
thresholds.

###  IV-C Anomaly detection

Following the architecture suggested in [26], we construct the feature
extractor $\phi$ to be the ResNet-101 convolutional neural network followed by
an additional l2 normalization layer. We augment the normal samples by
sequentially applying 224 by 224-pixel crop from a randomly resize image and
random horizontal flips. The model is initialized with the ImageNet pre-
trained model and then the two last blocks of the ResNet-101 are fine-tuned
for 5000 epochs with the loss function in Equation 7. The temperature $\tau$
in 2 is set to be 0.25. We use SGD optimizer with learning right of $1e-5$
weight decay of $5e-5$, and no momentum. For each minibatch size of 60, we
sample one frame from each video to prevent similar frames treated as negative
pairs. During training, we use the groundtruth bounding box to extract patches
from videos. We use KNN (k=2) in computing the frame-level anomaly score. We
conduct experiments with three random seeds for evaluations.

Results: We evaluate our methods on the anomaly detection task in the
previously introduced otoscopy video dataset. We adopt the Area Under the
Receiver Operating Characteristic curve (AUROC) and the Area Under the
Precision Recall Curve (AUPRC) as anomaly detection performance score. We
compare our approaches against the recently proposed top performing method MSC
[26]. The training and testing procedures are the same as ours except that it
does not utilize shift-transformations and uses the angular loss (3) instead
of shift angular loss (4) in computing the final loss function. In Table II,
we present a comparison between MSC and three variants of SCAD evaluated by
AUROC and AUPRC. We see that the SCAD-CJ-WF variation performs the best,
reaching an AUROC of 88.0% and AUPRC of 87.9% on average. This suggests that
i) color is indeed a very important feature in separating anomaly from normal
in otoscopy images, and that ii) color jitter on the whole image can produces
semantically more similar samples to the real anomaly. This could be due to
that both random cut and random region introduce artificial patterns on the
image that are utilized in the self-supervised training but are not seen in
the abnormal samples during testing. To intuitively visualize the feature
embedding produced by different methods, we plot the latent embedding for all
testing frames in Figure 5 and the corresponding kNN based frame-level anomaly
score in Figure 6. As we can see, SCAD-CJ-WF generates the embedding that is
most separable both qualitatively and quantitatively as measured by the
anomaly score.

Method | AUROC | AUPRC
---|---|---
MSC | $64.3\pm 4.0$ | $60.7\pm 2.4$
SCAD-CJ-RC | $78.0\pm 1.0$ | $69.1\pm 1.1$
SCAD-CJ-RR | $82.3\pm 1.2$ | $74.4\pm 1.1$
SCAD-CJ-WF | $\mathbf{88.0\pm 1.0}$ | $\mathbf{87.9\pm 1.8}$
TABLE II: Results evaluated by AUROC % and AUPRC %.

(a) MSC. (b) AMSC-CJ-RC. (c) AMSC-CJ-RR. (d) AMSC-CJ-WF.

Figure 5: Visualization of the embedding for all frames in the testing set of
normal and abnormal videos using t-SNE.

(a) MSC. (b) AMSC-CJ-RC. (c) AMSC-CJ-RR. (d) AMSC-CJ-WF.

Figure 6: Distribution of the frame-level anomaly score for all frames in the
testing set of normal and abnormal videos.

Training objective: The individual effect of each loss component is presented
in Table III. We note that neither $\mathcal{L}_{\emph{msc}}$ or
$\mathcal{L}_{\emph{angular}}$ individually performs well in our dataset.
$\mathcal{L}_{\emph{shift\\_angular}}$ outperforms all other individual
objective and combining it with $\mathcal{L}_{\emph{msc}}$ results in further
improvement.

Metric | $\mathcal{L}_{\emph{msc}}$ | $\mathcal{L}_{\emph{angular}}$ | $\mathcal{L}_{\emph{shift\\_angular}}$ |  $\mathcal{L}_{\emph{msc}}$ \+ $\mathcal{L}_{\emph{shift\\_angular}}$
---|---|---|---|---
AUROC | $59.3\pm 0.6$ | $57.7\pm 4.0$ | $78.0\pm 1.0$ | $\mathbf{88.0\pm 1.0}$
AUPRC | $55.1\pm 1.6$ | $54.5\pm 3.0$ | $69.8\pm 0.8$ | $\mathbf{87.9\pm 1.8}$
TABLE III: Training objective ablation study (AUROC % and AUPRC %).
$\mathcal{L}_{\emph{shift\\_angular}}$ is of the SCAD-CJ-WF variant. Note that
the MSC method uses $\mathcal{L}_{\emph{msc}}$ \+
$\mathcal{L}_{\emph{angular}}$ as loss function and the result is presented in
Table II

###  IV-D Comparison with human clinicians

To further evaluate our performance in a real-world setting, we performed a
cross-sectional study comparing the performance of clinicians versus our
model. The study was approved by the local institutional review board.

Study Setting: Clinicians who routinely evaluate eardrums by otoscopy in a
primary care, urgent care, or emergency medicine setting were invited to
participate in this study. Qualifying clinicians included physicians (MD/DO),
nurse practitioners (NP), and physician assistants (PAs). These clinician
group reflects frontline healthcare providers who typically diagnose and
manage ear disease. Exclusion criteria were practicing outside of the United
states, employment by the Johns Hopkins University or Johns Hopkins Hospital,
visual impairment limiting video analysis, or non-English speakers.
Recruitment took place through digital messaging via professional email
mailing lists and professional social media groups during November 2020.

Study Procedure: Ten normal and ten videos designated as acute otitis media
(AOM) by ground truth were selected as the testing set of images, as described
above. Clinicians completed an online servey in which they were instructed to
review the 20 eardrum videos without any further clinical information
provided. Following each video, participants were asked to indicate whether
the eardrums appeared normal or abnormal. They were also asked to rate their
level of confidence in the accuracy of their diagnosis using a Likert scale
(1= not at all confident to 4 = very confident). Participants were also asked
to report the number of times they viewed each video. Participant responses
were collected and managed using REDCap electronic data capture tools [43]
hosted at Johns Hopkins University.

Figure 7: Consort flow diagram of clinician enrollment. After removing ineligible samples, 25 surveys were included for comparative analysis. Characteristic | Sample Size | Accuracy
---|---|---
Sex |  |
Female | 15 (60%) | $73.0\pm 15.1$
Male | 10 (40%) | $60.0\pm 17.5$
Current Practice |  |
General Pediatrics | 12 (48%) | $73.0\pm 15.1$
Pediatric Specialty | 5 (20%) | $70.0\pm 22.4$
General Internal Medicine | 3 (12%) | $55.0\pm 8.7$
Others | 5 (20%) | $67.0\pm 8.4$
TABLE IV:  Clinician demographics and corresponding accuracy (mean % and
standard deviation %) by groups. Other practices include emergency medicine,
family medicine and urgent care.

Study Population: A total of 59 clinicians were recruited and assessed for
eligibility. After removing ineligible, declined, or incomplete surveys, 25
surveys were included for analysis. The process is shown in the consort flow
diagram in Figure 7. Participant demographics are summarized in Table IV. All
respondents are early-to-mid career physicians with the majority being female
and practicing pediatrics.

Results: Clinician confidence has a mean of 3 and a standard deviation of 0.5
and the view count per video has a mean of 1.5 and a standard deviation of
0.5. We use linear regression to analyze the correlation of both clinicians’
confidence and view count with their score. No statistically significant
correlation between clinician confidence and the accuracy score (P=.07) or
with view count and score (P=.32) were found (a P-value of $\leq$ .05 was
considered statistically significant).

In Table V, we compare anomaly detection methods against clinicians in
accuracy, sensitivity, specificity and precision. We compute the metrics for
all individual clinicians and then compute the mean and standard deviation of
the group. We choose the decision threshold of machine learning approaches of
which the sensitivity (TPR) on the validation set to be 90%, and then report
the metric values on the testing set. Compared to clinicians, our methods
perform better and achieve lower variances across all metrics. We further
notice that SCAD-CJ-WF outperforms both clinicians and all other methods in
most metrics.

Method | Accuracy | Sensitivity | Specificity | Precision
---|---|---|---|---
Clinicians | $67.8\pm 17.0$ | $63.2\pm 31.8$ | $72.4\pm 17.4$ | $64.5\pm 24.3$
MSC | $51.7\pm 7.6$ | $43.3\pm 15.3$ | $60.0\pm 0.0$ | $51.0\pm 8.6$
SCAD-CJ-RC | $73.3\pm 2.9$ | $70.0\pm 0.0$ | $76.7\pm 5.8$ | $75.2\pm 4.5$
SCAD-CJ-RR | $78.3\pm 7.6$ | $\mathbf{90.0\pm 0.0}$ | $66.7\pm 15.3$ | $73.7\pm 8.8$
SCAD-CJ-WF | $\mathbf{81.7\pm 2.9}$ | $83.3\pm 5.8$ | $\mathbf{80.0\pm 0.0}$ | $\mathbf{80.6\pm 1.0}$
TABLE V: Comparison of machine learning models against a group of clinicians
(n=25) in performances (mean and standard deviation of accuracy %, sensitivity
%, specificity % and precision %) on testing set. Machine learning models’
decision thresholds are chosen by setting the sensitivity (TPR) on the
validation set to be 90%. Figure 8: Overlaying clinicians’ TPR/FPR samples
(n=25) and group average on SCAD-CF-WF’s Receiver Operative Characteristic
(ROC) curve. Multiple clinicians’ responses may have the same TPR/FPR values.

In Figure 8, we provide further visualizations comparing our best performing
model SCAD-CJ-WF to clinicians. All seeds of our model produced an ROC curve
that exceeds the average clinician response. Controlling the same False
Positive Rate, our model outperforms the clinician performance in True
Positive Rate (Sensitivity) by a large margin from 61.5% to 90%. Besides the
improvement in average performance, our model also generates more consistent
predictions. We can see from the figure that different clinicians produce
results that are not consistent with each other. While the best clinician
included in our study did achieve a perfect True and False Positive Rate, the
majority of clinician evaluations are scattered with a large variation. The
wide variability and lower accuracy of clinician’s scores from this sample are
comparable to prior reports of pediatricians and general practitioners at
identifying AOM from otoscopy [44, 45]. This variability between clinician
scores reflects the clinical challenge of accurately diagnosing AOM.
Individual factors such as a clinicians’ training and experience, combined
with patient factors such as cooperativeness and a non-obstructed view
contribute to a successful diagnosis. In contrast, our model across different
seeds generated more consistent and more accurate results. Given the clinician
variability to diagnose AOM illustrated our comparative study and previous
studies, there is significant value to improve diagnosis in a clinical setting
by adopting deep learning based anomaly detection screening like that of the
present study.

Discussion: Our presented method shows promising performance in anomaly
detection on otoscopy video sequences. The feasibility of adopting computer-
aided diagnostics to detect and report anomaly may assist front-line
clinicians in primary care and eventually even patients at-home.

The eventual successful deployment of such technology necessitates more than
the successful development of an diagnostic AI algorithm with high accuracy.
For such technology to integrate within current healthcare models and work
flow there is a need to create a whole new technological infrastructure. The
implementation of such technology beyond it obvious need to be adopted by
healthcare providers, perhaps the more pressing questions where the required
financial investment for the creation of this infrastructure will come from.

Meanwhile, our study exhibits several limitations at the current stage. i) The
dataset is relatively small, skewed towards normal and collected through
convenience sampling. As such our data serve to provide some additional
support for the proof of concept but does not negate the need for larger and
more rigorously collected data-set before such technology becomes clinically
viable. The development of a large video database would likely further improve
the algorithm’s diagnostic performance and enable training of specific
diagnosis as apposed to being limited to detecting normal vs abnormal ear
exams. ii) Cellscope oto does not permit pneumotoscopy (blow air at eardrum to
evaluate for movement) which can also serve as a critical step for accurate
diagnosis of ear disease. iii) In our comparative study clinicians were not
provided with patient history which might increase their diagnostic accuracy
in a clinical setting.

##  V Conclusion

In this study, we present a two-stage anomaly detection method that is
designed to perform normal / abnormal classification for otoscopy videos and
that can be developed based on small dataset and highly skewed class
distribution. We demonstrate that our method outperforms baseline algorithms
as well as the average of 25 clinicians, constituting a promising step towards
computer-aided diagnosis of the middle and external ear pathology. The further
development of computer-aided diagnosis method, like ours, could contribute to
timely and high-quality management of ear disease.

## Acknowledgment

The authors gratefully acknowledge funding by the Leon Lowenstein Foundation.

## References

  * [1] H. Gurnaney, D. Spor, D. G. Johnson, and R. Propp, “Diagnostic accuracy and the observation option in acute otitis media: the capital region otitis project,” _International journal of pediatric otorhinolaryngology_ , vol. 68, no. 10, pp. 1315–1325, 2004.
  * [2] D. L. Brinker Jr, E. L. MacGeorge, and N. Hackman, “Diagnostic accuracy, prescription behavior, and watchful waiting efficacy for pediatric acute otitis media,” _Clinical pediatrics_ , vol. 58, no. 1, pp. 60–65, 2019.
  * [3] N. M. Poole, D. J. Shapiro, K. E. Fleming-Dutra, L. A. Hicks, A. L. Hersh, and M. P. Kronman, “Antibiotic prescribing for children in united states emergency departments: 2009–2014,” _Pediatrics_ , vol. 143, no. 2, 2019\.
  * [4] R. M. Rosenfeld, L. Culpepper, K. J. Doyle, K. M. Grundfast, A. Hoberman, M. A. Kenna, A. S. Lieberthal, M. Mahoney, R. A. Wahl, C. R. Woods Jr, _et al._ , “Clinical practice guideline: otitis media with effusion,” _Otolaryngology-Head and Neck Surgery_ , vol. 130, no. 5, pp. S95–S118, 2004\.
  * [5] T. Canares, W. Wang, M. Unberath, and J. Clark, “Artificial intelligence to diagnose ear disease using otoscopic image analysis: A review,” _Journal of Investigative Medicine_ , In production.
  * [6] M. Viscaino, J. Maass, P. Delano, M. Torrente, C. Stott, and F. A. Cheein., “Computer-aided diagnosis of external and middle ear conditions: A machine learning approach,” _PLoS_ , vol. 15, p. 3, 2020.
  * [7] C. Zafer, “Fusing fine-tuned deep features for recognizing different tympanic membranes,” _Biocybernetics and Biomedical Engineering_ , vol. 40(1), pp. 40–51, 2020.
  * [8] M. Khan, S. Kwon, J. Choo, S. Hong, S. Kang, I. Park, and et al, “Automatic detection of tympanic membrane and middle ear infection from oto-endoscopic images via convolutional neural networks,” _Neural Networks_ , 2020.
  * [9] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger, “Densely connected convolutional networks,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 4700–4708.
  * [10] R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, and D. Batra, “Grad-cam: Visual explanations from deep networks via gradient-based localization,” in _Proceedings of the IEEE international conference on computer vision_ , 2017, pp. 618–626.
  * [11] D. Cha, C. Pae, S. B. Seong, J. Choi, and H. Park, “Automated diagnosis of ear disease using ensemble deep learning with a big otoendoscopy image database,” _EBioMedicine_ , vol. 45, p. 606–614, 2019.
  * [12] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, “Imagenet: A large-scale hierarchical image database,” in _2009 IEEE conference on computer vision and pattern recognition_. Ieee, 2009, pp. 248–255.
  * [13] R. Chalapathy and S. Chawla, “Deep learning for anomaly detection: A survey,” _arXiv preprint arXiv:1901.03407_ , 2019.
  * [14] L. Ruff, J. R. Kauffmann, R. A. Vandermeulen, G. Montavon, W. Samek, M. Kloft, T. G. Dietterich, and K.-R. Müller, “A unifying review of deep and shallow anomaly detection,” _Proceedings of the IEEE_ , 2021.
  * [15] B. Zong, Q. Song, M. R. Min, W. Cheng, C. Lumezanu, D. Cho, and H. Chen, “Deep autoencoding gaussian mixture model for unsupervised anomaly detection,” in _International conference on learning representations_ , 2018.
  * [16] S. Zhai, Y. Cheng, W. Lu, and Z. Zhang, “Deep structured energy based models for anomaly detection,” in _International Conference on Machine Learning_. PMLR, 2016, pp. 1100–1109.
  * [17] L. Ruff, R. Vandermeulen, N. Goernitz, L. Deecke, S. A. Siddiqui, A. Binder, E. Müller, and M. Kloft, “Deep one-class classification,” in _Proceedings of the 35th International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, J. Dy and A. Krause, Eds., vol. 80. PMLR, 10–15 Jul 2018, pp. 4393–4402. [Online]. Available: http://proceedings.mlr.press/v80/ruff18a.html
  * [18] K. Sohn, C.-L. Li, J. Yoon, M. Jin, and T. Pfister, “Learning and evaluating representations for deep one-class classification,” _arXiv preprint arXiv:2011.02578_ , 2020.
  * [19] T. Schlegl, P. Seeböck, S. M. Waldstein, U. Schmidt-Erfurth, and G. Langs, “Unsupervised anomaly detection with generative adversarial networks to guide marker discovery,” in _International conference on information processing in medical imaging_. Springer, 2017, pp. 146–157.
  * [20] A. Kolesnikov, X. Zhai, and L. Beyer, “Revisiting self-supervised visual representation learning,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2019, pp. 1920–1929.
  * [21] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A simple framework for contrastive learning of visual representations,” in _International conference on machine learning_. PMLR, 2020, pp. 1597–1607.
  * [22] S. Choi and S.-Y. Chung, “Novelty detection via blurring,” _arXiv preprint arXiv:1911.11943_ , 2019.
  * [23] C.-L. Li, K. Sohn, J. Yoon, and T. Pfister, “Cutpaste: Self-supervised learning for anomaly detection and localization,” 2021.
  * [24] J. Tack, S. Mo, J. Jeong, and J. Shin, “Csi: Novelty detection via contrastive learning on distributionally shifted instances,” _Advances in neural information processing systems_ , 2020.
  * [25] V. Sehwag, M. Chiang, and P. Mittal, “Ssd: A unified framework for self-supervised outlier detection,” _arXiv preprint arXiv:2103.12051_ , 2021\.
  * [26] T. Reiss and Y. Hoshen, “Mean-shifted contrastive loss for anomaly detection,” _arXiv preprint arXiv:2106.03844_ , 2021.
  * [27] T. Fernando, H. Gammulle, S. Denman, S. Sridharan, and C. Fookes, “Deep learning for medical anomaly detection–a survey,” _arXiv preprint arXiv:2012.02364_ , 2020.
  * [28] D. Davletshina, V. Melnychuk, V. Tran, H. Singla, M. Berrendorf, E. Faerman, M. Fromm, and M. Schubert, “Unsupervised anomaly detection for x-ray images,” _arXiv preprint arXiv:2001.10883_ , 2020.
  * [29] B. Bozorgtabar, D. Mahapatra, G. Vray, and J.-P. Thiran, “Anomaly detection on x-rays using self-supervised aggregation learning,” _arXiv preprint arXiv:2010.09856_ , 2020.
  * [30] D. Sato, S. Hanaoka, Y. Nomura, T. Takenaga, S. Miki, T. Yoshikawa, N. Hayashi, and O. Abe, “A primitive study on unsupervised anomaly detection with an autoencoder in emergency head ct volumes,” in _Medical Imaging 2018: Computer-Aided Diagnosis_ , vol. 10575. International Society for Optics and Photonics, 2018, p. 105751P.
  * [31] N. Pawlowski, M. C. Lee, M. Rajchl, S. McDonagh, E. Ferrante, K. Kamnitsas, S. Cooke, S. Stevenson, A. Khetani, T. Newman, _et al._ , “Unsupervised lesion detection in brain ct using bayesian convolutional autoencoders,” 2018\.
  * [32] C. Baur, S. Denner, B. Wiestler, N. Navab, and S. Albarqouni, “Autoencoders for unsupervised anomaly segmentation in brain mr images: a comparative study,” _Medical Image Analysis_ , p. 101952, 2021.
  * [33] C. Han, L. Rundo, K. Murao, T. Noguchi, Y. Shimahara, Z. Á. Milacski, S. Koshino, E. Sala, H. Nakayama, and S. Satoh, “Madgan: unsupervised medical anomaly detection gan using multiple adjacent brain mri slice reconstruction,” _BMC bioinformatics_ , vol. 22, no. 2, pp. 1–20, 2021.
  * [34] C. Baur, R. Graf, B. Wiestler, S. Albarqouni, and N. Navab, “Steganomaly: Inhibiting cyclegan steganography for unsupervised anomaly detection in brain mri,” in _International Conference on Medical Image Computing and Computer-Assisted Intervention_. Springer, 2020, pp. 718–727.
  * [35] Y. Liu, Y. Tian, G. Maicas, L. Z. Pu, R. Singh, J. W. Verjans, and G. Carneiro, “Unsupervised dual adversarial learning for anomaly detection in colonoscopy video frames,” _arXiv preprint arXiv:1910.10345_ , 2019.
  * [36] D. Zimmerer, S. A. Kohl, J. Petersen, F. Isensee, and K. H. Maier-Hein, “Context-encoding variational autoencoder for unsupervised anomaly detection,” _arXiv preprint arXiv:1812.05941_ , 2018.
  * [37] J. Tan, B. Hou, J. Batten, H. Qiu, and B. Kainz, “Detecting outliers with foreign patch interpolation,” _arXiv preprint arXiv:2011.04197_ , 2020.
  * [38] J. Tan, B. Hou, T. Day, J. Simpson, D. Rueckert, and B. Kainz, “Detecting outliers with poisson image interpolation,” _arXiv preprint arXiv:2107.02622_ , 2021.
  * [39] K. He, H. Fan, Y. Wu, S. Xie, and R. Girshick, “Momentum contrast for unsupervised visual representation learning,” _Proceedings of the IEEE conference on Computer Vision and Pattern Recognition_ , 2020.
  * [40] A. Krizhevsky, G. Hinton, _et al._ , “Learning multiple layers of features from tiny images,” 2009.
  * [41] D. Livingstone, A. S. Talai, and J. Chau, “Building an otoscopic screening prototype tool using deep learning,” _J of Otolaryngol - Head & Neck Surg_, vol. 48, p. 66, 2019.
  * [42] H. Rezatofighi, N. Tsoi, J. Gwak, A. Sadeghian, I. Reid, and S. Savarese, “Generalized intersection over union: A metric and a loss for bounding box regression,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 658–666.
  * [43] P. Harris, R. Taylor, R. Thielke, J. Payne, N. Gonzalez, and J. Conde, “Research electronic data capture (redcap)–a metadata-driven methodology and workflow process for providing translational research informatics support,” _Journal of Biomedical Informatics_ , vol. 42, pp. 2377–381, 2009\.
  * [44] M. E. Pichichero and M. D. Poole, “Comparison of performance by otolaryngologists, pediatricians, and general practioners on an otoendoscopic diagnostic video examination,” _International journal of pediatric otorhinolaryngology_ , vol. 69, no. 3, pp. 361–366, 2005.
  * [45] ——, “Assessing diagnostic accuracy and tympanocentesis skills in the management of otitis media,” _Archives of pediatrics & adolescent medicine_, vol. 155, no. 10, pp. 1137–1142, 2001.


# Topological invariants of sorting networks.

Maxim Arnold and Christian Kondor<EMAIL_ADDRESS><EMAIL_ADDRESS>The University of Texas at Dallas. 800 W.
Campbell Rd., Richardson, TX 75080

###### Abstract.

In this note we investigate finite-type Vassiliev invariants of the pure
braids arising from signed sorting networks.

##  1\. Introduction

For a given $n$, a sorting network on $n$ elements is a representation of the
reverse permutation

$\begin{pmatrix}1&\dots&n\\\ n&\dots&1\end{pmatrix}$

as a product of elementary transpositions ${\bf j}:=(j,j+1)$,
$j=1,\dots,(n-1)$. This product must be of minimal length, being
$\binom{n}{2}$ transpositions in length. We will denote this “length of a
sorting network” by $N$.

Sorting networks have been extensively studied in the past few decades from
the viewpoints of both computer science and pure mathematics. These studies
have been motivated by sorting networks’ rich combinatorial structure and wide
range of applications (see e.g. [AHRV07], [ADHV19] and references therein).
One such application is to view a sorting network as a discretization of the
generalized $1$-dimensional Euler equation (see for example [Bre08]). In this
view, physical space is discretized to a finite number $n$ of particles, the
incompressibility condition is interpreted as the conservation of this $n$,
and the action to be minimized is the squared discrete velocity [RVV16]. Our
aim in this paper is to investigate asymptotic topological invariants of the
trajectories in this limited case. We hope that our study will add to the
discussion regarding the question posed by Arnold [Arn86].

##  2\. Basic Definitions

This paper concerns the symmetric group $\mathfrak{S}_{n}$ on $n$ elements,
since $\mathfrak{S}_{n}$ is generated by the set $\\{{\bf j}:=(j,j+1),\quad
j=1,\dots,n-1\\}$ of elementary transpositions. We will often refer to
$\mathfrak{S}_{n}$ in terms of its Cayley graph, the $n$-permutahedron. To
construct the $n$-permutahedron, consider
$\mathbf{x}=(x_{1},\cdots,x_{j},\cdots,x_{n})$ and a central arrangement of
$(n-1)$-dimensional hyperplanes $p_{j,k}:x_{j}=x_{k}$ in $\mathbb{R}^{n}$.
These hyperplanes split the whole of $\mathbb{R}^{n}$ into $n!$ disjoint
domains. If from one of these domains we choose the base point
$\mathbf{b}=(1,2,\cdots,n)$, we can associate each domain to the image of
$\mathbf{b}$ under reflections through the corresponding hyperplanes. Letting
these images be vertices and representing hyperplane reflections as edges
connecting the images, we construct the graph known as the $n$-permutahedron.
Labelling the vertices of the $n$-permutahedron with the inverse permutations
coming from the components of the vertex provide a labelling on the edges.
This edge labelling corresponds to the generators $\bf j$. For example, a
vertex with coordinates $(3,1,2)$ would be labelled with the permutation

$\begin{pmatrix}1&2&3\\\ 2&3&1\end{pmatrix}=(2,3)(1,2).$

Sorting networks correspond to paths on the $n$-permutahedron. In particular,
a sorting network on $n$ elements corresponds to a minimal-length path on the
$n$-permutahedron that connects the identity permutation to the so-called
“long element” of the symmetric group, as in figure 1.

Figure 1. Permutahedron for $\mathfrak{S}(3)$. Planes $p_{1,2}$, $p_{2,3}$ and
$p_{1,3}$ form central hyperplane arrangement and split $\mathbb{R}^{3}$ into
six chambers. Every chamber corresponds to one of the orderings of components
$(x_{1},x_{2},x_{3})$. The sorting network $\bf 121$ corresponds to the path
from the chamber $(+++):=\\{x_{2}>x_{1},\,x_{3}>x_{1},\,x_{3}>x_{2}\\}$ to the
chamber $(---):=\\{x_{2}<x_{1},\,x_{3}<x_{1},\,x_{3}<x_{2}\\}$.

This, however, is not the only geometric representation of sorting networks.
To see another representation consider the sorting network $S={\bf
j}_{1}\cdots{\bf j}_{N}$ acting on $n$ elements, where $N=\binom{n}{2}$ as
before. We can correspond this network to the trajectories of $n$ particles in
a $1$-dimensional dynamical system during finite time $N$. The trajectory of
the $i$-th particle, then, is the the sequence $s_{k}^{-1}(i)$ of locations of
the letter $i$ in the configuration $s_{k}={\bf j}_{1}\cdots{\bf j}_{k}$. The
set of all trajectories taken together forms the wiring diagram of the sorting
network111Throughout the paper we will denote the sorting network as a
sequence of elementary transpositions $\mathbf{j}_{k}$ written in the order of
wiring diagram, which is opposite to the order of the factors if we write the
corresponding product. For example, the product of transpositions
$(2,3)(1,2)$, providing the permutation $(2,3,1)$ would be denoted as $\bf
12$..

This wiring diagram lends itself to further geometric constructions. If we
assign a signature $+$ or $-$ to each transposition (or “crossing”) in a
wiring diagram, we can construct a braid from a sorting network. This braid
construction can be taken a step further: if we consider two sorting networks
$S={\bf j}_{1}\cdots{\bf j}_{N}$ and $T={\bf k}_{1}\cdots{\bf k}_{N}$, the
word $ST={\bf j}_{1}\cdots{\bf j}_{N}{\bf k}_{1}\cdots{\bf k}_{N}$ corresponds
to a closed loop on the permutahedron tracing a path from the identity to the
long element along the word $S$ and then returning back to the identity along
the word $T$. We will call such a loop a _sorting loop_. If as before we
assign signatures to each crossing in the wiring diagram for $TS$, we create a
pure braid on $n$ elements as in figure 3.

###  2.1. The Free Group on Strands, Finite-Type Invariants.

Constructing pure braids from sorting networks lends us a number of new tools
for analyzing sorting networks and distinguishing them from one another. One
subset of these tools is the finite-type invariants, which may be used to
classify braids, knots, and links (see [PV94] and [MK99]). An invariant of
particular interest is the Milnor number, sometimes referred to as the
generalized linking number of a braid.

We can take the closure of a pure braid – that is, connecting the “starting”
and “terminal” points of each strand – to form a closed link of $n$ loops
embedded in $\mathbb{R}^{3}$. If all of the loops in such a link can be
separated without any cuts, the link and its corresponding pure braid are
called trivial. Two loops, however, may be interlinked. This interlinking
results from two strands of a braid being wrapped around one another. How
interlinked two strands $i$ and $j$ are can be quantified by their linking
number $\mathsf{lk}(i,j)$, a well-known invariant of braids. The number
$\mathsf{lk}(i,j)$ is the number of times that $i$ crosses over $j$ minus the
number of times $j$ crosses over $i$.

Figure 2. The Borromean rings, a non-trivial 3-link with all pairwise linking
numbers equal to zero.

Intuition might lead us to believe that if the vector
$\\{\mathsf{lk}(i,j)\\}\in\mathbb{Z}^{N}$ of the linking numbers of each pair
of trajectories is the zero vector, then the braid is trivial, but this is
certainly not the case! The most famous example of a non-trivial $3$-link
having all pairwise linking numbers equal to zero is the Borromean rings,
shown in figure 2. Thankfully, we can generalize the notion of a linking
number to describe the entanglement of larger groups of strands and detect
finer structures. The generalization to $m$ strands is known as the $m$-th
Milnor number, defined as follows:

######  Definition 1 (see e.g. [MK99]).

The $m$-th Milnor number or _generalized linking number_ is defined for the
pure braid $\beta$ on $m$ strands $\\{p_{1},\dots,p_{m}\\}$. Let $w$ be a word
corresponding to the braid $\beta$ in the free group presentation. Then $w$
can be written as $S_{p_{m}}p_{m}S_{p_{m}}^{-1}$. For any subsequence
$d=t_{j_{1}}^{\epsilon_{j_{1}}}\cdots t_{j_{m-1}}^{\epsilon_{j_{m-1}}}$ inside
$S_{p_{m}}$ such that $j_{1}<\dots<j_{m-1}$ and $t_{j_{1}}=p_{1}$,
$\dots,t_{j_{m-1}}=p_{m-1}$, define the sign of the subsequence as
$\mathsf{sgn}(d)=\epsilon_{j_{1}}\cdots\epsilon_{j_{m-1}}$. Then the $m$-th
Milnor number for the braid $\beta$ is defined as the sum of $\mathsf{sgn}(d)$
over all sub-sequences of $S_{j_{m}}$.

In particular, $\mu_{2}(i,j)=\mathsf{lk}(i,j)$. Informally, $\mu_{m}$ detects
whether a set of loops forms a Brunnian link, a structure in which $m$ loops
are linked, such that any $(m-1)$ of them are unlinked. For example, for the
Borromean rings, $\mu_{2}(1,2)=\mu_{2}(2,3)=\mu_{2}(1,3)=0$ and
$\mu_{3}(1,2,3)=1$.

##  3\. The Case of $n=3$.

Consider the case of $n=3$ particles. The symmetric group $\mathfrak{S}_{3}$
has $6$ elements as illustrated in Fig. 1. The Cayley graph of
$\mathfrak{S}_{3}$ is therefore a regular hexagon contained by the plane
$x_{1}+x_{2}+x_{3}=1+2+3$. On this graph, the set of the reduced words for the
long element is generated by ${\bf 1}=(1,2)$ and ${\bf 2}=(2,3)$,
corresponding to reflections in the planes $p_{1,2}:\,x_{1}-x_{2}=0$ and
$p_{2,3}:\,x_{2}-x_{3}=0$, respectively. These two transpositions generate
exactly two different sorting networks: $S_{1}={\bf 121}$ and $S_{2}={\bf
212}$, connecting the identity element (corresponding to the vertex $(123)$ )
to the longest element (corresponding to the vertex $(321)$ ).

Up to the obvious symmetries, we can construct two different sorting loops:
$L_{1}=S_{1}S_{2}$ and $L_{2}=S_{1}S_{1}$, each having $2^{6}=64$ possible
signatures. In a given loop, denote the first and last crossings of the $i$-th
and $j$-th particles by $f(i,j)$ and $\ell(i,j)$. For the loop $L_{2}$ we have
$\ell(i,j)+f(i,j)=7$, while for the loop $L_{1}$ we have that
$\ell(i,j)-f(i,j)=3$. The linking number of the $i$-th and $j$-th strands only
depends on the signatures of $f(i,j)$ and $\ell(i,j)$; in particular,
$|\mathsf{lk}(i,j)|=\frac{1}{2}|\sigma_{f(i,j)}+\sigma_{\ell(i,j)}|$. It
follows that out of $64$ possible signatures, only $8$ provide a braid with
the vector $\vec{\mathsf{lk}}=\vec{0}$. For the case of $n=3$, the vector
$\vec{\mu}_{3}$ consists of one number: $\mu_{3}(1,2,3)$. We will now consider
the two loops $L_{1}$ and $L_{2}$ separately.

###  3.1. Sorting loop $L_{1}$

Out of $8$ possible signatures for the first half of the sorting loop,
precisely $6$ provide $\mu_{3}(1,2,3)=0$. The two signatures which provide a
non-trivial braid are $(+-+)$ and $(-+-)$ (see Fig. 3A. The explanation for
this lies in the inversion sets corresponding to each permutation.

Recall the definition of the inversion set:

######  Definition 2.

For any permutation $p\in\mathfrak{S}_{n}$ we can associate the inversion set,
a vector $\mathsf{inv}(p)\in\\{\pm 1\\}^{N}$ where

$\mathsf{inv}(p)_{i,j}=\begin{cases}\phantom{-}1&\mbox{ if }p(i)>p(j)\\\
-1&\mbox{ if }p(i)<p(j)\end{cases}.$

Each permutation is uniquely determined by its inversion set (see Fig.1). A
permutation $p$ can be thought of as a new total ordering on the set of $n$
elements. Because of this, not all of the $2^{N}$ possible signatures will
occur as an inversion set of some permutation, but only the $n!$ of these
signatures corresponding to the $p$-determined order relation, $>_{p}$.

A total order such as the one we have defined must satisfy transitivity
property: i.e., from $a>_{p}b>_{p}c$ it should follow that $a>_{p}c$. The
signatures $(+-+)$ and $(-+-)$ listed above are precisely those which violate
transitivity. Any other signature corresponds to a well-defined total order
relation on the strands; for example the signature $(+--)$ corresponds to the
ordering $3>1>2$ as shown in Fig. 1. Using this signature, we could draw the
wiring diagram of $L_{1}$ on three parallel planes, each containing exactly
one strand with no intersections.

In general for the signature given by the inversion set it is possible to
establish a total ordering on the set of strands. The total ordering implies
the parallel planes construction from our example, and so the resulting $n!$
sorting braids are trivial. The two signatures $(+-+)$ and $(-+-)$ do not
belong to the inversion set and so do not provide a complete ordering on the
strands. In fact, the closures of the pure braids resulting from these
signatures both create the Borromean rings (see Fig. 3A). Later on we will
call such loop to be of _type I_.

$3$$3$$2$$2$$1$$1$ A $\mathfrak{A}({\bf 121},{(+-+)})$

$3$$3$$2$$2$$1$$1$ B $\mathfrak{D}({\bf 121},{(+-+)})$

Figure 3. Sorting braids with zero vector of linking numbers.

###  3.2. Sorting loop $L_{2}$.

Our argument for $L_{1}$, it turns out, is independent of the sorting loop.
The six signatures we mentioned above for $L_{1}$ also provide trivial braids
for $L_{2}$. The two remaining signatures, however, act quite differently–they
lead to trivial braids as well (see Fig. 3B)! The reason for triviality is
much simpler: the word
$\mathbf{1\bar{2}1\bar{1}2\bar{1}}=\mathbf{1\bar{2}2\bar{1}}=\mathbf{1\bar{1}}=e$
reduces to the trivial one. Later we will refer to such loops as being of
_Type II._

##  4\. The General Case.

Now armed with some experience in the case of $n=3$, we proceed to examine the
general case for $n$ particles. Recall the notation
$N=\binom{n}{2}=\frac{n(n-1)}{2}$.

###  4.1. The Asymptotics of linking numbers.

We begin with a general statement on the signatures for a given sorting loop
$ST$:

######  Theorem 1.

Out of all $2^{2N}$ signatures for a sorting loop $ST$, there are exactly
$\binom{N}{k}2^{N}$ signatures for which $|\vec{\mathsf{lk}}|=k$.

###### Proof.

For any pair of particles $(p_{1},p_{2})$ there exist exactly two indices
corresponding to the crossing of their trajectories, namely
$f(p_{1},p_{2})\in[1,N]$ and $\ell(p_{1},p_{2})\in[N,2N]$. It follows that
$\mathsf{lk}(p_{1},p_{2})=\frac{1}{2}(\sigma_{f}+\sigma_{\ell})$. Furthermore,
$\mathsf{lk}(p_{1},p_{2})=0$ if $\sigma_{f}=-\sigma_{\ell}$ and has magnitude
$1$ otherwise. This implies that for any choice of the signatures
$\sigma_{1},\dots,\sigma_{N}$ of the first half of the loop, there exist
exactly $\binom{N}{k}$ ways we may choose signatures on the second part to
provide a non-zero vector of linking numbers. ∎

This theorem has the interesting effect of splitting the set of all signatures
with respect to the magnitude of the signature’s corresponding
$\vec{\mathsf{lk}}$: $2^{2N}=2^{N}\sum\limits_{k=0}^{N}\binom{N}{k}$.

We can also compute the average linking number for a sorting loop over all
possible signatures. The computation follows:

$\mathbb{E}(\mathsf{lk})=2^{-2N}\sum\limits_{k=0}^{N}k\binom{N}{k}2^{N}=2^{-N}\sum\limits_{k=0}^{N}k\binom{N}{k}=2^{-N}N2^{N-1}=\frac{N}{2}.$

###  4.2. Third order Milnor invariant.

As it turns out, we can make many general statements about the $\mu_{3}$
invariants corresponding to signatures of much larger sorting networks. We
begin with a definition:

######  Definition 3.

For the sorting network $S={\bf j}_{1}{\bf j}_{2}\cdots{\bf j}_{N}$ define the
conjugate sorting network as $S^{*}={\bf j}_{N}{\bf j}_{N-1}\cdots{\bf
j}_{1}$.

Conjugation of a sorting network is a well-defined involution. Intuitively,
$S^{*}$ corresponds to the sorting network that describes the path $S$ on the
permutahedron, but read back from the long element to the identity. We can
extend our definition of conjugation to apply to signed sorting networks as
well:

$({\bf j}_{1}^{\sigma_{1}}{\bf j}_{2}^{\sigma_{2}}\cdots{\bf
j}_{N}^{\sigma_{N}})^{*}:={\bf j}_{N}^{-\sigma_{N}}{\bf
j}_{N-1}^{-\sigma_{N-1}}\cdots{\bf j}_{1}^{-\sigma_{1}}.$

From our construction, we can immediately deduce the following lemma:

######  Lemma 1.

For any signed sorting network $S(\sigma)$, the sorting braid $SS^{*}$ is
trivial.

###### Proof.

The word ${\bf j}_{1}^{\sigma_{1}}{\bf j}_{2}^{\sigma_{2}}\cdots{\bf
j}_{N}^{\sigma_{N}}{\bf j}_{N}^{-\sigma_{N}}{\bf
j}_{N-1}^{-\sigma_{N-1}}\cdots{\bf j}_{1}^{-\sigma_{1}}$ is trivial since the
product ${\bf j}_{N}^{-\sigma_{N}}{\bf j}_{N-1}^{-\sigma_{N-1}}\cdots{\bf
j}_{1}^{-\sigma_{1}}$ is exactly the inverse of the product ${\bf
j}_{1}^{\sigma_{1}}{\bf j}_{2}^{\sigma_{2}}\cdots{\bf j}_{N}^{\sigma_{N}}$. ∎

###  4.3. $\mu_{3}$ invariants.

According to the theorem by Tits, any two sorting networks are connected in
$\mathcal{S}(n)$ by a sequence of commutations $\bf ij\sim ji$, $|i-j|>1$ and
braid moves $\bf jij\mapsto iji$, $|i-j|=1$ (see e.g. [RR13] and references
therein). From the sequence of commutations we can construct the following
definition:

######  Definition 4.

For two sorting networks $S,T\in\mathcal{S}(n)$, define $\langle S,T\rangle$
to be the number of braid moves in the minimal sequence connecting $S$ and
$T$. Note that $\langle S,T\rangle=\langle T,S\rangle$.

This definition allows us to connect the construction of signed braids to the
$n$-permutahedron. Similarly to before, let $\vec{\mu}_{3}$ be the vector of
$\mu_{3}$ invariants for a particular sorting braid.

######  Theorem 2.

For $S,T\in\mathcal{S}(n)$, the loop $ST$ corresponds to $2^{\langle
S,T\rangle}$ unlinked braids for which $|\vec{\mu}_{3}|\neq 0$.

###### Proof.

The value $\langle S,T\rangle$ is by definition the minimal number of
hexagonal faces bounded by the loop $ST$ on the $n$-permutahedron. Each
hexagon represents a sorting 3-subbraid of type I. As discussed in section 3,
exactly two signatures for every hexagon of type I provide a non-zero
$mu_{3}$. The theorem follows directly. ∎

Theorem 2 gives us a rough picture of the set of signatures providing non-
trivial $\vec{\mu}_{3}$; however, it does not detect the finer details of the
distribution of the magnitudes of $\vec{\mu}_{3}$ vectors. This is a
consequence of $\mu_{3}(p,r,s)$ and $\mu_{3}(p,q,r)$ not being independent.

###  4.4. Higher Milnor invariants.

As defined before, Milnor invariants can be generalized to an arbitrary number
of strands to detect entanglement (or more precisely, Brunnian arrangements in
the closure of a braid). Because of this, it is natural to ask what can be
said about higher linking invariants for our sorting braids. To this end, we
will prove the following theorem:

######  Theorem 3.

If a sorting braid ST has all $\mu_{3}$ invariants equal to zero, the braid is
trivial.

The proof relies on the following lemma:

######  Lemma 2.

If a sorting braid $ST$ has all $\mu_{3}$ invariants equal to zero, then it is
isotopic to the braid $T^{*}T$.

Figure 4. Braid move for the signed braid

###### Proof.

For our signed sorting network, define a commutation as ${\bf
j}_{1}^{\sigma_{1}}{\bf j}_{2}^{\sigma_{2}}\to{\bf j}_{2}^{\sigma_{2}}{\bf
j}_{1}^{\sigma_{1}}$ and a braid move as ${\bf j}_{1}^{\sigma_{1}}{\bf
j}_{2}^{\sigma_{2}}{\bf j}_{1}^{\sigma_{3}}\to{\bf j}_{2}^{\sigma_{3}}{\bf
j}_{1}^{\sigma_{2}}{\bf j}_{2}^{\sigma_{1}}$.

Let $ST$ represent a sorting loop. By definition, there exists a path
$S=S_{0}\to S_{1}\to\cdots\to S_{m}=T^{*}$ in $\mathcal{S}(n)$ such that each
step is either a braid move or a commutation. Since all $\mu_{3}$ are zero by
assumption, the signature for every triple involved in the braid move must
belong to the inversion set. Each of our braid moves, then, can be extended to
the signed sorting networks $S(\sigma)$ and $T^{*}(\sigma^{\prime})$ as well.
This directly implies that the same minimal path leads from $S(\sigma)$ to
$T^{*}(\sigma^{\prime})$. As shown in figure 4, every step in the sequence is
a braid isotopy. ∎

The claim of Theorem 3 now follows from lemmas 2 and 1.

##  5\. Two Particular Constructions.

For a signed sorting network, we are free to choose several rules for defining
the second half of a sorting loop. In this section, we will describe two
particular constructions that come from seemingly natural choices of the
second half.

###  5.1. Algebraic sorting braids.

When considering sorting networks in terms of the $n$-permutahedron, a natural
geometric approach to constructing a sorting loop is to approximate an
equator. We will call the pure braid arising from this construction an
algebraic sorting braid, which we will define as follows:

######  Definition 5.

An _Algebraic sorting braid_ (ASB) corresponding to the sorting network
$S={\bf j}_{1}^{\sigma_{1}}\cdots{\bf j}_{N}^{\sigma_{N}}\in\mathcal{S}(n)$ is
constructed as

$\mathfrak{A}(S,\sigma)={\bf j}_{1}^{\sigma_{1}}\cdots{\bf
j}_{N}^{\sigma_{N}}({\bf n-j}_{1})^{-\sigma_{1}}\cdots({\bf
n-j}_{N})^{-\sigma_{N}}.$

We can immediately prove that

######  Lemma 3.

Every algebraic sorting braid is unlinked.

###### Proof.

For any pair of particles $(p_{1},p_{2})$ such that $j_{k}=f(p_{1},p_{2})$ it
follows that ${\bf j}_{k}$ corresponds to the crossing of the $p_{1}$-st and
$p_{2}$-nd strands in the wiring diagram. It follows that $({\bf n-j}_{k})$ in
the product $({\bf n-j}_{1})\cdots({\bf n-j}_{N})$ corresponds to the crossing
of $(n+1-p_{1})$-st and $(n+1-p_{2})$-nd strands. These strands, though, are
exactly those corresponding to the trajectories of the particles at the
positions $(n+1-p_{1})$ and $(n+1-p_{2})$ after applying ${\bf
j}_{1}\cdots{\bf j}_{N}$. Since the latter product corresponds to the sorting
network, the particles are $p_{1}$ and $p_{2}$; therefore
$\ell(p_{1},p_{2})=N+k$. Since $\sigma_{N-k}=-\sigma_{k}$, we have that
$\mathsf{lk}(p_{1},p_{2})=\frac{1}{2}(\sigma_{k}-\sigma_{k})=0$. ∎

We can also easily demonstrate the following:

######  Theorem 4.

Out of all $2^{N}$ signatures of the sorting network $S$ there are exactly
$n!$ which provide a trivial ASB.

###### Proof.

Any $m$-subbraid of an algebraic sorting braid is algebraic sorting braid on
$m$ particles. Therefore, any $3$-subbraid is a sorting $3$-braid of type $I$
(see Fig. 3A). In order for the whole braid to be trivial, all $\mu_{3}$
invariants must to be zero. Therefore all the signatures corresponding to the
$3$-braids must to belong to the inversion set and therefore provide an
ordering on respective trajectories. Since any two trajectories intersect
exactly once in the first half of a sorting braid, the orderings on two
different triples $\\{i,j,k\\}$ and $\\{i,j,k^{\prime}\\}$ must match at the
intersection $\\{i,j\\}$. Hence, the order relation extends to a total order
on the whole set of trajectories. Since there are exactly $n!$ different total
orderings on the set of $n$ elements, the theorem follows. ∎

Since any triplet of strands in ASB form a sorting loop of type $I$ and
signatures providing inversion sets for the triplets $\\{i,j,k\\}$ and
$\\{i,j,k^{\prime}\\}$ conditioned to the relative order of $i$ and $j$ are
independent, the averaged magnitude of the third Milnor invariant of any ASB
equals to

$\mathbb{E}|\mu_{3}(\mathfrak{A}(S))|=\frac{1}{4}\binom{n}{3}.$

###  5.2. Dynamic sorting braid

Another natural choice of the second half of a sorting loop is motivated by a
dynamical interpretation of the word representation of the sorting network. We
can construct a loop from a word on elementary transpositions simply by
iterating the word a second time. Since this second part of the loop will act
on the reversed interval, is seems natural to reverse the signatures of each
transposition in the second part of the loop. We will define the resulting
pure braid as follows:

######  Definition 6.

A _Dynamic sorting braid_ (DSB) is obtained from the signed sorting network
$S(\sigma)={\bf j}_{1}^{\sigma_{1}}\cdots{\bf j}_{N}^{\sigma_{N}}$ by the
concatenation

$\mathfrak{D}(S,\sigma)=S(\sigma)S({-\sigma})={\bf
j}_{1}^{\sigma_{1}}\cdots{\bf j}_{N}^{\sigma_{N}}{\bf
j}_{1}^{-\sigma_{1}}\cdots{\bf j}_{N}^{-\sigma_{N}}$

(see Fig. 3B)).

Our definition of the DSB leads us to the following lemma:

######  Lemma 4.

Let $f(p_{1},p_{2})$ be the index of the transposition corresponding to the
first crossing of the trajectories of the particles $p_{1}$ and $p_{2}$. If
the signatures of the transpositions $f(p_{1},p_{2})$ and
$f(n+1-p_{1},n+1-p_{2})$ do not coincide, then $|\mathsf{lk}(p_{1},p_{2})|=1$.

###### Proof.

If $f(p_{1},p_{2})=k$, it follows that the transposition ${\bf j}_{k}$ in the
word $S={\bf j}_{1}\cdots{\bf j}_{N}$ corresponds to the intersection of the
trajectories of the particles $p_{1}$ and $p_{2}$. Hence, the $(N+k)$-th
transposition in the sorting loop $SS$ corresponds to the second intersection
of the trajectories of the particles $(n+1-p_{2})$ and $(n+1-p_{1})$. By the
definition of the dynamic sorting braid we have $\sigma_{N+k}=-\sigma_{k}$.
Therefore,
$|\mathsf{lk}(n+1-p_{1},n+1-p_{2})|=|\sigma_{f(n+1-p_{1},n+1-p_{2})}-\sigma_{k}|$.
Since $k=f(p_{1},p_{2})$ the lemma follows. ∎

This lemma implies

######  Theorem 5.

Given $S\in\mathcal{S}(n)$ there are $2^{\lfloor n/2\rfloor\lceil n/2\rceil}$
signatures $\sigma$ corresponding to the unlinked DSB.

###  5.3. Slim networks.

Unlike the case of the ASB, the amount of signatures providing a trivial DSB
depends on the particular sorting network. It is natural to ask if there
exists a network for which this amount attains its extremal values. This lead
us to the following definition:

######  Definition 7.

We will call a sorting network $S$ to be _slim_ if for any signature $\sigma$,
the braid $\mathfrak{D}(S,\sigma)$ is trivial.

The simplest example of a slim sorting network is the construction
$S_{0}=\bf(12\cdots(n-1))\,(1\cdot 2\cdots(n-2))\,\cdots(12)\,1$. In fact,
$S_{0}^{2}$ corresponds to the _full twist_ element in the braid group. We
think it instructive to show for this particular example that we can achieve
$S_{0}^{*}={\bf 1(21)(321)\cdots((n-1)\cdots 21)}$ from $S_{0}$ using only
commutations. The process is iterative:

  1. (1)

Start with the second block $\bf(12\cdots(n-2))$ in $S_{0}$. The generator
$\bf 1$ commutes with everything but the second letter in the block left to
it. Therefore, we obtain

$\bf(1(21)(34\cdots(n-1))(23\cdots(n-2))\cdots(12)1.$

Likewise, the generator $\bf 2$ in the second block commutes with everything
but the $\bf 3$ on the left. We get

$\bf(1(21)(32)(4\cdots(n-1)))(3\cdots(n-2))\cdots(12)1$

etc. Finally we obtain

$\bf(1(21)(32)(43)\cdots((n-1)(n-2)))(12\cdots(n-3))\cdots(12)1.$

  2. (2)

Now perform the same operations with the block $\bf(12\cdots(n-3))$: Starting
with the first generator, move it by commutations to the leftmost possible
position, then move to the second generator and move it to the leftmost
position achievable by commutations, etc. We get

$\bf(1(21)(321)(432)\cdots((n-1)(n-2)(n-3)))(12\cdots(n-4))\cdots(12)1.$

  3. (3)

Continuing this procedure with each block we obtain $S_{0}^{*}$.

This construction may be generalized to provide a family of slim sorting
networks. However, we were unable to find a defining property of a sorting
network to be slim obstructing us from the complete description of the set of
slim networks. As such, we are left with two lingering questions about the
nature of this set:

######  Question 1.

How big is the set of slim sorting networks?

One can also come up with the construction of somewhat opposite to the slim
network, i.e. the network providing the DSB with the most magnitude of the
vector of Milnor invariants (see Fig. 6). Therefore we can define the
_fatness_ of a sorting network as the maximal magnitude of $\vec{\mu}_{3}$ for
$\mathfrak{D}(S,\sigma)$ over all signatures $\sigma$. This measures the area
on the permutahedron bounded by the loop $SS^{*}$, leading us to our second
question:

######  Question 2.

How “fat” is the typical sorting network, and more generally what does the
distribution of fatness look like across all sorting networks?

###  5.4. Experimental results.

According to a theorem by Stanley, the cardinality of the set $\mathcal{S}(n)$
grows immensely fast.

######  Theorem 6 ([Sta84]).

Set $\mathcal{S}(n)$ has exactly

$\frac{\binom{n}{2}!}{1^{n-1}3^{n-2}5^{n-3}\cdots(2n-3)^{1}}$

elements.

Because of this, we were only able to directly analyze results for the first
few values of $n$. Below we examine the set $\mathcal{S}(4)$ in detail, and we
provide numerical results obtained for the set $\mathcal{S}(5)$.

####  5.4.1. The case of $n=4$

The set $\mathcal{S}(4)$ contains sixteen distinct sorting networks, which
form four equivalence classes up to obvious symmetries:

$\displaystyle S_{0}=\bf 123121,$ $\displaystyle S_{1}=\bf 123212,$
$\displaystyle S_{2}=\bf 132312,$ $\displaystyle S_{3}=\bf 132132.$

From theorem 5, it follows that sixteen signatures for each of these sorting
networks provides an unlinked DSB. Among them, $S_{0}$ is a slim network for
which no signatures provide a non-trivial linked braid (see Fig. 5A).

A $S_{0}$

B $S_{1}$

C $S_{2}$

D $S_{3}$

Figure 5. Loops on $4$-permutohedron, corresponding to the Dynamical Sorting
Braids.

The DSB obtained from the sorting network $S_{1}$ encircles two hexagons (see
Fig. 5B). In other words, $S_{1}^{*}$ can be obtained from $S_{1}$ by two
braid moves: $\bf 123212\mapsto 123121\sim 121321\mapsto 212321$, where $\sim$
denotes commutations and $\mapsto$ denotes braid moves. It follows from
theorem 2 that there are $4$ signatures which provide a nontrivial unlinked
DSB. Interestingly, each of these braids has $|\vec{\mu}_{3}|=2$. For example,
$\bf{1\bar{2}32\bar{1}2}$ provides the pure braid with
$|\mu_{3}(1,2,3)|=|\mu_{3}(2,3,4)|=1$ (see figure 6). The expected $\mu_{3}$
invariant for a signed DSB based on $S_{1}$ is thus $1/2$.

$4$$4$$3$$3$$2$$2$$1$$1$ Figure 6. $\mathfrak{D}(S_{3},(+-++-+))$ provides
$|\mu_{3}(1,2,3)|=|\mu_{3}(1,3,4)|=1$.

For both $S_{2}$ and $S_{3}$, the corresponding sorting loop encircles four
hexagons (see Figs 5C and 5D). Theorem 2 ensures the existence of $16$
different signatures leading to the non-trivial braid. However, not all these
signatures correspond to the dynamical sorting braid. Recall, that DSB
definition implies symmetry conditions on the signatures on the first and
second halves of the sorting loop. For instance, contrary to $64$ possible
signatures providing unlinked braid in general setting, only $16$ of them
provide unlinked DSB (see theorem 5). Out of these $16$ there are $8$
providing nontrivial unlinked DSB, with $|\vec{\mu}_{3}|=2$. The expected
$\mu_{3}$ invariant for these networks is $1$, so the expected invariant over
all unlinked DSB is $5/8$.

####  5.4.2. The case of $n=5$

$\mathcal{S}(5)$ contains $768$ distinct sorting networks, each having $64$
signatures providing an unlinked DSB. Below, we present a distribution of the
magnitudes of $\vec{\mu}_{3}$ averaged over all signatures for each network.
These averaged quantities happen to belong to the set $\\{0,1/2,1,3/2,2\\}$
with majority being $0$. The expected $\mu_{3}$ invariant for $\mathcal{S}(5)$
is $113/128$.

Figure 7. Distribution of the magnitudes of the third Milnor invariants among
all $768$ unlinked DSB on $5$ strands.

The experimental results lead us to numerous conjectures. In particular, it
seems plausible that:

  * •

Set of slim networks occupies non-zero measure in the set of all sorting
networks, i.e. $\frac{|\mbox{slim}|}{|\mathcal{S}(n)|}\to\varepsilon\neq 0$ as
$n\to\infty$.

  * •

Expected magnitude of the vector of Milnor invariants tends to some constant
$\mathbb{E}|\mu_{3}|\to\mu<\infty$ as $\nLeftarrow\to\infty$.

### Acknowledgements.

The authors are thankful to N. Stein for the beautiful COCALC environment
where most of our numerical studies were conducted, as well as to N. Williams
for deep and fruitful discussions. Part of work of C.K. was completed during
his summer spent participating in Budapest Semesters in Mathematics.

## References

  * [ADHV19] Omer Angel, Duncan Dauvergne, Alexander E. Holroyd, and Bálint Virág.  The local limit of random sorting networks.  Ann. Inst. H. Poincaré Probab. Statist., 55(1):412–440, 02 2019\.
  * [AHRV07] Omer Angel, Alexander E. Holroyd, Dan Romik, and Bálint Virág.  Random sorting networks.  Advances in Mathematics, 215(2):839 – 868, 2007.
  * [Arn86] V. I. Arnold.  The asymptotic Hopf invariant and its applications.  Selecta Math. Soviet., 5(4):327–345, 1986.  Selected translations.
  * [Bre08] Yann Brenier.  Generalized solutions and hydrostatic approximation of the euler equations.  Physica D: Nonlinear Phenomena, 237(14):1982 – 1988, 2008.  Euler Equations: 250 Years On.
  * [MK99] Kunio Murasugi and Bohdan I. Kurpita.  A Study of Braids.  Mathematics and Its Applications. Springer, 1999.
  * [PV94] Michael Polyak and Oleg Viro.  Gauss diagram formulas for vassiliev invariants.  International Mathematics Research Notices, 11, 1994.
  * [RR13] Victor Reiner and Yuval Roichman.  Diameter of graphs of reduced words and galleries.  Transactions of the American Mathematical Society, 365(5):2779 – 2802, 2013.
  * [RVV16] Mustazee Rahman, Balint Virag, and Mate Vizer.  Geometry of Permutation Limits.  arXiv e-prints, page arXiv:1609.03891, Sep 2016.
  * [Sta84] Richard Stanley.  On the number of reduced decompositions of elements of coxeter groups.  European Journal of Combinatorics, 5(4):359–372, 1984.