# SARS-CoV-2 Result Interpretation based on Image Analysis of Lateral Flow
Devices

Neeraj Vashistha

###### Abstract

The widely used gene quantisation technique, Lateral Flow Device (LFD), is now
commonly used to detect the presence of SARS-CoV-2. It is enabling the control
and prevention of the spread of the virus. Depending on the viral load, LFD
have different sensitivity and self-test for normal user present additional
challenge to interpret the result. With the evolution of machine learning
algorithms, image processing and analysis has seen unprecedented growth. In
this interdisciplinary study, we employ novel image analysis methods of
computer vision and machine learning field to study visual features of the
control region of LFD. Here, we automatically derive results for any image
containing LFD into positive, negative or inconclusive. This will reduce the
burden of human involvement of health workers and perception bias.

Keywords : SARS-CoV-2, Lateral Flow Device, Machine Learning Algorithms,
Computer Vision.

##  1 Introduction

The most common test to detect the presence of SARS-Cov2 are rapid antigen
tests using LFD and real-time quantitative polymerase chain reaction(RT-qPCR).
However, rapid antigen tests are the most cost-effective measure after the RT-
qPCR test. LFD using rapid antigen test provides users with a simple, binary
result (positive, negative) in as little as 5 to 20 minutes. These tests have
been widely proliferated since the beginning of SARS-Cov2 and are commonly
used to test pregnancy, malaria, H5N1, SARS, antibiotic resistance and
legionella [WQHX16] The LDF contain a strip that is coated with a specific
antigen on which human fluid (blood, plasma, mucus or urine) along with a
reagent is placed. This strip contains two bands, (control band and test band)
which appear after 5 mins to 20 mins. [vAVAK18]

Although, RT-qPCR is considered the more accurate but the cost associated with
it considering human efforts and equipment is very high. Rapid antigen tests
on the other hand can detect viruses (both asymptomatic and symptomatic) at
par with RT-qPCR tests[LRP+21]. Rapid antigen tests are simple to use,
relatively cheap and require minimal training. Both type of tests have helped
in the suppression and spread of viruses within the UK boroughs [PVZA+21]. Two
studies in 2020, showed LFD results when compared with PCR test showed the
specificity is almost 99.9% with sensitivity in the test to be around 48.9% to
78.9%[PAA+21] and LFD positivity is directly proportional with viral
load[GFHC+21].

A potential issue with the use of LFD is interpreting the test
results[MHJC20]. As previously mentioned, if the count of antibodies is less,
the band colours in LFD are faintly visible. This sometimes tempts users to
falsely read a test as negative or sometimes the eye acuity reads the test
result wrong. This task of interpretation becomes more complex in very elderly
or young individuals, which can waste a lot of effort of mass testing and
virus suppression.

The first objective of this study is to find a simple and intuitive way to
determine the presence of the control band and the test band in the LFD result
image using state-of-art machine learning models. The second objective is to
scale this system both vertically and horizontally, which would mean that it
would work for multiple retail vendors which are providing Lateral Flow kits
in the market as well as can sustain a multitude of verification -requests at
the same time. Third, reduce the human dependence for some borderline results
which are too small to detect by using active learning techniques.

To achieve the above laid out objectives, we propose a multi-model machine
learning strategy. These machine learning models are built in three phases. In
the first phase, the models are trained on different collections of images of
the Lateral Flow test results, these images are curated from labs and
different experiments to represent real-world data. In the second phase, we
use the result models to test on actual real-world images and tune the
parameter settings of each of these models. This enables the models to
understand the distribution of data in a real-world scenario. In the third
phase, we use these tested models and deploy them to infer in real-time to
classify a Lateral Flow test image result into positive, negative or
inconclusive.

##  2 Model Design

The proposed multi-model machine learning framework is illustrated in Figure
1. Here, multiple arbitrary sized images containing Lateral Flow results were
processed in the data preparation stage to extract the Region-of-Interest(ROI)
part using Detectron2[WKM+19] instance segmentation based Mask R-CNN
model[HGDG18]. This results in a uniform fixed-size image extraction of the
test strip which is only a selected portion of the test cassette.

On these uniformed test strip image datasets, we train various classification
models and segmentation models, such as vendor recognition image classifier, a
binary and multi-class classifier with active learning component enabled. The
binary classification models provide a probabilistic result for each test
strip image as positive, negative or inconclusive. The multi-class classifier
provides a more detailed probabilistic result for each test strip image giving
results in terms of different classes such as Positive IGG, Positive IGM,
Positive IGG and IGM, Negative or inconclusive, this multi-class classifier
model can be tuned according to different vendors.

Figure 1:  Machine Learning Model Architecture, the left section in blue
represents the data preparation stage, the left-centre section in black
represents Training and model development stage, the centre section in red
represents testing phase, the right section in green represents Testing with
validation, hyperparameter optimisation and inference stage.

The segmentation models[Yak19] are a set of validation models for the above
classifier models. In computer vision, image segmentation is the process of
partitioning a digital image into multiple segments (sets of pixels, also
known as image objects). The goal of segmentation is to simplify and/or change
the representation of an image into something more meaningful and easier to
analyze. Image segmentation is typically used to locate objects and boundaries
(lines, curves, etc.) in images. To support the classification models, we use
the result of the segmentation models and integrate it to provide a conclusive
result. Since the segmentation model has to work in coordination with the
classification model, the parameter tuning of segmentation model and loss
function tuning is performed only after classification model parameters are
tuned.

After hyperparameter tuning of classification model and segmentation model is
completed, these models are optimised for scalability which enables them to be
run in parallel and inference at real time. Sometimes classification models
are sufficient if the probabilistic scores for a class are above a certain
threshold, in such case segmentation models are not inferred to save the
memory and inference time. Thus, each model is configurable and independent of
each other with certain constraints of the overall architect design.

##  3 Data

Data used in this work is curated manually at Medusa19 and Avacta labs, it
consists of only image samples and associated test results. It is the only
information required for the development, training, testing and validation of
machine learning models. Table 1, describe the overall data on which this
study was performed.

Type | Medusa19 | Avacta
---|---|---
Positive images | 306, 77.2% | 18, 10.4%
Negative images | 104, 22.1% | 153, 88.4%
Inconclusive images | 3, 0.6% | 2, 1.15%
Table 1: Training Samples for model building Figure 2:  Training Sample
Images, left Medusa19 with QR code, center Medusa19 w/o QR code and right
Avacta w/o QR code.

The raw images can come in various sizes and shapes with different orientation
but the default configuration are, portrait mode, with acceptable angle of
rotation of the test cassette from 0 to +-18 degrees. Figure 2, shows some of
the curated raw images.

##  4 Model Development

There are four models which effectively work in this system. The Instant
Segmentation model, Classification Models, Image Segmentation Model and Active
Learning model. The instant Segmentation model and classification model are
one of the most critical models, although independent of each other but these
need to work one after another, and loaded parallelly into the memory as the
result of the first is consumed in the latter. Image Segmentation Model and
Active Learning models improve the results of the other two models. In the
below section, we will describe how each model is built, used and its
performance metrics.

###  4.1 Instant Segmentation Model

This model is built on an instance segmentation data-set and trained on the
Detectron2 model. The main aim of this model is to perform instance
segmentation to extract test strips from test cassettes. Figure 3, showcase
this aim. To build a robust model, we curated images using different mobile
cameras, in different backgrounds, varying lighting conditions as well as
random objects in the background. After taking raw images, we transformed them
into a set resolution, 2582 x 1937, this is the most common resolution of a
modern 5MP camera where images are in 4:3 aspect ratio.

Figure 3:  Left raw Medusa19 test cassette image (2722 x 3607) and right test
strip extract (300 x 875)

After gathering the data, we use the Labelme[Wad] tool for annotation, as it
supports both VOC and COCO formats. The annotation built using this tool on
the images will be used to perform the instance segmentation. We transform our
annotations into COCO format which are easily trainable.

Figure 4:  Detectron2 COCO Instance Segmentation with Mask R-CNN model
performance over Lateral flow test cassettes.

To build this model we are using COCO Instance Segmentation with Mask R-CNN
model[Abd17], this is done to build a robust model by using a pre-trained
model which scores very highly, 41.0 on box Average Precision and 37.2 on mask
Average Precision and thus reduce the dependence on training a very large
dataset and while maintaining the performance of the model. Figure 4,
demonstrates the learning process of the model, as the loss decreases the
model accuracy to determine the presence of test strips from the test
cassettes increases. The model returns an impressive benchmark inference speed
score of processing an image at an average of 0.19 with 5.18 fps.

###  4.2 Binary and Multi Class Classification Model

For these models, the required input image dataset is the test strips which
are obtained from the previous Instance segmentation model. It becomes
critical that the right training data is fed into each model. The effect of
training data on loss function guides us through this. Binary Classifiers will
be trained with all images partitioned into two classes, positive and
unknown(negative and inconclusive).

Multi-Class Classifiers will be trained with all the images being segregated
into different classes such as Positive IGG, Positive IGM, Positive IGG and
IGM, Negative or inconclusive. The defined architecture has 5 output classes
with one neuron which provides the probability of each class.

Each image is of 300 x 875 resolution and there are a total of 337 training,
85 validation and 17 test images for binary and multi-class classification
models. Since we have a small size of data, we are using the Xception[Cho17]
pre-trained model’s weights and fine-tune[C+20] the existing model with our
new data. Below figure 5 describe the configuration of the model.

Figure 5:  Xception Model Architecture, the input embedding layer and end
neuron layers are modified according to our data and expected result. The data
first goes through the entry flow, then through the middle flow which is
repeated eight times, and finally through the exit flow. Note that all
Convolution and Separable Convolution layers are followed by batch
normalization (not included in the diagram). All Separable Convolution layers
use a depth multiplier of 1 (no depth expansion).

Figure 6:  Binary Classification model Loss over 30 epochs

In figure 6, the binary cross-entropy loss of the model can be seen to have
subtle variations on the validation set. After performing hyperparameter
optimisation, in each epoch, loss measure can be seen decreasing with some
exceptions. This implies that the model is learning but sometimes it has a
tough time to generalise on unseen data sets while predicting the classes. The
best weights are found at the 19th epoch and shown in Table 2.

|  Train set |  |
---|---|---|---
|  Validation set |  |
|  Test set evaluation |  |
| evaluation scores |  |
Binary cross entropy | 0.202241 | 0.240638 | 0.194755
Accuracy | 0.923630 | 0.912064 | 0.926810
F1 score mean | 0.921999 | 0.912423 | 0.921435
Mean Precision | 0.949316 | 0.937087 | 0.955327
Mean Recall | 0.905966 | 0.898664 | 0.902135

Table 2: Binary Classification model performance measures

The binary classification model has good performance on train, validation and
test datasets. The values of loss and metrics can be seen to be similar in
these datasets. This tells that the model is not overfitting on the dataset.
The f1 score of 0.921 on the validation dataset is acceptable.

Figure 7:  Multi-Class Classification Loss over 30 epochs

In figure 7, the categorical cross-entropy loss of the model is initially
struggling to converge over the validation dataset for the multi-class
classifier as we see a large variation in the 0-10 epoch. But it later
improves and generalisation is overall good.

|  Train set |  |
---|---|---|---
|  Validation set |  |
|  Test set evaluation |  |
| evaluation scores |  |
Categorical cross entropy | 0.081054 | 0.092119 | 0.094178
Accuracy | 0.968118 | 0.962500 | 0.965517
F1 score mean | 0.940510 | 0.929417 | 0.936398
Mean Precision | 0.945815 | 0.929264 | 0.941134
Mean Recall | 0.937232 | 0.931588 | 0.93854

Table 3: Multi-Class Classification model performance measures

The best weights were found at the 18th epoch and Table 4.2, shows the
performance of the model. In conclusion, the multi-class classification model
is generalizing well on the unseen data, as the values of the evaluation test
set and validation sets are closer to the train set.

###  4.3 Segmentation Model

The Segmentation models are used to find the presence of a certain type of
class. To prepare the dataset for this model, we use the VGG Image Annotator
(VIA)[DZ19] tool and create masks for all the images. We use the pre-trained
model, for Legendary UNet[RFB15] with the EfficientNetB1[TL20] backbone
model’s weight. For this model we are using the Dice Coefficient[Dic45] as the
performance metric, to gauge similarity of two samples. The Dice coefficient
can be used to compare the pixel-wise agreement between a predicted
segmentation and its corresponding ground truth. The formula is given by:

$C=\frac{2*|X{\displaystyle\cap}Y|}{|X|+|Y|}$

where X is the predicted set of pixels and Y is the ground truth. The Dice
coefficient is defined to be 1 when both X and Y are empty. Our objective is
to maximize the Dice coefficient and identify/locate the type of defect
present in the image.

Figure 8:  Positive-IGG class Segmentation model dice coefficient over epochs

In figure 8, the dice coefficient stabilizes around 65% after 15 epochs for
class 1 which is to detect Positive IGG. This score is good in terms of the
number of images for this class.

|  Train set |  |
---|---|---|---
|  Validation set |  |
|  Test set evaluation |  |
| evaluation scores |  |
Dice loss | 0.285742 | 0.334879 | 0.388797
Dice coefficient | 0.714258 | 0.665121 | 0.611203

Table 4: Positive-IGG class Dice performance scores

The Dice coefficient can be considered as an F1 score to evaluate the
performance metrics of the segmentation models. In Table 4.3, for class
Positive-IGG, an acceptable score of 0.665 is obtained on the validation data
set.

Figure 9:  Positive-IGM class Segmentation model dice coefficient over epochs.

In figure 9, the dice coefficient stabilizes around 65% at 20 epoch for class
2 which is to detect Positive IGM. This score is good in terms of the number
of images for this class.

|  Train set |  |
---|---|---|---
|  Validation set |  |
|  Test set evaluation |  |
| evaluation scores |  |
Dice loss | 0.295981 | 0.308811 | 0.338572
Dice coefficient | 0.766948 | 0.678812 | 0.655394

Table 5: Positive-IGM class Dice performance scores

The Dice coefficient can be considered as an F1 score to evaluate the
performance metrics of the segmentation models. In Table 4.3, for class
Positive-IGG, an acceptable score of 0.665 is obtained on the validation data
set.

Similarly, we calculate the Dice coefficient for Positive-IGG and IGM class
and Negative class and inconclusive class. Below table 4.3. summarises the
Dice coefficient for all the classes.

Dataset |  Positive IGG |  |  |  |
---|---|---|---|---|---
|  Positive IGM |  |  |  |
|  Positive IGG |  |  |  |
|  Inconclusive |  |  |  |
|  Negative |  |  |  |
| class |  |  |  |
X Train | 0.714258 | 0.766948 | 0.735519 | 0.67421 | 0.821943
X Validation | 0.665121 | 0.678812 | 0.709641 | 0.63597 | 0.76066
X Test | 0.611203 | 0.655394 | 0.698548 | 0.619753 | 0.78822

Table 6: Summary of Dice Coefficient on Unet - EfficientNetB1 model for
different classes on data splits.

###  4.4 Active Learning

This is an addition to the existing classification models. The role of active
learning[RM01] is to improve the overall performance of binary and multiclass
classifiers. Using the Active Learning framework we allow increasing
classification performance by intelligently querying the label for the most
informative instances.

Figure 10:  Active Learning on Binary and Multi-class classifiers. Souce:
modal

We are using the modAL framework[DH] to implement Active Learning, the working
is described in figure 10. The key components are choosing the uncertain
measure and query strategy for the labels. We are handling the first component
by asserting the probability score of each of the classes of binary and multi-
class classifiers at a set threshold. This threshold is derived through
rigorous testing of the overall system. The second component is handled using
a human-in-the-loop response or a pre-built framework of modAL.

##  5 Inference and Results

From the training above, the best models are saved to make inferences on
images. We define the different dependencies for loading saved models as it is
important while we use custom metrics. We define a custom method for
generating predictions using Classifier models and another custom method for
generating predictions using Segmentation models. Different thresholds such as
Area thresholds and Classification thresholds are applied to the predictions
of the models. Figure 11, shows sample evaluation on a single image. The image
has a negative class(no defect) and the models have perfectly detected that
image has a negative class(no defect).

Figure 11:  Negative class evaluation on Multi-class classifier and
Segmentation Models

Dataset |  Class 1 |  |  |  |
---|---|---|---|---|---
|  Class 2 |  |  |  |
|  Class 3 |  |  |  |
|  Class 4 |  |  |  |
|  Class 5 |  |  |  |
| Negative |  |  |  |
X Train | 0.9524 | 0.9823 | 0.8203 | 0.8724 | 0.9787
X Validation | 0.9479 | 0.9806 | 0.8162 | 0.8611 | 0.9738
X Test | 0.9442 | 0.9826 | 0.8047 | 0.8341 | 0.9787

Table 7: Mean Dice Coefficient class wise.

From the table 5, we can understand that each metric can be used to see the
fact that the model did not overfit on the train set, validated well on
validation dataset and generalising well in unseen test set.

Figure 12:  Overall evaluation scores on all the Training images on
Segmentation Model, Note: Some images are added by transformations and noise
insertion techniques.

In figure 12, we observe the classification report for the multi-class
classifier on augmented images. Due to data augmentation, the F1 scores are
low but overall precision and recall are good. The mean Dice coefficient for
overall classes is perfect and is around 93% and the mean Dice coefficient for
each class is above 85%. Thus, if the Multi-class classifier produces an
inaccurate result, these can be well covered up by using segmentation models.

##  6 Conclusion

In this study we have carried out an analysis of accuracy on machine learning
model to evaluate the results of lateral flow devices tests. The image reader
demonstrate a novel AL algorithm can perform and scale well across different
vendors with 84% sensitivity and 87% specificity. The model built, is self-
adaptive to learn new images with human-in-the-loop active learning model
present and the current values of sensitivity and specificity increases with
each new inference when validated by human. This overall improves the false
positive rates and maintains true positive rate as the results are identical
with each iteration with active learning. Another advantage is the processing
time is below 2 seconds and it remains constant with increase in load, thus it
scales vertically and horizontally. Additionally, due to the presence of
segmentation model, the model excellently detects faint color bands which are
due to low vial load and this likely help users who are unclear about reading
LFDs and are prone to make mistake likely or otherwise. One disadvantage is
the quality of image, we have considered the current trend of available
smartphone collect images above 5 mega pixels if the images are below this
metric or have issues such as blur, brightness, reflection or darkness, the
model would likely reject such images. Apart from this limitation, the model
works at a level of consistent performance with the capacity of processing
millions of images rapidly with feature to improve the ineffectualness with
active learning, thus evolving better every time, and reducing the user error
and increasing overall accuracy of model.

## References

  * [Abd17] Waleed Abdulla.  Mask r-cnn for object detection and instance segmentation on keras and tensorflow.  https://github.com/matterport/Mask_RCNN, 2017.
  * [C+20] Francois Chollet et al.  Transfer learning and fine-tuning, 2020.
  * [Cho17] François Chollet.  Xception: Deep learning with depthwise separable convolutions, 2017.
  * [DH] Tivadar Danka and Peter Horvath.  modAL: A modular active learning framework for Python.  available on arXiv at https://arxiv.org/abs/1805.00979.
  * [Dic45] Lee R. Dice.  Measures of the amount of ecologic association between species.  Ecology, 26(3):297–302, 1945.
  * [DZ19] Abhishek Dutta and Andrew Zisserman.  The VIA annotation software for images, audio and video.  In Proceedings of the 27th ACM International Conference on Multimedia, MM ’19, New York, NY, USA, 2019. ACM.
  * [GFHC+21] Marta García-Fiñana, David M Hughes, Christopher P Cheyne, Girvan Burnside, Mark Stockbridge, Tom A Fowler, Veronica L Fowler, Mark H Wilcox, Malcolm G Semple, and Iain Buchan.  Performance of the innova sars-cov-2 antigen rapid lateral flow test in the liverpool asymptomatic testing pilot: population based cohort study.  BMJ, 374:n1637, 2021.
  * [HGDG18] Kaiming He, Georgia Gkioxari, Piotr Dollár, and Ross Girshick.  Mask r-cnn, 2018.
  * [LRP+21] Lennard Y W Lee, Stefan Rozmanowski, Matthew Pang, Andre Charlett, Charlotte Anderson, Gareth J Hughes, Matthew Barnard, Leon Peto, Richard Vipond, Alex Sienkiewicz, Susan Hopkins, John Bell, Derrick W Crook, Nick Gent, A Sarah Walker, Tim E A Peto, and David W Eyre.  Severe acute respiratory syndrome coronavirus 2 (sars-cov-2) infectivity by viral load, s gene variants and demographic factors, and the utility of lateral flow devices to prevent transmission.  Clinical infectious diseases : an official publication of the Infectious Diseases Society of America, 74(3):407–415, 2021.
  * [MHJC20] Isabela A. Mattioli, Ayaz Hassan, Osvaldo N. Oliveira Jr., and Frank N. Crespilho.  On the challenges for the diagnosis of sars-cov-2 based on a review of current methodologies.  ACS Sensors, 5(12):3655–3677, 2020.
  * [PAA+21] Tim Peto, Dominic Affron, Babak Afrough, Anita Agasu, Mark Ainsworth, Alison Allanson, Katherine Allen, Collette Allen, Lorraine Archer, Natasha Ashbridge, Iman Aurfan, Miriam Avery, Ellena Badenoch, Priya Bagga, Rishab Balaji, Ella Baldwin, and Sophie Barraclough et al.  Covid-19: Rapid antigen detection for sars-cov-2 by lateral flow assay: A national systematic evaluation of sensitivity and specificity for mass-testing.  EClinicalMedicine, 36(100924), 2021.
  * [PVZA+21] Martin Pavelka, Kevin Van-Zandvoort, Sam Abbott, Katharine Sherratt, Marek Majdan, CMMID COVID-19 working group, Inštitút Zdravotných Analýz, Pavol Jarčuška, Marek Krajčí, and Stefan Flascheand Sebastian Funk.  The impact of population-wide rapid antigen testing on sars-cov-2 prevalence in slovakia.  Science (New York, N.Y.), 372(6542):635–641, 2021.
  * [RFB15] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation, 2015\.
  * [RM01] Nicholas Roy and Andrew McCallum.  Toward optimal active learning through monte carlo estimation of error reduction.  ICML, Williamstown, 2:441–448, 2001.
  * [TL20] Mingxing Tan and Quoc V. Le.  Efficientnet: Rethinking model scaling for convolutional neural networks, 2020.
  * [vAVAK18] A. van Amerongen, J. Veen, H. A. Arends, and M. Koets.  Lateral flow immunoassays in handbook of immunoassay technologies.  Academic Press, page 157–182, 2018.
  * [Wad] Kentaro Wada.  Labelme: Image Polygonal Annotation with Python.
  * [WKM+19] Yuxin Wu, Alexander Kirillov, Francisco Massa, Wan-Yen Lo, and Ross Girshick.  Detectron2.  https://github.com/facebookresearch/detectron2, 2019.
  * [WQHX16] Kan Wang, Weijian Qin, Yafei Hou, and Kun Xiao.  The application of lateral flow immunoassay in point of care testing: A review.  Nano Biomed. Eng., 8(3):172–183, 2016.
  * [Yak19] Pavel Yakubovskiy.  Segmentation models.  https://github.com/qubvel/segmentation_models, 2019.


varioSection #1 variothmTheorem #1 variolemLemma #1 variocorCorollary #1
variodefDefinition #1 varioobsObservation #1 varioasmAssumption #1 varioFig.
#1 varioappAppendix #1 variopropProposition #1 varioalgAlgorithm #1
variotabTable #1 vario(#1)

# Privacy-Preserving Wireless Federated Learning Exploiting Inherent Hardware
Impairments

Sina Rezaei Aghdam${}^{\text{1}}$, Ehsan Amid${}^{\text{2}}$, Marija
Furdek${}^{\text{1}}$, and Alexandre Graell i Amat${}^{\text{1}}$
The research in this paper has been supported by Chalmers Artificial
Intelligence Research Centre (CHAIR). ${}^{\text{1}}$Chalmers University of
Technology, Gothenburg, Sweden
${}^{\text{2}}$Google Research, Mountain View, CA

###### Abstract

We consider a wireless federated learning system where multiple data holder
edge devices collaborate to train a global model via sharing their parameter
updates with an honest-but-curious parameter server. We demonstrate that the
inherent hardware-induced distortion perturbing the model updates of the edge
devices can be exploited as a privacy-preserving mechanism. In particular, we
model the distortion as power-dependent additive Gaussian noise and present a
power allocation strategy that provides privacy guarantees within the
framework of differential privacy. We conduct numerical experiments to
evaluate the performance of the proposed power allocation scheme under
different levels of hardware impairments.

##  I Introduction

Federated learning has recently emerged as a promising technique for
distributed machine learning in a variety of applications, including mobile
edge computing [1]. Unlike conventional centralized machine learning
techniques that require datasets to reside on a single server, federated
learning leverages the local processing capabilities of edge devices and
requires only a limited amount of communication between the devices and the
server. In particular, instead of offloading the entire data, each user uses
its own dataset to train a local model and sends a small update to the global
model maintained by the server. This offers significant gains in terms of
privacy and communication efficiency [2].

When deployed over wireless networks, the superposition nature of wireless
channels can be used as a natural data aggregator for federated learning. This
phenomenon is referred to as _over-the-air aggregation_ , where simultaneously
transmitted analog waves from different edge devices are weighted by channel
coefficients and superposed at the server[3]. Accurate and fair computation of
the global model requires that the gradient estimates received from different
devices have the same power at the parameter server. One way to address this
issue is to perform channel-inversion power control at the edge devices [4,
5]. Alternatively, the parameter server can mitigate channel fading by using
multiple receiver antennas [6, 7].

Despite its various attractive features, wireless federated learning in its
primary form does not provide sufficiently strong privacy guarantees. Although
edge devices do not explicitly share their data in its original format,
several recent works have shown that model parameters can themselves be
informative to an honest-but-curious server (see, e.g., [8, 9, 10, 11, 12]),
leading to privacy leakage. Accordingly, the study of privacy for wireless
federated learning has been the subject of much recent research [13, 14, 15,
16]. These works adopt the notion of differential privacy as a formal model
for quantifying information disclosure about individuals and introduce
strategies for limiting it. Specifically, in [13, 14], the authors show that
by adding artificial noise to the model parameters on the edge device and
properly adjusting the variance of this noise, different levels of
differential privacy protection can be achieved. Moreover, [15, 16] show that
the inherent channel noise can be harnessed via proper power allocation to
achieve privacy for free.

In this paper, we demonstrate that in addition to channel noise, the
distortion introduced by the imperfect hardware of edge devices can also be
used as a privacy-preserving mechanism. Motivated by the analytical analysis
and experimental validations in [17, 18, 19], we model distortion as an
additive Gaussian noise whose variance is proportional to the signal power.
Within the differential privacy framework, we derive an upper-bound on the
privacy violation probability and use it for formulating a privacy
preservation condition. We then propose a distortion-aware power allocation
scheme to satisfy this condition. We conduct numerical experiments to evaluate
the performance of the proposed scheme under different levels of hardware
impairments. Our numerical results reveal that in realistic scenarios with
imperfect hardware, taking into account the impact of distortion in power
allocation can considerably reduce the performance loss caused by enforcing
differential privacy constraints. Exploiting distortion in realistic hardware
can even result in improved performance compared to ideal distortion-free
hardware.

The remainder of this paper is organized as follows. In Section II, we
introduce our system model and the required technical background, including
the definition of differential privacy. The derivation of the privacy
preservation criterion is described in Section III. In Section IV, we present
a power allocation scheme that satisfies the differential privacy criterion.
In Section V, we evaluate the performance of the proposed solution using
numerical examples. Finally, we conclude the paper in Section VI.

Notation: Throughout the paper, $\mathbb{R}$ and $\mathbb{C}$ represent the
sets of real and complex values, respectively. We denote a zero-mean normal
distribution with variance $\sigma^{2}$ by $\mathcal{N}(0,\sigma^{2})$. Zero-
mean complex Gaussian random variables are represented by
$\mathcal{CN}(0,\sigma^{2})$. The cardinality of a set $\mathcal{A}$ is
denoted by $|\mathcal{A}|$. Furthermore, the difference between two sets is
defined as
$\mathcal{A}^{\prime}-\mathcal{A}^{\prime\prime}=\\{x~{}|~{}x\in\mathcal{A}^{\prime},x\not\in\mathcal{A}^{\prime\prime}\\}$
and $\|\mathcal{A}^{\prime}-\mathcal{A}^{\prime\prime}\|_{1}$ denotes its
cardinality. The union of sets $\mathcal{A}_{1},\dots,\mathcal{A}_{K}$ is
represented by $\mathcal{A}=\cup_{k=1}^{K}\mathcal{A}_{k}$. Vectors are
denoted by boldfaced letters. The transpose and the Euclidean norm of a vector
$\bm{a}$ are denoted by $\bm{a}^{\mathsf{T}}$ and $\|\bm{a}\|$, respectively.

##  II System Model

We consider a wireless federated learning system (as shown in Fig. 1)
consisting of $K$ single-antenna edge devices that attempt to collaboratively
learn a shared model with the aid of a single-antenna parameter server. Each
device $k$ stores a local dataset, denoted by $\mathcal{B}_{k}$, which
consists of

$\mathcal{B}_{k}=\Big{\\{}\big{(}\bm{u}_{k,i},\upsilon_{k,i}\big{)}\Big{\\}}_{i=1}^{|\mathcal{B}_{k}|},$
(1)

where $\bm{u}_{k,i}$ denotes the $i$-th training sample at the $k$-th device
and $\upsilon_{k,i}\in\mathbb{R}$ is its corresponding label. The federated
learning process finds the optimal model parameters
$\bm{\theta}\in\mathbb{R}^{d}$ that minimize the global loss function
$F(\bm{\theta})$:

$\underset{\bm{\theta}}{\text{minimize}}~{}~{}F(\bm{\theta})\overset{\Delta}{=}\frac{1}{K}\sum_{k=1}^{K}F_{k}(\bm{\theta}).$
(2)

Here, $F_{k}(\bm{\theta})$ is the local loss function given by

$F_{k}(\bm{\theta})=\frac{1}{|\mathcal{B}_{k}|}\sum_{(\bm{u}_{k,i},\upsilon_{k,i})\in\mathcal{B}_{k}}f(\bm{\theta},\bm{u}_{k,i},\upsilon_{k,i}),$
(3)

and $f(\bm{\theta},\bm{u}_{k,i},\upsilon_{k,i})$ is the sample loss that
quantifies the prediction error of the model $\bm{\theta}$ on the training
samples $\bm{u}_{k,i}$ with respect to the labels $\upsilon_{k,i}$.

###  II-A Learning Protocol

In order to solve (2), a distributed stochastic gradient descent (SGD) is
adopted by the parameter server and the edge devices via the following
procedure:

  * •

Step 1: The parameter server broadcasts an initial model, i.e.,
$\bm{\theta}^{(1)}$, to the edge devices allowing them to initialize their
learning model.

  * •

Step 2: Each device performs local training on its own dataset and transmits
the trained learning model parameters to the server. More specifically, at
each communication round $t=1,\dots,T$, each device computes the local
gradient

$\nabla
F_{k}(\bm{\theta}^{(t)})=\frac{1}{|\mathcal{B}_{k}|}\sum_{(\bm{u}_{k,i},\upsilon_{k,i})\in\mathcal{B}_{k}}\nabla
f(\bm{\theta}^{(t)},\bm{u}_{k,i},\upsilon_{k,i}),$ (4)

and sends it to the server. Here, $F_{k}(\bm{\theta}^{(t)})$ represents the
local loss function.

  * •

Step 3: The transmitted updates are aggregated over the air and received by
the server. In particular the received aggregated signal at the parameter
server, i.e., $\widehat{\nabla F}(\bm{\theta}^{(t)})$, yields an estimation of
the global gradient. This allows the server to update the global model using
gradient descent

$\bm{\theta}^{(t+1)}=\bm{\theta}^{(t)}-\eta\widehat{\nabla
F}(\bm{\theta}^{(t)}),$ (5)

where $\eta$ denotes the learning rate. The updated model is then broadcasted
back to the edge devices.

  * •

Step 4: Steps 2 and 3 are repeated until a convergence criterion is met.

Figure 1: Wireless federated learning with $K$ edge devices and an honest-but-
curious parameter server.

###  II-B Communication Model

All edge devices communicate with the edge server over a shared wireless
channel. We assume an ideal downlink channel where the server broadcasts the
initial and the aggregated models (Steps 1 and 3) through a distortion-free
channel. This assumption can be justified by the fact that the edge server can
communicate through a base station, which usually has less stringent
constraints in terms of power and hardware complexity compared to edge
devices. In the uplink transmission (Step 2), in order to facilitate the over-
the-air aggregation of model updates, the local gradients are amplitude-
modulated at the devices and sent to the server using an analog transmission
scheme [20]. The non-ideal hardware in the edge devices distorts the transmit
signal. We model this distortion as power-dependent additive Gaussian noise.
At communication round $t=1,\dots,T$, the transmitted signal from the $k$-th
device can be written as

$\bm{x}_{k}^{(t)}=\sqrt{\frac{\rho_{k}^{(t)}}{\|\bm{g}_{k}^{(t)}\|^{2}}}\bm{g}_{k}^{(t)}+\bm{e}_{k}^{(t)},$
(6)

where $\rho_{k}^{(t)}$ denotes the transmit power of device $k$ at
communication round $t$, and $\bm{g}_{k}^{(t)}=|\mathcal{B}_{k}|\nabla
F_{k}\big{(}\bm{\theta}^{(t)}\big{)}$ is the scaled version of the local
gradient in (4). The additive distortion noise $\bm{e}_{k}^{(t)}$ is
distributed according to
$\mathcal{N}\big{(}\textbf{0},\kappa_{k}\rho_{k}^{(t)}\textbf{I}\big{)}$111Gaussianity
can be justified analytically by the central limit theorem, since
$\bm{e}_{k}^{(t)}$ describes the aggregate effect of different residual
hardware impairments [19]. where $\kappa_{k}\in[0,1)$ is the proportionality
coefficient. The rationale behind this model is that, in many practical cases,
a fixed portion of the signal is turned into distortion [17]. This can happen,
for example, due to amplitude-to-amplitude (AM-AM) nonlinearities in the power
amplifier [18].

We further assume that the edge devices have perfect channel state information
and accordingly, they compensate for the phases of their channels in the
course of transmission. The received signal at the server can therefore be
written as

$\bm{y}^{(t)}=\sum_{k=1}^{K}\big{|}h_{k}^{(t)}\big{|}\bm{x}_{k}^{(t)}+\bm{w}^{(t)},$
(7)

where $h_{k}^{(t)}$ denotes the fading channel coefficient for the $k$-th
device at the $t$-th communication round, and
$\bm{w}^{(t)}\sim\mathcal{N}(0,N_{0})$ is independent and identically
distributed (i.i.d.) additive Gaussian noise. We consider the case of block
fading channels where the channel coefficients remain unchanged within each
time slot, and change independently from one communication round to the other.

###  II-C Threat Model and Privacy Mechanism

We consider an honest-but-curious server that follows the protocol
instructions correctly, but it may attempt to break privacy through observing
the signals received throughout the learning process. Specifically, the server
observes

$\bm{y}^{(t)}=\sum_{k=1}^{K}b_{k}^{(t)}\big{|}h_{k}^{(t)}\big{|}\bm{g}_{k}^{(t)}+\bm{w}_{\text{eff}}^{(t)},$
(8)

over $t=1,\dots,T$, and can infer information about the edge devices’ private
data. Here, $b_{k}^{(t)}$ is given by

$b_{k}^{(t)}=\sqrt{\frac{\rho_{k}^{(t)}}{\|\bm{g}_{k}^{(t)}\|^{2}}},$ (9)

and
$\bm{w}_{\text{eff}}^{(t)}=\bm{w}^{(t)}+\sum_{k=1}^{K}\big{|}h_{k}^{(t)}\big{|}\bm{e}_{k}^{(t)}$
is the effective noise distributed according to
$\mathcal{N}\big{(}0,(\sigma^{(t)})^{2}\big{)}$, with

$\sigma^{(t)}=\sqrt{N_{0}+\sum_{k=1}^{K}\kappa_{k}\rho_{k}^{(t)}|h_{k}^{(t)}|^{2}}.$
(10)

The effective noise $\bm{w}_{\text{eff}}^{(t)}$ will be exploited as a
privacy-preserving mechanism in the next section.

###  II-D Differential Privacy

Differential privacy is a strong and provable privacy guarantee for an
individual’s input to a randomized function, i.e., a privacy mechanism.
Informally, this guarantee implies that the behavior of the mechanism is
essentially unchanged over inputs from any pair of adjacent datasets. In our
case, adjacent datasets are the same, except for an example associated with a
single edge device. Mathematically, this can be expressed as two global
datasets $\mathcal{B}^{\prime}=\cup_{k=1}^{K}\mathcal{B}_{k}^{\prime}$ and
$\mathcal{B}^{\prime\prime}=\cup_{k=1}^{K}\mathcal{B}_{k}^{\prime\prime}$
where $\|\mathcal{B}_{i}^{\prime}-\mathcal{B}_{i}^{\prime\prime}\|_{1}=1$ for
some device $i$ and
$\|\mathcal{B}_{k}^{\prime}-\mathcal{B}_{k}^{\prime\prime}\|_{1}=0$ for all
$k=1,\dots,K$ except for $k=i$ [16]. Using this definition, we can now
formally define the notion of $(\epsilon,\delta)$ differential privacy.

######  Definition 1 (Differential Privacy [21, 16]).

Let $\epsilon>0$ and $0\leq\delta\leq 1$. For any two adjacent datasets
$\mathcal{B}^{\prime}$ and $\mathcal{B}^{\prime\prime}$, $(\epsilon,\delta)$
differential privacy requires that

$\text{P}(\bm{y}|\mathcal{B}^{\prime})\leq\exp{(\epsilon)}\text{P}(\bm{y}|\mathcal{B}^{\prime\prime})+\delta,$
(11)

where $\text{P}(\bm{y}|\mathcal{B}^{\prime})$ and
$\text{P}(\bm{y}|\mathcal{B}^{\prime\prime})$ stand for the distributions of
received signals $\bm{y}=\big{\\{}\bm{y}^{(t)}\big{\\}}_{t=1}^{T}$ conditioned
on the uses of either of the adjacent datasets.

The $(\epsilon,\delta)$ differential privacy condition in (11) can be
rewritten as [21, Lemma 3.17]

$\text{Pr}(|\mathcal{L}_{\mathcal{B}^{\prime},\mathcal{B}^{\prime\prime}}(\bm{y})|\leq\epsilon)\geq
1-\delta,$ (12)

where

$\mathcal{L}_{\mathcal{B}^{\prime},\mathcal{B}^{\prime\prime}}(\bm{y})=\ln\bigg{(}\frac{\text{P}(\bm{y}|\mathcal{B}^{\prime})}{\text{P}(\bm{y}|\mathcal{B}^{\prime\prime})}\bigg{)}$
(13)

is referred to as differential privacy loss. According to (12),
$(\epsilon,\delta)$ differential privacy ensures that the absolute value of
the privacy loss is bounded by $\epsilon$ with probability at least
$1-\delta$. Smaller values of $\epsilon$ and $\delta$ provide therefore
stronger privacy guarantees. More precisely, if (12) is satisfied for
sufficiently small $\epsilon$ and $\delta$, this implies that it is
statistically impossible for the parameter server to detect whether a
particular training sample is included in the training dataset, even in the
extreme case when all other training samples are known.

##  III Privacy Preservation Condition

In this section, we derive a privacy preservation condition by applying the
notion of differential privacy to our problem. To this end, we first derive an
expression for the privacy loss due to the disclosure of the received signal
$\bm{y}=\big{[}\bm{y}^{(1)},\dots,\bm{y}^{(T)}\big{]}$ at the parameter server
by exploiting inherent distortion as the privacy-preserving mechanism. This
can be written as

$\displaystyle\mathcal{L}_{\mathcal{B}^{\prime},\mathcal{B}^{\prime\prime}}(\bm{y})$
$\displaystyle=\ln\Bigg{(}\prod_{t=1}^{T}\frac{\text{P}(\bm{y}^{(t)}|\bm{y}^{(t-1)},\dots,\bm{y}^{(1)},\mathcal{B}^{\prime})}{\text{P}(\bm{y}^{(t)}|\bm{y}^{(t-1)},\dots,\bm{y}^{(1)},\mathcal{B}^{\prime\prime})}\Bigg{)}$
$\displaystyle=\sum_{t=1}^{T}\ln\Bigg{(}\frac{\text{P}(\bm{y}^{(t)}|\bm{y}^{(t-1)},\dots,\bm{y}^{(1)},\mathcal{B}^{\prime})}{\text{P}(\bm{y}^{(t)}|\bm{y}^{(t-1)},\dots,\bm{y}^{(1)},\mathcal{B}^{\prime\prime})}\Bigg{)}.$
(14)

From (8) and for a given dataset $\mathcal{B}=\cup_{k=1}^{K}\mathcal{B}_{k}$,
we have

P
$\displaystyle\big{(}\bm{y}^{(t)}|\bm{y}^{(t-1)},\dots,\bm{y}^{(1)},\mathcal{B}\big{)}=\frac{1}{\sigma^{(t)}\sqrt{2\pi}}$
$\displaystyle\cdot\exp\left(-\frac{\Big{\|}\bm{y}^{(t)}-\sum_{k=1}^{K}b_{k}^{(t)}(\mathcal{B}_{k})\big{|}h_{k}^{(t)}\big{|}\bm{g}_{k}^{(t)}(\mathcal{B}_{k})\Big{\|}^{2}}{2(\sigma^{(t)})^{2}}\right),$
(15)

where $b_{k}^{(t)}(\mathcal{B}_{k})$ can be calculated by replacing
$\bm{g}_{k}^{(t)}(\mathcal{B}_{k})$ in (9). Now, by using (III) in (III) we
obtain

$\displaystyle\mathcal{L}_{\mathcal{B}^{\prime},\mathcal{B}^{\prime\prime}}(\bm{y})$
$\displaystyle~{}~{}=\sum_{t=1}^{T}\ln\left(\frac{\exp\left(-\frac{\Big{\|}\bm{y}^{(t)}-\sum_{k=1}^{K}b_{k}^{(t)}(\mathcal{B}_{k}^{\prime})\left|h_{k}^{(t)}\right|\bm{g}_{k}^{(t)}(\mathcal{B}_{k}^{\prime})\Big{\|}^{2}}{2(\sigma^{(t)})^{2}}\right)}{\exp\left(-\frac{\Big{\|}\bm{y}^{(t)}-\sum_{k=1}^{K}b_{k}^{(t)}(\mathcal{B}_{k}^{\prime\prime})\left|h_{k}^{(t)}\right|\bm{g}_{k}^{(t)}(\mathcal{B}_{k}^{\prime\prime})\Big{\|}^{2}}{2(\sigma^{(t)})^{2}}\right)}\right)$
$\displaystyle~{}~{}=\sum_{t=1}^{T}\ln\left(\frac{\exp\left(-\frac{\|\bm{w}_{\text{eff}}^{(t)}\|^{2}}{2(\sigma^{(t)})^{2}}\right)}{\exp\left(-\frac{\|\bm{w}_{\text{eff}}^{(t)}+\bm{v}^{(t)}\|^{2}}{2(\sigma^{(t)})^{2}}\right)}\right)$
$\displaystyle~{}~{}=\sum_{t=1}^{T}\frac{\|\bm{v}^{(t)}\|^{2}+2(\bm{w}_{\text{eff}}^{(t)})^{\mathsf{T}}\bm{v}^{(t)}}{2(\sigma^{(t)})^{2}}\overset{\Delta}{=}\Gamma,$
(16)

where

$\bm{v}^{(t)}=\sum_{k=1}^{K}\big{|}h_{k}^{(t)}\big{|}\Big{(}b_{k}^{(t)}(\mathcal{B}_{k}^{\prime\prime})\bm{g}_{k}^{(t)}(\mathcal{B}_{k}^{\prime\prime})-b_{k}^{(t)}(\mathcal{B}_{k}^{\prime})\bm{g}_{k}^{(t)}(\mathcal{B}_{k}^{\prime})\Big{)}.$
(17)

By applying the triangle inequality, the norm of $\bm{v}^{(t)}$ can be bounded
as

$\|\bm{v}^{(t)}\|\leq
2\max_{k}\sqrt{\rho_{k}^{(t)}}\big{|}h_{k}^{(t)}\big{|}\overset{\Delta}{=}\Delta^{(t)}.$
(18)

Following similar steps as in [21, Appendix A] and [16, Eq. (52)], the
following upper-bound is obtained on the privacy violation probability

$\displaystyle\text{Pr}\left(\left|\Gamma\right|>\epsilon\right)$
$\displaystyle\overset{(\text{i})}{\leq}\text{Pr}\left(\left|\sum_{t=1}^{T}\frac{(\bm{w}_{\text{eff}}^{(t)})^{\mathsf{T}}\bm{v}^{(t)}}{(\sigma^{(t)})^{2}}\right|>\epsilon-\sum_{t=1}^{T}\frac{\|\bm{v}^{(t)}\|^{2}}{2(\sigma^{(t)})^{2}}\right)$
$\displaystyle=2\text{Pr}\left(\sum_{t=1}^{T}\frac{(\bm{w}_{\text{eff}}^{(t)})^{\mathsf{T}}\bm{v}^{(t)}}{(\sigma^{(t)})^{2}}>\epsilon-\sum_{t=1}^{T}\frac{\|\bm{v}^{(t)}\|^{2}}{2(\sigma^{(t)})^{2}}\right)$
$\displaystyle\overset{(\text{ii})}{\leq}2\text{Pr}\left(\Lambda>\epsilon-\sum_{t=1}^{T}\frac{1}{2}\left(\frac{\Delta^{(t)}}{\sigma^{(t)}}\right)^{2}\right),$
(19)

where (i) comes from the inequality
$\text{Pr}(|X+c|>\epsilon)\leq\text{Pr}(|X|+c>\epsilon)$ for any arbitrary
$c\geq 0$, and (ii) is obtained using (18) where

$\Lambda\sim\mathcal{CN}\left(0,\sum_{t=1}^{T}\left(\frac{\Delta^{(t)}}{\sigma^{(t)}}\right)^{2}\right).$
(20)

The upper-bound in (III) can therefore be written as

$\displaystyle\text{Pr}\left(\left|\Gamma\right|>\epsilon\right)$
$\displaystyle\leq\frac{2}{\sqrt{2\pi\sum_{t=1}^{T}\left(\frac{\Delta^{(t)}}{\sigma^{(t)}}\right)^{2}}}\int_{s}^{\infty}\exp\left(-\frac{z^{2}}{2\sum_{t=1}^{T}\left(\frac{\Delta^{(t)}}{\sigma^{(t)}}\right)^{2}}\right)\text{d}z,$
(21)

where

$s=\epsilon-\sum_{t=1}^{T}\frac{1}{2}\left(\frac{\Delta^{(t)}}{\sigma^{(t)}}\right)^{2}>0.$
(22)

Finally, by using (III)-(22) in (12), the privacy preservation condition is
expressed as

$\frac{2}{\sqrt{2\pi\nu}}\int_{\epsilon-\frac{\nu}{2}}^{\infty}\exp\left(-\frac{z^{2}}{2\nu}\right)\text{d}z<\delta,$
(23)

where

$\nu=\sum_{t=1}^{T}\left(\frac{\Delta^{(t)}}{\sigma^{(t)}}\right)^{2}=\sum_{t=1}^{T}\frac{\left(2\max_{k}\sqrt{\rho_{k}^{(t)}}\big{|}h_{k}^{(t)}\big{|}\right)^{2}}{N_{0}+\sum_{k=1}^{K}\kappa_{k}\rho_{k}^{(t)}|h_{k}^{(t)}|^{2}}.$
(24)

##  IV Privacy-Preserving Power Allocation Scheme

In this section, we present an offline-optimized privacy-preserving power
allocation scheme. In particular, we assume that the parameters
$\\{h_{k}^{(t)},\kappa_{k},N_{0}\\}$ are known and we seek the values
$\rho_{k}^{(t)}$ for $k=1,\dots,K$ and $t=1,\dots,T$ to minimize the
distortion of the recovered gradient from the received signal in (8) with
respect to the ground-truth global gradient, i.e.,

$\displaystyle\text{MSE}^{(t)}=\mathbb{E}\left[\left\|\widehat{\nabla
F}(\bm{\theta}^{(t)})-{\nabla F}(\bm{\theta}^{(t)})\right\|^{2}\right]$
$\displaystyle=\frac{1}{K^{2}}\mathbb{E}\left[\left\|\sum_{k=1}^{K}\frac{\sqrt{\rho_{k}^{(t)}}\big{|}h_{k}^{(t)}\big{|}\bm{g}_{k}^{(t)}}{\xi^{(t)}\big{\|}\bm{g}_{k}^{(t)}\big{\|}}+\frac{\bm{w}_{\text{eff}}^{(t)}}{\xi^{(t)}}-\sum_{k=1}^{K}\frac{\bm{g}_{k}^{(t)}}{\big{\|}\bm{g}_{k}^{(t)}\big{\|}}\right\|^{2}\right],$
(25)

while satisfying the differential privacy condition in (23). In (IV), the
expectation is over the distribution of the effective noise
$\bm{w}_{\text{eff}}$ and $\xi^{(t)}$ is the normalization factor applied by
the parameter server for recovering the gradient of interest. We assume that
the edge devices are subject to a peak power constraint, i.e.,

$(1+\kappa_{k})\rho_{k}^{(t)}\leq\rho_{\text{max}}.$ (26)

Accordingly, the relevant optimization problem can be defined as

$\displaystyle\underset{\rho_{k}^{(t)}}{\text{minimize}}~{}~{}$
$\displaystyle\text{MSE}^{(t)}$ (27) subject to
$\displaystyle\eqref{eq:DPcondition}~{}\text{and}~{}\eqref{eq:MaxPower}.$

To guarantee a perfectly unbiased global model, the power scaling factors
should be set such that the following gradient alignment condition is
satisfied,

$\sqrt{\rho_{k}^{(t)}}\big{|}h_{k}^{(t)}\big{|}=\lambda^{(t)}.$ (28)

In this case, by setting the normalization factor $\xi^{(t)}=\lambda^{(t)}$,
the mean squared error in (IV) simplifies to

$\text{MSE}^{(t)}=\mathbb{E}\left[\left\|\frac{\bm{w}_{\text{eff}}^{(t)}}{\lambda^{(t)}}\right\|^{2}\right]=\frac{N_{0}+\sum_{k=1}^{K}\kappa_{k}\lambda^{(t)}}{K^{2}(\lambda^{(t)})^{2}}.$
(29)

When no privacy constraint is in place, the optimal solution for minimizing
(29) under the power constraint in (26) is given by [22]

$\stackon[-1ex]{\rho}{\mathchar
12414\relax}_{k}^{(t)}=\frac{\rho_{\text{max}}\big{|}h_{j}^{(t)}\big{|}^{2}}{(1+\kappa_{j})\big{|}h_{k}^{(t)}\big{|}^{2}},$
(30)

where $j$ is the index of the weakest channel, i.e.,

$j=\operatorname*{arg\,min}_{k}\big{|}h_{k}^{(t)}\big{|}.$ (31)

In words, the user with the weakest channel performs a full power transmission
while the others apply power control with channel inversion to satisfy (28)
using the resulting $\stackon[-1ex]{\lambda}{\mathchar
12414\relax}^{(t)}=\sqrt{\rho_{\text{max}}\big{|}h_{j}^{(t)}\big{|}^{2}/(1+\kappa_{j})}$.
Note that when differential privacy is enforced, due to the constraint (23),
$\stackon[-1ex]{\lambda}{\mathchar 12414\relax}^{(t)}$ may no longer be a
feasible solution. In this case, noting that the privacy violation probability
is an increasing function of $\lambda^{(t)}$, the value of
$\stackon[-1ex]{\lambda}{\mathchar 12414\relax}^{(t)}$ should be decreased
until the constraint (23) is satisfied. Denote the maximum values of
$\lambda^{(t)}$ satisfying (23) by $\lambda_{\text{p}}^{(t)}$. The solution to
the optimization problem in (27) can then be stated as

$\breve{\rho}_{k}^{(t)}=\frac{\big{|}\breve{\lambda}^{(t)}\big{|}^{2}}{\big{|}h_{k}^{(t)}\big{|}^{2}},$
(32)

where

$\breve{\lambda}^{(t)}=\min\left\\{\stackon[-1ex]{\lambda}{\mathchar
12414\relax}^{(t)},\lambda_{\text{p}}^{(t)}\right\\},$ (33)

where $\lambda_{\text{p}}^{(t)}$ can be obtained using a simple line search.

##  V Numerical Results

In this section, we evaluate the performance of the proposed privacy-
preserving power allocation schemes under different hardware impairment
levels.222For the sake of reproducibility, the codes used to obtain the
simulation results will be made available after the review process is
completed. We consider the learning task of handwritten digit recognition
using the well-known MNIST dataset, which consists of $60000$ training samples
and $10000$ test samples. Each example is a $28\times 28$ pixel gray scale
image. The training samples are evenly distributed among the $K$ devices. As a
baseline classification model, we consider a feedforward neural network with a
single hidden layer containing $100$ hidden units followed by a softmax output
layer. We use Adam as the optimizer with an initial learning rate of
$\eta=0.001$; the optimizer automatically adjusts the learning rate as the
number of learning iterations increases to achieve faster convergence. We
assume that the number of local steps an edge device performs before
transmitting an update is equal to $30$. The batch size is set to $128$.

Figure 2: Test accuracy for different power allocation strategies under
several hardware impairment levels.

We consider a scenario with $K=50$ edge devices and set the number of
communication rounds to $T=10$. The channel coefficients $h_{k}^{(t)}$ are
drawn from a complex circularly-symmetric Gaussian distribution. The peak
power constraint is set to $\rho_{\text{max}}=10$ dBm for all devices. The
receiver noise variance is $N_{0}=-20$ dBm and the differential privacy
parameters are set to $\epsilon=25$ and $\delta=0.05$. Finally, we assume that
the proportionality coefficient is equal for all edge devices, i.e.,
$\kappa_{k}=\kappa$ for $k=1,\dots,K$. We consider three different hardware
impairment levels, $\kappa\in\\{0,0.01,0.1\\}$. Note that $\kappa=0$
corresponds to a scenario with ideal hardware, whereas $\kappa=0.01$ and
$\kappa=0.1$ represent cases with low-to-moderate and high levels of hardware
impairment, respectively.333The proportionality coefficient $\kappa$ is
related to the transmit error vector magnitude (EVM) as
$\text{EVM}=\sqrt{\kappa}$ [17]. Accordingly, $\kappa=0.01$ and $\kappa=0.1$
correspond to $10\%$ and $31.62\%$ EVM values, respectively. Typical EVM
values can be derived from the minimum requirements specified in established
standards. In each of these scenarios, we evaluate the learning performance
versus the communication round. As benchmarks, we also include the results for
the cases with distortion-unaware power allocation (where the effect of
hardware impairments is ignored and the power allocation is carried out
considering $\kappa=0$ regardless of its actual value) and the case with ideal
hardware and no privacy constraint.

Fig. 2 depicts the test accuracy, averaged over 50 independent trials with
random splitting of the dataset in each trial, for the cases mentioned above.
From the figure, it can be seen that for ideal hardware, enforcing the privacy
constraint leads to a significant degradation in learning performance. This is
because in this case, receiver noise is the only privacy-preserving mechanism,
and in order to satisfy the differential privacy constraint, devices must
scale down their transmit power. However, by adopting the proposed power
allocation schemes in realistic scenarios with imperfect hardware, the
hardware-induced distortion at the devices can be used together with the
receiver noise to satisfy the privacy requirements with a lower penalty in
terms of learning performance. The proposed power allocation strategy can,
therefore, achieve considerable performance gains compared to the conventional
strategies that ignore the effect of hardware-induced distortion. In certain
cases, the imperfect hardware devices can even achieve better performance than
the ideal hardware scenario. An example of such gain attained by utilizing
hardware-induced distortion can be seen for the case with $\kappa=0.01$ (red
curves with square markers). In particular, with $\kappa=0.01$, the devices
can transmit more power and yet guarantee the privacy preserving condition in
(23), resulting in improved learning performance despite the slightly
increased effective noise variance.

It should be noted that whether the imperfect hardware yields an improved
performance with respect to the ideal hardware depends on the level of
impairments as well as the stringency of the privacy requirements. For
instance, in the example given in Fig. 2, it can be seen that the case with
$\kappa=0.1$ (yellow curves with circle markers) yields a degraded performance
with respect to the ideal hardware. The reason behind this is that the amount
of distortion injected into the transmitted updates is more than enough for
satisfying the differential privacy requirements. As a result, the amount of
transmit power at the devices is determined by the peak power constraint in
(26) and the excessive distortion degrades the learning performance. Adopting
a truncation-based approach by excluding the edge devices that experience deep
fading channels similar to the solution in [4], can reduce the amount of
excessive distortion. A detailed investigation of this approach is left for
future research.

##  VI Conclusion

We studied differentially-private wireless federated learning in realistic
scenarios where edge devices are equipped with imperfect hardware. By modeling
the hardware-induced distortion as power-dependent additive white Gaussian
noise, we derived an expression for the privacy preservation condition and
proposed a power allocation scheme to satisfy it. Our numerical results
indicate that exploiting the inherent distortion considerably reduces the
performance loss caused by the enforcement of differential privacy
constraints. As a result, our proposed power allocation scheme achieves
significant gains in terms of learning performance compared to conventional
approaches that ignore the effect of hardware-induced distortion.

## References

  * [1] W. Y. B. Lim, N. C. Luong, D. T. Hoang, Y. Jiao, Y.-C. Liang, Q. Yang, D. Niyato, and C. Miao, “Federated learning in mobile edge networks: A comprehensive survey,” _IEEE Commun. Surv. Tut._ , vol. 22, no. 3, pp. 2031–2063, third quarter 2020.
  * [2] P. Kairouz, H. B. McMahan, B. Avent, A. Bellet, M. Bennis, A. N. Bhagoji, K. Bonawitz, Z. Charles, G. Cormode, R. Cummings _et al._ , “Advances and open problems in federated learning,” _arXiv preprint arXiv:1912.04977_ , 2019.
  * [3] K. Yang, T. Jiang, Y. Shi, and Z. Ding, “Federated learning via over-the-air computation,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 3, pp. 2022–2035, Jan. 2020.
  * [4] G. Zhu, Y. Wang, and K. Huang, “Broadband analog aggregation for low-latency federated edge learning,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 1, pp. 491–506, Jan. 2020.
  * [5] M. M. Amiri and D. Gündüz, “Federated learning over wireless fading channels,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 5, pp. 3546–3557, May 2020.
  * [6] M. M. Amiri, T. M. Duman, D. Gunduz, S. R. Kulkarni, and H. V. Poor, “Blind federated edge learning,” _arXiv preprint arXiv:2010.10030_ , 2020.
  * [7] B. Tegin and T. M. Duman, “Machine learning at wireless edge with ofdm and low resolution ADC and DAC,” _arXiv preprint arXiv:2010.00350_ , 2020.
  * [8] M. Fredrikson, S. Jha, and T. Ristenpart, “Model inversion attacks that exploit confidence information and basic countermeasures,” in _Proc. ACM SIGSAC Conf. Comput. Commun. Security_ , 2015, pp. 1322–1333.
  * [9] L. Melis, C. Song, E. De Cristofaro, and V. Shmatikov, “Exploiting unintended feature leakage in collaborative learning,” in _Proc. IEEE Symp. Secur. Privacy (SP)_ , May 2019, pp. 691–706.
  * [10] H. Yin, A. Mallya, A. Vahdat, J. M. Alvarez, J. Kautz, and P. Molchanov, “See through gradients: Image batch recovery via gradinversion,” in _Proc. IEEE/CVF Conf. Computer Vision Pattern Recognition_ , 2021, pp. 16 337–16 346.
  * [11] J. Geiping, H. Bauermeister, H. Dröge, and M. Moeller, “Inverting gradients–how easy is it to break privacy in federated learning?” _arXiv preprint arXiv:2003.14053_ , 2020.
  * [12] F. Mo, A. Borovykh, M. Malekzadeh, H. Haddadi, and S. Demetriou, “Quantifying information leakage from gradients,” _arXiv preprint arXiv:2105.13929_ , 2021\.
  * [13] M. Seif, R. Tandon, and M. Li, “Wireless federated learning with local differential privacy,” in _Proc. IEEE Int. Symp. Inf. Theory (ISIT)_ , Los Angeles, CA, USA, Jun. 2020, pp. 2604–2609.
  * [14] K. Wei, J. Li, M. Ding, C. Ma, H. H. Yang, F. Farokhi, S. Jin, T. Q. Quek, and H. V. Poor, “Federated learning with differential privacy: Algorithms and performance analysis,” _IEEE Trans. Inf. Forensics Security_ , vol. 15, pp. 3454–3469, Apr. 2020.
  * [15] Y. Koda, K. Yamamoto, T. Nishio, and M. Morikura, “Differentially private AirComp federated learning with power adaptation harnessing receiver noise,” in _Proc. IEEE Global Commun. Conf. (GLOBECOM)_ , Taipei, Taiwan, Dec. 2020.
  * [16] D. Liu and O. Simeone, “Privacy for free: Wireless federated learning via uncoded transmission with adaptive power control,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 1, pp. 170–185, Jan. 2021.
  * [17] E. Björnson, J. Hoydis, M. Kountouris, and M. Debbah, “Massive MIMO systems with non-ideal hardware: Energy efficiency, estimation, and capacity limits,” _IEEE Trans. Inf. Theory_ , vol. 11, no. 60, pp. 7112–7139, Nov. 2014.
  * [18] T. Schenk, _RF imperfections in high-rate wireless systems: Impact and digital compensation_. Springer Science & Business Media, 2008.
  * [19] C. Studer, M. Wenk, and A. Burg, “MIMO transmission with residual transmit-RF impairments,” in _Proc. Int. ITG Workshop on Smart Antennas (WSA)_ , Bremen, Germany, Feb. 2010, pp. 189–196.
  * [20] M. M. Amiri and D. Gündüz, “Machine learning at the wireless edge: Distributed stochastic gradient descent over-the-air,” _IEEE Trans. Signal Process._ , vol. 68, pp. 2155–2169, Mar. 2020.
  * [21] C. Dwork and A. Roth, “The algorithmic foundations of differential privacy,” _Foundations and Trends in Theoretical Computer Science_ , vol. 9, no. 3-4, pp. 211–407, 2014.
  * [22] X. Cao, G. Zhu, J. Xu, and K. Huang, “Optimized power control for over-the-air computation in fading channels,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 11, pp. 7498–7513, Aug. 2020.


∎

e1e-mail<EMAIL_ADDRESS>e2e-mail<EMAIL_ADDRESS>
11institutetext: National Centre for Nuclear Research, Department of
Fundamental Research, 02-093 Warsaw, Poland 22institutetext: Institute of
Physics, Jan Kochanowski University, 25-406 Kielce, Poland

# Beyond the relaxation time approximation

Grzegorz Wilke1,addr1  Zbigniew Włodarczyke2,addr2

(Received: date / Accepted: date)

###### Abstract

The relaxation time approximation (RTA) is a well known method of describing
the time evolution of a statistical ensemble by linking distributions of the
variables of interest at different stages of their temporal evolution. We show
that if all the distributions occurring in the RTA have the same functional
form of a quasi-power Tsallis distribution the time evolution of which depends
on the time evolution of its control parameter, nonextensivity $q(t)$, then it
is more convenient to consider only the time evolution of this control
parameter.

††journal: Eur. Phys. J. A

##  1 Introduction

Many problems in science involve understanding the time evolution of a
statistical ensemble. Here we will focus on a system of particles described by
the probability distribution $f(r,p,t)$ which depends on position $r$,
momentum $p$ and time $t$. In general, for an evolving physical system
operating irreversibly out of thermodynamical equilibrium $f(r,p,t)$ will
differ from that of a Boltzmannian ensemble and its evolution is usually
studied using the Boltzmann transport equation (BTE), the general form of
which is Huang ,

$\frac{df(r,p,t)}{dt}=\frac{\partial f}{\partial
t}+\vec{u}\cdot\nabla_{r}f+\vec{F}\cdot\nabla_{p}f=C[f],$ (1)

where $\vec{F}$ is the external force, $\vec{u}$ the velocity and $C[f]$ the
collision term. Narrowing our interest to situations where the considered
system is homogeneous (i.e. $\nabla_{r}f=0$) and when no external forces are
acting (i.e. when $\vec{F}=0$) Eq. (1) simplifies to the form

$\frac{df(r,p,t)}{dt}=\frac{\partial f}{\partial t}=C[f(t)].$ (2)

Because of the freedom in choosing the functional form of $C[f(t)]$ it is
still a very general equation that allows one to deal with a variety of
situations. However, usually in many applications it is further simplified to
a form called the relaxation time approximation (RTA) which consists in using
such a simple form of the collision term Huang ; BGK ; AW ; FR :

$C[f]=\frac{f_{eq}-f}{\tau},$ (3)

where $f_{eq}$ is the local equilibrium distribution and $\tau$ is the
relaxation time, understood as the time taken by the non-equilibrium system to
reach equilibrium. In this approximation the BTE simplifies further to

$\frac{\partial f}{\partial t}=\frac{f_{eq}-f}{\tau}.$ (4)

Solving this equation for the initial conditions such that at $t=0$ one has as
initial (assumed) distribution, $f=f_{in}$, and at freeze-out time, $t=t_{f}$
one has a final distribution, $f=f_{fo}$ (which we identify with the
distribution we are looking for that actually describes the distribution
obtained experimentally)111The statistical system produced in multiparticle
production processes quickly reaches an initial distribution (pre-equilibrium
state) which slowly evolves to equilibrium but becomes frozen at the freeze-
out time (usually the system at freeze-out is not in thermodynamic
equilibrium). The experimentally measured spectra of the produced particles
reflect the state of the system at freeze-out. one finds that

$f_{fo}=f_{eq}+\left(f_{in}-f_{eq}\right)\exp\left(-\frac{t_{f}}{\tau}\right).$
(5)

The continued popularity of such an approach to the analysis of various
particle production processes can be proved by the fact that recently the
Boltzmann transport equation in the RTA approximation was used to analyze
various observables in nucleus-nucleus collisions (in particular, to study the
time evolution of temperature fluctuations in a non-equilibrium system BGSS ,
to describe the elliptic flow YTTS , transverse momentum spectra ZLD , as well
as for the study of nuclear modification factors at RHIC and LHC energies TKTS
; TBGKSC ; QCGZZ ).

##  2 The two-component nature of the RTA approximation

In what follows we consider (for simplicity) massless objects in a one-
dimensional statistical ensemble and assume that all probability distributions
have the form $f(z,t)$, where $z$ is a scaled variable, $z=x/x_{0}$, for
observable $x$ with $x_{0}$ being the corresponding scale parameter. The
evolution of our system reaches asymptotically, for $t_{f}\to\infty$, a state
of local equilibrium in which $f_{fo}=f_{eq}$. In such a case

$f_{eq}\left(z\right)=exp\left(-z\right),$ (6)

which for $z=E/T$ (with $E$ being energy and the scale parameter $T$ being
temperature) considered here is simply the Boltzmann-Gibbs (BG)
distribution222The choice of the exponential distribution (6) as the final
equilibrium distribution results from the following premises: it describes the
spectra in the thermalized system, it is a memoryless probability
distribution, it describes the processes in which events occur independently
and it is the most probable distribution for the fixed expected value..

Now note that the final distribution given by Eq. (5) can be rewritten in a
two-component form,

$f_{fo}(z)=f_{in}(z)\exp\left(-\frac{t_{f}}{\tau}\right)+f_{eq}(z)\left[1-\exp\left(-\frac{t_{f}}{\tau}\right)\right].$
(7)

This emphasizes that the sought two-component final distribution $f_{fo}$ is
bistable and jumps between the BG distribution $f_{eq}$ and some assumed
initial distribution $f_{in}$ with probability described by the parameter
$t_{f}/\tau$. To illustrate this effect we use quasi-power like initial
distribution $f_{in}(z)$

$f_{in}\left(z\right)=C\cdot\left[1+\frac{z}{\gamma}\right]^{-\gamma},$ (8)

with $\gamma=4$ and $C=(\gamma-1)/\gamma$ being the normalization factor, and
the exponential distribution (6) for $f_{eq}$; for this choice we get $f_{fo}$
for different values of $t_{f}/\tau$ as presented in Fig. 1. The clearly
visible structure of $f_{fo}(z)$ for some choice of parameters clearly
illustrates the potential bi-stability of the final system being a consequence
of its composition of quasi-power and exponential distributions.

Figure 1: Schematic probability distributions $f_{fo}(z)$ resulting from the
relaxation time approximation scenario calculated using Eq.(7) for $f_{in}$
given by Eq.(8) with $\gamma=4$ and $f_{eq}$ given by Eq.(6), for $t/\tau=0$,
$5$, $10$, $\infty$ (curves from top down).

Let us summarize this part. In the RTA we always have at our disposal two
distributions: the initial $f_{in}$ and the exponential $f_{eq}$. In the case
considered here at the beginning, i.e. for $t=0$, all the objects of the
system are described by the power distribution. Over time, some objects
(decreasing with time) are still described by the power distribution, but the
remainder (increasing with time) become thermalized and are described by the
exponential distribution. So: for a given $t$, every object of the system is
either in the initial state or in the final state and the occupation numbers
of these states changes with time.

##  3 Beyond the RTA - time evolution of internal variable

Note however that a given non-equilibrium process, regardless of the specific
dynamics, evolves the probability distribution over the system configurations.
The system traverses a probability manifold, which is not a manifold of
equilibrium states because the distribution at each point on the manifold need
not correspond to a BG distribution NatPhys . This means that for every $t$
one has some distribution $f(z,t)$ and all objects in the system are described
by that (single) distribution. In our case this distribution should smoothly
evolve with time from power to exponential (equilibrium) form, i.e., it should
be a quasi-power like distribution. The most common distribution of this type
is the Tsallis distribution T ; T1 ; T3 ,

$f(z,t)=(2-q)\left[1+(q-1)z\right]^{\frac{1}{1-q}},$ (9)

characterized by the energy dependent parameter $q=q(t)$ (note that for $q\to
1$ the Tsallis distribution becomes the BG from Eq. (6)). To justify such a
choice let us recall a unique feature of distribution (9) which distinguishes
it from all other distributions used so far in this context. Namely, whereas
for $z\to\infty$ (or for $z>>1/(q-1)$) it becomes a power distribution, for
$z\to 0$ (or for $z<<1/(q-1)$) it goes into an exponential distribution. Thus,
for an appropriate selection of the parameter $q(t)$, it can describe all the
distributions occurring in the RTA formula: $f_{in}(z,t)$ for
$q\left(t_{f}=0\right)=q_{in}>1$ and $f_{eq}$ for
$q\left(t_{f}\to\infty\right)=1$.

The effectiveness of the Tsallis distribution is best evidenced by the results
of works WW1 ; WWCT ; WW2 . In particular, as shown in WWCT , it nicely
describes a wide range of the measured transverse momenta ($0.1<p_{T}<100$ GeV
which corresponds in this case to $1<z<700$) in which the cross section spans
a range of $\sim 14$ orders of magnitude333Recently, inspired by the Tsallis
statistics for a non-equilibrium distributions, the non-extensive BTE in the
relaxation time approximation was discussed in qBTE ..

The time-dependent parameter $q=q(t)$ (more specifically, its deviation from
$q=1$) represents the degree of the non-extensivity or, in other words, the
degree of deviation of the system from the thermalized or equilibrated system,
which is usually described by the well known BG statistical mechanics. It is
also a control parameter that fully defines the shape of the Tsallis
distribution, in particular its evolution over time through moments such as
the expected value, $\langle z(t)\rangle$ and variance, $Var[z(t)]$:

$\displaystyle\langle z(t)\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{3-2q(t)},$ (10) $\displaystyle Var(z)$ $\displaystyle=$
$\displaystyle\frac{2-q(t)}{(3-2q(t))^{2}(4-3q(t))}.$ (11)

In general, the moments $<z^{n}>$, limited to $n+1<1/(q-1)$, are related by
the recurrence relation

$\langle z^{n}\rangle=\frac{n}{1-(n+1)(q-1)}\langle z^{n-1}\rangle.$ (12)

As mentioned above, we assume that the dynamic evolution of a system over time
smoothly and monotonically transforms the probability distribution $f(z,t)$ in
$(z,t)$-space and that $f(z,t)$ is a Tsallis distribution fully described by
the time-dependent control parameter $q(t)$.

Now note that if we had completely formally used the Tsallis distributions for
all distributions in Eq. (5) defining the RTA and comparing $f(z=0)$ or $<z>$,
we would get such a relationship between the parameters $q$ appearing there
(remembering that $f_{eq}$ is assumed as a BG distribution which is equavalent
to a Tsallis distribution with $q_{eq}=1$):

$\left(2-q_{fo}\right)=1+\left[\left(2-q_{in}\right)-1\right]\cdot\exp\left(-\frac{t_{f}}{\tau}\right).$
(13)

However, if from the very beginning we decide to describe the entire process
using only quasi-power Tsallis distributions, the time evolution of which is
given only by the time evolution of their control parameters $q=q(t)$, which
means that $f(t)=f[q(t)]$, we should go back to Eq. (2), which is now

$\frac{\partial f(t)}{\partial t}=F[q(t)].$ (14)

This equation replaces Eq. (2). The form of the function $F$ from Eq. (14) can
be deduced by taking $f(t)$ given by the Tsallis distribution with $q=q(t)$
and calculating $df/dt$. As a result, we get that

$\displaystyle F\left[z,q(t)\right]$ $\displaystyle=$ $\displaystyle
f\left[z,q(t)\right]\cdot Q\cdot\frac{1}{[q(t)-1]^{2}}\frac{dq(t)}{dt},$
$\displaystyle Q$ $\displaystyle=$
$\displaystyle\ln\left[1+(q(t)-1)z\right]+\frac{1}{1+[q(t)-1]z}-$ (15)
$\displaystyle-\frac{[q(t)-1]^{2}}{2-q(t)}-1,$

and expressing the dependence of $Q(z)$ on the variable $z$ by the function
$f=f[z,q(t),t]$ we have

$\displaystyle\frac{\partial f[z,q(t),t]}{\partial t}=-\frac{f}{\tau}\cdot$
$\displaystyle~{}~{}~{}~{}~{}\cdot\left\\{\frac{1}{q-1}\left[\left(\frac{f}{2-q}\right)^{q-1}-\frac{1}{2-q}\right]-\ln\left(\frac{f}{2-q}\right)\right\\}.$
(16)

To go further, we need to set the time dependence of the parameter $q$ in some
way. Note that the nonextensivity parameter $q(t)$ describes deviations of the
state of a statistical system from equilibrium and in this sense it plays the
role of an internal variable discussed in Refs. IP ; VP . Therefore, following
such an approach VP  we assume that the equation of the dynamics describing
the control parameter $q(t)$ is of the form of the equation of a relaxation:

$\frac{dq}{dt}=\frac{q_{eq}-q}{\tau}.$ (17)

Remembering that we always assume that $q_{eq}=1$, the solution of (17) is

$q-1=\left(q_{in}-1\right)\exp\left(-\frac{t_{f}}{\tau}\right)$ (18)

which coincides with Eq. (13). Fig. 2 shows the resultant schematic
distributions $f_{fo}$ for different $t_{f}/\tau$; they all have the form of a
Tsallis distribution with $q=q\left(t=t_{f}\right)$ as given by Eq. (18). As
one can see the result is now different from that using the RTA approximation
shown in Fig. 1.

Figure 2: Schematic probability distributions $f_{fo}(z)$ calculated from Eq.
(14) for $t/\tau=0$, $0.5$, $1$, $2$, $\infty$ (curves from top down). The
distributions $f_{in}$ and $f_{eq}$ are the same as in Fig. 1 (given by the
Tsallis formula (9) with $q_{in}=1.25$ and $q_{eq}=1$, respectively).

Notice that the relaxation times $\tau=\tau_{f}$ in Eq.(3) and $\tau=\tau_{q}$
in Eq.(17) describe the relaxation of different quantities, respectively the
entire distribution and its control parameter. Comparing (for the same time)
mean values $\langle
z\rangle=1+\exp\left(-\frac{t}{\tau_{f}}\right)\left(\langle
z\rangle_{in}-1\right)$ evaluated from the RTA, Eq. (7), and its $q$ version,
$\langle z\rangle=\left[1+\exp\left(-\frac{t}{\tau_{q}}\right)\left(\langle
z\rangle_{in}^{-1}-1\right)\right]^{-1}$, Eq. (9) for $q(t)$ given by Eq.
(18), we have that

$\frac{t}{\tau_{f}}=\frac{t}{\tau_{q}}+\ln[<z>_{in}+(1-<z>_{in})\exp(-t/\tau_{q})]$
(19)

and the ratio of relaxation times in both approaches is

$\frac{\tau_{f}}{\tau_{q}}=\frac{t/\tau_{q}}{t/\tau_{q}+\ln[<z>_{in}+(1-<z>_{in})\exp(-t/\tau_{q})]}$
(20)

changing from $\tau_{f}/\tau_{q}=1/<z>_{in}$ for $t\rightarrow 0$ to
$\tau_{f}/\tau_{q}=1$ for $t\rightarrow\infty$ (cf. Fig 3).

Figure 3: Time dependence of the the relaxation times ratio
$\tau_{f}/\tau_{q}$ for an initial mean value $<z>_{in}=2$ (corresponding to
$q_{in}=1.25$).

Note that in a situation where in some isolated system we have a fixed number
of particles $N$ and a fixed total energy $U$, we have a constant average
energy $\langle E\rangle$. Therefore, in such a case, the variability in time
of $\langle z\rangle=\langle E\rangle/T$ must be caused by the appropriate
variability in time of the scale parameter $T$ (here the temperature). This
means that the scale parameter in our scaled variable $z=x/x_{0}$ also changes
with time:

$x_{0}(t)=x_{0}^{*}[1-2(q_{in}-1)\exp(-t/\tau)],$ (21)

where $x_{0}=x_{0}^{*}$ at equilibrium ($t\longrightarrow\infty$). Because, as
known from WW , fluctuations of the scale parameter $x_{0}$ are directly
connected with the parameter $q$:

$\frac{Var(1/x_{0})}{<1/x_{0}>^{2}}=q-1,$ (22)

we can write that

$Var(1/x_{0})=(q_{in}-1)\langle 1/x_{0}\rangle^{2}\exp\left(-t/\tau\right).$
(23)

This therefore means that the relaxation time $\tau$ now describes the
temporal evolution of the fluctuations of the scale parameter $x_{0}$ (in the
scaled variable $z=x/x_{0}$).

##  4 Conclusions

A few remarks that may inspire further research in this direction could be of
interest here NatPhys . Notice that in the language of information theory
based on Shannon entropy,

$S=-\int dzf(z)ln[f(z)],$ (24)

our time evolution Eq. (14) can be expressed as

$\frac{\partial f}{\partial t}=-f\frac{\partial S_{z}}{\partial t}$ (25)

where the $S_{z}=-\ln[f(z)]$ is suprisal, which measures the information
gained by observing the outcome $z$ in the system and the Shannon entropy is
its expectation value. Now note that for the Tsallis distribution (9) we have
that $S_{z}\sim z$ and $\partial S_{z}/\partial t\sim z$. The entropy rate is
given by

$\frac{\partial S}{\partial t}=-\left<S_{z}\frac{\partial S_{z}}{\partial
t}\right>.$ (26)

From NatPhys  we know that the linear relationship between $\partial
S_{z}/\partial t$ and $z$ guaranties that distribution (9) saturates the time-
information uncertainty bound:

$\Bigg{|}Cov\left(\frac{\partial S}{\partial
t},z\right)\Bigg{|}\leq\sqrt{Var\left(\frac{\partial S}{\partial
t}\right)Var(z)}.$ (27)

More precisely, for the Tsallis distribution

$Cov\left(\frac{\partial S_{z}}{\partial t},z\right)=\left<\frac{\partial
S_{z}}{\partial t}\cdot z\right>-\left<\frac{\partial S_{z}}{\partial
t}\right>\langle z\rangle=0$ (28)

because $\langle\partial S_{z}/\partial t\cdot z\rangle=0$ and (for any
distribution) $\langle\partial S_{z}/\partial t\rangle=0$. The distance of a
given distribution (in our case defined by $q(t)$) from the equilibrium
distribution (defined by $q=1$) is given by the difference of the
corresponding entropies:

$\displaystyle S[q(t)]\\!-\\!S(q=1)\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\frac{q(t)\\!-1\\!}{1\\!-\\![q(t)-1)]}-\ln\\{1\\!-\\![q(t)\\!-\\!1]\\}\simeq$
(29) $\displaystyle\simeq$ $\displaystyle\frac{2[q(t)-1]}{2-q(t)}.$

From equations (16), (17) and (29) we get that

$\frac{df(t)}{dt}=\frac{f(t)}{\tau}\cdot\frac{q(t)-1}{2-q(t)}=\frac{f(t)}{\tau}\cdot\frac{S[q(t)]-S(q=1)]}{2}.$
(30)

In conclusion, we propose a new, modified form for the relaxation time
approximation for the collision term in the Boltzmann equation (2) allowing a
smooth transition to the thermalized distribution. It consists of replacing
the simple form of this term, given by the Eq. (3), where the relaxation time
$\tau$ determines how fast the equilibrium state of the studied distribution
$f$ is reached, by Eq. (17) describing the time evolution of the most
important (control) parameter of the analyzed distribution. The relaxation
time $\tau$ would now control the rate of change of this parameter from some
value to one that corresponds to the equilibrium state. We argue that this is
possible if we use for the phenomenological description of the distributions
of interest the quasi-power law Tsallis distribution given by Eq. (9) which is
able to describe a given process at all stages of time evolution. Its control
parameter is the time-dependent non-extensivity parameter $q(t)$ and the
relaxation time parameter $\tau$ describes its time evolution, as shown in Eq.
(17). This single quasi-power law probability distribution (9) smoothly
evolving towards thermalization would then replace the two-component
distribution given by Eq. (7) which arises from the RTA. The proposed scheme
offers multiple applications in situations where one wants to study the time
evolution of an ensemble but one does not want to invoke the kinetic theory
with complicated collision integrals.

###### Acknowledgements.

This research was supported by the Polish National Science Centre grant
UMO-2016/22/M/ST2/00176 and by the Polish Ministry of Science and Higher
Education grant DIR/WK/2016/2018/17-1 (GW). We thank Dr Nicholas Keeley for
reading the manuscript.

## References

  * (1) K.Huang, Statistical Mechanics, John Wiley and Sons, New York, 1963.
  * (2) P.L.Bhatnagar, E.P. Gross, M.Krook, Phys. Rev. 94, 511 (1954).
  * (3) J.L. Anderson, H.R. Witting, Physica 14, 466 (1974).
  * (4) W. Florkowski, R. Ryblewski, Phys. Rev. C 93, 064903 (2016).
  * (5) T. Bhattacharyya, P. Garg, R. Sahoo, P. Samantray, Eur. Phys. J. A 52, 283 (2016).
  * (6) M. Younus, S. Tripathy, S.K. Towar, R.Sahoo, Adv. High. En. Phys. ID 4728649 (2020).
  * (7) H-F. Zhao, B-C. Li, H-W. Dong, Adv. High En. Phys. ID3724761 (2020).
  * (8) S. Tripathy, A. Khuntia, S.K. Tiwari, R. Sahoo, Eur. Phys. J. A 53, 99 (2017).
  * (9) S. Tripathy, T. Bhattacharyya, P. Garg, P. Kumar, R. Sahoo, J. Cleymans, Eur. Phys. J. A 52, 289 (2016).
  * (10) L. Qiao, G. Che, J. Gu, H. Zheng, W. Zhang, J. Phys. G 47, 075101 (2020).
  * (11) S.B.Nicholson, L.P. Garcia-Pintos, A. del Campo, J.R. Green, Nat. Phys. 16, 1211 (2020).
  * (12) C. Tsallis, J. Stat. Phys. 52, 479 (1988).
  * (13) C. Tsallis, Eur. Phys. J. A 40, 257 (2009).
  * (14) C. Tsallis, Introduction to Nonextensive Statistical Mechanics (Springer, Berlin, 2009). For an updated bibliography on this subject, see http://tsallis.cat.cbpf.br/biblio.htm.
  * (15) G. Wilk, Z. Włodarczyk, Eur. Phys. J. A 48, 161 (2012).
  * (16) C-Y. Wong, G. Wilk, L.J.L. Cirto, C. Tsallis, Phys. Rev. D 91, 114027 (2015).
  * (17) G. Wilk, Z. Włodarczyk, Int. J. Mod. Phys. A 33, 193008 (2018).
  * (18) T. Bhattacharyya, Non-extensive Boltzmann Transport Equation: the Relaxation Time Approximation and Beyond, arXiv: 2104.05949 (2021).
  * (19) I. Prigogine, Introduction To Thermodynamics of Irreversible Processes, Edited by C. C.Thomas, Sprinfild, Geneseo, Ill, USA, 1955.
  * (20) V. N. Pokrovski, ISRN Thermodynamics, ID 906136 (2013), http://dx.doi.org/10.1155/2013/906136.
  * (21) G. Wilk, Z. Włodarczyk, Phys. Rev. Lett. 84, 2770 (2000).


Further author information contact Holger Roth<EMAIL_ADDRESS>or Ronald
Summers<EMAIL_ADDRESS>
# Deep convolutional networks for automated detection of posterior-element
fractures on spine CT

Holger R. Roth  Imaging Biomarkers and Computer-Aided Diagnosis Laboratory,
Radiology and Imaging Sciences, National Institutes of Health Clinical Center,
Bethesda, MD 20892-1182, USA.  Yinong Wang  Imaging Biomarkers and Computer-
Aided Diagnosis Laboratory, Radiology and Imaging Sciences, National
Institutes of Health Clinical Center, Bethesda, MD 20892-1182, USA.  Jianhua
Yao  Imaging Biomarkers and Computer-Aided Diagnosis Laboratory, Radiology and
Imaging Sciences, National Institutes of Health Clinical Center, Bethesda, MD
20892-1182, USA.  Le Lu  Imaging Biomarkers and Computer-Aided Diagnosis
Laboratory, Radiology and Imaging Sciences, National Institutes of Health
Clinical Center, Bethesda, MD 20892-1182, USA.  Joseph E. Burns  Department of
Radiological Sciences, University of California-Irvine, Orange, CA 92868, USA.
Ronald M. Summers  Imaging Biomarkers and Computer-Aided Diagnosis Laboratory,
Radiology and Imaging Sciences, National Institutes of Health Clinical Center,
Bethesda, MD 20892-1182, USA.

###### Abstract

Injuries of the spine, and its posterior elements in particular, are a common
occurrence in trauma patients, with potentially devastating consequences.
Computer-aided detection (CADe) could assist in the detection and
classification of spine fractures. Furthermore, CAD could help assess the
stability and chronicity of fractures, as well as facilitate research into
optimization of treatment paradigms.

In this work, we apply deep convolutional networks (ConvNets) for the
automated detection of posterior element fractures of the spine. First, the
vertebra bodies of the spine with its posterior elements are segmented in
spine CT using multi-atlas label fusion. Then, edge maps of the posterior
elements are computed. These edge maps serve as candidate regions for
predicting a set of probabilities for fractures along the image edges using
ConvNets in a 2.5D fashion (three orthogonal patches in axial, coronal and
sagittal planes). We explore three different methods for training the ConvNet
using 2.5D patches along the edge maps of ‘positive’, i.e. fractured
posterior-elements and ‘negative’, i.e. non-fractured elements.

An experienced radiologist retrospectively marked the location of 55 displaced
posterior-element fractures in 18 trauma patients. We randomly split the data
into training and testing cases. In testing, we achieve an area-under-the-
curve of 0.857. This corresponds to 71% or 81% sensitivities at 5 or 10 false-
positives per patient, respectively. Analysis of our set of trauma patients
demonstrates the feasibility of detecting posterior-element fractures in spine
CT images using computer vision techniques such as deep convolutional
networks.

###### keywords:

spine CT, computer-aided detection, posterior-element fractures, deep learning

##  1 INTRODUCTION

Injuries of the spine and its posterior elements in particular (see Figure 1)
are a common occurrence in traumatic patients, with potentially devastating
consequences [1]. Spine fractures are detected using volumetric imaging such
as computed tomography (CT) in order to access the degree of injury. Spine
injuries are a critical concern in blunt trauma, particularly in cases of
motor vehicle collision and fall from significant heights. More than 140,000
vertebral fractures occur in the U.S. each year [2]. However, the traditional
method of qualitative visual assessment of images for diagnosis could miss
fractures, and is time-consuming, potentially causing delays in time-critical
situations such as the treatment of spine injuries. Computer-aided detection
(CADe) has the potential to expedite the assessment of trauma cases, reduce
the chance of misclassification of fractures of the spine, and decrease inter-
observer variability. Furthermore, CADe could help assess the stability and
chronicity of fractures, as well as facilitate research into optimization of
treatment paradigms.

---
Figure 1:  Example 3D rendering of a posterior element fracture of the spinous
process (left) and reformatted images (right) show a posterior element
fracture of the transverse process in axial and sagittal views. Manual gold
standard marks by radiologist are shown in red.

##  2 METHODS

Recently, the availability of large annotated training sets and the
accessibility of affordable parallel computing resources via GPUs has made it
feasible to train deep convolutional networks (ConvNets), also popularized
under the keyword “deep learning”, for computer vision classification tasks.
Great advances in classification of natural images have been achieved [3, 4].
Studies that have tried to apply deep learning and ConvNets to medical imaging
applications also showed promise [5, 6, 7, 8].

###  2.1 Convolutional networks

ConvNets are named for their convolutional filters which are used to compute
image features for classification [3]. In this work, we use a ConvNet similar
to [3]111https://code.google.com/archive/p/cuda-convnet2/ with 5 convolutional
layers, 3 fully connected layers, and a final softmax layer for classificaton.
However, we use input images of size $64\times 64$ and a stride of 1 in the
first convolutional layer. All convolutional filter kernel elements are
trained from the data in a fully-supervised fashion. This has major advantages
over more traditional CADe approaches that use hand-crafted features, designed
from human experience. ConvNets have a better chance of capturing the
“essence” of the imaging data set used for training than when using hand-
crafted features [3, 4]. Furthermore, we can train the similarly configured
ConvNet architectures from randomly initialized or pre-trained model
parameters for detecting different lesions or pathologies (with heterogeneous
appearances), with no manual intervention of system and feature design.
Examples of trained filters from the first convolutional layer for posterior-
element fracture detection are shown in Fig. 2. Dropout [9] is used during
training as a form of regularization that avoids overfitting to the training
data.

---
Figure 2:  Convolution filters trained on examples of posterior-element
fractures and non-fractures.

###  2.2 Application to spine CT

In this work, we apply ConvNets for the automated detection of posterior
element fractures of the spine. First, the vertebra bodies of the spine with
its posterior elements are segmented in spine CT using multi-atlas label
fusion [10, 11]. A set of atlases of the vertebra bodies are registered to the
target image with free-form deformation. Then, an edge map $E$ of the
posterior elements is computed using the Sobel operators

$S_{x}=\begin{Bmatrix}-1&0&+1\\\ -2&0&+2\\\ -1&0&+1\end{Bmatrix}\ \
\mathrm{and}\ \ S_{y}=\begin{Bmatrix}-1&-2&-1\\\ 0&0&0\\\
+1&+2&+1\end{Bmatrix}$ (1)

in order to find the horizontal ($G_{x}=S_{x}*I$) and vertical
($G_{y}=S_{x}*I$) approximations of the image derivative for each CT slice.
Here, $*$ denotes a convolutional operation. Edge points are located at the
maximum of the absolute gradient $\left|G\right|=\sqrt{G_{x}^{2}+G_{y}^{2}}$.
The edge maps $E$ serve as candidate regions for predicting a set of
probabilities $P$ for fractures along an image edge using ConvNets. An example
of posterior-element segmentation and edge map estimation is shown in Fig. 3.

---
Figure 3:  (a) Posterior element fracture; (b) posterior element segmentation
using multi-atlas label fusion[10]; (c) edge map estimation.

We explore three different methods for training the ConvNet using 2.5D patches
(three orthogonal patches in axial, coronal and sagittal planes) along the
edge maps of ‘positive’, i.e. fractured posterior-elements and ‘negative’,
non-fractured elements:

  1. 1.

$P_{\mathrm{original}}$: 2.5D patches along each edge voxel $e\in E$ are
aligned to original scanner coordinates in axial, coronal, and sagittal
planes.

  2. 2.

$P_{\mathrm{mirrored}}$: 2.5D patches at voxels $e$ near fractures are
mirrored along the $y$-axis of the image in order to increase ‘positive’
training examples.

  3. 3.

$P_{\mathrm{oriented}}$: 2.5D patches are oriented along the principal axis of
an edge at voxel $e$ (while ‘positive’ training examples are handled as in
case 2, i.e. they are also mirrored).

The alignment of a training patch along its edge orientation is illustrated in
Fig. 4. This step can be also seen as a form of data augmentation that
artificially increases the variation of training examples (as opposed to just
having slightly shifted versions of nearby patches as in case 1 and 2).
Orientating a patch along the anatomical structures of interest has also been
shown to improve performance in other applications such as in detection of
pulmonary embolism [12] where the image patch can be aligned along the
direction of a blood vessel. The orientation of an edge is estimated as in
Equation 2. In 2D, the eigenvector $v_{1}$ corresponding to the largest
eigenvalue $\lambda_{1}$ corresponds to the major axis of a local edge. Hence,
the orientation $\theta$ of the edge can be estimated as the angle between
this eigenvector and the origin $v_{0}$:

$\theta=\tan^{-1}\left(\frac{v_{1}}{v_{0}}\right)$ (2)

---
Figure 4:  Illustration of 2D patches aligned along the candidate edge
orientation with respect to the origin for training convolutional networks
(ConvNets).

##  3 RESULTS

###  3.1 Data

An experienced radiologist retrospectively marked the location of 55 displaced
posterior-element fractures in 18 trauma patients admitted for traumatic
emergencies to the University of California-Irvine Medical Center. Image sizes
range within $512\times 512\times[259-582]$ and resolutions range within
$[0.29-0.41]\times[0.29-0.41]\times[1.0-2.0]$ mm. We use a random subset of 12
of these patients with spine fractures for training ConvNets as described in
Section 2; 6 patients are reserved for testing. An additional set of 5 spine
CTs of healthy patients were added to the training set in order to increase
the number of non-fractured examples. A total of 800,000 2.5D patches are
randomly selected from the candidate edge maps $\left\\{E\right\\}$ of the
training set and used for learning the ConvNet parameters. After convergence,
the ConvNet is applied to a testing case edge map E in order to produce a
probability map $P$ for fractures. Figure 5 shows examples of probability maps
for posterior-element fracture detection along candidate edges. The
radiologist’s markings are indicated by crosshairs.

---
Figure 5:  Examples of probability maps for posterior-element fracture
detection along candidate edges. The posterior element fractures are marked by
crosshairs.

###  3.2 Performance

The classification performance is evaluated using voxel-wise ROC and per
posterior-element process (left, right and spinous process) FROC analysis.
Fig. 6 shows the voxel-wise ROC and processes-wise FROC performance using the
differently trained ConvNets. A clear advantage of using ‘Oriented’ patches
for training can be observed in testing, compared to ‘Original’ and ‘Mirrored’
patches for training. In testing, we achieve area-under-the-curve (AUC) values
of 0.761, 0.796 and 0.857 for ‘Original’, ‘Mirrored’, and ‘Oriented’ patch
training respectively. This corresponds to 71% or 81% sensitivities at 5 or 10
false-positives per patient, respectively, in the case of training and testing
with ‘Oriented’ patches.

---
Figure 6:  Voxelwise-ROC performance (left) and per-fracture FROC performance
(right) evaluating the classification accuracy.

##  4 CONCLUSION

While there has been work on automated detection of other spinal injuries,
such as vertebral compression fractures [2], to the best of our knowledge,
this is the first work to explore automated posterior-element fracture
detection. It demonstrates that deep convolutional networks (ConvNets) can be
useful for detection tasks, such as the detection of fractures in spine CT.

Results from analysis of our set of trauma patients demonstrate the
feasibility of detecting posterior-element fractures in spine CT images using
computer vision techniques such as ConvNets. We evaluated three methods for
training a ConvNet and show the advantages of applying an oriented patch
extraction method ($P_{\mathrm{oriented}}$) for better classification
performance.

###### Acknowledgements.

This work was supported by the Intramural Research Program of the NIH Clinical
Center.

## References

  * [1] Parizel, P., van der Zijden, T., Gaudino, S., Spaepen, M., Voormolen, M., Venstermans, C., De Belder, F., van den Hauwe, L., and Van Goethem, J., “Trauma of the spine and spinal cord: imaging strategies,” European Spine Journal 19(1), 8–17 (2010).
  * [2] Yao, J., Burns, J. E., Muñoz, H., and Summers, R. M., “Cortical shell unwrapping for vertebral body abnormality detection on computed tomography,” Computerized Medical Imaging and Graphics 38(7), 628–638 (2014).
  * [3] Krizhevsky, A., Sutskever, I., and Hinton, G. E., “Imagenet classification with deep convolutional neural networks,” in [Advances in neural information processing systems ], 1097–1105 (2012).
  * [4] Zeiler, M. D. and Fergus, R., “Visualizing and understanding convolutional networks,” in [Computer Vision–ECCV 2014 ], 818–833, Springer (2014).
  * [5] Yan, Z., Zhan, Y., Peng, Z., Liao, S., Shinagawa, Y., Metaxas, D. N., and Zhou, X. S., “Bodypart recognition using multi-stage deep learning,” in [Information Processing in Medical Imaging ], 449–461, Springer (2015).
  * [6] Cireşan, D. C., Giusti, A., Gambardella, L. M., and Schmidhuber, J., “Mitosis detection in breast cancer histology images with deep neural networks,” in [Medical Image Computing and Computer-Assisted Intervention–MICCAI 2013 ], 411–418, Springer (2013).
  * [7] Roth, H. R., Lu, L., Seff, A., Cherry, K. M., Hoffman, J., Wang, S., Liu, J., Turkbey, E., and Summers, R. M., “A new 2.5 d representation for lymph node detection using random sets of deep convolutional neural network observations,” in [Medical Image Computing and Computer-Assisted Intervention–MICCAI 2014 ], 520–527, Springer (2014).
  * [8] Roth, H. R., Lu, L., Liu, J., Yao, J., Seff, A., Kevin, C., Kim, L., and Summers, R. M., “Improving computer-aided detection using convolutional neural networks and random view aggregation,” arXiv preprint arXiv:1505.03046 (2015).
  * [9] Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R., “Dropout: A simple way to prevent neural networks from overfitting,” The Journal of Machine Learning Research 15(1), 1929–1958 (2014).
  * [10] Wang, Y., Yao, J., Roth, H. R., Burns, J. E., and Summers, R. M., “Multi-atlas segmentation with joint label fusion of osteoporotic vertebral compression fractures on ct,” Recent Advances in Computational Methods and Clinical Applications for Spine Imaging (2015).
  * [11] Wang, H., Suh, J. W., Das, S. R., Pluta, J. B., Craige, C., Yushkevich, P., et al., “Multi-atlas segmentation with joint label fusion,” Pattern Analysis and Machine Intelligence, IEEE Transactions on 35(3), 611–623 (2013).
  * [12] Tajbakhsh, N., Gotway, M. B., and Liang, J., “Computer-aided pulmonary embolism detection using a novel vessel-aligned multi-planar image representation and convolutional neural networks,” in [Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015 ], 62–69, Springer (2015).