# DeepOpht: Medical Report Generation for Retinal Images via Deep Models and
Visual Explanation

Jia-Hong Huang1†, C.-H. Huck Yang2‡,3, Fangyu Liu4, Meng Tian5, Yi-Chieh
Liu3,9, Ting-Wei Wu3,7,
I-Hung Lin6, Kang Wang8, Hiromasa Morikawa2, Hernghua Chang9, Jesper Tegner2,
Marcel Worring1
1University of Amsterdam, 2King Abdullah University of Science and Technology
(KAUST), 3Georgia Institute of Technology,
4University of Cambridge, 5Department of Ophthalmology, Bern University
Hospital, 6Tri-Service General Hospital,
7University of California, Berkeley, 8Beijing Friendship Hospital, and
9National Taiwan University
$\dagger$ corresponding author<EMAIL_ADDRESS>$\ddagger$ work conducted
during visiting KAUST.

###### Abstract

In this work, we propose an AI-based method that intends to improve the
conventional retinal disease treatment procedure and help ophthalmologists
increase diagnosis efficiency and accuracy. The proposed method is composed of
a deep neural networks-based (DNN-based) module, including a retinal disease
identifier and clinical description generator, and a DNN visual explanation
module. To train and validate the effectiveness of our DNN-based module, we
propose a large-scale retinal disease image dataset. Also, as ground truth, we
provide a retinal image dataset manually labeled by ophthalmologists to
qualitatively show, the proposed AI-based method is effective. With our
experimental results, we show that the proposed method is quantitatively and
qualitatively effective. Our method is capable of creating meaningful retinal
image descriptions and visual explanations that are clinically relevant.
DeepOpht Github.

##  1 Introduction

The World Health Organization (WHO) estimates that typical retinal diseases
such as Age-related Macular Degeneration (AMD) and Diabetic Retinopathy (DR)
are expected to affect over 500 million people worldwide shortly [49].
Besides, generally speaking, the traditional process of retinal disease
diagnosis and creating a medical report for a patient takes time in practice.
The above means that ophthalmologists will become busier and busier.

As we may know, the current state of the art in Artificial Intelligence (AI)
involves deep learning research, and we claim deep learning is one of the
promising ways to help ophthalmologists and improve the traditional retinal
disease treatment procedure. Deep learning based models such as convolutional
neural networks (CNN) or recurrent neural networks (RNN) for computer vision
or natural language processing tasks, respectively, have achieved, and, in
some cases, even exceeded human-level performance. There is no better time
than now to propose an AI-based medical diagnosis method to aid
ophthalmologists.

Figure 1: This figure shows the proposed AI-based medical diagnosis method in
the ophthalmology expert domain. It contains DNN-based and DNN Visual
explanation modules. The DNN-based module is composed of two sub-modules,
i.e., a retinal disease identifier and a clinical description generator
reinforced by our proposed keyword-driven method, referring to our Methodology
section. The input of our method is a retinal image, and the output is a
table-based [68] medical report. In Figure 2, we shows how to exploit this AI-
based method to improve the traditional retinal diseases treatment procedure.
Note that, in this figure, DNN indicates deep neural networks.

In this paper, we propose an AI-based method for automatic medical report
generation based on an input retinal image, as illustrated in Figure 1. The
proposed method intends to improve the traditional retinal disease diagnosis
procedure, referring to Figure 2, and help ophthalmologists increase diagnosis
efficiency and accuracy. The main idea of this method is to exploit the deep
learning based models, including an effective retinal disease identifier (RDI)
and an effective clinical description generator (CDG), to automate part of the
traditional treatment procedure. Then, the proposed method will make the
diagnosis more efficient.

To train our deep learning models and validate the effectiveness of our RDI
and CDG, we introduce a new large-scale retinal disease image dataset, called
DeepEyeNet (DEN). Besides, as ground truth, we provide a retinal image dataset
manually labeled by ophthalmologists to qualitatively show that the proposed
AI-based model is effective. The dataset helps us show the activation maps of
our deep models are aligned with image features that are clinically recognized
by ophthalmologists as linked with the identified disease. Our experimental
results show that the proposed AI-based method is effective and successfully
improves the traditional retinal disease treatment procedure. Our main
contributions are summarized as follows:

Contributions.

  * •

To improve the traditional retinal disease treatment procedure and help
ophthalmologists increase diagnosis efficiency and accuracy, we propose an AI-
based method to generate medical reports for retinal images. In this method,
we exploit the deep learning based models including an RDI and a CDG to
automate part of the conventional treatment procedure.

  * •

We propose a large-scale retinal disease image dataset, called DeepEyeNet
(DEN) dataset, with 15,709 images to train our deep models and validate the
effectiveness of the proposed RDI and CDG quantitatively.

  * •

We provide another dataset with 300 retinal images labeled by ophthalmologists
to qualitatively show our method is effective by visually confirming the
activation maps of our models are aligned with image features clinically
recognized by ophthalmologists.

Figure 2:  (a) is an existing traditional medical treatment procedure for
retinal diseases [56]. Typically, doctors have to handle most of the jobs in
the traditional procedure. In (b), we incorporate the AI-based medical
diagnosis method, referring to Figure 1, in the traditional treatment
procedure to improve the efficiency of (a), based on the point-of-care (POC)
[47] concept. In the proposed method, it mainly contains DNN-based and DNN
visual explanation modules. The outputs of the DNN-based module are “Disease
class “A”” and “Clinical description”. DNN visual explanation module visualize
the information from the DNN-based module on the classification tasks. Please
refer to our Methodology section for a more detailed explanation. Note that
DNN indicates deep neural networks in this figure.

##  2 Related Work

In this section, we divide the related works into retinal disease
classification, image captioning, neural networks visual explanation, and
retinal dataset comparison.

2.1 Retinal Disease Classification

Optical Coherence Tomography (OCT), Fluorescein Angiography (FA), and Color
Fundus Photography (CFP) are the three most commonly used and important
imaging methods for the diagnosis of retinal diseases [61]. Optical Coherence
Tomography (OCT) is a technology of emerging biomedical imaging, and it
provides high-resolution and non-invasive real-time imaging of highly
scattering tissues. That is, OCT images [11, 33, 17] usually are used to show
the structure of the retina. [5] have proposed an algorithm to segment and
detect six different retinal layers, including Nerve Fiber Layer (NFL),
Ganglion Cell Layer (GCL) + Inner Plexiform Layer (IPL), Inner Nuclear Layer
(INL), Outer Plexiform Layer (OPL), Outer Nuclear Layer (ONL) + Photoreceptor
Inner Segments (PIS), and Photoreceptor Outer Segments (POS), in OCT retinal
images. Fluorescein Angiography (FA) has been used to realize the
pathophysiologic course of Retinopathy of Prematurity (ROP) following
intravitreal anti-Vascular Endothelial Growth Factor (Anti-VEGF) [35]. Color
Fundus photography (CFP) is a simple and cost-effective technology for trained
medical professionals. Image preprocessing is one of the important issues in
the automated analysis of CFP. The authors of [67] have proposed a method to
reduce the vignetting effect caused by non-uniform illumination of a retinal
image. In this work, we mainly exploit the DNN-based methods [18, 53, 65] to
further investigate the retinal disease classification [62] toward a multi-
label task combining language information.

2.2 Image Captioning

Recently, computer vision researchers have proposed a new task, image
captioning, and [30, 59, 14] are early works. In [30], the proposed model can
embed visual and language information into a common multimodal space. The
authors of [14] exploit a natural language model to combine a set of possible
words, which are related to several small parts of the image, and then
generate the caption of the given image. The authors of [59] use CNN to
extract the image feature and use it as the input at the first time step of
the RNN to generate the caption of the input image. The authors of [16]
propose a new deliberate residual attention network for image captioning. The
layer of first-pass residual-based attention prepares the visual attention and
hidden states for generating a preliminary version of the captions, while the
layer of second-pass deliberate residual-based attention refines them. Since
the second-pass is based on the global features captured by the hidden layer
and visual attention in the first-pass, their method has the potentials to
generate better captions. In [40], the authors mention that existing image
captioning models are usually trained via maximum likelihood estimation.
However, the log-likelihood score of some captions cannot correlate well with
human assessments of quality. Standard syntactic text evaluation metrics, such
as METEOR [6], BLEU [48], and ROUGE [38], are also not well correlated. The
authors of [40] show how to use a policy gradient method to optimize a linear
combination of CIDEr [58] and SPICE [3]. In [19], the authors propose a method
that focuses on discriminating properties of the visible object, jointly
predicts a class label, and explains why the predicted label is proper for a
given image. Through a loss function based on reinforcement learning and
sampling, their model learns to generate captions. According to [59, 30, 16],
existing image captioning models are only able to generate the rough
description for a given image. So, in this work, we exploit keywords to make
our CDG have better reasoning ability.

Table 1: Summary of available retinal datasets. Based on this table, we find
our proposed DEN is much larger than the other retinal image datasets. It
contains three types of labels including the name of the disease, keywords,
and clinical description. Most of the retinal dataset only contains image
data, and the dataset size is not large. Note that “ Text*” denotes clinical
description and keywords, referring to our Dataset Introduction and Analysis
section. “ Text” denotes only clinical description. So, our DEN is unique.

Name of Dataset |  Field of View |  Resolution |  Data Type |  Number of Images
---|---|---|---|---
VICAVR [57] |  $45^{\circ}$ |  $768*584$ |  Image |  58
VARIA [46] |  $20^{\circ}$ |  $768*584$ |  Image |  233
STARE [21] |  $\approx 30^{\circ}-45^{\circ}$ |  $700*605$ |  Image + Text |  397
CHASE-DB1 [15] |  $\approx 25^{\circ}$ |  $999*960$ |  Image |  14
RODREP [1] |  $45^{\circ}$ |  $2000*1312$ |  Image |  1,120
HRF [45] |  $45^{\circ}$ |  3504*2336 |  Image |  45
e-ophtha [12] |  $\approx 45^{\circ}$ |  $2544*1696$ |  Image |  463
ROC [43] |  $\approx 30^{\circ}-45^{\circ}$ |  $768*576-1386*1391$ |  Image |  100
REVIEW [2] |  $\approx 45^{\circ}$ |  $1360*1024-3584*2438$ |  Image |  14
ONHSD [10] |  $45^{\circ}$ |  $640*480$ |  Image |  99
INSPIRE-AVR [44] |  $30^{\circ}$ |  $2392*2048$ |  Image |  40
DIARETDB1 [32] |  $50^{\circ}$ |  $1500*1152$ |  Image + Text |  89
DIARETDB0 [31] |  $50^{\circ}$ |  $1500*1152$ |  Image |  130
MESSIDOR [13] |  $45^{\circ}$ |  $1440*960-2304*1536$ |  Image + Text |  1,200
Drishti-GS [54] |  $\approx 25^{\circ}$ |  $2045*1752$ |  Image |  101
FIRE [20] |  $45^{\circ}$ |  $2912*2912$ |  Image |  129
DRIONS-DB [7] |  $\approx 30^{\circ}$ |  $600*400$ |  Image |  110
IDRiD [50] |  $50^{\circ}$ |  $4288*2848$ |  Image |  516
DRIVE [55] |  $45^{\circ}$ |  $565*584$ |  Image |  40
DeepEyeNet (DEN) |  $\approx\textbf{30}^{\circ}-\textbf{60}^{\circ}$ |  various |  Image + Text* |  15,709

2.3 Neural Networks Visual Explanation

There are some popular CNN visualization tools, [70, 52, 64]. The authors of
[70] have proposed a technique, called Class Activation Mapping (CAM), for
CNN. It makes classification-trained CNN learn how to perform the task of
object localization, without using a bounding box. In our previous works [42,
63], we exploit class activation maps to visualize the predicted class scores
on retinal images, highlighting the discriminative object parts which are
detected by the CNN. In [52], the authors have proposed the other similar
features visualization tool, called Gradient-weighted Class Activation Mapping
(Grad-CAM), for making a CNN-based model transparent by producing visual
explanations of features. The authors of [69] introduce a CNN visualization
technique that gives insight into the operation of the classifier and the
function of intermediate feature layers. These visualizations allow us to find
architectures of CNN models. The authors of [8] propose a generalized method,
Grad-CAM++, based on Grad-CAM. The Grad-CAM++ method provides better visual
explanations of CNN model predictions than Grad-CAM, in terms of better object
localization [66] and occurrences explanation of multiple object instances in
a single image. In [37], the authors propose another method different from the
above methods which are trying to explain the network. They build up an end-
to-end model to provide supervision directly on the visual explanations.
Furthermore, the authors validate that the supervision can guide the network
to focus on some expected regions. The aforementioned is more related to image
data only visualization. The authors of [68, 51, 28] have proposed some
methods for the multimedia data, such as text and images, visualization. In
[68], the authors introduce five popular multimedia visualization concepts,
including basic grid, similarity space, similarity-based, spreadsheet, and
thread-based concepts. In this work, we exploit CAM to visually show that the
activation maps of our deep models are aligned with image features that are
clinically recognized by ophthalmologists as linked with the identified
disease. In addition, we use a table-based concept, similar to the static
spreadsheet concept, to visualize our medical report.

2.4 Retinal Dataset Comparison

Retinal disease research already has long history and many retinal datasets
have been proposed, such as [55, 50, 7, 20, 54, 13, 31, 32, 44, 10, 2, 43, 12,
45, 1, 15, 21, 46, 57]. The DRIVE dataset [55] contains 40 retina images which
are obtained from a diabetic retinopathy screening program in the Netherlands.
These 40 images have been divided into a half training set and a half test
set. For the training images, a single manual segmentation of the vasculature
is available. For the test cases, two manual segmentations are available. The
IDRiD dataset [50] is a dataset for retinal fundus image consisting of 516
images. The authors of IDRiD dataset provide ground truths associated with the
signs of Diabetic Macular Edema (DME) and Diabetic Retinopathy (DR) and normal
retinal structures given below and described as follows: (i) Pixel level
labels of typical DR lesions and optic disc; (ii) Image level disease severity
grading of DR, and DME; (iii) Optic disc and fovea center coordinates. The
DRIONS-DB dataset [7] consists of 110 color digital retinal images, and it
contains several visual characteristics, such as cataract (severe or
moderate), light artifacts, some of the rim blurred or missing, moderate
peripapillary atrophy, concentric peripapillary atrophy/artifacts, and strong
pallor distractor. The FIRE dataset [20] consists of 129 retinal images
forming 134 image pairs, and image pairs are split into three different
categories depending on their characteristics. The Drishti-GS dataset [54]
contains 101 images, and it is divided into 50 training and 51 testing images.
The MESSIDOR dataset [13] has 1200 eye fundus color numerical images. Although
the dataset contains a medical diagnosis for each image, there is no manual
annotation, such as lesions contours or position, on the images. The DIARETDB0
dataset [31] consists of 130 color fundus images of which 20 are normal, and
110 contain signs of the DR. The DIARETDB1 dataset [32] consists of 89 color
fundus images of which 84 contain at least mild non-proliferative signs of the
DR, and five are considered as normal which do not contain any signs of the
DR. The INSPIRE-AVR dataset [44] has 40 colorful images of the vessels and
optic disc and an arterio-venous ratio reference standard. The ONHSD dataset
[10] has 99 retinal images and it is mainly used for the segmentation task.
The REVIEW dataset [2] consists of 14 images, and it is also mainly used for
the segmentation task. The ROC dataset [43] aims to help patients with
diabetes through improving computer-aided detection and diagnosis of DR. The
e-ophtha [12] is a dataset of color fundus images specially designed for
scientific research in DR. The HRF dataset [45] contains at the moment 15
images of healthy patients, 15 images of patients with DR and 15 images of
glaucomatous patients. Also, binary gold standard vessel segmentation images
are available for each image. The RODREP dataset [1] contains repeated 4-field
color fundus photos (1120 in total) of 70 patients in the DR screening program
of the Rotterdam Eye Hospital. The CHASE-DB1 dataset [15] is mainly used for
retinal vessel analysis, and it contains 14 images. The STARE dataset [21] has
397 images and it is used to develop an automatic system for diagnosing
diseases of the human eye. The VARIA [46] is a dataset of retinal images used
for authentication purposes, and it includes 233 images from 139 different
individuals. The VICAVR dataset [57] includes 58 images, and it is used for
the computation of the ratio of A/V, (Artery/Vein). In this work, we propose a
large-scale retinal images dataset, DeepEyeNet (DEN), to train our deep
learning based models and validate our RDI and CDG. For convenience, we
summarize the above retinal datasets in Table 1.

##  3 Dataset Introduction and Analysis

In this section, we start to describe our proposed DEN dataset in terms of
types of retinal images and labels and some statistics of the dataset. Note
that some of our group members are experienced ophthalmologists and they help
us build the proposed DEN dataset sorted by 265 unique retinal symptoms from
the clinical definition and their professional domain knowledge. In our
proposed DEN dataset, there are two types of retinal images, grey scale FA and
colorful CFP. The total amount of images is 15,709, including 1,811 FA and
13,898 CFP. As with most of the large-scale datasets for deep learning
research, we create standard splits, separating the whole dataset into
$60\%/20\%/20\%$, i.e., $9425/3142/3142$, for training/validation/testing,
respectively. Each retinal image has three corresponding labels including the
name of the disease, keywords, and clinical description. For the total number
of retinal diseases, the dataset contains 265 different retinal diseases
including the common and non-common. For the keyword and clinical description,
it contains 15,709 captions and 15,709 keywords labels. Keyword label denotes
important information in the diagnosis process. Clinical description label
represents the corresponding caption of a given retinal image. Note that all
the labels are defined by retina specialists or ophthalmologists. To better
understand our dataset, we show some data examples from the DEN dataset in
Figure 3. Also, in Figure 4, we show the word length distribution of the
keyword and clinical description labels. Based on Figure 4, we observe the
longest word length in our dataset is more than $15$ words for keywords and
$50$ words for clinical descriptions. Note that the longest word length of
existing datasets for natural image captioning or VQA [4, 9, 39, 26, 23] is
only around $10$ words. It implies that our proposed dataset is challenging.
Additionally, we provide the Venn-style word cloud visualization results
clinical description labels, referring to Figure 5. Based on Figure 5, in
clinical description labels, we see there are specific abstract concepts,
which makes the dataset more challenging.

Figure 3: Examples from our DEN dataset. Each image has three labels including
the name of the disease, keywords, and clinical description. Note that
ophthalmologists define all the labels. Figure 4: This figure shows the word
length distribution of the keyword and clinical description labels. Based on
the figure, the word length in our DEN dataset is mainly between $5$ and $10$
words. Figure 5: The figure represents Venn-style word cloud for clinical
description labels. Note that the word size indicates the normalized counts.
Based on this figure, we can see there are specific abstract concepts, which
makes image captioning algorithms more difficult to generate descriptions with
good quality. Figure 6: This figure conceptually depicts the clinical
description generator with our proposed keyword-driven method. In our clinical
description generator, we exploit a pre-trained CNN model to extract the
retinal image feature. So, the CNN model is a so-called image encoder. Then,
we use an LSTM model, i.e., recurrent neural networks (RNN), as a decoder to
generate a word at each time step. Finally, all of the collected words will
form a clinical description.

##  4 Methodology

In this section, we start to describe the proposed AI-based method for
automatic medical report generation. The proposed method is mainly composed of
the DNN-based module and DNN visual explanation module.

4.1 DNN-based Module

The DNN-based module contains two components, i.e., a retinal disease
identifier (RDI) and a clinical description generator (CDG). We introduce them
in the following subsections. Note that we hypothesize an effective RDI and
effective CDG help improve the conventional retinal disease treatment
procedure and help ophthalmologists increase diagnosis efficiency and
accuracy.

Retinal Disease Identifier (RDI). To identify retinal diseases, in our RDI
sub-module, we provide two types of deep learning models based on [18, 53],
pre-trained on ImageNet, and then trained on the proposed DEN dataset. From
the lower level feature perspective, such as color, most of the medical
images, e.g., radiology images of the chest, are mainly grey-scale [34] but
retinal images are mainly colorful in our dataset. Using the ImageNet pre-
trained at least helps extract the better lower level features information.
So, in this case, we expect that pre-training on ImageNet can improve model
performance.

Clinical Description Generator (CDG). To generate the clinical description for
an input retinal image, we use a pre-trained CNN-based model, such as
MobileNetV2, VGG16, VGG19, or InceptionV3, as our image feature encoder and a
Long Short-term Memory (LSTM) as our decoder to generate text, referring to
Figure 6. When we try to generate the clinical description by the LSTM unit,
we incorporate the beam-search mechanism to get the better final output
description. In ophthalmological practice, commonly existing keywords, with
unordered nature, help ophthalmologists create medical reports. Inspired by
this, we exploit keywords to reinforce our CDG sub-module. As shown in Figure
6, we use a keyword embedding module, such as bag of words, to encode our
keyword information. Note that when keywords are used to reinforce CDG, it
means we will have two types of input features, i.e., image and text features.
In this work, we use the average method to fuse these two types of features,
referring to Figure 6.

4.2 DNN Visual Explanation Modules

There are some existing DNN visual explanation methods, such as [70, 52, 22,
41]. The authors of [70] have proposed a technique, called Class Activation
Mapping (CAM), for CNN. It makes classification-trained CNN learn how to
perform the task of object localization, without using a bounding box.
Furthermore, they exploit class activation maps to visualize the predicted
class scores on a given image, highlighting the discriminative object parts
which are detected by the CNN. To improve the conventional retinal disease
treatment procedure, we incorporate the DNN visual explanation module in our
proposed AI-based method. Also, we exploit this module to help verify the
effectiveness of the method, referring to our Experiments section.

Table 2:  This table shows the quantitative results of different RDI models
based on our DEN. The RDI model based on [53] with ImageNet pre-training has
the best performance. “Pre-trained” indicates the model is initialized from
the pre-trained weights of ImageNet. “Random init” means the model’s weights
are initialized randomly. Prec@k indicates how often the ground truth label is
within the top $k$ ranked labels after the softmax layer. We investigate
Prec@1 and Prec@5 due to the need to shortlist candidates of diseases in real-
world scenarios. Note that since we have $265$ retinal disease candidates and
limited training data, it is hard to have good performance in the sense of
Prec@1. The situation of limited data is common in medicine.

Model | Precision
---|---
Pre-trained | Random init
Prec@1 | Prec@5 | Prec@1 | Prec@5
He, et al. [18] | 37.09 | 63.36 | 36.60 | 62.87
Simonyan, et al. [53] | 54.23 | 80.75 | 35.93 | 73.73
Jing, et al. [29] | 32.72 | 63.75 | 29.11 | 60.68

4.3 Medical Report Generation

According to [68, 51], proper multimedia data visualization helps people get
insight from the data efficiently. In some sense, we can say that multimedia
visualization is a way to visually arrange multimedia data, and it sometimes
even helps people get a deeper understanding and extra information from the
visualized data. In this work, it contains five multimedia data, including the
name of the disease, keyword, clinical description, retinal image, and CAM
result image. So, we exploit the table-based concept, which is similar to the
static spreadsheet-based concept [68], to visualize our medical report,
referring to Figure 7. The medical report visualization intends to help
ophthalmologists get the insight from the above image and text data
efficiently and also increase the diagnostic accuracy.

Figure 7: This figure shows the medical reports based on the table-based
concept [68]. Since retinal diseases may have some implicit common property or
relation, we can put the diseases with the common property or relation
together on the table. The table-based medical report intends to help
ophthalmologists get more insights. Figure 8: This figure shows some generated
results by our clinical description generator. Based on this figure, we know
that our models can generate meaningful clinical descriptions for
ophthalmologists. Note that, in practice, “age” and “gender” are hard to be
generated correctly by automatic algorithms. The first row with correct “age”
prediction is just a special case.

##  5 Experiments

In this section, we compare our proposed method to baselines for verifying the
effectiveness based on the assumption described in our Methodology section.

5.1 Retinal Disease Identifier (RDI) Verification

In our experiment, we try to show that the RDI model with ImageNet pre-
training is better than the RDI model without ImageNet pre-training, i.e., our
baseline. We exploit the ImageNet-pre-trained DNN-based deep model and non-
ImageNet-pre-trained DNN-based deep model with different architectures to do
fine-tuning on DEN. For empirical reasons, we use two recipes to train
different models. For the RDI model based on [18], we start with a learning
rate of $0.1$ and decay it $5$ times for every $50$ epoch. For the RDI model
based on [53], we start with a learning rate of $0.001$ and decay it $5$ times
for every $50$ epoch. According to the evaluation results in Table 2, we find
that the RDI model based on [53] with ImageNet pre-training has better
performance than others. We conjecture that RDI models based on [18, 29] may
be too complicated for the proposed DEN dataset. Although DEN is a large-scale
dataset from the retinal field perspective, the number of training images is
still not enough for very deep models. Note that our proposed DEN dataset has
$265$ classes, including common and non-common retinal diseases or symptoms,
and only $8512$ training images, so it is not easy to achieve high Prec@1
accuracy for human doctors and AI machines. That is one of the reasons why we
investigate both Prec@1 and Prec@5. Also, reporting Prec@5 accuracy is more
appropriate from the real-world scenario perspective.

Table 3: This table shows the evaluation results of our keyword-driven and
non-keyword-driven clinical description generators (CDGs). Note that we
highlight the best scores of keyword-driven and non-keyword-driven generators
in each column, respectively. “w/o” denotes non-keyword-driven baseline
generators, and “w/” denotes our proposed keyword-driven generators. “BLEU-
avg” denotes the average score of BLEU-1, BLEU2, BLEU-3, and BLEU-4. Note that
the model based on “Jing, et al. [29]” has the best performance among all the
non-keyword-driven models, and the keyword-driven model based on “Jing, et al.
[29]” also has the best performance among all the models. All the keyword-
driven models, with the average feature fusion method, are superior to the
non-keyword-driven models. So, using keywords to reinforce the CDGs is
effective.

Model | BLEU-1 | BLEU-2 | BLEU-3 | BLEU-4 | BLEU-avg | CIDEr | ROUGE
---|---|---|---|---|---|---|---
Karpathy, et al. [30] | w/o | 0.081 | 0.031 | 0.009 | 0.004 | 0.031 | 0.117 | 0.134
w/ | 0.169 | 0.103 | 0.060 | 0.017 | 0.087 | 0.120 | 0.202
Vinyals, et al. [59] | w/o | 0.054 | 0.018 | 0.002 | 0.001 | 0.019 | 0.056 | 0.083
w/ | 0.144 | 0.092 | 0.052 | 0.021 | 0.077 | 0.296 | 0.197
Jing, et al. [29] | w/o | 0.130 | 0.083 | 0.044 | 0.012 | 0.067 | 0.167 | 0.149
w/ | 0.184 | 0.114 | 0.068 | 0.032 | 0.100 | 0.361 | 0.232
Li, et al. [36] | w/o | 0.111 | 0.060 | 0.026 | 0.006 | 0.051 | 0.066 | 0.129
w/ | 0.181 | 0.107 | 0.062 | 0.032 | 0.096 | 0.453 | 0.230

5.2 Clinical Description Generator (CDG) Verification

In [24, 27, 25], the authors mention that the evaluation of image description
generators is very subjective and there is no such thing as the most proper
metric to evaluate the text-to-text similarity. Different text-to-text
similarity metrics have different properties, so we exploit six commonly used
metrics, including BLEU-1, BLEU-2, BLEU-3, BLEU-4 [48], ROUGE [38], and CIDEr
[58], to evaluate generated results by our CDG. Table 3 contains the
evaluation results of our CDGs based on the above six different text-to-text
similarity metrics. All CDG modules with the keyword-driven method have better
performance than the non-keyword-driven CDGs, i.e., our baseline. It implies
that using keywords to reinforce the CDGs is effective. Based on Table 3 and
[29, 60], we find that the evaluation score of the medical image captioning,
based on the above commonly used text evaluation metrics, is much lower than
the evaluation score of the natural image captioning. One reason is that,
typically, the length of the medical image caption is much longer than the
natural image caption. Also, the medical image caption has more abstract words
or concepts than the natural image caption. These abstract words/concepts will
make algorithms difficult to generate correct captions. The other possible
reason is that the innate property of the commonly used text-to-text
similarity metrics [24, 27, 25] makes this happen. In addition, in Figure 8,
we show some generated clinical description results. Based on Figure 8, we
find that although our CDG module cannot always generate correct “age” or
“gender”, the models are capable of generating correct descriptions to
important characteristics for retinal images.

Based on the assumption mentioned in our Methodology section, subsection 5.1,
and subsection 5.2, we have shown the proposed AI-based method is
quantitatively effective.

5.3 Evaluation by DNN Visual Explanation Module

The main idea of DNN visual explanation module evaluation is that if our DNN
visual explanation results generated by CAM [70] are accepted by
ophthalmologists, it implies that the proposed method is qualitatively
effective. To prove the claim, we build the other retinal image dataset with
$300$ retinal images labeled by ophthalmologists and exploit the CNN
visualization tool, CAM, to visualize the learned feature and compare it to
the ground truth retinal image. We show the qualitative results in Figure 9.
In Figure 9, row (a) shows the four different kinds of raw images of retina
diseases and each raw image has a yellow sketch labeled by the ophthalmologist
to highlight the lesion areas on the retina. The numbers from (1) to (4)
denote the four different diseases, including Optic Neuritis, Macular
Dystrophy, Albinotic Spots in Macula, and Stargardt Cone-Rod Dystrophy,
respectively. We exploit CAM to generate row (b) to demonstrate the
visualization results of our DNN-based model. Then, row (c) is produced by the
same method as the row (b). Note that both row (b) and row (c) use the same
pre-trained weights of ImageNet but row (b) has fine-tuning on DEN dataset and
row (c) has no fine-tuning on DEN. The comparison of row (b) and row (c) shows
that the DNN-based model successfully learns the robust features of retinal
images by training on our DEN dataset. Also, row (b) indicates that the
features learned by DNN agree with the domain knowledge of ophthalmologists.
That is to say, the activation maps of our deep models are aligned with image
features that are clinically recognized by ophthalmologists as linked with the
identified disease. The above experimental results show our proposed AI-based
method is qualitatively effective.

Figure 9:  This figure shows the randomly selected qualitative results of CAM.
For the detailed explanation, please refer to subsection 5.3.

##  6 Conclusion

To sum up, we propose an AI-based method to automatically generate medical
reports for retinal images to improve the traditional retinal diseases
treatment procedure. The proposed method is composed of a DNN-based module,
including RDI and CDG sub-modules, and DNN visual explanation module. To train
our deep models and validate the effectiveness of our RDI and CDG, we propose
a large-scale retinal disease image dataset, DEN. Also, we provide another
retinal image dataset manually labeled by ophthalmologists to qualitatively
evaluate the proposed method. Our experimental results show the proposed AI-
based method is effective and successfully improves the conventional treatment
procedure of retinal diseases.

##  7 Acknowledgments

This work is supported by competitive research funding from King Abdullah
University of Science and Technology (KAUST) and University of Amsterdam.

## References

  * [1] Kedir M Adal, Peter G van Etten, Jose P Martinez, Lucas J van Vliet, and Koenraad A Vermeer.  Accuracy assessment of intra-and intervisit fundus image registration for diabetic retinopathy screening.  Investigative ophthalmology & visual science, 56(3):1805–1812, 2015.
  * [2] Bashir Al-Diri, Andrew Hunter, David Steel, Maged Habib, Taghread Hudaib, and Simon Berry.  A reference data set for retinal vessel profiles.  In 2008 30th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, pages 2262–2265. IEEE, 2008.
  * [3] Peter Anderson, Basura Fernando, Mark Johnson, and Stephen Gould.  Spice: Semantic propositional image caption evaluation.  In European Conference on Computer Vision, pages 382–398. Springer, 2016.
  * [4] Stanislaw Antol, Aishwarya Agrawal, Jiasen Lu, Margaret Mitchell, Dhruv Batra, C Lawrence Zitnick, and Devi Parikh.  Vqa: Visual question answering.  In Proceedings of the ICCV, pages 2425–2433, 2015.
  * [5] Ahmet Murat Bagci, Mahnaz Shahidi, Rashid Ansari, Michael Blair, Norman Paul Blair, and Ruth Zelkha.  Thickness profiles of retinal layers by optical coherence tomography image segmentation.  American journal of ophthalmology, 2008.
  * [6] Satanjeev Banerjee and Alon Lavie.  Meteor: An automatic metric for mt evaluation with improved correlation with human judgments.  In Proceedings of the acl workshop on intrinsic and extrinsic evaluation measures for machine translation and/or summarization, pages 65–72, 2005.
  * [7] Enrique J Carmona, Mariano Rincón, Julián García-Feijoó, and José M Martínez-de-la Casa.  Identification of the optic nerve head with genetic algorithms.  Artificial Intelligence in Medicine, 43(3):243–259, 2008.
  * [8] Aditya Chattopadhay, Anirban Sarkar, Prantik Howlader, and Vineeth N Balasubramanian.  Grad-cam++: Generalized gradient-based visual explanations for deep convolutional networks.  In 2018 IEEE Winter Conference on Applications of Computer Vision (WACV), pages 839–847. IEEE, 2018.
  * [9] Xinlei Chen, Hao Fang, Tsung-Yi Lin, Ramakrishna Vedantam, Saurabh Gupta, Piotr Dollár, and C Lawrence Zitnick.  Microsoft coco captions: Data collection and evaluation server.  arXiv preprint arXiv:1504.00325, 2015.
  * [10] Retinal Image Computing.  Understanding,“onhsd-optic nerve head segmentation dataset,” university of lincoln, united kingdom, 2004, 2012.
  * [11] Jeffrey De Fauw, Joseph R Ledsam, Bernardino Romera-Paredes, Stanislav Nikolov, Nenad Tomasev, Sam Blackwell, Harry Askham, Xavier Glorot, Brendan O’Donoghue, Daniel Visentin, et al.  Clinically applicable deep learning for diagnosis and referral in retinal disease.  Nature Medicine, page 1, 2018.
  * [12] Etienne Decencière, Guy Cazuguel, Xiwei Zhang, Guillaume Thibault, J-C Klein, Fernand Meyer, Beatriz Marcotegui, Gwénolé Quellec, Mathieu Lamard, Ronan Danno, et al.  Teleophta: Machine learning and image processing methods for teleophthalmology.  Irbm, 34(2):196–203, 2013.
  * [13] Etienne Decencière, Xiwei Zhang, Guy Cazuguel, Bruno Lay, Béatrice Cochener, Caroline Trone, Philippe Gain, Richard Ordonez, Pascale Massin, Ali Erginay, et al.  Feedback on a publicly distributed image database: the messidor database.  Image Analysis & Stereology, 33(3):231–234, 2014.
  * [14] Hao Fang, Saurabh Gupta, Forrest Iandola, Rupesh K Srivastava, Li Deng, Piotr Dollár, Jianfeng Gao, Xiaodong He, Margaret Mitchell, John C Platt, et al.  From captions to visual concepts and back.  In CVPR, pages 1473–1482, 2015.
  * [15] Muhammad Moazam Fraz, Paolo Remagnino, Andreas Hoppe, Bunyarit Uyyanonvara, Alicja R Rudnicka, Christopher G Owen, and Sarah A Barman.  An ensemble classification-based approach applied to retinal blood vessel segmentation.  IEEE Transactions on Biomedical Engineering, 59(9):2538–2548, 2012\.
  * [16] Lianli Gao, Kaixuan Fan, Jingkuan Song, Xianglong Liu, Xing Xu, and Heng Tao Shen.  Deliberate attention networks for image captioning.  AAAI, 2019.
  * [17] Ulrich Gerckens, Lutz Buellesfeld, Edward McNamara, and Eberhard Grube.  Optical coherence tomography (oct).  Herz, 28(6):496–500, 2003.
  * [18] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In CVPR, pages 770–778, 2016.
  * [19] Lisa Anne Hendricks, Zeynep Akata, Marcus Rohrbach, Jeff Donahue, Bernt Schiele, and Trevor Darrell.  Generating visual explanations.  In European Conference on Computer Vision, pages 3–19. Springer, 2016.
  * [20] Carlos Hernandez-Matas, Xenophon Zabulis, Areti Triantafyllou, Panagiota Anyfanti, Stella Douma, and Antonis A Argyros.  Fire: fundus image registration dataset.  Journal for Modeling in Ophthalmology, 1(4):16–28, 2017.
  * [21] Adam Hoover and Michael Goldbaum.  Locating the optic nerve in a retinal image using the fuzzy convergence of the blood vessels.  IEEE transactions on medical imaging, 22(8):951–958, 2003.
  * [22] Tao Hu, Pascal Mettes, Jia-Hong Huang, and Cees GM Snoek.  Silco: Show a few images, localize the common object.  In Proceedings of the IEEE International Conference on Computer Vision, pages 5067–5076, 2019.
  * [23] Jia-Hong Huang.  Robustness analysis of visual question answering models by basic questions.  King Abdullah University of Science and Technology, Master Thesis, 2017.
  * [24] Jia-Hong Huang, Modar Alfadly, and Bernard Ghanem.  Vqabq: Visual question answering by basic questions.  VQA Challenge Workshop, CVPR, 2017.
  * [25] Jia-Hong Huang, Modar Alfadly, and Bernard Ghanem.  Robustness analysis of visual qa models by basic questions.  VQA Challenge and Visual Dialog Workshop, CVPR, 2018.
  * [26] Jia-Hong Huang, Modar Alfadly, Bernard Ghanem, and Marcel Worring.  Assessing the robustness of visual question answering.  arXiv preprint arXiv:1912.01452, 2019.
  * [27] Jia-Hong Huang, Cuong Duc Dao, Modar Alfadly, and Bernard Ghanem.  A novel framework for robustness analysis of visual qa models.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 8449–8456, 2019.
  * [28] Jia-Hong Huang and Marcel Worring.  Query-controllable video summarization.  In Proceedings of the 2020 International Conference on Multimedia Retrieval, pages 242–250, 2020.
  * [29] Baoyu Jing, Pengtao Xie, Eric Xing, Baoyu Jing, Pengtao Xie, and Eric Xing.  On the automatic generation of medical imaging reports.  ACL, 2018.
  * [30] Andrej Karpathy and Li Fei-Fei.  Deep visual-semantic alignments for generating image descriptions.  In CVPR, pages 3128–3137, 2015.
  * [31] T Kauppi, V Kalesnykiene, et al.  Diaretdb0-standard diabetic retinopathy database, calibration level 0\. imageret project 2007.
  * [32] Tomi Kauppi, Valentina Kalesnykiene, Joni-Kristian Kamarainen, Lasse Lensu, Iiris Sorri, A Raninen, R Voutilainen, J Pietilä, H Kälviäinen, and H Uusitalo.  Diaretdb1—standard diabetic retinopathy database calibration level 1, 2007.
  * [33] Andrew Lang, Aaron Carass, Matthew Hauser, Elias S Sotirchos, Peter A Calabresi, Howard S Ying, and Jerry L Prince.  Retinal layer segmentation of macular oct images using boundary classification.  Biomedical optics express, 2013.
  * [34] Jonathan Laserson, Christine Dan Lantsman, Michal Cohen-Sfady, Itamar Tamir, Eli Goz, Chen Brestel, Shir Bar, Maya Atar, and Eldad Elnekave.  Textray: Mining clinical reports to gain a broad understanding of chest x-rays.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 553–561. Springer, 2018.
  * [35] Domenico Lepore, Graham E Quinn, Fernando Molle, Lorenzo Orazi, Antonio Baldascino, Marco H Ji, Maria Sammartino, Fabio Sbaraglia, Daniela Ricci, and Eugenio Mercuri.  Follow-up to age 4 years of treatment of type 1 retinopathy of prematurity intravitreal bevacizumab injection versus laser: fluorescein angiographic findings.  Ophthalmology, 125(2):218–226, 2018.
  * [36] Christy Y Li, Xiaodan Liang, Zhiting Hu, and Eric P Xing.  Knowledge-driven encode, retrieve, paraphrase for medical image report generation.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 6666–6673, 2019.
  * [37] Kunpeng Li, Ziyan Wu, Kuan-Chuan Peng, Jan Ernst, and Yun Fu.  Tell me where to look: Guided attention inference network.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 9215–9223, 2018.
  * [38] Chin-Yew Lin.  Rouge: A package for automatic evaluation of summaries.  Text Summarization Branches Out, 2004.
  * [39] Tsung-Yi Lin, Michael Maire, Serge Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C Lawrence Zitnick.  Microsoft coco: Common objects in context.  In European conference on computer vision, pages 740–755. Springer, 2014.
  * [40] Siqi Liu, Zhenhai Zhu, Ning Ye, Sergio Guadarrama, and Kevin Murphy.  Improved image captioning via policy gradient optimization of spider.  In Proceedings of the IEEE international conference on computer vision, pages 873–881, 2017.
  * [41] Yi-Chieh Liu, Yung-An Hsieh, Min-Hung Chen, C-H Huck Yang, Jesper Tegner, and Y-C James Tsai.  Interpretable self-attention temporal reasoning for driving behavior understanding.  In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 2338–2342. IEEE, 2020.
  * [42] Yi-Chieh Liu, Hao-Hsiang Yang, Chao-Han Huck Yang, Jia-Hong Huang, Meng Tian, Hiromasa Morikawa, Yi-Chang James Tsai, and Jesper Tegner.  Synthesizing new retinal symptom images by multiple generative models.  In Asian Conference on Computer Vision, pages 235–250. Springer, 2018.
  * [43] Meindert Niemeijer, Bram Van Ginneken, Michael J Cree, Atsushi Mizutani, Gwénolé Quellec, Clara I Sánchez, Bob Zhang, Roberto Hornero, Mathieu Lamard, Chisako Muramatsu, et al.  Retinopathy online challenge: automatic detection of microaneurysms in digital color fundus photographs.  IEEE transactions on medical imaging, 29(1):185–195, 2010.
  * [44] M Niemeijer, X Xu, A Dumitrescu, P Gupta, B van Ginneken, J Folk, and M Abramoff.  Inspire-avr: Iowa normative set for processing images of the retina-artery vein ratio, 2011.
  * [45] J Odstrčilík, Jiri Jan, J Gazárek, and R Kolář.  Improvement of vessel segmentation by matched filtering in colour retinal images.  In World Congress on Medical Physics and Biomedical Engineering, Munich, Germany. Springer, 2009.
  * [46] Marcos Ortega, Manuel G Penedo, José Rouco, Noelia Barreira, and María J Carreira.  Retinal verification using a feature points-based biometric pattern.  EURASIP Journal on Advances in Signal Processing, 2009:2, 2009.
  * [47] Nitika Pant Pai, Caroline Vadnais, Claudia Denkinger, Nora Engel, and Madhukar Pai.  Point-of-care testing for infectious diseases: diversity, complexity, and barriers in low-and middle-income countries.  PLoS medicine, 9(9):e1001306, 2012.
  * [48] Kishore Papineni, Salim Roukos, Todd Ward, and Wei-Jing Zhu.  Bleu: a method for automatic evaluation of machine translation.  In Proceedings of ACL, pages 311–318. Association for Computational Linguistics, 2002.
  * [49] Louis Pizzarello, Adenike Abiose, Timothy Ffytche, Rainaldo Duerksen, R Thulasiraj, Hugh Taylor, Hannah Faal, Gullapali Rao, Ivo Kocur, and Serge Resnikoff.  Vision 2020: The right to sight: a global initiative to eliminate avoidable blindness.  Archives of ophthalmology, 122(4):615–620, 2004.
  * [50] Prasanna Porwal, Samiksha Pachade, Ravi Kamble, Manesh Kokare, Girish Deshmukh, Vivek Sahasrabuddhe, and Fabrice Meriaudeau.  Indian diabetic retinopathy image dataset (idrid): a database for diabetic retinopathy screening research.  Data, 3(3):25, 2018.
  * [51] Ork De Rooij and Marcel Worring.  Efficient targeted search using a focus and context video browser.  ACM Transactions on Multimedia Computing, Communications, and Applications (TOMM), 8(4):51, 2012.
  * [52] Ramprasaath R Selvaraju, Michael Cogswell, Abhishek Das, Ramakrishna Vedantam, Devi Parikh, and Dhruv Batra.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In ICCV, pages 618–626, 2017.
  * [53] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv:1409.1556, 2014.
  * [54] Jayanthi Sivaswamy, SR Krishnadas, Gopal Datt Joshi, Madhulika Jain, and A Ujjwaft Syed Tabish.  Drishti-gs: Retinal image dataset for optic nerve head (onh) segmentation.  In 2014 IEEE 11th International Symposium on Biomedical Imaging (ISBI), pages 53–56. IEEE, 2014.
  * [55] Joes Staal, Michael D Abràmoff, Meindert Niemeijer, Max A Viergever, and Bram Van Ginneken.  Ridge-based vessel segmentation in color images of the retina.  TMI, 23(4):501–509, 2004.
  * [56] Melissa H Tukey and Renda Soylemez Wiener.  The impact of a medical procedure service on patient safety, procedure quality and resident training opportunities.  Journal of general internal medicine, 29(3):485–490, 2014.
  * [57] SG Vázquez, Brais Cancela, Noelia Barreira, Manuel G Penedo, M Rodríguez-Blanco, M Pena Seijo, G Coll de Tuero, Maria Antònia Barceló, and Marc Saez.  Improving retinal artery and vein classification by means of a minimal path approach.  Machine vision and applications, 24(5):919–930, 2013.
  * [58] Ramakrishna Vedantam, C Lawrence Zitnick, and Devi Parikh.  Cider: Consensus-based image description evaluation.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4566–4575, 2015.
  * [59] Oriol Vinyals, Alexander Toshev, Samy Bengio, and Dumitru Erhan.  Show and tell: A neural image caption generator.  In CVPR, pages 3156–3164, 2015.
  * [60] Xiaosong Wang, Yifan Peng, Le Lu, Zhiyong Lu, and Ronald M Summers.  Tienet: Text-image embedding network for common thorax disease classification and reporting in chest x-rays.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 9049–9058, 2018.
  * [61] Ryan T Yanagihara, Cecilia S Lee, Daniel Shu Wei Ting, and Aaron Y Lee.  Methodological challenges of deep learning in optical coherence tomography for retinal diseases: a review.  Translational Vision Science & Technology, 9(2):11–11, 2020.
  * [62] Chao-Han Huck Yang, Jia-Hong Huang, Fangyu Liu, Fang-Yi Chiu, Mengya Gao, Weifeng Lyu, I-Hung Lin, and Jesper Tegner.  A novel hybrid machine learning model for auto-classification of retinal diseases.  Workshop on Computational Biology, ICML, 2018.
  * [63] Chao-Han Huck Yang, Fangyu Liu, Jia-Hong Huang, Meng Tian, MD I-Hung Lin, Yi Chieh Liu, Hiromasa Morikawa, Hao-Hsiang Yang, and Jesper Tegner.  Auto-classification of retinal diseases in the limit of sparse data using a two-streams machine learning model.  In Asian Conference on Computer Vision, pages 323–338. Springer, 2018.
  * [64] Chao-Han Huck Yang, Yi-Chieh Liu, Pin-Yu Chen, Xiaoli Ma, and Yi-Chang James Tsai.  When causal intervention meets adversarial examples and image masking for deep neural networks.  In 2019 IEEE International Conference on Image Processing (ICIP), pages 3811–3815. IEEE, 2019.
  * [65] Hao-Hsiang Yang, Chao-Han Huck Yang, and Yi-Chang James Tsai.  Y-net: Multi-scale feature aggregation network with wavelet structure similarity loss function for single image dehazing.  In ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 2628–2632. IEEE, 2020.
  * [66] Hao-Hsiang Yang, Chao-Han Huck Yang, and Yu-Chiang Frank Wang.  Wavelet channel attention module with a fusion network for single image deraining.  arXiv preprint arXiv:2007.09163, 2020.
  * [67] Aliaa AA Youssif, Atef Z Ghalwash, Amr S Ghoneim, et al.  Comparative study of contrast enhancement and illumination equalization methods for retinal vasculature segmentation.  International Biomedical Engineering Conferemce, 2006.
  * [68] Jan Zahálka and Marcel Worring.  Towards interactive, intelligent, and integrated multimedia analytics.  In 2014 IEEE Conference on Visual Analytics Science and Technology (VAST), pages 3–12. IEEE, 2014.
  * [69] Matthew D Zeiler and Rob Fergus.  Visualizing and understanding convolutional networks.  In ECCV, pages 818–833. Springer, 2014.
  * [70] B. Zhou, A. Khosla, Lapedriza. A., A. Oliva, and A. Torralba.  Learning Deep Features for Discriminative Localization.  CVPR, 2016.


# Thread-Scalable Evaluation of Multi-Jet Observables

Walter T. Giele111email<EMAIL_ADDRESS>Gerben C. Stavenga222email:
<EMAIL_ADDRESS>and Jan Winter333email<EMAIL_ADDRESS>
Fermi National Accelerator Laboratory, Batavia, IL 60510, USA

###### Abstract:

A leading-order, leading-color parton-level event generator is developed for
use on a multi-threaded GPU. Speed-up factors between 150 and 300 are obtained
compared to an unoptimized CPU-based implementation of the event generator. In
this first paper we study the feasibility of a GPU-based event generator with
an emphasis on the constraints imposed by the hardware. Some studies of Monte
Carlo convergence and accuracy are presented for $PP\rightarrow 2,\ldots,10$
jet observables using of the order of $10^{11}$ events.

QCD, LO Computations, Jets, Hadronic Colliders

††preprint: FERMILAB-PUB-10-025-T

##  1 Introduction

Leading order (LO) parton-level Monte Carlos (MCs) play a prominent role in
collider phenomenology [1, 2, 3, 4, 5, 6]. As one needs to average the
calculation of the observable over many events, the evaluation time for the
event generation is a crucial issue in the development of LO parton level MCs.
Furthermore, to make full use of the recent progress in the calculation of
virtual corrections [7, 8, 9], fast tree-level event generators are needed for
the calculation of the radiative contributions in a next-to-leading order MC.

One can use large-scale grids for the generation of the tree-level events.
Such grids are expensive and need a large infrastructure. A more preferable
solution would be to run the MC on a single, affordable workstation. As we
will show this is possible using a massively parallel GPU. The NVIDIA Tesla
chip is designed for numerical applications [10] and the CUDA C compiler [11]
provides a familiar development environment. We will use the Tesla GPU
throughout the paper.444We thank the LQCD collaboration for giving us access
to the Tesla GPU processors.

In this paper we will execute all steps needed for event generation on the
GPU. These steps include the implementation of the unit-weight phase-space
generator Rambo [12], the evaluation of the strong coupling and parton density
function using LHAPDF [13], the evaluation of the leading-color $gg\rightarrow
2,\ldots,10$ gluon matrix elements at LO and the calculation of some
observables. The CPU is tasked with calculating the distributions using the
event weight and observables provided by the GPU. By utilizing memory with a
fast access time only, considerable speed-ups are obtained in the event
generation time. This memory is limited in size, requiring some coding effort.
As the GPU chips are developing fast, we can enhance the capabilities of our
parton-level generator in accordance.

In Refs. [14, 15] methods have been developed to evaluate multi-jet cross
sections on GPUs within the framework of the Helas matrix-element evaluator
[16], which forms the basis of the Madgraph event generator [1]. The method is
based on individual Feynman diagram evaluations. As such the scaling with the
number of external particles of the scattering process is faster than
factorial. Such an algorithm will have limited scalability properties, which
cannot be compensated by deploying a large number of threads. Instead, an
algorithm of polynomial complexity will have excellent scaling properties; its
only limitation is the available fast-access memory size. Polynomial
algorithms for the evaluation of ordered LO multi-parton matrix elements have
been formulated in the form of Berends–Giele (BG) recursion relations [17].
For a leading-color generator, any Standard Model matrix element can be
evaluated with an algorithm of polynomial complexity of degree 4 [18] or, by
using more memory storage, of degree 3 [19]. For any fixed color expansion,
the complexity remains polynomial. Therefore, we will use ordered recursive
evaluations of the matrix elements instead of Feynman diagram evaluations.

In this paper we present a GPU-based implementation of all basic tools needed
for a LO generator. In Section 2 we discuss the GPU and its hardware
limitations. According to these limitations, we will determine the optimal
running configuration as a function of the number of gluons. The algorithmic
implementation of the recursion algorithm and other tools such as phase-space
generation, experimental cuts and parton density functions are discussed in
Section 3. Finally, in Section 4 we put all pieces together and construct the
leading-color LO parton-level generator capable of generating up to $PP\to 10$
jets with sufficient statistics for serious phenomenology. The conclusions and
outlook are given in Section 5.

##  2 Thread-Scalable Algorithms for Event Generators

Monte Carlo algorithms belong to a class of algorithms, which can be trivially
parallelized, by dividing the events over the threads. Optimized for graphics
processing, the GPU works by having many threads executing essentially the
same instructions over different data. For a given class of events, e.g.
$n$-gluon scattering, the only difference between the events is due to the
external sources, i.e. the momentum and polarization four-vectors of each
gluon defining the state of the external gluon. The recursive algorithm acts
on these input sources in an identical manner. That is, each thread can
execute the same processor instructions to calculate the matrix-element
weight.

However, because of hardware constraints such a straightforward approach is
limited by the amount of available fast access memory. The GPU memory is
independent from the CPU memory and divided into the off-chip global memory
and the on-chip memory. This distinction is important as the off-chip memory
is large (of the order of gigabytes) but slow to access by the threads.
Therefore, we want to limit the access to the global memory by using it only
for the transfer of results to the CPU memory. The on-chip memory is fast to
access, but limited in size (of the order of tens of kilobytes). The first on-
chip memory structures are the registers. Each thread has its own registers,
which cannot be accessed by other threads. These registers are used within the
algorithm for variable storage, function evaluations, etc. The other on-chip
memory structure is shared memory, which is accessible to all the threads on a
multi-processor (MP). The current GPUs are not yet optimize-able to one event
per thread due to these shared memory and register constraints. With the next
generation of GPUs the shared memory will increase significantly, and we will
reach the point at which we can evaluate one event per thread up to large
multiplicities of gluons.

From this discussion the limitations are clear as each event requires a
certain amount of the limited register and shared memory. For the optimal
solution, we put the maximum number of events on one MP, such that the
evaluation does not exceed the available on-chip memory. The resulting
multiple threads per event can be used to “unroll” do-loops etc., thereby help
speed up the evaluation. This optimal solution is dependent on the rapidly
evolving hardware structure of the GPU chips.

By lowering the number of events per MP below the optimal solution, the number
of available threads per event increases. However, this will not lead to an
effective speed-up of the overall event generation as the total number of
threads per GPU is fixed. Once the number of events to be used per MP has been
determined, the GPU evaluation becomes scalable. The MC generator now simply
scales with the number of available MPs on the GPU.

$n$ | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12
---|---|---|---|---|---|---|---|---|---
events/MP | 102 | 68 | 48 | 36 | 28 | 22 | 18 | 15 | 13
threads/event | 10 | 15 | 21 | 28 | 36 | 45 | 55 | 66 | 78
Table 1: The number of $n$-gluon events, which can be simultaneously executed
on one MP (and is equal to $2048/[n\times(n+1)]$) and the number of available
threads per event (equal to $n\times(n+1)/2$). The total number of events
evaluated in parallel on the Tesla chip is 30$\times$(events/MP).

We will use the current NVIDIA Tesla chip for the numerical evaluations. This
chip consists of 30 MPs each capable of running up to 1024 threads. Each MP
has 16,384 32-bit registers and an internal shared memory of 16,384 bytes.
Each thread is assigned its own registers from the pool. Compiling the current
MC implementation indicates that 35 registers per thread are needed. This
gives us an upper maximum based solely on the use of registers of
$16,384/35=468$ threads per MP (each of which could potentially be used to
evaluate one event). The momenta and current storage is of more concern. As we
will see in the next section, for the evaluation of the $n$-gluon matrix
element, we need to store $n\times(n+1)/2$ four-vectors in single (float)
precision. This requires $8\times n\times(n+1)$ bytes of shared memory per
event on the MP.555A bit of calculus shows that when we need to store
$n\times(n+1)/2$ real-valued four-vectors in single precision we require
$4\times 4\times n\times(n+1)/2$ bytes of shared memory. The resulting maximum
number of events per MP as a function of the number of gluons is given in
Table 1. Note that up to 44-gluon scattering can be evaluated on the MP
(albeit with only one event per MP). Beyond 44 gluons the shared memory is too
small to store all the required four-vectors.

##  3 The Implementation of the Thread-Scalable Algorithm

Now that we have determined the optimal running configuration, i.e. the number
of events per MP, we can implement the algorithm. We will describe the
implementation of the Threaded EventS Simulator MC, which we name Tess MC, for
the NVIDIA Tesla chip.666The Tess MC code can be downloaded from the website:
http://vircol.fnal.gov/TESS.html. As we have many threads available per event,
we will use these threads to speed up the MC. In Figure 1 we show the thread
usage during different stages of the event generator for a $2\rightarrow n$
gluon process. The fraction of the evaluation time spent in each stage depends
on the gluon multiplicity. For 4 gluons, we get 20%, 20%, 50% and 0% for the
Rambo, PS-weight, ME-weight and the epilogue phases, respectively. For 12
gluons, the time consumption divides up into 2%, 18%, 75% and 4% for the four
phases.

Figure 1:  The thread usage for an $n$-gluon event as the algorithm progresses
through the stages of the event generation: flat phase-space generation,
phase-space weight evaluation (including parton density functions and
$\alpha_{S}$), matrix-element evaluation and finalization phase.

The initialization phase (not shown in Figure 1) consists of starting up the
kernel on the GPU. This is taken care of by the CUDA runtime code and does
essentially not depend on the number of threads it has to spawn. However it is
a significant part when the total kernel time is small, as for the 4-gluon
case.

The kernel starts initiating the unit-weight phase-space generator Rambo. On
the CPU this algorithm grows linearly with $n$ as we have to construct the
$n-2$ outgoing momenta. On the GPU we can employ $n-2$ threads to
simultaneously generate the outgoing momenta, making the Rambo code in
practice independent of $n$.777The Rambo algorithm has some summation
operations, which grow linearly with $n$, but this time scaling is very small
compared to the overall evaluation time of the Rambo algorithm.

After the momenta are generated, we have to calculate the strong coupling
constant, the parton density functions and the observables. We also determine,
if the event passes the canonical cuts. If the event fails the cuts, it is
only flagged as such; the matrix-element weight will still be evaluated as
this has no effect on the overall evaluation time. This means one can deviate
from the chosen canonical cuts on the CPU during the histogramming phase if so
desired. Note that we could in principle generate more events, which pass the
cut before starting the calculation of the matrix-element weights. This should
increase the performance of the Monte Carlo, at the cost of additional
bookkeeping.

The evaluation of the strong coupling constant and parton density functions
requires special attention. As we have used all shared memory for the four-
vector storage of the gluon currents and momenta, we have to use the off-chip
global memory to store the parton density and strong coupling constant
information in the form of grids. Furthermore, interpolation is required
between the grid points. To facilitate this, we use a special type of memory,
the so-called texture memory. This off-chip memory was designed for graphics
applications and performs hardware interpolations of the grid. Specifically,
we set up a 1-dimensional grid for the strong coupling constant. The value of
the strong coupling constant is stored as a function of the renormalization
scale at integer values of the grid. For the 2-dimensional grid used by the
parton density functions, the two dimensions are given by the factorization
scale and the parton fractions. This parton grid is directly obtained from
LHAPDF [13]. After the grid initialization, the texture memory can be accessed
by the GPU and its hardware will perform the appropriate linear interpolation
between the grid points when accessing the grid using non-integer values. This
way we have a very fast evaluation of the strong coupling and parton density
functions taking only about $6$% and $0.6$% of the total GPU time for 4-gluon
and 12-gluon processes, respectively.

The four-momenta are generated and the phase-space weight is determined, hence
we have to evaluate the matrix-element weight next. This happens at the core
of the event generator where we use recursion relations to compute these
weights. For this proof-of-concept program, we decided to use the recursion
relation of Ref. [17] and restricted ourselves to the case of pure gluonic
cross sections; quarks can be easily added at a later stage without changing
the event generator in a fundamental way. The recursion relations we employ
are given by

$J_{\mu}[m,\ldots,n]\;=\;\frac{1}{K[m,\ldots,n]^{2}}\;\Bigg{(}\,\sum_{i=m}^{n-1}\Big{[}J[m,\ldots,i],J[i+1,\ldots,n]\Big{]}_{\mu}\\\
+\;\sum_{i=m}^{n-2}\sum_{j=i+1}^{n-1}\Big{\\{}J[m,\ldots,i],J[i+1,\ldots,j],J[j+1,\ldots,n]\Big{\\}}_{\mu}\Bigg{)}\
,$ (1)

where $J_{\mu}[m,\ldots,n]$ is a conserved four-vector current depending on
the external gluons $\\{m,\ldots,n\\}$. Furthermore, we have used the
shorthand notations

$\begin{split}K^{\mu}[m,\ldots,n]\;&=\;\sum_{i=m}^{n}k_{i}^{\mu}\
,\\\\[8.53581pt]
\Big{[}J[\\{a\\}],J[\\{b\\}]\Big{]}_{\mu}\,&=\;2\,\Big{(}J[\\{a\\}]\cdot
K[\\{b\\}]\Big{)}\,J_{\mu}[\\{b\\}]-2\,\Big{(}K[\\{a\\}]\cdot
J[\\{b\\}]\Big{)}\,J_{\mu}[\\{a\\}]\\\ &+\Big{(}J[\\{a\\}]\cdot
J[\\{b\\}]\Big{)}\,\Big{(}K_{\mu}[\\{a\\}]-K_{\mu}[\\{b\\}]\Big{)}\
,\\\\[8.53581pt]
\Big{\\{}J[\\{a\\}],J[\\{b\\}],J[\\{c\\}]\Big{\\}}_{\mu}\,&=\;2\,\Big{(}J[\\{a\\}]\cdot
J[\\{c\\}]\Big{)}J_{\mu}[\\{b\\}]\\\ &-\Big{(}J[\\{a\\}]\cdot
J[\\{b\\}]\Big{)}\,J_{\mu}[\\{c\\}]-\Big{(}J[\\{c\\}]\cdot
J[\\{b\\}]\Big{)}\,J_{\mu}[\\{a\\}]\ ,\end{split}$ (2)

where the external gluon labeled $i$ has momentum $k_{i}^{\mu}$ and
polarization state $J^{\mu}[i]$. These four-vectors form the initial
conditions for the recursion relation. In addition to the $n$ momenta, the
recursion relation requires $n\times(n-1)/2$ four-vector currents to be stored
giving a total storage of $n\times(n+1)/2$ four-currents per event.

The recursion relations have a polynomial complexity of order $n^{4}$ for
calculating the currents [18]. By exploiting the available threads for each
event, we can reduce the algorithmic complexity of the BG recursion relation.
The relation is easily thread-able, which enables us to lower the polynomial
scaling of the evaluation time of the recursion relation to $n^{3}$. A full
recursion for an $n$-gluon process is completed in $n-1$ steps. In the first
step, we use $n-1$ threads to calculate the polarization vectors
$\\{J[2],\ldots,J[n]\\}$ needed as a starting point in the recursion relation.
We choose each polarization vector as a random unit vector orthogonal to the
respective gluon momentum. By doing this, instead of employing the
conventional helicity vectors, we obtain real-valued currents. This avoids
complex multiplications and reduces the shared-memory usage, resulting in a
significant time gain.

After the 1-currents have been determined, we use $n-2$ threads to calculate
the 2-currents $\\{J[2,3],J[3,4],\ldots,J[n-1,n]\\}$. We continue with the
$n-1$ steps until we have determined $J[2,3,\ldots,n]$ at which point we can
calculate the ordered amplitude and, hence, the matrix-element weight. Note
that because we make use of the multiple threads we have reduced the
computational effort from ${\cal O}(n^{4})$ to ${\cal O}(n^{3})$ complexity.

In principle we may be able to improve even further. The initial ${\cal
O}(n^{4})$ growth of the one-threaded recursion relation to calculate the
$J[2,3,\ldots,n]$ current can be reduced by rewriting the recursion relation
as

$J_{\mu}[m,\ldots,n]\;=\;\sum_{i=m}^{n-1}\left[\Big{(}W[i+1,\ldots,n]\cdot
J[m,\ldots,i]\Big{)}_{\mu}-\Big{(}W[m,\ldots,i]\cdot
J[i+1,\ldots,n]\Big{)}_{\mu}\right]\ ,$ (3)

where the tensor $W_{\mu\nu}$ is defined as

$\begin{split}W_{\mu\nu}[m,\ldots,n]\;&=\;2\,J_{\mu}[m,\ldots,n]\,K_{\nu}[m,\ldots,n]-K_{\mu}[m,\ldots,n]\,J_{\nu}[m,\ldots,n]\\\
&+\sum_{i=m}^{n-1}\Big{(}J_{\mu}[m,\ldots,i]\,J_{\nu}[i+1,\ldots,n]-J_{\mu}[i+1,\ldots,n]\,J_{\nu}[m,\ldots,i]\Big{)}\
.\end{split}$ (4)

By undoing the nested summations in the second term of Eq. 1 we have lowered
the complexity of the algorithm to ${\cal O}(n^{3})$. However, this is only
achieved at the cost of using significantly more storage. For each event, one
would have to store $n\times(n-1)/2$ $4\times 4$-tensors in addition to the
$n\times(n+1)/2$ momenta and current four-vectors. Up to $n\approx 10$ the
extra work of doing matrix multiplications together with the fact that the
relative pre-factor of the $n^{4}$-algorithm is small, $1/4$, compared to the
$n^{3}$-algorithm actually make the $n^{3}$-algorithm slower than the
$n^{4}$-algorithm. Moreover, the extra storage demand does not make the
$n^{3}$-algorithm attractive for our GPU implementation.

From the current for $n-1$ gluons we then obtain the amplitude for the
$n$-gluon matrix element by putting the off-shell leg on-shell, contracting in
with the final polarization vector and symmetrizing over the gluons in the
current. Specifically,

$\mathcal{A}(1,\ldots,n)\;=\;\sum_{\pi}^{(n-1)!}\mbox{Tr}\Big{[}T^{a_{\pi_{1}}}\ldots
T^{a_{\pi_{n-1}}}T^{a_{n}}\Big{]}\,m(\pi_{1},\ldots,\pi_{n-1},n)\ ,$ (5)

with

$m(\pi_{1},\ldots,\pi_{n-1},n)\;=\;\left.\Big{(}J[\pi_{1},\ldots,\pi_{n-1}]\cdot
J[n]\Big{)}\times K^{2}[1,\ldots,n-1]\right\rfloor_{K[1,\ldots,n]=0}\ .$ (6)

Notice that for a given phase-space point, we have to perform the permutation
sum requiring $(n-1)!$ steps to arrive at the full amplitude. This would
immediately lead to a factorial growth in the computer time. We can circumvent
the super-exponential sum over permutations in Eq. 5. In the leading-color
approximation this is easily accomplished and the color-summed, squared
amplitude is given by

$\left|\mathcal{A}(1,\ldots,n)\right|^{2}\;\sim\;N_{\mathrm{C}}^{n-2}\left(N_{\mathrm{C}}^{2}-1\right)\left(\sum_{\pi}^{(n-1)!}\left|m(\pi_{1},\ldots,\pi_{n-1},n)\right|^{2}+{\cal
O}\left(\frac{1}{N_{\mathrm{C}}^{2}}\right)\right)\ .$ (7)

As we will use this matrix element in a $2\rightarrow n-2$ gluon-scattering
phase-space integration, we can use the symmetry of the final state to remove
the permutation sum over the ordered amplitudes. In detail,

$\displaystyle d\,\sigma(PP\rightarrow n-2\ \mbox{jets})$ (8) $\displaystyle=$
$\displaystyle\int d\,x_{1}\ d\,x_{2}\
\frac{F_{g}(x_{1})F_{g}(x_{2})}{4\,p_{1}\cdot p_{2}}\frac{1}{(n-2)!}\int
d\,\Phi(p_{1}p_{2}\rightarrow p_{3}\cdots
p_{n})\sum_{\pi}^{(n-1)!}\left|m(\pi_{1},\ldots,\pi_{n-1},n)\right|^{2}$
$\displaystyle=$ $\displaystyle\int d\,x_{1}\ d\,x_{2}\
\frac{F_{g}(x_{1})F_{g}(x_{2})}{4\,p_{1}\cdot p_{2}}\,(n-1)\int
d\,\Phi(p_{1}p_{2}\rightarrow p_{3}\cdots
p_{n})\left|m(1,\sigma_{2},\ldots,\sigma_{n})\right|^{2}\ ,$

where $p_{1}=x_{1}\,P_{1}$, $p_{2}=x_{2}\,P_{2}$, the parton density function
is given by $F_{g}(x)$, $d\,\Phi$ is the phase-space integration measure and
$\\{\sigma_{2},\ldots,\sigma_{n}\\}$ is a permutation of the list
$\\{2,\ldots,n\\}$ assigned randomly for each MC phase-space point evaluation.

Eventually, in the very last step of our threaded event simulation all the
results are put together and returned to the CPU for processing.

By using the Tess MC, we can evaluate the differential $n$-jet cross sections
in the leading-color approximation. The algorithm is of polynomial complexity
and scales as $n^{3}$ with the number $n$ of gluons.

##  4 A Numerical Study of the Threaded Events Simulator

Figure 2: The horizontal axis is the number of events per MP in a sweep,
giving a total number of 30$\times$(events per MP) evaluated events per sweep.
The red curves used together with the vertical axes on the right indicate the
total GPU time in seconds for 1,000,000 sweeps. The blue curves depict the
evaluation time of one event in seconds as labeled by the vertical axes on the
left.

The first issue to study is the timing behavior of the Tess MC. We show our
results in Figure 2 where we plot the GPU timing as a function of events per
MP for several gluon multiplicities. Each MP will evaluate in a sweep a number
of events in parallel using the threads. In principle the sweep time should be
independent of the number of events evaluated by each MP as long as the
shared-memory constraints are not exceeded, cf. Table 1. However, we have to
execute a substantial amount of transcendental function calls per event, which
induces some queuing at the special function units each MP uses for evaluating
these functions. This queuing effect will increase as the number of events per
MP rises and, hence, lead to a slower execution of the sweep. In Figure 2, one
can see this complicated timing behavior, which is controlled by the GPU
hardware. As discussed the overall evaluation time increases with the number
of events per MP, see the red curves in the plots. In fact, the increase of
the overall evaluation time is overcome by the gain we achieve in evaluating
more events per MP. The more relevant quantity therefore is the evaluation
time per event, defined as the GPU evaluation time divided by the total number
of generated events. As clearly indicated by the blue curves in Figure 2, the
time consumption per event steadily decreases as the number of events per MP
increases. The best performance will be achieved by using the maximal number
of events available per MP.

$n$ |  $T_{n}^{\mathrm{GPU}}$ (seconds) | $P_{n}(3)$ |  $T_{n}^{\mathrm{CPU}}$ (seconds) | $P_{n}(4)$ | $G_{n}$
---|---|---|---|---|---
4 | $2.975\times 10^{-8}$ |  | $8.753\times 10^{-6}$ |  | 294
5 | $4.438\times 10^{-8}$ | 0.91 | $1.247\times 10^{-5}$ | 0.87 | 281
6 | $8.551\times 10^{-8}$ | 1.03 | $1.966\times 10^{-5}$ | 0.93 | 230
7 | $2.304\times 10^{-7}$ | 1.19 | $3.047\times 10^{-5}$ | 0.96 | 132
8 | $3.546\times 10^{-7}$ | 1.01 | $4.736\times 10^{-5}$ | 0.98 | 133
9 | $4.274\times 10^{-7}$ | 0.94 | $7.263\times 10^{-5}$ | 0.99 | 170
10 | $6.817\times 10^{-7}$ | 1.05 | $1.044\times 10^{-4}$ | 0.99 | 153
11 | $9.750\times 10^{-7}$ | 1.02 | $1.529\times 10^{-4}$ | 1.00 | 157
12 | $1.356\times 10^{-6}$ | 1.02 | $2.129\times 10^{-4}$ | 1.00 | 158
Table 2:  The GPU and CPU evaluation times per event, $T_{n}^{\mathrm{GPU}}$
and $T_{n}^{\mathrm{CPU}}$, given as a function of the number $n$ of gluons
for $gg\to(n-2)\,g$ processes. The polynomial scaling measures are also shown,
for the GPU, $P_{n}(3)$, and for the CPU, $P_{n}(4)$. The $P_{n}(m)$ are
defined as $P_{n}(m)=[(n-1)/n]\times\sqrt[m]{T_{n}/T_{n-1}\;}$. The rightmost
column finally displays the gain
$G=T_{n}^{\mathrm{CPU}}/T_{n}^{\mathrm{GPU}}$.

Now that we have determined the optimal running conditions, we give in Table 2
the evaluation time per event on the GPU compared to the evaluation time of
the same algorithm when executed on the CPU. As can be seen the speed-up in
evaluation time is substantial, ranging from almost a factor of 300 for
4-gluon processes to a factor of around 150 for 12-gluon processes. Note that
the speed-up is completely due to the fact that we evaluate in parallel 3060
and 390 events for the 4- and 12-gluon case, respectively. We have also tested
that running the events sequentially on the GPU using only one event per sweep
results in an event evaluation time, which is slower than the CPU evaluation
time. In particular, we found factors of 10 and 2 for the 4- and 12-gluon
computations, respectively. Because of the substantial time gains, a single
GPU can replace a large grid of hundreds of CPUs.

Also of interest is the scaling behavior of the algorithm. As expected, on the
CPU it is simply polynomial scaling with a factor of 4 in the limit of a large
number of gluons. We see from the table that this scaling is setting in
quickly. The GPU algorithm scales with a factor of 3 as discussed in Section
3. However, as the number of gluons increases, the number of events per MP
decreases. This makes the timing more dependent on specific hardware issues.
As can be seen from Table 2 the polynomial scaling is trending towards a
factor of 3.

Figure 3:  The number of sweeps versus several $gg\to(n-2)\,g$ cross sections
normalized to their respective best cross section estimates as given in Table
3. The error is the mean standard deviation.

Given the fast evaluation of events, we can easily generate ${\cal
O}(10^{11})$ events for the calculation of the LO cross sections. With these
large numbers of generated events, one has to carefully consider the
performance of the random number generators. In our case this should cause no
issues, since the number of generated random numbers is of the order of the
square root of the generator’s cycle length. However, as we average over
${\cal O}(10^{11})$ numbers, care has to be taken concerning the loss of
precision, which would result in a systematic underestimation of the cross
section. This is demonstrated in the first graph of Figure 3 where we have
used a single-precision summation to calculate the 4-gluon cross section. As
can be seen the effect becomes dramatic as the number of sweeps is rising and
we end up with a totally wrong determination of the cross section. We avoid
this problem by using the Kahan summation algorithm [20]. All other graphs of
the figure are produced by following this procedure. These additional graphs
display the convergence of the cross section estimate including its mean
standard deviation as a function of the number of GPU sweeps. The vertical
axis has been normalized to the respective best estimate of the cross section;
all of which are listed in Table 3. For this study, we have used Rambo as the
momenta generator, therefore, a severe under-sampling of small phase-space
regions with large weights may occur especially for larger gluon
multiplicities. Because the Rambo phase-space generation is flat and does not
reflect the scattering amplitudes’ strong dipole structure, such under-
sampling effects are expected and cause the peaking behavior of our cross
section estimates. Even with ${\cal O}(10^{10})$ phase-space points an
estimate of the 12-gluon cross section using the Rambo event generator is
quite unreliable and the mean standard deviation error estimate does not fully
reflect the true uncertainty. In a further development step, one may implement
a phase-space generator like Sarge [21], which is capable of adapting to the
QCD antenna structures as occurring in the matrix elements. As pointed out in
Ref. [6], this would resolve the phase-space integration issues we have seen
here.

$n$ |  $\sigma_{n}$ (pb) | $N_{\mathrm{generated}}$ | $N_{\mathrm{accepted}}$ |  $\sigma_{n}^{\textsc{Comix}}$ (pb)
---|---|---|---|---
4 | $(2.32421\pm 0.00047)\times 10^{8}$ | $3.06\times 10^{11}$ | $1.96848\times 10^{11}$ | $(2.3283\pm 0.0023)\times 10^{8}$
5 | $(1.4353\pm 0.0011)\times 10^{7}$ | $2.04\times 10^{11}$ | $1.12939\times 10^{11}$ | $(1.4355\pm 0.0014)\times 10^{7}$
6 | $(2.84780\pm 0.00096)\times 10^{6}$ | $1.44\times 10^{11}$ | $6.98918\times 10^{10}$ | $(2.8560\pm 0.0030)\times 10^{6}$
7 | $(6.356\pm 0.012)\times 10^{5}$ | $1.08\times 10^{11}$ | $4.49985\times 10^{10}$ | $(6.408\pm 0.015)\times 10^{5}$
8 | $(1.608\pm 0.011)\times 10^{5}$ | $8.40\times 10^{10}$ | $2.93316\times 10^{10}$ |
9 | $(4.38\pm 0.11)\times 10^{4}$ | $6.60\times 10^{10}$ | $1.88182\times 10^{10}$ |
10 | $(1.193\pm 0.024)\times 10^{4}$ | $5.40\times 10^{10}$ | $1.22356\times 10^{10}$ |
11 | $(3.550\pm 0.020)\times 10^{4}$ | $4.50\times 10^{10}$ | $7.88017\times 10^{9}$ |
12 | $(9.64361\pm 0.74)\times 10^{3}$ | $3.90\times 10^{10}$ | $5.13041\times 10^{9}$ |
Table 3:  The cross sections $\sigma_{n}$ for $gg\to(n-2)\,g$ and their mean
standard deviations in pb as calculated by the Tess MC using $10^{9}$ sweeps.
The two center columns show the total number of generated events and the
number of events passing the jet cuts. For comparison, the cross sections
$\sigma_{n}^{\textsc{Comix}}$ in pb as computed by Comix considering the full-
color dependence are also given.

Figure 4:  Left panels: the profile plots of the relative gauge invariance as
a function of the decimal logarithm of the matrix-element weight,
$\log_{10}W_{\mathrm{ME}}$; center panels: the normalized $H_{T}$
distributions and right panels: the normalized minimum $R$-separations between
pairs of jets. All of which is shown for $gg\to(n-2)\,g$ scatterings at a 14
TeV center-of-mass energy for $n=5,7,9,11$.

The convergence issues reflected in Figure 3 should be taken into account when
interpreting the uncertainties of our best cross section estimates, which are
listed in Table 3. For these cross section calculations of $gg\to(n-2)\,g$
scattering processes at a center-of-mass energy of 14 TeV, we have used the
CTEQ6L1 parton density function set [22] as implemented in LHAPDF [13] with a
fixed renormalization and factorization scale taken at $M_{Z}=91.188$ GeV. For
the jet cuts, we have chosen $p_{T}^{\mathrm{jet}}>20$ GeV,
$|\eta^{\mathrm{jet}}|<2.5$ and $\Delta R_{\textrm{jet--jet}}>0.4$. The cut
efficiencies for different numbers $n$ of gluons can be read off Table 3.
Employing these cuts we were also able to verify the jet production cross
sections that we have produced using Comix with the results reported in Ref.
[6]. To have a stringent comparison, we ran Comix for pure gluon scatterings
yielding cross sections that take the full-color dependence into account.
These results are also listed in the table; for the 4-gluon and 5-gluon
processes, they can be directly compared to the cross section estimates
obtained with the Tess MC on the GPU, since the leading-color approximation
already gives the exact result. The agreement is found to be satisfactory.

We show differential distributions in Figure 4. To obtain them we again used
$10^{9}$ sweeps where, for a certain gluon multiplicity, the total number of
generated events can be read off Table 3. We kept most of the input parameters
unaltered except for the jet cuts, which we changed to
$p_{T}^{\mathrm{jet}}>60$ GeV, $|\eta^{\mathrm{jet}}|<2.0$ and $\Delta
R_{\textrm{jet--jet}}>0.4$, and the choice of the renormalization and
factorization scales, which we decided to set dynamically using $H_{T}$ as a
scale. On the right hand side of Figure 4 we show for 3, 5, 7, 9 gluon jets in
the final state the normalized distributions for the $H_{T}$ observable and
the minimum $R$-separation, $R_{\mathrm{min}}$, which we define through the
jet–jet pair being closest in $R$-space, $R_{\mathrm{min}}=\min\\{\Delta
R_{ij}\\}$. As can be seen smooth distributions are easily obtained using the
Rambo phase-space generator. They are normalized to the total cross sections,
which have been calculated by Tess as $\sigma_{5}=(6.97838\pm 0.00044)\times
10^{4}$ pb, $\sigma_{7}=(4.9761\pm 0.0043)\times 10^{2}$ pb,
$\sigma_{9}=(4.532\pm 0.044)$ pb and $\sigma_{11}=(4.51\pm 0.19)\times
10^{-2}$ pb. On the left hand side of Figure 4 we have added profile plots
displaying the relative gauge invariance versus the decimal logarithm of the
matrix-element weight. Specifically, we show the average
$\left|K[1]\cdot\,J[2,\ldots,n]\right|^{2}/\left|J[1]\cdot\,J[2,\ldots,n]\right|^{2}$
and its mean standard deviation as a function of the matrix-element weight
$W_{\mathrm{ME}}=\left|m(1,2,\ldots,n)\right|^{2}=\left|(J[1]\cdot\,J[2,\ldots,n])\times
K^{2}[2,\ldots,n]\right|^{2}$. The behavior is as expected; for large weights,
we see gauge cancellations up to float precision. For small weights, the gauge
cancellations are less precise. However, these small-weight events are not
important since they do not contribute to the calculation of the observables.

Figure 5:  The ratio $(\sigma_{5}^{\textsc{Tess}}\times
d\sigma_{5}^{\textsc{Comix}}/dX)/(\sigma_{5}^{\textsc{Comix}}\times
d\sigma_{5}^{\textsc{Tess}}/dX)-1$ for the 5-gluon $X=H_{T}$ (left panel) and
$X=R_{\mathrm{min}}$ distributions, see the text. The mean standard deviation
error bars of the Comix calculation are also shown.

Finally, in Figure 5 we compare our results to the results obtained with the
Comix event generator [6]. For this comparison, we use both the $H_{T}$ and
$R_{\mathrm{min}}$ 5-gluon distributions and we fix the renormalization and
factorization scale through $M_{Z}=91.118$ GeV to avoid any issues resulting
from slight differences in the evolution codes for running scales between the
two MCs. Furthermore, to have a sole shape comparison, we plot the ratio
$(\sigma_{5}^{\textsc{Tess}}\times
d\sigma_{5}^{\textsc{Comix}}/dX)/(\sigma_{5}^{\textsc{Comix}}\times
d\sigma_{5}^{\textsc{Tess}}/dX)-1$ with the results shown in Figure 5 and $X$
being the observable in consideration. Note that for the minimum
$R$-separation distribution, we have excellent agreement with Comix. For the
$H_{T}$ distribution, we have to realize that the cross section spans 28
orders of magnitude. As Comix relies on importance sampling, it only sparsely
populates the tail of the distribution. This leads to large uncertainties at
large values of $H_{T}$ and, in these regions, Comix will hence tend to
underestimate the value for the cross section.

##  5 Conclusions and Outlook

In our first exploration of the potential of using multi-threaded GPU-based
workstations for Monte Carlo programs, we obtained very encouraging results.
We implemented the entire Tess Monte Carlo on the GPU chip; the only off-chip
usage occurs through utilizing the texture memory for the evaluation of the
parton density function and the strong coupling constant. The GPU global
memory is solely used for transferring the Monte Carlo results to the CPU
memory. At this exploratory phase of the project, we limited ourselves to the
calculation of leading-color leading-order $n$-gluon matrix elements. With
respect to the CPU-based implementation of our Monte Carlo we have found
impressive speed-ups in the computations reaching from ${\cal O}(300)$ for
$PP\to 2$ jets to ${\cal O}(150)$ for $PP\to 10$ jets.

Given these results we are encouraged to further develop the Tess Monte Carlo
by including quarks, vector bosons and subleading color contributions. We are
also planning to implement on the GPU a dipole-based phase-space generator
like Sarge as an alternative to the unit-weight phase-space generator Rambo.
This will avoid the under-sampling issues in high jet-multiplicity final
states. These improvements will result in a full leading-order parton-level
event generator, which will be at least two orders of magnitude faster than
existing leading-order parton-level generators.

More importantly, a GPU-based Monte Carlo can be used as the generator for the
real corrections in an automated next-to-leading order parton-level MC
generator. The virtual corrections can be calculated by using a generalized-
unitarity based method. These methods themselves spend about 90% of their
evaluation time calculating leading-order matrix elements. By using the GPU-
based matrix-element evaluator of the Tess Monte Carlo, one can safely
estimate a speed-up factor of ${\cal O}(10)$ in the evaluation time of the
virtual corrections. This means that the GPU-based automated NLO parton-level
generator can be successfully implemented on the GPU-based workstation and
obtain accurate results on a timescale of a day without resorting to a large-
scale PC farm.

Finally, GPU chips for numerical evaluations are still evolving rapidly. This
will lead to additional significant speed-ups over CPU-based Monte Carlos in
the coming years.

###### Acknowledgments.

We want to thank Jim Simone for suggesting the Tesla GPU chips for use in
event generators. We thank the High-Performance Computing Department at
Fermilab for giving us support and access to the LQCD Tesla-based
workstations. We would also like to thank Patrick Fox, who came up with the
name “Tess” and suggested it to us as the name for our Monte Carlo program.
Fermilab is operated by Fermi Research Alliance, LLC, under contract DE-
AC02-07CH11359 with the United States Department of Energy.

## References

  * [1] T. Stelzer and W. F. Long, Comput. Phys. Commun. 81 (1994) 357 [arXiv:hep-ph/9401258].
  * [2] F. Krauss, R. Kuhn and G. Soff, JHEP 0202 (2002) 044 [arXiv:hep-ph/0109036].
  * [3] P. D. Draggiotis, R. H. P. Kleiss and C. G. Papadopoulos, Eur. Phys. J. C 24 (2002) 447 [arXiv:hep-ph/0202201].
  * [4] M. L. Mangano, M. Moretti, F. Piccinini, R. Pittau and A. D. Polosa, JHEP 0307 (2003) 001 [arXiv:hep-ph/0206293].
  * [5] E. Boos et al. [CompHEP Collaboration], Nucl. Instrum. Meth. A 534 (2004) 250 [arXiv:hep-ph/0403113].
  * [6] T. Gleisberg and S. Höche, JHEP 0812 (2008) 039 [arXiv:0808.3674 [hep-ph]].
  * [7] W. T. Giele, Z. Kunszt and K. Melnikov, JHEP 0804, 049 (2008) [arXiv:0801.2237 [hep-ph]].
  * [8] C. F. Berger et al., Phys. Rev. D 78, 036003 (2008) [arXiv:0803.4180 [hep-ph]].
  * [9] A. van Hameren, C. G. Papadopoulos and R. Pittau, arXiv:0903.4665 [hep-ph].
  * [10] See the website “http://www.nvidia.com/object/product_tesla_c1060_us.html” for further information.
  * [11] See the website “http://www.nvidia.com/object/cuda_home.html” for further information.
  * [12] R. Kleiss, W. J. Stirling and S. D. Ellis, Comput. Phys. Commun. 40 (1986) 359.
  * [13] M. R. Whalley, D. Bourilkov and R. C. Group, arXiv:hep-ph/0508110.
  * [14] K. Hagiwara, J. Kanzaki, N. Okamura, D. Rainwater and T. Stelzer, arXiv:0908.4403 [physics.comp-ph].
  * [15] K. Hagiwara, J. Kanzaki, N. Okamura, D. Rainwater and T. Stelzer, arXiv:0909.5257 [hep-ph].
  * [16] H. Murayama, I. Watanabe and K. Hagiwara, “HELAS: HELicity Amplitude Subroutines for Feynman diagram evaluations”, KEK-91-11
  * [17] F. A. Berends and W. T. Giele, Nucl. Phys. B 306, 759 (1988).
  * [18] R. Kleiss and H. Kuijf, Nucl. Phys. B 312 (1989) 616.
  * [19] P. Draggiotis, R. H. P. Kleiss and C. G. Papadopoulos, Phys. Lett. B 439 (1998) 157 [arXiv:hep-ph/9807207].
  * [20] W. Kahan, ”Further remarks on reducing truncation errors”, Communications of the ACM 8 (1): 40 (1965)
  * [21] P. D. Draggiotis, A. van Hameren and R. Kleiss, Phys. Lett. B 483 (2000) 124 [arXiv:hep-ph/0004047].
  * [22] J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. M. Nadolsky and W. K. Tung, JHEP 0207 (2002) 012 [arXiv:hep-ph/0201195].