11institutetext: University of North Carolina at Chapel Hill, USA
22institutetext: Normandie Univ, UNIROUEN, UNIHAVRE, INSA Rouen, LITIS, 76000
Rouen, France

# Medical Image Synthesis with Context-Aware Generative Adversarial Networks

Dong Nie The authors contributed equally to this work11 Roger
Trullo11footnotemark: 1 22 Caroline Petitjean 22 Su Ruan 22 Dinggang Shen
Corresponding author11

###### Abstract

Computed tomography (CT) is critical for various clinical applications, e.g.,
radiotherapy treatment planning and also PET attenuation correction. However,
CT exposes radiation during acquisition, which may cause side effects to
patients. Compared to CT, magnetic resonance imaging (MRI) is much safer and
does not involve any radiations. Therefore, recently, researchers are greatly
motivated to estimate CT image from its corresponding MR image of the same
subject for the case of radiotherapy planning. In this paper, we propose a
data-driven approach to address this challenging problem. Specifically, we
train a fully convolutional network to generate CT given an MR image. To
better model the nonlinear relationship from MRI to CT and to produce more
realistic images, we propose to use the adversarial training strategy and an
image gradient difference loss function. We further apply AutoContext Model to
implement a context-aware generative adversarial network. Experimental results
show that our method is accurate and robust for predicting CT images from MRI
images, and also outperforms three state-of-the-art methods under comparison.

###### Keywords:

Generative Models, GAN, Image Synthesis, Deep Learning, AutoContext

##  1 Introduction

Computed Tomography (CT) imaging is widely used for both diagnostic and
therapeutic purposes in various clinical applications. In the cancer radiation
therapy, CT image provides Hounsfield units, which is essential for dose
calculation in treatment planning. Besides, CT image is also of great
importance for attenuation correction of positron emission tomography (PET) in
the popular PET-CT scanner [10].

However, patients are exposed to radiation during CT imaging, which can damage
normal body cells and further increase potential risks of cancer. Brenner et
al.brenner2007computed reported that 0.4% of cancers were due to CT scanning
performed in the past, and this rate will increase to as high as 1.5 to 2% in
the future. Therefore, the use of CT scan should be done with great caution.
Magnetic Resonance Imaging (MRI) on the other hand, is a safe imaging protocol
which also provides more anatomical details than CT for diagnostic purposes,
but unfortunately cannot be used for either dose calculation or attenuation
correction. To reduce unnecessary imaging dose for patients, it is clinically
desired to estimate CT images from MR images in many applications.

Figure 1: A pair of corresponding brain MR (left) and CT (right) images from
the same subject

It is technically difficult to directly estimate CT image from MR image. As
shown in Fig. 1, CT and MR images have very different appearances. MR images
contain much richer texture information than CT images. Therefore, it is
challenging to directly estimate a mapping from MRI to CT.

Recently, many researches focus on estimating one modal image from another
modality image, e.g., estimating CT image using MRI data. The first category
of methods is image segmentation based methods. Zaidi et al.[20] developed a
fuzzy clustering technique to segment MR images into different tissue classes,
and then refined segmentation manually. Berker et al.[1] proposed four class
tissue segmentation technique with a combined ultrashort-echo-time/Dixon MRI
sequence to fulfill MRI-based attenuation correction. These methods segment MR
images into different tissue classes, and then assign each class with a known
attenuation property. This category of methods highly depends on the
segmentation accuracy and always needs manual work to get final accurate
results. The second category of methods is atlas-based methods. These methods
first register an atlas (with the attenuation map) to the new subject’s MR
image, and then warp the corresponding attenuation map of the atlas to the new
MR image as its estimated attenuation map [3]. However, the performance of
these atlas-based methods highly depends on the registration accuracy. The
third category is learning-based methods, in which a non-linear mapping model
is learnt from MRI to CT image. Jog et al.[9] learned nonlinear regression
using random forest to improve MR Resolution. Tri et al.[7] presented an
approach to predict CT image from MRI using structured random forest. Since
estimating MRI from CT is not a one-to-one mapping task, using only the
intensity values cannot effectively distinguish the structural details. Such
methods often have to first represent the input MR image by features and then
map them to output the CT image. Thus, the performance of these methods is
bound to the quality of the extracted features and how well they can represent
the natural properties of the MR image.

On the other hand, recently the convolutional neural network (CNN) [11] became
popular in both computer vision and medical imaging fields. As a multi-layer
and fully trainable model, CNN is able to capture the complex non-linear
mapping from the input space to the output space. For the case of 2D images,
2D CNN has been widely used in many applications. However, it is unreasonable
to apply 2D CNN to process 3D medical images because 2D CNN considers the
image appearance slice by slice, thus potentially causing discontinuous
prediction results across slices. To address this issue, 3D CNNs have been
proposed. Ji et al.[8] presented a 3D CNN model for action recognition in an
uncontrolled environment. Tran et al.[16] used 3D CNN to effectively learn
spatio-temporal features on a large-scale video dataset. Dong et al.[4]
proposed a deep learning method for single image super-resolution. Li et
al.[13] applied deep learning models to estimate the missing PET image from
the MR image of the same subject.

Typically an L2 distance is used as loss function with the assumption that the
data is drawn from a Gaussian distribution. This can be problematic in the
case of multimodal distributions and tends to produce blurry results in the
output images [14]. Moreover, recently in the field of image generation,
Generative Adversarial Networks (GAN) have achieved state of the art results
producing very realistic images in an unsupervised setting [6], [15].

In this paper, we propose to learn the non-linear mapping from MR to CT images
through a 3D fully convolutional neural network (FCN), which is a variation of
the conventional CNN. Compared to CNN, FCN generates the structured output,
which can better preserve the neighborhood information in the predicted CT
image. This 3D FCN is used as generator in a Generative Adversarial framework
where an adversarial loss term in addition to the conventional reconstruction
error with the objective of producing more realistic CT data. The network is
trained on patches which makes its view to be restricted to the patch itself
and thus cannot provide long-range information. We use Auto-Context Model
(ACM) where each stage is trained using the GAN framework for making it
context-aware. To the best of our knowledge, this is the first application of
the GAN framework in the field of synthetic medical image generation.

The proposed method is evaluated on two real CT/MR datasets. Experimental
results demonstrate that our method can effectively predict CT image from MR
image, and also outperforms three state-of-the-art methods under the
comparison.

##  2 Methods

Deep learning models can learn a hierarchy of features, i.e., high-level
features built upon low-level features. CNN [5, 11] is one popular type of
deep learning models, in which trainable filters and local neighborhood
pooling operations are applied in an alternating sequence starting with the
raw input images. When trained with appropriate regularization, CNN can
achieve superior performance on visual object recognition and image
classification tasks [12]. However, most of CNNs are designed for 2D natural
images. They are not well suited for medical image analysis, since most of
medical images are 3D volumetric images, such as MRI, CT and PET.
Additionally, the output of a conventional CNN is a single target value, which
means that for a full volume, it would be necessary to compute one output for
every voxel. This framework has the obvious disadvantage of requiring a large
computational time, but also is unable to preserve neighborhood information in
the output space. The task of image synthesis can be seen as a regression
problem, where for every voxel in the input, an estimated output is required.
Typically an Euclidean loss function is used during training for regression
networks, which could potentially generate blurry images.

To address the above mentioned problems, we propose a generative adversarial
network in which fully convolutional networks form the generator. First, we
propose a basic 3D FCN structure to estimate the CT from MRI images. 3D
operations can better model the 3D spatial information and thus could solve
the discontinuity problem across slices, which can be present in 2D CNN.
Second, we utilize the adversarial training strategy [6] for the designed
network where an additional discriminator network can urge the generator’s
output to look like real CT as much as possible. We add an image gradient
difference term to the loss function of the generator, with the aim of
retaining the sharpness of the generated CT. Finally, we employ the Auto-
Context model to iteratively refine the output of the generator.

At testing stage, an input MR image is first partitioned into overlapping
patches. For each patch, the generator is used to predict the corresponding CT
patch. Finally, all predicted CT patches are merged into a single CT image by
averaging the intensities of overlapping CT regions.

In the following paragraphs, we will describe the framework that we use in the
MRI-to-CT prediction.

###  2.1 Proposed Generative Adversarial Networks

Generative Adversarial Networks (GANs) have been successfully applied in the
context of generative models. They were proposed in [6], and have achieved
impressive results in the context of image generation. GANs work by training
two different networks: a generator network $G$, and a discriminator network
$D$. $G$ is typically an FCN which generates images $G(\bm{z})$ from a random
noise vector $\bm{z}$, and $D$ is a CNN which estimates the probability that
an input image $x$ is drawn from the distribution of real images; that is, it
can classify an input image as real or synthetic. Both networks are trained
simultaneously with $D$ trying to correctly discriminate between real and
synthetic data, while $G$ is trying to produce realistic images that will
confuse $D$. More formally, for $D$ we would like to find its parameters:

$\max_{D}~{}\log(D(x))+\log(1-D(G(\bm{z})))$ (1)

For $G$ we would like to optimize

$\max_{G}~{}\log(D(G(\bm{z})))$ (2)

Inspired by the work in [14], where the authors use the GAN framework for
video frame prediction using the input frames as input instead of a random
vector, we design the GAN framework shown in Fig. 2. Our network includes a
generator which aims to estimate the CT and a discriminator which aims to
distinguish the real CT from the generated one. Specifically, we minimize the
binary cross entropy (bce) between the decisions of $D$ and the correct label
(real or synthetic), while the network $G$ is trying to minimize the binary
cross entropy between the decision done by $D$ and the wrong label for the
generated images, in addition to the traditional reconstruction error. In this
way, $D$ is trying to distinguish between real CT data, and the CT data
generated by $G$. At the same time, $G$ is trying to produce more realistic CT
images such that $D$ gets completely confused and cannot perform better than
chance.

Concretely, the loss function for $D$ is defined as:

$L_{D}=L_{bce}(D(x),1)+L_{bce}(D(G(x)),0)$ (3)

where $x$ is the input MRI image and

$L_{bce}(\hat{y},y)=-\frac{1}{N}\sum_{i=1}^{N}{y_{i}~{}log\hat{y_{i}}+(1-y_{i})~{}log(1-\hat{y_{i}})}$
(4)

$N$ is the number of samples in the minibatch, $y\in\\{0,1\\}$ represents the
label of the input data (0 for the generated and 1 for the real CT), and
$\hat{y}\in[0,1]$ is the estimated label by the discriminator network.

In the case of $G$, as mentioned above, we use a loss that includes an
adversarial term and a reconstruction error. It is defined as:

$L_{G}=\lambda_{1}L_{bce}(D(G(x)),1)+\lambda_{2}\lVert Y-G(x)\rVert_{2}^{2}$
(5)

where $Y$ is the corresponding ground-truth CT.

Similar to [14], in order to deal with the inherently blurry predictions
obtained from the L2 loss function (Eq.5), we propose to embed an image
gradient difference loss function in the generator training. The gradient
difference loss (gdl) between the generated CT and the real CT is given by:

$L_{gdl}=\left|{\left|{\nabla{Y_{x}}}\right|-\left|{\nabla{{\widehat{Y}}_{x}}}\right|}\right|^{2}+\left|{\left|{\nabla{Y_{y}}}\right|-\left|{\nabla{{\widehat{Y}}_{y}}}\right|}\right|^{2}+\left|{\left|{\nabla{Y_{z}}}\right|-\left|{\nabla{{\widehat{Y}}_{z}}}\right|}\right|^{2}$
(6)

where $Y$ is the ground truth image data, and $\hat{Y}$ is the estimated by
the network. This loss function tries to minimize the difference of the
magnitude of the gradient between the estimated image and the ground-truth CT.
In this way, the estimated data will try to keep zones with strong gradients
(for example edges) for effectively compensating for the L2 term. This is
approximated as finite differences during the implementation. Finally, the
total loss used for training the generator $G$ is defined as the weighted sum
of all the terms as shown in Eq.7.

$L\left({X,Y}\right)={L_{G}}\left({X,Y}\right)+{\lambda_{3}}{L_{gdl}}\left({X,Y}\right)$
(7)

The training is performed in an alternating fashion; first, $D$ is updated by
taking a minibatch of real CT data and a mini-batch of generated (the output
of $G$) data. Then, $G$ is updated by using another mini-batch of samples
including MRI and their corresponding CT.

Regarding the architecture, we make use of batch normalization, and also we
avoid the use of pooling layers since they reduce the spatial resolution of
feature maps. Although this property is desirable for some tasks, such as
image classification where the pooling over local neighborhoods can enhance
the invariance to certain image distortions, it is not desired in the task of
image prediction, where subtle image distortions need to be precisely captured
in the prediction process. These properties agree with several works that have
extended GAN, showing that some architectural constraints are necessary in
order to achieve a stable training. The authors in [15], for example, were
able to obtain very realistic images by using FCN without max pooling and with
batch normalization across different layers in both $G$ and $D$ in a way
similar to what we propose.

In Fig. 2, we also show the architecture of our generator network $G$ which
has the constraints mentioned above, where the numbers indicate the filter
sizes. The network takes as input an MRI image, and tries to generate the
corresponding CT image. It has 8 stages containing convolutions, Batch
Normalization and ReLu operations with number of filters 32, 32, 32, 64, 64,
64, 32, 32, respectively. The last layer only includes 32 convolutional
filters, and its output is considered as the estimated CT. Finally, the
Discriminator is a typical CNN architecture including three stages of
convolutions+Batch Normalization+ReLu+Max Pooling and a combination of
convolution with three fully connected layers, where the first two use ReLu as
activation function, and the last one uses sigmoid, whose output represents
the likelihood that the input data was drawn from the distribution of real CT.
The filter sizes are $5\times 5\times 5$, the numbers of filters are 32, 64,
128 and 256 for the convolutional layers, and the numbers of output nodes in
the fully connected layers are 512, 128 and 1.

Figure 2: Architecture used in the Generative Adversarial setting used for
estimation of synthetic images.

###  2.2 Auto-Context Model (ACM) for Refinement

Since our work is patch-based, the context information available for each
training sample is limited inside of the patch. This affects the modeling
capacity of our network. One way to enlarge the context during training is by
using the Auto-Context model which is commonly used in medical image analysis
applications and it has been shown to be very effective. It was first
introduced by [17], where it was used in the context of segmentation. The idea
is to train several classifiers iteratively where each classifier is trained
not only with the feature data, but also with the probability map obtained
from the previous classifier, which gives to the classifier additional context
information. At testing time, the features will be processed for each
classifier one after the other, concatenating the probability map to the input
features. This technique has shown to work remarkably well in the task of
semantic segmentation of MRI data [19], where Random Forest was used as
classifier. In this work, we show that the ACM can be successfully applied
also in regression tasks. Since the context features are extracted from the
whole previously estimated CT, they can encode information that is not
available using only patch based features. In our work, instead of building
several classifiers, we use Generator networks, and also instead of
concatenating probability maps, we concatenate the output of the previous
generator network. Specifically, we train iteratively several GANs that take
as input MRI patches and estimate corresponding CT patches. These patches are
concatenated as a second channel in the MRI patches, and this new data is used
as input during the training of the next GAN. An illustration of this scheme
is shown in Fig.3.

Figure 3: Proposed architecture for ACM with GAN.

##  3 Datasets

We use two datasets to test our proposed methods.

  * •

The brain dataset was acquired from 16 subjects with both MRI and CT scans in
the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database (see www.adni-
info.org for details). The MR images were acquired using a Siemens Triotim
scanner, with voxel size $1.2\times 1.2\times 1~{}\text{mm}^{3}$ , TE 2.95 ms,
TR 2300 ms, and flip angle $9^{\circ}$; The CT images, with voxel size
$0.59\times 0.59\times 3~{}\text{mm}^{3}$, were acquired on a Siemens Somatom
scanner. A typical example of preprocessed CT and MR images is given in Fig.
1.

  * •

Our pelvic dataset consists of 22 subjects, each with MR and CT images. The
spacings of CT and MR images are $1.172\times 1.172\times 1~{}\text{mm}^{3}$
and $1\times 1\times 1~{}\text{mm}^{3}$, respectively. In the training stage,
CT image is manually aligned to MR image to build the voxel-level
correspondence. After alignment, CT and MR images of the same patient have the
same image size and spacing. Since only pelvic regions are concerned, we
further crop the aligned CT and MR images to reduce the computational burden.
Finally, each preprocessed image has a size of $153\times 193\times 50$ and a
spacing of $1\times 1\times 1~{}\text{mm}^{3}$.

##  4 Experimental Results

We extracted randomly MRI patches of size $32\times 32\times 32$ with
corresponding CT of size $16\times 16\times 16$ using the same center point as
training samples. The networks were trained using the Adam optimizer with a
learning rate of $10^{-6}$, $\beta_{1}=0.5$ as suggested in [15], and mini-
batch size of 10. The generator was trained using
$\lambda_{1}=0.5,\lambda_{2}=\lambda_{3}=1$.

The code is implemented using the TensorFlow library, and it will be publicly
released upon acceptance. To demonstrate the advantage of the proposed method
in terms of prediction accuracy, we compare it with three widely used
approaches: an atlas-based method, a sparse representation based method, and
structured random forest with auto-context model. We used our own
implementation of the first two methods, while for the the structured random
forest we show the results declared in [7]. And all experiments are done in a
leave-one-out fashion.

  * •

Atlas-based method (Atlas): Here, the MR image of each atlas is first aligned
[18] onto the target MR image, and the resulting deformation field is used to
warp the CT image of each atlas. The final prediction is obtained by averaging
all warped CT images of all atlases.

  * •

Sparse representation based method (SR): After warping the atlases to the
target image space as described above, a local sparse representation is then
performed.

  * •

Structured random forest and auto-context model (SRF+) based method: Besides
the structured random forest, auto-context model (ACM) [17] is further used to
iteratively refine the prediction of CT images.

To compare the performance of different methods, we utilize mean absolute
error (MAE) and peak signal-to-noise ratio (PSNR) as measurements.

###  4.1 Impact of Proposed GAN model

To show the contribution of the proposed GAN model, we conduct comparison
experiments between the traditional FCN (just the generator shown in Fig. 2)
and the proposed GAN model. The PSNR are 24.7 and 25.9 for the traditional FCN
and the proposed approach, respectively. These results do not include the ACM.
We can visualize the results in Fig. 4, where the leftmost image is the input
MRI, and the rightmost image is the ground-truth CT. We can clearly see that
the generated data using the GAN approach has less artifacts than the
traditional FCN, by estimating an image that is closer to the desired output
not only quantitatively but also qualitatively.

Figure 4: Visual comparison for impact of adversarial training. FCN means
without adversarial training, and GAN means with adversarial training.

###  4.2 Impact of Auto-Context Refinement

To show the contribution of ACM, we present the performance of the proposed
method with iterations of ACM in Fig. 5. We can quantitatively observe that
both MAE and PSNR are improved gradually and consistently with iterations,
especially in the first two iterations. This is because ACM could solve the
short-range dependency by providing long-range information. Considering the
trade-off between benefit and training time, we choose 2 iterations of ACM in
our experiments for both datasets.

Figure 5: Performance of ACM on brain dataset with iterations.

###  4.3 Experimental Results for Brain Dataset

To qualitatively compare the estimated CT by different methods, we visualize
the generated CT with the ground truth in Fig. 6. We can see that the proposed
algorithm can better preserve the continuity, coalition and smoothness in the
prediction results, since it uses image gradient difference constraints in the
image patch as discussed in Section 2.1. Furthermore, the generated CT looks
closer to the real CT compared to all others, and we argue that this is due to
the adversarial training strategy which constrains the generated images to be
so similar to the real ones that a even a complex discriminator cannot perform
better than chance.

Figure 6: Visual comparison of original MR images, the estimated CT images by
our method and other methods, and the ground-truth CT images on 1 typical
subject on the brain dataset

We also quantitatively compare the predicted results in Table 1 in terms of
PSNR and MAE. Our proposed method outperforms all other methods in both
metrics, and it further demonstrates the advantage of our architecture.

Table 1: Average of MAE and PSNR on 16 subjects from the brain dataset by 4
different methods: Atlas, SR, SRF+, and Proposed

MAE | Mean(std.) | Med.
---|---|---
Atlas | 171.5(35.7) | 170.2
SR | 159.8(37.4) | 161.1
SRF+ | 99.9(14.2) | 97.6
Proposed | 92.5(13.9) | 92.1

PSNR | Mean(std.) | Med.
---|---|---
Atlas | 20.8(1.6) | 20.6
SR | 21.3(1.7) | 21.2
SRF+ | 26.3(1.4) | 26.3
Proposed | 27.6(1.3) | 27.6

###  4.4 Experimental Results for Pelvic Dataset

The prediction results by the same methods used above but on the pelvic
dataset are shown in Fig. 7. It can be clearly seen that our results are
consistent with the ground-truth CT. The quantitative results based on the
same metrics used in the previous dataset are shown in Table 2.

Figure 7: Visual comparison of original MR images, the estimated CT images by
our method and other methods, and the ground-truth CT images on 1 typical
subject on the pelvic dataset Table 2: Average of MAE and PSNR on 22 subjects
from the pelvic dataset by 4 different methods: Atlas, SR, SRF+, and Proposed

MAE | Mean(std.) | Med.
---|---|---
Atlas | 66.1(6.9) | 66.7
SR | 52.1(9.8) | 52.3
SRF+ | 48.1(4.6) | 48.3
Proposed | 39.0(4.6) | 39.1

PSNR | Mean(std.) | Med.
---|---|---
Atlas | 29.0(2.1) | 29.6
SR | 30.3(2.6) | 31.1
SRF+ | 32.1(0.9) | 31.8
Proposed | 34.1(1.0) | 34.1

Quantitative results in Tables 2 show that our method outperforms the other
methods in terms of both MAE and PSNR. Specifically, our method gives an
average PSNR of 34.1, which is considerably higher than the 32.1 obtained by
the state-of-the-art SRF+ method.

##  5 Conclusions

We have developed a 3D GAN model for estimating CT images from MRI images by
directly taking MR image patches as input and CT patches as output. The
performance is improved by the use of ACM since the context of the GAN is
effectively enlarged during the training process, which makes it context
aware. We have applied this model to predict CT images from their
corresponding MR images where the experiments demonstrate that our proposed
method can significantly outperform three state-of-the-art methods, showing
its suitability in regression tasks. Although we only considered the task of
CT image prediction, our proposed model can also be applied to other related
tasks involving a generative process in medical image analysis such as super-
resolution, image denoising and so on.

## References

  * [1] Yannick Berker et al.  Mri-based attenuation correction for hybrid pet/mri systems: a 4-class tissue segmentation technique using a combined ultrashort-echo-time/dixon mri sequence.  J. Nucl. Med., 53(5):796–804, 2012.
  * [2] David J Brenner and Eric J Hall.  Computed tomography¡ªan increasing source of radiation exposure.  N. Engl. J. Med., 357(22):2277–2284, 2007.
  * [3] Ciprian Catana et al.  Toward implementing an mri-based pet attenuation-correction method for neurologic studies on the mr-pet brain prototype.  J. Nucl. Med., 51(9):1431–1438, 2010.
  * [4] Chao Dong et al.  Learning a deep convolutional network for image super-resolution.  In Computer Vision–ECCV 2014, pages 184–199. Springer, 2014.
  * [5] Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In International conference on artificial intelligence and statistics, pages 249–256, 2010.
  * [6] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In Advances in Neural Information Processing Systems, pages 2672–2680, 2014.
  * [7] Tri Huynh et al.  Estimating CT Image from MRI Data Using Structured Random Forest and Auto-context Model.  2015\.
  * [8] Shuiwang Ji et al.  3d convolutional neural networks for human action recognition.  Pattern Analysis and Machine Intelligence, IEEE Transactions on, 35(1):221–231, 2013.
  * [9] Amod Jog et al.  Improving magnetic resonance resolution with supervised learning.  In Biomedical Imaging (ISBI), 2014 IEEE 11th International Symposium on, pages 987–990. IEEE, 2014.
  * [10] PE Kinahan et al.  Attenuation correction for a combined 3d pet/ct scanner.  Med. Phys., 25(10):2046–2053, 1998.
  * [11] Y. LeCun et al.  Deep learning.  Nature, 521(7553):436–444, 2015.
  * [12] Yann LeCun et al.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [13] Rongjian Li et al.  Deep learning based imaging data completion for improved brain disease diagnosis.  In MICCAI 2014, pages 305–312. Springer, 2014.
  * [14] Michael Mathieu, Camille Couprie, and Yann LeCun.  Deep multi-scale video prediction beyond mean square error.  arXiv preprint arXiv:1511.05440, 2015.
  * [15] Alec Radford, Luke Metz, and Soumith Chintala.  Unsupervised representation learning with deep convolutional generative adversarial networks.  arXiv preprint arXiv:1511.06434, 2015.
  * [16] Du Tran et al.  Learning spatiotemporal features with 3d convolutional networks.  arXiv preprint arXiv:1412.0767, 2014.
  * [17] Zhuowen Tu and Xiang Bai.  Auto-context and its application to high-level vision tasks and 3d brain image segmentation.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(10):1744–1757, 2010.
  * [18] Tom Vercauteren et al.  Diffeomorphic demons: Efficient non-parametric image registration.  NeuroImage, 45(1):S61–S72, 2009.
  * [19] Li Wang, Yaozong Gao, Feng Shi, Gang Li, John H Gilmore, Weili Lin, and Dinggang Shen.  Links: Learning-based multi-source integration framework for segmentation of infant brain images.  NeuroImage, 108:160–172, 2015.
  * [20] Habib Zaidi et al.  Magnetic resonance imaging-guided attenuation and scatter corrections in three-dimensional brain positron emission tomography.  Med. Phys., 30(5):937–948, 2003.


# Parameterized Task Graph Scheduling Algorithm for Comparing Algorithmic
Components††thanks: This work was supported in part by Army Research
Laboratory under Cooperative Agreement W911NF-17-2-0196. The authors
acknowledge the Center for Advanced Research Computing (CARC) at the
University of Southern California for providing computing resources that have
contributed to the research results reported within this publication. URL:
https://carc.usc.edu.

1st Jared Coleman  Thomas Lord Department of Computer Science
University of Southern California
Los Angeles, California, USA
<EMAIL_ADDRESS>2nd Ravi Vivek Agrawal  Thomas Lord Department of Computer
Science
University of Southern California
Los Angeles, California, USA
<EMAIL_ADDRESS>3rd Ebrahim Hirani  Thomas Lord Department of Computer
Science
University of Southern California
Los Angeles, California, USA
<EMAIL_ADDRESS>4th Bhaskar Krishnamachari  Ming Hsieh Department of
Electrical and Computer Engineering
University of Southern California
Los Angeles, California, USA
<EMAIL_ADDRESS>
###### Abstract

Scheduling distributed applications modeled as directed, acyclic task graphs
to run on heterogeneous compute networks is a fundamental (NP-Hard) problem in
distributed computing for which many heuristic algorithms have been proposed
over the past decades. Many of these algorithms fall under the list-scheduling
paradigm, whereby the algorithm first computes priorities for the tasks and
then schedules them greedily to the compute node that minimizes some cost
function. Thus, many algorithms differ from each other only in a few key
components (e.g., the way they prioritize tasks, their cost functions, where
the algorithms consider inserting tasks into a partially complete schedule,
etc.). In this paper, we propose a generalized list-scheduling algorithm that
allows mixing and matching different task prioritization and greedy node
selection schemes to produce 72 unique algorithms. We benchmark these
algorithms on four datasets to study the individual effect of different
algorithmic components on performance and runtime.

###### Index Terms:

scheduling, task graph, workflow, benchmarking

##  I Introduction

Task graph scheduling is a fundamental problem in distributed computing.
Essentially, the goal is to assign computational tasks to different compute
nodes in such a way that minimizes/maximizes some performance metric (e.g.,
total execution time, energy consumption, throughput, etc.). In this paper, we
will focus on the task scheduling problem concerning heterogeneous task graphs
and compute networks with the objective of minimizing makespan (total
execution time) under the related machines model111In the related machines
model, if the same task executes faster one some compute node $n_{1}$ than on
node $n_{2}$, then $n_{1}$ must execute all tasks faster than $n_{2}$ ($n_{1}$
is strictly faster than $n_{2}$). The related machines model as it pertains to
the task scheduling problem we study in this paper is described further in
Section I-A.. As this problem is NP-hard [1] and also not polynomial-time
approximable within a constant factor [2], many heuristic algorithms have been
proposed. One of the most popular paradigms for heuristic scheduling
algorithms is the list-scheduling paradigm. Essentially, list-scheduling
algorithms involve the following steps:

  1. 1.

Compute a priority for each task such that every task has a higher priority
than its dependents.

  2. 2.

Greedily schedule tasks in order of their computed priorities (from highest to
lowest) to run on the node that minimizes/maximizes some predefined cost
function.

Thus, many list scheduling algorithms differ only in a few algorithmic
components (e.g., their prioritization functions, their cost functions, where
the algorithms consider inserting tasks into an existing schedule, etc.). In
this paper, we propose a general parametric scheduler (extending SAGA [3], an
open-source python library for comparing task scheduling algorithms) that
allows us to mix-and-match different algorithmic components and evaluate how
they individually contribute to a list-scheduling algorithm’s performance and
runtime. Interestingly, we find that many new algorithms (composed of
previously unstudied combinations of algorithmic components) are pareto-
optimal with respect to performance and runtime. We also report the average
effects that both individual components and combinations of different
components have on performance and runtime, presenting evidence that the way
algorithmic components interact with each other is problem-dependent (i.e.,
depends on the task graph structure, whether or not the application is
communication or computation heavy, etc.).

###  I-A Problem Definition

Let us denote the task graph as $G=(T,D)$, where $T$ is the set of tasks and
$D$ contains the directed edges or dependencies between these tasks. An edge
$(t,t^{\prime})\in D$ implies that the output from task $t$ is required input
for task $t^{\prime}$. Thus, task $t^{\prime}$ cannot start executing until it
has received the output of task $t$. This is often referred to as a precedence
constraint. For a given task $t\in T$, its compute cost is represented by
$c(t)\in\mathbb{R}^{+}$ and the size of the data exchanged between two
dependent tasks, $(t,t^{\prime})\in D$, is $c(t,t^{\prime})\in\mathbb{R}^{+}$.
Let $N=(V,E)$ denote the compute node network, where $N$ is a complete
undirected graph. $V$ is the set of nodes and $E$ is the set of edges. The
compute speed of a node $v\in V$ is $s(v)\in\mathbb{R}^{+}$ and the
communication strength between nodes $(v,v^{\prime})\in E$ is
$s(v,v^{\prime})\in\mathbb{R}^{+}$. Under the related machines model [4], the
execution time of a task $t\in T$ on a node $v\in V$ is $\frac{c(t)}{s(v)}$,
and the data communication time between tasks $(t,t^{\prime})\in D$ from node
$v$ to node $v^{\prime}$ (i.e., $t$ executes on $v$ and $t^{\prime}$ executes
on $v^{\prime}$) is $\frac{c(t,t^{\prime})}{s(v,v^{\prime})}$.

The goal is to schedule the tasks on different compute nodes in such a way
that minimizes the makespan (total execution time) of the task graph. Let
$\mathcal{A}$ denote a task scheduling algorithm. For a given a problem
instance $(N,G)$ which represents a network/task graph pair, a schedule is a
set of tuples of the form $(t,v,r,e)$ where $t\in T$, $v\in V$ is the node on
which task $t$ is scheduled to run, $r\in\mathbb{R}^{+}$ is the start time,
and $e\in\mathbb{R}^{+}$ is the end time. A valid schedule $S$ must satisfy
the following properties

  * •

All tasks must be scheduled exactly once:

$\displaystyle\forall t\in T,\exists(t,v,r,e)\in S$
$\displaystyle\forall(t,v,r,e),(t^{\prime},v^{\prime},r^{\prime},e^{\prime})\in
S,t=t^{\prime}$ $\displaystyle\hskip
80.00012pt\Rightarrow(v,r,e)=(v^{\prime},r^{\prime},e^{\prime})$

  * •

All tasks must have valid start and end times:

$\forall(t,v,r,e)\in S,e-r=\frac{c(t)}{s(v)}$

  * •

Only one task can be scheduled on a node at a time (i.e., their start/end
times cannot overlap):

$\forall(t,v,r,e),(t^{\prime},v^{\prime},r^{\prime},e^{\prime})\in S,e\leq
r^{\prime}\lor e^{\prime}\leq r$

  * •

A task cannot start executing until all of its dependencies have finished
executing and their outputs have been received at the node on which the task
is scheduled:

$\displaystyle\forall(t,v,r,e),(t^{\prime},v^{\prime},r^{\prime},e^{\prime})\in
S,(t,t^{\prime})\in D$ $\displaystyle\hskip 80.00012pt\Rightarrow
e+\frac{c(t,t^{\prime})}{s(v,v^{\prime})}\leq r^{\prime}$

(a) Task Graph

(b) Network

(c) Schedule

Figure 1: Example problem instance and schedule.

Figure 1 depicts an example problem instance (task graph and network) and
solution (schedule). We define the makespan of a schedule $S$ as the time at
which the last task finishes executing:

$m(S)=\max_{(t,v,r,e)\in S}e$

Because the problem of minimizing makespan is NP-hard for this model [1], many
heuristic algorithms have been proposed. Traditionally, these heuristic
algorithms are evaluated on a set of problem instances and compared against
other algorithms based on their makespan ratio, which for a given problem
instance $(N,G)$ is the makespan of the schedule produced by the algorithm
divided by the minimum makespan of the schedules produced by the baseline
algorithms. Let $S_{\mathcal{A},N,G}$ denote the schedule produced by
algorithm $\mathcal{A}$ on problem instance $(N,G)$. Then the makespan ratio
of an algorithm $\mathcal{A}$ against a set of baseline algorithms
$\mathcal{A}_{1},\mathcal{A}_{2},\ldots$ for a problem instance $(N,G)$ can be
written

$\displaystyle\frac{m\left(S_{\mathcal{A},N,G}\right)}{\min\left\\{m\left(S_{\mathcal{A}_{1},N,G}\right),m\left(S_{\mathcal{A}_{2},N,G}\right),m\left(S_{\mathcal{A}_{3},N,G}\right),\ldots\right\\}}.$

Makespan ratio is a commonly used concept in the task scheduling literature
[5, 6]. It is common to measure the makespan ratios of an algorithm (against a
set of baseline algorithms) for a dataset of problem instances. System
designers use this technique, called benchmarking, to decide which
algorithm(s) best support(s) their application. In this paper, we study how
algorithmic components affect both the makespan and runtime (amount of time an
algorithm takes to produce a schedule). Let $r_{\mathcal{A}}(N,G)$ denote the
amount of time it takes to produce a schedule (on a particular system) using
algorithm $\mathcal{A}$. Then, the runtime ratio of an algorithm $\mathcal{A}$
against a set of baseline algorithms $\mathcal{A}_{1},\mathcal{A}_{2},\ldots$
for a problem instance $(N,G)$ can be written

$\displaystyle\frac{r_{\mathcal{A}}\left(N,G\right)}{\min\left\\{r_{\mathcal{A}_{1}}\left(N,G\right),r_{\mathcal{A}_{2}}\left(N,G\right),r_{\mathcal{A}_{3}}\left(N,G\right),\ldots\right\\}}.$

In reality, the actual runtime of an algorithm even on a particular instance
is not deterministic (due, for example, to background processes running on the
system). Thus, the runtime ratios reported in this paper should be interpreted
as estimated runtime ratios, rather than absolute deterministic values (unlike
the makespan ratios).

###  I-B Contributions

The main contributions of this work are:

  * •

Proposes an open-source generalized parametric list scheduling algorithm that
allows mixing and matching of different algorithmic components.

  * •

Reports benchmarking results for 72 algorithms produced by combining five
different algorithmic components on 20 publicly available datasets.

  * •

Reports both individual and combined effects of different algorithmic
components on average performance and runtime across all datasets.

  * •

Reports both individual and combined effects of different algorithmic
components on performance and runtime for each individual dataset.

The rest of this paper is organized as follows. In Section II we survey
related work in comparing task scheduling algorithms and algorithmic
components of such algorithms. Then, in Section III we present our methodology
for comparing different algorithmic components and the generalized parametric
list scheduling algorithm for doing so. In Section IV we report benchmarking
results for 72 algorithms created by combining five different types of
algorithmic components. We study the effects that individual components have
on performance and runtime and their interactions with other components.
Finally, in Section V we conclude the paper with a summary of the
contributions and directions for future work.

##  II Related Work

Existing approaches to comparing task scheduling algorithms mostly involve
benchmarking, whereby a set of algorithms is evaluated on one or more datasets
and various performance metrics are reported for comparison. One of the most
common metrics (and the one used in this paper) is the makespan ratio (also
called schedule length ratio) which is the makespan, or total execution time,
of the schedule produced by an algorithm for a given problem instance
normalized by the minimum makespan produced by any of the algorithms being
evaluated. Other metrics often reported include speedup (how many times larger
the makespan would be if the algorithm scheduled all tasks to run on a single
node), efficiency (average speedup per node), frequency that the algorithm is
the best algorithm among those being evaluated, and slack (a measurement of
schedule robustness) [7].

Among the many notable benchmarking efforts in the literature, one study [8]
benchmarks 15 scheduling algorithms and even discusses many of the algorithmic
components that we consider in this paper, including Assigning Priorities to
Nodes (tasks), Insertion vs. Non-Insertion, and Critical-Path-Based vs. non
Critical-Path-Based. They evaluate on different task graphs: Peer Set Graphs
and Random Graphs without studying the effects of these components on
performance. Another evaluation [9] evaluates eleven algorithms for
independent non-communicating tasks with uniformly random costs and
heterogeneous compute nodes with uniformly random speeds and yet another [7]
evaluates six algorithms on randomly generated application graphs and real-
world application graphs (FFT, Gaussian Elimination, Montage, and Epigenomics
Scientific Workflow). Metaheuristic algorithms like simulated annealing and
genetic algorithms have also been evaluated using a similar methodology [10].
A comprehensive evaluation of 31 list scheduling algorithms produced by
combining different algorithmic components has been conducted [6] (the same
paper also presents results for cluster scheduling which group tasks into
clusters before scheduling them), reporting useful benchmarking results but do
nothing regarding effect that individual algorithmic components have on
performance. Unfortunately, many of the existing comparison efforts are also
difficult to reproduce since the datasets typically not publicly available
and/or the scheduler implementations are not open-sourced.

To the best of our knowledge, this work is the first to evaluate the
performance of all possible combinations of a set of algorithmic components
and use the results to study the individual and combined effects these
components have on the performance (makespan ratio) and runtime of list-
scheduling algorithms.

##  III A Generalized List-Scheduling Algorithm

To study the effects of individual algorithmic components in list-scheduling
algorithms, we extended SAGA [3] (a python library for comparing task graph
scheduling algorithms) with a generalized parametric scheduling algorithm that
allows users to specify the following algorithmic components:

  * •

Priority Function: Function used to determine the order in which to schedule
tasks.

  * •

Comparison Function: Function used to determine which node to schedule a task
to.

  * •

Insertion-based vs. append-only scheduling: Insertion-based algorithms insert
tasks into the earliest sufficiently large window (on a particular node) in an
existing schedule. Append-only algorithms only consider scheduling tasks after
the last currently scheduled task finishes on a particular node.

  * •

Critical path reservation vs. no reservation: Whether or not the algorithm
commits to scheduling all tasks on the critical path — the longest (with
respect to node/edge weights) chain of tasks in the task graph — on the
fastest compute node 222This is consistent with the original definition [5] of
critical path and critical path reservation for the task scheduling problem we
study in this paper (fully heterogeneous, related machines model) but may not
be for other task scheduling varaints (e.g., for the unrelated machines
model).

  * •

Sufferage vs. no-sufferage consideration: Sufferage schedulers [11] consider
the top 2 highest-priority task in each iteration, choosing to schedule the
one for which scheduling on the second-best node would cause the greatest
detriment (with respect to the comparison function).

The algorithm works first by using the configured priority function to
determine the order in which to schedule tasks. The priority function takes
the problem instance (network/task graph pair) as input and returns a sequence
of tasks. In each iteration, the algorithm uses the comparison function and
whether or not the algorithm is using insertion-based or append-only
scheduling to decide which node it should schedule a task on. In insertion-
based scheduling, the scheduler finds, for each node, the earliest window of
time for which the node is idle and the window is large enough for the task to
execute (taking into account communication delays for dependency data
transfer). In the append-only scheme, the scheduling algorithm only considers
scheduling tasks after the last currently scheduled task on a node finishes
executing. Once the algorithm computes the candidate start/end times for the
task on each candidate node, it uses the comparison method to choose which
node to select. The comparison function may consider the candidate task’s
start time, end time, and/or other information to make this decision. If the
algorithm is a sufferage scheduler, then the comparison function is used
identify the best and the second-best node on which to schedule the task and
then used again to quantify the candidate task’s sufferage: the difference in
quality between scheduling on the second-best and best node. The algorithm
then chooses to schedule the task the highest sufferage value (the idea being
that not scheduling this task on its preferred node is more detrimental) and
returns the other task to the queue.

The priority and comparison functions are very general algorithmic components
for which there are many possible implementations. In this paper, we consider
the following implementations of each algorithmic component:

  * •

Priority Function (initial_priority)

    * –

UpwardRank: the priority function used in the HEFT scheduling algorithm [5]

    * –

CPoPRank: the priority function used in the CPoP scheduling algorithm [5]

    * –

ArbitraryTopological: an arbitrary topological sort of the task graph

  * •

Comparison Function (compare)

    * –

EFT: (Earliest Finish Time): schedules tasks to the node on which the task can
finish the soonest

    * –

EST (Earliest Start Time): schedules tasks to the node on which the task can
start the soonest

    * –

Quickest: schedules tasks to the node on which it executes in the least amount
of time

  * •

Insertion-based vs. append-only scheduling (append_only = True or False)

  * •

Critical path reservation vs. no reservation (critical_path = True or False)

  * •

Sufferage vs. no-sufferage consideration (sufferage = True or False)

Algorithms 1,2, and 3 implement the EFT, EST, and Quickest comparison
functions, respectively. Algorithms 4 and 5 implement functionality for
append-only and insertion-based scheduling. Implementations for the priority
functions are not included. Details for the UpwardRank and CPoPRank priority
functions can be found in the original paper they were proposed in [5].
Finally, Algorithm 6 represents the full parametric scheduling algorithm and
uses the functions defined in Algorithms 1-5.

Algorithm 1 Earliest Finish Time Compare Algorithm

1:function EFT($(start,end)$, $(start^{\prime},end^{\prime})$)

2:inputs: $(start,end)\in\mathbb{R}^{2}$ $\triangleright$ start/end time of
first candidate task

3:inputs: $(start^{\prime},end^{\prime})\in\mathbb{R}^{2}$ $\triangleright$
start/end time of second candidate task

4: return $end-end^{\prime}$ $\triangleright$ $<0$ iff the first task is
better (finishes earlier) than the second task

5:end function

Algorithm 2 Earliest Start Time Compare Algorithm

1:function EST($(start,end)$, $(start^{\prime},end^{\prime})$)

2:inputs: $(start,end)\in\mathbb{R}^{2}$ $\triangleright$ start/end time of
first candidate task

3:inputs: $(start^{\prime},end^{\prime})\in\mathbb{R}^{2}$ $\triangleright$
start/end time of second candidate task

4: return $start-start^{\prime}$ $\triangleright$ $<0$ iff the first task is
better (starts earlier) than the second task

5:end function

Algorithm 3 Quickest Execution Time Compare Algorithm

1:function EST($(start,end)$, $(start^{\prime},end^{\prime})$)

2:inputs: $(start,end)\in\mathbb{R}^{2}$ $\triangleright$ start/end time of
first candidate task

3:inputs: $(start^{\prime},end^{\prime})\in\mathbb{R}^{2}$ $\triangleright$
start/end time of second candidate task

4: return $(end-start)-(end^{\prime}-start^{\prime})$ $\triangleright$ $<0$
iff the first task is better (executes faster) than the second task

5:end function

Algorithm 4 Append-only available window finding algorithm

1:function GetOpenWindowAppendOnly($u$, $t$)

2: Let $est$ be the finish time of the last task scheduled on node $u$ ($0$ if
no tasks scheduled on $u$)

3: Let $dat$ be the minimum time at which all data from task $t$ dependencies
can be available on node $u$

4: $start\leftarrow\max\\{est,dat\\}$

5: $end\leftarrow start+\frac{c(t)}{s(u)}$

6: return $(start,end)$

7:end function

Algorithm 5 Insertion-based available window finding algorithm

1:

2:function GetOpenWindowInsertionBased($u$, $t$)

3: $T\leftarrow$ list of tasks scheduled on node $u$ in order of start time

4: if $T$ is empty then return $\textsc{GetOpenWindowAppendOnly}(()u,t)$

5: end if

6: Let $dat$ be the minimum time at which all data from task $t$ dependencies
can be available on node $u$

7: for task $t^{\prime}$ scheduled on node $u$ do $\triangleright$ iterate in
order of task start time

8: Let $est$ be the time that task $t^{\prime}$ finishes

9: $start\leftarrow\max\\{est,dat\\}$

10: $end\leftarrow start+\frac{c(t)}{s(u)}$

11: if $end\leq$ start time of next task in $T$ or $t$ is the last task in $T$
then return $(start,end)$

12: end for

13:end function

Algorithm 6 Generalized Parametric Scheduling Algorithm

1:parameters:
$\textsc{GetPriority}\in\\{\textsc{UpwardRanking},\textsc{CPoPRanking},\textsc{ArbitraryTopological}\\}$

2:parameters:
$\textsc{Compare}\in\\{\textsc{EFT},\textsc{EST},\textsc{Quickest}\\}$

3:parameters: $append\\_only\in\\{0,1\\}$

4:parameters: $sufferage\in\\{0,1\\}$

5:parameters: $critical\\_path\in\\{0,1\\}$

6:Initialize empty schedule $S$

7:Compute priority for each task in $G$ using Priority Function

8:Sort tasks in descending order of priority

9:if $append\\_only$ then

10: Let $\textsc{GetWindow}\leftarrow\textsc{GetOpenWindowAppendOnly}$

11:else

12: Let $\textsc{GetWindow}\leftarrow\textsc{GetOpenWindowInsertionBased}$

13:end if

14:while unscheduled tasks exist do

15: initialize $best$ and $second\\_best$ to arbitrary nodes

16: Let $t$ be the unscheduled task with highest priority

17: for $u\in V$ do

18: if
$\textsc{Compare}(\textsc{GetWindow}(t,u),\textsc{GetWindow}(t,best))<0$ then

19: $second\\_best\leftarrow best$

20: $best\leftarrow u$

21: else if
$\textsc{Compare}(\textsc{GetWindow}(t,u),\textsc{GetWindow}(t,second\\_best))<0$
then

22: $second\\_best\leftarrow u$

23: end if

24: end for

25: if $sufferage$ and there are at least two unscheduled tasks then

26: initialize $best^{\prime}$ and $second\\_best^{\prime}$ to arbitrary nodes

27: Let $t^{\prime}$ be the unscheduled task with second-highest priority

28: for $u\in V$ do

29: if
$\textsc{Compare}(\textsc{GetWindow}(t^{\prime},u),\textsc{GetWindow}(t^{\prime},best^{\prime}))<0$
then

30: $second\\_best^{\prime}\leftarrow best^{\prime}$

31: $best^{\prime}\leftarrow u$

32: else if
$\textsc{Compare}(\textsc{GetWindow}(t^{\prime},u),\textsc{GetWindow}(t^{\prime},second\\_best^{\prime}))<0$
then

33: $second\\_best^{\prime}\leftarrow u$

34: end if

35: end for

36: $sufferage\\_t\leftarrow$Compare($(t,second\\_best),(t,best)$)

37:
$sufferage\\_t^{\prime}\leftarrow$Compare($(t^{\prime},second\\_best^{\prime}),(t^{\prime},best^{\prime})$)

38: if $sufferage\\_t^{\prime}>sufferage\\_t$ then

39: $t,best\leftarrow t^{\prime},best^{\prime}$

40: end if

41: end if

42: Add $(t,start,end,best)$ to $S$

43:end while

44:return $S$

The combinations of these algorithmic components allow for $72$ unique
scheduling algorithms. We evaluated each of these on twenty datasets based on
four task graph structures (chains, in_trees, and out_trees, cycles) and five
communication-to-computation ratios, or CCRs, (CCR $1/5$, $1/2$, $1$, $2$,
$5$) for each. Each in_trees, out_trees, and chains datasets consists of $100$
randomly generated network/task graph pairs following a common methodology
used in the literature [12]. In-trees and out-trees are generated with between
$2$ and $4$ levels (chosen uniformly at random), a branching factor of either
$2$ or $3$ (chosen uniformly at random), and node/edge-weights drawn from a
clipped Gaussian distribution (mean: $1$, standard deviation: $1/3$, min: $0$,
max: $2$). Parallel chains task graphs are generated with between $2$ and $5$
parallel chains (chosen uniformly at random) of length between $2$ and $5$
(chosen uniformly at random) and node/edge-weights drawn from the same clipped
Gaussian distribution. Randomly weighted networks are complete graphs with
between $3$ and $5$ nodes (chosen uniformly at random) and node/edge-weights
drawn from the same clipped Gaussian distribution. Finally, the network
communication strengths are scaled to achieve the desired CCR of either $1/5$,
$1/2$, $1$, $2$, or $5$. The cycles dataset is based on the cycles scientific
workflow [13] and represents a multi-crop, multi-year agro-ecosystem model for
simulating crop production. For each workflow, the runtime of each task,
input/output sizes in bytes, and speedup factor (compute speed) for each
machine are available from public execution trace
information333https://github.com/wfcommons/pegasus-instances,
https://github.com/wfcommons/makeflow-instances. We set communication
strengths to be homogeneous so that the average CCR is
$\frac{1}{5},\frac{1}{2},1,2,$ or $5$ (resulting in five different datasets).
Figure 2 shows an example in_trees, out_trees, chains, and cycles task graph.

(a) Example task graph structure for the in-trees datasets.

(b) Example task graph structure for the out-trees datasets.

(c) Example task graph structure for the chains datasets.

(d) Example task graph structure for the cycles datasets.

Figure 2: Example task graphs for each of the datasets evaluated.

##  IV Results

In this section we discuss the results of running each of the $72$ schedulers
on $20$ datasets ($4$ dataset types and $5$ CCRs). Each dataset consists of
$100$ problem instances. All schedulers that are pareto optimal for at least
one of the evaluated datasets are listed in Table I. There are $24$ such
schedulers ($48$ of the $72$ total schedulers evaluated were strictly
dominated by at least one other scheduler for every dataset).

component | initial_priority | append_only | compare | critical_path | sufferage
---|---|---|---|---|---
scheduler |  |  |  |  |
EFT_Ins_AT | ArbitraryTopological | False | EFT | False | False
EFT_Ins_AT_Suf | ArbitraryTopological | False | EFT | False | True
EST_Ins_AT | ArbitraryTopological | False | EST | False | False
EST_Ins_AT_Suf | ArbitraryTopological | False | EST | False | True
EST_Ins_CP_AT | ArbitraryTopological | False | EST | True | False
MCT [9] | ArbitraryTopological | True | EFT | False | False
Sufferage [11] | ArbitraryTopological | True | EFT | False | True
EFT_App_CP_AT | ArbitraryTopological | True | EFT | True | False
EST_App_AT | ArbitraryTopological | True | EST | False | False
EST_App_AT_Suf | ArbitraryTopological | True | EST | False | True
MET [9] | ArbitraryTopological | True | Quickest | False | False
EFT_Ins_CR | CPoPRanking | False | EFT | False | False
EFT_Ins_CR_Suf | CPoPRanking | False | EFT | False | True
EST_Ins_CR | CPoPRanking | False | EST | False | False
EST_Ins_CR_Suf | CPoPRanking | False | EST | False | True
EFT_App_CR_Suf | CPoPRanking | True | EFT | False | True
HEFT [5] | UpwardRanking | False | EFT | False | False
EFT_Ins_UR_Suf | UpwardRanking | False | EFT | False | True
EST_Ins_UR | UpwardRanking | False | EST | False | False
EST_Ins_UR_Suf | UpwardRanking | False | EST | False | True
EFT_App_UR | UpwardRanking | True | EFT | False | False
EFT_App_UR_Suf | UpwardRanking | True | EFT | False | True
EST_App_UR | UpwardRanking | True | EST | False | False
EST_App_UR_Suf | UpwardRanking | True | EST | False | True
TABLE I: All scheduling algorithms that were pareto-optimal for at least one
evaluated dataset ($24$ of $72$ schedulers). All other schedulers ($48$ of
$72$) were strictly dominated by another algorithm in every dataset.

A pareto-optimal scheduler is one such that no other scheduler has both a
lower average makespan ratio and a lower average runtime ratio on a given
dataset. Figure 3(a) depicts the pareto-optimal schedulers for each dataset.
Each subplot has $24$ markers (one for each pareto-optimal scheduler) but only
the schedulers that are pareto optimal for each dataset are colored blue. For
example, for the dataset in_trees_ccr_0.2 (in_trees with CCR=1/5), two
schedulers that are pareto-optimal for another dataset are strictly dominated
by another scheduler for this dataset (and so are colored red). Figure 3(b)
serves as a kind of legend for the scatter plot. Each cell represents the
pareto-optimal scheduler’s rank compared to other pareto-optimal schedulers
for the same dataset with respect to its runtime ratio. For each dataset, the
scheduler labeled $1$ is that with the least runtime ratio (furthest blue
marker to the left in the corresponding scatter plot in Figure 3(a)). Because
it is pareto-optimal, this also indicates it has the highest makespan ratio of
pareto-optimal schedulers for the dataset (and so is also the highest blue
marker in the corresponding scatter plot in Figure 3(a)). Blank cells indicate
that the scheduling algorithm is not pareto-optimal for the dataset. Thus,
schedulers that have consistently low order numbers (like Sufferage) are
pareto-optimal mostly (or entirely) due to their low runtimes and not due to
their low makespan ratios. Results indicate that these are scheduling
algorithms that are fast, but not very performant. On the other hand,
scheduling algorithms that have consistently high order numbers (like HEFT)
are pareto-optimal mostly (or entirely) due to their low makespan ratios
despite higher runtimes. These algorithms are performant, but slower in
generating schedules.

(a) Markers represent scheduling algorithms that are pareto-optimal for at
least one of the evaluated datasets.

(b) The chart number represents the pareto-optimal schedule’s position (from
left to right) for each dataset.

Figure 3: Pareto-optimal scheduling algorithms for each dataset.

###  IV-A Effects of Algorithmic Components

In this section, we study the effect that different algorithmic components
have on performance and runtime. Figure 4 suggests that, across all datasets,
the priority function has a small effect on makespan ratio with UpwardRanking
just slightly out-performing ArbitraryTopological and CPoPRanking. Similar
results for the append-only and sufferage components are shown in Figures 6
and 8. Figure 5 indicates that while the Quickest comparison function is
clearly the least performant, EFT and EST have roughly similar performance
across all datasets. Figure 7 shows that critical path reservation tends to
increase makespan ratios and, more dramatically, increase runtime ratios,
suggesting that critical path reservation, in general, is a poor strategy.
Recall, however, that Table I from the previous section though, does indicate
that there are datasets for which critical path reserving schedulers are
pareto-optimal. This suggests that the effects of these components might be
dataset-specific.

Figure 4: Effect of the initial priority function on the makespan and runtime
ratios over all datasets. UR stands for UpwardRanking, AT for
ArbitraryTopological, and CR for CPoPRanking.

Figure 5: Effect of the comparison function on the makespan and runtime ratios
over all datasets.

Figure 6: Effect of the insertion vs. append-only scheme on the makespan and
runtime ratios over all datasets.

Figure 7: Effect of the critical path reservation on the makespan and runtime
ratios over all datasets.

Figure 8: Effect of the sufferage selection scheme on the makespan and runtime
ratios over all datasets.

We do, in fact, observe in some cases quite different behavior for individual
datasets. For example, Figure 9 depicts the effect of the comparison function
on the makespan and runtime ratios for the cycles_ccr_5 dataset (cycles task
graphs with CCR=5). The Quickest comparison function, which performs generally
terribly compared to EFT and EST for other datasets, outperforms EFT and EST
by a large margin! The dataset isn’t the only thing that interacts with the
effects of algorithmic components, though. Figure 10 depicts a few of the more
interesting interactions between algorithmic components (averaged over all
datasets). Figure 10(a) depicts an interaction between the append_only and
initial_priority parameters, suggesting that the combination of
append_only=True and intial_priority=CPoPRanking has a more detrimental effect
on the makespan ratio that either of the parameters do by themselves. In other
words, the append-only strategy is particularly bad when using the CPoPRanking
priority function. Figure 10(b) shows that the Quickest comparison function is
generally bad but less so on communication-heavy applications (those with high
CCR). Figure 10(c) shows that it is also a particularly bad comparison
function for the out_trees datasets. Figure 10(d) suggests that the small
difference in the makespan ratio for schedulers that do critical path
reservation may be due almost entirely to how critical path reservation
increases makespan ratios for in_trees datasets. Plots for the individual
effects and interactions between algorithmic components can be found in the
appendix.

Figure 9: Effect of the comparison function on the makespan and runtime ratios
for the cycles_ccr_5 dataset (cycles task graphs with CCR=5).

(a) Interaction between append_only and priority parameters.

(b) Interaction between the compare paremeter and the task graph’s CCRs.

(c) Interaction between compare function parameter and the dataset type.

(d) Interaction between the critical_path parameter and the dataset type.

Figure 10: Interactions Between Algorithmic Components

##  V Conclusion

In this paper, we proposed a generalized parametric list scheduling approach
for studying the individual and combined effects of different algorithmic
components. We evaluated 72 algorithms produced by combining five different
algorithmic components on 20 datasets and presented results on the individual
and combined effects of these components on average performance and runtime
across all datasets. We also discuss how these results differ for individual
datasets, suggesting that the way algorithmic components interact with each
other is problem-dependent (i.e., depends on the task graph structure, whether
or not the application is communication or computation heavy, etc.). This work
opens many avenues for future research. First, this work can be extended by
considering new algorithmic components (e.g., k-depth lookahead), new
implementations for the five current algorithmic components, and other
datasets. In particular, it would be interesting to see more results for
application-specific datasets like cycles.

In this paper, we compared algorithms using the traditional benchmarking
approach, whereby we ran each scheduler on different datasets. It has been
shown that this approach, while certainly useful, can be misleading in some
instances [14]. An adversarial approach to comparing algorithms was recently
proposed to address this. It may be interesting to evaluate the scheduling
algorithms and algorithmic components using this approach.

## Acknowledgements

This work was supported in part by Army Research Laboratory under Cooperative
Agreement W911NF-17-2-0196. The authors acknowledge the Center for Advanced
Research Computing (CARC) at the University of Southern California for
providing computing resources that have contributed to the research results
reported within this publication. URL: https://carc.usc.edu.

## References

  * [1] M. R. Garey and D. S. Johnson, _Computers and Intractability: A Guide to the Theory of NP-Completeness_. W. H. Freeman, 1979.
  * [2] A. Bazzi and A. Norouzi-Fard, “Towards tight lower bounds for scheduling problems,” in _Algorithms - ESA 2015 - 23rd Annual European Symposium, Patras, Greece, September 14-16, 2015, Proceedings_ , ser. Lecture Notes in Computer Science, N. Bansal and I. Finocchi, Eds., vol. 9294. Springer, 2015, pp. 118–129. [Online]. Available: https://doi.org/10.1007/978-3-662-48350-3_11
  * [3] A. Authors, “Scheduling algorithms gathered,” Github, 2023. [Online]. Available: https://anonymous.4open.science/r/saga-1F6D/README.md
  * [4] R. L. Graham, “Bounds on multiprocessing timing anomalies,” _SIAM Journal on Applied Mathematics_ , vol. 17, no. 2, pp. 416–429, 1969. [Online]. Available: https://doi.org/10.1137/0117039
  * [5] H. Topcuoglu, S. Hariri, and M. Wu, “Task scheduling algorithms for heterogeneous processors,” in _8th Heterogeneous Computing Workshop, HCW 1999, San Juan, Puerto Rico, April12, 1999_. IEEE Computer Society, 1999, pp. 3–14. [Online]. Available: https://doi.org/10.1109/HCW.1999.765092
  * [6] H. Wang and O. Sinnen, “List-scheduling versus cluster-scheduling,” _IEEE Trans. Parallel Distributed Syst._ , vol. 29, no. 8, pp. 1736–1749, 2018. [Online]. Available: https://doi.org/10.1109/TPDS.2018.2808959
  * [7] A. K. Maurya and A. K. Tripathi, “On benchmarking task scheduling algorithms for heterogeneous computing systems,” _J. Supercomput._ , vol. 74, no. 7, pp. 3039–3070, 2018. [Online]. Available: https://doi.org/10.1007/s11227-018-2355-0
  * [8] Y.-K. Kwok and I. Ahmad, “Benchmarking the task graph scheduling algorithms,” in _Proceedings of the First Merged International Parallel Processing Symposium and Symposium on Parallel and Distributed Processing_ , 1998, pp. 531–537.
  * [9] T. D. Braun, H. J. Siegel, N. Beck, L. Bölöni, M. Maheswaran, A. I. Reuther, J. P. Robertson, M. D. Theys, B. Yao, D. A. Hensgen, and R. F. Freund, “A comparison of eleven static heuristics for mapping a class of independent tasks onto heterogeneous distributed computing systems,” _J. Parallel Distributed Comput._ , vol. 61, no. 6, pp. 810–837, 2001. [Online]. Available: https://doi.org/10.1006/jpdc.2000.1714
  * [10] T. Braun, H. Siegal, N. Beck, L. Boloni, M. Maheswaran, A. Reuther, J. Robertson, M. Theys, B. Yao, D. Hensgen, and R. Freund, “A comparison study of static mapping heuristics for a class of meta-tasks on heterogeneous computing systems,” in _Proceedings. Eighth Heterogeneous Computing Workshop (HCW’99)_ , 1999, pp. 15–29.
  * [11] T. N’Takpé and F. Suter, “Critical path and area based scheduling of parallel task graphs on heterogeneous platforms,” in _12th International Conference on Parallel and Distributed Systems, ICPADS 2006, Minneapolis, Minnesota, USA, July 12-15, 2006_. IEEE Computer Society, 2006, pp. 3–10. [Online]. Available: https://doi.org/10.1109/ICPADS.2006.32
  * [12] D. Cordeiro, G. Mouniê, S. Perarnau, D. Trystram, J.-M. Vincent, and F. Wagner, “Random graph generation for scheduling simulations.” ICST, 5 2010.
  * [13] R. F. da Silva, R. Mayani, Y. Shi, A. R. Kemanian, M. Rynge, and E. Deelman, “Empowering agroecosystem modeling with HTC scientific workflows: The cycles model use case,” in _2019 IEEE International Conference on Big Data (IEEE BigData), Los Angeles, CA, USA, December 9-12, 2019_ , C. K. Baru, J. Huan, L. Khan, X. Hu, R. Ak, Y. Tian, R. S. Barga, C. Zaniolo, K. Lee, and Y. F. Ye, Eds. IEEE, 2019, pp. 4545–4552. [Online]. Available: https://doi.org/10.1109/BigData47090.2019.9006107
  * [14] J. Coleman and B. Krishnamachari, “Comparing task graph scheduling algorithms: An adversarial approach,” under review.